[GiNaC-list] multi index storage of expressions

[Charlls Quarra]

--- Charlls Quarra <charlls_quarra at 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