Correct color definition

Add rectangle plot for MCMC chains

work/tools/python/camel.py

@@ -403,7 +403,7 @@ def posterior1d( chains, params, nbin=50, smooth=1,
 #        ax = plt.gca()
 #        color_cycle = ax._get_lines.prop_cycler #ax._get_lines.color_cycle
 #        colors = [next(color_cycle) for i in range(nchain)]
-        colors = rcParams['axes.color_cycle']
+        colors = plt.rcParams['axes.color_cycle']
 
     if len(linestyles) < nchain:
         linestyles = ['-']*nchain
@@ -569,7 +569,7 @@ def triangle( chains, params, nbin=50, parnames=None, names=[], colors=[], lines
 #        ax = plt.gca()
 #        color_cycle = ax._get_lines.prop_cycler  #ax._get_lines.color_cycle
 #        colors = [next(color_cycle) for i in range(nchain)]
-        colors = rcParams['axes.color_cycle']
+        colors = plt.rcParams['axes.color_cycle']
 
     if len(cmaps) < nchain:
         cmaps = ['Blues','Greens','Greys','Oranges','Purples']
@@ -673,6 +673,94 @@ def triangle( chains, params, nbin=50, parnames=None, names=[], colors=[], lines
 
 
 
+def rectangle( chains, par0, params, nbin=50, parnames=None, names=[], colors=[], linestyles=[], smooth=1., fontsize=12,
+               contour=True, fill=False, cmaps=[], Norm=True, weights=None, ylim=None):
+    import matplotlib.ticker as mtick
+    
+    fig=plt.figure()
+    plt.subplots_adjust(hspace=0.001,wspace=0.001)
+
+    #force tuple
+    if not( isinstance(chains, list) or isinstance(chains, tuple)):
+        chains = (chains,)
+
+    nchain = len(chains)
+    npar   = len(params)
+
+    if len(colors) < nchain:
+        colors = plt.rcParams['axes.color_cycle']
+
+    if len(cmaps) < nchain:
+        cmaps = ['Blues','Greens','Greys','Oranges','Purples']
+
+    if type(fill) == bool:
+        fill = [fill]*nchain
+    if type(contour) == bool:
+        contour = [contour]*nchain
+
+    if len(linestyles) < nchain:
+        linestyles = ['-']*nchain
+    
+    if weights == None:
+        weights = [np.ones(len(c['chi2'])) for c in chains]
+
+    if parnames == None:
+        parnames = [parname.get(p,p) for p in params]
+
+    if ylim == None:
+        tmp = []
+        for chain in chains:
+            tmp = tmp+list(chain[par0])
+        ylim = list(np.mean(tmp) + np.dot(np.std(tmp), [-4,4]))
+
+    rge = {}
+    for par in params:
+        tmp = []
+        for chain in chains:
+            tmp = tmp+list(chain[par])
+        rge[par] = list(np.mean(tmp) + np.dot(np.std(tmp), [-4,4]))
+
+
+    for par in params:
+        xpar = params.index(par)
+
+        ax=plt.subplot( 1, npar, xpar+1)
+
+        for c in range(nchain):
+            chain = chains[c]
+            weight = weights[c]
+            hmap,binX,binY = np.histogram2d( chain[par], chain[par0], bins=nbin, range=[rge[par],ylim], weights=weight)
+            X,Y = np.meshgrid((binX[1:]+binX[:-1])/2.,(binY[1:]+binY[:-1])/2.)
+            hmap = nd.gaussian_filter( hmap.T, 2)
+            if fill[c]:
+                plt.contourf( X,Y, hmap, levels=np.append(list(ctr_level( hmap, [0.68,0.95])[::-1]),hmap.max()),
+                              colors=None, extent=tuple(rge[par]+ylim), cmap=plt.cm.get_cmap(cmaps[c]))
+            if contour[c]:
+                plt.contour( X,Y, hmap, levels=ctr_level( hmap, [0.68,0.95]),
+                             colors=colors[c], linestyles=linestyles[c], extent=tuple(rge[par]+ylim))
+
+        ax.ticklabel_format(useOffset=False)
+        ax.locator_params(tight=True, nbins=5)
+        ax.set_xlabel(parname.get(parnames[xpar],parnames[xpar]))
+        ax.tick_params(axis='both', which='major', labelsize=fontsize)
+        ax.yaxis.set_visible(xpar==0)
+        if xpar == 0:
+            ax.set_ylabel(parname.get(par0,par0))
+        
+#        if 'PS' in parname.get(parnames[par],parnames[par]):
+#            ax.xaxis.set_major_formatter(mtick.FormatStrFormatter('%.0e'))
+        
+    #write legend
+##     if len(names) == nchain:
+##         ax=plt.subplot( npar, npar, int(0.25*npar)*npar+0.85*npar-1)
+##         for c in range(nchain):
+##             plt.text(0.5,1.-0.2*(c+1), names[c], color=colors[c])
+##         ax.axis('off')
+
+    return(fig)
+
+
+
 def hist1d(x, bins=50, weights=None, smooth=0., normmax=True, **kwargs):
 #    ax = kwargs.get("ax", plt.gca())
     extents = kwargs.get("extents",[x.min(),x.max()])
@@ -1635,13 +1723,13 @@ class prof:
             all plot arguments
         """
         ax = plt.gca()
-##         from matplotlib import __version__
-##         if __version__ >= '1.5':
-##             colordef = next(ax._get_lines.prop_cycler)
-##         else:
-##             colordef = next(ax._get_lines.color_cycle)
+        from matplotlib import __version__
+        if __version__ >= '1.5':
+            colordef = next(ax._get_lines.prop_cycler)
+        else:
+            colordef = next(ax._get_lines.color_cycle)
             
-        color     = kwargs.pop( "color",     None)
+        color     = kwargs.pop( "color",  colordef)
         linestyle = kwargs.pop( "linestyle",  '-')
         label     = kwargs.pop( "label",     None)
         marker    = kwargs.pop( "marker",    'o')