Commit 11a11fa5 authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files

Some fixes to branching expressions in automatic differentiation

parent 5aaf8930
...@@ -212,9 +212,12 @@ int main(int argc,char ** argv) ...@@ -212,9 +212,12 @@ int main(int argc,char ** argv)
NK(k,1) = nK(0,1); NK(k,1) = nK(0,1);
NK(k,2) = nK(0,2); NK(k,2) = nK(0,2);
} //end of loop over faces } //end of loop over faces
rMatrix SU,SS,SV;
nKGRAD = NK*(NK.Transpose()*R).Invert(true).first*NK.Transpose(); //stability part nKGRAD = NK*(NK.Transpose()*R).Invert(true).first*NK.Transpose(); //stability part
nKGRAD += (rMatrix::Unit(NF) - R*(R.Transpose()*R).Invert(true).first*R.Transpose())*(2.0/(static_cast<real>(NF)*vP)*(NK*K.Invert(true).first*NK.Transpose()).Trace());
/* /*
rMatrix SU,SS,SV;
std::cout << "W" << std::endl; std::cout << "W" << std::endl;
nKGRAD.Print(); nKGRAD.Print();
nKGRAD.SVD(SU,SS,SV); nKGRAD.SVD(SU,SS,SV);
...@@ -226,7 +229,7 @@ int main(int argc,char ** argv) ...@@ -226,7 +229,7 @@ int main(int argc,char ** argv)
SV.Print(); SV.Print();
std::cout << "Check " << (nKGRAD - SU*SS*SV.Transpose()).FrobeniusNorm() << std::endl; std::cout << "Check " << (nKGRAD - SU*SS*SV.Transpose()).FrobeniusNorm() << std::endl;
*/ */
/*
int rank = 0; //size of matrix U int rank = 0; //size of matrix U
{ //Retrive orthogonal to R matrix D { //Retrive orthogonal to R matrix D
...@@ -353,7 +356,7 @@ int main(int argc,char ** argv) ...@@ -353,7 +356,7 @@ int main(int argc,char ** argv)
//(U*D.Transpose()*R).Print(); //(U*D.Transpose()*R).Print();
nKGRAD += D*U*D.Transpose(); nKGRAD += D*U*D.Transpose();
*/
//std::cout << "W: " << std::endl; //std::cout << "W: " << std::endl;
//nKGRAD.Print(); //nKGRAD.Print();
bulk & isDMP = cell->Bulk(tag_DMP); bulk & isDMP = cell->Bulk(tag_DMP);
......
...@@ -1296,7 +1296,7 @@ namespace INMOST ...@@ -1296,7 +1296,7 @@ namespace INMOST
throw NotImplemented; throw NotImplemented;
} }
}; };
template<class A, class B, class C> template<class A, class B, class C>
class condition_expression : public shell_expression<condition_expression<A,B,C> > class condition_expression : public shell_expression<condition_expression<A,B,C> >
{ {
......
...@@ -454,76 +454,92 @@ namespace INMOST ...@@ -454,76 +454,92 @@ namespace INMOST
abstract_dynamic_variable * Copy() const {return static_cast<abstract_dynamic_variable *>(new table_variable(*this));} abstract_dynamic_variable * Copy() const {return static_cast<abstract_dynamic_variable *>(new table_variable(*this));}
}; };
/// This class makes possible to evaluate different expressions on different element types.
template<class A, class B> /// See etype_branch function.
class etype_branch_variable : public shell_dynamic_variable< binary_pool_expression< branch_expression<typename A::Var, typename B::Var>, typename A::Var, typename B::Var >, etype_branch_variable<A,B> > template<class A, class B>
{ class etype_branch_variable : public shell_dynamic_variable< multivar_expression, etype_branch_variable<A,B> >
private: {
A ArgA; private:
B ArgB; A ArgA; //< Variable expression to be evaluated when type of provided element matches selected types.
ElementType types_true; B ArgB; //< Variable expression to be evaluated when type of provided element does not match selected types.
public: ElementType types_true; //< Selected types of elements.
etype_branch_variable(ElementType _types_true, const A & _ArgA, const B & _ArgB) : types_true(_types_true), ArgA(_ArgA), ArgB(_ArgB) {} public:
etype_branch_variable(const etype_branch_variable & other) : types_true(other.types_true), ArgA(other.ArgA), ArgB(other.ArgB) {} /// Constructor. Used by etype_branch function.
etype_branch_variable & operator =(etype_branch_variable const & other) etype_branch_variable(ElementType _types_true, const A & _ArgA, const B & _ArgB) : types_true(_types_true), ArgA(_ArgA), ArgB(_ArgB) {}
{ /// Copy constructor.
types_true = other.types_true; etype_branch_variable(const etype_branch_variable & other) : types_true(other.types_true), ArgA(other.ArgA), ArgB(other.ArgB) {}
ArgA = other.ArgA; /// Assignment operator.
ArgB = other.ArgB; etype_branch_variable & operator =(etype_branch_variable const & other)
return *this; {
} types_true = other.types_true;
INMOST_DATA_REAL_TYPE Value(const Storage & e) const {return (*this)[e].GetValue();} ArgA = other.ArgA;
multivar_expression Variable(const Storage & e) const ArgB = other.ArgB;
{ return *this;
multivar_expression ret = (*this)[e]; }
return ret; /// Get value of variable expression on provided element e.
} INMOST_DATA_REAL_TYPE Value(const Storage & e) const {return (*this)[e].GetValue();}
binary_pool_expression<branch_expression<typename A::Var,typename B::Var>, typename A::Var, typename B::Var > operator [](const Storage & e) const /// Get value with derivatives of variable expression on provided element e.
{ /// This function collapses associated expression tree into multivar_expression.
binary_pool<branch_expression<typename A::Var,typename B::Var>,typename A::Var,typename B::Var> pool(ArgA[e],ArgB[e]); multivar_expression Variable(const Storage & e) const { return (*this)[e]; }
pool.get_op().SetCondition(e->GetElementType() & types_true ? true : false); /// Build an expression associated with variable expression on provided element e.
return binary_pool_expression<branch_expression<typename A::Var,typename B::Var>, typename A::Var, typename B::Var >(pool); multivar_expression operator [](const Storage & e) const
} {
void GetVariation(const Storage & e, Sparse::Row & r) const { (*this)[e].GetJacobian(1.0,r); } if( e->GetElementType() & types_true )
void GetVariation(const Storage & e, Sparse::RowMerger & r) const { (*this)[e].GetJacobian(1.0,r); } return ArgA[e];
abstract_dynamic_variable * Copy() const {return static_cast<abstract_dynamic_variable *>(new etype_branch_variable(*this));} else return ArgB[e];
}; }
/// Retrive first derivatives of variable expression on provided element e, default approach.
template<class A, class B> void GetVariation(const Storage & e, Sparse::Row & r) const { (*this)[e].GetJacobian(1.0,r); }
class marker_branch_variable : public shell_dynamic_variable< binary_pool_expression< branch_expression<typename A::Var, typename B::Var>, typename A::Var, typename B::Var >, marker_branch_variable<A,B> > /// Retrive first derivatives of variable expression on provided element e, with supplimentary structure Sparse::RowMerger.
{ void GetVariation(const Storage & e, Sparse::RowMerger & r) const { (*this)[e].GetJacobian(1.0,r); }
private: /// Make a copy of this class, used to reproduce and store a tree of variable expressions.
A ArgA; abstract_dynamic_variable * Copy() const {return static_cast<abstract_dynamic_variable *>(new etype_branch_variable(*this));}
B ArgB; };
MarkerType marker;
public: /// This class makes possible to evaluate different expressions depending on the markers.
marker_branch_variable(MarkerType _marker, const A & _ArgA, const B & _ArgB) : marker(_marker), ArgA(_ArgA), ArgB(_ArgB) {} /// Works similarly for shared and private markers.
marker_branch_variable(const marker_branch_variable & other) : marker(other.marker), ArgA(other.ArgA), ArgB(other.ArgB) {} /// See marker_branch function.
marker_branch_variable & operator =(marker_branch_variable const & other) template<class A, class B>
{ class marker_branch_variable : public shell_dynamic_variable< multivar_expression, marker_branch_variable<A,B> >
marker = other.marker; {
ArgA = other.ArgA; private:
ArgB = other.ArgB; A ArgA; //< Variable expression to be evaluated when marker is set on the element.
return *this; B ArgB; //< Variable expression to be evaluated when marker is not set on the element.
} MarkerType marker; //< Marker.
INMOST_DATA_REAL_TYPE Value(const Storage & e) const {return (*this)[e].GetValue();} public:
multivar_expression Variable(const Storage & e) const /// Constructor. Used by marker_branch function.
{ marker_branch_variable(MarkerType _marker, const A & _ArgA, const B & _ArgB) : marker(_marker), ArgA(_ArgA), ArgB(_ArgB) {}
multivar_expression ret = (*this)[e]; /// Copy constructor.
return ret; marker_branch_variable(const marker_branch_variable & other) : marker(other.marker), ArgA(other.ArgA), ArgB(other.ArgB) {}
} /// Assignment operator.
binary_pool_expression<branch_expression<typename A::Var,typename B::Var>, typename A::Var, typename B::Var > operator [](const Storage & e) const marker_branch_variable & operator =(marker_branch_variable const & other)
{ {
binary_pool<branch_expression<typename A::Var,typename B::Var>,typename A::Var,typename B::Var> pool(ArgA[e],ArgB[e]); marker = other.marker;
pool.get_op().SetCondition(isPrivateMarker(marker) ? e->GetPrivateMarker(marker) : e->GetMarker(marker)); ArgA = other.ArgA;
return binary_pool_expression<branch_expression<typename A::Var,typename B::Var>, typename A::Var, typename B::Var >(pool); ArgB = other.ArgB;
} return *this;
void GetVariation(const Storage & e, Sparse::Row & r) const { (*this)[e].GetJacobian(1.0,r); } }
void GetVariation(const Storage & e, Sparse::RowMerger & r) const { (*this)[e].GetJacobian(1.0,r); } /// Get value of variable expression on provided element e.
abstract_dynamic_variable * Copy() const {return static_cast<abstract_dynamic_variable *>(new marker_branch_variable(*this));} INMOST_DATA_REAL_TYPE Value(const Storage & e) const {return (*this)[e].GetValue();}
}; /// Get value with derivatives of variable expression on provided element e.
/// This function collapses associated expression tree into multivar_expression.
multivar_expression Variable(const Storage & e) const { return (*this)[e]; }
/// Build an expression associated with variable expression on provided element e.
multivar_expression operator [](const Storage & e) const
{
if( isPrivateMarker(marker) ? e->GetPrivateMarker(marker) : e->GetMarker(marker) )
return ArgA[e];
else return ArgB[e];
}
/// Retrive first derivatives of variable expression on provided element e, default approach.
void GetVariation(const Storage & e, Sparse::Row & r) const { (*this)[e].GetJacobian(1.0,r); }
/// Retrive first derivatives of variable expression on provided element e, with supplimentary structure Sparse::RowMerger.
void GetVariation(const Storage & e, Sparse::RowMerger & r) const { (*this)[e].GetJacobian(1.0,r); }
/// Make a copy of this class, used to reproduce and store a tree of variable expressions.
abstract_dynamic_variable * Copy() const {return static_cast<abstract_dynamic_variable *>(new marker_branch_variable(*this));}
};
template<class Expr, class A> template<class Expr, class A>
class unary_custom_variable : public shell_dynamic_variable< unary_pool_expression<Expr, typename A::Var >,unary_custom_variable<Expr,A> > class unary_custom_variable : public shell_dynamic_variable< unary_pool_expression<Expr, typename A::Var >,unary_custom_variable<Expr,A> >
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment