|
| 1 | + |
| 2 | +#include "plumed/tools/SwitchingFunction.h" |
| 3 | +#include <fstream> |
| 4 | +#include <iostream> |
| 5 | +#include <iomanip> |
| 6 | +#include <ostream> |
| 7 | +#include <utility> |
| 8 | + |
| 9 | + |
| 10 | +using PLMD::SwitchingFunction; |
| 11 | +//using PLMD::SwitchingFunctionAccelerable; |
| 12 | + |
| 13 | + |
| 14 | +struct callSQR { |
| 15 | + template<typename T> |
| 16 | + static std::pair <double,double> call(const T& sw, double point) { |
| 17 | + double deriv; |
| 18 | + point *= point; |
| 19 | + double value = sw.calculateSqr(point,deriv); |
| 20 | + return std::make_pair(value,deriv); |
| 21 | + } |
| 22 | +}; |
| 23 | + |
| 24 | +struct callplain { |
| 25 | + template<typename T> |
| 26 | + static std::pair <double,double> call(const T& sw, const double point) { |
| 27 | + double deriv; |
| 28 | + double value = sw.calculate(point,deriv); |
| 29 | + return std::make_pair(value,deriv); |
| 30 | + } |
| 31 | +}; |
| 32 | + |
| 33 | +template<typename caller,typename SW> |
| 34 | +std::ostream& calculateAndSubtract(const SW& sw, |
| 35 | + const double point, |
| 36 | + const double value, |
| 37 | + const double derivative, |
| 38 | + std::ostream& out) { |
| 39 | + auto [myval,mydev] = caller::call(sw,point); |
| 40 | + out << (myval -value)<< " " << (mydev - derivative)<<" " ; |
| 41 | + return out; |
| 42 | +} |
| 43 | + |
| 44 | + |
| 45 | +void test(const std::string& name, std::string definition,std::string custom,bool stretch=true) { |
| 46 | + if (!stretch) { |
| 47 | + definition += " NOSTRETCH "; |
| 48 | + custom += " NOSTRETCH "; |
| 49 | + } |
| 50 | + |
| 51 | + std::ofstream os("out_"+name+(stretch ? "" : "_nostretch")); |
| 52 | + os <<std::fixed<<std::setprecision(6); |
| 53 | + SwitchingFunction sw; |
| 54 | + std::string error; |
| 55 | + sw.set(definition,error); |
| 56 | + if (!error.empty()) { |
| 57 | + std::cerr<<error<<"\n"; |
| 58 | + } |
| 59 | + error.clear(); |
| 60 | + |
| 61 | + SwitchingFunction swc; |
| 62 | + swc.set(custom,error); |
| 63 | + if (!error.empty()) { |
| 64 | + std::cerr<<error<<"\n"; |
| 65 | + } |
| 66 | + error.clear(); |
| 67 | + /* |
| 68 | + SwitchingFunctionAccelerable swa; |
| 69 | + swa.set(definition,error); |
| 70 | + if (!error.empty()) { |
| 71 | + std::cerr<<error<<"\n"; |
| 72 | + } |
| 73 | + */ |
| 74 | + os << "point : value deriv vsq_delta dsq_delta " |
| 75 | + "vAcc_delta dAcc_delta " |
| 76 | + //"vAccsq_delta dAccsq_delta " |
| 77 | + "vcus_delta dcus_delta " |
| 78 | + "vsqcus_delta dsqcus_delta\n"; |
| 79 | + for (int i=0; i<10; i++) { |
| 80 | + double point=i/2.50; |
| 81 | + double derivative; |
| 82 | + double value; |
| 83 | + os<< point <<" : "; |
| 84 | + std::tie(value,derivative) = callplain::call(sw, point); |
| 85 | + os << value << " " << derivative <<" "; |
| 86 | + calculateAndSubtract<callSQR>(sw,point, value, derivative,os) << " "; |
| 87 | + //calculateAndSubtract<callplain>(swa,point, value, derivative,os) << " "; |
| 88 | + //calculateAndSubtract<callSQR>(swa,point, value, derivative,os) << " "; |
| 89 | + calculateAndSubtract<callplain>(swc,point, value, derivative,os) << " "; |
| 90 | + calculateAndSubtract<callSQR>(swc,point, value, derivative,os) << "\n"; |
| 91 | + } |
| 92 | +} |
| 93 | + |
| 94 | + |
| 95 | +int main() { |
| 96 | +//the test will output the value of the switch along with de difference of it and |
| 97 | +//the same calculation under different settings (with squared/non squared input, |
| 98 | +//with the "accelerable version and with custom) |
| 99 | + for(const auto &x: { |
| 100 | + std::tuple<std::string,std::string,std::string> { |
| 101 | + "COSINUS R_0=2.6", |
| 102 | + "CUSTOM FUNC=0.5*(cos(x*pi)+1) R_0=2.6 D_MAX=2.6", |
| 103 | + "cosinus" |
| 104 | +}, { |
| 105 | + "EXP R_0=0.8 D_0=0.5 D_MAX=2.6", |
| 106 | + "CUSTOM FUNC=exp(-x) R_0=0.8 D_0=0.5 D_MAX=2.6", |
| 107 | + "exp" |
| 108 | +}, { |
| 109 | + "GAUSSIAN R_0=1.0 D_0=0.0 D_MAX=2.6", |
| 110 | + "CUSTOM FUNC=exp(-x^2/2) R_0=1.0 D_0=0.0 D_MAX=2.6", |
| 111 | + "fastgaussian" |
| 112 | +}, { |
| 113 | + "GAUSSIAN R_0=1.0 D_0=0.3 D_MAX=2.6", |
| 114 | + "CUSTOM FUNC=exp(-x^2/2) R_0=1.0 D_0=0.3 D_MAX=2.6", |
| 115 | + "gaussian" |
| 116 | +}, { |
| 117 | + "RATIONAL R_0=1.3 NN=6 MM=10 D_MAX=2.6", |
| 118 | + "CUSTOM FUNC=(1-x^6)/(1-x^10) D_MAX=2.6 R_0=1.3", |
| 119 | + "fastrational" |
| 120 | +}, { |
| 121 | + "RATIONAL R_0=1.3 D_MAX=2.6", |
| 122 | + "CUSTOM FUNC=(1-x^6)/(1-x^12) D_MAX=2.6 R_0=1.3", |
| 123 | + "fastrational_NNeq2MM" |
| 124 | +}, { |
| 125 | + "RATIONAL R_0=1.3 NN=5 MM=11 D_MAX=2.6", |
| 126 | + "CUSTOM FUNC=(1-x^5)/(1-x^11) D_MAX=2.6 R_0=1.3", |
| 127 | + "rational" |
| 128 | +}, { |
| 129 | + "RATIONAL R_0=1.3 NN=5 D_MAX=2.6", |
| 130 | + "CUSTOM FUNC=(1-x^5)/(1-x^10) D_MAX=2.6 R_0=1.3", |
| 131 | + "rational_NNeq2MM" |
| 132 | +}, { |
| 133 | + "Q R_0=1.0 D_0=0.3 BETA=5.0 LAMBDA=1.0 REF=1.3 D_MAX=2.6", |
| 134 | + "CUSTOM FUNC=1/(1+exp(5.0*((x+0.3)-1.3))) R_0=1.0 D_0=0.3 D_MAX=2.6", |
| 135 | + "q" |
| 136 | +}, { |
| 137 | + "TANH R_0=1.3 D_MAX=2.6", |
| 138 | + "CUSTOM FUNC=1-tanh(x) R_0=1.3 D_MAX=2.6", |
| 139 | + "tanh" |
| 140 | +}, { |
| 141 | + "SMAP R_0=1.3 A=3 B=2 D_MAX=2.6", |
| 142 | + "CUSTOM FUNC=(1+(2^(3/2)-1)*x^3)^(-2/3) R_0=1.3 D_MAX=2.6", |
| 143 | + "smap" |
| 144 | +}, { |
| 145 | + "CUBIC D_MAX=2.6 D_0=0.6", |
| 146 | +//Here R_O = DMAX-D_0 |
| 147 | + "CUSTOM FUNC=(x-1)^2*(1+2*x) R_0=2.0 D_0=0.6 D_MAX=2.6", |
| 148 | + "cubic" |
| 149 | +} |
| 150 | + }) { |
| 151 | + auto [definition, custom, name] = x; |
| 152 | + test (name,definition,custom); |
| 153 | + test (name,definition,custom,false); |
| 154 | + |
| 155 | + } |
| 156 | + |
| 157 | + return 0; |
| 158 | +} |
0 commit comments