def objective(pars,Z,ycovar):
    """The objective function"""
    bcost= pars[0]
    t= pars[1]
    Pb= pars[2]
    Xb= pars[3]
    Yb= pars[4]
    Zb= sc.array([Xb,Yb])
    Vb1= sc.exp(pars[5])
    Vb2= sc.exp(pars[6])
    corr= pars[7]
    V= sc.array([[Vb1,sc.sqrt(Vb1*Vb2)*corr],[sc.sqrt(Vb1*Vb2)*corr,Vb2]])
    v= sc.array([-sc.sin(t),sc.cos(t)])
    if Pb < 0. or Pb > 1.:
        return -sc.finfo(sc.dtype(sc.float64)).max
    if corr < -1. or corr > 1.:
        return -sc.finfo(sc.dtype(sc.float64)).max
    delta= sc.dot(v,Z.T)-bcost
    sigma2= sc.dot(v,sc.dot(ycovar,v))

    ndata= Z.shape[0]
    detVycovar= sc.zeros(ndata)
    deltaOUT= sc.zeros(ndata)
    for ii in range(ndata):
        detVycovar[ii]= m.sqrt(linalg.det(V+ycovar[:,ii,:]))
        deltaOUT[ii]= sc.dot(Z[ii,:]-Zb,sc.dot(linalg.inv(V+ycovar[:,ii,:]),Z[ii,:]-Zb))
    return sc.sum(sc.log((1.-Pb)/sc.sqrt(2.*m.pi*sigma2)*
                         sc.exp(-0.5*delta**2./sigma2)
                         +Pb/2./m.pi/detVycovar
                         *sc.exp(-0.5*deltaOUT)))

def rotation_matrix_from_cross_prod(a,b):
    """
    Returns the rotation matrix which rotates the
    vector :samp:`a` onto the the vector :samp:`b`.
    
    :type a: 3 sequence of :obj:`float`
    :param a: Vector to be rotated on to :samp:`{b}`.
    :type b: 3 sequence of :obj:`float`
    :param b: Vector.
    :rtype: :obj:`numpy.array`
    :return: 3D rotation matrix.
    """

    crs = np.cross(a,b)
    dotProd = np.dot(a,b)
    crsNorm = sp.linalg.norm(crs)
    eps = sp.sqrt(sp.finfo(a.dtype).eps)
    r = sp.eye(a.size, a.size, dtype=a.dtype)
    if (crsNorm > eps):
        theta = sp.arctan2(crsNorm, dotProd)
        r = axis_angle_to_rotation_matrix(crs, theta)
    elif (dotProd < 0):
        r = -r
    
    return r

def pre_compute_E_Beta(x, y, sig, kernel="RBF"):
    """
    Function that pre computes the kernel eigenvalues/eigenfunctions during the cross-validation
    Input:
    x,y: the sample matrix and the label
    sig: the value of the kernel parameters
    Output:
    E_: a list of eigenvalues
    Beta_: a list of corresponding eigenvectors
    """
    C = int(y.max())
    eps = sp.finfo(sp.float64).eps
    E_ = []
    Beta_ = []

    for i in range(C):
        t = sp.where(y == (i + 1))[0]
        ni = t.size
        Ki = KERNEL()
        Ki.compute_kernel(x[t, :], kernel=kernel, sig=sig)
        Ki.center_kernel()
        Ki.scale_kernel(ni)

        E, Beta = linalg.eigh(Ki.K)
        idx = E.argsort()[::-1]
        E = E[idx]
        E[E < eps] = eps
        Beta = Beta[:, idx]
        E_.append(E)
        Beta_.append(Beta)

    del E, Beta, Ki
    return E_, Beta_

def edge_property_vs_depth(G,property,intervals,eIndices=None,function=None):
    """Generic function to compile and optionally process edge information of a
    vascular graph versus the cortical depth.
    INPUT: G: Vascular graph in iGraph format.
           property: Which edge property to operate on.
           intervals: Intervals of cortical depth in which the sample is split.
                      (Expected 
           eIndices: (Optional.) Indices of edges to consider. If not provided,
                     all edges are taken into account.
           function: (Optional.) Function which to perform on the compiled data
                     of each interval.
    OUTPUT: The compiled (and possibly processed) information as a list (one 
            entry per interval).
    """
    
    intervals[-1] = (intervals[-1][0], intervals[-1][1] + sp.finfo(float).eps)
    database = []
    for interval in intervals:
        if eIndices:
            data = G.es(eIndices,depth_ge=interval[0], 
                        depth_lt=interval[1])[property]
        else:                
            data = G.es(depth_ge=interval[0], 
                        depth_lt=interval[1])[property]
        if function:
            database.append(function(data))
        else:
            database.append(data)
    return database   

def implant_srxtm_cube(Ga, Gd, origin=None, crWidth=150):
    """Implants a cubic srXTM sample in an artificial vascular network
    consisting of pial vessels, arterioles, venoles and capillaries. At the 
    site of insertion, artificial vessels are removed to make space for the 
    srXTM sample. At the border between artificial and data networks,
    connections are made between the respective loose ends.
    INPUT: Ga: VascularGraph of the artificial vasculature.
           Gd: VascularGraph of the srXTM data. The graph is expected to have
               the attributes 'av' and 'vv' that denote the indices of 
               endpoints of the large penetrating arteries and veins 
               respectively.
           origin: The two-dimensional origin in the xy plane, where the srXTM
                   sample should be inserted. This will be the network's 
                   center of mass, if not provided.
           crWidth: Width of connection region. After computing the center and
                    radius of the SRXTM sample, loose ends of the SRXTM at
                    radius - 2*crWidth are connected to loose ends of the
                    artificial network at radius + 2*crWidth.
    OUTPUT: Ga: Modified input Ga, with Gd inserted at origin.       
    """
    
    # Create Physiology object with appropriate default units:
    P = vgm.Physiology(Ga['defaultUnits'])
    eps = finfo(float).eps * 1e4


    return Ga

    def train(self, x, y, mu=None, sig=None):
        # Initialization
        n = y.shape[0]
        C = int(y.max())
        eps = sp.finfo(sp.float64).eps

        if (mu is None) and (self.mu is None):
            mu = 10 ** (-7)
        elif self.mu is None:
            self.mu = mu

        if (sig is None) and (self.sig is None):
            self.sig = 0.5
        elif self.sig is None:
            self.sig = sig

        # Compute K and
        K = KERNEL()
        K.compute_kernel(x, sig=self.sig)
        G = KERNEL()
        G.K = self.mu * sp.eye(n)

        for i in range(C):
            t = sp.where(y == (i + 1))[0]
            self.ni.append(sp.size(t))
            self.prop.append(float(self.ni[i]) / n)

            # Compute K_k
            Ki = KERNEL()
            Ki.compute_kernel(x, z=x[t, :], sig=self.sig)
            T = sp.eye(self.ni[i]) - sp.ones((self.ni[i], self.ni[i]))
            Ki.K = sp.dot(Ki.K, T)
            del T
            G.K += sp.dot(Ki.K, Ki.K.T) / self.ni[i]
        G.scale_kernel(C)

        # Solve the generalized eigenvalue problem
        a, A = linalg.eigh(G.K, b=K.K)
        idx = a.argsort()[::-1]
        a = a[idx]
        A = A[:, idx]

        # Remove negative eigenvalue
        t = sp.where(a > eps)[0]
        a = a[t]
        A = A[:, t]

        # Normalize the eigenvalue
        for i in range(a.size):
            A[:, i] /= sp.sqrt(sp.dot(sp.dot(A[:, i].T, K.K), A[:, i]))

        # Update model
        self.a = a.copy()
        self.A = A.copy()
        self.S = sp.dot(sp.dot(self.A, sp.diag(self.a ** (-1))), self.A.T)

        # Free memory
        del G, K, a, A

def objective(pars,X,Y,yerr):
    """The objective function"""
    b= pars[0]
    s= pars[1]
    Pb= pars[2]
    Yb= pars[3]
    Vb= m.exp(pars[4])
    if Pb < 0. or Pb > 1.:
        return -sc.finfo(sc.dtype(sc.float64)).max
    return sc.sum(sc.log((1.-Pb)/sc.sqrt(2.*m.pi)/yerr*sc.exp(-0.5*(Y-s*X-b)**2./yerr**2.)+Pb/sc.sqrt(2.*m.pi*(Vb+yerr**2.))*sc.exp(-0.5*(Y-Yb)**2./(Vb+yerr**2.))))#+pars[4]

 def __init__(self, module, dataset, totalIterations = 100,
              xPrecision = finfo(float).eps, fPrecision = finfo(float).eps,
              init_scg=True, **kwargs):
     """Create a SCGTrainer to train the specified `module` on the
     specified `dataset`.
     """
     Trainer.__init__(self, module)
     self.setData(dataset)
     self.input_sequences = self.ds.getField('input')
     self.epoch = 0
     self.totalepochs = 0
     self.module = module
     #self.tmp_module = module.copy()
     if init_scg:
         self.scg = SCG(self.module.params, self.f, self.df, self,
                        totalIterations, xPrecision, fPrecision,
                        evalFunc = lambda x: str(x / self.ds.getLength()))
     else:
         print "Warning: SCG trainer not initialized!"

    def __init__(self, N, uni_ground):
        self.odr = 'C'
        self.typ = sp.complex128

        self.zero_tol = sp.finfo(self.typ).resolution
        """Tolerance for detecting zeros. This is used when (pseudo-) inverting 
           l and r."""

        self._sanity_checks = False

        self.N = N
        """The number of sites. Do not change after initializing."""
        
        self.N_centre = N / 2
        """The 'centre' site. This affects the gauge-fixing and canonical
           form. It is the site between the left-gauge parts and the 
           right-gauge parts."""
           
        self.D = sp.repeat(uni_ground.D, self.N + 2)
        """Vector containing the bond-dimensions. A[n] is a 
           q[n] x D[n - 1] x D[n] tensor."""
           
        self.q = sp.repeat(uni_ground.q, self.N + 2)
        """Vector containing the site Hilbert space dimensions. A[n] is a 
           q[n] x D[n - 1] x D[n] tensor."""
           
        self.S_hc = sp.repeat(sp.NaN, self.N + 1)
        """Vector containing the von Neumann entropy S_hc[n] corresponding to 
           splitting the state between sites n and n + 1. Available only
           after performing update(restore_CF=True) or restore_CF()."""   

        self.uni_l = copy.deepcopy(uni_ground)
        self.uni_l.symm_gauge = False
        self.uni_l.sanity_checks = self.sanity_checks
        self.uni_l.update()

        self.uni_r = copy.deepcopy(uni_ground)
        self.uni_r.sanity_checks = self.sanity_checks
        self.uni_r.symm_gauge = False
        self.uni_r.update()

        self.grown_left = 0
        self.grown_right = 0
        self.shrunk_left = 0
        self.shrunk_right = 0

        self._init_arrays()

        for n in xrange(self.N + 2):
            self.A[n][:] = self.uni_l.A

        self.r[self.N] = self.uni_r.r
        self.r[self.N + 1] = self.r[self.N]
        self.l[0] = self.uni_l.l

    def __init__(self, N, uni_ground):
        self.odr = 'C'
        self.typ = sp.complex128

        self.zero_tol = sp.finfo(self.typ).resolution
        """Tolerance for detecting zeros. This is used when (pseudo-) inverting 
           l and r."""

        self._sanity_checks = False

        self.N = N
        """The number of sites. Do not change after initializing."""
        
        self.N_centre = N / 2
        """The 'centre' site. This affects the gauge-fixing and canonical
           form. It is the site between the left-gauge parts and the 
           right-gauge parts."""
           
        self.D = sp.repeat(uni_ground.D, self.N + 2)
        """Vector containing the bond-dimensions. A[n] is a 
           q[n] x D[n - 1] x D[n] tensor."""
           
        self.q = sp.repeat(uni_ground.q, self.N + 2)
        """Vector containing the site Hilbert space dimensions. A[n] is a 
           q[n] x D[n - 1] x D[n] tensor."""
          
        self.uni_l = copy.deepcopy(uni_ground)
        self.uni_l.symm_gauge = True
        self.uni_l.sanity_checks = self.sanity_checks
        self.uni_l.update()
        
        if not N % self.uni_l.L == 0:
            print "Warning: Length of nonuniform window is not a multiple of the uniform block size."

        self.uni_r = copy.deepcopy(self.uni_l)

        self.grown_left = 0
        self.grown_right = 0
        self.shrunk_left = 0
        self.shrunk_right = 0

        self._init_arrays()

        for n in xrange(1, self.N + 1):
            self.A[n][:] = self.uni_l.A[(n - 1) % self.uni_l.L]
        
        for n in xrange(self.N + 2):
            self.r[n][:] = sp.asarray(self.uni_l.r[(n - 1) % self.uni_l.L])
            self.l[n][:] = sp.asarray(self.uni_l.l[(n - 1) % self.uni_l.L])

 def billard_setup_strip(self, R1, R2, d, N):
     # d = increment of radius
     # N = number of balls for the unit circle
     rads = sp.linspace(R1, R2, (R2-R1)/d)
     pts = []
     for r in rads:
         num = r * N
         base = sp.linspace(0, (2-sp.finfo(float).eps)*sp.pi, num) # split up the circumference to num pieces
         for arg in base: 
             x = r*(sp.cos(arg))
             y = r*(sp.sin(arg))
             pts.append((x,y))
     
     self.balls = sp.array(pts)
     self.corner_slopes()

def safe_logdet(cov):
    '''
    The function computes a secure version of the logdet of a covariance matrix and it returns the rcondition number of the matrix.
    Inputs:
        cov
    Outputs:
        the logdet
    '''
    eps = sp.finfo(sp.float64).eps
    e = linalg.eigvalsh(cov)
    if e.max()<eps:
        rcond = 0
    else:
        rcond = e.min()/e.max()    
    e = sp.where(e<eps,eps,e)    
    return sp.sum(sp.log(e)),rcond

    def __init__(self, numsites, uni_ground):
        self.u_gnd_l = uni.EvoMPS_TDVP_Uniform(uni_ground.D, uni_ground.q)
        self.u_gnd_l.sanity_checks = self.sanity_checks
        self.u_gnd_l.h_nn = uni_ground.h_nn
        self.u_gnd_l.h_nn_cptr = uni_ground.h_nn_cptr
        self.u_gnd_l.A = uni_ground.A.copy()
        self.u_gnd_l.l = uni_ground.l.copy()
        self.u_gnd_l.r = uni_ground.r.copy()

        self.u_gnd_l.symm_gauge = False
        self.u_gnd_l.update()
        self.u_gnd_l.calc_lr()
        self.u_gnd_l.calc_B()
        self.eta_uni = self.u_gnd_l.eta
        self.u_gnd_l_kmr = la.norm(self.u_gnd_l.r / la.norm(self.u_gnd_l.r) - 
                                   self.u_gnd_l.K / la.norm(self.u_gnd_l.K))

        self.u_gnd_r = uni.EvoMPS_TDVP_Uniform(uni_ground.D, uni_ground.q)
        self.u_gnd_r.sanity_checks = self.sanity_checks
        self.u_gnd_r.symm_gauge = False
        self.u_gnd_r.h_nn = uni_ground.h_nn
        self.u_gnd_r.h_nn_cptr = uni_ground.h_nn_cptr
        self.u_gnd_r.A = self.u_gnd_l.A.copy()
        self.u_gnd_r.l = self.u_gnd_l.l.copy()
        self.u_gnd_r.r = self.u_gnd_l.r.copy()
        
        self.grown_left = 0
        self.grown_right = 0
        self.shrunk_left = 0
        self.shrunk_right = 0

        self.h_nn = self.wrap_h
        self.h_nn_mat = None

        self.eps = sp.finfo(self.typ).eps

        self.N = numsites

        self._init_arrays()

        for n in xrange(self.N + 2):
            self.A[n][:] = self.u_gnd_l.A

        self.r[self.N] = self.u_gnd_r.r
        self.r[self.N + 1] = self.r[self.N]
        self.l[0] = self.u_gnd_l.l

def evaluate_hull(x,hull):
    """evaluate_hull: evaluate h_u(x) and (optional) h_l(x)

    Input:
       x     - abcissa
       hull  - the hull (see setup_hull for a definition)

    Output:
      hu(x) (optional), hl(x)

    History:
       2009-05-21 - Written - Bovy (NYU)
    """
    #Find in which [z_{i-1},z_i] interval x lies
    if x < hull[4][0]:
        #x lies in the first interval
        hux= hull[3][0]*(x-hull[1][0])+hull[2][0]
        indx= 0
    else:
        if len(hull[5]) == 1:
            #There are only two intervals
            indx= 1
        else:
            indx= 1
            while indx < len(hull[4]) and hull[4][indx] < x:
                indx= indx+1
            indx= indx-1
        hux= hull[3][indx]*(x-hull[1][indx])+hull[2][indx]
    #Now evaluate hlx
    neginf= sc.finfo(sc.dtype(sc.float64)).min
    if x < hull[1][0] or x > hull[1][-1]:
        hlx= neginf
    else:
        if indx == 0:
            hlx= ((hull[1][1]-x)*hull[2][0]+(x-hull[1][0])*hull[2][1])/(hull[1][1]-hull[1][0])
        elif indx == len(hull[4]):
            hlx= ((hull[1][-1]-x)*hull[2][-2]+(x-hull[1][-2])*hull[2][-1])/(hull[1][-1]-hull[1][-2])
        elif x < hull[1][indx+1]:
            hlx= ((hull[1][indx+1]-x)*hull[2][indx]+(x-hull[1][indx])*hull[2][indx+1])/(hull[1][indx+1]-hull[1][indx])
        else:
            hlx= ((hull[1][indx+2]-x)*hull[2][indx+1]+(x-hull[1][indx+1])*hull[2][indx+2])/(hull[1][indx+2]-hull[1][indx+1])
    return hux, hlx

    def learn(self,x,y):
        '''
        Function that learns the GMM with ridge regularizationb from training samples
        Input:
            x : the training samples
            y :  the labels
        Output:
            the mean, covariance and proportion of each class, as well as the spectral decomposition of the covariance matrix
        '''

        ## Get information from the data
        C = sp.unique(y).shape[0]
        #C = int(y.max(0))  # Number of classes
        n = x.shape[0]  # Number of samples
        d = x.shape[1]  # Number of variables
        eps = sp.finfo(sp.float64).eps

        ## Initialization
        self.ni = sp.empty((C,1))    # Vector of number of samples for each class
        self.prop = sp.empty((C,1))  # Vector of proportion
        self.mean = sp.empty((C,d))  # Vector of means
        self.cov = sp.empty((C,d,d)) # Matrix of covariance
        self.Q = sp.empty((C,d,d)) # Matrix of eigenvectors
        self.L = sp.empty((C,d)) # Vector of eigenvalues
        self.classnum = sp.empty(C).astype('uint8')
        
        ## Learn the parameter of the model for each class
        for c,cR in enumerate(sp.unique(y)):
            
            j = sp.where(y==(cR))[0]
            
            self.classnum[c] = cR # Save the right label
            self.ni[c] = float(j.size)    
            self.prop[c] = self.ni[c]/n
            self.mean[c,:] = sp.mean(x[j,:],axis=0)
            self.cov[c,:,:] = sp.cov(x[j,:],bias=1,rowvar=0)  # Normalize by ni to be consistent with the update formulae

            # Spectral decomposition
            L,Q = linalg.eigh(self.cov[c,:,:])
            idx = L.argsort()[::-1]
            self.L[c,:] = L[idx]
            self.Q[c,:,:]=Q[:,idx]

def mylstsq(a, b, rcond):
    '''
    Compute a safe/fast least square fitting. 
    For that particular case, the number of unknown parameters is equal to the number of equations. 
    However, a is a covariance matrix that might estimated with a number of samples less than the number of variables, leading to a badly conditionned covariance matrix.
    So, the rcond number is check: if it is ok, we use the fast linalg.solve function; otherwise, we use the slow, but safe, linalg.lstsq function.
    Inputs:
    a: a symmetric definite positive matrix d times d
    b: a d times n matrix
    Outputs:
    x: a d times n matrix
    '''
    eps = sp.finfo(sp.float64).eps
    if rcond>eps: # If the condition number is not too bad try linear system
        try:
            x = linalg.solve(a,b)
        except linalg.LinAlgError: # If error, use least square estimations
            x = linalg.lstsq(a,b)[0]
    else:# Use least square estimate
        x = linalg.lstsq(a,b)[0]
    return x

    def newton(self, y, t=1, tol=None, maxit=100, verbose=False):
        """ Newton iteration

        Parameters
        ------------
        y : ndarray, shape (155, )
           Initial parameter guess. Changed in place.
        tol : float
           Tolerance for convergence
        maxit : int
           Maximum number of iterations
        verbose : bool
            Whether to print messages during iterations

        Returns
        -------
        info: integer
           Did function converge? True if it did.
        relres : float
           Relative residual of last iteration
        i : int
           Number of iterations


        """
        if tol is None:
            tol = sp.sqrt(sp.finfo(float).eps)
        info = False
        for i in range(maxit):
            f, df, dft = self.func(y, t)
            dy = la.solve(df, f)
            y -= dy
            relres = la.norm(dy)
            if verbose:
                print("%d: %f" % (i, relres))
            if relres  < tol:
                info = True
                break
        return (info, relres, i)

def fixpoint(f, x, tol=None, maxit=100, **kwargs):
    """ Fixed point function iteration

    Parameters
    -------------
    f : function
      Function for which to find the fixed point.
    x : array 
      Initial guess.
    tol : float, optional
      Tolerance for convergence.
    maxit : int, optional
      Maximum number of iterations.

    Returns
    -----------
    info : integer
       Did function converge? -1 if it did not, and 0 if it did.
    relres : float
       Relative residual at last iteration.
    t : int
       Number of iterations
    gval : array
       Fixed point, where f(x) = x.
       
    """
    if tol is None:
        tol = sqrt(sp.finfo(float).eps)
    info = -1
    t = 0
    for it in range(maxit):
        t += 1
        gval = f(x, **kwargs)
        relres = la.norm(gval - x)
        if relres < tol:
            info = 0
            break
        x = gval
    return (info, relres, t, gval)

    def predict(self, xt, x, y, out_decision=None, out_proba=None):
        nt = xt.shape[0]
        C = int(y.max())
        D = sp.empty((nt, C))
        D += self.prop
        eps = sp.finfo(sp.float64).eps

        # Pre compute the Gramm kernel matrix
        Kt = KERNEL()
        Kt.compute_kernel(xt, z=x, sig=self.sig)
        Ki = KERNEL()

        for i in range(C):
            t = sp.where(y == (i + 1))[0]
            Ki.compute_kernel(x, z=x[t, :], sig=self.sig)
            T = Kt.K - sp.dot(Ki.K, sp.ones((self.ni[i])) / self.ni[i])
            temp = sp.dot(T, self.S)
            D[:, i] = sp.sum(T * temp, axis=1)

        # Check if negative value
        if D.min() < 0:
            D -= D.min()

        yp = D.argmin(1) + 1
        yp.shape = (nt, 1)

        # Format the output
        if out_proba is None:
            if out_decision is None:
                return yp
            else:
                return yp, D
        else:
            # Compute posterior !! Should be changed to a safe version
            P = sp.exp(-0.5 * D)
            P /= sp.sum(P, axis=1).reshape(nt, 1)
            P[P < eps] = 0
        return yp, D, P

    def decomposition(self, M):
        """
            Compute the decompostion of symmetric matrix
            Inputs:
                M:   matrix to decompose
            Outputs:
                vp:    eigenvalues
                Q:     eigenvectors
                rcond: conditioning
        """
        # Decomposition
        vp,Q = linalg.eigh(M)

        # Compute conditioning
        eps = sp.finfo(sp.float64).eps
        if vp.max()<eps:
            rcond = 0
        else:
            rcond = vp.min()/vp.max()

        vp[vp<eps] = eps

        return vp,Q,rcond