Commit e1d46236 authored by Antoine Berchet's avatar Antoine Berchet
Browse files

Plotting the article figures during test

parent 4e52c710
verbose: 2
logfile: pycif.logtest
workdir: /home/satellites13/aberchet/cif_paper/simulations_20210323/test_integration_inversion_dum0
workdir: /home/satellites13/aberchet/cif_paper/simulations_20210324/test_integration_inversion_dum0
datei: 2010-01-01
datef: 2010-01-05 00:00:00
mode:
......@@ -35,7 +35,7 @@ model:
name: dummy
version: std
obsvect:
dir_obsvect: /home/satellites13/aberchet/cif_paper/simulations_20210323/test_integration_inversion_dum0/../ref_obsvect/
dir_obsvect: /home/satellites13/aberchet/cif_paper/simulations_20210324/test_integration_inversion_dum0/../ref_obsvect/
dump_type: nc
plugin:
name: standard
......@@ -47,7 +47,7 @@ datavect:
CH4:
duration: 2H
frequency: 1H22min
nstations: 50
nstations: 5
plugin:
name: random
type: measurements
......@@ -77,7 +77,7 @@ datavect:
square: null
variable: zlon
hcorrelations:
dircorrel: /home/satellites13/aberchet/cif_paper/simulations_20210323/test_integration_inversion_dum0/datavect/
dircorrel: /home/satellites13/aberchet/cif_paper/simulations_20210324/test_integration_inversion_dum0/datavect/
dump_hcorr: true
evalmin: 0
landsea: false
......
verbose: 2
logfile: pycif.logtest
workdir: /home/satellites13/aberchet/cif_paper/simulations_20210323/test_integration_inversion_dum5
workdir: /home/satellites13/aberchet/cif_paper/simulations_20210324/test_integration_inversion_dum1
datei: 2010-01-01
datef: 2010-01-05 00:00:00
mode:
......@@ -35,7 +35,7 @@ model:
name: dummy
version: std
obsvect:
dir_obsvect: /home/satellites13/aberchet/cif_paper/simulations_20210323/test_integration_inversion_dum5/../ref_obsvect/
dir_obsvect: /home/satellites13/aberchet/cif_paper/simulations_20210324/test_integration_inversion_dum1/../ref_obsvect/
dump_type: nc
plugin:
name: standard
......@@ -47,7 +47,7 @@ datavect:
CH4:
duration: 2H
frequency: 1H22min
nstations: 50
nstations: 5
plugin:
name: random
type: measurements
......@@ -77,11 +77,11 @@ datavect:
square: null
variable: zlon
hcorrelations:
dircorrel: /home/satellites13/aberchet/cif_paper/simulations_20210323/test_integration_inversion_dum5/datavect/
dircorrel: /home/satellites13/aberchet/cif_paper/simulations_20210324/test_integration_inversion_dum1/datavect/
dump_hcorr: true
evalmin: 0
landsea: false
sigma: 200000.0
sigma: 500
hresol: hpixels
nlev: 1
period: 5D
......
verbose: 2
logfile: pycif.logtest
workdir: /home/satellites13/aberchet/cif_paper/simulations_20210324/test_integration_inversion_dum13
datei: 2010-01-01
datef: 2010-01-05 00:00:00
mode:
minimizer:
df1: 0.01
epsg: 0.02
maxiter: 1
nsim: 1
plugin:
name: congrad
version: std
simulator:
plugin:
name: gausscost
version: std
reload_from_previous: true
plugin:
name: 4dvar
version: std
save_out_netcdf: true
controlvect:
plugin:
name: standard
version: std
save_out_netcdf: true
model:
chemistry:
acspecies:
CH4: null
file_pg: /home/users/aberchet/CIF/model_sources/dummy_gauss/Pasquill-Gifford.txt
plugin:
name: dummy
version: std
obsvect:
dir_obsvect: /home/satellites13/aberchet/cif_paper/simulations_20210324/test_integration_inversion_dum13/../ref_obsvect/
dump_type: nc
plugin:
name: standard
version: std
datavect:
components:
concs:
parameters:
CH4:
duration: 2H
frequency: 1H22min
nstations: 5
plugin:
name: random
type: measurements
version: param
random_subperiod_shift: true
seed: true
zmax: 100
fluxes:
parameters:
CH4:
err: 1
errtype: max
flx_formula:
- product:
- sum:
- cos: null
period: 500
variable: zlat
- period: 1000
sin: null
variable: zlon
- sum:
- period: 1000
square: null
variable: zlat
- period: 1000
square: null
variable: zlon
hcorrelations:
dircorrel: /home/satellites13/aberchet/cif_paper/simulations_20210324/test_integration_inversion_dum13/datavect/
dump_hcorr: true
evalmin: 0
landsea: false
sigma: 500
hresol: global
nlev: 1
period: 5D
plugin:
name: dummy
type: fluxes
version: txt
tcorrelations:
sigma_t: 12
vresol: vpixels
meteo:
plugin:
name: dummy
type: meteo
version: csv
resolution: 1H
seed: true
plugin:
name: standard
version: std
domain:
nlat: 15
nlon: 30
plugin:
name: dummy
version: std
xmax: 2500
xmin: 0
ymax: 2000
ymin: 0
platform:
plugin:
name: LSCE
version: obelix
verbose: 2
logfile: pycif.logtest
workdir: /workdir/
datei: 2010-01-01
datef: 2010-01-05 00:00:00
mode:
obserror: 0.01
perturb_obsvect: true
plugin:
name: forward
version: std
controlvect:
plugin:
name: standard
version: std
save_out_netcdf: true
model:
chemistry:
acspecies:
CH4: null
MCF: null
file_pg: /tmp/CIF/model_sources/dummy_gauss/Pasquill-Gifford.txt
plugin:
name: dummy
version: std
obsvect:
dir_obsvect: /home/satellites13/aberchet/cif_paper/simulations_20210324/test_integration_inversion_dum2/../ref_obsvect/
dump_type: nc
plugin:
name: standard
version: std
datavect:
components:
concs:
parameters:
CH4:
duration: 2H
frequency: 1H22min
nstations: 5
plugin:
name: random
type: measurements
version: param
random_subperiod_shift: true
seed: true
zmax: 100
MCF:
duration: 4H
frequency: 3H33min
nstations: 10
plugin:
name: random
type: measurements
version: param
random_subperiod_shift: true
seed: true
seed_id: 100
zmax: 100
fluxes:
parameters:
CH4:
err: 1
errtype: max
flx_formula:
- product:
- sum:
- cos: null
period: 500
variable: zlat
- period: 1000
sin: null
variable: zlon
- sum:
- period: 1000
square: null
variable: zlat
- period: 1000
square: null
variable: zlon
hcorrelations:
dircorrel: /home/satellites13/aberchet/cif_paper/simulations_20210324/test_integration_inversion_dum2/datavect/
dump_hcorr: true
evalmin: 0
landsea: false
sigma: 500
hresol: hpixels
nlev: 1
period: 5D
plugin:
name: dummy
type: fluxes
version: txt
tcorrelations:
sigma_t: 12
vresol: vpixels
MCF:
err: 1
errtype: max
flx_formula:
- product:
- sum:
- cos: null
period: 100
variable: zlat
- period: 200
sin: null
variable: zlon
- sum:
- period: 500
square: null
variable: zlat
- period: 600
square: null
variable: zlon
hcorrelations:
dircorrel: /home/satellites13/aberchet/cif_paper/simulations_20210324/test_integration_inversion_dum2/datavect/
dump_hcorr: true
evalmin: 0
landsea: false
sigma: 500
hresol: hpixels
nlev: 1
period: 5D
plugin:
name: dummy
type: fluxes
version: txt
vresol: vpixels
meteo:
plugin:
name: dummy
type: meteo
version: csv
resolution: 1H
seed: true
plugin:
name: standard
version: std
domain:
nlat: 15
nlon: 30
plugin:
name: dummy
version: std
xmax: 2500
xmin: 0
ymax: 2000
ymin: 0
......@@ -35,7 +35,7 @@ model:
name: dummy
version: std
obsvect:
dir_obsvect: /home/satellites13/aberchet/cif_paper/simulations_20210323/test_integration_fwd_dummy_con0/../ref_obsvect/
dir_obsvect: /home/satellites13/aberchet/cif_paper/simulations_20210324/test_integration_fwd_dummy_con0/../ref_obsvect/
dump_type: nc
plugin:
name: standard
......@@ -47,7 +47,7 @@ datavect:
CH4:
duration: 2H
frequency: 1H22min
nstations: 50
nstations: 5
plugin:
name: random
type: measurements
......@@ -81,7 +81,7 @@ datavect:
sin: null
variable: zlon
hcorrelations: &id001
dircorrel: /home/satellites13/aberchet/cif_paper/simulations_20210323/test_integration_fwd_dummy_con0/datavect/
dircorrel: /home/satellites13/aberchet/cif_paper/simulations_20210324/test_integration_fwd_dummy_con0/datavect/
dump_hcorr: true
evalmin: 0
landsea: false
......@@ -143,7 +143,7 @@ datavect:
square: null
variable: zlon
hcorrelations:
dircorrel: /home/satellites13/aberchet/cif_paper/simulations_20210323/test_integration_fwd_dummy_con0/datavect/
dircorrel: /home/satellites13/aberchet/cif_paper/simulations_20210324/test_integration_fwd_dummy_con0/datavect/
dump_hcorr: true
evalmin: 0
landsea: false
......
......@@ -74,4 +74,4 @@ def execute(self, **kwargs):
to_netcdf=self.save_out_netcdf,
dir_netcdf="{}/controlvect/".format(workdir),
)
return chiopt
return controlvect, obsvect
......@@ -22,7 +22,7 @@ def ref_dummy(tmpdir):
"logfile": "pycif.logtest",
"workdir": tmpdir_str,
"datei": datetime.date(2010, 1, 1),
"datef": datetime.date(2010, 1, 5),
"datef": datetime.date(2010, 1, 2),
"mode": {
"plugin": {"name": "forward", "version": "std"},
"perturb_obsvect": True,
......@@ -151,7 +151,7 @@ def ref_dummy(tmpdir):
"type": "measurements"
},
"frequency": "1H22min",
"nstations": 50,
"nstations": 5,
"duration": "2H",
"random_subperiod_shift": True,
"zmax": 100,
......
......@@ -11,10 +11,10 @@ from pycif.utils.yml import ordered_dump
params=
[
{"resol": "full"},
pytest.param({"resol": "full", "correlations": 2e5},
marks=pytest.mark.article),
{"resol": "bands"},
{"resol": "global"},
# pytest.param({"resol": "full", "correlations": 2e5},
# marks=pytest.mark.article),
# {"resol": "bands"},
# {"resol": "global"},
])
def dummy_config_inversion(dummy_config_fwd, request):
"""
......
import os
import time
import shutil
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
import pytest
......@@ -16,9 +18,10 @@ from pycif.utils.path import init_dir
@pytest.mark.parametrize(
"settings", [
{"mode": "4dvar", "minimizer": "M1QN3"},
{"mode": "4dvar", "minimizer": "congrad"},
{"mode": "analytical"},
{"mode": "ensrf"}]
# {"mode": "4dvar", "minimizer": "congrad"},
# {"mode": "analytical"},
# {"mode": "ensrf"}]
]
)
def test_integration_inversion(dummy_config_inversion, settings):
"""
......@@ -141,35 +144,109 @@ def test_integration_inversion(dummy_config_inversion, settings):
with open(dummy_config_file, "w") as outfile:
ordered_dump(outfile, config)
# Compute the value of the cost function for a varying number of simulation
if os.path.isdir("{}/varying_cost/".format(tmpdir)):
shutil.rmtree("{}/varying_cost/".format(tmpdir))
init_dir("{}/varying_cost/".format(tmpdir))
# Loop with different number of simulations
for nsim in range(10, nsimmax + 1, 10):
file_config = "{}/dummy_config.yml".format(tmpdir)
inv_setup = Setup.from_yaml(file_config)
obs_root = "obsvect_posterior"
control_root = "controlvect"
title = "Analytical"
if settings.get("minimizer", None) == "congrad":
title = "CONGRAD"
obs_root = "obsoperator/fwd_-002/obsvect"
control_root = "obsoperator/fwd_-002/controlvect"
inv_setup["mode"]["minimizer"]["maxiter"] = int(nsim / 2)
if os.path.isdir("{}/obsoperator/fwd_-002".format(tmpdir)):
shutil.rmtree("{}/obsoperator/fwd_-002".format(tmpdir))
elif settings.get("minimizer", None) == "M1QN3":
title = "M1QN3"
obs_root = "obsoperator/fwd_-002/obsvect"
control_root = "obsoperator/fwd_-002/controlvect"
inv_setup["mode"]["minimizer"]["maxiter"] = int(nsim / 2)
inv_setup["mode"]["minimizer"]["nsim"] = int(nsim / 2)
if os.path.isdir("{}/obsoperator/fwd_-002".format(tmpdir)):
shutil.rmtree("{}/obsoperator/fwd_-002".format(tmpdir))
elif settings["mode"] == "EnSRF":
title = "EnSRF"
inv_setup["mode"]["nsample"] = nsim
# Execute a partial inversion with less simulations
inv_setup = Setup.from_dict(inv_setup, convert_none=True)
inv_setup = inv_setup.load_config(inv_setup)
result = inv_setup.mode.execute()
controlvect, obsvect = inv_setup.mode.execute()
# Compute the cost function
departures = obsvect.ysim - obsvect.yobs
j_o = 0.5 * (departures * obsvect.rinvprod(departures)).sum()
dx = controlvect.x - controlvect.xb
bfull = np.linalg.inv(
inv_setup.controlvect
.build_b(inv_setup.controlvect)
)
j_b = (bfull.dot(dx[:, np.newaxis]) * dx).sum()
# Dump to text file for later plot
with open("{}/varying_cost_function.txt".format(tmpdir), "a") as f:
f.write("{} {} {}\n".format(nsim, j_o, j_b))
# Domain limits
xmin = inv_setup.domain.xmin
xmax = inv_setup.domain.xmax
ymin = inv_setup.domain.ymin
ymax = inv_setup.domain.ymax
# Read observations
file_obs = "{}/obsvect/concs/CH4/monitor.nc".format(tmpdir)
monitor_ref = read_datastore(file_obs)
coords = monitor_ref.loc[:, ["lon", "lat", "alt"]].drop_duplicates()
shutil.copytree("{}/{}".format(tmpdir, obs_root),
"{}/varying_cost/{:04d}".format(tmpdir, nsim))
# Compute fluxes from control vector
file_flx = "{}/{}/controlvect/fluxes/" \
"controlvect_fluxes_CH4.nc".format(tmpdir, control_root)
ds = xr.open_dataset(file_flx)
dflx = ds["x_phys"].mean(axis=(0, 1)) - ds["xb_phys"].mean(axis=(0, 1))
# Fetch resolution for figure name
resol = ""
if config["datavect"]["components"][
"fluxes"]["parameters"]["CH4"]["hresol"] == "hpixels":
if config["datavect"]["components"][
"fluxes"]["parameters"]["CH4"]["hcorrelations"]["sigma"] == 500:
resol = "lowcorr"
else:
resol = "highcorr"
elif config["datavect"]["components"]["fluxes"] \
["parameters"]["CH4"]["hresol"] == "bands":
resol = "bands"
elif config["datavect"]["components"]["fluxes"] \
["parameters"]["CH4"]["hresol"] == "global":
resol = "global"
# Plot the figure
plt.figure(figsize=(21, 11))
ax0 = plt.axes([0.05, 0.05, 0.73, 0.87])
im = plt.imshow(dflx, extent=(xmin, xmax, ymin, ymax), cmap="YlOrRd",
vmin=-0.2, vmax=0.5)
sc = plt.scatter(coords["lon"], coords["lat"], c=coords["alt"],
cmap="Blues", linewidths=1, edgecolors="k", s=600)
plt.yticks(fontsize=25)
plt.xticks(fontsize=25)
ax1 = plt.axes([0.74, 0.05, 0.05, 0.87])
cb1 = plt.colorbar(im, cax=ax1)
plt.yticks(fontsize=25)
plt.ylabel("Fluxes (a.u.)", fontsize=30)
ax2 = plt.axes([0.86, 0.05, 0.05, 0.87])
cb2 = plt.colorbar(sc, cax=ax2)
plt.yticks(fontsize=25)
plt.ylabel("Station altitude (m a.g.l)", fontsize=30)
ax0.set_title(title, fontsize=45)
plt.savefig("{}/map_dx_{}_{}.pdf".format(tmpdir, title, resol))
plt.close()
Supports Markdown
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