def integral(A=None,dF=None,F=None,axis = <NUM_LIT:0>,trapez = False,cumulative = False):
ndim = max(v.ndim for v in (A,dF,F) if v is not None)<EOL>def broadcast(x):<EOL><INDENT>new_shape = [<NUM_LIT:1>]*ndim<EOL>new_shape[axis] = -<NUM_LIT:1><EOL>return np.reshape(x,new_shape)<EOL><DEDENT>if F is not None:<EOL><INDENT>assert(dF is None)<EOL>if F.ndim<ndim:<EOL><INDENT>F = broadcast(F)<EOL><DEDENT>N = F.shape[axis]<EOL>dF = F.take(indices = range(<NUM_LIT:1>,N),axis = axis)-F.take(indices = range(N-<NUM_LIT:1>),axis = axis)<EOL><DEDENT>elif dF is not None:<EOL><INDENT>if dF.ndim<ndim:<EOL><INDENT>dF = broadcast(dF)<EOL><DEDENT>N = dF.shape[axis]+<NUM_LIT:1><EOL><DEDENT>else:<EOL><INDENT>if A.ndim<ndim:<EOL><INDENT>A = broadcast(A)<EOL><DEDENT>N = A.shape[axis]<EOL><DEDENT>if A is not None:<EOL><INDENT>if trapez:<EOL><INDENT>midA = (A.take(indices = range(<NUM_LIT:1>,N),axis = axis)+A.take(indices = range(N-<NUM_LIT:1>),axis = axis))/<NUM_LIT:2><EOL><DEDENT>else:<EOL><INDENT>midA = A.take(indices=range(N-<NUM_LIT:1>),axis=axis)<EOL><DEDENT>if dF is not None:<EOL><INDENT>dY = midA*dF<EOL><DEDENT>else:<EOL><INDENT>dY = midA<EOL><DEDENT><DEDENT>else:<EOL><INDENT>dY = dF<EOL><DEDENT>pad_shape = list(dY.shape)<EOL>pad_shape[axis] = <NUM_LIT:1><EOL>pad = np.zeros(pad_shape)<EOL>if cumulative:<EOL><INDENT>return np.concatenate((pad,np.cumsum(dY,axis = axis)),axis = axis)<EOL><DEDENT>else:<EOL><INDENT>return np.sum(dY,axis = axis)<EOL><DEDENT>
Turns an array A of length N (the function values in N points)
and an array dF of length N-1 (the masses of the N-1 intervals)
into an array of length N (the integral \int A dF at N points, with first entry 0)

:param A: Integrand (optional, default ones, length N)
:param dF: Integrator (optional, default ones, length N-1)
:param F: Alternative to dF (optional, length N)
:param trapez: Use trapezoidal rule (else left point)