TabulatedFunction.cpp 13.2 KB
Newer Older
1
2
3
4
5
6
7
8
/* -------------------------------------------------------------------------- *
 *                                   OpenMM                                   *
 * -------------------------------------------------------------------------- *
 * This is part of the OpenMM molecular simulation toolkit originating from   *
 * Simbios, the NIH National Center for Physics-Based Simulation of           *
 * Biological Structures at Stanford, funded under the NIH Roadmap for        *
 * Medical Research, grant U54 GM072970. See https://simtk.org.               *
 *                                                                            *
9
 * Portions copyright (c) 2014-2022 Stanford University and the Authors.      *
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
 * Authors: Peter Eastman                                                     *
 * Contributors:                                                              *
 *                                                                            *
 * Permission is hereby granted, free of charge, to any person obtaining a    *
 * copy of this software and associated documentation files (the "Software"), *
 * to deal in the Software without restriction, including without limitation  *
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
 * and/or sell copies of the Software, and to permit persons to whom the      *
 * Software is furnished to do so, subject to the following conditions:       *
 *                                                                            *
 * The above copyright notice and this permission notice shall be included in *
 * all copies or substantial portions of the Software.                        *
 *                                                                            *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
 * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
 * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
 * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
 * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
 * -------------------------------------------------------------------------- */

#include "openmm/TabulatedFunction.h"
#include "openmm/OpenMMException.h"

using namespace OpenMM;
using namespace std;

38
39
40
41
bool TabulatedFunction::getPeriodic() const {
    return periodic;
}

42
43
44
45
int TabulatedFunction::getUpdateCount() const {
    return updateCount;
}

46
47
48
Continuous1DFunction::Continuous1DFunction(const vector<double>& values, double min, double max, bool periodic) {
    this->periodic = periodic;
    setFunctionParameters(values, min, max);
49
50
}

51
void Continuous1DFunction::getFunctionParameters(vector<double>& values, double& min, double& max) const {
52
53
54
55
56
    values = this->values;
    min = this->min;
    max = this->max;
}

57
void Continuous1DFunction::setFunctionParameters(const vector<double>& values, double min, double max) {
58
59
    if (max <= min)
        throw OpenMMException("Continuous1DFunction: max <= min for a tabulated function.");
60
61
62
63
    int n = values.size();
    if (periodic) {
        if (n < 3)
            throw OpenMMException("Continuous1DFunction: a periodic tabulated function must have at least three points");
64
65
66
       // Note: value-matching at boundary is eventually checked at spline creation.
       // if (values[0] != values[n-1])
       //     throw OpenMMException("Continuous1DFunction: with periodic=true, the first and last points must have the same value");
67
68
69
    }
    else if (n < 2)
        throw OpenMMException("Continuous1DFunction: a non-periodic tabulated function must have at least two points");
70
71
72
    this->values = values;
    this->min = min;
    this->max = max;
73
    updateCount++;
74
}
peastman's avatar
peastman committed
75

76
Continuous1DFunction* Continuous1DFunction::Copy() const {
77
78
79
    vector<double> new_vec(values.size());
    for (size_t i = 0; i < values.size(); i++)
        new_vec[i] = values[i];
80
    return new Continuous1DFunction(new_vec, min, max);
81
82
}

83
84
85
86
87
88
89
90
91
bool Continuous1DFunction::operator==(const TabulatedFunction& other) const {
    const Continuous1DFunction* fn = dynamic_cast<const Continuous1DFunction*>(&other);
    if (fn == NULL)
        return false;
    if (fn->min != min || fn->max != max)
        return false;
    return (fn->values == values);
}

92
93
94
Continuous2DFunction::Continuous2DFunction(int xsize, int ysize, const vector<double>& values, double xmin, double xmax, double ymin, double ymax, bool periodic) {
    this->periodic = periodic;
    setFunctionParameters(xsize, ysize, values, xmin, xmax, ymin, ymax);
95
96
97
98
99
100
101
102
103
104
105
106
107
}

void Continuous2DFunction::getFunctionParameters(int& xsize, int& ysize, vector<double>& values, double& xmin, double& xmax, double& ymin, double& ymax) const {
    values = this->values;
    xsize = this->xsize;
    ysize = this->ysize;
    xmin = this->xmin;
    xmax = this->xmax;
    ymin = this->ymin;
    ymax = this->ymax;
}

void Continuous2DFunction::setFunctionParameters(int xsize, int ysize, const vector<double>& values, double xmin, double xmax, double ymin, double ymax) {
108
109
110
    if (periodic) {
        if (xsize < 3 || ysize < 3)
           throw OpenMMException("Continuous2DFunction: must have at least three points along each axis if periodic");
111
        // Note: value-matching at boundary is eventually checked at 2D-spline creation.
112
113
    }
    else if (xsize < 2 || ysize < 2)
114
115
116
117
118
119
120
121
122
123
124
125
126
127
        throw OpenMMException("Continuous2DFunction: must have at least two points along each axis");
    if (values.size() != xsize*ysize)
        throw OpenMMException("Continuous2DFunction: incorrect number of values");
    if (xmax <= xmin)
        throw OpenMMException("Continuous2DFunction: xmax <= xmin for a tabulated function.");
    if (ymax <= ymin)
        throw OpenMMException("Continuous2DFunction: ymax <= ymin for a tabulated function.");
    this->values = values;
    this->xsize = xsize;
    this->ysize = ysize;
    this->xmin = xmin;
    this->xmax = xmax;
    this->ymin = ymin;
    this->ymax = ymax;
128
    updateCount++;
129
130
}

131
Continuous2DFunction* Continuous2DFunction::Copy() const {
132
133
134
    vector<double> new_vec(values.size());
    for (size_t i = 0; i < values.size(); i++)
        new_vec[i] = values[i];
135
    return new Continuous2DFunction(xsize, ysize, new_vec, xmin, xmax, ymin, ymax);
136
137
}

138
139
140
141
142
143
144
145
146
147
148
149
150
bool Continuous2DFunction::operator==(const TabulatedFunction& other) const {
    const Continuous2DFunction* fn = dynamic_cast<const Continuous2DFunction*>(&other);
    if (fn == NULL)
        return false;
    if (fn->xsize != xsize || fn->ysize != ysize)
        return false;
    if (fn->xmin != xmin || fn->xmax != xmax)
        return false;
    if (fn->ymin != ymin || fn->ymax != ymax)
        return false;
    return (fn->values == values);
}

151
152
153
Continuous3DFunction::Continuous3DFunction(int xsize, int ysize, int zsize, const vector<double>& values, double xmin, double xmax, double ymin, double ymax, double zmin, double zmax, bool periodic) {
    this->periodic = periodic;
    setFunctionParameters(xsize, ysize, zsize, values, xmin, xmax, ymin, ymax, zmin, zmax);
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
}

void Continuous3DFunction::getFunctionParameters(int& xsize, int& ysize, int& zsize, vector<double>& values, double& xmin, double& xmax, double& ymin, double& ymax, double& zmin, double& zmax) const {
    values = this->values;
    xsize = this->xsize;
    ysize = this->ysize;
    zsize = this->zsize;
    xmin = this->xmin;
    xmax = this->xmax;
    ymin = this->ymin;
    ymax = this->ymax;
    zmin = this->zmin;
    zmax = this->zmax;
}

void Continuous3DFunction::setFunctionParameters(int xsize, int ysize, int zsize, const vector<double>& values, double xmin, double xmax, double ymin, double ymax, double zmin, double zmax) {
170
171
172
    if (periodic) {
       if (xsize < 3 || ysize < 3 || zsize < 3)
           throw OpenMMException("Continuous3DFunction: must have at least three points along each axis if periodic");
173
       // Note: value-matching at boundary is eventually checked at 3D-spline creation.
174
175
    }
    else if (xsize < 2 || ysize < 2 || zsize < 2)
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
        throw OpenMMException("Continuous3DFunction: must have at least two points along each axis");
    if (values.size() != xsize*ysize*zsize)
        throw OpenMMException("Continuous3DFunction: incorrect number of values");
    if (xmax <= xmin)
        throw OpenMMException("Continuous3DFunction: xmax <= xmin for a tabulated function.");
    if (ymax <= ymin)
        throw OpenMMException("Continuous3DFunction: ymax <= ymin for a tabulated function.");
    if (zmax <= zmin)
        throw OpenMMException("Continuous3DFunction: zmax <= zmin for a tabulated function.");
    this->values = values;
    this->xsize = xsize;
    this->ysize = ysize;
    this->zsize = zsize;
    this->xmin = xmin;
    this->xmax = xmax;
    this->ymin = ymin;
    this->ymax = ymax;
    this->zmin = zmin;
    this->zmax = zmax;
195
    updateCount++;
196
197
}

198
Continuous3DFunction* Continuous3DFunction::Copy() const {
199
200
201
    vector<double> new_vec(values.size());
    for (size_t i = 0; i < values.size(); i++)
        new_vec[i] = values[i];
202
    return new Continuous3DFunction(xsize, ysize, zsize, new_vec, xmin, xmax, ymin, ymax, zmin, zmax);
203
204
}

205
206
207
208
209
210
211
212
213
214
215
216
217
218
bool Continuous3DFunction::operator==(const TabulatedFunction& other) const {
    const Continuous3DFunction* fn = dynamic_cast<const Continuous3DFunction*>(&other);
    if (fn == NULL)
        return false;
    if (fn->xsize != xsize || fn->ysize != ysize || fn->zsize != zsize)
        return false;
    if (fn->xmin != xmin || fn->xmax != xmax)
        return false;
    if (fn->ymin != ymin || fn->ymax != ymax)
        return false;
    if (fn->zmin != zmin || fn->zmax != zmax)
        return false;
    return (fn->values == values);
}
219

220
Discrete1DFunction::Discrete1DFunction(const vector<double>& values) {
peastman's avatar
peastman committed
221
222
223
    this->values = values;
}

224
void Discrete1DFunction::getFunctionParameters(vector<double>& values) const {
peastman's avatar
peastman committed
225
226
227
    values = this->values;
}

228
void Discrete1DFunction::setFunctionParameters(const vector<double>& values) {
peastman's avatar
peastman committed
229
    this->values = values;
230
    updateCount++;
peastman's avatar
peastman committed
231
}
232

233
Discrete1DFunction* Discrete1DFunction::Copy() const {
234
235
236
    vector<double> new_vec(values.size());
    for (size_t i = 0; i < values.size(); i++)
        new_vec[i] = values[i];
237
    return new Discrete1DFunction(new_vec);
238
239
}

240
241
242
243
244
245
246
bool Discrete1DFunction::operator==(const TabulatedFunction& other) const {
    const Discrete1DFunction* fn = dynamic_cast<const Discrete1DFunction*>(&other);
    if (fn == NULL)
        return false;
    return (fn->values == values);
}

247
Discrete2DFunction::Discrete2DFunction(int xsize, int ysize, const vector<double>& values) {
248
249
250
251
252
253
254
    if (values.size() != xsize*ysize)
        throw OpenMMException("Discrete2DFunction: incorrect number of values");
    this->xsize = xsize;
    this->ysize = ysize;
    this->values = values;
}

255
void Discrete2DFunction::getFunctionParameters(int& xsize, int& ysize, vector<double>& values) const {
256
257
258
259
260
    xsize = this->xsize;
    ysize = this->ysize;
    values = this->values;
}

261
void Discrete2DFunction::setFunctionParameters(int xsize, int ysize, const vector<double>& values) {
262
263
264
265
266
    if (values.size() != xsize*ysize)
        throw OpenMMException("Discrete2DFunction: incorrect number of values");
    this->xsize = xsize;
    this->ysize = ysize;
    this->values = values;
267
    updateCount++;
268
269
}

270
Discrete2DFunction* Discrete2DFunction::Copy() const {
271
272
273
    vector<double> new_vec(values.size());
    for (size_t i = 0; i < values.size(); i++)
        new_vec[i] = values[i];
274
    return new Discrete2DFunction(xsize, ysize, new_vec);
275
276
}

277
278
279
280
281
282
283
284
285
bool Discrete2DFunction::operator==(const TabulatedFunction& other) const {
    const Discrete2DFunction* fn = dynamic_cast<const Discrete2DFunction*>(&other);
    if (fn == NULL)
        return false;
    if (fn->xsize != xsize || fn->ysize != ysize)
        return false;
    return (fn->values == values);
}

286
Discrete3DFunction::Discrete3DFunction(int xsize, int ysize, int zsize, const vector<double>& values) {
287
288
289
290
291
292
293
294
    if (values.size() != xsize*ysize*zsize)
        throw OpenMMException("Discrete3DFunction: incorrect number of values");
    this->xsize = xsize;
    this->ysize = ysize;
    this->zsize = zsize;
    this->values = values;
}

295
void Discrete3DFunction::getFunctionParameters(int& xsize, int& ysize, int& zsize, vector<double>& values) const {
296
297
298
299
300
301
    xsize = this->xsize;
    ysize = this->ysize;
    zsize = this->zsize;
    values = this->values;
}

302
void Discrete3DFunction::setFunctionParameters(int xsize, int ysize, int zsize, const vector<double>& values) {
303
304
305
306
307
308
    if (values.size() != xsize*ysize*zsize)
        throw OpenMMException("Discrete3DFunction: incorrect number of values");
    this->xsize = xsize;
    this->ysize = ysize;
    this->zsize = zsize;
    this->values = values;
309
    updateCount++;
310
}
311

312
Discrete3DFunction* Discrete3DFunction::Copy() const {
313
314
315
    vector<double> new_vec(values.size());
    for (size_t i = 0; i < values.size(); i++)
        new_vec[i] = values[i];
316
    return new Discrete3DFunction(xsize, ysize, zsize, new_vec);
317
}
318
319
320
321
322
323
324
325
326

bool Discrete3DFunction::operator==(const TabulatedFunction& other) const {
    const Discrete3DFunction* fn = dynamic_cast<const Discrete3DFunction*>(&other);
    if (fn == NULL)
        return false;
    if (fn->xsize != xsize || fn->ysize != ysize || fn->zsize != zsize)
        return false;
    return (fn->values == values);
}