Skip to content

Commit

Permalink
Fix beta in simulated annealing
Browse files Browse the repository at this point in the history
  • Loading branch information
@JBjoernskov
JBjoernskov committed Jan 11, 2024
1 parent 500837c commit 63e5ed0
Show file tree
Hide file tree
Showing 4 changed files with 27 additions and 52 deletions.
2 changes: 2 additions & 0 deletions copy.sh
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
#!/bin/sh
scp -rp jabj@cei-datamining.sandbox.tek.sdu.dk:/home/jabj/Twin4Build/twin4build/estimator/tests/generated_files/model_parameters/chain_logs/ "C:/Users/jabj/OneDrive - Syddansk Universitet/PhD_Project_Jakob/Twin4build/python/BuildingEnergyModel/remote_results/"
28 changes: 8 additions & 20 deletions twin4build/estimator/estimator.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -119,6 +119,8 @@ def run_emcee_estimation(self,
assert np.all(np.abs(self.x0-self.lb)>self.tol), f"The difference between x0 and lb must be larger than {str(self.tol)}. {np.array(self.flat_attr_list)[(np.abs(self.x0-self.lb)>self.tol)==False]} violates this condition."
assert np.all(np.abs(self.x0-self.ub)>self.tol), f"The difference between x0 and ub must be larger than {str(self.tol)}. {np.array(self.flat_attr_list)[(np.abs(self.x0-self.ub)>self.tol)==False]} violates this condition."
self.use_simulated_annealing = use_simulated_annealing
if self.use_simulated_annealing:
assert n_temperature==1, "Simulated annealing can only be used if \"n_temperature\" is 1."
self.model.make_pickable()
self.model.cache(stepSize=self.stepSize,
startTime=self.startTime_train,
Expand Down Expand Up @@ -157,8 +159,7 @@ def run_emcee_estimation(self,
else:
loglike = self._loglike_wrapper

if self.use_simulated_annealing:
assert n_temperature==1, "Simulated can only be used if \"n_temperature\" is 1."


n_walkers = int(ndim*fac_walker) #*4 #Round up to nearest even number and multiply by 2

Expand Down Expand Up @@ -193,8 +194,8 @@ def run_emcee_estimation(self,
mapper=pool.imap)

burnin = 500
betas = np.linspace(0, 1, burnin)[1:]
self.beta = 0
betas = np.linspace(0, 1, burnin)[2:]
self.beta = betas[1]
chain = sampler.chain(x0_start)
n_save_checkpoint = 50 if n_sample>=50 else 1
result = {"integratedAutoCorrelatedTime": [],
Expand All @@ -219,8 +220,8 @@ def run_emcee_estimation(self,
}

for i, ensemble in tqdm(enumerate(chain.iterate(n_sample)), total=n_sample):
if i<burnin-1:
self.beta = float(betas[i])
if i<burnin-2:
self.beta = betas[i]
result["integratedAutoCorrelatedTime"].append(chain.get_acts())
result["chain.jumps_accepted"].append(chain.jumps_accepted.copy())
result["chain.jumps_proposed"].append(chain.jumps_proposed.copy())
Expand Down Expand Up @@ -313,21 +314,13 @@ def _loglike_wrapper(self, theta):
# Python returns a complex number for (-x)**y where x and y is positive. Python returns the real numbered root for -(x)**y where x and y is positive.
# Therefore abs is used to convert the negative loglike and then the sign is added after taking the power.
loglike = loglike*self.beta

return loglike

def _loglike(self, theta):
'''
This function calculates the log-likelihood. It takes in an array x representing the parameters to be optimized,
sets these parameter values in the model and simulates the model to obtain the predictions.
'''
# Set parameters for the model
# x = theta[:-n_sigma]
# sigma = theta[-n_sigma:]

outsideBounds = np.any(theta<self.lb) or np.any(theta>self.ub)
if outsideBounds: #####################################################h
return -1e+10

self.model.set_parameters_from_array(theta, self.flat_component_list, self.flat_attr_list)
self.simulator.simulate(self.model,
Expand Down Expand Up @@ -382,12 +375,7 @@ def _loglike_gaussian_process(self, theta):
This function calculates the log-likelihood. It takes in an array x representing the parameters to be optimized,
sets these parameter values in the model and simulates the model to obtain the predictions.
'''
# Set parameters for the model
# x = theta[:-n_sigma]
# sigma = theta[-n_sigma:]



theta_kernel = np.exp(theta[-self.n_par:])
theta = theta[:-self.n_par]

Expand Down Expand Up @@ -433,7 +421,7 @@ def _loglike_gaussian_process(self, theta):
print("=================")
print("")

return float(loglike)
return loglike

def uniform_logprior(self, theta):
outsideBounds = np.any(theta<self.lb) or np.any(theta>self.ub)
Expand Down
2 changes: 1 addition & 1 deletion twin4build/estimator/tests/test_estimator_wbypass.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -81,7 +81,7 @@ def test_estimator():
"fac_walker": 8, #Scaling factor for the number of ensemble walkers per chain. Minimum is 2.
"prior": "uniform", #Prior distribution - "gaussian" is also implemented
"walker_initialization": "uniform",#Initialization of parameters - "gaussian" is also implemented
"n_cores": 8,
# "n_cores": 1,
"assume_uncorrelated_noise": True,
"use_simulated_annealing": True
}
Expand Down
47 changes: 16 additions & 31 deletions twin4build/estimator/tests/test_load_emcee_chain.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,7 @@ def test_load_emcee_chain():
do_jump_plot = True
do_corner_plot = True
do_inference = False
gaussian_noise = False

# loaddir = os.path.join(uppath(os.path.abspath(__file__), 2), "chain_logs", "20230829_155706_chain_log.pickle")
# loaddir = os.path.join(uppath(os.path.abspath(__file__), 2), "chain_logs", "20230830_194210_chain_log.pickle")
Expand Down Expand Up @@ -114,10 +115,11 @@ def test_load_emcee_chain():
loaddir = os.path.join(uppath(os.path.abspath(__file__), 1), "generated_files", "model_parameters", "chain_logs", "model_20240107_224328_.pickle") #15 temps , 8*walkers, 30tau, test bypass valve, lower massflow and pressure, gaussian prior, GlycolEthanol, valve more parameters, lower UA, lower massflow, Kp
loaddir = os.path.join(uppath(os.path.abspath(__file__), 1), "generated_files", "model_parameters", "chain_logs", "model_20240107_224328_.pickle") #15 temps , 8*walkers, 30tau, test bypass valve, lower massflow and pressure, gaussian prior, GlycolEthanol, valve more parameters, lower UA, lower massflow, Kp
loaddir = os.path.join(uppath(os.path.abspath(__file__), 1), "generated_files", "model_parameters", "chain_logs", "model_20240108_175437_.pickle") #15 temps , 8*walkers, 30tau, test bypass valve, lower massflow and pressure, gaussian prior, GlycolEthanol, valve more parameters, lower UA, lower massflow, Kp
loaddir = os.path.join(uppath(os.path.abspath(__file__), 1), "generated_files", "model_parameters", "chain_logs", "model_20240109_110253_.pickle") #15 temps , 8*walkers, 30tau, test bypass valve, lower massflow and pressure, gaussian prior, GlycolEthanol, valve more parameters, lower UA, lower massflow, Kp
loaddir = os.path.join(uppath(os.path.abspath(__file__), 1), "generated_files", "model_parameters", "chain_logs", "model_20240109_143730_.pickle") #15 temps , 8*walkers, 30tau, test bypass valve, lower massflow and pressure, gaussian prior, GlycolEthanol, valve more parameters, lower UA, lower massflow, Kp
loaddir = os.path.join(uppath(os.path.abspath(__file__), 1), "generated_files", "model_parameters", "chain_logs", "model_20240110_093807_.pickle") #15 temps , 8*walkers, 30tau, test bypass valve, lower massflow and pressure, gaussian prior, GlycolEthanol, valve more parameters, lower UA, lower massflow, Kp
loaddir = os.path.join(uppath(os.path.abspath(__file__), 1), "generated_files", "model_parameters", "chain_logs", "model_20240110_141839_.pickle") #15 temps , 8*walkers, 30tau, test bypass valve, lower massflow and pressure, gaussian prior, GlycolEthanol, valve more parameters, lower UA, lower massflow, Kp
# loaddir = os.path.join(uppath(os.path.abspath(__file__), 1), "generated_files", "model_parameters", "chain_logs", "model_20240109_110253_.pickle") #15 temps , 8*walkers, 30tau, test bypass valve, lower massflow and pressure, gaussian prior, GlycolEthanol, valve more parameters, lower UA, lower massflow, Kp
# loaddir = os.path.join(uppath(os.path.abspath(__file__), 1), "generated_files", "model_parameters", "chain_logs", "model_20240109_143730_.pickle") #15 temps , 8*walkers, 30tau, test bypass valve, lower massflow and pressure, gaussian prior, GlycolEthanol, valve more parameters, lower UA, lower massflow, Kp
# loaddir = os.path.join(uppath(os.path.abspath(__file__), 1), "generated_files", "model_parameters", "chain_logs", "model_20240110_093807_.pickle") #15 temps , 8*walkers, 30tau, test bypass valve, lower massflow and pressure, gaussian prior, GlycolEthanol, valve more parameters, lower UA, lower massflow, Kp
# loaddir = os.path.join(uppath(os.path.abspath(__file__), 1), "generated_files", "model_parameters", "chain_logs", "model_20240110_141839_.pickle") #15 temps , 8*walkers, 30tau, test bypass valve, lower massflow and pressure, gaussian prior, GlycolEthanol, valve more parameters, lower UA, lower massflow, Kp
loaddir = os.path.join(uppath(os.path.abspath(__file__), 1), "generated_files", "model_parameters", "chain_logs", "model_20240110_174122_.pickle") #15 temps , 8*walkers, 30tau, test bypass valve, lower massflow and pressure, gaussian prior, GlycolEthanol, valve more parameters, lower UA, lower massflow, Kp


with open(loaddir, 'rb') as handle:
Expand Down Expand Up @@ -146,15 +148,6 @@ def test_load_emcee_chain():
nwalkers = result["chain.x"].shape[2] #Round up to nearest even number and multiply by 2


# for i in range(5):
# s = f"$a_{str(i)}$"
# s = r'{}'.format(s)
# flat_attr_list.append(s)
# for j in range(4):
# s = r'$l_{%.0f,%.0f}$' % (j,i, )
# flat_attr_list.append(s)


# startTime = datetime.datetime(year=2022, month=1, day=1, hour=0, minute=0, second=0, tzinfo=tz.gettz("Europe/Copenhagen"))
# endTime = datetime.datetime(year=2022, month=2, day=15, hour=0, minute=0, second=0, tzinfo=tz.gettz("Europe/Copenhagen"))
startTime = datetime.datetime(year=2022, month=2, day=3, hour=0, minute=0, second=0, tzinfo=tz.gettz("Europe/Copenhagen"))
Expand Down Expand Up @@ -249,16 +242,16 @@ def test_load_emcee_chain():
# ax.plot(simulator.dateTimeSteps, simulator.model.component_dict[])



for j, measuring_device in enumerate(targetMeasuringDevices):
for i in range(n_par_map[measuring_device.id]):
if i==0:
s = f"$a_{str(j)}$"
s = r'{}'.format(s)
flat_attr_list.append(s)
else:
s = r'$l_{%.0f,%.0f}$' % (j,i, )
flat_attr_list.append(s)
if gaussian_noise:
for j, measuring_device in enumerate(targetMeasuringDevices):
for i in range(n_par_map[measuring_device.id]):
if i==0:
s = f"$a_{str(j)}$"
s = r'{}'.format(s)
flat_attr_list.append(s)
else:
s = r'$l_{%.0f,%.0f}$' % (j,i, )
flat_attr_list.append(s)


assert len(flat_attr_list) == ndim, f"Number of parameters in flat_attr_list ({len(flat_attr_list)}) does not match number of parameters in chain.x ({ndim})"
Expand Down Expand Up @@ -349,14 +342,6 @@ def test_load_emcee_chain():
chain_logl[bool_] = np.nan
chain_logl[np.isnan(chain_logl)] = np.nanmin(chain_logl)

logl_min = np.min(chain_logl)
logl_max = np.max(chain_logl)
min_alpha = 0.1
max_alpha = 1
# vmin = np.min(result["chain.T"])
# vmax = np.max(result["chain.T"])
# vmin = np.min(result["chain.betas"])
# vmax = np.max(result["chain.betas"])
for nt in reversed(range(ntemps)):
for nw in range(nwalkers):
x = result["chain.x"][:, nt, nw, :]
Expand Down

0 comments on commit 63e5ed0

Please sign in to comment.