00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
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