def __init__(self, capacity=100, cost=100, number=None):

        Vehicle = namedtuple("Vehicle", ["index", "capacity", "cost"])

        if number is None:
            self.number = np.size(capacity)
        else:
            self.number = number
        idxs = np.array(range(0, self.number))

        if np.isscalar(capacity):
            capacities = capacity * np.ones_like(idxs)
        elif np.size(capacity) != np.size(capacity):
            print("capacity is neither scalar, nor the same size as num!")
        else:
            capacities = capacity

        if np.isscalar(cost):
            costs = cost * np.ones_like(idxs)
        elif np.size(cost) != self.number:
            print(np.size(cost))
            print("cost is neither scalar, nor the same size as num!")
        else:
            costs = cost

        self.vehicles = [Vehicle(idx, capacity, cost) for idx, capacity, cost in zip(idxs, capacities, costs)]

def get_dummy_particles():
    x, y = numpy.mgrid[-5 * dx : box_length + 5 * dx + 1e-10 : dx, -5 * dx : box_height + 5 * dx + 1e-10 : dx]

    xd, yd = x.ravel(), y.ravel()

    md = numpy.ones_like(xd) * m
    hd = numpy.ones_like(xd) * h

    rhod = numpy.ones_like(xd) * ro
    cd = numpy.ones_like(xd) * co
    pd = numpy.zeros_like(xd)

    dummy_fluid = base.get_particle_array(name="dummy_fluid", type=Fluid, x=xd, y=yd, h=hd, rho=rhod, c=cd, p=pd)

    # remove indices within the square

    indices = []

    np = dummy_fluid.get_number_of_particles()
    x, y = dummy_fluid.get("x", "y")

    for i in range(np):
        if -dx / 2 <= x[i] <= box_length + dx / 2:
            if -dx / 2 <= y[i] <= box_height + dx / 2:
                indices.append(i)

    to_remove = base.LongArray(len(indices))
    to_remove.set_data(numpy.array(indices))

    dummy_fluid.remove_particles(to_remove)

    return dummy_fluid

def test_arithmetic_overload_ccddata_operand(ccd_data):
    ccd_data.uncertainty = StdDevUncertainty(np.ones_like(ccd_data))
    operand = ccd_data.copy()
    result = ccd_data.add(operand)
    assert len(result.meta) == 0
    np.testing.assert_array_equal(result.data,
                                  2 * ccd_data.data)
    np.testing.assert_array_equal(result.uncertainty.array,
                                  np.sqrt(2) * ccd_data.uncertainty.array)

    result = ccd_data.subtract(operand)
    assert len(result.meta) == 0
    np.testing.assert_array_equal(result.data,
                                  0 * ccd_data.data)
    np.testing.assert_array_equal(result.uncertainty.array,
                                  np.sqrt(2) * ccd_data.uncertainty.array)

    result = ccd_data.multiply(operand)
    assert len(result.meta) == 0
    np.testing.assert_array_equal(result.data,
                                  ccd_data.data ** 2)
    expected_uncertainty = (np.sqrt(2) * np.abs(ccd_data.data) *
                            ccd_data.uncertainty.array)
    np.testing.assert_allclose(result.uncertainty.array,
                               expected_uncertainty)

    result = ccd_data.divide(operand)
    assert len(result.meta) == 0
    np.testing.assert_array_equal(result.data,
                                  np.ones_like(ccd_data.data))
    expected_uncertainty = (np.sqrt(2) / np.abs(ccd_data.data) *
                            ccd_data.uncertainty.array)
    np.testing.assert_allclose(result.uncertainty.array,
                               expected_uncertainty)

def fix_chip_wavelength(model_orders, data_orders, band_cutoff=1870):
    """ Adjust the wavelength in data_orders to be self-consistent
    """
    # H band
    model_orders_H = [o.copy() for o in model_orders if o.x[-1] < band_cutoff]
    data_orders_H = [o.copy() for o in data_orders if o.x[-1] < band_cutoff]
    ordernums_H = 121.0 - np.arange(len(model_orders_H))
    p_H = fit_wavelength(model_orders_H, ordernums_H, first_order=3, last_order=len(ordernums_H) - 4)

    # K band
    model_orders_K = [o.copy() for o in model_orders if o.x[-1] > band_cutoff]
    data_orders_K = [o.copy() for o in data_orders if o.x[-1] > band_cutoff]
    ordernums_K = 92.0 - np.arange(len(model_orders_K))
    p_K = fit_wavelength(model_orders_K, ordernums_K, first_order=7, last_order=len(ordernums_K) - 4)

    new_orders = []
    for i, order in enumerate(data_orders):
        pixels = np.arange(order.size(), dtype=np.float)
        if order.x[-1] < band_cutoff:
            # H band
            ordernum = ordernums_H[i] * np.ones_like(pixels)
            wave = p_H(pixels, ordernum) / ordernum
        else:
            # K band
            ordernum = ordernums_K[i-len(ordernums_H)] * np.ones_like(pixels)
            wave = p_K(pixels, ordernum) / ordernum
            
        new_orders.append(DataStructures.xypoint(x=wave, y=order.y, cont=order.cont, err=order.err))
    return new_orders

def get_fluid():
    """ Get the fluid particle array """

    x, y = numpy.mgrid[dx : box_length - 1e-10 : dx, dx : box_height - 1e-10 : dx]

    xf, yf = x.ravel(), y.ravel()

    mf = numpy.ones_like(xf) * m
    hf = numpy.ones_like(xf) * h

    rhof = numpy.ones_like(xf) * ro
    cf = numpy.ones_like(xf) * co
    pf = numpy.zeros_like(xf)

    fluid = base.get_particle_array(name="fluid", type=Fluid, x=xf, y=yf, h=hf, rho=rhof, c=cf, p=pf)

    # remove indices within the square

    indices = []

    np = fluid.get_number_of_particles()
    x, y = fluid.get("x", "y")

    for i in range(np):
        if 1.0 - dx / 2 <= x[i] <= 2.0 + dx / 2:
            if 2.0 - dx / 2 <= y[i] <= 3.0 + dx / 2:
                indices.append(i)

    to_remove = base.LongArray(len(indices))
    to_remove.set_data(numpy.array(indices))

    fluid.remove_particles(to_remove)

    return fluid

    def noise_variance_feedpairs(self, fi, fj, f_indices, nt_per_day, ndays=None):
        ndays = self.ndays if not ndays else ndays # Set to value if not set.
        t_int = ndays * units.t_sidereal / nt_per_day
        # bw = 1.0e6 * (self.freq_upper - self.freq_lower) / self.num_freq
        bw = np.abs(self.frequencies[1] - self.frequencies[0]) * 1e6

        return np.ones_like(fi) * np.ones_like(fj) * 2.0*self.tsys(f_indices)**2 / (t_int * bw) # 2.0 for two pol

    def prime_to_pixel(self, xprime, yprime,  color=0):
        color0 = self._get_ricut()
        g0, g1, g2, g3 = self._get_drow()
        h0, h1, h2, h3 = self._get_dcol()
        px, py, qx, qy = self._get_cscc()

        # #$(%*&^(%$%*& bad documentation.
        (px,py) = (py,px)
        (qx,qy) = (qy,qx)

        qx = qx * np.ones_like(xprime)
        qy = qy * np.ones_like(yprime)
        xprime -= np.where(color < color0, px * color, qx)
        yprime -= np.where(color < color0, py * color, qy)

        # Now invert:
        #   yprime = y + g0 + g1 * x + g2 * x**2 + g3 * x**3
        #   xprime = x + h0 + h1 * x + h2 * x**2 + h3 * x**3
        x = xprime - h0
        # dumb-ass Newton's method
        dx = 1.
        # FIXME -- should just update the ones that aren't zero
        # FIXME -- should put in some failsafe...
        while np.max(np.abs(np.atleast_1d(dx))) > 1e-10:
            xp    = x + h0 + h1 * x + h2 * x**2 + h3 * x**3
            dxpdx = 1 +      h1     + h2 * 2*x +  h3 * 3*x**2
            dx = (xprime - xp) / dxpdx
            x += dx
        y = yprime - (g0 + g1 * x + g2 * x**2 + g3 * x**3)
        return (x, y)

def histgram_3D(data):
    '''
    入力された二次元配列を3Dhistgramとして表示する
    '''
    from mpl_toolkits.mplot3d import Axes3D
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    
    x = data[:,0]
    y = data[:,1]

    hist, xedges, yedges = np.histogram2d(x, y, bins=30)
    X, Y = np.meshgrid(xedges[:-1] + 0.25, yedges[:-1] + 0.25)

    # bar3dでは行にする
    X = X.flatten()
    Y = Y.flatten()
    Z = np.zeros(len(X))

    # 表示するバーの太さ
    dx = (xedges[1] - xedges[0]) * np.ones_like(Z)
    dy = (yedges[1] - yedges[0]) * np.ones_like(Z)
    dz = hist.flatten() # これはそのままでok

    # 描画
    ax.bar3d(X, Y, Z, dx, dy, dz, color='b', zsort='average')

    def isclose(a, b, rtol=1e-05, atol=1e-08, equal_nan=False):
        def within_tol(x, y, atol, rtol):
            result = np.less_equal(np.abs(x-y), atol + rtol * np.abs(y))
            if np.isscalar(a) and np.isscalar(b):
                result = np.bool(result)
            return result

        x = np.array(a, copy=False, subok=True, ndmin=1)
        y = np.array(b, copy=False, subok=True, ndmin=1)
        xfin = np.isfinite(x)
        yfin = np.isfinite(y)
        if np.all(xfin) and np.all(yfin):
            return within_tol(x, y, atol, rtol)
        else:
            finite = xfin & yfin
            cond = np.zeros_like(finite, subok=True)
            # Because we're using boolean indexing, x & y must be the same shape.
            # Ideally, we'd just do x, y = broadcast_arrays(x, y). It's in
            # lib.stride_tricks, though, so we can't import it here.
            x = x * np.ones_like(cond)
            y = y * np.ones_like(cond)
            # Avoid subtraction with infinite/nan values...
            cond[finite] = within_tol(x[finite], y[finite], atol, rtol)
            # Check for equality of infinite values...
            cond[~finite] = (x[~finite] == y[~finite])
            if equal_nan:
                # Make NaN == NaN
                cond[np.isnan(x) & np.isnan(y)] = True
            return cond

    def _prepare_sw_arguments(self, ncol, nlay):
        aldif = _climlab_to_rrtm_sfc(self.aldif * np.ones_like(self.Ts))
        aldir = _climlab_to_rrtm_sfc(self.aldir * np.ones_like(self.Ts))
        asdif = _climlab_to_rrtm_sfc(self.asdif * np.ones_like(self.Ts))
        asdir = _climlab_to_rrtm_sfc(self.asdir * np.ones_like(self.Ts))
        coszen = _climlab_to_rrtm_sfc(self.coszen * np.ones_like(self.Ts))
        #  THE REST OF THESE ARGUMENTS ARE STILL BEING HARD CODED.
        #   NEED TO FIX THIS UP...

        #  These arrays have an extra dimension for number of bands
        dim_sw1 = [nbndsw,ncol,nlay]     # [nbndsw,ncol,nlay]
        dim_sw2 = [ncol,nlay,nbndsw]  # [ncol,nlay,nbndsw]
        tauc = np.zeros(dim_sw1) # In-cloud optical depth
        ssac = np.zeros(dim_sw1) # In-cloud single scattering albedo
        asmc = np.zeros(dim_sw1) # In-cloud asymmetry parameter
        fsfc = np.zeros(dim_sw1) # In-cloud forward scattering fraction (delta function pointing forward "forward peaked scattering")

        # AEROSOLS
        tauaer = np.zeros(dim_sw2)   # Aerosol optical depth (iaer=10 only), Dimensions,  (ncol,nlay,nbndsw)] #  (non-delta scaled)
        ssaaer = np.zeros(dim_sw2)   # Aerosol single scattering albedo (iaer=10 only), Dimensions,  (ncol,nlay,nbndsw)] #  (non-delta scaled)
        asmaer = np.zeros(dim_sw2)   # Aerosol asymmetry parameter (iaer=10 only), Dimensions,  (ncol,nlay,nbndsw)] #  (non-delta scaled)
        ecaer  = np.zeros([ncol,nlay,naerec])   # Aerosol optical depth at 0.55 micron (iaer=6 only), Dimensions,  (ncol,nlay,naerec)] #  (non-delta scaled)

        return (aldif,aldir,asdif,asdir,coszen,tauc,ssac,asmc,
                fsfc,tauaer,ssaaer,asmaer,ecaer)

 def _scalars_changed(self, s):
     self.dataset.point_data.scalars = s
     self.dataset.point_data.scalars.name = 'scalars'
     self.set(vectors=np.c_[np.ones_like(s),
                               np.ones_like(s),
                               s])
     self.update()

    def set_jds(self, val1, val2):
        self._check_scale(self._scale)  # Validate scale.

        sum12, err12 = two_sum(val1, val2)
        iy_start = np.trunc(sum12).astype(np.int)
        extra, y_frac = two_sum(sum12, -iy_start)
        y_frac += extra + err12

        val = (val1 + val2).astype(np.double)
        iy_start = np.trunc(val).astype(np.int)

        imon = np.ones_like(iy_start)
        iday = np.ones_like(iy_start)
        ihr = np.zeros_like(iy_start)
        imin = np.zeros_like(iy_start)
        isec = np.zeros_like(y_frac)

        # Possible enhancement: use np.unique to only compute start, stop
        # for unique values of iy_start.
        scale = self.scale.upper().encode('ascii')
        jd1_start, jd2_start = erfa.dtf2d(scale, iy_start, imon, iday,
                                          ihr, imin, isec)
        jd1_end, jd2_end = erfa.dtf2d(scale, iy_start + 1, imon, iday,
                                      ihr, imin, isec)

        t_start = Time(jd1_start, jd2_start, scale=self.scale, format='jd')
        t_end = Time(jd1_end, jd2_end, scale=self.scale, format='jd')
        t_frac = t_start + (t_end - t_start) * y_frac

        self.jd1, self.jd2 = day_frac(t_frac.jd1, t_frac.jd2)

def plot(y, title, t):
    
        a = y[0]
        b = y[1]
        
        print a
        
        if a and b:
        
            bins = numpy.linspace(min(a+b), max(a+b), 20)
            pyplot.clf()
            
            if a:
                w0 = numpy.ones_like(a)/float(len(a))
                pyplot.hist(a, bins, weights=w0,alpha=0.5, color='r', histtype='stepfilled', label='link')
            
            if b:
                w1 = numpy.ones_like(b)/float(len(b))
                pyplot.hist(b, bins,weights=w1, alpha=0.5, color='b', histtype='stepfilled', label='no link')
            
            pyplot.title(title)
            pyplot.ylabel("Fraction over population")
            pyplot.xlabel("Similarity")
            pyplot.legend();
            #plt.savefig("/Users/spoulson/Dropbox/my_papers/figs/"+title.replace(' ','_')+'_'+ str(t) +'.png')
            pyplot.show()

def get_circular_patch(name="", type=0, dx=0.05):
    
    x,y = numpy.mgrid[-1.05:1.05+1e-4:dx, -1.05:1.05+1e-4:dx]
    x = x.ravel()
    y = y.ravel()
 
    m = numpy.ones_like(x)*dx*dx
    h = numpy.ones_like(x)*2*dx
    rho = numpy.ones_like(x)

    p = 0.5*1.0*100*100*(1 - (x**2 + y**2))

    cs = numpy.ones_like(x) * 100.0

    u = 0*x
    v = 0*y

    indices = []

    for i in range(len(x)):
        if numpy.sqrt(x[i]*x[i] + y[i]*y[i]) - 1 > 1e-10:
            indices.append(i)
            
    pa = base.get_particle_array(x=x, y=y, m=m, rho=rho, h=h, p=p, u=u, v=v,
                                 cs=cs,name=name, type=type)

    la = base.LongArray(len(indices))
    la.set_data(numpy.array(indices))

    pa.remove_particles(la)

    pa.set(idx=numpy.arange(len(pa.x)))

    return pa

    def siggen_model(s, rad, phi, z, e, temp, num_1, num_2, num_3, den_1, den_2, den_3):
      out = np.zeros_like(data)
      
      detector.SetTemperature(temp)
      siggen_wf= detector.GetSiggenWaveform(rad, phi, z, energy=2600)

      if siggen_wf is None:
        return np.ones_like(data)*-1.
      if np.amax(siggen_wf) == 0:
        print "wtf is even happening here?"
        return np.ones_like(data)*-1.
      siggen_wf = np.pad(siggen_wf, (detector.zeroPadding,0), 'constant', constant_values=(0, 0))

      num = [num_1, num_2, num_3]
      den = [1,   den_1, den_2, den_3]
#      num = [-1.089e10,  5.863e17,  6.087e15]
#      den = [1,  3.009e07, 3.743e14,5.21e18]
      system = signal.lti(num, den)
      t = np.arange(0, len(siggen_wf)*10E-9, 10E-9)
      tout, siggen_wf, x = signal.lsim(system, siggen_wf, t)
      siggen_wf /= np.amax(siggen_wf)
      
      siggen_data = siggen_wf[detector.zeroPadding::]
      
      siggen_data = siggen_data*e
      
      out[s:] = siggen_data[0:(len(data) - s)]

      return out

    def test_ewma(self):
        ewma = EWMAVariance()

        sv = ewma.starting_values(self.resids)
        assert_equal(sv.shape[0], ewma.num_params)

        bounds = ewma.bounds(self.resids)
        assert_equal(len(bounds), 0)
        var_bounds = ewma.variance_bounds(self.resids)
        backcast = ewma.backcast(self.resids)
        parameters = np.array([])

        names = ewma.parameter_names()
        names_target = []
        assert_equal(names, names_target)

        ewma.compute_variance(parameters, self.resids, self.sigma2,
                              backcast, var_bounds)
        cond_var_direct = np.zeros_like(self.sigma2)
        parameters = np.array([0.0, 0.06, 0.94])
        rec.garch_recursion(parameters,
                            self.resids ** 2.0,
                            np.sign(self.resids),
                            cond_var_direct,
                            1, 0, 1, self.T, backcast, var_bounds)
        # sigma3 = np.zeros_like(self.sigma2)
        # sigma3[0] = backcast
        # for t in range(1,self.T):
        # sigma3[t] = 0.94 * sigma3[t-1] + 0.06 * self.resids[t-1]**2.0

        assert_allclose(self.sigma2 / cond_var_direct,
                        np.ones_like(self.sigma2))

        A, b = ewma.constraints()
        A_target = np.empty((0, 0))
        b_target = np.empty((0,))
        assert_array_equal(A, A_target)
        assert_array_equal(b, b_target)
        state = np.random.get_state()
        rng = Normal()
        sim_data = ewma.simulate(parameters, self.T, rng.simulate([]))
        np.random.set_state(state)
        e = np.random.standard_normal(self.T + 500)
        initial_value = 1.0

        sigma2 = np.zeros(self.T + 500)
        data = np.zeros(self.T + 500)
        sigma2[0] = initial_value
        data[0] = np.sqrt(initial_value)
        for t in range(1, self.T + 500):
            sigma2[t] = 0.94 * sigma2[t - 1] + 0.06 * data[t - 1] ** 2.0
            data[t] = e[t] * np.sqrt(sigma2[t])

        data = data[500:]
        sigma2 = sigma2[500:]
        assert_almost_equal(data - sim_data[0] + 1.0, np.ones_like(data))
        assert_almost_equal(sigma2 / sim_data[1], np.ones_like(sigma2))

        assert_equal(ewma.num_params, 0)
        assert_equal(ewma.name, 'EWMA/RiskMetrics')

def logit(prop, max_events=None):
    """Convert proportion (expressed in the range [0, 1]) to logit.

    Parameters
    ----------
    prop : float | array-like
        the occurrence proportion.
    max_events : int | array-like | None
        the number of events used to calculate ``prop``. Used in a correction
        factor for cases when ``prop`` is 0 or 1, to prevent returning ``inf``.
        If ``None``, no correction is done, and ``inf`` or ``-inf`` may result.

    Returns
    -------
    lgt : ``numpy.ndarray``, with shape matching ``numpy.array(prop).shape``.
    """
    prop = np.atleast_1d(prop).astype(float)
    if np.any([prop > 1, prop < 0]):
        raise ValueError('Proportions must be in the range [0, 1].')
    if max_events is not None:
        # add equivalent of half an event to 0s, and subtract same from 1s
        max_events = np.atleast_1d(max_events) * np.ones_like(prop)
        corr_factor = 0.5 / max_events
        for loc in zip(*np.where(prop == 0)):
            prop[loc] = corr_factor[loc]
        for loc in zip(*np.where(prop == 1)):
            prop[loc] = 1 - corr_factor[loc]
    return np.log(prop / (np.ones_like(prop) - prop))

def reflective_transformation(y, lb, ub):
    """Compute reflective transformation and its gradient."""
    if in_bounds(y, lb, ub):
        return y, np.ones_like(y)

    lb_finite = np.isfinite(lb)
    ub_finite = np.isfinite(ub)

    x = y.copy()
    g_negative = np.zeros_like(y, dtype=bool)

    mask = lb_finite & ~ub_finite
    x[mask] = np.maximum(y[mask], 2 * lb[mask] - y[mask])
    g_negative[mask] = y[mask] < lb[mask]

    mask = ~lb_finite & ub_finite
    x[mask] = np.minimum(y[mask], 2 * ub[mask] - y[mask])
    g_negative[mask] = y[mask] > ub[mask]

    mask = lb_finite & ub_finite
    d = ub - lb
    t = np.remainder(y[mask] - lb[mask], 2 * d[mask])
    x[mask] = lb[mask] + np.minimum(t, 2 * d[mask] - t)
    g_negative[mask] = t > d[mask]

    g = np.ones_like(y)
    g[g_negative] = -1

    return x, g

    def test_blocked(self):
        # test alignments offsets for simd instructions
        # alignments for vz + 2 * (vs - 1) + 1
        for dt, sz in [(np.float32, 11), (np.float64, 7)]:
            for out, inp1, inp2, msg in _gen_alignment_data(dtype=dt,
                                                            type='binary',
                                                            max_size=sz):
                exp1 = np.ones_like(inp1)
                inp1[...] = np.ones_like(inp1)
                inp2[...] = np.zeros_like(inp2)
                assert_almost_equal(np.add(inp1, inp2), exp1, err_msg=msg)
                assert_almost_equal(np.add(inp1, 1), exp1 + 1, err_msg=msg)
                assert_almost_equal(np.add(1, inp2), exp1, err_msg=msg)

                np.add(inp1, inp2, out=out)
                assert_almost_equal(out, exp1, err_msg=msg)

                inp2[...] += np.arange(inp2.size, dtype=dt) + 1
                assert_almost_equal(np.square(inp2),
                                    np.multiply(inp2, inp2),  err_msg=msg)
                assert_almost_equal(np.reciprocal(inp2),
                                    np.divide(1, inp2),  err_msg=msg)

                inp1[...] = np.ones_like(inp1)
                inp2[...] = np.zeros_like(inp2)
                np.add(inp1, 1, out=out)
                assert_almost_equal(out, exp1 + 1, err_msg=msg)
                np.add(1, inp2, out=out)
                assert_almost_equal(out, exp1, err_msg=msg)

    def __test_pass(axis, data, idx):
        # By default, mean aggregation
        dsync = librosa.util.sync(data, idx, axis=axis)
        if data.ndim == 1 or axis == -1:
            assert np.allclose(dsync, 2 * np.ones_like(dsync))
        else:
            assert np.allclose(dsync, data)

        # Explicit mean aggregation
        dsync = librosa.util.sync(data, idx, aggregate=np.mean, axis=axis)
        if data.ndim == 1 or axis == -1:
            assert np.allclose(dsync, 2 * np.ones_like(dsync))
        else:
            assert np.allclose(dsync, data)

        # Max aggregation
        dsync = librosa.util.sync(data, idx, aggregate=np.max, axis=axis)
        if data.ndim == 1 or axis == -1:
            assert np.allclose(dsync, 4 * np.ones_like(dsync))
        else:
            assert np.allclose(dsync, data)

        # Min aggregation
        dsync = librosa.util.sync(data, idx, aggregate=np.min, axis=axis)
        if data.ndim == 1 or axis == -1:
            assert np.allclose(dsync, np.zeros_like(dsync))
        else:
            assert np.allclose(dsync, data)

        # Test for dtype propagation
        assert dsync.dtype == data.dtype