In [1]:
from gp_emulator import MultivariateEmulator
angles = [ [15, 15, 0] ]
#emulator, samples, validate = create_emulators ( state, [""], angles=angles )
#for i,(s,v,r) in enumerate(angles):
# fname = "%02d_sza_%02d_vza_000_raa" % (s,v)
# emulator[i].dump_emulator(fname)
#emulators = {}
#emulators[(s,v)] = emulator[0]
emulators = {}
for i,(s,v,r) in enumerate(angles):
fname = "%02d_sza_%02d_vza_000_raa.npz" % (s,v)
emulators[(v,s)]= MultivariateEmulator ( dump=fname )
In [4]:
emulators.keys()
Out[4]:
In [5]:
emulator = emulators[(15,15)]
In [38]:
import prosail
import scipy.stats as ss
from gp_emulator import MultivariateEmulator, lhd
def fixnan(x):
'''
the RT model sometimes fails so we interpolate over nan
This method replaces the nans in the vector x by their interpolated values
'''
for i in xrange(x.shape[0]):
sample = x[i]
ww = np.where(np.isnan(sample))
nww = np.where(~np.isnan(sample))
sample[ww] = np.interp(ww[0],nww[0],sample[nww])
x[i] = sample
return x
def create_emulators ( state, fnames, random=False, v_size=300, n_size=175, angles=None, \
vzax = np.arange ( 5, 60, 5 ), szax = np.arange ( 5, 60, 5 ), \
raax = np.arange (-180, 180, 30 ) ):
"""
This function produces the sampling of parameter space based on the
state object that it receives. It then goes to generate the forward
models required, as well as a random sampling to validate against.
TODO completeme!!!!!
"""
distributions = []
for i, (param, minval ) in enumerate ( state.parameter_min.iteritems() ):
maxval = state.bounds[i][1]
minval = state.bounds[i][0]
distributions.append ( ss.uniform( loc=minval, \
scale=(maxval-minval)) )
if random:
samples = []
for d in distributions:
samples.append ( d.rvs( n_size ) )
samples = np.array ( samples ).T
else:
samples = lhd ( dist=distributions, size=n_size )
if angles is None:
angles = np.array ( [angle \
for angle in itertools.product ( szax, vzax, raax )] )
validate = []
for d in distributions:
validate.append ( d.rvs( v_size ) )
validate = np.array ( validate ).T
V = np.zeros ((v_size, 16))
V[:, 9:11] = validate[:, 8:]
V[:, 11] = 0.01
V[:, -1] = 2
S = np.zeros ((n_size, 16))
S[:, 9:11] = samples[:, 8:]
S[:, 11] = 0.01 # hotspot
S[:, -1] = 2
for i,k in enumerate ( state.parameter_max.keys()[:8]):
if state.invtransformation_dict.has_key(k):
S[:, i] = state.invtransformation_dict[k](samples[:, i] )
V[:, i] = state.invtransformation_dict[k](validate[:, i] )
else:
S[:, i] = samples[:, i]
V[:, i] = validate[:, i]
gps = []
for i, angle in enumerate ( angles ):
S[:, 12:15] = angle
V[:, 12:15] = angle
train_brf = np.array ( [ prosail.run_prosail ( *s ) for s in S ] )
validate_brf = np.array ( [ prosail.run_prosail ( *s ) for s in V ] )
train_brf = fixnan ( train_brf )
validate_brf = fixnan ( validate_brf )
emu, rmse = do_mv_emulation ( samples, validate, train_brf, validate_brf )
print "(RMSE:%g)" % ( rmse )
#emu.dump_emulator ( fnames[i] )
gps.append ( emu )
return gps, samples, validate
def do_mv_emulation ( xtrain, xvalidate, train_brf, validate_brf ):
wv = np.arange (400, 2501)
N = xvalidate.shape[0]
emu = MultivariateEmulator (y=xtrain, X=train_brf )
rmse = np.sqrt(np.mean([(validate_brf[i] - \
emu.predict ( np.atleast_2d(xvalidate[i,:]))[0])**2 \
for i in xrange(N)]))
for i in xrange(5):
plt.plot ( wv, validate_brf[i], '-', lw=2.1 )
plt.plot ( wv, emu.predict ( np.atleast_2d(xvalidate[i,:]))[0], '--k', lw=1)
plt.xlim(400,2500)
plt.ylabel(r'$\rho(\Omega_{i}, \Omega_{s})\;[-]$')
plt.xlabel(r'$\lambda\; [nm]$')
plt.title("RMSE: %g" % rmse )
return emu, rmse
In [35]:
from collections import OrderedDict
try:
from operators import *
from state import *
except ImportError:
import sys
sys.path.append ( "../eoldas_ng/")
from operators import *
from state import *
state_config = OrderedDict ()
state_config['n'] = FIXED
state_config['cab'] = VARIABLE
state_config['car'] = FIXED
state_config['cbrown'] = FIXED
state_config['cw'] = FIXED
state_config['cm'] = FIXED
state_config['lai'] = VARIABLE
state_config['ala'] = FIXED
state_config['bsoil'] = FIXED
state_config['psoil'] = FIXED
# Now define the default values
default_par = OrderedDict ()
default_par['n'] = 1.5
default_par['cab'] = 40.
default_par['car'] = 10.
default_par['cbrown'] = 0.01
default_par['cw'] = 0.018 # Say?
default_par['cm'] = 0.0065 # Say?
default_par['lai'] = 2
default_par['ala'] = 45.
default_par['bsoil'] = 1.
default_par['psoil'] = 0.1
parameter_min = OrderedDict()
parameter_max = OrderedDict()
min_vals = [ 0.8, 0.2, 0.0, 0.0, 0.0043, 0.0017,0.01, 40, 0., 0., 0.]
max_vals = [2.5, 77., 5., 1., 0.0713, 0.0331, 8., 50., 2., 1.]
for i, param in enumerate ( state_config.keys() ):
parameter_min[param] = min_vals[i]
parameter_max[param] = max_vals[i]
# Define parameter transformations
transformations = {
'lai': lambda x: np.exp ( -x/2. ), \
'cab': lambda x: np.exp ( -x/100. ), \
'car': lambda x: np.exp ( -x/100. ), \
'cw': lambda x: np.exp ( -50.*x ), \
'cm': lambda x: np.exp ( -100.*x ), \
'ala': lambda x: x/90. }
inv_transformations = {
'lai': lambda x: -2*np.log ( x ), \
'cab': lambda x: -100*np.log ( x ), \
'car': lambda x: -100*np.log( x ), \
'cw': lambda x: (-1/50.)*np.log ( x ), \
'cm': lambda x: (-1/100.)*np.log ( x ), \
'ala': lambda x: 90.*x }
# Define the state grid. In time in this case
state_grid = np.arange ( 1, 366 )
# Define the state
# L'etat, c'est moi
state = State ( state_config, state_grid, default_par, \
parameter_min, parameter_max )
# Set the transformations
state.set_transformations ( transformations, inv_transformations )
In [36]:
import json
import matplotlib
import urllib2
url = "http://github.com/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers/raw/master/styles/bmh_matplotlibrc.json"
s = json.load( urllib2.urlopen(url))
matplotlib.rcParams.update(s)
matplotlib.rcParams['text.usetex'] = True
In [39]:
emulator, samples, validate = create_emulators ( state, [""], angles=angles )
plt.savefig("/tmp/emulator_perf.pdf")
In [ ]: