Re: Help with Differentiation

Markus

Thank you for the reply. The new examples help a lot to understand 
how to use it.

I would like to be able to have expressions with constants e.g.

x + y + 1

but that will not compile.  What I have to do is to add a variable z

x + y + z

and give z a value, but then it has increased the size of the 
problem to 3.  I then found that I needed a third function and when I 
put in just z it would not compile, as that is a variable, not an 
expression.

Ideally, I would like to have a set of variables which are parameters 
and not differentiated. That could be done by changing the interior 
of the Solver object, which I have started to do.

I attach an example which I have built which solves a two 
dimensional newton problem with two variables and one parameter.

Just playing at the moment. I can see a lot of uses for this.

Daixtrose Versions
==============
So far I have only run Solver with the old 0.0.2 daixtrose which I 
use on one system which has the old boost.  That worked fine.

I will compare your new one with what I have done, which checks 
the version of boost and will run with either.  It also checks the 
version of Loki as well, as they have changed the name of the 
TYPELIST macros and also have a template version which 
replaces the macros.

I will then report back on that.

Thanks

John

// mysolver_2.cc adapted from Solver_2.C
//
// This has two variables x and y and two functions. These are differentiated wrt x and  y and
// the resulting jacobian problem is solved.

#include "daixtrose/Daixt.h"
#include "tiny/TinyMatAndVec.h" 

// FIXIT: find out why wstring did not work on some gcc
#define DISABLE_WIDE_CHAR_SUPPORT
#include <boost/lexical_cast.hpp>
#include <boost/mpl/int.hpp>

#include <iostream>
#include <string>
#include <vector>

// code for ML response December 2005

////////////////////////////////////////////////////////////////////////////////
// ************************************************************************** //
////////////////////////////////////////////////////////////////////////////////

// A named and compile-time-numbered variable which can be used for
// differentiation. 

struct DisambiguatedVariable {};

template <std::size_t Number> 
class Variable
{
public:
	inline Variable() {} // intel's C++-70 needs this one
	typedef DisambiguatedVariable Disambiguation;
	std::string GetName() const;// { return "undefined"; }
};


////////////////////////////////////////////////////////////////////////////
// convenience macros

#define DECLARE_VARIABLE(NAME, NUMBER, LATEX_NAME)	\
	template <>										\
	std::string Variable<NUMBER>::GetName() const	\
	{												\
		return LATEX_NAME;							\
	}												\
													\
	Variable<NUMBER> NAME

#define EXTERN_DECLARE_VARIABLE(NAME, NUMBER)	\
	extern Variable<NUMBER> NAME


////////////////////////////////////////////////////////////////////////////
// access to the expression name

template <std::size_t Number>
inline std::string GetName(const Variable<Number> & v) 
{ 
	return v.GetName(); 
} 
   
template <class ARG, class OP> 
inline std::string GetName(const Daixt::UnOp<ARG, OP>& UO) 
{ 
	return std::string("(") + OP::Symbol() + " " + GetName(UO.arg()) + ")";
} 
   
template <class LHS, class RHS, class OP> 
inline std::string GetName(const Daixt::BinOp<LHS, RHS, OP>& BO) 
{ 
	return 
		"(" + GetName(BO.lhs()) 
		+ " " + OP::Symbol() + " " 
		+ GetName(BO.rhs()) + ")"; 
} 

template<class T> 
inline std::string GetName(const Daixt::Expr<T>& E) 
{ 
	return GetName(E.content()); 
} 

inline std::string 
GetName(Daixt::IsNull<DisambiguatedVariable>) 
{
	return "0";
}
   
inline std::string 
GetName(Daixt::IsOne<DisambiguatedVariable>) 
{
	return "1";
}

inline std::string 
GetName(const Daixt::Scalar<DisambiguatedVariable> & v)
{
	return boost::lexical_cast<std::string>(v.Value());
}

template<typename ARG, int m, int n> 
inline std::string 
GetName(const Daixt::UnOp<ARG, Daixt::DefaultOps::RationalPower<m, n> >& UO)
{
	return "(" + GetName(UO.arg()) + ")^(" 
		+ boost::lexical_cast<std::string>(m) 
		+ "/" 
		+ boost::lexical_cast<std::string>(n) + ")"; 
}



////////////////////////////////////////////////////////////////////////////
// some features for our expressions 

template <int NumberOfVars>
struct AccessibleBaseClass 
{ 
	virtual std::string GetName() const = 0; 
	virtual double 
	GetValue(std::vector<double> const & V,
			 std::size_t index) const = 0; 
	virtual AccessibleBaseClass<NumberOfVars> * DeepCopy() const = 0; 
	virtual ~AccessibleBaseClass() {} 
};

 
// specializations which allows us to plug polymorphic behviour to our 
// expressions

namespace Daixt 
{ 

template <class T, int NumberOfVars> 
class FeaturesOfExpression<AccessibleBaseClass<NumberOfVars>, T> 
	: public AccessibleBaseClass<NumberOfVars> 
{ 
public: 
	std::string GetName() const 
	{ 
		return ::GetName(static_cast<const T&>(*this)); 
	} 
  
	////////////////////////////////////////////////////////////////////////
	// getting values 

private:
	////////////////////////////////////////////////////////////////////////
	// access to the expression values:

	template <std::size_t Number>
	inline double GetValue(const Variable<Number> & v,
						   std::vector<double> const & V,
						   std::size_t index) const 
	{ 
		// assuming memory layout 
		//
		// var_1[0]
		// var_2[0]
		// ...
		// var_N[0]
		// var_1[1]
		// var_2[1]
		// ...
		// var_N[1]
		// etc ...

		// please check this quick hack access formula, it might be wrong!
		// Number starts at index 1, but std::vector at index 0
		std::size_t const myindex = index * NumberOfVars + Number - 1;
#ifndef NDEBUG
		return V.at(myindex);
#else
		return V[myindex];
#endif
	}
   
	template <class ARG, class OP> 
	inline double GetValue(const Daixt::UnOp<ARG, OP>& UO,
						   std::vector<double> const & V,
						   std::size_t index) const 
	{ 
		return OP::Apply(this->GetValue(UO.arg(), V, index), 
						 Daixt::Hint<double>());
	} 
   
	template <class LHS, class RHS, class OP> 
	inline double GetValue(const Daixt::BinOp<LHS, RHS, OP>& BO,
						   std::vector<double> const & V,
						   std::size_t index) const 
	{ 
		return OP::Apply(this->GetValue(BO.lhs(), V, index), 
						 this->GetValue(BO.rhs(), V, index), 
						 Daixt::Hint<double>());
	} 

	template<class TT> 
	inline double GetValue(const Daixt::Expr<TT>& E,
						   std::vector<double> const & V,
						   std::size_t index) const 
	{ 
		return this->GetValue(E.content(), V, index); 
	} 

	inline double 
	GetValue(Daixt::IsNull<DisambiguatedVariable>,
			 std::vector<double> const & V,
			 std::size_t index) const 
	{
		return 0.0;
	}
   
	inline double 
	GetValue(Daixt::IsOne<DisambiguatedVariable>,
			 std::vector<double> const & V,
			 std::size_t index) const 
	{
		return 1.0;
	}

	inline double 
	GetValue(const Daixt::Scalar<DisambiguatedVariable> & v,
			 std::vector<double> const & V,
			 std::size_t index) const
	{
		return v.Value();
	}

public:
	////////////////////////////////////////////////////////////////////////

	double GetValue(std::vector<double> const & V,
					std::size_t index) const
	{
		return this->GetValue(static_cast<const T&>(*this), V, index);
	}

	FeaturesOfExpression<
		AccessibleBaseClass<NumberOfVars>, T>* DeepCopy() const 
	{ 
		return new T(static_cast<const T&>(*this)); 
	} 
  
	virtual ~FeaturesOfExpression() {}; 
}; 

} // namespace Daixt 
   

   

////////////////////////////////////////////////////////////////////////////
// ********************************************************************** //
////////////////////////////////////////////////////////////////////////////

template <int NumberOfVars, 
          typename VectorT, 
          typename MatrixT>
class Solver
{
	VectorT Expressions_;
	MatrixT Jacobian_;

	typedef AccessibleBaseClass<NumberOfVars> Accessor;

	typedef boost::mpl::int_<NumberOfVars> number_of_vars_t;

	// storing values outside the variable objects !!!!
	// this allows sharing with f77 or C legacy code
	std::vector<double> Values_;

public:
	void Resize(std::size_t size)
	{
		Values_.resize(NumberOfVars * size);

		// fill in some useful values, please replace by your own
		// initialize method
		for (int i = 0; i < size; ++i)
		{
			for (int j = 0; j < NumberOfVars; ++j)
			{
				std::cerr << "Values_[" << i * NumberOfVars + j << "] = " 
						  << j + 1.0 << std::endl;
				Values_[i * NumberOfVars + j] = j + 1.0;
			}
		}
	}
	
private:
	template<int Row, int Col = NumberOfVars>
	struct StoreDiff
	{
		template<typename T>
		static inline void 
		Apply(Daixt::Expr<T> const & e, MatrixT & Jacobian)
		{
			using namespace Daixt::ExprManip;
			using namespace Daixt::Differentiation;

			Variable<Col> Differ;

      
			Jacobian(Row, Col) = 
				Daixt::ChangeDisambiguation<Accessor>
				(Simplify(Diff(e, Differ))).DeepCopy(); 
      
			std::cerr << "Stored d(" << GetName(e) 
					  << ") / d(" << Differ.GetName()
					  << ") = " 
					  << Jacobian(Row, Col)->GetName()
					  << " at (" << Row << ", " << Col << ")" << std::endl;
      

			// recursive call
			StoreDiff<Row, Col-1>::Apply(e, Jacobian);
		}
	};

  
	// end of recursion
	template<int Row>
	struct StoreDiff<Row, 1>
	{
		template<typename T>
		static inline void 
		Apply(Daixt::Expr<T> const & e, MatrixT & Jacobian)
		{
			using namespace Daixt::ExprManip;
			using namespace Daixt::Differentiation;

			Variable<1> Differ;

      
			Jacobian(Row, 1) = 
				Daixt::ChangeDisambiguation<Accessor>
				(Simplify(Diff(e, Differ))).DeepCopy(); 

			std::cerr << "Stored d(" << GetName(e) 
					  << ") / d(" << GetName(Differ) 
					  << ") = " 
					  << Jacobian(Row, 1)->GetName()
					  << " at (" << Row << ", " << 1 << ")" << std::endl;


			// no more recursive call
		}
	};
  
public:

	template<int Row, typename T>
	void AddExpression(Daixt::Expr<T> const & e)
	{
		// you must delete this pointer in the destructor
		Expressions_(Row) = 
			Daixt::ChangeDisambiguation<Accessor>(e).
			DeepCopy(); 

		std::cerr << "registrating expression '" 
				  << Expressions_(Row)->GetName() << "'" << std::endl;

		// this stores all d(e)/d(var(i))
		StoreDiff<Row, NumberOfVars>::Apply(e, Jacobian_);
	}


public:
	template<typename T>
	inline void PrintValue(Daixt::Expr<T> const & e,
						   std::size_t index)
	{
		std::cout << "Value of "
				  << (Daixt::ChangeDisambiguation<Accessor>(e)
					  .GetName())
				  << "[" << index << "] = "
				  << (Daixt::ChangeDisambiguation<Accessor>(e)
					  .GetValue(Values_, index))
				  << std::endl;
	}

	inline void PrintJacobian(std::size_t index)
	{
		std::cout << "\nValues of Jacobian at " << index << ":\n";

		for (int j = 1; j < NumberOfVars + 1; ++j)
		{ 
			for (int i = 1; i < NumberOfVars + 1; ++i)
			{ 
				std::cout << Jacobian_(i, j)->GetValue(Values_, index) << " ";
			}
			std::cout << std::endl;
		}
	}

        inline void Step(std::size_t index)
        {
    	        std::size_t i = index + 1;
	        std::cout << "\nStep of solution at " << index << ":\n";
		// This is where to find the values of the variables
                double x = Values_[index * NumberOfVars];
                double y = Values_[index * NumberOfVars + 1];
		std::cout << "x  = " << x << std::endl;
		std::cout << "y  = " << y << std::endl;
		// This is where to find the values of the functions
                double f1 = Expressions_(1)->GetValue(Values_,index);
                double f2 = Expressions_(2)->GetValue(Values_,index);
		std::cout << "f1 = " << f1 << std::endl;
		std::cout << "f2 = " << f2 << std::endl;
		// This is where to find the values of the jacobian
                double g11 = Jacobian_(1, 1)->GetValue(Values_, index);
                double g12 = Jacobian_(1, 2)->GetValue(Values_, index);
                double g21 = Jacobian_(2, 1)->GetValue(Values_, index);
                double g22 = Jacobian_(2, 2)->GetValue(Values_, index);
		// compute determinant
                double det = g22*g11 - g12*g21;
		// compute corrections
                double delx = (f2*g12 - f1*g22)/det;
                double dely = (f2*g11 - f1*g21)/det;
                // adjust variables for the next step.
                Values_[i * NumberOfVars] = x + delx;
                Values_[i * NumberOfVars + 1] = y - dely;
		// need to check for convergence.
        }
};


////////////////////////////////////////////////////////////////////////////////
// ************************************************************************** //
////////////////////////////////////////////////////////////////////////////////


DECLARE_VARIABLE(x, 1, "x");
DECLARE_VARIABLE(y, 2, "y");
DECLARE_VARIABLE(z, 3, "z");



int main()
{
	using namespace Daixt::DefaultOps;
	using namespace TinyMatAndVec; 

	Solver<3, 
		TinyVector<AccessibleBaseClass<3> *, 3>, 
		TinyQuadraticMatrix<AccessibleBaseClass<3> *, 3> 
		> MySolver;

        MySolver.AddExpression<1>(x + y);
        MySolver.AddExpression<2>(x * x + y * y - z);
	MySolver.AddExpression<3>(z - z);

	int n = 10;
	MySolver.Resize(n+1);

	MySolver.PrintValue(Daixt::make_expr(x), 2);
	MySolver.PrintValue(Daixt::make_expr(y), 2);
	MySolver.PrintValue(Daixt::make_expr(z), 2);

	//MySolver.PrintJacobian(2);

	for (int i = 0; i < n; i++)
           MySolver.Step(i);


}