def test_subclass(self):
     a = array(1).view(PhysicalQuantity2)
     b = array(7).view(PhysicalQuantity2)
     gs = geomspace(a, b)
     assert type(gs) is PhysicalQuantity2
     assert_equal(gs, geomspace(1.0, 7.0))
     gs = geomspace(a, b, 1)
     assert type(gs) is PhysicalQuantity2
     assert_equal(gs, geomspace(1.0, 7.0, 1))

    def test_dtype(self):
        y = geomspace(1, 1e6, dtype='float32')
        assert_equal(y.dtype, dtype('float32'))
        y = geomspace(1, 1e6, dtype='float64')
        assert_equal(y.dtype, dtype('float64'))
        y = geomspace(1, 1e6, dtype='int32')
        assert_equal(y.dtype, dtype('int32'))

        # Native types
        y = geomspace(1, 1e6, dtype=float)
        assert_equal(y.dtype, dtype('float_'))
        y = geomspace(1, 1e6, dtype=complex)
        assert_equal(y.dtype, dtype('complex'))

 def test_start_stop_array(self):
     # Try to use all special cases.
     start = array([1.e0, 32., 1j, -4j, 1+1j, -1])
     stop = array([1.e4, 2., 16j, -324j, 10000+10000j, 1])
     t1 = geomspace(start, stop, 5)
     t2 = stack([geomspace(_start, _stop, 5)
                 for _start, _stop in zip(start, stop)], axis=1)
     assert_equal(t1, t2)
     t3 = geomspace(start, stop[0], 5)
     t4 = stack([geomspace(_start, stop[0], 5)
                 for _start in start], axis=1)
     assert_equal(t3, t4)
     t5 = geomspace(start, stop, 5, axis=-1)
     assert_equal(t5, t2.T)

def user_bkg_spec_a():
    """Create "user-background" data for testing.

    Returns
    -------
    m_bkg_spec : `~jwst.datamodels.MultiSpecModel`
    """

    # This data type is used for creating a MultiSpecModel.
    spec_dtype = datamodels.SpecModel().spec_table.dtype

    # m_bkg_spec doesn't have to be a MultiSpecModel, but that's an option.
    m_bkg_spec = datamodels.MultiSpecModel()
    wavelength = np.geomspace(1.5, 4.5, num=25, endpoint=True,
                              dtype=np.float64)
    flux = np.linspace(13., 25., num=25, endpoint=True, retstep=False,
                       dtype=np.float64)
    fl_error = np.ones_like(wavelength)
    dq = np.zeros(wavelength.shape, dtype=np.uint32)
    net = np.zeros_like(wavelength)
    nerror = np.ones_like(wavelength)
    background = np.ones_like(wavelength)
    berror = np.ones_like(wavelength)
    npixels = np.ones_like(wavelength)
    otab = np.array(list(zip(wavelength, flux, fl_error, dq,
                             net, nerror, background, berror,
                             npixels)),
                    dtype=spec_dtype)
    spec = datamodels.SpecModel(spec_table=otab)
    m_bkg_spec.spec.append(spec)

    return m_bkg_spec

def bin_dataframe(df, n_bins):
    """
    Assign a "bin" column to the dataframe to indicate which bin the true
    charges belong to.

    Bins are assigned in log space.

    Parameters
    ----------
    df : pd.DataFrame
    n_bins : int
        Number of bins to allow in range

    Returns
    -------
    pd.DataFrame
    """
    true = df['true'].values
    min_ = true.min()
    max_ = true.max()
    bins = np.geomspace(min_, max_, n_bins)
    bins = np.append(bins, 10**(np.log10(bins[-1]) +
                                np.diff(np.log10(bins))[0]))
    df['bin'] = np.digitize(true, bins, right=True) - 1

    return df

def user_bkg_spec_c():
    """Create "user-background" data for testing.

    Returns
    -------
    m_bkg_spec : `~jwst.datamodels.CombinedSpecModel`
    """

    # This is the data type of an output table from combine_1d.
    spec_table_dtype = datamodels.CombinedSpecModel().spec_table.dtype

    wavelength = np.geomspace(1.5, 4.5, num=25, endpoint=True,
                              dtype=np.float64)
    flux = np.linspace(13., 25., num=25, endpoint=True, retstep=False,
                       dtype=np.float64)
    fl_error = np.ones_like(wavelength)
    dq = np.zeros(wavelength.shape, dtype=np.uint32)
    net = np.zeros_like(wavelength)
    weight = np.ones_like(wavelength)
    n_input = np.ones_like(wavelength)                  # yes, float64
    data = np.array(list(zip(wavelength, flux, fl_error, net,
                             dq, weight, n_input)),
                    dtype=spec_table_dtype)
    m_bkg_spec = datamodels.CombinedSpecModel(spec_table=data)

    return m_bkg_spec

def geomspace_stepsize(start, stop, dlnx):
	""" return equally log spaced array from start to stop (sometimes over) with step size dlnx """
	num = np.ceil(np.absolute(np.log(stop) - np.log(start))/dlnx) + 1

	if stop > start:
		stop_real = np.exp(np.log(start)+(num-1)*dlnx)
	else: 
		stop_real = np.exp(np.log(start)-(num-1)*dlnx)

	return np.geomspace(start, stop_real, num=num, endpoint=True)

    def get_reynolds(self, velocity):
        Re = self.foil.Reynolds(velocity)
        re_space = np.round(np.geomspace(30000, 2e6, 20), -4)
        idx = np.argmin(abs(re_space - Re))
        reynolds = re_space[idx] # np.round(Re, -4)  # Round to nearest 1000

        if (reynolds < 30000.0):
            reynolds = 30000.0
            
        return reynolds

 def evaluate(self):
     var_reports = {}
     for variable in self.variables:
         initial_value = self.tf_session.run(variable)
         start_value = max(np.abs(initial_value/100), .01)
         base_space = np.geomspace(start_value, np.abs(initial_value), self.width // 2)
         space = np.hstack((initial_value + base_space, [initial_value], initial_value - base_space))
         reports = []
         for test_value in space:
             self._assign_value(variable.name, test_value)
             reports.append(self.evaluator.evaluate())
         self._assign_value(variable.name, initial_value)
         var_reports[variable.name] = reports
     return SensitivityReport(var_reports)

    def test_scott_vs_stone(self):
        """Verify that Scott's rule and Stone's rule converges for normally distributed data"""

        def nbins_ratio(seed, size):
            rng = np.random.RandomState(seed)
            x = rng.normal(loc=0, scale=2, size=size)
            a, b = len(np.histogram(x, 'stone')[0]), len(np.histogram(x, 'scott')[0])
            return a / (a + b)

        ll = [[nbins_ratio(seed, size) for size in np.geomspace(start=10, stop=100, num=4).round().astype(int)]
              for seed in range(10)]

        # the average difference between the two methods decreases as the dataset size increases.
        avg = abs(np.mean(ll, axis=0) - 0.5)
        assert_almost_equal(avg, [0.15, 0.09, 0.08, 0.03], decimal=2)

    def __init__(self, f_max, mu_min, population_distribution,
                 delta_fitness, mu_multiple, fraction_beneficial,
                 fraction_accurate, fraction_mu2mu, K):
        """
        The population is described both by its state and the parameters of the
        model. Population_distribution should be a single number or a numpy
        array.
        """

        self.population_distribution = \
            np.atleast_2d(np.array(population_distribution, dtype='int64'))
        if np.any(self.population_distribution) < 0:
            raise ValueError('the population distribution must be nonnegative')
        pop_shape = self.population_distribution.shape

        self.delta_fitness = delta_fitness
        if self.delta_fitness <= 0:
            raise ValueError('delta fitness must be positive')

        self.mu_multiple = mu_multiple
        if self.mu_multiple <= 1:
            raise ValueError('mu multiple must be greater than one')

        self.fraction_beneficial = fraction_beneficial
        if self.fraction_beneficial >= 1 or self.fraction_beneficial < 0:
            raise ValueError('fraction beneficial must be >= 0 and < 1')

        self.fraction_accurate = fraction_accurate
        if self.fraction_accurate >= 1 or self.fraction_accurate < 0:
            raise ValueError('fraction accurate must be >=0 and < 1')

        self.fraction_mu2mu = fraction_mu2mu
        if self.fraction_mu2mu >= 1 or self.fraction_mu2mu < 0:
            raise ValueError('fraction_mu2mu must be >=0 and < 1')

        self.pop_cap = K
        if self.pop_cap < 100:
            raise ValueError('pop_cap must be greater than or equal to 100')

        f_min = f_max - delta_fitness*(pop_shape[0]-1)
        self.fitness_list = np.transpose(np.atleast_2d(np.linspace(f_max,
                                         f_min, pop_shape[0])))

        self.mutation_list = np.geomspace(mu_min,
                                          mu_min*mu_multiple**(pop_shape[1]-1),
                                          pop_shape[1])

 def __init__(
         self,
         start_center_hz,
         stop_center_hz,
         bandwidth_ratio,
         n_bands,
         always_even=False):
     self.__bands = [
         FrequencyBand.from_center(cf, cf * bandwidth_ratio)
         for cf in np.geomspace(start_center_hz, stop_center_hz, num=n_bands)
         ]
     band = FrequencyBand(self.__bands[0].start_hz, self.__bands[-1].stop_hz)
     super(GeometricScale, self).__init__(
         band, n_bands, always_even=always_even)
     self.start_center_hz = start_center_hz
     self.stop_center_hz = stop_center_hz
     self.bandwidth_ratio = bandwidth_ratio

    def test_array_scalar(self):
        lim1 = array([120, 100], dtype="int8")
        lim2 = array([-120, -100], dtype="int8")
        lim3 = array([1200, 1000], dtype="uint16")
        t1 = geomspace(lim1[0], lim1[1], 5)
        t2 = geomspace(lim2[0], lim2[1], 5)
        t3 = geomspace(lim3[0], lim3[1], 5)
        t4 = geomspace(120.0, 100.0, 5)
        t5 = geomspace(-120.0, -100.0, 5)
        t6 = geomspace(1200.0, 1000.0, 5)

        # t3 uses float32, t6 uses float64
        assert_allclose(t1, t4, rtol=1e-2)
        assert_allclose(t2, t5, rtol=1e-2)
        assert_allclose(t3, t6, rtol=1e-5)

def user_bkg_spec_b():
    """Create "user-background" data for testing.

    `expand_to_2d` uses `np.interp` for interpolation, and the wavelength
    array that is passed to `np.interp` must be increasing.  `expand_to_2d`
    is supposed to handle the case that the wavelengths are decreasing.
    Create data for checking that the results are the same even if the
    wavelength array in `m_bkg_spec` is reversed so that the values are
    decreasing, and the corresponding flux array is also reversed to retain
    the original (wavelength, flux) relation.

    Returns
    -------
    m_bkg_spec : `~jwst.datamodels.MultiSpecModel`
    """

    # This data type is used for creating a MultiSpecModel.
    spec_dtype = datamodels.SpecModel().spec_table.dtype

    m_bkg_spec = datamodels.MultiSpecModel()
    wavelength = np.geomspace(1.5, 4.5, num=25, endpoint=True,
                              dtype=np.float64)[::-1]
    flux = np.linspace(13., 25., num=25, endpoint=True, retstep=False,
                       dtype=np.float64)[::-1]
    fl_error = np.ones_like(wavelength)
    dq = np.zeros(wavelength.shape, dtype=np.uint32)
    net = np.zeros_like(wavelength)
    nerror = np.ones_like(wavelength)
    background = np.ones_like(wavelength)
    berror = np.ones_like(wavelength)
    npixels = np.ones_like(wavelength)
    otab = np.array(list(zip(wavelength, flux, fl_error, dq,
                             net, nerror, background, berror,
                             npixels)),
                    dtype=spec_dtype)
    spec = datamodels.SpecModel(spec_table=otab)
    m_bkg_spec.spec.append(spec)

    return m_bkg_spec

    def overview(self):
        """
        Plots an overview figure.
        """
        is_nm = self.freq_unit is "nm"
        if is_nm:
            ph.vis_mode()
        else:
            ph.ir_mode()
        ds = self.dataset
        x = ds.wavelengths if is_nm else ds.wavenumbers
        fig, axs = plt.subplots(
            3, 1, figsize=(5, 12), gridspec_kw=dict(height_ratios=(2, 1, 1))
        )
        self.map(ax=axs[0])

        times = np.hstack((0, np.geomspace(0.1, ds.t.max(), 6)))
        sp = self.spec(times, ax=axs[1])
        freqs = np.unique(np.linspace(x.min(), x.max(), 6))
        tr = self.trans(freqs, ax=axs[2])
        OverviewPlot = namedtuple("OverviewPlot", "fig axs trans spec")
        return OverviewPlot(fig, axs, tr, sp)

    def test_basic(self):
        y = geomspace(1, 1e6)
        assert_(len(y) == 50)
        y = geomspace(1, 1e6, num=100)
        assert_(y[-1] == 10 ** 6)
        y = geomspace(1, 1e6, endpoint=False)
        assert_(y[-1] < 10 ** 6)
        y = geomspace(1, 1e6, num=7)
        assert_array_equal(y, [1, 10, 100, 1e3, 1e4, 1e5, 1e6])

        y = geomspace(8, 2, num=3)
        assert_allclose(y, [8, 4, 2])
        assert_array_equal(y.imag, 0)

        y = geomspace(-1, -100, num=3)
        assert_array_equal(y, [-1, -10, -100])
        assert_array_equal(y.imag, 0)

        y = geomspace(-100, -1, num=3)
        assert_array_equal(y, [-100, -10, -1])
        assert_array_equal(y.imag, 0)

 def test_physical_quantities(self):
     a = PhysicalQuantity(1.0)
     b = PhysicalQuantity(5.0)
     assert_equal(geomspace(a, b), geomspace(1.0, 5.0))

log = get_logger('clusters')

parser = argparse.ArgumentParser(description='Save clusters membership in mongodb')
parser.add_argument('--num', type=int, default=50, help='Specify the number of cluster levels [default: %(default)s]')
args = vars(parser.parse_args())

uri = uri_parser.parse_uri(MONGODB_URL)
client = MongoClient(uri['nodelist'][0][0], uri['nodelist'][0][1])
db = client.windmobile

now = datetime.now().timestamp()
all_stations = list(db.stations.find({
    'status': {'$ne': 'hidden'},
    'last._id': {'$gt': now - 30 * 24 * 3600}
}))
range_clusters = np.geomspace(20, len(all_stations), num=args['num'], dtype=np.int)

ids = np.array([station['_id'] for station in all_stations])

x = [station['loc']['coordinates'][0] for station in all_stations]
y = [station['loc']['coordinates'][1] for station in all_stations]
X = np.array((x, y))
X = X.T


try:
    mongo_bulk = db.stations.initialize_ordered_bulk_op()
    mongo_bulk.find({}).update({'$set': {'clusters': []}})

    for n_clusters in reversed(range_clusters):

    def __init__(self, f_max, mu_min, population_distribution,
                 delta_fitness, mu_multiple, fraction_beneficial,
                 fraction_accurate, fraction_mu2mu, K):
        """
        The population is described both by its state and the parameters of the
        model. Population_distribution should be a single number or a numpy
        array.
        """

        self.population_distribution = \
            np.atleast_2d(np.array(population_distribution, dtype='int64'))
        if np.any(self.population_distribution) < 0:
            raise ValueError('the population distribution must be nonnegative')
        pop_shape = self.population_distribution.shape

        self.delta_fitness = delta_fitness
        if self.delta_fitness <= 0:
            raise ValueError('delta fitness must be positive')

        self.mu_multiple = mu_multiple
        if self.mu_multiple <= 1:
            raise ValueError('mu multiple must be greater than one')

        self.fraction_beneficial = fraction_beneficial
        if self.fraction_beneficial >= 1 or self.fraction_beneficial < 0:
            raise ValueError('fraction beneficial must be >= 0 and < 1')

        self.fraction_accurate = fraction_accurate
        if self.fraction_accurate >= 1 or self.fraction_accurate < 0:
            raise ValueError('fraction accurate must be >=0 and < 1')

        self.fraction_mu2mu = fraction_mu2mu
        if self.fraction_mu2mu >= 1 or self.fraction_mu2mu < 0:
            raise ValueError('fraction_mu2mu must be >=0 and < 1')

        self.pmus = \
            np.array([fraction_mu2mu*fraction_accurate,
                     fraction_mu2mu*(1-fraction_accurate),
                     (1-fraction_mu2mu)*fraction_beneficial,
                     (1-fraction_mu2mu)*(1-fraction_beneficial)]).reshape(4,
                                                                          1, 1)

        self.pop_cap = K
        if self.pop_cap < 100:
            raise ValueError('pop_cap must be greater than or equal to 100')
        if K <= 10**9:
            self.wright_fisher = wf.wright_fisher_fitness_update
            self.multinomial = arrm.array_multinomial
        else:
            self.wright_fisher = wf.wright_fisher_fitness_update_bigN
            self.multinomial = arrm.array_multinomial_int64

        f_min = f_max - delta_fitness*(pop_shape[0]-1)
        self.fitness_list = np.transpose(np.atleast_2d(np.linspace(f_max,
                                         f_min, pop_shape[0])))

        self.mutation_list = np.minimum(np.geomspace(mu_min,
                                          mu_min*mu_multiple**(pop_shape[1]-1),
                                          pop_shape[1]),1)
        
        if self.mutation_list.size == 1:
            if self.mutation_list[0] > 1.0:
                raise ValueError('Your population distribution implies mutation'
                                 ' rates exceeding one. This is not possible in'
                                 ' this model.')
        #else:
        #    if self.mutation_list[-2] >= 1.0:
        #        raise ValueError('Your population distribution implies mutation'
        #                         ' rates exceeding one. This is not possible in'
        #                         ' this model.')

        self.stencil = np.array([[0, -1],
                                 [0, 1],
                                 [-1, 0],
                                 [1, 0],
                                 [0, 0]])
        self.summer = stsum.fixedStencilSum(3, 0, (5,), self.stencil)

    def test_complex(self):
        # Purely imaginary
        y = geomspace(1j, 16j, num=5)
        assert_allclose(y, [1j, 2j, 4j, 8j, 16j])
        assert_array_equal(y.real, 0)

        y = geomspace(-4j, -324j, num=5)
        assert_allclose(y, [-4j, -12j, -36j, -108j, -324j])
        assert_array_equal(y.real, 0)

        y = geomspace(1+1j, 1000+1000j, num=4)
        assert_allclose(y, [1+1j, 10+10j, 100+100j, 1000+1000j])

        y = geomspace(-1+1j, -1000+1000j, num=4)
        assert_allclose(y, [-1+1j, -10+10j, -100+100j, -1000+1000j])

        # Logarithmic spirals
        y = geomspace(-1, 1, num=3, dtype=complex)
        assert_allclose(y, [-1, 1j, +1])

        y = geomspace(0+3j, -3+0j, 3)
        assert_allclose(y, [0+3j, -3/sqrt(2)+3j/sqrt(2), -3+0j])
        y = geomspace(0+3j, 3+0j, 3)
        assert_allclose(y, [0+3j, 3/sqrt(2)+3j/sqrt(2), 3+0j])
        y = geomspace(-3+0j, 0-3j, 3)
        assert_allclose(y, [-3+0j, -3/sqrt(2)-3j/sqrt(2), 0-3j])
        y = geomspace(0+3j, -3+0j, 3)
        assert_allclose(y, [0+3j, -3/sqrt(2)+3j/sqrt(2), -3+0j])
        y = geomspace(-2-3j, 5+7j, 7)
        assert_allclose(y, [-2-3j, -0.29058977-4.15771027j,
                            2.08885354-4.34146838j, 4.58345529-3.16355218j,
                            6.41401745-0.55233457j, 6.75707386+3.11795092j,
                            5+7j])

        # Type promotion should prevent the -5 from becoming a NaN
        y = geomspace(3j, -5, 2)
        assert_allclose(y, [3j, -5])
        y = geomspace(-5, 3j, 2)
        assert_allclose(y, [-5, 3j])