--- Charlls Quarra <charlls_quarra@yahoo.com.ar> escribió:
Hi,
I wanted to define a simple matrix of derivatives (i.e: a Jacobian) first as an indexed object:
ex f = a*x + b*y + c*x*y + ... ex g = u*pow(x,2) + v*pow(y,2) + w*x,y + ....
Ok i'm happy to have something working that seems to be ok for my purposes now (the code to make this example work is at the bottom of the mail). I wish to eventually make this mindex_array a fully fledged ginac object, but i guess ill take small steps: mindex_array t3_idx( 3 , 2 ); // 3 is the number of indices, 2 the range mindex_array tink = symbolic_vector( 2 , "X" ); for (int i=0; i< 2; i++) cout << tink(i) << endl; ex fleet = sin(x*tink(0))*3 + pow( tink(1)/tink(0) , 3 ); cout << fleet.diff( ex_to<symbol>(tink(0)) ) << " ... " << fleet.diff( ex_to<symbol>(tink(1)) ) << endl; for (int qi=0 ; qi < 2 ; qi++) { for (int qj = 0; qj < 2 ; qj ++) { for (int qk = 0; qk < 2 ; qk++) { t3_idx( qi , qj , qk ) = fleet.diff( ex_to<symbol>(tink(qi)) ).diff( ex_to<symbol>(tink(qj)) ).diff( ex_to<symbol>(tink(qk) ) ); } } }; for (int qi=0 ; qi < 2 ; qi++) { for (int qj=0 ; qj < 2 ; qj++) { for (int qk=0 ; qk < 2 ; qk++) { cout << " t( " << qi <<" , " << qj << " , " << qk << " ) = " << t3_idx( qi , qj , qk ) << endl; } } } //////////////////////////////// o //////////////////////// //needed code is here! inline matrix& open_mt( ex& m ) { return ( *(matrix*) & ex_to<matrix>( m ) ); }; template <typename T> std::string stringify(const T& t) { std::ostringstream oss; oss.operator << (t); return oss.str(); } matrix multi_index_matrix( int n , int range ) { if (n == 1) { return matrix( 1 , range ); //return; } if (n == 2) { return matrix( range , range ); //return; } if (n % 2 == 1) { // *put = matrix(1 , range); matrix put( 1 , range ); for (int i=0; i < range; i++) { put( 0 , i ) = multi_index_matrix( n-1 , range ); } return put; } matrix put( range , range ); for (int i=0; i< range ; i++) { for (int j=0; j < range; j++) { put( i , j ) = multi_index_matrix( n - 2 , range ); } } return put; }; struct mindex_array { int nb_index; int nb_index_parity; int range; int* index_placeholder_array; matrix mb_array; mindex_array( int nb_indices , int range_ ) : nb_index( nb_indices ) , nb_index_parity( nb_index % 2 ) , range( range_ ) { index_placeholder_array = (int*)malloc( nb_indices * sizeof(int) ); mb_array = multi_index_matrix( nb_indices , range ); } // mindex_array& operator =(const ex& e) //{ //} ex& operator() (int i) { return mb_array( 0 , i ); } ex& operator() (int i , int j) { return mb_array(i , j); } ex& operator() (int i , int j , int k) { return open_mt( mb_array(0 , i) )( j , k ); } ex& operator() (int i , int j , int k , int m) { return open_mt( mb_array(i , j) )( k , m ); } ex& operator() (int i , int j , int k , int m , int n) { return open_mt( open_mt( mb_array(0 , i) )( j , k ) )( m , n ); } ex& operator() (int i , int j , int k , int m , int n , int p) { return open_mt( open_mt( mb_array(i , j) )( k , m ) )( n , p ); } ex& operator () () { matrix foo; foo = mb_array; if ( nb_index_parity == 1 ) { foo = ex_to<matrix>( foo(0, index_placeholder_array[0] ) ); } for (int i=nb_index_parity; i< nb_index-2 ; i+=2) { foo = ex_to<matrix>( foo( index_placeholder_array[ i ] , index_placeholder_array[ i+1 ] ) ); }; return foo( index_placeholder_array[ nb_index - 2 ] , index_placeholder_array[ nb_index - 1 ] ); } }; mindex_array symbolic_vector( int range , const string& base_name ) { mindex_array foo( 1 , range ); string full_name;// = base_name; for (int i=0; i<range ; i++) { full_name = base_name + "_"; full_name += stringify<int>(i); foo( 0 , i ) = symbol( full_name ); } return foo; }; __________________________________________________ Preguntá. Respondé. Descubrí. Todo lo que querías saber, y lo que ni imaginabas, está en Yahoo! Respuestas (Beta). ¡Probalo ya! http://www.yahoo.com.ar/respuestas