Problem with Solver

Marcus

Here is a copy of a program which fails when compiled with 

-DPROBLEM

to illustrate the problem I amhaving when I pass expressions to a 
function and then attempt to add them to the solver.

I have not tried this with the new release of Daixtrose.

John-- 
Dr John P. Fletcher   Tel: (44) 121 204 3389 (direct line)
Chemical Engineering and Applied Chemistry (CEAC),
School of Engineering and Applied Science (SEAS),
Aston University, Aston Triangle, BIRMINGHAM B4 7ET  U.K.               
CEAC Web site http://www.ceac.aston.ac.uk/
FAX: (44) 121 204 3679

// mysolver_6.cc adapted and extended 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.

// Modified to illustrate the problem of Expression<n>.

// Compile with -DPROBLEM to see the problem.

// John Fletcher ([email protected]) 2006

#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;

        Accessor *Determinant_;

	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<typename T1, typename T2>
        inline void StoreDeterminant(Daixt::Expr<T1> const & e1,Daixt::Expr<T2> const & e2)
        {
	      using namespace Daixt::ExprManip;
	      using namespace Daixt::Differentiation;
	      Variable<1> Diff1;
	      Variable<2> Diff2;
	      Determinant_ = Daixt::ChangeDisambiguation<Accessor>
                  (Simplify
		   (Diff(e1, Diff1)*Diff(e2, Diff2) - Diff(e2, Diff1)*Diff(e1, Diff2) )
                  ).DeepCopy(); 
              std::cerr << "Determinant = " << Determinant_->GetName() <<  std::endl;
        }



	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 << "registering expression '" 
				  << Expressions_(Row)->GetName() << "'" << std::endl;

		// this stores all d(e)/d(var(i))
		// This code has to be commented to save the call being compiled.
		//if (diffs) StoreDiff<Row, NumberOfVars>::Apply(e, Jacobian_);
	}


        template<typename T>
	void AddExpression1(Daixt::Expr<T> const & e)
        {
	  AddExpression<1>(e);
        }  
        template<typename T>
	void AddExpression2(Daixt::Expr<T> const & e)
        {
	  AddExpression<2>(e);
        }  
        template<typename T>
	void AddExpression3(Daixt::Expr<T> const & e)
        {
	  AddExpression<3>(e);
        }  
        template<typename T>
	void AddExpression4(Daixt::Expr<T> const & e)
        {
	  AddExpression<4>(e);
        }  

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 double 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 = g11*g22 - g12*g21;
                double det = Determinant_->GetValue(Values_, index);
		std::cout << "det  = " << det  << std::endl;
		//std::cout << "det2 = " << det2 << std::endl;
		// compute corrections
                double delx = (Expressions_(3)->GetValue(Values_,index))/det;
                //double delx = (f2*g12 - f1*g22)/det;
		std::cout << "delx  = " << delx  << std::endl;
		//std::cout << "delx2 = " << delx2  << std::endl;
                //double dely = (-f2*g11 + f1*g21)/det;
                double dely = (Expressions_(4)->GetValue(Values_,index))/det;
               // adjust variables for the next step.
 		std::cout << "dely  = " << dely  << std::endl;
		//std::cout << "dely2 = " << dely2  << std::endl;
                Values_[i * NumberOfVars] = x + delx;
                Values_[i * NumberOfVars + 1] = y + dely;
		// need to check for convergence.
                return f1*f1 + f2*f2;
        }
};


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


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

// Versions like this with the expressions e1 and e2 as parameters
// failed to compile. 
// The failures are at the calls to AddExpression<n> with the message
// error: invalid use of member (did you forget the '&' ?)
// The call to StoreDeterminant is O.K.
// In version 4 the same code is compiled O.K. and runs where the expression is not an argument.
// This looks like a problem with the compiler, something to do with the integer template value.
// The solution is to add new member functions AddExpressionN. The calls compile and so does the body
// containing the call which did not compile here. Reason???

// Advantage: Consistent code where each expression is entered once at the user's level.

template <class T1, class T2>
int problem(Daixt::Expr<T1> const & e1,Daixt::Expr<T2> const & e2)
{
	using namespace Daixt::DefaultOps;
	using namespace TinyMatAndVec; 
	using namespace Daixt::ExprManip;
	using namespace Daixt::Differentiation;

	typedef AccessibleBaseClass<2> Accessor_t;
        typedef TinyVector<Accessor_t *, 4> Vector_t;
	//typedef TinyQuadraticMatrix<Accessor_t *, 4> Matrix_t;
	typedef Solver<2, Vector_t> Solver_t;

	Solver_t MySolver;

#ifdef PROBLEM
        MySolver.AddExpression<1>(e1);
#else
        MySolver.AddExpression1(e1);
#endif
        MySolver.AddExpression2(e2);
	// This is how to make a variable into an expression.
	//MySolver.AddExpression<3>(Daixt::make_expr(z));
        MySolver.AddExpression3(
	    Simplify(e2*Diff(e1,y)-(e1*Diff(e2,y)))
        );
        MySolver.AddExpression4(
	    Simplify(e1*Diff(e2,x)-(e2*Diff(e1,x)))
        );
	// Compile time evaluation of 2 by 2 determinant from expressions.
	// I have not been able to do this from the stored versions.
	MySolver.StoreDeterminant(e1, e2);

	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++) {
            double residual = MySolver.Step(i);
	    std::cout << "Residual = " << residual << std::endl;
            if (residual < 1.e-16) {
	        std::cout << "Converged" << std::endl;
                break; 
	    }
        }
}


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

        problem(x+y,x*x + y*y -3.);

}