Skip to content
Snippets Groups Projects
Commit 0fae3c42 authored by Sumbel Shangareeva's avatar Sumbel Shangareeva
Browse files

spotpy script for optimization

parent 0a74b9dd
No related branches found
No related tags found
No related merge requests found
date,1
2014-01-01,4.75
2024-01-01,4.75
import os
import numpy as np
import spotpy
import f90nml as nml
from source.spotpy_optim.spot_setup import spot_setup
from plot import main_plot
parameters_nml_path = os.path.join('parameters_config.nml')
config_nml_path = os.path.join('ui1_config.nml')
res_path = 'results/spotpy_exp'
if not os.path.exists(res_path):
os.mkdir(res_path)
n_it = 100
setup = spot_setup(setup_path=config_nml_path, params_path=parameters_nml_path)
sampler = spotpy.algorithms.rope(setup, dbname=f'{res_path}/arkha_rope_{n_it}',
dbformat='csv',
save_sim=False
)
sampler.sample(n_it)
results = sampler.getdata()
kd_in_optim = spotpy.analyser.get_best_parameterset(results)[0]
params_nml = nml.read(parameters_nml_path)
kd_in = np.array(params_nml['parameters_config_namelist']['kd_in'])
kd_in[0, 5] = kd_in_optim
params_nml['parameters_config_namelist']['kd_in'] = kd_in.tolist()
params_nml.write(parameters_nml_path, force=True)
os.system(f'./run.exe')
main_plot()
&parameters_config_namelist
kd_in(1:6,1) = 0.03, 0.03, 0.0402, 0.08581933, 0.012, 0.07
kd_in(1:6,1) = 0.03, 0.03, 0.0402, 0.08581933, 0.012, 0.07087475
kd_in(1:6,2) = 0.03, 0.03, 0.05067, 0.08581933, 0.007, 0.007
kd_in(1:6,3) = 0.035, 0.035, 0.04185, 0.08581933, 0.007, 0.007
kd_in(1:6,4) = 0.03, 0.03, 0.0, 0.0, 0.0, 0.0
......
......@@ -194,7 +194,7 @@ def pools_sums_plots(file_names: List[str], obs_data: pd.DataFrame, out_path: st
plt.close()
def main():
def main_plot():
obs_data = pd.read_csv(obs_data_path)
obs_data['date'] = pd.to_datetime(obs_data['date'])
out_dir = ''
......@@ -213,4 +213,4 @@ def main():
if __name__ == '__main__':
main()
main_plot()
matplotlib==3.9.1
numpy==2.0.0
pandas==2.2.2
f90nml~=1.4.4
spotpy~=1.6.6
\ No newline at end of file
......@@ -35,25 +35,25 @@ program main
call carbon_pools_init()
print*
print('(1x,a,1x,10x,1x,a)'), 'модель углеродного цикла ', trim(carbon_model_type)
print('(1x,a,1x,10x,1x,a)'), 'данные об окружающей среде ', trim(environment_data_type)
if (environment_data_type == 'lsm_offline') then
print('(1x,a,1x,10x,1x,a)'), ' файл ', trim(lsm_datafile)
print('(1x,a,1x,f10.2)'), ' шаг по широте ', dlat_nc
print('(1x,a,1x,f10.2)'), ' шаг по долготе ', dlon_nc
print('(1x,a,1x,i10)'), ' шаг по времени, c ', dt_nc
endif
print*,'---'
print('(1x,a,1x,f10.2)'), 'шаг по широте ', dlat
print('(1x,a,1x,f10.2)'), 'шаг по долготе ', dlon
print('(1x,a,1x,i10)'), 'всего ячеек в области ', ncell_tot
print('(1x,a,1x,i10,1x,"(",i3,"%)")'), ' из них обсчитывается ', ncell_mask, nint(ncell_mask*100./ncell_tot)
print*,'---'
print('(1x,a,1x,i10)'), 'шаг по времени, с ', dt
print('(1x,a,1x,i10)'), 'число шагов по времени ', ntime
print('(1x,a,1x,10x,a)'), 'стартовая дата ', date_fst%timestamp
print('(1x,a,1x,10x,a)'), 'финальная дата ', date_lst%timestamp
! print*
! print('(1x,a,1x,10x,1x,a)'), 'модель углеродного цикла ', trim(carbon_model_type)
! print('(1x,a,1x,10x,1x,a)'), 'данные об окружающей среде ', trim(environment_data_type)
! if (environment_data_type == 'lsm_offline') then
! print('(1x,a,1x,10x,1x,a)'), ' файл ', trim(lsm_datafile)
! print('(1x,a,1x,f10.2)'), ' шаг по широте ', dlat_nc
! print('(1x,a,1x,f10.2)'), ' шаг по долготе ', dlon_nc
! print('(1x,a,1x,i10)'), ' шаг по времени, c ', dt_nc
! endif
! print*,'---'
! print('(1x,a,1x,f10.2)'), 'шаг по широте ', dlat
! print('(1x,a,1x,f10.2)'), 'шаг по долготе ', dlon
! print('(1x,a,1x,i10)'), 'всего ячеек в области ', ncell_tot
! print('(1x,a,1x,i10,1x,"(",i3,"%)")'), ' из них обсчитывается ', ncell_mask, nint(ncell_mask*100./ncell_tot)
! print*,'---'
! print('(1x,a,1x,i10)'), 'шаг по времени, с ', dt
! print('(1x,a,1x,i10)'), 'число шагов по времени ', ntime
! print('(1x,a,1x,10x,a)'), 'стартовая дата ', date_fst%timestamp
! print('(1x,a,1x,10x,a)'), 'финальная дата ', date_lst%timestamp
call environment_init()
......@@ -95,17 +95,17 @@ program main
if (firstcall) then
if (tt == cputime_averaging_period) then
call cpu_time(cputime1)
print*
print('(1x,a,4x,f7.2)'), 'примерная прод-ть расчета, мин ', ntime*(cputime1-cputime0)/(60.*cputime_averaging_period)
print*
print*, '----------------------------------------------------------------'
print*
! print*
! print('(1x,a,4x,f7.2)'), 'примерная прод-ть расчета, мин ', ntime*(cputime1-cputime0)/(60.*cputime_averaging_period)
! print*
! print*, '----------------------------------------------------------------'
! print*
firstcall = .false.
endif
else
progress_pc = floor(tt*100./ntime)
if (progress_pc /= progress_pc_mem) then
print'(1x,a,4x,i3,1x,"%")', date_c%timestamp, progress_pc
! print'(1x,a,4x,i3,1x,"%")', date_c%timestamp, progress_pc
progress_pc_mem = progress_pc
endif
endif
......
......@@ -2,17 +2,80 @@ import os
import f90nml as nml
import numpy as np
import pandas as pd
import spotpy
BASE_DIR = '../../'
parameters_nml_path = os.path.join(BASE_DIR, 'parameters_config.nml')
config_nml_path = os.path.join(BASE_DIR, 'ui1_config.nml')
class spot_setup(object):
# add start date
def __init__(self, setup_path: str, params_path: str,
obj_func=None):
self.station_opt = None
self.obj_func = obj_func
self.params = [
spotpy.parameter.Uniform('kd_in', low=0.01, high=0.5, optguess=0.07),
]
self.setup_path = setup_path
self.params_path = params_path
self.model_out_path = None
self.obs_data_path = None
self.read_setup_data()
self.obs_data = None
self.read_obs_data()
params_nml = nml.read(parameters_nml_path)
config_nml = nml.read(config_nml_path)['config_namelist']
kd_in = np.array(params_nml['parameters_config_namelist']['kd_in'])
kd_in[0, 5] = 0.07
def read_setup_data(self):
config_nml = nml.read(self.setup_path)['config_namelist']
self.station_opt = int(config_nml['station_opt'])
self.model_out_path = os.path.join('results',
config_nml['carbon_model_type'],
f"{config_nml['station_name']}_{self.station_opt}.txt")
self.obs_data_path = os.path.join('data', f"obs_data_{config_nml['station_name']}.csv")
def read_obs_data(self):
obs_data = pd.read_csv(self.obs_data_path)
self.obs_data = obs_data.iloc[0,1]
def read_model_output(self):
model_out = pd.read_csv(self.model_out_path, sep=';')
model_out = model_out.iloc[-1, 1:].sum()
return model_out
def simulation(self, x):
params_nml = nml.read(self.params_path)
kd_in = np.array(params_nml['parameters_config_namelist']['kd_in'])
kd_in[self.station_opt - 1, 5] = x[0]
params_nml['parameters_config_namelist']['kd_in'] = kd_in.tolist()
params_nml.write(parameters_nml_path, force=True)
params_nml.write(self.params_path, force=True)
os.system(f'./run.exe')
model_out = self.read_model_output()
return model_out
def evaluation(self):
return self.obs_data
def objectivefunction(self, simulation, evaluation, params=None):
# SPOTPY expects to get one or multiple values back,
# that define the performance of the model run.
if not self.obj_func:
# This is used if not overwritten by user
like = -abs(evaluation - simulation)
else:
# Way to ensure flexible spot setup class
like = self.obj_func(simulation, evaluation)
return like
def parameters(self):
return spotpy.parameter.generate(self.params)
......@@ -135,7 +135,7 @@
if_datafile_read_at_first = .false.
if_standard_print = .true.
if_standard_print = .false.
if_standard_output = .false.
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment