def xyshift(path):
      path1=path[24:70]
      path1=path1[0:29]
      f1=h5py.File(path1[0],'r')
      f2=f1["MAPS"]["XRF_roi"][...]


      ## element to use for crosscorrelation horizontally
      element1=15
      ## element to use for xcor vertically
      element2=18
      ## element for output
      element3=15

      matrix1=zeros([len(path1),24,38,151])
      matrix2=zeros([len(path1),24,38,151])
      matrix3=zeros([len(path1),24,38,151])
      theta=zeros(len(path1))
      
      for i in arange(len(path1)):
            fp=h5py.File(path1[i],'r')["MAPS"]
            f=fp["XRF_roi"][...]
            matrix2[i,:,:,:]=f
            matrix3[i,:,:,:]=f
            matrix1[i,:,:,:]=f
            temp=fp["extra_pvs_as_csv"][99]
            theta[i]=temp[temp.rfind(",")+2:]

      matrix2[:,element3,:,:]=np.roll(matrix2[:,element3,:,:],40,axis=2)
      matrix3[:,element3,:,:]=np.roll(matrix1[:,element3,:,:],40,axis=2)     

      shift=zeros([2,len(path1)],dtype=int)

      for i in arange(len(path1)-1):
            cor2d=correlate2d(matrix1[i,element1,:,:],matrix1[i+1,element1,:,:],fillvalue=np.average(matrix1[i,element1,:,:]))
            b1=np.where(cor2d==cor2d.max())
            cor2d=correlate2d(matrix1[i,element2,:,:],matrix1[i+1,element2,:,:],fillvalue=np.average(matrix1[i,element2,:,:]))
            b2=np.where(cor2d==cor2d.max())
            x1=b1[1][0]
            y1=b2[0][0]

            x=151
            y=38

            for j in arange(len(path1)-1-i):
                  matrix2[i+j+1,element3,:,:]=np.roll(matrix2[i+j+1,element3,:,:],x1-x+1,axis=1)
                  matrix2[i+j+1,element3,:,:]=np.roll(matrix2[i+j+1,element3,:,:],y1-y+1,axis=0)
                  
            print y1,x1
 
            shift[0,i+1]=y1
            shift[1,i+1]=x1
##      for i in arange(len(path1)):
##            j=Image.fromarray(matrix2[i,element2,8:30,:].astype(np.float32))
##            if i>=10:
##                  j.save("/Users/youngpyohong/Documents/Work/2014-2/hong/hong/projections/shifted"+str(i)+".tiff")
##            else:
##                  j.save("/Users/youngpyohong/Documents/Work/2014-2/hong/hong/projections/shifted0"+str(i)+".tiff")
##      ## wkatl
      return shift,matrix2[:,element3,8:30,:],theta

def xcor(data,element1,element2):
      xmatrix=data[element1,:,:,:]
      ymatrix=data[element2,:,:,:]

      finalmatrix=zeros(data.shape,dtype=float32)

      shift=zeros([2,data.shape[1]],dtype=int)
      for i in arange(len(data[0,:,0,0])-1):
            cor2dx=correlate2d(xmatrix[i,:,:],xmatrix[i+1,:,:],fillvalue=np.average(xmatrix[i,:,:]))
            position1=np.where(cor2dx==cor2dx.max())
            cor2dy=correlate2d(ymatrix[i,:,:],ymatrix[i+1,:,:],fillvalue=np.average(ymatrix[i,:,:]))
            position2=np.where(cor2dy==cor2dy.max())

            x1=position1[1][0]
            y1=position2[0][0]

            shift[0,i+1]=y1
            shift[1,i+1]=x1

            x=len(xmatrix[i,0,:])
            y=len(ymatrix[i,:,0])

            for j in arange(len(data[0,:,0,0])-1-i):
                  finalmatrix[i+j+1,:,:,:]=np.roll(data[i+j+1,:,:,:],x1-x+1,axis=2)
                  finalmatrix[i+j+1,:,:,:]=np.roll(data[i+j+1,:,:,:],y1-y+1,axis=1)
            print "Xcor", i
      return finalmatrix, shift

def gradients(image):
  deriv_filter = np.array([[-1.,0.,1.],[-2.,0,2.],[-1.,0.,1.]])

  #derivative in x direction
  deriv_x = correlate2d(image, deriv_filter, mode="same")
  #derivative in y direction
  deriv_y = correlate2d(image, deriv_filter.T, mode="same")

  gradients = np.sqrt(deriv_x**2 + deriv_y**2)
  return gradients

def matchIms():
	bmp=Image.open('C:/Copy/workspace/TyperSharkAI/Images/Play4.png').convert('L')
	shark=Image.open('basic_template.gif')
	shark.load()
	#shark.show()
	bmp=np.array(bmp)
	shark=np.array(shark)
	#Image.fromarray(shark).show()
	Image.fromarray(bmp).show()
	#bmp.show()
	print(bmp.shape)
	print(shark.shape)
	signal.correlate2d(shark,bmp)

def get_pssm_scores(encoded_sequences, pssm):
    encoded_sequences = np.squeeze(encoded_sequences, axis=1)
    num_samples, num_bases, seq_length = np.shape(encoded_sequences)
    scores = np.ones((num_samples, num_bases, seq_length))
    for base_indx in range(num_bases):
        base_pssm = pssm[base_indx].reshape(1, len(pssm[0]))
        fwd_scores = correlate2d(
            encoded_sequences[:, base_indx, :], base_pssm, mode='same')
        rc_base_pssm = pssm[-(base_indx + 1), ::-1].reshape(1, len(pssm[0]))
        rc_scores = correlate2d(
            encoded_sequences[:, base_indx, :], rc_base_pssm, mode='same')
        scores[:, base_indx, :] = np.maximum(fwd_scores, rc_scores)

    return scores.sum(axis=1)

def corr_NVs_no_subset(baseline_image, new_image):
    """
    # Tracks drift by correlating new and old images, and returns shift in pixels
    :param baseline_image: original image
    :param new_image: new (drifted) image. Should be same size as baseline_image in pixels
    :return: shift from baseline image to new image in pixels
    """
    # subtracts mean to sharpen each image and sharpen correlation
    baseline_image_sub = baseline_image - baseline_image.mean()
    new_image_sub = new_image - new_image.mean()

    #takes center part of baseline image
    x_len = len(baseline_image_sub[0])
    y_len = len(baseline_image_sub)
    old_image = baseline_image_sub[(x_len/4):(x_len*3/4),(y_len/4):(y_len*3/4)]

    # correlate with new image. mode='valid' ignores all correlation points where an image is out of bounds. if baseline
    # and new image are NxN, returns a (N/2)x(N/2) correlation
    corr = signal.correlate2d(new_image_sub, baseline_image_sub)
    y, x = np.unravel_index(np.argmax(corr), corr.shape)

    # finds shift by subtracting center of initial coordinates, x_shift = x + (x_len/4) - (x_len/2)
    x_shift = x - (x_len)
    y_shift = y - (y_len)

    #return (x_shift, y_shift) #, corr, old_image --- test outputs
    return (x_shift, y_shift, corr, old_image) # --- test outputs

 def forward(ctx, input, filter, bias):
     # detach so we can cast to NumPy
     input, filter, bias = input.detach(), filter.detach(), bias.detach()
     result = correlate2d(input.numpy(), filter.numpy(), mode='valid')
     result += bias.numpy()
     ctx.save_for_backward(input, filter, bias)
     return torch.from_numpy(result)

 def test_consistency_correlate_funcs(self):
     # Compare np.correlate, signal.correlate, signal.correlate2d
     a = np.arange(5)
     b = np.array([3.2, 1.4, 3])
     for mode in ["full", "valid", "same"]:
         assert_almost_equal(np.correlate(a, b, mode=mode), signal.correlate(a, b, mode=mode))
         assert_almost_equal(np.squeeze(signal.correlate2d([a], [b], mode=mode)), signal.correlate(a, b, mode=mode))

def autoCorr(matrix):
    averagePotential = np.mean(matrix)
    matrix = matrix - averagePotential
    matrix = correlate2d(matrix, matrix)
    shape = matrix.shape
    max_ele = matrix[shape[0] / 2, shape[1] / 2]
    return matrix[shape[0] / 2:, shape[1] / 2:] / max_ele

def main():
  scale = 100 #lower this value to make the correlation go faster
  image = imread2("./waldo.png")
  image = imresize(image, scale)
  template = imread2("./template.png")
  template = imresize(template, scale)
  # make grayscale
  image_gray = grayscale(image)
  template = grayscale(template)
  template_w, template_h = template.shape

  gradients_image = gradients(image_gray)
  gradients_image /= np.linalg.norm(gradients_image.flatten())
  gradients_template = gradients(template)
  gradients_template /= np.sum(gradients_template)

  # use cross correlation
  convolved_gradients = correlate2d(gradients_image, gradients_template, mode="same")

  position = np.argmax(convolved_gradients)
  position_x, position_y = np.unravel_index(position, gradients_image.shape)

  #put a big red dot in the middle of where we found our maxima
  dot_rad = 8
  image[position_x-dot_rad:position_x+dot_rad,position_y-dot_rad:position_y+dot_rad,0] = 255
  image[position_x-dot_rad:position_x+dot_rad,position_y-dot_rad:position_y+dot_rad,1:2] = 0

  imsave("./image_matched.png", image )

 def dynamics(self, data, verbose=0):
    dA = correlate2d(data,array([[0,1,0],[1,-4,1],[0,1,0]]),boundary='wrap')
    dA = dA[1:N+1,1:M+1]
    if verbose: print(dA[0:5,0:5])
 
    data = data + dt*((1+self.alpha*1j)*self.scale*dA + data - (1+1j*self.beta)*data*power(abs(data),2));
    return data

def correlation_retrieval(data,K=None,ARGS=False):
  f_data = np.fft.fft2(data)
  kx = np.fft.fftfreq(f_data.shape[0])
  ky = np.fft.fftfreq(f_data.shape[1])
  # get k-vector, if not give
  if K is None:
    K = get_wave(f_data)

  # get the frequency data information
  Kfreq = kx[K[0]],ky[K[1]]
  # make an X-Y grid
  pi = np.pi
  nx,ny = np.shape(data)
  Y,X = np.meshgrid(np.arange(ny),np.arange(nx),\
      sparse=False,indexing='xy')
  print Kfreq
  # produce a correlation function
  phplot.imageshow(wave_function(X,Y,Kfreq))
  zsum = np.zeros(data.shape)
  Nphi = 1.
  a_phi = np.linspace(0,2*np.pi,Nphi)
  for phi in a_phi:
    zsum += correlate2d(data,wave_function(X,Y,K,phi),mode='same')

  phplot.imageshow(zsum)
  return data

def TestLaws( mgnames, NJ = 100 ):

    # create laws filters
    filts = texture.BuildLawsFilters()
    # allocate for jets
    NI = len( mgnames ) # number of images
    jets = np.zeros( (NJ*NI, 25 ))
    # for each image
    for i in xrange( NI ):
        # load
        # correlate
        #corrs = BruteCorrelate( data, filts )
        data = mgnames[i]+0
        corrs = map( lambda x: correlate2d( data, x ), filts )
        for j in range( 25 ):
            corrs[i] = cspline2d( abs(corrs[i]), 200 )
        corrs = np.array( corrs )
        # extract random jets
        V,H = data.shape
        vs = range( V )
        hs = range( H )
        np.random.shuffle( vs ); np.random.shuffle( hs )
        for j in range( NJ ):
            jets[i*NJ + j] = corrs[:,vs[j], hs[j] ]
    # k-means clustering
    clust, mmb = kmeans.KMeans( NI, jets )
    #return jets
    cffs,evecs = pca.PCA(clust,3)
    cffs = pca.Map2PCA(clust,evecs)
    gnu.Save('Laws_results.txt',cffs)
    return clust,cffs

def test_corr():

    #imshp = (3,2,20,20) # num images, channels, szy, szx
    kshp = (10,2,5,10) # features, channels, szy, szx
    featshp = (3,10,11,11) # num images, features, szy, szx

    theano_correlate2d = get_theano_correlate2d(kshp=kshp,featshp=featshp)

    features = np.random.randn(*featshp)
    kernel = np.random.randn(*kshp)

    output_sz = (featshp[0], kshp[1], kshp[2] + featshp[2] - 1, kshp[3] + featshp[3] - 1)

    scipy_output = np.zeros(output_sz)
    for im_i in range(featshp[0]):
        for im_j in range(kshp[1]):
            for k_i in range(kshp[0]):
                scipy_output[im_i,im_j,:,:] += correlate2d(np.squeeze(features[im_i,k_i,:,:]),np.squeeze(kernel[k_i,im_j,:,:]),mode='full')

    theano_output = theano_correlate2d(features,kernel)

    print 'scipy:', scipy_output.shape
    print 'theano:', theano_output.shape

    np.testing.assert_allclose(scipy_output,theano_output)

    def _convolve(self, imgs, filters):
        assert(imgs.ndim == 3 and filters.ndim == 3)
        assert(imgs.shape[-2] >= filters.shape[-2] and imgs.shape[-1] >= filters.shape[-1])
        assert(filters.shape[-2] == filters.shape[-1] and filters.shape[-1] % 2 != 0)
        lx = filters.shape[-1]//2
        rx = imgs.shape[-1] - lx - 1
        ly = lx
        ry = imgs.shape[-2] - ly - 1
        #print "f " + str(filters.shape[0])
        output = np.zeros((filters.shape[0], rx - lx + 1, ry - ly + 1))
        for f in range(0, filters.shape[0]):
            filter = filters[f]
            filter_map = np.zeros((rx - lx + 1, ry - ly + 1))

            for i in range(0, imgs.shape[0]):
                img = imgs[i]
                convolved = np.zeros((rx - lx + 1, ry - ly + 1))
                #print "convolved shape " + str(convolved.shape)
                #print "lx " + str(lx) + " rx " + str(rx) + " ly " + str(ly) + " ry " + str(ry)
                for x in range(lx, rx + 1):
                    for y in range(ly, ry + 1):
                        subimg = img[y - ly:y + ly + 1:,x - lx:x + lx + 1]
                        convolved[y - ly, x - lx] = (subimg * filter).sum()
                if self.debug:
                    lib_convolved = correlate2d(img, filter, "valid")
                    if not np.all(np.abs(convolved - lib_convolved) < 0.000001):
                        print "Convolved:\n{}\nLib Convolved:\n{}\nFilter:\n{}".format(convolved, lib_convolved, filter)
                        assert(False)

                filter_map += convolved
            output[f]=filter_map
        return output

    def Conv_Forward(self, x, w, b):
        """
        - x: Input data of shape (N, C, H, W)
        - w: Filter weights of shape (F, C, HH, WW)
        - b: Biases, of shape (F,)
        Returns a tuple of:
        - out: Output data.
        - cache: (x, w, b, conv_param)
        """
        out = None
        pad = self.conv_param['pad']
        stride = self.conv_param['stride']
        H = x.shape[2]
        W = x.shape[3]
        HH = w.shape[2]
        WW = w.shape[3]
        N = x.shape[0]
        F = w.shape[0]

        out = np.zeros([N,F,H,W])

        for i in xrange(x.shape[0]): #For every data sample
            for j in xrange(x.shape[1]): #For every color
                temp = x[i,j,:,:] #store layer to make things simpler
                for f in xrange(w.shape[0]): #For every filter
                    filt = w[f,j,:,:]
                    out[i,f,:,:] += signal.correlate2d(temp, filt, mode='same', boundary='fill', fillvalue=0) + b[f] #/x.shape[1]
        cache = (x, w, b, self.conv_param)
        return out, cache

def get_shift_correlation(frame1, frame2):
    #frame2 = frame2[65:200, 75:220]
    frame2 = frame2[48:208, 48:208]
    #init 135, 145
    #second center 201, 219
    save_to_file("frame2.raw", frame2, np.float32)
    return signal.correlate2d(frame1, frame2, mode='same', boundary='symm')

 def __call__(self,x):
     y = np.empty_like(x)
     sh = x.shape
     xf = x.reshape((-1,)+sh[-2:])
     yf = y.reshape((-1,)+sh[-2:])
     for i in range(len(xf)):
         yf[i] = correlate2d(xf[i], np.array([[.0625, .125, .0625], [.125, .25, .125], [.0625, .125, .0625]]), mode='same')
     return y

    def xcorr( self, p=[0,0] ):
        from scipy.signal import correlate2d
        from scipy.ndimage import affine_transform

        im1 = self.im1.copy()
        im2 = self.im2.copy()
        im1_tf = affine_transform(im1, np.identity(2) ,mode='wrap',offset=p )
        result = correlate2d( (im1_tf - im1_tf.mean()), (im2 - im2.mean()), mode='full', boundary='wrap' )
        return np.sum( result ) / result.size

def dynamics(data,verbose=0):
   global beta;
   dA = correlate2d(data,array([[1,1,1],[1,0,1],[1,1,1]]),boundary='wrap');
   dA = dA[1:N+1,1:M+1];

   r = random.rand(N,M);

   data = 2*floor(r*(1+exp(beta*dA)))-1;
   return data;