powerint.cpp
Go to the documentation of this file.
00001 /*
00002  *    This file is part of ACADO Toolkit.
00003  *
00004  *    ACADO Toolkit -- A Toolkit for Automatic Control and Dynamic Optimization.
00005  *    Copyright (C) 2008-2014 by Boris Houska, Hans Joachim Ferreau,
00006  *    Milan Vukov, Rien Quirynen, KU Leuven.
00007  *    Developed within the Optimization in Engineering Center (OPTEC)
00008  *    under supervision of Moritz Diehl. All rights reserved.
00009  *
00010  *    ACADO Toolkit is free software; you can redistribute it and/or
00011  *    modify it under the terms of the GNU Lesser General Public
00012  *    License as published by the Free Software Foundation; either
00013  *    version 3 of the License, or (at your option) any later version.
00014  *
00015  *    ACADO Toolkit is distributed in the hope that it will be useful,
00016  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
00017  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018  *    Lesser General Public License for more details.
00019  *
00020  *    You should have received a copy of the GNU Lesser General Public
00021  *    License along with ACADO Toolkit; if not, write to the Free Software
00022  *    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
00023  *
00024  */
00025 
00026 
00027 
00035 #include <acado/utils/acado_utils.hpp>
00036 #include <acado/symbolic_operator/symbolic_operator.hpp>
00037 
00038 
00039 
00040 BEGIN_NAMESPACE_ACADO
00041 
00042 
00043 
00044 
00045 Power_Int::Power_Int() : SmoothOperator( )
00046 {
00047     nCount = 0;
00048     derivative = 0;
00049     derivative2 = 0;
00050 }
00051 
00052 Power_Int::Power_Int( Operator *_argument, int _exponent ) : SmoothOperator( )
00053 {
00054     argument          = _argument                        ;
00055     exponent          = _exponent                        ;
00056     dargument         = NULL                             ;
00057     argument_result   = (double*)calloc(1,sizeof(double));
00058     dargument_result  = (double*)calloc(1,sizeof(double));
00059     bufferSize        = 1                                ;
00060     curvature         = CT_UNKNOWN                       ;
00061     monotonicity      = MT_UNKNOWN                       ;
00062 
00063     nCount = 0;
00064     derivative = 0;
00065     derivative2 = 0;
00066 }
00067 
00068 
00069 Power_Int::Power_Int( const Power_Int &arg ){
00070 
00071     int run1;
00072 
00073     bufferSize       = arg.bufferSize;
00074     exponent         = arg.exponent;
00075 
00076         argument         = arg.argument->clone();
00077         
00078 //      argument         = arg.argument;
00079 //      argument->nCount++;
00080 
00081     if( arg.dargument == NULL ){
00082         dargument = NULL;
00083     }
00084     else{
00085         dargument = arg.dargument->clone();
00086     }
00087 
00088     argument_result  = (double*)calloc(bufferSize,sizeof(double));
00089     dargument_result = (double*)calloc(bufferSize,sizeof(double));
00090 
00091     for( run1 = 0; run1 < bufferSize; run1++ ){
00092 
00093         argument_result[run1] = arg.argument_result[run1];
00094        dargument_result[run1] = arg.dargument_result[run1];
00095 
00096     }
00097     curvature         = arg.curvature   ;
00098     monotonicity      = arg.monotonicity;
00099 
00100     nCount = 0;
00101     derivative = 0;
00102     derivative2 = 0;
00103     if( arg.derivative != 0 ) {
00104         derivative = arg.derivative->clone();
00105         derivative2 = arg.derivative2->clone();
00106     }
00107 }
00108 
00109 
00110 Power_Int::~Power_Int(){
00111 
00112  delete argument;
00113 //     if( argument != 0 ){
00114 // 
00115 //         if( argument->nCount == 0 ){
00116 //             delete argument;
00117 //             argument = 0;
00118 //         }
00119 //         else{
00120 //             argument->nCount--;
00121 //         }
00122 //     }
00123 
00124     if( dargument != NULL ){
00125         delete dargument;
00126     }
00127 
00128     free(  argument_result );
00129     free( dargument_result );
00130 
00131     if( derivative != 0 ) {
00132         delete derivative;
00133     }
00134     if( derivative2 != 0 ) {
00135         delete derivative2;
00136     }
00137     derivative = 0;
00138     derivative2 = 0;
00139 }
00140 
00141 Power_Int& Power_Int::operator=( const Power_Int &arg ){
00142 
00143     if( this != &arg ){
00144 
00145          
00146          delete argument;
00147 //         if( argument != 0 ){
00148 // 
00149 //             if( argument->nCount == 0 ){
00150 //                 delete argument;
00151 //                 argument = 0;
00152 //             }
00153 //             else{
00154 //                 argument->nCount--;
00155 //             }
00156 //         }
00157 
00158         if( dargument != NULL ){
00159             delete dargument;
00160         }
00161 
00162         free(  argument_result );
00163         free( dargument_result );
00164 
00165         argument = arg.argument->clone();
00166                 
00167 //         argument = arg.argument;
00168 //              argument->nCount++;
00169 
00170         exponent          = arg.exponent                       ;
00171         dargument         = NULL                               ;
00172         bufferSize        = arg.bufferSize                     ;
00173         argument_result   = (double*)calloc(bufferSize,sizeof(double))  ;
00174         dargument_result  = (double*)calloc(bufferSize,sizeof(double))  ;
00175 
00176         curvature         = arg.curvature   ;
00177         monotonicity      = arg.monotonicity;
00178 
00179         nCount = 0;
00180     }
00181 
00182     return *this;
00183 }
00184 
00185 
00186 returnValue Power_Int::evaluate( int number, double *x, double *result ){
00187 
00188     if( number >= bufferSize ){
00189         bufferSize += number;
00190         argument_result  = (double*)realloc( argument_result,bufferSize*sizeof(double));
00191         dargument_result = (double*)realloc(dargument_result,bufferSize*sizeof(double));
00192     }
00193     argument->evaluate( number, x , &argument_result[number] );
00194 
00195     result[0] = pow( argument_result[number], exponent );
00196 
00197     return SUCCESSFUL_RETURN;
00198 }
00199 
00200 
00201 returnValue Power_Int::evaluate( EvaluationBase *x ){
00202  
00203     x->powerInt(*argument,exponent);
00204     return SUCCESSFUL_RETURN;
00205 }
00206 
00207 
00208 Operator* Power_Int::differentiate( int index ){
00209 
00210   if( exponent == 0 ){
00211     return new DoubleConstant( 0.0 , NE_ZERO );
00212   }
00213   
00214   dargument = argument->differentiate( index );
00215 
00216   Operator *result = myProd( derivative, dargument );
00217 
00218   return result;
00219 }
00220 
00221 
00222 Operator* Power_Int::AD_forward( int dim,
00223                                    VariableType *varType,
00224                                    int *component,
00225                                    Operator **seed,
00226                                    int &nNewIS,
00227                                    TreeProjection ***newIS ){
00228 
00229         if( dargument != 0 )
00230                 delete dargument;
00231 
00232         dargument = argument->AD_forward(dim,varType,component,seed,nNewIS,newIS);
00233 
00234         Operator *result = myProd( derivative, dargument );
00235 
00236         return result;
00237 }
00238 
00239 
00240 
00241 returnValue Power_Int::AD_backward( int           dim      , 
00242                                         VariableType *varType  , 
00243                                         int          *component, 
00244                                         Operator     *seed     , 
00245                                         Operator    **df       , 
00246                                         int           &nNewIS  , 
00247                                         TreeProjection ***newIS   ){
00248 
00249     Operator *result = myProd( derivative, seed );
00250 
00251     argument->AD_backward( dim, varType, component, result, df, nNewIS, newIS );
00252 
00253     delete seed;
00254     return SUCCESSFUL_RETURN;
00255 }
00256 
00257 
00258 
00259 returnValue Power_Int::AD_symmetric( int            dim       , 
00260                                         VariableType  *varType   , 
00261                                         int           *component , 
00262                                         Operator      *l         , 
00263                                         Operator     **S         , 
00264                                         int            dimS      , 
00265                                         Operator     **dfS       , 
00266                                         Operator     **ldf       , 
00267                                         Operator     **H         , 
00268                                       int            &nNewLIS  , 
00269                                       TreeProjection ***newLIS , 
00270                                       int            &nNewSIS  , 
00271                                       TreeProjection ***newSIS , 
00272                                       int            &nNewHIS  , 
00273                                       TreeProjection ***newHIS    ){
00274     
00275     TreeProjection tmp( *derivative );
00276     TreeProjection tmp2( *derivative2 );
00277 
00278     return ADsymCommon( argument, tmp, tmp2, dim, varType, component, l, S, dimS, dfS,
00279                          ldf, H, nNewLIS, newLIS, nNewSIS, newSIS, nNewHIS, newHIS );
00280 }
00281 
00282 
00283 returnValue Power_Int::initDerivative() {
00284 
00285         if( initialized ) return SUCCESSFUL_RETURN;
00286         initialized = BT_TRUE;
00287 
00288         Operator *powerTmp = myPowerInt( argument, exponent-1 );
00289         Operator *expTmp = new DoubleConstant( (double) exponent, NE_NEITHER_ONE_NOR_ZERO );
00290 
00291         derivative = convert2TreeProjection(myProd( expTmp, powerTmp ));
00292 
00293         Operator *powerTmp2 = myPowerInt( argument, exponent-2 );
00294         Operator *expTmp2 = new DoubleConstant( (double) exponent-1, NE_NEITHER_ONE_NOR_ZERO );
00295         Operator *prodTmp = myProd( expTmp, expTmp2 );
00296 
00297         derivative2 = convert2TreeProjection(myProd( prodTmp, powerTmp2 ));
00298 
00299         delete powerTmp;
00300         delete expTmp;
00301         delete powerTmp2;
00302         delete expTmp2;
00303         delete prodTmp;
00304 
00305         return argument->initDerivative();
00306 }
00307 
00308 
00309 
00310 Operator* Power_Int::substitute( int index, const Operator *sub ){
00311 
00312     return new Power_Int( argument->substitute( index , sub ), exponent );
00313 
00314 }
00315 
00316 
00317 NeutralElement Power_Int::isOneOrZero() const{ return NE_NEITHER_ONE_NOR_ZERO; }
00318 
00319 
00320 BooleanType Power_Int::isDependingOn( VariableType var ) const{
00321 
00322     return argument->isDependingOn(var);
00323 }
00324 
00325 
00326 BooleanType Power_Int::isDependingOn( int dim,
00327                                         VariableType *varType,
00328                                         int *component,
00329                                         BooleanType   *implicit_dep ){
00330 
00331     if( exponent == 0 ){
00332         return BT_FALSE;
00333     }
00334 
00335     return argument->isDependingOn( dim, varType, component, implicit_dep );
00336 }
00337 
00338 
00339 BooleanType Power_Int::isLinearIn( int dim,
00340                                      VariableType *varType,
00341                                      int *component,
00342                                      BooleanType   *implicit_dep ){
00343 
00344     if( exponent == 0 ){
00345         return BT_TRUE;
00346     }
00347 
00348     if( exponent == 1 && argument->isLinearIn( dim, varType, component, implicit_dep ) == BT_TRUE ){
00349         return BT_TRUE;
00350     }
00351 
00352     return BT_FALSE;
00353 }
00354 
00355 
00356 BooleanType Power_Int::isPolynomialIn( int dim,
00357                                          VariableType *varType,
00358                                          int *component,
00359                                          BooleanType   *implicit_dep ){
00360 
00361     if(  argument->isPolynomialIn( dim, varType, component, implicit_dep )    == BT_TRUE &&
00362          exponent >= 0 ){
00363         return BT_TRUE;
00364     }
00365 
00366     return BT_FALSE;
00367 }
00368 
00369 
00370 BooleanType Power_Int::isRationalIn( int dim,
00371                                        VariableType *varType,
00372                                        int *component,
00373                                        BooleanType   *implicit_dep ){
00374 
00375     if(  argument->isRationalIn( dim, varType, component, implicit_dep ) == BT_TRUE ){
00376         return BT_TRUE;
00377     }
00378 
00379     return BT_FALSE;
00380 }
00381 
00382 
00383 MonotonicityType Power_Int::getMonotonicity( ){
00384 
00385     if( monotonicity != MT_UNKNOWN )  return monotonicity;
00386 
00387     const MonotonicityType m = argument->getMonotonicity();
00388 
00389     if( m == MT_CONSTANT )  return MT_CONSTANT;
00390 
00391     // if exponent is even:
00392     if( fabs( ceil( ((double) exponent)/2.0 - EPS ) - ((double) exponent)/2.0 ) < 10.0*EPS ){
00393 
00394         if( exponent == 0 ) return MT_CONSTANT;
00395 
00396         return MT_NONMONOTONIC;
00397     }
00398     else{
00399 
00400         if( exponent > 0  )  return m;
00401         return MT_NONMONOTONIC;
00402     }
00403 
00404     return MT_NONMONOTONIC;
00405 }
00406 
00407 
00408 CurvatureType Power_Int::getCurvature( ){
00409 
00410     if( curvature != CT_UNKNOWN )  return curvature;
00411 
00412     const CurvatureType cc = argument->getCurvature();
00413 
00414     if( cc == CT_CONSTANT )  return CT_CONSTANT;
00415 
00416     // if exponent is even:
00417     if( fabs( ceil( ((double) exponent)/2.0 - EPS ) - ((double) exponent)/2.0 ) < 10.0*EPS ){
00418 
00419         if( exponent  < 0   ) return CT_NEITHER_CONVEX_NOR_CONCAVE;
00420         if( exponent == 0   ) return CT_CONSTANT;
00421         if( cc == CT_AFFINE ) return CT_CONVEX  ;
00422 
00423         return CT_NEITHER_CONVEX_NOR_CONCAVE;
00424     }
00425     else{
00426 
00427         if( exponent == 1 ) return cc;
00428         return CT_NEITHER_CONVEX_NOR_CONCAVE;
00429     }
00430 
00431     return CT_NEITHER_CONVEX_NOR_CONCAVE;
00432 }
00433 
00434 
00435 returnValue Power_Int::setMonotonicity( MonotonicityType monotonicity_ ){
00436 
00437     monotonicity = monotonicity_;
00438     return SUCCESSFUL_RETURN;
00439 }
00440 
00441 
00442 returnValue Power_Int::setCurvature( CurvatureType curvature_ ){
00443 
00444     curvature = curvature_;
00445     return SUCCESSFUL_RETURN;
00446 }
00447 
00448 
00449 returnValue Power_Int::AD_forward( int number, double *x, double *seed,
00450                                    double *f, double *df ){
00451 
00452     if( number >= bufferSize ){
00453         bufferSize += number;
00454         argument_result  = (double*)realloc( argument_result,bufferSize*sizeof(double));
00455         dargument_result = (double*)realloc(dargument_result,bufferSize*sizeof(double));
00456     }
00457     argument->AD_forward( number, x, seed, &argument_result[number],
00458                                   &dargument_result[number] );
00459 
00460       f[0] =  pow( argument_result[number],exponent );
00461      df[0] =  exponent*pow( argument_result[number],exponent-1 )*dargument_result[number];
00462 
00463      return SUCCESSFUL_RETURN;
00464 }
00465 
00466 
00467 
00468 returnValue Power_Int::AD_forward( int number, double *seed, double *df ){
00469 
00470 
00471     argument->AD_forward( number, seed, &dargument_result[number] );
00472 
00473      df[0] =  exponent*pow( argument_result[number],exponent-1 )*dargument_result[number];
00474 
00475      return SUCCESSFUL_RETURN;
00476 }
00477 
00478 
00479 returnValue Power_Int::AD_backward( int number, double seed, double *df ){
00480     return argument->AD_backward( number, exponent*pow( argument_result[number],
00481                                   exponent-1 )*seed, df );
00482 }
00483 
00484 
00485 returnValue Power_Int::AD_forward2( int number, double *seed, double *dseed,
00486                                     double *df, double *ddf ){
00487 
00488     double      ddargument_result;
00489     double      dargument_result2;
00490 
00491     argument->AD_forward2( number, seed, dseed, &dargument_result2, &ddargument_result);
00492 
00493     const double nn = exponent*pow( argument_result[number],exponent-1 );
00494      df[0] = nn*dargument_result2;
00495     ddf[0] = nn*ddargument_result + exponent*(exponent-1)*dargument_result[number]
00496                                   * dargument_result2*
00497                                     pow( argument_result[number],exponent-2 )   ;
00498 
00499     return SUCCESSFUL_RETURN;
00500 }
00501 
00502 
00503 returnValue Power_Int::AD_backward2( int number, double seed1, double seed2,
00504                                double *df, double *ddf ){
00505 
00506     const double nn = exponent*pow(argument_result[number],exponent-1);
00507 
00508     argument->AD_backward2( number,
00509                             seed1*nn,
00510                             seed2*nn +
00511                             seed1*exponent*(exponent-1)
00512                                  *pow(argument_result[number],exponent-2)
00513                                  *dargument_result[number],
00514                             df, ddf );
00515 
00516     return SUCCESSFUL_RETURN;
00517 }
00518 
00519 
00520 std::ostream& Power_Int::print( std::ostream &stream ) const{
00521 
00522         if ( argument->getName() == ON_POWER ) 
00523         {
00524                 Power* dummy = (Power*)argument;
00525 
00526                 if ( acadoIsEqual( dummy->argument2->getValue(),1.0/((double)exponent) ) == BT_TRUE )
00527                 {
00528                         return stream << "(" << *(dummy->argument1) << ")";
00529                 }
00530         }
00531 
00532         if ( exponent == 1 )
00533         {
00534             return stream << "(" << *argument << ")";
00535         }
00536         else
00537         {
00538                 VariableType dummy1;
00539                 int dummy2;
00540 
00541                 if ( ( exponent == 2 ) && ( argument->isVariable( dummy1,dummy2 ) == BT_TRUE ) )
00542                 {
00543                         return stream << "((" << *argument << ")*(" << *argument << "))";
00544                 }
00545                 else
00546                 {
00547                         return stream << "(pow(" << *argument << "," << exponent << "))";
00548                 }
00549         }
00550 }
00551 
00552 
00553 
00554 Operator* Power_Int::clone() const{
00555 
00556     return new Power_Int(*this);
00557 }
00558 
00559 
00560 returnValue Power_Int::clearBuffer(){
00561 
00562     if( bufferSize > 1 ){
00563         bufferSize = 1;
00564         argument_result  = (double*)realloc( argument_result,bufferSize*sizeof(double));
00565         dargument_result = (double*)realloc(dargument_result,bufferSize*sizeof(double));
00566     }
00567 
00568     return SUCCESSFUL_RETURN;
00569 }
00570 
00571 
00572 returnValue Power_Int::enumerateVariables( SymbolicIndexList *indexList ){
00573 
00574     return argument->enumerateVariables( indexList );
00575 }
00576 
00577 
00578 //
00579 // PROTECTED MEMBER FUNCTIONS:
00580 // ---------------------------
00581 
00582 
00583 OperatorName Power_Int::getName(){
00584 
00585     return ON_POWER_INT;
00586 }
00587 
00588 BooleanType Power_Int::isVariable( VariableType &varType, int &component ) const
00589 {
00590     return BT_FALSE;
00591 }
00592 
00593 returnValue Power_Int::loadIndices( SymbolicIndexList *indexList ){
00594 
00595     return argument->loadIndices( indexList );
00596 }
00597 
00598 
00599 BooleanType Power_Int::isSymbolic() const{
00600 
00601     if( argument->isSymbolic() == BT_FALSE ) return BT_FALSE;
00602     return BT_TRUE;
00603 }
00604 
00605 returnValue Power_Int::setVariableExportName(   const VariableType &_type,
00606                                                                                                 const std::vector< std::string >& _name
00607                                                                                                 )
00608 {
00609         argument->setVariableExportName(_type, _name);
00610 
00611         return Operator::setVariableExportName(_type, _name);
00612 }
00613 
00614 CLOSE_NAMESPACE_ACADO
00615 
00616 // end of file.


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Sat Jun 8 2019 19:38:42