def setUp(self):
        SP.random.seed(1)
        self.n_dimensions=2
        self.n_samples = 100
        self.setXy()
        #self.X = SP.rand(self.n_samples,self.n_dimensions)
        covar  = limix.CCovSqexpARD(self.n_dimensions)
        ll  = limix.CLikNormalIso()
        covar_params = SP.array([1,1,1])
        lik_params   = SP.array([0.5])
        hyperparams0 = limix.CGPHyperParams()
        hyperparams0['covar'] = covar_params
        hyperparams0['lik'] = lik_params
        self.constrainU = limix.CGPHyperParams()
        self.constrainL = limix.CGPHyperParams()
        self.constrainU['covar'] = +10*SP.ones_like(covar_params);
        self.constrainL['covar'] = 0*SP.ones_like(covar_params);
        self.constrainU['lik'] = +5*SP.ones_like(lik_params);
        self.constrainL['lik'] = 0*SP.ones_like(lik_params);

        self.gp=limix.CGPbase(covar,ll)
        self.gp.setX(self.X)
        self.gp.setParams(hyperparams0)
        #self.genY()
        self.gp.setY(self.y)

 def run(self,phase,inlets,throat_prop='throat.capillary_pressure'):
     r'''
     Perform the algorithm
     
     Parameters
     ----------
     phase : OpenPNM Phase object
         The phase to be injected into the Network.  The Phase must have the
         capillary entry pressure values for the system.
         
     inlets : array_like
         The list of inlet pores from which the Phase can enter the Network
         
     throat_prop : string
         The name of the throat property containing the capillary entry
         pressure.  The default is 'throat.capillary_pressure'.
     
     '''
     import heapq as hq
     queue = []
     hq.heapify(queue)
     self._phase = phase
     net = self._net
     # Setup arrays and info
     t_entry = phase[throat_prop]
     t_sorted = sp.argsort(t_entry,axis=0)  # Indices into t_entry giving a sorted list
     t_order = sp.zeros_like(t_sorted)
     t_order[t_sorted] = sp.arange(0,net.Nt)  # Location in sorted list
     t_inv = -sp.ones_like(net.Ts)  # List for tracking throat invasion order
     p_inv = -sp.ones_like(net.Ps)  # List for tracking pore invasion order
     p_inv[inlets] = 0  # Set inlet pores to invaded
     # Perform initial analysis on input pores
     Ts = net.find_neighbor_throats(pores=inlets)
     [hq.heappush(queue,T) for T in t_order[Ts]]  # Push the new throats to the heap
     tcount = 1
     while len(queue) > 0:
         t = hq.heappop(queue)  # Find throat at the top of the queue
         t_next = t_sorted[t]  # Extract actual throat number
         t_inv[t_next] = tcount  # Note invasion sequence
         while (len(queue)>0) and (queue[0] == t):  # If throat is duplicated
             t = hq.heappop(queue)  # Note: Preventing duplicate entries below might save some time here
         Ps = net['throat.conns'][t_next]  # Find pores connected to newly invaded throat
         Ps = Ps[p_inv[Ps]<0]  # Remove already invaded pores from Ps
         if len(Ps)>0:
             p_inv[Ps] = tcount  # Note invasion sequence
             Ts = net.find_neighbor_throats(pores=Ps)  # Find connected throats
             Ts = Ts[t_inv[Ts]<0]  # Remove already invaded throats from Ts
             [hq.heappush(queue,T) for T in t_order[Ts]]  # Add new throats to queue
         tcount += 1
     self['throat.invasion_sequence'] = t_inv
     self['pore.invasion_sequence'] = p_inv

    def _set_info(self,element='',label='',locations='',mode='merge'):
        r'''
        This is the actual info setter method, but it should not be called directly.  
        Wrapper methods have been created.  Use set_pore_info and get_pore_info.
        
        See Also
        --------
        set_pore_info, set_throat_info
        '''
        if type(locations)==list: 
            try: locations = getattr(self,'get_'+element+'_indices')(locations)
            except: locations = sp.array(locations,ndmin=1)
        elif type(locations)==sp.ndarray:
            try: locations = getattr(self,'get_'+element+'_indices')(locations)
            except : pass 
        if locations!='':
            
            try: 
                locations = locations.name
                label = locations
            except: pass
            if type(locations)==str: locations = getattr(self,'get_'+element+'_indices')([locations])           
            locations=sp.array(locations,ndmin=1)
            if label:
                if label=='all':
                    try: 
                        old_label = getattr(self,'_'+element+'_info')[label]
                        if sp.shape(old_label)[0]<sp.shape(locations)[0]:
                            getattr(self,'_'+element+'_info')[label] = sp.ones_like(locations,dtype=bool)
                            self._logger.info('label=all for '+element+'has been updated to a bigger size!')
                            for info_labels in getattr(self,'_'+element+'_info').keys():
                                if info_labels!=label:
                                    temp = sp.zeros((getattr(self,'num_'+element+'s')(),),dtype=bool)
                                    temp[old_label] = getattr(self,'_'+element+'_info')[info_labels]
                                    getattr(self,'_'+element+'_info')[info_labels] = temp
                        elif sp.shape(old_label)[0]>sp.shape(locations)[0]: 
                            self._logger.error('To apply a new numbering label (label=all) to '+element+'s, size of the locations cannot be less than total number of '+element+'s!!')
                    except: getattr(self,'_'+element+'_info')[label] = sp.ones_like(locations,dtype=bool)
                else:    
                    try: getattr(self,'_'+element+'_info')[label]
                    except: getattr(self,'_'+element+'_info')[label] = sp.zeros((getattr(self,'num_'+element+'s')(),),dtype=bool)
                    if mode=='overwrite':
                        getattr(self,'_'+element+'_info')[label] = sp.zeros((getattr(self,'num_'+element+'s')(),),dtype=bool)
                        getattr(self,'_'+element+'_info')[label][locations] = True
                    elif mode=='remove':  getattr(self,'_'+element+'_info')[label][locations] = False                           
                    elif mode=='merge':  getattr(self,'_'+element+'_info')[label][locations] = True
            else: self._logger.error('No label has been defined for these locations')                

        elif mode=='remove':  del getattr(self,'_'+element+'_info')[label]
        else:  getattr(self,'_'+element+'_info')[label] = sp.zeros((getattr(self,'num_'+element+'s')(),),dtype=bool)

 def test_different_windows(self) :
     window1 = sp.ones_like(self.wave1)
     window1 += sp.sin(sp.arange(self.n)/40.0)
     window2 = 2*sp.ones_like(self.wave1)
     window2 += 2*sp.sin(sp.arange(self.n)/62.0)
     window2[window2<0.5] = 0
     wave1 = self.wave1*window1
     wave2 = self.wave1*window2
     power = np.windowed_power(wave1, window1, wave2, window2)
     self.assertAlmostEqual(power[self.mode1]/self.amp1**2/self.n*4, 1.0, 3)
     self.assertTrue(sp.allclose(power[:self.mode1], 0, 
                                 atol=self.amp1**2*self.n/1e3))
     self.assertTrue(sp.allclose(power[self.mode1+1:], 0, 
                                 atol=self.amp1**2*self.n/1e3))

def readMahalih5(filename,des_site):
    """ This function will read the mahali GPS data into a GeoData data structure.
        The user only has to give a filename and name of the desired site.
        Input
            filename - A string that holds the file name.
            des_site - The site name. Should be listed in the h5 file in the
                table sites.
    """
    h5fn = Path(filename).expanduser()

    with h5py.File(str(h5fn), "r", libver='latest') as f:

        despnts = sp.where(f['data']['site']==des_site)[0]
        # TODO: hard coded for now
        doy =  doy= f['data']['time'][despnts]
        year = 2015*sp.ones_like(doy,dtype=int)

        TEC = f['data']['los_tec'][despnts]

        nTEC = f['data']['err_los_tec'][despnts]

        vTEC = f['data']['vtec'][despnts]
        az2sat = f['data']['az'][despnts]
        el2sat = f['data']['az'][despnts]

        piercelat = f['data']['pplat'][despnts]
        piercelong = f['data']['pplon'][despnts]
        satnum= f['data']['prn'][despnts]
        recBias = f['data']['rec_bias'][despnts]
        nrecBias = f['data']['err_rec_bias'][despnts]

    # Make the integration time on the order of 15 seconds.
    if (year==year[1]).all():
        unixyear =(datetime(year[0],1,1,0,0,0,tzinfo=UTC) - EPOCH).total_seconds()
        uttime = unixyear + sp.round_(24*3600*sp.column_stack((doy,doy+15./24./3600.))) # Making the difference in time to be a minute
    else:
        (y_u,y_iv) = np.unique(year,return_inverse=True)
        unixyearu = sp.array([(datetime(iy,1,1,0,0,0,tzinfo=UTC) - EPOCH).total_seconds() for iy in y_u])
        unixyear = unixyearu[y_iv]
        uttime = unixyear + 24*3600*sp.column_stack((doy,doy+15./24./3600.))


    data = {'TEC':TEC,'nTEC':nTEC,'vTEC':vTEC,'recBias':recBias,'nrecBias':nrecBias,'satnum':satnum,'az2sat':az2sat,'el2sat':el2sat}
    coordnames = 'WGS84'
    sensorloc = sp.nan*sp.ones(3)
    dataloc = sp.column_stack((piercelat,piercelong, 350e3*sp.ones_like(piercelat)))

    return (data,coordnames,dataloc,sensorloc,uttime)

 def test_convolve_normalization(self) :
     window = sp.ones_like(self.wave1)
     power = npow.calculate_power(self.wave1)
     window_power = npow.calculate_power(window)
     # Convolving with the window function should do nothing for all ones.
     convolved_power = npow.convolve_power(power, window_power)
     self.assertTrue(sp.allclose(power, convolved_power))

 def _additionalInit(self):
     # default setting
     if self.slow_constant is None:
         self.slow_constant = max(1, int(self.paramdim / 10.))
     
     # get a few initial samples to work with
     tmp = self.batch_size
     self.batch_size = self.init_samples * tmp
     self._collectGradients()
     self._num_updates += self.init_samples
     self.batch_size = tmp
     
     # mean gradient vector
     self._gbar = mean(self._last_gradients, axis=0)
     # mean squared gradient
     self._vbar = (mean(self._last_gradients ** 2, axis=0) + self.epsilon) * self.slow_constant
     self._vpart = self._gbar ** 2 / self._vbar
     
     # mean diagonal Hessian
     #hs = clip(mean(self._last_diaghessians, axis=0), 1, 1 / self.epsilon)
     #self._hbar = hs * self.slow_constant
     self._hbar = mean(self._last_diaghessians, axis=0)
     
     # time constants
     self._taus = (ones_like(self.parameters) + self.epsilon) * 2#* self.init_samples
     
     # for debugging
     self._print_quantities = [('p', self.parameters),
                               ('tau', self._taus),
                               ('g', self._gbar),
                               ('v', self._vbar),
                               ('vpa', self._vpart),
                               ('h', self._hbar),
                               ]

 def _get_indices(self,element,labels,return_indices,mode):
     r'''
     This is the actual method for getting indices, but should not be called
     directly.  
     '''
     if mode == 'union':
         union = sp.zeros_like(self._get_info(element=element,label='all'),dtype=bool)
         for item in labels: #iterate over labels list and collect all indices
                 union = union + self._get_info(element=element,label=item)
         ind = union
     elif mode == 'intersection':
         intersect = sp.ones_like(self._get_info(element=element,label='all'),dtype=bool)
         for item in labels: #iterate over labels list and collect all indices
                 intersect = intersect*self._get_info(element=element,label=item)
         ind = intersect
     elif mode == 'not_intersection':
         not_intersect = sp.zeros_like(self._get_info(element=element,label='all'),dtype=int)
         for item in labels: #iterate over labels list and collect all indices
             info = self._get_info(element=element,label=item)
             not_intersect = not_intersect + sp.int8(info)
         ind = (not_intersect == 1)
     elif mode == 'none':
         none = sp.zeros_like(self._get_info(element=element,label='all'),dtype=int)
         for item in labels: #iterate over labels list and collect all indices
             info = self._get_info(element=element,label=item)
             none = none - sp.int8(info)
         ind = (none == 0)
     if return_indices: ind = sp.where(ind==True)[0]
     return ind

def steepest_descent(X,Y, step=.001, tol=1e-5, maxiter=5000, bounds=[-5,5], res=.1):
    w = betahat(X,Y)
    a = sp.exp(X.dot(w))
    yhat = (a.T/sp.sum(a, axis=1)).T
    grad = X.T.dot(yhat-Y)
    while la.norm(grad)>tol and maxiter > 0:
        w = w - grad*step
        a=sp.exp(X.dot(w))
        yhat = (a.T/sp.sum(a,axis=1)).T
        grad = X.T.dot(yhat-Y)
        maxiter -= 1

    rang = sp.arange(bounds[0],bounds[1]+res, res)
    Xg, Yg = sp.meshgrid(rang, rang)
    
    Xm = sp.c_[sp.ones_like(Xg.flatten()), Xg.flatten(), Yg.flatten()]
    Xmdot = Xm.dot(w)
    
    types = Xmdot.argmax(axis=1)
    print Xm.shape
    #plot the regions
    c = ['b','r','g']
    for i in range(Xmdot.shape[1]):
        plot_on = types==i
        plt.plot(Xm[plot_on,1], Xm[plot_on,2], c[i]+'.', X[Y[:,i]==1,1], X[Y[:,i]==1,2], c[i]+'o')

    #plot the data segmented
    #tmp = sp.where(Ydat[:,0]==True)[0]
    #plt.plot(Xdat[tmp,1], Xdat[tmp,2], 'o', Xdat[~tmp,1], Xdat[~tmp,2], 'o')
    plt.show()

def compactness(geometry, throat_perimeter='throat.perimeter',
                throat_area='throat.area', **kwargs):
    r"""
    Mortensen et al. have shown that the Hagen-Poiseuille hydraluic resistance is
    linearly dependent on the compactness. Defined as perimeter^2/area.
    The dependence is not universal as shapes with sharp corners provide more
    resistance than those that are more elliptical. Count the number of vertices
    and apply the right correction.
    """
    # Only apply to throats with an area
    ts = geometry.throats()[geometry[throat_area] > 0]
    P = geometry[throat_perimeter]
    A = geometry[throat_area]
    C = _sp.ones(geometry.num_throats())
    C[ts] = P[ts]**2/A[ts]
    verts = geometry['throat.offset_vertices']
    alpha = _sp.ones_like(C)
    for i in ts:
        if len(verts[i]) == 3:
            # Triangular Correction
            alpha[i] = C[i]*(25/17) + (40*_sp.sqrt(3)/17)
        elif len(verts[i]) == 4:
            # Rectangular Correction
            alpha[i] = C[i]*(22/7) - (65/3)
        elif len(verts[i]) > 4:
            # Approximate Elliptical Correction
            alpha[i] = C[i]*(8/3) - (8*_sp.pi/3)
    # For a perfect circle alpha = 8*pi so normalize by this
    alpha /= 8*_sp.pi
    # Very small throats could have values less than one
    alpha[alpha < 1.0] = 1.0

    return alpha

def compactness(target, throat_perimeter='throat.perimeter',
                throat_area='throat.area'):
    r"""
    Mortensen et al. have shown that the Hagen-Poiseuille hydraluic resistance
    is linearly dependent on the compactness. Defined as perimeter^2/area.
    The dependence is not universal as shapes with sharp corners provide more
    resistance than those that are more elliptical. Count the number of
    vertices and apply the right correction.

    Parameters
    ----------
    target : OpenPNM Object
        The object which this model is associated with. This controls the
        length of the calculated array, and also provides access to other
        necessary properties.

    throat_perimeter : string
        The dictionary key of the array containing the throat perimeter values.

    throat_area : string
        The dictionary key of the array containing the throat area values.

    Returns
    -------
    alpha : NumPy ndarray
        Array containing throat compactness values.

    References
    ----------
    Mortensen N.A, Okkels F., and Bruus H. Reexamination of Hagen-Poiseuille
    flow: Shape dependence of the hydraulic resistance in microchannels.
    Physical Review E, v.71, pp.057301 (2005).

    """
    # Only apply to throats with an area
    ts = target.throats()[target[throat_area] > 0]
    P = target[throat_perimeter]
    A = target[throat_area]
    C = _sp.ones(target.num_throats())
    C[ts] = P[ts]**2/A[ts]
    alpha = _sp.ones_like(C)*8*_sp.pi
    if 'throat.offset_vertices' in target.props():
        verts = target['throat.offset_vertices']
        for i in ts:
            if ~_sp.any(_sp.isnan(verts[i])):
                if len(verts[i]) == 3:
                    # Triangular Correction
                    alpha[i] = C[i]*(25/17) + (40*_sp.sqrt(3)/17)
                elif len(verts[i]) == 4:
                    # Rectangular Correction
                    alpha[i] = C[i]*(22/7) - (65/3)
                elif len(verts[i]) > 4:
                    # Approximate Elliptical Correction
                    alpha[i] = C[i]*(8/3) - (8*_sp.pi/3)
    # For a perfect circle alpha = 8*pi so normalize by this
    alpha /= 8*_sp.pi
    # Very small throats could have values less than one
    alpha[alpha < 1.0] = 1.0
    return alpha

 def test_find_path_with_weights(self):
     w = sp.ones_like(self.net.Ts)
     w[0] = self.net.Nt + 1
     a = misc.find_path(network=self.net,
                        pore_pairs=([0, 1]),
                        weights=w)
     assert len(a['pores'][0]) > 2
     assert len(a['throats'][0]) > 1

 def test_convolve(self) :
     window = sp.ones_like(self.wave1)
     window += sp.sin(sp.arange(self.n)/50.0)
     self.wave1 *= window
     power = npow.calculate_power(self.wave1)
     window_power = npow.calculate_power(window)
     deconvolved_power = npow.deconvolve_power(power, window_power)
     reconvolved_power = npow.convolve_power(deconvolved_power, window_power)
     self.assertTrue(sp.allclose(power, reconvolved_power))

 def predict(self, XTest=None, k=None, mean=None):
     # k is cross covariance KTestTrain
     self.mean = mean
     if self.mean is None:
         self.mean = SP.ones_like(self.yTrain)*self.yTrain.mean()
     if k is None:
         k = self.kernel
     Core = SP.dot(k, LA.inv((self.kernel + SP.eye(self.nTrain)*self.delta)))
     return SP.dot(Core, self.yTrain-self.mean)

def planeFromPoints(points):
    """ Produces a plane from N points using least squares, and returns of the form:
        Z = aX + bY + c
        Follows logic from http://www.velocityreviews.com/forums/t368189-re-linear-regression-in-3-dimensions.html
    """
    x, y, z = zip(*points)
    A = column_stack([x, y, ones_like(x)])
    abc, residuals, rank, s = lstsq(A,z)
    return abc

 def test_window_with_zeros(self) :
     window = sp.ones_like(self.wave1)
     window += sp.sin(sp.arange(self.n)/50.0)
     self.wave1 *= window
     power = np.windowed_power(self.wave1, window)
     self.assertAlmostEqual(power[self.mode1]/self.amp1**2/self.n*4, 1.0, 3)
     self.assertTrue(sp.allclose(power[:self.mode1], 0, 
                                 atol=self.amp1**2*self.n/1e3))
     self.assertTrue(sp.allclose(power[self.mode1+1:], 0, 
                                 atol=self.amp1**2*self.n/1e3))

def lower_confidence_bound(x_test, x_obs, y_obs, p=0.01, **kwargs):
    if len(x_obs) == 0:
        return sp.inf*sp.ones_like(x_test)
    
    mu_test, sigma_test = evaluate(x_test, x_obs, y_obs, **kwargs)
        
    std_test = sp.sqrt(sp.diag(sigma_test))
    lcb = mu_test.flatten() + sp.stats.norm.ppf(0.01)*std_test

    return lcb

def filt(k,kcut,beta,alph):

    #SPINS default parameters: 0.6, 2.0, 20.0

    knyq = max(k)
    kxcut = kcut*knyq
    filt = np.ones_like(k)
    filt = np.exp(-alph*((np.absolute(k) - kxcut)/(knyq - kxcut))**beta)*(np.absolute(k)>kxcut) + (np.absolute(k)<=kxcut)

    return filt

 def controlJacobian(self,v):
     eps = self.param['eps']
     # compute weights per bin for perturbation 
     # w = 1 + eps * v/rho
     w_bin = 1. + eps * v/norm(v)/self.u
     # transform to weights per particle
     w_part = w_bin[self.bin]
     total = scipy.sum(w_part)/self.N #divide by Nlarge
     # note that the perturbed state is no longer a probability distribution
     # so norming the histograms is not entirely correct
     u_eps_Dt = total*self.restriction.restrict(self.x_Dt,self.grid,self.domain,w=w_part)
     c = self.control
     # print "self.control", c,
     # print self.lambd
     if c==None:
         print "no control variable ..."
         control = scipy.zeros_like(v)
     else: 
         print "control variable ..."
         eps = self.param['eps']
         skip = self.skip
         Nlarge = self.param['Nlarge']
         Nsmall = self.param['Nsmall']
         w_part_large = scipy.ones_like(c.x)
         w_part_small = scipy.ones_like(self.x)
         w_bin_eps = 1. + eps * v/norm(v)/c.u
         w_eps_part_large = w_bin[c.bin]
         w_eps_part_small = w_bin[c.bin[skip/2:Nlarge:skip]]
         total_large = scipy.sum(w_part_large)/Nlarge
         total_small = scipy.sum(w_part_small)/Nsmall
         u_Dt_large = total_large*self.restriction.restrict(\
             c.x_Dt,c.grid,c.domain,w=w_part_large)
         u_eps_Dt_large = total_large*self.restriction.restrict(\
             c.x_Dt,c.grid,c.domain,w=w_eps_part_large)
         control_large = u_eps_Dt_large - u_Dt_large
         u_Dt_small = total_small*self.restriction.restrict(\
             c.x_Dt[skip/2:Nlarge:skip],c.grid,c.domain,w=w_part_small)
         u_eps_Dt_small = total_small*self.restriction.restrict(\
             c.x_Dt[skip/2:Nlarge:skip],c.grid,c.domain,w=w_eps_part_small)
         control_small = (u_eps_Dt_small - u_Dt_small)
         control = (control_small - control_large)/eps*norm(v)
     return control      

 def test_no_window(self) :
     window = sp.ones_like(self.wave1)
     power = np.windowed_power(self.wave1, window)
     self.assertAlmostEqual(power[self.mode1]/self.amp1**2/self.n*4, 1)
     self.assertTrue(sp.allclose(power[:self.mode1], 0, 
                                 atol=self.amp1**2*self.n/1e15))
     self.assertTrue(sp.allclose(power[self.mode1+1:], 0, 
                                 atol=self.amp1**2*self.n/1e15))
     # With no window, we have a quick way to the answer
     quick_power = (abs(fft.fft(self.wave1))**2)[:self.n//2]/self.n
     self.assertTrue(sp.allclose(power, quick_power))