import numpy as np 
import matplotlib.pyplot as plt #We might need this


#First, let us load the data
#Catalog from HSC 
cat_hsc = np.loadtxt('./Catalog_HSC.csv')
x_hsc = cat_hsc[:,0]
y_hsc = cat_hsc[:,1]
#Catalog from HST
cat_hst = np.loadtxt('./Catalog_HST.csv')
x_hst = cat_hst[:,0]
y_hst = cat_hst[:,1]

"""
Explanation: Exercise numpy
The ultimate goal of this exercise is to compare the position of stars in a patch of sky as measured in two different surveys. The main task at hand is to identify matching positions of stars between the surveys. For this, we will need to compare the positions of all stars in one survey to the position of all stars in the other survey. This task can be extremely time consuming if not implemented properly, we will therefore use this to compare different coding style and their impact on computation time. 
If time allows, we will move on to represent the results of our analysis in a meaningfull way.
End of explanation
"""


#First, check the number of stars in each survey:
ns_hst = #fill in
ns_hsc = #...
#Print the result
print()

#This is a graphic representation of our data content:
%matplotlib qt
plt.title('star catalogs in COSMOS')
plt.plot(x_hsc, y_hsc, 'or', label = 'hsc catalog')
plt.plot(x_hst, y_hst, 'ob', label = 'hst catalog')
plt.legend()
plt.xlabel('ra')
plt.ylabel('dec')
plt.show()

"""
Explanation: Check that the loaded data are consistent with what we expect: (ra, dec) coordinates of the same patch of sky
End of explanation
"""


def distance(point1, point2):
    ''' Returns the distance between two points with coordinates (x,y).
    
    Parameters
    ----------
    point1: list
        2D coordinates of a point 
    point2: list
        2D coordinates of a point 
        
    Returns
    -------
    d: float
        the distance between point1 and point2
    '''
    
    return

"""
Explanation: To begin with, let's write a function that returns the algebraic distance between two points
End of explanation
"""


point1 = [x_hst[0], y_hst[0]]
point2 = [x_hsc[0], y_hsc[0]]
print(distance(point1, point2))
# Answer should be 0.6648994838877168

"""
Explanation: Now let's test it by comparing the distance between the first point of each dataset.
End of explanation
"""


def point_to_points_distance(point, coordinates):
    ''' Returns the distance between one point and all the points in coordinates.
    
    Parameters
    ----------
    point: list
        2D coordinates of a point 
    coordinates: list
        set of N 2D coordinates stored in a list with shape Nx2
        
    Returns
    -------
    d: list
        the distance between point and each point in coordinates in an array with size N
    '''
    #Declaring an empty list
    d = []
    for c in coordinates:
        # for each point in coordinates, take the distance to point and concatenate it to d 
        d.append(distance(point, c))
    #make d a numpy array and return it
    return np.array(d)

"""
Explanation: Let's take it one step further and compare the distance between one point and a set of points
End of explanation
"""


coords = np.concatenate((x_hsc[:10,None], y_hsc[:10,None]), axis = 1)
print(point_to_points_distance(point1, coords))
# The answer should look like [0.66489948 0.4628197  0.39672485 0.43854084 0.32165335 0.30223269
# 0.65765909 0.65411548 0.6474303  0.79301678]

"""
Explanation: Let's test it on the first 10 points in the HSC catalog and the first point of the HST catalog
End of explanation
"""


def your_function(coord1, coord2): # Choose an adequate name for your function
    ''' Returns the distance between points in two sets of coordinates.
    
    Parameters
    coord1: array
        array of size Nx2 that contains the [x, y] positions of a catalog 
    coord2: array
        array of size Mx2 that contains the [x, y] positions of a catalog 
    
    Returns
    dist: array
        array of size NxM that contains the euclidean distances between points in the two datasets
    '''
    
    return


"""
Explanation: Now let's get to work. We would like to associate stars in one survey to their conterpart (if it exists) in the other survey. We will start by comparing the positions between each point of one survey to the position of each point in the other survey.
First, write a function that takes two sets of coordinates (hsc and hst) and returns the distance from each point of one survey to each point of the other, such that the output should be an array of size (n_hst x n_hsc) or (n_hsc x n_hst).
PS: if you have several (different) ideas about how to implement this, feel free to code them!
End of explanation
"""


# In order not to spend the whole evening here, let us reduce the dataset size:
#Select stars in hsc in the frame: 150.0<x<150.1 and 2.0<y<2.1
loc_hsc = #please fill these
x_hsc_exp = x_hsc[loc_hsc]
y_hsc_exp = y_hsc[loc_hsc]

loc_hst = #And that
x_hst_exp = x_hst[loc_hst]
y_hst_exp = y_hst[loc_hst]
#Once you are done with the exercise, feel free to try with larger selections to see how it impacts computation time

import distances as dt
# Insert the names of your functions in the following array:
methods = [your_function, dt.double_loop, dt.with_indices, dt.one_loop, dt.one_loop_reverse, dt.scipy_version, dt.newaxis_magic]
#An empty variable to store computation time
timers = []
# Making sets of coordinates of size Nx2 to feed your functions with the right format
c2 = np.concatenate((x_hst_exp[:,None], y_hst_exp[:,None]), axis = 1)
c1 = np.concatenate((x_hsc_exp[:,None], y_hsc_exp[:,None]), axis = 1)

for f in methods:
    print(f.__name__)
    r = %timeit -o f(c1, c2)
    timers.append(r)

#View the results:
plt.figure(figsize=(10,6))
plt.bar(np.arange(len(methods)), [r.best*1000 for r in timers], log=True)  # Set log to True for logarithmic scale
plt.xticks(np.arange(len(methods))+0.2, [f.__name__ for f in methods], rotation=30)
plt.xlabel('Method')
plt.ylabel('Time (ms)')
plt.yscale('log')
plt.show()

"""
Explanation: Now, let us take a look at the computation times:
End of explanation
"""


#Let us compute the distances as we did before, but this time, with the whole dataset.
#Of course, a fast method is to be prefered

c1 = #Please fill these. Same as before but with all the dataset
c2 = #


def get_match(coord_ref, coord2, rad):
    '''
    matches coordinates of stars between two datasets and computes the distance between the position of the stars in the 2 datasets

    Parameters
    coord_ref: numpy array (Nx2)
        coordinates (ra, dec) of stars in a FoV from a given dataset
    coord2: numpy array (Mx2)
        coordinates (ra dec) of stars in the same FoV in an other dataset
    rad: float
        radius (deg) around stars in coord_ref where to find a corresponding star in coord2
            
    Returns
    modulus:numpy array (N')
        containing the distance between matching stars
    v_coord: numpy array(N',2)
        coordinates in the coord_ref set of matching stars
            

    '''
    #Declare two empty arrays to store the coordinates and distances.
    #...
    s = np.size(coord_ref[:,0])#This is just for representation
    print('number of points in reference catalog: {0}'.format(s))
    #for each star in coord_ref
    for i,c in enumerate(coord_ref):

        #This is just here to keep track of the algorithm's progression
        if i % 3000 == 0:
            print('point number {0} out of {1}'.format(i, s))
        #compute the distance from c to all stars in coord2
        r = #...
        #Find the closest star from coord 2 to c
        loc = #...

        #Make sure that there is only one star matching (it can happen that 2 match)
        #Here I just arbitrarily pick one, but you can find a way to discard these stars
        if np.size(loc) > 1:
            loc = loc[0]

        #record the distance between matching stars
        rmin = #...
        
        #Check whether the closest distance is smaller than rad
        if #...:
            #if yes, place the coordinates and the distance in an array
            #... tip: use append()

    return #...

# Use your function 
coord , r = get_match(c1, c2, 0.3/3600.)

"""
Explanation: Identifying matching stars (optional)
Now that we know all the distances, let us find the stars in each datasets that correspond to one another.
This is done by finding, for each star, the minimum distance recorded between the two datasets.
One problem that arises with deriving an array that computes all the distances is that we end up with a very LARGE array, which becomes impractical for fast computations. Instead, we will modify one of the previous functions so that it returns the coordinates of stars that have a match in both datasets along with their distance.
Because all stars in a given set do not have a counter part in the other, we will only accept a match if the minimum distance between two points is smaller than 0.17 arcseconds (the size of an HSC pixel).
In other words, for each star in one dataset, find the star in the other dataset that is the closest (minimum distance), check wether that star is closer that 0.17 arcseconds, if yes, store its coordinates along with the computed distance. At the end of the function, return arrays with the matching star coordinates and their distance to their match in the other dataset.
End of explanation
"""


#Spatial distribution of distances
plt.title('distribution of distances accross the FoV')
#...

#Global representation
#...

"""
Explanation: Now I would like to have a representation for the work we have done that informs me about what is in my datasets. Namely, what is the error on star positions between the two datasets? We would like to have a global view of this error but also an impression of the error as a function of the position on the field. For the latter, I suggest you use the 'scatter' function from matplotlib.
End of explanation
"""


