1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16 package broadwick.montecarlo.markovchain;
17
18 import broadwick.BroadwickException;
19 import broadwick.io.FileOutput;
20 import broadwick.montecarlo.MonteCarlo;
21 import broadwick.montecarlo.MonteCarloResults;
22 import broadwick.montecarlo.MonteCarloScenario;
23 import broadwick.montecarlo.MonteCarloStep;
24 import broadwick.montecarlo.acceptor.MetropolisHastings;
25 import broadwick.montecarlo.acceptor.MonteCarloAcceptor;
26 import broadwick.montecarlo.markovchain.controller.MarkovChainController;
27 import broadwick.montecarlo.markovchain.controller.MarkovChainMaxNumStepController;
28 import broadwick.montecarlo.markovchain.observer.MarkovChainObserver;
29 import broadwick.rng.RNG;
30 import com.google.common.base.Throwables;
31 import java.util.HashSet;
32 import java.util.Set;
33 import lombok.Getter;
34 import lombok.Setter;
35 import lombok.extern.slf4j.Slf4j;
36
37
38
39
40 @Slf4j
41 public class MarkovChainMonteCarlo {
42
43
44
45
46
47
48
49
50 public MarkovChainMonteCarlo(final MonteCarloScenario model, final int numSimulations,
51 final MonteCarloResults consumer,
52 final MarkovStepGenerator generator) {
53 this(model, numSimulations, consumer, new MarkovChainMaxNumStepController(1000), generator);
54 }
55
56
57
58
59
60
61
62
63
64 public MarkovChainMonteCarlo(final MonteCarloScenario model, final int numSimulations,
65 final MonteCarloResults consumer, final MarkovChainController controller,
66 final MarkovStepGenerator generator) {
67 this.observers = new HashSet<>(1);
68 this.model = model;
69 this.numSimulations = numSimulations;
70 this.consumer = consumer;
71 this.mcController = controller;
72 this.pathGenerator = generator;
73 this.acceptor = new MetropolisHastings(GENERATOR.getInteger(Integer.MIN_VALUE, Integer.MAX_VALUE));
74 }
75
76
77
78
79 public final void run() {
80
81 try {
82 MonteCarloStep currentStep = pathGenerator.getInitialStep();
83
84 writer.write("# Steps taken [1]\n");
85 writer.write("# Current step accepted? [2]\n");
86 writer.write(String.format("# Current step coordinates [3-%d]%n", 2 + currentStep.getCoordinates().size()));
87 writer.write(String.format("# Results for current step [%d-n]%n", 3 + currentStep.getCoordinates().size()));
88
89 for (MarkovChainObserver observer : observers) {
90 observer.started();
91 }
92
93 log.info("Running Monte Carlo simulation with initial step {}", currentStep.toString());
94 model.setStep(currentStep);
95
96 MonteCarlo mc = new MonteCarlo(model, numSimulations);
97 mc.setResultsConsumer(consumer);
98 mc.run();
99 MonteCarloResults prevResults = mc.getResults();
100 for (MarkovChainObserver observer : observers) {
101 observer.step();
102 observer.takeMeasurements();
103 }
104 writer.write("%d,%d,%s,%s\n", numStepsTaken, 1, currentStep.toString(), prevResults.toCsv());
105 numStepsTaken++;
106
107 while (mcController.goOn(this)) {
108 int stepsSinceLastMeasurement = 0;
109 final MonteCarloStep proposedStep = pathGenerator.generateNextStep(currentStep);
110 model.setStep(proposedStep);
111
112 mc = new MonteCarlo(model, numSimulations);
113 mc.setResultsConsumer(consumer);
114 mc.run();
115 final MonteCarloResults currentResults = mc.getResults();
116
117 for (MarkovChainObserver observer : observers) {
118 observer.step();
119 }
120
121 final boolean accepted = acceptor.accept(prevResults, currentResults);
122
123 writer.write("%d,%d,%s,%s\n", numStepsTaken, accepted ? 1 : 0,
124 proposedStep.toString(), currentResults.toCsv());
125 if (accepted) {
126 log.info("Accepted Monte Carlo step ({})", proposedStep.toString());
127 numAcceptedSteps++;
128 if (numStepsTaken > burnIn
129 && (thinningInterval == 0 || stepsSinceLastMeasurement % thinningInterval == 0)) {
130 stepsSinceLastMeasurement++;
131 for (MarkovChainObserver observer : observers) {
132 observer.takeMeasurements();
133 }
134 }
135 currentStep = proposedStep;
136 prevResults = currentResults;
137 } else {
138 log.info("Rejected Monte Carlo step ({})", proposedStep.toString());
139 }
140 numStepsTaken++;
141 }
142
143 for (MarkovChainObserver observer : observers) {
144 observer.finished();
145 }
146 } catch (Exception e) {
147 log.error("Error running Monte Carlo simulation. {}", Throwables.getStackTraceAsString(e));
148 throw new BroadwickException("Error running Monte Carlo simulation." + Throwables.getStackTraceAsString(e));
149 }
150 }
151
152
153
154
155
156 public final void addObserver(final MarkovChainObserver observer) {
157 observers.add(observer);
158 }
159 @Getter
160 private int numStepsTaken = 0;
161 private int numSimulations;
162 @Getter
163 @SuppressWarnings("PMD.UnusedPrivateField")
164 private int numAcceptedSteps = 0;
165 @Setter
166 private int burnIn = 0;
167 @Setter
168 private int thinningInterval = 0;
169 @Setter
170 private MonteCarloAcceptor acceptor;
171 @Setter
172 private MarkovChainController mcController;
173 @Getter
174 private Set<MarkovChainObserver> observers;
175 @Setter
176 private FileOutput writer = new FileOutput();
177 @Getter
178 private MonteCarloScenario model;
179 private MarkovStepGenerator pathGenerator;
180 @Getter
181 private MonteCarloResults consumer;
182 private static final RNG GENERATOR = new RNG(RNG.Generator.Well19937c);
183
184 }