diff --git a/models.py b/models.py index 62c584d910e7403c826be3f85dc28cb8466bce74..9ffdbc3f5742dc307ba03d7bedaf4d98995a6496 100644 --- a/models.py +++ b/models.py @@ -2,13 +2,55 @@ # -*- coding: utf-8 -*- import numpy as np from scipy.interpolate import CubicSpline +from scipy.special import erfinv,erf verbose = True +def transform_none(x,kind): + if kind=='dir': + return x + if kind=='inv': + return x + +def transform_log(x,kind): + if kind=='dir': + return np.log(x) + if kind=='inv': + return np.exp(x) + +def transform_pit(x,kind): + # mu and sigma could be guessed from the data, but... + a = 0.5 + b = 4.5 + mu = 0.5*(a+b) + sigma = (b-a)/12**0.5 + # takes a uniform, returns a normal + if kind=='dir': + p = (x-a)/(b-a) + xt = mu+sigma*2**0.5*erfinv(2*p-1) + return xt + # takes a normal), returns a uniform + if kind=='inv': + p = 0.5*(1+erf((x-mu)/sigma/2**0.5)) + xt = a+p*(b-a) + return xt + class CBL: def __init__(self,settings): for k, v in settings.items(): setattr(self, k, v) + self.parameter_true = np.array([self.parameter_true]) + self.parameter_ensemble_min = np.array([self.parameter_ensemble_min]) + self.parameter_ensemble_max = np.array([self.parameter_ensemble_max]) + self.parameter_inflation_rtps_alpha = np.array([self.parameter_inflation_rtps_alpha]) + if self.parameter_transform == 'pit': + self.parameter_transform = [ transform_pit ] + elif self.parameter_transform == 'log': + self.parameter_transform = [ transform_log ] + elif self.parameter_transform == 'none': + self.parameter_transform = [ transform_none ] + if not hasattr(self,'do_parameter_estimation'): + self.do_parameter_estimation == False def initialize(self,argin,coordinates=None):