diff --git a/CHANGELOG.md b/CHANGELOG.md index a189e50..4f67a85 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,13 @@ All notable changes to `libcasm-clexulator` will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [2.0a4] - 2023-10-26 + +### Added + +- Added gradients of correlations functionality in C++, Python functions +- Added FADBAD library to support compiling of Clexulators capable of calculating gradients + ## [2.0a3] - 2023-10-25 ### Fixed diff --git a/CMakeLists.txt b/CMakeLists.txt index d3aa357..ab53021 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -78,48 +78,52 @@ ENDIF(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT) # create libcasm_clexulator set( libcasm_clexulator_HEADERS - ${PROJECT_SOURCE_DIR}/include/casm/clexulator/ConfigDoFValues.hh + ${PROJECT_SOURCE_DIR}/include/casm/clexulator/Correlations.hh + ${PROJECT_SOURCE_DIR}/include/casm/clexulator/BasicClexParamPack.hh + ${PROJECT_SOURCE_DIR}/include/casm/clexulator/SparseCoefficients.hh ${PROJECT_SOURCE_DIR}/include/casm/clexulator/NeighborList.hh - ${PROJECT_SOURCE_DIR}/include/casm/clexulator/LocalClusterExpansion.hh - ${PROJECT_SOURCE_DIR}/include/casm/clexulator/ConfigDoFValuesTools.hh + ${PROJECT_SOURCE_DIR}/include/casm/clexulator/ImpactTable.hh + ${PROJECT_SOURCE_DIR}/include/casm/clexulator/ConfigDoFValuesTools_impl.hh ${PROJECT_SOURCE_DIR}/include/casm/clexulator/ClusterExpansion.hh - ${PROJECT_SOURCE_DIR}/include/casm/clexulator/OrderParameter.hh - ${PROJECT_SOURCE_DIR}/include/casm/clexulator/DoFSpace.hh - ${PROJECT_SOURCE_DIR}/include/casm/clexulator/DoFSpaceAxisInfo.hh - ${PROJECT_SOURCE_DIR}/include/casm/clexulator/Clexulator.hh - ${PROJECT_SOURCE_DIR}/include/casm/clexulator/SparseCoefficients.hh + ${PROJECT_SOURCE_DIR}/include/casm/clexulator/ConfigDoFValues.hh ${PROJECT_SOURCE_DIR}/include/casm/clexulator/LocalCorrelations.hh - ${PROJECT_SOURCE_DIR}/include/casm/clexulator/ConfigDoFValuesTools_impl.hh - ${PROJECT_SOURCE_DIR}/include/casm/clexulator/ImpactTable.hh - ${PROJECT_SOURCE_DIR}/include/casm/clexulator/Correlations.hh + ${PROJECT_SOURCE_DIR}/include/casm/clexulator/BaseClexulator.hh + ${PROJECT_SOURCE_DIR}/include/casm/clexulator/DoFSpace.hh + ${PROJECT_SOURCE_DIR}/include/casm/clexulator/LocalClusterExpansion.hh ${PROJECT_SOURCE_DIR}/include/casm/clexulator/version.hh + ${PROJECT_SOURCE_DIR}/include/casm/clexulator/Clexulator.hh + ${PROJECT_SOURCE_DIR}/include/casm/clexulator/ConfigDoFValuesTools.hh + ${PROJECT_SOURCE_DIR}/include/casm/clexulator/DoFSpaceAxisInfo.hh ${PROJECT_SOURCE_DIR}/include/casm/clexulator/ClexParamPack.hh - ${PROJECT_SOURCE_DIR}/include/casm/clexulator/BasicClexParamPack.hh - ${PROJECT_SOURCE_DIR}/include/casm/clexulator/BaseClexulator.hh ${PROJECT_SOURCE_DIR}/include/casm/clexulator/DiffClexParamPack.hh + ${PROJECT_SOURCE_DIR}/include/casm/clexulator/OrderParameter.hh + ${PROJECT_SOURCE_DIR}/include/casm/clexulator/external/fadbad/fadbad.h + ${PROJECT_SOURCE_DIR}/include/casm/clexulator/external/fadbad/badiff.h + ${PROJECT_SOURCE_DIR}/include/casm/clexulator/external/fadbad/tadiff.h + ${PROJECT_SOURCE_DIR}/include/casm/clexulator/external/fadbad/fadiff.h + ${PROJECT_SOURCE_DIR}/include/casm/clexulator/io/json/ConfigDoFValues_json_io.hh ${PROJECT_SOURCE_DIR}/include/casm/clexulator/io/json/SparseCoefficients_json_io.hh ${PROJECT_SOURCE_DIR}/include/casm/clexulator/io/json/Clexulator_json_io.hh ${PROJECT_SOURCE_DIR}/include/casm/clexulator/io/json/DoFSpace_json_io.hh - ${PROJECT_SOURCE_DIR}/include/casm/clexulator/io/json/ConfigDoFValues_json_io.hh ) set( libcasm_clexulator_SOURCES - ${PROJECT_SOURCE_DIR}/src/casm/clexulator/LocalCorrelations.cc - ${PROJECT_SOURCE_DIR}/src/casm/clexulator/ImpactTable.cc - ${PROJECT_SOURCE_DIR}/src/casm/clexulator/Correlations.cc - ${PROJECT_SOURCE_DIR}/src/casm/clexulator/ClusterExpansion.cc + ${PROJECT_SOURCE_DIR}/src/casm/clexulator/LocalClusterExpansion.cc ${PROJECT_SOURCE_DIR}/src/casm/clexulator/OrderParameter.cc - ${PROJECT_SOURCE_DIR}/src/casm/clexulator/DoFSpace.cc ${PROJECT_SOURCE_DIR}/src/casm/clexulator/DoFSpaceAxisInfo.cc + ${PROJECT_SOURCE_DIR}/src/casm/clexulator/ClusterExpansion.cc ${PROJECT_SOURCE_DIR}/src/casm/clexulator/Clexulator.cc - ${PROJECT_SOURCE_DIR}/src/casm/clexulator/LocalClusterExpansion.cc - ${PROJECT_SOURCE_DIR}/src/casm/clexulator/NeighborList.cc ${PROJECT_SOURCE_DIR}/src/casm/clexulator/BaseClexulator.cc + ${PROJECT_SOURCE_DIR}/src/casm/clexulator/NeighborList.cc + ${PROJECT_SOURCE_DIR}/src/casm/clexulator/ImpactTable.cc + ${PROJECT_SOURCE_DIR}/src/casm/clexulator/DoFSpace.cc + ${PROJECT_SOURCE_DIR}/src/casm/clexulator/LocalCorrelations.cc + ${PROJECT_SOURCE_DIR}/src/casm/clexulator/Correlations.cc ${PROJECT_SOURCE_DIR}/src/casm/clexulator/version.cc - ${PROJECT_SOURCE_DIR}/src/casm/clexulator/io/json/SparseCoefficients_json_io.cc - ${PROJECT_SOURCE_DIR}/src/casm/clexulator/io/json/ConfigDoFValues_json_io.cc ${PROJECT_SOURCE_DIR}/src/casm/clexulator/io/json/DoFSpace_json_io.cc ${PROJECT_SOURCE_DIR}/src/casm/clexulator/io/json/Clexulator_json_io.cc + ${PROJECT_SOURCE_DIR}/src/casm/clexulator/io/json/ConfigDoFValues_json_io.cc + ${PROJECT_SOURCE_DIR}/src/casm/clexulator/io/json/SparseCoefficients_json_io.cc ) add_library(casm_clexulator SHARED ${libcasm_clexulator_SOURCES}) target_include_directories(casm_clexulator diff --git a/include/casm/clexulator/Correlations.hh b/include/casm/clexulator/Correlations.hh index 1c55947..4acca34 100644 --- a/include/casm/clexulator/Correlations.hh +++ b/include/casm/clexulator/Correlations.hh @@ -128,6 +128,18 @@ class Correlations { DoFKey const &key, Eigen::VectorXd const &new_value, Eigen::VectorXd const &reference_extensive_correlations); + // --- Gradients of correlations + /// \brief Calculates and returns gradients of correlations with respect + /// to DoF corresponding to key + /// Returns a M x N matrix, where M represents the number of correlations, + /// and N represents the dimensions of the degrees of freedom corresponding to + /// key. For example, if key is Hstrain, N will be 6. If key is disp, and + /// there are 3 sites in the configuration, N will be 9 (x, y, z displacements + /// for each site). Each row of the matrix corresponds to gradients of all M + /// correlations with respect to one of the dimensions of degrees of freedom + /// corresponding to key + Eigen::MatrixXd const grad_correlations(DoFKey const &key); + private: /// Holds all correlation indices std::vector m_correlation_indices; diff --git a/include/casm/clexulator/DiffClexParamPack.hh b/include/casm/clexulator/DiffClexParamPack.hh index dedeb51..46f9fab 100644 --- a/include/casm/clexulator/DiffClexParamPack.hh +++ b/include/casm/clexulator/DiffClexParamPack.hh @@ -9,8 +9,8 @@ #include "casm/casm_io/enum/json_io.hh" #include "casm/casm_io/enum/stream_io.hh" #include "casm/clexulator/ClexParamPack.hh" -#include "casm/external/fadbad/badiff.h" -#include "casm/external/fadbad/fadiff.h" +#include "casm/clexulator/external/fadbad/badiff.h" +#include "casm/clexulator/external/fadbad/fadiff.h" #include "casm/global/definitions.hh" namespace CASM { diff --git a/include/casm/clexulator/external/fadbad/COPYRIGHT b/include/casm/clexulator/external/fadbad/COPYRIGHT new file mode 100755 index 0000000..9831b36 --- /dev/null +++ b/include/casm/clexulator/external/fadbad/COPYRIGHT @@ -0,0 +1,28 @@ +Copyright (C) 1996-2006 Ole Stauning (fadbad@uning.dk) +All rights reserved. + +This code is provided "as is", without any warranty of any kind, +either expressed or implied, including but not limited to, any implied +warranty of merchantibility or fitness for any purpose. In no event +will any party who distributed the code be liable for damages or for +any claim(s) by any other party, including but not limited to, any +lost profits, lost monies, lost data or data rendered inaccurate, +losses sustained by third parties, or any other special, incidental or +consequential damages arising out of the use or inability to use the +program, even if the possibility of such damages has been advised +against. The entire risk as to the quality, the performance, and the +fitness of the program for any particular purpose lies with the party +using the code. + +This code, and any derivative of this code, may not be used in a +commercial package without the prior explicit written permission of +the authors. Verbatim copies of this code may be made and distributed +in any medium, provided that this copyright notice is not removed or +altered in any way. No fees may be charged for distribution of the +codes, other than a fee to cover the cost of the media and a +reasonable handling fee. + +*************************************************************** +ANY USE OF THIS CODE CONSTITUTES ACCEPTANCE OF THE TERMS OF THE + COPYRIGHT NOTICE +*************************************************************** diff --git a/include/casm/clexulator/external/fadbad/README b/include/casm/clexulator/external/fadbad/README new file mode 100755 index 0000000..65a7ce3 --- /dev/null +++ b/include/casm/clexulator/external/fadbad/README @@ -0,0 +1,76 @@ + + FADBAD++ Templates for Automatic Differentiation + ================================================ + + FADBAD++ defines C++ templates for performing automatic differentiation +of functions implemented as C/C++ programs. Three types of automatic +differentiation has been implemented. The forward and the backward methods +are implemented to compute first order derivatives and the Taylor expansion +method is implemented to compute Taylor series expansion. + All templates are flexible in the way that the arithmetic used in as +template-arguments can be chosen by the user. This way differentiation of +programs, based on any arithmetic such as: double, interval, multiprecision +etc. is possible. This flexibility also makes it possible to perform automatic +differentiation on a program which itself uses automatic differentiation. +All three methods can be mixed in an application and the user can obtain +derivatives using the most optimal combination of methods for the application. + +APPLICATIONS: + + * Forward automatic differentiation on a function f : R^n->R^n, evaluated + in interval arithmetics, using the BIAS/PROFIL package, to obtain function + values and derivatives. Used with the interval Newton method to obtain + guarenteed enlosures of all solutions to the nonlinear equation f(x)=0. + + * Forward-Backward automatic differentiation on a function f : R^n->R, + evaluated in interval arithmetics, using the BIAS/PROFIL package and the + backward method to obtain first order derivatives (the gradient) and the + forward method to differentiate the first order derivatives to obtain the + second order derivatives (the Hessian). Used with the interval Krawczyk + method to perform global optimization obtaining a guarenteed enclosure of + the global minimum. + + * Numerical integration of a function f : RxR^n->R, using a + three-point-two-derivative formula. The Backward method has been used to + differentiate the Numerical integration program obtaining the n partial + derivatives of the integral with respect to the n parameters in the + function. + + * Taylor expansion of the solution to an ordinary differential equation, + used to solve initial value problems. The Forward method has been used to + differentiate the initial value problem solver to obtain the solution of + the variational problem. + + * Taylor expansion of the solution to an ordinary differential equation + using interval arithmetics and using the Forward method to obtain + derivatives of the Taylor coefficients with respect to the point of + expansion, which are the values of the Taylor coefficients for the + solution of the variational problem. Used in a method which solves initial + value problems with guaranteed enclosures. + +LICENSING: + +FADBAD++ is distributed under the dual licensing business model similer to the +model used by MySQL, Trolltech and Sleepycat. In return for the advantages you +realize from using FADBAD++ in your application, we require that you do one of +the following: + + * Either: Contribute to the continued development of the product by + purchasing commercial licenses from the authors. This option secures you + the right to distribute your application under the license terms of your + choice. + + * Or: Contribute to the Open Source community by placing your application + under an Open Source license (e.g. the GPL). This option secures all users + the rights to obtain the application's full source code, modify it, and + redistribute it. + +Contact the authors if you want to obtain a commercial licence or if you are +unsure about what license to use . + +More information about automatic differentiation in general and source code +for FADBAD++ with documentation can be obtained from the FADBAD++ homepage: + + http://www.fadbad.com + + diff --git a/include/casm/clexulator/external/fadbad/badiff.h b/include/casm/clexulator/external/fadbad/badiff.h new file mode 100644 index 0000000..e1fbbcd --- /dev/null +++ b/include/casm/clexulator/external/fadbad/badiff.h @@ -0,0 +1,1125 @@ +// Copyright (C) 1996-2007 Ole Stauning & Claus Bendtsen (fadbad@uning.dk) +// All rights reserved. + +// This code is provided "as is", without any warranty of any kind, +// either expressed or implied, including but not limited to, any implied +// warranty of merchantibility or fitness for any purpose. In no event +// will any party who distributed the code be liable for damages or for +// any claim(s) by any other party, including but not limited to, any +// lost profits, lost monies, lost data or data rendered inaccurate, +// losses sustained by third parties, or any other special, incidental or +// consequential damages arising out of the use or inability to use the +// program, even if the possibility of such damages has been advised +// against. The entire risk as to the quality, the performance, and the +// fitness of the program for any particular purpose lies with the party +// using the code. + +// This code, and any derivative of this code, may not be used in a +// commercial package without the prior explicit written permission of +// the authors. Verbatim copies of this code may be made and distributed +// in any medium, provided that this copyright notice is not removed or +// altered in any way. No fees may be charged for distribution of the +// codes, other than a fee to cover the cost of the media and a +// reasonable handling fee. + +// *************************************************************** +// ANY USE OF THIS CODE CONSTITUTES ACCEPTANCE OF THE TERMS OF THE +// COPYRIGHT NOTICE +// *************************************************************** + +#ifndef _BADIFF_H +#define _BADIFF_H + +#include "fadbad.h" + +#include +#include + +namespace fadbad +{ + +template +class Derivatives +{ +public: + class RecycleBin + { + std::stack< std::vector* > m_recycle; + RecycleBin(const RecycleBin&){/*illegal*/} + public: + RecycleBin(){} + std::vector* popRecycle(const unsigned int n) + { + if (m_recycle.empty()) return new std::vector(n); + std::vector* elm=m_recycle.top(); + m_recycle.pop(); + USER_ASSERT(elm->size()==n,"Size mismatch "<size()<<"!="<* elm) + { + m_recycle.push(elm); + } + ~RecycleBin() + { + while(!m_recycle.empty()) + { + delete m_recycle.top(); + m_recycle.pop(); + } + } + }; +private: + std::vector* m_values; + + unsigned int size() const { return (unsigned int) m_values->size(); } + +public: + Derivatives():m_values(0){} + + void recycle(RecycleBin& bin) + { + USER_ASSERT(m_values!=0,"Nothing to recycle") + bin.pushRecycle(m_values); + m_values=0; + } + + bool haveValues() const { return m_values!=0; } + + U& diff(RecycleBin& bin, const unsigned int i, const unsigned int n) + { + USER_ASSERT(i::myZero(); + } + USER_ASSERT(m_values->size()==n,"Size mismatch "<size()<<"!="<::myOne(); + } + void add(RecycleBin& bin, const Derivatives& d) + { + USER_ASSERT(d.size()>0,"Propagating node with no derivatives") + if (m_values==0) + { + m_values=bin.popRecycle(d.size()); + for(unsigned int i=0;isize();++i) (*m_values)[i]=(*d.m_values)[i]; + } + else + { + USER_ASSERT(m_values->size()==d.size(),"Size mismatch "<size()<<"!="<size();++i) Op::myCadd((*m_values)[i],(*d.m_values)[i]); + } + } + void sub(RecycleBin& bin, const Derivatives& d) + { + USER_ASSERT(d.size()>0,"Propagating node with no derivatives") + if (m_values==0) + { + m_values=bin.popRecycle(d.size()); + for(unsigned int i=0;isize();++i) (*m_values)[i]=Op::myNeg((*d.m_values)[i]); + } + else + { + USER_ASSERT(m_values->size()==d.size(),"Size mismatch "<size()<<"!="<size();++i) Op::myCsub((*m_values)[i],(*d.m_values)[i]); + } + } + template + void add(RecycleBin& bin, const V& a, const Derivatives& d) + { + USER_ASSERT(d.size()>0,"Propagating node with no derivatives") + if (m_values==0) + { + m_values=bin.popRecycle(d.size()); + for(unsigned int i=0;isize();++i) (*m_values)[i]=a*(*d.m_values)[i]; + } + else + { + USER_ASSERT(m_values->size()==d.size(),"Size mismatch "<size()<<"!="<size();++i) Op::myCadd((*m_values)[i],a*(*d.m_values)[i]); + } + } + template + void sub(RecycleBin& bin, const V& a, const Derivatives& d) + { + USER_ASSERT(d.size()>0,"Propagating node with no derivatives") + if (m_values==0) + { + m_values=bin.popRecycle(d.size()); + for(unsigned int i=0;isize();++i) (*m_values)[i]=Op::myNeg(a*(*d.m_values)[i]); + } + else + { + USER_ASSERT(m_values->size()==d.size(),"Size mismatch "<size()<<"!="<size();++i) Op::myCsub((*m_values)[i],a*(*d.m_values)[i]); + } + } + + U& operator[](const unsigned int i) + { + if (m_values!=0) + { + USER_ASSERT(isize(),"Index "<size()<<"]") + return (*m_values)[i]; + } + else + { + static U zero; + zero=Op::myZero(); + return zero; + } + } + const U& operator[](const unsigned int i) const + { + if (m_values!=0) + { + USER_ASSERT(isize(),"Index "<size()<<"]") + return (*m_values)[i]; + } + else + { + static U zero; + zero=Op::myZero(); + return zero; + } + } +}; + + +template +class BTypeNameHV // Heap Value +{ + U m_val; + mutable unsigned int m_rc; + +protected: + mutable Derivatives m_derivatives; + virtual ~BTypeNameHV(){} +public: + BTypeNameHV():m_rc(0),m_derivatives(){} + template explicit BTypeNameHV(const V& val):m_val(val),m_rc(0),m_derivatives(){} + const U& val() const { return m_val; } + U& val() { return m_val; } + void decRef(BTypeNameHV*& pBTypeNameHV) + { + INTERNAL_ASSERT(m_rc>0,"Resource counter negative"); + if (--m_rc==0) + { + if (m_derivatives.haveValues()) + { + typename Derivatives::RecycleBin bin; + propagate(bin); + m_derivatives.recycle(bin); + propagateChildren(bin); + } + delete this; + } + pBTypeNameHV=0; + } + void decRef(typename Derivatives::RecycleBin& bin, BTypeNameHV*& pBTypeNameHV) + { + INTERNAL_ASSERT(m_rc>0,"Resource counter negative"); + if (--m_rc==0) + { + if (m_derivatives.haveValues()) + { + propagate(bin); + m_derivatives.recycle(bin); + propagateChildren(bin); + } + delete this; + } + pBTypeNameHV=0; + } + void incRef() const {++m_rc;} + + U& diff(typename Derivatives::RecycleBin& bin, const unsigned int idx, const unsigned int size) + { + return m_derivatives.diff(bin,idx,size); + } + virtual void propagate(typename Derivatives::RecycleBin&) {} + virtual void propagateChildren(typename Derivatives::RecycleBin&) {} + void add(typename Derivatives::RecycleBin& bin, const Derivatives& d) { m_derivatives.add(bin,d); } + void sub(typename Derivatives::RecycleBin& bin, const Derivatives& d) { m_derivatives.sub(bin,d); } + void add(typename Derivatives::RecycleBin& bin, const U& a, const Derivatives& d) { m_derivatives.add(bin,a,d); } + void sub(typename Derivatives::RecycleBin& bin, const U& a, const Derivatives& d) { m_derivatives.sub(bin,a,d); } + U& deriv(const unsigned int i) + { + USER_ASSERT(m_rc==1,"Still non-propagated dependencies ("< +class BTypeName +{ + struct SV // Stack Value refers to reference-counted Heap Value: + { + mutable BTypeNameHV* m_pBTypeNameHV; + SV(BTypeNameHV* pBTypeNameHV):m_pBTypeNameHV(pBTypeNameHV){ m_pBTypeNameHV->incRef(); } + SV(const typename BTypeName::SV& sv):m_pBTypeNameHV(sv.m_pBTypeNameHV){ m_pBTypeNameHV->incRef(); } + ~SV(){ m_pBTypeNameHV->decRef(m_pBTypeNameHV); } + BTypeNameHV* getBTypeNameHV() const { return m_pBTypeNameHV; } + void setBTypeNameHV(BTypeNameHV* pBTypeNameHV) + { + if (m_pBTypeNameHV!=pBTypeNameHV) + { + m_pBTypeNameHV->decRef(m_pBTypeNameHV); + m_pBTypeNameHV=pBTypeNameHV; + m_pBTypeNameHV->incRef(); + } + } + const U& val() const { return m_pBTypeNameHV->val(); } + U& val() { return m_pBTypeNameHV->val(); } + const U& deriv(const unsigned int i) const { return m_pBTypeNameHV->deriv(i); } + U& deriv(const unsigned int i) { return m_pBTypeNameHV->deriv(i); } + + U& diff(const unsigned int idx, const unsigned int size) + { + typename Derivatives::RecycleBin bin; + U& res(m_pBTypeNameHV->diff(bin,idx,size)); + BTypeNameHV* pHV=new BTypeNameHV(this->val()); + m_pBTypeNameHV->decRef(bin,m_pBTypeNameHV); + m_pBTypeNameHV=pHV; + m_pBTypeNameHV->incRef(); + return res; + } + } m_sv; +public: + typedef U UnderlyingType; + BTypeName():m_sv(new BTypeNameHV()){} + BTypeName(BTypeNameHV* pBTypeNameHV):m_sv(pBTypeNameHV){} + explicit BTypeName(const typename BTypeName::SV& sv):m_sv(sv){} + template /*explicit*/ BTypeName(const V& val):m_sv(new BTypeNameHV(val)){} + BTypeName& operator=(const BTypeName& val) + { + if (this==&val) return *this; + m_sv.setBTypeNameHV(val.m_sv.getBTypeNameHV()); + return *this; + } + template BTypeName& operator=(const V& val) { m_sv.setBTypeNameHV(new BTypeNameHV(val)); return *this; } + BTypeNameHV* getBTypeNameHV() const { return m_sv.getBTypeNameHV(); } + void setBTypeNameHV(const BTypeNameHV* pBTypeNameHV) { m_sv.setBTypeNameHV(pBTypeNameHV); } + const U& val() const { return m_sv.val(); } + U& val() { return m_sv.val(); } + U& x() { return m_sv.val(); } + const U& deriv(const unsigned int i) const { return m_sv.deriv(i); } + U& d(const unsigned int i) { return m_sv.deriv(i); } + const U& d(const unsigned int i) const { return m_sv.deriv(i); } + U& diff(const unsigned int idx, const unsigned int size) { return m_sv.diff(idx,size); } + + BTypeName& operator+=(const BTypeName& val); + BTypeName& operator-=(const BTypeName& val); + BTypeName& operator*=(const BTypeName& val); + BTypeName& operator/=(const BTypeName& val); + template BTypeName& operator+=(const V& val); + template BTypeName& operator-=(const V& val); + template BTypeName& operator*=(const V& val); + template BTypeName& operator/=(const V& val); +}; + +template bool operator==(const BTypeName& val1, const BTypeName& val2) { return Op::myEq(val1.val(),val2.val()); } +template bool operator!=(const BTypeName& val1, const BTypeName& val2) { return Op::myNe(val1.val(),val2.val()); } +template bool operator<(const BTypeName& val1, const BTypeName& val2) { return Op::myLt(val1.val(),val2.val()); } +template bool operator<=(const BTypeName& val1, const BTypeName& val2) { return Op::myLe(val1.val(),val2.val()); } +template bool operator>(const BTypeName& val1, const BTypeName& val2) { return Op::myGt(val1.val(),val2.val()); } +template bool operator>=(const BTypeName& val1, const BTypeName& val2) { return Op::myGe(val1.val(),val2.val()); } +template bool operator==(const BTypeName& val1, const V& val2) { return Op::myEq(val1.val(),val2); } +template bool operator==(const V& val1, const BTypeName& val2) { return Op::myEq(val1,val2.val()); } +template bool operator!=(const BTypeName& val1, const V& val2) { return Op::myNe(val1.val(),val2); } +template bool operator!=(const V& val1, const BTypeName& val2) { return Op::myNe(val1,val2.val()); } +template bool operator<(const BTypeName& val1, const V& val2) { return Op::myLt(val1.val(),val2); } +template bool operator<(const V& val1, const BTypeName& val2) { return Op::myLt(val1,val2.val()); } +template bool operator<=(const BTypeName& val1, const V& val2) { return Op::myLe(val1.val(),val2); } +template bool operator<=(const V& val1, const BTypeName& val2) { return Op::myLe(val1,val2.val()); } +template bool operator>(const BTypeName& val1, const V& val2) { return Op::myGt(val1.val(),val2); } +template bool operator>(const V& val1, const BTypeName& val2) { return Op::myGt(val1,val2.val()); } +template bool operator>=(const BTypeName& val1, const V& val2) { return Op::myGe(val1.val(),val2); } +template bool operator>=(const V& val1, const BTypeName& val2) { return Op::myGe(val1,val2.val()); } + +// Binary operator base class: + +template +class BinBTypeNameHV : public BTypeNameHV +{ + BTypeNameHV* m_pOp1; + BTypeNameHV* m_pOp2; +public: + BinBTypeNameHV(const U& val, BTypeNameHV* pOp1, BTypeNameHV* pOp2):BTypeNameHV(val),m_pOp1(pOp1),m_pOp2(pOp2) + { + m_pOp1->incRef(); + m_pOp2->incRef(); + } + virtual void propagateChildren(typename Derivatives::RecycleBin& bin) + { + m_pOp1->decRef(bin,m_pOp1); + m_pOp2->decRef(bin,m_pOp2); + } + virtual ~BinBTypeNameHV() + { + if (m_pOp1) m_pOp1->decRef(m_pOp1); + if (m_pOp2) m_pOp2->decRef(m_pOp2); + } + BTypeNameHV* op1() { return m_pOp1; } + BTypeNameHV* op2() { return m_pOp2; } +}; + +// Unary operator base class: + +template +class UnBTypeNameHV : public BTypeNameHV +{ + BTypeNameHV* m_pOp; +public: + UnBTypeNameHV(const U& val, BTypeNameHV* pOp):BTypeNameHV(val),m_pOp(pOp) + { + m_pOp->incRef(); + } + virtual void propagateChildren(typename Derivatives::RecycleBin& bin) + { + m_pOp->decRef(bin,m_pOp); + } + virtual ~UnBTypeNameHV() + { + if (m_pOp) m_pOp->decRef(m_pOp); + } + BTypeNameHV* op() { return m_pOp; } +}; + +// ADDITION: + +template +struct BTypeNameADD : public BinBTypeNameHV +{ + BTypeNameADD(const U& val, BTypeNameHV* pOp1, BTypeNameHV* pOp2):BinBTypeNameHV(val,pOp1,pOp2){} + virtual void propagate(typename Derivatives::RecycleBin& bin) + { + this->op1()->add(bin,this->m_derivatives); + this->op2()->add(bin,this->m_derivatives); + } +private: + void operator=(const BTypeNameADD&){} // not allowed +}; +template +struct BTypeNameADD1 : public UnBTypeNameHV +{ + const V m_a; + BTypeNameADD1(const U& val, const V& a, BTypeNameHV* pOp2):UnBTypeNameHV(val,pOp2),m_a(a){} + virtual void propagate(typename Derivatives::RecycleBin& bin) + { + this->op()->add(bin,this->m_derivatives); + } +private: + void operator=(const BTypeNameADD1&){} // not allowed +}; +template +struct BTypeNameADD2 : public UnBTypeNameHV +{ + const V m_b; + BTypeNameADD2(const U& val, BTypeNameHV* pOp1, const V& b):UnBTypeNameHV(val,pOp1),m_b(b){} + virtual void propagate(typename Derivatives::RecycleBin& bin) + { + this->op()->add(bin,this->m_derivatives); + } +private: + void operator=(const BTypeNameADD2&){} // not allowed +}; +template +BTypeName operator+(const BTypeName& val1, const BTypeName& val2) +{ + return BTypeName(static_cast*>(new BTypeNameADD(val1.val()+val2.val(),val1.getBTypeNameHV(),val2.getBTypeNameHV()))); +} +/* +template +BTypeName operator+(const typename Op::Underlying& a, const BTypeName& val2) +{ + return BTypeName(static_cast*>( + new BTypeNameADD1::Underlying>(a+val2.val(), a, val2.getBTypeNameHV()) + )); +} +template +BTypeName operator+(const BTypeName& val1, const typename Op::Underlying& b) +{ + return BTypeName(static_cast*>( + new BTypeNameADD2::Underlying>(val1.val()+b, val1.getBTypeNameHV(), b) + )); +} +template +BTypeName operator+(const typename Op::Base& a, const BTypeName& val2) +{ + return BTypeName(static_cast*>( + new BTypeNameADD1::Base>(a+val2.val(), a, val2.getBTypeNameHV()) + )); +} +template +BTypeName operator+(const BTypeName& val1, const typename Op::Base& b) +{ + return BTypeName(static_cast*>( + new BTypeNameADD2::Base>(val1.val()+b, val1.getBTypeNameHV(), b) + )); +} +*/ +template +BTypeName operator+(const V& a, const BTypeName& val2) +{ + return BTypeName(static_cast*>( + new BTypeNameADD1(a+val2.val(), a, val2.getBTypeNameHV()) + )); +} +template +BTypeName operator+(const BTypeName& val1, const V& b) +{ + return BTypeName(static_cast*>( + new BTypeNameADD2(val1.val()+b, val1.getBTypeNameHV(), b) + )); +} + +// SUBTRACTION: + +template +struct BTypeNameSUB : public BinBTypeNameHV +{ + BTypeNameSUB(const U& val, BTypeNameHV* pOp1, BTypeNameHV* pOp2):BinBTypeNameHV(val,pOp1,pOp2){} + virtual void propagate(typename Derivatives::RecycleBin& bin) + { + this->op1()->add(bin,this->m_derivatives); + this->op2()->sub(bin,this->m_derivatives); + } +private: + void operator=(const BTypeNameSUB&){} // not allowed +}; +template +struct BTypeNameSUB1 : public UnBTypeNameHV +{ + const V m_a; + BTypeNameSUB1(const U& val, const V& a, BTypeNameHV* pOp2):UnBTypeNameHV(val,pOp2),m_a(a){} + virtual void propagate(typename Derivatives::RecycleBin& bin) + { + this->op()->sub(bin,this->m_derivatives); + } +private: + void operator=(const BTypeNameSUB1&){} // not allowed +}; +template +struct BTypeNameSUB2 : public UnBTypeNameHV +{ + const V m_b; + BTypeNameSUB2(const U& val, BTypeNameHV* pOp1, const V& b):UnBTypeNameHV(val,pOp1),m_b(b){} + virtual void propagate(typename Derivatives::RecycleBin& bin) + { + this->op()->add(bin,this->m_derivatives); + } +private: + void operator=(const BTypeNameSUB2&){} // not allowed +}; +template +BTypeName operator-(const BTypeName& val1, const BTypeName& val2) +{ + return BTypeName(static_cast*>(new BTypeNameSUB(val1.val()-val2.val(),val1.getBTypeNameHV(),val2.getBTypeNameHV()))); +} +/* +template +BTypeName operator-(const typename Op::Underlying& a, const BTypeName& val2) +{ + return BTypeName(static_cast*>( + new BTypeNameSUB1::Underlying>(a-val2.val(), a, val2.getBTypeNameHV()) + )); +} +template +BTypeName operator-(const BTypeName& val1, const typename Op::Underlying& b) +{ + return BTypeName(static_cast*>( + new BTypeNameSUB2::Underlying>(val1.val()-b, val1.getBTypeNameHV(), b) + )); +} +template +BTypeName operator-(const typename Op::Base& a, const BTypeName& val2) +{ + return BTypeName(static_cast*>( + new BTypeNameSUB1::Base>(a-val2.val(), a, val2.getBTypeNameHV()) + )); +} +template +BTypeName operator-(const BTypeName& val1, const typename Op::Base& b) +{ + return BTypeName(static_cast*>( + new BTypeNameSUB2::Base>(val1.val()-b, val1.getBTypeNameHV(), b) + )); +} +*/ +template +BTypeName operator-(const V& a, const BTypeName& val2) +{ + return BTypeName(static_cast*>( + new BTypeNameSUB1(a-val2.val(), a, val2.getBTypeNameHV()) + )); +} +template +BTypeName operator-(const BTypeName& val1, const V& b) +{ + return BTypeName(static_cast*>( + new BTypeNameSUB2(val1.val()-b, val1.getBTypeNameHV(), b) + )); +} + +// MULTIPLICATION: + +template +struct BTypeNameMUL : public BinBTypeNameHV +{ + BTypeNameMUL(const U& val, BTypeNameHV* pOp1, BTypeNameHV* pOp2):BinBTypeNameHV(val,pOp1,pOp2){} + virtual void propagate(typename Derivatives::RecycleBin& bin) + { + this->op1()->add(bin,this->op2()->val(),this->m_derivatives); + this->op2()->add(bin,this->op1()->val(),this->m_derivatives); + } +private: + void operator=(const BTypeNameMUL&){} // not allowed +}; +template +struct BTypeNameMUL1 : public UnBTypeNameHV +{ + const V m_a; + BTypeNameMUL1(const U& val, const V& a, BTypeNameHV* pOp2):UnBTypeNameHV(val,pOp2),m_a(a){} + virtual void propagate(typename Derivatives::RecycleBin& bin) + { + this->op()->add(bin,m_a,this->m_derivatives); + } +private: + void operator=(const BTypeNameMUL1&){} // not allowed +}; +template +struct BTypeNameMUL2 : public UnBTypeNameHV +{ + const V m_b; + BTypeNameMUL2(const U& val, BTypeNameHV* pOp1, const V& b):UnBTypeNameHV(val,pOp1),m_b(b){} + virtual void propagate(typename Derivatives::RecycleBin& bin) + { + this->op()->add(bin,m_b,this->m_derivatives); + } +private: + void operator=(const BTypeNameMUL2&){} // not allowed +}; +template +BTypeName operator*(const BTypeName& val1, const BTypeName& val2) +{ + return BTypeName(static_cast*>(new BTypeNameMUL(val1.val()*val2.val(),val1.getBTypeNameHV(),val2.getBTypeNameHV()))); +} +/* +template +BTypeName operator*(const typename Op::Underlying& a, const BTypeName& val2) +{ + return BTypeName(static_cast*>( + new BTypeNameMUL1::Underlying>(a*val2.val(), a, val2.getBTypeNameHV()) + )); +} +template +BTypeName operator*(const BTypeName& val1, const typename Op::Underlying& b) +{ + return BTypeName(static_cast*>( + new BTypeNameMUL2::Underlying>(val1.val()*b, val1.getBTypeNameHV(), b) + )); +} +template +BTypeName operator*(const typename Op::Base& a, const BTypeName& val2) +{ + return BTypeName(static_cast*>( + new BTypeNameMUL1::Base>(a*val2.val(), a, val2.getBTypeNameHV()) + )); +} +template +BTypeName operator*(const BTypeName& val1, const typename Op::Base& b) +{ + return BTypeName(static_cast*>( + new BTypeNameMUL2::Base>(val1.val()*b, val1.getBTypeNameHV(), b) + )); +} +*/ +template +BTypeName operator*(const V& a, const BTypeName& val2) +{ + return BTypeName(static_cast*>( + new BTypeNameMUL1(a*val2.val(), a, val2.getBTypeNameHV()) + )); +} +template +BTypeName operator*(const BTypeName& val1, const V& b) +{ + return BTypeName(static_cast*>( + new BTypeNameMUL2(val1.val()*b, val1.getBTypeNameHV(), b) + )); +} + +// DIVISION: + +template +struct BTypeNameDIV : public BinBTypeNameHV +{ + BTypeNameDIV(const U& val, BTypeNameHV* pOp1, BTypeNameHV* pOp2):BinBTypeNameHV(val,pOp1,pOp2){} + virtual void propagate(typename Derivatives::RecycleBin& bin) + { + U tmp=Op::myInv(this->op2()->val()); + this->op1()->add(bin,tmp,this->m_derivatives); + this->op2()->sub(bin,tmp*this->val(),this->m_derivatives); + } +private: + void operator=(const BTypeNameDIV&){} // not allowed +}; +template +struct BTypeNameDIV1 : public UnBTypeNameHV +{ + const V m_a; + BTypeNameDIV1(const U& val, const V& a, BTypeNameHV* pOp2):UnBTypeNameHV(val,pOp2),m_a(a){} + virtual void propagate(typename Derivatives::RecycleBin& bin) + { + this->op()->sub(bin,Op::myInv(this->op()->val())*this->val(),this->m_derivatives); + } +private: + void operator=(const BTypeNameDIV1&){} // not allowed +}; +template +struct BTypeNameDIV2 : public UnBTypeNameHV +{ + const V m_b; + BTypeNameDIV2(const U& val, BTypeNameHV* pOp1, const V& b):UnBTypeNameHV(val,pOp1),m_b(b){} + virtual void propagate(typename Derivatives::RecycleBin& bin) + { + this->op()->add(bin,Op::myInv(m_b),this->m_derivatives); + } +private: + void operator=(const BTypeNameDIV2&){} // not allowed +}; +template +BTypeName operator/(const BTypeName& val1, const BTypeName& val2) +{ + return BTypeName(static_cast*>( + new BTypeNameDIV(val1.val()/val2.val(),val1.getBTypeNameHV(),val2.getBTypeNameHV()) + )); +} +/* +template +BTypeName operator/(const typename Op::Underlying& a, const BTypeName& val2) +{ + return BTypeName(static_cast*>( + new BTypeNameDIV1::Underlying>(a/val2.val(), a, val2.getBTypeNameHV()) + )); +} +template +BTypeName operator/(const BTypeName& val1, const typename Op::Underlying& b) +{ + return BTypeName(static_cast*>( + new BTypeNameDIV2::Underlying>(val1.val()/b, val1.getBTypeNameHV(), b) + )); +} +template +BTypeName operator/(const typename Op::Base& a, const BTypeName& val2) +{ + return BTypeName(static_cast*>( + new BTypeNameDIV1::Base>(a/val2.val(), a, val2.getBTypeNameHV()) + )); +} +template +BTypeName operator/(const BTypeName& val1, const typename Op::Base& b) +{ + return BTypeName(static_cast*>( + new BTypeNameDIV2::Base>(val1.val()/b, val1.getBTypeNameHV(), b) + )); +} +*/ +template +BTypeName operator/(const V& a, const BTypeName& val2) +{ + return BTypeName(static_cast*>( + new BTypeNameDIV1(a/val2.val(), a, val2.getBTypeNameHV()) + )); +} +template +BTypeName operator/(const BTypeName& val1, const V& b) +{ + return BTypeName(static_cast*>( + new BTypeNameDIV2(val1.val()/b, val1.getBTypeNameHV(), b) + )); +} + +// COMPOUND ASSIGNMENTS: + +template BTypeName& BTypeName::operator+=(const BTypeName& val) { return (*this)=(*this)+val; } +template BTypeName& BTypeName::operator-=(const BTypeName& val) { return (*this)=(*this)-val; } +template BTypeName& BTypeName::operator*=(const BTypeName& val) { return (*this)=(*this)*val; } +template BTypeName& BTypeName::operator/=(const BTypeName& val) { return (*this)=(*this)/val; } +template template BTypeName& BTypeName::operator+=(const V& val) { return (*this)=(*this)+val; } +template template BTypeName& BTypeName::operator-=(const V& val) { return (*this)=(*this)-val; } +template template BTypeName& BTypeName::operator*=(const V& val) { return (*this)=(*this)*val; } +template template BTypeName& BTypeName::operator/=(const V& val) { return (*this)=(*this)/val; } + +// UNARY MINUS + +template +struct BTypeNameUMINUS : public UnBTypeNameHV +{ + BTypeNameUMINUS(const U& val, BTypeNameHV* pOp):UnBTypeNameHV(val,pOp){} + virtual void propagate(typename Derivatives::RecycleBin& bin) + { + this->op()->sub(bin,this->m_derivatives); + } +private: + void operator=(const BTypeNameUMINUS&){} // not allowed +}; + +template +BTypeName operator-(const BTypeName& val) +{ + return BTypeName(static_cast*>(new BTypeNameUMINUS(Op::myNeg(val.val()),val.getBTypeNameHV()))); +} + +// UNARY PLUS + +template +struct BTypeNameUPLUS : public UnBTypeNameHV +{ + BTypeNameUPLUS(const U& val, BTypeNameHV* pOp):UnBTypeNameHV(val,pOp){} + virtual void propagate(typename Derivatives::RecycleBin& bin) + { + this->op()->add(bin,this->m_derivatives); + } +private: + void operator=(const BTypeNameUPLUS&){} // not allowed +}; + +template +BTypeName operator+(const BTypeName& val) +{ + return BTypeName(static_cast*>(new BTypeNameUPLUS(Op::myPos(val.val()),val.getBTypeNameHV()))); +} + +// POWER + +template +struct BTypeNamePOW : public BinBTypeNameHV +{ + BTypeNamePOW(const U& val, BTypeNameHV* pOp1, BTypeNameHV* pOp2):BinBTypeNameHV(val,pOp1,pOp2){} + virtual void propagate(typename Derivatives::RecycleBin& bin) + { + U tmp1(this->op2()->val() * Op::myPow(this->op1()->val(),this->op2()->val()-Op::myOne())); + U tmp2(this->val() * Op::myLog(this->op1()->val())); + this->op1()->add(bin,tmp1,this->m_derivatives); + this->op2()->add(bin,tmp2,this->m_derivatives); + } +private: + void operator=(const BTypeNamePOW&){} // not allowed +}; +template +struct BTypeNamePOW1 : public UnBTypeNameHV +{ + const V m_a; + BTypeNamePOW1(const U& val, const V& a, BTypeNameHV* pOp2):UnBTypeNameHV(val,pOp2),m_a(a){} + virtual void propagate(typename Derivatives::RecycleBin& bin) + { + U tmp2(this->val() * Op::myLog(m_a)); + this->op()->add(bin,tmp2,this->m_derivatives); + } +private: + void operator=(const BTypeNamePOW1&){} // not allowed +}; +template +struct BTypeNamePOW2 : public UnBTypeNameHV +{ + const V m_b; + BTypeNamePOW2(const U& val, BTypeNameHV* pOp1, const V& b):UnBTypeNameHV(val,pOp1),m_b(b){} + virtual void propagate(typename Derivatives::RecycleBin& bin) + { + U tmp1(m_b * Op::myPow(this->op()->val(),m_b-Op::myOne())); + this->op()->add(bin,tmp1,this->m_derivatives); + } +private: + void operator=(const BTypeNamePOW2&){} // not allowed +}; +template +BTypeName pow(const BTypeName& val1, const BTypeName& val2) +{ + return BTypeName(static_cast*>(new BTypeNamePOW(Op::myPow(val1.val(),val2.val()),val1.getBTypeNameHV(),val2.getBTypeNameHV()))); +} +/* +template +BTypeName pow(const typename Op::Underlying& a, const BTypeName& val2) +{ + return BTypeName(static_cast*>( + new BTypeNamePOW1::Underlying>(Op::myPow(a,val2.val()), a, val2.getBTypeNameHV()) + )); +} +template +BTypeName pow(const BTypeName& val1, const typename Op::Underlying& b) +{ + return BTypeName(static_cast*>( + new BTypeNamePOW2::Underlying>(Op::myPow(val1.val(),b), val1.getBTypeNameHV(), b) + )); +} +template +BTypeName pow(const typename Op::Base& a, const BTypeName& val2) +{ + return BTypeName(static_cast*>( + new BTypeNamePOW1::Base>(Op::myPow(a,val2.val()), a, val2.getBTypeNameHV()) + )); +} +template +BTypeName pow(const BTypeName& val1, const typename Op::Base& b) +{ + return BTypeName(static_cast*>( + new BTypeNamePOW2::Base>(Op::myPow(val1.val(),b), val1.getBTypeNameHV(), b) + )); +} +*/ +template +BTypeName pow(const V& a, const BTypeName& val2) +{ + return BTypeName(static_cast*>( + new BTypeNamePOW1(Op::myPow(a,val2.val()), a, val2.getBTypeNameHV()) + )); +} +template +BTypeName pow(const BTypeName& val1, const V& b) +{ + return BTypeName(static_cast*>( + new BTypeNamePOW2(Op::myPow(val1.val(),b), val1.getBTypeNameHV(), b) + )); +} + +// SQR + +template +struct BTypeNameSQR : public UnBTypeNameHV +{ + BTypeNameSQR(const U& val, BTypeNameHV* pOp):UnBTypeNameHV(val,pOp){} + virtual void propagate(typename Derivatives::RecycleBin& bin) + { + U tmp(Op::myTwo() * this->op()->val()); + this->op()->add(bin,tmp,this->m_derivatives); + } +private: + void operator=(const BTypeNameSQR&){} // not allowed +}; +template +BTypeName sqr(const BTypeName& val) +{ + return BTypeName(static_cast*>(new BTypeNameSQR(Op::mySqr(val.val()), val.getBTypeNameHV()))); +} + +// SQRT + +template +struct BTypeNameSQRT : public UnBTypeNameHV +{ + BTypeNameSQRT(const U& val, BTypeNameHV* pOp):UnBTypeNameHV(val,pOp){} + virtual void propagate(typename Derivatives::RecycleBin& bin) + { + U tmp(Op::myInv(this->val()*Op::myTwo())); + this->op()->add(bin,tmp,this->m_derivatives); + } +private: + void operator=(const BTypeNameSQRT&){} // not allowed +}; +template +BTypeName sqrt(const BTypeName& val) +{ + return BTypeName(static_cast*>(new BTypeNameSQRT(Op::mySqrt(val.val()), val.getBTypeNameHV()))); +} + +// EXP + +template +struct BTypeNameEXP : public UnBTypeNameHV +{ + BTypeNameEXP(const U& val, BTypeNameHV* pOp):UnBTypeNameHV(val,pOp){} + virtual void propagate(typename Derivatives::RecycleBin& bin) + { + this->op()->add(bin,this->val(),this->m_derivatives); + } +private: + void operator=(const BTypeNameEXP&){} // not allowed +}; +template +BTypeName exp(const BTypeName& val) +{ + return BTypeName(static_cast*>(new BTypeNameEXP(Op::myExp(val.val()), val.getBTypeNameHV()))); +} + +// LOG + +template +struct BTypeNameLOG : public UnBTypeNameHV +{ + BTypeNameLOG(const U& val, BTypeNameHV* pOp):UnBTypeNameHV(val,pOp){} + virtual void propagate(typename Derivatives::RecycleBin& bin) + { + this->op()->add(bin,Op::myInv(this->op()->val()),this->m_derivatives); + } +private: + void operator=(const BTypeNameLOG&){} // not allowed +}; +template +BTypeName log(const BTypeName& val) +{ + return BTypeName(static_cast*>(new BTypeNameLOG(Op::myLog(val.val()), val.getBTypeNameHV()))); +} + +// SIN + +template +struct BTypeNameSIN : public UnBTypeNameHV +{ + BTypeNameSIN(const U& val, BTypeNameHV* pOp):UnBTypeNameHV(val,pOp){} + virtual void propagate(typename Derivatives::RecycleBin& bin) + { + U tmp(Op::myCos(this->op()->val())); + this->op()->add(bin,tmp,this->m_derivatives); + } +private: + void operator=(const BTypeNameSIN&){} // not allowed +}; +template +BTypeName sin(const BTypeName& val) +{ + return BTypeName(static_cast*>(new BTypeNameSIN(Op::mySin(val.val()), val.getBTypeNameHV()))); +} + +// COS + +template +struct BTypeNameCOS : public UnBTypeNameHV +{ + BTypeNameCOS(const U& val, BTypeNameHV* pOp):UnBTypeNameHV(val,pOp){} + virtual void propagate(typename Derivatives::RecycleBin& bin) + { + U tmp(Op::mySin(this->op()->val())); + this->op()->sub(bin,tmp,this->m_derivatives); + } +private: + void operator=(const BTypeNameCOS&){} // not allowed +}; +template +BTypeName cos(const BTypeName& val) +{ + return BTypeName(static_cast*>(new BTypeNameCOS(Op::myCos(val.val()), val.getBTypeNameHV()))); +} + +// TAN + +template +struct BTypeNameTAN : public UnBTypeNameHV +{ + BTypeNameTAN(const U& val, BTypeNameHV* pOp):UnBTypeNameHV(val,pOp){} + virtual void propagate(typename Derivatives::RecycleBin& bin) + { + U tmp(Op::mySqr(this->val())+Op::myOne()); + this->op()->add(bin,tmp,this->m_derivatives); + } +private: + void operator=(const BTypeNameTAN&){} // not allowed +}; +template +BTypeName tan(const BTypeName& val) +{ + return BTypeName(static_cast*>(new BTypeNameTAN(Op::myTan(val.val()), val.getBTypeNameHV()))); +} + +// ASIN + +template +struct BTypeNameASIN : public UnBTypeNameHV +{ + BTypeNameASIN(const U& val, BTypeNameHV* pOp):UnBTypeNameHV(val,pOp){} + virtual void propagate(typename Derivatives::RecycleBin& bin) + { + U tmp(Op::myInv(Op::mySqrt(Op::myOne()-Op::mySqr(this->op()->val())))); + this->op()->add(bin,tmp,this->m_derivatives); + } +private: + void operator=(const BTypeNameASIN&){} // not allowed +}; +template +BTypeName asin(const BTypeName& val) +{ + return BTypeName(static_cast*>(new BTypeNameASIN(Op::myAsin(val.val()), val.getBTypeNameHV()))); +} + +// ACOS + +template +struct BTypeNameACOS : public UnBTypeNameHV +{ + BTypeNameACOS(const U& val, BTypeNameHV* pOp):UnBTypeNameHV(val,pOp){} + virtual void propagate(typename Derivatives::RecycleBin& bin) + { + U tmp(Op::myInv(Op::mySqrt(Op::myOne()-Op::mySqr(this->op()->val())))); + this->op()->sub(bin,tmp,this->m_derivatives); + } +private: + void operator=(const BTypeNameACOS&){} // not allowed +}; +template +BTypeName acos(const BTypeName& val) +{ + return BTypeName(static_cast*>(new BTypeNameACOS(Op::myAcos(val.val()), val.getBTypeNameHV()))); +} + +// ATAN + +template +struct BTypeNameATAN : public UnBTypeNameHV +{ + BTypeNameATAN(const U& val, BTypeNameHV* pOp):UnBTypeNameHV(val,pOp){} + virtual void propagate(typename Derivatives::RecycleBin& bin) + { + U tmp(Op::myInv(Op::mySqr(this->op()->val())+Op::myOne())); + this->op()->add(bin,tmp,this->m_derivatives); + } +private: + void operator=(const BTypeNameATAN&){} // not allowed +}; +template +BTypeName atan(const BTypeName& val) +{ + return BTypeName(static_cast*>(new BTypeNameATAN(Op::myAtan(val.val()), val.getBTypeNameHV()))); +} + +template struct Op< BTypeName > +{ + typedef BTypeName T; + typedef BTypeName Underlying; + typedef typename Op::Base Base; + static Base myInteger(const int i) { return Base(i); } + static Base myZero() { return myInteger(0); } + static Base myOne() { return myInteger(1);} + static Base myTwo() { return myInteger(2); } + static Base myPI() { return Op::myPI(); } + static T myPos(const T& x) { return +x; } + static T myNeg(const T& x) { return -x; } + template static T& myCadd(T& x, const V& y) { return x+=y; } + template static T& myCsub(T& x, const V& y) { return x-=y; } + template static T& myCmul(T& x, const V& y) { return x*=y; } + template static T& myCdiv(T& x, const V& y) { return x/=y; } + static T myInv(const T& x) { return myOne()/x; } + static T mySqr(const T& x) { return fadbad::sqr(x); } + template + static T myPow(const X& x, const Y& y) { return fadbad::pow(x,y); } + static T mySqrt(const T& x) { return fadbad::sqrt(x); } + static T myLog(const T& x) { return fadbad::log(x); } + static T myExp(const T& x) { return fadbad::exp(x); } + static T mySin(const T& x) { return fadbad::sin(x); } + static T myCos(const T& x) { return fadbad::cos(x); } + static T myTan(const T& x) { return fadbad::tan(x); } + static T myAsin(const T& x) { return fadbad::asin(x); } + static T myAcos(const T& x) { return fadbad::acos(x); } + static T myAtan(const T& x) { return fadbad::atan(x); } + static bool myEq(const T& x, const T& y) { return x==y; } + static bool myNe(const T& x, const T& y) { return x!=y; } + static bool myLt(const T& x, const T& y) { return xy; } + static bool myGe(const T& x, const T& y) { return x>=y; } +}; + +} //namespace fadbad + +#endif diff --git a/include/casm/clexulator/external/fadbad/fadbad.h b/include/casm/clexulator/external/fadbad/fadbad.h new file mode 100644 index 0000000..7c18c81 --- /dev/null +++ b/include/casm/clexulator/external/fadbad/fadbad.h @@ -0,0 +1,157 @@ +// Copyright (C) 1996-2007 Ole Stauning & Claus Bendtsen (fadbad@uning.dk) +// All rights reserved. + +// This code is provided "as is", without any warranty of any kind, +// either expressed or implied, including but not limited to, any implied +// warranty of merchantibility or fitness for any purpose. In no event +// will any party who distributed the code be liable for damages or for +// any claim(s) by any other party, including but not limited to, any +// lost profits, lost monies, lost data or data rendered inaccurate, +// losses sustained by third parties, or any other special, incidental or +// consequential damages arising out of the use or inability to use the +// program, even if the possibility of such damages has been advised +// against. The entire risk as to the quality, the performance, and the +// fitness of the program for any particular purpose lies with the party +// using the code. + +// This code, and any derivative of this code, may not be used in a +// commercial package without the prior explicit written permission of +// the authors. Verbatim copies of this code may be made and distributed +// in any medium, provided that this copyright notice is not removed or +// altered in any way. No fees may be charged for distribution of the +// codes, other than a fee to cover the cost of the media and a +// reasonable handling fee. + +// *************************************************************** +// ANY USE OF THIS CODE CONSTITUTES ACCEPTANCE OF THE TERMS OF THE +// COPYRIGHT NOTICE +// *************************************************************** + +#ifndef _FADBAD_H +#define _FADBAD_H + +#include + +namespace fadbad +{ + // NOTE: + // The following template allows the user to change the operations that + // are used in FADBAD++ for computing the derivatives. This is useful + // for example for specializing with non-standard types such as interval + // arithmetic types. + template struct Op // YOU MIGHT NEED TO SPECIALIZE THIS TEMPLATE: + { + typedef T Base; + static Base myInteger(const int i) { return Base(i); } + static Base myZero() { return myInteger(0); } + static Base myOne() { return myInteger(1);} + static Base myTwo() { return myInteger(2); } + static Base myPI() { return 3.14159265358979323846; } + static T myPos(const T& x) { return +x; } + static T myNeg(const T& x) { return -x; } + template static T& myCadd(T& x, const U& y) { return x+=y; } + template static T& myCsub(T& x, const U& y) { return x-=y; } + template static T& myCmul(T& x, const U& y) { return x*=y; } + template static T& myCdiv(T& x, const U& y) { return x/=y; } + static T myInv(const T& x) { return myOne()/x; } + static T mySqr(const T& x) { return x*x; } + template + static T myPow(const X& x, const Y& y) { return ::pow(x,y); } + static T mySqrt(const T& x) { return ::sqrt(x); } + static T myLog(const T& x) { return ::log(x); } + static T myExp(const T& x) { return ::exp(x); } + static T mySin(const T& x) { return ::sin(x); } + static T myCos(const T& x) { return ::cos(x); } + static T myTan(const T& x) { return ::tan(x); } + static T myAsin(const T& x) { return ::asin(x); } + static T myAcos(const T& x) { return ::acos(x); } + static T myAtan(const T& x) { return ::atan(x); } + static bool myEq(const T& x, const T& y) { return x==y; } + static bool myNe(const T& x, const T& y) { return x!=y; } + static bool myLt(const T& x, const T& y) { return xy; } + static bool myGe(const T& x, const T& y) { return x>=y; } + }; +} //namespace fadbad + +// Name for backward AD type: +#define BTypeName B + +// Name for forward AD type: +#define FTypeName F + +// Name for taylor AD type: +#define TTypeName T + +// Should always be inline: +#define INLINE0 inline + +// Methods with only one line: +#define INLINE1 inline + +// Methods with more than one line: +#define INLINE2 inline + +#ifdef __SUNPRO_CC +// FOR SOME REASON SOME INLINES CAUSES +// UNRESOLVED SMBOLS ON SUN. +#undef INLINE0 +#undef INLINE1 +#undef INLINE2 +#define INLINE0 +#define INLINE1 +#define INLINE2 +#endif + +// Define this if you want assertions, etc.. +#ifdef _DEBUG + +#include +#include + +inline void ReportError(const char* errmsg) +{ + std::cout< +class FTypeName // STACK-BASED +{ + T m_val; + T m_diff[N]; + bool m_depend; +public: + typedef T UnderlyingType; + FTypeName():m_depend(false){} + FTypeName(const FTypeName& val):m_val(val.m_val),m_depend(val.m_depend) + { + if (m_depend) for(unsigned int i=0;i /*explicit*/ FTypeName(const U& val):m_val(val),m_depend(false) + { + } + template FTypeName& operator=(const U& val) + { + m_val=val; + m_depend=false; + return *this; + } + FTypeName& operator=(const FTypeName& val) + { + if (this==&val) return *this; + m_val=val.m_val; + m_depend=val.m_depend; + if (m_depend) for(unsigned int i=0;i::myZero(); + return zero; + } + const T& d(const unsigned int i) const + { + USER_ASSERT(i::myZero(); + return zero; + } + + T& d(const unsigned int i) + { + USER_ASSERT(i::myZero(); + return zero; + } + + T& diff(unsigned int idx) + { + USER_ASSERT(idx::myZero(); + m_diff[i++]=Op::myOne(); + for( ;i::myZero(); + m_depend=true; + return m_diff[idx]; + } + bool depend() const { return m_depend; } + void setDepend(const FTypeName&) { m_depend=true; } + void setDepend(const FTypeName&, const FTypeName&) { m_depend=true; } + + FTypeName& operator+=(const FTypeName& val); + FTypeName& operator-=(const FTypeName& val); + FTypeName& operator*=(const FTypeName& val); + FTypeName& operator/=(const FTypeName& val); + template FTypeName& operator+=(const V& val); + template FTypeName& operator-=(const V& val); + template FTypeName& operator*=(const V& val); + template FTypeName& operator/=(const V& val); + +}; + +template +class FTypeName // HEAP-BASED +{ + T m_val; + unsigned int m_size; + T* m_diff; +public: + typedef T UnderlyingType; + FTypeName():m_val(),m_size(0),m_diff(0){} + FTypeName(const FTypeName& val):m_val(val.m_val),m_size(val.m_size),m_diff(m_size==0?0:new T[m_size]) + { + for(unsigned int i=0;i /*explicit*/ FTypeName(const U& val):m_val(val),m_size(0),m_diff(0){} + ~FTypeName(){ delete[] m_diff; } + template FTypeName& operator=(const U& val) + { + m_val=val; + m_size=0; + delete[] m_diff; + m_diff=0; + return *this; + } + FTypeName& operator=(const FTypeName& val) + { + if (this==&val) return *this; + m_val=val.m_val; + if (val.m_size>0) + { + if (m_size==0) + { + m_size = val.m_size; + m_diff = new T[m_size]; + } + USER_ASSERT(m_size==val.m_size,"derivative vectors not of same size"); + for(unsigned int i=0;i0) + { + for(unsigned int i=0;i::myZero(); + return zero; + } + + const T& d(const unsigned int i) const + { + if (i::myZero(); + return zero; + } + + T& d(const unsigned int i) + { + if (i::myZero(); + return zero; + } + + T& diff(unsigned int idx, unsigned int N) + { + USER_ASSERT(idx::myZero(); + m_diff[i++]=Op::myOne(); + for( ;i::myZero(); + return m_diff[idx]; + } + bool depend() const { return m_size!=0; } + void setDepend(const FTypeName& val) + { + INTERNAL_ASSERT(val.m_size>0,"input is not a dependent variable") + if (m_size==0) + { + m_size = val.m_size; + m_diff = new T[m_size]; + } + else + { + USER_ASSERT(m_size==val.m_size,"derivative vectors not of same size "<& val1, const FTypeName& val2) + { + USER_ASSERT(val1.m_size==val2.m_size,"derivative vectors not of same size "<0,"lhs-input is not a dependent variable") + INTERNAL_ASSERT(val2.m_size>0,"rhs-input is not a dependent variable") + if (m_size==0) + { + m_size=val1.m_size; + m_diff = new T[m_size]; + } + else + { + USER_ASSERT(m_size==val1.m_size,"derivative vectors not of same size "<& operator+=(const FTypeName& val); + FTypeName& operator-=(const FTypeName& val); + FTypeName& operator*=(const FTypeName& val); + FTypeName& operator/=(const FTypeName& val); + template FTypeName& operator+=(const V& val); + template FTypeName& operator-=(const V& val); + template FTypeName& operator*=(const V& val); + template FTypeName& operator/=(const V& val); + +}; + +template bool operator==(const FTypeName& val1, const FTypeName& val2) { return Op::myEq(val1.val(),val2.val()); } +template bool operator!=(const FTypeName& val1, const FTypeName& val2) { return Op::myNe(val1.val(),val2.val()); } +template bool operator<(const FTypeName& val1, const FTypeName& val2) { return Op::myLt(val1.val(),val2.val()); } +template bool operator<=(const FTypeName& val1, const FTypeName& val2) { return Op::myLe(val1.val(),val2.val()); } +template bool operator>(const FTypeName& val1, const FTypeName& val2) { return Op::myGt(val1.val(),val2.val()); } +template bool operator>=(const FTypeName& val1, const FTypeName& val2) { return Op::myGe(val1.val(),val2.val()); } +template bool operator==(const FTypeName& val1, const U& val2) { return Op::myEq(val1.val(),val2); } +template bool operator==(const U& val1, const FTypeName& val2) { return Op::myEq(val1,val2.val()); } +template bool operator!=(const FTypeName& val1, const U& val2) { return Op::myNe(val1.val(),val2); } +template bool operator!=(const U& val1, const FTypeName& val2) { return Op::myNe(val1,val2.val()); } +template bool operator<(const FTypeName& val1, const U& val2) { return Op::myLt(val1.val(),val2); } +template bool operator<(const U& val1, const FTypeName& val2) { return Op::myLt(val1,val2.val()); } +template bool operator<=(const FTypeName& val1, const U& val2) { return Op::myLe(val1.val(),val2); } +template bool operator<=(const U& val1, const FTypeName& val2) { return Op::myLe(val1,val2.val()); } +template bool operator>(const FTypeName& val1, const U& val2) { return Op::myGt(val1.val(),val2); } +template bool operator>(const U& val1, const FTypeName& val2) { return Op::myGt(val1,val2.val()); } +template bool operator>=(const FTypeName& val1, const U& val2) { return Op::myGe(val1.val(),val2); } +template bool operator>=(const U& val1, const FTypeName& val2) { return Op::myGe(val1,val2.val()); } + +template +INLINE2 FTypeName add1(const U& a, const FTypeName& b) +{ + FTypeName c(a+b.val()); + if (!b.depend()) return c; + c.setDepend(b); + for(unsigned int i=0;i +INLINE2 FTypeName add1(const U& a, const FTypeName& b) +{ + FTypeName c(a+b.val()); + if (!b.depend()) return c; + c.setDepend(b); + for(unsigned int i=0;i +INLINE2 FTypeName add2(const FTypeName& a, const U& b) +{ + FTypeName c(a.val()+b); + if (!a.depend()) return c; + c.setDepend(a); + for(unsigned int i=0;i +INLINE2 FTypeName add2(const FTypeName& a, const U& b) +{ + FTypeName c(a.val()+b); + if (!a.depend()) return c; + c.setDepend(a); + for(unsigned int i=0;i +INLINE2 FTypeName operator+ (const U& a, const FTypeName& b) +{ + return add1(a,b); +} + +template +INLINE2 FTypeName operator+ (const FTypeName& a, const U& b) +{ + return add2(a,b); +} + +template +INLINE2 FTypeName add3(const FTypeName& a, const FTypeName& b) +{ + FTypeName c(a.val()+b.val()); + c.setDepend(a,b); + for(unsigned int i=0;i +INLINE2 FTypeName add3(const FTypeName& a, const FTypeName& b) +{ + FTypeName c(a.val()+b.val()); + c.setDepend(a,b); + for(unsigned int i=0;i +INLINE2 FTypeName operator+ (const FTypeName& a, const FTypeName& b) +{ + switch ((a.depend()?1:0)|(b.depend()?2:0)) + { + case 0: return FTypeName(a.val()+b.val()); + case 1: return add2(a,b.val()); + case 2: return add1(a.val(),b); + } + return add3(a,b); +} + +template +INLINE2 FTypeName sub1(const U& a, const FTypeName& b) +{ + FTypeName c(a-b.val()); + if (!b.depend()) return c; + c.setDepend(b); + for(unsigned int i=0;i::myNeg(b[i]); + return c; +} +template +INLINE2 FTypeName sub1(const U& a, const FTypeName& b) +{ + FTypeName c(a-b.val()); + if (!b.depend()) return c; + c.setDepend(b); + for(unsigned int i=0;i::myNeg(b[i]); + return c; +} + +template +INLINE2 FTypeName sub2(const FTypeName& a, const U& b) +{ + FTypeName c(a.val()-b); + if (!a.depend()) return c; + c.setDepend(a); + for(unsigned int i=0;i +INLINE2 FTypeName sub2(const FTypeName& a, const U& b) +{ + FTypeName c(a.val()-b); + if (!a.depend()) return c; + c.setDepend(a); + for(unsigned int i=0;i +INLINE2 FTypeName operator- (const U& a, const FTypeName& b) +{ + return sub1(a,b); +} + +template +INLINE2 FTypeName operator- (const FTypeName& a, const U& b) +{ + return sub2(a,b); +} + +template +INLINE2 FTypeName sub3(const FTypeName& a, const FTypeName& b) +{ + FTypeName c(a.val()-b.val()); + c.setDepend(a,b); + for(unsigned int i=0;i +INLINE2 FTypeName sub3(const FTypeName& a, const FTypeName& b) +{ + FTypeName c(a.val()-b.val()); + c.setDepend(a,b); + for(unsigned int i=0;i +INLINE2 FTypeName operator- (const FTypeName& a, const FTypeName& b) +{ + switch ((a.depend()?1:0)|(b.depend()?2:0)) + { + case 0: return FTypeName(a.val()-b.val()); + case 1: return sub2(a,b.val()); + case 2: return sub1(a.val(),b); + } + return sub3(a,b); +} + +template +INLINE2 FTypeName mul1(const U& a, const FTypeName& b) +{ + FTypeName c(a*b.val()); + if (!b.depend()) return c; + c.setDepend(b); + for(unsigned int i=0;i +INLINE2 FTypeName mul1(const U& a, const FTypeName& b) +{ + FTypeName c(a*b.val()); + if (!b.depend()) return c; + c.setDepend(b); + for(unsigned int i=0;i +INLINE2 FTypeName mul2(const FTypeName& a, const U& b) +{ + FTypeName c(a.val()*b); + if (!a.depend()) return c; + c.setDepend(a); + for(unsigned int i=0;i +INLINE2 FTypeName mul2(const FTypeName& a, const U& b) +{ + FTypeName c(a.val()*b); + if (!a.depend()) return c; + c.setDepend(a); + for(unsigned int i=0;i +INLINE2 FTypeName operator* (const U& a, const FTypeName& b) +{ + return mul1(a,b); +} + +template +INLINE2 FTypeName operator* (const FTypeName& a, const U& b) +{ + return mul2(a,b); +} + +template +INLINE2 FTypeName mul3 (const FTypeName& a, const FTypeName& b) +{ + const T& aval(a.val()); + const T& bval(b.val()); + FTypeName c(aval*bval); + c.setDepend(a,b); + for(unsigned int i=0;i +INLINE2 FTypeName mul3 (const FTypeName& a, const FTypeName& b) +{ + const T& aval(a.val()); + const T& bval(b.val()); + FTypeName c(aval*bval); + c.setDepend(a,b); + for(unsigned int i=0;i +INLINE2 FTypeName operator* (const FTypeName& a, const FTypeName& b) +{ + switch ((a.depend()?1:0)|(b.depend()?2:0)) + { + case 0: return FTypeName(a.val()*b.val()); + case 1: return mul2(a,b.val()); + case 2: return mul1(a.val(),b); + } + return mul3(a,b); +} + +template +INLINE2 FTypeName div1(const U& a, const FTypeName& b) +{ + FTypeName c(a/b.val()); + if (!b.depend()) return c; + T tmp(Op::myNeg(c.val()/b.val())); + c.setDepend(b); + for(unsigned int i=0;i +INLINE2 FTypeName div1(const U& a, const FTypeName& b) +{ + FTypeName c(a/b.val()); + if (!b.depend()) return c; + T tmp(Op::myNeg(c.val()/b.val())); + c.setDepend(b); + for(unsigned int i=0;i +INLINE2 FTypeName div2(const FTypeName& a, const U& b) +{ + FTypeName c(a.val()/b); + if (!a.depend()) return c; + c.setDepend(a); + for(unsigned int i=0;i +INLINE2 FTypeName div2(const FTypeName& a, const U& b) +{ + FTypeName c(a.val()/b); + if (!a.depend()) return c; + c.setDepend(a); + for(unsigned int i=0;i +INLINE2 FTypeName operator/ (const U& a, const FTypeName& b) +{ + return div1(a,b); +} + +template +INLINE2 FTypeName operator/ (const FTypeName& a, const U& b) +{ + return div2(a,b); +} + +template +INLINE2 FTypeName div3(const FTypeName& a, const FTypeName& b) +{ + const T& bval(b.val()); + FTypeName c(a.val()/bval); + c.setDepend(a,b); + const T& cval(c.val()); + for(unsigned int i=0;i +INLINE2 FTypeName div3(const FTypeName& a, const FTypeName& b) +{ + const T& bval(b.val()); + FTypeName c(a.val()/bval); + c.setDepend(a,b); + const T& cval(c.val()); + for(unsigned int i=0;i +INLINE2 FTypeName operator/ (const FTypeName& a, const FTypeName& b) +{ + switch ((a.depend()?1:0)|(b.depend()?2:0)) + { + case 0: return FTypeName(a.val()/b.val()); + case 1: return div2(a,b.val()); + case 2: return div1(a.val(),b); + } + return div3(a,b); +} + +template +FTypeName& FTypeName::operator+=(const FTypeName& val) +{ + Op::myCadd(m_val,val.m_val); + if (!val.depend()) return *this; + if (this->depend()) + { + for(unsigned int i=0;i::myCadd(m_diff[i],val[i]); + } + else + { + this->setDepend(val); + for(unsigned int i=0;i +FTypeName& FTypeName::operator-=(const FTypeName& val) +{ + Op::myCsub(m_val,val.m_val); + if (!val.depend()) return *this; + if (this->depend()) + { + for(unsigned int i=0;i::myCsub(m_diff[i],val[i]); + } + else + { + this->setDepend(val); + for(unsigned int i=0;i::myNeg(val[i]); + } + return *this; +} + +template +FTypeName& FTypeName::operator*=(const FTypeName& val) +{ + if (this->depend() && val.depend()) + { + for(unsigned int i=0;idepend()) + { + for(unsigned int i=0;i::myCmul(m_diff[i],val.m_val); + } + else // (val.depend()) + { + this->setDepend(val); + for(unsigned int i=0;i::myCmul(m_val,val.m_val); + return *this; +} + +template +FTypeName& FTypeName::operator/=(const FTypeName& val) +{ + Op::myCdiv(m_val,val.m_val); + if (this->depend() && val.depend()) + { + for(unsigned int i=0;idepend()) + { + for(unsigned int i=0;i::myCdiv(m_diff[i],val.m_val); + } + else // (val.depend()) + { + this->setDepend(val); + for(unsigned int i=0;i::myNeg(m_val*val.m_diff[i]/val.m_val); + } + return *this; +} + +template +template FTypeName& FTypeName::operator+=(const V& val) +{ + Op::myCadd(m_val,val); + return *this; +} +template +template FTypeName& FTypeName::operator-=(const V& val) +{ + Op::myCsub(m_val,val); + return *this; +} + +template +template FTypeName& FTypeName::operator*=(const V& val) +{ + Op::myCmul(m_val,val); + if (!this->depend()) return *this; + for(unsigned int i=0;i::myCmul(m_diff[i],val); + return *this; +} + +template +template FTypeName& FTypeName::operator/=(const V& val) +{ + Op::myCdiv(m_val,val); + if (!this->depend()) return *this; + for(unsigned int i=0;i::myCdiv(m_diff[i],val); + return *this; +} + +template +FTypeName& FTypeName::operator+=(const FTypeName& val) +{ + Op::myCadd(m_val,val.m_val); + if (!val.depend()) return *this; + if (this->depend()) + { + for(unsigned int i=0;isize();++i) Op::myCadd(m_diff[i],val[i]); + } + else + { + this->setDepend(val); + for(unsigned int i=0;isize();++i) m_diff[i]=val[i]; + } + return *this; +} +template +FTypeName& FTypeName::operator-=(const FTypeName& val) +{ + Op::myCsub(m_val,val.m_val); + if (!val.depend()) return *this; + if (this->depend()) + { + for(unsigned int i=0;isize();++i) Op::myCsub(m_diff[i],val[i]); + } + else + { + this->setDepend(val); + for(unsigned int i=0;isize();++i) m_diff[i]=Op::myNeg(val[i]); + } + return *this; +} +template +FTypeName& FTypeName::operator*=(const FTypeName& val) +{ + if (this->depend() && val.depend()) + { + for(unsigned int i=0;isize();++i) m_diff[i]=m_diff[i]*val.m_val+val.m_diff[i]*m_val; + } + else if (this->depend()) + { + for(unsigned int i=0;isize();++i) Op::myCmul(m_diff[i],val.m_val); + } + else // (val.depend()) + { + this->setDepend(val); + for(unsigned int i=0;isize();++i) m_diff[i]=val.m_diff[i]*m_val; + } + Op::myCmul(m_val,val.m_val); + return *this; +} +template +FTypeName& FTypeName::operator/=(const FTypeName& val) +{ + Op::myCdiv(m_val,val.m_val); + if (this->depend() && val.depend()) + { + for(unsigned int i=0;isize();++i) m_diff[i]=(m_diff[i]-m_val*val.m_diff[i])/val.m_val; + } + else if (this->depend()) + { + for(unsigned int i=0;isize();++i) Op::myCdiv(m_diff[i],val.m_val); + } + else // (val.depend()) + { + this->setDepend(val); + for(unsigned int i=0;isize();++i) m_diff[i]=Op::myNeg(m_val*val.m_diff[i]/val.m_val); + } + return *this; +} +template +template FTypeName& FTypeName::operator+=(const V& val) +{ + Op::myCadd(m_val,val); + return *this; +} +template +template FTypeName& FTypeName::operator-=(const V& val) +{ + Op::myCsub(m_val,val); + return *this; +} +template +template FTypeName& FTypeName::operator*=(const V& val) +{ + Op::myCmul(m_val,val); + if (!this->depend()) return *this; + for(unsigned int i=0;isize();++i) Op::myCmul(m_diff[i],val); + return *this; +} +template +template FTypeName& FTypeName::operator/=(const V& val) +{ + Op::myCdiv(m_val,val); + if (!this->depend()) return *this; + for(unsigned int i=0;isize();++i) Op::myCdiv(m_diff[i],val); + return *this; +} + +template +INLINE2 FTypeName pow1(const U& a, const FTypeName& b) +{ + FTypeName c(Op::myPow(a,b.val())); + if (!b.depend()) return c; + T tmp(c.val()*Op::myLog(a)); + c.setDepend(b); + for(unsigned int i=0;i +INLINE2 FTypeName pow1(const U& a, const FTypeName& b) +{ + FTypeName c(Op::myPow(a,b.val())); + if (!b.depend()) return c; + T tmp(c.val()*Op::myLog(a)); + c.setDepend(b); + for(unsigned int i=0;i +INLINE2 FTypeName pow2(const FTypeName& a, const U& b) +{ + FTypeName c(Op::myPow(a.val(),b)); + if (!a.depend()) return c; + T tmp(b*Op::myPow(a.val(),b-Op::myOne())); + c.setDepend(a); + for(unsigned int i=0;i +INLINE2 FTypeName pow2(const FTypeName& a, const U& b) +{ + FTypeName c(Op::myPow(a.val(),b)); + if (!a.depend()) return c; + T tmp(b*Op::myPow(a.val(),b-Op::myOne())); + c.setDepend(a); + for(unsigned int i=0;i +INLINE2 FTypeName pow (const U& a, const FTypeName& b) +{ + return pow1(a,b); +} + +template +INLINE2 FTypeName pow (const FTypeName& a,const U& b) +{ + return pow2(a,b); +} + +template +INLINE2 FTypeName pow3(const FTypeName& a, const FTypeName& b) +{ + FTypeName c(Op::myPow(a.val(),b.val())); + T tmp(b.val()*Op::myPow(a.val(),b.val()-Op::myOne())),tmp1(c.val()*Op::myLog(a.val())); + c.setDepend(a,b); + for(unsigned int i=0;i +INLINE2 FTypeName pow3(const FTypeName& a, const FTypeName& b) +{ + FTypeName c(Op::myPow(a.val(),b.val())); + T tmp(b.val()*Op::myPow(a.val(),b.val()-Op::myOne())),tmp1(c.val()*Op::myLog(a.val())); + c.setDepend(a,b); + for(unsigned int i=0;i +INLINE2 FTypeName pow (const FTypeName& a, const FTypeName& b) +{ + switch ((a.depend()?1:0)|(b.depend()?2:0)) + { + case 0: return FTypeName(Op::myPow(a.val(),b.val())); + case 1: return pow2(a,b.val()); + case 2: return pow1(a.val(),b); + } + return pow3(a,b); +} + +/* Unary operators */ +template +INLINE2 FTypeName operator+ (const FTypeName& a) +{ + FTypeName c(a.val()); + if (!a.depend()) return c; + c.setDepend(a); + for(unsigned int i=0;i +INLINE2 FTypeName operator+ (const FTypeName& a) +{ + FTypeName c(a.val()); + if (!a.depend()) return c; + c.setDepend(a); + for(unsigned int i=0;i +INLINE2 FTypeName operator- (const FTypeName& a) +{ + FTypeName c(Op::myNeg(a.val())); + if (!a.depend()) return c; + c.setDepend(a); + for(unsigned int i=0;i::myNeg(a[i]); + return c; +} +template +INLINE2 FTypeName operator- (const FTypeName& a) +{ + FTypeName c(Op::myNeg(a.val())); + if (!a.depend()) return c; + c.setDepend(a); + for(unsigned int i=0;i::myNeg(a[i]); + return c; +} + +template +INLINE2 FTypeName sqr (const FTypeName& a) +{ + FTypeName c(Op::mySqr(a.val())); + if (!a.depend()) return c; + T tmp(Op::myTwo()*a.val()); + c.setDepend(a); + for(unsigned int i=0;i +INLINE2 FTypeName sqr (const FTypeName& a) +{ + FTypeName c(Op::mySqr(a.val())); + if (!a.depend()) return c; + T tmp(Op::myTwo()*a.val()); + c.setDepend(a); + for(unsigned int i=0;i +INLINE2 FTypeName exp (const FTypeName& a) +{ + FTypeName c(Op::myExp(a.val())); + if (!a.depend()) return c; + c.setDepend(a); + const T& cval(c.val()); + for(unsigned int i=0;i +INLINE2 FTypeName exp (const FTypeName& a) +{ + FTypeName c(Op::myExp(a.val())); + if (!a.depend()) return c; + c.setDepend(a); + const T& cval(c.val()); + for(unsigned int i=0;i +INLINE2 FTypeName log (const FTypeName& a) +{ + FTypeName c(Op::myLog(a.val())); + if (!a.depend()) return c; + c.setDepend(a); + const T& aval(a.val()); + for(unsigned int i=0;i +INLINE2 FTypeName log (const FTypeName& a) +{ + FTypeName c(Op::myLog(a.val())); + if (!a.depend()) return c; + c.setDepend(a); + const T& aval(a.val()); + for(unsigned int i=0;i +INLINE2 FTypeName sqrt (const FTypeName& a) +{ + FTypeName c(Op::mySqrt(a.val())); + if (!a.depend()) return c; + T tmp(c.val()*Op::myTwo()); + c.setDepend(a); + for(unsigned int i=0;i +INLINE2 FTypeName sqrt (const FTypeName& a) +{ + FTypeName c(Op::mySqrt(a.val())); + if (!a.depend()) return c; + T tmp(c.val()*Op::myTwo()); + c.setDepend(a); + for(unsigned int i=0;i +INLINE2 FTypeName sin (const FTypeName& a) +{ + FTypeName c(Op::mySin(a.val())); + if (!a.depend()) return c; + T tmp(Op::myCos(a.val())); + c.setDepend(a); + for(unsigned int i=0;i +INLINE2 FTypeName sin (const FTypeName& a) +{ + FTypeName c(Op::mySin(a.val())); + if (!a.depend()) return c; + T tmp(Op::myCos(a.val())); + c.setDepend(a); + for(unsigned int i=0;i +INLINE2 FTypeName cos (const FTypeName& a) +{ + FTypeName c(Op::myCos(a.val())); + if (!a.depend()) return c; + T tmp(-Op::mySin(a.val())); + c.setDepend(a); + for(unsigned int i=0;i +INLINE2 FTypeName cos (const FTypeName& a) +{ + FTypeName c(Op::myCos(a.val())); + if (!a.depend()) return c; + T tmp(-Op::mySin(a.val())); + c.setDepend(a); + for(unsigned int i=0;i +INLINE2 FTypeName tan (const FTypeName& a) +{ + FTypeName c(Op::myTan(a.val())); + if (!a.depend()) return c; + T tmp(Op::myOne()+Op::mySqr(c.val())); + c.setDepend(a); + for(unsigned int i=0;i +INLINE2 FTypeName tan (const FTypeName& a) +{ + FTypeName c(Op::myTan(a.val())); + if (!a.depend()) return c; + T tmp(Op::myOne()+Op::mySqr(c.val())); + c.setDepend(a); + for(unsigned int i=0;i +INLINE2 FTypeName asin (const FTypeName& a) +{ + FTypeName c(Op::myAsin(a.val())); + if (!a.depend()) return c; + T tmp(Op::myInv(Op::mySqrt(Op::myOne()-Op::mySqr(a.val())))); + c.setDepend(a); + for(unsigned int i=0;i +INLINE2 FTypeName asin (const FTypeName& a) +{ + FTypeName c(Op::myAsin(a.val())); + if (!a.depend()) return c; + T tmp(Op::myInv(Op::mySqrt(Op::myOne()-Op::mySqr(a.val())))); + c.setDepend(a); + for(unsigned int i=0;i +INLINE2 FTypeName acos (const FTypeName& a) +{ + FTypeName c(Op::myAcos(a.val())); + if (!a.depend()) return c; + T tmp(Op::myNeg(Op::myInv(Op::mySqrt(Op::myOne()-Op::mySqr(a.val()))))); + c.setDepend(a); + for(unsigned int i=0;i +INLINE2 FTypeName acos (const FTypeName& a) +{ + FTypeName c(Op::myAcos(a.val())); + if (!a.depend()) return c; + T tmp(Op::myNeg(Op::myInv(Op::mySqrt(Op::myOne()-Op::mySqr(a.val()))))); + c.setDepend(a); + for(unsigned int i=0;i +INLINE2 FTypeName atan (const FTypeName& a) +{ + FTypeName c(Op::myAtan(a.val())); + if (!a.depend()) return c; + T tmp(Op::myInv(Op::myOne()+Op::mySqr(a.val()))); + c.setDepend(a); + for(unsigned int i=0;i +INLINE2 FTypeName atan (const FTypeName& a) +{ + FTypeName c(Op::myAtan(a.val())); + if (!a.depend()) return c; + T tmp(Op::myInv(Op::myOne()+Op::mySqr(a.val()))); + c.setDepend(a); + for(unsigned int i=0;i struct Op< FTypeName > +{ + typedef FTypeName T; + typedef FTypeName Underlying; + typedef typename Op::Base Base; + static Base myInteger(const int i) { return Base(i); } + static Base myZero() { return myInteger(0); } + static Base myOne() { return myInteger(1);} + static Base myTwo() { return myInteger(2); } + static Base myPI() { return Op::myPI(); } + static T myPos(const T& x) { return +x; } + static T myNeg(const T& x) { return -x; } + template static T& myCadd(T& x, const V& y) { return x+=y; } + template static T& myCsub(T& x, const V& y) { return x-=y; } + template static T& myCmul(T& x, const V& y) { return x*=y; } + template static T& myCdiv(T& x, const V& y) { return x/=y; } + static T myInv(const T& x) { return myOne()/x; } + static T mySqr(const T& x) { return fadbad::sqr(x); } + template + static T myPow(const X& x, const Y& y) { return fadbad::pow(x,y); } + static T mySqrt(const T& x) { return fadbad::sqrt(x); } + static T myLog(const T& x) { return fadbad::log(x); } + static T myExp(const T& x) { return fadbad::exp(x); } + static T mySin(const T& x) { return fadbad::sin(x); } + static T myCos(const T& x) { return fadbad::cos(x); } + static T myTan(const T& x) { return fadbad::tan(x); } + static T myAsin(const T& x) { return fadbad::asin(x); } + static T myAcos(const T& x) { return fadbad::acos(x); } + static T myAtan(const T& x) { return fadbad::atan(x); } + static bool myEq(const T& x, const T& y) { return x==y; } + static bool myNe(const T& x, const T& y) { return x!=y; } + static bool myLt(const T& x, const T& y) { return xy; } + static bool myGe(const T& x, const T& y) { return x>=y; } +}; + +} // namespace fadbad + +#endif diff --git a/include/casm/clexulator/external/fadbad/tadiff.h b/include/casm/clexulator/external/fadbad/tadiff.h new file mode 100644 index 0000000..e1bf8d0 --- /dev/null +++ b/include/casm/clexulator/external/fadbad/tadiff.h @@ -0,0 +1,1283 @@ +// Copyright (C) 1996-2007 Ole Stauning & Claus Bendtsen (fadbad@uning.dk) +// All rights reserved. + +// This code is provided "as is", without any warranty of any kind, +// either expressed or implied, including but not limited to, any implied +// warranty of merchantibility or fitness for any purpose. In no event +// will any party who distributed the code be liable for damages or for +// any claim(s) by any other party, including but not limited to, any +// lost profits, lost monies, lost data or data rendered inaccurate, +// losses sustained by third parties, or any other special, incidental or +// consequential damages arising out of the use or inability to use the +// program, even if the possibility of such damages has been advised +// against. The entire risk as to the quality, the performance, and the +// fitness of the program for any particular purpose lies with the party +// using the code. + +// This code, and any derivative of this code, may not be used in a +// commercial package without the prior explicit written permission of +// the authors. Verbatim copies of this code may be made and distributed +// in any medium, provided that this copyright notice is not removed or +// altered in any way. No fees may be charged for distribution of the +// codes, other than a fee to cover the cost of the media and a +// reasonable handling fee. + +// *************************************************************** +// ANY USE OF THIS CODE CONSTITUTES ACCEPTANCE OF THE TERMS OF THE +// COPYRIGHT NOTICE +// *************************************************************** + +#ifndef _TADIFF_H +#define _TADIFF_H + +#include + +#ifndef MaxLength +#define MaxLength 40 +#endif + +#include "fadbad.h" + +namespace fadbad +{ + +template +class TValues +{ + unsigned int m_n; + U m_val[N]; +public: + TValues():m_n(0){std::fill(m_val,m_val+N,Op::myZero());} + template explicit TValues(const V& val):m_n(1){m_val[0]=val;std::fill(m_val+1,m_val+N,Op::myZero());} + U& operator[](const unsigned int i) + { + USER_ASSERT(i +class TTypeNameHV // Heap Value +{ + TValues m_val; + mutable unsigned int m_rc; +protected: + virtual ~TTypeNameHV(){} +public: + TTypeNameHV():m_rc(0){} + template explicit TTypeNameHV(const V& val):m_val(val),m_rc(0){} + const U& val(const unsigned int i) const { return m_val[i]; } + U& val(const unsigned int i) { return m_val[i]; } + unsigned int length() const { return m_val.length(); } + unsigned int& length() { return m_val.length(); } + void decRef(TTypeNameHV*& pTTypeNameHV) const { if (--m_rc==0) { delete this; pTTypeNameHV=0;} } + void incRef() const {++m_rc;} + + virtual void reset(){m_val.reset();} + virtual unsigned int eval(const unsigned int k){return k+1;} +}; + +template +class TTypeName +{ +private: + struct SV // Stack Value refers to reference-counted Heap Value: + { + mutable TTypeNameHV* m_pTTypeNameHV; + SV(TTypeNameHV* pTTypeNameHV):m_pTTypeNameHV(pTTypeNameHV){ m_pTTypeNameHV->incRef(); } + SV(const typename TTypeName::SV& sv):m_pTTypeNameHV(sv.m_pTTypeNameHV){ m_pTTypeNameHV->incRef(); } + ~SV(){ m_pTTypeNameHV->decRef(m_pTTypeNameHV); } + TTypeNameHV* getTTypeNameHV() const { return m_pTTypeNameHV; } + void setTTypeNameHV(TTypeNameHV* pTTypeNameHV) + { + if (m_pTTypeNameHV!=pTTypeNameHV) + { + m_pTTypeNameHV->decRef(m_pTTypeNameHV); + m_pTTypeNameHV=pTTypeNameHV; + m_pTTypeNameHV->incRef(); + } + } + const U& val() const { return m_pTTypeNameHV->val(0); } + const U& val(const unsigned int i) const { return m_pTTypeNameHV->val(i); } + unsigned int length() const { return m_pTTypeNameHV->length(); } + unsigned int& length() { return m_pTTypeNameHV->length(); } + U& val(const unsigned int i) { return m_pTTypeNameHV->val(i); } + + void reset(){m_pTTypeNameHV->reset();} + unsigned int eval(const unsigned int i){return m_pTTypeNameHV->eval(i);} + } m_sv; +public: + typedef U UnderlyingType; + TTypeName():m_sv(new TTypeNameHV()){} + TTypeName(TTypeNameHV* pTTypeNameHV):m_sv(pTTypeNameHV){} + explicit TTypeName(const typename TTypeName::SV& sv):m_sv(sv){} + template /*explicit*/ TTypeName(const V& val):m_sv(new TTypeNameHV(val)){m_sv.length()=N;} + TTypeName& operator=(const TTypeName& val) + { + if (this==&val) return *this; + m_sv.setTTypeNameHV(val.m_sv.getTTypeNameHV()); + return *this; + } + template TTypeName& operator=(const V& val) + { + m_sv.setTTypeNameHV(new TTypeNameHV(val)); + m_sv.length()=N; + return *this; + } + TTypeNameHV* getTTypeNameHV() const { return m_sv.getTTypeNameHV(); } + void setTTypeNameHV(const TTypeNameHV* pTTypeNameHV) { m_sv.setTTypeNameHV(pTTypeNameHV); } + const U& val() const { return m_sv.val(); } + unsigned int length() const { return m_sv.length(); } + const U& operator[](const unsigned int i) const { return m_sv.val(i); } + U& operator[](const unsigned int i) { if (i>=m_sv.length()) m_sv.length()=i+1; return m_sv.val(i);} + + TTypeName& operator+=(const TTypeName& val); + TTypeName& operator-=(const TTypeName& val); + TTypeName& operator*=(const TTypeName& val); + TTypeName& operator/=(const TTypeName& val); + template TTypeName& operator+=(const V& val); + template TTypeName& operator-=(const V& val); + template TTypeName& operator*=(const V& val); + template TTypeName& operator/=(const V& val); + + void reset(){m_sv.reset();} + unsigned int eval(const unsigned int i){return m_sv.eval(i);} +}; + +template bool operator==(const TTypeName& val1, const TTypeName& val2) { return Op::myEq(val1.val(),val2.val()); } +template bool operator!=(const TTypeName& val1, const TTypeName& val2) { return Op::myNe(val1.val(),val2.val()); } +template bool operator<(const TTypeName& val1, const TTypeName& val2) { return Op::myLt(val1.val(),val2.val()); } +template bool operator<=(const TTypeName& val1, const TTypeName& val2) { return Op::myLe(val1.val(),val2.val()); } +template bool operator>(const TTypeName& val1, const TTypeName& val2) { return Op::myGt(val1.val(),val2.val()); } +template bool operator>=(const TTypeName& val1, const TTypeName& val2) { return Op::myGe(val1.val(),val2.val()); } +template bool operator==(const TTypeName& val1, const V& val2) { return Op::myEq(val1.val(),val2); } +template bool operator==(const V& val1, const TTypeName& val2) { return Op::myEq(val1,val2.val()); } +template bool operator!=(const TTypeName& val1, const V& val2) { return Op::myNe(val1.val(),val2); } +template bool operator!=(const V& val1, const TTypeName& val2) { return Op::myNe(val1,val2.val()); } +template bool operator<(const TTypeName& val1, const V& val2) { return Op::myLt(val1.val(),val2); } +template bool operator<(const V& val1, const TTypeName& val2) { return Op::myLt(val1,val2.val()); } +template bool operator<=(const TTypeName& val1, const V& val2) { return Op::myLe(val1.val(),val2); } +template bool operator<=(const V& val1, const TTypeName& val2) { return Op::myLe(val1,val2.val()); } +template bool operator>(const TTypeName& val1, const V& val2) { return Op::myGt(val1.val(),val2); } +template bool operator>(const V& val1, const TTypeName& val2) { return Op::myGt(val1,val2.val()); } +template bool operator>=(const TTypeName& val1, const V& val2) { return Op::myGe(val1.val(),val2); } +template bool operator>=(const V& val1, const TTypeName& val2) { return Op::myGe(val1,val2.val()); } + +// Binary operator base class: + +template +class BinTTypeNameHV : public TTypeNameHV +{ + TTypeNameHV* m_pOp1; + TTypeNameHV* m_pOp2; +public: + BinTTypeNameHV(const U& val, TTypeNameHV* pOp1, TTypeNameHV* pOp2):TTypeNameHV(val),m_pOp1(pOp1),m_pOp2(pOp2) + { + m_pOp1->incRef();m_pOp2->incRef(); + } + BinTTypeNameHV(TTypeNameHV* pOp1, TTypeNameHV* pOp2):TTypeNameHV(),m_pOp1(pOp1),m_pOp2(pOp2) + { + m_pOp1->incRef();m_pOp2->incRef(); + } + virtual ~BinTTypeNameHV() + { + m_pOp1->decRef(m_pOp1);m_pOp2->decRef(m_pOp2); + } + TTypeNameHV* op1() { return m_pOp1; } + TTypeNameHV* op2() { return m_pOp2; } + + unsigned int op1Eval(const unsigned int k){return this->op1()->eval(k);} + unsigned int op2Eval(const unsigned int k){return this->op2()->eval(k);} + const U& op1Val(const unsigned int k) {return this->op1()->val(k);} + const U& op2Val(const unsigned int k) {return this->op2()->val(k);} + void reset(){op1()->reset();op2()->reset();TTypeNameHV::reset();} +}; + +// Unary operator base class: + +template +class UnTTypeNameHV : public TTypeNameHV +{ + TTypeNameHV* m_pOp; +public: + UnTTypeNameHV(const U& val, TTypeNameHV* pOp):TTypeNameHV(val),m_pOp(pOp) + { + m_pOp->incRef(); + } + UnTTypeNameHV(TTypeNameHV* pOp):TTypeNameHV(),m_pOp(pOp) + { + m_pOp->incRef(); + } + virtual ~UnTTypeNameHV() + { + m_pOp->decRef(m_pOp); + } + TTypeNameHV* op() { return m_pOp; } + + unsigned int opEval(const unsigned int k){return this->op()->eval(k);} + const U& opVal(const unsigned int k) {return this->op()->val(k);} + void reset(){op()->reset();TTypeNameHV::reset();} +}; + +// ADDITION: + +template +struct TTypeNameADD : public BinTTypeNameHV +{ + TTypeNameADD(const U& val, TTypeNameHV* pOp1, TTypeNameHV* pOp2):BinTTypeNameHV(val,pOp1,pOp2){} + TTypeNameADD(TTypeNameHV* pOp1, TTypeNameHV* pOp2):BinTTypeNameHV(pOp1,pOp2){} + unsigned int eval(const unsigned int k) + { + unsigned int l=std::min(this->op1Eval(k),this->op2Eval(k)); + for(unsigned int i=this->length();ival(i)=this->op1Val(i)+this->op2Val(i); + return this->length()=l; + } +private: + void operator=(const TTypeNameADD&){} // not allowed +}; +template +struct TTypeNameADD1 : public UnTTypeNameHV +{ + const V m_a; + TTypeNameADD1(const U& val, const V& a, TTypeNameHV* pOp2):UnTTypeNameHV(val,pOp2),m_a(a){} + TTypeNameADD1(const V& a, TTypeNameHV* pOp2):UnTTypeNameHV(pOp2),m_a(a){} + unsigned int eval(const unsigned int k) + { + unsigned int l=this->opEval(k); + if (0==this->length()) { this->val(0)=m_a+this->opVal(0); this->length()=1; } + for(unsigned int i=this->length();ival(i)=this->opVal(i); + return this->length()=l; + } +private: + void operator=(const TTypeNameADD1&){} // not allowed +}; +template +struct TTypeNameADD2 : public UnTTypeNameHV +{ + const V m_b; + TTypeNameADD2(const U& val, TTypeNameHV* pOp1, const V& b):UnTTypeNameHV(val,pOp1),m_b(b){} + TTypeNameADD2(TTypeNameHV* pOp1, const V& b):UnTTypeNameHV(pOp1),m_b(b){} + unsigned int eval(const unsigned int k) + { + unsigned int l=this->opEval(k); + if (0==this->length()) { this->val(0)=this->opVal(0)+m_b; this->length()=1; } + for(unsigned int i=this->length();ival(i)=this->opVal(i); + return this->length()=l; + } +private: + void operator=(const TTypeNameADD2&){} // not allowed +}; +template +TTypeName operator+(const TTypeName& val1, const TTypeName& val2) +{ + TTypeNameHV* pHV=val1.length()>0 && val2.length()>0 ? + new TTypeNameADD(val1.val()+val2.val(),val1.getTTypeNameHV(),val2.getTTypeNameHV()): + new TTypeNameADD(val1.getTTypeNameHV(),val2.getTTypeNameHV()); + return TTypeName(pHV); +} +/* +template +TTypeName operator+(const typename Op::Underlying& a, const TTypeName& val2) +{ + TTypeNameHV* pHV=val2.length()>0 ? + new TTypeNameADD1::Underlying>(a+val2.val(), a, val2.getTTypeNameHV()): + new TTypeNameADD1::Underlying>(a, val2.getTTypeNameHV()); + return TTypeName(pHV); +} +template +TTypeName operator+(const TTypeName& val1, const typename Op::Underlying& b) +{ + TTypeNameHV* pHV=val1.length()>0? + new TTypeNameADD2::Underlying>(val1.val()+b, val1.getTTypeNameHV(), b): + new TTypeNameADD2::Underlying>(val1.getTTypeNameHV(), b); + return TTypeName(pHV); +} +template +TTypeName operator+(const typename Op::Base& a, const TTypeName& val2) +{ + TTypeNameHV* pHV=val2.length()>0 ? + new TTypeNameADD1::Base>(a+val2.val(), a, val2.getTTypeNameHV()): + new TTypeNameADD1::Base>(a, val2.getTTypeNameHV()); + return TTypeName(pHV); +} +template +TTypeName operator+(const TTypeName& val1, const typename Op::Base& b) +{ + TTypeNameHV* pHV=val1.length()>0? + new TTypeNameADD2::Base>(val1.val()+b, val1.getTTypeNameHV(), b): + new TTypeNameADD2::Base>(val1.getTTypeNameHV(), b); + return TTypeName(pHV); +} +*/ +template +TTypeName operator+(const V& a, const TTypeName& val2) +{ + TTypeNameHV* pHV=val2.length()>0 ? + new TTypeNameADD1(a+val2.val(), a, val2.getTTypeNameHV()): + new TTypeNameADD1(a, val2.getTTypeNameHV()); + return TTypeName(pHV); +} +template +TTypeName operator+(const TTypeName& val1, const V& b) +{ + TTypeNameHV* pHV=val1.length()>0? + new TTypeNameADD2(val1.val()+b, val1.getTTypeNameHV(), b): + new TTypeNameADD2(val1.getTTypeNameHV(), b); + return TTypeName(pHV); +} + +// SUBTRACTION: + +template +struct TTypeNameSUB : public BinTTypeNameHV +{ + TTypeNameSUB(const U& val, TTypeNameHV* pOp1, TTypeNameHV* pOp2):BinTTypeNameHV(val,pOp1,pOp2){} + TTypeNameSUB(TTypeNameHV* pOp1, TTypeNameHV* pOp2):BinTTypeNameHV(pOp1,pOp2){} + unsigned int eval(const unsigned int k) + { + unsigned int l=std::min(this->op1Eval(k),this->op2Eval(k)); + for(unsigned int i=this->length();ival(i)=this->op1Val(i)-this->op2Val(i); + return this->length()=l; + } +private: + void operator=(const TTypeNameSUB&){} // not allowed +}; +template +struct TTypeNameSUB1 : public UnTTypeNameHV +{ + const V m_a; + TTypeNameSUB1(const U& val, const V& a, TTypeNameHV* pOp2):UnTTypeNameHV(val,pOp2),m_a(a){} + TTypeNameSUB1(const V& a, TTypeNameHV* pOp2):UnTTypeNameHV(pOp2),m_a(a){} + unsigned int eval(const unsigned int k) + { + unsigned int l=this->opEval(k); + if (0==this->length()) { this->val(0)=m_a-this->opVal(0); this->length()=1; } + for(unsigned int i=this->length();ival(i)=Op::myNeg(this->opVal(i)); + return this->length()=l; + } +private: + void operator=(const TTypeNameSUB1&){} // not allowed +}; +template +struct TTypeNameSUB2 : public UnTTypeNameHV +{ + const V m_b; + TTypeNameSUB2(const U& val, TTypeNameHV* pOp1, const V& b):UnTTypeNameHV(val,pOp1),m_b(b){} + TTypeNameSUB2(TTypeNameHV* pOp1, const V& b):UnTTypeNameHV(pOp1),m_b(b){} + unsigned int eval(const unsigned int k) + { + unsigned int l=this->opEval(k); + if (0==this->length()) { this->val(0)=this->opVal(0)-m_b; this->length()=1; } + for(unsigned int i=this->length();ival(i)=this->opVal(i); + return this->length()=l; + } +private: + void operator=(const TTypeNameSUB2&){} // not allowed +}; +template +TTypeName operator-(const TTypeName& val1, const TTypeName& val2) +{ + TTypeNameHV* pHV=val1.length()>0 && val2.length()>0 ? + new TTypeNameSUB(val1.val()-val2.val(),val1.getTTypeNameHV(),val2.getTTypeNameHV()): + new TTypeNameSUB(val1.getTTypeNameHV(),val2.getTTypeNameHV()); + return TTypeName(pHV); +} +/* +template +TTypeName operator-(const typename Op::Underlying& a, const TTypeName& val2) +{ + TTypeNameHV* pHV=val2.length()>0 ? + new TTypeNameSUB1::Underlying>(a-val2.val(), a, val2.getTTypeNameHV()): + new TTypeNameSUB1::Underlying>(a, val2.getTTypeNameHV()); + return TTypeName(pHV); +} +template +TTypeName operator-(const TTypeName& val1, const typename Op::Underlying& b) +{ + TTypeNameHV* pHV=val1.length()>0 ? + new TTypeNameSUB2::Underlying>(val1.val()-b, val1.getTTypeNameHV(), b): + new TTypeNameSUB2::Underlying>(val1.getTTypeNameHV(), b); + return TTypeName(pHV); +} +template +TTypeName operator-(const typename Op::Base& a, const TTypeName& val2) +{ + TTypeNameHV* pHV=val2.length()>0 ? + new TTypeNameSUB1::Base>(a-val2.val(), a, val2.getTTypeNameHV()): + new TTypeNameSUB1::Base>(a, val2.getTTypeNameHV()); + return TTypeName(pHV); +} +template +TTypeName operator-(const TTypeName& val1, const typename Op::Base& b) +{ + TTypeNameHV* pHV=val1.length()>0 ? + new TTypeNameSUB2::Base>(val1.val()-b, val1.getTTypeNameHV(), b): + new TTypeNameSUB2::Base>(val1.getTTypeNameHV(), b); + return TTypeName(pHV); +} +*/ +template +TTypeName operator-(const V& a, const TTypeName& val2) +{ + TTypeNameHV* pHV=val2.length()>0 ? + new TTypeNameSUB1(a-val2.val(), a, val2.getTTypeNameHV()): + new TTypeNameSUB1(a, val2.getTTypeNameHV()); + return TTypeName(pHV); +} +template +TTypeName operator-(const TTypeName& val1, const V& b) +{ + TTypeNameHV* pHV=val1.length()>0 ? + new TTypeNameSUB2(val1.val()-b, val1.getTTypeNameHV(), b): + new TTypeNameSUB2(val1.getTTypeNameHV(), b); + return TTypeName(pHV); +} + +// MULTIPLICATION: + +template +struct TTypeNameMUL : public BinTTypeNameHV +{ + TTypeNameMUL(const U& val, TTypeNameHV* pOp1, TTypeNameHV* pOp2):BinTTypeNameHV(val,pOp1,pOp2){} + TTypeNameMUL(TTypeNameHV* pOp1, TTypeNameHV* pOp2):BinTTypeNameHV(pOp1,pOp2){} + unsigned int eval(const unsigned int k) + { + unsigned int l=std::min(this->op1Eval(k),this->op2Eval(k)); + for(unsigned int i=this->length();ival(i)=Op::myZero(); + for(unsigned int j=0;j<=i;++j) + Op::myCadd(this->val(i),this->op1Val(j)*this->op2Val(i-j)); + } + return this->length()=l; + } +private: + void operator=(const TTypeNameMUL&){} // not allowed +}; +template +struct TTypeNameMUL1 : public UnTTypeNameHV +{ + const V m_a; + TTypeNameMUL1(const U& val, const V& a, TTypeNameHV* pOp2):UnTTypeNameHV(val,pOp2),m_a(a){} + TTypeNameMUL1(const V& a, TTypeNameHV* pOp2):UnTTypeNameHV(pOp2),m_a(a){} + unsigned int eval(const unsigned int k) + { + unsigned int l=this->opEval(k); + for(unsigned int i=this->length();ival(i)=m_a*this->opVal(i); + return this->length()=l; + } +private: + void operator=(const TTypeNameMUL1&){} // not allowed +}; +template +struct TTypeNameMUL2 : public UnTTypeNameHV +{ + const V m_b; + TTypeNameMUL2(const U& val, TTypeNameHV* pOp1, const V& b):UnTTypeNameHV(val,pOp1),m_b(b){} + TTypeNameMUL2(TTypeNameHV* pOp1, const V& b):UnTTypeNameHV(pOp1),m_b(b){} + unsigned int eval(const unsigned int k) + { + unsigned int l=this->opEval(k); + for(unsigned int i=this->length();ival(i)=this->opVal(i)*m_b; + return this->length()=l; + } +private: + void operator=(const TTypeNameMUL2&){} // not allowed +}; +template +TTypeName operator*(const TTypeName& val1, const TTypeName& val2) +{ + TTypeNameHV* pHV=val1.length()>0 && val2.length()>0 ? + new TTypeNameMUL(val1.val()*val2.val(),val1.getTTypeNameHV(),val2.getTTypeNameHV()): + new TTypeNameMUL(val1.getTTypeNameHV(),val2.getTTypeNameHV()); + return TTypeName(pHV); +} +/* +template +TTypeName operator*(const typename Op::Underlying& a, const TTypeName& val2) +{ + TTypeNameHV* pHV=val2.length()>0 ? + new TTypeNameMUL1::Underlying>(a*val2.val(), a, val2.getTTypeNameHV()): + new TTypeNameMUL1::Underlying>(a, val2.getTTypeNameHV()); + return TTypeName(pHV); +} +template +TTypeName operator*(const TTypeName& val1, const typename Op::Underlying& b) +{ + TTypeNameHV* pHV=val1.length()>0 ? + new TTypeNameMUL2::Underlying>(val1.val()*b, val1.getTTypeNameHV(), b): + new TTypeNameMUL2::Underlying>(val1.getTTypeNameHV(), b); + return TTypeName(pHV); +} +template +TTypeName operator*(const typename Op::Base& a, const TTypeName& val2) +{ + TTypeNameHV* pHV=val2.length()>0 ? + new TTypeNameMUL1::Base>(a*val2.val(), a, val2.getTTypeNameHV()): + new TTypeNameMUL1::Base>(a, val2.getTTypeNameHV()); + return TTypeName(pHV); +} +template +TTypeName operator*(const TTypeName& val1, const typename Op::Base& b) +{ + TTypeNameHV* pHV=val1.length()>0 ? + new TTypeNameMUL2::Base>(val1.val()*b, val1.getTTypeNameHV(), b): + new TTypeNameMUL2::Base>(val1.getTTypeNameHV(), b); + return TTypeName(pHV); +} +*/ +template +TTypeName operator*(const V& a, const TTypeName& val2) +{ + TTypeNameHV* pHV=val2.length()>0 ? + new TTypeNameMUL1(a*val2.val(), a, val2.getTTypeNameHV()): + new TTypeNameMUL1(a, val2.getTTypeNameHV()); + return TTypeName(pHV); +} +template +TTypeName operator*(const TTypeName& val1, const V& b) +{ + TTypeNameHV* pHV=val1.length()>0 ? + new TTypeNameMUL2(val1.val()*b, val1.getTTypeNameHV(), b): + new TTypeNameMUL2(val1.getTTypeNameHV(), b); + return TTypeName(pHV); +} + +// DIVISION: + +template +struct TTypeNameDIV : public BinTTypeNameHV +{ + TTypeNameDIV(const U& val, TTypeNameHV* pOp1, TTypeNameHV* pOp2):BinTTypeNameHV(val,pOp1,pOp2){} + TTypeNameDIV(TTypeNameHV* pOp1, TTypeNameHV* pOp2):BinTTypeNameHV(pOp1,pOp2){} + unsigned int eval(const unsigned int k) + { + unsigned int l=std::min(this->op1Eval(k),this->op2Eval(k)); + for(unsigned int i=this->length();ival(i)=this->op1Val(i); + for(unsigned int j=1;j<=i;++j) Op::myCsub(this->val(i),this->op2Val(j)*this->val(i-j)); + Op::myCdiv(this->val(i), this->op2Val(0)); + } + return this->length()=l; + } +private: + void operator=(const TTypeNameDIV&){} // not allowed +}; +template +struct TTypeNameDIV1 : public UnTTypeNameHV +{ + const V m_a; + TTypeNameDIV1(const U& val, const V& a, TTypeNameHV* pOp2):UnTTypeNameHV(val,pOp2),m_a(a){} + TTypeNameDIV1(const V& a, TTypeNameHV* pOp2):UnTTypeNameHV(pOp2),m_a(a){} + unsigned int eval(const unsigned int k) + { + unsigned int l=this->opEval(k); + if (0==this->length()) { this->val(0)=m_a/this->opVal(0); this->length()=1; } + for(unsigned int i=this->length();ival(i)=Op::myZero(); + for(unsigned int j=1;j<=i;++j) Op::myCsub(this->val(i),this->opVal(j)*this->val(i-j)); + Op::myCdiv(this->val(i), this->opVal(0)); + } + return this->length()=l; + } +private: + void operator=(const TTypeNameDIV1&){} // not allowed +}; +template +struct TTypeNameDIV2 : public UnTTypeNameHV +{ + const V m_b; + TTypeNameDIV2(const U& val, TTypeNameHV* pOp1, const V& b):UnTTypeNameHV(val,pOp1),m_b(b){} + TTypeNameDIV2(TTypeNameHV* pOp1, const V& b):UnTTypeNameHV(pOp1),m_b(b){} + unsigned int eval(const unsigned int k) + { + unsigned int l=this->opEval(k); + for(unsigned int i=this->length();ival(i)=this->opVal(i)/m_b; + return this->length()=l; + } +private: + void operator=(const TTypeNameDIV2&){} // not allowed +}; +template +TTypeName operator/(const TTypeName& val1, const TTypeName& val2) +{ + TTypeNameHV* pHV=val1.length()>0 && val2.length()>0 ? + new TTypeNameDIV(val1.val()/val2.val(),val1.getTTypeNameHV(),val2.getTTypeNameHV()): + new TTypeNameDIV(val1.getTTypeNameHV(),val2.getTTypeNameHV()); + return TTypeName(pHV); +} +/* +template +TTypeName operator/(const typename Op::Underlying& a, const TTypeName& val2) +{ + TTypeNameHV* pHV=val2.length()>0 ? + new TTypeNameDIV1::Underlying>(a/val2.val(), a, val2.getTTypeNameHV()): + new TTypeNameDIV1::Underlying>(a, val2.getTTypeNameHV()); + return TTypeName(pHV); +} +template +TTypeName operator/(const TTypeName& val1, const typename Op::Underlying& b) +{ + TTypeNameHV* pHV=val1.length()>0 ? + new TTypeNameDIV2::Underlying>(val1.val()/b, val1.getTTypeNameHV(), b): + new TTypeNameDIV2::Underlying>(val1.getTTypeNameHV(), b); + return TTypeName(pHV); +} +template +TTypeName operator/(const typename Op::Base& a, const TTypeName& val2) +{ + TTypeNameHV* pHV=val2.length()>0 ? + new TTypeNameDIV1::Base>(a/val2.val(), a, val2.getTTypeNameHV()): + new TTypeNameDIV1::Base>(a, val2.getTTypeNameHV()); + return TTypeName(pHV); +} +template +TTypeName operator/(const TTypeName& val1, const typename Op::Base& b) +{ + TTypeNameHV* pHV=val1.length()>0 ? + new TTypeNameDIV2::Base>(val1.val()/b, val1.getTTypeNameHV(), b): + new TTypeNameDIV2::Base>(val1.getTTypeNameHV(), b); + return TTypeName(pHV); +} +*/ +template +TTypeName operator/(const V& a, const TTypeName& val2) +{ + TTypeNameHV* pHV=val2.length()>0 ? + new TTypeNameDIV1(a/val2.val(), a, val2.getTTypeNameHV()): + new TTypeNameDIV1(a, val2.getTTypeNameHV()); + return TTypeName(pHV); +} +template +TTypeName operator/(const TTypeName& val1, const V& b) +{ + TTypeNameHV* pHV=val1.length()>0 ? + new TTypeNameDIV2(val1.val()/b, val1.getTTypeNameHV(), b): + new TTypeNameDIV2(val1.getTTypeNameHV(), b); + return TTypeName(pHV); +} + +// COMPOUND ASSIGNMENTS: + +template TTypeName& TTypeName::operator+=(const TTypeName& val) { return (*this)=(*this)+val; } +template TTypeName& TTypeName::operator-=(const TTypeName& val) { return (*this)=(*this)-val; } +template TTypeName& TTypeName::operator*=(const TTypeName& val) { return (*this)=(*this)*val; } +template TTypeName& TTypeName::operator/=(const TTypeName& val) { return (*this)=(*this)/val; } +template template TTypeName& TTypeName::operator+=(const V& val) { return (*this)=(*this)+val; } +template template TTypeName& TTypeName::operator-=(const V& val) { return (*this)=(*this)-val; } +template template TTypeName& TTypeName::operator*=(const V& val) { return (*this)=(*this)*val; } +template template TTypeName& TTypeName::operator/=(const V& val) { return (*this)=(*this)/val; } + +// UNARY MINUS + +template +struct TTypeNameUMINUS : public UnTTypeNameHV +{ + TTypeNameUMINUS(const U& val, TTypeNameHV* pOp):UnTTypeNameHV(val,pOp){} + TTypeNameUMINUS(TTypeNameHV* pOp):UnTTypeNameHV(pOp){} + unsigned int eval(const unsigned int k) + { + unsigned int l=this->opEval(k); + for(unsigned int i=this->length();ival(i)=Op::myNeg(this->opVal(i)); + return this->length()=l; + } +private: + void operator=(const TTypeNameUMINUS&){} // not allowed +}; + +template +TTypeName operator-(const TTypeName& val) +{ + TTypeNameHV* pHV=val.length()>0 ? + new TTypeNameUMINUS(Op::myNeg(val.val()),val.getTTypeNameHV()) : + new TTypeNameUMINUS(val.getTTypeNameHV()); + return TTypeName(pHV); +} + +// UNARY PLUS + +template +struct TTypeNameUPLUS : public UnTTypeNameHV +{ + TTypeNameUPLUS(const U& val, TTypeNameHV* pOp):UnTTypeNameHV(val,pOp){} + TTypeNameUPLUS(TTypeNameHV* pOp):UnTTypeNameHV(pOp){} + unsigned int eval(const unsigned int k) + { + unsigned int l=this->opEval(k); + for(unsigned int i=this->length();ival(i)=+this->opVal(i); + return this->length()=l; + } +private: + void operator=(const TTypeNameUPLUS&){} // not allowed +}; + +template +TTypeName operator+(const TTypeName& val) +{ + TTypeNameHV* pHV=val.length()>0 ? + new TTypeNameUPLUS(+val.val(),val.getTTypeNameHV()) : + new TTypeNameUPLUS(val.getTTypeNameHV()); + return TTypeName(pHV); +} + +// POWER + +template +struct TTypeNamePOW : public UnTTypeNameHV +{ + TTypeNamePOW(const U& val, TTypeNameHV* pOp):UnTTypeNameHV(val,pOp){} + TTypeNamePOW(TTypeNameHV* pOp):UnTTypeNameHV(pOp){} + unsigned int eval(const unsigned int k) + { + unsigned int l=this->opEval(k); + for(unsigned int i=this->length();ival(i)=this->opVal(i); + return this->length()=l; + } +private: + void operator=(const TTypeNamePOW&){} // not allowed +}; +template +struct TTypeNamePOW1 : public UnTTypeNameHV +{ + TTypeNamePOW1(const U& val, TTypeNameHV* pOp2):UnTTypeNameHV(val,pOp2){} + TTypeNamePOW1(TTypeNameHV* pOp2):UnTTypeNameHV(pOp2){} + unsigned int eval(const unsigned int k) + { + unsigned int l=this->opEval(k); + for(unsigned int i=this->length();ival(i)=this->opVal(i); + return this->length()=l; + } +private: + void operator=(const TTypeNamePOW1&){} // not allowed +}; +template +struct TTypeNamePOW2 : public UnTTypeNameHV +{ + TTypeNamePOW2(const U& val, TTypeNameHV* pOp1):UnTTypeNameHV(val,pOp1){} + TTypeNamePOW2(TTypeNameHV* pOp1):UnTTypeNameHV(pOp1){} + unsigned int eval(const unsigned int k) + { + unsigned int l=this->opEval(k); + for(unsigned int i=this->length();ival(i)=this->opVal(i); + return this->length()=l; + } +private: + void operator=(const TTypeNamePOW2&){} // not allowed +}; +template +TTypeName pow(const TTypeName& val1, const TTypeName& val2) +{ + TTypeName tmp(exp(val2*log(val1))); + TTypeNameHV* pHV=val1.length()>0 && val2.length()>0 ? + new TTypeNamePOW(Op::myPow(val1.val(),val2.val()),tmp.getTTypeNameHV()) : + new TTypeNamePOW(tmp.getTTypeNameHV()); + return TTypeName(pHV); +} +/* +template +TTypeName pow(const typename Op::Underlying& a, const TTypeName& val2) +{ + TTypeName tmp(exp(val2*Op::myLog(a))); + TTypeNameHV* pHV=val2.length()>0 ? + new TTypeNamePOW1::Underlying>(Op::myPow(a,val2.val()), tmp.getTTypeNameHV()) : + new TTypeNamePOW1::Underlying>(tmp.getTTypeNameHV()); + return TTypeName(pHV); +} +template +TTypeName pow(const TTypeName& val1, const typename Op::Underlying& b) +{ + TTypeName tmp(exp(b*log(val1))); + TTypeNameHV* pHV=val1.length()>0 ? + new TTypeNamePOW2::Underlying>(Op::myPow(val1.val(),b), tmp.getTTypeNameHV()) : + new TTypeNamePOW2::Underlying>(tmp.getTTypeNameHV()); + return TTypeName(pHV); +} +template +TTypeName pow(const typename Op::Base& a, const TTypeName& val2) +{ + TTypeName tmp(exp(val2*Op::Base>::myLog(a))); + TTypeNameHV* pHV=val2.length()>0 ? + new TTypeNamePOW1::Base>(Op::myPow(a,val2.val()), tmp.getTTypeNameHV()) : + new TTypeNamePOW1::Base>(tmp.getTTypeNameHV()); + return TTypeName(pHV); +} +template +TTypeName pow(const TTypeName& val1, const typename Op::Base& b) +{ + TTypeName tmp(exp(b*log(val1))); + TTypeNameHV* pHV=val1.length()>0 ? + new TTypeNamePOW2::Base>(Op::myPow(val1.val(),b), tmp.getTTypeNameHV()) : + new TTypeNamePOW2::Base>(tmp.getTTypeNameHV()); + return TTypeName(pHV); +} +*/ +template +TTypeName pow(const V& a, const TTypeName& val2) +{ + TTypeName tmp(exp(val2*Op::myLog(a))); + TTypeNameHV* pHV=val2.length()>0 ? + new TTypeNamePOW1(Op::myPow(a,val2.val()), tmp.getTTypeNameHV()) : + new TTypeNamePOW1(tmp.getTTypeNameHV()); + return TTypeName(pHV); +} +template +TTypeName pow(const TTypeName& val1, const V& b) +{ + TTypeName tmp(exp(b*log(val1))); + TTypeNameHV* pHV=val1.length()>0 ? + new TTypeNamePOW2(Op::myPow(val1.val(),b), tmp.getTTypeNameHV()) : + new TTypeNamePOW2(tmp.getTTypeNameHV()); + return TTypeName(pHV); +} + +// SQR + +template +struct TTypeNameSQR : public UnTTypeNameHV +{ + TTypeNameSQR(const U& val, TTypeNameHV* pOp):UnTTypeNameHV(val,pOp){} + TTypeNameSQR(TTypeNameHV* pOp):UnTTypeNameHV(pOp){} + unsigned int eval(const unsigned int k) + { + unsigned int l=this->opEval(k); + if (0==this->length()) { this->val(0)=Op::mySqr(this->opVal(0)); this->length()=1; } + for(unsigned int i=this->length();ival(i)=Op::myZero(); + unsigned int m=(i+1)/2; + for(unsigned int j=0;j::myCadd(this->val(i), this->opVal(i-j)*this->opVal(j)); + Op::myCmul(this->val(i), Op::myTwo()); + if (0==i%2) Op::myCadd(this->val(i), Op::mySqr(this->opVal(m))); + } + return this->length()=l; + } +private: + void operator=(const TTypeNameSQR&){} // not allowed +}; +template +TTypeName sqr(const TTypeName& val) +{ + TTypeNameHV* pHV=val.length()>0 ? + new TTypeNameSQR(Op::mySqr(val.val()), val.getTTypeNameHV()) : + new TTypeNameSQR(val.getTTypeNameHV()); + return TTypeName(pHV); +} + +// SQRT + +template +struct TTypeNameSQRT : public UnTTypeNameHV +{ + TTypeNameSQRT(const U& val, TTypeNameHV* pOp):UnTTypeNameHV(val,pOp){} + TTypeNameSQRT(TTypeNameHV* pOp):UnTTypeNameHV(pOp){} + unsigned int eval(const unsigned int k) + { + unsigned int l=this->opEval(k); + if (0==this->length()) { this->val(0)=Op::mySqrt(this->opVal(0)); this->length()=1; } + for(unsigned int i=this->length();ival(i)=Op::myZero(); + unsigned int m=(i+1)/2; + for(unsigned int j=1;j::myCadd(this->val(i), this->val(i-j)*this->val(j)); + Op::myCmul(this->val(i), Op::myTwo()); + if (0==i%2) Op::myCadd(this->val(i), Op::mySqr(this->val(m))); + this->val(i)=(this->opVal(i)-this->val(i))/(Op::myTwo()*this->val(0)); + } + return this->length()=l; + } +private: + void operator=(const TTypeNameSQRT&){} // not allowed +}; +template +TTypeName sqrt(const TTypeName& val) +{ + TTypeNameHV* pHV=val.length()>0 ? + new TTypeNameSQRT(Op::mySqrt(val.val()), val.getTTypeNameHV()): + new TTypeNameSQRT(val.getTTypeNameHV()); + return TTypeName(pHV); +} + +// EXP + +template +struct TTypeNameEXP : public UnTTypeNameHV +{ + TTypeNameEXP(const U& val, TTypeNameHV* pOp):UnTTypeNameHV(val,pOp){} + TTypeNameEXP(TTypeNameHV* pOp):UnTTypeNameHV(pOp){} + unsigned int eval(const unsigned int k) + { + unsigned int l=this->opEval(k); + if (0==this->length()) { this->val(0)=Op::myExp(this->opVal(0)); this->length()=1; } + for(unsigned int i=this->length();ival(i)=Op::myZero(); + for(unsigned int j=0;j::myCadd(this->val(i), (Op::myOne()-Op::myInteger(j) / + Op::myInteger(i))*this->opVal(i-j)*this->val(j)); + } + return this->length()=l; + } +private: + void operator=(const TTypeNameEXP&){} // not allowed +}; +template +TTypeName exp(const TTypeName& val) +{ + TTypeNameHV* pHV=val.length()>0 ? + new TTypeNameEXP(Op::myExp(val.val()), val.getTTypeNameHV()): + new TTypeNameEXP(val.getTTypeNameHV()); + return TTypeName(pHV); +} + +// LOG + +template +struct TTypeNameLOG : public UnTTypeNameHV +{ + TTypeNameLOG(const U& val, TTypeNameHV* pOp):UnTTypeNameHV(val,pOp){} + TTypeNameLOG(TTypeNameHV* pOp):UnTTypeNameHV(pOp){} + unsigned int eval(const unsigned int k) + { + unsigned int l=this->opEval(k); + if (0==this->length()) { this->val(0)=Op::myLog(this->opVal(0)); this->length()=1; } + for(unsigned int i=this->length();ival(i)=this->opVal(i); + for(unsigned int j=1;j::myCsub(this->val(i), (Op::myOne()-Op::myInteger(j) / + Op::myInteger(i))*this->opVal(j)*this->val(i-j)); + Op::myCdiv(this->val(i), this->opVal(0)); + } + return this->length()=l; + } +private: + void operator=(const TTypeNameLOG&){} // not allowed +}; +template +TTypeName log(const TTypeName& val) +{ + TTypeNameHV* pHV=val.length()>0 ? + new TTypeNameLOG(Op::myLog(val.val()), val.getTTypeNameHV()): + new TTypeNameLOG(val.getTTypeNameHV()); + return TTypeName(pHV); +} + +// SIN + +template +struct TTypeNameSIN : public UnTTypeNameHV +{ + U m_COS[N]; + TTypeNameSIN(const U& val, TTypeNameHV* pOp):UnTTypeNameHV(val,pOp){m_COS[0]=Op::myCos(this->opVal(0));} + TTypeNameSIN(TTypeNameHV* pOp):UnTTypeNameHV(pOp){} + unsigned int eval(const unsigned int k) + { + unsigned int l=this->opEval(k); + if (0==this->length()) + { + this->val(0)=Op::mySin(this->opVal(0)); + m_COS[0]=Op::myCos(this->opVal(0)); + this->length()=1; + } + for(unsigned int i=this->length();ival(i)=Op::myZero(); + for(unsigned int j=0;j::myCadd(this->val(i), Op::myInteger(j+1)*m_COS[i-1-j]*this->opVal(j+1)); + Op::myCdiv(this->val(i), Op::myInteger(i)); + m_COS[i]=Op::myZero(); + for(unsigned int j=0;j::myCsub(m_COS[i], Op::myInteger(j+1)*this->val(i-1-j)*this->opVal(j+1)); + Op::myCdiv(m_COS[i], Op::myInteger(i)); + } + return this->length()=l; + } +private: + void operator=(const TTypeNameSIN&){} // not allowed +}; +template +TTypeName sin(const TTypeName& val) +{ + TTypeNameHV* pHV=val.length()>0 ? + new TTypeNameSIN(Op::mySin(val.val()), val.getTTypeNameHV()): + new TTypeNameSIN(val.getTTypeNameHV()); + return TTypeName(pHV); +} + +// COS + +template +struct TTypeNameCOS : public UnTTypeNameHV +{ + U m_SIN[N]; + TTypeNameCOS(const U& val, TTypeNameHV* pOp):UnTTypeNameHV(val,pOp){m_SIN[0]=Op::mySin(this->opVal(0));} + TTypeNameCOS(TTypeNameHV* pOp):UnTTypeNameHV(pOp){} + unsigned int eval(const unsigned int k) + { + unsigned int l=this->opEval(k); + if (0==this->length()) + { + this->val(0)=Op::myCos(this->opVal(0)); + m_SIN[0]=Op::mySin(this->opVal(0)); + this->length()=1; + } + for(unsigned int i=this->length();ival(i)=Op::myZero(); + for(unsigned int j=0;j::myCsub(this->val(i), Op::myInteger(j+1)*m_SIN[i-1-j]*this->opVal(j+1)); + Op::myCdiv(this->val(i), Op::myInteger(i)); + m_SIN[i]=Op::myZero(); + for(unsigned int j=0;j::myCadd(m_SIN[i], Op::myInteger(j+1)*this->val(i-1-j)*this->opVal(j+1)); + Op::myCdiv(m_SIN[i], Op::myInteger(i)); + } + return this->length()=l; + } +private: + void operator=(const TTypeNameCOS&){} // not allowed +}; +template +TTypeName cos(const TTypeName& val) +{ + TTypeNameHV* pHV=val.length()>0 ? + new TTypeNameCOS(Op::myCos(val.val()), val.getTTypeNameHV()): + new TTypeNameCOS(val.getTTypeNameHV()); + return TTypeName(pHV); +} + +// TAN + +template +struct TTypeNameTAN : public BinTTypeNameHV +{ + TTypeNameTAN(const U& val, TTypeNameHV* pOp, TTypeNameHV* pSqrCos):BinTTypeNameHV(val,pOp,pSqrCos){} + TTypeNameTAN(TTypeNameHV* pOp, TTypeNameHV* pSqrCos):BinTTypeNameHV(pOp,pSqrCos){} + unsigned int eval(const unsigned int k) + { + unsigned int l=std::min(this->op1Eval(k),this->op2Eval(k)); + if (0==this->length()) { this->val(0)=Op::myTan(this->op1Val(0)); this->length()=1; } + for(unsigned int i=this->length();ival(i)=Op::myZero(); + for(unsigned int j=1;j::myCadd(this->val(i), Op::myInteger(j)*this->val(j)*this->op2Val(i-j)); + this->val(i)=(this->op1Val(i)-this->val(i)/Op::myInteger(i))/this->op2Val(0); + } + return this->length()=l; + } +private: + void operator=(const TTypeNameTAN&){} // not allowed +}; +template +TTypeName tan(const TTypeName& val) +{ + TTypeName tmp(sqr(cos(val))); + TTypeNameHV* pHV=val.length()>0 ? + new TTypeNameTAN(Op::myTan(val.val()), val.getTTypeNameHV(), tmp.getTTypeNameHV()): + new TTypeNameTAN(val.getTTypeNameHV(), tmp.getTTypeNameHV()); + return TTypeName(pHV); +} + +// ASIN + +template +struct TTypeNameASIN : public BinTTypeNameHV +{ + TTypeNameASIN(const U& val, TTypeNameHV* pOp, TTypeNameHV* pSqrt):BinTTypeNameHV(val,pOp,pSqrt){} + TTypeNameASIN(TTypeNameHV* pOp, TTypeNameHV* pSqrt):BinTTypeNameHV(pOp,pSqrt){} + unsigned int eval(const unsigned int k) + { + unsigned int l=std::min(this->op1Eval(k),this->op2Eval(k)); + if (0==this->length()) { this->val(0)=Op::myAsin(this->op1Val(0)); this->length()=1; } + for(unsigned int i=this->length();ival(i)=Op::myZero(); + for(unsigned int j=1;j::myCadd(this->val(i), Op::myInteger(j)*this->val(j)*this->op2Val(i-j)); + this->val(i)=(this->op1Val(i)-this->val(i)/Op::myInteger(i))/this->op2Val(0); + } + return this->length()=l; + } +private: + void operator=(const TTypeNameASIN&){} // not allowed +}; +template +TTypeName asin(const TTypeName& val) +{ + TTypeName tmp(sqrt(Op::myOne()-sqr(val))); + TTypeNameHV* pHV=val.length()>0 ? + new TTypeNameASIN(Op::myAsin(val.val()), val.getTTypeNameHV(), tmp.getTTypeNameHV()): + new TTypeNameASIN(val.getTTypeNameHV(), tmp.getTTypeNameHV()); + return TTypeName(pHV); +} + +// ACOS + +template +struct TTypeNameACOS : public BinTTypeNameHV +{ + TTypeNameACOS(const U& val, TTypeNameHV* pOp, TTypeNameHV* pSqrt):BinTTypeNameHV(val,pOp,pSqrt){} + TTypeNameACOS(TTypeNameHV* pOp, TTypeNameHV* pSqrt):BinTTypeNameHV(pOp,pSqrt){} + unsigned int eval(const unsigned int k) + { + unsigned int l=std::min(this->op1Eval(k),this->op2Eval(k)); + if (0==this->length()) { this->val(0)=Op::myAcos(this->op1Val(0)); this->length()=1; } + for(unsigned int i=this->length();ival(i)=Op::myZero(); + for(unsigned int j=1;j::myCadd(this->val(i), Op::myInteger(j)*this->val(j)*this->op2Val(i-j)); + this->val(i)=Op::myNeg((this->op1Val(i)+this->val(i)/Op::myInteger(i))/this->op2Val(0)); + } + return this->length()=l; + } +private: + void operator=(const TTypeNameACOS&){} // not allowed +}; +template +TTypeName acos(const TTypeName& val) +{ + TTypeName tmp(sqrt(Op::myOne()-sqr(val))); + TTypeNameHV* pHV=val.length()>0 ? + new TTypeNameACOS(Op::myAcos(val.val()), val.getTTypeNameHV(), tmp.getTTypeNameHV()): + new TTypeNameACOS(val.getTTypeNameHV(), tmp.getTTypeNameHV()); + return TTypeName(pHV); +} + +// ATAN + +template +struct TTypeNameATAN : public BinTTypeNameHV +{ + TTypeNameATAN(const U& val, TTypeNameHV* pOp, TTypeNameHV* p1pSqr):BinTTypeNameHV(val,pOp,p1pSqr){} + TTypeNameATAN(TTypeNameHV* pOp, TTypeNameHV* p1pSqr):BinTTypeNameHV(pOp,p1pSqr){} + unsigned int eval(const unsigned int k) + { + unsigned int l=std::min(this->op1Eval(k),this->op2Eval(k)); + if (0==this->length()) { this->val(0)=Op::myAtan(this->op1Val(0)); this->length()=1; } + for(unsigned int i=this->length();ival(i)=Op::myZero(); + for(unsigned int j=1;j::myCadd(this->val(i), Op::myInteger(j)*this->val(j)*this->op2Val(i-j)); + this->val(i)=(this->op1Val(i)-this->val(i)/Op::myInteger(i))/this->op2Val(0); + } + return this->length()=l; + } +private: + void operator=(const TTypeNameATAN&){} // not allowed +}; +template +TTypeName atan(const TTypeName& val) +{ + TTypeName tmp(Op::myOne()+sqr(val)); + TTypeNameHV* pHV=val.length()>0 ? + new TTypeNameATAN(Op::myAtan(val.val()), val.getTTypeNameHV(), tmp.getTTypeNameHV()): + new TTypeNameATAN(val.getTTypeNameHV(), tmp.getTTypeNameHV()); + return TTypeName(pHV); +} + +// Ned's diff operator + +template +struct DIFF : public UnTTypeNameHV +{ + int m_b; + DIFF(const U& val, TTypeNameHV* pOp, const int b):UnTTypeNameHV(val,pOp),m_b(b){} + DIFF(TTypeNameHV* pOp, const int b):UnTTypeNameHV(pOp),m_b(b){} + unsigned int eval(const unsigned int k) + { +// IN ORDER TO COMPUTE i'th ORDER COEFFICIENTS OF diff(m_o1,b) +// WE NEED (i+b)'th ORDER COEFFICIENTS OF TaylorOp. + unsigned int l=this->opEval(k+m_b); +// NOW WE SHOULD HAVE op1 EXPANDED TO DEGREE (i+b)'th ORDER, +// WE CAN PROCEED TO COMPUTE UP TO i'TH ORDER COEFFICIENTS +// OF diff(op1,b). + if (this->length()+m_blength();ii;--j){ fact*=j; } + this->val(i)=this->opVal(i+m_b)*fact; + } + this->length()=l-m_b; + } + return this->length(); + } +private: + void operator=(const DIFF&){} // not allowed +}; +template +TTypeName diff(const TTypeName& val, const int b) +{ +// IF THE ARGUMENT TAYLOR EXPANSION TaylorOp HAS BEEN EVALUATED TO +// DEGREE i THEN WE CAN EVALUATE THE ZERO ORDER VALUE OF +// diff(TaylorOp,i). +// THIS FUNCTION EVALUATES THE 0.ORDER COEFFICIENT + TTypeNameHV* pHV=0; + if (val.length()>b) + { + unsigned int fact=1; + for(unsigned int j=b;j>1;--j){ fact*=j; } + pHV=new DIFF(val[b]*fact, val.getTTypeNameHV(), b); + } + else + { + pHV=new DIFF(val.getTTypeNameHV(), b); + } + return TTypeName(pHV); +} + + +template struct Op< TTypeName > +{ + typedef TTypeName V; + typedef TTypeName Underlying; + typedef typename Op::Base Base; + static Base myInteger(const int i) { return Base(i); } + static Base myZero() { return myInteger(0); } + static Base myOne() { return myInteger(1);} + static Base myTwo() { return myInteger(2); } + static Base myPI() { return Op::myPI(); } + static V myPos(const V& x) { return +x; } + static V myNeg(const V& x) { return -x; } + template static V& myCadd(V& x, const Y& y) { return x+=y; } + template static V& myCsub(V& x, const Y& y) { return x-=y; } + template static V& myCmul(V& x, const Y& y) { return x*=y; } + template static V& myCdiv(V& x, const Y& y) { return x/=y; } + static V myInv(const V& x) { return myOne()/x; } + static V mySqr(const V& x) { return fadbad::sqr(x); } + template + static V myPow(const X& x, const Y& y) { return fadbad::pow(x,y); } + static V mySqrt(const V& x) { return fadbad::sqrt(x); } + static V myLog(const V& x) { return fadbad::log(x); } + static V myExp(const V& x) { return fadbad::exp(x); } + static V mySin(const V& x) { return fadbad::sin(x); } + static V myCos(const V& x) { return fadbad::cos(x); } + static V myTan(const V& x) { return fadbad::tan(x); } + static V myAsin(const V& x) { return fadbad::asin(x); } + static V myAcos(const V& x) { return fadbad::acos(x); } + static V myAtan(const V& x) { return fadbad::atan(x); } + static bool myEq(const V& x, const V& y) { return x==y; } + static bool myNe(const V& x, const V& y) { return x!=y; } + static bool myLt(const V& x, const V& y) { return xy; } + static bool myGe(const V& x, const V& y) { return x>=y; } +}; + +} // namespace fadbad + +#endif diff --git a/python/src/clexulator.cpp b/python/src/clexulator.cpp index 3fbeb4b..e4040ba 100644 --- a/python/src/clexulator.cpp +++ b/python/src/clexulator.cpp @@ -1161,7 +1161,41 @@ PYBIND11_MODULE(_clexulator, m) { Change in extensive correlations, relative to `reference_extensive_correlations`. )pbdoc", py::arg("key"), py::arg("new_value"), - py::arg("`reference_extensive_correlations")); + py::arg("reference_extensive_correlations")) + .def("grad_correlations", &clexulator::Correlations::grad_correlations, + py::return_value_policy::reference_internal, R"pbdoc( + Calculates and returns gradients of correlations with respect to + degrees of freedom corresponding to `dof_key` + + Returns a :math:`\\mathbf{M} \times \\mathbf{N}` matrix, + where :math:`\\mathbf{M}` represents the number of correlations, + and :math:`\\mathbf{N} represents the dimensions of the degrees + of freedom corresponding to `dof_key` + + + For example, if `dof_key` is `Hstrain`, N will be 6. + if `dof_key` is `disp`, and there are 3 sites in the configuration, + N will be 9 (:math:`x, y, z` displacements for each site) + + Each row of the matrix corresponds to gradients of all :math:`\\mathbf{M}` + correlations with respect to one of the dimensions of degrees of freedom + corresponding to `dof_key` + + Parameters + ---------- + dof_key: str + Specifies the type of DoF + Can be `occ`, `disp`, `Hstrain`, etc. + DoF that is set in the `Prim` + + Returns + ------- + gradients_of_correlations: np.ndarray[np.float64[m,n]] + gradients_of_correlations + where m represents number of correlations and + n represents number of degrees of freedom + )pbdoc", + py::arg("dof_key")); py::class_>( diff --git a/python/tests/clexulator/data/DiffClexulatorSquareLatticeTest/basis_sets/bset.default/DiffClexulatorSquareLatticeTest_Clexulator_default.cc b/python/tests/clexulator/data/DiffClexulatorSquareLatticeTest/basis_sets/bset.default/DiffClexulatorSquareLatticeTest_Clexulator_default.cc new file mode 100644 index 0000000..99a16a0 --- /dev/null +++ b/python/tests/clexulator/data/DiffClexulatorSquareLatticeTest/basis_sets/bset.default/DiffClexulatorSquareLatticeTest_Clexulator_default.cc @@ -0,0 +1,633 @@ +#include +#include "casm/clexulator/BaseClexulator.hh" +#include "casm/global/eigen.hh" +#include "casm/clexulator/DiffClexParamPack.hh" + + + + +/****** PROJECT SPECIFICATIONS ****** + + ****** prim.json ****** + +{ + "basis" : [ + { + "coordinate" : [ 0.000000000000, 0.000000000000, 0.000000000000 ], + "dofs" : { + "disp" : { + "axis_names" : [ "dx", "dy", "dz" ], + "basis" : [ + [ 1.000000000000, 0.000000000000, 0.000000000000 ], + [ 0.000000000000, 1.000000000000, 0.000000000000 ], + [ 0.000000000000, 0.000000000000, 1.000000000000 ] + ] + } + }, + "occupants" : [ "A" ] + } + ], + "coordinate_mode" : "Fractional", + "lattice_vectors" : [ + [ 3.000000000000, 0.000000000000, 0.000000000000 ], + [ 0.000000000000, 3.000000000000, 0.000000000000 ], + [ 0.000000000000, 0.000000000000, 10.000000000000 ] + ], + "title" : "DiffClexulatorSquareLatticeTest" +} + + ****** bspecs.json ****** + +{ + "basis_function_specs" : { + "dofs" : [ "disp" ], + "global_max_poly_order" : 5, + "param_pack_type" : "DIFF" + }, + "cluster_specs" : { + "method" : "periodic_max_length", + "params" : { + "generating_group" : [ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15 ], + "orbit_branch_specs" : { + "0" : { + "max_length" : 0.000000000000 + }, + "1" : { + "max_length" : 3.010000000000 + } + } + } + } +} + +**/ + + +/// \brief Returns a clexulator::BaseClexulator* owning a DiffClexulatorSquareLatticeTest_Clexulator_default +extern "C" CASM::clexulator::BaseClexulator *make_DiffClexulatorSquareLatticeTest_Clexulator_default(); + +namespace CASM { +namespace clexulator { + + /****** GENERATED CLEXPARAMPACK DEFINITION ******/ + + + typedef DiffClexParamPack ParamPack; + + + /****** GENERATED CLEXULATOR DEFINITION ******/ + + class DiffClexulatorSquareLatticeTest_Clexulator_default : public clexulator::BaseClexulator { + + public: + + DiffClexulatorSquareLatticeTest_Clexulator_default(); + + ~DiffClexulatorSquareLatticeTest_Clexulator_default(); + + ClexParamPack const ¶m_pack() const override { + return m_params; + } + + ClexParamPack ¶m_pack() override { + return m_params; + } + + + template + Scalar eval_bfunc_0_0() const; + + template + Scalar eval_bfunc_1_0() const; + template + Scalar eval_bfunc_1_1() const; + template + Scalar eval_bfunc_1_2() const; + template + Scalar eval_bfunc_1_3() const; + template + Scalar eval_bfunc_1_4() const; + template + Scalar eval_bfunc_1_5() const; + + template + Scalar site_eval_bfunc_1_0_at_0() const; + template + Scalar site_eval_bfunc_1_1_at_0() const; + template + Scalar site_eval_bfunc_1_2_at_0() const; + template + Scalar site_eval_bfunc_1_3_at_0() const; + template + Scalar site_eval_bfunc_1_4_at_0() const; + template + Scalar site_eval_bfunc_1_5_at_0() const; + + + private: + + // ParamPack object, which stores temporary data for calculations + mutable ParamPack m_params; + + // typedef for method pointers of scalar type double + typedef double (DiffClexulatorSquareLatticeTest_Clexulator_default::*BasisFuncPtr_0)() const; + + // typedef for method pointers + typedef double (DiffClexulatorSquareLatticeTest_Clexulator_default::*DeltaBasisFuncPtr_0)(int, int) const; + + // array of pointers to member functions for calculating basis functions of scalar type double + BasisFuncPtr_0 m_orbit_func_table_0[7]; + + // array of pointers to member functions for calculating flower functions of scalar type double + BasisFuncPtr_0 m_flower_func_table_0[1][7]; + + // array of pointers to member functions for calculating DELTA flower functions of scalar type double + DeltaBasisFuncPtr_0 m_delta_func_table_0[1][7]; + + // typedef for method pointers of scalar type ParamPack::DiffScalar + typedef ParamPack::DiffScalar (DiffClexulatorSquareLatticeTest_Clexulator_default::*BasisFuncPtr_1)() const; + + // typedef for method pointers + typedef ParamPack::DiffScalar (DiffClexulatorSquareLatticeTest_Clexulator_default::*DeltaBasisFuncPtr_1)(int, int) const; + + // array of pointers to member functions for calculating basis functions of scalar type ParamPack::DiffScalar + BasisFuncPtr_1 m_orbit_func_table_1[7]; + + // array of pointers to member functions for calculating flower functions of scalar type ParamPack::DiffScalar + BasisFuncPtr_1 m_flower_func_table_1[1][7]; + + // array of pointers to member functions for calculating DELTA flower functions of scalar type ParamPack::DiffScalar + DeltaBasisFuncPtr_1 m_delta_func_table_1[1][7]; + + //ClexParamPack allocation for evaluated correlations + ParamPack::Key m_corr_param_key; + //ClexParamPack allocation for DoF disp + ParamPack::Key m_disp_var_param_key; + + /// \brief Clone the DiffClexulatorSquareLatticeTest_Clexulator_default + BaseClexulator *_clone() const override { + return new DiffClexulatorSquareLatticeTest_Clexulator_default(*this); + } + + /// \brief Calculate contribution to global correlations from one unit cell + /// Result is recorded in ClexParamPack + void _calc_global_corr_contribution() const override; + + /// \brief Calculate contribution to global correlations from one unit cell /// Result is recorded in double array starting at corr_begin + void _calc_global_corr_contribution(double *corr_begin) const override; + + /// \brief Calculate contribution to select global correlations from one unit cell into ClexParamPack + /// Result is recorded in ClexParamPack + void _calc_restricted_global_corr_contribution(size_type const *ind_list_begin, size_type const *ind_list_end) const override; + + /// \brief Calculate contribution to select global correlations from one unit cell + /// Result is recorded in double array starting at corr_begin + void _calc_restricted_global_corr_contribution(double *corr_begin, size_type const *ind_list_begin, size_type const *ind_list_end) const override; + + /// \brief Calculate point correlations about neighbor site 'nlist_ind' + /// For global clexulators, 'nlist_ind' only ranges over sites in the cell + /// For local clexulators, 'nlist_ind' ranges over all sites in the neighborhood + /// Result is recorded in ClexParamPack + void _calc_point_corr(int nlist_ind) const override; + + /// \brief Calculate point correlations about neighbor site 'nlist_ind' + /// For global clexulators, 'nlist_ind' only ranges over sites in the cell + /// For local clexulators, 'nlist_ind' ranges over all sites in the neighborhood + /// Result is recorded in double array starting at corr_begin + void _calc_point_corr(int nlist_ind, double *corr_begin) const override; + + /// \brief Calculate select point correlations about neighbor site 'nlist_ind' + /// For global clexulators, 'nlist_ind' only ranges over sites in the cell + /// For local clexulators, 'nlist_ind' ranges over all sites in the neighborhood + /// Result is recorded in ClexParamPack + void _calc_restricted_point_corr(int nlist_ind, size_type const *ind_list_begin, size_type const *ind_list_end) const override; + + /// \brief Calculate select point correlations about neighbor site 'nlist_ind' + /// For global clexulators, 'nlist_ind' only ranges over sites in the cell + /// For local clexulators, 'nlist_ind' ranges over all sites in the neighborhood + /// Result is recorded in double array starting at corr_begin + void _calc_restricted_point_corr(int nlist_ind, double *corr_begin, size_type const *ind_list_begin, size_type const *ind_list_end) const override; + + /// \brief Calculate the change in point correlations due to changing an occupant at neighbor site 'nlist_ind' + /// For global clexulators, 'nlist_ind' only ranges over sites in the cell + /// For local clexulators, 'nlist_ind' ranges over all sites in the neighborhood + /// Result is recorded in ClexParamPack + void _calc_delta_point_corr(int nlist_ind, int occ_i, int occ_f) const override; + + /// \brief Calculate the change in point correlations due to changing an occupant at neighbor site 'nlist_ind' + /// For global clexulators, 'nlist_ind' only ranges over sites in the cell + /// For local clexulators, 'nlist_ind' ranges over all sites in the neighborhood + /// Result is recorded in double array starting at corr_begin + void _calc_delta_point_corr(int nlist_ind, int occ_i, int occ_f, double *corr_begin) const override; + + /// \brief Calculate the change in select point correlations due to changing an occupant at neighbor site 'nlist_ind' + /// For global clexulators, 'nlist_ind' only ranges over sites in the cell + /// For local clexulators, 'nlist_ind' ranges over all sites in the neighborhood + /// Result is recorded in ClexParamPack + void _calc_restricted_delta_point_corr(int nlist_ind, int occ_i, int occ_f, size_type const *ind_list_begin, size_type const *ind_list_end) const override; + + /// \brief Calculate the change in select point correlations due to changing an occupant at neighbor site 'nlist_ind' + /// For global clexulators, 'nlist_ind' only ranges over sites in the cell + /// For local clexulators, 'nlist_ind' ranges over all sites in the neighborhood + /// Result is recorded in double array starting at corr_begin + void _calc_restricted_delta_point_corr(int nlist_ind, int occ_i, int occ_f, double *corr_begin, size_type const *ind_list_begin, size_type const *ind_list_end) const override; + + template + void _global_prepare() const; + + template + void _point_prepare(int nlist_ind) const; + + // disp evaluators and accessors for basis site 0: + double eval_disp_var_0_0(const int &nlist_ind) const { + return m_local_dof_ptrs[m_disp_var_param_key.index()]->col(_l(nlist_ind))[0]; + } + + double eval_disp_var_0_1(const int &nlist_ind) const { + return m_local_dof_ptrs[m_disp_var_param_key.index()]->col(_l(nlist_ind))[1]; + } + + double eval_disp_var_0_2(const int &nlist_ind) const { + return m_local_dof_ptrs[m_disp_var_param_key.index()]->col(_l(nlist_ind))[2]; + } + + template + Scalar const &disp_var_0(const int &nlist_ind) const { + return ParamPack::Val::get(m_params, m_disp_var_param_key, 0, nlist_ind); + } + template + Scalar const &disp_var_1(const int &nlist_ind) const { + return ParamPack::Val::get(m_params, m_disp_var_param_key, 1, nlist_ind); + } + template + Scalar const &disp_var_2(const int &nlist_ind) const { + return ParamPack::Val::get(m_params, m_disp_var_param_key, 2, nlist_ind); + } + //default functions for basis function evaluation + template + Scalar zero_func() const { + return Scalar(0.0); + } + + template + Scalar zero_func(int, int) const { + return Scalar(0.0); + } + + + }; + + //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + DiffClexulatorSquareLatticeTest_Clexulator_default::DiffClexulatorSquareLatticeTest_Clexulator_default() : + BaseClexulator(1, 7, 1) { + m_disp_var_param_key = m_params.allocate("disp_var", 3, 1, true); + _register_local_dof("disp", m_disp_var_param_key.index()); + + + m_corr_param_key = m_params.allocate("corr", corr_size(), 1, false); + + m_orbit_func_table_0[0] = &DiffClexulatorSquareLatticeTest_Clexulator_default::eval_bfunc_0_0; + m_orbit_func_table_0[1] = &DiffClexulatorSquareLatticeTest_Clexulator_default::eval_bfunc_1_0; + m_orbit_func_table_0[2] = &DiffClexulatorSquareLatticeTest_Clexulator_default::eval_bfunc_1_1; + m_orbit_func_table_0[3] = &DiffClexulatorSquareLatticeTest_Clexulator_default::eval_bfunc_1_2; + m_orbit_func_table_0[4] = &DiffClexulatorSquareLatticeTest_Clexulator_default::eval_bfunc_1_3; + m_orbit_func_table_0[5] = &DiffClexulatorSquareLatticeTest_Clexulator_default::eval_bfunc_1_4; + m_orbit_func_table_0[6] = &DiffClexulatorSquareLatticeTest_Clexulator_default::eval_bfunc_1_5; + + + m_flower_func_table_0[0][0] = &DiffClexulatorSquareLatticeTest_Clexulator_default::zero_func; + m_flower_func_table_0[0][1] = &DiffClexulatorSquareLatticeTest_Clexulator_default::site_eval_bfunc_1_0_at_0; + m_flower_func_table_0[0][2] = &DiffClexulatorSquareLatticeTest_Clexulator_default::site_eval_bfunc_1_1_at_0; + m_flower_func_table_0[0][3] = &DiffClexulatorSquareLatticeTest_Clexulator_default::site_eval_bfunc_1_2_at_0; + m_flower_func_table_0[0][4] = &DiffClexulatorSquareLatticeTest_Clexulator_default::site_eval_bfunc_1_3_at_0; + m_flower_func_table_0[0][5] = &DiffClexulatorSquareLatticeTest_Clexulator_default::site_eval_bfunc_1_4_at_0; + m_flower_func_table_0[0][6] = &DiffClexulatorSquareLatticeTest_Clexulator_default::site_eval_bfunc_1_5_at_0; + + + m_delta_func_table_0[0][0] = &DiffClexulatorSquareLatticeTest_Clexulator_default::zero_func; + m_delta_func_table_0[0][1] = &DiffClexulatorSquareLatticeTest_Clexulator_default::zero_func; + m_delta_func_table_0[0][2] = &DiffClexulatorSquareLatticeTest_Clexulator_default::zero_func; + m_delta_func_table_0[0][3] = &DiffClexulatorSquareLatticeTest_Clexulator_default::zero_func; + m_delta_func_table_0[0][4] = &DiffClexulatorSquareLatticeTest_Clexulator_default::zero_func; + m_delta_func_table_0[0][5] = &DiffClexulatorSquareLatticeTest_Clexulator_default::zero_func; + m_delta_func_table_0[0][6] = &DiffClexulatorSquareLatticeTest_Clexulator_default::zero_func; + + + m_orbit_func_table_1[0] = &DiffClexulatorSquareLatticeTest_Clexulator_default::eval_bfunc_0_0; + m_orbit_func_table_1[1] = &DiffClexulatorSquareLatticeTest_Clexulator_default::eval_bfunc_1_0; + m_orbit_func_table_1[2] = &DiffClexulatorSquareLatticeTest_Clexulator_default::eval_bfunc_1_1; + m_orbit_func_table_1[3] = &DiffClexulatorSquareLatticeTest_Clexulator_default::eval_bfunc_1_2; + m_orbit_func_table_1[4] = &DiffClexulatorSquareLatticeTest_Clexulator_default::eval_bfunc_1_3; + m_orbit_func_table_1[5] = &DiffClexulatorSquareLatticeTest_Clexulator_default::eval_bfunc_1_4; + m_orbit_func_table_1[6] = &DiffClexulatorSquareLatticeTest_Clexulator_default::eval_bfunc_1_5; + + + m_flower_func_table_1[0][0] = &DiffClexulatorSquareLatticeTest_Clexulator_default::zero_func; + m_flower_func_table_1[0][1] = &DiffClexulatorSquareLatticeTest_Clexulator_default::site_eval_bfunc_1_0_at_0; + m_flower_func_table_1[0][2] = &DiffClexulatorSquareLatticeTest_Clexulator_default::site_eval_bfunc_1_1_at_0; + m_flower_func_table_1[0][3] = &DiffClexulatorSquareLatticeTest_Clexulator_default::site_eval_bfunc_1_2_at_0; + m_flower_func_table_1[0][4] = &DiffClexulatorSquareLatticeTest_Clexulator_default::site_eval_bfunc_1_3_at_0; + m_flower_func_table_1[0][5] = &DiffClexulatorSquareLatticeTest_Clexulator_default::site_eval_bfunc_1_4_at_0; + m_flower_func_table_1[0][6] = &DiffClexulatorSquareLatticeTest_Clexulator_default::site_eval_bfunc_1_5_at_0; + + + m_delta_func_table_1[0][0] = &DiffClexulatorSquareLatticeTest_Clexulator_default::zero_func; + m_delta_func_table_1[0][1] = &DiffClexulatorSquareLatticeTest_Clexulator_default::zero_func; + m_delta_func_table_1[0][2] = &DiffClexulatorSquareLatticeTest_Clexulator_default::zero_func; + m_delta_func_table_1[0][3] = &DiffClexulatorSquareLatticeTest_Clexulator_default::zero_func; + m_delta_func_table_1[0][4] = &DiffClexulatorSquareLatticeTest_Clexulator_default::zero_func; + m_delta_func_table_1[0][5] = &DiffClexulatorSquareLatticeTest_Clexulator_default::zero_func; + m_delta_func_table_1[0][6] = &DiffClexulatorSquareLatticeTest_Clexulator_default::zero_func; + + + m_weight_matrix.row(0) << 1, 0, 0; + m_weight_matrix.row(1) << 0, 1, 0; + m_weight_matrix.row(2) << 0, 0, 11; + + m_sublat_indices = std::set{0}; + + m_n_sublattices = 1; + + m_neighborhood = std::set { + xtal::UnitCell(0, 0, 0) + }; + + + m_orbit_neighborhood.resize(corr_size()); + m_orbit_site_neighborhood.resize(corr_size()); + m_orbit_neighborhood[1] = std::set { + xtal::UnitCell(0, 0, 0) + }; + m_orbit_neighborhood[2] = m_orbit_neighborhood[1]; + m_orbit_neighborhood[3] = m_orbit_neighborhood[1]; + m_orbit_neighborhood[4] = m_orbit_neighborhood[1]; + m_orbit_neighborhood[5] = m_orbit_neighborhood[1]; + m_orbit_neighborhood[6] = m_orbit_neighborhood[1]; + + m_orbit_site_neighborhood[1] = std::set { + xtal::UnitCellCoord(0, 0, 0, 0) + }; + m_orbit_site_neighborhood[2] = m_orbit_site_neighborhood[1]; + m_orbit_site_neighborhood[3] = m_orbit_site_neighborhood[1]; + m_orbit_site_neighborhood[4] = m_orbit_site_neighborhood[1]; + m_orbit_site_neighborhood[5] = m_orbit_site_neighborhood[1]; + m_orbit_site_neighborhood[6] = m_orbit_site_neighborhood[1]; + + } + + + DiffClexulatorSquareLatticeTest_Clexulator_default::~DiffClexulatorSquareLatticeTest_Clexulator_default() { + //nothing here for now + } + + /// \brief Calculate contribution to global correlations from one unit cell + void DiffClexulatorSquareLatticeTest_Clexulator_default::_calc_global_corr_contribution(double *corr_begin) const { + _calc_global_corr_contribution(); + for(size_type i = 0; i < corr_size(); i++) { + *(corr_begin + i) = ParamPack::Val::get(m_params, m_corr_param_key, i); + } + } + + /// \brief Calculate contribution to global correlations from one unit cell + void DiffClexulatorSquareLatticeTest_Clexulator_default::_calc_global_corr_contribution() const { + m_params.pre_eval(); + if(m_params.eval_mode() == ParamPack::DEFAULT) { + _global_prepare(); + for(size_type i = 0; i < corr_size(); i++) { + ParamPack::Val::set(m_params, m_corr_param_key, i, (this->*m_orbit_func_table_0[i])()); + } + } + else if(m_params.eval_mode() == ParamPack::DIFF) { + _global_prepare(); + for(size_type i = 0; i < corr_size(); i++) { + ParamPack::Val::set(m_params, m_corr_param_key, i, (this->*m_orbit_func_table_1[i])()); + } + } + m_params.post_eval(); + } + + /// \brief Calculate contribution to select global correlations from one unit cell + void DiffClexulatorSquareLatticeTest_Clexulator_default::_calc_restricted_global_corr_contribution(double *corr_begin, size_type const *ind_list_begin, size_type const *ind_list_end) const { + _calc_restricted_global_corr_contribution(ind_list_begin, ind_list_end); + for(; ind_list_begin < ind_list_end; ind_list_begin++) { + *(corr_begin + *ind_list_begin) = ParamPack::Val::get(m_params, m_corr_param_key, *ind_list_begin); + } + } + + /// \brief Calculate contribution to select global correlations from one unit cell + void DiffClexulatorSquareLatticeTest_Clexulator_default::_calc_restricted_global_corr_contribution(size_type const *ind_list_begin, size_type const *ind_list_end) const { + m_params.pre_eval(); + if(m_params.eval_mode() == ParamPack::DEFAULT) { + _global_prepare(); + for(; ind_list_begin < ind_list_end; ind_list_begin++) { + ParamPack::Val::set(m_params, m_corr_param_key, *ind_list_begin, (this->*m_orbit_func_table_0[*ind_list_begin])()); + } + } + else if(m_params.eval_mode() == ParamPack::DIFF) { + _global_prepare(); + for(; ind_list_begin < ind_list_end; ind_list_begin++) { + ParamPack::Val::set(m_params, m_corr_param_key, *ind_list_begin, (this->*m_orbit_func_table_1[*ind_list_begin])()); + } + } + m_params.post_eval(); + } + + /// \brief Calculate point correlations about basis site 'nlist_ind' + void DiffClexulatorSquareLatticeTest_Clexulator_default::_calc_point_corr(int nlist_ind, double *corr_begin) const { + _calc_point_corr(nlist_ind); + for(size_type i = 0; i < corr_size(); i++) { + *(corr_begin + i) = ParamPack::Val::get(m_params, m_corr_param_key, i); + } + } + + /// \brief Calculate point correlations about basis site 'nlist_ind' + void DiffClexulatorSquareLatticeTest_Clexulator_default::_calc_point_corr(int nlist_ind) const { + m_params.pre_eval(); + if(m_params.eval_mode() == ParamPack::DEFAULT) { + _point_prepare(nlist_ind); + for(size_type i = 0; i < corr_size(); i++) { + ParamPack::Val::set(m_params, m_corr_param_key, i, (this->*m_flower_func_table_0[nlist_ind][i])()); + } + } + else if(m_params.eval_mode() == ParamPack::DIFF) { + _point_prepare(nlist_ind); + for(size_type i = 0; i < corr_size(); i++) { + ParamPack::Val::set(m_params, m_corr_param_key, i, (this->*m_flower_func_table_1[nlist_ind][i])()); + } + } + m_params.post_eval(); + } + + /// \brief Calculate select point correlations about basis site 'nlist_ind' + void DiffClexulatorSquareLatticeTest_Clexulator_default::_calc_restricted_point_corr(int nlist_ind, double *corr_begin, size_type const *ind_list_begin, size_type const *ind_list_end) const { + _calc_restricted_point_corr(nlist_ind, ind_list_begin, ind_list_end); + for(; ind_list_begin < ind_list_end; ind_list_begin++) { + *(corr_begin + *ind_list_begin) = ParamPack::Val::get(m_params, m_corr_param_key, *ind_list_begin); + } + } + + /// \brief Calculate select point correlations about basis site 'nlist_ind' + void DiffClexulatorSquareLatticeTest_Clexulator_default::_calc_restricted_point_corr(int nlist_ind, size_type const *ind_list_begin, size_type const *ind_list_end) const { + m_params.pre_eval(); + if(m_params.eval_mode() == ParamPack::DEFAULT) { + _point_prepare(nlist_ind); + for(; ind_list_begin < ind_list_end; ind_list_begin++) { + ParamPack::Val::set(m_params, m_corr_param_key, *ind_list_begin, (this->*m_flower_func_table_0[nlist_ind][*ind_list_begin])()); + } + } + else if(m_params.eval_mode() == ParamPack::DIFF) { + _point_prepare(nlist_ind); + for(; ind_list_begin < ind_list_end; ind_list_begin++) { + ParamPack::Val::set(m_params, m_corr_param_key, *ind_list_begin, (this->*m_flower_func_table_1[nlist_ind][*ind_list_begin])()); + } + } + m_params.post_eval(); + } + + /// \brief Calculate the change in point correlations due to changing an occupant + void DiffClexulatorSquareLatticeTest_Clexulator_default::_calc_delta_point_corr(int nlist_ind, int occ_i, int occ_f, double *corr_begin) const { + _calc_delta_point_corr(nlist_ind, occ_i, occ_f); + for(size_type i = 0; i < corr_size(); i++) { + *(corr_begin + i) = ParamPack::Val::get(m_params, m_corr_param_key, i); + } + } + + /// \brief Calculate the change in point correlations due to changing an occupant + void DiffClexulatorSquareLatticeTest_Clexulator_default::_calc_delta_point_corr(int nlist_ind, int occ_i, int occ_f) const { + m_params.pre_eval(); + if(m_params.eval_mode() == ParamPack::DEFAULT) { + _point_prepare(nlist_ind); + for(size_type i = 0; i < corr_size(); i++) { + ParamPack::Val::set(m_params, m_corr_param_key, i, (this->*m_delta_func_table_0[nlist_ind][i])(occ_i, occ_f)); + } + } + else if(m_params.eval_mode() == ParamPack::DIFF) { + _point_prepare(nlist_ind); + for(size_type i = 0; i < corr_size(); i++) { + ParamPack::Val::set(m_params, m_corr_param_key, i, (this->*m_delta_func_table_1[nlist_ind][i])(occ_i, occ_f)); + } + } + m_params.post_eval(); + } + + /// \brief Calculate the change in select point correlations due to changing an occupant + void DiffClexulatorSquareLatticeTest_Clexulator_default::_calc_restricted_delta_point_corr(int nlist_ind, int occ_i, int occ_f, double *corr_begin, size_type const *ind_list_begin, size_type const *ind_list_end) const { + _calc_restricted_delta_point_corr(nlist_ind, occ_i, occ_f, ind_list_begin, ind_list_end); + for(; ind_list_begin < ind_list_end; ind_list_begin++) { + *(corr_begin + *ind_list_begin) = ParamPack::Val::get(m_params, m_corr_param_key, *ind_list_begin); + } + } + + /// \brief Calculate the change in select point correlations due to changing an occupant + void DiffClexulatorSquareLatticeTest_Clexulator_default::_calc_restricted_delta_point_corr(int nlist_ind, int occ_i, int occ_f, size_type const *ind_list_begin, size_type const *ind_list_end) const { + m_params.pre_eval(); + if(m_params.eval_mode() == ParamPack::DEFAULT) { + _point_prepare(nlist_ind); + for(; ind_list_begin < ind_list_end; ind_list_begin++) { + ParamPack::Val::set(m_params, m_corr_param_key, *ind_list_begin, (this->*m_delta_func_table_0[nlist_ind][*ind_list_begin])(occ_i, occ_f)); + } + } + else if(m_params.eval_mode() == ParamPack::DIFF) { + _point_prepare(nlist_ind); + for(; ind_list_begin < ind_list_end; ind_list_begin++) { + ParamPack::Val::set(m_params, m_corr_param_key, *ind_list_begin, (this->*m_delta_func_table_1[nlist_ind][*ind_list_begin])(occ_i, occ_f)); + } + } + m_params.post_eval(); + } + + + template + void DiffClexulatorSquareLatticeTest_Clexulator_default::_point_prepare(int nlist_ind) const { + switch(nlist_ind) { + case 0: + if(m_params.eval_mode(m_disp_var_param_key) != ParamPack::READ) { + ParamPack::Val::set(m_params, m_disp_var_param_key, 0, 0, eval_disp_var_0_0(0)); + ParamPack::Val::set(m_params, m_disp_var_param_key, 1, 0, eval_disp_var_0_1(0)); + ParamPack::Val::set(m_params, m_disp_var_param_key, 2, 0, eval_disp_var_0_2(0)); + } + break; + } + } + template + void DiffClexulatorSquareLatticeTest_Clexulator_default::_global_prepare() const { + if(m_params.eval_mode(m_disp_var_param_key) != ParamPack::READ) { + ParamPack::Val::set(m_params, m_disp_var_param_key, 0, 0, eval_disp_var_0_0(0)); + ParamPack::Val::set(m_params, m_disp_var_param_key, 1, 0, eval_disp_var_0_1(0)); + ParamPack::Val::set(m_params, m_disp_var_param_key, 2, 0, eval_disp_var_0_2(0)); + } + } + + // Basis functions for empty cluster: + template + Scalar DiffClexulatorSquareLatticeTest_Clexulator_default::eval_bfunc_0_0() const { + return 1; + } + + /**** Basis functions for orbit 1**** +0.0000000 0.0000000 0.0000000 A + + ****/ + template + Scalar DiffClexulatorSquareLatticeTest_Clexulator_default::eval_bfunc_1_0() const { + return (0.707107 * pow(disp_var_0(0), 2)+0.707107 * pow(disp_var_1(0), 2)); + } + template + Scalar DiffClexulatorSquareLatticeTest_Clexulator_default::eval_bfunc_1_1() const { + return pow(disp_var_2(0), 2); + } + template + Scalar DiffClexulatorSquareLatticeTest_Clexulator_default::eval_bfunc_1_2() const { + return (0.707107 * pow(disp_var_0(0), 4)+0.707107 * pow(disp_var_1(0), 4)); + } + template + Scalar DiffClexulatorSquareLatticeTest_Clexulator_default::eval_bfunc_1_3() const { + return 2.44949 * pow(disp_var_0(0), 2) * pow(disp_var_1(0), 2); + } + template + Scalar DiffClexulatorSquareLatticeTest_Clexulator_default::eval_bfunc_1_4() const { + return (1.73205 * pow(disp_var_0(0), 2) * pow(disp_var_2(0), 2)+1.73205 * pow(disp_var_1(0), 2) * pow(disp_var_2(0), 2)); + } + template + Scalar DiffClexulatorSquareLatticeTest_Clexulator_default::eval_bfunc_1_5() const { + return pow(disp_var_2(0), 4); + } + + template + Scalar DiffClexulatorSquareLatticeTest_Clexulator_default::site_eval_bfunc_1_0_at_0() const { + return (0.707107 * pow(disp_var_0(0), 2)+0.707107 * pow(disp_var_1(0), 2)); + } + template + Scalar DiffClexulatorSquareLatticeTest_Clexulator_default::site_eval_bfunc_1_1_at_0() const { + return pow(disp_var_2(0), 2); + } + template + Scalar DiffClexulatorSquareLatticeTest_Clexulator_default::site_eval_bfunc_1_2_at_0() const { + return (0.707107 * pow(disp_var_0(0), 4)+0.707107 * pow(disp_var_1(0), 4)); + } + template + Scalar DiffClexulatorSquareLatticeTest_Clexulator_default::site_eval_bfunc_1_3_at_0() const { + return 2.44949 * pow(disp_var_0(0), 2) * pow(disp_var_1(0), 2); + } + template + Scalar DiffClexulatorSquareLatticeTest_Clexulator_default::site_eval_bfunc_1_4_at_0() const { + return (1.73205 * pow(disp_var_0(0), 2) * pow(disp_var_2(0), 2)+1.73205 * pow(disp_var_1(0), 2) * pow(disp_var_2(0), 2)); + } + template + Scalar DiffClexulatorSquareLatticeTest_Clexulator_default::site_eval_bfunc_1_5_at_0() const { + return pow(disp_var_2(0), 4); + } + +} // namespace clexulator +} // namespace CASM + + +extern "C" { + /// \brief Returns a clexulator::BaseClexulator* owning a DiffClexulatorSquareLatticeTest_Clexulator_default + CASM::clexulator::BaseClexulator *make_DiffClexulatorSquareLatticeTest_Clexulator_default() { + return new CASM::clexulator::DiffClexulatorSquareLatticeTest_Clexulator_default(); + } + +} + diff --git a/python/tests/clexulator/data/square_lattice_prim.json b/python/tests/clexulator/data/square_lattice_prim.json new file mode 100644 index 0000000..190258c --- /dev/null +++ b/python/tests/clexulator/data/square_lattice_prim.json @@ -0,0 +1,14 @@ +{ + "basis": [ + { + "coordinate": [0.000000000000, 0.000000000000, 0.000000000000], + "occupants": ["A"] + }], + "coordinate_mode": "Fractional", + "lattice_vectors": [ + [3.000000000000, 0.00000000000, 0.00000000000], + [0.00000000000, 3.000000000000, 0.00000000000], + [0.000000000000, 0.000000000000, 10.000000000000] + ], + "title": "A" +} diff --git a/python/tests/clexulator/test_diff_clex.py b/python/tests/clexulator/test_diff_clex.py new file mode 100644 index 0000000..975aed4 --- /dev/null +++ b/python/tests/clexulator/test_diff_clex.py @@ -0,0 +1,159 @@ +import pytest +import numpy as np +import libcasm.xtal +import libcasm.clexulator +from .functions import make_source + + +@pytest.fixture +def diff_clexulator(session_shared_datadir): + source = make_source( + session_shared_datadir, "DiffClexulatorSquareLatticeTest", "default" + ) + + prim_neighbor_list = libcasm.clexulator.PrimNeighborList() + clexulator = libcasm.clexulator.make_clexulator(source, prim_neighbor_list) + return clexulator, prim_neighbor_list + + +def test_make_diff_clexulator(diff_clexulator): + """Tests if diff clexulator made had the + correct number of basis functions + + There should be 7 basis functions + + Generated using casm1.X + + .. math:: + \\Phi_{0} = 1 + \\Phi_{1} = \\sqrt{1/2}(dx_{0}^{2} +dy_{0}^{2} ) + \\Phi_{2} = dz_{0}^{2} + \\Phi_{3} = \\sqrt{1/2}(dx_{0}^{4} +dy_{0}^{4} ) + \\Phi_{4} = \\sqrt{6}dx_{0}^{2} dy_{0}^{2} + \\Phi_{5} = \\sqrt{3}(dx_{0}^{2} dz_{0}^{2} +dy_{0}^{2} dz_{0}^{2} ) + \\Phi_{6} = dz_{0}^{4} + + Parameters + ---------- + diff_clexulator : libcasm.clexulator.Clexulator + + """ + clexulator, _ = diff_clexulator + assert isinstance(clexulator, libcasm.clexulator.Clexulator) + assert clexulator.n_functions() == 7 + + +def test_gradients_of_correlations(diff_clexulator): + """Tests if the gradients of correlations being made + matches the analytical gradients + + Generated using casm1.X + + .. math:: + \\Phi_{0} = 1 + \\Phi_{1} = \\sqrt{1/2}(dx_{0}^{2} +dy_{0}^{2} ) + \\Phi_{2} = dz_{0}^{2} + \\Phi_{3} = \\sqrt{1/2}(dx_{0}^{4} +dy_{0}^{4} ) + \\Phi_{4} = \\sqrt{6}dx_{0}^{2} dy_{0}^{2} + \\Phi_{5} = \\sqrt{3}(dx_{0}^{2} dz_{0}^{2} +dy_{0}^{2} dz_{0}^{2} ) + \\Phi_{6} = dz_{0}^{4} + + Parameters + ---------- + diff_clexulator : libcasm.clexulator.Clexulator + + """ + + clexulator, prim_neighbor_list = diff_clexulator + + # make a configuration of supercell size 1 + T = np.eye(3) + + # assign disp site dofs to the configuration + # There is only one site in volume 1 supercell + # dof_values correspond to dx, dy, dz of that site + dof_values = np.array([1.0, -0.5, 0.5]) + config_dof_values = libcasm.clexulator.ConfigDoFValues() + config_dof_values.set_occupation([0]) + config_dof_values.insert_or_assign_local_dof_values("disp", dof_values) + + # make SuperNeighborList + supercell_neighbor_list = libcasm.clexulator.SuperNeighborList( + T, prim_neighbor_list + ) + + # make correlations + correlations = libcasm.clexulator.Correlations( + supercell_neighbor_list, clexulator, config_dof_values + ) + + # query gradients of correlations for this config + gradients_of_correlations = correlations.grad_correlations("disp") + + # the shape of the gradients of correlations matrix should be 3x7 + # because there are 3 degrees of freedom (dx, dy, dz) and 7 correlations + # you will have derivatives with each of the correlations with each of the degrees of freedom + + assert gradients_of_correlations.shape == (3, 7) + + # comparing the values with analytical derivates + dx = dof_values[0] + dy = dof_values[1] + dz = dof_values[2] + + # first row of gradients_of_correlations will be derviatives of all the 7 correlations with dx + derviatives_with_dx = gradients_of_correlations[0, :] + + # derivative of \\Phi_{0} with dx will be 0.0 + assert np.allclose(derviatives_with_dx[0], 0.0) + # derivative of \\Phi_{1} with dx will be \\sqrt{2} * dx + assert np.allclose(derviatives_with_dx[1], np.sqrt(2) * dx) + # derivative of \\Phi_{2} with dx will be 0.0 + assert np.allclose(derviatives_with_dx[2], 0.0) + # derivative of \\Phi_{3} with dx will be 2 * \\sqrt(2) * dx**3 + assert np.allclose(derviatives_with_dx[3], 2 * np.sqrt(2) * np.power(dx, 3)) + # derivative of \\Phi_{4} with dx will be 2 * \\sqrt(6) * dx * dy**2 + assert np.allclose(derviatives_with_dx[4], 2 * np.sqrt(6) * dx * np.power(dy, 2)) + # derivative of \\Phi_{5} with dx will be 2 * \\sqrt(3) * dx * dz**2 + assert np.allclose(derviatives_with_dx[5], 2 * np.sqrt(3) * dx * np.power(dz, 2)) + # derivative of \\Phi_{6} with dx will be 0.0 + assert np.allclose(derviatives_with_dx[6], 0.0) + + # second row of gradients_of_correlations will be derviatives of all the 7 correlations with dy + derviatives_with_dy = gradients_of_correlations[1, :] + + # derivative of \\Phi_{0} with dy will be 0.0 + assert np.allclose(derviatives_with_dy[0], 0.0) + # derivative of \\Phi_{1} with dy will be \\sqrt{2} * dy + assert np.allclose(derviatives_with_dy[1], np.sqrt(2) * dy) + # derivative of \\Phi_{2} with dy will be 0.0 + assert np.allclose(derviatives_with_dy[2], 0.0) + # derivative of \\Phi_{3} with dy will be 2 * \\sqrt(2) * dy**3 + assert np.allclose(derviatives_with_dy[3], 2 * np.sqrt(2) * np.power(dy, 3)) + # derivative of \\Phi_{4} with dy will be 2 * \\sqrt(6) * dx**2 * dy + assert np.allclose(derviatives_with_dy[4], 2 * np.sqrt(6) * np.power(dx, 2) * dy) + # derivative of \\Phi_{5} with dy will be 2 * \\sqrt(3) * dy * dz**2 + assert np.allclose(derviatives_with_dy[5], 2 * np.sqrt(3) * dy * np.power(dz, 2)) + # derivative of \\Phi_{6} with dy will be 0.0 + assert np.allclose(derviatives_with_dy[6], 0.0) + + # third row of gradients_of_correlations will be derviatives of all the 7 correlations with dz + derviatives_with_dz = gradients_of_correlations[2, :] + + # derivative of \\Phi_{0} with dz will be 0.0 + assert np.allclose(derviatives_with_dz[0], 0.0) + # derivative of \\Phi_{1} with dz will be 0.0 + assert np.allclose(derviatives_with_dz[1], 0.0) + # derivative of \\Phi_{2} with dz will be 2 * dz + assert np.allclose(derviatives_with_dz[2], 2 * dz) + # derivative of \\Phi_{3} with dz will be 0.0 + assert np.allclose(derviatives_with_dz[3], 0.0) + # derivative of \\Phi_{4} with dz will be 0.0 + assert np.allclose(derviatives_with_dz[4], 0.0) + # derivative of \\Phi_{5} with dz will be 2 * \\sqrt(3) * dz (dx**2 + dy**2) + assert np.allclose( + derviatives_with_dz[5], + 2 * np.sqrt(3) * dz * (np.power(dx, 2) + np.power(dy, 2)), + ) + # derivative of \\Phi_{6} with dz will be 4 * dz**3 + assert np.allclose(derviatives_with_dz[6], 4 * np.power(dz, 3)) diff --git a/src/casm/clexulator/Correlations.cc b/src/casm/clexulator/Correlations.cc index 71ad9db..9bfd418 100644 --- a/src/casm/clexulator/Correlations.cc +++ b/src/casm/clexulator/Correlations.cc @@ -1,7 +1,9 @@ #include "casm/clexulator/Correlations.hh" +#include "casm/clexulator/ClexParamPack.hh" #include "casm/clexulator/Clexulator.hh" #include "casm/clexulator/NeighborList.hh" +#include "casm/crystallography/AnisoValTraits.hh" namespace CASM { @@ -613,5 +615,89 @@ std::vector all_correlation_indices(Index n) { return corr_indices; } +/// Gradients of Correlations + +Eigen::MatrixXd const Correlations::grad_correlations(DoFKey const &key) { + ClexParamKey paramkey; + clexulator::ClexParamKey corr_key(m_clexulator->param_pack().key("corr")); + clexulator::ClexParamKey dof_key; + + if (key == "occ") { + paramkey = + m_clexulator->param_pack().key("diff/corr/" + key + "_site_func"); + dof_key = m_clexulator->param_pack().key("occ_site_func"); + } else { + paramkey = m_clexulator->param_pack().key("diff/corr/" + key + "_var"); + dof_key = m_clexulator->param_pack().key(key + "_var"); + } + + std::string em_corr, em_dof; + em_corr = m_clexulator->param_pack().eval_mode(corr_key); + em_dof = m_clexulator->param_pack().eval_mode(dof_key); + + // this const_cast is not great... + // but it seems like the only place passing const Clexulator is a problem and + // it is not actually changing clexulator before/after this function + + const_cast(*m_clexulator) + .param_pack() + .set_eval_mode(corr_key, "DIFF"); + const_cast(*m_clexulator) + .param_pack() + .set_eval_mode(dof_key, "DIFF"); + + Eigen::MatrixXd gcorr; + Index scel_vol = m_supercell_neighbor_list->n_unitcells(); + if (AnisoValTraits(key).global()) { + Eigen::MatrixXd gcorr_func = m_dof_values->global_dof_values.at(key); + gcorr.setZero(gcorr_func.size(), m_clexulator->corr_size()); + // Holds contribution to global correlations from a particular Neighborhood + + // std::vector corr(clexulator.corr_size(), 0.0); + for (int v = 0; v < scel_vol; v++) { + // Fill up contributions + m_clexulator->calc_global_corr_contribution( + *m_dof_values, m_supercell_neighbor_list->sites(v).data()); + + for (Index c = 0; c < m_clexulator->corr_size(); ++c) + gcorr.col(c) += m_clexulator->param_pack().read(paramkey(c)); + } + } else { + Eigen::MatrixXd gcorr_func; + gcorr.setZero(m_dof_values->local_dof_values.at(key).size(), + m_clexulator->corr_size()); + // Holds contribution to global correlations from a particular Neighborhood + Index l; + for (int v = 0; v < scel_vol; v++) { + // Fill up contributions + m_clexulator->calc_global_corr_contribution( + *m_dof_values, m_supercell_neighbor_list->sites(v).data()); + + for (Index c = 0; c < m_clexulator->corr_size(); ++c) { + gcorr_func = m_clexulator->param_pack().read(paramkey(c)); + + for (Index n = 0; n < m_clexulator->nlist_size(); ++n) { + l = m_supercell_neighbor_list->sites(v)[n]; + // for(Index i=0; i(*m_clexulator) + .param_pack() + .set_eval_mode(corr_key, em_corr); + const_cast(*m_clexulator) + .param_pack() + .set_eval_mode(dof_key, em_dof); + + return gcorr; +} + } // namespace clexulator } // namespace CASM diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 3d88e8b..41be781 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -100,8 +100,8 @@ endif() ### libcasm_testing ### set( libcasm_testing_SOURCES - ${PROJECT_SOURCE_DIR}/unit/autotools.cc ${PROJECT_SOURCE_DIR}/unit/testdir.cc + ${PROJECT_SOURCE_DIR}/unit/autotools.cc ) add_library(casm_testing SHARED ${libcasm_testing_SOURCES}) target_include_directories(casm_testing @@ -127,19 +127,19 @@ enable_testing() # casm_unit_clexulator add_executable(casm_unit_clexulator ${PROJECT_SOURCE_DIR}/unit/gtest_main_run_all.cpp + ${PROJECT_SOURCE_DIR}/unit/clexulator/Clexulator_json_io_test.cpp ${PROJECT_SOURCE_DIR}/unit/clexulator/ConfigDoFValuesTools_test.cpp + ${PROJECT_SOURCE_DIR}/unit/clexulator/OccClexulator_test.cpp + ${PROJECT_SOURCE_DIR}/unit/clexulator/version_test.cpp + ${PROJECT_SOURCE_DIR}/unit/clexulator/DoFSpace_test.cpp ${PROJECT_SOURCE_DIR}/unit/clexulator/ConfigDoFValues_json_io_test.cpp - ${PROJECT_SOURCE_DIR}/unit/clexulator/OrderParameter_test.cpp ${PROJECT_SOURCE_DIR}/unit/clexulator/LocalClusterExpansion_test.cpp - ${PROJECT_SOURCE_DIR}/unit/clexulator/ClusterExpansion_test.cpp ${PROJECT_SOURCE_DIR}/unit/clexulator/LocalOccClexulator_test.cpp - ${PROJECT_SOURCE_DIR}/unit/clexulator/ConfigDoFValues_test.cpp - ${PROJECT_SOURCE_DIR}/unit/clexulator/DoFSpace_test.cpp ${PROJECT_SOURCE_DIR}/unit/clexulator/DoFSpace_values_test.cpp - ${PROJECT_SOURCE_DIR}/unit/clexulator/version_test.cpp - ${PROJECT_SOURCE_DIR}/unit/clexulator/OccClexulator_test.cpp - ${PROJECT_SOURCE_DIR}/unit/clexulator/Clexulator_json_io_test.cpp ${PROJECT_SOURCE_DIR}/unit/clexulator/DoFSpace_json_io_test.cpp + ${PROJECT_SOURCE_DIR}/unit/clexulator/ClusterExpansion_test.cpp + ${PROJECT_SOURCE_DIR}/unit/clexulator/OrderParameter_test.cpp + ${PROJECT_SOURCE_DIR}/unit/clexulator/ConfigDoFValues_test.cpp ${PROJECT_SOURCE_DIR}/unit/clexulator/SparseCoefficients_json_io_test.cpp ) target_link_libraries(casm_unit_clexulator