def test_reduceND(self):
        from numpy import add, arange

        a = arange(12).reshape(3, 4)
        assert (add.reduce(a, 0) == [12, 15, 18, 21]).all()
        assert (add.reduce(a, 1) == [6.0, 22.0, 38.0]).all()
        raises(ValueError, add.reduce, a, 2)

 def test_reduce_keepdims(self):
     from numpy import add, arange
     a = arange(12).reshape(3, 4)
     b = add.reduce(a, 0, keepdims=True)
     assert b.shape == (1, 4)
     assert (add.reduce(a, 0, keepdims=True) == [12, 15, 18, 21]).all()
     assert (add.reduce(a, 0, None, None, True) == [12, 15, 18, 21]).all()

def make_F_matrix(E_matrix):
    """takes an E matrix and returns an F matrix

    input is output of make_E_matrix

    for each element in matrix subtract mean of corresponding row and 
    column and add the mean of all elements in the matrix
    """
    num_rows, num_cols = shape(E_matrix)
    #make a vector of the means for each row and column
    #column_means = (add.reduce(E_matrix) / num_rows)
    column_means = (add.reduce(E_matrix) / num_rows)[:,newaxis]
    trans_matrix = transpose(E_matrix)
    row_sums = add.reduce(trans_matrix)
    row_means = row_sums / num_cols
    #calculate the mean of the whole matrix
    matrix_mean = sum(row_sums) / (num_rows * num_cols)
    #adjust each element in the E matrix to make the F matrix

    E_matrix -= row_means
    E_matrix -= column_means
    E_matrix += matrix_mean

    #for i, row in enumerate(E_matrix):
    #    for j, val in enumerate(row):
    #        E_matrix[i,j] = E_matrix[i,j] - row_means[i] - \
    #                column_means[j] + matrix_mean
    return E_matrix

def _basic_simps(y,start,stop,x,dx,axis):
    nd = len(y.shape)
    if start is None:
        start = 0
    step = 2
    all = (slice(None),)*nd
    slice0 = tupleset(all, axis, slice(start, stop, step))
    slice1 = tupleset(all, axis, slice(start+1, stop+1, step))
    slice2 = tupleset(all, axis, slice(start+2, stop+2, step))

    if x is None:  # Even spaced Simpson's rule.
        result = add.reduce(dx/3.0* (y[slice0]+4*y[slice1]+y[slice2]),
                                    axis)
    else:
        # Account for possibly different spacings.
        #    Simpson's rule changes a bit.
        h = diff(x,axis=axis)
        sl0 = tupleset(all, axis, slice(start, stop, step))
        sl1 = tupleset(all, axis, slice(start+1, stop+1, step))
        h0 = h[sl0]
        h1 = h[sl1]
        hsum = h0 + h1
        hprod = h0 * h1
        h0divh1 = h0 / h1
        result = add.reduce(hsum/6.0*(y[slice0]*(2-1.0/h0divh1) + \
                                              y[slice1]*hsum*hsum/hprod + \
                                              y[slice2]*(2-h0divh1)),axis)
    return result

def gaussian_convolution(data, ijk_linewidths):

  from numpy import float32, zeros, add, divide, outer, reshape
  if data.dtype.type != float32:
    data = data.astype(float32)

  from math import exp
  gaussians = []
  for a in range(3):
    size = data.shape[a]
    gaussian = zeros((size,), float32)
    hw = ijk_linewidths[2-a] / 2.0
    for i in range(size):
      u = min(i,size-i) / hw
      p = min(u*u/2, 100)               # avoid OverflowError with exp()
      gaussian[i] = exp(-p)
    area = add.reduce(gaussian)
    divide(gaussian, area, gaussian)
    gaussians.append(gaussian)

  g01 = outer(gaussians[0], gaussians[1])
  g012 = outer(g01, gaussians[2])
  g012 = reshape(g012, data.shape)
  
  cdata = zeros(data.shape, float32)

  from numpy.fft import fftn, ifftn
  # TODO: Fourier transform Gaussian analytically to reduce computation time
  #       about 30% (one of three fft calculations).
  ftg = fftn(g012)
  ftd = fftn(data)
  gd = ifftn(ftg * ftd)
  gd = gd.astype(float32)
  return gd

    def plane(self, matrix):

        from numpy import ravel, minimum, maximum, add, multiply, array, float32
        matrix_1d = matrix.ravel()
        dmin = minimum.reduce(matrix_1d)
        if self.min == None or dmin < self.min:
            self.min = dmin
        dmax = maximum.reduce(matrix_1d)
        if self.max == None or dmax > self.max:
            self.max = dmax
        self.sum += add.reduce(matrix_1d)
        # TODO: Don't copy array to get standard deviation.
        # Avoid overflow when squaring integral types
        m2 = array(matrix_1d, float32)
        multiply(m2, m2, m2)
        self.sum2 += add.reduce(m2)

    def test_reduce_errors(self):
        from numpy import sin, add, maximum, zeros

        raises(ValueError, sin.reduce, [1, 2, 3])
        assert add.reduce(1) == 1

        assert list(maximum.reduce(zeros((2, 0)), axis=0)) == []
        exc = raises(ValueError, maximum.reduce, zeros((2, 0)), axis=None)
        assert exc.value[0] == ("zero-size array to reduction operation " "maximum which has no identity")
        exc = raises(ValueError, maximum.reduce, zeros((2, 0)), axis=1)
        assert exc.value[0] == ("zero-size array to reduction operation " "maximum which has no identity")

        a = zeros((2, 2)) + 1
        assert (add.reduce(a, axis=1) == [2, 2]).all()
        assert (add.reduce(a, axis=(1,)) == [2, 2]).all()
        exc = raises(ValueError, add.reduce, a, axis=2)
        assert exc.value[0] == "'axis' entry is out of bounds"

def sample_from_histogram(p, n_samples=1):
    """
    returns the indice of bin according to the histogram p

    @param p: histogram
    @type p: numpy.array
    @param n_samples: number of samples to generate
    @type n_samples: integer
    """
    
    from numpy import add, less, argsort, take, arange
    from numpy.random import random

    indices = argsort(p)
    indices = take(indices, arange(len(p) - 1, -1, -1))

    c = add.accumulate(take(p, indices)) / add.reduce(p)

    return indices[add.reduce(less.outer(c, random(n_samples)), 0)]

 def getStats(self, histo):
     m = array(histo[1:])
     sum = 0.0
     sum2 = 0.0
     n = float(add.reduce(m))
     for j in range(len(m)):
         sum = sum + j * m[j]
         sum2 = sum2 + (j ** 2) * float(m[j])
     var = (sum2-(sum**2.0)/n)/n
     return sum/n,sqrt(var)

def gaussian(sdev, size):

  from math import exp
  from numpy import empty, single as floatc, add, divide

  g = empty((size,), floatc)
  for i in range(size):
    u = min(i,size-i) / sdev
    p = min(u*u/2, 100)               # avoid OverflowError with exp()
    g[i] = exp(-p)
  area = add.reduce(g)
  divide(g, area, g)
  return g

def _make_f_matrix(matrix):
    """It takes an E matrix and returns an F matrix

    The input is the output of make_E_matrix

    For each element in matrix subtract mean of corresponding row and
    column and add the mean of all elements in the matrix
    """
    num_rows, num_cols = matrix.shape
    # make a vector of the means for each row and column
    # column_means = (add.reduce(E_matrix) / num_rows)
    column_means = (add.reduce(matrix) / num_rows)[:, newaxis]
    trans_matrix = transpose(matrix)
    row_sums = add.reduce(trans_matrix)
    row_means = row_sums / num_cols
    # calculate the mean of the whole matrix
    matrix_mean = nsum(row_sums) / (num_rows * num_cols)
    # adjust each element in the E matrix to make the F matrix

    matrix -= row_means
    matrix -= column_means
    matrix += matrix_mean

    return matrix

def RombergMethod(y, dx, show=False):

    axis=-1
    y = asarray(y)
    nd = len(y.shape)
    Nsamps = y.shape[axis]
    Ninterv = Nsamps-1
    n = 1
    k = 0

    while n < Ninterv:
        n <<= 1
        k += 1

    R = {}
    all = (slice(None),) * nd
    slice0 = tupleset(all, axis, 0)
    slicem1 = tupleset(all, axis, -1)
    h = Ninterv*asarray(dx)*1.0
    R[(1,1)] = (y[slice0] + y[slicem1])/2.0*h
    slice_R = all
    start = stop = step = Ninterv
    for i in range(2,k+1):
        start >>= 1
        slice_R = tupleset(slice_R, axis, slice(start,stop,step))
        step >>= 1
        R[(i,1)] = 0.5*(R[(i-1,1)] + h*add.reduce(y[slice_R],axis))
        for j in range(2,i+1):
            R[(i,j)] = R[(i,j-1)] + \
                       (R[(i,j-1)]-R[(i-1,j-1)]) / ((1 << (2*(j-1)))-1)
        h = h / 2.0

    if show:
        precis = 5
        width = 8
        formstr = "%" + str(width) + '.' + str(precis)+'f'

        print('\nMétodo de Romberg')
        print('----------------------------------')
        for i in range(1,k+1):
            for j in range(1,i+1):
                print(formstr % R[(i,j)], end=' ')
            print()
        print('----------------------------------')

    return R[(k,k)]

def notes_roc (la, lb, eps):
    from numpy import transpose, add, resize 
    """ creates a matrix of size len(la)*len(lb) then look for hit and miss
    in it within eps tolerance windows """
    gdn,fpw,fpg,fpa,fdo,fdp = 0,0,0,0,0,0
    m = len(la)
    n = len(lb)
    x =           resize(la[:][0],(n,m))
    y = transpose(resize(lb[:][0],(m,n)))
    teps =  (abs(x-y) <= eps[0]) 
    x =           resize(la[:][1],(n,m))
    y = transpose(resize(lb[:][1],(m,n)))
    tpitc = (abs(x-y) <= eps[1]) 
    res = teps * tpitc
    res = add.reduce(res,axis=0)
    for i in range(len(res)) :
        if res[i] > 1:
            gdn+=1
            fdo+=res[i]-1
        elif res [i] == 1:
            gdn+=1
    fpa = n - gdn - fpa
    return gdn,fpw,fpg,fpa,fdo,fdp

def distanceA2AEuclideanSquared(x, std=[], w=[]):
    """
    This function calcule the Euclidean Squared distance between
    two or more variables.
    """
    if std:
        x = nparray(x)
        x = stdobs(x)  #  standardize
        x = x.tolist()
    if w:
        x = nparray(x)
        w = w / float(npadd.reduce(w))
        x = x * w  #  weights
        x = x.tolist()

    numrows = len(x)
    distance = [0]*(numrows-1)

    for row in xrange(numrows - 1):
        npsublist = npsubtract(x[row], x[row + 1])
        sublist = npsublist.tolist()
        distance[row] = [square_double(sublist)]

    return distance

    def test_reduce_1d(self):
        import numpy as np
        from numpy import array, add, maximum, less, float16, complex64

        assert less.reduce([5, 4, 3, 2, 1])
        assert add.reduce([1, 2, 3]) == 6
        assert maximum.reduce([1]) == 1
        assert maximum.reduce([1, 2, 3]) == 3
        raises(ValueError, maximum.reduce, [])

        assert add.reduce(array([True, False] * 200)) == 200
        assert add.reduce(array([True, False] * 200, dtype="int8")) == 200
        assert add.reduce(array([True, False] * 200), dtype="int8") == -56
        assert type(add.reduce(array([True, False] * 200, dtype="float16"))) is float16
        assert type(add.reduce(array([True, False] * 200, dtype="complex64"))) is complex64

        for dtype in ["bool", "int"]:
            assert np.equal.reduce([1, 2], dtype=dtype) == True
            assert np.equal.reduce([1, 2, 0], dtype=dtype) == False

def romb(y, dx=1.0, axis=-1, show=False):
    """
    Romberg integration using samples of a function.

    Parameters
    -----------
    y : array_like
        A vector of ``2**k + 1`` equally-spaced samples of a function.
    dx : array_like, optional
        The sample spacing. Default is 1.
    axis : array_like?, optional
        The axis along which to integrate. Default is -1 (last axis).
    show : bool, optional
           When y is a single 1-D array, then if this argument is True
           print the table showing Richardson extrapolation from the
           samples. Default is False.

    Returns
    -------
    ret : array_like?
        The integrated result for each axis.

    See also
    --------
    quad - adaptive quadrature using QUADPACK
    romberg - adaptive Romberg quadrature
    quadrature - adaptive Gaussian quadrature
    fixed_quad - fixed-order Gaussian quadrature
    dblquad, tplquad - double and triple integrals
    simps, trapz - integrators for sampled data
    cumtrapz - cumulative integration for sampled data
    ode, odeint - ODE integrators

    """
    y = asarray(y)
    nd = len(y.shape)
    Nsamps = y.shape[axis]
    Ninterv = Nsamps-1
    n = 1
    k = 0
    while n < Ninterv:
        n <<= 1
        k += 1
    if n != Ninterv:
        raise ValueError("Number of samples must be one plus a "
                "non-negative power of 2.")

    R = {}
    all = (slice(None),) * nd
    slice0 = tupleset(all, axis, 0)
    slicem1 = tupleset(all, axis, -1)
    h = Ninterv*asarray(dx)*1.0
    R[(1,1)] = (y[slice0] + y[slicem1])/2.0*h
    slice_R = all
    start = stop = step = Ninterv
    for i in range(2,k+1):
        start >>= 1
        slice_R = tupleset(slice_R, axis, slice(start,stop,step))
        step >>= 1
        R[(i,1)] = 0.5*(R[(i-1,1)] + h*add.reduce(y[slice_R],axis))
        for j in range(2,i+1):
            R[(i,j)] = R[(i,j-1)] + \
                       (R[(i,j-1)]-R[(i-1,j-1)]) / ((1 << (2*(j-1)))-1)
        h = h / 2.0

    if show:
        if not isscalar(R[(1,1)]):
            print("*** Printing table only supported for integrals" + \
                  " of a single data set.")
        else:
            try:
                precis = show[0]
            except (TypeError, IndexError):
                precis = 5
            try:
                width = show[1]
            except (TypeError, IndexError):
                width = 8
            formstr = "%" + str(width) + '.' + str(precis)+'f'

            print("\n       Richardson Extrapolation Table for Romberg Integration       ")
            print("====================================================================")
            for i in range(1,k+1):
                for j in range(1,i+1):
                    print(formstr % R[(i,j)], end=' ')
                print()
            print("====================================================================\n")

    return R[(k,k)]

def model(SWnet, LWnet, Hs, Hl):
    N = SWnet.shape[0]
    Rnet_array = EnergyNetFluxBalance(SWnet, LWnet, Hs, Hl)
    Rnet = add.reduce(Rnet_array, axis=0) / N
    return Rnet

def Sum_j(f):
    return add.reduce(f,1)

def Sum_i(f):
    return add.reduce(f,0)

    def Intersectar(self,tipo='real'):
        '''!
        @brief: Método que cálcula la intersección entre las dos líneas.
        @param tipo str: opción que permite establecer si las intersecciones calculadas son reales o virtuales.
        @note tipo: Valores validos: real,virtual 
        '''
        self.__checkTipo(tipo)
        #Coordenadas de cada punto.
        xi1=self.__l1.getPuntoInicial().getX()
        yi1=self.__l1.getPuntoInicial().getY()
        
        xf1=self.__l1.getPuntoFinal().getX()
        yf1=self.__l1.getPuntoFinal().getY()
        
        xi2=self.__l2.getPuntoInicial().getX()
        yi2=self.__l2.getPuntoInicial().getY()
        
        xf2=self.__l2.getPuntoFinal().getX()
        yf2=self.__l2.getPuntoFinal().getY()
        
        v1=[xf1-xi1,yf1-yi1]
        v2=[xf2-xi2,yf2-yi2]
        A=matrix([[v1[1],-v1[0]],
                  [v2[1],-v2[0]]])
        B=array([xf1*v1[1]-yf1*v1[0],xf2*v2[1]-yf2*v2[0]])
        #print(A,B)
        if det(A)==0:
            A=add.reduce(A, 0)
            B=add.reduce(B, 0)
            #print(A,B)
            if B==0:
                # sol ay=bx
                return None
            if B!=0:
                if A.item((0, 0))==0:
                    x=xi1
                    y=B/A.item((0, 1))
                elif A.item((0, 1))==0:
                    x=B/A.item((0, 0))
                    y=yi1
        else:
            sol=solve(A,B)
            x=sol[0]
            y=sol[1]
            #print(x,y)
        #Cálculo de la interseccion #y=ax+b
#         a1=None
#         a2=None
#         b1=None
#         b2=None
#             
#         if abs(yf1-yi1)==0.0:
#             a1=0.0
#         elif abs(xf1-xi1)==0.0:
#             a1=None
#         else:
#             a1=(yf1-yi1)/(xf1-xi1)
#             
#         if abs(yf2-yi2)==0.0:
#             a2=0.0
#         elif abs(xf2-xi2)==0.0:
#             a2=None
#         else:
#             a2=(yf2-yi2)/(xf2-xi2)
#             
#         if a1==None:
#             b1=xi1
#         else:
#             b1=yi1-a1*xi1
#             
#         if a2==None:
#             b2=xi2
#         else:
#             b2=yi2-a2*xi2
#         
#             
# #         print(a1,a2,b1,b2)
#         if a1==None and a2==None:
#             #Paralelas.
#             return None
#         elif a1==None and a2==0.0:
#             x=xi2
#             y=a2
#             
#         elif a1==0.0 and a2==None:
#             x=b2
#             y=b1
#         elif a1==None:
#             x=b1
#             y=(b2-b1)/(-a2)
#         elif a2==None:
#             x=b2
#             y=(b2-b1)/(a1)
#         elif a1==None and a2==None:
#             x=(b2-b1)
#         else:
#             try:
#                 x=(b2-b1)/(a1-a2)
#                 y=a1*x+b1
#             except:
#.........这里部分代码省略.........