Source code for pycif.plugins.models.TM5.io.native2inputs_adj
import pandas as pd
import datetime
import os
import glob
from netCDF4 import Dataset
import xarray as xr
import numpy as np
# JvP 20210802: added import statements, MODULE_NAME module level logger
import logging
import sys
from .....utils.check.errclass import CifRuntimeError
from .....utils.hdf5 import _hdf5_lock
MODULE_NAME = __name__[__name__.index('TM5'):] if 'TM5' in __name__ else __name__
logger = logging.getLogger(MODULE_NAME)
# Original native2inputs_adj function by A. Berchet
[docs]
def native2inputs_adj_AB(self, datastore, input_type, datei, datef, runsubdir, mode="fwd", **kwargs):
"""Read sensitivity to input and converts them to pyCIF data
Args:
self: the model Plugin
input_type (str): one of 'fluxes'
datastore: data to convert
if input_type == 'fluxes',
datei, datef: date interval of the sub-simulation
mode (str): running mode: one of 'fwd', 'adj' and 'tl'
runsubdir (str): sub-directory for the current simulation
workdir (str): the directory of the whole pyCIF simulation
Notes:
- CHIMERE expects hourly inputs;
"""
if datastore == {}:
return datastore
ddi = min(datei, datef)
# List of TM5 dates
dref = datetime.datetime.strptime(
os.path.basename(os.path.normpath(runsubdir)), "%Y-%m-%d_%H-%M"
)
list_dates = self.input_dates[ddi]
# Reading only output files related to given input_type
ref_names = {
"inicond": "corresponding file name",
"flux": "corresponding file name",
}
# If looking for sensitivity to input not provided by TM5, do nothing
if input_type not in ref_names:
return datastore
# Otherwise, loop over component/parameters in the datastore
for trid in datastore:
# Here you need to read the correspondig sensitivity
# for that tracer and input_type
ndate = len(list_dates)
nmerid = self.domain.nlat
nzonal = self.domain.nlon
if trid[0] == "flux":
nvert = 1
adj_file = f"{runsubdir}/output/adj_emissions.nc4"
with _hdf5_lock:
with Dataset(adj_file) as f:
group = f.groups["glb600x400"].groups["total"]
data = group.variables["adj_emis"][:][:, np.newaxis]
# Then you put in "adj_out" as a xr.DataArray
datastore[trid]["adj_out"] = xr.DataArray(
data,
coords={"time": list_dates},
dims=("time", "lev", "lat", "lon"),
)
elif trid[0] == "inicond":
nvert = 25
adj_file = dref.strftime(
f"{runsubdir}/output/save/adj_iniconc_%Y%m%d0000.nc")
with _hdf5_lock:
with Dataset(adj_file) as f:
data = f.groups["glb600x400"].variables[trid[1]][:]
data = data[np.newaxis]
# Then you put in "adj_out" as a xr.DataArray
datastore[trid]["adj_out"] = xr.DataArray(
data,
coords={"time": [dref]},
dims=("time", "lev", "lat", "lon"),
)
return datastore
# New native2inputs_adj function by J.C.A. van Peet
[docs]
def native2inputs_adj(
self, datastore, input_type, datei, datef, runsubdir,
check_transforms=False,
mode="fwd", **kwargs):
"""
PURPOSE
Read sensitivity to input and converts them to pyCIF data.
ARGS
self = the model Plugin
datastore = data to convert if input_type == 'fluxes'
input_type = (str) one of 'fluxes'
datei = initial date of the sub-simulation interval
datef = final date of the sub-simulation interval
runsubdir = (str) sub-directory for the current simulation
KWARGS
mode = (str) running mode, one of 'fwd', 'adj' and 'tl'
workdir = (str) the directory of the whole pyCIF simulation
NOTE
None
VERSION HISTORY
2.1 DD-MM-YYYY by J.C.A. van Peet
Continued development
2.0 02-08-2021 by J.C.A. van Peet
Adapted for TM5.
1.0 ??-??-???? by A. Berchet
See original code above.
"""
# Set the name of this function
PROG_NAME = MODULE_NAME+".native2inputs_adj"
# Local logger
logger = logging.getLogger(PROG_NAME)
logger.setLevel(logging.DEBUG)
# Return if datastore is empty
if( datastore == {} ):
return datastore
# end if
# Return if input_type is something else than "flux".
if( input_type != 'flux' ):
return datastore
# end if
# Some debug statements...
logger.debug("")
logger.debug("*"*30)
logger.debug(PROG_NAME+" => DEBUG:")
#logger.debug(" self = %s", str(self))
#logger.debug(" dir(self) = %s", str(("\n"+" "*27).join( [ ', '.join(dir(self)[i:i+5]) for i in range(0,len(dir(self)),5) ] ) ))
logger.debug(f" datastore = {datastore}" )
logger.debug(f" input_type = {input_type}")
logger.debug(f" datei = {datei}" )
logger.debug(f" datef = {datef}" )
logger.debug(f" runsubdir = {runsubdir}" )
logger.debug(f" mode = {mode}" )
logger.debug(f" kwargs = {kwargs}" )
# Get a list of dates
# Not sure if dref is required, see function by Antoine above.
ddi = min(datei, datef)
list_dates = self.input_dates[ddi]
logger.debug(f" list_dates = {list_dates}" )
# Loop over component / parameters in the datastore
for trid in datastore:
# Log the "trid"...
logger.debug
# Get the adjoint fluxes
if trid[0] == "flux":
# As far as I can follow the TM5 code, the output in adj_emissions.nc4
# contains the updates of the emission factors relative to the aprior
# emissions. These factors are copied in the g2gc job to a file called
# g_obs.nc (see jobs in TMVar.py). The long_name attribute changes from
# "adjoint emission factors" to "gradient towards relative deviation in
# emission of CH4 for category XXX", and the units remain 1. These
# adjoint emission factors are then preconditioned and written to file
# gc_obs. The long_name attribute is changed to "pre-conditioned gradient
# towards relative deviation in emission of CH4 for category XXX".
# Then something similar is done for something that is called gc_x, which
# I assume is the gradient of the background term. Next, gc_obs and gc_x
# are combined into a file called gc.nc. The long_name attribute in gc.nc
# is still "pre-conditioned gradient towards relative deviation in
# emission of CH4 for category XXX". Then some norms are calculated and
# written to csv file. Finally, the contents of gc.nc appear to be flattened
# and written to a file called dJ_dxc.nc which serves as input for the next
# step.
# The results from the g2gc job (dJ_dxc.nc) and the cost function file
# from the dep job (J.nc) are both used in the opt job (optimizer) to take
# a step in the right direction to minimise the cost function. This step
# creates two files called optimizer-input.nc and optimizer-output.nc
# respectively.
# The files from the opt job are probably read by the next job, xc2x.
# This job saves the preconditioned state to a file called xc.nc, where
# the emissions have a long_name attribute "pre-conditioned relative
# deviation in emission of CH4 for category XXX". The preconditioning is
# removed, and the result is saved into a file called x.nc, where the
# emission long_name attribute is "relative deviation in emission of CH4
# for category XXX".
# Finally the x.nc file is converted to actual emissions by the function
# State2Model (TMVar.py, around line 2140), which converts these relative
# deviations to actual emissions as follows:
#
# * Emissions input is created be expanding the factors 'x' in the state
# and applying them to the base emissions 'em0':
#
# em = { em0 exp( x) , x < 0
# { em0 (1 + x) , x >= 0
#
# With em0 the prior emissions and x the values derived from the
# adj_emissions.nc4 file by jobs g2gc and opt.
#
# Throughout all these processing steps, the "adjoint emissions" remain
# relative deviations and only in the last step are they converted into
# actual emissions that can be used in the next iteration. But I guess I
# have to pass the actual sensitivities back to CIF, not the realtive ones,
# so I have to convert the contents of adj_emissions.nc4 to something like
# "gradient towards ABSOLUTE deviation in emission of CH4 for category XXX".
# So I will assume that I can do that by inserting the contents of
# adj_emissions.nc4 into the "em" equation above, and pass the result back
# to CIF...
# The first step is to get the prior emissions.
#
# "self" is the TM5 model plugin, which contains a variable called fluxes,
# which points to the fluxes plugin, which contains the read function you
# also elsewhere to read the original TM5 fluxes. You can check that by
# printing the following:
# logger.debug(" self.flux.read.__doc__ = %s" % self.flux.read.__doc__)
# The directory .../pycif/plugin/fluxes/tm5/ also seems like the logical
# location to store a reading function for the adjoint fluxes...
#
# Anyway, a link to the file with prior emissions is created in
# .../plugins/models/TM5/io/inputs/make_fluxes.py
# with the name prior_emissions.nc, so you can just read that file here.
#
#logger.debug("")
#logger.debug(" self.flux = %s" % str(self.flux))
#logger.debug(" dir(self.flux) = %s" % str(("\n"+" "*27).join( [ ', '.join(dir(self.flux)[i:i+5]) for i in range(0,len(dir(self)),5) ] ) ))
#logger.debug(" self.flux.read.__doc__ = %s" % self.flux.read.__doc__)
#logger.debug(" self.input_dates = %s" % self.input_dates )
#logger.debug(" os.getcwd() = %s" % os.getcwd() )
#logger.debug("")
#
name = 'CH4'
tracdir = runsubdir+'/'
tracfile = 'prior_emissions.nc4'
varnames = None
#
# The function fluxes.read is apparently also called from line 97 in
# .../pycif/plugins/transforms/fromcontrol/forward.py
# which refers to the read function in
# .../pycif/plugins/fluxes/tm5/read.py
# If you can believe the traceback from that latter function following
# an induces RuntimeError. In the function forward.py, the dates passed
# to the read function are determined by:
# ddi = min(di, df)
# So here, I've replaced di with self.datei and df with self.datef.
#2021-10-26 14:08:26: self.input_dates = {datetime.datetime(2018, 1, 1, 0, 0): DatetimeIndex(['2018-01-01'], dtype='datetime64[ns]', freq='MS')}
#ddi = min(self.datei, self.datef)
#dates = self.input_dates[ddi]
#
prior_emissions = self.flux.read(
name, varnames,
datastore[trid]["input_dates"][ddi],
datastore[trid]["input_files"][ddi])
# Now read the file with posterior emissions. Whatever the unit is...
# The function read_adj can be found in
# .../pycif/plugins/fluxes/tm5/read_adj.py
# Note that you've also updated the __init__.py in the same directory.
filename = tracdir+'output/adj_emissions.nc4'
adj_emissions = self.flux.read_adj(filename)
# The prior_emissions will have a shape of [time, lev, lat, lon], where lev
# will always be 1. There is no time info in the adj_emissions.nc4 file,
# so it will have a shape of [time, lat, lon]. In order to easily combine
# both arrays, add an axis in adj_emissions, so that the shapes are the same.
adj_emissions = np.expand_dims(adj_emissions, 1)
# Check if the shape of prior_emissions and adj_emissions match.
if prior_emissions.shape != adj_emissions.shape:
logger.critical("")
logger.critical("*"*30)
logger.critical(PROG_NAME+" => ERROR: The shapes of prior_emissions and adj_emissions do not comply!")
logger.critical(f" prior_emissions.shape = {prior_emissions.shape}" )
logger.critical(f" adj_emissions.shape = {adj_emissions.shape}" )
logger.critical("*"*30)
logger.critical("")
raise CifRuntimeError
# end if
# Convert the contents of adj_emissions to actual emissions.
# em = { em0 exp( x) , x < 0
# { em0 (1 + x) , x >= 0
# See note above for more information.
adj_emissions_orig = np.copy( adj_emissions )
# mask = adj_emissions < 0
# np.exp (adj_emissions, out=adj_emissions, where=mask )
# np.add (adj_emissions, 1.0, out=adj_emissions, where=~mask)
# np.multiply(adj_emissions, prior_emissions, out=adj_emissions )
logger.debug(f" adj_emissions_orig mean = {np.mean(adj_emissions_orig)}")
logger.debug(f" adj_emissions_orig sdev = {np.std(adj_emissions_orig)}")
logger.debug(f" adj_emissions_orig min = {np.amin(adj_emissions_orig)}")
logger.debug(f" adj_emissions_orig max = {np.amax(adj_emissions_orig)}")
logger.debug(f" adj_emissions mean = {np.mean(adj_emissions)}")
logger.debug(f" adj_emissions sdev = {np.std(adj_emissions)}")
logger.debug(f" adj_emissions min = {np.amin(adj_emissions)}")
logger.debug(f" adj_emissions max = {np.amax(adj_emissions)}")
# Put the adjoint emissions in datastore as an xr.DataArray
datastore[trid]["data"][ddi]["adj_out"] = xr.DataArray(
adj_emissions,
coords={"time": prior_emissions.time},
dims=("time", "lev", "lat", "lon"))
else:
logger.critical("")
logger.critical("*"*30)
logger.critical(PROG_NAME+" => ERROR: Unknown value for trid[0]!")
logger.critical(f" trid = {trid}" )
logger.critical("*"*30)
logger.critical("")
raise CifRuntimeError
# end if trid[0]==fluxes
# end for trid
# Some debug statements...
logger.debug("*"*30)
logger.debug("")
#logger.debug("")
#logger.debug("*"*30)
#logger.debug(PROG_NAME+" => Computer says no!")
#logger.debug("*"*30)
#logger.debug("")
#try:
# raise RuntimeError
#except RuntimeError as e:
# #logger.exception("OOPS!")
# logger.critical(e, exc_info=True)
# #raise # => Will display the traceback on screen a second time
# sys.exit() # => Just exit.
## end try
# Return the datastore
return datastore
# end function native2inputs_adj