Python numpy.piecewise函数代码示例

        def _arc_cosine_hypothesis_test(ac_output, X, W, b, s):
            """
            Runs hypothesis test for Arc Cosine layer.

            Inputs:
                ac_output -- output of net after running arc cosine layer
                X -- input data
                W -- weight parameter from train_init_net
                b -- bias parameter from train_init_net
                s -- degree parameter
            """
            # Get output from net
            net_output = workspace.FetchBlob(ac_output)

            # Computing output directly
            x_rand = np.matmul(X, np.transpose(W)) + b
            x_pow = np.power(x_rand, s)
            if s > 0:
                h_rand_features = np.piecewise(x_rand,
                                               [x_rand <= 0, x_rand > 0],
                                               [0, 1])
            else:
                h_rand_features = np.piecewise(x_rand,
                                               [x_rand <= 0, x_rand > 0],
                                               [0, lambda x: x / (1 + x)])
            output_ref = np.multiply(x_pow, h_rand_features)

            # Comparing net output and computed output
            npt.assert_allclose(net_output, output_ref, rtol=1e-3, atol=1e-3)

        def EpEvo(t,eta,j,Gamma0,q):

            
            #
            # For now I'm going to set a few params at good values
            # Adjust later

            z=2.332

            #q=power(10.,q)
            #q=1.E-3
            #E0 = 1.E54
            E0 = 1.97E55        
            #E0 = power(10.,E0)
            
            g = (j+1.-eta)*0.5

            n0 = 1.E2

            #xd = ((3.-eta)*E0 / ( 4.*pi*n0*Gamma0**2. * mp  ) )**(1./3.)
            xd = 1.8E16*((1.+j-eta)*(E0/1.E54)/((n0/100.)*(Gamma0/300.)))**(1./3.)
            #td = (1.+z)*xd / (Gamma0**2. * c)
            td = 6.7*(1.+z)*((1.+j-eta)*(E0/1.E54)/((n0/100.)*(Gamma0/300.)**8.))**(1./3.)



            ### Calculate X(t)  ###

            test = td/(2. * g + 1.) *( Gamma0**(2.+1./g) + 2.*g)
            
            condition1 = t<td
            condition2 = logical_and(td<=t, t<=test)


            X = piecewise(t, [condition1, condition2],\
            [lambda t: t/td, \
            lambda t: ((2.*g+1.)*(t/td) - 2.*g)**(1./(2.*g+1.)) ])

            ### Calculate X(t)  ###



            ### Calculate Gamma(X)  ###

            condition3 = X<1.
            condition4 = logical_and(1.<=X, X<=Gamma0**(1./g))

            Gamma = piecewise(X, [condition3, condition4],\
            [lambda X: Gamma0, \
            lambda X: Gamma0*X**(-g) ])

            ### Calculate Gamma(X)  ###


            eE0 = 3.E-8 * (1E-3)  * n0**(.5)*q*Gamma0**4. /(1.+z)


            result = 511.8 *  eE0*(Gamma/Gamma0)**4. * (X)**(-eta/2.)
            
            return result

def makewdnoise(f,wdstyle,obs='X'):
    if wdstyle == 'mldc':
        x = 2.0 * math.pi * model.lisaL * f
        t = 4 * x**2 * N.sin(x)**2 * (1.0 if obs == 'X' else 1.5)
        
        return t * ( N.piecewise(f,(f >= 1.0e-4  ) & (f < 1.0e-3  ),[lambda f: 10**-44.62 * f**-2.3, 0]) + \
                     N.piecewise(f,(f >= 1.0e-3  ) & (f < 10**-2.7),[lambda f: 10**-50.92 * f**-4.4, 0]) + \
                     N.piecewise(f,(f >= 10**-2.7) & (f < 10**-2.4),[lambda f: 10**-62.8  * f**-8.8, 0]) + \
                     N.piecewise(f,(f >= 10**-2.4) & (f < 10**-2.0),[lambda f: 10**-89.68 * f**-20.0,0])     )
    elif wdstyle in wdnoise:
        mod, p = wdnoise[wdstyle]
        p = p[0] if obs == 'X' else p[1] # assume AE if not X
        y = N.log10(f)
    
        if mod == 'rat42':
            return 10.0**( (p[0]*y**4+p[1]*y**3+p[2]*y**2+p[3]*y+p[4])/(y**2+p[5]*y+p[6]) )
        elif mod == 'poly4':
            return 10.0**( p[0]*y**4+p[1]*y**3+p[2]*y**2+p[3]*y+p[4] )
        else:
            raise NotImplementedError
    else:
        if '.txt' in wdstyle:
            conf = N.loadtxt(wdstyle)
            conf[N.isnan(conf[:,1]),1] = 0
        
            return N.interp(f,conf[:,0],conf[:,1])
        else:
            raise NotImplementedError

        def EpEvo(t,eta,g,Gamma0,q):


            #
            # For now I'm going to set a few params at good values
            # Adjust later

            z=2.332
            
            #            q=1.E-4
#            q=1E-3


            
            E0 = 9.1E53        
            
            n0 = 1.E2

            #n0 = power(10,n0)
            
            #xd = ((3.-eta)*E0 / ( 4.*pi*n0*Gamma0**2. * mp  ) )**(1./3.)
            xd = 2.6E16*((1.-eta/3.)*(E0/1.E54)/((n0/100.)*(Gamma0/300.)))**(1./3.)
            #td = (1.+z)*xd / (Gamma0**2. * c)
            td = 9.7*(1.+z)*((1.-eta/3.)*(E0/1.E54)/((n0/100.)*(Gamma0/300.)**8.))**(1./3.)

            ### Calculate X(t)  ###
            test = td/(2. * g + 1.) *( Gamma0**(2.+1./g) + 2.*g)
            #frac = t/td

            condition1 = t<td
            condition2 = logical_and(td<=t, t<=test)


            X = piecewise(t, [condition1, condition2],\
            [lambda t: t/td, \
            lambda t: ((2.*g+1.)*(t/td) - 2.*g)**(1./(2.*g+1.)) ])

            ### Calculate X(t)  ###



            ### Calculate Gamma(X)  ###



            condition3 = X<1.
            condition4 = logical_and(1.<=X, X<=Gamma0**(1./g))

            Gamma = piecewise(X, [condition3, condition4],\
            [lambda X: Gamma0, \
            lambda X: Gamma0*X**(-g) ])

            ### Calculate Gamma(X)  ###
    
            eE0 = 3.E-8 * (1E-3)  * n0**(.5)*q  *Gamma0**4. /(1.+z)

#            eE0 = 3.E-8 * n0**(.5)*q*Gamma0**4. /(1.+z)

            return 511. *  eE0*(Gamma/Gamma0)**4. * (X/xd)**(-eta/2.)

 def __init__(self, Ymax, Y0, mumax, lag, x):
         # x is x data      
          
     a1 = np.piecewise(x, [x < lag, x >= lag], [0.0, 1.0])
     tc = (Ymax-Y0)/mumax+lag
     a2 = np.piecewise(x, [x < tc, x >= tc], [0.0, 1.0])
     self.Buchanan = Y0*(1.0-a1)+(Y0 + mumax*(x-lag))*a1*(1.0-a2) +  Ymax*a2
     
     print x, self.Buchanan

    def _get_solution(self):
        """Get the solution from the solver output

        Fills the dictionary self.sol with the information:
            * 's': The optimal s as a function of time
            * 't': The time vector
            * 'states': Numerical values of the states defined in self.sys

        TODO: perform accurate integration to determine time
        TODO: Do exact interpolation
        """
        solver = self.prob['solver']
        N = self.options['N']
        x_opt = np.array(solver.getOutput("x")).ravel()
        b_opt = np.reshape(x_opt, (N + 1, -1), order='F')
        self.sol['b'] = b_opt

        # Determine time on a sufficiently fine spatial grid
        s0 = np.linspace(0, 1, 1001)
        delta = s0[1] - s0[0]

        pieces = [lambda s, b=b_opt, ss=ss, j=j:
                  sum([bb * (s - ss) ** i / fact[i]
                       for i, bb in enumerate(b[j])])
                  for j, ss in enumerate(self.prob['s'][:-1])]
        conds = lambda s0: [np.logical_and(self.prob['s'][i] <= s0,
                                           s0 <= self.prob['s'][i+1])
                            for i in range(N)]
        b0_opt = np.piecewise(s0, conds(s0), pieces)
        b0_opt[b0_opt < 0] = 0
        time = np.cumsum(np.hstack([0, 2 * delta / (np.sqrt(b0_opt[:-1]) +
                                    np.sqrt(b0_opt[1:]))]))

        # Resample to constant time-grid
        t = np.arange(time[0], time[-1], self.options['Ts'])
        st = np.interp(t, time, s0)
        # Evaluate solution on equidistant time grid
        b_opt = np.c_[[np.piecewise(st, conds(st), pieces, b=b_opt[:, i:])
                       for i in range(self.sys.order)]].T
        st = np.matrix(st)

        # Determine s and derivatives from b_opt
        b, Ds = self._make_path()[1:]
        Ds_f = cas.SXFunction([b], [Ds])  # derivatives of s wrt b
        Ds_f.init()
        s_opt = np.hstack((st.T, np.array([evalf(Ds_f, bb).toArray().ravel()
                                          for bb in b_opt])))
        self.sol['s'] = np.asarray(s_opt)
        self.sol['t'] = t
        # Evaluate the states
        f = cas.SXFunction([self.s], [cas.substitute(cas.vertcat(self.sys.x.values()),
                                           self.sys.y, self.path)])
        f_val = np.array([evalf(f, s.T).toArray().ravel() for s in s_opt])
        self.sol['states'] = dict([(k, f_val[:, i]) for i, k in
                          enumerate(self.sys.x.keys())])

def I(w, wa, wb, wc):
    wa,wb,wc = float(wa), float(wb), float(wc)
    f = (wb - wa) / (wc - wa)
    eps = (wc-wa)/1e3
    if (wb-wa)>eps and (wc-wb)>eps:
        w1 = wa+eps
        w2 = wb-eps
        w3 = wc-eps
        return np.piecewise(w, 
                            [(w <= wa-eps),
                            (wa-eps < w) & (w <= wa+eps),
                            (wa+eps < w) & (w < wb-eps),
                            (wb-eps <= w) & (w <= wb+eps),
                            (wb+eps < w) & (w < wc-eps),
                            (wc-eps <= w) & (w < wc+eps),
                            (wc+eps <= w)],
                            [0,
                             4/np.sqrt(f*(1-((w1 - wa) / (wc - wa)))) * K([(1-f)*((w1 - wa) / (wc - wa)) / (f*(1-((w1 - wa) / (wc - wa))))]),
                             lambda w: 4/np.sqrt(f*(1-((w - wa) / (wc - wa)))) * K((1-f)*((w - wa) / (wc - wa)) / (f*(1-((w - wa) / (wc - wa))))),
                             4/np.sqrt(f*(1-((w2 - wa) / (wc - wa)))) * K([(1-f)*((w2 - wa) / (wc - wa)) / (f*(1-((w2 - wa) / (wc - wa))))]),
                             lambda w: 4/np.sqrt(((w - wa) / (wc - wa))*(1-f)) * K((1-((w - wa) / (wc - wa)))*f / (((w - wa) / (wc - wa))*(1-f))),
                             4/np.sqrt(((w3 - wa) / (wc - wa))*(1-f)) * K([(1-((w3 - wa) / (wc - wa)))*f / (((w3 - wa) / (wc - wa))*(1-f))]),
                             0])
    elif (wb-wa)<=eps:
        w2 = wb+eps
        w3 = wc-eps
        return np.piecewise(w, 
                            [(w <= wa-eps),
                            (wa-eps < w) & (w <= wb+eps),
                            (wb+eps < w) & (w < wc-eps),
                            (wc-eps <= w) & (w < wc+eps),
                            (wc+eps <= w)],
                            [0,
                             4/np.sqrt(((w2 - wa) / (wc - wa))*(1-f)) * K([(1-((w2 - wa) / (wc - wa)))*f / (((w2 - wa) / (wc - wa))*(1-f))]),
                             lambda w: 4/np.sqrt(((w - wa) / (wc - wa))*(1-f)) * K((1-((w - wa) / (wc - wa)))*f / (((w - wa) / (wc - wa))*(1-f))),
                             4/np.sqrt(((w3 - wa) / (wc - wa))*(1-f)) * K([(1-((w3 - wa) / (wc - wa)))*f / (((w3 - wa) / (wc - wa))*(1-f))]),
                             0])
    elif (wc-wb)<=eps:
        w1 = wa+eps
        w2 = wb-eps
        return np.piecewise(w, 
                            [(w <= wa-eps),
                            (wa-eps < w) & (w <= wa+eps),
                            (wa+eps < w) & (w < wb-eps),
                            (wb-eps <= w) & (w <= wc+eps),
                            (wc+eps <= w)],
                            [0,
                             4/np.sqrt(f*(1-((w1 - wa) / (wc - wa)))) * K([(1-f)*((w1 - wa) / (wc - wa)) / (f*(1-((w1 - wa) / (wc - wa))))]),
                             lambda w: 4/np.sqrt(f*(1-((w - wa) / (wc - wa)))) * K((1-f)*((w - wa) / (wc - wa)) / (f*(1-((w - wa) / (wc - wa))))),
                             4/np.sqrt(f*(1-((w2 - wa) / (wc - wa)))) * K([(1-f)*((w2 - wa) / (wc - wa)) / (f*(1-((w2 - wa) / (wc - wa))))]),
                             0])

def EpEvo(t,A,eta,g,E0,Gamma0,n0,q):

        #c = 2.99E10 #cm/s
    c=1.
    mp = 1.67E-26   # keV ??

    z=1.
    #q=1.E-3
    #Gamma0 = 300.

    #g = (3.-eta)/2.
    #n0 = 1.E2

    #xd = ((3.-eta)*E0 / ( 4.*pi*n0*Gamma0**2. * mp  ) )**(1./3.)
    xd = 2.6E16*((1.-eta/3.)*(E0/1.E54)/((n0/100.)*(Gamma0/300.)))**(1./3.)
    #td = (1.+z)*xd / (Gamma0**2. * c)
    td = 9.7*(1.+z)*((1.-eta/3.)*(E0/1.E54)/((n0/100.)*(Gamma0/300.)**8.))**(1./3.)
        
    ### Calculate X(t)  ###
    test = (td/(2. * g + 1.) * Gamma0**(2.+1./g) + 2.*g)
    #frac = t/td

    condition1 = t<td
    condition2 = logical_and(td<=t, t<=test)

        
    X = piecewise(t, [condition1, condition2],\
    [lambda t: t/td, \
    lambda t: ((2.*g+1.)*(t/td) - 2.*g)**(1./(2.*g+1.)) ])

    ### Calculate X(t)  ###



    ### Calculate Gamma(X)  ###



    condition3 = X<1.
    condition4 = logical_and(1.<=X, X<=Gamma0**(1./g))

    Gamma = piecewise(X, [condition3, condition4],\
    [lambda X: Gamma0, \
    lambda X: Gamma0*X**(-g) ])

    ### Calculate Gamma(X)  ###

        
    eE0 = 3.E-8 * n0**(.5)*q*Gamma0**4. /(1.+z)

    return A*eE0*(Gamma/Gamma0)**4. * (X/xd)**(-eta/2.)

 def _evaluate(self, t):
     if self.direction == 1:
         kernel = np.piecewise(
             t, [t < 0, t >= 0], [
                 lambda t: 0,
                 lambda t: (1.0 / self._sigma_scaled.magnitude) * np.exp(
                     (-t / self._sigma_scaled).magnitude)]) / t.units
     elif self.direction == -1:
         kernel = np.piecewise(
             t, [t < 0, t >= 0], [
                 lambda t: (1.0 / self._sigma_scaled.magnitude) * np.exp(
                     (t / self._sigma_scaled).magnitude),
                 lambda t: 0]) / t.units
     return kernel

 def loadExampleData():
     data = pd.DataFrame({'Seconds': np.arange(0,15,1/200)})
     t = data.loc[:,'Seconds'].values
     data.loc[:,'T ACP'] = \
     np.piecewise(t, [t < 4, 4 <= t], [1,0]) * (1800/4*t) + \
     np.piecewise(t, [t < 4,(4<=t)*(t<=12)], [0,1]) * (1800 + (1200-1800)/8*(t-4)) + \
     np.piecewise(t, [12 < t, 20 < t], [1,0]) * (1500 + -400/8*(t-12))
      
     data.loc[:,'time'] = t           
     self.data = data
     self.initExampleData()
      
     self.camera = 'rad' 
     self.times = [4,7,12]          

def phi(xs, Z):
    """Evaluate phi(x) on a numpy list.
    
    This uses invibro.phi.phi0_cache to interpolate values."""
    Z0, x_bound = phi0_cache['Z0'], phi0_cache['xs'][0]
    c1 = pi ** 2 /6; c2 = 7 * pi ** 4 / 60.0
    neg = lambda x: log((Z - x) / (Z + x))
    interp = lambda x: phi0_cache['interp'](x) + log((Z + x) / (Z0 + x))
    large = lambda x: log(1 + Z / x) + c1 / x ** 2 + c2 / x ** 4
    return (
        piecewise(xs, [xs < 0], [neg, 0.0]) +
        piecewise(abs(xs), [abs(xs) < x_bound], [interp, large]) +
        complex(0, -0.5) * n(xs)
    )

def reposition(coords, box_size):
    com = compute_com_pbc(coords, box_size)
    cob = box_size / 2.0
    coords_recenter = coords - com + cob
    coords_recenter_x = coords_recenter[:,0]
    coords_recenter_y = coords_recenter[:,1]
    coords_recenter_z = coords_recenter[:,2]
    #print coords_recenter
    coords_recenter_x = np.piecewise(coords_recenter_x, [coords_recenter_x < 0.0, (coords_recenter_x >= 0.0) * (coords_recenter_x <= box_size[0]), coords_recenter_x > box_size[0]], \
        [lambda coords_recenter_x: coords_recenter_x + box_size[0], lambda coords_recenter_x: coords_recenter_x, lambda coords_recenter_x: coords_recenter_x - box_size[0]])
    coords_recenter_y = np.piecewise(coords_recenter_y, [coords_recenter_y < 0.0, (coords_recenter_y >= 0.0) * (coords_recenter_y <= box_size[1]), coords_recenter_y > box_size[1]], \
        [lambda coords_recenter_y: coords_recenter_y + box_size[1], lambda coords_recenter_y: coords_recenter_y, lambda coords_recenter_y: coords_recenter_y - box_size[1]])
    coords_recenter_z = np.piecewise(coords_recenter_z, [coords_recenter_z < 0.0, (coords_recenter_z >= 0.0) * (coords_recenter_z <= box_size[2]), coords_recenter_z > box_size[2]], \
        [lambda coords_recenter_z: coords_recenter_z + box_size[2], lambda coords_recenter_z: coords_recenter_z, lambda coords_recenter_z: coords_recenter_z - box_size[2]])
    return np.array(zip(coords_recenter_x,coords_recenter_y,coords_recenter_z))

def I_axial(w, wa, wc):
    if wc > wa:
        return np.piecewise(w, [(w <= wa),
                                (wa < w) & (w < wc),
                                (wc <= w)],
                                [0, 
                                 lambda w: 1 / (2*np.sqrt((w - wa) / (wc - wa))),
                                 0])
    elif wa > wc:
        return np.piecewise(w, [(w <= wc),
                            (wc < w) & (w < wa),
                            (wa <= w)],
                            [0, 
                             lambda w: 1 / (2*np.sqrt((w - wa) / (wc - wa))),
                             0])

def f(x, a1, b1, c1, d1, b2, c2, v0):
    a2_calculated = a2(a1, b1, c1, d1, b2, c2, v0)
    d2_calculated = d2(a1, b1, c1, d1, b2, c2, v0)
    return np.piecewise(
            x,
            [x < v0, x >= v0],
            [lambda x: p(x, a1, b1, c1, d1), lambda x: p(x, a2_calculated, b2, c2, d2_calculated)])

def intgl_simp38(f,a,b,steps=-1,h=1):
    if steps>0:
        xis = np.linspace(a,b,steps+1)
        h  = xis[1]-xis[0]
    fxis = f(xis)
    wis = np.zeros(steps+1)
    pcs = []; fpcs = []
    for i in xrange(0,steps-2,3):
        wis[i:i+4] += [1,3,3,1]
        pcs.append(xis[i:i+4])
        fpcs.append(fxis[i:i+4])
    wis *= 3*h/8
    if steps%3==2:
        wis[-3:] += [h/3,4*h/3,h/3]
        pcs.append(xis[-3:])
        fpcs.append(fxis[-3:])
    elif steps%3==1:
        wis[-2:] += [h/2,h/2]
        pcs.append(xis[-2:])
        fpcs.append(fxis[-2:])
    fapprox = lambda x: np.piecewise(x,
                                     [np.logical_and(p[0]<=x,x<=p[-1]) for p in pcs],
                                     [lagrange(pcs[i],fpcs[i]) for i in xrange(len(pcs))])# np.interp(x,xis,fxis)
    # fapprox = lambda x: np.interp(x,xis,fxis)
    return (sum(fxis*wis),xis,fxis,wis,fapprox) # h/2 * sum(np.array([f(x) for x in xs]) * np.array([1]+[2]*(len(xs)-2)+[1]))