Update mm_rework to match mm_rework_lingfeng for comparison

flystar/analysis.py

@@ -1,19 +1,12 @@
 import numpy as np
 import pylab as plt
-from flystar import starlists
-from flystar import startables
-from flystar import align
-from flystar import match
-from flystar import transforms
+from flystar import starlists, match
 from astropy import table
 from astropy.table import Table, Column
 from astropy.coordinates import SkyCoord
 from astropy import units as u
-from astropy.wcs import WCS
 from astroquery.gaia import Gaia
-from astroquery.mast import Observations, Catalogs
 import pdb, copy
-import math
 from scipy.stats import f
 
 ##################################################
@@ -49,7 +42,7 @@ def query_gaia(ra, dec, search_radius=30.0, table_name='gaiadr3'):
     search_radius *= u.arcsec
 
     Gaia.ROW_LIMIT = 50000
-    gaia_job = Gaia.cone_search_async(target_coords, search_radius, table_name = table_name + '.gaia_source')
+    gaia_job = Gaia.cone_search_async(target_coords, radius=search_radius, table_name=table_name + '.gaia_source')
     gaia = gaia_job.get_results()
 
     #Change new 'SOURCE_ID' column header back to lowercase 'source_id' so all subsequent functions still work:
@@ -475,85 +468,6 @@ def startable_subset(tab, idx, mag_trans=True, mag_trans_orig=False):
 # Old codes.
 ##################################################
 
-def calc_chi2(ref_mat, starlist_mat, transform, errs='both'):
-    """
-    calculate the chi2 and reduced chi2 of the position 
-    between two matched starlists.
-    Input:
-    ref_mat: astropy table
-        Reference starlist only containing matched stars that were used in the
-        transformation. Standard column headers are assumed.
-
-    starlist_mat: astropy table
-        Transformed starlist only containing the matched stars used in
-        the transformation. Standard column headers are assumed.
-
-    transform: transformation object
-        Transformation object of final transform. Used in chi-square
-        determination
-
-    errs: string; 'both', 'reference', or 'starlist'
-        If both, add starlist errors in quadrature with reference errors.
-
-        If reference, only consider reference errors. This should be used if the starlist
-        does not have valid errors
-
-        If starlist, only consider starlist errors. This should be used if the reference
-        does not have valid errors
-
-    Output:
-    chi_sq: float
-        chi2 = sum (diff_x**2 / xerr**2 + diff_y**2 /yerr**2)
-    chi_sq_red: float
-        reduced chi2 = chi2/ degree of freedom
-    deg_freedom: int
-        degree of freedom
-
-    """
-    diff_x = ref_mat['x'] - starlist_mat['x']
-    diff_y = ref_mat['y'] - starlist_mat['y']
-
-    # Set errors as per user input
-    if errs == 'both':
-        xerr = np.hypot(ref_mat['xe'], starlist_mat['xe'])
-        yerr = np.hypot(ref_mat['ye'], starlist_mat['ye'])
-    elif errs == 'reference':
-        xerr = ref_mat['xe']
-        yerr = ref_mat['ye']
-    elif errs == 'starlist':
-        xerr = starlist_mat['xe']
-        yerr = starlist_mat['ye']
-
-
-    # For both X and Y, calculate chi-square. Combine arrays to get combined
-    # chi-square
-    chi_sq_x = diff_x**2. / xerr**2.
-    chi_sq_y = diff_y**2. / yerr**2.
-
-    chi_sq = np.append(chi_sq_x, chi_sq_y)
-
-    # Calculate degrees of freedom in transformation
-    num_mod_params = calc_nparam(transform)
-    deg_freedom = len(chi_sq) - num_mod_params
-
-    # Calculate reduced chi-square
-    chi_sq = np.sum(chi_sq)
-    chi_sq_red = chi_sq / deg_freedom
-
-    return chi_sq, chi_sq_red, deg_freedom
-
-
-def calc_nparam(transformation):
-    """
-    calculate the degree of freedom for a transformation
-    """
-    # Read transformation: Extract X, Y coefficients from transform
-    if transformation.__class__.__name__ == 'four_paramNW':
-        nparam = 4
-    elif transformation.__class__.__name__ == 'PolyTransform':
-        order = transformation.order
-        nparam = (order+1) * (order+2) 
-    return nparam
 
 def calc_F(red_chi2_1, red_chi2_2, v1, v2):
     """

flystar/examples.py

@@ -1,8 +1,4 @@
-from flystar import transforms
-from flystar import match
-from flystar import align
-from flystar import starlists
-from flystar import plots
+from . import transforms, match, align, starlists, plots
 import numpy as np
 import copy
 import pdb

flystar/match.py

@@ -1,5 +1,4 @@
 import numpy as np
-from flystar import starlists, transforms, startables, align
 from collections import Counter
 from scipy.spatial import cKDTree as KDT
 from astropy.table import Column, Table
@@ -462,210 +461,3 @@ def add_votes(votes, match1, match2):
     votes.flat[unique_idx] += deltas
 
     return
-
-
-def generic_match(sl1, sl2, init_mode='triangle',
-                  model=transforms.PolyTransform, order_dr=(1, 1.0),
-                  dr_final=1.0,
-                  xy_match=(None, None, None, None, None, None, None, None),
-                  m_match=(None, None, None, None), sigma_match=None,
-                  n_bright=100, verbose=True, **kwargs):
-    """
-    Finds the transformation between two starlists using the first one
-    as reference frame. Different matching methods can be used. If no
-    transformation is found, it returns an error message.
-
-
-    Parameters
-    sl1 : StarList
-        starlist used for reference frame
-    sl2 : StarList
-        starlist transformed
-    init_mode : str
-        Initial matching method.
-        If 'triangle', uses the blind triangle method.
-        If 'match_name', uses match by name
-        If 'load', uses the transformation from a loaded file
-    model : str
-        Transformation model to be used with the 'triangle' initial mode
-    poly_order : int
-        Order of the transformation model
-    order_dr : int, float [n, 2]
-        Combinations of polinomial order (first column) and search radius
-        (second column) to refine the transformation. Rows are executed in
-        orders
-    dr_final: float
-        Search radius used for the final matching
-    n_bright : int
-        Number of bright stars used in the initial blind triangles matching
-    xy_match : array
-        Area of the images to remove in the matching [reference catalog min x,
-        reference catalog max x, reference catalog min y, reference catalog max y,
-        transformed catalog min x, transformed catalog max x,
-        transformed catalog min y, transformed catalog max y]. Use None for values not used.
-    m_match : array
-        Magnitude limits of matching stars used to find transformations
-        [reference catalog min mag, reference catalog max mag, transformed
-        catalog min mag, transformed catalog max mag]. Use None for values not
-        used
-    sigma_match : array
-        Number of Deltap movement sigmas [0] used for sigma-cutting matched
-        stars for a number of times [1]. Use None for no sigma-cut. The last
-        polynomial order and search radius in 'order_dr' are used
-    transf_file : str
-        File name and path of the transformation file used with the 'load'
-        init_mode
-    verbose : bool, optional
-        Prints on screen information on the matching
-
-    Returns
-    -------
-    transf : Transform2D
-        Transformation of the second starlist respect to the first
-    st : StarTable
-        Startable of the two matched catalogs
-
-    """
-
-    #  Check the input StarLists and transform them into astropy Tables
-    if not isinstance(sl1, starlists.StarList):
-        raise TypeError("The first catalog has to be a StarList")
-    if not isinstance(sl2, starlists.StarList):
-        raise TypeError("The second catalog has to be a StarList")
-
-    #  Find the initial transformation
-    if init_mode == 'triangle': #  Blind triangles method
-
-        #  Prepare the reduced starlists for matching
-        sl1_cut = copy.deepcopy(sl1)
-        sl2_cut = copy.deepcopy(sl2)
-        sl1_cut.restrict_by_value(x_min=xy_match[0], x_max=xy_match[1],
-                                  y_min=xy_match[2], y_max=xy_match[3])
-        sl2_cut.restrict_by_value(x_min=xy_match[4], x_max=xy_match[5],
-                                  y_min=xy_match[6], y_max=xy_match[7])
-        sl1_cut.restrict_by_value(m_min=m_match[0], m_max=m_match[1])
-        sl2_cut.restrict_by_value(m_min=m_match[2], m_max=m_match[3])
-
-        # Find the transformation
-        # TODO: test 'initial_align' with StarList input
-        transf = align.initial_align(sl1_cut, sl2_cut, briteN=n_bright,
-                                     transformModel=model, order=order_dr[0]) #order_dr[i_loop][0] ?
-
-    elif init_mode == 'match_name': #  Name match
-        sl1_idx_init, sl2_idx_init, _ = starlists.restrict_by_name(sl1, sl2)
-        transf = model(sl2['x'][sl2_idx_init], sl2['y'][sl2_idx_init],
-                       sl1['x'][sl1_idx_init], sl1['y'][sl1_idx_init],
-                       order=int(order_dr[0][0]))
-
-    elif init_mode == 'load': #  Load a transformation file
-        transf = transforms.Transform2D.from_file(kwargs['transf_file'])
-
-    else: #  None of the above
-        raise TypeError("Unrecognized initial matching method")
-
-    # Restrict the matching catalogs
-    sl1_match = copy.deepcopy(sl1)
-    sl2_match = copy.deepcopy(sl2)
-    sl1_match.restrict_by_value(m_min=m_match[0], m_max=m_match[1])
-    sl2_match.restrict_by_value(m_min=m_match[2], m_max=m_match[3])
-
-    #  Refine the transformation
-    if sigma_match:
-        order_dr_len = len(order_dr)
-
-        for i_loop in range(sigma_match[1]):
-            order_dr = np.vstack((np.array(order_dr), np.array(order_dr[-1])))
-
-    for i_loop in range(len(order_dr)):
-
-        #  Transform and match the catalog to the reference frame
-#        sl2_idx, sl1_idx = align.transform_and_match(sl2_match, sl1_match, transf,
-#                                                     dr_tol=order_dr[i_loop][1],
-#                                                     verbose=verbose)
-
-        sl2_idx, sl1_idx = align.transform_and_match(sl2_match, sl1_match, transf,
-                                                     dr_tol=order_dr[1],
-                                                     verbose=verbose)
-
-        #  Transform the catalog to the reference frame
-        sl2_transf_match = align.transform_from_object(sl2_match, transf)
-
-        # Sigma-rejection
-        if sigma_match and (i_loop >= order_dr_len):
-            resid = np.sqrt((sl1_match['x'][sl1_idx] -
-                            sl2_transf_match['x'][sl2_idx])**2 +
-                            (sl1_match['y'][sl1_idx] -
-                            sl2_transf_match['y'][sl2_idx])**2)
-            sl1_idx = sl1_idx[resid <= (sigma_match[0] * np.std(resid))]
-            sl2_idx = sl2_idx[resid <= (sigma_match[0] * np.std(resid))]
-
-        # Test section to observe the matching catalogs before refining the transformation
-        """
-        from matplotlib import pyplot
-
-        _, axarr = pyplot.subplots(nrows=1, ncols=1, figsize=(10,10))
-        axarr.scatter(sl1_match['x'][sl1_idx], sl1_match['y'][sl1_idx])
-        xlim = axarr.get_xlim()
-        ylim = axarr.get_ylim()
-
-        _, axarr = pyplot.subplots(nrows=1, ncols=1, figsize=(10, 10))
-        axarr.scatter(sl2_transf_match['x'][sl2_idx], sl2_transf_match['y'][sl2_idx])
-        axarr.set_xlim(xlim)
-        axarr.set_ylim(ylim)
-        """
-
-        #  Find a better transformation
-        transf, _ = align.find_transform(sl2_match[sl2_idx],
-                                         sl2_transf_match[sl2_idx],
-                                         sl1_match[sl1_idx], transModel=model,
-                                         order=order_dr[0], verbose=verbose)
-#                                         order=int(order_dr[i_loop][0]), verbose=verbose)
-
-        # This section was used for testing transformations with normalized
-        # coordinates. Only several catalogs had reduced residuals when using
-        # high order polynomials (>3), some of them became unstable
-        """sl1_match_norm = sl1_match[sl1_idx]
-        sl2_match_norm = sl2_match[sl2_idx]
-        sl2_transf_match_norm = sl2_transf_match[sl2_idx]
-        mm = max(max(sl1_match_norm['x']), max(sl1_match_norm['y']),
-                 max(sl2_transf_match_norm['x']), max(sl2_transf_match_norm['y']))
-        sl1_match_norm['x'] = sl1_match_norm['x'] / mm
-        sl1_match_norm['y'] = sl1_match_norm['y'] / mm
-        sl2_match_norm['x'] = sl2_match_norm['x'] / mm
-        sl2_match_norm['y'] = sl2_match_norm['y'] / mm
-        sl2_transf_match_norm['x'] = sl2_transf_match_norm['x'] / mm
-        sl2_transf_match_norm['y'] = sl2_transf_match_norm['y'] / mm
-        transf, _ = align.find_transform(sl2_match_norm, sl2_transf_match_norm,
-                                         sl1_match_norm, transModel=model,
-                                         order=poly_order, verbose=verbose)
-        c_exp = np.zeros(len(transf.px._parameters))
-
-        for i_c in range(len(transf.px._parameters)):
-            c_exp[i_c] = int(transf.px._param_names[i_c][1:].split('_')[0]) +\
-                         int(transf.px._param_names[i_c][1:].split('_')[1])
-
-        c_corr = mm ** (1 - c_exp)
-        transf.px._parameters = transf.px._parameters * c_corr
-        transf.py._parameters = transf.py._parameters * c_corr"""
-
-    # Do the final transformation and matching using
-    sl2_idx, sl1_idx = align.transform_and_match(sl2, sl1, transf, dr_tol=dr_final,
-                                                 verbose=verbose)
-    #  StarTable output
-    sl2_transf = align.transform_from_object(sl2, transf)
-    unames = np.array(range(len(sl1_idx)))
-    st = startables.StarTable(name=unames,
-         x=np.column_stack((np.array(sl1['x'][sl1_idx]), np.array(sl2_transf['x'][sl2_idx]))),
-         y=np.column_stack((np.array(sl1['y'][sl1_idx]), np.array(sl2_transf['y'][sl2_idx]))),
-         m=np.column_stack((np.array(sl1['m'][sl1_idx]), np.array(sl2_transf['m'][sl2_idx]))),
-         ep_name=np.column_stack((np.array(sl1['name'][sl1_idx]), np.array(sl2_transf['name'][sl2_idx]))))
-#         ep_name=np.column_stack((np.array(sl1['name'][sl1_idx]), np.array(sl2_transf['name'][sl2_idx]))),
-#         list_times=[sl1.meta['list_time'], sl2.meta['list_time']],
-#         list_names=[sl1.meta['list_name'], sl2.meta['list_name']])
-
-    for col in sl1.colnames:
-        if col in sl2.colnames:
-            if col not in ['name', 'x', 'y', 'm']:
-                st.add_column(Column(np.column_stack((np.array(sl1[col][sl1_idx]),np.array(sl2_transf[col][sl2_idx]))), name=col))
-
-    return transf, st