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/function/function_evaluation_tree.hpp>
00036 #include <acado/symbolic_expression/symbolic_expression.hpp>
00037
00038 using namespace std;
00039
00040 BEGIN_NAMESPACE_ACADO
00041
00042
00043
00044
00045
00046
00047
00048 FunctionEvaluationTree::FunctionEvaluationTree( ){
00049
00050 f = NULL;
00051 sub = NULL;
00052 lhs_comp = NULL;
00053 indexList = new SymbolicIndexList();
00054 dim = 0;
00055 n = 0;
00056
00057 globalExportVariableName = "acado_aux";
00058 }
00059
00060 FunctionEvaluationTree::FunctionEvaluationTree( const FunctionEvaluationTree& arg ){
00061
00062 int run1;
00063
00064 dim = arg.dim;
00065 n = arg.n ;
00066
00067 globalExportVariableName = arg.globalExportVariableName;
00068
00069 if( arg.f == NULL ){
00070 f = NULL;
00071 }
00072 else{
00073 f = (Operator**)calloc(dim,sizeof(Operator*));
00074
00075 for( run1 = 0; run1 < dim; run1++ ){
00076 f[run1] = arg.f[run1]->clone();
00077 }
00078 }
00079
00080 indexList = new SymbolicIndexList(*arg.indexList);
00081
00082 if( arg.sub == NULL ){
00083 sub = NULL;
00084 lhs_comp = NULL;
00085 }
00086 else{
00087 sub = (Operator**)calloc(n,sizeof(Operator*));
00088 lhs_comp = (int*)calloc(n,sizeof(int));
00089
00090 for( run1 = 0; run1 < n; run1++ ){
00091 sub[run1] = arg.sub[run1]->clone();
00092 lhs_comp[run1] = arg.lhs_comp[run1];
00093 }
00094 }
00095
00096 safeCopy = arg.safeCopy;
00097 }
00098
00099
00100 FunctionEvaluationTree::~FunctionEvaluationTree( ){
00101
00102 int run1;
00103
00104 for( run1 = 0; run1 < dim; run1++ ){
00105 delete f[run1];
00106 }
00107 if( f != NULL){
00108 free(f);
00109 }
00110
00111 for( run1 = 0; run1 < n; run1++ ){
00112
00113 delete sub[run1];
00114 }
00115
00116 if( sub != NULL ){
00117 free(sub);
00118 }
00119
00120 if( lhs_comp != NULL ){
00121 free(lhs_comp);
00122 }
00123
00124 delete indexList;
00125 }
00126
00127
00128 FunctionEvaluationTree& FunctionEvaluationTree::operator=( const FunctionEvaluationTree& arg ){
00129
00130 int run1;
00131
00132 if( this != &arg ){
00133
00134 for( run1 = 0; run1 < dim; run1++ ){
00135 delete f[run1];
00136 }
00137 if( f != NULL ){
00138 free(f);
00139 }
00140
00141 for( run1 = 0; run1 < n; run1++ ){
00142 delete sub[run1];
00143 }
00144 if( sub != 0 ){
00145 free(sub);
00146 }
00147
00148 if( lhs_comp != NULL ){
00149 free(lhs_comp);
00150 }
00151
00152 delete indexList;
00153
00154 dim = arg.dim;
00155 n = arg.n ;
00156
00157 globalExportVariableName = arg.globalExportVariableName;
00158
00159 if( arg.f == NULL ){
00160 f = NULL;
00161 }
00162 else{
00163 f = (Operator**)calloc(dim,sizeof(Operator*));
00164
00165 for( run1 = 0; run1 < dim; run1++ ){
00166 f[run1] = arg.f[run1]->clone();
00167 }
00168 }
00169
00170 indexList = new SymbolicIndexList(*arg.indexList);
00171
00172 if( arg.sub == NULL ){
00173 sub = NULL;
00174 lhs_comp = NULL;
00175 }
00176 else{
00177 sub = (Operator**)calloc(n,sizeof(Operator*));
00178 lhs_comp = (int*)calloc(n,sizeof(int));
00179
00180 for( run1 = 0; run1 < n; run1++ ){
00181 sub[run1] = arg.sub[run1]->clone();
00182 lhs_comp[run1] = arg.lhs_comp[run1];
00183 }
00184 }
00185 safeCopy = arg.safeCopy;
00186 }
00187
00188 return *this;
00189 }
00190
00191
00192 returnValue FunctionEvaluationTree::getExpression( Expression& expression ) const{
00193
00194 expression = safeCopy;
00195 return SUCCESSFUL_RETURN;
00196 }
00197
00198
00199 returnValue FunctionEvaluationTree::operator<<( const Expression& arg ){
00200
00201 safeCopy << arg;
00202
00203 uint run1;
00204
00205 for( run1 = 0; run1 < arg.getDim(); run1++ ){
00206
00207 int nn;
00208
00209 nn = n;
00210
00211 f = (Operator**)realloc(f,(dim+1)*sizeof(Operator*));
00212
00213 f[dim] = arg.element[run1]->clone();
00214 f[dim]-> loadIndices ( indexList );
00215
00216 sub = (Operator**)realloc(sub,
00217 (indexList->getNumberOfOperators())*sizeof(Operator*));
00218 lhs_comp = (int*)realloc(lhs_comp,
00219 (indexList->getNumberOfOperators())*sizeof(int));
00220
00221 indexList->getOperators( sub, lhs_comp, &n );
00222
00223 while( nn < n ){
00224
00225 sub[nn]-> enumerateVariables( indexList );
00226 nn++;
00227 }
00228
00229 f[dim]-> enumerateVariables( indexList );
00230
00231 dim++;
00232 }
00233 return SUCCESSFUL_RETURN;
00234 }
00235
00236
00237
00238 returnValue FunctionEvaluationTree::evaluate( double *x, double *result ){
00239
00240 int run1;
00241
00242 for( run1 = 0; run1 < n; run1++ ){
00243
00244 sub[run1]->evaluate( 0, x, &x[ indexList->index(VT_INTERMEDIATE_STATE,
00245 lhs_comp[run1] ) ] );
00246 }
00247 for( run1 = 0; run1 < dim; run1++ ){
00248 f[run1]->evaluate( 0, x, &result[run1] );
00249 }
00250
00251 return SUCCESSFUL_RETURN;
00252 }
00253
00254
00255 returnValue FunctionEvaluationTree::evaluate( double *x, double *result, PrintLevel printL ){
00256
00257 int run1;
00258
00259 if (printL == MEDIUM || printL == HIGH)
00260 cout << "Symbolic expression evaluation:" << endl;
00261
00262 for( run1 = 0; run1 < n; run1++ ){
00263
00264 sub[run1]->evaluate( 0, x, &x[ indexList->index(VT_INTERMEDIATE_STATE,
00265 lhs_comp[run1] ) ] );
00266 if( printL == HIGH )
00267 cout << "sub[" << lhs_comp[ run1 ] << "] = "
00268 << scientific << x[ indexList->index(VT_INTERMEDIATE_STATE, lhs_comp[run1]) ]
00269 << endl;
00270 }
00271
00272 for (run1 = 0; run1 < dim; run1++) {
00273 f[run1]->evaluate(0, x, &result[run1]);
00274
00275 if (printL == HIGH || printL == MEDIUM)
00276 cout << "f[" << run1 << "] = " << scientific << result[run1] << endl;
00277 }
00278
00279 return SUCCESSFUL_RETURN;
00280 }
00281
00282
00283
00284 returnValue FunctionEvaluationTree::evaluate( int number, double *x, double *result ){
00285
00286 int run1;
00287
00288 for( run1 = 0; run1 < n; run1++ ){
00289 sub[run1]->evaluate( number, x, &x[ indexList->index(VT_INTERMEDIATE_STATE,
00290 lhs_comp[run1] ) ] );
00291 }
00292 for( run1 = 0; run1 < dim; run1++ ){
00293 f[run1]->evaluate( number, x, &result[run1] );
00294 }
00295
00296 return SUCCESSFUL_RETURN;
00297 }
00298
00299
00300
00301 FunctionEvaluationTree* FunctionEvaluationTree::differentiate( int index_ ){
00302
00303 ACADOERROR(RET_NOT_IMPLEMENTED_YET);
00304
00305 return this;
00306 }
00307
00308
00309 FunctionEvaluationTree FunctionEvaluationTree::substitute( VariableType variableType_, int index_, double sub_ ){
00310
00311 int run1;
00312 FunctionEvaluationTree tmp;
00313
00314 int sub_index = index(variableType_, index_ );
00315
00316 tmp.dim = dim;
00317 tmp.n = n ;
00318
00319 if( f == NULL ){
00320 tmp.f = NULL;
00321 }
00322 else{
00323 tmp.f = (Operator**)calloc(dim,sizeof(Operator*));
00324
00325 Operator *temp;
00326
00327 if( fabs( sub_ ) > 10.0*EPS ){
00328 if( fabs( 1.0 - sub_ ) < 10.0*EPS ){
00329 temp = new DoubleConstant(1.0,NE_ONE);
00330 for( run1 = 0; run1 < dim; run1++ ){
00331 tmp.f[run1] = f[run1]->substitute(sub_index,temp);
00332 }
00333 delete temp;
00334 }
00335 else{
00336 temp = new DoubleConstant(sub_,NE_NEITHER_ONE_NOR_ZERO);
00337 for( run1 = 0; run1 < dim; run1++ ){
00338 tmp.f[run1] = f[run1]->substitute(sub_index,temp);
00339 }
00340 delete temp;
00341 }
00342 }
00343 else{
00344 temp = new DoubleConstant(0.0,NE_ZERO);
00345 for( run1 = 0; run1 < dim; run1++ ){
00346 tmp.f[run1] = f[run1]->substitute(sub_index,temp);
00347 }
00348 delete temp;
00349 }
00350 }
00351
00352 if( sub == NULL ){
00353 tmp.sub = NULL;
00354 tmp.lhs_comp = NULL;
00355 }
00356 else{
00357 tmp.sub = (Operator**)calloc(n,sizeof(Operator*));
00358 tmp.lhs_comp = (int*)calloc(n,sizeof(int));
00359
00360 Operator *temp;
00361
00362 if( fabs(sub_) > 10.0*EPS ){
00363 if( fabs( 1.0-sub_ ) < 10.0*EPS ){
00364 temp = new DoubleConstant(1.0,NE_ONE);
00365 for( run1 = 0; run1 < n; run1++ ){
00366 tmp.sub[run1] = sub[run1]->substitute(sub_index,temp);
00367 tmp.lhs_comp[run1] = lhs_comp[run1];
00368 }
00369 delete temp;
00370 }
00371 else{
00372 temp = new DoubleConstant(sub_,NE_NEITHER_ONE_NOR_ZERO);
00373 for( run1 = 0; run1 < n; run1++ ){
00374 tmp.sub[run1] = sub[run1]->substitute(sub_index,temp);
00375 tmp.lhs_comp[run1] = lhs_comp[run1];
00376 }
00377 delete temp;
00378 }
00379 }
00380 else{
00381 temp = new DoubleConstant(0.0,NE_ZERO);
00382 for( run1 = 0; run1 < n; run1++ ){
00383 tmp.sub[run1] = sub[run1]->substitute(sub_index,temp);
00384 tmp.lhs_comp[run1] = lhs_comp[run1];
00385 }
00386 delete temp;
00387 }
00388 }
00389
00390 delete tmp.indexList;
00391 tmp.indexList = indexList->substitute(variableType_, index_);
00392
00393 return tmp;
00394 }
00395
00396
00397
00398 NeutralElement FunctionEvaluationTree::isOneOrZero(){
00399
00400 int run1;
00401
00402 NeutralElement e = NE_NEITHER_ONE_NOR_ZERO;
00403
00404 if( dim > 0 ){
00405 if( f[0]->isOneOrZero() == NE_ZERO ){
00406 e = NE_ZERO;
00407 }
00408 if( f[0]->isOneOrZero() == NE_ONE ){
00409 e = NE_ONE;
00410 }
00411 }
00412
00413 if( e == NE_NEITHER_ONE_NOR_ZERO ){
00414 return e;
00415 }
00416
00417 for( run1 = 1; run1 < dim; run1++ ){
00418
00419 if( e == NE_ONE ){
00420 if( f[run1]->isOneOrZero() != NE_ONE ){
00421 return NE_NEITHER_ONE_NOR_ZERO;
00422 }
00423 }
00424 if( e == NE_ZERO ){
00425 if( f[run1]->isOneOrZero() != NE_ZERO ){
00426 return NE_NEITHER_ONE_NOR_ZERO;
00427 }
00428 }
00429 }
00430
00431 return e;
00432 }
00433
00434
00435
00436 BooleanType FunctionEvaluationTree::isDependingOn( const Expression &variable ){
00437
00438 int nn = variable.getDim();
00439
00440 int run1;
00441 BooleanType *implicit_dep = new BooleanType [n ];
00442 VariableType *varType = new VariableType[nn];
00443 int *component = new int [nn];
00444
00445 for( run1 = 0; run1 < nn; run1++ ){
00446
00447 Operator *tmp2 = (variable.element[run1])->clone();
00448
00449 if( tmp2->isVariable( varType[run1], component[run1] ) == BT_FALSE ){
00450
00451 ACADOERROR(RET_INVALID_ARGUMENTS);
00452 delete tmp2 ;
00453 delete[] varType ;
00454 delete[] component;
00455 return BT_TRUE ;
00456 }
00457
00458 if( varType[run1] == VT_INTERMEDIATE_STATE ){
00459
00460 ACADOERROR(RET_INVALID_ARGUMENTS);
00461 delete tmp2 ;
00462 delete[] varType ;
00463 delete[] component;
00464 return BT_TRUE ;
00465 }
00466 delete tmp2;
00467 }
00468
00469
00470 for( run1 = 0; run1 < n; run1++ ){
00471 implicit_dep[run1] = sub[run1]->isDependingOn( 1, varType, component, implicit_dep );
00472 }
00473 for( run1 = 0; run1 < dim; run1++ ){
00474 if( f[run1]->isDependingOn( 1, varType, component, implicit_dep ) == BT_TRUE ){
00475 delete[] implicit_dep;
00476 delete[] varType ;
00477 delete[] component;
00478 return BT_TRUE;
00479 }
00480 }
00481
00482 delete[] implicit_dep;
00483 delete[] varType ;
00484 delete[] component;
00485 return BT_FALSE;
00486 }
00487
00488
00489 BooleanType FunctionEvaluationTree::isLinearIn( const Expression &variable ){
00490
00491 int nn = variable.getDim();
00492
00493 int run1;
00494 BooleanType *implicit_dep = new BooleanType [n ];
00495 VariableType *varType = new VariableType[nn];
00496 int *component = new int [nn];
00497
00498 for( run1 = 0; run1 < nn; run1++ ){
00499
00500 Operator *tmp2 = (variable.element[run1])->clone();
00501
00502 if( tmp2->isVariable( varType[run1], component[run1] ) == BT_FALSE ){
00503
00504 ACADOERROR(RET_INVALID_ARGUMENTS);
00505 delete tmp2 ;
00506 delete[] varType ;
00507 delete[] component;
00508 return BT_TRUE ;
00509 }
00510
00511 if( varType[run1] == VT_INTERMEDIATE_STATE ){
00512
00513 ACADOERROR(RET_INVALID_ARGUMENTS);
00514 delete tmp2 ;
00515 delete[] varType ;
00516 delete[] component;
00517 return BT_TRUE ;
00518 }
00519 delete tmp2;
00520 }
00521
00522
00523 for( run1 = 0; run1 < n; run1++ ){
00524 implicit_dep[run1] = sub[run1]->isLinearIn( 1, varType, component, implicit_dep );
00525 }
00526 for( run1 = 0; run1 < dim; run1++ ){
00527 if( f[run1]->isLinearIn( 1, varType, component, implicit_dep ) == BT_FALSE ){
00528 delete[] implicit_dep;
00529 delete[] varType ;
00530 delete[] component;
00531 return BT_FALSE;
00532 }
00533 }
00534
00535 delete[] implicit_dep;
00536 delete[] varType ;
00537 delete[] component;
00538 return BT_TRUE;
00539 }
00540
00541
00542 BooleanType FunctionEvaluationTree::isPolynomialIn( const Expression &variable ){
00543
00544 int nn = variable.getDim();
00545
00546 int run1;
00547 BooleanType *implicit_dep = new BooleanType [n ];
00548 VariableType *varType = new VariableType[nn];
00549 int *component = new int [nn];
00550
00551 for( run1 = 0; run1 < nn; run1++ ){
00552
00553 Operator *tmp2 = (variable.element[run1])->clone();
00554
00555 if( tmp2->isVariable( varType[run1], component[run1] ) == BT_FALSE ){
00556
00557 ACADOERROR(RET_INVALID_ARGUMENTS);
00558 delete tmp2 ;
00559 delete[] varType ;
00560 delete[] component;
00561 return BT_TRUE ;
00562 }
00563
00564 if( varType[run1] == VT_INTERMEDIATE_STATE ){
00565
00566 ACADOERROR(RET_INVALID_ARGUMENTS);
00567 delete tmp2 ;
00568 delete[] varType ;
00569 delete[] component;
00570 return BT_TRUE ;
00571 }
00572 delete tmp2;
00573 }
00574
00575
00576 for( run1 = 0; run1 < n; run1++ ){
00577 implicit_dep[run1] = sub[run1]->isPolynomialIn( 1, varType, component, implicit_dep );
00578 }
00579 for( run1 = 0; run1 < dim; run1++ ){
00580 if( f[run1]->isPolynomialIn( 1, varType, component, implicit_dep ) == BT_FALSE ){
00581 delete[] implicit_dep;
00582 delete[] varType ;
00583 delete[] component;
00584 return BT_FALSE;
00585 }
00586 }
00587
00588 delete[] implicit_dep;
00589 delete[] varType ;
00590 delete[] component;
00591 return BT_TRUE;
00592 }
00593
00594
00595
00596 BooleanType FunctionEvaluationTree::isRationalIn( const Expression &variable ){
00597
00598 int nn = variable.getDim();
00599
00600 int run1;
00601 BooleanType *implicit_dep = new BooleanType [n ];
00602 VariableType *varType = new VariableType[nn];
00603 int *component = new int [nn];
00604
00605 for( run1 = 0; run1 < nn; run1++ ){
00606
00607 Operator *tmp2 = (variable.element[run1])->clone();
00608
00609 if( tmp2->isVariable( varType[run1], component[run1] ) == BT_FALSE ){
00610
00611 ACADOERROR(RET_INVALID_ARGUMENTS);
00612 delete tmp2 ;
00613 delete[] varType ;
00614 delete[] component;
00615 return BT_TRUE ;
00616 }
00617
00618 if( varType[run1] == VT_INTERMEDIATE_STATE ){
00619
00620 ACADOERROR(RET_INVALID_ARGUMENTS);
00621 delete tmp2 ;
00622 delete[] varType ;
00623 delete[] component;
00624 return BT_TRUE ;
00625 }
00626 delete tmp2;
00627 }
00628
00629
00630 for( run1 = 0; run1 < n; run1++ ){
00631 implicit_dep[run1] = sub[run1]->isRationalIn( 1, varType, component, implicit_dep );
00632 }
00633 for( run1 = 0; run1 < dim; run1++ ){
00634 if( f[run1]->isRationalIn( 1, varType, component, implicit_dep ) == BT_FALSE ){
00635 delete[] implicit_dep;
00636 delete[] varType ;
00637 delete[] component;
00638 return BT_FALSE;
00639 }
00640 }
00641
00642 delete[] implicit_dep;
00643 delete[] varType ;
00644 delete[] component;
00645 return BT_TRUE;
00646 }
00647
00648
00649 MonotonicityType FunctionEvaluationTree::getMonotonicity()
00650 {
00651 int run1;
00652 MonotonicityType m = MT_CONSTANT;
00653
00654 for (run1 = 0; run1 < dim; run1++)
00655 {
00656 int mf = f[run1]->getMonotonicity();
00657
00658 switch (mf)
00659 {
00660 case MT_NONDECREASING:
00661 if (m == MT_CONSTANT) m = MT_NONDECREASING;
00662 if (m == MT_NONINCREASING) return MT_NONMONOTONIC;
00663
00664 break;
00665
00666 case MT_NONINCREASING:
00667 if (m == MT_CONSTANT) m = MT_NONINCREASING;
00668 if (m == MT_NONDECREASING) return MT_NONMONOTONIC;
00669
00670 break;
00671
00672 case MT_NONMONOTONIC:
00673 return MT_NONMONOTONIC;
00674
00675 case MT_UNKNOWN:
00676 return MT_UNKNOWN;
00677 }
00678 }
00679 return m;
00680 }
00681
00682
00683
00684 CurvatureType FunctionEvaluationTree::getCurvature( ){
00685
00686 int run1;
00687 CurvatureType m = CT_CONSTANT;
00688
00689 for( run1 = 0; run1 < dim; run1++ ){
00690
00691 CurvatureType mf = f[run1]->getCurvature();
00692
00693 if( mf == CT_AFFINE ){
00694
00695 if( m == CT_CONSTANT ) m = CT_AFFINE;
00696 }
00697
00698 if( mf == CT_CONVEX ){
00699
00700 if( m == CT_CONSTANT ) m = CT_CONVEX;
00701 if( m == CT_AFFINE ) m = CT_CONVEX;
00702 if( m == CT_CONCAVE ) return CT_NEITHER_CONVEX_NOR_CONCAVE;
00703 }
00704
00705 if( mf == CT_CONCAVE ){
00706
00707 if( m == CT_CONSTANT ) m = CT_CONCAVE;
00708 if( m == CT_AFFINE ) m = CT_CONCAVE;
00709 if( m == CT_CONVEX ) return CT_NEITHER_CONVEX_NOR_CONCAVE;
00710 }
00711
00712 if( mf == CT_NEITHER_CONVEX_NOR_CONCAVE ) return CT_NEITHER_CONVEX_NOR_CONCAVE;
00713 }
00714 return m;
00715 }
00716
00717
00718
00719 returnValue FunctionEvaluationTree::AD_forward( double *x, double *seed, double *ff,
00720 double *df ){
00721
00722 int run1;
00723
00724 for( run1 = 0; run1 < n; run1++ ){
00725 sub[run1]->AD_forward( 0, x, seed,
00726 &x [ indexList->index(VT_INTERMEDIATE_STATE, lhs_comp[run1])],
00727 &seed[ indexList->index(VT_INTERMEDIATE_STATE, lhs_comp[run1])] );
00728 }
00729 for( run1 = 0; run1 < dim; run1++ ){
00730 f[run1]->AD_forward( 0, x, seed, &ff[run1], &df[run1] );
00731 }
00732
00733 return SUCCESSFUL_RETURN;
00734 }
00735
00736
00737
00738 returnValue FunctionEvaluationTree::AD_forward( int number, double *x, double *seed,
00739 double *ff, double *df ){
00740
00741 int run1;
00742
00743 for( run1 = 0; run1 < n; run1++ ){
00744 sub[run1]->AD_forward( number, x, seed,
00745 &x [ indexList->index(VT_INTERMEDIATE_STATE, lhs_comp[run1])],
00746 &seed[ indexList->index(VT_INTERMEDIATE_STATE, lhs_comp[run1])] );
00747 }
00748 for( run1 = 0; run1 < dim; run1++ ){
00749 f[run1]->AD_forward( number, x, seed, &ff[run1], &df[run1] );
00750 }
00751
00752 return SUCCESSFUL_RETURN;
00753 }
00754
00755
00756 returnValue FunctionEvaluationTree::AD_forward( int number, double *seed, double *df ){
00757
00758 int run1;
00759
00760 for( run1 = 0; run1 < n; run1++ ){
00761 sub[run1]->AD_forward( number, seed,
00762 &seed[ indexList->index(VT_INTERMEDIATE_STATE, lhs_comp[run1])] );
00763 }
00764 for( run1 = 0; run1 < dim; run1++ ){
00765 f[run1]->AD_forward( number, seed, &df[run1] );
00766 }
00767
00768 return SUCCESSFUL_RETURN;
00769 }
00770
00771
00772 returnValue FunctionEvaluationTree::AD_backward( double *seed, double *df ){
00773
00774 int run1;
00775
00776 for( run1 = dim-1; run1 >= 0; run1-- ){
00777 f[run1]->AD_backward( 0, seed[run1], df );
00778 }
00779
00780
00781 for( run1 = n-1; run1 >= 0; run1-- ){
00782
00783 sub[run1]->AD_backward( 0, df[ indexList->index(VT_INTERMEDIATE_STATE, lhs_comp[run1])],
00784 df
00785 );
00786 }
00787
00788 return SUCCESSFUL_RETURN;
00789 }
00790
00791
00792
00793 returnValue FunctionEvaluationTree::AD_backward( int number, double *seed, double *df ){
00794
00795 int run1;
00796
00797 for( run1 = dim-1; run1 >= 0; run1-- ){
00798 f[run1]->AD_backward( number, seed[run1], df );
00799 }
00800
00801 for( run1 = n-1; run1 >= 0; run1-- ){
00802 sub[run1]->AD_backward( number,
00803 df[ indexList->index(VT_INTERMEDIATE_STATE, lhs_comp[run1])],
00804 df );
00805 }
00806
00807 return SUCCESSFUL_RETURN;
00808 }
00809
00810
00811 returnValue FunctionEvaluationTree::AD_forward2( int number, double *seed,
00812 double *dseed, double *df,
00813 double *ddf ){
00814
00815 int run1;
00816
00817 for( run1 = 0; run1 < n; run1++ ){
00818 sub[run1]->AD_forward2( number, seed, dseed,
00819 &seed [ indexList->index(VT_INTERMEDIATE_STATE, lhs_comp[run1])],
00820 &dseed[ indexList->index(VT_INTERMEDIATE_STATE, lhs_comp[run1])] );
00821 }
00822 for( run1 = 0; run1 < dim; run1++ ){
00823 f[run1]->AD_forward2( number, seed, dseed, &df[run1], &ddf[run1] );
00824 }
00825
00826 return SUCCESSFUL_RETURN;
00827 }
00828
00829
00830 returnValue FunctionEvaluationTree::AD_backward2( int number, double *seed1, double *seed2,
00831 double *df, double *ddf ){
00832
00833 int run1;
00834
00835 for( run1 = dim-1; run1 >= 0; run1-- ){
00836 f[run1]->AD_backward2( number, seed1[run1], seed2[run1], df, ddf );
00837 }
00838
00839
00840 for( run1 = n-1; run1 >= 0; run1-- ){
00841 sub[run1]->AD_backward2( number,
00842 df [ indexList->index(VT_INTERMEDIATE_STATE, lhs_comp[run1])],
00843 ddf [ indexList->index(VT_INTERMEDIATE_STATE, lhs_comp[run1])],
00844 df,
00845 ddf
00846 );
00847 }
00848
00849 return SUCCESSFUL_RETURN;
00850 }
00851
00852
00853
00854 returnValue FunctionEvaluationTree::C_print( std::ostream& stream,
00855 const char *fcnName ,
00856 const char *realString
00857 ) const
00858 {
00859 stream << "/* This file was auto-generated by ACADO Toolkit. */" << endl << endl;
00860
00861 exportForwardDeclarations(stream, fcnName, realString);
00862 stream << endl;
00863 exportCode(stream, fcnName, realString);
00864
00865 return SUCCESSFUL_RETURN;
00866 }
00867
00868 returnValue FunctionEvaluationTree::exportForwardDeclarations( std::ostream& stream,
00869 const char *fcnName,
00870 const char *realString
00871 ) const
00872 {
00873 stream <<
00874 "\n"
00875 "/** Export of an ACADO symbolic function.\n"
00876 " *\n"
00877 " * \\param in Input to the exported function.\n"
00878 " * \\param out Output of the exported function.\n"
00879 " */\n"
00880 << "void " << fcnName << "(const " << realString << "* in, "
00881 << realString << "* out);" << endl;
00882
00883 return SUCCESSFUL_RETURN;
00884 }
00885
00886
00887 returnValue FunctionEvaluationTree::exportCode( std::ostream& stream,
00888 const char *fcnName,
00889 const char *realString,
00890 uint _numX,
00891 uint _numXA,
00892 uint _numU,
00893 uint _numP,
00894 uint _numDX,
00895 uint _numOD,
00896 bool allocateMemory,
00897 bool staticMemory
00898 ) const{
00899
00900 int run1;
00901 int nni = 0;
00902
00903 for (run1 = 0; run1 < n; run1++)
00904 if (lhs_comp[run1] + 1 > nni)
00905 nni = lhs_comp[run1] + 1;
00906
00907 unsigned numX = _numX > 0 ? _numX : getNX();
00908 unsigned numXA = _numXA > 0 ? _numXA : getNXA();
00909 unsigned numU = _numU > 0 ? _numU : getNU();
00910 unsigned numP = _numP > 0 ? _numP : getNP();
00911 unsigned numDX = _numDX > 0 ? _numDX : getNDX();
00912 unsigned numOD = _numOD > 0 ? _numOD : getNOD();
00913
00914 unsigned offset = 0;
00915
00916 stream << "void " << fcnName << "(const " << realString << "* in, " << realString << "* out)\n{\n";
00917
00918 if (numX > 0)
00919 stream << "const " << realString << "* xd = in;" << endl;
00920 offset += numX;
00921
00922 if (numXA > 0)
00923 stream << "const " << realString << "* xa = in + " << offset << ";" << endl;
00924 offset += numXA;
00925
00926 if (getNU() > 0)
00927 stream << "const " << realString << "* u = in + " << offset << ";" << endl;
00928 offset += numU;
00929
00930 if (getNUI() > 0)
00931 stream << "const " << realString << "* v = in + " << offset << ";" << endl;
00932 offset += getNUI();
00933
00934 if (numP > 0)
00935 stream << "const " << realString << "* p = in + " << offset << ";" << endl;
00936 offset += numP;
00937
00938 if (getNOD() > 0)
00939 stream << "const " << realString << "* od = in + " << offset << ";" << endl;
00940 offset += numOD;
00941
00942 if (getNPI() > 0)
00943 stream << "const " << realString << "* q = in + " << offset << ";" << endl;
00944 offset += getNPI();
00945
00946 if (getNW() > 0)
00947 stream << "const " << realString << "* w = in + " << offset << ";" << endl;
00948 offset += getNW();
00949
00950 if (numDX > 0)
00951 stream << "const " << realString << "* dx = in + " << offset << ";" << endl;
00952 offset += numDX;
00953
00954 if (getNT() > 0)
00955 stream << "const " << realString << "* t = in + " << offset << ";" << endl;
00956 offset += getNT();
00957
00958 if (n > 0)
00959 {
00960 stream << "/* Vector of auxiliary variables; number of elements: " << n << ". */" << endl;
00961
00962 if ( allocateMemory )
00963 {
00964 if ( staticMemory )
00965 {
00966 stream << "static ";
00967 }
00968 stream << realString << " a[" << n << "];";
00969 }
00970 else
00971 stream << realString << "* a = " << globalExportVariableName << ";";
00972 stream << endl << endl;
00973
00974 stream << "/* Compute intermediate quantities: */" << endl;
00975 }
00976
00977 vector< string > auxVarIndividualNames;
00978 auxVarIndividualNames.resize( nni );
00979 for (run1 = 0; run1 < n; run1++)
00980 {
00981 stringstream ss;
00982 ss << "a" << "[" << run1 << "]";
00983 auxVarIndividualNames[ lhs_comp[ run1 ] ] = ss.str();
00984 }
00985
00986 IoFormatter iof( stream );
00987 iof.set(16, iof.width, ios::scientific);
00988
00989
00990 for (run1 = 0; run1 < n; run1++)
00991 {
00992
00993 sub[run1]->setVariableExportName(VT_INTERMEDIATE_STATE, auxVarIndividualNames);
00994
00995 stream << "a[" << run1 << "] = " << *sub[ run1 ] << ";" << endl;
00996 }
00997
00998
00999 stream << endl << "/* Compute outputs: */" << endl;
01000 for (run1 = 0; run1 < dim; run1++)
01001 {
01002
01003 f[run1]->setVariableExportName(VT_INTERMEDIATE_STATE, auxVarIndividualNames);
01004
01005 stream << "out[" << run1 << "] = " << *f[ run1 ] << ";" << endl;
01006 }
01007
01008 iof.reset();
01009
01010 stream << "}" << endl << endl;
01011
01012 return SUCCESSFUL_RETURN;
01013 }
01014
01015
01016 returnValue FunctionEvaluationTree::clearBuffer(){
01017
01018 int run1;
01019 returnValue returnvalue;
01020
01021 for( run1 = 0; run1 < n; run1++ ){
01022 returnvalue = sub[run1]->clearBuffer();
01023 if( returnvalue != SUCCESSFUL_RETURN ){
01024 return returnvalue;
01025 }
01026 }
01027 for( run1 = 0; run1 < dim; run1++ ){
01028 returnvalue = f[run1]->clearBuffer();
01029 if( returnvalue != SUCCESSFUL_RETURN ){
01030 return returnvalue;
01031 }
01032 }
01033
01034 return SUCCESSFUL_RETURN;
01035 }
01036
01037
01038 returnValue FunctionEvaluationTree::makeImplicit(){
01039
01040 return makeImplicit(dim);
01041 }
01042
01043
01044 returnValue FunctionEvaluationTree::makeImplicit( int dim_ ){
01045
01046 int run1;
01047 int var_counter = indexList->makeImplicit(dim_);
01048
01049 for( run1 = 0; run1 < dim_; run1++ ){
01050
01051 Operator *tmp = f[run1]->clone();
01052 delete f[run1];
01053 Projection pp;
01054 pp.variableType = VT_DDIFFERENTIAL_STATE;
01055 pp.vIndex = run1 ;
01056 pp.variableIndex = var_counter-dim_+run1 ;
01057
01058 f[run1] = new Subtraction( pp.clone(), tmp->clone() );
01059 delete tmp;
01060 }
01061
01062 return SUCCESSFUL_RETURN;
01063 }
01064
01065
01066 int FunctionEvaluationTree::getNX () const{
01067
01068 return indexList->getNX();
01069 }
01070
01071 int FunctionEvaluationTree::getNXA () const{
01072
01073 return indexList->getNXA();
01074 }
01075
01076 int FunctionEvaluationTree::getNDX () const{
01077
01078 return indexList->getNDX();
01079 }
01080
01081 int FunctionEvaluationTree::getNU () const{
01082
01083 return indexList->getNU();
01084 }
01085
01086 int FunctionEvaluationTree::getNUI () const{
01087
01088 return indexList->getNUI();
01089 }
01090
01091 int FunctionEvaluationTree::getNP () const{
01092
01093 return indexList->getNP();
01094 }
01095
01096 int FunctionEvaluationTree::getNPI () const{
01097
01098 return indexList->getNPI();
01099 }
01100
01101 int FunctionEvaluationTree::getNW () const{
01102
01103 return indexList->getNW();
01104 }
01105
01106 int FunctionEvaluationTree::getNT () const{
01107
01108 return indexList->getNT();
01109 }
01110
01111 int FunctionEvaluationTree::getNOD () const{
01112
01113 return indexList->getOD();
01114 }
01115
01116 int FunctionEvaluationTree::index( VariableType variableType_, int index_ ) const{
01117
01118 return indexList->index( variableType_, index_ );
01119 }
01120
01121 double FunctionEvaluationTree::scale( VariableType variableType_, int index_ ) const{
01122
01123 return indexList->scale( variableType_, index_ );
01124 }
01125
01126
01127 int FunctionEvaluationTree::getNumberOfVariables() const{
01128
01129 return indexList->getNumberOfVariables();
01130 }
01131
01132
01133 Operator* FunctionEvaluationTree::getExpression( uint componentIdx
01134 ) const
01135 {
01136 if ( (int)componentIdx < getDim( ) )
01137 return f[componentIdx]->clone( );
01138 else
01139 return new DoubleConstant( 0.0,NE_ZERO );
01140 }
01141
01142
01143 BooleanType FunctionEvaluationTree::isSymbolic() const{
01144
01145 int run1;
01146 for( run1 = 0; run1 < n; run1++ ){
01147 if( sub[run1]->isSymbolic() == BT_FALSE ) return BT_FALSE;
01148 }
01149 for( run1 = 0; run1 < dim; run1++ ){
01150 if( f[run1]->isSymbolic() == BT_FALSE ) return BT_FALSE;
01151 }
01152 return BT_TRUE;
01153 }
01154
01155
01156 returnValue FunctionEvaluationTree::setScale( double *scale_ )
01157 {
01158 return ACADOERROR(RET_INVALID_USE_OF_FUNCTION);
01159 }
01160
01161 returnValue FunctionEvaluationTree::setGlobalExportVariableName(const std::string& _name)
01162 {
01163 if (_name.size() == 0)
01164 return ACADOERROR( RET_INVALID_ARGUMENTS );
01165
01166 globalExportVariableName = _name;
01167
01168 return SUCCESSFUL_RETURN;
01169 }
01170
01171 std::string FunctionEvaluationTree::getGlobalExportVariableName() const
01172 {
01173 return globalExportVariableName;
01174 }
01175
01176 unsigned FunctionEvaluationTree::getGlobalExportVariableSize() const
01177 {
01178 return n;
01179 }
01180
01181 CLOSE_NAMESPACE_ACADO
01182
01183