Commit da012a02 authored by Pablo Iturrieta Rebolledo's avatar Pablo Iturrieta Rebolledo
Browse files

Fixed error with numyp seed

parent b6f33441
...@@ -2,7 +2,7 @@ ...@@ -2,7 +2,7 @@
<module type="PYTHON_MODULE" version="4"> <module type="PYTHON_MODULE" version="4">
<component name="NewModuleRootManager"> <component name="NewModuleRootManager">
<content url="file://$MODULE_DIR$" /> <content url="file://$MODULE_DIR$" />
<orderEntry type="jdk" jdkName="Python 3.7 (csep2)" jdkType="Python SDK" /> <orderEntry type="jdk" jdkName="Python 3.8 (venv)" jdkType="Python SDK" />
<orderEntry type="sourceFolder" forTests="false" /> <orderEntry type="sourceFolder" forTests="false" />
</component> </component>
<component name="TestRunnerService"> <component name="TestRunnerService">
......
<?xml version="1.0" encoding="UTF-8"?> <?xml version="1.0" encoding="UTF-8"?>
<project version="4"> <project version="4">
<component name="ProjectRootManager" version="2" project-jdk-name="Python 3.7 (csep2)" project-jdk-type="Python SDK" /> <component name="ProjectRootManager" version="2" project-jdk-name="Python 3.8 (venv)" project-jdk-type="Python SDK" />
</project> </project>
\ No newline at end of file
...@@ -26,14 +26,12 @@ import pandas as pd ...@@ -26,14 +26,12 @@ import pandas as pd
def get_simulations(model, catalog, lons, lats): def get_simulations(model, catalog, lons, lats):
seed = 23
forecast_data = model.spatial_counts() forecast_data = model.spatial_counts()
n_obs = catalog.get_number_of_events() n_obs = catalog.get_number_of_events()
num_events_to_simulate = n_obs num_events_to_simulate = n_obs
sampling_weights = numpy.cumsum(forecast_data) / numpy.sum(forecast_data) sampling_weights = numpy.cumsum(forecast_data) / numpy.sum(forecast_data)
sim_fore = numpy.zeros(sampling_weights.shape) sim_fore = numpy.zeros(sampling_weights.shape)
numpy.random.seed(seed)
sim_func = csep.core.poisson_evaluations._simulate_catalog sim_func = csep.core.poisson_evaluations._simulate_catalog
sim = sim_func(num_events_to_simulate, sampling_weights, sim_fore).astype(int) sim = sim_func(num_events_to_simulate, sampling_weights, sim_fore).astype(int)
ind_sim = np.where(sim != 0)[0] ind_sim = np.where(sim != 0)[0]
...@@ -279,18 +277,18 @@ def plot_results(Results, alpha=0.05): ...@@ -279,18 +277,18 @@ def plot_results(Results, alpha=0.05):
(L_cat<L_down)*(L_cat-L_down) (L_cat<L_down)*(L_cat-L_down)
L_over[np.isnan(L_over)] = 0. L_over[np.isnan(L_over)] = 0.
L_diff[key] = L_over L_diff[key] = L_over
# for i in L_sims: for i in L_sims:
# sns.lineplot(x=r_sims, y=i, lw=0.05, color='blue', ls='--') sns.lineplot(x=r_sims, y=i, lw=0.05, color='blue', ls='--')
# sns.lineplot(x=r_cat, y=L_cat, color='r', label='Observed catalog') sns.lineplot(x=r_cat, y=L_cat, color='r', label='Observed catalog')
# sns.lineplot(x=r_sims, y=L_avg, label='Sim. average') sns.lineplot(x=r_sims, y=L_avg, label='Sim. average')
# plt.fill_between(r_sims, L_down, L_up, color='gray', alpha=0.4, label=r'Sim. envelope ($\alpha=%.2f$)' % (1-alpha)) plt.fill_between(r_sims, L_down, L_up, color='gray', alpha=0.4, label=r'Sim. envelope ($\alpha=%.2f$)' % (1-alpha))
# plt.title("Model - %s" % key, fontsize=16) plt.title("Model - %s" % key, fontsize=16)
# plt.xlabel(r"$r~~\mathrm{[km]}$") plt.xlabel(r"$r~~\mathrm{[km]}$")
# plt.ylabel(r"$L(r) = \sqrt{\frac{\hat{K}(r)}{\pi}}$") plt.ylabel(r"$L(r) = \sqrt{\frac{\hat{K}(r)}{\pi}}$")
# plt.legend(loc=4) plt.legend(loc=4)
# plt.savefig(paths.get_ripleyk_figpath('L', 10, key)) plt.savefig(paths.get_kripley_figpath('L', 10, key))
# plt.show() plt.show()
#
# for i in L_sims: # for i in L_sims:
# sns.lineplot(x=r_sims, y=i-r_sims, lw=0.05, color='blue', ls='--') # sns.lineplot(x=r_sims, y=i-r_sims, lw=0.05, color='blue', ls='--')
# sns.lineplot(x=r_cat, y=L_cat-r_sims, color='r', label='Observed catalog') # sns.lineplot(x=r_cat, y=L_cat-r_sims, color='r', label='Observed catalog')
...@@ -300,7 +298,7 @@ def plot_results(Results, alpha=0.05): ...@@ -300,7 +298,7 @@ def plot_results(Results, alpha=0.05):
# plt.xlabel(r"$r~~\mathrm{[km]}$") # plt.xlabel(r"$r~~\mathrm{[km]}$")
# plt.ylabel(r"$L(r)-r = \sqrt{\frac{\hat{K}(r)}{\pi}}-r$") # plt.ylabel(r"$L(r)-r = \sqrt{\frac{\hat{K}(r)}{\pi}}-r$")
# plt.legend(loc=4) # plt.legend(loc=4)
# plt.savefig(paths.get_ripleyk_figpath('Lr', 10, key)) # plt.savefig(paths.get_kripley_figpath('Lr', 10, key))
# plt.show() # plt.show()
# #
pcf_sims = value['PCF_sims'] pcf_sims = value['PCF_sims']
...@@ -328,8 +326,8 @@ def plot_results(Results, alpha=0.05): ...@@ -328,8 +326,8 @@ def plot_results(Results, alpha=0.05):
# plt.show() # plt.show()
# env_diff['r'] = r_sims # env_diff['r'] = r_sims
b = plot_combined(L_diff, pcf_diff, r_sims, pcfr_sims) # b = plot_combined(L_diff, pcf_diff, r_sims, pcfr_sims)
return b # return b
plt.show() plt.show()
...@@ -440,7 +438,7 @@ def plot_combined(A,B, x, xx, order='inc'): ...@@ -440,7 +438,7 @@ def plot_combined(A,B, x, xx, order='inc'):
plt.show() plt.show()
return b return b
def run_parallel(): def run():
models = get_models.ten_years() models = get_models.ten_years()
catalogs = [get_catalogs.filter_10yr_01_2010(get_catalogs.bollettino()), catalogs = [get_catalogs.filter_10yr_01_2010(get_catalogs.bollettino()),
...@@ -495,14 +493,14 @@ def run_parallel(): ...@@ -495,14 +493,14 @@ def run_parallel():
assert numpy.allclose(np.mean(Rk_sims, axis=0), np.array(Rk_sims[0])) assert numpy.allclose(np.mean(Rk_sims, axis=0), np.array(Rk_sims[0]))
assert numpy.allclose(np.mean(Rpcf_sims, axis=0), np.array(Rpcf_sims[0])) assert numpy.allclose(np.mean(Rpcf_sims, axis=0), np.array(Rpcf_sims[0]))
L = np.divide(np.sqrt(K_sims),np.sqrt(np.pi)) L_sims = np.divide(np.sqrt(K_sims),np.sqrt(np.pi))
L_cat = np.sqrt(np.array(K_cat) / np.pi) L_cat = np.sqrt(np.array(K_cat) / np.pi)
Results[model.name] = {'K_sims': K, Results[model.name] = {'K_sims': K_sims,
'Rk_sims': Rk_sims, 'Rk_sims': Rk_sims,
'K_cat': K_cat, 'K_cat': K_cat,
'rk_cat': rk_cat, 'rk_cat': rk_cat,
'L_sims': L, 'L_sims': L_sims,
'L_cat':L_cat, 'L_cat':L_cat,
'PCF_sims': PCF_sims, 'PCF_sims': PCF_sims,
'Rpcf_sims':Rpcf_sims, 'Rpcf_sims':Rpcf_sims,
...@@ -515,78 +513,12 @@ def run_parallel(): ...@@ -515,78 +513,12 @@ def run_parallel():
return Results return Results
def run():
models = get_models.ten_years()
catalog = get_catalogs.filter_10yr_01_2010(get_catalogs.bollettino())
Results = dict.fromkeys([i.name for i in models])
for model in models:
print(model.name)
points = np.stack((catalog.get_longitudes(), catalog.get_latitudes())).T
r = np.linspace(0,800,100)
K = K_r_spatstat(points, model, r=r)
rk_cat = K[0]
K_cat = K[2]
PCF = pcf_spatstat(K)
rpcf_cat = PCF[0]
pcf_cat = PCF[2]
lons = model.get_longitudes()
lats = model.get_latitudes()
nsim = 200
sim_catalogs = iter([get_simulations(model, catalog, lons, lats) for i in range(nsim)])
Rk_sims = []
K_sims = []
Rpcf_sims = []
PCF_sims = []
start = time.process_time()
for cat, inds in sim_catalogs:
K = K_r_spatstat(cat, model, r=r)
rk_i = list(K[0])
K_i = list(K[2])
PCF = pcf_spatstat(K)
rpcf_i = list(PCF[0])
pcf_i = list(PCF[2])
Rk_sims.append(rk_i)
K_sims.append(K_i)
Rpcf_sims.append(rpcf_i)
PCF_sims.append(pcf_i)
print(time.process_time() - start)
assert numpy.allclose(np.mean(Rk_sims, axis=0), np.array(Rk_sims[0]))
assert numpy.allclose(np.mean(Rpcf_sims, axis=0), np.array(Rpcf_sims[0]))
L = np.divide(np.sqrt(K_sims),np.sqrt(np.pi))
L_cat = np.sqrt(np.array(K_cat) / np.pi)
Results[model.name] = {'K_sims': K,
'Rk_sims': Rk_sims,
'K_cat': K_cat,
'rk_cat': rk_cat,
'L_sims': L,
'L_cat':L_cat,
'PCF_sims': PCF_sims,
'Rpcf_sims':Rpcf_sims,
'pcf_cat': pcf_cat,
'rpcf_cat': rpcf_cat}
with open(paths.get_kripley_result_path('KRipley_bolletino', 10), 'wb') as file_:
pickle.dump(Results, file_)
return Results
if __name__ == "__main__": if __name__ == "__main__":
Results = run_parallel() Results = run()
# # for cat in ['bollettino', 'horus', 'emrcmt'][:1]:
with open(paths.get_kripley_result_path('KRipley_testpcf_bollettino', 10), 'rb') as file_: # with open(paths.get_kripley_result_path('K_%s' % cat, 10), 'rb') as file_:
Results = pickle.load(file_) # Results = pickle.load(file_)
a = plot_results(Results) # a = plot_results(Results)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment