def test_pce_f2py(self):
        """ test solve_pce_f2py """ 

        N = 128; Yp = 0.24; Nrho = NT = N; z=3.0

        T = np.linspace( 1.0e4, 1.0e5, NT ) * ra.U.K
        fcA_H2 = 1.0; fcA_He2 = 1.0; fcA_He3 = 1.0

        hmr = ra.uv_bgnd.HM12_Photorates_Table()
        H1i = np.ones(N) * hmr.H1i(z)
        He1i = np.ones(N) * hmr.He1i(z)
        He2i = np.ones(N) * hmr.He2i(z)

        kchem = ra.atomic.ChemistryRates( T, fcA_H2, fcA_He2, fcA_He3,
                                          H1i=H1i, He1i=He1i, He2i=He2i )

        nH = np.linspace( 1.0e-4, 1.0e-3, Nrho ) / ra.U.cm**3
        nHe = nH * 0.25 * Yp / (1-Yp)

        x_pce_1D = ra.f2py.Solve_PCE( nH, nHe, kchem )

        y = ( x_pce_1D.He2 + 2 * x_pce_1D.He3 ) * nHe / nH
        xH1 = H.analytic_soltn_xH1( nH, y, kchem )

        err = np.abs( (x_pce_1D.H1 - xH1) / xH1 )
        ok = not np.any( err > TOL )

        self.assertTrue( ok )        

def hex_cube(x0, x1, y0, y1, z0, z1, Mx, My, Mz):
    """Creates a uniform hexahedral Q1 mesh covering [x0,x1]*[y0,y1]*[z0,z1]."""

    x = np.linspace(x0, x1, Mx)
    y = np.linspace(y0, y1, My)
    z = np.linspace(z0, z1, Mz)

    def ii(i, j, k):
        return (Mx*My)*k + (j*Mx + i)

    pts = np.zeros((Mx*My*Mz, 3), dtype='f8')
    for k in xrange(Mz):
        for j in xrange(My):
            for i in xrange(Mx):
                pts[ii(i,j,k), :] = (x[i], y[j], z[k])

    hexes = np.zeros(((Mx-1)*(My-1)*(Mz-1), 8), dtype='i')
    n = 0
    for k in xrange(Mz-1):
        for j in xrange(My-1):
            for i in xrange(Mx-1):
                hexes[n] = [ii(i,j,k), ii(i+1,j,k), ii(i+1,j+1,k), ii(i,j+1,k),
                            ii(i,j,k+1), ii(i+1,j,k+1), ii(i+1,j+1,k+1), ii(i,j+1,k+1)]
                n += 1

    return x, y, z, pts, hexes

    def _setup_rank3(self):
        a = np.linspace(0, 39, 40).reshape((2, 4, 5), order='F').astype(self.dt)
        b = np.linspace(0, 23, 24).reshape((2, 3, 4), order='F').astype(self.dt)

        y_r = array([[[    0.,   184.,   504.,   912.,  1360.,   888.,   472.,   160.,],
            [   46.,   432.,  1062.,  1840.,  2672.,  1698.,   864.,   266.,],
            [  134.,   736.,  1662.,  2768.,  3920.,  2418.,  1168.,   314.,],
            [  260.,   952.,  1932.,  3056.,  4208.,  2580.,  1240.,   332.,] ,
            [  202.,   664.,  1290.,  1984.,  2688.,  1590.,   712.,   150.,] ,
            [  114.,   344.,   642.,   960.,  1280.,   726.,   296.,    38.,]],

            [[   23.,   400.,  1035.,  1832.,  2696.,  1737.,   904.,   293.,],
             [  134.,   920.,  2166.,  3680.,  5280.,  3306.,  1640.,   474.,],
             [  325.,  1544.,  3369.,  5512.,  7720.,  4683.,  2192.,   535.,],
             [  571.,  1964.,  3891.,  6064.,  8272.,  4989.,  2324.,   565.,],
             [  434.,  1360.,  2586.,  3920.,  5264.,  3054.,  1312.,   230.,],
             [  241.,   700.,  1281.,  1888.,  2496.,  1383.,   532.,    39.,]],

            [[   22.,   214.,   528.,   916.,  1332.,   846.,   430.,   132.,],
             [   86.,   484.,  1098.,  1832.,  2600.,  1602.,   772.,   206.,],
             [  188.,   802.,  1698.,  2732.,  3788.,  2256.,  1018.,   218.,],
             [  308.,  1006.,  1950.,  2996.,  4052.,  2400.,  1078.,   230.,],
             [  230.,   692.,  1290.,  1928.,  2568.,  1458.,   596.,    78.,],
             [  126.,   354.,   636.,   924.,  1212.,   654.,   234.,     0.,]]],
            dtype=self.dt)

        return a, b, y_r

def mpl_palette(name, n_colors=6):
    """Return discrete colors from a matplotlib palette.

    Note that this handles the qualitative colorbrewer palettes
    properly, although if you ask for more colors than a particular
    qualitative palette can provide you will fewer than you are
    expecting.

    Parameters
    ----------
    name : string
        name of the palette
    n_colors : int
        number of colors in the palette

    Returns
    -------
    palette : list of tuples
        palette colors in r, g, b format

    """
    brewer_qual_pals = {"Accent": 8, "Dark2": 8, "Paired": 12,
                        "Pastel1": 9, "Pastel2": 8,
                        "Set1": 9, "Set2": 8, "Set3": 12}

    cmap = getattr(mpl.cm, name)
    if name in brewer_qual_pals:
        bins = np.linspace(0, 1, brewer_qual_pals[name])[:n_colors]
    else:
        bins = np.linspace(0, 1, n_colors + 2)[1:-1]
    palette = list(map(tuple, cmap(bins)[:, :3]))

    return palette

def plotISVar():
    plt.figure()
    plt.title('Variance minimization problem (call).\nVertical lines mark the minima.')
    for K in [0.6, 0.8, 1.0, 1.2]:
        theta = np.linspace(-0.6, 2)
        var = [BS.exactCallVar(K*s0, theta) for theta in theta]
        minth = theta[np.argmin(var)]
        line, = plt.plot(theta, var, label=str(K))
        plt.axvline(minth, color=line.get_color())

    plt.xlabel(r'$\theta$')
    plt.ylabel('call variance')
    plt.legend(title=r'$K/s_0$', loc='upper left')
    plt.autoscale(tight=True)

    plt.figure()
    plt.title('Variance minimization problem (put).\nVertical lines mark the minima.')
    for K in [0.8, 1.0, 1.2, 1.4]:
        theta = np.linspace(-2, 0.5)
        var = [BS.exactPutVar(K*s0, theta) for theta in theta]
        minth = theta[np.argmin(var)]
        line, = plt.plot(theta, var, label=str(K))
        plt.axvline(minth, color=line.get_color())

    plt.xlabel(r'$\theta$')
    plt.ylabel('put variance')
    plt.legend(title=r'$K/s_0$', loc='upper left')
    plt.autoscale(tight=True)

def test_BadBCDerivativesNoParam():
    problem = scikits.bvp_solver.ProblemDefinition(num_ODE = 1,
                                                   num_parameters =0,
                                        num_left_boundary_conditions = 1,
                                        boundary_points = (-numpy.pi/2.0, numpy.pi/2.0),
                                        function = functionNoParamGood,
                                        boundary_conditions = boundary_conditionsNoParamGood,
                                        function_derivative = dfunctionNoParamGood,
                                        boundary_conditions_derivative = dbconditionsNoParamBad1)

    nose.tools.assert_raises(ValueError, scikits.bvp_solver.solve, problem,
                                solution_guess = 0,
                                initial_mesh = numpy.linspace(problem.boundary_points[0],problem.boundary_points[1], 21))

    problem2 = scikits.bvp_solver.ProblemDefinition(num_ODE = 1,
                                                   num_parameters =0,
                                        num_left_boundary_conditions = 0,
                                        boundary_points = (-numpy.pi/2.0, numpy.pi/2.0),
                                        function = functionNoParamGood,
                                        boundary_conditions = boundary_conditionsNoParamGood2,
                                        function_derivative = dfunctionNoParamGood,
                                        boundary_conditions_derivative = dbconditionsNoParamBad2)

    nose.tools.assert_raises(ValueError, scikits.bvp_solver.solve, problem2,
                                solution_guess = 0,
                                initial_mesh = numpy.linspace(problem2.boundary_points[0],problem2.boundary_points[1], 21))

def test_more_known_parametrization_together():
    R = 1
    P = 1
    toll = 7.e-3
    intervals = 5
    vs_order = 2
    n = (intervals*(vs_order)+1-1)

    #n = 18
    ii = np.linspace(0,1,n+1)
    n_1 = 2
    n_2 = 4
    control_points_3d = np.asarray(np.zeros([n+1,n_1,n_2,3]))#[np.array([R*np.cos(5*i * np.pi / (n + 1)), R*np.sin(5*i * np.pi / (n + 1)), P * i]) for i in range(0, n+1)]
    for k in range(n_1):
        for j in range(n_2):
            control_points_3d[:,k,j,0] = np.array([R*np.cos(5*i * np.pi / (n + 1))for i in ii])
            control_points_3d[:,k,j,1] = np.array([R*np.sin(5*i * np.pi / (n + 1))for i in ii])
            control_points_3d[:,k,j,2] = np.array([(k+j+1)*P*i for i in range(n+1)])
    #vsl = IteratedVectorSpace(UniformLagrangeVectorSpace(vs_order+1), np.linspace(0,1,intervals+1))
    vsl = AffineVectorSpace(UniformLagrangeVectorSpace(n+1),0,1)
    arky = ArcLengthParametrizer(vsl, control_points_3d)
    new_control_points_3d = arky.reparametrize()

    #print control_points_3d.shape, new_control_points_3d.shape
    tt = np.linspace(0,1,128)
    for k in range(n_1):
        for j in range(n_2):
            vals = vsl.element(control_points_3d)(tt)
            new_vals = vsl.element(new_control_points_3d)(tt)
            print (np.amax(np.abs(vals-new_vals))/(k+j+1)/P, (k+j+1))
            assert np.amax(np.abs(vals-new_vals))/(k+j+1)/P < toll

    def test_cmac(self):
        input_train = np.reshape(np.linspace(0, 2 * np.pi, 100), (100, 1))
        input_train_before = input_train.copy()
        input_test = np.reshape(np.linspace(np.pi, 2 * np.pi, 50), (50, 1))
        input_test_before = input_test.copy()

        target_train = np.sin(input_train)
        target_train_before = target_train.copy()
        target_test = np.sin(input_test)

        cmac = algorithms.CMAC(
            quantization=100,
            associative_unit_size=32,
            step=0.2,
            verbose=False,
        )
        cmac.train(input_train, target_train, epochs=100)

        predicted_test = cmac.predict(input_test)
        predicted_test = predicted_test.reshape((len(predicted_test), 1))
        error = metrics.mean_absolute_error(target_test, predicted_test)

        self.assertAlmostEqual(error, 0.0024, places=4)

        # Test that algorithm didn't modify data samples
        np.testing.assert_array_equal(input_train, input_train_before)
        np.testing.assert_array_equal(input_train, input_train_before)
        np.testing.assert_array_equal(target_train, target_train_before)

def plot_mle_graph(function,
                   mle_params,
                   x_start=eps, x_end=1 - eps,
                   y_start=eps, y_end=1 - eps, resolution=100,
                   x_label="x", y_label="y",
                   show_constraint=False,
                   show_optimum=False):
    x = np.linspace(x_start, x_end, resolution)
    y = np.linspace(y_start, y_end, resolution)
    xx, yy = np.meshgrid(x, y)
    np_func = np.vectorize(lambda x, y: function(x, y))
    z = np_func(xx, yy)

    optimal_loss = function(*mle_params)
    levels_before = np.arange(optimal_loss - 3.0, optimal_loss, 0.25)
    levels_after = np.arange(optimal_loss, min(optimal_loss + 2.0, -0.1), 0.25)

    fig = plt.figure()
    contour = plt.contour(x, y, z, levels=np.concatenate([levels_before, levels_after]))
    plt.xlabel(x_label)
    plt.ylabel(y_label)

    if show_constraint:
        plt.plot(x, 1 - x)
    if show_optimum:
        plt.plot(mle_params[0], mle_params[1], 'ro')
    plt.clabel(contour)
    return mpld3.display(fig)

def allowed_region( V_nj, ave_j ):

    # read PCs
    PC1 = V_nj[0]
    PC2 = V_nj[1]
    n_band = len( PC1 )
    band_ticks = np.arange( n_band )

    x_ticks = np.linspace(-0.4,0.2,RESOLUTION)
    y_ticks = np.linspace(-0.2,0.4,RESOLUTION)
    x_mesh, y_mesh, band_mesh = np.meshgrid( x_ticks, y_ticks, band_ticks, indexing='ij' )
    vec_mesh = x_mesh * PC1[ band_mesh ] + y_mesh * PC2[ band_mesh ] + ave_j[ band_mesh ]

    x_grid, y_grid = np.meshgrid( x_ticks, y_ticks, indexing='ij' )
    prohibited_grid = np.zeros_like( x_grid )

    for ii in xrange( len( x_ticks ) ) :
        for jj in xrange( len( y_ticks ) ) :

            if np.any( vec_mesh[ii][jj] < 0. ) :
                prohibited_grid[ii][jj] = 1
                if np.any( vec_mesh[ii][jj] > 1. ) :
                    prohibited_grid[ii][jj] = 3
            elif np.any( vec_mesh[ii][jj] > 1. ) :
                prohibited_grid[ii][jj] = 2
            else :
                prohibited_grid[ii][jj] = 0

    return x_grid, y_grid, prohibited_grid

def make_let_im(let_file, dim = 16, y_lo = 70, y_hi = 220, x_lo = 10, x_hi = 200, edge_pix = 150, plot_let = False):

    letter = mpimg.imread(let_file)

    letter = letter[y_lo:y_hi, x_lo:x_hi, 0]
    for i in range(letter.shape[1]):
        if letter[0:edge_pix, i].any() == 0:   # here is to remove the edge
            letter[0:edge_pix, i] = 1
    
    plt.imshow(letter, cmap='gray')
    plt.grid('off')
    plt.show()
        
    x = np.arange(letter.shape[1])
    y = np.arange(letter.shape[0])

    f2d = interp2d(x, y, letter)

    x_new = np.linspace(0, letter.shape[1], dim)    # dim = 16
    y_new = np.linspace(0, letter.shape[0], dim)

    letter_new = f2d(x_new, y_new)
    letter_new -= np.mean(letter_new)
    
    if plot_let: 
        plt.imshow(letter_new, cmap = 'gray')
        plt.grid('off')
        plt.show()
        
    letter_flat = letter_new.flatten()   # letter_flat is a 1-dimensional array containing 256 elements
    
    return letter_new, letter_flat

def normalize_input(params):
    if pc_id == 0:
        print 'normalize_input'
        dt = params['dt_rate'] # [ms] time step for the non-homogenous Poisson process 
        L_input = np.zeros((params['n_exc'], params['t_stimulus']/dt))

        v_max = params['v_max']
        if params['log_scale']==1:
            v_rho = np.linspace(v_max/params['N_V'], v_max, num=params['N_V'], endpoint=True)
        else:
            v_rho = np.logspace(np.log(v_max/params['N_V'])/np.log(params['log_scale']),
                            np.log(v_max)/np.log(params['log_scale']), num=params['N_V'],
                            endpoint=True, base=params['log_scale'])
        v_theta = np.linspace(0, 2*np.pi, params['N_theta'], endpoint=False)
        index = 0
        for i_RF in xrange(params['N_RF_X']*params['N_RF_Y']):
            index_start = index
            for i_v_rho, rho in enumerate(v_rho):
                for i_theta, theta in enumerate(v_theta):
                    fn = params['input_rate_fn_base'] + str(index) + '.dat'
                    L_input[index, :] = np.loadtxt(fn)
                    print 'debug', fn
                    index += 1
            index_stop = index
            print 'before', i_RF, L_input[index_start:index_stop, :].sum()
            if (L_input[index_start:index_stop, :].sum() > 1):
                L_input[index_start:index_stop, :] /= L_input[index_start:index_stop, :].sum()
            print 'after', i_RF, L_input[index_start:index_stop, :].sum()

        for i in xrange(params['n_exc']):
            output_fn = params['input_rate_fn_base'] + str(i) + '.dat'
            print 'output_fn:', output_fn
            np.savetxt(output_fn, L_input[i, :])
    if comm != None:
        comm.barrier()

 def all_GL(self, q, maxpiv=None):
     """return (piv, f_binodal_gas, f_binodal_liquid, f_spinodal_gas, f_spinodal_liquid) at insersion works piv sampled between the critical point and maxpiv (default to 2.2*critical pressure)"""
     fc, pivc = self.critical_point(q)
     Fc = np.log(fc)
     #start sensibly above the critical point
     startp = pivc*1.1
     fm = fminbound(self.mu, fc, self.maxf(), args=(startp, q))
     fM = fminbound(lambda f: -self.pv(f, startp, q), 0, fc)
     initial_guess = np.log([0.5*fM, 0.5*(fm+self.maxf())])
     #construct the top of the GL binodal
     if maxpiv is None:
         maxpiv = startp*2
     topp = 1./np.linspace(1./startp, 1./maxpiv)
     topGL = [initial_guess]
     for piv in topp:
         topGL.append(self.binodalGL(piv, q, topGL[-1]))
     #construct the GL binodal between the starting piv and the critical point
     botp = np.linspace(startp, pivc)[:-1]
     botGL = [initial_guess]
     for piv in botp:
         botGL.append(self.binodalGL(piv, q, botGL[-1]))
     #join the two results and convert back from log
     binodal = np.vstack((
         [[pivc, fc, fc]],
         np.column_stack((botp, np.exp(botGL[1:])))[::-1],
         np.column_stack((topp, np.exp(topGL[1:])))[1:]
         ))
     #spinodal at the same pivs
     spinodal = self.spinodalGL(q, binodal[:,0])
     #join everything
     return np.column_stack((binodal, spinodal[:,1:]))

def simulate():
    # Plotting the PDF of Unif(0, 1)
    pyplot.subplot(211)
    x = np.linspace(stats.uniform.ppf(0), stats.uniform.ppf(1), 100)
    pyplot.title('PDF of Unif(0, 1)')
    pyplot.plot(x, stats.uniform.pdf(x))

    print("Xn is for n=1000: ", get_xn(1000))
    E_Xn = 0.5  # As we know, E(Xn) is equal to mu which is 0.5
    print("E(Xn) is : ", E_Xn)

    n = np.linspace(1, 1000, 1000)
    X_ns = []
    for i in range(1, 1001):
        X_ns.append(get_xn(i))
    pyplot.subplot(212)
    pyplot.title('f(n,Xn)')
    pyplot.plot(n, X_ns, '-g')

    print("Xn for n=1", get_xn(1))
    print("Xn for n=5", get_xn(5))
    print("Xn for n=25", get_xn(25))
    print("Xn for n=100", get_xn(100))

    pyplot.show()

def mesh(npts=(101, 101), closed=False):
    """Generate a meshgrid on the unit sphere.

    Uniformly sample the polar angle, theta, and the azimuthal angle, phi.

    Parameters
    ----------
    npts : int or tuple
        Number of angle points sampled.
    closed : bool
        Whether to generate an open mesh (like `np.ogrid`)
        or a closed mesh like `np.mgrid` or `np.meshgrid`.
        By default, an open grid is generated.

    Returns
    -------
    theta, phi : (N,) or (N,M) ndarray
        Sampling of the polar angle.  Shape depends on ``open``.

    """
    if np.isscalar(npts):
        npts = (npts, npts)

    theta = np.linspace(0, np.pi, npts[0])[:, None]
    phi = np.linspace(0, 2 * np.pi, npts[1])

    if open:
        return theta, phi
    else:
        mg_phi, mg_theta = np.meshgrid(phi_grid, theta_grid)
        return mg_theta, mg_phi

def disc_norm():
    x = np.linspace(-3,3,100)
    y = st.norm.pdf(x,0,1)
    fig, ax = plt.subplots()
    fig.canvas.draw()
    
    ax.plot(x,y)
    
    fill1_x = np.linspace(-2,-1.5,100)
    fill1_y = st.norm.pdf(fill1_x,0,1)
    fill2_x = np.linspace(-1.5,-1,100)
    fill2_y = st.norm.pdf(fill2_x,0,1)
    ax.fill_between(fill1_x,0,fill1_y,facecolor = 'blue', edgecolor = 'k',alpha = 0.75)
    ax.fill_between(fill2_x,0,fill2_y,facecolor = 'blue', edgecolor = 'k',alpha = 0.75)
    for label in ax.get_yticklabels():
        label.set_visible(False)
    for tick in ax.get_xticklines():
        tick.set_visible(False)
    for tick in ax.get_yticklines():
        tick.set_visible(False)
    
    plt.rc("font", size = 16)
    plt.xticks([-2,-1.5,-1])
    labels = [item.get_text() for item in ax.get_xticklabels()]
    labels[0] = r"$v_k$"
    labels[1] = r"$\varepsilon_k$"
    labels[2] = r"$v_{k+1}$"
    ax.set_xticklabels(labels)
    plt.ylim([0, .45])

    
    plt.savefig('discnorm.pdf')
    plt.clf()

def color_cycle_example(ax):
    L = 6
    x = np.linspace(0, L)
    ncolors = len(plt.rcParams['axes.prop_cycle'])
    shift = np.linspace(0, L, ncolors, endpoint=False)
    for s in shift:
        ax.plot(x, np.sin(x + s), 'o-')

def get_nxnyr_cd():
    box = cfg.pms['shape']
    boxMpc = np.array([cfg.pms['xyMpc'],cfg.pms['xyMpc'],cfg.pms['zMpc']])
    lx = boxMpc[0]/2.; ly = boxMpc[1]/2.; lz = boxMpc[2]
    z,d = get_z_d(cfg.pms['zi'],cfg.pms['zf'])
    
    # front of box -- don't use bc grid will extend
    #                 outside the standard box
    # nx_max = lx / np.sqrt(lx*lx+d[0]*d[0]) # nx_min = - nx_max
    # ny_max = ly / np.sqrt(ly*ly+d[0]*d[0]) # ny_min = - ny_max
    # r_max = np.sqrt(lx*lx+ly*ly+(d[0]+lz)*(d[0]+lz)) # r_min = d[0]

    # back of box -- throws away half the box but whatever
    df = d[0]+lz
    nx_max = lx / np.sqrt(lx*lx+df*df) # nx_min = - nx_max
    ny_max = ly / np.sqrt(ly*ly+df*df) # ny_min = - ny_max
    r_max = np.sqrt(lx*lx+ly*ly+df*df) # r_min = d[0]

    print nx_max, ny_max

    nxcd = np.linspace(-nx_max,nx_max,box[0])
    nycd = np.linspace(-ny_max,ny_max,box[1])
    print 2*nx_max/box[0], 2*ny_max/box[1]
    rcd = np.linspace(d[0],r_max,box[2])
    return nxcd,nycd,rcd

def test_known_parametrization():
    R = 1
    P = 1
    toll = 2.e-3

    n = 10
    ii = np.linspace(0,1,n+1)
    control_points_3d = np.asarray(np.zeros([n+1,3]))#[np.array([R*np.cos(5*i * np.pi / (n + 1)), R*np.sin(5*i * np.pi / (n + 1)), P * i]) for i in range(0, n+1)]
    print (control_points_3d.shape)
    control_points_3d[:,0] = np.array([R*np.cos(5*i * np.pi / (n + 1))for i in ii])
    control_points_3d[:,1] = np.array([R*np.sin(5*i * np.pi / (n + 1))for i in ii])
    control_points_3d[:,2] = np.array([P*i for i in range(n+1)])
    vsl = AffineVectorSpace(UniformLagrangeVectorSpace(n+1),0,1)
    arky = ArcLengthParametrizer(vsl, control_points_3d)
    new_control_points_3d = arky.reparametrize()

    #new_arky = ArcLengthParametrizer(vsl, new_control_points_3d)
    #new_new_control_points_3d = arky.reparametrize()
    tt = np.linspace(0, 1, 128)

    vals = vsl.element(control_points_3d)(tt)
    #print vals
    new_vals = vsl.element(new_control_points_3d)(tt)
    #print vals.shape, new_vals.shape
    print (np.amax((np.abs(vals-new_vals))))
    assert np.amax(np.abs(control_points_3d-new_control_points_3d))/P < toll

    def test_parse_spot_request(self):

        tests = [('spot:20S,154W|4,3',
                  {'location': {'latitude': -20., 'longitude':-154.},
                   'model': 'gfs',
                   'type': 'spot',
                   'send-image': False,
                   'hours': np.linspace(0, 96, 33).astype('int'),
                   'vars': ['wind'],
                   'warnings': ['No variable requested, defaulting to WIND']}),
                 ('spot: 20S,154W |4 ,3',
                  {'location': {'latitude': -20., 'longitude':-154.},
                   'model': 'gfs',
                   'send-image': False,
                   'type': 'spot',
                   'hours': np.linspace(0, 96, 33).astype('int'),
                   'vars': ['wind'],
                   'warnings': ['No variable requested, defaulting to WIND']}),
                 ('spot:20S,154W|4,3|wind|image',
                  {'location': {'latitude': -20., 'longitude':-154.},
                   'model': 'gfs',
                   'type': 'spot',
                   'send-image': True,
                   'hours': np.linspace(0, 96, 33).astype('int'),
                   'vars': ['wind'],
                   'warnings': []}),
                 ]

        for request, expected in tests:
            self.maxDiff = 3000
            actual = saildocs.parse_spot_request(request)
            np.testing.assert_array_equal(actual.pop('hours'),
                                          expected.pop('hours'))
            self.assertDictEqual(actual, expected)