Tool for simple MC validation studies -- specifically to create multiple realizations of
    a satellite given an CompositeIsochrone object, Kernel object, stellar mass (M_sol) for normalization,
    and distance_modulus.
168418
def validateSatellite(config, isochrone, kernel, stellar_mass, distance_modulus, trials=1, debug=False, seed=0):
    
    logger.info('=== Validate Satellite ===')

    config.params['kernel']['params'] = [kernel.r_h] # TODO: Need better solution to update size??
    logger.debug('Using Plummer profile spatial model with half-light radius %.2f deg'%(config.params['kernel']['params'][0]))
    roi = ugali.observation.roi.ROI(config, kernel.lon, kernel.lat)
    simulator = ugali.simulation.simulator.Simulator(config, roi=roi)
    catalog_base = ugali.observation.catalog.Catalog(config,roi=roi)
    mask = ugali.observation.mask.Mask(config, roi)

    coords = (kernel.lon, kernel.lat)
    
    results = {'mc_lon': [],
               'mc_lat': [],
               'mc_distance_modulus': [],
               'mc_stellar_mass': [],
               'mc_radius': [],
               'richness': [],
               'log_likelihood': [],
               'richness_lower': [],
               'richness_upper': [],
               'richness_limit': [],
               'f': [],
               'stellar_mass': []}

    numpy.random.seed(seed)

    for ii in range(0, trials):
        logger.info('=== Running Satellite %i ==='%ii)

        # Simulate
        catalog_satellite = simulator.satellite(isochrone, kernel, stellar_mass, distance_modulus, mc_source_id=1)
        #catalog_bootstrap = catalog_base.bootstrap()
        #catalog_merge = ugali.observation.catalog.mergeCatalogs([catalog_bootstrap, catalog_satellite])
        catalog_background = simulator.satellite(mc_source_id=1)
        catalog_merge = ugali.observation.catalog.mergeCatalogs([catalog_background, catalog_satellite])

        # Analyze
        likelihood = ugali.analysis.likelihood.Likelihood(config, roi, mask, catalog_merge, isochrone, kernel)
                                                               
        likelihood.precomputeGridSearch([distance_modulus])
        richness, log_likelihood, richness_lower, richness_upper, richness_upper_limit, richness_raw, log_likelihood_raw, p, f = likelihood.gridSearch(coords=coords, distance_modulus_index=0)

        results['mc_lon'].append(kernel.lon)
        results['mc_lat'].append(kernel.lat)
        results['mc_distance_modulus'].append(distance_modulus)
        results['mc_stellar_mass'].append(stellar_mass)
        results['mc_radius'].append(kernel.r_h)
        results['richness'].append(richness)
        results['log_likelihood'].append(log_likelihood)
        results['richness_lower'].append(richness_lower)
        results['richness_upper'].append(richness_upper)
        results['richness_limit'].append(richness_upper_limit)
        results['f'].append(f)
        results['stellar_mass'].append(richness * isochrone.stellarMass())

        logger.info('MC Stellar Mass = %.2f, Measured Stellar Mass = %.2f'%(stellar_mass,richness * isochrone.stellarMass()))
        if debug:
            return likelihood, richness, log_likelihood, richness_lower, richness_upper, richness_upper_limit, richness_raw, log_likelihood_raw, p, f

    return results
13690
def build_landmark_op(self):
        
        tasklogger.log_start("landmark operator")
        is_sparse = sparse.issparse(self.kernel)
        # spectral clustering
        tasklogger.log_start("SVD")
        _, _, VT = randomized_svd(self.diff_aff,
                                  n_components=self.n_svd,
                                  random_state=self.random_state)
        tasklogger.log_complete("SVD")
        tasklogger.log_start("KMeans")
        kmeans = MiniBatchKMeans(
            self.n_landmark,
            init_size=3 * self.n_landmark,
            batch_size=10000,
            random_state=self.random_state)
        self._clusters = kmeans.fit_predict(
            self.diff_op.dot(VT.T))
        # some clusters are not assigned
        tasklogger.log_complete("KMeans")

        # transition matrices
        pmn = self._landmarks_to_data()

        # row normalize
        pnm = pmn.transpose()
        pmn = normalize(pmn, norm='l1', axis=1)
        pnm = normalize(pnm, norm='l1', axis=1)
        landmark_op = pmn.dot(pnm)  # sparsity agnostic matrix multiplication
        if is_sparse:
            # no need to have a sparse landmark operator
            landmark_op = landmark_op.toarray()
        # store output
        self._landmark_op = landmark_op
        self._transitions = pnm
        tasklogger.log_complete("landmark operator")
1112285
def calibration(date, satellite):
    
    counts_shift = CountsShift()
    space_measurement = SpaceMeasurement()
    prelaunch = PreLaunch()
    postlaunch = PostLaunch()
    return {
        'counts_shift': counts_shift.coefficient(satellite),
        'space_measurement': space_measurement.coefficient(satellite),
        'prelaunch': prelaunch.coefficient(satellite),
        'postlaunch': postlaunch.coefficient(date, satellite)
    }
315114
def sigma_r2(self, r, kwargs_profile, kwargs_anisotropy, kwargs_light):
        
        if self._mass_profile == 'power_law':
            if self._anisotropy_type == 'r_ani':
                if self._light_profile == 'Hernquist':
                    sigma_r = self.power_law_anisotropy(r, kwargs_profile, kwargs_anisotropy, kwargs_light)
                else:
                    raise ValueError('light profile %s not supported for Jeans solver' % self._light_profile)
            else:
                raise ValueError('anisotropy type %s not implemented in Jeans equation modelling' % self._anisotropy_type)
        else:
            raise ValueError('mass profile type %s not implemented in Jeans solver' % self._mass_profile)
        return sigma_r
877320
def _spherical_kmeans_single_lloyd(
    X,
    n_clusters,
    sample_weight=None,
    max_iter=300,
    init="k-means++",
    verbose=False,
    x_squared_norms=None,
    random_state=None,
    tol=1e-4,
    precompute_distances=True,
):
    
    random_state = check_random_state(random_state)

    sample_weight = _check_sample_weight(X, sample_weight)

    best_labels, best_inertia, best_centers = None, None, None

    # init
    centers = _init_centroids(
        X, n_clusters, init, random_state=random_state, x_squared_norms=x_squared_norms
    )
    if verbose:
        print("Initialization complete")

    # Allocate memory to store the distances for each sample to its
    # closer center for reallocation in case of ties
    distances = np.zeros(shape=(X.shape[0],), dtype=X.dtype)

    # iterations
    for i in range(max_iter):
        centers_old = centers.copy()

        # labels assignment
        # TODO: _labels_inertia should be done with cosine distance
        #       since ||a - b|| = 2(1 - cos(a,b)) when a,b are unit normalized
        #       this doesn't really matter.
        labels, inertia = _labels_inertia(
            X,
            sample_weight,
            x_squared_norms,
            centers,
            precompute_distances=precompute_distances,
            distances=distances,
        )

        # computation of the means
        if sp.issparse(X):
            centers = _k_means._centers_sparse(
                X, sample_weight, labels, n_clusters, distances
            )
        else:
            centers = _k_means._centers_dense(
                X, sample_weight, labels, n_clusters, distances
            )

        # l2-normalize centers (this is the main contibution here)
        centers = normalize(centers)

        if verbose:
            print("Iteration %2d, inertia %.3f" % (i, inertia))

        if best_inertia is None or inertia < best_inertia:
            best_labels = labels.copy()
            best_centers = centers.copy()
            best_inertia = inertia

        center_shift_total = squared_norm(centers_old - centers)
        if center_shift_total <= tol:
            if verbose:
                print(
                    "Converged at iteration %d: "
                    "center shift %e within tolerance %e" % (i, center_shift_total, tol)
                )
            break

    if center_shift_total > 0:
        # rerun E-step in case of non-convergence so that predicted labels
        # match cluster centers
        best_labels, best_inertia = _labels_inertia(
            X,
            sample_weight,
            x_squared_norms,
            best_centers,
            precompute_distances=precompute_distances,
            distances=distances,
        )

    return best_labels, best_inertia, best_centers, i + 1
762228
def validateMembership(likelihood, p, mc_source_id=1):
    

    cut_mc_source_id = (likelihood.catalog.mc_source_id == mc_source_id)

    # Spatial

    projector = ugali.utils.projector.Projector(likelihood.kernel.lon, likelihood.kernel.lat)
    x, y = projector.sphereToImage(likelihood.catalog.lon, likelihood.catalog.lat)
    
    pylab.figure()
    pylab.scatter(x[cut_mc_source_id], y[cut_mc_source_id], c='gray', s=100, edgecolors='none')
    pylab.scatter(x, y, c=p, edgecolors='none')
    pylab.colorbar()
    pylab.title('Membership Probability')
    pylab.xlabel(r'$\delta$ Lon (deg)')
    pylab.ylabel(r'$\delta$ Lat (deg)')
    
    # Spectral
    
    pylab.figure()
    pylab.scatter(likelihood.catalog.color[cut_mc_source_id], likelihood.catalog.mag[cut_mc_source_id], c='gray', s=100, edgecolors='none')
    pylab.scatter(likelihood.catalog.color, likelihood.catalog.mag, c=p, edgecolors='none')
    pylab.colorbar()
    pylab.title('Membership Probability')
    pylab.xlim(likelihood.roi.bins_color[0], likelihood.roi.bins_color[-1])
    pylab.ylim(likelihood.roi.bins_mag[-1], likelihood.roi.bins_mag[0])
    pylab.xlabel('Color (mag)')
    pylab.ylabel('Magnitude (mag)')
    
    # Membership accuracy

    p_array = []
    frac_members_array = []
    frac_members_array_low = []
    frac_members_array_high = []

    probability_bins = numpy.arange(0., 1. + 1.e-10, 0.1)
    for ii in range(0, len(probability_bins) - 1):

        if not numpy.any(numpy.logical_and(p >= probability_bins[ii], p < probability_bins[ii + 1])):
            continue
        
        #p_array.append(0.5 * (probability_bins[ii] + probability_bins[ii + 1]))
        cut = numpy.logical_and(p >= probability_bins[ii], p < probability_bins[ii + 1])
        p_array.append(numpy.median(p[cut]))

        n_members = numpy.sum(numpy.all([likelihood.catalog.mc_source_id == mc_source_id,
                                         p >= probability_bins[ii],
                                         p < probability_bins[ii + 1]], axis=0))
        n_not_members = numpy.sum(numpy.all([likelihood.catalog.mc_source_id != mc_source_id,
                                             p >= probability_bins[ii],
                                             p < probability_bins[ii + 1]], axis=0))

        frac_members, (frac_members_low, frac_members_high) = ugali.utils.bayesian_efficiency.confidenceInterval(n_members + n_not_members,
                                                                                                                 n_members,
                                                                                                                 0.68)
        frac_members_array.append(frac_members)
        frac_members_array_low.append(frac_members - frac_members_low)
        frac_members_array_high.append(frac_members_high - frac_members)
        
    pylab.figure()
    pylab.plot([0,1], [0,1], c='red')
    pylab.errorbar(p_array, frac_members_array, yerr=[frac_members_array_low, frac_members_array_high],
                   marker='o', color='blue', linestyle='none')
    pylab.xlim(0, 1)
    pylab.ylim(0, 1)
    pylab.xlabel('Membership Probability')
    pylab.ylabel('Fraction of True Satellite Members')

    # Where does richness come from?
    
    p_array = []
    p_sum_array = []
    n_members_array = []

    for ii in range(0, len(probability_bins) - 1):
        p_array.append(0.5 * (probability_bins[ii] + probability_bins[ii + 1]))
        p_sum_array.append(numpy.sum(p[numpy.logical_and(p >= probability_bins[ii],
                                                         p < probability_bins[ii + 1])]))
        n_members_array.append(numpy.sum(numpy.all([likelihood.catalog.mc_source_id == mc_source_id,
                                                    p >= probability_bins[ii],
                                                    p < probability_bins[ii + 1]], axis=0)))
        
    pylab.figure()
    pylab.scatter(p_array, p_sum_array, marker='o', c='blue')
    pylab.scatter(p_array, n_members_array, marker='o', c='red')
    pylab.xlim(0, 1)
    #pylab.ylim(0, 1)
    pylab.xlabel('Membership Probability')
    pylab.ylabel('Membership Probability Sum')

    # Purity and completeness

    x = numpy.linspace(0., 1., 1001) 
    purity = []
    completeness = []

    purity_reconstructed = []
    completeness_reconstructed = []

    for ii in range(0, len(x)):
        cut = p > (1 - x[ii])

        if numpy.sum(cut) < 1:
            purity.append(1.)
            completeness.append(0.)
            purity_reconstructed.append(1.)
            completeness_reconstructed.append(0.)
            continue
        
        purity_cut = numpy.logical_and(cut, likelihood.catalog.mc_source_id == mc_source_id)
        completeness_cut = likelihood.catalog.mc_source_id == mc_source_id
        
        purity.append(float(numpy.sum(purity_cut)) / numpy.sum(cut))
        completeness.append(float(numpy.sum(purity_cut)) / numpy.sum(completeness_cut))

        purity_reconstructed.append(numpy.mean(p[cut]))
        completeness_reconstructed.append(numpy.sum(p[cut]) / numpy.sum(p))
        
    pylab.figure()
    pylab.plot(x, purity, c='red', label='Purity')
    pylab.plot(x, completeness, c='blue', label='Completeness')
    pylab.plot(x, purity_reconstructed, c='red', linestyle='--', label='Purity Reconstructed')
    pylab.plot(x, completeness_reconstructed, c='blue', linestyle='--', label='Completeness Reconstructed')
    pylab.xlim(0, 1)
    pylab.ylim(0, 1)
    pylab.xlabel('1 - Membership Probability')
    #pylab.ylabel('Purity or Completeness')
    pylab.legend(loc='lower center')
508264
def cv_precompute(self, mask, b):
        '''
        

        '''

        # Get current chunk and mask outliers
        m1 = self.get_masked_chunk(b)
        flux = self.fraw[m1]
        K = GetCovariance(self.kernel, self.kernel_params,
                          self.time[m1], self.fraw_err[m1])
        med = np.nanmedian(flux)

        # Now mask the validation set
        M = lambda x, axis = 0: np.delete(x, mask, axis=axis)
        m2 = M(m1)
        mK = M(M(K, axis=0), axis=1)
        f = M(flux) - med

        # Pre-compute the matrices
        A = [None for i in range(self.pld_order)]
        B = [None for i in range(self.pld_order)]
        for n in range(self.pld_order):
            # Only compute up to the current PLD order
            if self.lam_idx >= n:
                X2 = self.X(n, m2)
                X1 = self.X(n, m1)
                A[n] = np.dot(X2, X2.T)
                B[n] = np.dot(X1, X2.T)
                del X1, X2

        if self.transit_model is None:
            C = 0
        else:
            C = np.zeros((len(m2), len(m2)))
            mean_transit_model = med * \
                np.sum([tm.depth * tm(self.time[m2])
                        for tm in self.transit_model], axis=0)
            f -= mean_transit_model
            for tm in self.transit_model:
                X2 = tm(self.time[m2]).reshape(-1, 1)
                C += tm.var_depth * np.dot(X2, X2.T)
                del X2

        return A, B, C, mK, f, m1, m2
188936
def fitlin(imgarr,refarr):
    
    # Initialize variables
    _mat = np.zeros((3,3),dtype=np.float64)
    _xorg = imgarr[0][0]
    _yorg = imgarr[0][1]
    _xoorg = refarr[0][0]
    _yoorg = refarr[0][1]
    _sigxox = 0.
    _sigxoy = 0.
    _sigxo = 0.
    _sigyox = 0.
    _sigyoy = 0.
    _sigyo = 0.

    _npos = len(imgarr)
    # Populate matrices
    for i in range(_npos):
        _mat[0][0] += np.power((imgarr[i][0] - _xorg),2)
        _mat[0][1] += (imgarr[i][0] - _xorg) * (imgarr[i][1] - _yorg)
        _mat[0][2] += (imgarr[i][0] - _xorg)
        _mat[1][1] += np.power((imgarr[i][1] - _yorg),2)
        _mat[1][2] += imgarr[i][1] - _yorg

        _sigxox += (refarr[i][0] - _xoorg)*(imgarr[i][0] - _xorg)
        _sigxoy += (refarr[i][0] - _xoorg)*(imgarr[i][1] - _yorg)
        _sigxo += refarr[i][0] - _xoorg
        _sigyox += (refarr[i][1] - _yoorg)*(imgarr[i][0] -_xorg)
        _sigyoy += (refarr[i][1] - _yoorg)*(imgarr[i][1] - _yorg)
        _sigyo += refarr[i][1] - _yoorg

    _mat[2][2] = _npos
    _mat[1][0] = _mat[0][1]
    _mat[2][0] = _mat[0][2]
    _mat[2][1] = _mat[1][2]

    # Now invert this matrix
    _mat = linalg.inv(_mat)

    _a  = _sigxox*_mat[0][0]+_sigxoy*_mat[0][1]+_sigxo*_mat[0][2]
    _b  = -1*(_sigxox*_mat[1][0]+_sigxoy*_mat[1][1]+_sigxo*_mat[1][2])
    #_x0 = _sigxox*_mat[2][0]+_sigxoy*_mat[2][1]+_sigxo*_mat[2][2]

    _c  = _sigyox*_mat[1][0]+_sigyoy*_mat[1][1]+_sigyo*_mat[1][2]
    _d  = _sigyox*_mat[0][0]+_sigyoy*_mat[0][1]+_sigyo*_mat[0][2]
    #_y0 = _sigyox*_mat[2][0]+_sigyoy*_mat[2][1]+_sigyo*_mat[2][2]

    _xt = _xoorg - _a*_xorg+_b*_yorg
    _yt = _yoorg - _d*_xorg-_c*_yorg

    return [_a,_b,_xt],[_c,_d,_yt]
477447
def _create_wcs (fitsheader):
    """"""
    wcsmodule = _load_wcs_module ()
    is_pywcs = hasattr (wcsmodule, 'UnitConverter')

    wcs = wcsmodule.WCS (fitsheader)
    wcs.wcs.set ()
    wcs.wcs.fix () # I'm interested in MJD computation via datfix()

    if hasattr (wcs, 'wcs_pix2sky'):
        wcs.wcs_pix2world = wcs.wcs_pix2sky
        wcs.wcs_world2pix = wcs.wcs_sky2pix

    return wcs
703558
def smpar(C):
    """"""
    m2 = C['m2'].real
    Lambda = C['Lambda'].real
    v = (sqrt(2 * m2 / Lambda) + 3 * m2**(3 / 2) /
         (sqrt(2) * Lambda**(5 / 2)) * C['phi'])
    GF = 1 / (sqrt(2) * v**2)  # TODO
    Mh2 = 2 * m2 * (1 - m2 / Lambda * (3 * C['phi'] - 4 * Lambda * C['phiBox'] +
                                       Lambda * C['phiD']))
    eps = C['phiWB'] * (v**2)
    gb = (C['g'] / (1 - C['phiW'] * (v**2))).real
    gpb = (C['gp'] / (1 - C['phiB'] * (v**2))).real
    gsb = (C['gs'] / (1 - C['phiG'] * (v**2))).real
    MW = gb * v / 2
    ZG0 = 1 + C['phiD'] * (v**2) / 4
    MZ = (sqrt(gb**2 + gpb**2) / 2 * v
          * (1 + eps * gb * gpb / (gb**2 + gpb**2)) * ZG0)
    Mnup = -(v**2) * C['llphiphi']
    Mep = v / sqrt(2) * (C['Ge'] - C['ephi'] * (v**2) / 2)
    Mup = v / sqrt(2) * (C['Gu'] - C['uphi'] * (v**2) / 2)
    Mdp = v / sqrt(2) * (C['Gd'] - C['dphi'] * (v**2) / 2)
    UeL, Me, UeR = ckmutil.diag.msvd(Mep)
    UuL, Mu, UuR = ckmutil.diag.msvd(Mup)
    UdL, Md, UdR = ckmutil.diag.msvd(Mdp)
    UnuL, Mnu = ckmutil.diag.mtakfac(Mnup)
    eb = (gb * gpb / sqrt(gb**2 + gpb**2) *
          (1 - eps * gb * gpb / (gb**2 + gpb**2)))
    K = UuL.conj().T @ UdL
    # U = UeL.conj().T @ UnuL
    sm = {}
    sm['GF'] = GF
    sm['alpha_e'] = eb**2 / (4 * pi)
    sm['alpha_s'] = gsb**2 / (4 * pi)
    sm['Vub'] = abs(K[0, 2])
    sm['Vcb'] = abs(K[1, 2])
    sm['Vus'] = abs(K[0, 1])
    sm['delta'] = phase(-K[0, 0] * K[0, 2].conj()
                        / (K[1, 0] * K[1, 2].conj()))
    # sm['U'] = Uu
    sm['m_W'] = MW
    sm['m_Z'] = MZ
    sm['m_h'] = sqrt(abs(Mh2))
    sm['m_u'] = Mu[0]
    sm['m_c'] = Mu[1]
    sm['m_t'] = Mu[2]
    sm['m_d'] = Md[0]
    sm['m_s'] = Md[1]
    sm['m_b'] = Md[2]
    sm['m_e'] = Me[0]
    sm['m_mu'] = Me[1]
    sm['m_tau'] = Me[2]
    return  {k: v.real for k, v in sm.items()}
529293
def cross_validate(self, ax):
        '''
        

        '''

        # The CDPP to beat
        cdpp_opt = self.get_cdpp_arr()

        # Loop over all chunks
        for b, brkpt in enumerate(self.breakpoints):

            log.info("Cross-validating chunk %d/%d..." %
                     (b + 1, len(self.breakpoints)))

            # Mask for current chunk
            m = self.get_masked_chunk(b)

            # Mask transits and outliers
            time = self.time[m]
            flux = self.fraw[m]
            ferr = self.fraw_err[m]
            med = np.nanmedian(self.fraw)

            # Setup the GP
            gp = GP(self.kernel, self.kernel_params, white=False)
            gp.compute(time, ferr)

            # The masks
            masks = list(Chunks(np.arange(0, len(time)),
                                len(time) // self.cdivs))

            # The pre-computed matrices
            pre_v = [self.cv_precompute(mask, b) for mask in masks]

            # Initialize with the nPLD solution
            log_lam_opt = np.log10(self.lam[b])
            scatter_opt = self.validation_scatter(
                log_lam_opt, b, masks, pre_v, gp, flux, time, med)
            log.info("Iter 0/%d: " % (self.piter) +
                     "logL = (%s), s = %.3f" %
                     (", ".join(["%.3f" % l for l in log_lam_opt]),
                      scatter_opt))

            # Do `piter` iterations
            for p in range(self.piter):

                # Perturb the initial condition a bit
                log_lam = np.array(
                    np.log10(self.lam[b])) * \
                    (1 + self.ppert * np.random.randn(len(self.lam[b])))
                scatter = self.validation_scatter(
                    log_lam, b, masks, pre_v, gp, flux, time, med)
                log.info("Initializing at: " +
                         "logL = (%s), s = %.3f" %
                         (", ".join(["%.3f" % l for l in log_lam]), scatter))

                # Call the minimizer
                log_lam, scatter, _, _, _, _ = \
                    fmin_powell(self.validation_scatter, log_lam,
                                args=(b, masks, pre_v, gp, flux, time, med),
                                maxfun=self.pmaxf, disp=False,
                                full_output=True)

                # Did it improve the CDPP?
                tmp = np.array(self.lam[b])
                self.lam[b] = 10 ** log_lam
                self.compute()
                cdpp = self.get_cdpp_arr()[b]
                self.lam[b] = tmp
                if cdpp < cdpp_opt[b]:
                    cdpp_opt[b] = cdpp
                    log_lam_opt = log_lam

                # Log it
                log.info("Iter %d/%d: " % (p + 1, self.piter) +
                         "logL = (%s), s = %.3f" %
                         (", ".join(["%.3f" % l for l in log_lam]), scatter))

            # The best solution
            log.info("Found minimum: logL = (%s), s = %.3f" %
                     (", ".join(["%.3f" % l for l in log_lam_opt]),
                      scatter_opt))
            self.lam[b] = 10 ** log_lam_opt

        # We're just going to plot lambda as a function of chunk number
        bs = np.arange(len(self.breakpoints))
        color = ['k', 'b', 'r', 'g', 'y']
        for n in range(self.pld_order):
            ax[0].plot(bs + 1, [np.log10(self.lam[b][n])
                                for b in bs], '.', color=color[n])
            ax[0].plot(bs + 1, [np.log10(self.lam[b][n])
                                for b in bs], '-', color=color[n], alpha=0.25)
        ax[0].set_ylabel(r'$\log\Lambda$', fontsize=5)
        ax[0].margins(0.1, 0.1)
        ax[0].set_xticks(np.arange(1, len(self.breakpoints) + 1))
        ax[0].set_xticklabels([])

        # Now plot the CDPP
        cdpp_arr = self.get_cdpp_arr()
        ax[1].plot(bs + 1, cdpp_arr, 'b.')
        ax[1].plot(bs + 1, cdpp_arr, 'b-', alpha=0.25)
        ax[1].margins(0.1, 0.1)
        ax[1].set_ylabel(r'Scatter (ppm)', fontsize=5)
        ax[1].set_xlabel(r'Chunk', fontsize=5)
        ax[1].set_xticks(np.arange(1, len(self.breakpoints) + 1))
87366
def generate(self, trilegal_filename, ra=None, dec=None,
                 n=2e4, ichrone='mist', MAfn=None,
                 mags=None, maxrad=None, f_binary=0.4, **kwargs):
        
        n = int(n)

        #generate/load BG primary stars from TRILEGAL simulation
        bgpop = BGStarPopulation_TRILEGAL(trilegal_filename,
                                        ra=ra, dec=dec, mags=mags,
                                        maxrad=maxrad, **kwargs)

        # Make sure that
        # properties of stars are within allowable range for isochrone.
        # This is a bit hacky, admitted.
        mass = bgpop.stars['m_ini'].values
        age = bgpop.stars['logAge'].values
        feh = bgpop.stars['[M/H]'].values

        ichrone = get_ichrone(ichrone)

        pct = 0.05 #pct distance from "edges" of ichrone interpolation
        mass[mass < ichrone.minmass*(1+pct)] = ichrone.minmass*(1+pct)
        mass[mass > ichrone.maxmass*(1-pct)] = ichrone.maxmass*(1-pct)
        age[age < ichrone.minage*(1+pct)] = ichrone.minage*(1+pct)
        age[age > ichrone.maxage*(1-pct)] = ichrone.maxage*(1-pct)
        feh[feh < ichrone.minfeh+0.05] = ichrone.minfeh+0.05
        feh[feh > ichrone.maxfeh-0.05] = ichrone.maxfeh-0.05

        distance = bgpop.stars['distance'].values

        #Generate binary population to draw eclipses from
        pop = MultipleStarPopulation(mA=mass, age=age, feh=feh,
                                            f_triple=0, f_binary=1,
                                            distance=distance,
                                            ichrone=ichrone)

        all_stars = pop.stars.dropna(subset=['mass_A'])
        all_stars.reset_index(inplace=True)

        #generate eclipses
        stars = pd.DataFrame()
        df_orbpop = pd.DataFrame()
        tot_prob = None; tot_dprob=None; prob_norm=None

        n_adapt = n
        while len(stars) < n:
            n_adapt = int(n_adapt)
            inds = np.random.randint(len(all_stars), size=n_adapt)

            s = all_stars.iloc[inds]

            #calculate limb-darkening coefficients
            u1A, u2A = ldcoeffs(s['Teff_A'], s['logg_A'])
            u1B, u2B = ldcoeffs(s['Teff_B'], s['logg_B'])

            inds, df, (prob,dprob) = calculate_eclipses(s['mass_A'], s['mass_B'],
                                                        s['radius_A'], s['radius_B'],
                                                        s['{}_mag_A'.format(self.band)],
                                                        s['{}_mag_B'.format(self.band)],
                                                        u11s=u1A, u21s=u2A,
                                                        u12s=u1B, u22s=u2B,
                                                        band=self.band,
                                                        period=self.period,
                                                        calc_mininc=True,
                                                        return_indices=True,
                                                        MAfn=MAfn)
            s = s.iloc[inds].copy()
            s.reset_index(inplace=True)
            for col in df.columns:
                s[col] = df[col]
            stars = pd.concat((stars, s))

            #new_df_orbpop = pop.orbpop.orbpop_long.dataframe.iloc[inds].copy()
            #new_df_orbpop.reset_index(inplace=True)

            #df_orbpop = pd.concat((df_orbpop, new_df_orbpop))

            logging.info('{} BEB systems generated (target {})'.format(len(stars),n))
            #logging.debug('{} nans in stars[dpri]'.format(np.isnan(stars['dpri']).sum()))
            #logging.debug('{} nans in df[dpri]'.format(np.isnan(df['dpri']).sum()))

            if tot_prob is None:
                prob_norm = (1/dprob**2)
                tot_prob = prob
                tot_dprob = dprob
            else:
                prob_norm = (1/tot_dprob**2 + 1/dprob**2)
                tot_prob = (tot_prob/tot_dprob**2 + prob/dprob**2)/prob_norm
                tot_dprob = 1/np.sqrt(prob_norm)

            n_adapt = min(int(1.2*(n-len(stars)) * n_adapt//len(s)), 5e5)
            #logging.debug('n_adapt = {}'.format(n_adapt))
            n_adapt = max(n_adapt, 100)
            n_adapt = int(n_adapt)

        stars = stars.iloc[:n]

        if 'level_0' in stars:
            stars.drop('level_0', axis=1, inplace=True) #dunno where this came from
        stars = stars.reset_index()
        stars.drop('index', axis=1, inplace=True)

        stars['mass_1'] = stars['mass_A']
        stars['radius_1'] = stars['radius_A']
        stars['mass_2'] = stars['mass_B']
        stars['radius_2'] = stars['radius_B']

        MultipleStarPopulation.__init__(self, stars=stars,
                                        #orbpop=orbpop,
                                        f_triple=0, f_binary=f_binary,
                                        period_long=self.period)

        priorfactors = {'f_binary':f_binary}

        #attributes needed for BGStarPopulation
        self.density = bgpop.density
        self.trilegal_args = bgpop.trilegal_args
        self._maxrad = bgpop._maxrad

        #create an OrbitPopulation here?

        EclipsePopulation.__init__(self, stars=stars, #orbpop=orbpop,
                                   period=self.period, cadence=self.cadence,
                                   model=self.model,
                                   lhoodcachefile=self.lhoodcachefile,
                                   priorfactors=priorfactors, prob=tot_prob)

        #add Rsky property
        self.stars['Rsky'] = randpos_in_circle(len(self.stars),
                                               self._maxrad, return_rad=True)
1184815
def imrphenomc_tmplt(**kwds):
    
    # Pull out the input arguments
    f_min = float128(kwds['f_lower'])
    f_max = float128(kwds['f_final'])
    delta_f = float128(kwds['delta_f'])
    distance = float128(kwds['distance'])
    mass1 = float128(kwds['mass1'])
    mass2 = float128(kwds['mass2'])
    spin1z = float128(kwds['spin1z'])
    spin2z = float128(kwds['spin2z'])

    if 'out' in kwds:
        out = kwds['out']
    else:
        out = None

    # Calculate binary parameters
    M = mass1 + mass2
    eta = mass1 * mass2 / (M * M)
    Xi = (mass1 * spin1z / M) + (mass2 * spin2z / M)
    Xisum = 2.*Xi
    Xiprod = Xi*Xi
    Xi2 = Xi*Xi

    m_sec = M * lal.MTSUN_SI;
    piM = lal.PI * m_sec;

    ## The units of distance given as input is taken to pe Mpc. Converting to SI
    distance *= (1.0e6 * lal.PC_SI / (2. * sqrt(5. / (64.*lal.PI)) * M * lal.MRSUN_SI * M * lal.MTSUN_SI))

    # Check if the value of f_max is correctly given, else replace with the fCut
    # used in the PhenomB code in lalsimulation. The various coefficients come
    # from Eq.(4.18) of http://arxiv.org/pdf/0710.2335 and
    # Table I of http://arxiv.org/pdf/0712.0343
    if not f_max:
        f_max = (1.7086 * eta * eta - 0.26592 * eta + 0.28236) / piM

    # Transform the eta, chi to Lambda parameters, using Eq 5.14, Table II of Main
    # paper.
    z101 = -2.417e-03
    z102 = -1.093e-03
    z111 = -1.917e-02
    z110 = 7.267e-02
    z120 = -2.504e-01

    z201 = 5.962e-01
    z202 = -5.600e-02
    z211 = 1.520e-01
    z210 = -2.970e+00
    z220 = 1.312e+01

    z301 = -3.283e+01
    z302 = 8.859e+00
    z311 = 2.931e+01
    z310 = 7.954e+01
    z320 = -4.349e+02

    z401 = 1.619e+02
    z402 = -4.702e+01
    z411 = -1.751e+02
    z410 = -3.225e+02
    z420 = 1.587e+03

    z501 = -6.320e+02
    z502 = 2.463e+02
    z511 = 1.048e+03
    z510 = 3.355e+02
    z520 = -5.115e+03

    z601 = -4.809e+01
    z602 = -3.643e+02
    z611 = -5.215e+02
    z610 = 1.870e+03
    z620 = 7.354e+02

    z701 = 4.149e+00
    z702 = -4.070e+00
    z711 = -8.752e+01
    z710 = -4.897e+01
    z720 = 6.665e+02

    z801 = -5.472e-02
    z802 = 2.094e-02
    z811 = 3.554e-01
    z810 = 1.151e-01
    z820 = 9.640e-01

    z901 = -1.235e+00
    z902 = 3.423e-01
    z911 = 6.062e+00
    z910 = 5.949e+00
    z920 = -1.069e+01

    eta2 = eta*eta
    Xi2 = Xiprod

    # Calculate alphas, gamma, deltas from Table II and Eq 5.14 of Main paper
    a1 = z101 * Xi + z102 * Xi2 + z111 * eta * Xi + z110 * eta + z120 * eta2
    a2 = z201 * Xi + z202 * Xi2 + z211 * eta * Xi + z210 * eta + z220 * eta2
    a3 = z301 * Xi + z302 * Xi2 + z311 * eta * Xi + z310 * eta + z320 * eta2
    a4 = z401 * Xi + z402 * Xi2 + z411 * eta * Xi + z410 * eta + z420 * eta2
    a5 = z501 * Xi + z502 * Xi2 + z511 * eta * Xi + z510 * eta + z520 * eta2
    a6 = z601 * Xi + z602 * Xi2 + z611 * eta * Xi + z610 * eta + z620 * eta2

    g1 = z701 * Xi + z702 * Xi2 + z711 * eta * Xi + z710 * eta + z720 * eta2

    del1 = z801 * Xi + z802 * Xi2 + z811 * eta * Xi + z810 * eta + z820 * eta2
    del2 = z901 * Xi + z902 * Xi2 + z911 * eta * Xi + z910 * eta + z920 * eta2

    # Get the spin of the final BH
    afin = FinalSpin( Xi, eta )
    Q = Qa( abs(afin) )

    # Get the fRD
    frd = fRD( abs(afin), M)
    Mfrd = frd * m_sec

    # Define the frequencies where SPA->PM->RD
    f1 = 0.1 * frd
    Mf1 = m_sec * f1
    f2 = frd
    Mf2 = m_sec * f2
    d1 = 0.005
    d2 = 0.005
    f0 = 0.98 * frd
    Mf0 = m_sec * f0
    d0 = 0.015

    # Now use this frequency for calculation of betas
    # calculate beta1 and beta2, that appear in Eq 5.7 in the main paper.
    b2 = ((-5./3.)* a1 * pow(Mfrd,(-8./3.)) - a2/(Mfrd*Mfrd) - \
      (a3/3.)*pow(Mfrd,(-4./3.)) + (2./3.)* a5 * pow(Mfrd,(-1./3.)) + a6)/eta

    psiPMrd = (a1 * pow(Mfrd,(-5./3.)) + a2/Mfrd + a3 * pow(Mfrd,(-1./3.)) + \
      a4 + a5 * pow(Mfrd,(2./3.)) + a6 * Mfrd)/eta
    b1 = psiPMrd - (b2 * Mfrd)

    ### Calculate the PN coefficients, Eq A3 - A5 of main paper ###
    pfaN = 3.0/(128.0 * eta)
    pfa2 = (3715./756.) + (55.*eta/9.0)
    pfa3 = -16.0*lal.PI + (113./3.)*Xi - 38.*eta*Xisum/3.
    pfa4 = (152.93365/5.08032) - 50.*Xi2 + eta*(271.45/5.04 + 1.25*Xiprod) + \
        3085.*eta2/72.
    pfa5 = lal.PI*(386.45/7.56 - 65.*eta/9.) - \
        Xi*(735.505/2.268 + 130.*eta/9.) + Xisum*(1285.0*eta/8.1 + 170.*eta2/9.) - \
        10.*Xi2*Xi/3. + 10.*eta*Xi*Xiprod
    pfa6 = 11583.231236531/4.694215680 - 640.0*lal.PI*lal.PI/3. - \
        6848.0*lal.GAMMA/21. - 684.8*log(64.)/6.3 + \
        eta*(2255.*lal.PI*lal.PI/12. - 15737.765635/3.048192) + \
        76.055*eta2/1.728 - (127.825*eta2*eta/1.296) + \
        2920.*lal.PI*Xi/3. - (175. - 1490.*eta)*Xi2/3. - \
        (1120.*lal.PI/3. - 1085.*Xi/3.)*eta*Xisum + \
        (269.45*eta/3.36 - 2365.*eta2/6.)*Xiprod

    pfa6log = -6848./63.

    pfa7 = lal.PI*(770.96675/2.54016 + 378.515*eta/1.512 - 740.45*eta2/7.56) - \
        Xi*(20373.952415/3.048192 + 1509.35*eta/2.24 - 5786.95*eta2/4.32) + \
        Xisum*(4862.041225*eta/1.524096 + 1189.775*eta2/1.008 - 717.05*eta2*eta/2.16 - 830.*eta*Xi2/3. + 35.*eta2*Xiprod/3.) - \
        560.*lal.PI*Xi2 + 20.*lal.PI*eta*Xiprod + \
        Xi2*Xi*(945.55/1.68 - 85.*eta) + Xi*Xiprod*(396.65*eta/1.68 + 255.*eta2)


    xdotaN = 64.*eta/5.
    xdota2 = -7.43/3.36 - 11.*eta/4.
    xdota3 = 4.*lal.PI - 11.3*Xi/1.2 + 19.*eta*Xisum/6.
    xdota4 = 3.4103/1.8144 + 5*Xi2 + eta*(13.661/2.016 - Xiprod/8.) + 5.9*eta2/1.8
    xdota5 = -lal.PI*(41.59/6.72 + 189.*eta/8.) - Xi*(31.571/1.008 - 116.5*eta/2.4) + \
          Xisum*(21.863*eta/1.008 - 79.*eta2/6.) - 3*Xi*Xi2/4. + \
          9.*eta*Xi*Xiprod/4.
    xdota6 = 164.47322263/1.39708800 - 17.12*lal.GAMMA/1.05 + \
          16.*lal.PI*lal.PI/3 - 8.56*log(16.)/1.05 + \
          eta*(45.1*lal.PI*lal.PI/4.8 - 561.98689/2.17728) + \
          5.41*eta2/8.96 - 5.605*eta*eta2/2.592 - 80.*lal.PI*Xi/3. + \
          eta*Xisum*(20.*lal.PI/3. - 113.5*Xi/3.6) + \
          Xi2*(64.153/1.008 - 45.7*eta/3.6) - \
          Xiprod*(7.87*eta/1.44 - 30.37*eta2/1.44)

    xdota6log = -856./105.

    xdota7 = -lal.PI*(4.415/4.032 - 358.675*eta/6.048 - 91.495*eta2/1.512) - \
          Xi*(252.9407/2.7216 - 845.827*eta/6.048 + 415.51*eta2/8.64) + \
          Xisum*(158.0239*eta/5.4432 - 451.597*eta2/6.048 + 20.45*eta2*eta/4.32 + 107.*eta*Xi2/6. - 5.*eta2*Xiprod/24.) + \
          12.*lal.PI*Xi2 - Xi2*Xi*(150.5/2.4 + eta/8.) + \
          Xi*Xiprod*(10.1*eta/2.4 + 3.*eta2/8.)


    AN = 8.*eta*sqrt(lal.PI/5.)
    A2 = (-107. + 55.*eta)/42.
    A3 = 2.*lal.PI - 4.*Xi/3. + 2.*eta*Xisum/3.
    A4 = -2.173/1.512 - eta*(10.69/2.16 - 2.*Xiprod) + 2.047*eta2/1.512
    A5 = -10.7*lal.PI/2.1 + eta*(3.4*lal.PI/2.1)

    A5imag = -24.*eta

    A6 = 270.27409/6.46800 - 8.56*lal.GAMMA/1.05 + \
      2.*lal.PI*lal.PI/3. + \
      eta*(4.1*lal.PI*lal.PI/9.6 - 27.8185/3.3264) - \
      20.261*eta2/2.772 + 11.4635*eta*eta2/9.9792 - \
      4.28*log(16.)/1.05

    A6log = -428./105.

    A6imag = 4.28*lal.PI/1.05

    ### Define other parameters needed by waveform generation ###
    kmin = int(f_min / delta_f)
    kmax = int(f_max / delta_f)
    n = kmax + 1;

    if not out:
        htilde = FrequencySeries(zeros(n,dtype=numpy.complex128), delta_f=delta_f, copy=False)
    else:
        if type(out) is not Array:
            raise TypeError("Output must be an instance of Array")
        if len(out) < kmax:
            raise TypeError("Output array is too small")
        if out.dtype != complex64:
            raise TypeError("Output array is the wrong dtype")
        htilde = FrequencySeries(out, delta_f=delta_f, copy=False)

    phenomC_kernel(htilde.data[kmin:kmax], kmin, delta_f, eta, Xi, distance,
                                       m_sec,  piM,  Mfrd,
                                       pfaN,  pfa2,  pfa3,  pfa4, pfa5,  pfa6,  pfa6log,  pfa7,
                                       a1,  a2,  a3,  a4, a5,  a6,  b1,  b2,
                                       Mf1,  Mf2,  Mf0, d1,  d2,  d0,
                                       xdota2,  xdota3,  xdota4, xdota5,  xdota6,  xdota6log,
                                       xdota7,  xdotaN,  AN, A2,  A3,  A4,  A5,
                                       A5imag,  A6,  A6log,  A6imag,
                                       g1,  del1,  del2,  Q )
    hp = htilde
    hc = htilde * 1j
    return hp, hc
738816
def lensNMF(data, k, ks=1):
    
    # TODO: why is this not working
    n_rounds = k/ks
    R_i = data.copy()
    nmf = NMF(ks)
    nmf2 = NMF(ks, init='custom')
    w_is = []
    h_is = []
    rs = []
    w_i = np.zeros((data.shape[0], ks))
    h_i = np.zeros((ks, data.shape[1]))
    for i in range(n_rounds):
        R_i = R_i - w_i.dot(h_i)
        R_i[R_i < 0] = 0
        """
        P_r = R_i.sum(1)/R_i.sum()
        print(P_r.shape)
        P_c = R_i.sum(0)/R_i.sum()
        print(P_c.shape)
        row_choice = np.random.choice(range(len(P_r)), p=P_r)
        print(row_choice)
        col_choice = np.random.choice(range(len(P_c)), p=P_c)
        print(col_choice)
        D_r = cosine_similarity(data[row_choice:row_choice+1,:], data)
        D_c = cosine_similarity(data[:,col_choice:col_choice+1].T, data.T)
        D_r = np.diag(D_r.flatten())
        D_c = np.diag(D_c.flatten())
        R_L = D_r.dot(R_i).dot(D_c)
        w_i = nmf.fit_transform(R_L)
        """
        w_i = nmf.fit_transform(R_i)
        h_i = nmf.components_
        #nmf2.fit_transform(R_i, W=w_i, H=nmf.components_)
        #h_i = nmf2.components_
        #h_i[h_i < 0] = 0
        w_is.append(w_i)
        h_is.append(h_i)
        rs.append(R_i)
    return np.hstack(w_is), np.vstack(h_is), rs
285440
def do_objs(kbos):
    """"""

    import orbfit, ephem, math
    import re
    re_string=w.FilterVar.get()
    vlist=[]
    for name in kbos:
        if not re.search(re_string,name):
            continue
        vlist.append(name)
        if type(kbos[name])==type(ephem.EllipticalBody()):
	    kbos[name].compute(w.date.get())
	    ra=kbos[name].ra
	    dec=kbos[name].dec
	    a=math.radians(10.0/3600.0)
	    b=a
	    ang=0.0
	    color='blue'
            yoffset=+10
	    xoffset=+10
        else:
            yoffset=-10
            xoffset=-10
            file=kbos[name]
	    jdate=ephem.julian_date(w.date.get())
	    obs=568
 	    try:
	        position=orbfit.predict(file,jdate,obs)
	    except:
		continue
	    ra=math.radians(position[0])
	    dec=math.radians(position[1])
	    a=math.radians(position[2]/3600.0)
	    b=math.radians(position[3]/3600.0)
	    ang=math.radians(position[4])
	    if ( a> math.radians(1.0) ):
	        color='green'
            else:
	        color='black'
	if w.show_ellipse.get()==1 :
            if ( a < math.radians(5.0) ): 
	        w.create_ellipse(ra,dec,a,b,ang)
        if ( a < math.radians(1.0) ): 
            w.create_point(ra,dec,size=2,color=color)
	if w.show_labels.get()==1:
	    w.label(ra,dec,name,offset=[xoffset,yoffset])
    vlist.sort()
    for v in vlist:
        w.objList.insert(END,v)
    w.plot_pointings()
152796
function pm(m) {
    if (m[4] === 0 && m[5] === 0) {
      if (m[1] === 0 && m[2] === 0) {
        if (m[0] === 1 && m[3] === 1) {
          return '';
        }
        return 'scale(' + pf(m[0]) + ' ' + pf(m[3]) + ')';
      }
      if (m[0] === m[3] && m[1] === -m[2]) {
        var a = Math.acos(m[0]) * 180 / Math.PI;
        return 'rotate(' + pf(a) + ')';
      }
    } else {
      if (m[0] === 1 && m[1] === 0 && m[2] === 0 && m[3] === 1) {
        return 'translate(' + pf(m[4]) + ' ' + pf(m[5]) + ')';
      }
    }
    return 'matrix(' + pf(m[0]) + ' ' + pf(m[1]) + ' ' + pf(m[2]) + ' ' +
      pf(m[3]) + ' ' + pf(m[4]) + ' ' + pf(m[5]) + ')';
  }
1357769
def generate(self,rprs=None, mass=None, radius=None,
                n=2e4, fp_specific=0.01, u1=None, u2=None,
                 starmodel=None,
                Teff=None, logg=None, rbin_width=0.3,
                MAfn=None, lhoodcachefile=None):
        

        n = int(n)

        if starmodel is None:
            if type(mass) is type((1,)):
                mass = dists.Gaussian_Distribution(*mass)
            if isinstance(mass, dists.Distribution):
                mdist = mass
                mass = mdist.rvs(1e5)

            if type(radius) is type((1,)):
                radius = dists.Gaussian_Distribution(*radius)
            if isinstance(radius, dists.Distribution):
                rdist = radius
                radius = rdist.rvs(1e5)
        else:
            samples = starmodel.random_samples(1e5)
            mass = samples['mass_0_0'].values
            radius = samples['radius_0_0'].values
            Teff = samples['Teff_0_0'].mean()
            logg = samples['logg_0_0'].mean()

        logging.debug('star mass: {}'.format(mass))
        logging.debug('star radius: {}'.format(radius))
        logging.debug('Teff: {}'.format(Teff))
        logging.debug('logg: {}'.format(logg))

        if u1 is None or u2 is None:
            if Teff is None or logg is None:
                logging.warning('Teff, logg not provided; using solar limb darkening')
                u1 = 0.394; u2=0.296
            else:
                u1,u2 = ldcoeffs(Teff, logg)

        #use point estimate of rprs to construct planets in radius bin
        #rp = self.rprs*np.median(radius)
        #rbin_min = (1-rbin_width)*rp
        #rbin_max = (1+rbin_width)*rp

        rprs_bin_min = (1-rbin_width)*self.rprs
        rprs_bin_max = (1+rbin_width)*self.rprs

        radius_p = radius * (np.random.random(int(1e5))*(rprs_bin_max - rprs_bin_min) + rprs_bin_min)
        mass_p = (radius_p*RSUN/REARTH)**2.06 * MEARTH/MSUN #hokey, but doesn't matter

        logging.debug('planet radius: {}'.format(radius_p))

        stars = pd.DataFrame()
        #df_orbpop = pd.DataFrame() #for orbit population

        tot_prob = None; tot_dprob = None; prob_norm = None
        n_adapt = n
        while len(stars) < n:
            n_adapt = int(n_adapt)
            inds = np.random.randint(len(mass), size=n_adapt)

            #calculate eclipses.
            ecl_inds, df, (prob,dprob) = calculate_eclipses(mass[inds], mass_p[inds],
                                                        radius[inds], radius_p[inds],
                                                        15, np.inf, #arbitrary
                                                        u11s=u1, u21s=u2,
                                                        band=self.band,
                                                        period=self.period,
                                                        calc_mininc=True,
                                                        return_indices=True,
                                                        MAfn=MAfn)

            df['mass_A'] = mass[inds][ecl_inds]
            df['mass_B'] = mass_p[inds][ecl_inds]
            df['radius_A'] = radius[inds][ecl_inds]
            df['radius_B'] = radius_p[inds][ecl_inds]
            df['u1'] = u1 * np.ones_like(df['mass_A'])
            df['u2'] = u2 * np.ones_like(df['mass_A'])
            df['P'] = self.period * np.ones_like(df['mass_A'])

            ok = (df['dpri']>0) & (df['T14_pri'] > 0)

            stars = pd.concat((stars, df[ok]))

            logging.info('{} Transiting planet systems generated (target {})'.format(len(stars),n))
            logging.debug('{} nans in stars[dpri]'.format(np.isnan(stars['dpri']).sum()))

            if tot_prob is None:
                prob_norm = (1/dprob**2)
                tot_prob = prob
                tot_dprob = dprob
            else:
                prob_norm = (1/tot_dprob**2 + 1/dprob**2)
                tot_prob = (tot_prob/tot_dprob**2 + prob/dprob**2)/prob_norm
                tot_dprob = 1/np.sqrt(prob_norm)

            n_adapt = min(int(1.2*(n-len(stars)) * n_adapt//len(df)), 5e4)
            n_adapt = max(n_adapt, 100)

        stars = stars.reset_index()
        stars.drop('index', axis=1, inplace=True)
        stars = stars.iloc[:n]

        stars['mass_1'] = stars['mass_A']
        stars['radius_1'] = stars['radius_A']
        stars['mass_2'] = stars['mass_B']
        stars['radius_2'] = stars['radius_B']

        #make OrbitPopulation?

        #finish below.

        if fp_specific is None:
            rp = stars['radius_2'].mean() * RSUN/REARTH
            fp_specific = fp_fressin(rp)

        priorfactors = {'fp_specific':fp_specific}

        self._starmodel = starmodel

        EclipsePopulation.__init__(self, stars=stars,
                                   period=self.period, cadence=self.cadence,
                                   model=self.model,
                                   priorfactors=priorfactors, prob=tot_prob,
                                   lhoodcachefile=lhoodcachefile)