def fourier_interp(sig_objs):
  nws = len(sig_objs[0].aspec)
  nqs = len(sig_objs) 
  pads = 8*nqs
  aspec = np.zeros([nws, nqs+pads])
  aux = np.zeros(nqs)
  f = open('aspecinterp.dat', 'w')
  for iw in range(nws):
    for i, sig in enumerate(sig_objs):
      tmp = map(float, sig.aspec[iw].split())
#trace over valence manifold:
      aux[i] = sum(tmp[2:9])
#  now have fourier coefficients interpolate back on to dense grid
    aux[:] = dct(aux[:], 2, norm='ortho')
#dct type 3 is the 
    auxd = np.pad(aux, (0, pads), 'constant')
    auxd = dct(auxd[:], 3, norm='ortho')
    for iq in range(len(auxd)):
      aspec[iw][iq] = auxd[iq]
  for iq in range(nqs+pads)[::-1]:
#pads with zeros
    for iw in range(nws):
      print >>f,  -iq, (sig.aspec[iw].split())[0], aspec[iw][iq]
    print >>f, ''
  for iq in range(nqs+pads):
#pads with zeros
    for iw in range(nws):
      print >>f, iq, (sig.aspec[iw].split())[0], aspec[iw][iq]
    print >>f, ''
  f.close()

def lpf(image, sigma, mode=2):
    (Mx, My) = image.shape
    if mode == 1:
        kernel = matlab_style_gauss2D(image.shape, sigma)
        kernel /= numpy.max(kernel)

        if Mx == 1 or My == 1:
            fft = numpy.fft.fft(image)
            fft = numpy.fft.fftshift(fft)
            fft *= kernel
            result = numpy.real(numpy.fft.ifft(numpy.fft.fftshift(fft)))
        else:
            fft = numpy.fft.fftshift(numpy.fft.fft2(image))
            fft *= kernel
            result = numpy.real(numpy.fft.ifft2(numpy.fft.fftshift(fft)))
    elif mode == 2:
        new_dim = 2 * array(image.shape)
        kernel = matlab_style_gauss2D((new_dim[0], new_dim[1]), sigma * 2)
        kernel /= numpy.max(kernel)
        kernel = kernel[Mx:, My:]

        image = image.astype(numpy.double)
        if Mx == 1 or My == 1:
            dct = fftpack.dct(image, type=1)
            dct *= kernel
            result = numpy.real(fftpack.idct(dct, type=1))
        else:
            dct = fftpack.dct(fftpack.dct(image.T, type=2, norm='ortho').T, type=2, norm='ortho')
            dct *= kernel
            result = numpy.real(fftpack.idct(fftpack.idct(dct.T, type=2, norm='ortho').T, type=2, norm='ortho'))
    return result

def energy_99 ():

	files = []
	# directory for read files.
	directory = "/Users/cyan/Desktop/color_hist_py/cropImage_large/"
	for infile in glob.glob(os.path.join(directory,'*.jpg')):
		files.append(infile)
		# print "current file is " + infile

	result_idx = []
	for i in range(100):
		print i
		img = cv2.imread(files[i],0)
		img = cv2.resize(img, (200,200))
		imf = np.float32(img)
		dct_ = dct(dct(imf.T, norm='ortho').T, norm='ortho')
		# dct_ = np.matrix([[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]])

		# read as dct, return as a list
		list_ = zigzagRead(dct_)
		total_energy = np.linalg.norm(list_)
		print "total enery" + str(total_energy)
		print np.linalg.norm(list_[:1])

		for j in range(1, len(list_)+1):
			if np.linalg.norm(list_[:j]) > 0.99 * total_energy:
				print "J :" + str(j)
				result_idx.append(j)
				break

	return result_idx

def doDCT (filename):
	DCTSlot = 30
	img = cv2.imread(filename,0)
	img = cv2.resize(img, (200,200))
	imf = np.float32(img)
	dct_ = dct(dct(imf.T, norm='ortho').T, norm='ortho')
	# dct_ = np.matrix([[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]])

	# read as dct, return as a list
	list_ = zigzagRead_half(dct_)
	list_ = np.abs(list_)
	# print dct_
	# print list_

	# first 10 slots x, next 10 slots 2x, the third 10 slots 4x
	slotLen = len(list_) / 70 + 1
	# slotLen = len(list_) / DCTSlot + 1
	# print slotLen
	dct_result = np.zeros((DCTSlot,1))

	# get into each slot
	for i in range(len(list_)):
		if i+1 < 10*slotLen:
			dct_result[(i+1) / slotLen] += list_[i]
		elif i+1 >=10*slotLen and i+1 < 30 * slotLen:
			dct_result[10 + (i+1-10*slotLen) / (2*slotLen)] += list_[i]
		else:
			dct_result[20 + (i+1-30*slotLen) / (4*slotLen)] += list_[i]

	# normalization
	base = np.linalg.norm(dct_result)
	dct_result /= base
	dct_result = np.transpose(dct_result)

	return dct_result

def A_dct2(x, n, omega=None):
    """
    Take the 2-dimensional type II DCT of the flattened (n,n) image
    contained in vector x. Works across columns.

    Parameters
    ----------

    x : ndarray, shape (n*n,n_col)
      flattened image vector
    n : int
      image column/row size
    omega : ndarray
      support of the output (ie, indices at which output vector is sampled)

    Returns
    -------

    y = dct2(x.reshape(n,n))[omega]
    """
    col_shapes = x.shape[1:] if len(x.shape) > 1 else ()
    x.shape = (n,n) + col_shapes
    y = fftpack.dct(x, type=2, axis=0, norm='ortho')
    y = fftpack.dct(y, type=2, axis=1, norm='ortho')
    x.shape = (n*n,) + col_shapes
    # this syntax needed because y is discontiguous in memory
    y = np.reshape(y, x.shape)
    if omega:
        return np.take(y, omega, axis=0)
    return y

def dct_distance(rep_one, rep_two, norm=True, num_coefficients=3):
    if not isinstance(rep_one, np.ndarray):
        rep_one = rep_one.to_array()
    if not isinstance(rep_two, np.ndarray):
        rep_two = rep_two.to_array()
    assert (rep_one.shape[1] == rep_two.shape[1])
    num_bands = rep_one.shape[1]
    distance = 0
    for i in range(num_bands):
        try:
            source_dct = dct(rep_one[:, i], norm='ortho')
        except ValueError:
            print(rep_one)
            raise
        if norm:
            source_dct = source_dct[1:]
        source_dct = source_dct[0:num_coefficients]
        target_dct = dct(rep_two[:, i], norm='ortho')
        if norm:
            target_dct = target_dct[1:]
        target_dct = target_dct[0:num_coefficients]
        if len(target_dct) < num_coefficients:
            source_dct = source_dct[:len(target_dct)]
        if len(source_dct) < num_coefficients:
            target_dct = target_dct[:len(source_dct)]
        distance += euclidean(source_dct, target_dct)
    return distance / num_bands

def dctFromPixel(luminance_image, x, y):
	# get the square of size sampling_side centered at (x, y)
	neighbors = find_neighbors(luminance_image, x, y)
	# compute the DCT for the square and reshape it into a 1D array
	discrete_cosine_transform = dct(dct(neighbors.T, norm='ortho').T, norm='ortho')
	feature = np.reshape(discrete_cosine_transform, -1).tolist()
	return feature

def A_dct2(x, n, omega):
    """
    Take the 2-dimensional type II DCT of the flattened (n,n) image
    contained in vector x.

    Parameters
    ----------

    x : ndarray, shape (n*n,)
      flattened image vector
    n : int
      image column/row size
    omega : ndarray
      support of the output (ie, indices at which output vector is sampled)

    Returns
    -------

    y = dct2(x.reshape(n,n))[omega]
    """
    x.shape = (n,n)
    y = fftpack.dct(x, type=2, axis=0, norm='ortho')
    y = fftpack.dct(y, type=2, axis=1, norm='ortho')
    x.shape = (n*n,)
    return y.flat[omega]

def extract_dct_features(time_windows, class_attr=None, n_comps=48):
    X_matrix = []
    y_vect = None
    if class_attr is not None:
        y_vect = []

    for tw in time_windows:        
        x = tw['x'].values
        y = tw['y'].values
        z = tw['z'].values
        m = mag(x,y,z)
        
        dct_x = np.abs(fftpack.dct(x))
        dct_y = np.abs(fftpack.dct(y))
        dct_z = np.abs(fftpack.dct(z))
        dct_m = np.abs(fftpack.dct(m))
        
        v = np.array([])       
        v = np.concatenate((v, dct_x[:n_comps]))            
        v = np.concatenate((v, dct_y[:n_comps]))
        v = np.concatenate((v, dct_z[:n_comps]))
        v = np.concatenate((v, dct_m[:n_comps]))       
        X_matrix.append(v)
        if y_vect is not None:
            y_vect.append(tw[class_attr].iloc[0])       
    
    X_matrix = np.array(X_matrix) 
           
    if y_vect is None:
        return X_matrix
    else:
        return X_matrix, y_vect

def process_cube(img, weight, quality):
    # TODO: check to make sure that size of img, and q_tables are consistent

    img = img.copy()

    # print('process_cube input: {}'.format(img))

    this_quality = np.round(np.max(weight)*quality)

    if this_quality < 0:
        this_quality = 0
    if this_quality > quality - 1:
        this_quality = quality - 1

    for i in range(img.shape[3]):
        img[:, :, :, i] = cv2.cvtColor(img[:, :, :, i], cv2.COLOR_BGR2LAB)
    img = np.float32(img)

    # print('process_cube pre DCT: {}'.format(img))

    # img_dct = dct(dct(dct(img, axis=0)/4, axis=1)/4, axis=3)/4
    img_dct = dct(dct(img, axis=0)/4, axis=1)/4

    Q_luma = luminance_tables[:, :, :, this_quality].astype(np.float32)
    Q_chroma = chrominance_tables[:, :, :, this_quality].astype(np.float32)

    # Q_luma[:, :, :] = .01
    # Q_chroma[:, :, :] = .01

    # print('Q_luma: {}'.format(Q_luma))
    # print('Q_chroma: {}'.format(Q_chroma))

    # print('dct, pre rounding: {}'.format(img_dct))
    img_dct[:, :, 0, :] /= Q_luma
    img_dct[:, :, 1, :] /= Q_chroma
    img_dct[:, :, 2, :] /= Q_chroma

    img_dct = np.round(img_dct)

    img_dct[:, :, 0, :] *= Q_luma
    img_dct[:, :, 1, :] *= Q_chroma
    img_dct[:, :, 2, :] *= Q_chroma
    # print('dct, post rounding: {}'.format(img_dct))

    # img_processed = idct(idct(idct(img_dct, axis=0)/4, axis=1)/4, axis=3)/4
    img_processed = idct(idct(img_dct, axis=0)/4, axis=1)/4

    # print('process_cube post DCT: {}'.format(img_processed))

    img_processed = np.clip(img_processed, 0, 255)
    img_processed = np.uint8(img_processed)

    for i in range(img.shape[3]):
        img_processed[:,:,:,i] = cv2.cvtColor(img_processed[:,:,:,i], cv2.COLOR_LAB2BGR)

    # print('process_cube output: {}'.format(img))

    # print('pre dct / post_dct: {}'.format(pre_dct / post_dct))
    return img_processed

def discrete_cosine(reference, target):
    
    reference_transform = fftpack.dct( reference )
    target_transform = fftpack.dct( target )
    reference_curve = reference_transform.mean(axis = 1)
    target_curve = target_transform.mean(axis = 1)

    return reference_transform, target_transform, reference_curve, target_curve

 def fastChebScalar(self, fj, fk):
     """Fast Chebyshev scalar product."""
     N = fj.shape[0]
     if self.quad == "GL":
         fk = dct(fj, 2, axis=0)*np.pi/(2*N)
     elif self.quad == "GC":
         fk = dct(fj, 1, axis=0)*np.pi/(2*(N-1))
     return fk

def dct2(arr):
    """
    @todo: 2D-dct, constructed from scipy's dct2
    @params { np.ndarray } arr
    @return { np.ndarray }
    """
    array = np.float64(arr)
    result = dct(dct(array, axis=0), axis=1)
    return result

def laplacian_pca_TV(res, x, f0, lam, gam, iter = 10):
    '''
    TV version of Laplacian embedding
    :param res: resolution of the grid
    :param x: numpy array of data in rows
    :param f0: initial embedding matrix
    :param lam: sparsity parameter
    :param gam: fidelity parameter
    :param iter: number of iterations to carry out
    :return: returns embedding matrix
    '''
    # f0 is an initial projection
    n = res ** 2
    num_data = x.shape[0]

    D = sparse_discrete_diff(res)
    M = 1/(lam*laplacian_eigenvalues(res).reshape(n)+gam)

    f = f0
    y = x .dot(f)
    z = shrink(y .dot(D.T), lam)

    for i in range(iter):

        # Update z
        z_old = z
        z = shrink(y .dot (D.T), lam)

        # Update f
        f_old = f
        u, s, v = la.svd(x.T .dot (y), full_matrices=False)
        f = u .dot(v)

        # Update y
        y_old = y
        q = lam * z .dot (D) + gam * x .dot(f)
        # print('norm of y before is %f' % np.sum(q ** 2))
        y = fftpack.dct(q.reshape((num_data, res, res)), norm='ortho') # Images unraveled as rows
        y = fftpack.dct(np.swapaxes(y,1,2), norm='ortho') # Swap and apply dct on the other side
        # print('norm of y after is %f' % np.sum(y ** 2))
        y = np.apply_along_axis(lambda v: M * v, 1, y.reshape((num_data, n)))
        y = fftpack.idct(y.reshape((num_data, res, res)), norm='ortho')
        y = fftpack.idct(np.swapaxes(y,1,2), norm='ortho')
        y = y.reshape((num_data, n))

        zres = np.sqrt(np.sum((z - z_old) ** 2))
        znorm = np.sqrt(np.sum((z)**2))
        yres = np.sqrt(np.sum((y - y_old) ** 2))
        ynorm = np.sqrt(np.sum((y)**2))
        fres = np.sqrt(np.sum((f - f_old) ** 2))

        value = np.sum(abs(z)) + 0.5*lam*np.sum((z-y .dot(D.T))**2) + 0.5*gam*np.sum((y- x .dot(f) )**2)
        print('Iter %d Val %f Z norm %f Z res %f Ynorm %f Y res %f F res %f' % (i, value, znorm, zres, ynorm, yres, fres))

    return f

 def ifct(self, fk, cj):
     """Inverse fast Chebyshev transform."""
     if self.quad == "GL":
         cj = 0.5*dct(fk, 3, axis=0)
         cj += 0.5*fk[0]
     elif self.quad == "GC":
         cj = 0.5*dct(fk, 1, axis=0)
         cj += 0.5*fk[0]
         cj[::2] += 0.5*fk[-1]
         cj[1::2] -= 0.5*fk[-1]
     return cj

 def fct(self, fj, cj):
     """Fast Chebyshev transform."""
     N = fj.shape[0]
     if self.quad == "GL":
         cj = dct(fj, 2, axis=0)
         cj /= N
         cj[0] /= 2        
     elif self.quad == "GC":
         cj = dct(fj, 1, axis=0)/(N-1)
         cj[0] /= 2
         cj[-1] /= 2
     return cj

def At_dct2(y, n, omega):
    """
    Take the 2-dimensional type III DCT of the flattened (n,n) matrix
    contained in vector y. This is the adjoint operator to the A_dct2
    operator defined above
    """

    y2 = np.zeros((n,n), 'd')
    y2.flat[omega] = y
    w = fftpack.dct(y2, type=3, axis=0, norm='ortho')
    w = fftpack.dct(w, type=3, axis=1, norm='ortho')
    return w.flatten()

def dct2(y):
    M = y.shape[0]
    N = y.shape[1]
    a = empty([M,N],float)
    b = empty([M,N],float)

    for i in range(M):
        a[i,:] = dct(y[i,:],norm='ortho')
    for j in range(N):
        b[:,j] = dct(a[:,j],norm='ortho')

    return b

def calculate_DCTII_2D(matrix):
    """
    Calculates the 2D transform of the DCT II algorithm.
    Assumes a square matrix.

    See:
        http://en.wikipedia.org/wiki/Discrete_cosine_transform#DCT-II

    We are using the plain version, which seems to work better.
    """
    a = numpy.reshape(numpy.array(matrix), (32, 32))
    return fftpack.dct(fftpack.dct(a.T).T)

def multiwavelet_from_rgb(rgb):
    from scipy.fftpack import dct
    from pywt import wavedec2

    r = rgb[:, :, 0].astype(np.float)
    g = rgb[:, :, 1].astype(np.float)

    dctr = dct(r, norm='ortho').ravel()
    dctg = dct(g, norm='ortho').ravel()
    daubr = _unpack(wavedec2(r, 'db4'))
    daubg = _unpack(wavedec2(g, 'db4'))
    return np.hstack([dctr, dctg, daubr, daubg])