def _sanitize_SchedCpuCapacity(self):
        """
        Add more columns to cpu_capacity data frame if the energy model is
        available and the platform is big.LITTLE.
        """
        if not self.hasEvents('cpu_capacity') \
           or 'nrg_model' not in self.platform \
           or not self.has_big_little:
            return

        df = self._dfg_trace_event('cpu_capacity')

        # Add column with LITTLE and big CPUs max capacities
        nrg_model = self.platform['nrg_model']
        max_lcap = nrg_model['little']['cpu']['cap_max']
        max_bcap = nrg_model['big']['cpu']['cap_max']
        df['max_capacity'] = np.select(
                [df.cpu.isin(self.platform['clusters']['little'])],
                [max_lcap], max_bcap)
        # Add LITTLE and big CPUs "tipping point" threshold
        tip_lcap = 0.8 * max_lcap
        tip_bcap = 0.8 * max_bcap
        df['tip_capacity'] = np.select(
                [df.cpu.isin(self.platform['clusters']['little'])],
                [tip_lcap], tip_bcap)

def soft_hard_burst(playercard):
    ## cards have 1
    if sum(np.select([playercard == 1],[playercard]))>0:
        criterion = sum(np.select([playercard == 1],[playercard])) - 1 + sum(np.select([playercard != 1],[playercard]))
        if criterion < 11:
            string="Soft"
            count =make_score(playercard)
            #count = playercard[0]+playercard[1]+10
            string = string + str(count)
            if count < 22:
                return string
            else:
                return "burst"
        elif criterion >10:
            string = "Hard"
            count =make_score(playercard)
            string = string + str(count)
            if count < 22:
                return string
            else:
                return "burst"
    ## cards have no 1
    elif sum(np.select([playercard == 1],[playercard]))==0:
        string="Hard"
        count = sum(playercard)
        string = string+str(count)
        if count < 22:
            return string
        else:
            return "burst"

def zpeaki(source,order=1,fpeak=fhigh):
    '''
        寻找n阶高/低点
        返回值为高点数据序列，以及该高点最大跨度的坐标(即计算该高/低点所需用到的最远的未来数据的坐标)
        order默认为1,小于1当作1
        返回值中第一个是高/低点非0,其余为0的序列 sh
                第二个是该高低点的最远未来数据的坐标序列 si
                其中 sh[np.nonzero(sh)]为高点序列, si[np.nonzero(sh)]为坐标序列,sif.time[si[np.nonzero(sh)]]为坐标的影响时间序列
    '''
    tsx1 = fpeak(source)
    sx1 = np.select([tsx1!=0],[source],0)
    icovered = rollx(np.arange(len(source)),-1)
    if order <= 1:
        return sx1,np.select([tsx1],[icovered],0)
    icursx = np.nonzero(tsx1)[0]
    for i in xrange(1,order):   #必然进入循环
        sxx = source[icursx]
        tsxx = fpeak(sxx)
        icovered[icursx] = rollx(icovered[icursx],-1)   #当前高/低点的计算范围,即之前顶点的范围左转一位(排除掉不是顶点的)
        icursx = icursx[np.nonzero(tsxx)[0]]
    osx = np.zeros_like(source)
    osx[icursx] = source[icursx]
    iz = np.zeros_like(source)
    iz[icursx] = icovered[icursx]   #去掉icovered之中不必要的那些数字
    return osx,iz

    def bin_indexes(self, arr, edgemode=0) :

        indmin, indmax = self._set_limit_indexes(edgemode)

        if self._equalbins :
            factor = float(self._nbins)/(self._edges[-1]-self._edges[0])
            nbins1 = self._nbins-1
            nparr = (np.array(arr, dtype=self._vtype)-self._edges[0])*factor
            ind = np.array(np.floor(nparr), dtype=np.int32)
            return np.select((ind<0, ind>nbins1), (indmin, indmax), default=ind)

        else :
            conds = None
            if self._ascending :
                conds = np.array([arr<edge for edge in self.binedges()], dtype=np.bool)
            else :            
                conds = np.array([arr>edge for edge in self.binedges()], dtype=np.bool)

            inds1d = range(-1, self._nbins)
            inds1d[0] = indmin # re-define index for underflow
            inds = np.array(len(arr)*inds1d, dtype=np.int32)
            inds.shape = (len(arr),self._nbins+1)
            inds = inds.transpose()
            #print 'indmin, indmax = ', indmin, indmax
            #print 'XXX conds:\n', conds
            #print 'XXX inds:\n', inds
            return np.select(conds, inds, default=indmax)

def vec_dam_break(x, t, h0=1.0, h1=10.0):
    import math
    import numpy
    from anuga import g

    msg = "Argument x should be a numpy array"
    assert isinstance(x, numpy.ndarray), msg

    h2 = calc_h2(h0, h1)
    u2 = 2.0 * (math.sqrt(g * h1) - math.sqrt(g * h2))

    try:
        s = u2 * h2 / (h2 - h0)
    except ZeroDivisionError:
        s = math.sqrt(g * h2)

    c1 = math.sqrt(g * h1)
    c2 = math.sqrt(g * h2)

    condlist = [x < -t * c1, x < t * (u2 - c2), x < s * t]
    hchoicelist = [h1, 1.0 / g * (2.0 / 3.0 * c1 - 1.0 / 3.0 * x / t) ** 2, h2]
    uchoicelist = [0.0, 2.0 / 3.0 * (c1 + x / t), u2]

    h = numpy.select(condlist, hchoicelist, default=h0)
    u = numpy.select(condlist, uchoicelist, default=0.0)

    return h, u

    def __getitem__(self, key):
        # If the key is a string, just get the subnode
        if isinstance(key, str):
            return self.__getattr__(key)
        # If the key is a vector, e.g. ['zone_1', 'zone_2', 'zone_1']
        elif isinstance(key, np.ndarray):
            if not np.issubdtype(key.dtype, np.str_):
                # In case the key is not a string vector, stringify it
                if key.dtype == object and issubclass(type(key[0]), Enum):
                    enum = type(key[0])
                    key = np.select([key == item for item in enum], [item.name for item in enum])
                elif isinstance(key, EnumArray):
                    enum = key.possible_values
                    key = np.select([key == item.index for item in enum], [item.name for item in enum])
                else:
                    key = key.astype('str')
            names = list(self.dtype.names)  # Get all the names of the subnodes, e.g. ['zone_1', 'zone_2']
            default = np.full_like(self.vector[key[0]], np.nan)  # In case of unexpected key, we will set the corresponding value to NaN.
            conditions = [key == name for name in names]
            values = [self.vector[name] for name in names]
            result = np.select(conditions, values, default)
            if contains_nan(result):
                unexpected_key = set(key).difference(self.vector.dtype.names).pop()
                raise ParameterNotFound('.'.join([self._name, unexpected_key]), self._instant_str)

            # If the result is not a leaf, wrap the result in a vectorial node.
            if np.issubdtype(result.dtype, np.record):
                return VectorialParameterNodeAtInstant(self._name, result.view(np.recarray), self._instant_str)

            return result

def getLargVal_man(*inA):
    inputlen = len(inA) 
    
    if inputlen == 2:
        condlist = [ inA[0] > inA[1] ]
        choicelist  = [ inA[0] ]
        result = np.select(condlist, choicelist, inA[1])
    
    elif inputlen == 3:
        condlist = [ np.logical_and(inA[0]>inA[1],inA[0]>inA[2]), 
                                    inA[1]>inA[2] ]
        choicelist  = [ inA[0], inA[1] ]
        result = np.select(condlist, choicelist, inA[2])        
        
        
    elif inputlen == 4: 
        condlist = [ np.logical_and(inA[0]>inA[1], 
                     np.logical_and(inA[0]>inA[2], inA[0]>inA[3])),
                     np.logical_and(inA[1]>inA[2], inA[1]>inA[3]),
                     inA[2]>inA[3] ]
        choicelist  = [ inA[0], inA[1], inA[2] ]
        result = np.select(condlist, choicelist, inA[3])
    
    else:
        print("Only up to 4 arrays supported")
    
    
    return result

    def update_particles(self, tick):
        mtick = tick/1000.0

        norm = np.sqrt(self.particle_pos[:,0]**2 + self.particle_pos[:,1]**2)
        norm = np.select([norm==0], [0.0000001], default=norm)
        posx = self.particle_pos[:,0]/norm
        posy = self.particle_pos[:,1]/norm

        radial = np.array([posx, posy])
        tangential = np.array([-posy, posx])

        radial = np.swapaxes(radial, 0, 1)
        radial *= self.particle_rad
        tangential = np.swapaxes(tangential, 0, 1)
        tangential *= self.particle_tan

        self.particle_dir +=  (tangential + radial + self.particle_grav)*mtick

        self.particle_pos += self.particle_dir*mtick

        self.particle_life -= mtick

        if self.position_type == POSITION_FREE:
            tuple = np.array(self.origin)
            tmp = tuple - self.particle_start_pos
            self.particle_pos -= tmp

        self.particle_color += self.particle_delta_color*mtick
        self.particle_color[:,3] = np.select([self.particle_life[:,0] < 0], [0], default=self.particle_color[:,3])

def up_seller(stock,buy_signal,xstop=25,ret=50,**kwargs):
    '''
        如果买入日为阴线，则开盘卖出
        如果价格小于最近5日高点5%，则卖出
        xstop为根据买入价的止损
        ret为从高点向下的回退值
    '''
    t = stock.transaction

    #阴线处理
    sol = rollx(gand(buy_signal,t[CLOSE] < t[OPEN]),1)

    #从顶下落处理,前5天的收盘/开盘的高者和今天的开盘的高者 回落ret之后

    #hhret = gmax(rollx(tmax(gmax(t[OPEN],t[CLOSE]),5),1),t[OPEN])* (1000-ret)/1000
    hhret = gmax(rollx(tmax(t[HIGH],5),1),t[OPEN])* (1000-ret)/1000
   
    sdl = t[LOW] < hhret
 
    #止损处理2.5%
    stop_price = extend2next(rollx(stock.buyprice,1) * (1000-xstop)/1000)   #要求buyprice只有在buyer日才有数据,否则extend2next无意义
    stopl = t[LOW] < stop_price

    cut_price = gmin(gmax(hhret,stop_price),t[HIGH])    #首先，止损线和退回线高者先被触及，同时，穿越时可能跳低，所以找与t[HIGH]的低点
    cut_signal = gor(sdl,stopl)

    cut_signal = select([t[VOLUME]>0],[cut_signal]) #默认为0，即未交易的日子卖出信号不能发出，否则会合并到下一交易日


    ssignal = gor(sol,cut_signal)
    stock.sellprice = select([cut_signal],[cut_price],default=t[OPEN])
    #止损和退回用cut_price, 阴线出局和停牌平移都用开盘价

    return ssignal

def follow_seller(stock,buy_signal,xstop=25,ret=50,**kwargs):
    '''
        如果价格小于最近5日高点5%，则卖出
        xstop为根据买入价的止损
        ret为从高点向下的回退值
    '''
    t = stock.transaction

    #从顶下落处理,前5天的收盘/开盘的高者和今天的开盘的高者 回落ret之后

    #hhret = gmax(rollx(tmax(gmax(t[OPEN],t[CLOSE]),5),1),t[OPEN])* (1000-ret)/1000
    hhret = gmax(rollx(tmax(t[HIGH],5),1),t[OPEN])* (1000-ret)/1000
    
    #hhret = rollx(tmax(t[HIGH],5),1) * (1000-ret)/1000
    sdl = t[LOW] < hhret
 
    #止损处理2.5%
    stop_price = extend2next(rollx(stock.buyprice,1) * (1000-xstop)/1000)
    stopl = t[LOW] < stop_price

    cut_price = gmin(gmax(hhret,stop_price),t[HIGH])    #首先，止损线和退回线高者先被触及，同时，穿越时可能跳低，所以找与t[HIGH]的低点
    cut_signal = gor(sdl,stopl)

    cut_signal = select([t[VOLUME]>0],[cut_signal]) #默认为0，即未交易的日子卖出信号不能发出，否则会合并到下一交易日

    bs = gand(buy_signal,cut_signal)
    rbs = rollx(bs)

    sell_signal = select([bs],[0],default=cut_signal) + rbs #如果当日冲销，则后推一日，但如果前一日也是当日，则信号被屏蔽

    stock.sellprice = select([cut_signal],[cut_price],default=t[OPEN])
    #止损和退回用cut_price, 当日卖出信号平移用开盘价，停牌平移用开盘价

    return cut_signal

 def pseudo_colr_amount(self):
     """
     Calculate pseudo Cost-of-Living Refund amount.
     Note this is simply meant to illustrate a Python programming technique;
     this function does NOT calculate an exact Cost-of-Living Refund amount.
     See setting of parameters above in specify_pseudo_COLR_policy method.
     """
     recs = self.__records
     # create MARS-specific policy parameter arrays
     mars_indicators = [recs.MARS == 1, recs.MARS == 2, recs.MARS == 3,
                        recs.MARS == 4, recs.MARS == 5]
     colr_c = np.select(mars_indicators, self.colr_param['COLR_c'])
     colr_ps = np.select(mars_indicators, self.colr_param['COLR_ps'])
     colr_rt = self.colr_param['COLR_rt']
     colr_prt = self.colr_param['COLR_prt']
     # compute colr_amt
     amt_pre_phaseout = np.minimum(recs.e00200 * colr_rt, colr_c)
     phaseout = np.maximum((recs.c00100 - colr_ps) * colr_prt, 0.)
     colr_amt = np.maximum(amt_pre_phaseout - phaseout, 0.)
     setattr(recs, 'colr_amount', colr_amt)
     # reduce income and combined taxes because COLR is a refundable credit
     recs.iitax -= colr_amt
     recs.combined -= colr_amt
     # delete local arrays used only in this method
     del mars_indicators
     del colr_c
     del colr_ps
     del amt_pre_phaseout
     del phaseout
     del colr_amt

        def loyer_retenu():

            # loyer mensuel réel, multiplié par 2/3 pour les meublés
            L1 = round_(loyer * where(statut_occupation == 5, 2 / 3, 1))

            zone_apl = simulation.calculate('zone_apl_famille', period)

            # Paramètres contenant les plafonds de loyer pour cette zone
            plafonds_by_zone = [[0] + [al.loyers_plafond[ 'zone' + str(zone) ][ 'L' + str(i) ] for zone in range(1, 4)] for i in range(1, 5)]
            L2_personne_seule = take(plafonds_by_zone[0], zone_apl)
            L2_couple = take(plafonds_by_zone[1], zone_apl)
            L2_famille = take(plafonds_by_zone[2], zone_apl) + (al_pac > 1) * (al_pac - 1) * take(plafonds_by_zone[3], zone_apl)

            L2 = select(
                [personne_seule * (al_pac == 0) + chambre, al_pac > 0],
                [L2_personne_seule, L2_famille],
                default = L2_couple
                )

            # taux à appliquer sur le loyer plafond
            coeff_chambre_colloc = select(
                [chambre, coloc],
                [al.loyers_plafond.chambre, al.loyers_plafond.colocation],
                default = 1)

            L2 = round_(L2 * coeff_chambre_colloc, 2)

            # loyer retenu
            L = min_(L1, L2)

            return L

def divideArraysSafely(num, den) :
    """Per evement divides numpy arrays result = num/den. Protected for 0 values. Arrays should have the same size."""
    if num.shape != den.shape :
        print 'divideArraysSafely: non-equal array shapes for numerator and denumerator: ', num.shape, den.shape
    num_corr =  np.select([den<1], [0], default=num)    
    den_corr =  np.select([den<1], [1], default=den)    
    return num_corr/den_corr

def split_x(x, split_pos):
    # NOTE: do not support multiple sentence tensors
    # sequence input , non-sequence input, and no non-sequence input
    # sequence input:
    if type(x) is not list:
        x=[x]

    if len(x) == 1:
        # sec1,                 sec2, sec3,...
        # sent1, sent2, sent5
        x01, x02 = tuple(np.split(x[0],[split_pos]))
        cond_list=[x02>=0,x02<0]
        offset = x02[0][0]
        choice_list=[x02-offset, x02 ]
        x02 = np.select(cond_list, choice_list)
        return ([x01],[x02])

    # doc1 doc2 doc3
    # sec1 sec2 ...

    # sec1, sec2, ...
    # sent1, sent2, ...

    x01, x02 = tuple(np.split(x[0], [split_pos]))
    offset = x02[0][0]
    x1, x2 = split_x(x[1:], offset)
    cond_list = [x02 >= 0, x02 < 0]
    choice_list = [x02 - offset, x02]
    x02 = np.select(cond_list, choice_list)
    return ([x01] + x1, [x02]+x2)

    def update_particles(self, delta):
        # radial: posx + posy
        norm = numpy.sqrt(self.particle_pos[:, 0] ** 2 + self.particle_pos[:, 1] ** 2)
        # XXX prevent div by 0
        norm = numpy.select([norm == 0], [0.0000001], default=norm)
        posx = self.particle_pos[:, 0] / norm
        posy = self.particle_pos[:, 1] / norm

        radial = numpy.array([posx, posy])
        tangential = numpy.array([-posy, posx])

        # update dir
        radial = numpy.swapaxes(radial, 0, 1)
        radial *= self.particle_rad
        tangential = numpy.swapaxes(tangential, 0, 1)
        tangential *= self.particle_tan

        self.particle_dir += (tangential + radial + self.particle_grav) * delta

        # update pos with updated dir
        self.particle_pos += self.particle_dir * delta

        # life
        self.particle_life -= delta

        # color
        self.particle_color += self.particle_delta_color * delta

        # if life < 0, set alpha in 0
        self.particle_color[:, 3] = numpy.select([self.particle_life[:, 0] < 0], [0], default=self.particle_color[:, 3])

    def _sanitize_SchedEnergyDiff(self):
        if not self.hasEvents('sched_energy_diff') \
           or 'nrg_model' not in self.platform:
            return
        nrg_model = self.platform['nrg_model']
        em_lcluster = nrg_model['little']['cluster']
        em_bcluster = nrg_model['big']['cluster']
        em_lcpu = nrg_model['little']['cpu']
        em_bcpu = nrg_model['big']['cpu']
        lcpus = len(self.platform['clusters']['little'])
        bcpus = len(self.platform['clusters']['big'])
        SCHED_LOAD_SCALE = 1024

        power_max = em_lcpu['nrg_max'] * lcpus + em_bcpu['nrg_max'] * bcpus + \
                    em_lcluster['nrg_max'] + em_bcluster['nrg_max']
        print "Maximum estimated system energy: {0:d}".format(power_max)

        df = self.df('sched_energy_diff')
        df['nrg_diff_pct'] = SCHED_LOAD_SCALE * df.nrg_diff / power_max

        # Tag columns by usage_delta
        ccol = df.usage_delta
        df['usage_delta_group'] = np.select(
            [ccol < 150, ccol < 400, ccol < 600],
            ['< 150', '< 400', '< 600'], '>= 600')

        # Tag columns by nrg_payoff
        ccol = df.nrg_payoff
        df['nrg_payoff_group'] = np.select(
            [ccol > 2e9, ccol > 0, ccol > -2e9],
            ['Optimal Accept', 'SchedTune Accept', 'SchedTune Reject'], 'Suboptimal Reject')

def supdown2(sopen,sclose,shigh,slow):
    ''' 计算每日的上升行程和下降行程
        以距离开盘价距离近的方向为运行方向
        则若最低近,运行轨迹为 开盘-->最低-->最高-->收盘
          若最高近,运行轨迹为 开盘-->最高-->最低-->收盘
          平开往低走
          另，如果开盘大于昨日收盘，则上升段 +　开盘－昨收盘
            　小于昨日收盘，则下降段　+ 昨收盘 - 开盘
    '''
    if len(sopen) == 0:
        return np.array([],int),np.array([],int)
    sc1 = rollx(sclose)
    sc1[0] = sopen[0]   #前一日收盘价视同首日开盘价
    u_hlc = shigh-sopen+sclose-slow
    u_lhc = shigh - slow
    d_hlc = shigh - slow
    d_lhc = sopen-slow+shigh-sclose
    ou = np.where(sopen > sc1)
    od = np.where(sopen < sc1)
    doc = sopen-sc1
    u_hlc[ou] = u_hlc[ou] + doc[ou]
    u_lhc[ou] = u_lhc[ou] + doc[ou]
    d_hlc[od] = d_hlc[od] - doc[od] #doc[od]<0
    d_lhc[od] = d_lhc[od] - doc[od] #doc[od]<0
    is_up = shigh-sopen < sopen-slow #True为向上，False为向下
    u = np.select([is_up],[u_hlc],default=u_lhc)
    d = np.select([is_up],[d_hlc],default=d_lhc)
    return u,d

 def get_intens_for_stat_q_bins(sp, intens_map) :
     q_map_stat = sp.get_q_map_for_stat_bins()
     counts = sp.get_counts_for_stat_q_bins()
     # print 'counts = ', counts
     intens = bincount(q_map_stat, intens_map, sp.ana_stat_part_q+1)
     counts_prot = np.select([counts<=0.], [-1.], default=counts)
     intens_aver = np.select([counts_prot<=0.], [0.], default=intens/counts_prot)
     return intens_aver

 def get_g2_for_dyna_bins_itau(sp, itau) :
     q_phi_map_dyna = sp.get_q_phi_map_for_dyna_bins()
     g2_map         = sp.get_g2_map_for_itau(itau)
     intens_dyna    = sp.bincount(q_phi_map_dyna, g2_map, sp.npart_dyna)
     counts         = sp.get_counts_for_dyna_bins()
     counts_prot    = np.select([counts==0], [-1], default=counts) 
     sp.g2_for_dyna_bins = np.select([counts_prot<0], [0], default=intens_dyna/counts_prot)
     return sp.g2_for_dyna_bins

    def get_q_average_for_dyna_bins(sp) :
        if sp.q_average_dyna != None : return sp.q_average_dyna

        q_map_masked      = sp.get_q_map() * sp.get_mask_total()
        sum_q_dyna        = sp.bincount(sp.get_q_phi_map_for_dyna_bins(), q_map_masked, length=sp.npart_dyna)
        counts_dyna       = sp.get_counts_for_dyna_bins()
        counts_dyna_prot  = np.select([counts_dyna<=0], [-1], counts_dyna)
        sp.q_average_dyna = np.select([counts_dyna_prot<0], [0], default=sum_q_dyna/counts_dyna_prot)
        print 'get_q_average_for_dyna_bins():\n', sp.q_average_dyna
        return sp.q_average_dyna