Source code for hydrac.inversion.inversion

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Inversion schemes (:mod:`hydrac.inversion.inversion`)
=====================================================

.. autoclass:: Inversion
   :members:
   :private-members:

"""

from hydrac.model.scattering.scattering import Scattering
from hydrac.model.particle import Particle
# from .run_inversion_attenuation import Run_inversion_attenuation
# from .run_inversion_implicit import Run_inversion_implicit
from .run_inversion_NNLS import Run_inversion_NNLS
from hydrac.util.parameters import Parameters
from .run_inversion_implicit import Run_inversion_implicit
import numpy as np


[docs]class Inversion(object): """Inversion class This class assimilates hydroacoustic objects and invert them according to the specifyed method and desired scattering model. The inversion consists in determining both the size and the concentration of suspended particles. There are as many unknowns as size classes in case the particles follow a size distribution. Two kinds of inversions are thus possible : Singlefrequency inversions : to invert the hydroacoustic data, the user is prompted with a mean diameter and a distribution width to constrain the backscattering model by making it function of the mean diameter. This reduces the problem to two unknowns, the mean size and the concentration. The user is thus prompted with a mean diameter which makes it possible to invert the data. Multifrequency inversions : The ratio of the backscattering instensities of the same volume at different frequencies is function of the mean diameter. Hence there is no need to indicate a mean diameter, neither a distribution width, especially for the methods returning the concentration volume distribution (NNLS, Regularization...). However, the user has still the possibility to indicate a distribution width, which is necessary for the multifrequency implicit inversion. The inversion results are stored in the HAC object and contain various information that can vary from one method to the other, but always contain the concentration values. Parameters ---------- o_hac : Hydroacoustic object derived from the :class:`hydrac.instruments.acoustique.Acoustic` class. Contains all the data and acquisition settings of the acoustic data to be inverted particletype : str, {'TFS', 'Thorne2008', 'copepDWBA', ...} Determines the scattering model to be used to invert the data. This str is used as input of the :class:`hydrac.model.particle.Particle class, itself input to the class :class:`hydrac.model.scattering.scattering.Scattering` (see :mod:`hydrac.model.scattering.scattering`) method = None : str, {'NNLS', 'implicit', ...} Inversion method to be used. If none, the user is prompted with an answer. The possible inversion methods differ for the single and multi frequency instruments : .. list-table:: Inversion methods :widths: 25 25 :header-rows: 1 * - Singlefrequency - Multifrequency * - implicit (mean size prompted) - implicit * - (-) - Non-negetive Least Square Algorithm (NNLS) * - (-) - Attenuation * - (-) - Regularization am, aM : float Minimum, maximum size of the particle to be considered (sphere radius for spherical particles, semi-minor axis for prolate spheroids, section radius for cylindical shaped particles). spacing : str {lin, log} Spacing of the size vector to be considered (linear or log-space). This vector is mostly used for the regularization or NNLS methods. dsigma = None : float Size distribution width, for averaging of the scattering model over the size distribution to make the form-function or backscattering corss-section function of the mean radius. See https://doi.org/10.1121/1.3242374. Set to None by default and prompted in case of singlefrequency inversion. a0 = None : float Mean diameter of the particles if known. Set to None by default for multifrequency inversion. Prompted for single-frequency inversions. """ def __init__(self, o_hac, particletype, method=None, am=3e-6, aM=500e-6, spacing='log', dsigma=None, a0=None, **kwargs): __acceptable_keys_list = ['depr', 'display'] method_pool = {'multi': ['NNLS', 'Implicit', 'Attenuation'], 'mono': [ 'Implicit']} [self.__setattr__(k, kwargs[k]) if k in kwargs.keys( ) else self.__setattr__(k, None ) for k in __acceptable_keys_list] self.o_hac = o_hac self.particletype = particletype self.am = am self.aM = aM self.a0 = a0 self.spacing = spacing self.dsigma = dsigma self.particle = Particle(self.particletype) if method is None: self.method = input('Please enter a method amonsgt the ' + 'following list' + str(method_pool[ self.o_hac.tag_M]) + ':') else: self.method = method if self.method == 'NNLS': Run_inversion_NNLS(self, depr=self.depr) self._reshape_solution() self._bv_to_conc() elif self.method == 'Implicit': if self.dsigma is None: self.dsigma = np.float(input( 'Please enter a distribution width: ')) if self.o_hac.tag_M == 'mono': if self.a0 is None: self.a0 = np.float(input( 'Please enter a median radius: ')) self.am = self.a0-.1*self.a0 self.aM = self.a0+.1*self.a0 Run_inversion_implicit(self, a0=self.a0) else: if self.a0 is not None: # temporary... self.am = self.a0 - .1 * self.a0 self.aM = self.a0 + .1 * self.a0 Run_inversion_implicit(self, a0=self.a0) elif self.method == 'Attenuation': pass self.o_hac.method_inv = self.method def calc_inv(self, func, obj): if hasattr(self.o_hac, 'single_model')\ and self.o_hac.single_model is True: single_scattering = Scattering( Particle(self.particletype), water=obj.P0.w_obj, freq=obj.P0.FrequencyRx, tag_M=self.o_hac.tag_M, mode_depl=self.o_hac.mode_depl, am=self.am, aM=self.aM, spacing=self.spacing, dsigma=self.dsigma) [obj['P' + str(i)].update({'inv': func(obj['P' + str(i)].asv, single_scattering), 'S':single_scattering.scattering}) for i in range(0, len(obj))] else: for i in range(0, len(obj)): s = Scattering(Particle(self.particletype), water=obj['P' + str(i)].w_obj, freq=obj['P' + str(i)].FrequencyRx, tag_M=self.o_hac.tag_M, mode_depl=self.o_hac.mode_depl, am=self.am, aM=self.aM, spacing=self.spacing, dsigma=self.dsigma) obj['P' + str(i)].update({'inv': func(obj['P' + str(i)].asv, s), 'S': s.scattering}) self.o_hac.inverted = 1 def _reshape_solution(self): if hasattr(self.o_hac, 'inverted'): [self.o_hac.param['P' + str(i)].update({ 'solution': np.squeeze(np.array( [h.solution for h in self.o_hac.param['P' + str(i)].inv])), 'bv':np.squeeze(np.array([ h.bv for h in self.o_hac.param['P' + str(i)].inv])), 'norm':np.squeeze(np.array([ h.norm for h in self.o_hac.param['P' + str(i)].inv])), 'rnorm':np.squeeze(np.array([ h.rnorm for h in self.o_hac.param['P' + str(i)].inv]))} ) for i in range(0, len( self.o_hac.param))] self.reshape_sol = 1 [self.o_hac.param['P' + str(i)].update({ 'inv': Parameters({'bv': self.o_hac.param['P' + str(i)].bv, 'norm': self.o_hac.param['P' + str(i)].norm, 'rnorm': self.o_hac.param['P' + str(i) ].rnorm})} ) for i in range(0, len( self.o_hac.param))] [[self.o_hac.param['P' + str(i)].__delattr__(k) for k in ['bv', 'rnorm', 'norm']] for i in range(0, len(self.o_hac.param))] def _bv_to_conc(self): if hasattr(self.o_hac, 'inverted') and hasattr(self, 'reshape_sol'): # deep copy bv to manipulate it into conc ?? for i in range(0, len(self.o_hac.param)): if len(np.shape(self.o_hac.param['P' + str(i)].inv.bv)) == 3\ and np.size(self.o_hac.param['P' + str(i)].inv.bv, 2) != len(self.o_hac.param['P' + str(i)].S.a_UND): self.o_hac.param['P' + str(i)].inv.bv = np.transpose( self.o_hac.param['P' + str(i)].inv.bv, (2, 0, 1)) self.o_hac.param['P' + str(i)].inv.norm = np.transpose( self.o_hac.param['P' + str(i)].inv.norm, (2, 0, 1)) self.o_hac.param['P' + str(i)].inv.update({'conc_i': np.zeros(np.shape(self.o_hac.param['P' + str(i)].inv.bv))}) if self.particle.nature != 'HYB': self.o_hac.param['P' + str(i)].inv.conc_i =\ self.o_hac.param['P' + str(i)].inv.bv *\ self.particle.param.rhos else: if len(np.shape(self.o_hac.param['P' + str(i)].inv.bv)) == 3: a_HR = self.o_hac.param['P' + str(i)].S.a_HR a_UND = self.o_hac.param['P' + str(i)].S.a_UND for j in range(0, np.size(self.o_hac.param['P' + str(i)].inv.bv, 2)): ind = np.where(np.abs(a_HR-a_UND[j]) == np.min(np.abs(a_HR-a_UND[j])) )[0][0] self.o_hac.param['P' + str(i)].inv.conc_i[:, :, j] =\ self.o_hac.param['P' + str(i)].inv.bv[:, :, j] *\ self.o_hac.param['P' + str(i)].S.rhoex[ind] self.o_hac.param['P' + str(i) ].inv.update({'conc': np.nansum(self.o_hac.param['P' + str(i)].inv.conc_i, len(np.shape(self.o_hac.param['P' + str(i)].inv.bv))-1)})