Sensitivity and its normalization, optional max sens value or layer thickness

Hello Paul,

This is Luca Peruzzo, from Padova.

How are you?

I wanted to ask you about the sensitivity and its normalization.

In particular,

1. I would like to make the normalization by the maximum value optional, keeping it as default for backward compatibility.

2. I would like to add an alternative normalization by the thickness of the layer, so that the sens profile does not depend on the thickness of the layers.

I made the changes that should reflect these two points. I don't think I can submit a pull request, but the changes are so limited that I can paste them here as a first preview.

Thank you in advance for your time and for EMagPy.

Best, Luca

    def computeSens(self, forwardModel='CS', coils=None, models=[], depths=[], normalise='max'):
        """Compute normalised local sensitivity using perturbation method.
        
        Parameters
        ----------
        forwardModel : str, optional
            Type of forward model:
                - CS : Cumulative sensitivity (default)
                - FS : Full Maxwell solution with low-induction number (LIN) approximation
                - FSeq : Full Maxwell solution without LIN approximation (see Andrade 2016)
        coils : list of str, optional
            If `None`, then the default attribute of the object will be used (foward
            mode on inverted solution).
            If specified, the coil spacing, orientation and height above the ground
            will be set. In this case you need to assign at models and depths (full forward mode).
            The ECa values generated will be incorporated as a new Survey object.
        models : list of numpy.array of float
            List of array of shape Nsample x Nlayer with conductiivty in mS/m. If empty,
            `self.models` will be used.
        depths : list of numpy.array of float
            List of array of shape Nsample x (Nlayer - 1) with the depth (positive number)
            of the bottom of the layer in meters relative to the surface.
            If empty `self.depths` will be used.
        normalise : string, 'max' or 'thickness'
            Define which type of normalisation should be applied to the calculated sensitivities.
            Default is 'max' for back maintain the previous behavior
            'max' sensitivities are divided by the maximum sensitivity values, i.e., range 0 - 1.
            'thickness' each sensitivity is divided by the thickness of the associated layer.
            Values larger than 1000 mS/m / S/m can be expected with the thickness normalisation.

        Returns
        -------
        senss : list of numpy.array of float
            List of matrix of size Nsample x Ncoils x Nlayers containing the normalised
            local sensitivity.
        """
        if len(models) == 0:
            models = self.models
        if len(depths) == 0:
            depths = self.depths
            
        senss = []
        for model, depth in zip(models, depths): # for each model
            npos = model.shape[0]  
            nlayer = depth.shape[1] + 1 # number of layer (last layer is infinite)
            nprofile = nlayer + 1 # number of 1D profile (last profile is reference)
            smodels = np.dstack([model]*nprofile) # Nsample x Nlayer x Nprofile
            ix = np.arange(nlayer)
            smodels[:,ix,ix] = smodels[:,ix,ix] + 1 # perturbation
            sdepths = np.dstack([depth]*nprofile)
            lmodels = [smodels[i,:,:].T for i in range(npos)]
            ldepths = [sdepths[i,:,:].T for i in range(npos)]
            dfs = self.forward(forwardModel=forwardModel, coils=coils,
                               models=lmodels, depths=ldepths,
                               noise=0.0)
            eca = np.dstack([df.values for df in dfs]) # Nsample x Ncoils x Nprofiles
            # sens = eca[:-1,:,:] / eca[-1,:,:][None,:,:] - 1 # dividing by ref (undisturbed ECa)
            sens = eca[:-1,:,:] - eca[-1,:,:][None,:,:] # subtracting the ref ECa (slighly better up to 1e-16)

            if normalise == 'max':
                sens = sens/np.max(sens, axis=0)[None,:,:] # normalising

            elif normalise == 'thickness':
                thick_nofirst = np.diff(depth, axis=1)
                thick_withfirst = np.hstack([
                    depth[:, 0].reshape(-1, 1),  # add first depth as first thickness
                    thick_nofirst,  # all the thicknesses from the depth differences
                    np.ones((1, 1))  # add one for the last (infinite) layer
                ]).reshape(-1, 1, sens.shape[2])  # reshape to broadcast with sens (nlayer, ncoils, npos)
                sens = sens / thick_withfirst  # normalize the sens of each layer by its thickness

            senss.append(sens)
        return senss

The changes are limited to the last part of the method, the if and elif statements, and the normalise argument.