code: # -*- coding: utf-8 -*-
  #pylint: disable-msg=E0611, E1101, C0103, R0901, R0902, R0903, R0904, W0232
  #------------------------------------------------------------------------------
  # Copyright (c) 2007-2021, Acoular Development Team.
  #------------------------------------------------------------------------------
  """Implements beamformers in the frequency domain.
  
  .. autosummary::
      :toctree: generated/
      
      SteeringVector
  
      
      BeamformerBase
      BeamformerFunctional
      BeamformerCapon
      BeamformerEig
      BeamformerMusic
      BeamformerClean
      BeamformerDamas
      BeamformerDamasPlus
      BeamformerOrth
      BeamformerCleansc
      BeamformerCMF
      BeamformerSODIX
      BeamformerGIB
  
      PointSpreadFunction
      L_p
      integrate
  
  """
  
  # imports from other packages
  from __future__ import print_function, division
  
  from numpy import array, ones, full, hanning, hamming, bartlett, blackman, \
  invert, dot, newaxis, zeros, empty, fft, float32, float64, complex64, linalg, \
  where, searchsorted, pi, multiply, sign, diag, arange, sqrt, exp, log10, int,\
  reshape, hstack, vstack, eye, tril, size, clip, tile, round, delete, \
  absolute, argsort, sort, sum, hsplit, fill_diagonal, zeros_like, isclose, \
  vdot, flatnonzero, einsum, ndarray, isscalar, inf, real
  
  from numpy.linalg import norm
  
  from sklearn.linear_model import LassoLars, LassoLarsCV, LassoLarsIC,\
  OrthogonalMatchingPursuit, ElasticNet, OrthogonalMatchingPursuitCV, Lasso
  
  from scipy.optimize import nnls, linprog, fmin_l_bfgs_b
  from scipy.linalg import inv, eigh, eigvals, fractional_matrix_power
  from warnings import warn
  
  from traits.api import HasPrivateTraits, Float, Int, ListInt, ListFloat, \
  CArray, Property, Instance, Trait, Bool, Range, Delegate, Enum, Any, \
  cached_property, on_trait_change, property_depends_on
  from traits.trait_errors import TraitError
  
  from .fastFuncs import beamformerFreq, calcTransfer, calcPointSpreadFunction, \
  damasSolverGaussSeidel
  
  from .h5cache import H5cache
  from .h5files import H5CacheFileBase
  from .internal import digest
  from .grids import Grid, Sector
  from .microphones import MicGeom
  from .configuration import config
  from .environments import Environment
  from .spectra import PowerSpectra
  
  class SteeringVector( HasPrivateTraits ):
      """ 
      Basic class for implementing steering vectors with monopole source transfer models
      """
      
      #: :class:`~acoular.grids.Grid`-derived object that provides the grid locations.
      grid = Trait(Grid, 
          desc="beamforming grid")
      
      #: :class:`~acoular.microphones.MicGeom` object that provides the microphone locations.
      mics = Trait(MicGeom, 
          desc="microphone geometry")
          
      #: Type of steering vectors, see also :ref:`Sarradj, 2012<Sarradj2012>`. Defaults to 'true level'.
      steer_type = Trait('true level', 'true location', 'classic', 'inverse',
                    desc="type of steering vectors used")
      
      #: :class:`~acoular.environments.Environment` or derived object, 
      #: which provides information about the sound propagation in the medium.
      #: Defaults to standard :class:`~acoular.environments.Environment` object.
      env = Instance(Environment(), Environment)
      
      # TODO: add caching capability for transfer function
      # Flag, if "True" (not default), the transfer function is 
      # cached in h5 files and does not have to be recomputed during subsequent 
      # program runs. 
      # Be aware that setting this to "True" may result in high memory usage.
      #cached = Bool(False, 
      #              desc="cache flag for transfer function")    
      
      
      # Sound travel distances from microphone array center to grid 
      # points or reference position (readonly). Feature may change.
      r0 = Property(desc="array center to grid distances")
  
      # Sound travel distances from array microphones to grid 
      # points (readonly). Feature may change.
      rm = Property(desc="all array mics to grid distances")
      
      # mirror trait for ref
      _ref = Any(array([0.,0.,0.]),
                 desc="reference position or distance")
      
      #: Reference position or distance at which to evaluate the sound pressure 
      #: of a grid point. 
      #: If set to a scalar, this is used as reference distance to the grid points.
      #: If set to a vector, this is interpreted as x,y,z coordinates of the reference position.
      #: Defaults to [0.,0.,0.].
      ref = Property(desc="reference position or distance")
      
      def _set_ref (self, ref):
          if isscalar(ref):
              try:
                  self._ref = absolute(float(ref))
              except:
                  raise TraitError(args=self,
                                   name='ref', 
                                   info='Float or CArray(3,)',
                                   value=ref) 
          elif len(ref) == 3:
              self._ref = array(ref, dtype=float)
          else:
              raise TraitError(args=self,
                               name='ref', 
                               info='Float or CArray(3,)',
                               value=ref)
        
      def _get_ref (self):
          return self._ref
      
      
      # internal identifier
      digest = Property( 
          depends_on = ['steer_type', 'env.digest', 'grid.digest', 'mics.digest', '_ref'])
      
      # internal identifier, use for inverse methods, excluding steering vector type
      inv_digest = Property( 
          depends_on = ['env.digest', 'grid.digest', 'mics.digest', '_ref'])
          
      @property_depends_on('grid.digest, env.digest, _ref')
      def _get_r0 ( self ):
          if isscalar(self.ref) and self.ref > 0:
              return full((self.grid.size,), self.ref)
          else:
              return self.env._r(self.grid.pos())
  
      @property_depends_on('grid.digest, mics.digest, env.digest')
      def _get_rm ( self ):
          return self.env._r(self.grid.pos(), self.mics.mpos)
   
      @cached_property
      def _get_digest( self ):
          return digest( self )
      
      @cached_property
      def _get_inv_digest( self ):
          return digest( self )
      
      def transfer(self, f, ind=None):
          """
          Calculates the transfer matrix for one frequency. 
          
          Parameters
          ----------
          f   : float
              Frequency for which to calculate the transfer matrix
          ind : (optional) array of ints
              If set, only the transfer function of the gridpoints addressed by 
              the given indices will be calculated. Useful for algorithms like CLEAN-SC,
              where not the full transfer matrix is needed
          
          Returns
          -------
          array of complex128
              array of shape (ngridpts, nmics) containing the transfer matrix for the given frequency
          """
          #if self.cached:
          #    warn('Caching of transfer function is not yet supported!', Warning)
          #    self.cached = False
          
          if ind is None:
              trans = calcTransfer(self.r0, self.rm, array(2*pi*f/self.env.c))
          elif not isinstance(ind,ndarray):
              trans = calcTransfer(self.r0[ind], self.rm[ind, :][newaxis], array(2*pi*f/self.env.c))#[0, :]
          else:
              trans = calcTransfer(self.r0[ind], self.rm[ind, :], array(2*pi*f/self.env.c))
          return trans
      
      def steer_vector(self, f, ind=None):
          """
          Calculates the steering vectors based on the transfer function
          See also :ref:`Sarradj, 2012<Sarradj2012>`.
          
          Parameters
          ----------
          f   : float
              Frequency for which to calculate the transfer matrix
          ind : (optional) array of ints
              If set, only the steering vectors of the gridpoints addressed by 
              the given indices will be calculated. Useful for algorithms like CLEAN-SC,
              where not the full transfer matrix is needed
          
          Returns
          -------
          array of complex128
              array of shape (ngridpts, nmics) containing the steering vectors for the given frequency
          """
          func = {'classic' : lambda x: x / absolute(x) / x.shape[-1],
                  'inverse' : lambda x: 1. / x.conj() / x.shape[-1],
                  'true level' : lambda x: x / einsum('ij,ij->i',x,x.conj())[:,newaxis],
                  'true location' : lambda x: x / sqrt(einsum('ij,ij->i',x,x.conj()) * x.shape[-1])[:,newaxis]
                  }[self.steer_type]
          return func(self.transfer(f, ind))
      
      
  class BeamformerBase( HasPrivateTraits ):
      """
      Beamforming using the basic delay-and-sum algorithm in the frequency domain.
      """
      
      
      # Instance of :class:`~acoular.fbeamform.SteeringVector` or its derived classes
      # that contains information about the steering vector. This is a private trait.
      # Do not set this directly, use `steer` trait instead.
      _steer_obj = Instance(SteeringVector(), SteeringVector)   
      
      #: :class:`~acoular.fbeamform.SteeringVector` or derived object. 
      #: Defaults to :class:`~acoular.fbeamform.SteeringVector` object.
      steer = Property(desc="steering vector object")  
      
      def _get_steer(self):
          return self._steer_obj
      
      def _set_steer(self, steer):
          if isinstance(steer, SteeringVector):
              self._steer_obj = steer
          elif steer in ('true level', 'true location', 'classic', 'inverse'):
              # Type of steering vectors, see also :ref:`Sarradj, 2012<Sarradj2012>`.
              warn("Deprecated use of 'steer' trait. "
                   "Please use object of class 'SteeringVector' in the future.", 
                   Warning, stacklevel = 2)
              self._steer_obj.steer_type = steer
          else:
              raise(TraitError(args=self,
                               name='steer', 
                               info='SteeringVector',
                               value=steer))
  
      # --- List of backwards compatibility traits and their setters/getters -----------
      
      # :class:`~acoular.environments.Environment` or derived object. 
      # Deprecated! Only kept for backwards compatibility. 
      # Now governed by :attr:`steer` trait.
      env = Property()
      
      def _get_env(self):
          return self._steer_obj.env    
      
      def _set_env(self, env):
          warn("Deprecated use of 'env' trait. ", Warning, stacklevel = 2)
          self._steer_obj.env = env
      
      # The speed of sound.
      # Deprecated! Only kept for backwards compatibility. 
      # Now governed by :attr:`steer` trait.
      c = Property()
      
      def _get_c(self):
          return self._steer_obj.env.c
      
      def _set_c(self, c):
          warn("Deprecated use of 'c' trait. ", Warning, stacklevel = 2)
          self._steer_obj.env.c = c
     
      # :class:`~acoular.grids.Grid`-derived object that provides the grid locations.
      # Deprecated! Only kept for backwards compatibility. 
      # Now governed by :attr:`steer` trait.
      grid = Property()
  
      def _get_grid(self):
          return self._steer_obj.grid
      
      def _set_grid(self, grid):
          warn("Deprecated use of 'grid' trait. ", Warning, stacklevel = 2)
          self._steer_obj.grid = grid
      
      # :class:`~acoular.microphones.MicGeom` object that provides the microphone locations.
      # Deprecated! Only kept for backwards compatibility. 
      # Now governed by :attr:`steer` trait
      mpos = Property()
      
      def _get_mpos(self):
          return self._steer_obj.mics
      
      def _set_mpos(self, mpos):
          warn("Deprecated use of 'mpos' trait. ", Warning, stacklevel = 2)
          self._steer_obj.mics = mpos
      
      
      # Sound travel distances from microphone array center to grid points (r0)
      # and all array mics to grid points (rm). Readonly.
      # Deprecated! Only kept for backwards compatibility. 
      # Now governed by :attr:`steer` trait
      r0 = Property()
      def _get_r0(self):
          return self._steer_obj.r0
      
      rm = Property()
      def _get_rm(self):
          return self._steer_obj.rm
      
      # --- End of backwards compatibility traits --------------------------------------
  
      #: :class:`~acoular.spectra.PowerSpectra` object that provides the 
      #: cross spectral matrix and eigenvalues
      freq_data = Trait(PowerSpectra, 
                        desc="freq data object")
  
      #: Boolean flag, if 'True' (default), the main diagonal is removed before beamforming.
      r_diag = Bool(True, 
                    desc="removal of diagonal")
      
      #: If r_diag==True: if r_diag_norm==0.0, the standard  
      #: normalization = num_mics/(num_mics-1) is used. 
      #: If r_diag_norm !=0.0, the user input is used instead.  
      #: If r_diag==False, the normalization is 1.0 either way. 
      r_diag_norm = Float(0.0, 
                          desc="If diagonal of the csm is removed, some signal energy is lost." 
                          "This is handled via this normalization factor." 
                          "Internally, the default is: num_mics / (num_mics - 1).") 
      
      #: Floating point precision of property result. Corresponding to numpy dtypes. Default = 64 Bit.
      precision = Trait('float64', 'float32',
              desc="precision (32/64 Bit) of result, corresponding to numpy dtypes")
      
      #: Boolean flag, if 'True' (default), the result is cached in h5 files.
      cached = Bool(True, 
          desc="cached flag")
                    
      # hdf5 cache file
      h5f = Instance( H5CacheFileBase, transient = True )
      
      #: The beamforming result as squared sound pressure values 
      #: at all grid point locations (readonly).
      #: Returns a (number of frequencies, number of gridpoints) array of floats.
      result = Property(
          desc="beamforming result")
      
      # internal identifier
      digest = Property( 
          depends_on = ['freq_data.digest', 'r_diag', 'r_diag_norm', 'precision', '_steer_obj.digest'])
  
      # internal identifier
      ext_digest = Property( 
          depends_on = ['digest', 'freq_data.ind_low', 'freq_data.ind_high'], 
          )
  
      @cached_property
      def _get_digest( self ):
          return digest( self )
      
      @cached_property
      def _get_ext_digest( self ):
          return digest( self, 'ext_digest' )
      
      def _get_filecache( self ):
          """
          function collects cached results from file depending on 
          global/local caching behaviour. Returns (None, None) if no cachefile/data 
          exist and global caching mode is 'readonly'.
          """
  #        print("get cachefile:", self.freq_data.basename)
          H5cache.get_cache_file( self, self.freq_data.basename ) 
          if not self.h5f: 
  #            print("no cachefile:", self.freq_data.basename)
              return (None, None)# only happens in case of global caching readonly
  
          nodename = self.__class__.__name__ + self.digest
  #        print("collect filecache for nodename:",nodename)
          if config.global_caching == 'overwrite' and self.h5f.is_cached(nodename):
  #            print("remove existing data for nodename",nodename)
              self.h5f.remove_data(nodename) # remove old data before writing in overwrite mode
          
          if not self.h5f.is_cached(nodename):
  #            print("no data existent for nodename:", nodename)
              if config.global_caching == 'readonly': 
                  return (None, None)
              else:
  #                print("initialize data.")
                  numfreq = self.freq_data.fftfreq().shape[0]# block_size/2 + 1steer_obj
                  group = self.h5f.create_new_group(nodename)
                  self.h5f.create_compressible_array('result',
                                        (numfreq, self.steer.grid.size),
                                        self.precision,
                                        group)
                  self.h5f.create_compressible_array('freqs',
                                        (numfreq, ),
                                        'int8',#'bool', 
                                        group)
          ac = self.h5f.get_data_by_reference('result','/'+nodename)
          fr = self.h5f.get_data_by_reference('freqs','/'+nodename)
          return (ac,fr)        
  
      def _assert_equal_channels(self):
          numchannels = self.freq_data.numchannels
          if  numchannels != self.steer.mics.num_mics or numchannels == 0:
              raise ValueError("%i channels do not fit %i mics" % (numchannels, self.steer.mics.num_mics))        
  
      @property_depends_on('ext_digest')
      def _get_result ( self ):
          """
          This is the :attr:`result` getter routine.
          The beamforming result is either loaded or calculated.
          """
          f = self.freq_data
          numfreq = f.fftfreq().shape[0]# block_size/2 + 1steer_obj
          _digest = ''
          while self.digest != _digest:
              _digest = self.digest
              self._assert_equal_channels()
              if not ( # if result caching is active
                      config.global_caching == 'none' or 
                      (config.global_caching == 'individual' and self.cached == False)
                  ):
  #                print("get filecache..")
                  (ac,fr) = self._get_filecache() 
                  if ac and fr: 
  #                    print("cached data existent")
                      if not fr[f.ind_low:f.ind_high].all():
  #                        print("calculate missing results")                            
                          if config.global_caching == 'readonly': 
                              (ac, fr) = (ac[:], fr[:])
                          self.calc(ac,fr)
                          self.h5f.flush()
  #                    else:
  #                        print("cached results are complete! return.")
                  else:
  #                    print("no caching, calculate result")
                      ac = zeros((numfreq, self.steer.grid.size), dtype=self.precision)
                      fr = zeros(numfreq, dtype='int8')
                      self.calc(ac,fr)
              else:
  #                print("no caching activated, calculate result")
                  ac = zeros((numfreq, self.steer.grid.size), dtype=self.precision)
                  fr = zeros(numfreq, dtype='int8')
                  self.calc(ac,fr)
          return ac
        
      def sig_loss_norm(self):
          """ 
          If the diagonal of the CSM is removed one has to handle the loss 
          of signal energy --> Done via a normalization factor.
          """
          if not self.r_diag:  # Full CSM --> no normalization needed 
              normFactor = 1.0 
          elif self.r_diag_norm == 0.0:  # Removed diag: standard normalization factor 
              nMics = float(self.freq_data.numchannels) 
              normFactor = nMics / (nMics - 1) 
          elif self.r_diag_norm != 0.0:  # Removed diag: user defined normalization factor 
              normFactor = self.r_diag_norm 
          return normFactor
  
  
      def _beamformer_params(self):
          """
          Manages the parameters for calling of the core beamformer functionality.
          This is a workaround to allow faster calculation and may change in the
          future.
          
          Returns
          -------
              - String containing the steering vector type
              - Function for frequency-dependent steering vector calculation
                  
          """
          if type(self.steer) == SteeringVector: # for simple steering vector, use faster method
              param_type = self.steer.steer_type
              def param_steer_func(f): return (self.steer.r0, self.steer.rm, 2*pi*f/self.steer.env.c )
          else:
              param_type = 'custom'
              param_steer_func = self.steer.steer_vector
          return param_type, param_steer_func
  
      def calc(self, ac, fr):
          """
          Calculates the delay-and-sum beamforming result for the frequencies 
          defined by :attr:`freq_data`
          
          This is an internal helper function that is automatically called when 
          accessing the beamformer's :attr:`result` or calling
          its :meth:`synthetic` method.        
          
          Parameters
          ----------
          ac : array of floats
              This array of dimension ([number of frequencies]x[number of gridpoints])
              is used as call-by-reference parameter and contains the calculated
              value after calling this method. 
          fr : array of booleans
              The entries of this [number of frequencies]-sized array are either 
              'True' (if the result for this frequency has already been calculated)
              or 'False' (for the frequencies where the result has yet to be calculated).
              After the calculation at a certain frequency the value will be set
              to 'True'
          
          Returns
          -------
          This method only returns values through the *ac* and *fr* parameters
          """
          f = self.freq_data.fftfreq()#[inds]
          param_steer_type, steer_vector = self._beamformer_params()
          for i in self.freq_data.indices:
              if not fr[i]:
                  csm = array(self.freq_data.csm[i], dtype='complex128')
                  beamformerOutput = beamformerFreq(param_steer_type, 
                                                    self.r_diag, 
                                                    self.sig_loss_norm(), 
                                                    steer_vector(f[i]), 
                                                    csm)[0]
                  if self.r_diag:  # set (unphysical) negative output values to 0
                      indNegSign = sign(beamformerOutput) < 0
                      beamformerOutput[indNegSign] = 0.0
                  ac[i] = beamformerOutput
                  fr[i] = 1
      
      def synthetic( self, f, num=0):
          """
          Evaluates the beamforming result for an arbitrary frequency band.
          
          Parameters
          ----------
          f: float
              Band center frequency. 
          num : integer
              Controls the width of the frequency bands considered; defaults to
              0 (single frequency line).
              
              ===  =====================
              num  frequency band width
              ===  =====================
              0    single frequency line
              1    octave band
              3    third-octave band
              n    1/n-octave band
              ===  =====================
                
          Returns
          -------
          array of floats
              The synthesized frequency band values of the beamforming result at 
              each grid point .
              Note that the frequency resolution and therefore the bandwidth 
              represented by a single frequency line depends on 
              the :attr:`sampling frequency<acoular.tprocess.SamplesGenerator.sample_freq>` and 
              used :attr:`FFT block size<acoular.spectra.PowerSpectra.block_size>`.
          """
          res = self.result # trigger calculation
          freq = self.freq_data.fftfreq()
          if len(freq) == 0:
              return None
          
          indices = self.freq_data.indices
          
          if num == 0:
              # single frequency line
              ind = searchsorted(freq, f)
              if ind >= len(freq):
                  warn('Queried frequency (%g Hz) not in resolved '
                                'frequency range. Returning zeros.' % f, 
                                Warning, stacklevel = 2)
                  h = zeros_like(res[0])
              else:
                  if freq[ind] != f:
                      warn('Queried frequency (%g Hz) not in set of '
                           'discrete FFT sample frequencies. '
                           'Using frequency %g Hz instead.' % (f,freq[ind]), 
                           Warning, stacklevel = 2)
                  if not (ind in indices):
                      warn('Beamforming result may not have been calculated '
                           'for queried frequency. Check '
                           'freq_data.ind_low and freq_data.ind_high!',
                            Warning, stacklevel = 2)
                  h = res[ind]
          else:
              # fractional octave band
              if isinstance(num,list):
                  f1=num[0]
                  f2=num[-1]
              else:
                  f1 = f*2.**(-0.5/num)
                  f2 = f*2.**(+0.5/num)
              ind1 = searchsorted(freq, f1)
              ind2 = searchsorted(freq, f2)
              if ind1 == ind2:
                  warn('Queried frequency band (%g to %g Hz) does not '
                       'include any discrete FFT sample frequencies. '
                       'Returning zeros.' % (f1,f2), 
                       Warning, stacklevel = 2)
                  h = zeros_like(res[0])
              else:
                  h = sum(res[ind1:ind2], 0)
                  if not ((ind1 in indices) and (ind2 in indices)):
                      warn('Beamforming result may not have been calculated '
                           'for all queried frequencies. Check '
                           'freq_data.ind_low and freq_data.ind_high!',
                            Warning, stacklevel = 2)
          return h.reshape(self.steer.grid.shape)
  
  
      def integrate(self, sector):
          """
          Integrates result map over a given sector.
          
          Parameters
          ----------
          sector: array of floats
                Tuple with arguments for the 'indices' method 
                of a :class:`~acoular.grids.Grid`-derived class 
                (e.g. :meth:`RectGrid.indices<acoular.grids.RectGrid.indices>` 
                or :meth:`RectGrid3D.indices<acoular.grids.RectGrid3D.indices>`).
                Possible sectors would be *array([xmin, ymin, xmax, ymax])* 
                or *array([x, y, radius])*.
                
          Returns
          -------
          array of floats
              The spectrum (all calculated frequency bands) for the integrated sector.
          """
          #resp. array([rmin, phimin, rmax, phimax]), array([r, phi, radius]).
          
  #        ind = self.grid.indices(*sector)
  #        gshape = self.grid.shape
  #        r = self.result
  #        rshape = r.shape
  #        mapshape = (rshape[0], ) + gshape
  #        h = r[:].reshape(mapshape)[ (s_[:], ) + ind ]
  #        return h.reshape(h.shape[0], prod(h.shape[1:])).sum(axis=1)
          if isinstance(sector, Sector):
              ind = self.steer.grid.subdomain(sector)
          else:
              ind = self.steer.grid.indices(*sector)
          gshape = self.steer.grid.shape
          r = self.result
          h = zeros(r.shape[0])
          for i in range(r.shape[0]):
              h[i] = r[i].reshape(gshape)[ind].sum()
          return h
  
  class BeamformerFunctional( BeamformerBase ):
      """
      Functional beamforming after :ref:`Dougherty, 2014<Dougherty2014>`.
      """
  
      #: Functional exponent, defaults to 1 (= Classic Beamforming).
      gamma = Float(1, 
          desc="functional exponent")
  
      # internal identifier
      digest = Property(depends_on = ['freq_data.digest', '_steer_obj.digest', 'r_diag', 'gamma'])
      
      #: Functional Beamforming is only well defined for full CSM
      r_diag = Enum(False, 
                    desc="False, as Functional Beamformer is only well defined for the full CSM")
  
      @cached_property
      def _get_digest( self ):
          return digest( self )
  
      def calc(self, ac, fr):
          """
          Calculates the Functional Beamformer result for the frequencies defined by :attr:`freq_data`
          
          This is an internal helper function that is automatically called when 
          accessing the beamformer's :attr:`~BeamformerBase.result` or calling
          its :meth:`~BeamformerBase.synthetic` method.        
          
          Parameters
          ----------
          ac : array of floats
              This array of dimension ([number of frequencies]x[number of gridpoints])
              is used as call-by-reference parameter and contains the calculated
              value after calling this method. 
          fr : array of booleans
              The entries of this [number of frequencies]-sized array are either 
              'True' (if the result for this frequency has already been calculated)
              or 'False' (for the frequencies where the result has yet to be calculated).
              After the calculation at a certain frequency the value will be set
              to 'True'
          
          Returns
          -------
          This method only returns values through the *ac* and *fr* parameters
          """
          f = self.freq_data.fftfreq()
          normFactor = self.sig_loss_norm()
          param_steer_type, steer_vector = self._beamformer_params()
          for i in self.freq_data.indices:
              if not fr[i]:
                  if self.r_diag:
                      # This case is not used at the moment (see Trait r_diag)  
                      # It would need some testing as structural changes were not tested...
  #==============================================================================
  #                     One cannot use spectral decomposition when diagonal of csm is removed,
  #                     as the resulting modified eigenvectors are not orthogonal to each other anymore.
  #                     Therefor potentiating cannot be applied only to the eigenvalues.
  #                     --> To avoid this the root of the csm (removed diag) is calculated directly.
  #                     WATCH OUT: This doesn't really produce good results.
  #==============================================================================
                      csm = self.freq_data.csm[i]
                      fill_diagonal(csm, 0)
                      csmRoot = fractional_matrix_power(csm, 1.0 / self.gamma)
                      beamformerOutput, steerNorm = beamformerFreq(param_steer_type, 
                                                                   self.r_diag, 
                                                                   1.0, 
                                                                   steer_vector(f[i]), 
                                                                   csmRoot)
                      beamformerOutput /= steerNorm  # take normalized steering vec
                      
                      # set (unphysical) negative output values to 0
                      indNegSign = sign(beamformerOutput) < 0
                      beamformerOutput[indNegSign] = 0.0
                  else:
                      eva = array(self.freq_data.eva[i], dtype='float64') ** (1.0 / self.gamma)
                      eve = array(self.freq_data.eve[i], dtype='complex128')
                      beamformerOutput, steerNorm = beamformerFreq(param_steer_type, 
                                                                   self.r_diag, 
                                                                   1.0, 
                                                                   steer_vector(f[i]), 
                                                                   (eva, eve))
                      beamformerOutput /= steerNorm  # take normalized steering vec
                  ac[i] = (beamformerOutput ** self.gamma) * steerNorm * normFactor  # the normalization must be done outside the beamformer
                  fr[i] = 1
              
  class BeamformerCapon( BeamformerBase ):
      """
      Beamforming using the Capon (Mininimum Variance) algorithm, 
      see :ref:`Capon, 1969<Capon1969>`.
      """
      # Boolean flag, if 'True', the main diagonal is removed before beamforming;
      # for Capon beamforming r_diag is set to 'False'.
      r_diag = Enum(False, 
          desc="removal of diagonal")
  
      def calc(self, ac, fr):
          """
          Calculates the Capon result for the frequencies defined by :attr:`freq_data`
          
          This is an internal helper function that is automatically called when 
          accessing the beamformer's :attr:`~BeamformerBase.result` or calling
          its :meth:`~BeamformerBase.synthetic` method.        
          
          Parameters
          ----------
          ac : array of floats
              This array of dimension ([number of frequencies]x[number of gridpoints])
              is used as call-by-reference parameter and contains the calculated
              value after calling this method. 
          fr : array of booleans
              The entries of this [number of frequencies]-sized array are either 
              'True' (if the result for this frequency has already been calculated)
              or 'False' (for the frequencies where the result has yet to be calculated).
              After the calculation at a certain frequency the value will be set
              to 'True'
          
          Returns
          -------
          This method only returns values through the *ac* and *fr* parameters
          """        
          f = self.freq_data.fftfreq()
          nMics = self.freq_data.numchannels
          normFactor = self.sig_loss_norm() * nMics**2
          param_steer_type, steer_vector = self._beamformer_params()
          for i in self.freq_data.indices:
              if not fr[i]:
                  csm = array(linalg.inv(array(self.freq_data.csm[i], dtype='complex128')), order='C')
                  beamformerOutput = beamformerFreq(param_steer_type, 
                                                    self.r_diag, 
                                                    normFactor, 
                                                    steer_vector(f[i]), 
                                                    csm)[0]
                  ac[i] = 1.0 / beamformerOutput
                  fr[i] = 1
  
  class BeamformerEig( BeamformerBase ):
      """
      Beamforming using eigenvalue and eigenvector techniques,
      see :ref:`Sarradj et al., 2005<Sarradj2005>`.
      """
      #: Number of component to calculate: 
      #: 0 (smallest) ... :attr:`~acoular.tprocess.SamplesGenerator.numchannels`-1;
      #: defaults to -1, i.e. numchannels-1
      n = Int(-1, 
          desc="No. of eigenvalue")
  
      # Actual component to calculate, internal, readonly.
      na = Property(
          desc="No. of eigenvalue")
  
      # internal identifier
      digest = Property( 
          depends_on = ['freq_data.digest', '_steer_obj.digest', 'r_diag', 'n'])
  
      @cached_property
      def _get_digest( self ):
          return digest( self )
      
      @property_depends_on('steer.mics, n')
      def _get_na( self ):
          na = self.n
          nm = self.steer.mics.num_mics
          if na < 0:
              na = max(nm + na, 0)
          return min(nm - 1, na)
  
      def calc(self, ac, fr):
          """
          Calculates the result for the frequencies defined by :attr:`freq_data`
          
          This is an internal helper function that is automatically called when 
          accessing the beamformer's :attr:`~BeamformerBase.result` or calling
          its :meth:`~BeamformerBase.synthetic` method.        
          
          Parameters
          ----------
          ac : array of floats
              This array of dimension ([number of frequencies]x[number of gridpoints])
              is used as call-by-reference parameter and contains the calculated
              value after calling this method. 
          fr : array of booleans
              The entries of this [number of frequencies]-sized array are either 
              'True' (if the result for this frequency has already been calculated)
              or 'False' (for the frequencies where the result has yet to be calculated).
              After the calculation at a certain frequency the value will be set
              to 'True'
          
          Returns
          -------
          This method only returns values through the *ac* and *fr* parameters
          """
          f = self.freq_data.fftfreq()
          na = int(self.na)  # eigenvalue taken into account
          normFactor = self.sig_loss_norm()
          param_steer_type, steer_vector = self._beamformer_params()
          for i in self.freq_data.indices:
              if not fr[i]:
                  eva = array(self.freq_data.eva[i], dtype='float64')
                  eve = array(self.freq_data.eve[i], dtype='complex128')
                  beamformerOutput = beamformerFreq(param_steer_type, 
                                                    self.r_diag, 
                                                    normFactor, 
                                                    steer_vector(f[i]), 
                                                    (eva[na:na+1], eve[:, na:na+1]))[0]
                  if self.r_diag:  # set (unphysical) negative output values to 0
                      indNegSign = sign(beamformerOutput) < 0
                      beamformerOutput[indNegSign] = 0
                  ac[i] = beamformerOutput
                  fr[i] = 1
  
  class BeamformerMusic( BeamformerEig ):
      """
      Beamforming using the MUSIC algorithm, see :ref:`Schmidt, 1986<Schmidt1986>`.
      """
  
      # Boolean flag, if 'True', the main diagonal is removed before beamforming;
      # for MUSIC beamforming r_diag is set to 'False'.
      r_diag = Enum(False, 
          desc="removal of diagonal")
  
      # assumed number of sources, should be set to a value not too small
      # defaults to 1
      n = Int(1, 
          desc="assumed number of sources")
  
      def calc(self, ac, fr):
          """
          Calculates the MUSIC result for the frequencies defined by :attr:`freq_data`
          
          This is an internal helper function that is automatically called when 
          accessing the beamformer's :attr:`~BeamformerBase.result` or calling
          its :meth:`~BeamformerBase.synthetic` method.        
          
          Parameters
          ----------
          ac : array of floats
              This array of dimension ([number of frequencies]x[number of gridpoints])
              is used as call-by-reference parameter and contains the calculated
              value after calling this method. 
          fr : array of booleans
              The entries of this [number of frequencies]-sized array are either 
              'True' (if the result for this frequency has already been calculated)
              or 'False' (for the frequencies where the result has yet to be calculated).
              After the calculation at a certain frequency the value will be set
              to 'True'
          
          Returns
          -------
          This method only returns values through the *ac* and *fr* parameters
          """
          f = self.freq_data.fftfreq()
          nMics = self.freq_data.numchannels
          n = int(self.steer.mics.num_mics-self.na)
          normFactor = self.sig_loss_norm() * nMics**2
          param_steer_type, steer_vector = self._beamformer_params()
          for i in self.freq_data.indices:
              if not fr[i]:
                  eva = array(self.freq_data.eva[i], dtype='float64')
                  eve = array(self.freq_data.eve[i], dtype='complex128')
                  beamformerOutput = beamformerFreq(param_steer_type, 
                                                    self.r_diag, 
                                                    normFactor, 
                                                    steer_vector(f[i]), 
                                                    (eva[:n], eve[:, :n]))[0]
                  ac[i] = 4e-10*beamformerOutput.min() / beamformerOutput
                  fr[i] = 1
  
  class PointSpreadFunction (HasPrivateTraits):
      """
      The point spread function.
      
      This class provides tools to calculate the PSF depending on the used 
      microphone geometry, focus grid, flow environment, etc.
      The PSF is needed by several deconvolution algorithms to correct
      the aberrations when using simple delay-and-sum beamforming.
      """
      
          # Instance of :class:`~acoular.fbeamform.SteeringVector` or its derived classes
      # that contains information about the steering vector. This is a private trait.
      # Do not set this directly, use `steer` trait instead.
      _steer_obj = Instance(SteeringVector(), SteeringVector)   
      
      #: :class:`~acoular.fbeamform.SteeringVector` or derived object. 
      #: Defaults to :class:`~acoular.fbeamform.SteeringVector` object.
      steer = Property(desc="steering vector object")  
      
      def _get_steer(self):
          return self._steer_obj
      
      def _set_steer(self, steer):
          if isinstance(steer, SteeringVector):
              self._steer_obj = steer
          elif steer in ('true level', 'true location', 'classic', 'inverse'):
              # Type of steering vectors, see also :ref:`Sarradj, 2012<Sarradj2012>`.
              warn("Deprecated use of 'steer' trait. "
                   "Please use object of class 'SteeringVector' in the future.", 
                   Warning, stacklevel = 2)
              self._steer_obj = SteeringVector(steer_type = steer)
          else:
              raise(TraitError(args=self,
                               name='steer', 
                               info='SteeringVector',
                               value=steer))
  
      # --- List of backwards compatibility traits and their setters/getters -----------
      
      # :class:`~acoular.environments.Environment` or derived object. 
      # Deprecated! Only kept for backwards compatibility. 
      # Now governed by :attr:`steer` trait.
      env = Property()
      
      def _get_env(self):
          return self._steer_obj.env    
      
      def _set_env(self, env):
          warn("Deprecated use of 'env' trait. ", Warning, stacklevel = 2)
          self._steer_obj.env = env
      
      # The speed of sound.
      # Deprecated! Only kept for backwards compatibility. 
      # Now governed by :attr:`steer` trait.
      c = Property()
      
      def _get_c(self):
          return self._steer_obj.env.c
      
      def _set_c(self, c):
          warn("Deprecated use of 'c' trait. ", Warning, stacklevel = 2)
          self._steer_obj.env.c = c
     
      # :class:`~acoular.grids.Grid`-derived object that provides the grid locations.
      # Deprecated! Only kept for backwards compatibility. 
      # Now governed by :attr:`steer` trait.
      grid = Property()
  
      def _get_grid(self):
          return self._steer_obj.grid
      
      def _set_grid(self, grid):
          warn("Deprecated use of 'grid' trait. ", Warning, stacklevel = 2)
          self._steer_obj.grid = grid
      
      # :class:`~acoular.microphones.MicGeom` object that provides the microphone locations.
      # Deprecated! Only kept for backwards compatibility. 
      # Now governed by :attr:`steer` trait
      mpos = Property()
      
      def _get_mpos(self):
          return self._steer_obj.mics
      
      def _set_mpos(self, mpos):
          warn("Deprecated use of 'mpos' trait. ", Warning, stacklevel = 2)
          self._steer_obj.mics = mpos
      
      
      # Sound travel distances from microphone array center to grid points (r0)
      # and all array mics to grid points (rm). Readonly.
      # Deprecated! Only kept for backwards compatibility. 
      # Now governed by :attr:`steer` trait
      r0 = Property()
      def _get_r0(self):
          return self._steer_obj.r0
      
      rm = Property()
      def _get_rm(self):
          return self._steer_obj.rm
      
      # --- End of backwards compatibility traits --------------------------------------
      
      
      #: Indices of grid points to calculate the PSF for.
      grid_indices = CArray( dtype=int, value=array([]), 
                       desc="indices of grid points for psf") #value=array([]), value=self.grid.pos(),
      
      #: Flag that defines how to calculate and store the point spread function
      #: defaults to 'single'.
      #:
      #: * 'full': Calculate the full PSF (for all grid points) in one go (should be used if the PSF at all grid points is needed, as with :class:`DAMAS<BeamformerDamas>`)
      #: * 'single': Calculate the PSF for the grid points defined by :attr:`grid_indices`, one by one (useful if not all PSFs are needed, as with :class:`CLEAN<BeamformerClean>`)
      #: * 'block': Calculate the PSF for the grid points defined by :attr:`grid_indices`, in one go (useful if not all PSFs are needed, as with :class:`CLEAN<BeamformerClean>`)
      #: * 'readonly': Do not attempt to calculate the PSF since it should already be cached (useful if multiple processes have to access the cache file)
      calcmode = Trait('single', 'block', 'full', 'readonly',
                       desc="mode of calculation / storage")
                
      #: Floating point precision of property psf. Corresponding to numpy dtypes. Default = 64 Bit.
      precision = Trait('float64', 'float32',
              desc="precision (32/64 Bit) of result, corresponding to numpy dtypes")
      
      #: The actual point spread function.
      psf = Property(desc="point spread function")
      
      #: Frequency to evaluate the PSF for; defaults to 1.0. 
      freq = Float(1.0, desc="frequency")
  
      # hdf5 cache file
      h5f = Instance( H5CacheFileBase, transient = True )
      
      # internal identifier
      digest = Property( depends_on = ['_steer_obj.digest', 'precision'], cached = True)
  
      @cached_property
      def _get_digest( self ):
          return digest( self )
      
      def _get_filecache( self ):
          """
          function collects cached results from file depending on 
          global/local caching behaviour. Returns (None, None) if no cachefile/data 
          exist and global caching mode is 'readonly'.
          """
          filename = 'psf' + self.digest
          nodename = ('Hz_%.2f' % self.freq).replace('.', '_')
  #        print("get cachefile:", filename)
          H5cache.get_cache_file( self, filename ) 
          if not self.h5f: # only happens in case of global caching readonly
  #            print("no cachefile:", filename)
              return (None, None)# only happens in case of global caching readonly
                      
          if config.global_caching == 'overwrite' and self.h5f.is_cached(nodename):
  #            print("remove existing data for nodename",nodename)
              self.h5f.remove_data(nodename) # remove old data before writing in overwrite mode
          
          if not self.h5f.is_cached(nodename):
  #            print("no data existent for nodename:", nodename)
              if config.global_caching == 'readonly':
                  return (None, None)
              else:
  #                print("initialize data.")
                  gs = self.steer.grid.size
                  group = self.h5f.create_new_group(nodename)
                  self.h5f.create_compressible_array('result',
                                        (gs, gs),
                                        self.precision,
                                        group)
                  self.h5f.create_compressible_array('gridpts',
                                        (gs,),
                                        'int8',#'bool', 
                                        group)
          ac = self.h5f.get_data_by_reference('result','/'+nodename)
          gp = self.h5f.get_data_by_reference('gridpts','/'+nodename)
          return (ac,gp)        
  
      def _get_psf ( self ):
          """
          This is the :attr:`psf` getter routine.
          The point spread function is either loaded or calculated.
          """
          gs = self.steer.grid.size
          if not self.grid_indices.size: 
              self.grid_indices = arange(gs)
  
          if not config.global_caching == 'none':
  #            print("get filecache..")
              (ac,gp) = self._get_filecache()
              if ac and gp: 
  #                print("cached data existent")
                  if not gp[:][self.grid_indices].all():
  #                    print("calculate missing results")                            
                      if self.calcmode == 'readonly':
                          raise ValueError('Cannot calculate missing PSF (points) in \'readonly\' mode.')
                      if config.global_caching == 'readonly':
                          (ac, gp) = (ac[:], gp[:])
                          self.calc_psf(ac,gp)
                          return ac[:,self.grid_indices]
                      else:
                          self.calc_psf(ac,gp)
                          self.h5f.flush()
                          return ac[:,self.grid_indices]
  #                else:
  #                    print("cached results are complete! return.")
                  return ac[:,self.grid_indices]
              else: # no cached data/file
  #                print("no caching, calculate result")
                  ac = zeros((gs, gs), dtype=self.precision)
                  gp = zeros((gs,), dtype='int8')
                  self.calc_psf(ac,gp)
          else: # no caching activated
  #            print("no caching activated, calculate result")
              ac = zeros((gs, gs), dtype=self.precision)
              gp = zeros((gs,), dtype='int8')
              self.calc_psf(ac,gp)
          return ac[:,self.grid_indices] 
  
      def calc_psf( self, ac, gp ):
          """
          point-spread function calculation
          """
          if self.calcmode != 'full':
              # calc_ind has the form [True, True, False, True], except
              # when it has only 1 entry (value True/1 would be ambiguous)
              if self.grid_indices.size == 1:
                  calc_ind = [0]
              else:
                  calc_ind = invert(gp[:][self.grid_indices])
          
          # get indices which have the value True = not yet calculated
              g_ind_calc = self.grid_indices[calc_ind]
          
          if self.calcmode == 'single': # calculate selected psfs one-by-one
              for ind in g_ind_calc:
                  ac[:,ind] = self._psfCall([ind])[:,0]
                  gp[ind] = 1
          elif self.calcmode == 'full': # calculate all psfs in one go
              gp[:] = 1
              ac[:] = self._psfCall(arange(self.steer.grid.size))
          else: # 'block' # calculate selected psfs in one go
              hh = self._psfCall(g_ind_calc)
              indh = 0
              for ind in g_ind_calc:
                  gp[ind] = 1
                  ac[:,ind] = hh[:,indh]
                  indh += 1
  
      def _psfCall(self, ind):
          """
          Manages the calling of the core psf functionality.
          
          Parameters
          ----------
          ind : list of int
              Indices of gridpoints which are assumed to be sources.
              Normalization factor for the beamforming result (e.g. removal of diag is compensated with this.)
  
          Returns
          -------
          The psf [1, nGridPoints, len(ind)]
          """
          if type(self.steer) == SteeringVector: # for simple steering vector, use faster method
              result = calcPointSpreadFunction(self.steer.steer_type, 
                                               self.steer.r0, 
                                               self.steer.rm, 
                                               2*pi*self.freq/self.env.c, 
                                               ind, self.precision)
          else: # for arbitrary steering sectors, use general calculation
              # there is a version of this in fastFuncs, may be used later after runtime testing and debugging
              product = dot(self.steer.steer_vector(self.freq).conj(), self.steer.transfer(self.freq,ind).T)
              result = (product * product.conj()).real
          return result
  
  class BeamformerDamas (BeamformerBase):
      """
      DAMAS deconvolution, see :ref:`Brooks and Humphreys, 2006<BrooksHumphreys2006>`.
      Needs a-priori delay-and-sum beamforming (:class:`BeamformerBase`).
      """
  
      #: :class:`BeamformerBase` object that provides data for deconvolution.
      beamformer = Trait(BeamformerBase)
  
      #: :class:`~acoular.spectra.PowerSpectra` object that provides the cross spectral matrix; 
      #: is set automatically.
      freq_data = Delegate('beamformer')
  
      #: Boolean flag, if 'True' (default), the main diagonal is removed before beamforming; 
      #: is set automatically.
      r_diag =  Delegate('beamformer')
      
      #: instance of :class:`~acoular.fbeamform.SteeringVector` or its derived classes,
      #: that contains information about the steering vector. Is set automatically.
      steer = Delegate('beamformer')
      
      #: Floating point precision of result, is set automatically.
      precision = Delegate('beamformer')
      
      #: The floating-number-precision of the PSFs. Default is 64 bit.
      psf_precision = Trait('float64', 'float32', 
                            desc="precision of PSF")
  
      #: Number of iterations, defaults to 100.
      n_iter = Int(100, 
          desc="number of iterations")
      
      #: Damping factor in modified gauss-seidel
      damp = Float(1.0,
                            desc="damping factor in modified gauss-seidel-DAMAS-approach")
  
      #: Flag that defines how to calculate and store the point spread function, 
      #: defaults to 'full'. See :attr:`PointSpreadFunction.calcmode` for details.
      calcmode = Trait('full', 'single', 'block', 'readonly',
                       desc="mode of psf calculation / storage")
      
      # internal identifier
      digest = Property( 
          depends_on = ['beamformer.digest', 'n_iter', 'damp', 'psf_precision'], 
          )
  
      # internal identifier
      ext_digest = Property( 
          depends_on = ['digest', 'beamformer.ext_digest'], 
          )
      
      @cached_property
      def _get_digest( self ):
          return digest( self )
        
      @cached_property
      def _get_ext_digest( self ):
          return digest( self, 'ext_digest' )
      
      def calc(self, ac, fr):
          """
          Calculates the DAMAS result for the frequencies defined by :attr:`freq_data`
          
          This is an internal helper function that is automatically called when 
          accessing the beamformer's :attr:`~BeamformerBase.result` or calling
          its :meth:`~BeamformerBase.synthetic` method.        
          A Gauss-Seidel algorithm implemented in C is used for computing the result.
          
          Parameters
          ----------
          ac : array of floats
              This array of dimension ([number of frequencies]x[number of gridpoints])
              is used as call-by-reference parameter and contains the calculated
              value after calling this method. 
          fr : array of booleans
              The entries of this [number of frequencies]-sized array are either 
              'True' (if the result for this frequency has already been calculated)
              or 'False' (for the frequencies where the result has yet to be calculated).
              After the calculation at a certain frequency the value will be set
              to 'True'
          
          Returns
          -------
          This method only returns values through the *ac* and *fr* parameters
          """
          f = self.freq_data.fftfreq()
          p = PointSpreadFunction(steer=self.steer, calcmode=self.calcmode, precision=self.psf_precision)
          for i in self.freq_data.indices:
              if not fr[i]:
                  y = array(self.beamformer.result[i])
                  x = y.copy()
                  p.freq = f[i]
                  psf = p.psf[:]
                  damasSolverGaussSeidel(psf, y, self.n_iter, self.damp, x)
                  ac[i] = x
                  fr[i] = 1
  
  class BeamformerDamasPlus (BeamformerDamas):
      """
      DAMAS deconvolution, see :ref:`Brooks and Humphreys, 2006<BrooksHumphreys2006>`,
      for solving the system of equations, instead of the original Gauss-Seidel 
      iterations, this class employs the NNLS or linear programming solvers from 
      scipy.optimize or one  of several optimization algorithms from the scikit-learn module.
      Needs a-priori delay-and-sum beamforming (:class:`BeamformerBase`).
      """
      
      #: Type of fit method to be used ('LassoLars', 
      #: 'OMPCV', 'LP', or 'NNLS', defaults to 'NNLS').
      #: These methods are implemented in 
      #: the `scikit-learn <http://scikit-learn.org/stable/user_guide.html>`_ 
      #: module or within scipy.optimize respectively.
      method = Trait('NNLS','LP','LassoLars', 'OMPCV',  
                     desc="method used for solving deconvolution problem")
      
      #: Weight factor for LassoLars method,
      #: defaults to 0.0.
      # (Values in the order of 10^⁻9 should produce good results.)
      alpha = Range(0.0, 1.0, 0.0,
                    desc="Lasso weight factor")
      
      #: Maximum number of iterations,
      #: tradeoff between speed and precision;
      #: defaults to 500
      max_iter = Int(500,
                     desc="maximum number of iterations")
      
      #: Unit multiplier for evaluating, e.g., nPa instead of Pa. 
      #: Values are converted back before returning. 
      #: Temporary conversion may be necessary to not reach machine epsilon
      #: within fitting method algorithms. Defaults to 1e9.
      unit_mult = Float(1e9,
                        desc = "unit multiplier")
      
      # internal identifier
      digest = Property( 
          depends_on = ['beamformer.digest','alpha', 'method', 
                        'max_iter', 'unit_mult'], 
          )
  
      # internal identifier
      ext_digest = Property( 
          depends_on = ['digest', 'beamformer.ext_digest'], 
          )
      
      @cached_property
      def _get_digest( self ):
          return digest( self )
        
      @cached_property
      def _get_ext_digest( self ):
          return digest( self, 'ext_digest' )
      
      def calc(self, ac, fr):
          """
          Calculates the DAMAS result for the frequencies defined by :attr:`freq_data`
          
          This is an internal helper function that is automatically called when 
          accessing the beamformer's :attr:`~BeamformerBase.result` or calling
          its :meth:`~BeamformerBase.synthetic` method.        
          
          Parameters
          ----------
          ac : array of floats
              This array of dimension ([number of frequencies]x[number of gridpoints])
              is used as call-by-reference parameter and contains the calculated
              value after calling this method. 
          fr : array of booleans
              The entries of this [number of frequencies]-sized array are either 
              'True' (if the result for this frequency has already been calculated)
              or 'False' (for the frequencies where the result has yet to be calculated).
              After the calculation at a certain frequency the value will be set
              to 'True'
          
          Returns
          -------
          This method only returns values through the *ac* and *fr* parameters
          """
          f = self.freq_data.fftfreq()
          p = PointSpreadFunction(steer=self.steer, calcmode=self.calcmode, precision=self.psf_precision)
          unit = self.unit_mult
          for i in self.freq_data.indices:
              if not fr[i]:
                  y = self.beamformer.result[i] * unit
                  p.freq = f[i]
                  psf = p.psf[:]
  
                  if self.method == 'NNLS':
                      resopt = nnls(psf,y)[0]
                  elif self.method == 'LP': # linear programming (Dougherty)
                      if self.r_diag:
                          warn('Linear programming solver may fail when CSM main '
                                'diagonal is removed for delay-and-sum beamforming.', 
                                Warning, stacklevel = 5)
                      cT = -1*psf.sum(1) # turn the minimization into a maximization
                      resopt = linprog(c=cT, A_ub=psf, b_ub=y).x # defaults to simplex method and non-negative x
                  elif self.method == 'LassoLars':
                      model = LassoLars(alpha = self.alpha * unit, 
                                        max_iter = self.max_iter)
                  else: # self.method == 'OMPCV':
                      model = OrthogonalMatchingPursuitCV()
                  
                  
                  if self.method in ('NNLS','LP'):
                      ac[i] = resopt / unit
                  else: # sklearn models
                      model.fit(psf,y)
                      ac[i] = model.coef_[:] / unit
                  
                  fr[i] = 1
  
  class BeamformerOrth (BeamformerBase):
      """
      Orthogonal beamforming, see :ref:`Sarradj, 2010<Sarradj2010>`.
      Needs a-priori beamforming with eigenvalue decomposition (:class:`BeamformerEig`).
      """
  
      #: :class:`BeamformerEig` object that provides data for deconvolution.
      beamformer = Trait(BeamformerEig)
  
      #: :class:`~acoular.spectra.PowerSpectra` object that provides the cross spectral matrix 
      #: and eigenvalues, is set automatically.    
      freq_data = Delegate('beamformer')
  
      #: Flag, if 'True' (default), the main diagonal is removed before beamforming, 
      #: is set automatically.
      r_diag =  Delegate('beamformer')
      
      #: instance of :class:`~acoular.fbeamform.SteeringVector` or its derived classes,
      #: that contains information about the steering vector. Is set automatically.
      steer = Delegate('beamformer')
  
      #: List of components to consider, use this to directly set the eigenvalues
      #: used in the beamformer. Alternatively, set :attr:`n`.
      eva_list = CArray(dtype=int,
          desc="components")
          
      #: Number of components to consider, defaults to 1. If set, 
      #: :attr:`eva_list` will contain
      #: the indices of the n largest eigenvalues. Setting :attr:`eva_list` 
      #: afterwards will override this value.
      n = Int(1)
  
      # internal identifier
      digest = Property( 
          depends_on = ['beamformer.digest', 'eva_list'], 
          )
  
      # internal identifier
      ext_digest = Property( 
          depends_on = ['digest', 'beamformer.ext_digest'], 
          )
      
      @cached_property
      def _get_digest( self ):
          return digest( self )
  
      @cached_property
      def _get_ext_digest( self ):
          return digest( self, 'ext_digest' )
      
      @on_trait_change('n')
      def set_eva_list(self):
          """ sets the list of eigenvalues to consider """
          self.eva_list = arange(-1, -1-self.n, -1)
  
      def calc(self, ac, fr):
          """
          Calculates the Orthogonal Beamforming result for the frequencies 
          defined by :attr:`freq_data`.
          
          This is an internal helper function that is automatically called when 
          accessing the beamformer's :attr:`~BeamformerBase.result` or calling
          its :meth:`~BeamformerBase.synthetic` method.        
          
          Parameters
          ----------
          ac : array of floats
              This array of dimension ([number of frequencies]x[number of gridpoints])
              is used as call-by-reference parameter and contains the calculated
              value after calling this method. 
          fr : array of booleans
              The entries of this [number of frequencies]-sized array are either 
              'True' (if the result for this frequency has already been calculated)
              or 'False' (for the frequencies where the result has yet to be calculated).
              After the calculation at a certain frequency the value will be set
              to 'True'
          
          Returns
          -------
          This method only returns values through the *ac* and *fr* parameters
          """
          # prepare calculation
          ii = []
          for i in self.freq_data.indices:        
              if not fr[i]:
                  ii.append(i)
          numchannels = self.freq_data.numchannels
          e = self.beamformer
          for n in self.eva_list:
              e.n = n
              for i in ii:
                  ac[i, e.result[i].argmax()]+=e.freq_data.eva[i, n]/numchannels
          for i in ii:
              fr[i] = 1
      
  class BeamformerCleansc( BeamformerBase ):
      """
      CLEAN-SC deconvolution, see :ref:`Sijtsma, 2007<Sijtsma2007>`.
      Classic delay-and-sum beamforming is already included.
      """
  
      #: no of CLEAN-SC iterations
      #: defaults to 0, i.e. automatic (max 2*numchannels)
      n = Int(0, 
          desc="no of iterations")
  
      #: iteration damping factor
      #: defaults to 0.6
      damp = Range(0.01, 1.0, 0.6, 
          desc="damping factor")
  
      #: iteration stop criterion for automatic detection
      #: iteration stops if power[i]>power[i-stopn]
      #: defaults to 3
      stopn = Int(3, 
          desc="stop criterion index")
  
      # internal identifier
      digest = Property( 
          depends_on = ['freq_data.digest', '_steer_obj.digest', 'r_diag', 'n', 'damp', 'stopn'])
  
      @cached_property
      def _get_digest( self ):
          return digest( self )
  
      def calc(self, ac, fr):
          """
          Calculates the CLEAN-SC result for the frequencies defined by :attr:`freq_data`
          
          This is an internal helper function that is automatically called when 
          accessing the beamformer's :attr:`~BeamformerBase.result` or calling
          its :meth:`~BeamformerBase.synthetic` method.        
          
          Parameters
          ----------
          ac : array of floats
              This array of dimension ([number of frequencies]x[number of gridpoints])
              is used as call-by-reference parameter and contains the calculated
              value after calling this method. 
          fr : array of booleans
              The entries of this [number of frequencies]-sized array are either 
              'True' (if the result for this frequency has already been calculated)
              or 'False' (for the frequencies where the result has yet to be calculated).
              After the calculation at a certain frequency the value will be set
              to 'True'
          
          Returns
          -------
          This method only returns values through the *ac* and *fr* parameters
          """
  
          # prepare calculation
          normFactor = self.sig_loss_norm()
          numchannels = self.freq_data.numchannels
          f = self.freq_data.fftfreq()
          result = zeros((self.steer.grid.size), 'f') 
          normFac = self.sig_loss_norm()
          if not self.n:
              J = numchannels*2
          else:
              J = self.n
          powers = zeros(J, 'd')
          
          param_steer_type, steer_vector = self._beamformer_params()
          for i in self.freq_data.indices:
              if not fr[i]:
                  csm = array(self.freq_data.csm[i], dtype='complex128', copy=1)
                  #h = self.steer._beamformerCall(f[i], self.r_diag, normFactor, (csm,))[0]
                  h = beamformerFreq(param_steer_type, 
                                     self.r_diag, 
                                     normFactor, 
                                     steer_vector(f[i]), 
                                     csm)[0]
                  # CLEANSC Iteration
                  result *= 0.0
                  for j in range(J):
                      xi_max = h.argmax() #index of maximum
                      powers[j] = hmax = h[xi_max] #maximum
                      result[xi_max] += self.damp * hmax
                      if  j > self.stopn and hmax > powers[j-self.stopn]:
                          break
                      wmax = self.steer.steer_vector(f[i],xi_max) * sqrt(normFac)
                      wmax = wmax[0].conj()  # as old code worked with conjugated csm..should be updated
                      hh = wmax.copy()
                      D1 = dot(csm.T - diag(diag(csm)), wmax)/hmax
                      ww = wmax.conj()*wmax
                      for m in range(20):
                          H = hh.conj()*hh
                          hh = (D1+H*wmax)/sqrt(1+dot(ww, H))
                      hh = hh[:, newaxis]
                      csm1 = hmax*(hh*hh.conj().T)
                      
                      #h1 = self.steer._beamformerCall(f[i], self.r_diag, normFactor, (array((hmax, ))[newaxis, :], hh[newaxis, :].conjugate()))[0]
                      h1 = beamformerFreq(param_steer_type, 
                                          self.r_diag, 
                                          normFactor, 
                                          steer_vector(f[i]), 
                                          (array((hmax, )), hh.conj()))[0]
                      h -= self.damp * h1
                      csm -= self.damp * csm1.T#transpose(0,2,1)
                  ac[i] = result
                  fr[i] = 1
  
  class BeamformerClean (BeamformerBase):
      """
      CLEAN deconvolution, see :ref:`Hoegbom, 1974<Hoegbom1974>`.
      Needs a-priori delay-and-sum beamforming (:class:`BeamformerBase`).
      """
  
      # BeamformerBase object that provides data for deconvolution
      beamformer = Trait(BeamformerBase)
  
      # PowerSpectra object that provides the cross spectral matrix
      freq_data = Delegate('beamformer')
  
      # flag, if true (default), the main diagonal is removed before beamforming
      #r_diag =  Delegate('beamformer')
      
      #: instance of :class:`~acoular.fbeamform.SteeringVector` or its derived classes,
      #: that contains information about the steering vector. Is set automatically.
      steer = Delegate('beamformer')
      
      #: Floating point precision of result, is set automatically.
      precision = Delegate('beamformer')
      
      #: The floating-number-precision of the PSFs. Default is 64 bit.
      psf_precision = Trait('float64', 'float32', 
                       desc="precision of PSF.")
      
      # iteration damping factor
      # defaults to 0.6
      damp = Range(0.01, 1.0, 0.6, 
          desc="damping factor")
          
      # max number of iterations
      n_iter = Int(100, 
          desc="maximum number of iterations")
  
      # how to calculate and store the psf
      calcmode = Trait('block', 'full', 'single', 'readonly',
                       desc="mode of psf calculation / storage")
                       
      # internal identifier
      digest = Property( 
          depends_on = ['beamformer.digest', 'n_iter', 'damp', 'psf_precision'], 
          )
  
      # internal identifier
      ext_digest = Property( 
          depends_on = ['digest', 'beamformer.ext_digest'], 
          )
      
      @cached_property
      def _get_digest( self ):
          return digest( self )
        
      @cached_property
      def _get_ext_digest( self ):
          return digest( self, 'ext_digest' )
      
      def calc(self, ac, fr):
          """
          Calculates the CLEAN result for the frequencies defined by :attr:`freq_data`
          
          This is an internal helper function that is automatically called when 
          accessing the beamformer's :attr:`~BeamformerBase.result` or calling
          its :meth:`~BeamformerBase.synthetic` method.        
          
          Parameters
          ----------
          ac : array of floats
              This array of dimension ([number of frequencies]x[number of gridpoints])
              is used as call-by-reference parameter and contains the calculated
              value after calling this method. 
          fr : array of booleans
              The entries of this [number of frequencies]-sized array are either 
              'True' (if the result for this frequency has already been calculated)
              or 'False' (for the frequencies where the result has yet to be calculated).
              After the calculation at a certain frequency the value will be set
              to 'True'
          
          Returns
          -------
          This method only returns values through the *ac* and *fr* parameters
          """
          f = self.freq_data.fftfreq()
          gs = self.steer.grid.size
          
          if self.calcmode == 'full':
              warn("calcmode = 'full', possibly slow CLEAN performance. "
                   "Better use 'block' or 'single'.", Warning, stacklevel = 2)
          p = PointSpreadFunction(steer=self.steer, calcmode=self.calcmode, precision=self.psf_precision)
          for i in self.freq_data.indices:
              if not fr[i]:
                  p.freq = f[i]
                  dirty = self.beamformer.result[i].copy()
                  clean = zeros(gs, dtype=dirty.dtype)
                  
                  i_iter = 0
                  flag = True
                  while flag:
                      # TODO: negative werte!!!
                      dirty_sum = abs(dirty).sum(0)
                      next_max = dirty.argmax(0)
                      p.grid_indices = array([next_max])
                      psf = p.psf.reshape(gs,)
                      new_amp = self.damp * dirty[next_max] #/ psf[next_max]
                      clean[next_max] += new_amp
                      dirty -= psf * new_amp
                      i_iter += 1
                      flag = (dirty_sum > abs(dirty).sum(0) \
                              and i_iter < self.n_iter \
                              and max(dirty) > 0)
                  
                  ac[i] = clean            
                  fr[i] = 1
  
  class BeamformerCMF ( BeamformerBase ):
      """
      Covariance Matrix Fitting, see :ref:`Yardibi et al., 2008<Yardibi2008>`.
      This is not really a beamformer, but an inverse method.
      """
  
      #: Type of fit method to be used ('LassoLars', 'LassoLarsBIC', 
      #: 'OMPCV' or 'NNLS', defaults to 'LassoLars').
      #: These methods are implemented in 
      #: the `scikit-learn <http://scikit-learn.org/stable/user_guide.html>`_ 
      #: module.
      method = Trait('LassoLars', 'LassoLarsBIC',  \
          'OMPCV', 'NNLS','fmin_l_bfgs_b', desc="fit method used")
          
      #: Weight factor for LassoLars method,
      #: defaults to 0.0.
      #: (Use values in the order of 10^⁻9 for good results.)
      alpha = Range(0.0, 1.0, 0.0, 
          desc="Lasso weight factor")
      
      #: Maximum number of iterations,
      #: tradeoff between speed and precision;
      #: defaults to 500
      max_iter = Int(500, 
          desc="maximum number of iterations")
  
      
      #: Unit multiplier for evaluating, e.g., nPa instead of Pa. 
      #: Values are converted back before returning. 
      #: Temporary conversion may be necessary to not reach machine epsilon
      #: within fitting method algorithms. Defaults to 1e9.
      unit_mult = Float(1e9,
                        desc = "unit multiplier")
  
      # internal identifier
      digest = Property( 
          depends_on = ['freq_data.digest', 'alpha', 'method', 'max_iter', 'unit_mult', 'r_diag', 'steer.inv_digest'], 
          )
  
      @cached_property
      def _get_digest( self ):
          return digest( self )
     
  
      def calc(self, ac, fr):
          """
          Calculates the CMF result for the frequencies defined by :attr:`freq_data`
          
          This is an internal helper function that is automatically called when 
          accessing the beamformer's :attr:`~BeamformerBase.result` or calling
          its :meth:`~BeamformerBase.synthetic` method.        
          
          Parameters
          ----------
          ac : array of floats
              This array of dimension ([number of frequencies]x[number of gridpoints])
              is used as call-by-reference parameter and contains the calculated
              value after calling this method. 
          fr : array of booleans
              The entries of this [number of frequencies]-sized array are either 
              'True' (if the result for this frequency has already been calculated)
              or 'False' (for the frequencies where the result has yet to be calculated).
              After the calculation at a certain frequency the value will be set
              to 'True'
          
          Returns
          -------
          This method only returns values through the *ac* and *fr* parameters
          """
          
          # function to repack complex matrices to deal with them in real number space
          def realify(M):
              return vstack([M.real,M.imag])
  
              
          # prepare calculation
          i = self.freq_data.indices
          f = self.freq_data.fftfreq()
          nc = self.freq_data.numchannels
          numpoints = self.steer.grid.size
          unit = self.unit_mult
  
          for i in self.freq_data.indices:        
              if not fr[i]:
                  csm = array(self.freq_data.csm[i], dtype='complex128',copy=1)
  
                  h = self.steer.transfer(f[i]).T
                  
                  # reduced Kronecker product (only where solution matrix != 0)
                  Bc = ( h[:,:,newaxis] * \
                         h.conjugate().T[newaxis,:,:] )\
                           .transpose(2,0,1)
                  Ac = Bc.reshape(nc*nc,numpoints)
                  
                  # get indices for upper triangular matrices (use tril b/c transposed)
                  ind = reshape(tril(ones((nc,nc))), (nc*nc,)) > 0
                  
                  ind_im0 = (reshape(eye(nc),(nc*nc,)) == 0)[ind]
                  if self.r_diag:
                      # omit main diagonal for noise reduction
                      ind_reim = hstack([ind_im0, ind_im0])
                  else:
                      # take all real parts -- also main diagonal
                      ind_reim = hstack([ones(size(ind_im0),)>0,ind_im0])
                      ind_reim[0]=True # TODO: warum hier extra definiert??
  #                    if sigma2:
  #                        # identity matrix, needed when noise term sigma is used
  #                        I  = eye(nc).reshape(nc*nc,1)                
  #                        A = realify( hstack([Ac, I])[ind,:] )[ind_reim,:]
  #                        # ... ac[i] = model.coef_[:-1]
  #                    else:
  
                  A = realify( Ac [ind,:] )[ind_reim,:]
                  # use csm.T for column stacking reshape!
                  R = realify( reshape(csm.T, (nc*nc,1))[ind,:] )[ind_reim,:] * unit
                  # choose method
                  if self.method == 'LassoLars':
                      model = LassoLars(alpha = self.alpha * unit,
                                        max_iter = self.max_iter)
                  elif self.method == 'LassoLarsBIC':
                      model = LassoLarsIC(criterion = 'bic',
                                          max_iter = self.max_iter)
                  elif self.method == 'OMPCV':
                      model = OrthogonalMatchingPursuitCV()
  
                  # nnls is not in sklearn
                  if self.method == 'NNLS':
                      ac[i] , x = nnls(A,R.flat)
                      ac[i] /= unit
                      
                  elif self.method == 'fmin_l_bfgs_b':
                      #function to minimize
                      def function(x):
                          #function
                          func = x.T@A.T@A@x - 2*R.T@A@x + R.T@R                
                          #derivitaive
                          der = 2*A.T@A@x.T[:, newaxis] - 2*A.T@R 
                          return  func[0].T, der[:,0]
                      
                      # initial guess
                      x0 = ones([numpoints])   
                      #boundarys - set to non negative
                      boundarys = tile((0, +inf), (len(x0),1))
                      
                      #optimize
                      ac[i], yval, dicts =  fmin_l_bfgs_b(function, x0, fprime=None, args=(),  
                                                            approx_grad=0, bounds=boundarys, m=10,
                                                            factr=10000000.0, pgtol=1e-05, epsilon=1e-08,
                                                            iprint=-1, maxfun=15000, maxiter=self.max_iter,
                                                            disp=None, callback=None, maxls=20)
                      
                      ac[i] /= unit
                  else:
                      model.fit(A,R[:,0])
                      ac[i] = model.coef_[:] / unit
                  fr[i] = 1
                  
  
  
  
  class BeamformerSODIX( BeamformerBase ):
      """
      SODIX, see Funke, Ein Mikrofonarray-Verfahren zur Untersuchung der
      Schallabstrahlung von Turbofantriebwerken, 2017. and 
      Oertwig, Advancements in the source localization method SODIX and
      application to short cowl engine data, 2019
      
      Source directivity modeling in the cross-spectral matrix
      """
      #: Type of fit method to be used ('fmin_l_bfgs_b').
      #: These methods is implemented in 
      #: the `scipy module.
      method = Trait('fmin_l_bfgs_b', desc="fit method used")
          
      #: Maximum number of iterations,
      #: tradeoff between speed and precision;
      #: defaults to 200
      max_iter = Int(200, 
          desc="maximum number of iterations")
      
      #: Norm to consider for the regularization 
      #: defaults to L-1 Norm
      pnorm= Float(1,desc="Norm for regularization")
      
      #: Weight factor for regularization,
      #: defaults to 0.0.
      alpha = Range(0.0, 1.0, 0.0, 
          desc="regularization factor")
  
      #: Unit multiplier for evaluating, e.g., nPa instead of Pa. 
      #: Values are converted back before returning. 
      #: Temporary conversion may be necessary to not reach machine epsilon
      #: within fitting method algorithms. Defaults to 1e9.
      unit_mult = Float(1e9,
                        desc = "unit multiplier")
      
      #: The beamforming result as squared sound pressure values 
      #: at all grid point locations (readonly).
      #: Returns a (number of frequencies, number of gridpoints) array of floats.
      sodix_result = Property(
          desc="beamforming result")
  
      # internal identifier
      digest = Property( 
          depends_on = ['freq_data.digest', 'alpha', 'method', 'max_iter', 'unit_mult', 'r_diag', 'steer.inv_digest'], 
          )
  
      @cached_property
      def _get_digest( self ):
          return digest( self )
      
      def _get_filecache( self ):
              """
              function collects cached results from file depending on 
              global/local caching behaviour. Returns (None, None) if no cachefile/data 
              exist and global caching mode is 'readonly'.
              """
              H5cache.get_cache_file( self, self.freq_data.basename ) 
              if not self.h5f: 
                  return (None, None)# only happens in case of global caching readonly
      
              nodename = self.__class__.__name__ + self.digest
              if config.global_caching == 'overwrite' and self.h5f.is_cached(nodename):
                  self.h5f.remove_data(nodename) # remove old data before writing in overwrite mode
              
              if not self.h5f.is_cached(nodename):
                  if config.global_caching == 'readonly': 
                      return (None, None)
                  else:
      #                print("initialize data.")
                      numfreq = self.freq_data.fftfreq().shape[0]# block_size/2 + 1steer_obj
                      group = self.h5f.create_new_group(nodename)
                      self.h5f.create_compressible_array('result',
                                            (numfreq, self.steer.grid.size*self.steer.mics.num_mics),
                                            self.precision,
                                            group)
                      self.h5f.create_compressible_array('freqs',
                                            (numfreq, ),
                                            'int8',#'bool', 
                                            group)
              ac = self.h5f.get_data_by_reference('result','/'+nodename)
              fr = self.h5f.get_data_by_reference('freqs','/'+nodename)
              return (ac,fr)    
      
      @property_depends_on('ext_digest')
      def _get_sodix_result ( self ):
          """
          This is the :attr:`result` getter routine.
          The sodix beamforming result is either loaded or calculated.
          """
          f = self.freq_data
          numfreq = f.fftfreq().shape[0]# block_size/2 + 1steer_obj
          _digest = ''
          while self.digest != _digest:
              _digest = self.digest
              self._assert_equal_channels()
              if not ( # if result caching is active
                      config.global_caching == 'none' or 
                      (config.global_caching == 'individual' and self.cached == False)
                  ):
                  (ac,fr) = self._get_filecache() 
                  if ac and fr: 
                      if not fr[f.ind_low:f.ind_high].all():                       
                          if config.global_caching == 'readonly': 
                              (ac, fr) = (ac[:], fr[:])
                          self.calc(ac,fr)
                          self.h5f.flush()
  
                  else:
                      ac = zeros((numfreq, self.steer.grid.size*self.steer.mics.num_mics), dtype=self.precision)
                      fr = zeros(numfreq, dtype='int8')
                      self.calc(ac,fr)
              else:
                  ac = zeros((numfreq, self.steer.grid.size*self.steer.mics.num_mics), dtype=self.precision)
                  fr = zeros(numfreq, dtype='int8')
                  self.calc(ac,fr)
          return ac
      
      def synthetic( self, f, num=0):
          """
          Evaluates the beamforming result for an arbitrary frequency band.
          
          Parameters
          ----------
          f: float
              Band center frequency. 
          num : integer
              Controls the width of the frequency bands considered; defaults to
              0 (single frequency line).
              
              ===  =====================
              num  frequency band width
              ===  =====================
              0    single frequency line
              1    octave band
              3    third-octave band
              n    1/n-octave band
              ===  =====================
                
          Returns
          -------
          array of floats
              The synthesized frequency band values of the beamforming result at 
              each grid point and each microphone .
              Note that the frequency resolution and therefore the bandwidth 
              represented by a single frequency line depends on 
              the :attr:`sampling frequency<acoular.sources.SamplesGenerator.sample_freq>` and conjugate
              used :attr:`FFT block size<acoular.spectra.PowerSpectra.block_size>`.
          """
          res = self.sodix_result # trigger calculation
          freq = self.freq_data.fftfreq()
          if len(freq) == 0:
              return None
          
          indices = self.freq_data.indices
          
          if num == 0:
              # single frequency line
              ind = searchsorted(freq, f)
              if ind >= len(freq):
                  warn('Queried frequency (%g Hz) not in resolved '
                                'frequency range. Returning zeros.' % f, 
                                Warning, stacklevel = 2)
                  h = zeros_like(res[0])
              else:
                  if freq[ind] != f:
                      warn('Queried frequency (%g Hz) not in set of '
                           'discrete FFT sample frequencies. '
                           'Using frequency %g Hz instead.' % (f,freq[ind]), 
                           Warning, stacklevel = 2)
                  if not (ind in indices):
                      warn('Beamforming result may not have been calculated '
                           'for queried frequency. Check '
                           'freq_data.ind_low and freq_data.ind_high!',
                            Warning, stacklevel = 2)
                  h = res[ind]
          else:
              # fractional octave band
              f1 = f*2.**(-0.5/num)
              f2 = f*2.**(+0.5/num)
              ind1 = searchsorted(freq, f1)
              ind2 = searchsorted(freq, f2)
              if ind1 == ind2:
                  warn('Queried frequency band (%g to %g Hz) does not '
                       'include any discrete FFT sample frequencies. '
                       'Returning zeros.' % (f1,f2), 
                       Warning, stacklevel = 2)
                  h = zeros_like(res[0])
              else:
                  h = sum(res[ind1:ind2], 0)
                  if not ((ind1 in indices) and (ind2 in indices)):
                      warn('Beamforming result may not have been calculated '
                           'for all queried frequencies. Check '
                           'freq_data.ind_low and freq_data.ind_high!',
                            Warning, stacklevel = 2)
          return h.reshape([self.steer.grid.size,self.steer.mics.num_mics])
     
  
      def calc(self, ac, fr):
          """
          Calculates the SODIX result for the frequencies defined by :attr:`freq_data`
          
          This is an internal helper function that is automatically called when 
          accessing the beamformer's :attr:`~Beamformer.sodix_result` or calling
          its :meth:`~BeamformerSODIX.synthetic` method.        
          
          Parameters
          ----------
          ac : array of floats
              This array of dimension ([number of frequencies]x[number of gridpoints]x[number of microphones])
              is used as call-by-reference parameter and contains the calculated
            #/= unit   value after calling this method. 
          fr : array of booleans
              The entries of this [number of frequencies]-sized array are either 
              'True' (if the result for this frequency has already been calculated)
              or 'False' (for the frequencies where the result has yet to be calculated).
              After the calculation at a certain frequency the value will be set
              to 'True'
          
          Returns
          -------
          This method only returns values through the *ac* and *fr* parameters
          """
          
          # prepare calculation
          i = self.freq_data.indices
          f = self.freq_data.fftfreq()
          numpoints = self.steer.grid.size
          #unit = self.unit_mult
          num_mics = self.steer.mics.num_mics
  
          for i in self.freq_data.indices:        
              if not fr[i]:
                  
                  #measured csm
                  csm = array(self.freq_data.csm[i], dtype='complex128',copy=1) 
                  #transfer function
                  h = self.steer.transfer(f[i]).T           
                     
                  if self.method == 'fmin_l_bfgs_b':
                      #function to minimize
                      def function(D): 
                          '''
                          Parameters
                          ----------
                          D 
                          [numpoints*num_mics]
                          
                          Returns
                          -------
                          func - Sodix function to optimize
                               [1]
                          derdrl - derivitaives in direction of D
                              [num_mics*numpoints].
  
                          '''           
                          #### the sodix function ####
                          Djm = D.reshape([numpoints,num_mics])                           
                          csmmod = einsum('jm,jm,jn,jn->mn',h.T,Djm,Djm,h.T.conj() )        
                          func = sum(absolute((csm - csmmod)))**2 + self.alpha*norm(Djm,self.pnorm)
                          ####the sodix  derivitaive ####
                          inner = csm - einsum('jl,jl,jm,jm->lm',h.T,Djm,Djm,h.T.conj() )      
                          derdrl = -4 *  Djm * real(einsum('rm,rl,lm->rl',h.T,h.T.conj(),inner)) + \
                              self.alpha * (abs(Djm)/norm(Djm,self.pnorm))**(1-self.pnorm)*sign(Djm)
  
                          return  func, derdrl[:].flatten()  #func[0]
                      
                      ##### initial guess #### 
                      if all(ac[(i-1)]==0):
                           D0 = ones([numpoints,num_mics])
                      else:
                           D0 = ac[(i-1)]
                      
                      #boundarys - set to non negative [2*(numpoints*num_mics)]
                      boundarys = tile((0, +inf), (numpoints*num_mics,1))
  
                      #optimize with gradient solver
                      # see https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fmin_l_bfgs_b.html                    
  
                      qi = ones([numpoints,num_mics])
                      qi, yval, dicts =  fmin_l_bfgs_b(function, D0, fprime=None, args=(),  #None  
                                                           approx_grad=0, bounds=boundarys, #approx_grad 0 or True
                                                           factr=100.0, pgtol=1e-09, epsilon=1e-08,
                                                            iprint=0, maxfun=1500000, maxiter=self.max_iter,
                                                            disp=None, callback=None, maxls=20)
                      #squared pressure
                      ac[i]=qi**2
                  else:
                      pass
                  fr[i] = 1
  
                  
                  
                  
                  
  class BeamformerGIB(BeamformerEig):  #BeamformerEig #BeamformerBase
      """
      Beamforming GIB methods with different normalizations,
      """
      
      #: Unit multiplier for evaluating, e.g., nPa instead of Pa. 
      #: Values are converted back before returning. 
      #: Temporary conversion may be necessary to not reach machine epsilon
      #: within fitting method algorithms. Defaults to 1e9.
      unit_mult = Float(1e9,
                        desc = "unit multiplier")
  
      #: Maximum number of iterations,
      #: tradeoff between speed and precision;
      #: defaults to 10
      max_iter = Int(10, 
                     desc="maximum number of iterations")
  
      #: Type of fit method to be used ('Suzuki', 'LassoLars', 'LassoLarsCV', 'LassoLarsBIC', 
      #: 'OMPCV' or 'NNLS', defaults to 'Suzuki').
      #: These methods are implemented in 
      #: the `scikit-learn <http://scikit-learn.org/stable/user_guide.html>`_ 
      #: module.
      method = Trait('Suzuki', 'InverseIRLS', 'LassoLars', 'LassoLarsBIC','LassoLarsCV',  \
          'OMPCV', 'NNLS', desc="fit method used")
  
      #: Weight factor for LassoLars method,
      #: defaults to 0.0.
      alpha = Range(0.0, 1.0, 0.0, 
          desc="Lasso weight factor")
      # (use values in the order of 10^⁻9 for good results) 
      
      #: Norm to consider for the regularization in InverseIRLS and Suzuki methods 
      #: defaults to L-1 Norm
      pnorm= Float(1,desc="Norm for regularization")
  
      #: Beta - Fraction of sources maintained after each iteration
      #: defaults to 0.9 
      beta =  Float(0.9,desc="fraction of sources maintained")
      
      #: eps - Regularization parameter for Suzuki algorithm
      #: defaults to 0.05. 
      eps_perc =  Float(0.05,desc="regularization parameter")
  
      # This feature is not fully supported may be changed in the next release 
      # First eigenvalue to consider. Defaults to 0.
      m = Int(0,
                        desc = "First eigenvalue to consider")
      
      
      # internal identifier++++++++++++++++++++++++++++++++++++++++++++++++++
      digest = Property( 
          depends_on = ['steer.inv_digest', 'freq_data.digest', \
              'alpha', 'method', 'max_iter', 'unit_mult', 'eps_perc',\
              'pnorm', 'beta','n', 'm'], 
          )
  
      @cached_property
      def _get_digest( self ):
          return digest( self )
      
      @property_depends_on('n')
      def _get_na( self ):
          na = self.n
          nm = self.steer.mics.num_mics
          if na < 0:
              na = max(nm + na, 0)
          return min(nm - 1, na)
  
      def calc(self, ac, fr):
          
          """
          Calculates the result for the frequencies defined by :attr:`freq_data`
          
          This is an internal helper function that is automatically called when 
          accessing the beamformer's :attr:`~BeamformerBase.result` or calling
          its :meth:`~BeamformerBase.synthetic` method.        
          
          Parameters
          ----------
          ac : array of floats
              This array of dimension ([number of frequencies]x[number of gridpoints])
              is used as call-by-reference parameter and contains the calculated
              value after calling this method. 
          fr : array of booleans
              The entries of this [number of frequencies]-sized array are either 
              'True' (if the result for this frequency has already been calculated)
              or 'False' (for the frequencies where the result has yet to be calculated).
              After the calculation at a certain frequency the value will be set
              to 'True'
          
          Returns
          -------
          This method only returns values through the *ac* and *fr* parameters
  
          """        
          # prepare calculation
          f = self.freq_data.fftfreq()
          n = int(self.na)   #number of eigenvalues
          m = int(self.m)    #number of first eigenvalue
          numchannels = self.freq_data.numchannels   #number of channels
          numpoints = self.steer.grid.size
          hh = zeros((1, numpoints, numchannels), dtype='D')
          
          #Generate a cross spectral matrix, and perform the eigenvalue decomposition
          for i in self.freq_data.indices:        
              if not fr[i]:
                  #for monopole and source strenght Q needs to define density
                  #calculate a transfer matrix A 
                  hh = self.steer.transfer(f[i])
                  A=hh.T                 
                  #eigenvalues and vectors               
                  csm = array(self.freq_data.csm[i], dtype='complex128',copy=1)
                  eva,eve=eigh(csm)
                  eva = eva[::-1]
                  eve = eve[:, ::-1] 
                  eva[eva < max(eva)/1e12] = 0 #set small values zo 0, lowers numerical errors in simulated data
                  #init sources    
                  qi=zeros([n+m,numpoints], dtype='complex128')
                  #Select the number of coherent modes to be processed referring to the eigenvalue distribution.
                  #for s in arange(n):  
                  for s in list(range(m,n+m)):
                      if eva[s] > 0:                    
                          #Generate the corresponding eigenmodes
                          emode=array(sqrt(eva[s])*eve[:,s], dtype='complex128')
                          # choose method for computation
                          if self.method == 'Suzuki':
                              leftpoints=numpoints
                              locpoints=arange(numpoints)         
                              weights=diag(ones(numpoints))             
                              epsilon=arange(self.max_iter)              
                              for it in arange(self.max_iter): 
                                  if numchannels<=leftpoints:
                                      AWA= dot(dot(A[:,locpoints],weights),A[:,locpoints].conj().T)
                                      epsilon[it] = max(absolute(eigvals(AWA)))*self.eps_perc
                                      qi[s,locpoints]=dot(dot(dot(weights,A[:,locpoints].conj().T),inv(AWA+eye(numchannels)*epsilon[it])),emode)
                                  elif numchannels>leftpoints:
                                      AA=dot(A[:,locpoints].conj().T,A[:,locpoints])
                                      epsilon[it] = max(absolute(eigvals(AA)))*self.eps_perc
                                      qi[s,locpoints]=dot(dot(inv(AA+inv(weights)*epsilon[it]),A[:,locpoints].conj().T),emode)                                                       
                                  if self.beta < 1 and it > 1:   
                                      #Reorder from the greatest to smallest magnitude to define a reduced-point source distribution , and reform a reduced transfer matrix 
                                      leftpoints=int(round(numpoints*self.beta**(it+1)))                                                                                          
                                      idx = argsort(abs(qi[s,locpoints]))[::-1]   
                                      #print(it, leftpoints, locpoints, idx )
                                      locpoints= delete(locpoints,[idx[leftpoints::]])             
                                      qix=zeros([n+m,leftpoints], dtype='complex128')                      
                                      qix[s,:]=qi[s,locpoints]
                                      #calc weights for next iteration 
                                      weights=diag(absolute(qix[s,:])**(2-self.pnorm))    
                                  else:                          
                                      weights=diag((absolute(qi[s,:])**(2-self.pnorm)))    
                               
                          elif self.method == 'InverseIRLS':                         
                              weights=eye(numpoints)
                              locpoints=arange(numpoints)
                              for it in arange(self.max_iter): 
                                  if numchannels<=numpoints: 
                                      wtwi=inv(dot(weights.T,weights))  
                                      aH=A.conj().T                       
                                      qi[s,:]=dot(dot(wtwi,aH),dot(inv(dot(A,dot(wtwi,aH))),emode))                            
                                      weights=diag(absolute(qi[s,:])**((2-self.pnorm)/2))
                                      weights=weights/sum(absolute(weights))                                 
                                  elif numchannels>numpoints:
                                      wtw=dot(weights.T,weights)
                                      qi[s,:]= dot(dot(inv(dot(dot(A.conj.T,wtw),A)),dot( A.conj().T,wtw)) ,emode)
                                      weights=diag(absolute(qi[s,:])**((2-self.pnorm)/2))
                                      weights=weights/sum(absolute(weights))                  
                          else:
                              locpoints=arange(numpoints) 
                              unit = self.unit_mult
                              AB = vstack([hstack([A.real,-A.imag]),hstack([A.imag,A.real])])
                              R  = hstack([emode.real.T,emode.imag.T]) * unit
                              if self.method == 'LassoLars':
                                  model = LassoLars(alpha=self.alpha * unit,max_iter=self.max_iter)
                              elif self.method == 'LassoLarsBIC':
                                  model = LassoLarsIC(criterion='bic',max_iter=self.max_iter)
                              elif self.method == 'OMPCV':
                                  model = OrthogonalMatchingPursuitCV()
                              elif self.method == 'LassoLarsCV':
                                  model = LassoLarsCV()                        
                              if self.method == 'NNLS':
                                  x , zz = nnls(AB,R)
                                  qi_real,qi_imag = hsplit(x/unit, 2) 
                              else:
                                  model.fit(AB,R)
                                  qi_real,qi_imag = hsplit(model.coef_[:]/unit, 2)
                              #print(s,qi.size)    
                              qi[s,locpoints] = qi_real+qi_imag*1j
                      else:
                          warn('Eigenvalue %g <= 0 for frequency index %g. Will not be calculated!' % (s, i),Warning, stacklevel = 2)
                      #Generate source maps of all selected eigenmodes, and superpose source intensity for each source type.
                  temp = zeros(numpoints)
                  temp[locpoints] = sum(absolute(qi[:,locpoints])**2,axis=0)
                  ac[i] = temp
                  fr[i] = 1    
  
  def L_p ( x ):
      """
      Calculates the sound pressure level from the squared sound pressure.
      
      :math:`L_p = 10 \lg ( x / 4\cdot 10^{-10})`
      
      Parameters
      ----------
      x: array of floats
          The squared sound pressure values
          
      Returns
      -------
      array of floats
          The corresponding sound pressure levels in dB. 
          If `x<0`, -350.0 dB is returned.
      """
      # new version to prevent division by zero warning for float32 arguments
      return 10*log10(clip(x/4e-10,1e-35,None))
  #    return where(x>0, 10*log10(x/4e-10), -1000.)
  
  def integrate(data, grid, sector):
      """
      Integrates a sound pressure map over a given sector.
      
      This function can be applied on beamforming results to
      quantitatively analyze the sound pressure in a given sector.
      If used with :meth:`Beamformer.result()<acoular.fbeamform.BeamformerBase.result>`, 
      the output is identical to the result of the intrinsic 
      :meth:`Beamformer.integrate<acoular.fbeamform.BeamformerBase.integrate>` method.
      It can, however, also be used with the 
      :meth:`Beamformer.synthetic<acoular.fbeamform.BeamformerBase.synthetic>`
      output.
      
      Parameters
      ----------
      data: array of floats
          Contains the calculated sound pressures in Pa.        
          If data has the same number of entries than the number of grid points
          only one value is returned.
          In case of a 2-D array with the second dimension identical 
          to the number of grid points an array containing as many entries as
          the first dimension is returned.
      grid: Grid object 
          Object of a :class:`~acoular.grids.Grid`-derived class 
          that provides the grid locations.        
      sector: array of floats
          Tuple with arguments for the `indices` method 
          of a :class:`~acoular.grids.Grid`-derived class 
          (e.g. :meth:`RectGrid.indices<acoular.grids.RectGrid.indices>` 
          or :meth:`RectGrid3D.indices<acoular.grids.RectGrid3D.indices>`).
          Possible sectors would be `array([xmin, ymin, xmax, ymax])`
          or `array([x, y, radius])`.
            
      Returns
      -------
      array of floats
          The spectrum (all calculated frequency bands) for the integrated sector.
      """
      
      if isinstance(sector, Sector):
          ind = grid.subdomain(sector)
      else:
          ind = grid.indices(*sector)
      
      gshape = grid.shape
      gsize = grid.size
      if size(data) == gsize: # one value per grid point
          h = data.reshape(gshape)[ind].sum()
      elif data.ndim == 2 and data.shape[1] == gsize:
          h = zeros(data.shape[0])
          for i in range(data.shape[0]):
              h[i] = data[i].reshape(gshape)[ind].sum()
      return h
  
apis: sklearn.linear_model.LassoLarsCV sklearn.linear_model.OrthogonalMatchingPursuitCV sklearn.linear_model.LassoLars sklearn.linear_model.LassoLarsIC