def plot_stoch_value():    
    #Compute Solution==========================================================
    sigma = .5
    mu = 4*sigma
    K = 7
    Gamma, eps = discretenorm.discretenorm(K,mu,sigma)
    
    N = 100
    W = sp.linspace(0,1,N)
    V = sp.zeros((N,K))
    
    u = lambda c: sp.sqrt(c)
    beta = 0.99
    
    X,Y= sp.meshgrid(W,W)
    Wdiff = Y-X
    index = Wdiff < 0
    Wdiff[index] = 0
    
    util_grid = u(Wdiff)
    
    util3 = sp.tile(util_grid[:,:,sp.newaxis],(1,1,K))
    eps_grid = eps[sp.newaxis,sp.newaxis,:]
    eps_util = eps_grid*util3
    
    Gamma_grid = Gamma[sp.newaxis,:]
    
    delta = 1
    Vprime = V
    z = 0
    while (delta > 10**-9):
        z= z+1
        V = Vprime
        gamV = Gamma_grid*V
        Expval = sp.sum(gamV,1)
        Exp_grid = sp.tile(Expval[sp.newaxis,:,sp.newaxis],(N,1,K))
        arg = eps_util+beta*Exp_grid
        arg[index] = -10^10
        Vprime = sp.amax(arg,1)
        psi_ind = sp.argmax(arg,1)
        psi = W[psi_ind]
        delta = sp.linalg.norm(Vprime - V)
    
    #============================================================    
    #Plot 3D    
    x=sp.arange(0,N)
    y=sp.arange(0,K)
    X,Y=sp.meshgrid(x,y)
    fig1 = plt.figure()
    ax1= Axes3D(fig1)
    ax1.set_xlabel(r'$W$')
    ax1.set_ylabel(r'$\varepsilon$')
    ax1.set_zlabel(r'$V$')
    ax1.plot_surface(W[X],Y,sp.transpose(Vprime), cmap=cm.coolwarm)
    plt.savefig('stoch_value.pdf')

def stepfunction(s):
    ''' eddy diffusivity with discontinuity for testing of mixing scheme'''
    Nparticles = range(s.elements.z.shape[0])
    K = s.environment_profiles['z'] * 0 + 0.1
    K[s.environment_profiles['z'] < -20] = K[s.environment_profiles['z'] < -20] * 0 + 0.02
    N, Kprofiles = sp.meshgrid(Nparticles, K)
    return Kprofiles

def windspeed_Sundby1983(s):
    depths = s.environment_profiles['z']
    windspeed_squared = s.environment.x_wind**2 + s.environment.y_wind**2
    K = 76.1e-4 + 2.26e-4 * windspeed_squared
    #valid = windspeed_squared < 13**2
    Kprofiles, N = sp.meshgrid(K, depths)
    return Kprofiles

 def setUp(self):
     self.n = 100
     self.m = 10
     self.kappa1=0.0
     self.errLevelPos = 6
     self.errLevelCurv= 5
     self.x = numpy.linspace(-1.5,1.5,self.n)
     self.lamb = numpy.linspace(-1.0,1.0,self.m)
     X,LAMB = scipy.meshgrid(self.x,self.lamb)
     lambField = DC.FieldContainer(LAMB,
                                   unit = '1 V / m**3',
                                   longname='parameter',
                                   shortname='\lambda')
     xField = DC.FieldContainer(X[0],
                                unit = '1 m',
                                longname = 'position',
                                shortname = 'x')
     V=[]
     for i in xrange(len(lambField.data)):
         u = xField.data
         V.append(-lambField.data[i]/2* u**2 + u**4/4-u*self.kappa1)
     self.I = DC.FieldContainer(numpy.array(V),
                                longname = 'intensity',
                                shortname='I')
     self.I.dimensions[-1]=xField
     self.I0 = DC.FieldContainer(numpy.ones(self.I.data.shape,'float'),
                                longname = 'white reference',
                                shortname='I_0')
     self.I0.dimensions[-1]=xField
     self.Id = DC.FieldContainer(numpy.zeros(self.I.data.shape,'float'),
                                longname = 'darf reference',
                                shortname='I_d')
     self.Id.dimensions[-1]=xField
     self.sampleC = DC.SampleContainer([self.I,self.I0,self.Id])

def plot_mixture(fig, mixture, data):
  plt.subplot(fig)

  # Plot data.
  plt.scatter(data[:, 0], data[:, 1], color='r', marker='+', alpha=0.5)

  delta = 1.0
  coords = []
  # TODO: Should be configurable...
  xs = scipy.arange(-30, 30, delta)
  ys = scipy.arange(-30, 30, delta)
  for x in xs:
    for y in ys:
      coords.append((x, y))

  coords = scipy.asarray(coords)

  Z = None
  for i in range(mixture.degree):
    gaussian = MultivariateGaussianDensity(mixture.means[i], mixture.covs[i])
    z = mixture.mixcoeffs[i] * gaussian.multpdf(coords)

    if Z is None:
      Z = z
    else:
      Z += z

  X, Y = scipy.meshgrid(xs, ys)
  Z.shape = X.shape

  CS = plt.contour(X, Y, Z, 30)
  plt.plot(mixture.means[:, 1], mixture.means[:, 0], 'o', color='green')

def steepest_descent(X,Y, step=.001, tol=1e-5, maxiter=5000, bounds=[-5,5], res=.1):
    w = betahat(X,Y)
    a = sp.exp(X.dot(w))
    yhat = (a.T/sp.sum(a, axis=1)).T
    grad = X.T.dot(yhat-Y)
    while la.norm(grad)>tol and maxiter > 0:
        w = w - grad*step
        a=sp.exp(X.dot(w))
        yhat = (a.T/sp.sum(a,axis=1)).T
        grad = X.T.dot(yhat-Y)
        maxiter -= 1

    rang = sp.arange(bounds[0],bounds[1]+res, res)
    Xg, Yg = sp.meshgrid(rang, rang)
    
    Xm = sp.c_[sp.ones_like(Xg.flatten()), Xg.flatten(), Yg.flatten()]
    Xmdot = Xm.dot(w)
    
    types = Xmdot.argmax(axis=1)
    print Xm.shape
    #plot the regions
    c = ['b','r','g']
    for i in range(Xmdot.shape[1]):
        plot_on = types==i
        plt.plot(Xm[plot_on,1], Xm[plot_on,2], c[i]+'.', X[Y[:,i]==1,1], X[Y[:,i]==1,2], c[i]+'o')

    #plot the data segmented
    #tmp = sp.where(Ydat[:,0]==True)[0]
    #plt.plot(Xdat[tmp,1], Xdat[tmp,2], 'o', Xdat[~tmp,1], Xdat[~tmp,2], 'o')
    plt.show()

def test_mayaxes():

    from mayaxes import mayaxes
    from scipy import sqrt,sin,meshgrid,linspace,pi
    import mayavi.mlab as mlab
        
    resolution = 200
    lambda_var = 3
    theta = linspace(-lambda_var*2*pi,lambda_var*2*pi,resolution)
    
    x, y = meshgrid(theta, theta)
    r = sqrt(x**2 + y**2)
    z = sin(r)/r
    
    fig = mlab.figure(size=(1024,768))
    surf = mlab.surf(theta,theta,z,colormap='jet',opacity=1.0,warp_scale='auto') 
    mayaxes(title_string='Figure 1: Diminishing polar cosine series', \
        xlabel='X data',ylabel='Y data',zlabel='Z data',handle=surf)
    
    fig.scene.camera.position = [435.4093863309094, 434.1268937227623, 315.90311468125287]
    fig.scene.camera.focal_point = [94.434632665253829, 93.152140057106593, -25.071638984402856]
    fig.scene.camera.view_angle = 30.0
    fig.scene.camera.view_up = [0.0, 0.0, 1.0]
    fig.scene.camera.clipping_range = [287.45231734040635, 973.59247058049255]
    fig.scene.camera.compute_view_plane_normal()
    fig.scene.render()   
    
    mlab.show() 

def GetSq(filename,Nsamp=1):
    if type(filename)==str:
        filename=[filename]
    attrs=GetAttr(filename[0])
    #hfile=h5py.File(filename,'r')
    dpath,args=GetStat(filename,Nsamp)
    filetype=''
    try:
        filetype=attrs['type']
    except KeyError:
        mo=re.match('.*/?[0-9]+-(StatSpinStruct)\.h5',filename[0])
        filetype=mo.groups()[0]
    if filetype!='StatSpinStruct':
        raise InputFileError('\"{0}\" is not a static structure factor file'.format(filename))
    N=pow(attrs['L'],2)
    Sq=sc.zeros((Nsamp,5,N),complex)
    for sample,b in enumerate(args):
        for d in b:
            hfile=h5py.File(dpath[d][0],'r')
            dim=hfile[dpath[d][1]].shape[0]
            Sq[sample,:dim,:]+=hfile[dpath[d][1]][0:,0::2]+1j*hfile[dpath[d][1]][0:,1::2]
            hfile.close()
        Sq[sample,:,:]/=len(b)
    qx,qy=sc.meshgrid(np.arange(attrs['L']),np.arange(attrs['L']))
    return qx.flatten(),qy.flatten(),Sq

def plot_problem(X, y, h=None, surfaces=True) :
    '''
    Plots a two-dimensional labeled dataset (X,y) and, if function h(x) is given, 
    the decision boundaries (surfaces=False) or decision surfaces (surfaces=True)
    '''
    assert X.shape[1] == 2, "Dataset is not two-dimensional"
    if h!=None : 
        # Create a mesh to plot in
        r = 0.02  # mesh resolution
        x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1
        y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1
        xx, yy = sp.meshgrid(sp.arange(x_min, x_max, r),
                             sp.arange(y_min, y_max, r))
        XX = sp.c_[xx.ravel(), yy.ravel()]
        try:
            #Z_test = h(XX)
            #if sp.shape(Z_test) == () :
            #    # h returns a scalar when applied to a matrix; map explicitly
            #    Z = sp.array(map(h,XX))
            #else :
            #    Z = Z_test
            Z = sp.array(map(h,XX))
        except ValueError:
            # can't apply to a matrix; map explicitly
            Z = sp.array(map(h,XX))
        # Put the result into a color plot
        Z = Z.reshape(xx.shape)
        if surfaces :
            plt.contourf(xx, yy, Z, cmap=plt.cm.Pastel1)
        else :
            plt.contour(xx, yy, Z)
    # Plot the dataset
    plt.scatter(X[:,0],X[:,1],c=y, cmap=plt.cm.Paired,marker='o',s=50);

    def split(self, sagi, meri):
        """ utilizes geometry.grid to change the rectangle into a generalized surface,
        it is specified with a single set of basis vectors to describe the meridonial,
        normal, and sagittal planes."""
        ins = float((sagi - 1))/sagi
        inm = float((meri - 1))/meri
        stemp = self.norm.s/sagi
        mtemp = self.meri.s/meri

        z,theta = scipy.meshgrid(scipy.linspace(-self.norm.s*ins,
                                                self.norm.s*ins,
                                                sagi),
                                 scipy.linspace(-self.meri.s*inm,
                                                self.meri.s*inm,
                                                meri))

        vecin =geometry.Vecr((self.sagi.s*scipy.ones(theta.shape),
                              theta+scipy.pi/2,
                              scipy.zeros(theta.shape))) #this produces an artificial
        # meri vector, which is in the 'y_hat' direction in the space of the cylinder
        # This is a definite patch over the larger problem, where norm is not normal
        # to the cylinder surface, but is instead the axis of rotation.  This was
        # done to match the Vecr input, which works better with norm in the z direction
               
        pt1 = geometry.Point(geometry.Vecr((scipy.zeros(theta.shape),
                                            theta,
                                            z)),
                             self)

        pt1.redefine(self._origin)

        vecin = vecin.split()

        x_hat = self + pt1 #creates a vector which includes all the centers of the subsurface

        out = []
        #this for loop makes me cringe super hard
        for i in xrange(meri):
            try:
                temp = []
                for j in xrange(sagi):
                    inp = self.rot(vecin[i][j])
                    temp += [Cyl(geometry.Vecx(x_hat.x()[:,i,j]),
                                 self._origin,
                                 [2*stemp,2*mtemp],
                                 self.sagi.s,
                                 vec=[inp, self.norm.copy()],
                                 flag=self.flag)]
                out += [temp]
            except IndexError:
                inp = self.rot(vecin[i])
                out += [Cyl(geometry.Vecx(x_hat.x()[:,i]),
                            self._origin,
                            [2*stemp,2*mtemp],
                            self.norm.s,
                            vec=[inp, self.norm.copy()],
                            flag=self.flag)]
                

        return out

def Problem1Real():
    beta = 0.9;
    T = 10;
    N = 100;
    u = lambda c: sp.sqrt(c);
    W = sp.linspace(0,1,N);
    X, Y = sp.meshgrid(W,W);
    Wdiff = Y-X
    index = Wdiff <0;
    Wdiff[index] = 0;
    util_grid = u(Wdiff);
    util_grid[index] = -10**10;
    V = sp.zeros((N,T+2));
    psi = sp.zeros((N,T+1));


    for k in xrange(T,-1,-1):
        val = util_grid + beta*sp.tile(sp.transpose(V[:,k+1]),(N,1));
        vt = sp.amax(val, axis = 1);
        psi_ind = sp.argmax(val,axis = 1)
        V[:,k]    = vt;
        psi[:,k]    = W[psi_ind];

    
    return V,psi

def m_circles(mags, phase_min=-359.75, phase_max=-0.25):
    """Constant-magnitude contours of the function Gcl = Gol/(1+Gol), where
    Gol is an open-loop transfer function, and Gcl is a corresponding
    closed-loop transfer function.

    Usage
    =====
    contours = m_circles(mags, phase_min, phase_max)

    Parameters
    ----------
    mags : array-like
        Array of magnitudes in dB of the M-circles
    phase_min : degrees
        Minimum phase in degrees of the N-circles
    phase_max : degrees
        Maximum phase in degrees of the N-circles

    Return values
    -------------
    contours : complex array
        Array of complex numbers corresponding to the contours.
    """
    # Convert magnitudes and phase range into a grid suitable for
    # building contours
    phases = sp.radians(sp.linspace(phase_min, phase_max, 2000))
    Gcl_mags, Gcl_phases = sp.meshgrid(10.0**(mags/20.0), phases)
    return closed_loop_contours(Gcl_mags, Gcl_phases)

def plotDecisionBoundary(data, X, theta):
    """
    Plots the data points X and y into a new figure with the
    decision boundary defined by theta
    """
    plotData(data)
    if X.shape[1] <= 3:
        plot_x = sp.array([min(X[:, 1]), max(X[:, 1])])
        plot_y = (-1.0 / theta[2]) * (theta[1] * plot_x + theta[0])
        plt.plot(plot_x, plot_y)
        plt.axis([30, 100, 30, 100])
    else:
        u = sp.linspace(-1, 1.5, 50)
        v = sp.linspace(-1, 1.5, 50)
        z = sp.zeros((len(u), len(v)))

        for i in range(0, len(u)):
            for j in range(0, len(v)):
                z[i, j] = (mapFeature(sp.array([u[i]]), sp.array([v[j]]))).dot(theta)

        z = z.T # important to transpose z before calling contour

        u, v = sp.meshgrid(u, v)

        plt.contour(u, v, z, [0.0, 0.0])

def getMeshGrid(t, z, minT=None, maxT=None, export=False):
    if minT is None:
        minT = t.min()
    if maxT is None:
        maxT = t.max()

    nt = settings.XGDS_PLOT_PROFILE_TIME_GRID_RESOLUTION
    intervalStart = minT + TIME_OFFSET_DAYS
    if export:
        dt = EXPORT_TIME_RESOLUTION
        intervalStart = int(float(intervalStart) / dt) * dt
    else:
        dt = float(maxT - minT) / nt
    tvals = np.arange(intervalStart,
                      maxT + TIME_OFFSET_DAYS,
                      dt)

    if 0:
        minZ = z.min()
        maxZ = z.max()
        nz = settings.XGDS_PLOT_PROFILE_Z_GRID_RESOLUTION
        dz = float(maxZ - minZ) / nz
        zvals = np.arange(minZ, maxZ, dz)

    zvals = np.arange(*settings.XGDS_PLOT_PROFILE_Z_RANGE)

    return scipy.meshgrid(tvals, zvals)

def plotOrthogonalField(sh, b):
    center=(np.array(sh)-1)/2.0
    C,R=sp.meshgrid(np.array(range(sh[1]), dtype=np.float64), np.array(range(sh[0]), dtype=np.float64))
    R=R-center[0]+b[0]
    C=C-center[1]+b[1]
    plt.figure()
    plt.quiver(R, -C)

 def make_grid_flip_y(self):
     """A way to make a basic grid from x,y data, inverting y.
     """
     self.y_grid, self.x_grid = scipy.meshgrid(self.y[::-1], self.x[:], indexing='ij')
     self.dims = (self.ny, self.nx)
     self.dy = self.y[0]-self.y[1]
     self.dx = self.x[1]-self.x[0]

def n_circles(phases, mag_min=-40.0, mag_max=12.0):
    """Constant-phase contours of the function Gcl = Gol/(1+Gol), where
    Gol is an open-loop transfer function, and Gcl is a corresponding
    closed-loop transfer function.

    Usage
    ===== 
    contours = n_circles(phases, mag_min, mag_max)

    Parameters
    ----------
    phases : array-like
        Array of phases in degrees of the N-circles
    mag_min : dB
        Minimum magnitude in dB of the N-circles
    mag_max : dB
        Maximum magnitude in dB of the N-circles

    Return values
    -------------
    contours : complex array
        Array of complex numbers corresponding to the contours.
    """
    # Convert phases and magnitude range into a grid suitable for
    # building contours    
    mags = sp.linspace(10**(mag_min/20.0), 10**(mag_max/20.0), 2000)
    Gcl_phases, Gcl_mags = sp.meshgrid(sp.radians(phases), mags)
    return closed_loop_contours(Gcl_mags, Gcl_phases)

 def initiatespatial(self):
     self.spatial=True
     self.pMap = Basemap(projection='stere',lat_0=68.0,lon_0=15.0,llcrnrlon=3.,llcrnrlat=60.3,urcrnrlon=47.0,urcrnrlat=71.,resolution='c')
     self.x, self.y = self.pMap(self.lon, self.lat)
     XI,ETA = sp.meshgrid(range(self.y.shape[1]),range(self.y.shape[0]))
     self.xip, self.etap = XI.ravel(), ETA.ravel()
     self.Tree = spatial.KDTree(zip(self.x.ravel(), self.y.ravel()))

def Problem3Real():
    beta = 0.9
    N = 1000
    u = lambda c: sp.sqrt(c)
    W = sp.linspace(0,1,N)
    X, Y = sp.meshgrid(W,W)
    Wdiff = sp.transpose(X-Y)
    index = Wdiff <0
    Wdiff[index] = 0
    util_grid = u(Wdiff)
    util_grid[index] = -10**10
    
    Vprime = sp.zeros((N,1))
    psi = sp.zeros((N,1))
    delta = 1.0
    tol = 10**-9
    it = 0
    max_iter = 500
    
    while (delta >= tol) and (it < max_iter):
        V = Vprime
        it += 1;
        #print(it)
        val = util_grid + beta*sp.transpose(V)
        Vprime = sp.amax(val, axis = 1)
        Vprime = Vprime.reshape((N,1))
        psi_ind = sp.argmax(val,axis = 1)
        psi    = W[psi_ind]
        delta = sp.dot(sp.transpose(Vprime - V),Vprime-V)
    
    return psi

def plot_delta():     
    beta = 0.99
    N = 1000
    u = lambda c: sp.sqrt(c)
    W = sp.linspace(0,1,N)
    X, Y = sp.meshgrid(W,W)
    Wdiff = sp.transpose(X-Y)
    index = Wdiff <0
    Wdiff[index] = 0
    util_grid = u(Wdiff)
    util_grid[index] = -10**10
    
    Vprime = sp.zeros((N,1))
    delta = sp.ones(1)
    tol = 10**-9
    it = 0
    max_iter = 500
    
    while (delta[-1] >= tol) and (it < max_iter):
        V = Vprime
        it += 1;
        print(it)
        val = util_grid + beta*sp.transpose(V)
        Vprime = sp.amax(val, axis = 1)
        Vprime = Vprime.reshape((N,1))
        delta = sp.append(delta,sp.dot(sp.transpose(Vprime - V),Vprime-V))
        
    plt.figure()
    plt.plot(delta[1:])
    plt.ylabel(r'$\delta_k$')
    plt.xlabel('iteration')
    plt.savefig('convergence.pdf')