|
#!/usr/bin/python # -*- coding: utf-8 -*- except ImportError: from scipy import fft
trimMatrix, derampMatrix, squareMatrix)
(1./2000, num.inf), (1./2000, 1./500), (1./500, 1./10), (0, num.inf)]
"""Exponential function model to approximate a positive-definite covariance
We assume the following simple covariance model to describe the empirical noise observations:
.. math::
cov(d) = c \\cdot e^{\\frac{-d}{b}}
:param distance: Distance between :type distance: float or :class:`numpy.ndarray` :param a: Linear model coefficient :type a: float :param b: Exponential model coefficient :type b: float :returns: Covariance at ``distance`` :rtype: :class:`numpy.ndarray` """
r"""Exponential function model to approximate a positive-definite covariance
We assume the following simple covariance model to describe the empirical noise observations:
.. math::
cov(d) = c \\cdot e^{\\frac{-d}{b}} \\cdot \cos{\\frac{d-c}{d}}
:param distance: Distance between :type distance: float or :class:`numpy.ndarray` :param a: Linear model coefficient :type a: float :param b: Exponential model coefficient :type b: float :param c: Cosinus distance correction :type c: float :param c: Cosinus coefficient :type c: float :returns: Covariance at ``distance`` :rtype: :class:`numpy.ndarray` """ return a * num.exp(-distance/b) * num.cos((distance-c)/d)
"""Exponential linear model to estimate a log-linear power spectrum
We assume the following log-linear model for a measured power spectrum:
.. math::
pow(k) = \\frac{k^\\beta}{D}
:param k: Wavenumber :type k: float or :class:`numpy.ndarray` :param a: Exponential model factor :type a: float :param b: Fractional model factor :type b: float """ return (k**beta)/D
shape=(None,), dtype=num.float, serialize_as='list', optional=True, help='Noise patch coordinates and size,') optional=True, help='Covariance model coefficients. Either two (exponential) ' 'or three (exponential and cosine term) coefficients.' 'See also :func:`~kite.covariance.modelCovariance`.') choices=['exponential', 'exponential_cosine'], default='exponential', help='Covariance approximation function.') choices=['spectral', 'spatial'], default='spatial', help='Method for estimating the covariance and structure function.') default=75, help='Number of distance bins for spatial covariance sampling.') default=200000, help='Number of random pairs for spatial covariance sampling.') optional=True, help='Variance of the model.') default=True, help='Adaptive subsampling flag for full covariance calculation.') dtype=num.float, optional=True, serialize_as='base64', help='Cached covariance matrix, ' 'see :attr:`~kite.Covariance.covariance_matrix`.')
kwargs['model_coefficients'] = ( kwargs.pop('a'), kwargs.pop('b'))
"""Construct the variance-covariance matrix of quadtree subsampled data.
Variance and covariance estimates are used to construct the weighting matrix to be used later in an optimization.
Two different methods exist to propagate full-resolution data variances and covariances of :class:`kite.Scene.displacement` to the covariance matrix of the subsampled dataset:
1. The distance between :py:class:`kite.quadtree.QuadNode` leaf focal points, :py:class:`kite.covariance.Covariance.matrix_focal` defines the approximate covariance of the quadtree leaf pair. 2. The _accurate_ propagation of covariances by taking the mean of every node pair pixel covariances. This process is computational very expensive and can take a few minutes. :py:class:`kite.covariance.Covariance.matrix_focal`
:param quadtree: Quadtree to work on :type quadtree: :class:`~kite.Quadtree` :param config: Config object :type config: :class:`~kite.covariance.CovarianceConfig` """
return self.getLeafCovariance(*args, **kwargs)
""" Sets and updated the config of the instance
:param config: New config instance, defaults to configuration provided by parent :class:`~kite.Scene` :type config: :class:`~kite.covariance.CovarianceConfig`, optional """ config = self.scene.config.covariance
self._log.warning('Old format - resetting noise patch coordinates') config.covariance_matrix = None config.noise_coord = None
and (config.model_coefficients is not None or config.variance is not None): self.noise_data # init data array self.config.model_coefficients = config.model_coefficients self.config.variance = config.variance
def finished_combinations(self): return covariance_ext.get_finished_combinations()
def noise_coord(self): """ Coordinates of the noise patch in local coordinates.
:setter: Set the noise coordinates :getter: Get the noise coordinates :type: :class:`numpy.ndarray`, ``[llE, llN, sizeE, sizeN]`` """ self.noise_data
def noise_coord(self, values):
def noise_patch_size_km2(self): """ :getter: Noise patch size in :math:`km^2`. :type: float """ if self.noise_coord is None: return 0. size = (self.noise_coord[2] * self.noise_coord[3])*1e-6 if self.noise_data.size < self.NOISE_PATCH_MIN_PX: self._log.warning('Defined noise patch is instably small') return size
def noise_data(self, data): """ Noise data we process to estimate the covariance
:setter: Set the noise patch to analyze the covariance. :getter: If the noise data has not been set manually, we grab data through :func:`~kite.Covariance.selectNoiseNode`. :type: :class:`numpy.ndarray` """ return self._noise_data
def noise_data(self): *self.config.noise_coord[:2]) *self.config.noise_coord[2:])
else: node.sizeE, node.sizeN]
def noise_data(self, data):
def noise_data_gridE(self):
def noise_data_gridN(self):
""" Choose noise node from quadtree the biggest :class:`~kite.quadtree.QuadNode` from :class:`~kite.Quadtree`.
:returns: A quadnode with the least signal. :rtype: :class:`~kite.quadtree.QuadNode` """
if n.npixel > NOISE_PATCH_MIN_PX and n.nan_fraction < NOISE_PATCH_MAX_NAN]
% (time.time() - t0))
""" Helper function returning appropriate :class:`~kite.quadtree.QuadNode` and for maintaining the internal mapping with the matrices.
:param nx: matrix x position :type nx: int :param ny: matrix y position :type ny: int :returns: tuple of :class:`~kite.quadtree.QuadNode` s for ``nx`` and ``ny`` :rtype: tuple """
if self.config.covariance_matrix is None: return False return True
def covariance_matrix(self): """ Covariance matrix calculated from mean of all pixel pairs inside the node pairs (full and accurate propagation).
:type: :class:`numpy.ndarray`, size (:class:`~kite.Quadtree.nleaves` x :class:`~kite.Quadtree.nleaves`) """ self._calcCovarianceMatrix(method='full') self.config.covariance_matrix.reshape(nl, nl) except ValueError: self.config.covariance_matrix = None return self.covariance_matrix
def covariance_matrix_focal(self): """ Approximate Covariance matrix from quadtree leaf pair distance only. Fast, use for intermediate steps only and finallly use approach :attr:`~kite.Covariance.covariance_matrix`.
:type: :class:`numpy.ndarray`, size (:class:`~kite.Quadtree.nleaves` x :class:`~kite.Quadtree.nleaves`) """ return self._calcCovarianceMatrix(method='focal')
def weight_matrix(self): """ Weight matrix from full covariance :math:`cov^{-1}`.
:type: :class:`numpy.ndarray`, size (:class:`~kite.Quadtree.nleaves` x :class:`~kite.Quadtree.nleaves`) """ return num.linalg.inv(self.covariance_matrix)
def weight_matrix_L2(self): """ Weight matrix from full covariance :math:`\\sqrt{cov^{-1}}`.
:type: :class:`numpy.ndarray`, size (:class:`~kite.Quadtree.nleaves` x :class:`~kite.Quadtree.nleaves`) """ incov = num.linalg.inv(self.covariance_matrix) return sp.linalg.sqrtm(incov)
def weight_matrix_focal(self): """ Approximated weight matrix from fast focal method :math:`cov_{focal}^{-1}`.
:type: :class:`numpy.ndarray`, size (:class:`~kite.Quadtree.nleaves` x :class:`~kite.Quadtree.nleaves`) """ try: return num.linalg.inv(self.covariance_matrix_focal) except num.linalg.LinAlgError as e: self._log.exception(e) return num.eye(self.covariance_matrix_focal.shape[0])
def weight_vector(self): """ Weight vector from full covariance :math:`cov^{-1}`. :type: :class:`numpy.ndarray`, size (:class:`~kite.Quadtree.nleaves`) """ return num.sum(self.weight_matrix, axis=1)
def weight_vector_focal(self): """ Weight vector from fast focal method :math:`\\sqrt{cov_{focal}^{-1}}`. :type: :class:`numpy.ndarray`, size (:class:`~kite.Quadtree.nleaves`) """ return num.sum(self.weight_matrix_focal, axis=1)
"""Constructs the covariance matrix.
:param method: Either ``focal`` point distances are used - this is quick but only an approximation. Or ``full``, where the full quadtree pixel distances matrices are calculated , defaults to ``focal`` :type method: str, optional :returns: Covariance matrix :rtype: thon:numpy.ndarray """
model = self.getModelFunction()
coords = self.quadtree.leaf_focal_points_meter dist_matrix = num.sqrt( (coords[:, 0] - coords[:, 0, num.newaxis])**2 + (coords[:, 1] - coords[:, 1, num.newaxis])**2) cov_matrix = model(dist_matrix, *self.covariance_model)
# adding variance if self.variance < cov_matrix.max(): variance = cov_matrix.max() else: variance = self.variance if self.quadtree.leaf_mean_px_var is not None: self._log.debug( 'Adding variance from scene.displacement_px_var') variance += self.quadtree.leaf_mean_px_var num.fill_diagonal(cov_matrix, variance)
for nx, ny in num.nditer(num.triu_indices_from(dist_matrix)): self._mapLeaves(nx, ny)
dtype=num.uint32) leaf._slice_rows.stop) leaf._slice_cols.stop)
self.scene.frame.gridEmeter.filled(), self.scene.frame.gridNmeter.filled(), leaf_map, self.covariance_model, self.variance, nthreads, self.config.adaptive_subsampling)\ .reshape(nleaves, nleaves)
self._log.debug( 'Adding variance from scene.displacement_px_var') cov_matrix[num.diag_indices_from(cov_matrix)] +=\ self.quadtree.leaf_mean_px_var
else: raise TypeError('Covariance calculation %s method not defined!' % method)
(method, time.time()-t0))
self._log.debug('Checking whether matrix is positive-definite') if full: matrix = self.covariance_matrix else: matrix = self.covariance_matrix_focal
try: chol_decomp = num.linalg.cholesky(matrix) except num.linalg.linalg.LinAlgError: pos_def = False else: pos_def = ~num.all(num.iscomplex(chol_decomp)) finally: if not pos_def: self._log.warning('Covariance matrix is not positiv definite!') return pos_def
def _leafFocalDistance(leaf1, leaf2): return num.sqrt((leaf1.focal_point[0] - leaf2.focal_point[0])**2 + (leaf1.focal_point[1] - leaf2.focal_point[1])**2)
if not isinstance(leaf1, str): leaf1 = leaf1.id if not isinstance(leaf2, str): leaf2 = leaf2.id if not self._initialized: self.covariance_matrix_focal try: return self._leaf_mapping[leaf1], self._leaf_mapping[leaf2] except KeyError as e: raise KeyError('Unknown quadtree leaf with id %s' % e)
"""Get the covariance between ``leaf1`` and ``leaf2`` from distances.
:param leaf1: Leaf one :type leaf1: str of `leaf.id` or :class:`~kite.quadtree.QuadNode` :param leaf2: Leaf two :type leaf2: str of `leaf.id` or :class:`~kite.quadtree.QuadNode` :returns: Covariance between ``leaf1`` and ``leaf2`` :rtype: float """ return self.covariance_matrix[self._leafMapping(leaf1, leaf2)]
""" Get the total weight of ``leaf``, which is the summation of all single pair weights of :attr:`kite.Covariance.weight_matrix`.
.. math ::
w_{x} = \\sum_i W_{x,i}
:param model: ``Focal`` or ``full``, default ``focal`` :type model: str :param leaf: A leaf from :class:`~kite.Quadtree` :type leaf: :class:`~kite.quadtree.QuadNode`
:returns: Weight of the leaf :rtype: float """ (nl, _) = self._leafMapping(leaf, leaf) weight_mat = self.weight_matrix_focal return num.mean(weight_mat, axis=0)[nl]
anisotropic=False, rstate=None): """Create random synthetic noise from data noise power spectrum.
This function uses the power spectrum of the data noise (:attr:`noise_data`) (:func:`powerspecNoise`) to create synthetic noise, e.g. to use it for data pertubation in optinmizations. The default sampling distances are taken from :attr:`kite.scene.Frame.dE` and :attr:`kite.scene.Frame.dN`. They can be overwritten.
:param shape: shape of the desired noise patch. Pixels in northing and easting (`nE`, `nN`), defaults to `(1024, 1024)`. :type shape: tuple, optional :param dEdN: The sampling distance in east and north [m], defaults to (:attr:`kite.scene.Frame.dEmeter`, :attr:`kite.scene.Frame.dNmeter`). :type dEdN: tuple, floats :returns: synthetic noise patch :rtype: :class:`numpy.ndarray` """ # self._log.warning('Patch dimensions must be even, ' # 'ceiling dimensions!') pass
rstate = num.random.RandomState()
elif anisotropic: interp_pspec, _, _, _, skE, skN = self.powerspecNoise3D() kE = num.fft.fftshift(kE) kN = num.fft.fftshift(kN) mkE = num.logical_and(kE >= skE.min(), kE <= skE.max()) mkN = num.logical_and(kN >= skN.min(), kN <= skN.max()) mkRad = num.where( # noqa k_rad < num.sqrt(kN[mkN].max()**2 + kE[mkE].max()**2)) res = interp_pspec(kN[mkN, num.newaxis], kE[num.newaxis, mkE], grid=True) amp = res amp = num.fft.fftshift(amp)
# remove shape % 2 padding
'''Create noise for a :class:`~kite.quadtree.Quadtree`
Use :meth:`~kite.covariance.Covariance.getSyntheticNoise` to create data-driven noise on each quadtree leaf, summarized by
:param gather: Function gathering leaf's noise realisation, defaults to num.median. :type normalisation: callable, optional :returns: Array of noise level at each quadtree leaf. :rtype: :class:`numpy.ndarray` '''
shape=self.scene.displacement.shape, rstate=rstate)
syn_noise[lv._slice_rows, lv._slice_cols])
if self._powerspec1d_cached is None: self._powerspec1d_cached = self._powerspecNoise( data, norm='1d', ndeg=ndeg, nk=nk) return self._powerspec1d_cached
data, norm='2d', ndeg=ndeg, nk=nk)
if self._powerspec3d_cached is None: self._powerspec3d_cached = self._powerspecNoise( data, norm='3d') return self._powerspec3d_cached
"""Get the noise power spectrum from :attr:`kite.Covariance.noise_data`.
:param data: Overwrite Covariance.noise_data, defaults to `None` :type data: :class:`numpy.ndarray`, optional :returns: `(power_spec, k, f_spectrum, kN, kE)` :rtype: tuple """ else: noise = data.copy() raise AttributeError('norm must be 1d, 2d or 3d')
# noise = squareMatrix(noise)
d=self.quadtree.frame.dEmeter)) d=self.quadtree.frame.dNmeter))
# def power1d(k): # theta = num.linspace(-num.pi, num.pi, ndeg, False) # power = num.empty_like(k) # for i in range(k.size): # kE = num.cos(theta) * k[i] # kN = num.sin(theta) * k[i] # power[i] = num.median(power_interp.ev(kN, kE)) * k[i]\ # * num.pi * 4 # return power
theta = num.linspace(-num.pi, num.pi, ndeg, False) power = num.empty_like(k)
cos_theta = num.cos(theta) sin_theta = num.sin(theta) for i in range(k.size): kE = cos_theta * k[i] kN = sin_theta * k[i] power[i] = num.mean(power_interp.ev(kN, kE))
power *= 2 * num.pi return power
""" Mean 2D Power works! """
# Median is more stable than the mean here
return power_interp
elif norm == '3d': power = power3d
k_rad.max(), nk)
"""Calculating the cosine transform of the power spectrum.
The cosine transform of the power spectrum is an estimate of the data covariance (see Hanssen, 2001).""" cos = fft.idct(p_spec, type=3) return cos
""" Set the sampling method """ assert method in CovarianceConfig.sampling_method.choices
self.config.sampling_method = method self._clear(config=True, spectrum=False) self.evChanged.notify() self._log.debug('Changed sampling method to %s' % method)
""" Set number of spatial bins """ self.config.spatial_bins = nbins self._clear(config=True, spectrum=False) self.evChanged.notify() self._log.debug('Changed spatial distance bins to %s' % nbins)
""" Set number of random spatial pairs """ self.config.spatial_pairs = npairs self._clear(config=True, spectrum=False) self.evChanged.notify() self._log.debug('Changed random pairs to %s' % npairs)
assert model in CovarianceConfig.model_function.choices self.config.model_function = model self._clear(config=True, spectrum=True) self.evChanged.notify() self._log.debug('Changed model function to %s' % model)
if self.config.model_function == 'exponential_cosine': return modelCovarianceExponentialCosine
def covariance_spectral(self): """ Covariance function estimated directly from the power spectrum of displacement noise patch using the cosine transform.
:type: tuple, :class:`numpy.ndarray` (covariance, distance) """ power_spec, k, dk, _, _, _ = self.powerspecNoise1D() # power_spec -= self.variance
d = num.arange(1, power_spec.size+1) * dk cov = self._powerCosineTransform(power_spec)
return cov, d
def covariance_spatial(self):
abs(grdN.min() - grdN.max()))
# Select random coordinates
(grdE[idx0] - grdE[idx1])**2)
bin_distances[~num.isnan(variance)]) bin_distances[~num.isnan(covariance)])
""" Calculate the covariance function
:return: The covariance and distance :rtype: tuple """ elif self.config.sampling_method == 'spectral': return self.covariance_spectral
""" Covariance model parameters for :func:`~kite.covariance.modelCovariance` retrieved from :attr:`~kite.Covariance.getCovarianceFunction`.
.. note:: using this function implies several model fits: (1) fit of the spectrum and (2) the cosine transform. Not sure about the consequences, if this is useful and/or meaningful.
:getter: Get the parameters. :type: tuple, ``a`` and ``b`` """
elif self.config.model_function == 'exponential_cosine': coeff = (num.mean(covariance), num.mean(distance), num.mean(distance)*-.1, .1)
func = self.getModelFunction()
# Testing penalty function def model(*args): distance, a, b, c, d = args res = func(*args)
penalty = 0. if distance[-1]/b > (distance[-1]+c)/d: penalty = (b-d) * coeff[0] self._log.warning('Penalty %f' % penalty)
return res + penalty
# Overwrite with pure model function model = self.getModelFunction() # noqa
model, distance, covariance, p0=coeff) except (RuntimeError, TypeError) as e: self._log.exception(e) self._log.warning('Could not fit the %s covariance model', self.config.model_function) finally:
def covariance_model_rms(self): """ :getter: RMS missfit between :class:`~kite.Covariance.covariance_model` and :class:`~kite.Covariance.getCovarianceFunction` :type: float """ cov, d = self.getCovariance() model = self.getModelFunction() cov_mod = model(d, *self.covariance_model)
return num.sqrt(num.mean((cov - cov_mod)**2))
def structure_spatial(self):
def structure_spectral(self): """ Structure function derived from ``noise_patch`` :type: tuple, :class:`numpy.ndarray` (structure_spectral, distance)
Adapted from http://clouds.eos.ubc.ca/~phil/courses/atsc500/docs/strfun.pdf """ power_spec, k, dk, _, _, _ = self.powerspecNoise1D() d = num.arange(1, power_spec.size+1) * dk
def structure_spectral(power_spec, d, k): struc_func = num.zeros_like(k) for i, d in enumerate(d): for ik, tk in enumerate(k): # struc_func[i] += (1. - num.cos(tk*d))*power_spec[ik] struc_func[i] += (1. - sp.special.j0(tk*d))*power_spec[ik] struc_func *= 2./1 return struc_func
struc_func = structure_spectral(power_spec, d, k) return struc_func, d
""" Get the structure function
:param method: Either `spatial` or `spectral`, if `None` the method is taken from config :type method: str (optional)
:return: (variance, distance) :rtype: tuple """ if method is None: method = self.config.sampling_method if method == 'spatial': return self.structure_spatial elif method == 'spectral': return self.structure_spectral
def variance(self): """ Variance of data noise estimated from the high-frequency end of power spectrum.
:setter: Set the variance manually :getter: Retrieve the variance :type: float """ return self.config.variance
def variance(self, value): self.config.variance = float(value) # self._clear(config=False, spectrum=False, spatial=False) self.evChanged.notify()
def variance(self):
self.config.sampling_method == 'spatial':
num.mean(structure_spatial[(last_20p):]))
self.config.sampling_method == 'spectral'): power_spec, k, dk, spectrum, _, _ = self.powerspecNoise1D() cov, _ = self.covariance_spectral ma = self.covariance_model[0] # print(cov[1]) ps = power_spec * spectrum.size # print(spectrum.size) # print(num.mean(ps[-int(ps.size/9.):-1])) var = num.median(ps[-int(ps.size/9.):]) + ma self.config.variance = float(var)
""" Export the full :attr:`~kite.Covariance.weight_matrix` to an ASCII file. The data can be loaded through :func:`numpy.loadtxt`.
:param filename: path to export to :type filename: str """ self._log.debug('Exporting Covariance.weight_matrix to %s' % filename) header = 'Exported kite.Covariance.weight_matrix, '\ 'for more information visit https://pyrocko.org\n'\ '\nThe matrix is symmetric and ordered by QuadNode.id:\n' header += ', '.join([l.id for l in self.quadtree.leaves]) num.savetxt(filename, self.weight_matrix, header=header)
sha = sha1() sha.update(str(self.config).encode()) return sha.digest().hex()
def plot(self): """ Simple overview plot to summarize the covariance estimations. """ from kite.plot2d import CovariancePlot return CovariancePlot(self)
def plot_syntheticNoise(self): """ Simple overview plot to summarize the covariance estimations. """ from kite.plot2d import SyntheticNoisePlot return SyntheticNoisePlot(self) |