def grad_EVzxVzxT_by_hyper_exact(self, EVzxVzxT_list_this, Z, A, B, hyperno):

        P = Z.shape[0]
        R = Z.shape[1]
        N = A.shape[0]

        if hyperno != 0:
            return EVzxVzxT_list_this * 0

        alpha = self.length_scale * self.length_scale

        I = np.identity(R)
        S = np.diag(B[0, :] * B[0, :])
        Sinv = np.diag(1 / B[0, :] * B[0, :])
        C = I * alpha
        Cinv = I * (1 / alpha)
        CinvSinv = 2 * Cinv + Sinv
        CinvSinv_inv = np.diag(1 / CinvSinv.diagonal())

        dC = self.length_scale * I
        dCinv = -Cinv.dot(dC).dot(Cinv)
        dCinvSinv = 2 * dCinv
        dCinvSinv_inv = -CinvSinv_inv.dot(dCinvSinv).dot(CinvSinv_inv)

        S1 = (
            dCinv
            - dCinv.dot(CinvSinv_inv).dot(Cinv)
            - Cinv.dot(dCinvSinv_inv).dot(Cinv)
            - Cinv.dot(CinvSinv_inv).dot(dCinv)
        )
        S2 = -Sinv.dot(dCinvSinv_inv).dot(Sinv)
        S3 = Sinv.dot(dCinvSinv_inv).dot(Cinv) + Sinv.dot(CinvSinv_inv).dot(dCinv)
        S4 = dCinv.dot(CinvSinv_inv).dot(Cinv) + Cinv.dot(dCinvSinv_inv).dot(Cinv) + Cinv.dot(CinvSinv_inv).dot(dCinv)

        T1s = np.tile(Z.dot(S1).dot(Z.T).diagonal(), [P, 1])
        T1 = np.tile(T1s, [N, 1, 1])
        T2s = T1s.T
        T2 = np.tile(T2s, [N, 1, 1])
        T3 = np.tile(Z.dot(S4).dot(Z.T), [N, 1, 1])
        T4 = np.tile(A.dot(S2).dot(A.T).diagonal(), [P, 1]).T
        T4 = np.expand_dims(T4, axis=2)
        T4 = np.repeat(T4, P, axis=2)
        T5 = A.dot(S3).dot(Z.T)
        T5 = np.expand_dims(T5, axis=2)
        T5 = np.repeat(T5, P, axis=2)
        T6 = np.swapaxes(T5, 1, 2)

        SCinvI = 2 * Cinv.dot(S) + I
        SCinvI_inv = np.diag(1 / SCinvI.diagonal())
        (temp, logDetSCinvI) = np.linalg.slogdet(SCinvI)
        detSCinvI = np.exp(logDetSCinvI)
        dDetSCinvI = -0.5 * np.power(detSCinvI, -0.5) * SCinvI_inv.dot(2 * dCinv).dot(S).trace()

        expTerm = EVzxVzxT_list_this / np.power(detSCinvI, -0.5)

        res = EVzxVzxT_list_this * (-0.5 * T1 - 0.5 * T2 + T3 - 0.5 * T4 + T5 + T6) + dDetSCinvI * expTerm

        res = np.sum(res, axis=0)

        return res

    def test_convolve2d_king(self):
        
        gfn = lambda r, s: np.power(2*np.pi*s**2,-1)*np.exp(-r**2/(2*s**2))
        kfn = lambda r, s, g: np.power(2*np.pi*s**2,-1)*(1.-1./g)* \
            np.power(1+0.5/g*(r/s)**2,-g)

        kfn0 = lambda x, y, mux, muy, s, g: kfn(np.sqrt((x-mux)**2+(y-muy)**2),s,g)

        xaxis = Axis.create(-3,3,501)
        yaxis = Axis.create(-3,3,501)

        x, y = np.meshgrid(xaxis.center,yaxis.center)
        xbin, ybin = np.meshgrid(xaxis.width,yaxis.width)

        r = np.sqrt(x**2+y**2)

        # Scalar Input

        mux = 0.5
        muy = -0.2
        mur = (mux**2+muy**2)**0.5

        gsig = 0.1
        ksig = 0.2
        kgam = 4.0

        fval0 = np.sum(kfn0(x,y,mux,muy,ksig,kgam)*gfn(r,gsig)*xbin*ybin)
        fval1 = convolve2d_king(lambda t: gfn(t,gsig),mur,ksig,kgam,3.0,
                                nstep=10000)
#        fval2 = convolve2d_gauss(lambda t: kfn(t,ksig,kgam),mur,gsig,3.0,
#                                 nstep=1000)
#        print fval0, fval1, fval2, fval1/fval0

        assert_almost_equal(fval0,fval1,4)

    def test_power_zero(self):
        # ticket #1271
        zero = np.array([0j])
        one = np.array([1 + 0j])
        cinf = np.array([complex(np.inf, 0)])
        cnan = np.array([complex(np.nan, np.nan)])

        def assert_complex_equal(x, y):
            x, y = np.asarray(x), np.asarray(y)
            assert_array_equal(x.real, y.real)
            assert_array_equal(x.imag, y.imag)

        # positive powers
        for p in [0.33, 0.5, 1, 1.5, 2, 3, 4, 5, 6.6]:
            assert_complex_equal(np.power(zero, p), zero)

        # zero power
        assert_complex_equal(np.power(zero, 0), one)
        assert_complex_equal(np.power(zero, 0 + 1j), cnan)

        # negative power
        for p in [0.33, 0.5, 1, 1.5, 2, 3, 4, 5, 6.6]:
            assert_complex_equal(np.power(zero, -p), cnan)
        assert_complex_equal(np.power(zero, -1 + 0.2j), cnan)

        def test_fast_power(self):
            x = np.array([1, 2, 3], np.int16)
            assert (x ** 2.00001).dtype is (x ** 2.0).dtype

    def reg_score_function(X, y, mean, scale, shape, skewness):
        """ GAS Skew t Regression Update term using gradient only - native Python function

        Parameters
        ----------
        X : float
            datapoint for the right hand side variable
    
        y : float
            datapoint for the time series

        mean : float
            location parameter for the Skew t distribution

        scale : float
            scale parameter for the Skew t distribution

        shape : float
            tail thickness parameter for the Skew t distribution

        skewness : float
            skewness parameter for the Skew t distribution

        Returns
        ----------
        - Score of the Skew t family
        """
        m1 = (np.sqrt(shape)*sp.gamma((shape-1.0)/2.0))/(np.sqrt(np.pi)*sp.gamma(shape/2.0))
        mean = mean + (skewness - (1.0/skewness))*scale*m1
        if (y-mean)>=0:
            return ((shape+1)/shape)*((y-mean)*X)/(np.power(skewness*scale,2) + (np.power(y-mean,2)/shape))
        else:
            return ((shape+1)/shape)*((y-mean)*X)/(np.power(scale,2) + (np.power(skewness*(y-mean),2)/shape))

def check_skew_expect(distfn, arg, m, v, s, msg):
    if np.isfinite(s):
        m3e = distfn.expect(lambda x: np.power(x-m, 3), arg)
        npt.assert_almost_equal(m3e, s * np.power(v, 1.5),
                decimal=5, err_msg=msg + ' - skew')
    else:
        npt.assert_(np.isnan(s))

    def molar_heat_capacity_p(self, pressure, temperature, volume, params):
        """
        Returns the heat capacity [J/K/mol] as a function of pressure [Pa]
        and temperature [K].
        """
        a, b, c = mt.tait_constants(params)
        T = temperature
        T_e = params['T_einstein']
        n = params['n']
        Pth = self.__relative_thermal_pressure(T, params)

        ksi_over_ksi_0 = einstein.molar_heat_capacity_v(T, T_e, n) \
                         / einstein.molar_heat_capacity_v(params['T_0'], T_e, n)

        dintVdpdT = (params['V_0'] * params['a_0'] * params['K_0'] * a * ksi_over_ksi_0) * (
            np.power((1. + b * (pressure - params['P_0'] - Pth)), 0. - c) - np.power((1. - b * Pth), 0. - c))

        dSdT0 = params['V_0'] * params['K_0'] * np.power((ksi_over_ksi_0 * params['a_0']), 2.0) * \
                (np.power((1. + b * (pressure - params['P_0'] - Pth)), -1. - c) -
                 np.power((1. + b * (-Pth)), -1. - c))

        x = T_e/T
        dCv_einstdT = -(einstein.molar_heat_capacity_v(T, T_e, n) *
                        ( 1 - 2./x + 2./(np.exp(x) - 1.) ) * x/T)

        dSdT1 = -dintVdpdT * dCv_einstdT \
                / einstein.molar_heat_capacity_v(T, T_e, n)

        dSdT = dSdT0 + dSdT1
        return self.molar_heat_capacity_p0(temperature, params) + temperature * dSdT

def power_sutera_reweighing(Y, f_pow=lambda x: 1):
    """Re-weights the time series giving more value to the values of the time
    serie when there are a low global activity.

    References
    ---------
    .. [1] Antonio Sutera et al. Simple connectome inference from partial
    correlation statistics in calcium imaging

    """

    ## 0. Prepare variables needed
    m = Y.shape[1]
    global_y = np.sum(Y, axis=1)

    ## 1. Transformation
    Yt = np.zeros(Y.shape)
    for j in range(m):
        Yt[:, j] = np.power((Y[:, j] + 1.),
                            np.power((1.+np.divide(1., global_y)),
                                     f_pow(global_y)))
    # Correct global 0
    Yt[global_y == 0, :] = 1.

    return Yt

    def velocidadlimpio(self):
        #sino no llega al 2, asi funciona        
        v = np.arange(-3, 3,1)
        [x,y] = np.meshgrid(v,v)

        z=np.multiply(x,np.exp( -np.power(x,2) - np.power(y,2) ))

        #Matplotlib t invierte el orden de las matrices a diferenciade matlab
        [py,px] = np.gradient(z,1,1)

        print 'x '+str(x)
        print 'y '+str(y)
        print 'z '+str(z)
        print 'px '+str(px)
        print 'py '+str(py)
        
        #q = plt.quiver(X, Y, u, v, angles='xy', scale=40, color=['r'])
        #p = plt.quiverkey(q,1,16.5,50,"50 m/s",coordinates='data',color='r')        

        #primero los rangos, despues los valores q contiene
        q = plt.quiver(x,y, px, py)
        p = plt.quiverkey(q,1,16.5,50,"50 m/s",coordinates='data',color='r')
        plt.title('Velocidad')
        plt.show()
        print 'Fourth plot loaded...'  

    def p_corr(self):
        """
        calculate pearson correlation between users
        """
        data = self.data
        rows = self.rows
        (nrows, ncols) = data.shape
        p_corr_dict = {}

        for row_i in range(nrows):
            for row_j in range(nrows):
                valide_data_i = [data[row_i, :][n] \
                            for n in range(ncols) if \
                            data[row_i, n] != 0 and data[row_j, n] != 0]
                valide_data_j = [data[row_j, :][n] \
                                for n in range(ncols) if \
                                data[row_i, n] != 0 and data[row_j, n] != 0] 
                valide_data_i = np.array(valide_data_i) - np.mean(valide_data_i)
                valide_data_j = np.array(valide_data_j) - np.mean(valide_data_j)
                #print np.dot(data[row_i, :], data[row_j, :])
                p_corr = np.dot(valide_data_i, valide_data_j)*1.0/\
                         np.sqrt(sum(np.power(valide_data_i,2))*\
                                 sum(np.power(valide_data_j,2)))
                p_corr_dict[(rows[row_i], rows[row_j])] = p_corr  
        return p_corr_dict

def lee_parallel_betachEq11(x,n):
    beta=x.copy()
    beta_ch=11.6
    rmsV=5.
    cs=1.
    mach_bh= lee_mach_bh(rmsV,cs)
    return np.log( mach_bh**-2*np.power(mach_bh**n+np.power(beta_ch/beta,n/2.), -1./n) )

def read_excel(file='slab++.xlsx'):
    wb = open_workbook(file)
    sheet = wb.sheets()[0]
    number_of_rows = sheet.nrows
    number_of_columns = sheet.ncols
        
    items = []
    
    rows = []
    slabs= Catalog()
    for row in range(4, number_of_rows):
        slab= Slab()
        values = []
        for col in range(number_of_columns):
            par=str(sheet.cell(3,col).value)
            if par:
                print(par)
            value  = sheet.cell(row,col).value
            print(par,value)
            try:
                slab.params[par]= float(value)
            except:
                slab.params[par] = str(value)
                
        # save temperature at 600 km using thermal parameter
        phi = slab.params['thermalPar']
        z = 6.0
        Ta = 1338.0
        Tz = Ta * (1.0 - ( (2.0/np.pi)*np.exp(-1.0*( (np.power(np.pi,2.0)*z)/ (np.power(2.32,2.0) * phi) )) ) )
        slab.params['Temp600'] = Tz
        print(Tz)
        print(slab.params)
        slabs.append(slab)
    return slabs

def lee_parallel_neq8(x,beta_ch):
    beta=x.copy()
    n=8.
    rmsV=5.
    cs=1.
    mach_bh= lee_mach_bh(rmsV,cs)
    return np.log( mach_bh**-2*np.power(mach_bh**n+np.power(beta_ch/beta,n/2.), -1./n) )

def mdot_magnetic_neq8(x,beta_ch):
    beta=x.copy()
    n=8.
    rmsV=5.
    cs=1.
    mach_bh= lee_mach_bh(rmsV,cs)
    return mach_bh**-2*np.power(mach_bh**n+np.power(beta_ch/beta,n/2.), -1./n)

def evaluation_10_fold(root='./result/pytorch_result.mat'):
    ACCs = np.zeros(10)
    result = scipy.io.loadmat(root)
    for i in range(10):
        fold = result['fold']
        flags = result['flag']
        featureLs = result['fl']
        featureRs = result['fr']

        valFold = fold != i
        testFold = fold == i
        flags = np.squeeze(flags)

        mu = np.mean(np.concatenate((featureLs[valFold[0], :], featureRs[valFold[0], :]), 0), 0)
        mu = np.expand_dims(mu, 0)
        featureLs = featureLs - mu
        featureRs = featureRs - mu
        featureLs = featureLs / np.expand_dims(np.sqrt(np.sum(np.power(featureLs, 2), 1)), 1)
        featureRs = featureRs / np.expand_dims(np.sqrt(np.sum(np.power(featureRs, 2), 1)), 1)

        scores = np.sum(np.multiply(featureLs, featureRs), 1)
        threshold = getThreshold(scores[valFold[0]], flags[valFold[0]], 10000)
        ACCs[i] = getAccuracy(scores[testFold[0]], flags[testFold[0]], threshold)
    #     print('{}    {:.2f}'.format(i+1, ACCs[i] * 100))
    # print('--------')
    # print('AVE    {:.2f}'.format(np.mean(ACCs) * 100))
    return ACCs

def p_jds(Z, L, U, C, M):
    """
    Input:
         coefficient matrix Z
         desired numbers of dynamic active sets L
         dictionary atom label vector U
         the number of classes C
         the number of views M
    Output:
         Index matrix I for top-L dynamic active sets
    """   
    #Initialize
    I = np.zeros((L, M))
    V = np.zeros((C, M))
    _I = np.zeros((C, M))
    S = np.zeors(C)
    for l in xrange(L):
        for i in xrange(C):
            #c代表索引值
            c = find(U, i)
            for m in xrange(M):
                v, t = Max(Z[c, m])
                V[i, m] = v
                _I[i, m] = c[t]
            tmp = np.cumsum(np.power(V[i], 2))[-1]
            S[i] = np.power(tmp, 0.5)
        _v, _t = Max(S)
        I[l, :] = _I[_t,:]
        Z[_I[_t,:]] = 0
    return I 

def find_nearby(targets, catalog, max_dist=10):
  '''
  Finds all the sources within a max distance of each target.  
  '''

  # Determines the squared separation between targets
  sep = lambda x0,y0,x,y: np.power((x-x0)*np.cos(y0*pi/180),2) + np.power((y-y0),2)
  
  # Convert max distance from arcseconds to squarded degrees
  max_dist = max_dist / 3600.**2

  for each in targets:

    rvec = sep(each['degra'], each['degdec'], catalog['degra'], catalog['degdec'])


    # Ensure that target is recovered in catalog   
    nearby = np.where(rvec < max_dist)

    if each['hstid'] in catalog['hstid'][nearby]:
      write_cmd_file(catalog[nearby], each)

    else: 
      continue    
      print 'ERROR: Failed to find %s.' % each['hstid']

def _fit_quad_to_peak(x, y):
    """
    Fits a quadratic to the data points handed in
    to the from y = b[0](x-b[1])**2 + b[2]

    x -- locations
    y -- values

    returns (b, R2)

    """

    lenx = len(x)

    # some sanity checks
    if lenx < 3:
        raise Exception('insufficient points handed in ')
    # set up fitting array
    X = np.vstack((x ** 2, x, np.ones(lenx))).T
    # use linear least squares fitting
    beta, _, _, _ = np.linalg.lstsq(X, y)

    SSerr = np.sum(np.power(np.polyval(beta, x) - y, 2))
    SStot = np.sum(np.power(y - np.mean(y), 2))
    # re-map the returned value to match the form we want
    ret_beta = (beta[0],
                -beta[1] / (2 * beta[0]),
                beta[2] - beta[0] * (beta[1] / (2 * beta[0])) ** 2)

    return ret_beta, 1 - SSerr / SStot

def calc_v(pos, t, x_table, c_attract, dc_attract, c_repel, dc_repel):
    '''Calculate drift velocity as a function of c(pos,t) and dc(pos,t)'''
    v = 22. # µm/s
    chi = 50000. # µm^2/s
    chi2 = 50000. # µm^2/s
    k = 0.125 # mM
    k2 = 0.125 # mM
    if pos < x_table[0]:
        pos = x_table[0]
        c, dc = int_tables(pos, t, x_table, c_attract, dc_attract)
        c2, dc2 = int_tables(pos, t, x_table, c_repel, dc_repel)
    elif pos > x_table[-1]:
        pos = x_table[-1]
        c, dc = int_tables(pos, t, x_table, c_attract, dc_attract)
        c2, dc2 = int_tables(pos, t, x_table, c_repel, dc_repel)
    else:
        c, dc = int_tables(pos, t, x_table, c_attract, dc_attract)
        c2, dc2 = int_tables(pos, t, x_table, c_repel, dc_repel)
    v1 = ((8*v)/(3*np.pi))*np.tanh(((chi*np.pi)/(8*v))*(k/np.power((k+c),2))*dc)
    v2 = ((8*v)/(3*np.pi))*np.tanh(((chi2*np.pi)/(8*v))*(k2/np.power((k+c2),2))*dc2)
    #print v1, v2
    if np.abs(v1) > np.abs(v2):
        return v1
    else:
        return v2

def EN_CID(y):
    """
    CID measure from Batista, G. E. A. P. A., Keogh, E. J., Tataw, O. M. & de
    Souza, V. M. A. CID: an efficient complexity-invariant distance for time
    series. Data Min Knowl. Disc. 28, 634-669 (2014).
    
    Arguments
    ---------

    y: a nitime time-series object, or numpy vector

    """

    # Make the input a row vector of numbers:
    y = makeRowVector(vectorize(y))

    # Prepare the output dictionary
    out = {}
    
     # Original definition (in Table 2 of paper cited above)
    out['CE1'] = np.sqrt(np.mean(np.power(np.diff(y),2))); # sum -> mean to deal with non-equal time-series lengths

    # Definition corresponding to the line segment example in Fig. 9 of the paper
    # cited above (using Pythagoras's theorum):
    out['CE2'] = np.mean(np.sqrt(1 + np.power(np.diff(y),2)));

    return out

def mjt_discrete_movement(amp=1.0, dur=1.0, loc=0.,
                          time=np.arange(-0.5, 0.5, 0.01)):
    """
    Generate a discrete Minumum Jerk Trajectory (MJT) movement speed profile
    of the given amplitude, duration and time location for the given time span.

    Parameters
    ----------
    amp         : float
                  Amplitude of the MJT discrete movement.
    dur         : float
                  Duration of the MJT discrete movement.
    loc         : float
                  The temporal location of the center of the MJT speed
                  profile.
    time        : np.array
                  The time values for which the speed profile values are to be
                  returned.
    Returns
    -------
    movement    : np.array
                  The movement speed profile of the MJT discrete movement.

    Notes
    -----

    Examples
    --------
    """
    t = np.array([np.min([np.max([(_t + 0.5 * dur - loc) / dur, 0.]), 1.])
                  for _t in time])
    return amp * (30 * np.power(t, 4) -
                  60 * np.power(t, 3) +
                  30 * np.power(t, 2))