def slice_xp(current_data):
     from pylab import find,ravel
     x = current_data.x
     y = current_data.y
     dy = current_data.dy
     q = current_data.q
     ij = find((y <= latitude+dy+dy/2.) & (y > latitude+dy-dy/2.))
     x_slice = ravel(x)[ij]
     eta_slice = ravel(q[3,:,:])[ij]
     return x_slice, eta_slice

 def slice_y(current_data):
     from pylab import find,ravel
     x = current_data.x
     y = current_data.y
     dx = current_data.dx
     q = current_data.q
     ij = find((x <= longitude+dx/2.) & (x > longitude-dx/2.))
     y_slice = ravel(y)[ij]
     eta_slice = ravel(q[3,:,:])[ij]
     return y_slice, eta_slice

    def plot_topo_xsec(current_data):
        from pylab import find,ravel
        x = current_data.x
        y = current_data.y
        dy = current_data.dy

        ij = find((y <= dy/2.) & (y > -dy/2.))
        x_slice = ravel(x)[ij]
        b_slice = ravel(b(current_data))[ij]
        return x_slice, b_slice

 def slice_x(current_data):
     from pylab import find,ravel
     x = current_data.x
     y = current_data.y
     dy = current_data.dy
     q = current_data.q
     ij = find((y <= latitude+dy/2.) & (y > latitude-dy/2.))
     x_slice = ravel(x)[ij]
     eta_slice = ravel(q[3,:,:])[ij]
     print "+++ min eta = ",eta_slice.min()
     return x_slice, eta_slice

    def xsec(current_data):
        # Return x value and surface eta at this point, along y=0
        from pylab import find,ravel
        x = current_data.x
        y = current_data.y
        dy = current_data.dy

        ij = find((y <= dy/2.) & (y > -dy/2.))
        x_slice = ravel(x)[ij]
        eta_slice = ravel(plot_eta(current_data))[ij]
        return x_slice, eta_slice

    def xsec(current_data):
        # Return x value and surface eta at this point, along y=0
        from pylab import find,ravel
        x = current_data.x
        y = current_data.y
        dy = current_data.dy
        q = current_data.q

        ij = find((y <= y0+dy/2.+1e-8) & (y > y0-dy/2.))
        x_slice = ravel(x)[ij]
        eta_slice = ravel(q[3,:,:])[ij]
        return x_slice, eta_slice

 def xsec_hu(current_data):
     # Return x value and discharge at this point, along y=0
     from pylab import find,ravel,where,sqrt
     x = current_data.x
     y = current_data.y
     dy = current_data.dy
     q = current_data.q
     hu = q[1,:]
     ij = find((y <= 0.76+dy/2.) & (y > 0.76-dy/2.))
     x_slice = ravel(x)[ij]
     hu_slice = ravel(hu)[ij]
     return x_slice, hu_slice

 def B_x(current_data):
     from pylab import find,ravel
     x = current_data.x
     y = current_data.y
     dy = current_data.dy
     q = current_data.q
     ij = find((y <= latitude+dy+dy/2.) & (y > latitude+dy-dy/2.))
     x_slice = ravel(x)[ij]
     eta_slice = ravel(q[3,:,:])[ij]
     h_slice = ravel(q[0,:,:])[ij]
     B_slice = eta_slice - h_slice
     return x_slice, B_slice

 def B_y(current_data):
     from pylab import find,ravel
     x = current_data.x
     y = current_data.y
     dy = current_data.dy
     q = current_data.q
     dx = current_data.dx
     q = current_data.q
     ij = find((x <= longitude+dx/2.) & (x > longitude-dx/2.))
     y_slice = ravel(y)[ij]
     eta_slice = ravel(q[3,:,:])[ij]
     h_slice = ravel(q[0,:,:])[ij]
     B_slice = eta_slice - h_slice
     return y_slice, B_slice

 def ysec(current_data):
     # Return x value and surface eta at this point, along y=0
     from pylab import find,ravel
     x = current_data.x
     y = current_data.y
     dx = current_data.dx
     q = current_data.q
     t = current_data.t
     xlower = current_data.xlower
     xupper = current_data.xupper
     
     ij = find((x <= dx/2.) & (x > -dx/2.))
     y_slice = ravel(y)[ij]
     eta_slice = ravel(q[3,:,:])[ij]
     #print 'data'
     #print y_slice
     #fname = 'surface%d'%int(t)+'.txt'
     fname = r'/Users/huitang/clawpack-5.2.0/geoclaw/examples/tsunami/bowl-radial/_output/'+'surface_t=%d'%int(t)+'_ylower=%d'%int(y[1,1])+'.txt'
     fname2 = r'/Users/huitang/clawpack-5.2.0/geoclaw/examples/tsunami/bowl-radial/_output/location'+'_ylower=%d'%int(y[1,1])+'.txt'
     fname3 = r'/Users/huitang/clawpack-5.2.0/geoclaw/examples/tsunami/bowl-radial/_output/'+'momentums_t=%d'%int(t)+'_ylower=%d'%int(y[1,1])+'.txt'
     fname4 = r'/Users/huitang/clawpack-5.2.0/geoclaw/examples/tsunami/bowl-radial/_output/'+'depth_t=%d'%int(t)+'_ylower=%d'%int(y[1,1])+'.txt'
     savetxt(fname3,ravel(q[2,:,:])[ij])
     savetxt(fname,ravel(q[3,:,:])[ij])
     savetxt(fname4,ravel(q[0,:,:])[ij])
     savetxt(fname2,ravel(y)[ij])
     print '***********************************************'
     #print y[1]
     print ravel(q[0,:,:])[ij]
     #print y_slice
     #print y[1,1]
     #print xlower
     #print eta_slice
     return y_slice, eta_slice

    def xsec_B(current_data):
        # Return x value and B at this point, along y=0
        from pylab import find,ravel,where,sqrt
        x = current_data.x
        y = current_data.y
        dy = current_data.dy
        q = current_data.q
        h = q[0,:]
        eta = q[3,:]
        B = eta - h

        ij = find((y <= dy/2.) & (y > -dy/2.))
        x_slice = ravel(x)[ij]
        B_slice = ravel(B)[ij]
        return x_slice, B_slice

def Boxplot ( x, ax=None, offset=0.005, **kwargs ):
    """Creates a Tufte boxplot

    :Parameters:
        *x*
            values to be summarized
        *ax*
            target axes
        *offset*
            offset to mark central part and median gap

    :Optional keyword arguments:
        See pylab.boxplot
    """
    if ax is None:
        ax = pl.gca()
    # call a boxplot and manipulate it
    # how to get the offset right?
    if hasattr (x,'shape'):
        if len(x.shape)==1:
            if hasattr(x[0], 'shape'):
                x = list(x)
            else:
                x = [x,]
        elif len(x.shape) == 2:
            nr, nc = x.shape
            if nr==1:
                x = [x]
            elif nc==1:
                x = [x.ravel()]
            else:
                x = [x[:,i] for i in xrange(nc)]
        else:
            raise ValueError, "input x can have no more than 2 dimensions"
    if not hasattr(x[0],'__len__'):
        x = [x]

    positions = kwargs.setdefault ( 'positions', range(1,len(x)+1) )
    pdist = max(positions)-min(positions)
    offs = 0.5*float(pdist)/len(positions)

    if kwargs.setdefault('vert', 0 ) == 1:
        ax.set_ylim ( positions[0]-offs, positions[-1]+offs )
        ax.set_yticks ( positions )
    else:
        ax.set_xlim ( positions[0]-offs, positions[-1]+offs )
        ax.set_xticks ( positions )

    c = kwargs.setdefault ( 'color', [0,0,0] )
    c = dvis.color.colorsequence ( c )

    for i,pos in enumerate ( positions ):
        d = pl.ravel(x[i])
        kwargs['color'] = c[i%len(c)]
        ax.add_artist (
                BoxplotArtist ( pos,
                    calculate_boxplot_stats ( d, **kwargs ),
                    offset,
                    **kwargs ) )

    def xsec_s(current_data):
        # Return x value and speed at this point, along y=0
        from pylab import find,ravel,where,sqrt
        x = current_data.x
        y = current_data.y
        dy = current_data.dy
        q = current_data.q
        h = q[0,:]
        dry_tol = 0.001
        u = where(h>dry_tol, q[1,:]/h, 0.)
        v = where(h>dry_tol, q[2,:]/h, 0.)
        s = sqrt(u**2 + v**2)
        s = s / sqrt(9.81/0.97)  # so comparable to eta

        ij = find((y <= dy/2.) & (y > -dy/2.))
        x_slice = ravel(x)[ij]
        s_slice = ravel(s)[ij]
        return x_slice, s_slice

def pcolor_matrix_pylab(A, fname='pcolor_matrix_matplotlib'):
    """
    pcolor_matrix_pylab() implements a matlab-like 'pcolor' function to
    display the large elements of a matrix in pseudocolor using the Python Imaging
    Library.
    """
    try:
        import matplotlib
        matplotlib.use('Agg')
        import pylab
    except ImportError:
        if pedobj.kw['messages'] == 'verbose':
            print '[ERROR]: pyp_graphics/pcolor_matrix_pylab() was unable to import the matplotlib module!'
        logging.error('pyp_graphics/pcolor_matrix_pylab() was unable to import the matplotlib module!')
        return 0

    try:
        import numpy
        pylab.clf()
        x = pylab.arange(A.shape[0])
        X, Y = pylab.meshgrid(x,x)

        xmin = min(pylab.ravel(X))
        xmax = max(pylab.ravel(X))
        pylab.xlim(xmin, xmax)
        ymin = min(pylab.ravel(Y))
        ymax = max(pylab.ravel(Y))
        pylab.ylim(ymin, ymax)
        pylab.axis('off')

        pylab.pcolor(X, Y, pylab.transpose(A))#, shading='flat')
        pylab.clim(0.0, 1.0)
        plotfile = '%s.png' % (fname)
        myplotfile = open(plotfile,'w')
        pylab.savefig(myplotfile)
        myplotfile.close()
        return 1
    except:
        if pedobj.kw['messages'] == 'verbose':
            print '[ERROR]: pyp_graphics/pcolor_matrix_pylab() was unable to create the plot %s.' % (plotfile)
        logging.error('pyp_graphics/pcolor_matrix_pylab() was unable to create the plot %s.', (plotfile))
        return 0

def apply_mask(x):
    """
    Gets arrays with NaN from MAT files and applies python masked_where
    """
    f = pl.find(pl.isnan(x) == 1)
    l1, l2 = x.shape 
    x = pl.ravel(x)
    x[f] = 0
    x.shape = (l1,l2)
    x = pl.ma.masked_where(x == 0, x)
    return x

def plot_contour_lines(plot_data, xy_data=None, show_regression=False, \
        show_equation=False, smooth_steps=0, num_contour_lines=10, \
        label_contours=False, line_cmap=cm.hot, fill_cmap=cm.gray,dark=True,
        graph_shape='rect', prob_axes=False, **kwargs):
    """helper plots one or more series of contour line data
    calls the initializing functions, doesn't output figure
    
    takes: plot_data, xy_data, show_regression, show_equation, smooth,
        num_contour_lines, label_contours, line_cmap, fill_cmap, graph_shape,
        and **kwargs passed on to init_graph_display.
                 
           plot_data = (x_bin, y_bin, data_matrix dot_colors)
    """
    if prob_axes:
        extent = (0,1,0,1)
    else:
        a = gca()
        extent = a.get_xlim()+a.get_ylim()
    #init graph display, rectangular due to needed colorbar space
    init_graph_display(graph_shape=graph_shape,
        dark=dark, **kwargs)
    #plots the contour data
    for x_bin,y_bin,data in plot_data:
        orig_max = max(ravel(data))
        scaled_data = (data/orig_max*255).astype('b')
        if smooth_steps and (Image is not None):
            orig_shape = data.shape
            im = Image.fromstring('L', data.shape, scaled_data)
            for i in range(smooth_steps):
                im = im.filter(ImageFilter.BLUR)
            new_data = fromstring(im.tostring(), 'b')
            data = reshape(new_data.astype('i')/255.0 * orig_max, orig_shape)
        
        if fill_cmap is not None:
            im = imshow(data, interpolation='bicubic', extent=extent, \
                origin='lower', cmap=fill_cmap)
        result=contour(x_bin,y_bin,data, num_contour_lines,
                              origin='lower',linewidth=2,
                              extent=extent, cmap=line_cmap)
        if label_contours:
            clabel(result, fmt='%1.1g')

    #add the colorbar legend to the side
    cb = colorbar()
    cb.ax.axisbg = 'black'

    if show_regression:
        equation=plot_regression_line(xy_data[0],xy_data[1],prob_axes=prob_axes)
        if show_equation:
            add_regression_equations([equation])

def describe_vars(d):
    m = mc.Model(d)

    df = pandas.DataFrame(columns=["type", "value", "logp"], index=[n.__name__ for n in m.nodes], dtype=object)
    for n in m.nodes:
        k = n.__name__
        df.ix[k, "type"] = type(n).__name__

        if hasattr(n, "value"):
            rav = pl.ravel(n.value)
            if len(rav) == 1:
                df.ix[k, "value"] = n.value
            elif len(rav) > 1:
                df.ix[k, "value"] = "%.1f, ..." % rav[0]

        df.ix[k, "logp"] = getattr(n, "logp", pl.nan)

    return df.sort("logp")

def describe_vars(d):
    m = mc.Model(d)

    df = pd.DataFrame(columns=['type', 'value', 'logp'],
                          index=[n.__name__ for n in m.nodes],
                          dtype=object)
    for n in m.nodes:
        k = n.__name__
        df.ix[k, 'type'] = type(n).__name__

        if hasattr(n, 'value'):
            rav = pl.ravel(n.value)
            if len(rav) == 1:
                df.ix[k, 'value'] = n.value
            elif len(rav) > 1:
                df.ix[k, 'value'] = '%.1f, ...' % rav[0]

        df.ix[k, 'logp'] = getattr(n, 'logp', pl.nan)

    return df.sort('logp')

def norm_hist_bins(y, bins=10, normed='height'):
    """Just like the matplotlib mlab.hist, but can normalize by height.

    normed can be 'area' (produces matplotlib behavior, area is 1), 
    any False value (no normalization), or any True value (normalization).

    Original docs from matplotlib:

    Return the histogram of y with bins equally sized bins.  If bins
    is an array, use the bins.  Return value is
    (n,x) where n is the count for each bin in x

    If normed is False, return the counts in the first element of the
    return tuple.  If normed is True, return the probability density
    n/(len(y)*dbin)
    
    If y has rank>1, it will be raveled
    Credits: the Numeric 22 documentation
    """
    y = asarray(y)
    if len(y.shape)>1: y = ravel(y)

    if not iterable(bins):
        ymin, ymax = min(y), max(y)
        if ymin==ymax:
            ymin -= 0.5
            ymax += 0.5

        if bins==1: bins=ymax
        dy = (ymax-ymin)/bins
        bins = ymin + dy*arange(bins)
    n = searchsorted(sort(y), bins)
    n = diff(concatenate([n, [len(y)]]))
    if normed:
        if normed == 'area':
            db = bins[1]-bins[0]
        else:
            db = 1.0
        return 1/(len(y)*db)*n, bins
    else:
        return n, bins

    def xsec(current_data):
        # Return x value and surface eta at this point, along y=0
        from pylab import find,ravel
        x = current_data.x
        y = current_data.y
        dy = current_data.dy
        q = current_data.q
        aux = current_data.aux

        ij = find((y <= dy/2.) & (y > -dy/2.))
        x_slice = ravel(x)[ij]
        gamma_slice = ravel(aux[0,:,:])[ij]
        pinf_slice = ravel(aux[1,:,:])[ij]
        rho_slice = ravel(q[0,:,:])[ij]
        momx_slice = ravel(q[1,:,:])[ij]
        momy_slice = ravel(q[2,:,:])[ij]
        ene_slice = ravel(q[3,:,:])[ij]
        P_slice = (gamma_slice - 1.0)*(ene_slice - 0.5*(momx_slice**2 + momy_slice**2)/rho_slice)
        P_slice = P_slice - gamma_slice*pinf_slice
        # Convert to Psi and centimeters
        P_slice = P_slice*0.000145038 - 14.6959488
        x_slice = 100*x_slice
        return x_slice, P_slice