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

 def most_normal_transformation(self,trans_types=SUPPORTED_TRANSFORMATIONS,
             perform_trans=True, verbose=False):
     """
     Performs the transformation which results in most normal looking data, according to Shapiro-Wilk's test
     """
     from scipy import stats
     shapiro_pvals = []
     for trans_type in trans_types:
         if trans_type == 'most_normal':
             continue
         if trans_type != 'none':
             if not self.transform(trans_type=trans_type):
                 continue
         phen_vals = self.values
         #print 'sp.inf in phen_vals:', sp.inf in phen_vals
         if sp.inf in phen_vals:
             pval = 0.0
         else:
             r = stats.shapiro(phen_vals)
             if sp.isfinite(r[0]):
                 pval = r[1]
             else:
                 pval = 0.0
         shapiro_pvals.append(pval)
         if trans_type != 'none':
             self.revert_to_raw_values()
     argmin_i = sp.argmax(shapiro_pvals)
     trans_type = trans_types[argmin_i]
     shapiro_pval = shapiro_pvals[argmin_i]
     if perform_trans:
         self.transform(trans_type=trans_type)
     log.info("The most normal-looking transformation was %s, with a Shapiro-Wilk's p-value of %.2E" % \
             (trans_type, shapiro_pval))
     return trans_type, shapiro_pval

	def leap_prob_recurse(self, Z_chain, C, active_idx):
		"""
		Recursively compute to cumulative probability of transitioning from
		the beginning of the chain Z_chain to the end of the chain Z_chain.
		"""

		if sp.isfinite(C[0,-1,0]):
			## we've already visited this leaf
			cumu = C[0,-1,:].reshape((1,-1))
			return cumu, C

		if len(Z_chain) == 2:
			## the two states are one apart
			p_acc = self.leap_prob(Z_chain[0], Z_chain[1])
			p_acc = p_acc[:,active_idx]
			C[0,-1,:] = p_acc.ravel()
			return p_acc, C

		cum_forward, Cl = self.leap_prob_recurse(Z_chain[:-1], C[:-1,:-1,:], active_idx)
		C[:-1,:-1,:] = Cl
		cum_reverse, Cl = self.leap_prob_recurse(Z_chain[:0:-1], C[:0:-1,:0:-1,:], active_idx)
		C[:0:-1,:0:-1,:] = Cl

		H0 = self.Eval_H(Z_chain[0])
		H1 = self.Eval_H(Z_chain[-1])
		Ediff = H0 - H1
		Ediff = Ediff[:,active_idx]
		start_state_ratio = sp.exp(Ediff)
		prob =((sp.vstack  ((1. - cum_forward, start_state_ratio*(1. - cum_reverse)))).min(axis=0)).reshape((1,-1))
		cumu = cum_forward + prob
		C[0,-1,:] = cumu.ravel()
		return cumu, C

def EBTransitPhase(tset, kid_x):
    ebp = atpy.Table('%s/eb_pars.txt' %dir, type = 'ascii')
    
    lc = tset.tables[1]   
    time = lc.TIME
    flux = lc.PDCSAP_FLUX
    nobs = len(time)
    lg = scipy.isfinite(time)
    pylab.figure(52)
    pylab.clf()
    pylab.plot(time[lg], flux[lg])
    npl = 2
    phase = scipy.zeros((npl, nobs))
    inTr = scipy.zeros((npl, nobs), 'int')
    period = ebp.P[ebp.KID == kid_x]
    for ipl in scipy.arange(npl):
        if ipl == 0: t0 = ebp.Ep1[ebp.KID == kid_x]
        if ipl == 1: t0 = ebp.Ep2[ebp.KID == kid_x]
        if ipl == 0: dur = ebp.Dur1[ebp.KID == kid_x]
        if ipl == 1: dur = ebp.Dur2[ebp.KID == kid_x]
        dur /= 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 _box_cox_transform(self, verbose=False, method='standard'):
     """
     Performs the Box-Cox transformation, over different ranges, picking the optimal one w. respect to normality.
     """
     from scipy import stats
     a = sp.array(self.values)
     if method == 'standard':
         vals = (a - min(a)) + 0.1 * sp.var(a)
     else:
         vals = a
     sw_pvals = []
     lambdas = sp.arange(-2.0, 2.1, 0.1)
     for l in lambdas:
         if l == 0:
             vs = sp.log(vals)
         else:
             vs = ((vals ** l) - 1) / l
         r = stats.shapiro(vs)
         if sp.isfinite(r[0]):
             pval = r[1]
         else:
             pval = 0.0
         sw_pvals.append(pval)
     i = sp.argmax(sw_pvals)
     l = lambdas[i]
     if l == 0:
         vs = sp.log(vals)
     else:
         vs = ((vals ** l) - 1) / l
     self._perform_transform(vs,"box_cox")
     log.debug('optimal lambda was %0.1f' % l)
     return True

 def most_normal_transformation(self, pid, trans_types=['none', 'sqrt', 'log', 'sqr', 'exp', 'arcsin_sqrt'],
             perform_trans=True, verbose=False):
     """
     Performs the transformation which results in most normal looking data, according to Shapiro-Wilk's test
     """
     #raw_values = self.phen_dict[pid]['values']
     from scipy import stats
     shapiro_pvals = []
     for trans_type in trans_types:
         if trans_type != 'none':
             if not self.transform(pid, trans_type=trans_type):
                 continue
         phen_vals = self.get_values(pid)
         #print 'sp.inf in phen_vals:', sp.inf in phen_vals
         if sp.inf in phen_vals:
             pval = 0.0
         else:
             r = stats.shapiro(phen_vals)
             if sp.isfinite(r[0]):
                 pval = r[1]
             else:
                 pval = 0.0
         shapiro_pvals.append(pval)
         #self.phen_dict[pid]['values'] = raw_values
         if trans_type != 'none':
             self.revert_to_raw_values(pid)
     argmin_i = sp.argmax(shapiro_pvals)
     trans_type = trans_types[argmin_i]
     shapiro_pval = shapiro_pvals[argmin_i]
     if perform_trans:
         self.transform(pid, trans_type=trans_type)
     if verbose:
         print "The most normal-looking transformation was %s, with a Shapiro-Wilk's p-value of %0.6f" % \
             (trans_type, shapiro_pval)
     return trans_type, shapiro_pval