Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2019 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "core/ActionRegister.h"
23 : #include "Function.h"
24 :
25 : #include <cmath>
26 :
27 : using namespace std;
28 :
29 : namespace PLMD {
30 : namespace function {
31 :
32 : //+PLUMEDOC FUNCTION PIECEWISE
33 : /*
34 : Compute a piecewise straight line through its arguments that passes through
35 : a set of ordered control points.
36 :
37 : For variables less than the first
38 : (greater than the last) point, the value of the first (last) point is used.
39 :
40 : \f[
41 : \frac{y_{i+1}-y_i}{x_{i+1}-x_i}(s-x_i)+y_i ; if x_i<s<x_{i+1}
42 : \f]
43 : \f[
44 : y_N ; if x>x_{N-1}
45 : \f]
46 : \f[
47 : y_1 ; if x<x_0
48 : \f]
49 :
50 : Control points are passed using the POINT0=... POINT1=... syntax as in the example below
51 :
52 : If one argument is supplied, it results in a scalar quantity.
53 : If multiple arguments are supplied, it results
54 : in a vector of values. Each value will be named as the name of the original
55 : argument with suffix _pfunc.
56 :
57 : \par Examples
58 :
59 : \plumedfile
60 : dist1: DISTANCE ATOMS=1,10
61 : dist2: DISTANCE ATOMS=2,11
62 :
63 : pw: PIECEWISE POINT0=1,10 POINT1=1,PI POINT2=3,10 ARG=dist1
64 : ppww: PIECEWISE POINT0=1,10 POINT1=1,PI POINT2=3,10 ARG=dist1,dist2
65 : PRINT ARG=pw,ppww.dist1_pfunc,ppww.dist2_pfunc
66 : \endplumedfile
67 :
68 :
69 : */
70 : //+ENDPLUMEDOC
71 :
72 :
73 6 : class Piecewise :
74 : public Function
75 : {
76 : std::vector<std::pair<double,double> > points;
77 : public:
78 : explicit Piecewise(const ActionOptions&);
79 : void calculate();
80 : static void registerKeywords(Keywords& keys);
81 : };
82 :
83 :
84 6455 : PLUMED_REGISTER_ACTION(Piecewise,"PIECEWISE")
85 :
86 4 : void Piecewise::registerKeywords(Keywords& keys) {
87 4 : Function::registerKeywords(keys);
88 8 : keys.use("ARG");
89 16 : keys.add("numbered","POINT","This keyword is used to specify the various points in the function above.");
90 12 : keys.reset_style("POINT","compulsory");
91 4 : componentsAreNotOptional(keys);
92 16 : keys.addOutputComponent("_pfunc","default","one or multiple instances of this quantity will be referenceable elsewhere "
93 : "in the input file. These quantities will be named with the arguments of the "
94 : "function followed by the character string _pfunc. These quantities tell the "
95 : "user the values of the piecewise functions of each of the arguments.");
96 4 : }
97 :
98 3 : Piecewise::Piecewise(const ActionOptions&ao):
99 : Action(ao),
100 4 : Function(ao)
101 : {
102 9 : for(int i=0;; i++) {
103 : std::vector<double> pp;
104 24 : if(!parseNumberedVector("POINT",i,pp) ) break;
105 9 : if(pp.size()!=2) error("points should be in x,y format");
106 9 : points.push_back(std::pair<double,double>(pp[0],pp[1]));
107 21 : if(i>0 && points[i].first<=points[i-1].first) error("points abscissas should be monotonously increasing");
108 9 : }
109 :
110 9 : for(unsigned i=0; i<getNumberOfArguments(); i++)
111 4 : if(getPntrToArgument(i)->isPeriodic())
112 2 : error("Cannot use PIECEWISE on periodic arguments");
113 :
114 2 : if(getNumberOfArguments()==1) {
115 1 : addValueWithDerivatives();
116 1 : setNotPeriodic();
117 : } else {
118 5 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
119 5 : addComponentWithDerivatives( getPntrToArgument(i)->getName()+"_pfunc" );
120 2 : getPntrToComponent(i)->setNotPeriodic();
121 : }
122 : }
123 2 : checkRead();
124 :
125 2 : log.printf(" on points:");
126 22 : for(unsigned i=0; i<points.size(); i++) log.printf(" (%f,%f)",points[i].first,points[i].second);
127 2 : log.printf("\n");
128 2 : }
129 :
130 10 : void Piecewise::calculate() {
131 40 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
132 : double val=getArgument(i);
133 : unsigned p=0;
134 96 : for(; p<points.size(); p++) {
135 33 : if(val<points[p].first) break;
136 : }
137 : double f,d;
138 15 : if(p==0) {
139 5 : f=points[0].second;
140 : d=0.0;
141 10 : } else if(p==points.size()) {
142 8 : f=points[points.size()-1].second;
143 : d=0.0;
144 : } else {
145 12 : double m=(points[p].second-points[p-1].second) / (points[p].first-points[p-1].first);
146 6 : f=m*(val-points[p-1].first)+points[p-1].second;
147 : d=m;
148 : }
149 15 : if(getNumberOfArguments()==1) {
150 5 : setValue(f);
151 : setDerivative(i,d);
152 : } else {
153 10 : Value* v=getPntrToComponent(i);
154 10 : v->set(f);
155 : v->addDerivative(i,d);
156 : }
157 : }
158 10 : }
159 :
160 : }
161 4839 : }
162 :
163 :
|
