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 org.apache.commons.math3.util.FastMath.PI;
43 import static org.apache.commons.math3.util.FastMath.sqrt;
44
45 import ffx.crystal.ReflectionList;
46 import ffx.numerics.math.ComplexNumber;
47 import java.util.Random;
48 import java.util.logging.Level;
49 import java.util.logging.Logger;
50 import org.apache.commons.configuration2.CompositeConfiguration;
51
52
53
54
55
56
57
58 public class DiffractionRefinementData {
59
60 private static final Logger logger = Logger.getLogger(DiffractionRefinementData.class.getName());
61
62 public final int n;
63
64 public final int[] freeR;
65
66 public final double[][] fc;
67
68 public final double[][] fs;
69
70 public final double[][] fomPhi;
71
72 public final double[][] foFc1;
73
74 public final double fSigFCutoff;
75
76 final int nScale;
77
78 final double[][] fSigF;
79
80 final double[][] fcTot;
81
82 final double[][] foFc2;
83
84 final double[][] dFc;
85
86 final double[][] dFs;
87
88 final int nBins;
89
90 public double[] spline;
91
92 public double[] sigmaA;
93
94 public double[] sigmaW;
95
96
97
98
99 public int rFreeFlag;
100
101 double llkR, llkF;
102
103 CrystalReciprocalSpace crystalReciprocalSpaceFc;
104
105
106
107
108 CrystalReciprocalSpace crystalReciprocalSpaceFs;
109
110 double[] esqFc;
111
112 double[] esqFo;
113
114 double solventA, solventB;
115
116 double bulkSolventK, bulkSolventUeq;
117
118 double modelScaleK;
119
120 double[] modelAnisoB = new double[6];
121
122
123
124
125
126
127
128 public DiffractionRefinementData(
129 CompositeConfiguration properties, ReflectionList reflectionList) {
130
131 int rflag = -1;
132 if (properties != null) {
133 rflag = properties.getInt("rfree-flag", -1);
134 fSigFCutoff = properties.getDouble("f-sigf-cutoff", -1.0);
135 } else {
136 fSigFCutoff = -1.0;
137 }
138
139 this.n = reflectionList.hklList.size();
140 this.nScale = reflectionList.crystal.scaleN;
141 this.rFreeFlag = rflag;
142 this.nBins = reflectionList.nBins;
143 fSigF = new double[n][2];
144 freeR = new int[n];
145 fc = new double[n][2];
146 fs = new double[n][2];
147 fcTot = new double[n][2];
148 fomPhi = new double[n][2];
149 foFc2 = new double[n][2];
150 foFc1 = new double[n][2];
151 dFc = new double[n][2];
152 dFs = new double[n][2];
153
154 for (int i = 0; i < n; i++) {
155 fSigF[i][0] = fSigF[i][1] = Double.NaN;
156 fcTot[i][0] = fcTot[i][1] = Double.NaN;
157 }
158
159 spline = new double[nBins * 2];
160 sigmaA = new double[nBins];
161 sigmaW = new double[nBins];
162 esqFc = new double[nBins];
163 esqFo = new double[nBins];
164 for (int i = 0; i < nBins; i++) {
165 spline[i] = spline[i + nBins] = sigmaA[i] = esqFc[i] = esqFo[i] = 1.0;
166 sigmaW[i] = 0.05;
167 }
168
169
170 bulkSolventK = 0.33;
171 bulkSolventUeq = 50.0 / (8.0 * PI * PI);
172 modelScaleK = 0.0;
173 }
174
175
176
177
178
179
180
181 public double FoFc2F(int i) {
182 ComplexNumber c = new ComplexNumber(foFc2[i][0], foFc2[i][1]);
183 return c.abs();
184 }
185
186
187
188
189
190
191
192 public double FoFc2Phi(int i) {
193 ComplexNumber c = new ComplexNumber(foFc2[i][0], foFc2[i][1]);
194 return c.phase();
195 }
196
197
198
199
200
201
202
203 public double fcF(int i) {
204 ComplexNumber c = new ComplexNumber(fc[i][0], fc[i][1]);
205 return c.abs();
206 }
207
208
209
210
211
212
213
214 public double fcPhi(int i) {
215 ComplexNumber c = new ComplexNumber(fc[i][0], fc[i][1]);
216 return c.phase();
217 }
218
219
220
221
222
223
224
225 public double fcTotF(int i) {
226 ComplexNumber c = new ComplexNumber(fcTot[i][0], fcTot[i][1]);
227 return c.abs();
228 }
229
230
231
232
233
234
235
236 public double fcTotPhi(int i) {
237 ComplexNumber c = new ComplexNumber(fcTot[i][0], fcTot[i][1]);
238 return c.phase();
239 }
240
241
242
243
244
245
246
247 public double foFc1F(int i) {
248 ComplexNumber c = new ComplexNumber(foFc1[i][0], foFc1[i][1]);
249 return c.abs();
250 }
251
252
253
254
255
256
257
258 public double foFc1Phi(int i) {
259 ComplexNumber c = new ComplexNumber(foFc1[i][0], foFc1[i][1]);
260 return c.phase();
261 }
262
263
264
265
266
267
268
269 public double fsF(int i) {
270 ComplexNumber c = new ComplexNumber(fs[i][0], fs[i][1]);
271 return c.abs();
272 }
273
274
275
276
277
278
279
280 public double fsPhi(int i) {
281 ComplexNumber c = new ComplexNumber(fs[i][0], fs[i][1]);
282 return c.phase();
283 }
284
285
286
287
288
289
290 public void generateFsigFfromAnomalousFsigF(double[][] anomalousFsigF) {
291 for (int i = 0; i < n; i++) {
292 if (isNaN(anomalousFsigF[i][0]) && isNaN(anomalousFsigF[i][2])) {
293
294 } else if (isNaN(anomalousFsigF[i][0])) {
295 fSigF[i][0] = anomalousFsigF[i][2];
296 fSigF[i][1] = anomalousFsigF[i][3];
297 } else if (isNaN(anomalousFsigF[i][2])) {
298 fSigF[i][0] = anomalousFsigF[i][0];
299 fSigF[i][1] = anomalousFsigF[i][1];
300 } else {
301 fSigF[i][0] = (anomalousFsigF[i][0] + anomalousFsigF[i][2]) / 2.0;
302 fSigF[i][1] =
303 sqrt(
304 anomalousFsigF[i][1] * anomalousFsigF[i][1]
305 + anomalousFsigF[i][3] * anomalousFsigF[i][3]);
306 }
307 }
308 }
309
310
311 public void generateRFree() {
312 Random generator = new Random();
313 int free;
314 int nonfree;
315 int nfree = 0;
316
317 if (rFreeFlag == 0) {
318 free = 0;
319 nonfree = 1;
320 } else {
321 free = 1;
322 nonfree = 0;
323 }
324
325 for (int i = 0; i < n; i++) {
326 if (isNaN(fSigF[i][0])) {
327 freeR[i] = nonfree;
328 continue;
329 }
330
331 int randomi = generator.nextInt(100);
332 if (randomi < 5) {
333 freeR[i] = free;
334 nfree++;
335 } else {
336 freeR[i] = nonfree;
337 }
338 }
339
340 if (logger.isLoggable(Level.INFO)) {
341 logger.info(format(" Assigned %d of %d reflections to the R free set (5%%).\n", nfree, n));
342 }
343 }
344
345
346
347
348
349
350
351 public double getF(int i) {
352 return fSigF[i][0];
353 }
354
355
356
357
358
359
360
361 public double[] getFSigF(int i) {
362 return fSigF[i];
363 }
364
365
366
367
368
369
370
371 public ComplexNumber getFcTot(int i) {
372 return new ComplexNumber(fcTot[i][0], fcTot[i][1]);
373 }
374
375
376
377
378
379
380
381 public ComplexNumber getFoFc1(int i) {
382 return new ComplexNumber(foFc1[i][0], foFc1[i][1]);
383 }
384
385
386
387
388
389
390
391 public void getFoFc1IP(int i, ComplexNumber c) {
392 c.re(foFc1[i][0]);
393 c.im(foFc1[i][1]);
394 }
395
396
397
398
399
400
401
402 public ComplexNumber getFoFc2(int i) {
403 return new ComplexNumber(foFc2[i][0], foFc2[i][1]);
404 }
405
406
407
408
409
410
411
412 public void getFoFc2IP(int i, ComplexNumber c) {
413 c.re(foFc2[i][0]);
414 c.im(foFc2[i][1]);
415 }
416
417
418
419
420
421
422
423 public int getFreeR(int i) {
424 return freeR[i];
425 }
426
427
428
429
430
431
432
433 public ComplexNumber getFs(int i) {
434 return new ComplexNumber(fs[i][0], fs[i][1]);
435 }
436
437
438
439
440
441
442
443 public double getSigF(int i) {
444 return fSigF[i][1];
445 }
446
447
448
449
450
451
452 public void intensitiesToAmplitudes() {
453 double tmp;
454
455 for (int i = 0; i < n; i++) {
456 if (fSigF[i][0] > 0.0) {
457 tmp = fSigF[i][0];
458 fSigF[i][0] = sqrt(tmp);
459 if (fSigF[i][1] < tmp) {
460 fSigF[i][1] = fSigF[i][0] - sqrt(tmp - fSigF[i][1]);
461 } else {
462 fSigF[i][1] = fSigF[i][0];
463 }
464 } else if (!isNaN(fSigF[i][0])) {
465 fSigF[i][0] = 0.0;
466 fSigF[i][1] = 0.0;
467 }
468 }
469
470 if (logger.isLoggable(Level.WARNING)) {
471 StringBuilder sb = new StringBuilder();
472 sb.append("\n Internally converting intensities to amplitudes.\n");
473 sb.append(" This does not use French & Wilson scaling.\n");
474 logger.info(sb.toString());
475 }
476 }
477
478
479
480
481
482
483
484
485 public boolean isFreeR(int i, int f) {
486 return (freeR[i] == f);
487 }
488
489
490
491
492
493
494 public double likelihoodFree() {
495 return llkF;
496 }
497
498
499
500
501
502
503 public double likelihoodWork() {
504 return llkR;
505 }
506
507
508
509
510
511
512
513 public void setF(int i, double f) {
514 fSigF[i][0] = f;
515 }
516
517
518
519
520
521
522
523
524 public void setFSigF(int i, double f, double sigf) {
525 fSigF[i][0] = f;
526 fSigF[i][1] = sigf;
527 }
528
529
530
531
532
533
534
535 public void setFc(int i, ComplexNumber c) {
536 fc[i][0] = c.re();
537 fc[i][1] = c.im();
538 }
539
540
541
542
543
544
545
546 public void setFcTot(int i, ComplexNumber c) {
547 fcTot[i][0] = c.re();
548 fcTot[i][1] = c.im();
549 }
550
551
552
553
554
555
556
557 public void setFoFc1(int i, ComplexNumber c) {
558 foFc1[i][0] = c.re();
559 foFc1[i][1] = c.im();
560 }
561
562
563
564
565
566
567
568 public void setFoFc2(int i, ComplexNumber c) {
569 foFc2[i][0] = c.re();
570 foFc2[i][1] = c.im();
571 }
572
573
574
575
576
577
578
579 public void setFreeR(int i, int f) {
580 freeR[i] = f;
581 }
582
583
584
585
586
587
588 public void setFreeRFlag(int i) {
589 rFreeFlag = i;
590 }
591
592
593
594
595
596
597
598 public void setFs(int i, ComplexNumber c) {
599 fs[i][0] = c.re();
600 fs[i][1] = c.im();
601 }
602
603
604
605
606
607
608
609 public void setSigF(int i, double sigF) {
610 fSigF[i][1] = sigF;
611 }
612
613
614
615
616
617
618 void setCrystalReciprocalSpaceFc(CrystalReciprocalSpace crystalReciprocalSpace) {
619 this.crystalReciprocalSpaceFc = crystalReciprocalSpace;
620 }
621
622
623
624
625
626
627 void setCrystalReciprocalSpaceFs(CrystalReciprocalSpace crystalReciprocalSpace) {
628 this.crystalReciprocalSpaceFs = crystalReciprocalSpace;
629 }
630
631
632
633
634
635
636
637 boolean isFreeR(int i) {
638 return (freeR[i] == rFreeFlag);
639 }
640
641
642
643
644
645
646
647 ComplexNumber getFc(int i) {
648 return new ComplexNumber(fc[i][0], fc[i][1]);
649 }
650
651
652
653
654
655
656
657 void getFcIP(int i, ComplexNumber c) {
658 c.re(fc[i][0]);
659 c.im(fc[i][1]);
660 }
661
662
663
664
665
666
667
668 void getFsIP(int i, ComplexNumber c) {
669 c.re(fs[i][0]);
670 c.im(fs[i][1]);
671 }
672
673
674
675
676
677
678
679 void getFcTotIP(int i, ComplexNumber c) {
680 c.re(fcTot[i][0]);
681 c.im(fcTot[i][1]);
682 }
683 }