Python numpy.gradient() Examples

def test_log_interp_diff(self):
    """ Test the differentiation facility from the class log_interp_c """
    #import matplotlib.pyplot as plt
    rr,pp = funct_log_mesh(1024, 0.001, 20.0)
    logi = log_interp_c(rr)
    gc = 1.2030
    ff = np.array([np.exp(-gc*r**2) for r in rr])
    ffd_ref = np.array([np.exp(-gc*r**2)*(-2.0*gc*r) for r in rr])
    ffd = logi.diff(ff)
    ffd_np = np.gradient(ff, rr)
    s = 3
    for r,d,dref,dnp in zip(rr[s:],ffd[s:],ffd_ref[s:],ffd_np[s:]):
      self.assertAlmostEqual(d,dref)
      
    #plt.plot(rr, ff, '-', label='ff')
    #plt.plot(rr, ffd, '--', label='ffd')
    #plt.legend()
    #plt.show() 

def _ecg_delineator_peak_P_onset(rpeak, heartbeat, R, P):
    if P is None:
        return np.nan

    segment = heartbeat.iloc[:P]  # Select left of P wave
    try:
        signal = signal_smooth(segment["Signal"].values, size=R / 10)
    except TypeError:
        signal = segment["Signal"]

    if len(signal) < 2:
        return np.nan

    signal = np.gradient(np.gradient(signal))
    P_onset = np.argmax(signal)

    from_R = R - P_onset  # Relative to R
    return rpeak - from_R 

def _ecg_delineator_peak_T_offset(rpeak, heartbeat, R, T):
    if T is None:
        return np.nan

    segment = heartbeat.iloc[R + T :]  # Select left of P wave
    try:
        signal = signal_smooth(segment["Signal"].values, size=R / 10)
    except TypeError:
        signal = segment["Signal"]

    if len(signal) < 2:
        return np.nan

    signal = np.gradient(np.gradient(signal))
    T_offset = np.argmax(signal)

    return rpeak + T + T_offset


# =============================================================================
# Internals
# ============================================================================= 

def check_segmentation(fn_subject):
    from scipy import ndimage
    import matplotlib.pylab as pl
    from matplotlib.colors import ListedColormap
    files = simnibs.SubjectFiles(fn_subject + '.msh')
    T1 = nib.load(files.T1)
    masks = nib.load(files.final_contr).get_data()
    lines = np.linalg.norm(np.gradient(masks), axis=0) > 0
    print(lines.shape)
    viewer = NiftiViewer(T1.get_data(), T1.affine)
    cmap = pl.cm.jet
    my_cmap = cmap(np.arange(cmap.N))
    my_cmap[:,-1] = np.linspace(0, 1, cmap.N)
    my_cmap = ListedColormap(my_cmap)
    viewer.add_overlay(lines, cmap=my_cmap)
    viewer.show() 

def test_args(self):
        dx = np.cumsum(np.ones(5))
        dx_uneven = [1., 2., 5., 9., 11.]
        f_2d = np.arange(25).reshape(5, 5)

        # distances must be scalars or have size equal to gradient[axis]
        gradient(np.arange(5), 3.)
        gradient(np.arange(5), np.array(3.))
        gradient(np.arange(5), dx)
        # dy is set equal to dx because scalar
        gradient(f_2d, 1.5)
        gradient(f_2d, np.array(1.5))

        gradient(f_2d, dx_uneven, dx_uneven)
        # mix between even and uneven spaces and
        # mix between scalar and vector
        gradient(f_2d, dx, 2)

        # 2D but axis specified
        gradient(f_2d, dx, axis=1)

        # 2d coordinate arguments are not yet allowed
        assert_raises_regex(ValueError, '.*scalars or 1d',
            gradient, f_2d, np.stack([dx]*2, axis=-1), 1) 

def test_specific_axes(self):
        # Testing that gradient can work on a given axis only
        v = [[1, 1], [3, 4]]
        x = np.array(v)
        dx = [np.array([[2., 3.], [2., 3.]]),
              np.array([[0., 0.], [1., 1.]])]
        assert_array_equal(gradient(x, axis=0), dx[0])
        assert_array_equal(gradient(x, axis=1), dx[1])
        assert_array_equal(gradient(x, axis=-1), dx[1])
        assert_array_equal(gradient(x, axis=(1, 0)), [dx[1], dx[0]])

        # test axis=None which means all axes
        assert_almost_equal(gradient(x, axis=None), [dx[0], dx[1]])
        # and is the same as no axis keyword given
        assert_almost_equal(gradient(x, axis=None), gradient(x))

        # test vararg order
        assert_array_equal(gradient(x, 2, 3, axis=(1, 0)),
                           [dx[1]/2.0, dx[0]/3.0])
        # test maximal number of varargs
        assert_raises(TypeError, gradient, x, 1, 2, axis=1)

        assert_raises(np.AxisError, gradient, x, axis=3)
        assert_raises(np.AxisError, gradient, x, axis=-3)
        # assert_raises(TypeError, gradient, x, axis=[1,]) 

def bullen(self, depth):
        """
        Returns the Bullen parameter only for significant arrays
        """
        assert(len(depth) > 3)
        v_phi = self.v_phi(depth)
        density = self.density(depth)
        phi = v_phi * v_phi
        kappa = phi * density
        try:
            dkappadP = np.gradient(kappa, edge_order=2) / \
                       np.gradient(self.pressure(depth), edge_order=2)
            dphidz = np.gradient(phi,
                                 edge_order=2) / np.gradient(depth,
                                                             edge_order=2) / self.gravity(depth)
        except:
            dkappadP = np.gradient(kappa) / np.gradient(self.pressure(depth))
            dphidz = np.gradient(phi) / np.gradient(depth) / self.gravity(depth)
        bullen = dkappadP - dphidz
        return bullen 

def bullen(self):
        """
        Returns the Bullen parameter across the layer. 
        The Bullen parameter assess if compression as a function of pressure is 
        like homogeneous, adiabatic compression. 
        Bullen parameter =1  , homogeneous, adiabatic compression
        Bullen parameter > 1 , more compressed with pressure, e.g. across phase transitions
        Bullen parameter < 1, less compressed with pressure, e.g. across a boundary layer
        """
        kappa = self.bulk_sound_velocity * self.bulk_sound_velocity * self.density
        phi = self.bulk_sound_velocity * self.bulk_sound_velocity
        try:
            dkappadP = np.gradient(kappa, edge_order=2) / \
                       np.gradient(self.pressures, edge_order=2)
            dphidr = np.gradient(phi,edge_order=2) / np.gradient(self.radii,edge_order=2) / self.gravity
        except:
            dkappadP = np.gradient(kappa) / \
                       np.gradient(self.pressures)
            dphidr = np.gradient(phi) / np.gradient(self.radii) / self.gravity
        bullen = dkappadP + dphidr
        return bullen 

def Hillshade(raster_file, azimuth, angle_altitude): 
    
    array = ReadRasterArrayBlocks(raster_file,raster_band=1)    
    
    x, y = np.gradient(array)
    slope = np.pi/2. - np.arctan(np.sqrt(x*x + y*y))
    aspect = np.arctan2(-x, y)
    azimuthrad = np.azimuth*np.pi / 180.
    altituderad = np.angle_altitude*np.pi / 180.
     
 
    shaded = np.sin(altituderad) * np.sin(slope)\
     + np.cos(altituderad) * np.cos(slope)\
     * np.cos(azimuthrad - aspect)
    return 255*(shaded + 1)/2
#============================================================================== 

def test_specific_axes(self):
        # Testing that gradient can work on a given axis only
        v = [[1, 1], [3, 4]]
        x = np.array(v)
        dx = [np.array([[2., 3.], [2., 3.]]),
              np.array([[0., 0.], [1., 1.]])]
        assert_array_equal(gradient(x, axis=0), dx[0])
        assert_array_equal(gradient(x, axis=1), dx[1])
        assert_array_equal(gradient(x, axis=-1), dx[1])
        assert_array_equal(gradient(x, axis=(1, 0)), [dx[1], dx[0]])

        # test axis=None which means all axes
        assert_almost_equal(gradient(x, axis=None), [dx[0], dx[1]])
        # and is the same as no axis keyword given
        assert_almost_equal(gradient(x, axis=None), gradient(x))

        # test vararg order
        assert_array_equal(gradient(x, 2, 3, axis=(1, 0)), [dx[1]/2.0, dx[0]/3.0])
        # test maximal number of varargs
        assert_raises(SyntaxError, gradient, x, 1, 2, axis=1)

        assert_raises(ValueError, gradient, x, axis=3)
        assert_raises(ValueError, gradient, x, axis=-3)
        assert_raises(TypeError, gradient, x, axis=[1,]) 

def test_specific_axes(self):
        # Testing that gradient can work on a given axis only
        v = [[1, 1], [3, 4]]
        x = np.array(v)
        dx = [np.array([[2., 3.], [2., 3.]]),
              np.array([[0., 0.], [1., 1.]])]
        assert_array_equal(gradient(x, axis=0), dx[0])
        assert_array_equal(gradient(x, axis=1), dx[1])
        assert_array_equal(gradient(x, axis=-1), dx[1])
        assert_array_equal(gradient(x, axis=(1, 0)), [dx[1], dx[0]])

        # test axis=None which means all axes
        assert_almost_equal(gradient(x, axis=None), [dx[0], dx[1]])
        # and is the same as no axis keyword given
        assert_almost_equal(gradient(x, axis=None), gradient(x))

        # test vararg order
        assert_array_equal(gradient(x, 2, 3, axis=(1, 0)), [dx[1]/2.0, dx[0]/3.0])
        # test maximal number of varargs
        assert_raises(SyntaxError, gradient, x, 1, 2, axis=1)

        assert_raises(ValueError, gradient, x, axis=3)
        assert_raises(ValueError, gradient, x, axis=-3)
        assert_raises(TypeError, gradient, x, axis=[1,]) 

def test_args(self):
        dx = np.cumsum(np.ones(5))
        dx_uneven = [1., 2., 5., 9., 11.]
        f_2d = np.arange(25).reshape(5, 5)

        # distances must be scalars or have size equal to gradient[axis]
        gradient(np.arange(5), 3.)
        gradient(np.arange(5), np.array(3.))
        gradient(np.arange(5), dx)
        # dy is set equal to dx because scalar
        gradient(f_2d, 1.5)
        gradient(f_2d, np.array(1.5))

        gradient(f_2d, dx_uneven, dx_uneven)
        # mix between even and uneven spaces and
        # mix between scalar and vector
        gradient(f_2d, dx, 2)

        # 2D but axis specified
        gradient(f_2d, dx, axis=1)

        # 2d coordinate arguments are not yet allowed
        assert_raises_regex(ValueError, '.*scalars or 1d',
            gradient, f_2d, np.stack([dx]*2, axis=-1), 1) 

def test_specific_axes(self):
        # Testing that gradient can work on a given axis only
        v = [[1, 1], [3, 4]]
        x = np.array(v)
        dx = [np.array([[2., 3.], [2., 3.]]),
              np.array([[0., 0.], [1., 1.]])]
        assert_array_equal(gradient(x, axis=0), dx[0])
        assert_array_equal(gradient(x, axis=1), dx[1])
        assert_array_equal(gradient(x, axis=-1), dx[1])
        assert_array_equal(gradient(x, axis=(1, 0)), [dx[1], dx[0]])

        # test axis=None which means all axes
        assert_almost_equal(gradient(x, axis=None), [dx[0], dx[1]])
        # and is the same as no axis keyword given
        assert_almost_equal(gradient(x, axis=None), gradient(x))

        # test vararg order
        assert_array_equal(gradient(x, 2, 3, axis=(1, 0)),
                           [dx[1]/2.0, dx[0]/3.0])
        # test maximal number of varargs
        assert_raises(TypeError, gradient, x, 1, 2, axis=1)

        assert_raises(np.AxisError, gradient, x, axis=3)
        assert_raises(np.AxisError, gradient, x, axis=-3)
        # assert_raises(TypeError, gradient, x, axis=[1,]) 

def get_external_state(self):
        """
        Attempt to include avg decomp. of original normalised spread
        """
        x_sma = np.stack(
            [
                feature.get(size=self.p.time_dim) for feature in self.data.features
            ],
            axis=-1
        )
        scale = 1 / np.clip(self.data.std[0], 1e-10, None)
        x_sma *= scale  # <-- more or less ok

        # Gradient along features axis:
        dx = np.gradient(x_sma, axis=-1)

        # TODO: different conv. encoders for these two types of features:
        x = np.concatenate([x_sma, dx], axis=-1)

        # Crop outliers:
        x = np.clip(x, -10, 10)
        return x[:, None, :] 

def get_data_model_state(self):
        """
         Spread stochastic model parameters.
        """
        state = self.data_model.s.process.get_state()
        cross_corr = cov2corr(state.filtered.covariance)[[0, 0, 1], [1, 2, 2]]
        update = np.concatenate(
            [
                state.filtered.mean.flatten(),
                state.filtered.variance.flatten(),
                cross_corr,
            ]
        )
        self.external_model_state = np.concatenate(
            [
                self.external_model_state[1:, :, :],
                update[None, None, :]
            ],
            axis=0
        )
        # self.external_model_state = np.gradient(self.external_model_state, axis=-1)
        return self.external_model_state 

def get_external_state(self):

        x = np.stack(
            [
                np.frombuffer(self.data.open.get(size=self.time_dim)),
                np.frombuffer(self.data.sma_4.get(size=self.time_dim)),
                np.frombuffer(self.data.sma_8.get(size=self.time_dim)),
                np.frombuffer(self.data.sma_16.get(size=self.time_dim)),
                np.frombuffer(self.data.sma_32.get(size=self.time_dim)),
                np.frombuffer(self.data.sma_64.get(size=self.time_dim)),
                np.frombuffer(self.data.sma_128.get(size=self.time_dim)),
                np.frombuffer(self.data.sma_256.get(size=self.time_dim)),
            ],
            axis=-1
        )
        # Gradient along features axis:
        x = np.gradient(x, axis=1) * self.p.state_ext_scale

        # Log-scale:
        x = log_transform(x)
        return x[:, None, :] 

def get_external_state(self):

        x_sma = np.stack(
            [
                np.frombuffer(self.data.sma_16.get(size=self.time_dim)),
                np.frombuffer(self.data.sma_32.get(size=self.time_dim)),
                np.frombuffer(self.data.sma_64.get(size=self.time_dim)),
                np.frombuffer(self.data.sma_128.get(size=self.time_dim)),
                np.frombuffer(self.data.sma_256.get(size=self.time_dim)),
            ],
            axis=-1
        )
        # Gradient along features axis:
        dx = np.gradient(x_sma, axis=-1) * self.p.state_ext_scale

        x = tanh(dx)

        return x[:, None, :] 

def get_external_state(self):
        x_sma = np.stack(
            [
                np.frombuffer(self.data.sma_16.get(size=self.time_dim)),
                np.frombuffer(self.data.sma_32.get(size=self.time_dim)),
                np.frombuffer(self.data.sma_64.get(size=self.time_dim)),
                np.frombuffer(self.data.sma_128.get(size=self.time_dim)),
                np.frombuffer(self.data.sma_256.get(size=self.time_dim)),
            ],
            axis=-1
        )
        # Gradient along features axis:
        diff = np.gradient(x_sma, axis=-1) * self.p.state_ext_scale
        diff = tanh(diff)
        avg = np.gradient(x_sma, axis=0) * self.p.state_ext_scale
        avg = tanh(avg)

        return {'avg': avg[:, None, :], 'diff': diff[:, None, :]} 

def get_external_state(self):
        # Use Hi-Low median as signal:
        x = (
            np.frombuffer(self.data.high.get(size=self.time_dim)) +
            np.frombuffer(self.data.low.get(size=self.time_dim))
        ) / 2

        # Differences along time dimension:
        d_x = np.gradient(x, axis=0) * self.p.cwt_signal_scale

        # Compute continuous wavelet transform using Ricker wavelet:
        cwt_x = signal.cwt(d_x, signal.ricker, self.cwt_width).T

        # Note: differences taken once again along channels axis,
        # apply weighted scaling to normalize channels
        norm_x = np.gradient(cwt_x, axis=-1)
        norm_x = zscore(norm_x, axis=0) * self.p.state_ext_scale
        #out_x = tanh(norm_x)
        out_x = np.clip(norm_x, -10, 10)

        return out_x[:, None, :] 

def get_hessian(ccm, hes_norm=True, hes_smth=False, **kwargs):
    """ Find Hessian of the input cross correlation matrix <ccm>

    Parameters
    ----------
    ccm : 2D numpy array, cross-correlation matrix
    hes_norm : bool, normalize Hessian by AVG and STD?
    hes_smth : bool, smooth Hessian?

    """
    if hes_smth:
        ccm2 = nd.filters.gaussian_filter(ccm, 1)
    else:
        ccm2 = ccm
    # Jacobian components
    dcc_dy, dcc_dx = np.gradient(ccm2)
    # Hessian components
    d2cc_dx2 = np.gradient(dcc_dx)[1]
    d2cc_dy2 = np.gradient(dcc_dy)[0]
    hes = np.hypot(d2cc_dx2, d2cc_dy2)
    if hes_norm:
        hes = (hes - np.median(hes)) / np.std(hes)

    return hes 

def RdS_gradient(Frame_List, Prior_List, embed_mask_List, entropy="simple", norm_reg=True, **kwargs):
    #The Jacobian_Factor is already part of the entropy gradient that this function calls
    N_frame = Frame_List.shape[0]
    S_List = [static_regularizer(np.array([Frame_List[j]]), np.array([Prior_List[j]]), np.array([embed_mask_List[j]]), Prior_List[0].total_flux(), Prior_List[0].psize, entropy=entropy, norm_reg=norm_reg, **kwargs) for j in range(len(Frame_List))]
    S_grad_List = np.array([static_regularizer_gradient(np.array([Frame_List[j]]), np.array([Prior_List[j]]), np.array([embed_mask_List[j]]), Prior_List[0].total_flux(), Prior_List[0].psize, entropy=entropy, norm_reg=norm_reg, **kwargs) for j in range(len(Frame_List))])

    grad = np.copy(S_grad_List)*0.0

    for i in range(1,N_frame):
        grad[i] = grad[i] + 2.0*(S_List[i] - S_List[i-1])*S_grad_List[i]

    for i in range(N_frame-1):
        grad[i] = grad[i] + 2.0*(S_List[i] - S_List[i+1])*S_grad_List[i]

    return np.concatenate([grad[i] for i in range(N_frame)])


######## Rdt, RdI, and Rflow master functions ######## 

def test_args(self):
        dx = np.cumsum(np.ones(5))
        dx_uneven = [1., 2., 5., 9., 11.]
        f_2d = np.arange(25).reshape(5, 5)

        # distances must be scalars or have size equal to gradient[axis]
        gradient(np.arange(5), 3.)
        gradient(np.arange(5), np.array(3.))
        gradient(np.arange(5), dx)
        # dy is set equal to dx because scalar
        gradient(f_2d, 1.5)
        gradient(f_2d, np.array(1.5))

        gradient(f_2d, dx_uneven, dx_uneven)
        # mix between even and uneven spaces and
        # mix between scalar and vector
        gradient(f_2d, dx, 2)

        # 2D but axis specified
        gradient(f_2d, dx, axis=1)

        # 2d coordinate arguments are not yet allowed
        assert_raises_regex(ValueError, '.*scalars or 1d',
            gradient, f_2d, np.stack([dx]*2, axis=-1), 1) 

def test_specific_axes(self):
        # Testing that gradient can work on a given axis only
        v = [[1, 1], [3, 4]]
        x = np.array(v)
        dx = [np.array([[2., 3.], [2., 3.]]),
              np.array([[0., 0.], [1., 1.]])]
        assert_array_equal(gradient(x, axis=0), dx[0])
        assert_array_equal(gradient(x, axis=1), dx[1])
        assert_array_equal(gradient(x, axis=-1), dx[1])
        assert_array_equal(gradient(x, axis=(1, 0)), [dx[1], dx[0]])

        # test axis=None which means all axes
        assert_almost_equal(gradient(x, axis=None), [dx[0], dx[1]])
        # and is the same as no axis keyword given
        assert_almost_equal(gradient(x, axis=None), gradient(x))

        # test vararg order
        assert_array_equal(gradient(x, 2, 3, axis=(1, 0)),
                           [dx[1]/2.0, dx[0]/3.0])
        # test maximal number of varargs
        assert_raises(TypeError, gradient, x, 1, 2, axis=1)

        assert_raises(np.AxisError, gradient, x, axis=3)
        assert_raises(np.AxisError, gradient, x, axis=-3)
        # assert_raises(TypeError, gradient, x, axis=[1,]) 

def test_basic(self):
        v = [[1, 1], [3, 4]]
        x = np.array(v)
        dx = [np.array([[2., 3.], [2., 3.]]),
              np.array([[0., 0.], [1., 1.]])]
        assert_array_equal(gradient(x), dx)
        assert_array_equal(gradient(v), dx) 

def test_badargs(self):
        f_2d = np.arange(25).reshape(5, 5)
        x = np.cumsum(np.ones(5))

        # wrong sizes
        assert_raises(ValueError, gradient, f_2d, x, np.ones(2))
        assert_raises(ValueError, gradient, f_2d, 1, np.ones(2))
        assert_raises(ValueError, gradient, f_2d, np.ones(2), np.ones(2))
        # wrong number of arguments
        assert_raises(TypeError, gradient, f_2d, x)
        assert_raises(TypeError, gradient, f_2d, x, axis=(0,1))
        assert_raises(TypeError, gradient, f_2d, x, x, x)
        assert_raises(TypeError, gradient, f_2d, 1, 1, 1)
        assert_raises(TypeError, gradient, f_2d, x, x, axis=1)
        assert_raises(TypeError, gradient, f_2d, 1, 1, axis=1) 

def test_datetime64(self):
        # Make sure gradient() can handle special types like datetime64
        x = np.array(
            ['1910-08-16', '1910-08-11', '1910-08-10', '1910-08-12',
             '1910-10-12', '1910-12-12', '1912-12-12'],
            dtype='datetime64[D]')
        dx = np.array(
            [-5, -3, 0, 31, 61, 396, 731],
            dtype='timedelta64[D]')
        assert_array_equal(gradient(x), dx)
        assert_(dx.dtype == np.dtype('timedelta64[D]')) 

def test_masked(self):
        # Make sure that gradient supports subclasses like masked arrays
        x = np.ma.array([[1, 1], [3, 4]],
                        mask=[[False, False], [False, False]])
        out = gradient(x)[0]
        assert_equal(type(out), type(x))
        # And make sure that the output and input don't have aliased mask
        # arrays
        assert_(x.mask is not out.mask)
        # Also check that edge_order=2 doesn't alter the original mask
        x2 = np.ma.arange(5)
        x2[2] = np.ma.masked
        np.gradient(x2, edge_order=2)
        assert_array_equal(x2.mask, [False, False, True, False, False]) 

def test_timedelta64(self):
        # Make sure gradient() can handle special types like timedelta64
        x = np.array(
            [-5, -3, 10, 12, 61, 321, 300],
            dtype='timedelta64[D]')
        dx = np.array(
            [2, 7, 7, 25, 154, 119, -21],
            dtype='timedelta64[D]')
        assert_array_equal(gradient(x), dx)
        assert_(dx.dtype == np.dtype('timedelta64[D]')) 

def test_inexact_dtypes(self):
        for dt in [np.float16, np.float32, np.float64]:
            # dtypes should not be promoted in a different way to what diff does
            x = np.array([1, 2, 3], dtype=dt)
            assert_equal(gradient(x).dtype, np.diff(x).dtype) 

def test_values(self):
        # needs at least 2 points for edge_order ==1
        gradient(np.arange(2), edge_order=1)
        # needs at least 3 points for edge_order ==1
        gradient(np.arange(3), edge_order=2)

        assert_raises(ValueError, gradient, np.arange(0), edge_order=1)
        assert_raises(ValueError, gradient, np.arange(0), edge_order=2)
        assert_raises(ValueError, gradient, np.arange(1), edge_order=1)
        assert_raises(ValueError, gradient, np.arange(1), edge_order=2)
        assert_raises(ValueError, gradient, np.arange(2), edge_order=2)