numpy.gradient的逆 *** 作

numpy.gradient的逆 *** 作 numpy.gradient的逆 *** 作

numpy.gradient可以对函数进行求导，在设置edge_order=1时其逆 *** 作可以由以下代码实现。改写自stackoverflow的回答。该程序仅适用于一维导数且dx为常数的情况，但是可以指定axis：

import numpy as np

def integrate(dydx, y0, dx=1, axis=0):
    ''' 
    Inverse of numpy.gradient, only for 1d, evenly spaced x and edge_order=1. 
    dydx: array_like, (N1, N2, ..., N_axis, ...), output of gradient
    y0: array_like, (N1, N2, ...), the initial value, one dimension less than dydx
    dx: float, the dx used in np.gradient
    axis: int, the axis used in np.gradient
    return: ndarray, the same as the data before np.gradient
    '''
    dydx = np.asarray(dydx)
    dydx = np.moveaxis(dydx, axis, 0)
    y1 = np.append(dydx[:1]*0., dydx[1:-1:2].cumsum(axis=0), axis=0)
    y2 = dydx[::2].cumsum(axis=0) - dydx[:1] / 2 
    y1 = np.expand_dims(y1, 1)
    y2 = np.expand_dims(y2, 1)
    y3 = np.hstack([y1, y2])
    y3 = y3.reshape(-1, *y1.shape[2:])
    out = y0 + 2 * y3[:dydx.shape[0]] * dx
    out = np.moveaxis(out, 0, axis)
    return out 

if __name__ == '__main__':
    y = np.random.rand(2, 10, 3)
    dydx = np.gradient(y, 0.143, edge_order=1, axis=1)
    y1 = integrate(dydx, y[:,0,:], 0.143, axis=1)
    print(y)
    print(y1)
    print('=============')
    print(np.abs(y1-y).max())