text: # -*- coding: utf-8 -*-
  
  import numpy as np
  import scipy.io as sio
  from scipy.spatial.distance import cdist
  
  
  def ITQ(V, n_iter):
      # Main function for  ITQ which finds a rotation of the PCA embedded data
      # Input:
      #     V: nxc PCA embedded data, n is the number of images and c is the code length
      #     n_iter: max number of iterations, 50 is usually enough
      # Output:
      #     B: nxc binary matrix
      #     R: the ccc rotation matrix found by ITQ
      # Publications:
      #     Yunchao Gong and Svetlana Lazebnik. Iterative Quantization: A
      #     Procrustes Approach to Learning Binary Codes. In CVPR 2011.
      # Initialize with a orthogonal random rotation initialize with a orthogonal random rotation
  
      bit = V.shape[1]
      np.random.seed(n_iter)
      R = np.random.randn(bit, bit)
      U11, S2, V2 = np.linalg.svd(R)
      R = U11[:, :bit]
  
      # ITQ to find optimal rotation
      for iter in range(n_iter):
          Z = np.matmul(V, R)
          UX = np.ones((Z.shape[0], Z.shape[1])) * -1
          UX[Z >= 0] = 1
          C = np.matmul(np.transpose(UX), V)
          UB, sigma, UA = np.linalg.svd(C)
          R = np.matmul(UA, np.transpose(UB))
  
      # Make B binary
      B = UX
      B[B < 0] = 0
  
      return B, R
  
  
  def compressITQ(Xtrain, Xtest, n_iter=50):
  
      # compressITQ runs ITQ
      # Center the data, VERY IMPORTANT
      Xtrain = Xtrain - np.mean(Xtrain, axis=0, keepdims=True)
      Xtest = Xtest - np.mean(Xtest, axis=0, keepdims=True)
  
      # PCA
      C = np.cov(Xtrain, rowvar=False)
      l, pc = np.linalg.eigh(C, 'U')
      idx = l.argsort()[::-1]
      pc = pc[:, idx]
      XXtrain = np.matmul(Xtrain, pc)
      XXtest = np.matmul(Xtest, pc)
  
      # ITQ
      _, R = ITQ(XXtrain, n_iter)
  
      Ctrain = np.matmul(XXtrain, R)
      Ctest = np.matmul(XXtest, R)
  
      Ctrain = Ctrain > 0
      Ctest = Ctest > 0
  
      return Ctrain, Ctest
  
  
  if __name__ == '__main__':
      # out_vec = np.random.randint(10, size=(3000, 64))
      # out_vec = np.random.rand(3000, 64)
      # Y, pc, R = runITQ(out_vec, 64)
      mat_content = sio.loadmat('/home/adutta/test.mat')
      X = mat_content['X']
      C = compressITQ(X, 10)
      D = cdist(C, C, metric='hamming')
      print(D)
  