1
2
3
4
5
6
7
8
9
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
38 package ffx.xray;
39
40 import static java.lang.Double.isNaN;
41 import static java.lang.String.format;
42 import static java.util.Arrays.fill;
43 import static org.apache.commons.math3.util.FastMath.abs;
44 import static org.apache.commons.math3.util.FastMath.pow;
45 import static org.apache.commons.math3.util.FastMath.sqrt;
46
47 import ffx.crystal.Crystal;
48 import ffx.crystal.HKL;
49 import ffx.crystal.ReflectionList;
50 import ffx.crystal.ReflectionSpline;
51 import ffx.numerics.OptimizationInterface;
52 import ffx.numerics.math.ComplexNumber;
53
54 import java.util.logging.Logger;
55
56
57
58
59
60
61
62
63
64 public class SplineEnergy implements OptimizationInterface {
65
66 private static final Logger logger = Logger.getLogger(SplineEnergy.class.getName());
67 private final ReflectionList reflectionList;
68 private final ReflectionSpline spline;
69 private final int nParams;
70 private final int type;
71 private final Crystal crystal;
72 private final DiffractionRefinementData refinementData;
73 private final double[][] fc;
74 private final double[][] fo;
75 private final ComplexNumber fct = new ComplexNumber();
76 private double[] optimizationScaling = null;
77 private double totalEnergy;
78
79
80
81
82
83
84
85
86
87 SplineEnergy(
88 ReflectionList reflectionList,
89 DiffractionRefinementData refinementData,
90 int nParams,
91 int type) {
92 this.reflectionList = reflectionList;
93 this.crystal = reflectionList.crystal;
94 this.refinementData = refinementData;
95 this.type = type;
96 this.fc = refinementData.fc;
97 this.fo = refinementData.fSigF;
98
99
100 this.spline = new ReflectionSpline(reflectionList, nParams);
101 this.nParams = nParams;
102 }
103
104
105
106
107 @Override
108 public boolean destroy() {
109
110 return true;
111 }
112
113
114
115
116 @Override
117 public double energy(double[] x) {
118 unscaleCoordinates(x);
119 double sum = target(x, null, false, false);
120 scaleCoordinates(x);
121 return sum;
122 }
123
124
125
126
127 @Override
128 public double energyAndGradient(double[] x, double[] g) {
129 unscaleCoordinates(x);
130 double sum = target(x, g, true, false);
131 scaleCoordinatesAndGradient(x, g);
132 return sum;
133 }
134
135
136
137
138 @Override
139 public double[] getCoordinates(double[] parameters) {
140 throw new UnsupportedOperationException("Not supported yet.");
141 }
142
143
144
145
146 @Override
147 public int getNumberOfVariables() {
148 throw new UnsupportedOperationException("Not supported yet.");
149 }
150
151
152
153
154 @Override
155 public double[] getScaling() {
156 return optimizationScaling;
157 }
158
159
160
161
162 @Override
163 public void setScaling(double[] scaling) {
164 if (scaling != null && scaling.length == nParams) {
165 optimizationScaling = scaling;
166 } else {
167 optimizationScaling = null;
168 }
169 }
170
171
172
173
174 @Override
175 public double getTotalEnergy() {
176 return totalEnergy;
177 }
178
179
180
181
182
183
184
185
186
187
188 public double target(double[] x, double[] g, boolean gradient, boolean print) {
189
190 double r = 0.0;
191 double rf = 0.0;
192 double rfree = 0.0;
193 double rfreef = 0.0;
194 double sum = 0.0;
195 double sumfo = 1.0;
196
197
198 if (gradient) {
199 fill(g, 0.0);
200 }
201
202 for (HKL ih : reflectionList.hklList) {
203 int i = ih.getIndex();
204 if (isNaN(fc[i][0]) || isNaN(fo[i][0]) || fo[i][1] <= 0.0) {
205 continue;
206 }
207
208 if (type == Type.FOTOESQ && fo[i][0] <= 0.0) {
209 continue;
210 }
211
212 double eps = ih.getEpsilon();
213 double s = crystal.invressq(ih);
214
215 double fh = spline.f(s, x);
216 refinementData.getFcTotIP(i, fct);
217
218 double d2 = 0.0;
219 double dr = 0.0;
220 double w = 0.0;
221 switch (type) {
222 case Type.FOFC:
223 w = 1.0;
224 double f1 = refinementData.getF(i);
225 double f2 = fct.abs();
226 double d = f1 - fh * f2;
227 d2 = d * d;
228 dr = -2.0 * f2 * d;
229 sumfo += f1 * f1;
230 break;
231 case Type.F1F2:
232 w = 2.0 / ih.epsilonc();
233 double ieps = 1.0 / eps;
234 f1 = pow(fct.abs(), 2.0) * ieps;
235 f2 = pow(refinementData.getF(i), 2) * ieps;
236 d = fh * f1 - f2;
237 d2 = d * d / f1;
238 dr = 2.0 * d;
239 break;
240 case Type.FCTOESQ:
241 w = 2.0 / ih.epsilonc();
242 f1 = pow(fct.abs() / sqrt(eps), 2);
243 d = f1 * fh - 1.0;
244 d2 = d * d / f1;
245 dr = 2.0 * d;
246 break;
247 case Type.FOTOESQ:
248 w = 2.0 / ih.epsilonc();
249 f1 = pow(refinementData.getF(i) / sqrt(eps), 2);
250 d = f1 * fh - 1.0;
251 d2 = d * d / f1;
252 dr = 2.0 * d;
253 break;
254 }
255
256 sum += w * d2;
257
258 double afo = abs(fo[i][0]);
259 double afh = abs(fh * fct.abs());
260 if (refinementData.isFreeR(i)) {
261 rfree += abs(afo - afh);
262 rfreef += afo;
263 } else {
264 r += abs(afo - afh);
265 rf += afo;
266 }
267
268 if (gradient) {
269 int i0 = spline.i0();
270 int i1 = spline.i1();
271 int i2 = spline.i2();
272 double g0 = spline.dfi0();
273 double g1 = spline.dfi1();
274 double g2 = spline.dfi2();
275
276 g[i0] += w * dr * g0;
277 g[i1] += w * dr * g1;
278 g[i2] += w * dr * g2;
279 }
280 }
281
282 if (gradient && type == Type.FOFC) {
283 double isumfo = 1.0 / sumfo;
284 for (int i = 0; i < g.length; i++) {
285 g[i] *= isumfo;
286 }
287 }
288
289 if (print) {
290 StringBuilder sb = new StringBuilder("\n");
291 sb.append(" Computed Potential Energy\n");
292 sb.append(format(" residual: %8.3f\n", sum / sumfo));
293 if (type == Type.FOFC || type == Type.F1F2) {
294 sb.append(
295 format(
296 " R: %8.3f Rfree: %8.3f\n", (r / rf) * 100.0, (rfree / rfreef) * 100.0));
297 }
298 sb.append("x: ");
299 for (double x1 : x) {
300 sb.append(format("%8g ", x1));
301 }
302 sb.append("\ng: ");
303 for (double v : g) {
304 sb.append(format("%8g ", v));
305 }
306 sb.append("\n");
307 logger.info(sb.toString());
308 }
309
310 totalEnergy = sum / sumfo;
311 return sum / sumfo;
312 }
313
314 public interface Type {
315
316 int FOFC = 1;
317 int F1F2 = 2;
318 int FCTOESQ = 3;
319 int FOTOESQ = 4;
320 }
321 }