-
Notifications
You must be signed in to change notification settings - Fork 0
/
inference_nobias.py
193 lines (157 loc) · 7.66 KB
/
inference_nobias.py
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
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
# third party imports
import numpy as np
import numpy.random as rnd
import matplotlib.pyplot as plt
import arviz as az
# local imports (problem definition)
from probeye.definition.inverse_problem import InverseProblem
from probeye.definition.forward_model import ForwardModelBase
from probeye.definition.distribution import Normal, Uniform
from probeye.definition.sensor import Sensor
from probeye.definition.likelihood_model import GaussianLikelihoodModel
# local imports (problem solving)
from probeye.inference.scipy.solver import MaxLikelihoodSolver
from probeye.inference.emcee.solver import EmceeSolver
from probeye.inference.dynesty.solver import DynestySolver
# local imports (inference data post-processing)
from probeye.postprocessing.sampling_plots import create_pair_plot
from probeye.postprocessing.sampling_plots import create_posterior_plot
from probeye.postprocessing.sampling_plots import create_trace_plot
from TestBridgeTemperatureDT.dummy_model import DummyModel
from TestBridgeTemperatureDT.utils import save_results_to_h5
import pandas as pd
import json
import h5py
from pathlib import Path
class DummyForwardModel(ForwardModelBase):
def __init__(self, name: str, *args, **kwargs):
super().__init__(name)
self.args = args
self.kwargs = kwargs
self.model = DummyModel(0.0, 0.0, 0.0, 0.0, 0.0, 11, 11)
self.data_path = Path(kwargs["parameters"]["data_path"])
self.figures_path = Path(kwargs["parameters"]["figures_path"])
self.file_name = kwargs["parameters"]["file_name"]
def interface(self):
self.parameters = ['T', 'WS', 'Zen', 'Az', 'Elev']
self.input_sensors = [Sensor("temperatures"), Sensor("wind_speeds"), Sensor("sun_zeniths"), Sensor("sun_azimuths"), Sensor("sun_elevations") ]
self.output_sensors = Sensor("y", std_model="sigma")
def response(self, inp: dict) -> dict:
temperatures = inp["temperatures"]
wind_speeds = inp["wind_speeds"]
sun_zeniths = inp["sun_zeniths"]
sun_azimuths = inp["sun_azimuths"]
sun_elevations = inp["sun_elevations"]
input_dict = {"temp_const":0.0, "wind_const":0.0, "zenith_constant":0.0, "azimuth_constant":0.0, "elevation_constant":0.0}
equivalence_map = {"T":"temp_const", "WS":"wind_const", "Zen":"zenith_constant", "Az":"azimuth_constant", "Elev":"elevation_constant"}
for parameter in self.parameters:
if parameter in inp.keys():
input_dict[equivalence_map[parameter]] = inp[parameter]
self.model.update_parameters(**input_dict)
self.model.calculate_grid(temperatures, wind_speeds, sun_zeniths, sun_azimuths, sun_elevations)
return {"y": self.model.grid[5,0,:].flatten()}
def inference_nobias(parameters_file, dataset_input_file, dataset_output_file, show=False, steps=1000, burn=100,
output_data_path="data/", output_figures_path="figures/", output_filename="no_bias"):
## Generate dataset for inference
# Load the HDF5 file as a dataframe
training_input = pd.read_hdf(dataset_input_file, key='data')
temperatures = training_input['Temperature'].to_numpy()
wind_speeds = training_input['Wind Speed'].to_numpy()
sun_zeniths = training_input['Sun Zenith'].to_numpy()
sun_azimuths = training_input['Sun Azimuth'].to_numpy()
sun_elevations = training_input['Sun Elevation'].to_numpy()
with h5py.File(dataset_output_file, 'r') as f:
training_output = f['Model output'][:]
with open(parameters_file) as f:
sensitivity_analysis_results = json.load(f)
parameters_dict = sensitivity_analysis_results["S1"].keys()
# Initialize forward model
parameters_forward_model = {"data_path": output_data_path, "figures_path": output_figures_path, "file_name": output_filename}
model = DummyForwardModel("Dummy Forward Model", parameters = parameters_forward_model)
# initialize the problem (the print_header=False is only set to avoid the printout of
# the probeye header which is not helpful here)
problem = InverseProblem("Dummy problem", print_header=False)
# add the problem's parameters
possible_parameters = ['T', 'WS', 'Zen', 'Az', 'Elev']
for parameter in possible_parameters:
if parameter in parameters_dict:
problem.add_parameter(
parameter,
tex=parameter,
info=f"Parameter {parameter}",
prior=Normal(mean=0.0, std=1.0),
)
else:
problem.add_parameter(
parameter,
tex=parameter,
info=f"Parameter {parameter}",
value=0.0,
)
problem.add_parameter(
"sigma",
tex=r"$\sigma$",
info="Standard deviation, of zero-mean Gaussian noise model",
value=2E-2,
)
input_dict = {"temperatures": temperatures, "wind_speeds": wind_speeds, "sun_zeniths": sun_zeniths, "sun_azimuths": sun_azimuths, "sun_elevations": sun_elevations, "y": training_output.flatten()}
# experimental data
problem.add_experiment(
name="TestSeries_1",
sensor_data=input_dict,
)
# forward model
problem.add_forward_model(model, experiments="TestSeries_1")
# likelihood model
problem.add_likelihood_model(
GaussianLikelihoodModel(experiment_name="TestSeries_1", model_error="additive")
)
# print problem summary
problem.info(print_header=True)
# this is for using the emcee-solver (MCMC sampling)
solver = EmceeSolver(problem, show_progress=True)
inference_data = solver.run(n_steps=steps, n_initial_steps=burn)
az.to_netcdf(inference_data, output_data_path + output_filename + "_inference.nc")
# save the results to HDF5
parameters_results_dict = {"parameters_mean_dict":{}, "parameters_sd_dict":{}, "parameters_median_dict":{}, "bias_mean_dict":None, "bias_sd_dict":None}
for key_parameter, value_parameter in inference_data.posterior.items():
parameters_results_dict["parameters_mean_dict"][key_parameter] = np.mean(value_parameter).values
parameters_results_dict["parameters_sd_dict"][key_parameter] = np.std(value_parameter).values
parameters_results_dict["parameters_median_dict"][key_parameter] = np.median(value_parameter)
save_results_to_h5(model, parameters_results_dict, input_dict=input_dict, sensors_flag=False)
#############################################################################
# Compare the datasets
true_values = {"T": 1.0, "WS":0.01}
# this is an overview plot that allows to visualize correlations
pair_plot_array = create_pair_plot(
inference_data,
solver.problem,
focus_on_posterior=True,
true_values = true_values,
show_legends=True,
title="Sampling results from emcee-Solver (pair plot)",
)
# figure = pair_plot_array.figure
figure = pair_plot_array[0][0].figure
figure.savefig(output_figures_path + output_filename + "_pair.png")
# this is a posterior plot, without including priors
post_plot_array = create_posterior_plot(
inference_data,
solver.problem,
# title="Simple 1D-Case",
round_to=3,
show=show
)
# figure = post_plot_array.figure
figure = post_plot_array[0].figure
figure.savefig(output_figures_path + output_filename + "_posterior.png")
# trace plots are used to check for "healthy" sampling
trace_plot_array = create_trace_plot(
inference_data,
solver.problem,
# title="Simple 1D-Case",
show=show
)
# figure = trace_plot_array.ravel()[0].figure
figure = trace_plot_array[0][0].figure
figure.savefig(output_figures_path + output_filename + "_trace.png")