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 edu.rit.pj.ParallelTeam;
41 import ffx.algorithms.Terminatable;
42 import ffx.crystal.Crystal;
43 import ffx.crystal.HKL;
44 import ffx.crystal.ReflectionList;
45 import ffx.crystal.ReflectionSpline;
46 import ffx.numerics.math.ComplexNumber;
47 import ffx.numerics.optimization.LBFGS;
48 import ffx.numerics.optimization.LineSearch.LineSearchResult;
49 import ffx.numerics.optimization.OptimizationListener;
50
51 import javax.annotation.Nullable;
52 import java.util.logging.Level;
53 import java.util.logging.Logger;
54
55 import static ffx.utilities.Constants.NS2SEC;
56 import static java.lang.Double.isNaN;
57 import static java.lang.String.format;
58 import static java.lang.System.arraycopy;
59 import static org.apache.commons.math3.util.FastMath.pow;
60 import static org.apache.commons.math3.util.FastMath.sqrt;
61
62
63
64
65
66
67
68 public class SigmaAMinimize implements OptimizationListener, Terminatable {
69
70 private static final Logger logger = Logger.getLogger(SigmaAMinimize.class.getName());
71 private static final double toSeconds = 1.0e-9;
72 protected final DiffractionRefinementData refinementData;
73 private final ReflectionList reflectionList;
74 private final Crystal crystal;
75 private final SigmaAEnergy sigmaAEnergy;
76 private final int n;
77 private final double[] x;
78 private final double[] grad;
79 private final double[] scaling;
80 private boolean done = false;
81 private boolean terminate = false;
82 private long time;
83 private double grms;
84 private int nSteps;
85
86
87
88
89
90
91
92
93 SigmaAMinimize(
94 ReflectionList reflectionList,
95 DiffractionRefinementData refinementData,
96 ParallelTeam parallelTeam) {
97 this.reflectionList = reflectionList;
98 this.refinementData = refinementData;
99 this.crystal = reflectionList.crystal;
100
101 n = refinementData.nBins * 2;
102 sigmaAEnergy = new SigmaAEnergy(reflectionList, refinementData, parallelTeam);
103 x = new double[n];
104 grad = new double[n];
105 scaling = new double[n];
106
107
108 SplineEnergy.SplineType type = SplineEnergy.SplineType.FCTOESQ;
109 SplineMinimize splineMinimize =
110 new SplineMinimize(reflectionList, refinementData, refinementData.esqFc, type);
111 splineMinimize.minimize(7, 1.0e-1);
112
113
114 type = SplineEnergy.SplineType.FOTOESQ;
115 splineMinimize = new SplineMinimize(reflectionList, refinementData, refinementData.esqFo, type);
116 splineMinimize.minimize(7, 1.0e-1);
117
118 int offset = refinementData.nBins;
119 for (int i = 0; i < refinementData.nBins; i++) {
120
121 x[i] = 1.0;
122 scaling[i] = 1.0;
123
124 x[i + offset] = 0.05;
125 scaling[i + offset] = 20.0;
126 }
127 sigmaAEnergy.setScaling(scaling);
128 sigmaAEnergy.scaleCoordinates(x);
129
130 logger.fine(" S=1.0, W=0.05");
131 logger.fine(format(" Likelihood Target: %16.8f", calculateLikelihood()));
132 logger.fine(format(" Likelihood Target (Free): %16.8f", calculateLikelihoodFree()));
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149 }
150
151 private void logSigmA() {
152 int offset = refinementData.nBins;
153 StringBuilder sigA = new StringBuilder(" Sigma A: ");
154 StringBuilder sigAScale = new StringBuilder(" Scale A: ");
155 StringBuilder sigW = new StringBuilder(" Sigma W: ");
156 StringBuilder sigWScale = new StringBuilder(" Scale W: ");
157 for (int i = 0; i < refinementData.nBins; i++) {
158 sigA.append(format("%5.3f ", x[i]));
159 sigW.append(format("%5.3f ", x[i + offset]));
160 sigAScale.append(format("%5.3f ", scaling[i]));
161 sigWScale.append(format("%5.3f ", scaling[i + offset]));
162 }
163 logger.info(sigA.toString());
164 logger.info(sigAScale.toString());
165 logger.info(sigW.toString());
166 logger.info(sigWScale.toString());
167 }
168
169
170
171
172
173
174
175 public double[] getCoordinates(@Nullable double[] x) {
176 if (x == null) {
177 x = new double[this.x.length];
178 }
179 arraycopy(this.x, 0, x, 0, this.x.length);
180 return x;
181 }
182
183
184
185
186
187
188 public int getNumberOfVariables() {
189 return x.length;
190 }
191
192
193
194
195
196
197 public SigmaAEnergy minimize() {
198 return minimize(0.5);
199 }
200
201
202
203
204
205
206
207 public SigmaAEnergy minimize(double eps) {
208 return minimize(7, eps);
209 }
210
211
212
213
214
215
216
217
218 public SigmaAEnergy minimize(int m, double eps) {
219
220 sigmaAEnergy.setScaling(scaling);
221
222 double e = sigmaAEnergy.energyAndGradient(x, grad);
223
224 long mTime = -System.nanoTime();
225 time = -System.nanoTime();
226 done = false;
227 int status = LBFGS.minimize(n, m, x, e, grad, eps, sigmaAEnergy, this);
228 done = true;
229
230 switch (status) {
231 case 0:
232 logger.info(format("\n Optimization achieved convergence criteria: %8.5f\n", grms));
233 break;
234 case 1:
235 logger.info(format("\n Optimization terminated at step %d.\n", nSteps));
236 break;
237 default:
238 logger.warning("\n Optimization failed.\n");
239 }
240
241 sigmaAEnergy.unscaleCoordinates(x);
242 int offset = refinementData.nBins;
243 for (int i = 0; i < refinementData.nBins; i++) {
244 refinementData.sigmaA[i] = x[i];
245 refinementData.sigmaW[i] = x[i + offset];
246 }
247
248 if (logger.isLoggable(Level.INFO)) {
249 StringBuilder sb = new StringBuilder();
250 mTime += System.nanoTime();
251 sb.append(format(" Optimization time: %8.3f (sec)\n", mTime * toSeconds));
252 logger.info(sb.toString());
253 }
254
255 if (logger.isLoggable(Level.FINE)) {
256 logger.fine(" Final Result: ");
257 logSigmA();
258 sigmaAEnergy.setScaling(scaling);
259 sigmaAEnergy.scaleCoordinates(x);
260 logger.fine(format(" Likelihood Target: %16.8f", calculateLikelihood()));
261 logger.fine(format(" Likelihood Target (Free): %16.8f", calculateLikelihoodFree()));
262 sigmaAEnergy.setScaling(null);
263 }
264
265 return sigmaAEnergy;
266 }
267
268
269
270
271 @Override
272 public boolean optimizationUpdate(int iter, int nBFGS, int nfun, double grms, double xrms,
273 double f, double df, double angle, LineSearchResult info) {
274 long currentTime = System.nanoTime();
275 Double seconds = (currentTime - time) * NS2SEC;
276 time = currentTime;
277 this.grms = grms;
278 this.nSteps = iter;
279
280 if (iter == 0) {
281 if (nBFGS > 0) {
282 logger.info("\n Limited Memory BFGS Quasi-Newton Optimization of SigmaA Parameters\n");
283 } else {
284 logger.info("\n Steepest Decent Optimization of SigmaA Parameters\n");
285 }
286 logger.info(" Cycle Energy G RMS Delta E Delta X Angle Evals Time");
287 }
288 if (info == null) {
289 logger.info(format("%6d %12.2f %10.2f", iter, f, grms));
290 } else {
291 if (info == LineSearchResult.Success) {
292 logger.info(format("%6d %12.2f %10.2f %10.5f %9.5f %8.2f %6d %8.3f",
293 iter, f, grms, df, xrms, angle, nfun, seconds));
294 } else {
295 logger.info(format("%6d %12.2f %10.2f %10.5f %9.5f %8.2f %6d %8s",
296 iter, f, grms, df, xrms, angle, nfun, info));
297 }
298 }
299 if (terminate) {
300 logger.info(" The optimization recieved a termination request.");
301
302 return false;
303 }
304 return true;
305 }
306
307
308
309
310 @Override
311 public void terminate() {
312 terminate = true;
313 while (!done) {
314 synchronized (this) {
315 try {
316 wait(1);
317 } catch (Exception e) {
318 logger.log(Level.WARNING, "Exception terminating minimization.\n", e);
319 }
320 }
321 }
322 }
323
324 private void setWEstimate() {
325
326 ReflectionSpline spline = new ReflectionSpline(reflectionList, refinementData.nBins);
327 int[] nMean = new int[refinementData.nBins];
328 for (int i = 0; i < refinementData.nBins; i++) {
329 nMean[i] = 0;
330 }
331 double mean = 0.0;
332 double tot = 0.0;
333 double[][] fcTot = refinementData.fcTot;
334 double[][] fSigF = refinementData.fSigF;
335 for (HKL ih : reflectionList.hklList) {
336 int i = ih.getIndex();
337 if (ih.getAllowed() == 0.0 || isNaN(fcTot[i][0]) || isNaN(fSigF[i][0])) {
338 continue;
339 }
340
341 double s2 = crystal.invressq(ih);
342 double epsc = ih.epsilonc();
343 ComplexNumber fct = new ComplexNumber(fcTot[i][0], fcTot[i][1]);
344 double ecscale = spline.f(s2, refinementData.esqFc);
345 double eoscale = spline.f(s2, refinementData.esqFo);
346 double ec = fct.times(sqrt(ecscale)).abs();
347 double eo = fSigF[i][0] * sqrt(eoscale);
348 double wi = pow(eo - ec, 2.0) / epsc;
349
350 nMean[spline.i1()]++;
351 tot++;
352
353 x[spline.i1() + refinementData.nBins] +=
354 (wi - x[spline.i1() + refinementData.nBins]) / nMean[spline.i1()];
355
356 mean += (wi - mean) / tot;
357 }
358 logger.fine(format(" Starting mean w: %8.3f", mean));
359 double initialScale = 0.01;
360 if (mean > 0.0) {
361 initialScale = 1.0 / mean;
362 }
363 logger.fine(format(" Starting w scaling: %8.3f", initialScale));
364 for (int i = 0; i < refinementData.nBins; i++) {
365 x[i] -= x[i + refinementData.nBins];
366 x[i] *= scaling[i];
367 scaling[i + refinementData.nBins] = initialScale;
368 x[i + refinementData.nBins] *= scaling[i + refinementData.nBins];
369 }
370 }
371
372
373
374
375
376
377 SigmaAEnergy getSigmaAEnergy() {
378 return sigmaAEnergy;
379 }
380
381
382
383
384
385
386 double calculateLikelihood() {
387 sigmaAEnergy.setScaling(scaling);
388 sigmaAEnergy.energyAndGradient(x, grad);
389 sigmaAEnergy.setScaling(null);
390 return refinementData.llkR;
391 }
392
393
394
395
396
397
398 public double calculateLikelihoodFree() {
399 sigmaAEnergy.setScaling(scaling);
400 double energy = sigmaAEnergy.energyAndGradient(x, grad);
401 sigmaAEnergy.setScaling(null);
402 return energy;
403 }
404 }