Make sure to import GaussianProcess and prosail_fortran from:
from pylab import *
from prosail_fortran import run_prosail
import numpy as np
from GaussianProcess import GaussianProcess
def xify(p,params = ['n','xcab','xcar','cbrown','xcw','xcm','xala','bsoil','psoil','hspot','xlai']):
'''
Fix the state key names by putting x on where appropriate
'''
retval = {}
for i in p.keys():
if 'x' + i in params:
retval['x' + i] = p[i]
else:
retval[i] = p[i]
return retval
def invtransform(params,logvals = {'ala': 90., 'lai':-2.0, 'cab':-100., 'car':-100., 'cw':-1/50., 'cm':-1./100}):
retval = {}
for i in xify(params).keys():
if i[0] != 'x':
retval[i] = params[i]
else:
rest = i[1:]
try:
if logvals[rest] < 0:
retval[rest] = logvals[rest] * np.log(params[i])
else:
retval[rest] = params[i] * logvals[rest]
except:
retval[rest] = params[i]
return retval
def transform(params,logvals = {'ala': 90., 'lai':-2.0, 'cab':-100., 'car':-100., 'cw':-1/50., 'cm':-1./100}):
retval = {}
for i in xify(params).keys():
if i[0] != 'x':
retval[i] = params[i]
else:
rest = i[1:]
try:
if logvals[rest] < 0:
retval[i] = np.exp(params[rest] / logvals[rest])
else:
retval[i] = params[rest] / logvals[rest]
except:
retval[i] = params[i]
return retval
def limits():
'''
Set limits for parameters
'''
params = ['n','xcab','xcar','cbrown','xcw','xcm','xala','bsoil','psoil','hspot','xlai']
pmin = {}
pmax = {}
# set up defaults
for i in params:
pmin[i] = 0.001
pmax[i] = 1. - 0.001
# now specific limits
# These from Feret et al.
pmin['xlai'] = transform({'lai':15.})['xlai']
pmin['xcab'] = transform({'cab':0.2})['xcab']
pmax['xcab'] = transform({'cab':76.8})['xcab']
pmin['xcar'] = transform({'car':25.3})['xcar']
pmax['xcar'] = transform({'car':0.})['xcar']
pmin['cbrown'] = 1.
pmax['cbrown'] = 0.
pmin['xcm'] = transform({'cm':0.0017})['xcm']
pmax['xcm'] = transform({'cm':0.0331})['xcm']
pmin['xcw'] = transform({'cw':0.0043})['xcw']
pmax['xcw'] = transform({'cw':0.0713})['xcw']
pmin['n'] = 0.8
pmax['n'] = 2.5
pmin['bsoil'] = 0.0
pmax['bsoil'] = 2.0
pmin['psoil'] = 0.0
pmax['psoil'] = 1.0
pmin['xala'] = transform({'ala':90.})['xala']
pmax['xala'] = transform({'ala':0.})['xala']
for i in params:
if pmin[i] > pmax[i]:
tmp = pmin[i]
pmin[i] = pmax[i]
pmax[i] = tmp
return (pmin,pmax)
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 addDict(input,new):
'''
add dictionary items
'''
for i in new.keys():
try:
input[i] = np.hstack((np.atleast_1d(input[i]),np.atleast_1d(new[i])))
except:
input[i] = np.array([new[i]]).flatten()
return input
def unpack(params):
'''
Input a dictionary and output keys and array
'''
inputs = []
keys = np.sort(params.keys())
for i,k in enumerate(keys):
inputs.append(params[k])
inputs=np.array(inputs).T
return inputs,keys
def pack(inputs,keys):
'''
Input keys and array and output dict
'''
params = {}
for i,k in enumerate(keys):
params[k] = inputs[i]
return params
def sampler(pmin,pmax):
'''
Random sample
'''
params = {}
for i in pmin.keys():
params.update(invtransform({i:pmin[i] + np.random.rand()*(pmax[i]-pmin[i])}))
return params
def samples(n=2):
'''
Random samples over parameter space
'''
(pmin,pmax) = limits()
s = {}
for i in xrange(n):
s = addDict(s,sampler(pmin,pmax))
return s
def reconstruct(gps,SB,inputs):
'''
spectral reconstruction from GPs and spectral basis functions
'''
fwd = 0
inputs = np.atleast_2d(np.array(inputs))
for i in xrange(SB.shape[0]):
pred_mu, pred_var, grad = gps[i].predict(inputs)
fwd += np.matrix(pred_mu).T * np.matrix(SB[i])
return np.array(fwd)
def reconstructD(gps,SB,inputs):
'''
spectral reconstruction from GPs and spectral basis functions
'''
fwd = 0
inputs = np.atleast_2d(np.array(inputs))
deriv = np.zeros((inputs.shape[1],SB[0].shape[0]))
derivt = np.zeros((inputs.shape[1],SB[0].shape[0]))
for i in xrange(SB.shape[0]):
pred_mu, pred_var, grad = gps[i].predict(inputs)
fwd += np.matrix(pred_mu).T * np.matrix(SB[i])
deriv += np.matrix(grad).T * np.matrix(SB[i])
# check code , but its the same
#for j in xrange(inputs.shape[1]):
# for i in xrange(SB.shape[0]):
# pred_mu, pred_var, grad = gps[i].predict(inputs)
# derivt[j,:] += grad[:,j] * SB[i]
return np.array(fwd)[0],deriv
def brdfModel(s,vza,sza,raa):
'''
Run the full BRDF model for parameters in s
'''
try:
brdf = []
for i in xrange(len(s['n'])):
brdf.append(run_prosail(s['n'][i],s['cab'][i],s['car'][i],s['cbrown'][i],s['cw'][i],s['cm'][i],\
s['lai'][i],s['ala'][i],0.0,s['bsoil'][i],s['psoil'][i],s['hspot'][i],\
vza[i],sza[i],raa[i],2))
except:
brdf = []
for i in xrange(len(s['n'])):
brdf.append(run_prosail(s['n'][i],s['cab'][i],s['car'][i],s['cbrown'][i],s['cw'][i],s['cm'][i],\
s['lai'][i],s['ala'][i],0.0,s['bsoil'][i],s['psoil'][i],s['hspot'][i],\
vza,sza,raa,2))
return (vza,sza,raa),transform(s),fixnan(np.array(brdf))
def lut(n=2,vza=0.,sza=45.,raa=0.,brdf=None,bands=None):
'''
Generate a random LUT for given viewing and illumination geometry
'''
# get the parameter limits
(pmin,pmax) = limits()
# generate some random samples over the parameter space
s = samples(n=n)
return brdfModel(s,vza,sza,raa)
def funcd(params,brf,SB,gps):
fwd,derivs = reconstructD(gps,SB,params)
obs = brf
d = (fwd - obs)
Jprime = np.matrix(d) * np.matrix(derivs).T
e = np.sum(d*d)
return 0.5*e,Jprime.T
from scipy.optimize import fmin_l_bfgs_b
# function to do the minimisation in EOF space
def funcd2(params,projectedBrf,gps):
params = np.atleast_2d(np.array(params))
e = 0.
Jprime = np.zeros_like(params)
for i in xrange(len(gps)):
fwd,sd,deriv = gps[i].predict(params)
d = (fwd[0] - projectedBrf[i]) # /sd[0]
Jprime += d*deriv
e += d*d
return 0.5*e,Jprime
def solver(i,noshow=True):
params = fmin_l_bfgs_b(funcd,x,bounds=b,args=(test_brf[i],SB,gps),iprint=0)[0]
fwd = reconstruct(gps,SB,params)[0]
plt.figure(i)
plt.clf()
plt.plot(wl,test_brf[i],'k-')
plt.plot(wl,fwd,'r+')
if not noshow:
plt.show()
print 'parameters1: ',invtransform(pack(params,keys))
if not noshow:
print 'true : ',invtransform(pack(inputs_v[i],keys))
def solver2(i,noshow=True):
params = fmin_l_bfgs_b(funcd2,x,bounds=b,args=(test_SB[i],gps),iprint=0)[0]
fwd = reconstruct(gps,SB,params)[0]
if not noshow:
plt.figure(i)
plt.clf()
plt.plot(wl,test_brf[i],'k-')
plt.plot(wl,fwd,'g+')
plt.show()
print 'parameters2: ',invtransform(pack(params,keys))
print 'true : ',invtransform(pack(inputs_v[i],keys))
import pickle
# full sprectum
bands = None
thresh = 0.98
n = 150
lutfile = 'lutTest.dat.npz'
pickleFile = 'lutTest.dat.pkl'
vza = 0.
sza = 45.
raa = 0.
try:
gps = pickle.load(open(pickleFile,'rb'))
npzfile = np.load(lutfile)
angles,train_params,train_brf = npzfile['angles'],npzfile['train_params'],npzfile['train_brf']
theta_min0,theta_min1 = npzfile['theta_min0'],npzfile['theta_min1']
s,SB,train_SB = npzfile['s'],npzfile['SB'],npzfile['train_SB']
readLUT = True
print 'Read lutfile %s and %s'%(lutfile,pickleFile)
except:
readLUT = False
print 'Failed to read lutfile %s'%lutfile
angles,train_params,train_brf = lut(n,bands=bands,vza=vza,sza=sza,raa=raa)
# spectral reduction
U, s, V = np.linalg.svd(train_brf, full_matrices=True)
wt = s.cumsum()/s.sum()
ww = np.where(wt<=thresh)
SB = V[ww[0],:]
# now project the training data onto these axes
train_SB = np.dot(train_brf,SB.T).T
# unpack params
inputs_t,keys = unpack(train_params)
# set some arrays to collect information
gps = ['']*train_SB.shape[0]
theta_min0 = ['']*train_SB.shape[0]
theta_min1 = ['']*train_SB.shape[0]
# loop over each dimenstion and train a gp
for i in xrange(train_SB.shape[0]):
print 'doing PC %d'%i
yields_t = train_SB[i]
gps[i] = GaussianProcess ( inputs_t, yields_t )
theta_min0[i], theta_min1[i] = gps[i].learn_hyperparameters (n_tries=2)
pickle.dump(gps,open(pickleFile,'wb'))
np.savez(lutfile,angles=angles,train_params=train_params,\
train_brf=train_brf,\
theta_min0=theta_min0,theta_min1=theta_min1,\
s=s,SB=SB,train_SB=train_SB)
wl = np.arange(400,2501).astype(float)
Read lutfile lutTest.dat.npz and lutTest.dat.pkl
plt.clf()
for i in xrange(train_SB.shape[0]):
plt.plot(wl,SB[i],label='PC %d (%.2f)'%(i,(s/s.sum())[i]))
plt.legend()
plt.xlim(wl.min(),wl.max())
plt.show()
# now look at e.g. reflectance as a function of LAI
(pmin,pmax) = limits()
xlais = np.arange(pmin['xlai'],pmax['xlai'],0.1)
ones = np.ones_like(xlais)
params = {}
for i in pmin.keys():
params.update(invtransform({i:ones*(pmin[i]+pmax[i])*0.5}))
params.update(invtransform({'xlai':xlais}))
#import pdb;pdb.set_trace()
# now run true model
brdf = brdfModel(params,vza,sza,raa)[-1]
# now emulate (dont forget to transform!)
inputs,keys = unpack(transform(params))
brdfe = reconstruct(gps,SB,inputs)
plt.clf()
plt.plot(np.array([wl]*10).T,brdf.T,'k-')
plt.plot(np.array([wl]*10).T,brdfe.T,'r--')
plt.xlim(wl[0],wl[-1])
plt.xlabel('wavelength (nm)')
plt.ylabel('reflectance')
plt.title('true (black) and emulated (red) reflectance for varying LAI')
plt.show()
##### We now generate an independent dataset and use that for testing the GP
angles,test_params,test_brf = lut(n,bands=bands)
test_SB = np.dot(test_brf,SB.T).T
# unpack params
inputs_v,keys = unpack(test_params)
plt.clf()
# loop over each dimenstion and test a gp
for i in xrange(train_SB.shape[0]):
# This is the set of true projected parameters
yields_v = test_SB[i]
# This is what we predict from the emulator
pred_mu, pred_var, deriv = gps[i].predict ( inputs_v )
plt.plot(yields_v,pred_mu,'+',label='PC %d'%i)
plt.legend()
plt.show()
# Now look at some example reconstructions
fwd = reconstruct(gps,SB,inputs_v)
i = 0
true = test_brf[i]
em = fwd[i]
plt.clf()
plt.plot(wl,true,'k-')
plt.plot(wl,em,'r+')
plt.xlabel('wavelength (nm)')
plt.ylabel('reflectance')
plt.show()
i = 50
true = test_brf[i]
em = fwd[i]
plt.clf()
plt.plot(wl,true,'k-')
plt.plot(wl,em,'r+')
plt.xlabel('wavelength (nm)')
plt.ylabel('reflectance')
plt.show()
i = 100
true = test_brf[i]
em = fwd[i]
plt.clf()
plt.plot(wl,true,'k-')
plt.plot(wl,em,'r+')
plt.xlabel('wavelength (nm)')
plt.ylabel('reflectance')
plt.show()
# check derivatives on training set of all PCs
inputs_t,keys = unpack(train_params.tolist())
xf = lambda x,j:gps[j].predict ( np.atleast_2d(x) )[0]
import scipy.optimize
for j in xrange(len(gps)):
trueDerive = []
for i in xrange(inputs_t.shape[0]):
trueDerive.append(scipy.optimize.approx_fprime(inputs_t[i],xf,1e-10,(j)))
pred_mu, pred_var, deriv = gps[j].predict ( np.atleast_2d(inputs_t) )
# plot
plt.scatter(np.array(trueDerive).flatten(),deriv.flatten())
# check derivatives on test set of all PCs
for j in xrange(len(gps)):
trueDerive = []
for i in xrange(inputs_v.shape[0]):
trueDerive.append(scipy.optimize.approx_fprime(inputs_v[i],xf,1e-10,(j)))
pred_mu, pred_var, deriv = gps[j].predict ( np.atleast_2d(inputs_v) )
# plot
plt.scatter(np.array(trueDerive).flatten(),deriv.flatten())
# We now attempt to use the emulator to solve inverse problems (ultimately to allow fast data assimilation).
from scipy.optimize import fmin_l_bfgs_b
# form the state vector
(pmin,pmax) = limits()
keys = np.sort(pmin.keys())
ppmin = unpack(pmin)[0]
ppmax = unpack(pmax)[0]
b = [(ppmin[i],ppmax[i]) for i in xrange(len(keys))]
x = []
for k in keys:
x.append((pmin[k]+pmax[k])*0.5)
# project the test_brf data
test_SB = np.dot(test_brf,SB.T)
solver(0)
solver2(0)
parameters1: {'bsoil': 2.0, 'lai': 1.0288464414236769, 'cm': 0.0031670864300728618, 'car': 5.7927940429927025, 'n': 2.3018279471599672, 'cbrown': 0.0, 'psoil': 0.30291664950480746, 'cab': 17.433484976099766, 'ala': 19.282968992297693, 'cw': 0.0053248298367763727, 'hspot': 0.001}
parameters2: {'bsoil': 2.0, 'lai': 1.0289427087074099, 'cm': 0.0031684180346170376, 'car': 5.7959665441108328, 'n': 2.3015616933495635, 'cbrown': 0.0, 'psoil': 0.30333894461870015, 'cab': 17.432121475503337, 'ala': 19.28509961974304, 'cw': 0.0053240434061770569, 'hspot': 0.001} true : {'bsoil': 0.38463980331250447, 'lai': 1.3476339679025948, 'cm': 0.0018694991805296669, 'car': 5.0461058317899825, 'n': 1.8086975152880664, 'cbrown': 0.10538012583267453, 'psoil': 0.39197245292595684, 'cab': 15.890878922221249, 'ala': 38.928290861773945, 'cw': 0.005060167923528305, 'hspot': 0.35844017613866785}
solver(25)
solver2(25)
parameters1: {'bsoil': 0.97217154312287568, 'lai': 2.2069833210503735, 'cm': 0.011324946622961208, 'car': 25.300000000000001, 'n': 1.4057256465574939, 'cbrown': 0.98362925612801955, 'psoil': 0.49760623182034142, 'cab': 76.799999999999997, 'ala': 72.408641106983154, 'cw': 0.042471411951321071, 'hspot': 0.35266196411618583}
parameters2: {'bsoil': 0.97264612023037844, 'lai': 2.2051100238656676, 'cm': 0.011312888006075853, 'car': 25.300000000000001, 'n': 1.4093879133029454, 'cbrown': 0.98515896247270873, 'psoil': 0.49761773656842345, 'cab': 76.799999999999997, 'ala': 72.404073175563497, 'cw': 0.042508636445312921, 'hspot': 0.35017754298293852} true : {'bsoil': 0.25487836292536348, 'lai': 3.9627723199185283, 'cm': 0.017327762664483434, 'car': 15.388576663653247, 'n': 1.3272295140676365, 'cbrown': 0.68349705693840679, 'psoil': 0.15716298813882867, 'cab': 54.333613321248407, 'ala': 77.716271537510281, 'cw': 0.03206577152357229, 'hspot': 0.17619898982445925}
solver(50)
solver2(50)
parameters1: {'bsoil': 1.0054504416517878, 'lai': 5.9168285182593108, 'cm': 0.0016999999999999999, 'car': 12.948806675871637, 'n': 1.89674035103391, 'cbrown': 0.77249768507157135, 'psoil': 0.0, 'cab': 11.801865087990372, 'ala': 4.7562011744021513, 'cw': 0.0043, 'hspot': 0.089699283007800015}
parameters2: {'bsoil': 1.0131976031658587, 'lai': 5.9176865227451332, 'cm': 0.0016999999999999999, 'car': 13.006412146290069, 'n': 1.896656679139388, 'cbrown': 0.77249268937948146, 'psoil': 0.49296194130045989, 'cab': 11.801732091172669, 'ala': 4.7456956113471591, 'cw': 0.0043, 'hspot': 0.089651594507894422} true : {'bsoil': 0.51189615973810287, 'lai': 5.0832725379024639, 'cm': 0.0017861463740213518, 'car': 9.1133616003277425, 'n': 2.2197603658359366, 'cbrown': 0.75116203429465855, 'psoil': 0.31066402708010743, 'cab': 15.688470356749168, 'ala': 42.410592628282764, 'cw': 0.0050291627104636803, 'hspot': 0.77396264559812589}
solver(75)
solver2(75)
parameters1: {'bsoil': 1.0029812878205573, 'lai': 1.6010272611799012, 'cm': 0.012823468392655787, 'car': 25.300000000000001, 'n': 0.99111517038066521, 'cbrown': 1.0, 'psoil': 0.49552221075552938, 'cab': 16.743778815193263, 'ala': 89.413127252244948, 'cw': 0.052509469699587745, 'hspot': 0.29562942635741307}
parameters2: {'bsoil': 1.0029820326427359, 'lai': 1.6010048244370632, 'cm': 0.012824273258428614, 'car': 25.300000000000001, 'n': 0.99109798356424939, 'cbrown': 1.0, 'psoil': 0.49552180404812601, 'cab': 16.743801178709493, 'ala': 89.413208961658526, 'cw': 0.052509369226395521, 'hspot': 0.29563276046450326} true : {'bsoil': 0.77275384433249816, 'lai': 0.890427582841266, 'cm': 0.0056417932967290085, 'car': 9.8050112678704746, 'n': 2.3741778873167343, 'cbrown': 0.89837380876364747, 'psoil': 0.95861246365923425, 'cab': 40.735068228418086, 'ala': 87.678028671371152, 'cw': 0.010467970536657878, 'hspot': 0.20863625028845817}
solver(100)
solver2(100)
parameters1: {'bsoil': 1.0586424655022717, 'lai': 0.41023410712496428, 'cm': 0.0062442099760283863, 'car': 7.5650698517916757, 'n': 1.8859641680740722, 'cbrown': 0.40301561507121642, 'psoil': 0.50058524617816758, 'cab': 27.017551798987753, 'ala': 21.877126526793031, 'cw': 0.042443775617885988, 'hspot': 0.94100289723669839}
parameters2: {'bsoil': 1.0577437134336234, 'lai': 0.41593972010550473, 'cm': 0.0063818551119253228, 'car': 7.5015910699754507, 'n': 1.8884284712351815, 'cbrown': 0.39753652105903498, 'psoil': 0.50057998842365614, 'cab': 27.071074395566104, 'ala': 23.035798404389649, 'cw': 0.042180607576328473, 'hspot': 0.93795169929015987} true : {'bsoil': 0.027293706911639015, 'lai': 1.0504474180392158, 'cm': 0.0059607066850681581, 'car': 2.228060582783757, 'n': 1.579203466293738, 'cbrown': 0.18124965581191665, 'psoil': 0.15295568171553908, 'cab': 25.031166442921666, 'ala': 62.259709487107884, 'cw': 0.031637813894730978, 'hspot': 0.38813771702386185}
%timeit solver2(100)
parameters2: {'bsoil': 1.0577437134336234, 'lai': 0.41593972010550473, 'cm': 0.0063818551119253228, 'car': 7.5015910699754507, 'n': 1.8884284712351815, 'cbrown': 0.39753652105903498, 'psoil': 0.50057998842365614, 'cab': 27.071074395566104, 'ala': 23.035798404389649, 'cw': 0.042180607576328473, 'hspot': 0.93795169929015987} true : {'bsoil': 0.027293706911639015, 'lai': 1.0504474180392158, 'cm': 0.0059607066850681581, 'car': 2.228060582783757, 'n': 1.579203466293738, 'cbrown': 0.18124965581191665, 'psoil': 0.15295568171553908, 'cab': 25.031166442921666, 'ala': 62.259709487107884, 'cw': 0.031637813894730978, 'hspot': 0.38813771702386185}
parameters2: {'bsoil': 1.0577437134336234, 'lai': 0.41593972010550473, 'cm': 0.0063818551119253228, 'car': 7.5015910699754507, 'n': 1.8884284712351815, 'cbrown': 0.39753652105903498, 'psoil': 0.50057998842365614, 'cab': 27.071074395566104, 'ala': 23.035798404389649, 'cw': 0.042180607576328473, 'hspot': 0.93795169929015987} true : {'bsoil': 0.027293706911639015, 'lai': 1.0504474180392158, 'cm': 0.0059607066850681581, 'car': 2.228060582783757, 'n': 1.579203466293738, 'cbrown': 0.18124965581191665, 'psoil': 0.15295568171553908, 'cab': 25.031166442921666, 'ala': 62.259709487107884, 'cw': 0.031637813894730978, 'hspot': 0.38813771702386185}
parameters2: {'bsoil': 1.0577437134336234, 'lai': 0.41593972010550473, 'cm': 0.0063818551119253228, 'car': 7.5015910699754507, 'n': 1.8884284712351815, 'cbrown': 0.39753652105903498, 'psoil': 0.50057998842365614, 'cab': 27.071074395566104, 'ala': 23.035798404389649, 'cw': 0.042180607576328473, 'hspot': 0.93795169929015987} true : {'bsoil': 0.027293706911639015, 'lai': 1.0504474180392158, 'cm': 0.0059607066850681581, 'car': 2.228060582783757, 'n': 1.579203466293738, 'cbrown': 0.18124965581191665, 'psoil': 0.15295568171553908, 'cab': 25.031166442921666, 'ala': 62.259709487107884, 'cw': 0.031637813894730978, 'hspot': 0.38813771702386185}
parameters2: {'bsoil': 1.0577437134336234, 'lai': 0.41593972010550473, 'cm': 0.0063818551119253228, 'car': 7.5015910699754507, 'n': 1.8884284712351815, 'cbrown': 0.39753652105903498, 'psoil': 0.50057998842365614, 'cab': 27.071074395566104, 'ala': 23.035798404389649, 'cw': 0.042180607576328473, 'hspot': 0.93795169929015987} true : {'bsoil': 0.027293706911639015, 'lai': 1.0504474180392158, 'cm': 0.0059607066850681581, 'car': 2.228060582783757, 'n': 1.579203466293738, 'cbrown': 0.18124965581191665, 'psoil': 0.15295568171553908, 'cab': 25.031166442921666, 'ala': 62.259709487107884, 'cw': 0.031637813894730978, 'hspot': 0.38813771702386185} 1 loops, best of 3: 1.02 s per loop