Fit an orbit to data in a given potential
344053
def _fit_orbit(orb,vxvv,vxvv_err,pot,radec=False,lb=False,
               customsky=False,lb_to_customsky=None,
               pmllpmbb_to_customsky=None,
               tintJ=100,ntintJ=1000,integrate_method='dopr54_c',
               ro=None,vo=None,obs=None,disp=False):
    """"""
    # Need to turn this off for speed
    coords._APY_COORDS_ORIG= coords._APY_COORDS
    coords._APY_COORDS= False
    #Import here, because otherwise there is an infinite loop of imports
    from galpy.actionAngle import actionAngleIsochroneApprox, actionAngle
    #Mock this up, bc we want to use its orbit-integration routines
    class mockActionAngleIsochroneApprox(actionAngleIsochroneApprox):
        def __init__(self,tintJ,ntintJ,pot,integrate_method='dopr54_c'):
            actionAngle.__init__(self)
            self._tintJ= tintJ
            self._ntintJ=ntintJ
            self._tsJ= nu.linspace(0.,self._tintJ,self._ntintJ)
            self._integrate_dt= None
            self._pot= pot
            self._integrate_method= integrate_method
            return None
    tmockAA= mockActionAngleIsochroneApprox(tintJ,ntintJ,pot,
                                            integrate_method=integrate_method)
    opt_vxvv= optimize.fmin_powell(_fit_orbit_mlogl,orb.vxvv,
                                   args=(vxvv,vxvv_err,pot,radec,lb,
                                         customsky,lb_to_customsky,
                                         pmllpmbb_to_customsky,
                                         tmockAA,
                                         ro,vo,obs),
                                   disp=disp)
    maxLogL= -_fit_orbit_mlogl(opt_vxvv,vxvv,vxvv_err,pot,radec,lb,
                               customsky,lb_to_customsky,pmllpmbb_to_customsky,
                               tmockAA,
                               ro,vo,obs)
    coords._APY_COORDS= coords._APY_COORDS_ORIG
    return (opt_vxvv,maxLogL)
65347
def get_least_orbits(atom_index, cell, site_symmetry, symprec=1e-5):
    """"""
    orbits = _get_orbits(atom_index, cell, site_symmetry, symprec)
    mapping = np.arange(cell.get_number_of_atoms())

    for i, orb in enumerate(orbits):
        for num in np.unique(orb):
            if mapping[num] > mapping[i]:
                mapping[num] = mapping[i]

    return np.unique(mapping)
375495
def _coulomb(n1, n2, k, r):
    """"""
    # Get relevant positional data
    delta = [x2 - x1 for x1, x2 in zip(n1['velocity'], n2['velocity'])]
    distance = sqrt(sum(d ** 2 for d in delta))

    # If the deltas are too small, use random values to keep things moving
    if distance < 0.1:
        delta = [uniform(0.1, 0.2) for _ in repeat(None, 3)]
        distance = sqrt(sum(d ** 2 for d in delta))

    # If the distance isn't huge (ie. Coulomb is negligible), calculate
    if distance < r:
        force = (k / distance) ** 2
        n1['force'] = [f - force * d for f, d in zip(n1['force'], delta)]
        n2['force'] = [f + force * d for f, d in zip(n2['force'], delta)]
195776
def _fit_strel_to_im_2d(im, strel, r, x, y):
    
    """
    elem = strel.copy()
    x_dim, y_dim = im.shape
    x_min = x-r
    x_max = x+r+1
    y_min = y-r
    y_max = y+r+1
    if x_min < 0:
        x_adj = -x_min
        elem = elem[x_adj:, :]
        x_min = 0
    elif x_max > x_dim:
        x_adj = x_max - x_dim
        elem = elem[:-x_adj, :]
    if y_min < 0:
        y_adj = -y_min
        elem = elem[:, y_adj:]
        y_min = 0
    elif y_max > y_dim:
        y_adj = y_max - y_dim
        elem = elem[:, :-y_adj]
    ex, ey = elem.shape
    im[x_min:x_min+ex, y_min:y_min+ey] += elem
    return im
144679
def check(self, feature):
        """"""
        mapper = feature.as_dataframe_mapper()
        mapper.fit_transform(self.X, y=self.y)
185038
def _integrate_direct_python(self,t,pot,**kwargs):
        """"""
        #Prepare input for direct_nbody
        q= []
        p= []
        nq= len(self.orbits)
        dim= self.orbits[0].dim()
        if pot is None:
            thispot= None
        elif dim == 2:
            thispot= RZToplanarPotential(pot)
        else:
            thispot= pot
        for ii in range(nq):
            #Transform to rectangular frame
            if dim == 1:
                thisq= nu.array([self.orbits[ii].x()]).flatten()
                thisp= nu.array([self.orbits[ii].vx()]).flatten()
            elif dim == 2:
                thisq= nu.array([self.orbits[ii].x(),
                                 self.orbits[ii].y()]).flatten()
                thisp= nu.array([self.orbits[ii].vx(),
                                 self.orbits[ii].vy()]).flatten()
            elif dim == 3:
                thisq= nu.array([self.orbits[ii].x(),
                                 self.orbits[ii].y(),
                                 self.orbits[ii].z()]).flatten()
                thisp= nu.array([self.orbits[ii].vx(),
                                 self.orbits[ii].vy(),
                                 self.orbits[ii].vz()]).flatten()
            q.append(thisq)
            p.append(thisp)
        #Run simulation
        nbody_out= direct_nbody(q,p,self.masses,t,pot=thispot,**kwargs)
        #Post-process output
        nt= len(nbody_out)
        out= []
        for ii in range(nt):
            snap_orbits= []
            for jj in range(nq):
                if dim == 3:
                    #go back to the cylindrical frame
                    R= nu.sqrt(nbody_out[ii][0][jj][0]**2.
                               +nbody_out[ii][0][jj][1]**2.)
                    phi= nu.arccos(nbody_out[ii][0][jj][0]/R)
                    if nbody_out[ii][0][jj][1] < 0.: phi= 2.*nu.pi-phi
                    vR= nbody_out[ii][1][jj][0]*nu.cos(phi)\
                        +nbody_out[ii][1][jj][1]*nu.sin(phi)
                    vT= nbody_out[ii][1][jj][1]*nu.cos(phi)\
                        -nbody_out[ii][1][jj][0]*nu.sin(phi)
                    vxvv= nu.zeros(dim*2)
                    vxvv[3]= nbody_out[ii][0][jj][2]
                    vxvv[4]= nbody_out[ii][1][jj][2]
                    vxvv[0]= R
                    vxvv[1]= vR
                    vxvv[2]= vT
                    vxvv[5]= phi
                if dim == 2:
                    #go back to the cylindrical frame
                    R= nu.sqrt(nbody_out[ii][0][jj][0]**2.
                               +nbody_out[ii][0][jj][1]**2.)
                    phi= nu.arccos(nbody_out[ii][0][jj][0]/R)
                    if nbody_out[ii][0][jj][1] < 0.: phi= 2.*nu.pi-phi
                    vR= nbody_out[ii][1][jj][0]*nu.cos(phi)\
                        +nbody_out[ii][1][jj][1]*nu.sin(phi)
                    vT= nbody_out[ii][1][jj][1]*nu.cos(phi)\
                        -nbody_out[ii][1][jj][0]*nu.sin(phi)
                    vxvv= nu.zeros(dim*2)
                    vxvv[0]= R
                    vxvv[1]= vR
                    vxvv[2]= vT
                    vxvv[3]= phi
                if dim == 1:
                    vxvv= [nbody_out[ii][0][jj],nbody_out[ii][1][jj]]
                snap_orbits.append(Orbit(vxvv))
            out.append(Snapshot(snap_orbits,self.masses))
        return out
197406
def fit(self, x0=None, distribution='lognormal', n=None, **kwargs):
        '''
        '''
        dist = {'lognormal': PSDLognormal, 
                'GGS': PSDGatesGaudinSchuhman, 
                'RR': PSDRosinRammler}[distribution]
        
        if distribution == 'lognormal':
            if x0 is None:
                d_characteristic = sum([fi*di for fi, di in zip(self.fractions, self.Dis)])
                s = 0.4
                x0 = [d_characteristic, s]
        elif distribution == 'GGS':
            if x0 is None:
                d_characteristic = sum([fi*di for fi, di in zip(self.fractions, self.Dis)])
                m = 1.5
                x0 = [d_characteristic, m]
        elif distribution == 'RR':
            if x0 is None:
                x0 = [5E-6, 1e-2]
        from scipy.optimize import minimize
        return minimize(self._fit_obj_function, x0, args=(dist, n), **kwargs)
823752
def fit(self, X, y=None):
        
        X = array2d(X)
        self.n_features = X.shape[1]
        self.n_bins = self.n_bins_per_feature ** self.n_features

        if self.min is None:
            min = np.min(X, axis=0)
        elif isinstance(self.min, numbers.Number):
            min = self.min * np.ones(self.n_features)
        else:
            min = np.asarray(self.min)
            if not min.shape == (self.n_features,):
                raise ValueError('min shape error')

        if self.max is None:
            max = np.max(X, axis=0)
        elif isinstance(self.max, numbers.Number):
            max = self.max * np.ones(self.n_features)
        else:
            max = np.asarray(self.max)
            if not max.shape == (self.n_features,):
                raise ValueError('max shape error')

        self.grid = np.array(
            [np.linspace(min[i] - EPS, max[i] + EPS, self.n_bins_per_feature + 1)
             for i in range(self.n_features)])

        return self
354300
def orbit_y(self, angle):
        '''

        '''
        
        # Subtract pivot point
        self.position -= self.pivot
        
        # Rotate
        rot = rotation_matrix(-angle, self.b)[:3,:3]
        self.position = np.dot(rot, self.position)
        
        # Add again the pivot point
        self.position += self.pivot
        
        self.a = np.dot(rot, self.a)
        self.b = np.dot(rot, self.b)
        self.c = np.dot(rot, self.c)
744517
def fit(self, X, y):
        """"""

        self._data = X
        self._classes = np.unique(y)
        self._labels = y
        self._is_fitted = True
1024317
def orbit(self, orbit):
        

        self._orbit = orbit
        tle = Tle.from_orbit(orbit)
        lines = tle.text.splitlines()

        if len(lines) == 3:
            _, line1, line2 = lines
        else:
            line1, line2 = lines

        self.tle = twoline2rv(line1, line2, wgs72)
742312
def proper(self,x):
        """"""
        x[x < 0] = self.max_fit
        x[np.isnan(x)] = self.max_fit
        x[np.isinf(x)] = self.max_fit
        return x
529721
def _keplerian_to_cartesian(cls, coord, center):
        

        a, e, i, Ω, ω, ν = coord

        p = a * (1 - e ** 2)
        r = p / (1 + e * cos(ν))
        h = sqrt(center.µ * p)
        x = r * (cos(Ω) * cos(ω + ν) - sin(Ω) * sin(ω + ν) * cos(i))
        y = r * (sin(Ω) * cos(ω + ν) + cos(Ω) * sin(ω + ν) * cos(i))
        z = r * sin(i) * sin(ω + ν)
        vx = x * h * e / (r * p) * sin(ν) - h / r * (cos(Ω) * sin(ω + ν) + sin(Ω) * cos(ω + ν) * cos(i))
        vy = y * h * e / (r * p) * sin(ν) - h / r * (sin(Ω) * sin(ω + ν) - cos(Ω) * cos(ω + ν) * cos(i))
        vz = z * h * e / (r * p) * sin(ν) + h / r * sin(i) * cos(ω + ν)

        return np.array([x, y, z, vx, vy, vz], dtype=float)
862692
def spherical_galaxy_orbit(
    orbit_x,
    orbit_y,
    orbit_z,
    N_stars=100,
    sigma_r=1,
    orbit_visible=False,
    orbit_line_interpolate=5,
    N_star_orbits=10,
    color=[255, 220, 200],
    size_star=1,
    scatter_kwargs={},
):
    """"""
    if orbit_line_interpolate > 1:
        x = np.linspace(0, 1, len(orbit_x))
        x_smooth = np.linspace(0, 1, len(orbit_x) * orbit_line_interpolate)
        kind = 'quadratic'
        orbit_x_line = scipy.interpolate.interp1d(x, orbit_x, kind)(x_smooth)
        orbit_y_line = scipy.interpolate.interp1d(x, orbit_y, kind)(x_smooth)
        orbit_z_line = scipy.interpolate.interp1d(x, orbit_z, kind)(x_smooth)
    else:
        orbit_x_line = orbit_x
        orbit_y_line = orbit_y
        orbit_z_line = orbit_z
    line = ipv.plot(orbit_x_line, orbit_y_line, orbit_z_line, visible=orbit_visible)
    x = np.repeat(orbit_x, N_stars).reshape((-1, N_stars))
    y = np.repeat(orbit_y, N_stars).reshape((-1, N_stars))
    z = np.repeat(orbit_z, N_stars).reshape((-1, N_stars))
    xr, yr, zr = np.random.normal(0, scale=sigma_r, size=(3, N_stars))  # +
    r = np.sqrt(xr ** 2 + yr ** 2 + zr ** 2)

    for i in range(N_stars):
        a = np.linspace(0, 1, x.shape[0]) * 2 * np.pi * N_star_orbits
        xo = r[i] * np.sin(a)
        yo = r[i] * np.cos(a)
        zo = a * 0
        xo, yo, zo = np.dot(_randomSO3(), [xo, yo, zo])
        # print(x.shape, xo.shape)
        x[:, i] += xo
        y[:, i] += yo
        z[:, i] += zo

    sprite = ipv.scatter(
        x, y, z, texture=radial_sprite((64, 64), color), marker='square_2d', size=size_star, **scatter_kwargs
    )
    with sprite.material.hold_sync():
        sprite.material.blending = pythreejs.BlendingMode.CustomBlending  # pylint: disable=no-member
        sprite.material.blendSrc = pythreejs.BlendFactors.SrcColorFactor  # pylint: disable=no-member
        sprite.material.blendDst = pythreejs.BlendFactors.OneFactor  # pylint: disable=no-member
        sprite.material.blendEquation = 'AddEquation'
        sprite.material.transparent = True
        sprite.material.depthWrite = False
        sprite.material.alphaTest = 0.1
    return sprite, line
548403
def get_position(self, dt):
        """"""
        return self.sx + self.dx * dt, self.sy + self.dy * dt
965000
def _fit_s(D, a0, logp, tol=1e-7, maxiter=1000):
    ''''''
    N, K = D.shape
    s1 = a0.sum()
    m = a0 / s1
    mlogp = (m*logp).sum()
    for i in xrange(maxiter):
        s0 = s1
        g = psi(s1) - (m*psi(s1*m)).sum() + mlogp
        h = _trigamma(s1) - ((m**2)*_trigamma(s1*m)).sum()

        if g + s1 * h < 0:
            s1 = 1/(1/s0 + g/h/(s0**2))
        if s1 <= 0:
            s1 = s0 * exp(-g/(s0*h + g)) # Newton on log s
        if s1 <= 0:
            s1 = 1/(1/s0 + g/((s0**2)*h + 2*s0*g)) # Newton on 1/s
        if s1 <= 0:
            s1 = s0 - g/h # Newton
        if s1 <= 0:
            raise Exception('Unable to update s from {}'.format(s0))

        a = s1 * m
        if abs(s1 - s0) < tol:
            return a

    raise Exception('Failed to converge after {} iterations, s is {}'
            .format(maxiter, s1))
687319
def fitcircle(n, x, y): 
# n points, x points, y points
    """"""
    # adding in normalize vectors step
    #x = numpy.array(x) / max(x)
    #y = numpy.array(y) / max(y)
    #
    
    sx, sx2, sx3, sy, sy2, sy3, sxy, sxy2, syx2 = (0,) * 9
    print(type(sx), sx)
    for i in range(n):
        sx = sx + x[i]
        sx2 = sx2 + x[i]**2
        sx3 = sx3 + x[i]**3
        sy = sy + y[i]
        sy2 = sy2 + y[i]**2
        sy3 = sy3 + y[i]**3
        sxy = sxy + x[i] * y[i]
        sxy2 = sxy2 + x[i] * y[i]**2
        syx2 = syx2 + y[i] * x[i]**2

    A = n * sx2 - sx**2
    B = n * sxy - sx*sy
    C = n * sy2 - sy**2
    D = 0.5 * (n * sxy2 - sx * sy2 + n * sx3 - sx * sx2)
    E = 0.5 * (n * syx2 - sy * sx2 + n * sy3 - sy * sy2)
    # values check out up to here

    xo = old_div((D * C - B * E), (A * C - B**2))
    yo = old_div((A * E - B * D), (A * C - B**2))
    print("xo", xo)
    print("yo", yo)

    r = 0
    for z in range(n):
        r = r + old_div(numpy.sqrt( (x[z]-xo)**2 + (y[z]-yo)**2 ), n)

    if xo <= numpy.mean(x) and yo <= numpy.mean(y):
        k = old_div(-1.,r)
    else:
        k = old_div(1.,r)

    SSE = lib_k.get_SSE(xo, yo, r, x, y)
    print("r", r)
    return k, xo, yo, SSE