def test_frompyfunc_2d_sig(self):
        import sys
        from numpy import frompyfunc, dtype, arange

        if "__pypy__" not in sys.builtin_module_names:
            skip("PyPy only frompyfunc extension")

        def times_2(in_array, out_array):
            assert len(in_array.shape) == 2
            assert in_array.shape == out_array.shape
            out_array[:] = in_array * 2

        ufunc = frompyfunc(
            [times_2], 1, 1, signature="(m,n)->(n,m)", dtypes=[dtype(int), dtype(int)], stack_inputs=True
        )
        ai = arange(18, dtype=int).reshape(2, 3, 3)
        ai3 = ufunc(ai[0, :, :])
        ai2 = ufunc(ai)
        assert (ai2 == ai * 2).all()

        ufunc = frompyfunc(
            [times_2], 1, 1, signature="(m,m)->(m,m)", dtypes=[dtype(int), dtype(int)], stack_inputs=True
        )
        ai = arange(12 * 3 * 3, dtype="int32").reshape(12, 3, 3)
        exc = raises(ValueError, ufunc, ai[:, :, 0])
        assert "perand 0 has a mismatch in its core dimension 1" in exc.value.message
        ai3 = ufunc(ai[0, :, :])
        ai2 = ufunc(ai)
        assert (ai2 == ai * 2).all()
        # view
        aiV = ai[::-2, :, :]
        assert aiV.strides == (-72, 12, 4)
        ai2 = ufunc(aiV)
        assert (ai2 == aiV * 2).all()

def read_file ( filename ):
    """
    Lit un fichier USPS et renvoie un tableau de tableaux d'images.
    Chaque image est un tableau de nombres réels.
    Chaque tableau d'images contient des images de la même classe.
    Ainsi, T = read_file ( "fichier" ) est tel que T[0] est le tableau
    des images de la classe 0, T[1] contient celui des images de la classe 1,
    et ainsi de suite.
    """
    # lecture de l'en-tête
    infile = open ( filename, "r" )    
    nb_classes, nb_features = [ int( x ) for x in infile.readline().split() ]

    # creation de la structure de données pour sauver les images :
    # c'est un tableau de listes (1 par classe)
    data = np.empty ( 10, dtype=object )   
    filler = np.frompyfunc(lambda x: list(), 1, 1)
    filler( data, data )

    # lecture des images du fichier et tri, classe par classe
    for ligne in infile:
        champs = ligne.split ()
        if len ( champs ) == nb_features + 1:
            classe = int ( champs.pop ( 0 ) )
            data[classe].append ( map ( lambda x: float(x), champs ) )     
    infile.close ()

    # transformation des list en array
    output  = np.empty ( 10, dtype=object )
    filler2 = np.frompyfunc(lambda x: np.asarray (x), 1, 1)
    filler2 ( data, output )

    return output

    def test_frompyfunc_sig_broadcast(self):
        import sys
        from numpy import frompyfunc, dtype, arange

        if "__pypy__" not in sys.builtin_module_names:
            skip("PyPy only frompyfunc extension")

        def sum_along_0(in_array, out_array):
            out_array[...] = in_array.sum(axis=0)

        def add_two(in0, in1, out):
            out[...] = in0 + in1

        ufunc_add = frompyfunc(
            add_two,
            2,
            1,
            signature="(m,n),(m,n)->(m,n)",
            dtypes=[dtype(int), dtype(int), dtype(int)],
            stack_inputs=True,
        )
        ufunc_sum = frompyfunc(
            [sum_along_0], 1, 1, signature="(m,n)->(n)", dtypes=[dtype(int), dtype(int)], stack_inputs=True
        )
        ai = arange(18, dtype=int).reshape(3, 2, 3)
        aout = ufunc_add(ai, ai[0, :, :])
        assert aout.shape == (3, 2, 3)
        aout = ufunc_sum(ai)
        assert aout.shape == (3, 3)

    def plot_time_function(self, p):
        """Plot the time function.
        """
        n_steps = self.n_steps
        mats = self.mats
        step_size = self.step_size

        ls_t = linspace(0, step_size * n_steps, n_steps + 1)
        ls_fn = frompyfunc(self.time_function, 1, 1)
        ls_v = ls_fn(ls_t)

        p.subplot(321)
        p.plot(ls_t, ls_v, "ro-")

        final_epsilon = self.final_displ / self.length

        kappa = linspace(mats.epsilon_0, final_epsilon, 10)
        omega_fn = frompyfunc(lambda kappa: mats._get_omega(None, kappa), 1, 1)
        omega = omega_fn(kappa)
        kappa_scaled = step_size + (1 - step_size) * (kappa - mats.epsilon_0) / (final_epsilon - mats.epsilon_0)
        xdata = hstack([array([0.0], dtype=float), kappa_scaled])
        ydata = hstack([array([0.0], dtype=float), omega])
        p.plot(xdata, ydata, "g")
        p.xlabel("regular time [-]")
        p.ylabel("scaled time [-]")

def mov_average_expw(data, span, tol=1e-6):
    """Calculates the exponentially weighted moving average of a series.

:Parameters:
    $$data$$
    span : int 
        Time periods. The smoothing factor is 2/(span + 1)
    tol : float, *[1e-6]*
        Tolerance for the definition of the mask. When data contains masked 
        values, this parameter determinea what points in the result should be masked.
        Values in the result that would not be "significantly" impacted (as 
        determined by this parameter) by the masked values are left unmasked."""

    data = marray(data, copy=True, subok=True)
    ismasked = (data._mask is not nomask)
    data._mask = N.zeros(data.shape, bool_)
    _data = data._data
    #
    k = 2./float(span + 1)
    def expmave_sub(a, b):
        return a + k * (b - a)
    #
    data._data.flat = N.frompyfunc(expmave_sub, 2, 1).accumulate(_data)
    if ismasked:
        _unmasked = N.logical_not(data._mask).astype(float_)
        marker = 1. - N.frompyfunc(expmave_sub, 2, 1).accumulate(_unmasked)
        data._mask[marker > tol] = True
    data._mask[0] = True
    #
    return data

def __get_function__(astr, width):
    f = np.cos
    if astr == "cosine":
        f = np.cos
    elif astr == "rampup":
        f = np.frompyfunc(
            lambda x: w.ramp_up(x, width),
            1, 1
            )
    elif astr == "rampdown":
        f = np.frompyfunc(
            lambda x: w.ramp_down(x, width),
            1, 1
            )
    elif astr == "impulse":
        f = np.frompyfunc(
            lambda x: w.impulse(x, width),
            1, 1
            )
    elif astr == "step":
        f = np.frompyfunc(
            lambda x: w.step(x, width),
            1, 1
            )
    else:
        f = np.cos
    return f

    def __init__(self, rng, nin, nout, activation=logistic, activation_prime=dlogistic, W=None, b=None, inputs=None, learningRate=.9):
        self.ActivationFn = np.frompyfunc(activation, 1, 1)
        self.DActivationFn = np.frompyfunc(activation_prime, 1, 1)
        self.inputs = inputs
        self.activations = None
        self.activationHistory = None
        self.outputs = None
        self.outputHistory = None
        self.learningRate = learningRate
        self.momentumFactor = .7
        self.previousDelta = None
        self.previousbDelta = None

        if not W:
            self.W = np.asarray(
                rng.uniform(
                    low=4 * np.sqrt(6.0 / (nin + nout)),     # generic range of values
                    high=-4 * np.sqrt(6.0 / (nin + nout)),
                    size=(nin, nout)
                )
                , dtype=float
            )
        else:
            self.W = W
        if not b:
            self.b = np.zeros(nout)
        else:
            self.b = b

    def polar_workspace_init(radial_bins=256, angular_bins=256,
                             max_radius=None, centre=None): 
        #if (centre == None) and self.centre == None:
            #pass # Raise an exception

        xdim = self.image.shape[0]
        ydim = self.image.shape[1]

        if centre == None:
            xc = xdim * 0.5
            yc = ydim * 0.5 
        else:
            xc = centre[0]
            yc = centre[1]

        # Calculate minimum distance from centre to edge of image - this
        # determines the maximum radius in the polar image
        xsize = min (xdim + 0.5 - xc, xc)
        ysize = min (ydim + 0.5 - yc, yc)
        max_rad = m.sqrt(xsize**2 + ysize**2)

        if max_radius == None:
            max_radius = max_rad
        elif max_radius > max_rad:
            raise ValueError
        
        # Set up interpolation - cubic spline with no smoothing by default 
        x = numpy.indices((xdim,)) + 0.5 - centre[0]
        y = numpy.indices((ydim,)) + 0.5 - centre[1]
        interp = spint.RectBivariateSpline(x, y, self.image)

        # Polar image bin widths
        theta_bin_width = (2.0 * math.pi) / (theta_bins - 1.0)
        radial_bin_width = max_radius / (radial_bins - 1.0)

        # Calculate polar image values - use vectorization for efficiency
        # Because we broadcast when using a ufunc (created by frompyfunc
        # below), we could get away with an ogrid here to save time and space?
        r, theta = numpy.mgrid[0:radial_bins, 0:angular_bins]
        theta = (theta + 0.5) * theta_bin_width
        r = (r + 0.5) * radial_bin_width

        def polar_pix_val(r, theta):
            # Should we use the numpy.sin/cos functions here for more
            # efficiency ?
            return interp.ev(r * m.sin(theta), r * m.cos(theta))

        numpy.frompyfunc(polar_pix_val, 2, 1)
        self.pimage = polar_pix_val(r, theta)

        # Calculate polar image values - non-vectorized version
        self.pimage = numpy.empty(radial_bins, angular_bins)
        for r in radial_bins:
            R = (r + 0.5) * radial_bin_width;
            for t in theta_bins:
                theta = (t + 0.5) * theta_bin_width
                x = R * sin(theta)
                y = R * cos(theta)
                self.pimage[r, t] = interp.ev(x, y)

def Model2a():
    """
    A version of Model3 where the activation functions are generated
    using a different numpy function. For exploratory reasons only.
    """
    m2a=Model2()
    sn = m2a.GetGroupByName("SN")
    sp = m2a.GetGroupByName("SP")
    sn.SetActivationFunction(np.frompyfunc(lambda x: neural.STanh_plus(x, gain=2)), 1,1)
    sp.SetActivationFunction(np.frompyfunc(lambda x: neural.STanh_plus(x, gain=2)), 1,1) 

    def redraw(self, e = None):
        if ((self.idx_x < 0 and len(self.idx_x_arr) == 0) or
             (self.idx_y < 0 and len(self.idx_y_arr) == 0) or
             self._xdata == [] or
             self._ydata == []):
            return
        #
        if len(self.idx_x_arr) > 0:
            print 'x: summation for', self.idx_x_arr
            xarray = np.array(self._xdata)[:, self.idx_x_arr].sum(1)
        else:
            xarray = np.array(self._xdata)[:, self.idx_x]

        if len(self.idx_y_arr) > 0:
            print 'y: summation for', self.idx_y_arr
            yarray = np.array(self._ydata)[:, self.idx_y_arr].sum(1)

#            print 'yarray', yarray
#            yarray_arr = array( self._ydata )[:, self.idx_y_arr]
#            sym_weigth_arr = 2. * ones_like( yarray_arr[1] )
#            sym_weigth_arr[0] = 4.
#            print 'yarray_arr', yarray_arr
#            print 'sym_weigth_arr', sym_weigth_arr
#            yarray = dot( yarray_arr, sym_weigth_arr )
#            print 'yarray', yarray


        else:
            yarray = np.array(self._ydata)[:, self.idx_y]

        if self.transform_x:
            def transform_x_fn(x):
                '''makes a callable function out of the Str-attribute
                "transform_x". The vectorised version of this function is 
                then used to transform the values in "xarray". Note that
                the function defined in "transform_x" must be defined in
                terms of a lower case variable "x".
                '''
                return eval(self.transform_x)
            xarray = np.frompyfunc(transform_x_fn, 1, 1)(xarray)

        if self.transform_y:
            def transform_y_fn(y):
                '''makes a callable function out of the Str-attribute
                "transform_y". The vectorised version of this function is 
                then used to transform the values in "yarray". Note that
                the function defined in "transform_y" must be defined in
                terms of a lower case variable "y".
                '''
                return eval(self.transform_y)
            yarray = np.frompyfunc(transform_y_fn, 1, 1)(yarray)

        self.trace.xdata = np.array(xarray)
        self.trace.ydata = np.array(yarray)
        self.trace.data_changed = True

 def _refresh_fired(self):
     xdata = linspace(0,10,10)
     fneval1 = frompyfunc( lambda x: eval( self.expression1 ), 1, 1 )
     fneval2 = frompyfunc( lambda x: eval( self.expression2 ), 1, 1 )
     fneval3 = frompyfunc( lambda x: eval( self.expression3 ), 1, 1 )
     y1 = fneval1( xdata )
     y2 = fneval2( xdata )
     y3 = fneval3( xdata )
     ydata = column_stack((y1,y2,y3))
     self.mfn.set( xdata = xdata, ydata = ydata )
     self.mfn.data_changed = True

 def __init__(self, mu, alpha):
     self.map_pdf = {}
     self.map_logpdf = {}
     self.bins = []
     mu = float(mu)
     
     self.alpha = alpha
     self.mu = mu
     
     self.nbin = np.frompyfunc(self._get_value, 3, 1)
     self.nbin_log = np.frompyfunc(self._get_value_log, 3, 1)

    def _get_values(self):
        l_rho = self.l_rho
        n_points = self.n_points
        gl_b = self._get_gbundle_props()[0]
        gmu_b = self._get_gbundle_props()[1]
        m_f = self.m_f
        mu_r = self.mu_r
        l_r = self.l_r

        # for Gaussian bundle strength distribution 
        if self.l_plot <= gl_b:
            gl_arr = logspace( log( self.min_plot_length,10 ), log(gl_b,10), n_points )
            gstrength_arr = self.fl(gl_arr)/self.fl(self.l_r)*self.mu_r
        elif self.l_plot > gl_b:
            gl_1 = logspace( log( self.min_plot_length,10), log(gl_b,10), n_points )
            gl_2 = logspace( log( gl_b, 10 ), log( self.l_plot, 10 ), n_points )
            gl_arr = hstack( (gl_1, gl_2) )
            gstrength_1 = self.fl( gl_1 ) / self.fl( self.l_r ) * self.mu_r
            gstrength_22 = frompyfunc( self._get_gstrength, 1, 1 )
            gstrength_2 = array( gstrength_22( gl_2 ), dtype = 'float64' )
            gstrength_arr = hstack( ( gstrength_1, gstrength_2 ) )
        # Mirek's mean approximation
            strength_22 = frompyfunc( self.mean_approx, 1, 3 )
            strength_2 = array( strength_22( gl_2 )[0], dtype = 'float64' )
            mean_gumb = array( strength_22( gl_2 )[1], dtype = 'float64' )
            med_gumb = array( strength_22( gl_2 )[2], dtype = 'float64' )

        #asymptotes for the first two branches
        if self.l_plot <= l_rho:
            al_arr = array([self.min_plot_length, self.l_plot])
            astrength_arr = array([mu_r / self.fl(l_r), mu_r / self.fl(l_r)])
        elif l_rho < self.l_plot:
            al_arr = array([self.min_plot_length, l_rho, 10. * gl_b])
            astrength_1 = mu_r / self.fl(l_r)
            astrength_2 = (l_rho/al_arr[2])**(1/m_f) * astrength_1
            astrength_arr = hstack((astrength_1,astrength_1,astrength_2))
        
        # left asymptote
        self.mu_sigma_0 = astrength_arr[0]
            
        # standard deviation for the first branch = before fragmentation
        if self.l_plot <= gl_b:
            stl_arr = logspace( log( self.min_plot_length,10) , log(self.l_plot,10) , 
                                n_points / 2. )
            stdev_arr_plus  = self.fl( stl_arr ) / self.fl( l_r ) * mu_r * (1 + self.cov)
            stdev_arr_minus = self.fl( stl_arr ) / self.fl( l_r ) * mu_r * (1 - self.cov)
        else:
            stl_arr = logspace( log( self.min_plot_length,10), log(gl_b,10), n_points)
            stdev_arr_plus = self.fl( stl_arr )/self.fl( l_r ) * mu_r * (1 + self.cov)
            stdev_arr_minus = self.fl( stl_arr )/self.fl( l_r ) * mu_r * (1 - self.cov)        
        
        return gl_arr, al_arr, gstrength_arr, astrength_arr,\
               stl_arr, stdev_arr_plus, stdev_arr_minus, gl_2,\
               strength_2, mean_gumb, med_gumb

 def add_dates(df):
     fun_em = np.frompyfunc(get_date_em, 2, 1)
     fun_lig = np.frompyfunc(get_date_lig, 2, 1)
     df['date_emergence_leaf'] = fun_em(df['num_leaf_bottom'], df['fnl'])
     df['date_ligulation_leaf'] = fun_lig(df['num_leaf_bottom'], df['fnl'])
     if force_mean_fnl==False:
         df['date_emergence_flag_leaf'] = map(lambda fnl: hs_fit.TTemleaf(fnl, nff=fnl), df['fnl'])
         df['date_ligulation_flag_leaf'] = map(lambda fnl: hs_fit.TTligleaf(fnl, nff=fnl), df['fnl'])
     else:
         df['date_emergence_flag_leaf'] = hs_fit.TTemleaf(hs_fit.mean_nff, nff=None)[0]
         df['date_ligulation_flag_leaf'] = hs_fit.TTligleaf(hs_fit.mean_nff, nff=None)[0]
     return df

 def _get_mfn_plot(self):
     n_points = 100
     sigma_max = self.sigma_fu * self.rho
     
     sigma_arr = linspace( 0, sigma_max, n_points )
     
     get_epsilon_f = frompyfunc( lambda sigma: sigma / self.E_f, 1, 1 ) 
     epsilon_f_arr  = get_epsilon_f( sigma_arr )
     
     get_epsilon_c = frompyfunc( self._get_epsilon_c, 1, 1 )
     epsilon_c_arr = get_epsilon_c( sigma_arr )
     
     return MFnLineArray( xdata = epsilon_c_arr, ydata = sigma_arr )

 def _refresh_fired(self):
     
     #creates an empty plot data container as a list of MFnLineArray classes
     self.mfn.lines = self.No_of_all_curves * [MFnLineArray()]
     
     xdata = linspace(0,10,100)
     fneval1 = frompyfunc( lambda x: eval( self.expression1 ), 1, 1 )
     fneval2 = frompyfunc( lambda x: eval( self.expression2 ), 1, 1 )
     fneval3 = frompyfunc( lambda x: eval( self.expression3 ), 1, 1 )
     self.mfn.lines[0] = MFnLineArray(xdata = xdata, ydata = fneval1( xdata ))
     self.mfn.lines[1] = MFnLineArray(xdata = xdata, ydata = fneval2( xdata ))
     self.mfn.lines[2] = MFnLineArray(xdata = xdata, ydata = fneval3( xdata ))
     self.mfn.data_changed = True

    def std_of_flow_time_to_stream(streamlink, flow_dir_file, slope, radius, velocity, delta_s_file,
                                   flow_dir_code='TauDEM'):
        """Generate standard deviation of t0_s (flow time to the workflow channel from each cell).
        """
        strlk_r = RasterUtilClass.read_raster(streamlink)
        strlk_data = strlk_r.data
        rad_data = RasterUtilClass.read_raster(radius).data
        slo_data = RasterUtilClass.read_raster(slope).data

        vel_r = RasterUtilClass.read_raster(velocity)
        vel_data = vel_r.data
        xsize = vel_r.nCols
        ysize = vel_r.nRows
        nodata_value = vel_r.noDataValue

        def initial_variables(vel, strlk, slp, rad):
            """initial variables"""
            if abs(vel - nodata_value) < UTIL_ZERO:
                return DEFAULT_NODATA
            if strlk <= 0:
                tmp_weight = 1
            else:
                tmp_weight = 0
            # 0 is river
            if slp < 0.0005:
                slp = 0.0005
            # dampGrid = vel * rad / (slp / 100. * 2.) # No need to divide 100
            # in my view. By LJ
            damp_grid = vel * rad / (slp * 2.)
            celerity = vel * 5. / 3.
            tmp_weight *= damp_grid * 2. / numpy.power(celerity, 3.)
            return tmp_weight

        initial_variables_numpy = numpy.frompyfunc(initial_variables, 4, 1)
        weight = initial_variables_numpy(vel_data, strlk_data, slo_data, rad_data)

        delta_s_sqr = TerrainUtilClass.calculate_flow_length(flow_dir_file, weight, flow_dir_code)

        def cal_delta_s(vel, sqr):
            """Calculate delta s"""
            if abs(vel - nodata_value) < UTIL_ZERO:
                return nodata_value
            else:
                return sqrt(sqr) / 3600.

        cal_delta_s_numpy = numpy.frompyfunc(cal_delta_s, 2, 1)
        delta_s = cal_delta_s_numpy(vel_data, delta_s_sqr)

        RasterUtilClass.write_gtiff_file(delta_s_file, ysize, xsize, delta_s, strlk_r.geotrans,
                                         strlk_r.srs, DEFAULT_NODATA, GDT_Float32)

 def __init__(self, func, nin, nout):
     self._ufunc = np.frompyfunc(func, nin, nout)
     self._func = func
     self.nin = nin
     self.nout = nout
     self._name = funcname(func)
     self.__name__ = 'frompyfunc-%s' % self._name

    def IC_spec_gamma(self,EdNdE_ph,Eph,EdNdE_e,Ee,en):
       
        dlogE_e = np.log(Ee[1] / Ee[0]) * np.ones_like(Ee)#assumes log spacing of energy bins
        dNe = EdNdE_e * dlogE_e
        #print 'TEST: dlogE_e in IC: ',dlogE_e
        #print 'TEST: Ee[3]/Ee[2]: ',np.log(Ee[3]/Ee[2])
        #print 'TEST: Ee[6]/Ee[5]: ',np.log(Ee[6]/Ee[5])

        dLogE_ph = np.log(Eph[1] / Eph[0]) * np.ones_like(Eph)
        dN_ph = EdNdE_ph * dLogE_ph #changed
        '''print 'dNe: ',dNe
        print 'TEST: Eph: ',Eph
        print 'TEST: Ee: ',Ee
        print 'TEST: en: ',en
        print 'TEST: dLogE_ph: ',dLogE_ph
        print 'TEST: Eph[1]/Eph[0]: ',np.log(Eph[1]/Eph[0])
        print 'TEST: Eph[3]/Eph[2]: ',np.log(Eph[3]/Eph[2])
        print 'TEST: Eph[10]/Eph[9]: ',np.log(Eph[10]/Eph[9])
        print 'TEST: dN_ph: ',dN_ph'''
        
        def EdNdE_gamma(x):
            sigma=crs.sigma_ic(x,Eph,Ee)
            return const.C*const.m2cm*np.dot(dN_ph,np.dot(sigma,dNe))
        
        csec_vec = np.frompyfunc(EdNdE_gamma, 1, 1)
        return csec_vec(en)

    def pi0_spectrum(self,dNdp_p, p_p,E_g,n_H=1.):
        #print 'dNdp: ',dNdp_p
        #print 'E_g : ',E_g
        kappa_pi = 0.17
        dp_p = p_p[1:] - p_p[:-1]
        dNdp_p = np.sqrt(dNdp_p[1:] * dNdp_p[:-1])

        step = lambda x: (1. + np.sign(x))/2.
        positive = lambda x: x * step(x)

        E_p = np.sqrt(p_p**2 + const.MP_GeV**2)
        E_p = np.sqrt(E_p[1:] * E_p[:-1])
        E_kin = E_p - const.MP_GeV
        kin_mask = step(E_kin**2 / kappa_pi**2 - const.MPi0_GeV**2)

        p_pi = np.sqrt(positive(E_kin**2 / kappa_pi**2 - const.MPi0_GeV**2)) + const.epsilon

        def EdNdE_gamma(E_g):
            E_pi_min = E_g + const.MPi0_GeV**2 / (4 * E_g)
            E_p_min =  const.MP_GeV+ E_pi_min / kappa_pi
            dNpi = const.C*const.m2cm* n_H * crs.sigma_pp(E_kin) * dNdp_p * dp_p
            mask = kin_mask * step(E_p - E_p_min)
            return E_g * 2 * np.sum(dNpi / p_pi * mask)
        EdNdE_gamma_vec = np.frompyfunc(EdNdE_gamma, 1, 1)
        return EdNdE_gamma_vec(E_g)