def trainNB0(trainMatrix, trainClassMatrix):
    '''
    朴素贝叶斯分类器训练函数
    :param trainMatrix: 训练集矩阵
    :param trainClassMatrix: 包含训练集的类别矩阵
    :return: p0Vect: 第0类的一个向量，其中每个维度的值表示在第0类的所有样本的所有词组中，该维度的词所占的比例, p1Vect: 意思和p0Vect类同, pAbusive: 第1类样本的数量在训练集所占的比例
    '''
    # 获取训练集有多少个样本
    numOfTrainDocs = len(trainMatrix)
    # 获取词条向量所包含的属性个数
    numOfWords = len(trainMatrix[0])
    pAbusive = sum(trainClassMatrix) / numOfTrainDocs
    p0Num = ones(numOfWords)
    p1Num = ones(numOfWords)
    # 分母，表示在第0类下训练集所包含的词的总数量
    p0Denom = 2.0
    p1Denom = 2.0
    for i in range(numOfTrainDocs):
        if trainClassMatrix[i] == 1:
            p1Num += trainMatrix[i]
            p1Denom += sum(trainMatrix[i])
        else:
            p0Num += trainMatrix[i]
            p0Denom += sum(trainMatrix[i])
    # 这里用log是防止下溢
    p0Vect = log(p0Num / p0Denom)
    p1Vect = log(p1Num / p1Denom)
    return p0Vect, p1Vect, pAbusive

def classifyNB(inputVec, p0Vect, p1Vect, pAbsive):
    p0 = sum(inputVec * p0Vect) + log(1 - pAbsive)
    p1 = sum(inputVec * p1Vect) + log(pAbsive)
    if p1 > p0:
        return 1
    else:
        return 0

def costFunction_Regular(theta, *args):
    print "现在在调用正则化的cost函数"
    # print "@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@",shape(theta)
    dataArr = asarray(args[0])  # args0 为特征数据  args1 为类别数据
    labelArr = asarray(args[1])
    m, n = shape(dataArr)  # dataArr 已经加入bias
    # print "记录条数:",m
    theta = reshape(theta, (n, 1))
    #print "costFunction_Regular中theta的值为：", theta
    hx = sigmodFunction(dataArr, theta)  #计算预测的类别概率值
    loghx = ma.log(hx)
    #print "log hx:",loghx
    yhx = dot(loghx.transpose(), labelArr)
    #print "yhx:",yhx
    log1_hx = ma.log(1 - hx)
    #print "log1_hx:",log1_hx
    #print (-yhx-dot(log1_hx.transpose(),(1-labelArr)))*1.0/m
    #print "lameda 值为:", args[2]  # args[2]为传入的lameda参数
    jtheta = (-yhx - dot(log1_hx.transpose(), (1 - labelArr))) * 1.0 / m + args[2] * 1.0 / (2 * m) * (
        dot(theta.transpose(), theta) - theta[0, 0] ** 2)
    #gra=getGradient(dataArr,labelArr,theta,m)
    #print type(jtheta),type(gra.flatten())
    #print gra
    #print "###################",type(array(jtheta)[0]),array(jtheta)[0]
    #print "costFunction_Regular得到的jtheta值为:", type(jtheta.flatten()[0]), jtheta, jtheta.flatten()[0]
    print "&&&&&7jtheta.flatten()[0]", jtheta.flatten()[0]
    return jtheta.flatten()[0]

def entropy(array, dim=None):
    if dim is None:
        array = array.ravel()
        dim = 0
    n = ma.sum(array, dim)
    array = ma.log(array) * array
    sum = ma.sum(array, dim)
    return (ma.log(n) - sum / n) / ma.log(2.0)

def average_in_flux(mag, dmag, axis=None):
    flux = 10**(mag / -2.5)
    dflux = np.log(10) / 2.5 * flux * dmag
    avg_dflux = np.power(np.sum(np.power(dflux, -2), axis), -0.5)
    avg_flux = np.sum(flux * np.power(dflux, -2), axis) * avg_dflux**2
    avg_mag = -2.5 * np.log10(avg_flux)
    avg_dmag = 2.5 / np.log(10) * np.divide(avg_dflux, avg_flux)
    return avg_mag, avg_dmag

def computeCost(theta, X, y, lamda):
    m = np.shape(X)[0]
    hypo = sigmoid(X.dot(theta))
    term1 = log(hypo).dot(-y)
    term2 = log(1.0 - hypo).dot(1 - y)
    left_hand = (term1 - term2) / m
    right_hand = theta.transpose().dot(theta) * lamda / (2 * m)
    return left_hand + right_hand

def log_linear_vinterp(T,P,levs):
    '''
    # Author Charles Doutriaux
    # Version 1.1
    # Expect 2D field here so there''s no reorder which I suspect to do a memory leak
    # email: [email protected]
    # Converts a field from sigma levels to pressure levels
    # Log linear interpolation


    # Input
    # T :    temperature on sigma levels
    # P :    pressure field from TOP (level 0) to BOTTOM (last level)
    # levs : pressure levels to interplate to (same units as P)

    # Output
    # t :    temperature on pressure levels (levs)

    # External: Numeric'''
    import numpy.ma as MA
##     from numpy.oldnumeric.ma import ones,Float,greater,less,logical_and,where,equal,log,asarray,Float16
    sh=P.shape
    nsigma=sh[0] # Number of sigma levels
    try:
        nlev=len(levs)  # Number of pressure levels
    except:
        nlev=1  # if only one level len(levs) would breaks
    t=[]
    for ilv in range(nlev): # loop through pressure levels
        try:
            lev=levs[ilv] # get value for the level
        except:
            lev=levs  # only 1 level passed
#       print '          ......... level:',lev
        Pabv=MA.ones(P[0].shape,Numeric.Float)
        Tabv=-Pabv # Temperature on sigma level Above
        Tbel=-Pabv # Temperature on sigma level Below
        Pbel=-Pabv # Pressure on sigma level Below
        Pabv=-Pabv # Pressure on sigma level Above
        for isg in range(1,nsigma): # loop from second sigma level to last one
##             print 'Sigma level #',isg
            a = MA.greater(P[isg],  lev) # Where is the pressure greater than lev
            b = MA.less(P[isg-1],lev)    # Where is the pressure less than lev

            # Now looks if the pressure level is in between the 2 sigma levels
            # If yes, sets Pabv, Pbel and Tabv, Tbel
            Pabv=MA.where(MA.logical_and(a,b),P[isg],Pabv) # Pressure on sigma level Above
            Tabv=MA.where(MA.logical_and(a,b),T[isg],Tabv) # Temperature on sigma level Above
            Pbel=MA.where(MA.logical_and(a,b),P[isg-1],Pbel) # Pressure on sigma level Below
            Tbel=MA.where(MA.logical_and(a,b),T[isg-1],Tbel) # Temperature on sigma level Below
        # end of for isg in range(1,nsigma)
#       val=where(equal(Pbel,-1.),Pbel.missing_value,lev) # set to missing value if no data below lev if there is
        
        tl=MA.masked_where(MA.equal(Pbel,-1.),MA.log(lev/MA.absolute(Pbel))/MA.log(Pabv/Pbel)*(Tabv-Tbel)+Tbel) # Interpolation
        t.append(tl) # add a level to the output
    # end of for ilv in range(nlev)
    return asMA(t).astype(Numeric.Float32) # convert t to an array

def KL_Measure(i, j):
    '''
    计算KL散度
    :return:
    '''
    KL1 = sum(i*(log(i/j).data))
    KL2 = sum(j*(log(j/i).data))
    D = (KL1 + KL2)/2
    return 1/(1+ math.e ** D )

def _zfromp_MA(P, lapse_rate, P_bott, T_bott, z_bott):
    """Altitude given pressure in a constant lapse rate layer.

    The dry gas constant is used in calculations requiring the gas
    constant.  See the docstring for press2alt for references.

    Input Arguments:
    * P:  Pressure [hPa].
    * lapse_rate:  -dT/dz [K/m] over the layer.
    * P_bott:  Pressure [hPa] at the base of the layer.
    * T_bott:  Temperature [K] at the base of the layer.
    * z_bott:  Geopotential altitude [m] of the base of the layer.

    Output:
    * Altitude [m] for each element given in the input arguments.

    All input arguments can be either a scalar or an MA array.  All 
    arguments that are MA arrays, however, are of the same size and 
    shape.  If every input argument is a scalar, the output is a scalar.
    If any of the input arguments is an MA array, the output is an MA 
    array of the same size and shape.
    """
    import numpy as N
    #jfp was import Numeric as N
    import numpy.ma as MA
    #jfp was import MA
    from atmconst import AtmConst

    const = AtmConst()

    if MA.size(lapse_rate) == 1:
        if MA.array(lapse_rate)[0] == 0.0:
            return ( (-const.R_d * T_bott / const.g) * MA.log(P/P_bott) ) + \
                   z_bott
        else:
            exponent = (const.R_d * lapse_rate) / const.g
            return ((T_bott / lapse_rate) * (1. - (P/P_bott)**exponent)) + \
                   z_bott
    else:
        exponent = (const.R_d * lapse_rate) / const.g
        z = ((T_bott / lapse_rate) * (1. - (P/P_bott)**exponent)) + z_bott
        z_at_0 = ( (-const.R_d * T_bott / const.g) * MA.log(P/P_bott) ) + \
                 z_bott

        zero_lapse_mask = MA.filled(MA.where(lapse_rate == 0., 1, 0), 0)
        zero_lapse_mask_indices_flat = N.nonzero(N.ravel(zero_lapse_mask))
        z_flat = MA.ravel(z)
        MA.put( z_flat, zero_lapse_mask_indices_flat \
              , MA.take(MA.ravel(z_at_0), zero_lapse_mask_indices_flat) )
        return MA.reshape(z_flat, z.shape)

def classifyNB(vecToClassify, p0Vec, p1Vec, pClass1):
    '''
    利用训练好的模型进行分类
    :param param:
    :param p0Vec:
    :param p1Vec:
    :param pClass1:
    :return:
    '''
    p1 = sum(vecToClassify * p0Vec) + log(pClass1)
    p0 = sum(vecToClassify * p1Vec) + log(1 - pClass1)
    if p1 > p0:
        return 1
    else:
        return 0

    def __call__(self, value, clip=None):
        if clip is None:
            clip = self.clip

        if cbook.iterable(value):
            vtype = 'array'
            val = ma.asarray(value).astype(np.float)
        else:
            vtype = 'scalar'
            val = ma.array([value]).astype(np.float)

        self.autoscale_None(val)
        vmin, vmax = self.vmin, self.vmax
        if vmin > vmax:
            raise ValueError("minvalue must be less than or equal to maxvalue")
        elif vmin<=0:
            raise ValueError("values must all be positive")
        elif vmin==vmax:
            return 0.0 * val
        else:
            if clip:
                mask = ma.getmask(val)
                val = ma.array(np.clip(val.filled(vmax), vmin, vmax),
                                mask=mask)
            result = (ma.log(val)-np.log(vmin))/(np.log(vmax)-np.log(vmin))
        if vtype == 'scalar':
            result = result[0]
        return result

 def test_testUfuncs1(self):
     # Test various functions such as sin, cos.
     (x, y, a10, m1, m2, xm, ym, z, zm, xf, s) = self.d
     assert_(eq(np.cos(x), cos(xm)))
     assert_(eq(np.cosh(x), cosh(xm)))
     assert_(eq(np.sin(x), sin(xm)))
     assert_(eq(np.sinh(x), sinh(xm)))
     assert_(eq(np.tan(x), tan(xm)))
     assert_(eq(np.tanh(x), tanh(xm)))
     with np.errstate(divide='ignore', invalid='ignore'):
         assert_(eq(np.sqrt(abs(x)), sqrt(xm)))
         assert_(eq(np.log(abs(x)), log(xm)))
         assert_(eq(np.log10(abs(x)), log10(xm)))
     assert_(eq(np.exp(x), exp(xm)))
     assert_(eq(np.arcsin(z), arcsin(zm)))
     assert_(eq(np.arccos(z), arccos(zm)))
     assert_(eq(np.arctan(z), arctan(zm)))
     assert_(eq(np.arctan2(x, y), arctan2(xm, ym)))
     assert_(eq(np.absolute(x), absolute(xm)))
     assert_(eq(np.equal(x, y), equal(xm, ym)))
     assert_(eq(np.not_equal(x, y), not_equal(xm, ym)))
     assert_(eq(np.less(x, y), less(xm, ym)))
     assert_(eq(np.greater(x, y), greater(xm, ym)))
     assert_(eq(np.less_equal(x, y), less_equal(xm, ym)))
     assert_(eq(np.greater_equal(x, y), greater_equal(xm, ym)))
     assert_(eq(np.conjugate(x), conjugate(xm)))
     assert_(eq(np.concatenate((x, y)), concatenate((xm, ym))))
     assert_(eq(np.concatenate((x, y)), concatenate((x, y))))
     assert_(eq(np.concatenate((x, y)), concatenate((xm, y))))
     assert_(eq(np.concatenate((x, y, x)), concatenate((x, ym, x))))

def dewpoint(e):
    r'''Calculate the ambient dewpoint given the vapor pressure.

    Parameters
    ----------
    e : array_like
        Water vapor partial pressure in mb

    Returns
    -------
    array_like
        Dew point temperature in degrees Celsius.

    See Also
    --------
    dewpoint_rh, saturation_vapor_pressure, vapor_pressure

    Notes
    -----
    This function inverts the Bolton 1980 [3] formula for saturation vapor
    pressure to instead calculate the temperature. This yield the following
    formula for dewpoint in degrees Celsius:

    .. math:: T = \frac{243.5 log(e / 6.112)}{17.67 - log(e / 6.112)}

    References
    ----------
    .. [3] Bolton, D., 1980: The Computation of Equivalent Potential
           Temperature. Mon. Wea. Rev., 108, 1046-1053.
    '''

    val = log(e / sat_pressure_0c)
    return 243.5 * val / (17.67 - val)

def geometric_mean(array, axis=0):
    '''return the geometric mean of an array removing all zero-values but
    retaining total length
    '''
    non_zero = ma.masked_values(array, 0)
    log_a = ma.log(non_zero)
    return ma.exp(log_a.mean(axis=axis))

 def transform(self, a):
     sign = np.sign(a)
     masked = ma.masked_inside(a, -self.linthresh, self.linthresh, copy=False)
     log = sign * self.linthresh * (1 + ma.log(np.abs(masked) / self.linthresh))
     if masked.mask.any():
         return ma.where(masked.mask, a, log)
     else:
         return log

 def transform_non_affine(self, a):
     sign = np.sign(a)
     masked = ma.masked_inside(a, -self.linthresh, self.linthresh, copy=False)
     log = sign * self.linthresh * (self._linscale_adj + ma.log(np.abs(masked) / self.linthresh) / self._log_base)
     if masked.mask.any():
         return ma.where(masked.mask, a * self._linscale_adj, log)
     else:
         return log

 def transform(self, a):
     sign = np.sign(np.asarray(a))
     masked = ma.masked_inside(a, -self.linthresh, self.linthresh, copy=False)
     log = sign * ma.log(np.abs(masked)) / self._log_base
     if masked.mask.any():
         return np.asarray(ma.where(masked.mask, a * self._linadjust, log))
     else:
         return np.asarray(log)

def ln_shifted_auto(v):
    """If 'v' has values <= 0, it is shifted in a way that min(v)=1 before doing log. 
    Otherwise the log is done on the original 'v'."""
    vmin = ma.minimum(v)
    if vmin <= 0:
        values = v - vmin + 1
    else:
        values = v
    return ma.log(values)

def trainNB1(trainMatrix, trainCategory):
    numTrainDocs = len(trainMatrix)
    numWord = len(trainMatrix[0])
    pAbusive = sum(trainCategory) / float(numTrainDocs)
    p0Num = ones(numWord)
    p1Num = ones(numWord)
    p0Denom = 2.0
    p1Denom = 2.0
    for i in range(numTrainDocs):
        if trainCategory[i] == 1:
            p1Num += trainMatrix[i]
            p1Denom += sum(trainMatrix[i])
        else:
            p0Num += trainMatrix[i]
            p0Denom += sum(trainMatrix[i])
    p0Vect = log(p0Num / p0Denom)
    p1Vect = log(p1Num / p1Denom)
    return p0Vect, p1Vect, pAbusive

 def transform(self, a):
     a = np.asarray(a)
     sign = np.sign(a)
     masked = ma.masked_inside(a, -self.linthresh, self.linthresh, copy=False)
     if masked.mask.any():
         log = sign * (ma.log(np.abs(masked)) / self._log_base + self._linadjust)
         return np.asarray(ma.where(masked.mask, a * self._linscale, log))
     else:
         return sign * (np.log(np.abs(a)) / self._log_base + self._linadjust)