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.lang.System.arraycopy;
43 import static java.util.Arrays.fill;
44 import static org.apache.commons.math3.util.FastMath.log;
45
46 import edu.rit.pj.ParallelTeam;
47 import ffx.algorithms.Terminatable;
48 import ffx.crystal.Crystal;
49 import ffx.crystal.HKL;
50 import ffx.crystal.ReflectionList;
51 import ffx.numerics.optimization.LBFGS;
52 import ffx.numerics.optimization.LineSearch.LineSearchResult;
53 import ffx.numerics.optimization.OptimizationListener;
54 import ffx.xray.CrystalReciprocalSpace.SolventModel;
55
56 import javax.annotation.Nullable;
57 import java.util.logging.Level;
58 import java.util.logging.Logger;
59
60
61
62
63
64
65
66 public class ScaleBulkMinimize implements OptimizationListener, Terminatable {
67
68 private static final Logger logger = Logger.getLogger(ScaleBulkMinimize.class.getName());
69 private static final double toSeconds = 1.0e-9;
70 private final ReflectionList reflectionlist;
71 private final DiffractionRefinementData refinementData;
72 private final Crystal crystal;
73 private final CrystalReciprocalSpace crystalReciprocalSpace;
74 private final ScaleBulkEnergy bulkSolventEnergy;
75 private final int solventN;
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
94 public ScaleBulkMinimize(
95 ReflectionList reflectionList,
96 DiffractionRefinementData refinementData,
97 CrystalReciprocalSpace crystalReciprocalSpace,
98 ParallelTeam parallelTeam) {
99 this.reflectionlist = reflectionList;
100 this.refinementData = refinementData;
101 this.crystal = reflectionList.crystal;
102 this.crystalReciprocalSpace = crystalReciprocalSpace;
103
104 if (crystalReciprocalSpace.solventModel == SolventModel.NONE) {
105 solventN = 1;
106 } else {
107 solventN = 3;
108 }
109 n = solventN + refinementData.nScale;
110 bulkSolventEnergy = new ScaleBulkEnergy(reflectionList, refinementData, n, parallelTeam);
111
112 x = new double[n];
113 grad = new double[n];
114 scaling = new double[n];
115 fill(scaling, 1.0);
116
117 x[0] = refinementData.modelScaleK;
118 if (solventN > 1) {
119 x[1] = refinementData.bulkSolventK;
120 x[2] = refinementData.bulkSolventUeq;
121 }
122 for (int i = 0; i < 6; i++) {
123 if (crystal.scaleB[i] >= 0) {
124 x[solventN + crystal.scaleB[i]] = refinementData.modelAnisoB[i];
125 }
126 }
127
128 setInitialScale();
129 }
130
131
132
133
134
135
136
137 public double[] getCoordinates(@Nullable double[] x) {
138 if (x == null) {
139 x = new double[this.x.length];
140 }
141 arraycopy(this.x, 0, x, 0, this.x.length);
142 return x;
143 }
144
145
146
147
148
149
150 public int getNumberOfVariables() {
151 return x.length;
152 }
153
154
155 public void ksbsGridOptimize() {
156 if (solventN < 3) {
157 return;
158 }
159
160 bulkSolventEnergy.setScaling(scaling);
161
162 double min = Double.POSITIVE_INFINITY;
163 double bulkSolventK = refinementData.bulkSolventK;
164 double bulkSolventUeq = refinementData.bulkSolventUeq;
165 double kmin = 0.05;
166 double kmax = 0.9;
167 double kstep = 0.05;
168 double bmin = 10.0;
169 double bmax = 200.0;
170 double bstep = 10.0;
171 for (double i = kmin; i <= kmax; i += kstep) {
172 for (double j = bmin; j <= bmax; j += bstep) {
173 x[1] = i;
174 x[2] = j;
175 double sum = bulkSolventEnergy.energyAndGradient(x, grad);
176 logger.info(" ks: " + i + " bs: " + j + " sum: " + sum);
177 if (sum < min) {
178 min = sum;
179 bulkSolventK = i;
180 bulkSolventUeq = j;
181 }
182 }
183 }
184 logger.info(" minks: " + bulkSolventK + " minbs: " + bulkSolventUeq + " min: " + min);
185 refinementData.bulkSolventK = bulkSolventK;
186 refinementData.bulkSolventUeq = bulkSolventUeq;
187 bulkSolventEnergy.setScaling(null);
188 }
189
190
191
192
193
194
195 public ScaleBulkEnergy minimize() {
196 return minimize(0.5);
197 }
198
199
200
201
202
203
204
205 public ScaleBulkEnergy minimize(double eps) {
206 return minimize(5, eps);
207 }
208
209
210
211
212
213
214
215
216 public ScaleBulkEnergy minimize(int m, double eps) {
217 bulkSolventEnergy.setScaling(scaling);
218 double e = bulkSolventEnergy.energyAndGradient(x, grad);
219
220 long mtime = -System.nanoTime();
221 time = -System.nanoTime();
222 done = false;
223 int status = LBFGS.minimize(n, m, x, e, grad, eps, bulkSolventEnergy, this);
224 done = true;
225 switch (status) {
226 case 0:
227 logger.info(format("\n Optimization achieved convergence criteria: %10.5f\n", grms));
228 break;
229 case 1:
230 logger.info(format("\n Optimization terminated at step %d.\n", nSteps));
231 break;
232 default:
233 logger.warning("\n Optimization failed.\n");
234 }
235 refinementData.modelScaleK = x[0] / scaling[0];
236 if (solventN > 1) {
237 refinementData.bulkSolventK = x[1] / scaling[1];
238 refinementData.bulkSolventUeq = x[2] / scaling[2];
239 if (crystalReciprocalSpace != null) {
240 refinementData.solventA = crystalReciprocalSpace.solventA;
241 refinementData.solventB = crystalReciprocalSpace.solventB;
242 }
243 }
244 for (int i = 0; i < 6; i++) {
245 int offset = crystal.scaleB[i];
246 if (offset >= 0) {
247 int index = solventN + offset;
248 refinementData.modelAnisoB[i] = x[index] / scaling[index];
249 }
250 }
251
252 if (logger.isLoggable(Level.INFO)) {
253 mtime += System.nanoTime();
254 logger.info(format(" Optimization time: %8.3f (sec)\n", mtime * toSeconds));
255 }
256
257 bulkSolventEnergy.setScaling(null);
258
259 return bulkSolventEnergy;
260 }
261
262
263 @Override
264 public boolean optimizationUpdate(
265 int iter,
266 int nBFGS,
267 int nfun,
268 double grms,
269 double xrms,
270 double f,
271 double df,
272 double angle,
273 LineSearchResult info) {
274 long currentTime = System.nanoTime();
275 Double seconds = (currentTime - time) * 1.0e-9;
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 Fc to Fo Scaling Parameters\n");
283 } else {
284 logger.info("\n Steepest Decent Optimization of Fc to Fo Scaling 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.5f %10.6f", iter, f, grms));
290 } else {
291 if (info == LineSearchResult.Success) {
292 logger.info(format("%6d %12.5f %10.6f %10.6f %9.5f %8.2f %6d %8.3f",
293 iter, f, grms, df, xrms, angle, nfun, seconds));
294 } else {
295 logger.info(format("%6d %12.5f %10.6f %10.6f %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 @Override
309 public void terminate() {
310 terminate = true;
311 while (!done) {
312 synchronized (this) {
313 try {
314 wait(1);
315 } catch (Exception e) {
316 logger.log(Level.WARNING, "Exception terminating minimization.\n", e);
317 }
318 }
319 }
320 }
321
322
323
324
325
326
327 ScaleBulkEnergy getScaleBulkEnergy() {
328 return bulkSolventEnergy;
329 }
330
331 private void setInitialScale() {
332 double[][] fc = refinementData.fc;
333 double[][] fSigf = refinementData.fSigF;
334
335 bulkSolventEnergy.setScaling(scaling);
336 double e = bulkSolventEnergy.energyAndGradient(x, grad);
337 bulkSolventEnergy.setScaling(null);
338
339 double sumfofc = 0.0;
340 double sumfc = 0.0;
341 for (HKL ih : reflectionlist.hklList) {
342 int i = ih.getIndex();
343 if (isNaN(fc[i][0]) || isNaN(fSigf[i][0]) || fSigf[i][1] <= 0.0) {
344 continue;
345 }
346
347 double fcTotF = refinementData.fcTotF(i);
348 sumfofc += fSigf[i][0] * fcTotF;
349 sumfc += fcTotF * fcTotF;
350 }
351
352 x[0] = log(4.0 * sumfofc / sumfc);
353 }
354
355
356 void gridOptimize() {
357 if (crystalReciprocalSpace == null) {
358 return;
359 }
360
361 bulkSolventEnergy.setScaling(scaling);
362
363 double min = Double.POSITIVE_INFINITY;
364 double solventA = crystalReciprocalSpace.solventA;
365 double solventB = crystalReciprocalSpace.solventB;
366 double amin = solventA - 1.0;
367 double amax = (solventA + 1.0) / 0.9999;
368 double astep = 0.25;
369 double bmin = solventB - 0.2;
370 double bmax = (solventB + 0.2) / 0.9999;
371 double bstep = 0.05;
372 if (crystalReciprocalSpace.solventModel == SolventModel.BINARY) {
373 amin = solventA - 0.2;
374 amax = (solventA + 0.2) / 0.9999;
375 astep = 0.05;
376 }
377
378 logger.info(" Bulk Solvent Grid Search");
379 for (double i = amin; i <= amax; i += astep) {
380 for (double j = bmin; j <= bmax; j += bstep) {
381 crystalReciprocalSpace.setSolventAB(i, j);
382 crystalReciprocalSpace.computeDensity(refinementData.fs);
383 double sum = bulkSolventEnergy.energy(x);
384 logger.info(format(" A: %6.3f B: %6.3f Sum: %12.8f", i, j, sum));
385 if (sum < min) {
386 min = sum;
387 solventA = i;
388 solventB = j;
389 }
390 }
391 }
392
393 logger.info(format("\n Minimum at\n A: %6.3f B: %6.3f Sum: %12.8f", solventA, solventB, min));
394 crystalReciprocalSpace.setSolventAB(solventA, solventB);
395 refinementData.solventA = solventA;
396 refinementData.solventB = solventB;
397 crystalReciprocalSpace.computeDensity(refinementData.fs);
398 bulkSolventEnergy.setScaling(null);
399 }
400 }