power.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::Power( ):BinaryOperator( ){
00046 
00047         derivative01 = 0;
00048         derivative02 = 0;
00049         derivative12 = 0;
00050         derivative21 = 0;
00051         derivative22 = 0;
00052         derivative23 = 0;
00053 }
00054 
00055 Power::Power( Operator *_argument1, Operator *_argument2 )
00056       :BinaryOperator( _argument1, _argument2 ){
00057 
00058         derivative01 = 0;
00059         derivative02 = 0;
00060         derivative12 = 0;
00061         derivative21 = 0;
00062         derivative22 = 0;
00063         derivative23 = 0;
00064 }
00065 
00066 Power::Power( const Power &arg ):BinaryOperator( arg ){
00067 
00068         derivative01 = 0;
00069         derivative02 = 0;
00070         derivative12 = 0;
00071         derivative21 = 0;
00072         derivative22 = 0;
00073         derivative23 = 0;
00074         if( arg.derivative01 != 0 ) {
00075                 derivative01 = arg.derivative01->clone();
00076                 derivative02 = arg.derivative02->clone();
00077                 derivative12 = arg.derivative12->clone();
00078                 derivative21 = arg.derivative21->clone();
00079                 derivative22 = arg.derivative22->clone();
00080                 derivative23 = arg.derivative23->clone();
00081         }
00082 }
00083 
00084 
00085 Power::~Power(){
00086 
00087         if( derivative01 != 0 ) {
00088                 delete derivative01;
00089         }
00090         if( derivative02 != 0 ) {
00091                 delete derivative02;
00092         }
00093         if( derivative12 != 0 ) {
00094                 delete derivative12;
00095         }
00096         if( derivative21 != 0 ) {
00097                 delete derivative21;
00098         }
00099         if( derivative22 != 0 ) {
00100                 delete derivative22;
00101         }
00102         if( derivative23 != 0 ) {
00103                 delete derivative23;
00104         }
00105         derivative01 = 0;
00106         derivative02 = 0;
00107         derivative12 = 0;
00108         derivative21 = 0;
00109         derivative22 = 0;
00110         derivative23 = 0;
00111 }
00112 
00113 
00114 Power& Power::operator=( const Power &arg ){
00115 
00116     if( this != &arg ){
00117 
00118         BinaryOperator::operator=( arg );
00119     }
00120     return *this;
00121 }
00122 
00123 
00124 returnValue Power::evaluate( int number, double *x, double *result ){
00125 
00126     if( number >= bufferSize ){
00127         bufferSize += number;
00128         argument1_result  = (double*)realloc( argument1_result,bufferSize*sizeof(double));
00129         argument2_result  = (double*)realloc( argument2_result,bufferSize*sizeof(double));
00130         dargument1_result = (double*)realloc(dargument1_result,bufferSize*sizeof(double));
00131         dargument2_result = (double*)realloc(dargument2_result,bufferSize*sizeof(double));
00132     }
00133 
00134     argument1->evaluate( number, x , &argument1_result[number] );
00135     argument2->evaluate( number, x , &argument2_result[number] );
00136 
00137     result[0] = pow( argument1_result[number], argument2_result[number] );
00138 
00139     return SUCCESSFUL_RETURN;
00140 }
00141 
00142 
00143 returnValue Power::evaluate( EvaluationBase *x ){
00144 
00145     x->power(*argument1,*argument2);
00146     return SUCCESSFUL_RETURN;
00147 }
00148 
00149 
00150 Operator* Power::differentiate( int index ){
00151 
00152         dargument1 = argument1->differentiate( index );
00153         dargument2 = argument2->differentiate( index );
00154 
00155         Operator *prodTmp1 = myProd( this, derivative02);
00156         Operator *prodTmp2 = myProd( dargument2, prodTmp1 );
00157         Operator *prodTmp4 = myProd( dargument1, derivative12 );
00158 
00159         Operator *result = myAdd( prodTmp2, prodTmp4 );
00160 
00161         delete prodTmp1;
00162         delete prodTmp2;
00163         delete prodTmp4;
00164 
00165         return result;
00166 
00167 }
00168 
00169 
00170 Operator* Power::AD_forward( int dim,
00171                                VariableType *varType,
00172                                int *component,
00173                                Operator **seed,
00174                                int &nNewIS,
00175                                TreeProjection ***newIS ){
00176 
00177     if( dargument1 != 0 )
00178         delete dargument1;
00179 
00180     if( dargument2 != 0 )
00181         delete dargument2;
00182 
00183     dargument1 = argument1->AD_forward(dim,varType,component,seed,nNewIS,newIS);
00184     dargument2 = argument2->AD_forward(dim,varType,component,seed,nNewIS,newIS);
00185 
00186         Operator *prodTmp1 = myProd( this, derivative02);
00187     Operator *prodTmp2 = myProd( dargument2, prodTmp1 );
00188     Operator *prodTmp4 = myProd( dargument1, derivative12 );
00189 
00190     Operator *result = myAdd( prodTmp2, prodTmp4 );
00191 
00192     delete prodTmp1;
00193     delete prodTmp2;
00194     delete prodTmp4;
00195 
00196     return result;
00197 }
00198 
00199 
00200 returnValue Power::AD_backward( int           dim      , 
00201                                         VariableType *varType  , 
00202                                         int          *component, 
00203                                         Operator     *seed     , 
00204                                         Operator    **df       , 
00205                                         int           &nNewIS  , 
00206                                         TreeProjection ***newIS   ){
00207 
00208         if( seed->isOneOrZero() != NE_ZERO ){
00209 
00210                 TreeProjection tmp;
00211                 tmp = *seed;
00212 
00213                 Operator *prodTmp4 = myProd( &tmp, derivative12 );
00214 
00215                 argument1->AD_backward( dim, varType, component, prodTmp4->clone(), df, nNewIS, newIS );
00216 
00217                 Operator *prodTmp1 = myProd( this, derivative02);
00218                 Operator *prodTmp2 = myProd( &tmp, prodTmp1 );
00219 
00220                 argument2->AD_backward( dim, varType, component, prodTmp2->clone(), df, nNewIS, newIS );
00221 
00222                 delete prodTmp1;
00223                 delete prodTmp2;
00224                 delete prodTmp4;
00225         }
00226 
00227     delete seed;
00228     return SUCCESSFUL_RETURN;
00229 }
00230 
00231 
00232 returnValue Power::AD_symmetric( int            dim       , 
00233                                         VariableType  *varType   , 
00234                                         int           *component , 
00235                                         Operator      *l         , 
00236                                         Operator     **S         , 
00237                                         int            dimS      , 
00238                                         Operator     **dfS       , 
00239                                         Operator     **ldf       , 
00240                                         Operator     **H         , 
00241                                       int            &nNewLIS  , 
00242                                       TreeProjection ***newLIS , 
00243                                       int            &nNewSIS  , 
00244                                       TreeProjection ***newSIS , 
00245                                       int            &nNewHIS  , 
00246                                       TreeProjection ***newHIS    ){
00247 
00248     TreeProjection dyy, dy;
00249     
00250     TreeProjection dx( *derivative12 );
00251     dy = Product( clone(), derivative02->clone());
00252     TreeProjection dxx( *derivative22 );
00253     TreeProjection dxy( *derivative23 );
00254     dyy = Product( dy.clone(), derivative02->clone() );
00255     
00256     return ADsymCommon2( argument1,argument2,dx,dy,dxx,dxy,dyy, dim, varType, component, l, S, dimS, dfS,
00257                           ldf, H, nNewLIS, newLIS, nNewSIS, newSIS, nNewHIS, newHIS );
00258 }
00259 
00260 
00261 returnValue Power::initDerivative() {
00262 
00263         if( initialized ) return SUCCESSFUL_RETURN;
00264         initialized = BT_TRUE;
00265 
00266         Operator *oneTmp = new DoubleConstant(1.0, NE_ONE);
00267         Operator *subTmp = mySubtract( argument2, oneTmp );
00268 
00269         derivative01 = convert2TreeProjection(myPower( argument1, subTmp));
00270         derivative02 = convert2TreeProjection(myLogarithm( argument1 ));
00271 
00272         derivative12 = convert2TreeProjection(myProd( derivative01, argument2 ));
00273 
00274         Operator *twoTmp = new DoubleConstant(2.0,NE_NEITHER_ONE_NOR_ZERO);
00275         Operator *subTmp2 = mySubtract( argument2, twoTmp );
00276         Operator *prodTmp = myProd( argument2, subTmp );
00277         Operator *prodTmp2 = myProd( argument2, derivative02 );
00278         Operator *addTmp = myAdd( oneTmp, prodTmp2 );
00279 
00280         derivative21 = convert2TreeProjection(myPower( argument1, subTmp2));
00281         derivative22 = convert2TreeProjection(myProd( prodTmp, derivative21 ));
00282         derivative23 = convert2TreeProjection(myProd( derivative01, addTmp ));
00283 
00284         delete oneTmp;
00285         delete subTmp;
00286         delete twoTmp;
00287         delete subTmp2;
00288         delete prodTmp;
00289         delete prodTmp2;
00290         delete addTmp;
00291 
00292         argument1->initDerivative();
00293         return argument2->initDerivative();
00294 }
00295 
00296 
00297 Operator* Power::substitute( int index, const Operator *sub ){
00298 
00299     return new Power( argument1->substitute( index , sub ),
00300                       argument2->substitute( index , sub ) );
00301 
00302 }
00303 
00304 
00305 BooleanType Power::isLinearIn( int dim,
00306                                  VariableType *varType,
00307                                  int *component,
00308                                  BooleanType   *implicit_dep ){
00309 
00310     if( argument2->isOneOrZero() == NE_ZERO ){
00311         return BT_TRUE;
00312     }
00313 
00314     if( argument1->isOneOrZero() == NE_ONE ){
00315         return BT_TRUE;
00316     }
00317 
00318     if( argument1->isOneOrZero() == NE_ZERO ){
00319         return BT_TRUE;
00320     }
00321 
00322     if( argument1->isLinearIn( dim, varType, component, implicit_dep ) == BT_TRUE &&
00323         argument2->isOneOrZero() == NE_ONE ){
00324         return BT_TRUE;
00325     }
00326 
00327     if( isDependingOn( dim, varType, component, implicit_dep) == BT_FALSE ){
00328         return BT_TRUE;
00329     }
00330 
00331     return BT_FALSE;
00332 }
00333 
00334 
00335 BooleanType Power::isPolynomialIn( int dim,
00336                                      VariableType *varType,
00337                                      int *component,
00338                                      BooleanType   *implicit_dep ){
00339 
00340     if(  argument1->isPolynomialIn( dim, varType, component, implicit_dep ) == BT_TRUE  &&
00341          argument2->isDependingOn( dim, varType, component, implicit_dep )  == BT_FALSE ){
00342         return BT_TRUE;
00343     }
00344 
00345     return BT_FALSE;
00346 }
00347 
00348 
00349 BooleanType Power::isRationalIn( int dim,
00350                                    VariableType *varType,
00351                                    int *component,
00352                                    BooleanType   *implicit_dep ){
00353 
00354     if(  argument1->isRationalIn( dim, varType, component, implicit_dep )   == BT_TRUE  &&
00355          argument2->isDependingOn( dim, varType, component, implicit_dep )  == BT_FALSE ){
00356         return BT_TRUE;
00357     }
00358 
00359     return BT_FALSE;
00360 }
00361 
00362 
00363 MonotonicityType Power::getMonotonicity( ){
00364 
00365     if( monotonicity != MT_UNKNOWN )  return monotonicity;
00366 
00367     const MonotonicityType m1 = argument1->getMonotonicity();
00368     const MonotonicityType m2 = argument2->getMonotonicity();
00369 
00370     if( m1 == MT_CONSTANT ){
00371 
00372         if( m2 == MT_CONSTANT      )  return MT_CONSTANT;
00373 
00374         double res;
00375         argument1->evaluate(0,0,&res);
00376 
00377         if( res >= 1.0 ) return m2;
00378 
00379         if( m2 == MT_NONDECREASING )  return MT_NONINCREASING;
00380         if( m2 == MT_NONINCREASING )  return MT_NONDECREASING;
00381     }
00382 
00383 
00384     if( m2 == MT_CONSTANT ){
00385 
00386         double res;
00387         argument2->evaluate(0,0,&res);
00388 
00389         if( res >= 0.0 ) return m1;
00390 
00391         if( res < 0.0 ){
00392 
00393             if( m1 == MT_NONDECREASING )  return MT_NONINCREASING;
00394             if( m1 == MT_NONINCREASING )  return MT_NONDECREASING;
00395         }
00396     }
00397 
00398     return MT_NONMONOTONIC;
00399 }
00400 
00401 
00402 CurvatureType Power::getCurvature( ){
00403 
00404     if( curvature != CT_UNKNOWN )  return curvature;
00405 
00406     const CurvatureType c1 = argument1->getCurvature();
00407     const CurvatureType c2 = argument2->getCurvature();
00408 
00409     if( c1 == CT_CONSTANT ){
00410 
00411         if( c2 == CT_CONSTANT )  return CT_CONSTANT;
00412         if( c2 == CT_AFFINE   )  return CT_CONVEX  ;
00413 
00414         if( c2 == CT_CONVEX  ){
00415 
00416             double res;
00417             argument1->evaluate(0,0,&res);
00418 
00419             if( res >= 1.0 )  return CT_CONVEX;
00420         }
00421         if( c2 == CT_CONCAVE  ){
00422 
00423             double res;
00424             argument1->evaluate(0,0,&res);
00425 
00426             if( res < 1.0 )  return CT_CONVEX;
00427         }
00428 
00429         return CT_NEITHER_CONVEX_NOR_CONCAVE;
00430     }
00431 
00432     if( c2 == CT_CONSTANT ){
00433 
00434         double res;
00435         argument2->evaluate(0,0,&res);
00436 
00437         if( c1 == CT_AFFINE ){
00438 
00439             if( res >= 1.0 || res <  0.0 ) return CT_CONVEX ;
00440             if( res <  1.0 && res >= 0.0 ) return CT_CONCAVE;
00441         }
00442 
00443         if( c1 == CT_CONVEX ){
00444 
00445             if( res >= 1.0 ) return CT_CONVEX;
00446         }
00447 
00448         if( c1 == CT_CONCAVE ){
00449 
00450             if( res <  0.0               ) return CT_CONVEX ;
00451             if( res <= 1.0 && res >= 0.0 ) return CT_CONCAVE;
00452         }
00453     }
00454 
00455     return CT_NEITHER_CONVEX_NOR_CONCAVE;
00456 }
00457 
00458 
00459 returnValue Power::AD_forward( int number, double *x, double *seed,
00460                                double *f, double *df ){
00461 
00462     if( number >= bufferSize ){
00463         bufferSize += number;
00464         argument1_result  = (double*)realloc( argument1_result,bufferSize*sizeof(double));
00465         argument2_result  = (double*)realloc( argument2_result,bufferSize*sizeof(double));
00466         dargument1_result = (double*)realloc(dargument1_result,bufferSize*sizeof(double));
00467         dargument2_result = (double*)realloc(dargument2_result,bufferSize*sizeof(double));
00468     }
00469 
00470     argument1->AD_forward( number, x, seed, &argument1_result[number],
00471                            &dargument1_result[number] );
00472     argument2->AD_forward( number, x, seed, &argument2_result[number],
00473                            &dargument2_result[number] );
00474 
00475       f[0] = pow(argument1_result[number],argument2_result[number]);
00476      df[0] = argument2_result[number]*pow(argument1_result[number],argument2_result[number]
00477              -1.0)*dargument1_result[number]
00478              +f[0]*log(argument1_result[number])*dargument2_result[number];
00479 
00480      return SUCCESSFUL_RETURN;
00481 }
00482 
00483 
00484 returnValue Power::AD_forward( int number, double *seed, double *df ){
00485 
00486     argument1->AD_forward( number, seed, &dargument1_result[number] );
00487     argument2->AD_forward( number, seed, &dargument2_result[number] );
00488 
00489      df[0] = argument2_result[number]*pow(argument1_result[number],argument2_result[number]
00490              -1.0)*dargument1_result[number]
00491              +pow(argument1_result[number],argument2_result[number])
00492              *log(argument1_result[number])*dargument2_result[number];
00493 
00494      return SUCCESSFUL_RETURN;
00495 }
00496 
00497 
00498 returnValue Power::AD_backward( int number, double seed, double *df ){
00499 
00500     argument1->AD_backward(number, argument2_result[number]*pow(argument1_result[number],
00501                            argument2_result[number]-1.0)*seed, df );
00502 
00503     argument2->AD_backward(number, pow(argument1_result[number],argument2_result[number])
00504                            *log(argument1_result[number])*seed, df );
00505 
00506     return SUCCESSFUL_RETURN;
00507 }
00508 
00509 
00510 
00511 returnValue Power::AD_forward2( int number, double *seed, double *dseed,
00512                                 double *df, double *ddf ){
00513 
00514     double      ddargument1_result;
00515     double      ddargument2_result;
00516     double      dargument_result1 ;
00517     double      dargument_result2 ;
00518 
00519     argument1->AD_forward2( number, seed, dseed,
00520                             &dargument_result1, &ddargument1_result);
00521     argument2->AD_forward2( number, seed, dseed,
00522                             &dargument_result2, &ddargument2_result);
00523 
00524     const double nn1 = pow( argument1_result[number], argument2_result[number]     );
00525     const double nn2 = pow( argument1_result[number], argument2_result[number]-1.0 );
00526     const double nn3 = log( argument1_result[number] );
00527     const double nn4 = nn2*argument2_result[number];
00528     const double nn5 = nn1*nn3;
00529     const double nn6 = nn2*(argument2_result[number]*nn3 + 1.0);
00530 
00531      df[0] = nn4*dargument_result1
00532             +nn5*dargument_result2;
00533     ddf[0] = nn4*ddargument1_result
00534             +nn5*ddargument2_result
00535             +argument2_result[number]*(argument2_result[number]-1.0)
00536                  *pow(argument1_result[number],argument2_result[number]-2.0)
00537                  *dargument1_result[number]*dargument_result1
00538             +nn5*nn3*dargument2_result[number]*dargument_result2
00539             +nn6*dargument_result1*dargument2_result[number]
00540             +nn6*dargument1_result[number]*dargument_result2;
00541 
00542     return SUCCESSFUL_RETURN;
00543 }
00544 
00545 returnValue Power::AD_backward2( int number, double seed1, double seed2,
00546                                        double *df, double *ddf ){
00547 
00548     const double nn1 = pow( argument1_result[number], argument2_result[number]     );
00549     const double nn2 = pow( argument1_result[number], argument2_result[number]-1.0 );
00550     const double nn3 = log( argument1_result[number] );
00551     const double nn4 = nn2*argument2_result[number];
00552     const double nn5 = nn1*nn3;
00553     const double nn6 = nn2*(argument2_result[number]*nn3 + 1.0);
00554 
00555     argument1->AD_backward2(  number,
00556                               seed1*nn4,
00557                               seed2*nn4 + seed1*(
00558                                   argument2_result[number]*(argument2_result[number]-1.0)*
00559                                   pow( argument1_result[number],
00560                                        argument2_result[number]-2.0 )*
00561                                   dargument1_result[number]
00562                                 + nn6*dargument2_result[number] ),
00563                               df, ddf );
00564 
00565     argument2->AD_backward2(  number,
00566                               seed1*nn5,
00567                               seed2*nn5 + seed1*(
00568                                   nn5*nn3*dargument2_result[number]
00569                                 + nn6*dargument1_result[number]   ),
00570                               df, ddf );
00571 
00572     return SUCCESSFUL_RETURN;
00573 }
00574 
00575 
00576 std::ostream& Power::print( std::ostream &stream ) const{
00577 
00578         if ( acadoIsEqual( argument2->getValue(),0.5 ) == BT_TRUE )
00579         {
00580                 return stream << "(sqrt(" << *argument1 << "))";
00581         }
00582         else
00583         {
00584                 if ( acadoIsEqual( argument2->getValue(),-0.5 ) == BT_TRUE )
00585                 {
00586                         return stream << "(1.0/sqrt(" << *argument1 << "))";
00587                 }
00588                 else
00589                 {
00590                         return stream << "(pow(" << *argument1 << "," << *argument2 << "))";
00591                 }
00592         }
00593 }
00594 
00595 
00596 Operator* Power::clone() const{
00597 
00598     return new Power(*this);
00599 }
00600 
00601 
00602 OperatorName Power::getName(){
00603 
00604     return ON_POWER;
00605 
00606 }
00607 
00608 
00609 CLOSE_NAMESPACE_ACADO
00610 
00611 // end of file.


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