Add files via upload

Emergent constraint code needed to generate figures for the following paper: 
"Combining local model calibration with the emergent constraint approach to reduce uncertainty in the tropical land carbon cycle feedback"

Emergent_constraint/Booth_2012.csv

@@ -0,0 +1 @@
+x,y28.34,748.60729.012,852.62827.031,824.98429.684,726.90527.689,623.09430.325,580.74331.015,596.47631.691,604.74132.355,622.29733.676,604.43837.039,540.51135.022,588.28433.012,499.40334.358,513.23535.679,399.05436.371,379.53

Emergent_constraint/EC_pdf_basic.py

@@ -0,0 +1,167 @@
+def EC_pdf_basic(x,y,x_obs,dx_obs,xtitle,ytitle) :
+    import numpy as np
+    from pylab import plot, legend, xlabel, ylabel, figure, savefig, title
+    from pylab import mean, exp, sqrt, arange, ylim
+    from lin_reg_DP import lin_reg_DP
+    import matplotlib.mlab as mlab
+    import matplotlib.pyplot as plt
+    import matplotlib
+    import matplotlib.cm as cm
+    from plot_scat_basic import plot_scat_basic
+    import seaborn as sns
+    from scipy.stats import norm
+
+# Calculate mean and stdev of (equal model weight) prior
+    mn_pr=mean(y)
+    std_pr=np.std(y)
+    print(mn_pr)
+    print(std_pr)
+
+    mn=x_obs
+    std=dx_obs
+
+# Calculate best-fit straight-line between x & y
+    yf,a,b,da,db,xfit,yfit,yband=lin_reg_DP(x,y)
+
+# Plot contours of probability fro linear regression
+    mult=[-1,1]
+    jconts=len(mult)
+    figure()
+    plot_scat_basic(x,y,xtitle,ytitle)
+    plot(xfit,yfit,'k-.')
+    plot(32,600,'b')
+    for j in range (0,jconts):
+      y1=yfit+mult[j]*yband
+      plot(xfit,y1,'k--')
+    pass
+
+# Plot observational constraint
+    xbest=x_obs
+    xlo=x_obs-dx_obs
+    xhi=x_obs+dx_obs
+    plot([xbest,xbest],[300,900],'b-.')
+    plot([xlo,xlo],[300,900],'b--',label='adJULES Constraint')
+    plot([xhi,xhi],[300,900],'b--')
+    legend(loc=3,prop={'size':11})
+    savefig("Emergent_Relationship_Contours.pdf")
+
+# Calculate PDF for IAV constraints
+    x2=xfit
+    nfitx=len(xfit)
+    dx=x2[1]-x2[0]
+    Px=x2
+    Pi=3.142
+    Px=1/sqrt(2*Pi*std**2) * exp(-((x2-mn)/(sqrt(2)*std))**2)
+
+    miny=0.1*min(y)
+    maxy=1.1*max(y)
+    mfity=2000
+    dy=(maxy-miny)/float(mfity)
+    #y2=miny+dy*range(0,mfity)
+    y2=miny+[dy* i for i in range(0,mfity)]
+
+# Calculate "prior"
+    Py_pr=y2
+    Py_pr=1/sqrt(2*Pi*std_pr**2)*exp(-((y2-mn_pr)/(sqrt(2)*std_pr))**2)
+
+# Calculate contours of probability in (x,y) space
+    Pxy=np.zeros((nfitx,mfity))
+    Pyx=np.zeros((mfity,nfitx))
+    Py=np.zeros(mfity)
+    Py_norm=0.0
+    for m in range(0, mfity):
+        Py[m]=0.0
+        for n in range(0,nfitx):
+            Py_given_x=1/sqrt(2*Pi*yband[n]**2) \
+               * exp(-((y2[m]-yfit[n])/(sqrt(2)*yband[n]))**2)
+            Pxy[n,m]=Px[n]*Py_given_x
+            Pyx[m,n]=Pxy[n,m]
+# Integrate over x to get Py
+            Py[m]=Py[m]+Pxy[n,m]*dx
+        pass
+        Py_norm=Py_norm+Py[m]*dy
+    pass
+
+    print('NORMALISATION REQUIRED FOR POSTERIOR = ',Py_norm)
+# Normalise Py
+    for m in range(0, mfity):
+        Py[m]=Py[m]/Py_norm
+    pass
+
+ #   figure()
+ #   plot(x,y,'ro')
+    z=100*Pyx
+    CS = plt.contour(x2,y2,z)
+    ylim(300,900) #ylim(min(y),max(y))
+    plt.clabel(CS)
+    savefig("Emergent_Relationship_contours_full.pdf")
+
+# Plot PDF
+    figure()
+    plot(y2,Py,'k-',label='Emergent Constraint',linewidth=3)
+    xlabel(ytitle,size=14)
+    ylabel('Probablity Density',size=14)  
+    dum=np.argmax(Py)
+    ybest=y2[dum]
+    a=np.column_stack((y2,Py))
+    print(np.std(a))
+    print(ybest-np.std(a))
+    #plot(y2,norm.pdf(y2,ybest,91))
+
+    dum_pr=np.argmax(Py_pr)
+    ybest_pr=y2[dum_pr]
+    #binny=min(y2)+(max(y2)-min(y2))*arange(11)/10.0
+    binny=[150,250,350,450,550,650,750,850,950]
+    binny=[100,200,300,400,500,600,700,800,900]
+    n, bins, patches = plt.hist(y, bins=binny,\
+                        normed=1, facecolor='orange',label='Equal-weight prior',edgecolor = "none")
+    #sns.kdeplot(np.array(y),linewidth=3)
+    legend(prop={'size':11})
+    title("(a) PDF of Emergent Constraint",size=15,loc='left')
+    savefig("Emergent_Constraint_PDF.pdf")  
+
+# Calculate CDFs
+    CDF=np.zeros(mfity)
+    CDF_pr=np.zeros(mfity)
+    CDF[0]=Py[0]*dy
+    CDF_pr[0]=Py_pr[0]*dy
+    for m in range(1,mfity):
+        CDF[m]=CDF[m-1]+Py[m]*dy
+        CDF_pr[m]=CDF_pr[m-1]+Py_pr[m]*dy
+    pass
+
+
+# Plot CDF   
+    figure()
+    plot(y2,CDF,'k-',label='Emergent Constraint',linewidth=3)
+
+# Find 95% confidence limits    
+    dum_up=CDF-0.975
+    dum_975=dum_up**2
+    dum_lo=CDF-0.025
+    dum_025=dum_lo**2     
+    val, n_lo = min((val, idx) for (idx, val) in enumerate(dum_025))
+    val, n_hi = min((val, idx) for (idx, val) in enumerate(dum_975))
+    val, n_best = max((val, idx) for (idx, val) in enumerate(Py))
+    y_best=y2[n_best]
+    y_lo=y2[n_lo]
+    y_hi=y2[n_hi]
+
+    plot([min(y2),max(y2)],[0.025,0.025],'k-.')
+    plot([min(y2),max(y2)],[0.975,0.975],'k-.')
+    xlabel(ytitle,size=14)
+    ylabel('CDF (Probability<x)',size=14)
+    n, bins, patches = plt.hist(y, bins=binny,cumulative=1,\
+                                normed=1, facecolor='orange',\
+                                label='Equal-weight prior',edgecolor = "none")
+    legend(prop={'size':11})
+    title("(b) CDF of Emergent Constraint",size=15,loc='left')
+    savefig("Emergent_Constraint_CDF.pdf")
+
+    print('EC on Y             = ',y_best,\
+                          ',[',y_lo,'-',y_hi,']')
+
+
+    return ybest, ybest_pr;
+
+

Emergent_constraint/FigurePanelF_2012paper.txt

@@ -0,0 +1,17 @@
+      Topt        CO2 change (2100-1900)
+      31.0000      589.773
+      33.6667      597.978
+      36.3333      373.826
+      35.6667      392.767
+      27.6667      616.082
+      31.6667      599.082
+      27.0000      817.798
+      29.0000      845.664
+      29.6667      720.893
+      33.0000      492.605
+      32.3333      616.539
+      37.0000      534.170
+      34.3333      506.310
+      28.3333      742.913
+      30.3333      573.196
+      35.0000      581.707

Emergent_constraint/MAIN_EC_basic.py

@@ -0,0 +1,67 @@
+#!/usr/bin/env python
+import numpy as np
+from pylab import plot, show, bar, legend, colors, axes, xlabel, ylabel, hist
+from pylab import title, savefig, axis, figure, semilogx, mean, exp, sqrt
+from pylab import log, arctan
+from scipy import stats
+from plot_scat_basic import plot_scat_basic
+from EC_pdf_basic import EC_pdf_basic
+
+# Read x and y variables
+fname='Booth_2012.csv'
+xy_data=np.loadtxt(fname, \
+                 dtype={'names': ('xcol','ycol'),'formats' : ('f8','f8')},\
+                 comments='#',delimiter=',', converters=None, skiprows=1, \
+                 usecols=None, unpack=False, ndmin=0)
+xlab='T$_\mathregular{opt}$ (K)'
+ylab='CO$_2$ Change, 2100-1900 (ppmv)'
+x=xy_data['xcol']
+y=xy_data['ycol']
+npts=len(x)
+
+# Input obsrvational estimate of x
+x_obs=35.392 #35
+dx_obs=0.922
+
+#x_obs=43.027
+#dx_obs=3.376
+
+#x_obs=32
+#dx_obs = sqrt(25/sqrt(0.5))
+
+#x_obs=33.649
+#dx_obs=0.628
+
+#x_obs=37.17861#
+#dx_obs= 2.897104
+
+#x_obs=37.35398
+#dx_obs=1.24051
+# Plot Emergent Relationship
+figure()
+plot_scat_basic(x,y,xlab,ylab)
+deg=1
+p=np.polyfit(x,y,deg)
+yfit=np.zeros(npts)
+for n in range(0,npts):
+    yfit[n]=0.0
+    for j in range(0,deg+1):
+        yfit[n]=yfit[n]+p[deg-j]*x[n]**j
+    pass
+pass
+plot(x,yfit,'k-')
+xbest=x_obs
+xlo=x_obs-dx_obs
+xhi=x_obs+dx_obs
+plot([xbest,xbest],[min(y),max(y)],'b-.')
+plot([xlo,xlo],[min(y),max(y)],'b--',label='ADJULES Constraint')
+plot([xhi,xhi],[min(y),max(y)],'b--')
+legend(loc=3,prop={'size':11})
+savefig("Emergent_Relationship.pdf")
+
+# Calculate PDF from Emergent Constraint
+ytit=ylab
+xtit=xlab
+ybest,ybest_pr=EC_pdf_basic(x,y,x_obs,dx_obs,xlab,ylab)
+print(' Equal-weight prior Estimate of Y = ',ybest_pr)
+print('Emergent Constraint Estimate of Y = ',ybest)

Emergent_constraint/lin_reg_DP.py

@@ -0,0 +1,64 @@
+def lin_reg_DP(x,y) :
+    import numpy as np
+    from pylab import plot, legend, xlabel, ylabel
+    from pylab import  figure, mean, exp, sqrt, sum
+    from scipy import stats
+    from lin_reg_DP import lin_reg_DP
+
+# Based on "Least Squares fitting" equations from Math World website.
+# This version checked against data on the Wikipedia "Simple Linear Regression" pages.
+# It also calculates the +/- 1 sigma confidence limits in the regression [xfit,yband]
+#
+# IN THIS VERSION THE YBAND PREDICTION ERRORS ARE CALCULATED 
+# ACCORDING TO DAVID PEARSON (AS USED IN COX ET AL., 2013)
+
+    nx=len(x)
+    ny=len(y)
+
+    xm=mean(x)
+    ym=mean(y)
+
+    x2=x*x
+    y2=y*y
+    xy=x*y
+
+    ssxx=sum(x2)-nx*xm*xm
+    ssyy=sum(y2)-ny*ym*ym
+    ssxy=sum(xy)-ny*xm*ym
+
+    b=ssxy/ssxx
+    a=ym-b*xm
+
+    yf=a+b*x
+
+    r2=ssxy**2/(ssxx*ssyy)
+
+    e2=(y-yf)**2
+    s2=sum(e2)/(nx-2)
+
+    s=sqrt(s2)
+
+    da=s*sqrt(1.0/nx+xm*xm/ssxx)
+    db=s/sqrt(ssxx)
+
+
+# Calculate confidence limits on fit (see Wikipedia page on "Simple Linear Regression")
+#    minx=min(x)-1.1*(max(x)-min(x))
+#    maxx=max(x)+1.1*(max(x)-min(x))
+    minx=min(x)
+    maxx=max(x)
+    nfit=100
+    dx=(maxx-minx)/nfit
+    #xfit=minx+dx*range(0,nfit)
+    xfit=minx+[dx*i for i in range(0,100)]
+    #xfit=minx+dx*np.linspace(0,1,nfit-1)
+    yfit=a+b*xfit
+    yband=np.zeros(nfit)
+
+# David Pearson's formula for "Prediction Error"
+    for n in range (0,nfit):
+        yband[n]=sqrt(s2*(1.0+1.0/nx+(xfit[n]-xm)**2/ssxx))
+    pass
+
+    return yf,a,b,da,db,xfit,yfit,yband;
+

Emergent_constraint/plot_scat_basic.py

@@ -0,0 +1,17 @@
+def plot_scat_basic(x,y,xlab,ylab) :
+    import numpy as np
+    from pylab import plot, show, bar, legend, colors, axes, xlabel, ylabel
+    from pylab import title, savefig, axis, figure, semilogx, mean, exp, sqrt
+    from pylab import log, arctan
+    import scipy
+    import matplotlib.pyplot as plt
+
+
+    fig = plt.figure()
+    ax = fig.add_subplot(111)
+    plot(x,y,'ro')
+    xlabel(xlab,size=14)
+    ylabel(ylab,size=14)
+
+    return;
+