'''
i'm investigating different techniques and subroutines for performing 
association on the celestial sphere

Jamie Roberts
'''

# imports
import numpy as np
from matplotlib import pyplot as plt

for i in range(100):
    
    # initialize
    n      = 10
    nmatch = np.int(n*0.8)
    thresh = 0.01
    shift  = np.random.rand(1,2)*0.1
    noise  = np.random.randn(n,2)*0.001
    a      = np.random.rand(n,2)
    b      = a + shift
    na,_   = a.shape
    nb,_   = b.shape
    
    
    # scramble b
    if True:
        idx = np.random.permutation(n)
        b   = b[idx,:]
    
    
    
    # remove some points from b
    if False:
        nremove = 2
        b      = b[:-2,:]
        idx    = idx[:-2]
        nb    -= nremove
    
    
    # add some extra points to b
    if False:
        nadd        = 3
        newb        = np.random.rand(nb+nadd,2)
        newb[:nb,:] = b
        b           = newb + 0.0
        nb,_        = b.shape
    
    
    
    # plot
    if False:
        plt.figure(1)
        plt.clf()
        plt.plot(a[:,0],a[:,1],'rx')
        plt.plot(b[:,0],b[:,1],'b+')
    
    
    
    
    # form distances from each group to every other member of their group
    disa = np.zeros((na,na))
    for i1,a1 in enumerate(a):
        for i2,a2 in enumerate(a):
            disa[i1,i2] = np.sum((a1-a2)**2)
    disb = np.zeros((nb,nb))
    for i1,b1 in enumerate(b):
        for i2,b2 in enumerate(b):
            disb[i1,i2] = np.sum((b1-b2)**2)
    
    
    
    # compare the feature sets to find the best association
    idxa = np.zeros(na)
    maxcount = np.zeros(na)
    for ia1,a1 in enumerate(disa):
        count = np.zeros(nb)
        for ib1,b1 in enumerate(disb):
            tmpdiff = np.zeros((na,nb))
            for ia2 in range(na):
                for ib2 in range(nb):
                    tmpdiff[ia2,ib2] = np.abs(a1[ia2]-b1[ib2])
            
            count[ib1] = np.sum(np.min(tmpdiff,axis=1)<thresh)
            
        idxa[ia1] = np.argmax(count)
    
    
    
    #print idx
    #print np.argsort(idxa)
    print i, idx - np.argsort(idxa)
    if np.sum(np.abs(idx - np.argsort(idxa))) > 0:
        stop





