expression.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 
00034 #include <acado/symbolic_expression/expression.hpp>
00035 #include <acado/symbolic_expression/symbolic_expression.hpp>
00036 #include <acado/symbolic_expression/acado_syntax.hpp>
00037 #include <acado/function/function.hpp>
00038 #include <acado/symbolic_operator/symbolic_operator.hpp>
00039 
00040 BEGIN_NAMESPACE_ACADO
00041 
00042 
00043 // ---------------------------------------------------------------------------------------------------
00044 //                                 IMPLEMENTATION OF THE CONSTRUCTORS:
00045 // ---------------------------------------------------------------------------------------------------
00046 
00047 Expression::Expression()
00048 {
00049         construct(VT_UNKNOWN, 0, 0, 0, "");
00050 }
00051 
00052 Expression::Expression(const std::string& name_, uint nRows_, uint nCols_, VariableType variableType_, uint globalTypeID)
00053 {
00054         construct(variableType_, globalTypeID, nRows_, nCols_, name_);
00055 }
00056 
00057 Expression::Expression(const std::string &name_)
00058 {
00059         construct(VT_UNKNOWN, 0, 0, 0, name_);
00060 }
00061 
00062 Expression::Expression(uint nRows_, uint nCols_, VariableType variableType_, uint globalTypeID)
00063 {
00064         construct(variableType_, globalTypeID, nRows_, nCols_, "");
00065 }
00066 
00067 Expression::Expression(int nRows_, int nCols_, VariableType variableType_, int globalTypeID)
00068 {
00069         construct(variableType_, globalTypeID, nRows_, nCols_, "");
00070 }
00071 
00072 Expression::Expression(const Operator &tree_)
00073 {
00074         VariableType tmpType;
00075         int tmpComp;
00076 
00077         if (tree_.isVariable(tmpType, tmpComp) == BT_TRUE)
00078         {
00079                 construct(tmpType, tmpComp, 1, 1, "");
00080         }
00081         else
00082         {
00083                 construct(VT_UNKNOWN, 0, 1, 1, "");
00084         }
00085 
00086         delete element[0];
00087         element[0] = tree_.clone();
00088 }
00089 
00090 Expression::Expression( const double& rhs )
00091 {
00092         construct(VT_UNKNOWN, 0, 1, 1, "");
00093         delete element[ 0 ];
00094         element[ 0 ] = new DoubleConstant(rhs, NE_NEITHER_ONE_NOR_ZERO);
00095 }
00096 
00097 Expression::Expression( const DVector& rhs )
00098 {
00099         construct(VT_UNKNOWN, 0, rhs.getDim(), 1, "");
00100         for(unsigned el = 0; el < rhs.getDim(); el++ )
00101         {
00102                 delete element[ el ];
00103                 element[ el ] = new DoubleConstant(rhs( el ), NE_NEITHER_ONE_NOR_ZERO);
00104         }
00105 }
00106 
00107 Expression::Expression( const DMatrix& rhs )
00108 {
00109         construct(VT_UNKNOWN, 0, rhs.getNumRows(), rhs.getNumCols(), "");
00110         for(unsigned run1 = 0; run1 < rhs.getNumRows(); run1++ )
00111         {
00112                 for(unsigned run2 = 0; run2 < rhs.getNumCols(); run2++ )
00113                 {
00114                         delete element[rhs.getNumCols() * run1 + run2];
00115                         element[rhs.getNumCols() * run1 + run2] = 
00116                                 new DoubleConstant(rhs(run1, run2), NE_NEITHER_ONE_NOR_ZERO);
00117                 }
00118         }
00119 }
00120 
00121 Expression::Expression(const Expression& rhs)
00122 {
00123         copy(rhs);
00124 }
00125 
00126 Expression::~Expression()
00127 {
00128         deleteAll();
00129 }
00130 
00131 Expression& Expression::operator=( const Expression& rhs )
00132 {
00133         if (this != &rhs)
00134         {
00135                 deleteAll();
00136                 copy(rhs);
00137         }
00138         return *this;
00139 }
00140 
00141 Expression&  Expression::appendRows(const Expression& arg) {
00142         if (getDim()==0) {operator=(arg);return *this;}
00143         ASSERT(arg.getNumCols() == getNumCols());
00144         
00145         uint run1;
00146         uint oldDim = dim;
00147         dim     += arg.dim;
00148         nRows += arg.getNumRows();
00149 
00150         if( arg.variableType != variableType )
00151                 variableType = VT_UNKNOWN;
00152         element = (Operator**)realloc(element, dim*sizeof(Operator*) );
00153             
00154         for( run1 = oldDim; run1 < dim; run1++ )
00155                 element[run1] = arg.element[run1-oldDim]->clone();
00156         
00157         return *this;
00158 }
00159 
00160 // this is still very unefficient code
00161 
00162 Expression& Expression::appendCols(const Expression& arg) {
00163         if (getDim()==0) {operator=(arg);return *this;}
00164         ASSERT(arg.getNumRows() == getNumRows());
00165         
00166         
00167         Expression E(transpose());
00168         E.appendRows(arg.transpose());
00169         
00170         operator=(E.transpose());
00171         
00172         return *this;
00173 }
00174 
00175 Expression& Expression::operator<<( const Expression& arg ){
00176   
00177     if( dim == 0 ) return operator=(arg);
00178 
00179     uint run1;
00180     uint oldDim = dim;
00181 
00182     dim   += arg.dim  ;
00183     nRows += arg.dim  ;
00184 
00185     variableType = VT_VARIABLE;
00186 
00187     if( arg.isVariable() == BT_FALSE ) variableType = VT_UNKNOWN;
00188 
00189     element = (Operator**)realloc(element, dim*sizeof(Operator*) );
00190 
00191     for( run1 = oldDim; run1 < dim; run1++ )
00192         element[run1] = arg.element[run1-oldDim]->clone();
00193 
00194     return *this;
00195 }
00196 
00197 
00198 std::ostream& Expression::print(std::ostream& stream) const
00199 {
00200         uint run1;
00201         stream << "[ ";
00202         if (dim) {
00203                 for (run1 = 0; run1 < dim - 1; run1++)
00204                         stream << *element[run1] << " , ";
00205                 stream << *element[dim - 1];
00206         }
00207         stream << "]";
00208 
00209         return stream;
00210 }
00211 
00212 
00213 std::ostream& operator<<( std::ostream& stream, const Expression &arg )
00214 {
00215     return arg.print(stream);
00216 }
00217 
00218 Expression Expression::operator()( uint idx ) const{
00219 
00220     ASSERT( idx < getDim( ) );
00221 
00222     Expression tmp(1);
00223 
00224     delete tmp.element[0];
00225     tmp.element[0] = element[idx]->clone();
00226 
00227     tmp.component    = component + idx;
00228     tmp.variableType = variableType;
00229 
00230     return tmp;
00231 }
00232 
00233 Expression Expression::operator()( uint rowIdx, uint colIdx ) const{
00234 
00235     ASSERT( rowIdx < getNumRows( ) );
00236     ASSERT( colIdx < getNumCols( ) );
00237 
00238     Expression tmp(1);
00239 
00240     delete tmp.element[0];
00241     tmp.element[0] = element[rowIdx*getNumCols()+colIdx]->clone();
00242 
00243     tmp.component    = component + rowIdx*getNumCols() + colIdx;
00244     tmp.variableType = variableType;
00245 
00246     return tmp;
00247 }
00248 
00249 
00250 Operator& Expression::operator()( uint idx ){
00251 
00252     switch( variableType ){
00253 
00254         case  VT_INTERMEDIATE_STATE:
00255               ASSERT( idx < getDim( ) );
00256               return *element[idx];
00257 
00258         case VT_UNKNOWN:
00259               ASSERT( idx < getDim( ) );
00260               delete  element[idx];
00261               element[idx] = new TreeProjection();
00262               return *element[idx];
00263 
00264         default:
00265               ASSERT( idx < getDim( ) );
00266               return *element[idx];
00267     }
00268     ASSERT( 1 == 0 );
00269     return *element[0];
00270 }
00271 
00272 
00273 Operator& Expression::operator()( uint rowIdx, uint colIdx ){
00274 
00275     switch( variableType ){
00276 
00277         case  VT_INTERMEDIATE_STATE:
00278               ASSERT( rowIdx < getNumRows( ) );
00279               ASSERT( colIdx < getNumCols( ) );
00280               return *element[rowIdx*getNumCols()+colIdx];
00281 
00282 //        case  VT_UNKNOWN:
00283         default:
00284               ASSERT( rowIdx < getNumRows( ) );
00285               ASSERT( colIdx < getNumCols( ) );
00286               delete  element[rowIdx*getNumCols()+colIdx];
00287               element[rowIdx*getNumCols()+colIdx] = new TreeProjection();
00288               return *element[rowIdx*getNumCols()+colIdx];
00289 
00290 //               ASSERT( 1 == 0 );
00291 //               return *element[0];
00292     }
00293     ASSERT( 1 == 0 );
00294     return *element[0];
00295 }
00296 
00297 
00298 Expression operator+( const Expression  & arg1, const Expression  & arg2 )
00299 {return arg1.add(arg2);}
00300 Expression operator-( const Expression  & arg1, const Expression  & arg2 )
00301 {return arg1.sub(arg2);}
00302 Expression operator*( const Expression  & arg1, const Expression  & arg2 )
00303 {return arg1.mul(arg2);}
00304 Expression operator/( const Expression  & arg1, const Expression  & arg2 )
00305 {return arg1.div(arg2);}
00306 
00307 Expression& Expression::operator+=(const Expression &arg)
00308 {return static_cast<Expression&>(*this) = static_cast<Expression*>(this)->add(arg);}
00309 
00310 Expression& Expression::operator-=(const Expression &arg)
00311 {return static_cast<Expression&>(*this) = static_cast<Expression*>(this)->sub(arg);}
00312 
00313 Expression& Expression::operator*=(const Expression &arg)
00314 {return static_cast<Expression&>(*this) = static_cast<Expression*>(this)->mul(arg);}
00315 
00316 Expression& Expression::operator/=(const Expression &arg)
00317 {return static_cast<Expression&>(*this) = static_cast<Expression*>(this)->div(arg);}
00318 
00319 Expression Expression::add( const Expression& arg ) const{
00320 
00321     ASSERT( getNumRows() == arg.getNumRows() );
00322     ASSERT( getNumCols() == arg.getNumCols() );
00323 
00324     uint i,j;
00325 
00326     Expression tmp("", getNumRows(), getNumCols() );
00327 
00328     for( i=0; i<getNumRows(); ++i ){
00329         for( j=0; j<getNumCols(); ++j ){
00330 
00331             delete tmp.element[i*getNumCols()+j];
00332             if( element[i*getNumCols()+j]->isOneOrZero() != NE_ZERO ){
00333                 if( arg.element[i*getNumCols()+j]->isOneOrZero() != NE_ZERO )
00334                     tmp.element[i*getNumCols()+j] = new Addition( element[i*getNumCols()+j]->clone(),
00335                                                     arg.element[i*getNumCols()+j]->clone() );
00336                 else
00337                     tmp.element[i*getNumCols()+j] = element[i*getNumCols()+j]->clone();
00338             }
00339             else{
00340                 if( arg.element[i*getNumCols()+j]->isOneOrZero() != NE_ZERO )
00341                      tmp.element[i*getNumCols()+j] = arg.element[i*getNumCols()+j]->clone();
00342                 else tmp.element[i*getNumCols()+j] = new DoubleConstant(0.0,NE_ZERO);
00343             }
00344         }
00345     }
00346     return tmp;
00347 }
00348 
00349 
00350 Expression Expression::sub( const Expression& arg ) const{
00351 
00352     ASSERT( getNumRows() == arg.getNumRows() );
00353     ASSERT( getNumCols() == arg.getNumCols() );
00354 
00355     uint i,j;
00356 
00357     Expression tmp("", getNumRows(), getNumCols() );
00358 
00359     for( i=0; i<getNumRows(); ++i ){
00360         for( j=0; j<getNumCols(); ++j ){
00361 
00362             delete tmp.element[i*getNumCols()+j];
00363             if( element[i*getNumCols()+j]->isOneOrZero() != NE_ZERO ){
00364                 if( arg.element[i*getNumCols()+j]->isOneOrZero() != NE_ZERO )
00365                     tmp.element[i*getNumCols()+j] = new Subtraction( element[i*getNumCols()+j]->clone(),
00366                                                     arg.element[i*getNumCols()+j]->clone() );
00367                 else
00368                     tmp.element[i*getNumCols()+j] = element[i*getNumCols()+j]->clone();
00369             }
00370             else{
00371                 if( arg.element[i*getNumCols()+j]->isOneOrZero() != NE_ZERO )
00372                      tmp.element[i*getNumCols()+j] = new Subtraction( new DoubleConstant(0.0,NE_ZERO),
00373                                                                       arg.element[i*getNumCols()+j]->clone() );
00374                 else tmp.element[i*getNumCols()+j] = new DoubleConstant(0.0,NE_ZERO);
00375             }
00376         }
00377     }
00378     return tmp;
00379 }
00380 
00381 
00382 Operator* Expression::product( const Operator *a, const Operator *b ) const{
00383 
00384 
00385     switch( a->isOneOrZero() ){
00386 
00387         case NE_ZERO:
00388              return new DoubleConstant( 0.0, NE_ZERO );
00389 
00390         case NE_ONE:
00391              return b->clone();
00392 
00393         default:
00394 
00395              switch( b->isOneOrZero() ){
00396 
00397                  case NE_ZERO:
00398                       return new DoubleConstant( 0.0, NE_ZERO );
00399 
00400                  case NE_ONE:
00401                       return a->clone();
00402 
00403                  default:
00404                       return new Product( a->clone(), b->clone() );
00405              }
00406     }
00407     return 0;
00408 }
00409 
00410 
00411 
00412 Expression Expression::mul( const Expression& arg ) const{
00413     if (getDim()==0 || arg.getDim()==0) return Expression();
00414     // uninitialized expressions yield uninitialized results
00415 
00416     uint i,j,k;
00417 
00418     if( getNumRows() == 1 && getNumCols( ) == 1 ){
00419 
00420         Expression tmp("", arg.getNumRows(), arg.getNumCols() );
00421 
00422         for( i = 0; i< arg.getDim(); i++ ){
00423 
00424              delete tmp.element[i];
00425              Operator *prod = product( element[0], arg.element[i] );
00426              tmp.element[i] = prod->clone();
00427 
00428              delete prod;
00429         }
00430         return tmp;
00431     }
00432 
00433 
00434     if( arg.getNumRows() == 1 && arg.getNumCols( ) == 1 ){
00435 
00436         Expression tmp("", getNumRows(), getNumCols() );
00437 
00438         for( i = 0; i< getDim(); i++ ){
00439 
00440              delete tmp.element[i];
00441              Operator *prod = product( arg.element[0], element[i] );
00442              tmp.element[i] = prod->clone();
00443 
00444              delete prod;
00445         }
00446         return tmp;
00447     }
00448 
00449 
00450     ASSERT( getNumCols( ) == arg.getNumRows( ) );
00451 
00452     uint newNumRows = getNumRows( );
00453     uint newNumCols = arg.getNumCols( );
00454     IntermediateState tmp = zeros<double>( newNumRows, newNumCols );
00455 
00456     for( i=0; i<newNumRows; ++i ){
00457         for( j=0; j<newNumCols; ++j ){
00458             for( k=0; k<getNumCols( ); ++k ){
00459 
00460                  Operator *tmpO = product( element[i*getNumCols()+k],
00461                                            arg.element[k*arg.getNumCols()+j] );
00462 
00463                  if( tmpO->isOneOrZero() != NE_ZERO )
00464                      tmp(i,j) += *tmpO;
00465 
00466                  delete tmpO;
00467             }
00468         }
00469     }
00470     return tmp;
00471 }
00472 
00473 
00474 Expression Expression::div( const Expression& arg ) const{
00475 
00476     ASSERT( arg.getNumRows() == 1 );
00477     ASSERT( arg.getNumCols() == 1 );
00478 
00479     uint i;
00480 
00481     Expression tmp("", getNumRows(), getNumCols() );
00482 
00483     for( i = 0; i< getDim(); i++ ){
00484          delete tmp.element[i];
00485          tmp.element[i] = new Quotient( element[i]->clone(), arg.element[0]->clone() );
00486     }
00487     return tmp;
00488 }
00489 
00490 Expression Expression::getInverse() const{
00491 
00492     ASSERT( getNumRows() == getNumCols() );
00493 
00494     int i,j,k;                 // must really be int, never use uint here.
00495     int M = getNumRows();      // must really be int, never use uint here.
00496 
00497         IntermediateState tmp("", M, 2 * M);
00498 
00499     for( i = 0; i < M; i++ ){
00500         for( j = 0; j < 2*M; j++ ){
00501             if( j <  M   ) tmp(i,j) = operator()(i,j);
00502             else           tmp(i,j) = 0.0        ;
00503             if( j == M+i ) tmp(i,j) = 1.0        ;
00504         }
00505     }
00506 
00507     for( i = 0; i < M-1; i++ )
00508          for( k = i+1; k < M; k++ )
00509              for( j = 2*M-1; j >= i; j-- )
00510                   tmp(k,j) -= ( tmp(i,j)*tmp(k,i) )/tmp(i,i);
00511 
00512     for( i = M-1; i > 0; i-- )
00513         for( k = i-1; k >= 0; k-- )
00514             for( j = 2*M-1; j >= i; j-- )
00515                 tmp(k,j) -= (tmp(i,j)*tmp(k,i))/tmp(i,i);
00516 
00517     for( i = 0; i < M; i++ )
00518         for( j = 2*M-1; j >= 0; j-- )
00519              tmp(i,j) = tmp(i,j)/tmp(i,i);
00520 
00521     Expression I("", M,M);
00522     for( i = 0; i < M; i++ ){
00523         for( j = 0; j < M; j++ ){
00524             delete I.element[i*M+j];
00525             I.element[i*M+j] = tmp.element[i*2*M+j+M]->clone();
00526         }
00527     }
00528 
00529     return I;
00530 }
00531 
00532 
00533 Expression Expression::getRow( const uint& rowIdx ) const{
00534 
00535     uint run1;
00536     ASSERT( rowIdx < getNumRows() );
00537 
00538     Expression tmp("", 1, (int) getNumCols() );
00539 
00540     for( run1 = 0; run1 < getNumCols(); run1++ ){
00541         delete tmp.element[run1];
00542         tmp.element[run1] = element[rowIdx*getNumCols()+run1]->clone();
00543     }
00544     return tmp;
00545 }
00546 
00547 
00548 Expression Expression::getRows( const uint& rowIdx1, const uint& rowIdx2 ) const{
00549 
00550     uint run1, run2, _nRows;
00551     ASSERT( rowIdx1 < getNumRows() );
00552     ASSERT( rowIdx2 <= getNumRows() );
00553     _nRows = rowIdx2 - rowIdx1;
00554 
00555     Expression tmp("", _nRows, (int) getNumCols() );
00556 
00557     for( run1 = 0; run1 < _nRows; run1++ ){
00558         for( run2 = 0; run2 < getNumCols(); run2++ ){
00559                 delete tmp.element[run1*getNumCols()+run2];
00560                 tmp.element[run1*getNumCols()+run2] = element[(rowIdx1+run1)*getNumCols()+run2]->clone();
00561         }
00562     }
00563     return tmp;
00564 }
00565 
00566 
00567 Expression Expression::getCol( const uint& colIdx ) const{
00568 
00569     uint run1;
00570     ASSERT( colIdx < getNumCols() );
00571 
00572     Expression tmp("", (int) getNumRows(), 1 );
00573 
00574     for( run1 = 0; run1 < getNumRows(); run1++ ){
00575         delete tmp.element[run1];
00576         tmp.element[run1] = element[run1*getNumCols()+colIdx]->clone();
00577     }
00578     return tmp;
00579 }
00580 
00581 
00582 Expression Expression::getCols( const uint& colIdx1, const uint& colIdx2 ) const{
00583 
00584     uint run1, run2, _nCols;
00585     ASSERT( colIdx1 < getNumCols() );
00586     ASSERT( colIdx2 <= getNumCols() );
00587     _nCols = colIdx2 - colIdx1;
00588 
00589     Expression tmp("", (int) getNumRows(), _nCols );
00590 
00591     for( run1 = 0; run1 < getNumRows(); run1++ ){
00592         for( run2 = 0; run2 < _nCols; run2++ ){
00593                 delete tmp.element[run1*_nCols+run2];
00594                 tmp.element[run1*_nCols+run2] = element[run1*getNumCols()+colIdx1+run2]->clone();
00595         }
00596     }
00597     return tmp;
00598 }
00599 
00600 
00601 Expression Expression::getSubMatrix( const uint& rowIdx1, const uint& rowIdx2, const uint& colIdx1, const uint& colIdx2 ) const{
00602 
00603     uint run1, run2, _nRows, _nCols;
00604     ASSERT( colIdx1 < getNumCols() );
00605     ASSERT( colIdx2 <= getNumCols() );
00606     _nCols = colIdx2 - colIdx1;
00607     ASSERT( rowIdx1 < getNumRows() );
00608     ASSERT( rowIdx2 <= getNumRows() );
00609     _nRows = rowIdx2 - rowIdx1;
00610 
00611     Expression tmp("", _nRows, _nCols );
00612 
00613     for( run1 = 0; run1 < _nRows; run1++ ){
00614         for( run2 = 0; run2 < _nCols; run2++ ){
00615                 delete tmp.element[run1*_nCols+run2];
00616                 tmp.element[run1*_nCols+run2] = element[(rowIdx1+run1)*getNumCols()+colIdx1+run2]->clone();
00617         }
00618     }
00619     return tmp;
00620 }
00621 
00622 
00623 DMatrix Expression::getDependencyPattern( const Expression& arg ) const{
00624 
00625         DMatrix tmp;
00626         if( arg.getDim() == 0 ) return tmp;
00627 
00628     Function f;
00629     f << backwardDerivative( *this, arg );
00630 
00631 //    FILE *test=fopen("test-getdep.txt","w");
00632 //    test << f;
00633 //    fclose(test);
00634     
00635     double *x      = new double[ f.getNumberOfVariables()+1 ];
00636     double *result = new double[ f.getDim()                 ];
00637 
00638     // TODO: DANGEROUS CODE --> TALK TO MILAN ABOUT THIS
00639     int run1;
00640     srand(1.0);
00641     for( run1 = 0; run1 < f.getNumberOfVariables()+1; run1++ )
00642         x[run1] = 1.0 + (double)rand() / RAND_MAX;
00643 
00644     // EVALUATE f AT THE POINT  (tt,xx):
00645     // ---------------------------------
00646     f.evaluate( 0, x, result );
00647 
00648     tmp = DMatrix( getDim(), arg.getDim(), result );
00649 
00650     delete[] result;
00651     delete[] x;
00652 
00653     return tmp;
00654 }
00655 
00656 DMatrix Expression::getSparsityPattern() const
00657 {
00658         DMatrix res = zeros<double>(getNumRows(), getNumCols());
00659 
00660         for (unsigned el = 0; el < getDim(); ++el)
00661         {
00662                 Operator* foo = getOperatorClone(el);
00663                 if (foo->isOneOrZero() != NE_ZERO)
00664                         res(el) = 1.0;
00665                 delete foo;
00666         }
00667 
00668         return res;
00669 }
00670 
00671 
00672 
00673 Expression Expression::getSin( ) const{
00674 
00675         Expression tmp("", nRows, nCols);
00676     uint run1;
00677 
00678     for( run1 = 0; run1 < dim; run1++ ){
00679         delete tmp.element[run1];
00680         tmp.element[run1] = new Sin( element[run1]->clone() ); 
00681     }
00682     return tmp;
00683 }
00684 
00685 Expression Expression::getCos( ) const{
00686 
00687         Expression tmp("", nRows, nCols);
00688     uint run1;
00689 
00690     for( run1 = 0; run1 < dim; run1++ ){
00691         delete tmp.element[run1];
00692         tmp.element[run1] = new Cos( element[run1]->clone() ); 
00693     }
00694     return tmp;
00695 }
00696 
00697 Expression Expression::getTan( ) const{
00698 
00699         Expression tmp("", nRows, nCols);
00700     uint run1;
00701 
00702     for( run1 = 0; run1 < dim; run1++ ){
00703         delete tmp.element[run1];
00704         tmp.element[run1] = new Tan( element[run1]->clone() ); 
00705     }
00706     return tmp;
00707 }
00708 
00709 Expression Expression::getAsin( ) const{
00710 
00711         Expression tmp("", nRows, nCols);
00712     uint run1;
00713 
00714     for( run1 = 0; run1 < dim; run1++ ){
00715         delete tmp.element[run1];
00716         tmp.element[run1] = new Asin( element[run1]->clone() ); 
00717     }
00718     return tmp;
00719 }
00720 
00721 Expression Expression::getAcos( ) const{
00722 
00723         Expression tmp("", nRows, nCols);
00724     uint run1;
00725 
00726     for( run1 = 0; run1 < dim; run1++ ){
00727         delete tmp.element[run1];
00728         tmp.element[run1] = new Acos( element[run1]->clone() ); 
00729     }
00730     return tmp;
00731 }
00732 
00733 Expression Expression::getAtan( ) const{
00734 
00735         Expression tmp("", nRows, nCols);
00736     uint run1;
00737 
00738     for( run1 = 0; run1 < dim; run1++ ){
00739         delete tmp.element[run1];
00740         tmp.element[run1] = new Atan( element[run1]->clone() ); 
00741     }
00742     return tmp;
00743 }
00744 
00745 Expression Expression::getExp( ) const{
00746 
00747         Expression tmp("", nRows, nCols);
00748     uint run1;
00749 
00750     for( run1 = 0; run1 < dim; run1++ ){
00751         delete tmp.element[run1];
00752         tmp.element[run1] = new Exp( element[run1]->clone() ); 
00753     }
00754     return tmp;
00755 }
00756 
00757 Expression Expression::getSqrt( ) const{
00758 
00759         Expression tmp("", nRows, nCols);
00760     uint run1;
00761 
00762     for( run1 = 0; run1 < dim; run1++ ){
00763         delete tmp.element[run1];
00764         tmp.element[run1] = new Power( element[run1]->clone(), new DoubleConstant( 0.5, NE_NEITHER_ONE_NOR_ZERO ) ); 
00765     }
00766     return tmp;
00767 }
00768 
00769 Expression Expression::getLn( ) const{
00770 
00771         Expression tmp("", nRows, nCols);
00772     uint run1;
00773 
00774     for( run1 = 0; run1 < dim; run1++ ){
00775         delete tmp.element[run1];
00776         tmp.element[run1] = new Logarithm( element[run1]->clone() ); 
00777     }
00778     return tmp;
00779 }
00780 
00781 
00782 Expression Expression::getPow( const Expression& arg ) const{
00783 
00784     ASSERT( arg.getDim() == 1 );
00785 
00786     Expression tmp("", nRows, nCols);
00787     uint run1;
00788 
00789     for( run1 = 0; run1 < dim; run1++ ){
00790         delete tmp.element[run1];
00791         tmp.element[run1] = new Power( element[run1]->clone(), arg.element[0]->clone() ); 
00792     }
00793     return tmp;
00794 }
00795 
00796 
00797 Expression Expression::getPowInt( const int &arg ) const{
00798 
00799         Expression tmp("", nRows, nCols);
00800     uint run1;
00801 
00802     for( run1 = 0; run1 < dim; run1++ ){
00803         delete tmp.element[run1];
00804         tmp.element[run1] = new Power_Int( element[run1]->clone(), arg );
00805     }
00806     return tmp;
00807 }
00808 
00809 
00810 
00811 Expression Expression::transpose( ) const{
00812 
00813         Expression tmp("", getNumCols(), getNumRows());
00814 
00815     uint run1, run2;
00816 
00817     for( run1 = 0; run1 < getNumRows(); run1++ ){
00818         for( run2 = 0; run2 < getNumCols(); run2++ ){
00819              delete tmp.element[run2*getNumRows()+run1];
00820              tmp.element[run2*getNumRows()+run1] = element[run1*getNumCols()+run2]->clone();
00821         }
00822     }
00823     return tmp;
00824 }
00825 
00826 
00827 
00828 Expression Expression::getSumSquare( ) const{
00829 
00830     uint run1;
00831 
00832     Expression result = operator*(transpose(), *this);
00833 
00834     CurvatureType c = CT_CONSTANT;
00835     CurvatureType cc;
00836 
00837     for( run1 = 0; run1 < getDim(); run1++ ){
00838 
00839         cc = element[run1]->getCurvature();
00840 
00841         if( cc != CT_CONSTANT ){
00842             if( cc == CT_AFFINE &&
00843                 (c == CT_CONSTANT || c == CT_AFFINE) )
00844                  c = CT_AFFINE;
00845             else c = CT_NEITHER_CONVEX_NOR_CONCAVE;
00846         }
00847     }
00848 
00849     if( c == CT_CONSTANT )
00850         result.element[0]->setCurvature(CT_CONSTANT);
00851 
00852     if( c == CT_AFFINE )
00853         result.element[0]->setCurvature(CT_CONVEX);
00854 
00855     return result;
00856 
00857 }
00858 
00859 
00860 Expression Expression::getLogSumExp( ) const{
00861 
00862     uint run1;
00863 
00864     Expression result;
00865     result = exp( operator()(0) );
00866 
00867     for( run1 = 1; run1 < getDim(); run1++ )
00868         result = result + exp( operator()(run1) );
00869 
00870     result = ln( result );
00871 
00872     CurvatureType c = CT_CONSTANT;
00873     CurvatureType cc;
00874 
00875     for( run1 = 0; run1 < getDim(); run1++ ){
00876 
00877         cc = element[run1]->getCurvature();
00878 
00879         if( cc != CT_CONSTANT ){
00880             if( cc == CT_AFFINE ){
00881                 if( c == CT_CONSTANT || c == CT_AFFINE )
00882                     c = CT_AFFINE;
00883                 else{
00884                     if( c == CT_CONVEX ){
00885                         c = CT_CONVEX;
00886                     }
00887                     else{
00888                         c = CT_NEITHER_CONVEX_NOR_CONCAVE;
00889                     }
00890                 }
00891             }
00892             else{
00893                 if( cc == CT_CONVEX &&
00894                     (c == CT_CONSTANT || c == CT_AFFINE || c == CT_CONVEX ) )
00895                      c = CT_CONVEX;
00896                 else
00897                      c = CT_NEITHER_CONVEX_NOR_CONCAVE;
00898             }
00899         }
00900     }
00901 
00902     if( c == CT_CONSTANT )
00903         result.element[0]->setCurvature(CT_CONSTANT);
00904 
00905     if( c == CT_AFFINE || c == CT_CONVEX )
00906         result.element[0]->setCurvature(CT_CONVEX);
00907 
00908     return result;
00909 }
00910 
00911 
00912 Expression Expression::getEuclideanNorm( ) const{
00913 
00914     uint run1;
00915 
00916     Expression result = sqrt( getSumSquare() );
00917 
00918     CurvatureType c = CT_CONSTANT;
00919     CurvatureType cc;
00920 
00921     for( run1 = 0; run1 < getDim(); run1++ ){
00922 
00923         cc = element[run1]->getCurvature();
00924 
00925         if( cc != CT_CONSTANT ){
00926             if( cc == CT_AFFINE &&
00927                 (c == CT_CONSTANT || c == CT_AFFINE) )
00928                  c = CT_AFFINE;
00929             else c = CT_NEITHER_CONVEX_NOR_CONCAVE;
00930         }
00931     }
00932 
00933     if( c == CT_CONSTANT )
00934         result.element[0]->setCurvature(CT_CONSTANT);
00935 
00936     if( c == CT_AFFINE )
00937         result.element[0]->setCurvature(CT_CONVEX);
00938 
00939     return result;
00940 }
00941 
00942 
00943 
00944 
00945 Expression Expression::getEntropy( ) const{
00946 
00947     ACADOWARNING( RET_NOT_IMPLEMENTED_YET );
00948     return *this;
00949 }
00950 
00951 
00952 Expression Expression::getDot( ) const{
00953 
00954     if( variableType != VT_DIFFERENTIAL_STATE ){
00955 
00956         ACADOERROR( RET_INVALID_ARGUMENTS );
00957         ASSERT( 1 == 0 );
00958     }
00959 
00960     Expression tmp("", getNumRows(), getNumCols(), VT_DDIFFERENTIAL_STATE, component );
00961     return tmp;
00962 }
00963 
00964 
00965 Expression Expression::getNext( ) const{
00966 
00967     return getDot();
00968 }
00969 
00970 
00971 Expression Expression::ADforward ( const Expression &arg ) const{
00972 
00973     ASSERT(     getNumCols() == 1 );
00974     ASSERT( arg.getNumCols() == 1 );
00975 
00976         Expression result("", getNumRows(), arg.getNumRows());
00977 
00978     uint run1, run2;
00979 
00980     Expression seed( arg.getNumRows() );
00981 
00982     for( run1 = 0; run1 < arg.getNumRows(); run1++ ){
00983 
00984         delete seed.element[run1];
00985         seed.element[run1] = new DoubleConstant( 1.0, NE_ONE );
00986 
00987         Expression tmp = ADforward( arg, seed );
00988 
00989         delete seed.element[run1];
00990         seed.element[run1] = new DoubleConstant( 0.0, NE_ZERO );
00991 
00992         for( run2 = 0; run2 < getNumRows(); run2++ ){
00993             delete result.element[run2*arg.getNumRows()+run1];
00994             result.element[run2*arg.getNumRows()+run1] = tmp.element[run2]->clone();
00995         }
00996     }
00997 
00998     return result;
00999 }
01000 
01001 
01002 Expression Expression::ADforward ( const VariableType &varType_, const int *arg, int nV ) const{
01003  
01004     ASSERT( getNumCols() == 1 );
01005 
01006         Expression result("", (int) getNumRows(), nV);
01007 
01008     int run1, run2;
01009 
01010     IntermediateState seed( nV );
01011 
01012     for( run1 = 0; run1 < nV; run1++ ) seed(run1) = 0;
01013 
01014     for( run1 = 0; run1 < nV; run1++ ){
01015 
01016         seed(run1) = 1.0;
01017 
01018         Expression tmp = ADforward( varType_, arg, seed );
01019 
01020         seed(run1) = 0.0;
01021 
01022         for( run2 = 0; run2 < (int) getNumRows(); run2++ ){
01023             delete result.element[run2*nV+run1];
01024             result.element[run2*nV+run1] = tmp.element[run2]->clone();
01025         }
01026     }
01027     return result;
01028 }
01029 
01030 
01031 
01032 Expression Expression::ADbackward ( const Expression &arg ) const{
01033 
01034     ASSERT(     getNumCols() == 1 );
01035     ASSERT( arg.getNumCols() == 1 );
01036 
01037         Expression result("", getNumRows(), arg.getNumRows());
01038 
01039     uint run1, run2;
01040 
01041     Expression seed( getNumRows() );
01042 
01043     for( run1 = 0; run1 < getNumRows(); run1++ ){
01044 
01045         delete seed.element[run1];
01046         seed.element[run1] = new DoubleConstant( 1.0, NE_ONE );
01047 
01048         Expression tmp = ADbackward( arg, seed );
01049 
01050         delete seed.element[run1];
01051         seed.element[run1] = new DoubleConstant( 0.0, NE_ZERO );
01052 
01053         for( run2 = 0; run2 < arg.getNumRows(); run2++ ){
01054             delete result.element[run1*arg.getNumRows()+run2];
01055             result.element[run1*arg.getNumRows()+run2] = tmp.element[run2]->clone();
01056         }
01057     }
01058 
01059     return result;
01060 }
01061 
01062 
01063 Expression Expression::ADforward ( const Expression &arg, const Expression &seed ) const{
01064 
01065     unsigned int run1;
01066     const unsigned int n = arg.getDim();
01067 
01068     ASSERT( arg .isVariable() == BT_TRUE );
01069     ASSERT( seed.getDim    () == n       );
01070 
01071     VariableType *varType   = new VariableType[n];
01072     int          *Component = new int[n];
01073 
01074     for( run1 = 0; run1 < n; run1++ ){
01075         arg.element[run1]->isVariable(varType[run1],Component[run1]);
01076     }
01077 
01078         Expression result = ADforward( varType, Component, seed );
01079     delete[] Component;
01080     delete[] varType;
01081 
01082     return result;
01083 }
01084 
01085 
01086 Expression Expression::ADforward (  const VariableType &varType_,
01087                                                                            const int          *arg     ,
01088                                                                            const Expression   &seed      ) const{
01089 
01090         VariableType *varType = new VariableType[seed.getDim()];
01091         for( uint run1 = 0; run1 < seed.getDim(); run1++ ) varType[run1] = varType_;
01092         Expression tmp = ADforward( varType, arg, seed );
01093         delete[] varType;
01094         return tmp;
01095 }
01096 
01097 
01098 Expression Expression::ADforward ( const VariableType *varType_,
01099                                                                    const int          *arg     ,
01100                                                                    const Expression   &seed      ) const{
01101 
01102     unsigned int run1, run2;
01103     const unsigned int n = seed.getDim();
01104 
01105         Expression result("", getNumRows(), getNumCols());
01106 
01107     VariableType  *varType   = new VariableType[n];
01108     int           *Component = new int         [n];
01109     Operator     **seed1     = new Operator*   [n];
01110 
01111     for( run1 = 0; run1 < n; run1++ ){
01112         varType  [run1] = varType_[run1];
01113         Component[run1] = arg[run1];
01114         seed1    [run1] = seed.element[run1]->clone();
01115     }
01116 
01117     for( run1 = 0; run1 < getDim(); run1++ ){
01118 
01119         delete result.element[run1];
01120 
01121         int Dim = n;
01122         int nIS = 0;
01123         TreeProjection **IS = 0;
01124 
01125         element[run1]->initDerivative();
01126         result.element[run1] = element[run1]->AD_forward( Dim, varType, Component, seed1, nIS, &IS );
01127 
01128         for( run2 = 0; (int) run2 < nIS; run2++ ){
01129             if( IS[run2] != 0 ){
01130                 delete IS[run2];
01131             }
01132         }
01133         if( IS != 0 )
01134             free(IS);
01135     }
01136 
01137     delete[] varType  ;
01138     delete[] Component;
01139 
01140     for( run1 = 0; run1 < n; run1++ )
01141         delete seed1[run1];
01142 
01143     delete[] seed1;
01144 
01145     return result;
01146 }
01147 
01148 
01149 Expression Expression::getODEexpansion( const int &order, const int *arg ) const{
01150  
01151         IntermediateState coeff("", (int) dim, order+2 );
01152         
01153     VariableType  *vType = new VariableType[dim*(order+1)+1];
01154     int           *Comp  = new int         [dim*(order+1)+1];
01155     Operator     **seed  = new Operator*   [dim*(order+1)+1];
01156         
01157         Operator **der = new Operator*[dim*(order+1)];
01158         
01159         vType[0] = VT_TIME;
01160         Comp [0] = 0      ;
01161         seed [0] = new DoubleConstant( 1.0 , NE_ONE );
01162         
01163         for( uint i=0; i<dim; i++ ){
01164                 coeff(i,0)   = Expression("",1,1,VT_DIFFERENTIAL_STATE,arg[i]);
01165                 coeff(i,1)   = operator()(i);
01166                 vType[i+1]   = VT_DIFFERENTIAL_STATE;
01167                 Comp [i+1]   = arg[i];
01168                 seed [i+1]   = coeff.element[(order+2)*i+1];
01169                 der[i]       = element[i]->clone();
01170         }
01171         
01172         int nIS = 0;
01173         TreeProjection **IS = 0;
01174         
01175         for( int j=0; j<order; j++ ){
01176                 for( uint i=0; i<dim; i++ ){
01177                         der[dim*j+i]->initDerivative();
01178                         der[dim*(j+1)+i] = der[dim*j+i]->AD_forward( (j+1)*dim+1, vType, Comp, seed, nIS, &IS );
01179                 }
01180                 for( uint i=0; i<dim; i++ ){
01181                         coeff(i,j+2) = *der[dim*(j+1)+i];
01182                         vType[dim*(j+1)+i+1] = VT_INTERMEDIATE_STATE;
01183                         Comp [dim*(j+1)+i+1] = coeff.element[(order+2)*i+j+1]->getGlobalIndex();
01184                         seed [dim*(j+1)+i+1] = coeff.element[(order+2)*i+j+2];
01185                 }
01186         }
01187         
01188         for( int run = 0; run < nIS; run++ ){
01189                 
01190                 if( IS[run] != 0 ) delete IS[run];
01191         }
01192         if( IS != 0 ) free(IS);
01193         
01194         delete[] vType;
01195         delete[] Comp;
01196         for( uint i=0; i<dim*(order+1); i++ ) delete der[i];
01197         delete[] der;
01198         delete seed[0];
01199         delete[] seed;
01200         
01201         return coeff;
01202 }
01203 
01204 
01205 
01206 Expression Expression::ADbackward( const Expression &arg, const Expression &seed ) const{
01207 
01208     int Dim = arg.getDim();
01209     int run1, run2;
01210 
01211     ASSERT( arg .isVariable() == BT_TRUE  );
01212     ASSERT( seed.getDim    () == getDim() );
01213 
01214         Expression result("", arg.getNumRows(), arg.getNumCols());
01215 
01216     VariableType *varType   = new VariableType[Dim];
01217     int          *Component = new int         [Dim];
01218     Operator    **iresult   = new Operator*   [Dim];
01219 
01220     for( run1 = 0; run1 < Dim; run1++ ){
01221         arg.element[run1]->isVariable(varType[run1],Component[run1]);
01222     }
01223 
01224         int nIS = 0;
01225         TreeProjection **IS = 0;
01226         
01227     for( run1 = 0; run1 < (int) getDim(); run1++ ){
01228 
01229         Operator *seed1 = seed.element[run1]->clone();
01230 
01231         for( run2 = 0; run2 < Dim; run2++ )
01232              iresult[run2] = new DoubleConstant(0.0,NE_ZERO);
01233 
01234         element[run1]->initDerivative();
01235         element[run1]->AD_backward( Dim, varType, Component, seed1, iresult, nIS, &IS );
01236 
01237         for( run2 = 0; run2 < Dim; run2++ ){
01238             Operator *sum = result.element[run2]->clone();
01239             delete result.element[run2];
01240             result.element[run2] = new Addition( sum->clone(), iresult[run2]->clone() );
01241             delete sum;
01242             delete iresult[run2];
01243         }
01244     }
01245 
01246 
01247         for( int run = 0; run < nIS; run++ ){
01248                 
01249                 if( IS[run] != 0 ) delete IS[run];
01250         }
01251         if( IS != 0 ) free(IS);
01252     
01253     delete[] iresult   ;
01254     delete[] varType   ;
01255     delete[] Component ;
01256 
01257     return result;
01258 }
01259 
01260 Expression Expression::ADsymmetric(     const Expression &arg, 
01261                                         const Expression &S  , 
01262                                         const Expression &l  , 
01263                                         Expression *dfS,    
01264                                         Expression *ldf    ) const{
01265 
01266         int Dim = arg.getDim();
01267         ASSERT( (int) S.getNumRows() == Dim      );
01268         ASSERT( (int) S.getNumRows() == (int) S.getNumCols() );
01269         
01270         IntermediateState H    = ADsymmetric( arg, l, dfS, ldf );
01271         IntermediateState sum  = zeros<double>(Dim,Dim);
01272         IntermediateState sum2 = zeros<double>(Dim,Dim);
01273         
01274         int i,j,k;
01275         for( i=0; i<Dim; i++ )
01276                 for( j=0; j<Dim; j++ ){
01277                         for( k=0; k<=i; k++ ){
01278                                 sum(i,j) += H(i,k)*S(k,j);
01279                         }
01280                         for( k=i+1; k<Dim; k++ ){
01281                                 sum(i,j) += H(k,i)*S(k,j);
01282                         }
01283                 }
01284 
01285         for( i=0; i<Dim; i++ )
01286                 for( j=0; j<=i; j++ )
01287                         for( k=0; k<Dim; k++ )
01288                                 sum2(i,j) += S(k,i)*sum(k,j);
01289 
01290         for( i=0; i<Dim; i++ )
01291                 for( j=0; j<i; j++ )
01292                         sum2(j,i) = sum2(i,j);
01293 
01294         if( dfS != 0 ) {
01295                 // multiplication with the proper forward seeds:
01296                 *dfS = *dfS*S;
01297         }
01298           
01299         return sum2;
01300 }
01301 
01302 
01303 Expression Expression::ADsymmetric(     const Expression &arg, 
01304                                         const Expression &l  , 
01305                                         Expression *dfS,    
01306                                         Expression *ldf    ) const{
01307 
01308         int Dim = arg.getDim();
01309         int nS  = Dim;
01310         int run1, run2;
01311         
01312         Expression S = eye<double>(Dim);
01313         
01314         ASSERT( arg .isVariable()    == BT_TRUE  );
01315         ASSERT( l.getDim()           == getDim() );
01316         
01317         Expression result( nS, nS );
01318 
01319         VariableType *varType   = new VariableType[Dim];
01320         int          *Component = new int         [Dim];
01321         Operator    **dS        = new Operator*   [nS];
01322         Operator    **ld        = new Operator*   [Dim];
01323         Operator    **H         = new Operator*   [nS*nS];
01324 
01325         for( run1 = 0; run1 < Dim; run1++ ){
01326                 arg.element[run1]->isVariable(varType[run1],Component[run1]);
01327         }
01328 
01329         int nLIS = 0;
01330         int nSIS = 0;
01331         int nHIS = 0;
01332         TreeProjection **LIS = 0;
01333         TreeProjection **SIS = 0;
01334         TreeProjection **HIS = 0;
01335 
01336         Expression tmp((int) getDim(),Dim);
01337         Expression tmp2(Dim);
01338         
01339         for( run1 = 0; run1 < (int) getDim(); run1++ ){
01340 
01341                 Operator *l1 = l.element[run1]->clone();
01342                 Operator **S1 = new Operator*[Dim*nS];
01343 
01344                 for( run2 = 0; run2 < Dim*nS; run2++ )
01345                         S1[run2] = S.element[run2]->clone();
01346 
01347                 for( run2 = 0; run2 < nS; run2++ )
01348                         dS[run2] = new DoubleConstant(0.0,NE_ZERO);
01349 
01350                 for( run2 = 0; run2 < Dim; run2++ )
01351                         ld[run2] = new DoubleConstant(0.0,NE_ZERO);
01352 
01353                 for( run2 = 0; run2 < nS*nS; run2++ )
01354                         H[run2] = new DoubleConstant(0.0,NE_ZERO);
01355 
01356                 element[run1]->initDerivative();
01357                 element[run1]->AD_symmetric( Dim, varType, Component, l1, S1, nS, dS, ld, H, nLIS, &LIS, nSIS, &SIS, nHIS, &HIS );
01358 
01359                 int run3 = 0;
01360 
01361                 for( run2 = 0; run2 < nS; run2++ ){
01362                         for( run3 = 0; run3 < run2; run3++ ){
01363                                 Operator *sum = result.element[run2*nS+run3]->clone();
01364                                 delete result.element[run2*nS+run3];
01365                                 delete result.element[run3*nS+run2];
01366                                 result.element[run2*nS+run3] = sum->myAdd( sum, H[run2*nS+run3] );
01367                                 result.element[run3*nS+run2] = sum->myAdd( sum, H[run2*nS+run3] );
01368                                 delete sum;
01369                         }
01370                         Operator *sum = result.element[run2*nS+run2]->clone();
01371                         delete result.element[run2*nS+run2];
01372                         result.element[run2*nS+run2] = sum->myAdd( sum, H[run2*nS+run2] );
01373                         delete sum;
01374                 }
01375 
01376                 if( dfS != 0 ){
01377                         for( run2 = 0; run2 < nS; run2++ ){
01378                                 delete tmp.element[run1*nS+run2];
01379                                 tmp.element[run1*nS+run2] = dS[run2]->clone();
01380                         }
01381                 }
01382 
01383                 if( ldf != 0 ){
01384                    for( run2 = 0; run2 < nS; run2++ ){
01385                         Operator *sum = tmp2.element[run2]->clone();
01386                         delete tmp2.element[run2];
01387                         tmp2.element[run2] = sum->myAdd( sum, ld[run2] );
01388                         delete sum;
01389                   }
01390                 }
01391 
01392                 for( run2 = 0; run2 < Dim*nS; run2++ )
01393                         delete S1[run2];
01394 
01395                 for( run2 = 0; run2 < nS; run2++ )
01396                         delete dS[run2];
01397 
01398                 for( run2 = 0; run2 < Dim; run2++ )
01399                         delete ld[run2];
01400 
01401                 for( run2 = 0; run2 < nS*nS; run2++ )
01402                         delete H[run2];
01403 
01404                 delete[] S1;
01405         }
01406 
01407         if( dfS != 0 ) *dfS = tmp ;
01408         if( ldf != 0 ) *ldf = tmp2;
01409         
01410         
01411         for( int run = 0; run < nLIS; run++ ){
01412                 if( LIS[run] != 0 ) delete LIS[run];
01413         }
01414         if( LIS != 0 ) free(LIS);
01415 
01416         for( int run = 0; run < nSIS; run++ ){
01417                 if( SIS[run] != 0 ) delete SIS[run];
01418         }
01419         if( SIS != 0 ) free(SIS);
01420 
01421         for( int run = 0; run < nHIS; run++ ){
01422                 if( HIS[run] != 0 ) delete HIS[run];
01423         }
01424         if( HIS != 0 ) free(HIS);
01425 
01426         delete[] dS        ;
01427         delete[] ld        ;
01428         delete[] H         ;
01429         delete[] varType   ;
01430         delete[] Component ;
01431 
01432         return result;
01433 }
01434 
01435 
01436 returnValue Expression::substitute( int idx, const Expression &arg ) const{
01437 
01438     ASSERT( arg.getDim() == 1 );
01439 
01440     uint i;
01441     for( i = 0; i < getDim(); i++ )
01442         if( element[i] != 0 ) element[i]->substitute( idx, arg.element[0] );
01443 
01444     return SUCCESSFUL_RETURN;
01445 }
01446 
01447 
01448 Expression Expression::operator-() const{
01449 
01450     uint run1, run2;
01451         Expression tmp("", getNumRows(), getNumCols());
01452 
01453     for( run1 = 0; run1 < getNumRows(); run1++ ){
01454         for( run2 = 0; run2 < getNumCols(); run2++ ){
01455              delete tmp.element[run1*getNumCols()+run2];
01456              tmp.element[run1*getNumCols()+run2] = new Subtraction( new DoubleConstant(0.0,NE_ZERO),
01457                                                                     element[run1*getNumCols()+run2]->clone() );
01458         }
01459     }
01460     return tmp;
01461 }
01462 
01463 
01464 
01465 
01466 
01467 //
01468 // PROTECTED MEMBER FUNCTIONS:
01469 //
01470 
01471 void Expression::construct( VariableType variableType_, uint globalTypeID, uint nRows_, uint nCols_, const std::string &name_ ){
01472 
01473     nRows        = nRows_       ;
01474     nCols        = nCols_       ;
01475     dim          = nRows*nCols  ;
01476     variableType = variableType_;
01477     component    = globalTypeID ;
01478     name         = name_        ;
01479 
01480     uint i;
01481     element = (Operator**)calloc(dim,sizeof(Operator*));
01482 
01483     for( i = 0; i < dim; i++ ){
01484 
01485         switch( variableType ){
01486             case VT_UNKNOWN           : element[i] = new DoubleConstant( 0.0, NE_ZERO ); break;
01487             case VT_INTERMEDIATE_STATE:
01488                                         element[i] = new TreeProjection( "" );
01489                                         break;
01490             default                   : element[i] = new Projection( variableType_, globalTypeID+i, "" ); break;
01491         }
01492     }
01493 }
01494 
01495 
01496 void Expression::copy(const Expression &rhs)
01497 {
01498         nRows = rhs.nRows;
01499         nCols = rhs.nCols;
01500         dim = rhs.dim;
01501         variableType = rhs.variableType;
01502         component = rhs.component;
01503 
01504         uint i;
01505         element = (Operator**) calloc(dim, sizeof(Operator*));
01506 
01507         for (i = 0; i < dim; i++)
01508         {
01509                 if (rhs.element[i] != 0)
01510                         element[i] = rhs.element[i]->clone();
01511                 else
01512                         element[i] = 0;
01513         }
01514 
01515         // Name not copied?
01516 }
01517 
01518 Expression& Expression::assignmentSetup(const Expression &arg)
01519 {
01520         deleteAll();
01521         nRows = arg.getNumRows();
01522         nCols = arg.getNumCols();
01523         dim = nRows * nCols;
01524         variableType = VT_INTERMEDIATE_STATE;
01525 
01526         element = (Operator**) calloc(dim, sizeof(Operator*));
01527 
01528         VariableType tt = VT_UNKNOWN;
01529         int comp = 0;
01530 
01531         for (uint i = 0; i < dim; i++)
01532         {
01533                 arg.element[i]->isVariable(tt, comp);
01534                 if (tt == VT_INTERMEDIATE_STATE)
01535                         element[i] = arg.element[i]->clone();
01536                 else
01537                 {
01538                         std::stringstream tmpName;
01539                         if (name.empty() == false)
01540                         {
01541                                 if (dim > 1)
01542                                         tmpName << name << "[" << i << "]";
01543                                 else
01544                                         tmpName << name;
01545                         }
01546                         TreeProjection tmp(tmpName.str());
01547                         tmp.operator=(*(arg.element[i]));
01548                         element[i] = tmp.clone();
01549                 }
01550         }
01551         return *this;
01552 }
01553 
01554 
01555 void Expression::deleteAll( ){
01556 
01557     uint i;
01558 
01559     for( i = 0; i < dim; i++ )
01560         if( element[i] != 0 )
01561             delete element[i];
01562 
01563     if( element != 0 ) free(element);
01564 }
01565 
01566 
01567 
01568 
01569 BooleanType Expression::isDependingOn( VariableType type ) const{
01570 
01571     uint run1;
01572 
01573     for( run1 = 0; run1 < getDim(); run1++ )
01574         if( element[run1]->isDependingOn(type) == BT_TRUE )
01575             return BT_TRUE;
01576 
01577     return BT_FALSE;
01578 }
01579 
01580 
01581 
01582 BooleanType Expression::isDependingOn( const Expression &e ) const{
01583     ASSERT( e.getDim() ==1 );
01584     DVector sum=getDependencyPattern(e).sumRow();
01585 
01586     if( fabs(sum(0) - EPS) > 0 ) return BT_TRUE;
01587     return BT_FALSE;
01588 }
01589 
01590 Operator* Expression::getOperatorClone( uint idx ) const{
01591 
01592     ASSERT( idx < getDim() );
01593 
01594     Operator *tmp = element[idx]->passArgument();
01595     if( tmp == 0 ) tmp = element[idx];
01596 
01597     return tmp->clone();
01598 }
01599 
01600 
01601 CLOSE_NAMESPACE_ACADO
01602 
01603 // end of file.


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