quotient.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 Quotient::Quotient():BinaryOperator(){
00046         derivative0 = 0;
00047         derivative1 = 0;
00048         derivative2 = 0;
00049 }
00050 
00051 Quotient::Quotient( Operator *_argument1, Operator *_argument2 )
00052          :BinaryOperator( _argument1, _argument2 ){
00053 
00054         derivative0 = 0;
00055         derivative1 = 0;
00056         derivative2 = 0;
00057 }
00058 
00059 
00060 Quotient::Quotient( const Quotient &arg ):BinaryOperator( arg ){
00061 
00062         derivative0 = 0;
00063         derivative1 = 0;
00064         derivative2 = 0;
00065         if( arg.derivative0 != 0 && arg.derivative1 != 0 && arg.derivative2 != 0 ) {
00066                 derivative0 = arg.derivative0->clone();
00067                 derivative1 = arg.derivative1->clone();
00068                 derivative2 = arg.derivative2->clone();
00069         }
00070 }
00071 
00072 
00073 Quotient::~Quotient(){
00074 
00075         if( derivative2 != 0 ) {
00076                 delete derivative2;
00077         }
00078         if( derivative1 != 0 ) {
00079                 delete derivative1;
00080         }
00081         if( derivative0 != 0 ) {
00082                 delete derivative0;
00083         }
00084         derivative0 = 0;
00085         derivative1 = 0;
00086         derivative2 = 0;
00087 }
00088 
00089 Quotient& Quotient::operator=( const Quotient &arg ){
00090 
00091     if( this != &arg ){
00092 
00093         BinaryOperator::operator=( arg );
00094     }
00095     return *this;
00096 }
00097 
00098 
00099 
00100 returnValue Quotient::evaluate( int number, double *x, double *result ){
00101 
00102     if( number >= bufferSize ){
00103         bufferSize += number;
00104         argument1_result  = (double*)realloc( argument1_result,bufferSize*sizeof(double));
00105         argument2_result  = (double*)realloc( argument2_result,bufferSize*sizeof(double));
00106         dargument1_result = (double*)realloc(dargument1_result,bufferSize*sizeof(double));
00107         dargument2_result = (double*)realloc(dargument2_result,bufferSize*sizeof(double));
00108     }
00109 
00110     argument1->evaluate( number, x , &argument1_result[number] );
00111     argument2->evaluate( number, x , &argument2_result[number] );
00112 
00113     result[0] = argument1_result[number] / argument2_result[number];
00114 
00115     return SUCCESSFUL_RETURN;
00116 }
00117 
00118 
00119 returnValue Quotient::evaluate( EvaluationBase *x ){
00120 
00121     x->quotient(*argument1,*argument2);
00122     return SUCCESSFUL_RETURN;
00123 }
00124 
00125 
00126 Operator* Quotient::differentiate( int index ){
00127 
00128         dargument1 = argument1->differentiate( index );
00129         dargument2 = argument2->differentiate( index );
00130 
00131         Operator* prodTmp = myProd( dargument1, derivative0 );
00132         Operator* prodTmp2 = myProd( argument1, dargument2 );
00133         Operator* prodTmp3 = myProd( prodTmp2, derivative1 );
00134         Operator* result = mySubtract( prodTmp, prodTmp3 );
00135 
00136         delete prodTmp;
00137         delete prodTmp2;
00138         delete prodTmp3;
00139 
00140         return result;
00141 }
00142 
00143 
00144 
00145 
00146 Operator* Quotient::AD_forward( int dim,
00147                                   VariableType *varType,
00148                                   int *component,
00149                                   Operator **seed,
00150                                   int &nNewIS,
00151                                   TreeProjection ***newIS ){
00152 
00153     if( dargument1 != 0 )
00154         delete dargument1;
00155 
00156     if( dargument2 != 0 )
00157         delete dargument2;
00158 
00159     dargument1 = argument1->AD_forward(dim,varType,component,seed,nNewIS,newIS);
00160     dargument2 = argument2->AD_forward(dim,varType,component,seed,nNewIS,newIS);
00161 
00162     Operator* prodTmp = myProd( dargument1, derivative0 );
00163     Operator* prodTmp2 = myProd( argument1, dargument2 );
00164     Operator* prodTmp3 = myProd( prodTmp2, derivative1 );
00165     Operator* result = mySubtract( prodTmp, prodTmp3 );
00166 
00167     delete prodTmp;
00168     delete prodTmp2;
00169     delete prodTmp3;
00170 
00171     return result;
00172 }
00173 
00174 
00175 returnValue Quotient::AD_backward( int           dim      , 
00176                                         VariableType *varType  , 
00177                                         int          *component, 
00178                                         Operator     *seed     , 
00179                                         Operator    **df       , 
00180                                         int           &nNewIS  , 
00181                                         TreeProjection ***newIS   ){
00182 
00183         if( seed->isOneOrZero() != NE_ZERO ){
00184 
00185                 TreeProjection tmp;
00186                 tmp = *seed;
00187 
00188                 Operator *prodTmp = myProd( &tmp, derivative0 );
00189 
00190                 argument1->AD_backward( dim, varType, component, prodTmp->clone(), df, nNewIS, newIS );
00191 
00192                 Operator *prodTmp2 = myProd( argument1, &tmp );
00193                 Operator *prodTmp3 = myProd( prodTmp2, derivative1 );
00194                 Operator *zeroTmp = new DoubleConstant( 0.0, NE_ZERO );
00195                 Operator *subTmp = mySubtract( zeroTmp, prodTmp3 );
00196 
00197                 argument2->AD_backward( dim, varType, component, subTmp->clone(), df, nNewIS, newIS );
00198 
00199                 delete zeroTmp;
00200                 delete prodTmp;
00201                 delete prodTmp2;
00202                 delete prodTmp3;
00203                 delete subTmp;
00204         }
00205 
00206         delete seed;
00207         return SUCCESSFUL_RETURN;
00208 }
00209 
00210 
00211 returnValue Quotient::AD_symmetric( int            dim       , 
00212                                         VariableType  *varType   , 
00213                                         int           *component , 
00214                                         Operator      *l         , 
00215                                         Operator     **S         , 
00216                                         int            dimS      , 
00217                                         Operator     **dfS       , 
00218                                         Operator     **ldf       , 
00219                                         Operator     **H         , 
00220                                       int            &nNewLIS  , 
00221                                       TreeProjection ***newLIS , 
00222                                       int            &nNewSIS  , 
00223                                       TreeProjection ***newSIS , 
00224                                       int            &nNewHIS  , 
00225                                       TreeProjection ***newHIS    ){
00226 
00227     TreeProjection dy,dxx,dxy,dyy;
00228     
00229     TreeProjection dx(*derivative0);
00230     dxy = Product( new DoubleConstant(-1.0,NE_NEITHER_ONE_NOR_ZERO), derivative1->clone() );
00231     dxx = DoubleConstant(0.0,NE_ZERO);
00232     dy  = Product( dxy.clone(), argument1->clone() );
00233     dyy = Product( new Product( new DoubleConstant(2.0,NE_NEITHER_ONE_NOR_ZERO), derivative2->clone() ),
00234                     argument1->clone() );
00235     
00236     return ADsymCommon2( argument1,argument2,dx,dy,dxx,dxy,dyy, dim, varType, component, l, S, dimS, dfS,
00237                           ldf, H, nNewLIS, newLIS, nNewSIS, newSIS, nNewHIS, newHIS );
00238 }
00239 
00240 
00241 returnValue Quotient::initDerivative() {
00242 
00243         if( initialized ) return SUCCESSFUL_RETURN;
00244         initialized = BT_TRUE;
00245 
00246         derivative0 = convert2TreeProjection(new Quotient( new DoubleConstant(1.0,NE_ONE), argument2->clone() ));
00247         derivative1 = convert2TreeProjection(new Product( derivative0->clone(), derivative0->clone() ));
00248         derivative2 = convert2TreeProjection(new Product( derivative0->clone(), derivative1->clone() ));
00249 
00250         argument1->initDerivative();
00251         return argument2->initDerivative();
00252 }
00253 
00254 
00255 Operator* Quotient::substitute( int index, const Operator *sub ){
00256 
00257     return new Quotient( argument1->substitute( index , sub ),
00258                          argument2->substitute( index , sub ) );
00259 
00260 }
00261 
00262 
00263 BooleanType Quotient::isLinearIn( int dim,
00264                                     VariableType *varType,
00265                                     int *component,
00266                                     BooleanType   *implicit_dep ){
00267 
00268     if(  argument1->isLinearIn( dim, varType, component, implicit_dep )    == BT_TRUE &&
00269          argument2->isDependingOn( dim, varType, component, implicit_dep ) == BT_FALSE ){
00270         return BT_TRUE;
00271     }
00272 
00273     return BT_FALSE;
00274 }
00275 
00276 
00277 BooleanType Quotient::isPolynomialIn( int dim,
00278                                         VariableType *varType,
00279                                         int *component,
00280                                         BooleanType   *implicit_dep ){
00281 
00282     if(  argument1->isPolynomialIn( dim, varType, component, implicit_dep )    == BT_TRUE  &&
00283          argument2->isDependingOn( dim, varType, component, implicit_dep )     == BT_FALSE ){
00284         return BT_TRUE;
00285     }
00286 
00287     return BT_FALSE;
00288 }
00289 
00290 
00291 BooleanType Quotient::isRationalIn( int dim,
00292                                       VariableType *varType,
00293                                       int *component,
00294                                       BooleanType   *implicit_dep ){
00295 
00296     if(  argument1->isRationalIn( dim, varType, component, implicit_dep )    == BT_TRUE  &&
00297          argument2->isRationalIn( dim, varType, component, implicit_dep )    == BT_TRUE ){
00298         return BT_TRUE;
00299     }
00300 
00301     return BT_FALSE;
00302 }
00303 
00304 
00305 MonotonicityType Quotient::getMonotonicity( ){
00306 
00307     if( monotonicity != MT_UNKNOWN )  return monotonicity;
00308 
00309     MonotonicityType m1, m2;
00310 
00311     m1 = argument1->getMonotonicity();
00312     m2 = argument2->getMonotonicity();
00313 
00314     if( m2 == MT_CONSTANT ){
00315 
00316         if( m1 == MT_CONSTANT )  return MT_CONSTANT;
00317 
00318         double res;
00319         argument2->evaluate(0,0,&res);
00320 
00321         if( res >= 0.0 ) return m1;
00322 
00323         if( m1 == MT_NONDECREASING ) return MT_NONINCREASING;
00324         if( m1 == MT_NONINCREASING ) return MT_NONDECREASING;
00325 
00326         return MT_NONMONOTONIC;
00327     }
00328 
00329     return MT_NONMONOTONIC;
00330 }
00331 
00332 
00333 CurvatureType Quotient::getCurvature( ){
00334 
00335     if( curvature != CT_UNKNOWN )  return curvature;
00336 
00337     CurvatureType c1, c2;
00338 
00339     c1 = argument1->getCurvature();
00340     c2 = argument2->getCurvature();
00341 
00342     if( c2 == CT_CONSTANT ){
00343 
00344         if( c1 == CT_CONSTANT )  return CT_CONSTANT;
00345 
00346         double res;
00347         argument2->evaluate(0,0,&res);
00348 
00349         if( res >= 0.0 ) return c1;
00350 
00351         if( c1 == CT_AFFINE  ) return CT_AFFINE ;
00352         if( c1 == CT_CONVEX  ) return CT_CONCAVE;
00353         if( c1 == CT_CONCAVE ) return CT_CONVEX ;
00354 
00355         return CT_NEITHER_CONVEX_NOR_CONCAVE;
00356     }
00357 
00358     return CT_NEITHER_CONVEX_NOR_CONCAVE;
00359 }
00360 
00361 
00362 double Quotient::getValue() const
00363 { 
00364         if ( ( argument1 == 0 ) || ( argument2 == 0 ) )
00365                 return INFTY;
00366                 
00367         if ( ( acadoIsEqual( argument1->getValue(),INFTY ) == BT_TRUE ) ||
00368                  ( acadoIsEqual( argument2->getValue(),INFTY ) == BT_TRUE ) )
00369                 return INFTY;
00370 
00371         if (acadoIsZero( argument2->getValue() ) == BT_TRUE)
00372                 ACADOFATAL(RET_DIV_BY_ZERO);
00373 
00374         return (argument1->getValue() / argument2->getValue());
00375 }
00376 
00377 returnValue Quotient::AD_forward( int number, double *x, double *seed,
00378                                  double *f, double *df ){
00379 
00380     if( number >= bufferSize ){
00381         bufferSize += number;
00382         argument1_result  = (double*)realloc( argument1_result,bufferSize*sizeof(double));
00383         argument2_result  = (double*)realloc( argument2_result,bufferSize*sizeof(double));
00384         dargument1_result = (double*)realloc(dargument1_result,bufferSize*sizeof(double));
00385         dargument2_result = (double*)realloc(dargument2_result,bufferSize*sizeof(double));
00386     }
00387 
00388     argument1->AD_forward( number, x, seed, &argument1_result[number],
00389                            &dargument1_result[number] );
00390     argument2->AD_forward( number,
00391                            x, seed, &argument2_result[number], &dargument2_result[number] );
00392 
00393       f[0] =  argument1_result[number]/argument2_result[number];
00394      df[0] =  dargument1_result[number]/argument2_result[number]
00395              -(argument1_result[number]*dargument2_result[number])/
00396               (argument2_result[number]*argument2_result[number] );
00397 
00398      return SUCCESSFUL_RETURN;
00399 }
00400 
00401 
00402 
00403 returnValue Quotient::AD_forward( int number, double *seed, double *df ){
00404 
00405     argument1->AD_forward( number, seed, &dargument1_result[number] );
00406     argument2->AD_forward( number, seed, &dargument2_result[number] );
00407 
00408      df[0] =  dargument1_result[number]/argument2_result[number]
00409              -(argument1_result[number]*dargument2_result[number])/
00410               (argument2_result[number]*argument2_result[number] );
00411 
00412      return SUCCESSFUL_RETURN;
00413 }
00414 
00415 
00416 returnValue Quotient::AD_backward( int number, double seed, double *df ){
00417 
00418     argument1->AD_backward( number, seed/argument2_result[number], df );
00419     argument2->AD_backward( number, -argument1_result[number]*seed/
00420                             (argument2_result[number]*argument2_result[number]), df );
00421 
00422     return SUCCESSFUL_RETURN;
00423 }
00424 
00425 
00426 returnValue Quotient::AD_forward2( int number, double *seed, double *dseed,
00427                                    double *df, double *ddf ){
00428 
00429     double      ddargument1_result;
00430     double      ddargument2_result;
00431     double      dargument_result1;
00432     double      dargument_result2;
00433 
00434     argument1->AD_forward2( number, seed, dseed,
00435                             &dargument_result1, &ddargument1_result);
00436     argument2->AD_forward2( number, seed, dseed,
00437                             &dargument_result2, &ddargument2_result);
00438 
00439     const double gg  =   argument2_result[number]*argument2_result[number];
00440     const double ggg =   dargument_result2/gg;
00441 
00442      df[0] =   dargument_result1/argument2_result[number]
00443               -argument1_result[number]*ggg;
00444 
00445     ddf[0] =   ddargument1_result/argument2_result[number]
00446               -argument1_result[number]*ddargument2_result/gg
00447               -dargument_result2*dargument1_result[number]/gg
00448               -dargument_result1*dargument2_result[number]/gg
00449               +2.0*argument1_result[number]/(gg*argument2_result[number])
00450               *dargument_result2*dargument2_result[number];
00451 
00452     return SUCCESSFUL_RETURN;
00453 }
00454 
00455 
00456 returnValue Quotient::AD_backward2( int number, double seed1, double seed2,
00457                                        double *df, double *ddf ){
00458 
00459     const double gg  =   argument2_result[number]*argument2_result[number];
00460     const double ggg =   argument1_result[number]/gg;
00461 
00462     argument1->AD_backward2(  number, seed1/argument2_result[number],
00463                                       seed2/argument2_result[number] -
00464                                       seed1*dargument2_result[number]/gg, df, ddf );
00465 
00466     argument2->AD_backward2( number, -seed1*ggg,
00467                                      -seed2*ggg
00468                                      -seed1*dargument1_result[number]/gg
00469                                      +2.0*ggg/argument2_result[number]
00470                                       *seed1*dargument2_result[number],
00471                              df, ddf );
00472 
00473     return SUCCESSFUL_RETURN;
00474 }
00475 
00476 
00477 std::ostream& Quotient::print( std::ostream &stream ) const{
00478 
00479     return stream << "(" << *argument1 << "/" << *argument2 << ")";
00480 }
00481 
00482 
00483 Operator* Quotient::clone() const{
00484 
00485     return new Quotient(*this);
00486 }
00487 
00488 
00489 OperatorName Quotient::getName(){
00490 
00491     return ON_QUOTIENT;
00492 }
00493 
00494 
00495 CLOSE_NAMESPACE_ACADO
00496 
00497 // end of file.


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