def objective_gradient(self, X, J=None, return_K=False):
     """
     Computes MPF objective gradient on input data X given coupling
     strengths J.
     
     Parameters
     ----------
     X : numpy array
         (M, N)-dim array of binary input patterns of length N,
         where N is the number of nodes in the network
     J : numpy array, optional
         Coupling matrix of size N x N, where N denotes the number
         of nodes in the network (default None)
     return_K : bool, optional
         Flag wether to return K (default False)
     
     Returns
     -------
     dJ [, K] : numpy array [, numpy array]
         Update to coupling matrix J [and K if return_K is True]
     """
     if J is None:
         J = self._J
         J[np.eye(self._N, dtype=bool)] = -2 * self._theta
     X = np.atleast_2d(X)
     M, N = X.shape
     S = 2 * X - 1
     Kfull = np.exp(-S * np.dot(X, J.T) + .5 * np.diag(J)[None, :])
     dJ = -np.dot(X.T, Kfull * S) + .5 * np.diag(Kfull.sum(0))
     if self._symmetric is True:
         dJ = .5 * (dJ + dJ.T)
     if return_K:
         return Kfull.sum() / M, dJ / M
     else:
         return dJ / M

def phiHJC(eGAF, eGFA, kA):
    """
    Calculate initial HJC vector for openings by solving
    phi*(I-eGAF*eGFA)=0 (Eq. 10, HJC92)
    For shuttings exhange A by F and F by A in function call.

    Parameters
    ----------
    eGAF : array_like, shape (kA, kF)
    eGFA : array_like, shape (kF, kA)
    kA : int
        A number of open states in kinetic scheme.
    kF : int
        A number of shut states in kinetic scheme.

    Returns
    -------
    phi : array_like, shape (kA)
    """

    if kA == 1:
        phi = np.array([1])

    else:
        Qsub = np.eye(kA) - np.dot(eGAF, eGFA)
        u = np.ones((kA, 1))
        S = np.concatenate((Qsub, u), 1)
        phi = np.dot(u.transpose(), nplin.inv(np.dot(S, S.transpose())))[0]

    return phi

def CHSvec(roots, tres, tcrit, QFA, kA, expQAA, phiF, R):
    """
    Calculate initial and final CHS vectors for HJC likelihood function
    (Eqs. 5.5 or 5.7, CHS96).

    Parameters
    ----------
    roots : array_like, shape (1, kA)
        Roots of the asymptotic pdf.
    tres : float
        Time resolution (dead time).
    tcrit : float
        Critical time.
    QFA : array_like, shape(kF, kA)
    kA : int
    expQAA : array_like, shape(kA, kA)
    phiF : array_like, shape(1, kF)
    R : array_like, shape(kF, kF, kF)

    Returns
    -------
    start : ndarray, shape (1, kA)
        CHS start vector (Eq. 5.11, CHS96).
    end : ndarray, shape (kF, 1)
        CHS end vector (Eq. 5.8, CHS96).
    """

    H = HAF(roots, tres, tcrit, QFA, expQAA, R)
    u = np.ones((kA, 1))
    start = np.dot(phiF, H) / np.dot(np.dot(phiF, H), u)
    end = np.dot(H, u)

    return start, end

def read_abinit(filename):
    with open(filename) as f:
        abinit_in = AbinitIn(f.readlines())
    tags = abinit_in.get_variables()
    acell = tags['acell']
    rprim = tags['rprim'].T
    scalecart = tags['scalecart']
    lattice = rprim * acell
    if scalecart is not None:
        for i in range(3):
            lattice[i] *= scalecart[i]

    if tags['xcart'] is not None:
        pos_bohr = np.transpose(tags['xcart'])
        positions = np.dot(np.linalg.inv(lattice), pos_bohr).T
    elif tags['xangst'] is not None:
        pos_bohr = np.transpose(tags['xangst']) / Bohr
        positions = np.dot(np.linalg.inv(lattice), pos_bohr).T
    elif tags['xred'] is not None:
        positions = tags['xred']

    numbers = [tags['znucl'][x - 1] for x in tags['typat']]

    return Atoms(numbers=numbers,
                 cell=lattice.T,
                 scaled_positions=positions)

def eGs(GAF, GFA, kA, kF, expQFF):
    """
    Calculate eGAF, probabilities from transitions from apparently open to
    shut states regardles of when the transition occurs. Thease are Laplace
    transform of eGAF(t) when s=0. Used to calculat initial HJC vectors (HJC92).
    eGAF*(s=0) = (I - GAF * (I - expQFF) * GFA)^-1 * GAF * expQFF
    To caculate eGFA exhange A by F and F by A in function call.

    Parameters
    ----------
    GAF : array_like, shape (kA, kF)
    GFA : array_like, shape (kF, kA)
    kA : int
        A number of open states in kinetic scheme.
    kF : int
        A number of shut states in kinetic scheme.
    
    Returns
    -------
    eGAF : array_like, shape (kA, kF)
    """

    temp = np.eye(kA) - np.dot(np.dot(GAF, np.eye(kF) - expQFF), GFA)
    eGAF = np.dot(np.dot(nplin.inv(temp), GAF), expQFF)
    return eGAF

    def fitToData(self, data):
        '''
        param data: numpy array where [:,0] is x and [:,1] is y
        '''
        x = data[:, 0][:, np.newaxis]
        y = data[:, 1][:, np.newaxis]
        D = np.hstack((x*x, x*y, y*y, x, y, np.ones_like(x)))
        S = np.dot(D.T, D)
        C = np.zeros([6, 6])
        C[0, 2] = C[2, 0] = 2; C[1, 1] = -1
        E, V = eig(np.dot(inv(S), C))
        n = np.argmax(np.abs(E))
        self.parameters = V[:, n]

        axes = self.ellipse_axis_length()
        self.a = axes[0]
        self.b = axes[1]
        self.angle = self.ellipse_angle_of_rotation()

        if not self.a or not self.b or self.parameters == None or np.iscomplexobj(self.parameters) or \
           math.isnan(self.a) or math.isnan(self.b) or math.isnan(self.ellipse_center()[0]) or \
           np.iscomplex(self.ellipse_center()[0]) or np.iscomplex(self.a) or np.iscomplex(self.b) or \
           np.iscomplexobj(self.angle):
            self.a = 0
            self.b = 0
            self.parameters = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
            self.angle = 0
            self.error = True

def test_grid_search_precomputed_kernel():
    """Test that grid search works when the input features are given in the
    form of a precomputed kernel matrix """
    X_, y_ = make_classification(n_samples=200, n_features=100, random_state=0)

    # compute the training kernel matrix corresponding to the linear kernel
    K_train = np.dot(X_[:180], X_[:180].T)
    y_train = y_[:180]

    clf = SVC(kernel='precomputed')
    cv = GridSearchCV(clf, {'C': [0.1, 1.0]})
    cv.fit(K_train, y_train)

    assert_true(cv.best_score_ >= 0)

    # compute the test kernel matrix
    K_test = np.dot(X_[180:], X_[:180].T)
    y_test = y_[180:]

    y_pred = cv.predict(K_test)

    assert_true(np.mean(y_pred == y_test) >= 0)

    # test error is raised when the precomputed kernel is not array-like
    # or sparse
    assert_raises(ValueError, cv.fit, K_train.tolist(), y_train)

    def smooth_objective(self, x, mode='both', check_feasibility=False):
        """
        Evaluate a smooth function and/or its gradient

        if mode == 'both', return both function value and gradient
        if mode == 'grad', return only the gradient
        if mode == 'func', return only the function value
        """
        x = self.apply_offset(x)
        exp_x = np.exp(x)

        #TODO: Using transposes to scale the rows of a 2d array - should we use an affine_transform to do this?
        #JT: should be able to do this with np.newaxis

        if mode == 'both':
            ratio = ((self.trials/(1. + np.sum(exp_x, axis=1))) * exp_x.T).T
            f, g = -2. * self.scale(np.sum(self.firstcounts * x) -  np.dot(self.trials, np.log(1. + np.sum(exp_x, axis=1)))), - 2 * self.scale(self.firstcounts - ratio) 
            return f, g
        elif mode == 'grad':
            ratio = ((self.trials/(1. + np.sum(exp_x, axis=1))) * exp_x.T).T
            f, g = None, - 2 * self.scale(self.firstcounts - ratio) 
            return g
        elif mode == 'func':
            f, g = -2. * self.scale(np.sum(self.firstcounts * x) -  np.dot(self.trials, np.log(1. + np.sum(exp_x, axis=1)))), None
            return f
        else:
            raise ValueError("mode incorrectly specified")

def test_dot():
    # Test normal dot
    a = np.random.uniform(-3, 3, (3, 4))
    b = np.random.uniform(-3, 3, (4, 5))
    c = np.dot(a, b)
    A = mx.nd.array(a)
    B = mx.nd.array(b)
    C = mx.nd.dot(A, B)
    assert reldiff(c, C.asnumpy()) < 1e-5
    # Test dot with transpose kargs
    a = np.random.uniform(-3, 3, (3, 4))
    b = np.random.uniform(-3, 3, (3, 5))
    c = np.dot(a.T, b)
    A = mx.nd.array(a)
    B = mx.nd.array(b)
    C = mx.nd.dot(A, B, transpose_a=True)
    assert reldiff(c, C.asnumpy()) < 1e-5
    # Test dot with transpose kargs
    a = np.random.uniform(-3, 3, (3, 4))
    b = np.random.uniform(-3, 3, (5, 4))
    c = np.dot(a, b.T)
    A = mx.nd.array(a)
    B = mx.nd.array(b)
    C = mx.nd.dot(A, B, transpose_b=True)
    assert reldiff(c, C.asnumpy()) < 1e-5
    # Test dot with transpose kargs
    a = np.random.uniform(-3, 3, (4, 3))
    b = np.random.uniform(-3, 3, (5, 4))
    c = np.dot(a.T, b.T)
    A = mx.nd.array(a)
    B = mx.nd.array(b)
    C = mx.nd.dot(A, B, transpose_a=True, transpose_b=True)
    assert reldiff(c, C.asnumpy()) < 1e-5

    def InitialAlignment(self, scale = 0.15):
        """ Compute SVD and align object to be in a certain coordinate frame.
        
        Usage: model.InitialAlignment(scale)

        Input:
            scale - Desired scale for object. Scale is defined as the length
            along the leading eigenvector, in meters.
        """


        pts3D = self.pts3D

        # Compute eigenvecs and rotate according to them
        pc, evals, mean = utils.pca(pts3D, remove_mean = True)
        pts3D_rot = np.dot(pc.T, pts3D)

        # Find length according to max eigenvector
        mins = np.min(pts3D_rot, axis=1)
        maxs = np.max(pts3D_rot, axis=1)
        max_length = maxs[0] - mins[0]
        
        # Rotation matrix is the covariance matrix, but we want Z as the leading
        # eigenvector:
        rot = np.c_[-pc[2], pc[1], pc[0]]

        # Transform model to have zero mean, reasonable scale and rotation.
        self.transform(rot, np.dot(rot, -mean), float(scale) / max_length)

    def get_corr_pred( self, sctx, u, du, tn, tn1,
                      u_avg = None,
                      B_mtx_grid = None,
                      J_det_grid = None,
                      ip_coords = None,
                      ip_weights = None ):
        '''
        Corrector and predictor evaluation.

        @param u current element displacement vector
        '''
        if J_det_grid == None or B_mtx_grid == None:
            X_mtx = sctx.X

        show_comparison = True
        if ip_coords == None:
            ip_coords = self.ip_coords
            show_comparison = False
        if ip_weights == None:
            ip_weights = self.ip_weights

        ### Use for Jacobi Transformation

        n_e_dofs = self.n_e_dofs
        K = zeros( ( n_e_dofs, n_e_dofs ) )
        F = zeros( n_e_dofs )
        sctx.fets_eval = self
        ip = 0

        for r_pnt, wt in zip( ip_coords, ip_weights ):
            #r_pnt = gp[0]
            sctx.r_pnt = r_pnt
#caching cannot be switched off in the moment
#            if J_det_grid == None:
#                J_det = self._get_J_det( r_pnt, X_mtx )
#            else:
#                J_det = J_det_grid[ip, ... ]
#            if B_mtx_grid == None:
#                B_mtx = self.get_B_mtx( r_pnt, X_mtx )
#            else:
#                B_mtx = B_mtx_grid[ip, ... ]
            J_det = J_det_grid[ip, ... ]
            B_mtx = B_mtx_grid[ip, ... ]

            eps_mtx = dot( B_mtx, u )
            d_eps_mtx = dot( B_mtx, du )
            sctx.mats_state_array = sctx.elem_state_array[ip * self.m_arr_size: ( ip + 1 ) * self.m_arr_size]
            #print 'elem state ', sctx.elem_state_array
            #print 'mats state ', sctx.mats_state_array
            sctx.r_ls = sctx.ls_val[ip]
            sig_mtx, D_mtx = self.get_mtrl_corr_pred( sctx, eps_mtx, d_eps_mtx, tn, tn1 )
            k = dot( B_mtx.T, dot( D_mtx, B_mtx ) )
            k *= ( wt * J_det )
            K += k
            f = dot( B_mtx.T, sig_mtx )
            f *= ( wt * J_det )
            F += f
            ip += 1

        return F, K

 def train(self, inp, out, training_weight=1.):
     inp = np.mat(inp).T
     out = np.mat(out).T
     deriv = []
     val = inp
     vals = [val]
     # forward calculation of activations and derivatives
     for weight,bias in self.__weights:
         val = weight*val
         val += bias
         deriv.append(self.__derivative(val))
         vals.append(self.__activation(val))
     deriv = iter(reversed(deriv))
     weights = iter(reversed(self.__weights))
     errs = []
     errs.append(np.multiply(vals[-1]-out, next(deriv)))
     # backwards propagation of errors
     for (w,b),d in zip(weights, deriv):
         errs.append(np.multiply(np.dot(w.T, errs[-1]), d))
     weights = iter(self.__weights)
     for (w,b),v,e in zip(\
             self.__weights,\
             vals, reversed(errs)):
         e *= self.__learning_rate*training_weight
         w -= e*v.T
         b -= e
     tmp = vals[-1]-out
     return np.dot(tmp[0].T,tmp[0])*.5*training_weight

def test_grtm():
    l = language(1000)
    n_iter = 1000
    KL_thresh = 0.3

    mu = 0.
    nu2 = 1.
    np.random.seed(l['seed'])
    H = np.random.normal(loc=mu, scale=nu2, size=(l['K'], l['K']))
    zeta = pd.DataFrame([(i, j, np.dot(np.dot(l['thetas'][i], H),
                                       l['thetas'][j]))
                         for i, j in product(range(l['D']), repeat=2)],
                        columns=('tail', 'head', 'zeta'))
    zeta['y'] = (zeta.zeta >= 0).astype(int)
    y = zeta[['tail', 'head', 'y']].values
    skf = StratifiedKFold(y[:, 2], n_folds=100)
    _, train_idx = next(iter(skf))
    _K = l['K']
    _alpha = l['alpha'][:_K]
    _beta = np.repeat(0.01, l['V'])
    _b = 1.
    grtm = GRTM(_K, _alpha, _beta, mu, nu2, _b, n_iter, seed=l['seed'],
                n_report_iter=l['n_report_iters'])
    grtm.fit(l['doc_term_matrix'], y[train_idx])

    assert_probablity_distribution(grtm.phi)
    check_KL_divergence(l['topics'], grtm.phi, KL_thresh)

def EvaluatePolicy(s, w_pi, useRBFKernel = False):
  
    # the value of the improved policy
    value = np.zeros((len(s),1))

    # the new policy
    policy = [False] * len(s)

    # iterate through every state, 
    for idx in range(len(s)):

        # State-Action value function for actions 0.0 and 1.0
        if useRBFKernel == True:
            q0 = np.dot(computePhiRBF(s[idx], 0.0).T, w_pi)
            q1 = np.dot(computePhiRBF(s[idx], 1.0).T, w_pi)
        else:
            q0 = np.dot(np.append(s[idx, 0],0.0), w_pi)
            q1 = np.dot(np.append(s[idx, 0],1.0), w_pi)

        # update the value
        value[idx] = max(q0, q1)

        # update the policy
        policy[idx] = True if q1 > q0 else False
        
    return (policy, value)

 def forward(self, X):
     #Propogate inputs though network
     self.z2 = np.dot(X, self.W1)
     self.a2 = self.sigmoid(self.z2)
     self.z3 = np.dot(self.a2, self.W2)
     yHat = self.sigmoid(self.z3) 
     return yHat

	def updateParameters(self, articlePicked, click,  userID):	
		self.counter +=1
		self.Wlong = vectorize(self.W)
		featureDimension = len(articlePicked.featureVector)
		T_X = vectorize(np.outer(articlePicked.featureVector, self.W.T[userID])) 
		self.A += np.outer(T_X, T_X)	
		self.b += click*T_X
		self.AInv = np.linalg.inv(self.A)
		self.UserTheta = matrixize(np.dot(self.AInv, self.b), len(articlePicked.featureVector)) 

		Xi_Matirx = np.zeros(shape = (featureDimension, self.userNum))
		Xi_Matirx.T[userID] = articlePicked.featureVector
		W_X = vectorize( np.dot(np.transpose(self.UserTheta), Xi_Matirx))
		self.batchGradient +=evaluateGradient(W_X, click, self.Wlong, self.lambda_, self.regu  )

		if self.counter%self.windowSize ==0:
			self.Wlong -= 1/(float(self.counter/self.windowSize)+1)*self.batchGradient
			self.W = matrixize(self.Wlong, self.userNum)
			self.W = normalize(self.W, axis=0, norm='l1')
			#print 'SVD', self.W
			self.batchGradient = np.zeros(self.userNum*self.userNum)
			# Use Ridge regression to fit W
		'''
		plt.pcolor(self.W_b)
		plt.colorbar
		plt.show()
		'''
		if self.W.T[userID].any() <0 or self.W.T[userID].any()>1:
			print self.W.T[userID]

		self.CoTheta = np.dot(self.UserTheta, self.W)
		self.BigW = np.kron(np.transpose(self.W), np.identity(n=len(articlePicked.featureVector)))
		self.CCA = np.dot(np.dot(self.BigW , self.AInv), np.transpose(self.BigW))
		self.BigTheta = np.kron(np.identity(n=self.userNum) , self.UserTheta)

def ehist_equalize_melhist(d, sr, refMelHist, edges):
    """ Modify a signal in the Mel domain 
        by equalizing the Mel-subband histograms to match
        the passed-in ones """
    # Calculate the (Mel) spectrograms, and histogram, and axes
    melHist, edges, D, DmeldB, melmx, freqs = mel_hist(d, sr, edges=edges)
    # Build mapping & modify mel spectrogram
    histmaps = make_hist_maps(melHist, refMelHist, edges)
    # for some reason, extrapolating madly below bottom edge - clip it
    DmeldBmapped = np.maximum(edges[0], 
                              np.minimum(edges[-1], 
                                         apply_hist_maps(DmeldB, histmaps)))
    # Reconstruct audio based on mapped envelope
    # We map both original and modified Mel envelopes to FFT domain
    # then scale original STFT magnitudes by their ratio
    DmelInFFT = np.dot(melmx.T, idB(DmeldB))
    DmappedInFFT = np.dot(melmx.T, idB(DmeldBmapped))
    # Zero values in denominator will match to zeros in numerator, 
    # so it's OK to drop them
    Dmask = DmappedInFFT / (DmelInFFT + (DmelInFFT==0))
    # Median filter to remove short blips in gain
    medfiltwin = 7
    DmaskF = median_filter(Dmask, size=(1, medfiltwin))
    # Now scale each FFT val by their ratio
    Dmod = D * DmaskF
    # and resynthesize
    nfft = 2*(np.size(D, axis=0)-1)
    win = nfft
    hop = win/4
    dout = istft(Dmod.T, win, hop)
    return dout

def vertex_transform1(vertex):
    """
    This transform was applied on the original surface.
    """
    return np.dot(rotation_matrix(np.array([0.0, 0.0, 1.0]), math.pi),
                  np.dot(rotation_matrix(np.array([1.0, 0.0, 0.0]), -math.pi / 1.6),
                         np.array([float(x) / 1.5 for x in vertex[:3]]) + np.array([0.0, -40.0, 20.0])))

def runLabeling(file_path, gps_filename, output_name, frames_to_skip, final_frame, lp, rp, pickle_loc):
    video_reader = WarpedVideoReader(file_path)
    #video_reader.setSubsample(True)
    video_reader.setPerspectives(pickle_loc)
    gps_reader = GPSReader(gps_filename)
    gps_dat = gps_reader.getNumericData()

    cam = getCameraParams()
    cam_to_use = cam[int(output_name[-1]) - 1]

    lp = pixelTo3d(lp, cam_to_use)
    rp = pixelTo3d(rp, cam_to_use)
    tr = GPSTransforms(gps_dat, cam_to_use)
    pitch = -cam_to_use['rot_x']
    height = 1.106
    R_camera_pitch = euler_matrix(cam_to_use['rot_x'], cam_to_use['rot_y'], cam_to_use['rot_z'], 'sxyz')[0:3, 0:3]
    Tc = np.eye(4)
    Tc[0:3, 0:3] = R_camera_pitch.transpose()
    Tc[0:3, 3] = [-0.2, -height, -0.5]
    lpts = np.zeros((lp.shape[0], 4))
    rpts = np.zeros((rp.shape[0], 4))
    for t in range(min(tr.shape[0], lp.shape[0])):
        lpts[t, :] = np.dot(tr[t, :, :], np.linalg.solve(Tc, np.array([lp[t, 0], lp[t, 1], lp[t, 2], 1])))
        rpts[t, :] = np.dot(tr[t, :, :], np.linalg.solve(Tc, np.array([rp[t, 0], rp[t, 1], rp[t, 2], 1])))

    ldist = np.apply_along_axis(np.linalg.norm, 1, np.concatenate((np.array([[0, 0, 0, 0]]), lpts[1:] - lpts[0:-1])))
    rdist = np.apply_along_axis(np.linalg.norm, 1, np.concatenate((np.array([[0, 0, 0, 0]]), rpts[1:] - rpts[0:-1])))
    start_frame = frames_to_skip
    runBatch(video_reader, gps_dat, cam_to_use, output_name, start_frame, final_frame, lpts, rpts, ldist, rdist, tr)

    print "Done with %s" % output_name

    def backprop(self, x, y):
        activation = x
        activations = [x]
        zs = []

        for weight, bias in zip(self.weights, self.biases):
            z = np.dot(activation, weight)+bias
            zs.append(z)
            activation = sigmoid(z)
            activations.append(activation)

        delta = (activation-y)*sigmoid_prime(zs[-1])

        nabla_weights = [np.zeros(w.shape) for w in self.weights]
        nabla_biases = [np.zeros(b.shape) for b in self.biases]

        nabla_weights[-1] = np.dot(activations[-2].transpose(), delta)
        nabla_biases[-1] = delta

        for l in xrange(2, len(self.layers)):
            delta = np.dot(delta, self.weights[-l+1].transpose())*sigmoid_prime(zs[-l])
            nabla_weights[-l] = np.dot(activations[-l-1].transpose(), delta)
            nabla_biases[-l] = delta

        return (nabla_weights, nabla_biases)