Python scipy.isfinite函数代码示例

def normalize(series):
    '''
    Returns the series demeaned and divided by its standard deviation.
    '''
    mean = series[sp.isfinite(series)].mean()
    sdev = series[sp.isfinite(series)].std()
    return (series - mean) / sdev

def obs_d22(station, day, numdays=1, varlist=[]):
    '''
    station: offshore station name as string
    day: datetime object (should be with 00 hours)
    returns dictionary with hourly time, Hs, Tp, Tm, FF, DD
    each of the returned variable is a 2d array with time and the number of available sensors
    '''
    time, WMlist, WIlist = d22.read_d22(station,start=day, end=day+dt.timedelta(numdays-1))
    datadict = {'time':time['1hr']}

    # organize each variable as 2d-arrays with time and sensor number as dimensions
    Hs, Tp, Tm02 ,Tm01, DDP, DDM, WMnames, FF, DD, WInames = [],[],[],[], [], [], [], [], [], []
    for WM in WMlist:
        if sp.isfinite(WM['Hs']).sum()>1:
            Hs.append(WM['Hs'])
            Tp.append(WM['Tp'])
            Tm02.append(WM['Tm02'])
            Tm01.append(WM['Tm01'])
            DDP.append(WM['DDP'])
            DDM.append(WM['DDM'])
            WMnames.append(WM['name'])

    for WI in WIlist:
        if sp.isfinite(WI['FF']).sum()>1:
            FF.append(WI['FF'])
            DD.append(WI['DD'])
            WInames.append(WI['name'])

    datadict.update({'Hs':sp.array(Hs), 'Tp':sp.array(Tp), 'Tm02':sp.array(Tm02) , 'Tm01':sp.array(Tm01) , 
                     'FF':sp.array(FF), 'DD':sp.array(DD),'DDM':sp.array(DDM), 'DDP':sp.array(DDP)})
    datadict.update({'WMnames':WMnames, 'WInames':WInames}) # save sensor names

    return datadict

def bias(a,b):
    '''
    bias
    '''
    a,b = sp.array(a),sp.array(b)
    mask = sp.logical_and(sp.isfinite(a),sp.isfinite(b))
    a, b = a[mask], b[mask]
    return a.mean()-b.mean()

    def intercept(self, ray):
        """Solves for intersection point of surface and a ray or Beam
    
        Args:
            ray: Ray or Beam object
                It must be in the same coordinate space as the surface object.
            
        Returns:
            s: value of s [meters] which intercepts along norm, otherwise an
            empty tuple (for no intersection).
        
        Examples:
            Accepts all point and point-derived object inputs, though all data 
            is stored as a python object.

            Generate an y direction Ray in cartesian coords using a Vec from (0,0,1)::
            
                    cen = geometry.Center(flag=True)
                    ydir = geometry.Vecx((0,1,0))
                    zpt = geometry.Point((0,0,1),cen)

        """


        # Proceedure will be to generate 
        if self._origin is ray._origin:
            try:
                rcopy = ray.copy()
                rcopy.redefine(self)
                
                intersect = _beam.interceptCyl(scipy.atleast_2d(rcopy.x()[:,-1]), 
                                               scipy.atleast_2d(rcopy.norm.unit), 
                                               scipy.array([self.sagi.s,self.sagi.s]),
                                               scipy.array([-self.norm.s,self.norm.s])) + rcopy.norm.s[-1]
                
                if not scipy.isfinite(intersect):
                    #relies on r1 using arctan2 so that it sets the branch cut properly (-pi,pi]
                    return None
                elif self.edgetest(intersect, (rcopy(intersect)).r1()):
                        return intersect
                else:
                    rcopy.norm.s[-1] = intersect
                    intersect = _beam.interceptCyl(scipy.atleast_2d(rcopy.x()[:,-1]), 
                                                   scipy.atleast_2d(rcopy.norm.unit), 
                                                   scipy.array([self.sagi.s,self.sagi.s]),
                                                   scipy.array([-self.norm.s,self.norm.s])) + rcopy.norm.s[-1]
                    if not scipy.isfinite(intersect):
                        #relies on r1 using arctan2 so that it sets the branch cut properly (-pi,pi]
                        return None
                    elif self.edgetest(intersect, (rcopy(intersect)).r1()):
                        return None
                    else:
                        return None

            except AttributeError:
                raise ValueError('not a surface object')
        else:           
            raise ValueError('not in same coordinate system, use redefine and try again')

def d22mean(sample):
    ''' return the mean if at least 3 of 6  10min values are finite
     Returns NaN only if more than 3 values missing, or if all are the same (std=0)
     '''
    if (sum(sp.isfinite(sample)) > 2 and sp.std(sample) > 0.):
        average = sp.mean(sample[sp.isfinite(sample)])
    else:
        average = sp.nan
    return average

 def df(x):
     x_ = X0
     x_[Ifilter_x] = x
     rv = gpr.LMLgrad(param_list_to_dict(x_,param_struct,skeys),*args,**kw_args)
     rv = param_dict_to_list(rv,skeys)
     if (~SP.isfinite(rv)).any():
         idx = (~SP.isfinite(rv))
         rv[idx] = 1E6
     return rv[Ifilter_x]

 def df(x):
     x_ = X0
     x_[Ifilter_x] = x
     rv = gpr.LMLgrad(param_list_to_dict(x_, param_struct, skeys), *args, **kw_args)
     rv = param_dict_to_list(rv, skeys)
     # LG.debug("dL("+str(x_)+")=="+str(rv))
     if not SP.isfinite(rv).all():  # SP.isnan(rv).any():
         In = ~SP.isfinite(rv)#SP.isnan(rv)
         rv[In] = 1E6
     return rv[Ifilter_x]

 def test_remove_boundary_conditions(self):
     alg = op.algorithms.GenericTransport(network=self.net,
                                          phase=self.phase)
     alg.set_value_BC(pores=self.net.pores('top'), values=1)
     alg.set_value_BC(pores=self.net.pores('bottom'), values=0)
     assert sp.sum(sp.isfinite(alg['pore.bc_value'])) > 0
     alg.remove_BC(pores=self.net.pores('top'))
     assert sp.sum(sp.isfinite(alg['pore.bc_value'])) > 0
     alg.remove_BC(pores=self.net.pores('bottom'))
     assert sp.sum(sp.isfinite(alg['pore.bc_value'])) == 0

 def f(x):
     x_ = X0
     x_[Ifilter_x] = x
     lml = gpr.LML(param_list_to_dict(x_,param_struct,skeys))
     if SP.isnan(lml):
         lml=1E6
     lml_grad = gpr.LMLgrad()
     lml_grad = param_dict_to_list(lml_grad,skeys)
     if (~SP.isfinite(lml_grad)).any():
         idx = (~SP.isfinite(lml_grad))
         lml_grad[idx] = 1E6
     return lml, lml_grad[Ifilter_x]

def lognpdf(x,mean,sig):
    if x<0 or not scipy.isfinite(x):
        pdf = 0
    else:
        a   = 1./(x*sig*scipy.sqrt(2*scipy.pi))
        pdf =  a*scipy.exp(-(scipy.log(x)-mean)**2/(2.*sig**2)) 
    return pdf

 def set(self, arg):
     if type(arg) in [list, tuple, Array, NArray] and len(arg)==2:
         if arg[0] == arg[1]:
             self.set(arg[0])
         else:
             self.issingleton = False
             loval = arg[0]
             hival = arg[1]
             assert loval is not NaN and hival is not NaN, \
                    ("Cannot specify NaN as interval endpoint")
             if not loval < hival:
                 print "set() was passed loval = ", loval, \
                       " and hival = ", hival
                 raise PyDSTool_ValueError, ('Interval endpoints must be '
                                 'given in order of increasing size')
             self._intervalstr = '['+str(loval)+',' \
                                 +str(hival)+']'
             self._loval = loval
             self._hival = hival            
             self.defined = True
     elif type(arg) in [int, float]:
         assert isfinite(arg), \
                "Singleton interval domain value must be finite"
         self.issingleton = True
         self._intervalstr = str(arg)
         self._loval = arg
         self._hival = arg
         self.defined = True
     else:
         print "Error in argument: ", arg, "of type", type(arg)
         raise PyDSTool_TypeError, \
               'Interval spec must be a numeric or a length-2 sequence type'

    def atEndPoint(self, val, bdcode):
        """val, bdcode -> Bool

        Determines whether val is at the endpoint specified by bdcode,
        to the precision of the interval's _abseps tolerance.
        bdcode can be one of 'lo', 'low', 0, 'hi', 'high', 1"""

        assert self.defined, 'Interval undefined'
        assert isinstance(val, float) or isinstance(val, int), \
               'Invalid value type'
        assert isfinite(val), "Can only test finite argument values"
        if bdcode in ['lo', 'low', 0]:
            if self.type == Float:
                return abs(val - self._loval) < self._abseps
            elif self.type == Int:
                return val == self._loval
            else:
                raise TypeError, "Unsupported value type"
        elif bdcode in ['hi', 'high', 1]:
            if self.type == Float:
                return abs(val - self._hival) < self._abseps
            elif self.type == Int:
                return val == self._hival
            else:
                raise TypeError, "Unsupported value type"
        else:
            raise ValueError, 'Invalid boundary spec code'

def lnprob(theta, time, rv, err):
    lp = lnprior(theta)
    if not sp.isfinite(lp):
        return -sp.inf
    #print(lp)
    #print(lp + lnlike(theta, time, rv, err))
    return lp + lnlike(theta, time, rv, err)

def TransitPhase(tset):
    lc = copy.deepcopy(tset.tables[1])
    time = lc.TIME
    nobs = len(time)
    lg = scipy.isfinite(time)
    pl = tset.tables[0]
    npl = len(pl.Period)
    phase = scipy.zeros((npl, nobs))
    inTr = scipy.zeros((npl, nobs), "int")
    for ipl in scipy.arange(npl):
        period = pl.Period[ipl]
        t0 = pl.t0[ipl] + BJDREF_t0 - BJDREF_lc
        dur = pl.Dur[ipl] / 24.0 / period
        counter = 0
        while (time[lg] - t0).min() < 0:
            t0 -= period
            counter += 1
            if counter > 1000:
                break
        ph = ((time - t0) % period) / period
        ph[ph < -0.5] += 1
        ph[ph > 0.5] -= 1
        phase[ipl, :] = ph
        inTr[ipl, :] = abs(ph) <= dur / 1.5
    return phase, inTr

def PlotLc(id=None, dir=None, quarter=None, tset=None):
    if id != None:
        tset, status = GetLc(id=id, dir=dir, tr_out=True)
        if tset is None:
            return
    elif tset == None:
        print "no tset"
        return
    if quarter != None:
        tset.tables[1] = tset.tables[1].where(tset.tables[1].Q == quarter)
        if len(tset.tables[1].TIME) == 0:
            print "No data for Q%d" % quarter
            return
    time = tset.tables[1].TIME
    phase, inTr = TransitPhase(tset)
    col = ["r", "g", "b", "y", "c", "m", "grey"]
    npl, nobs = inTr.shape
    pylab.figure(1)
    pylab.clf()
    pylab.plot(time, tset.tables[1].PDCSAP_FLUX, "k-")
    for ipl in scipy.arange(npl):
        list = inTr[ipl, :].astype(bool)
        pylab.plot(time[list], tset.tables[1].PDCSAP_FLUX[list], ".", c=col[ipl])
    l = scipy.isfinite(time)
    pylab.xlim(time[l].min(), time[l].max())
    ttl = "KIC %d, P=" % (tset.tables[0].KID[0])
    for i in scipy.arange(npl):
        ttl = "%s%.5f " % (ttl, tset.tables[0].Period[i])
    if quarter != None:
        ttl += "Q%d" % quarter
    pylab.title(ttl)

    return