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.algorithms.commands.test;
39
40 import ffx.algorithms.cli.AlgorithmsCommand;
41 import ffx.numerics.math.RunningStatistics;
42 import ffx.numerics.math.SummaryStatistics;
43 import ffx.potential.ForceFieldEnergy;
44 import ffx.potential.bonded.Residue;
45 import ffx.potential.extended.ExtendedSystem;
46 import ffx.potential.parsers.XPHFilter;
47 import ffx.utilities.FFXBinding;
48 import org.apache.commons.lang3.ArrayUtils;
49 import picocli.CommandLine.Command;
50 import picocli.CommandLine.Option;
51 import picocli.CommandLine.Parameters;
52
53 import java.io.File;
54 import java.util.ArrayList;
55 import java.util.Random;
56
57 import static ffx.numerics.estimator.EstimateBootstrapper.getBootstrapIndices;
58 import static java.lang.String.format;
59
60
61
62
63
64
65
66
67 @Command(description = " Use the Rao-Blackwell estimator to get a free energy difference for residues in a CpHMD system.", name = "test.RaoBlackwellEstimator")
68 public class RaoBlackwellEstimator extends AlgorithmsCommand {
69
70 @Option(names = {"--aFi", "--arcFile"}, paramLabel = "traj",
71 description = "A file containing the the PDB from which to build the ExtendedSystem. There is currently no default.")
72 private String arcFileName = null;
73
74 @Option(names = {"--numSnaps"}, paramLabel = "-1", defaultValue = "-1",
75 description = "Number of snapshots to use from an archive file. Default is all.")
76 private int numSnaps;
77
78 @Option(names = {"--specifiedResidues", "--sR"}, paramLabel = "<selection>", defaultValue = "",
79 description = "Specified residues to do analysis.")
80 private String specified;
81
82 @Option(names = {"--startSnap"}, paramLabel = "-1", defaultValue = "-1",
83 description = "Start energy evaluations at a snap other than 2.")
84 private int startSnap;
85 @Option(names = {"--discreteStateResidues", "--dsR"}, paramLabel = "<selection>", defaultValue = "",
86 description = "Specified residues are set to discrete state when doing analysis.")
87 private String discreteStateResidues = "";
88
89 @Option(names = {"--discreteState", "--dS"}, paramLabel = "<selection>", defaultValue = "0.0",
90 description = "Specified residues are set to discrete state when doing analysis.")
91 private double discreteState = 0.0;
92
93 @Option(names = {"--chain", "--ch"}, paramLabel = " ", defaultValue = "A",
94 description = "Single character chain name (default is 'A').")
95 private Character chain = 'A';
96
97 @Option(names = {"--bootstrapIter"}, paramLabel = "100000", defaultValue = "100000",
98 description = "Number of bootstrap iterations. Set -1 for no bootstrapping.")
99 private int bootstrapIter;
100
101 @Option(names = {"--skip"}, paramLabel = "-1", defaultValue = "-1",
102 description = "Calculate energies on snaps with this interval.")
103 private int skip;
104
105 @Option(names = {"--writeFrequency"}, paramLabel = "100", defaultValue = "100",
106 description = "Calculate the RBE and print at this snapshot read frequency.")
107 private int writeFrequency;
108
109 @Parameters(arity = "1..*", paramLabel = "files",
110 description = "PDB input file in the same directory as the ARC file.")
111 private String filename;
112
113 private ForceFieldEnergy forceFieldEnergy;
114 private ArrayList<Double>[][] oneZeroDeltaLists;
115 private ArrayList<Double>[][] tautomerOneZeroDeltaList;
116 private int numESVs;
117 private int numTautomerESVs;
118
119
120
121
122 public RaoBlackwellEstimator() {
123 super();
124 }
125
126
127
128
129
130 public RaoBlackwellEstimator(FFXBinding binding) {
131 super(binding);
132 }
133
134
135
136
137
138 public RaoBlackwellEstimator(String[] args) {
139 super(args);
140 }
141
142 @Override
143 @SuppressWarnings("unchecked")
144 public RaoBlackwellEstimator run() {
145
146 if (!init()) {
147 return this;
148 }
149
150
151 File arcFile = new File(arcFileName);
152 if(!arcFile.exists()){
153 logger.severe(format(" ARC file %s does not exist.", arcFile));
154 }
155 else{
156 logger.info(format("Using ARC file %s.", arcFile));
157 }
158
159 boolean bootstrap = false;
160 if(bootstrapIter >= 50)
161 {
162 bootstrap = true;
163 } else if (bootstrapIter != -1){
164 logger.severe("Too few bootstrap iterations specified. Must be at least 50.");
165 }
166
167 activeAssembly = getActiveAssembly(filename);
168 if (activeAssembly == null) {
169 logger.info(helpString());
170 return this;
171 }
172 forceFieldEnergy = activeAssembly.getPotentialEnergy();
173
174
175 String filename = activeAssembly.getFile().getAbsolutePath();
176
177
178 ExtendedSystem esvSystem = new ExtendedSystem(activeAssembly, 7.0, null);
179
180
181 Residue specialResidue = null;
182 int numberOfStates = 1;
183 int[][] states = null;
184 if(esvSystem.getSpecialResidueList().size() > 1){
185 logger.severe(" Multiple special residues were identified in the key file. " +
186 "Only one can be specified with this algorithm.");
187 } else if (esvSystem.getSpecialResidueList().size() == 1) {
188 int specialResidueNumber = esvSystem.getSpecialResidueList().getFirst().intValue();
189 for (Residue residue : esvSystem.getTitratingResidueList()) {
190 if (residue.getResidueNumber() == specialResidueNumber && residue.getChainID() == chain) {
191 specialResidue = residue;
192 }
193 }
194 if(specialResidue != null){
195 numberOfStates = !esvSystem.isTautomer(specialResidue) ? 3 : 4;
196 switch (specialResidue.getName()) {
197
198 case "ASD":
199 case "GLD":
200 states = new int[3][2];
201 states[0][0] = 0;
202 states[0][1] = 0;
203 states[1][0] = 1;
204 states[1][1] = 0;
205 states[2][0] = 1;
206 states[2][1] = 1;
207 break;
208
209 case "HIS":
210 states = new int[3][2];
211 states[0][0] = 0;
212 states[0][1] = 0;
213 states[1][0] = 0;
214 states[1][1] = 1;
215 states[2][0] = 1;
216 states[2][1] = 0;
217 break;
218
219 case "LYS":
220 case "CYS":
221 states = new int[2][2];
222 states[0][0] = 0;
223 states[0][1] = 0;
224 states[1][0] = 1;
225 states[1][1] = 0;
226 break;
227 }
228
229 } else {
230 logger.severe(" The special residue specified in the key file was not found in the titrating residue list.");
231 }
232 }
233
234
235
236 int[] specifiedResidues = null;
237 if(specified != null && !specified.isEmpty()){
238 String[] specifiedResiduesString = specified.split(",");
239 specifiedResidues = new int[specifiedResiduesString.length];
240 for (int i = 0; i < specifiedResiduesString.length; i++) {
241 specifiedResidues[i] = Integer.parseInt(specifiedResiduesString[i].trim());
242 }
243 }
244
245 int[] dsResidueNums;
246 Residue[] dsResidues = null;
247 if(discreteStateResidues != null && !discreteStateResidues.isEmpty()){
248 String[] discreteStateResiduesString = discreteStateResidues.split(",");
249 dsResidueNums = new int[discreteStateResiduesString.length];
250 dsResidues = new Residue[dsResidueNums.length];
251 for (int i = 0; i < discreteStateResiduesString.length; i++) {
252 dsResidueNums[i] = Integer.parseInt(discreteStateResiduesString[i].trim());
253 }
254 int index=0;
255 for (Residue residue : esvSystem.getTitratingResidueList()) {
256 if (ArrayUtils.contains(dsResidueNums, residue.getResidueNumber()) && residue.getChainID() == chain) {
257 dsResidues[index] = residue;
258 index++;
259 }
260 }
261 }
262 ArrayList<Residue> onlyResidues = new ArrayList<>();
263 ArrayList<Integer> onlyResidueIndices = new ArrayList<>();
264 if(specifiedResidues != null){
265 for (int i = 0; i < esvSystem.getTitratingResidueList().size(); i++) {
266 Residue residue = esvSystem.getTitratingResidueList().get(i);
267 if (ArrayUtils.contains(specifiedResidues, residue.getResidueNumber()) && residue.getChainID() == chain) {
268 onlyResidues.add(residue);
269 onlyResidueIndices.add(i);
270 }
271 }
272 if(onlyResidues.size() != specifiedResidues.length){
273 logger.severe("Could not find all residues from --specifiedResidues input.");
274 }
275 }
276 else{
277 for (int i = 0; i < esvSystem.getTitratingResidueList().size(); i++) {
278 onlyResidueIndices.add(i);
279 }
280 }
281
282
283
284 numESVs = esvSystem.getTitratingResidueList().size();
285 oneZeroDeltaLists = new ArrayList[numESVs][numberOfStates + 1];
286 for (int i = 0; i < numESVs; i++) {
287 for (int j = 0; j < numberOfStates + 1; j++) {
288 oneZeroDeltaLists[i][j] = new ArrayList<>();
289 }
290 }
291
292 numTautomerESVs = esvSystem.getTautomerizingResidueList().size();
293 tautomerOneZeroDeltaList = new ArrayList[numTautomerESVs][numberOfStates + 1];
294 for (int i = 0; i < numTautomerESVs; i++) {
295 for (int j = 0; j < numberOfStates + 1; j++) {
296 tautomerOneZeroDeltaList[i][j] = new ArrayList<>();
297 }
298 }
299
300
301 activeAssembly.setFile(arcFile);
302 XPHFilter xphFilter = new XPHFilter(
303 arcFile,
304 activeAssembly,
305 activeAssembly.getForceField(),
306 activeAssembly.getProperties(),
307 esvSystem);
308 xphFilter.readFile();
309
310
311 esvSystem.setFixedTitrationState(true);
312 esvSystem.setFixedTautomerState(true);
313 forceFieldEnergy.attachExtendedSystem(esvSystem);
314 logger.info(format(" Attached extended system with %d residues.", numESVs));
315
316
317 double[] x = new double[forceFieldEnergy.getNumberOfVariables()];
318 forceFieldEnergy.getCoordinates(x);
319 forceFieldEnergy.energy(x, true);
320
321
322 double pH = 0.0;
323 String[] parts = xphFilter.getRemarkLines()[0].split(" ");
324 for(int i = 0; i < parts.length; i++) {
325 if (parts[i].contains("pH")) {
326 pH = Double.parseDouble(parts[i+1]);
327 }
328 }
329 logger.info("\n Setting constant pH to " + pH + ".");
330 esvSystem.setConstantPh(pH);
331
332
333 int evals = 0;
334 if(numSnaps != -1) {
335 logger.info(format(" Using %d snapshots.", numSnaps));
336 }
337 else {
338 logger.info(format(" Using all %d snapshots.", xphFilter.countNumModels()));
339 }
340
341
342 while(xphFilter.readNext()) {
343
344 if(startSnap != -1 && startSnap > 2 && evals == 0){
345 for(int i = 0; i < startSnap - 2; i++){
346 xphFilter.readNext();
347 }
348 }
349 for(Residue dsResidue: dsResidues){
350 esvSystem.setTitrationLambda(dsResidue, discreteState, false);
351 }
352
353
354 forceFieldEnergy.getCoordinates(x);
355 for (int i : onlyResidueIndices) {
356 double titrationState = 0;
357 double tautomerState = 0;
358 if(specialResidue != null){
359 titrationState = esvSystem.getTitrationLambda(specialResidue);
360 tautomerState = esvSystem.getTautomerLambda(specialResidue);
361 }
362
363
364 for (int j = 0; j < numberOfStates; j++) {
365 if (j != 0) {
366 esvSystem.setTitrationLambda(specialResidue, states[j-1][0], false);
367 if(numberOfStates != 3) {
368 esvSystem.setTautomerLambda(specialResidue, states[j-1][1], false);
369 }
370 }
371
372
373 ArrayList<Double> results = getZeroOneDeltas(i, esvSystem, forceFieldEnergy, x);
374
375 Residue res = esvSystem.getTitratingResidueList().get(i);
376 if(esvSystem.isTautomer(res)){
377 tautomerOneZeroDeltaList[esvSystem.getTautomerizingResidueList().indexOf(res)][j].add(results.get(0));
378 oneZeroDeltaLists[i][j].add(results.get(1));
379 } else{
380 oneZeroDeltaLists[i][j].add(results.getFirst());
381 }
382 }
383
384 if(specialResidue == esvSystem.getExtendedResidueList().get(i)) {
385 break;
386 }
387
388
389 if(specialResidue != null) {
390 esvSystem.setTitrationLambda(specialResidue, titrationState, false);
391 esvSystem.setTautomerLambda(specialResidue, tautomerState, false);
392 }
393 }
394
395 evals++;
396 if(evals % writeFrequency == 0 || evals == numSnaps) {
397
398 int tautomerCount = 0;
399 double[][] energyLists = new double[numESVs][numberOfStates];
400 double[][] energyStdLists = new double[numESVs][numberOfStates];
401 double[][] tautomerEnergyLists = new double[numTautomerESVs][numberOfStates];
402 double[][] tautomerEnergyStdLists = new double[numTautomerESVs][numberOfStates];
403 for (int i : onlyResidueIndices) {
404 Residue res = esvSystem.getExtendedResidueList().get(i);
405 logger.info("\n Performing Rao-Blackwell Estimator on " + res.getAminoAcid3() + ".");
406 if (bootstrap) {
407 logger.info(" Performing bootstrap with " + bootstrapIter + " iterations.");
408 } else {
409 logger.info(" Performing RBE without bootstrap. Ignore standard deviation values.");
410 }
411 for (int j = 0; j < numberOfStates; j++) {
412 double[] bootstrapMeanStd = RBE(oneZeroDeltaLists[i][j], bootstrap, bootstrapIter);
413 energyLists[i][j] = bootstrapMeanStd[0];
414 if (bootstrap) {
415 energyStdLists[i][j] = bootstrapMeanStd[1];
416 }
417
418 if (esvSystem.getTautomerizingResidueList().contains(res)) {
419 bootstrapMeanStd = RBE(tautomerOneZeroDeltaList[esvSystem.getTautomerizingResidueList().indexOf(res)][j],
420 bootstrap, bootstrapIter);
421 tautomerEnergyLists[tautomerCount][j] = bootstrapMeanStd[0];
422 if (bootstrap) {
423 tautomerEnergyStdLists[tautomerCount][j] = bootstrapMeanStd[1];
424 }
425 }
426 }
427 if (esvSystem.isTautomer(res)) {
428 tautomerCount++;
429 }
430 if (specialResidue == res) {
431 break;
432 }
433 }
434
435
436 printResults(specialResidue, esvSystem, energyLists, energyStdLists, tautomerEnergyLists, tautomerEnergyStdLists,
437 states, numberOfStates, numESVs, onlyResidueIndices);
438 }
439 if (numSnaps != -1 && evals >= numSnaps) {
440 break;
441 }
442
443 if(skip != -1){
444 for(int i = 0; i < skip-1; i++){
445 xphFilter.readNext();
446 }
447 }
448 }
449 return this;
450 }
451
452
453
454 private static double[] RBE(ArrayList<Double> deltaUList, boolean bootstrap, int bootstrapIter){
455 ArrayList<Double> deltaU = deltaUList;
456 double temperature = 298.0;
457 double boltzmann = 0.001985875;
458 double beta = 1.0 / (temperature * boltzmann);
459
460 ArrayList<Double> deltaExp = exp(mult(-beta, deltaU));
461 ArrayList<Double> numerator = div(mult(beta, mult(deltaU, deltaExp)), subtract(1.0, deltaExp));
462 ArrayList<Double> denominator = div(mult(beta, deltaU), subtract(1.0, deltaExp));
463
464 double[] deltaGRBE = bootstrap ?
465 bootStrap(numerator, denominator, bootstrapIter) :
466 new double[] {-(1.0 / beta) * Math.log(average(numerator) / average(denominator))};
467 return deltaGRBE;
468 }
469
470
471 private static double[] bootStrap(ArrayList<Double> numerator, ArrayList<Double> denominator, int iter) {
472 RunningStatistics estimates = new RunningStatistics();
473 for (int k = 0; k < iter; k++) {
474 Random rng = new Random();
475 int[] trial = getBootstrapIndices(numerator.size(), rng);
476 double estimate = estimateDg(numerator, denominator, trial);
477 estimates.addValue(estimate);
478 }
479 SummaryStatistics stats = new SummaryStatistics(estimates);
480 return new double[] {stats.mean, stats.getSd()};
481 }
482
483
484 private static double estimateDg(ArrayList<Double> num, ArrayList<Double> denom, int[] index) {
485 double temperature = 298.0;
486 double boltzmann = 0.001985875;
487 double beta = 1.0 / (temperature * boltzmann);
488 ArrayList<Double> numerator = new ArrayList<>();
489 numerator.ensureCapacity(index.length);
490 ArrayList<Double> denominator = new ArrayList<>();
491 denominator.ensureCapacity(index.length);
492
493 for (int j : index) {
494 numerator.add(num.get(j));
495 denominator.add(denom.get(j));
496 }
497
498 return -(1.0 / beta) * Math.log(average(numerator) / average(denominator));
499 }
500
501
502 private static void printResults(Residue specialResidue, ExtendedSystem esvSystem, double[][] energyLists, double[][] energyStdLists,
503 double[][] tautomerEnergyLists, double[][]tautomerStdLists, int[][] states,
504 int numberOfStates, int numESVs, ArrayList<Integer> onlyResidueIndex)
505 {
506 logger.info("\n Rao-Blackwell Estimator Results: ");
507 ArrayList<String> line = new ArrayList<>();
508 if(specialResidue != null){
509 logger.info(" Special Residue: " + specialResidue);
510 if(esvSystem.isTautomer(specialResidue)){
511 logger.info(format(" %-10s %-10s %-23s %-28s %-28s %-28s", "Residue", "Tautomer", "DeltaGTitr", "DeltaG-SpecialRes=(" + states[0][0] + "," + states[0][1] + ")", "DeltaG-SpecialRes=(" + states[1][0] + "," + states[1][1] + ")", "DeltaG-SpecialRes=(" + states[2][0] + "," + states[2][1] + ")"));
512 } else{
513 logger.info(format(" %-10s %-10s %-23s %-28s %-28s", "Residue", "Tautomer", "DeltaGTitr","DeltaG-SpecialRes=(" + states[0][0] + "," + states[0][1] + ")", "DeltaG-SpecialRes=(" + states[1][0] + "," + states[1][1] + ")"));
514 }
515 } else
516 {
517 logger.info(format(" %-10s %-10s %-23s", "Residue", "Tautomer", "DeltaGTitr"));
518 }
519 int tautomerCount = 0;
520 for(int i : onlyResidueIndex){
521 Residue res = esvSystem.getExtendedResidueList().get(i);
522 line.add(res.toString());
523 line.add("0");
524 line.add(Double.toString(energyLists[i][0]));
525 line.add(Double.toString(energyStdLists[i][0]));
526 for(int j = 1; j < numberOfStates; j++){
527 line.add(Double.toString(energyLists[i][j]));
528 line.add(Double.toString(energyStdLists[i][j]));
529 }
530 if(specialResidue != null && esvSystem.isTautomer(specialResidue)) {
531 logger.info(format(" %-10s %-10s %-10.5f +/- %-5.3f %-10.5f +/- %-5.3f %-10.5f +/- %-5.3f %-10.5f +/- %-5.3f",
532 line.get(0),
533 line.get(1),
534 Double.parseDouble(line.get(2)),
535 Double.parseDouble(line.get(3)),
536 Double.parseDouble(line.get(4)),
537 Double.parseDouble(line.get(5)),
538 Double.parseDouble(line.get(6)),
539 Double.parseDouble(line.get(7)),
540 Double.parseDouble(line.get(8)),
541 Double.parseDouble(line.get(9))));
542 } else if (specialResidue != null) {
543 logger.info(format(" %-10s %-10s %-10.5f +/- %-5.3f %-10.5f +/- %-5.3f %-10.5f +/- %-5.3f",
544 line.get(0),
545 line.get(1),
546 Double.parseDouble(line.get(2)),
547 Double.parseDouble(line.get(3)),
548 Double.parseDouble(line.get(4)),
549 Double.parseDouble(line.get(5)),
550 Double.parseDouble(line.get(6)),
551 Double.parseDouble(line.get(7))));
552 } else {
553 logger.info(format(" %-10s %-10s %-10.5f +/- %-5.3f",
554 line.get(0),
555 line.get(1),
556 Double.parseDouble(line.get(2)),
557 Double.parseDouble(line.get(3))));
558 }
559 line.clear();
560
561 if(esvSystem.isTautomer(res)) {
562 line.add(res.toString());
563 line.add("1");
564 line.add(Double.toString(tautomerEnergyLists[tautomerCount][0]));
565 line.add(Double.toString(tautomerStdLists[tautomerCount][0]));
566 for(int j = 1; j < numberOfStates; j++){
567 line.add(Double.toString(tautomerEnergyLists[tautomerCount][j]));
568 line.add(Double.toString(tautomerStdLists[tautomerCount][j]));
569 }
570 tautomerCount++;
571 if(specialResidue != null && esvSystem.isTautomer(specialResidue)) {
572 logger.info(format(" %-10s %-10s %-10.5f +/- %-5.3f %-10.5f +/- %-5.3f %-10.5f +/- %-5.3f %-10.5f +/- %-5.3f",
573 line.get(0),
574 line.get(1),
575 Double.parseDouble(line.get(2)),
576 Double.parseDouble(line.get(3)),
577 Double.parseDouble(line.get(4)),
578 Double.parseDouble(line.get(5)),
579 Double.parseDouble(line.get(6)),
580 Double.parseDouble(line.get(7)),
581 Double.parseDouble(line.get(8)),
582 Double.parseDouble(line.get(9))));
583 } else if (specialResidue != null) {
584 logger.info(format(" %-10s %-10s %-10.5f +/- %-5.3f %-10.5f +/- %-5.3f %-10.5f +/- %-5.3f",
585 line.get(0),
586 line.get(1),
587 Double.parseDouble(line.get(2)),
588 Double.parseDouble(line.get(3)),
589 Double.parseDouble(line.get(4)),
590 Double.parseDouble(line.get(5)),
591 Double.parseDouble(line.get(6)),
592 Double.parseDouble(line.get(7))));
593 } else {
594 logger.info(format(" %-10s %-10s %-10.5f +/- %-5.3f",
595 line.get(0),
596 line.get(1),
597 Double.parseDouble(line.get(2)),
598 Double.parseDouble(line.get(3))));
599 }
600 line.clear();
601 }
602 }
603 }
604
605 private static ArrayList<Double> getZeroOneDeltas(int i, ExtendedSystem esv,
606 ForceFieldEnergy forceFieldEnergy, double[] x)
607 {
608 ArrayList<Double> deltaU = new ArrayList<>();
609 Residue res = esv.getExtendedResidueList().get(i);
610 double titrationState = esv.getTitrationLambda(res);
611 double tautomerState = esv.getTautomerLambda(res);
612
613 if (esv.getTautomerizingResidueList().contains(res)) {
614 esv.setTautomerLambda(res, 1, false);
615
616 esv.setTitrationLambda(res, 0, false);
617 double zeroEnergy = forceFieldEnergy.energy(x, false);
618
619 esv.setTitrationLambda(res, 1, false);
620 double oneEnergy = forceFieldEnergy.energy(x, false);
621
622 esv.setTitrationLambda(res, titrationState, false);
623
624 deltaU.add(oneEnergy - zeroEnergy);
625 esv.setTautomerLambda(res, 0, false);
626 }
627
628 esv.setTitrationLambda(res, 0, false);
629 double zeroEnergy = forceFieldEnergy.energy(x, false);
630
631 esv.setTitrationLambda(res, 1, false);
632 double oneEnergy = forceFieldEnergy.energy(x, false);
633
634 esv.setTitrationLambda(res, titrationState, false);
635
636 deltaU.add(oneEnergy - zeroEnergy);
637 if (esv.getTautomerizingResidueList().contains(res)) {
638 esv.setTautomerLambda(res, tautomerState, false);
639 }
640
641 return deltaU;
642 }
643
644 private static double average(ArrayList<Double> list) {
645 double sum = 0.0;
646 for (Double d : list) {
647 sum += d;
648 }
649 return sum / list.size();
650 }
651
652 private static ArrayList<Double> mult(double a, ArrayList<Double> u) {
653 ArrayList<Double> result = new ArrayList<>();
654 for (Double d : u) {
655 result.add(a * d);
656 }
657 return result;
658 }
659
660 private static ArrayList<Double> mult(ArrayList<Double> v, ArrayList<Double> u) {
661 if (v.size() != u.size()) {
662 throw new IllegalArgumentException("Vector sizes must be equal.");
663 }
664 ArrayList<Double> result = new ArrayList<>();
665 for (int i = 0; i < v.size(); i++) {
666 result.add(v.get(i) * u.get(i));
667 }
668 return result;
669 }
670
671 private static ArrayList<Double> subtract(double a, ArrayList<Double> u) {
672 ArrayList<Double> result = new ArrayList<>();
673 for (Double d : u) {
674 result.add(a - d);
675 }
676 return result;
677 }
678
679 private static ArrayList<Double> exp(ArrayList<Double> u) {
680 ArrayList<Double> result = new ArrayList<>();
681 for (Double d : u) {
682 result.add(Math.exp(d));
683 }
684 return result;
685 }
686
687 private static ArrayList<Double> div(ArrayList<Double> a, ArrayList<Double> b){
688 if (a.size() != b.size()) {
689 throw new IllegalArgumentException("Vector sizes must be equal.");
690 }
691 ArrayList<Double> result = new ArrayList<>();
692 for (int i = 0; i < a.size(); i++) {
693 result.add(a.get(i) / b.get(i));
694 }
695 return result;
696 }
697 }