def calculate_risk(vector):
    init_dollars = 100          # change this if you want
    DOB = days_of_bettin = 82  # change this if you want
    num_runs = 4096            # change this if you want
    risks = np.array(np.linspace(0.01,0.99,99)).T # change this too, but be careful
    risks = risks.reshape((risks.size, 1))
    prob, odds = vector[0], vector[1]
    dollars = np.full_like(risks, init_dollars)
    talley = np.full_like(risks, 0)
    multiplier = odds_to_pct(odds)
    if prob*multiplier - (1-prob) < 0:
        return 0
    else:
        for _ in xrange(num_runs):
            for day in xrange(DOB):
                outcome = np.random.random()
                if outcome < prob:
                    dollars += dollars * multiplier * risks
                else:
                    dollars -= dollars * risks
            if talley.all() == np.full_like(risks, 0).all():
                talley = dollars
            else:
                talley = np.c_[talley, dollars]
            dollars = np.full_like(risks, init_dollars)
        expectation = np.median(talley, axis=1)
        index = np.argmax(expectation)
        return risks[index][0]

def Compare_Dihedral(Dihedral):
    dx = 0.0001
    N = len(Dihedral)
    print N
    Dih = np.arange(0, 2*3.1415, dx)
    U1 = np.full_like(Dih, 0.0)
    U2 = np.full_like(Dih, 0.0)
    U3 = np.full_like(Dih, 0.0)
    
    for i in range(N):
        U1 += Dihedral[i]*np.cos(Dih)**i
    
    for i in range(5):
        print i
        U2 += Dihedral[i]*np.cos(Dih)**i
    
    Popt, Pcov = curve_fit(Multi, Dih, U1)

    for i in range(5):
        U3 += Popt[i]*np.cos(Dih)**i
    
    print Popt
    plt.figure()
    plt.plot(Dih, U1, label = 'Full Potential')
    plt.plot(Dih, U2, label = 'Truncated Approximation')
    plt.plot(Dih, U3, label = 'Optimized Approximation')
    plt.ylim((U2.min(), U2.max()))
    plt.xlim((0.0, 2*3.1415))
    plt.title('P2P1P1P2', fontsize = 30)
    plt.xlabel('Dihedral Angle (radians)', fontsize = 20)
    plt.ylabel('Energy (Kcal/mol)', fontsize = 20)
    plt.legend()
    plt.show()
    return

def Gen_PDF_CDF_OPLS( V, Beta ):
	"""
	This function takes in a numpy array V that containes the energetic coefficients for the OPLS style dihedral potential of the form:
			U = (1/2)V1(1+cos(phi)) + (1/2)V2(1-cos(2phi)) + (1/2)V3(1+cos(3phi)) + ....
	It then uses Boltzmann statistics along with the inverse temperature Beta to generate a PDF and CDF of the dihedral angle
	
	The output is two numpy arrays that represent the PDF and CDF associated with this potential energy function
	"""
	dx = 0.0001
	x = np.arange(0, 6.28, dx) # Generate discretized values for x (phi)
	U = np.full_like(x, 0.0) # Initialize potential energy array
	PDF = np.full_like(x, 0.0) # Initialize PDF array
	CDF_NN = np.full_like(x, 0.0) # Initialize non-normalized CDF array
	CDF = np.full_like(x, 0.0) # Initialize normalized CDF array
	norm = 0
	L = len(x.tolist()) 
	U = 0.5*(V[0]*(1 + np.cos(x)) + V[1]*(1 - np.cos(2*x)) + V[2]*(1 + np.cos(3*x)) + V[3]*(1 - np.cos(4*x)))
	PDF = np.exp(-U*Beta)
	
	for i in range(L-1):
		CDF_NN[i+1] = CDF_NN[i] + PDF[i]*dx
	
	for i in range(L):
		PDF[i] = PDF[i]/CDF_NN[-1]
		norm += PDF[i]*dx
	
	for i in range(L-1):
		CDF[i+1] = CDF[i] + PDF[i]*dx
		
	return PDF, CDF

def Compare_Angle( Angle):
    dx = 0.0001
    M = [0, 2 , 3, 4 , 5 , 6, 8, 10, 12, 14, 16, 18, 20, 22, 24]
    N = len(Angle) - 1
    Th0 = Angle[0]*(3.1415/180.)
    dTh = 2.0
    Th = np.arange(Th0-dTh, Th0 + dTh, dx)
    U1 = np.full_like(Th, 0.0)
    U2 = np.full_like(Th, 0.0)
    
    for i in range(1, N+1):
        U1 += Angle[i]*(Th - Th0)**M[i-1]
    
    #U2 = Angle[2]*(Th - Th0)**2 + Angle[3]*(Th - Th0)**3 + Angle[4]*(Th-Th0)**4
    U2 =1000.*(Th - Th0)**2
    plt.figure()
    plt.plot(Th, U1, label = 'Full Potential')
    plt.plot(Th, U2, label = 'Harmonic Approximation')
    plt.ylim((U1.min(), U2.max() ))
    plt.xlim((Th0-dTh, Th0 + dTh))
    plt.title('P1P1P2')
    plt.xlabel('Angle (Radians)', fontsize = 20)
    plt.ylabel('Potential Energy (Kcal/mol)', fontsize = 20)
    plt.legend()
    plt.show()

    return  

def Compare_Bond( Bond , title):
    dx = 0.0001
    N = len(Bond) - 1
    DR = .25
    r = np.arange(Bond[0]-DR, Bond[0]+DR, dx)
    U1 = np.full_like(r, 0.0)
    U2 = np.full_like(r, 0.0)
    
    for i in range(1,N+1):
        U1 += Bond[i]*(r - Bond[0])**(i+1)
        print i+1
        
    U2 = Bond[1]*(r - Bond[0])**2 + Bond[2]*(r-Bond[0])**3 + Bond[3]*(r-Bond[0])**4
    
    plt.figure()
    plt.plot(r, U1, label = 'Full Potential')
    plt.plot(r, U2, label = 'Class2 Approximation')
    plt.xlim((Bond[0]-DR,Bond[0]+DR))
    plt.ylim((U1.min(), U1.max()))
    plt.title('%s' % title, fontsize = 30)
    plt.xlabel('Bond Length (Angstrom)', fontsize = 20)
    plt.ylabel('Potential Energy (Kcal/mol)', fontsize = 20)
    plt.legend()
    plt.show()
    return

def assignReads(anchor, highestPeak, clusterSize, blockCount):
    global tagCount
    global clusterStart
    readMeans = np.empty(tagCount)
    readHeights = np.empty(tagCount)
    readMeans = np.full_like(readMeans, -1, dtype=np.double)
    readHeights = np.full_like(readHeights, -1, dtype=np.double)
    
    meanCounter = 0

    counterNew = 0
    counterOld = -1

    while counterOld != counterNew:
        dev = stddev(readMeans, readHeights, tagCount)
        counterOld = counterNew
        for start in anchor:
            if start.block == -1:
                mean = ((start.start + start.end) / 2) - clusterStart
                variance = args.sizescale * (abs(start.end - start.start) / 2)

                if (((mean - variance - dev) <= highestPeak and (mean + variance + dev) >= highestPeak) or (mean >= (highestPeak - args.merge) and mean <= (highestPeak + args.merge))):
                    readMeans[meanCounter] = mean
                    readHeights[meanCounter] = start.height
                    meanCounter += 1
                    start.block = blockCount
                    counterNew += 1

    return counterNew

    def resample(self, bin_count=120):
        start_i = int(self.t0 * self.sample_rate)
        end_i = util.clip(start_i + int(self.dt * self.sample_rate),
                          start_i, sys.maxsize)
        bin_size = (end_i - start_i) // bin_count
        if bin_size < 1:
            bin_size = 1
        bin_count = len(np.arange(start_i, end_i, bin_size))

        data = np.empty((self.data.shape[1], 2*bin_count, 4), dtype=np.float32)

        for i, column in enumerate(self.data):
            v = mea.min_max_bin(self.data[column].values[start_i:end_i],
                                bin_size, bin_count+1)
            col, row = mea.coordinates_for_electrode(column)
            row = 12 - row - 1
            x = np.full_like(v, col, dtype=np.float32)
            y = np.full_like(v, row, dtype=np.float32)
            t = np.arange(0, bin_count, 0.5, dtype=np.float32)
            data[i] = np.column_stack((x, y, t, v))

        # Update shader
        self.program['a_position'] = data.reshape(
            2*self.data.shape[1]*bin_count, 4)
        self.program['u_width'] = bin_count

 def test_manual_bounds(self, cuda=False):
     device = torch.device("cuda") if cuda else torch.device("cpu")
     for dtype in (torch.float, torch.double):
         # get a test module
         train_x = torch.tensor([[1.0, 2.0, 3.0]], device=device, dtype=dtype)
         train_y = torch.tensor([4.0], device=device, dtype=dtype)
         likelihood = GaussianLikelihood()
         model = ExactGP(train_x, train_y, likelihood)
         model.covar_module = RBFKernel(ard_num_dims=3)
         model.mean_module = ConstantMean()
         model.to(device=device, dtype=dtype)
         mll = ExactMarginalLogLikelihood(likelihood, model)
         # test the basic case
         x, pdict, bounds = module_to_array(
             module=mll, bounds={"model.covar_module.raw_lengthscale": (0.1, None)}
         )
         self.assertTrue(np.array_equal(x, np.zeros(5)))
         expected_sizes = {
             "likelihood.noise_covar.raw_noise": torch.Size([1]),
             "model.covar_module.raw_lengthscale": torch.Size([1, 3]),
             "model.mean_module.constant": torch.Size([1]),
         }
         self.assertEqual(set(pdict.keys()), set(expected_sizes.keys()))
         for pname, val in pdict.items():
             self.assertEqual(val.dtype, dtype)
             self.assertEqual(val.shape, expected_sizes[pname])
             self.assertEqual(val.device.type, device.type)
         lower_exp = np.full_like(x, 0.1)
         for p in ("likelihood.noise_covar.raw_noise", "model.mean_module.constant"):
             lower_exp[_get_index(pdict, p)] = -np.inf
         self.assertTrue(np.equal(bounds[0], lower_exp).all())
         self.assertTrue(np.equal(bounds[1], np.full_like(x, np.inf)).all())

def test_hpxgeom_coord_to_idx(nside, nested, coordsys, region, axes):
    import healpy as hp

    geom = HpxGeom(nside, nested, coordsys, region=region, axes=axes)
    lon = np.array([112.5, 135.0, 105.0])
    lat = np.array([75.3, 75.3, 74.6])
    coords = make_test_coords(geom, lon, lat)
    zidx = tuple([ax.coord_to_idx(t) for t, ax in zip(coords[2:], geom.axes)])

    if geom.nside.size > 1:
        nside = geom.nside[zidx]
    else:
        nside = geom.nside

    phi, theta = coords.phi, coords.theta
    idx = geom.coord_to_idx(coords)
    assert_allclose(hp.ang2pix(nside, theta, phi), idx[0])
    for i, z in enumerate(zidx):
        assert_allclose(z, idx[i + 1])

    # Test w/ coords outside the geometry
    lon = np.array([0.0, 5.0, 10.0])
    lat = np.array([75.3, 75.3, 74.6])
    coords = make_test_coords(geom, lon, lat)
    zidx = [ax.coord_to_idx(t) for t, ax in zip(coords[2:], geom.axes)]

    idx = geom.coord_to_idx(coords)
    if geom.region is not None:
        assert_allclose(np.full_like(coords[0], -1, dtype=int), idx[0])

    idx = geom.coord_to_idx(coords, clip=True)
    assert np.all(np.not_equal(np.full_like(coords[0], -1, dtype=int), idx[0]))

def _extract_current_results(data, curr, data_time):
    grid = data['models']['simulationGrid']
    plate_spacing = _meters(grid['plate_spacing'])
    zmesh = np.linspace(0, plate_spacing, grid['num_z'] + 1) #holds the z-axis grid points in an array
    beam = data['models']['beam']
    if data.models.simulationGrid.simulation_mode == '3d':
        cathode_area = _meters(grid['channel_width']) * _meters(grid['channel_height'])
    else:
        cathode_area = _meters(grid['channel_width'])
    RD_ideal = sources.j_rd(beam['cathode_temperature'], beam['cathode_work_function']) * cathode_area
    JCL_ideal = sources.cl_limit(beam['cathode_work_function'], beam['anode_work_function'], beam['anode_voltage'], plate_spacing) * cathode_area

    if beam['currentMode'] == '2' or (beam['currentMode'] == '1' and beam['beam_current'] >= JCL_ideal):
        curr2 = np.full_like(zmesh, JCL_ideal)
        y2_title = 'Child-Langmuir cold limit'
    else:
        curr2 = np.full_like(zmesh, RD_ideal)
        y2_title = 'Richardson-Dushman'
    return {
        'title': 'Current for Time: {:.4e}s'.format(data_time),
        'x_range': [0, plate_spacing],
        'y_label': 'Current [A]',
        'x_label': 'Z [m]',
        'points': [
            curr.tolist(),
            curr2.tolist(),
        ],
        'x_points': zmesh.tolist(),
        'y_range': [min(np.min(curr), np.min(curr2)), max(np.max(curr), np.max(curr2))],
        'y1_title': 'Current',
        'y2_title': y2_title,
    }

def watershed(image):
    hsv_image = color.rgb2hsv(image)

    low_res_image = rescale(hsv_image[:, :, 0], SCALE)
    local_mean = mean(low_res_image, disk(50))
    local_minimum_flat = np.argmin(local_mean)
    local_minimum = np.multiply(np.unravel_index(local_minimum_flat, low_res_image.shape), round(1 / SCALE))

    certain_bone_pixels = np.full_like(hsv_image[:, :, 0], False, bool)
    certain_bone_pixels[
    local_minimum[0] - INITIAL_WINDOW_SIZE/2:local_minimum[0]+INITIAL_WINDOW_SIZE/2,
    local_minimum[1] - INITIAL_WINDOW_SIZE/2:local_minimum[1]+INITIAL_WINDOW_SIZE/2
    ] = True

    certain_non_bone_pixels = np.full_like(hsv_image[:, :, 0], False, bool)
    certain_non_bone_pixels[0:BORDER_SIZE, :] = True
    certain_non_bone_pixels[-BORDER_SIZE:-1, :] = True
    certain_non_bone_pixels[:, 0:BORDER_SIZE] = True
    certain_non_bone_pixels[:, -BORDER_SIZE:-1] = True

    smoothed_hsv = median(hsv_image[:, :, 0], disk(50))
    threshold = MU * np.median(smoothed_hsv[certain_bone_pixels])

    possible_bones = np.zeros_like(hsv_image[:, :, 0])
    possible_bones[smoothed_hsv < threshold] = 1

    markers = np.zeros_like(possible_bones)
    markers[certain_bone_pixels] = 1
    markers[certain_non_bone_pixels] = 2

    labels = morphology.watershed(-possible_bones, markers)

    return labels

def fill_array( var1, var2 ):
    """
    fix fill_array such that it returns two numpy arrays of equal size

    use numpy.full_like

    """
    var1_a = np.asarray( var1 )
    var2_a = np.asarray( var2 )

    if var1_a.shape==():
        var1_a = np.asarray( [var1] )
    if var2_a.shape==():
        var2_a = np.asarray( [var2] )

    # Begin try/except block to handle all cases for filling an array
    while True:
        try:
            assert var1_a.shape == var2_a.shape
            break
        except: pass
        try:
            var1_a = np.full_like( var2_a, var1_a )
            break
        except: pass
        try:
            var2_a = np.full_like( var1_a, var2_a )
            break
        except: pass

        # If none of the cases properly handle it, throw error
        assert False, 'var1 and var2 must both be equal shape or size=1'

    return var1_a, var2_a

def optimize_img(init_img, solver_type, solver_param, max_iter, display, root_dir, net,
                 all_target_blob_names, targets, target_data_list):
    ensuredir(root_dir)

    solver_param.update({
        'maxiter': max_iter,
        'disp': True,
    })

    # Set initial value and reshape net
    set_data(net, init_img)
    x0 = np.ravel(init_img).astype(np.float64)

    mins = np.full_like(x0, -128)
    maxs = np.full_like(x0, 128)

    bounds = zip(mins, maxs)
    display_func = DisplayFunctor(net, root_dir, display)

    opt_res = optimize.minimize(
        objective_func,
        x0,
        args=(net, all_target_blob_names, targets, target_data_list),
        bounds=bounds,
        method=solver_type,
        jac=True,
        callback=display_func,
        options=solver_param,
    )
    print opt_res

def march(x,u_e,nu):
    dx = numpy.diff(x)
    du_e = numpy.gradient(u_e,numpy.gradient(x))
    delta = numpy.full_like(x,0.)
    lam = numpy.full_like(x,lam0)

    # Initial conditions must be a stagnation point. If u_e[0]>0
    # assume stagnation is at x=0 and integrate from x=0..x[0].
    if u_e[0]<0.01:                     # stagnation point
        delta[0] = numpy.sqrt(lam0*nu/du_e[0])
    elif x[0]>0:                        # just downstream
        delta[0] = numpy.sqrt(lam0*nu*x[0]/u_e[0])
        delta[0] += 0.5*x[0]*g_pohl(delta[0],0,u_e,du_e,nu)
        lam[0] = delta[0]**2*du_e[0]/nu
    else:
        raise ValueError('x=0 must be stagnation point')

    # march!
    for i in range(len(x)-1):
        delta[i+1] = heun(g_pohl,delta[i],i,dx[i],
                          u_e,du_e,nu)  # ...additional arguments
        lam[i+1] = delta[i+1]**2*du_e[i+1]/nu

        if lam[i+1] < -12: i-=1; break  # separation condition

    return delta,lam,i+1                # return with separation index

def test_flux_unit_conversion():
    # By default the flux units should be set to Jy
    s = Spectrum1D(flux=np.array([26.0, 44.5]), spectral_axis=np.array([400, 500]) * u.nm)
    assert np.all(s.flux == np.array([26.0, 44.5]) * u.Jy)
    assert s.flux.unit == u.Jy

    # Simple Unit Conversion
    s = Spectrum1D(flux=np.array([26.0, 44.5]) * u.Jy, spectral_axis=np.array([400, 500])*u.nm)
    converted_value = s.to_flux(unit=u.uJy)[0]
    assert ((26.0 * u.Jy).to(u.uJy) == converted_value)

    # Make sure incompatible units raise UnitConversionError
    with pytest.raises(u.UnitConversionError):
        converted_value = s.to_flux(unit=u.m)

    # Pass custom equivalencies
    s = Spectrum1D(flux=np.array([26.0, 44.5]) * u.Jy, spectral_axis=np.array([400, 500]) * u.nm)
    eq = [[u.Jy, u.m,
          lambda x: np.full_like(np.array(x), 1000.0, dtype=np.double),
          lambda x: np.full_like(np.array(x), 0.001, dtype=np.double)]]
    converted_value = s.to_flux(unit=u.m, equivalencies=eq)[0]
    assert 1000.0 * u.m == converted_value

    # Check if suppressing the unit conversion works
    s = Spectrum1D(flux=np.array([26.0, 44.5]) * u.Jy, spectral_axis=np.array([400, 500]) * u.nm)
    s.to_flux("uJy", suppress_conversion=True)
    assert s.flux[0] == 26.0 * u.uJy

def plot_energy(run_summary, x_axis):
            
    plt.figure()
    
    if x_axis == "time":
        x_variable = run_summary.times / yr
        xlabel = "Time [yr]"
        xscale = "linear"
        plt.xscale(xscale)
        xfmt = plt.gca().get_xaxis().get_major_formatter() # needs to be set AFTER plt.xscale()
        if xscale == "log":
            mask = x_variable > 1
        elif xscale == "linear":
            mask = np.full_like(x_variable, True, dtype=bool) 
            xfmt.set_powerlimits((-2, 2)) # force scientific notation outside this range

    elif x_axis == "checkpoints":
        x_variable = np.arange(len(run_summary.times))
        xlabel = "Checkpoint"
        xscale = "linear"
        mask = np.full_like(x_variable, True, dtype=bool) 

        plt.xscale(xscale)
        xfmt = plt.gca().get_xaxis().get_major_formatter() # needs to be set AFTER plt.xscale()

    else:
        raise NotImplementedError("can't recognize x_axis value: " + x_axis)

    E_err = (run_summary.E_tot - run_summary.E_tot[0]) / run_summary.E_tot[0]
    plt.plot(x_variable[mask], E_err[mask])
    plt.xscale(xscale)
    plt.xlabel(xlabel)   
    plt.gca().xaxis.set_major_formatter(xfmt)
    plt.ylabel("Fractional Change (Energy)")
    SNe_distplot(run_summary, x_axis)

    plt.figure()
    plt.plot(x_variable[mask], run_summary.E_tot[mask], label="E_tot" )
    plt.plot(x_variable[mask], run_summary.E_kin[mask], label="E_kin" )
    plt.plot(x_variable[mask], run_summary.E_int[mask], label="E_int" )
    plt.legend(loc="best")
    plt.xscale(xscale)
    plt.xlabel(xlabel) 
    plt.gca().xaxis.set_major_formatter(xfmt)
    plt.ylabel("Energy [erg]")
    SNe_distplot(run_summary, x_axis)

    
    plt.figure()
    plt.plot(x_variable[mask], run_summary.E_R_tot[mask], label="E_Remnant" )
    plt.legend(loc="best")
    plt.xscale(xscale)
    plt.xlabel(xlabel)  
    plt.gca().xaxis.set_major_formatter(xfmt)
    plt.ylabel("Energy [erg]")
    SNe_distplot(run_summary, x_axis)

    if x_axis == "checkpoints":
        plt.xlim(xmin=0)

    def plot(self):
        self.graph.clearPlot()
        self.validindices = numpy.empty((0,), dtype=int)
        self.current_selection = []
        group, target_indices = self.selected_split()
        self.warning([0, 1])
        self.error(1)

        if self.data and group is not None and target_indices:
            X = self.data.X
            I1 = grouputils.group_selection_mask(
                self.data, group, target_indices)
            I2 = ~I1
            if isinstance(group, grouputils.RowGroup):
                X = X.T

            N1, N2 = numpy.count_nonzero(I1), numpy.count_nonzero(I2)

            if not N1 or not N2:
                self.error(
                    1, "Target labels most exclude/include at least one value."
                )

            if N1 < 2 and N2 < 2:
                self.warning(
                    0, "Insufficient data to compute statistics. "
                       "More than one measurement per class should be provided"
                )

            X1, X2 = X[:, I1], X[:, I2]
            if numpy.any(X1 < 0.0) or numpy.any(X2 < 0):
                self.error(
                    "Negative values in the input. The inputs cannot be in "
                    "ratio scale."
                )
                X1 = numpy.full_like(X1, numpy.nan)
                X2 = numpy.full_like(X2, numpy.nan)

            with numpy.errstate(divide="ignore", invalid="ignore"):
                fold = numpy.log2(numpy.mean(X1, axis=1) /
                                  numpy.mean(X2, axis=1))
                # TODO: handle missing values better (mstats)
                _, P = scipy.stats.ttest_ind(X1, X2, axis=1, equal_var=True)
                logP = numpy.log10(P)
                if numpy.isscalar(logP):
                    # ttest_ind does not preserve output shape if either
                    # a or b is empty
                    logP = numpy.full(fold.shape, numpy.nan)

            mask = numpy.isfinite(fold) & numpy.isfinite(logP)
            self.validindices = numpy.flatnonzero(mask)
            self.graph.setPlotData(numpy.array([fold[mask], -logP[mask]]).T)

            self.infoLabel.setText("%i genes on input" % len(fold))
            # ("{displayed} displayed, {undef} with undefined ratio "
            #  "or t-statistics.")

            if not len(numpy.flatnonzero(mask)):
                self.warning(1, "Could not compute statistics for any genes!")

def add_location_data(ds, lat, lon):
    lat = lat if lat else LAT_FILL
    lon = lon if lon else LON_FILL
    lat_array = np.full_like(ds.time.values, lat)
    lon_array = np.full_like(ds.time.values, lon)

    ds['lat'] = ('obs', lat_array, {'axis': 'Y', 'units': 'degrees_north', 'standard_name': 'latitude'})
    ds['lon'] = ('obs', lon_array, {'axis': 'X', 'units': 'degrees_east', 'standard_name': 'longitude'})

    def distance_to(self, source):
        src_lats = num.full_like(self.lats, fill_value=source.lat)
        src_lons = num.full_like(self.lons, fill_value=source.lon)

        target_coords = self.get_latlon()
        target_lats = target_coords[:, 0]
        target_lons = target_coords[:, 1]
        return distance_accurate50m_numpy(
            src_lats, src_lons, target_lats, target_lons)

    def test_basic(self):
        # Check derivative at endpoints
        n1_10 = np.arange(1, 10)
        dataset0 = np.column_stack([n1_10, np.full_like(n1_10, 0), np.full_like(n1_10, -1)])
        FuncData(_smirnovp, dataset0, (0, 1), 2, rtol=_rtol).check(dtypes=[int, float, float])

        n2_10 = np.arange(2, 10)
        dataset1 = np.column_stack([n2_10, np.full_like(n2_10, 1.0), np.full_like(n2_10, 0)])
        FuncData(_smirnovp, dataset1, (0, 1), 2, rtol=_rtol).check(dtypes=[int, float, float])