"""
Part of comet/snmr/modelling
Nuclear magnetic resonance base manager as used by MRS and MRT manager classes
"""
# COMET - COupled Magnetic resonance Electrical resistivity Tomography
# Copyright (C) 2019 Nico Skibbe
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
# general modules to import according to standards
import os
import numpy as np
import pygimli as pg
from pygimli.utils.base import gmat2numpy
from comet.snmr import kernel as k
from comet.snmr.survey import Survey
from comet.pyhed import log
from comet.pyhed.plot import quantile
# local functions in package
from . errors import InputError
from . snmrModelling import MRS1dBlockQTModelling, SNMRModelling
from scipy.io import loadmat # loading Matlab mat files
import matplotlib.pyplot as plt
from comet.pyhed.misc import temporal_printoptions
[docs]class SNMRBase():
""" Manager base class for MRS and MRT manager classes.
"""
def __init__(self,
survey=None,
fid=0,
kernel=None,
mtype='block',
dtype='rotatedAmplitudes',
update_kernel=False,
dim=1,
**kwargs):
"""
Parameters
----------
kernel: pg.Matrix or ndarray or class-instance or string [None]
Basically a simple precalculated kernel matrix in the .mrsi, .mrsk
or .npz format is sufficient for this class to work. If a
recalcualtion of the kernelmatrix is intended, an instance of a
kernelclass forward operator can be loaded.
data: string or comet.survey.Survey
Comet survey class instance or path to survey file or mrsd file.
mtype: string [ 'block' ]
Model type for the inversion.
The inversion can be done unsing 'smooth' discretization (1D, 2D)
or 'block' discretization (1D only).
dtype: string [ 'rotatedAmplitudes' ]
Data type of the inversion.
Possible options for the used inversion scheme are 'complex',
'rotatedAmplitudes' or 'amplitudes'.
update_kernel: boolean [ False ]
If providedwith an initialized kernel class instance, a
recalculation of the kernel based on for example new resistivity
information is possible by setting this switch to True.
dim: integer [ 1 ]
1D or 2D inversion is supported at this state.
verbose: boolean [False]
Switch for console output.
"""
self.inv = None
self.fop = None
self.model = None
self._dtype = None
self.setDataType(dtype)
self._mtype = None
self.setModelType(mtype)
self.dimension = dim
self.complex = self.dtype == 'complex'
self.num_cpu = kwargs.pop('num_cpu', 12)
self.kernel = None
self.survey = None
self.setSurvey(survey, fid=fid)
if self.fid is None:
raise Exception(
'Dis not found corresponding sounding in survey '
'({}) with index {}'.format(self.survey, self.fid_index))
# kernel import
if kernel is not None:
self.setKernel(kernel, **kwargs)
else:
self.setKernel(self.survey.createKernel(
self.fid_index, dimension=self.dimension))
# sets self.kernel: kernel matrix
# sets self.ZVec: model: coordinates of the model [m]
# sets self.fid.pulses: pulse moments [As]
# initialising attributes for forward operator
self.setBoundsAndTrans()
# sets Boundaries for Transforamtion (start, min, max)
# self.thkBounds = [10., 1., 30.] # thickness
# self.wcBounds = [0.3, 0.0, 0.45] # water content
# self.t2Bounds = [0.2, 0.02, 1.0] # relaxation times
# self.trans = ['log', 'cot', 'log']
self.update_kernel = update_kernel
self.modelL = None
self.modelU = None
self._numPara = 2
def __str__(self): # string represantation of the class
'''
Print all information about the attributes: Kernel Matrix, ZVec,
pulses, gates, data and error.
Also gives information about initialised parts of the class, such as
forward operators, inversion shemes and kernle class initialisation
(is recalcualtion of kernel is needed/wanted)
'''
if self._hasSound and self.fid.pulses is not None:
with temporal_printoptions():
Ps = 'Pulsemoments: {0} ({0.size})'.format(self.fid.pulses)
else:
Ps = 'Pulsemoments not loaded yet.'
if self._hasSound and self.fid.gates is not None:
with temporal_printoptions():
Gs = 'Gates: {0} ({0.size})'.format(self.fid.gates)
else:
Gs = 'Gates not loaded yet.'
if not self._hasData:
Ds = 'No data loaded/calculated yet.'
Es = 'No error loaded yet.'
else:
absd = np.absolute(self.data)
Ds = 'data: shape, min, max: {}, {}, {}'.format(
np.shape(absd), np.min(absd), np.max(absd))
abse = np.absolute(self.error)
Es = 'error: shape, min, max: {}, {}, {}'.format(
np.shape(abse), np.min(abse), np.max(abse))
string = '{!r}\n{}\n{}\n{}\n{}'.format(self, Ps, Gs, Ds, Es)
string += f'\ndata type: {self.dtype}'
string += f'\nmodel type: {self.mtype}'
string += '\n### state of initialization ###'
if self.K is None:
string += '\nKernel matrix not loaded/calculated yet.'
else:
string += f'K.shape: {self.K.shape} , {self.K.dtype}'
fopis = False
if self.fop is not None:
string += '\nFOP initialised: %r' % (fopis)
foptype = self.fop.__class__.__name__
fopc = self.fop.iscomplex
string += ', {}, complex={}'.format(foptype, fopc)
string += '\nKernel ready for calculation: %r' % \
(self.kernel.updatable)
string += '\nInversion prepared: %r' % (self.inv is not None)
return string
def __repr__(self):
string = f'NMR base, {self.mtype}, {self.dtype}'
return string
@property
def dtype(self):
return self._dtype
[docs] def setDataType(self, dtype):
dtypes = ['complex', 'rotatedAmplitudes', 'amplitudes']
assert dtype in dtypes, f'dtype has to be in {dtypes}'
self._dtype = dtype
if self.fop is not None:
self.fop.dtype = self.dtype
@property
def mtype(self):
return self._mtype
[docs] def setModelType(self, mtype):
mtypes = ['block', 'smooth']
assert mtype in mtypes, f'mtype has to be in {mtypes}'
self._mtype = mtype
@property
def fid(self):
""" Reference to sounding (FID) class instance in survey."""
try:
fid = self.survey.fids[self.fid_index]
except IndexError:
# no fid as of yet
fid = None
return fid
@property
def _hasSound(self):
return self.fid is not None
@property
def data(self):
""" Data Vector representation with respect to self.dtype.
Returns None if no sounding or data_cube in sounding is found.
"""
if self._hasData:
if self.fid.rotated and self.dtype != 'rotatedAmplitudes':
# case 1/2: need to rotate data back
if self.fid.raw_rotated:
# if not solvable: raise error
raise Exception(
f'Cannot return data with dtype: "amplitudes". '
'Raw data are already rotated. '
'Please provide unrotated data or set to '
'dtype to "rotatedAmplitudes"')
else:
self.fid._regating()
elif not self.fid.rotated and self.dtype == 'rotatedAmplitudes':
# case 2/2: need to rotate data
self.fid.rotateAmplitudes()
# case 1/2: data are alright: prepare shape with respect to dtype
if self.dtype == 'rotatedAmplitudes':
# take real part of data
data, _ = self.fid.getRotatedAmplitudes()
elif self.dtype == 'amplitudes':
# take absolute of data
data = np.absolute(self.fid.data_gated).ravel()
elif self.dtype == 'complex':
# stich together real and imag part of data
# = also raveled, but twice the length!
# including rotation for fixed data phase (only complex)
rot_data = self.fid.getComplexData()
data = np.append(rot_data.real,
rot_data.imag)
else:
# case 2/2: no data at all
data = None
return data # either real valued array or None
@property
def _hasData(self):
if not self._hasSound:
has_data = False
elif self.fid.data_gated is None:
has_data = False
else:
has_data = True
return has_data
@property
def error(self):
if self._hasData:
# data come with error...
if self.dtype == 'complex':
error = self.fid.error_gated.ravel()
if len(error) == len(self.data)/2:
if np.iscomplexobj(error):
# take complex error
error = np.append(error.real, error.imag)
else:
# take real error for both, although this is bad if
# given error amplitudes
error = np.tile(error, 2)
elif self.dtype == 'amplitudes':
error = self.fid.error_gated.ravel()
if np.iscomplexobj(error):
error = np.absolute(error)
elif self.dtype == 'rotatedAmplitudes':
error = self.fid.error_gated.ravel()
if np.iscomplexobj(error):
error = np.real(error)
else:
# ... or they don't
error = None
return error
@property
def K(self):
"""
"""
return self.kernel.getKernel()
[docs] def loadSurvey(self, dataname):
"""
"""
self.survey = Survey()
self.survey.load(dataname)
[docs] def setSurvey(self, survey, fid=0):
"""
"""
self.survey = survey
if isinstance(survey, str):
self.importSurvey(survey)
elif survey is None:
self.survey = Survey()
else:
self.survey = survey
self.fid_index = fid
if self.kernel is not None:
# this is only True during the init of SNMRBase
# the kernel creation is then handled explicitly
# so the next line will be done one way or the other
self.kernel.setSurvey(self.survey, fid=self.fid_index)
def _createFOP_block(self, nlay):
"""
Create block forward operator instance.
Please use createFOP('block', nlay=...) instead of a direct call.
"""
self.nlay = nlay
# self.fop = SNMRModelling(
# self.survey,
# self.kernel,
# dtype=self.dtype,
# fid=self.fid_index,
# num_cpu=self.num_cpu,
# update_kernel=self.update_kernel,
# verbose=verbose)
# edit: 22.01.: Fehler in MRS1dBlockQTModelling gefunden. Erstmal weg,
# im Prinzip müsste SNMRModelling das sowieso alles können.
self.fop = MRS1dBlockQTModelling(self.survey,
fid=self.fid_index,
kernel=self.kernel,
nlay=nlay,
dtype=self.dtype)
def _createFOP_smooth(self):
"""
Create block forward operator instance.
Please use createFOP('smooth') instead of a direct call.
"""
self.fop = SNMRModelling(
self.survey,
self.kernel,
dtype=self.dtype,
fid=self.fid_index,
num_cpu=self.num_cpu,
update_kernel=self.update_kernel)
def _createINV_block(self, verbose=True, debug=False, **kwargs):
"""
Create inversion instance following a Marquardt inversion scheme
(block). Please use createINV instead of this method directly.
Use the attribute Type to 'block' (default) to set the
inversion and forward operator type to block.
"""
self.inv = pg.core.RInversion(pg.Vector(self.data), self.fop,
verbose=verbose,
dosave=debug)
# self.inv.dataVals = self.data # only for framework
self.inv.setMarquardtScheme(kwargs.pop('lambdaFactor', 0.8))
self.inv.setRobustData(kwargs.pop('robust', False))
def _createINV_smooth(self, verbose=True, debug=False, **kwargs):
"""
Create inversion instance with smoothness constraint.
Please use createINV instead of this method directly.
Set the attribute Type to 'smooth' (not default) to set the
inversion and forward operator type to smooth.
"""
self.inv = pg.core.RInversion(pg.Vector(self.data), self.fop,
verbose=verbose, dosave=debug)
if kwargs.pop('logTrans', True) is True:
self.transmodel = pg.core.TransLog() # has to be "self"
self.inv.setTransModel(self.transmodel)
if kwargs.pop('blockyModel', False) is True:
self.inv.setBlockyModel(True) # L1 Norm
[docs] def createFOP(self, nlay=3):
"""
Creates the forward operator (FOP). Two possibilities are supported:
block and smooth. The choice affects the inversion process and
therefore its results.
possible keyword argument for block FOP is 'nlay' to define the number
of layers, the FOP is calculating (default = 3).
"""
if not self._hasSound:
raise Exception(
'Error: Cannot create FOP without sounding.'
'Please create sounding in corresponing survey class')
if self.mtype == 'smooth':
self._createFOP_smooth()
elif self.mtype == 'block':
self._createFOP_block(nlay=nlay)
[docs] def createINV(self, lam=1000, **kwargs):
"""
Create inversion instance (and fop if necessary with nlay).
arguments
---------
lam: float [100]
Lambda factor for inversion.
verbose: bool [True]
Additional verbose decission, can be True, even if the rest of the
Manager should remain silent. Most information of the different
iterations is printed in the console. It's recommended to set
verbose in this case to True (default).
special kwargs for Marquardt scheme (block)
-------------------------------------------
lambdaFactor, float [0.8]
Sets lambda factor for Marquardt scheme.
robust, bool [False]
Sets the robust flag for the data.
See pg.RInversion for more details
special kwargs for smooth scheme
---------------------------------
logTrans, bool [True]
Applies a logarithmic transformation to the data. Its recommended
to do so (default), due to the dealing with water contents,
which can't be negative. Logarithmic transformation is the easiest
way to archieve that.
blockyModel, bool [False]
Instead of the standard L2-Norm a L1 Norm can be used to allow for
more blocky models. Heavy changes in watercontent and relaxation
times can sometimes be fitted better this way.
"""
if self.fop is None:
self.createFOP(nlay=kwargs.pop('nlay', 3))
if self.mtype == 'block':
self._createINV_block(**kwargs)
elif self.mtype == 'smooth':
self._createINV_smooth(**kwargs)
self.applyBoundsAndTrans()
self.inv.setLambda(lam)
self.inv.stopAtChi1(kwargs.pop('stopAtChi1', True))
self.inv.setDeltaPhiAbortPercent(
kwargs.pop('deltaPhiAbortPercent', 1))
self.inv.setAbsoluteError(kwargs.pop('absoluteError', self.error))
self.inv.setMaxIter(kwargs.pop('maxIter', 15))
[docs] def setKernel(self, kernelfile=None, load_loopmesh=True,
load_kernelmesh=True, use_order_refinement=True):
"""
Load or initialize a new Kernel class instance for calcualting the
NMR kernels.
"""
if isinstance(kernelfile, str): # load a given file
if kernelfile.endswith(('.mrsk')):
# case 1/5: Use survey class to import mrsk and set Kernel
log.debug('setKernel: loadMRSK')
self.kernel = self.survey.loadMRSK(
kernelfile, tx=self.fid.tx_index,
rx=self.fid.rx_index)
else:
# case 2/5: load kernel class from numpy archive
log.debug('setKernel: load kernel')
self.kernel = k.Kernel(self.survey, fid=self.fid_index)
self.kernel.load(kernelfile, load_loopmesh=load_loopmesh,
load_kernelmesh=load_kernelmesh,
use_order_refinement=use_order_refinement)
elif isinstance(kernelfile, k.Kernel):
log.debug('setKernel: set kernel')
# case 3/5: already given an initialized kernel class
self.kernel = kernelfile
elif isinstance(kernelfile, np.ndarray):
log.debug('setKernel: set numpy matrix')
# case 4/5: given simple matrix
self.kernel = k.Kernel(self.survey, fid=self.fid_index)
self.kernel.K = kernelfile
elif kernelfile is None:
# case 5/5: new empty instance
log.debug('setKernel: create new Kernel')
self.kernel = self.survey.createKernel(fid=self.fid_index,
dimension=self.dimension)
else:
raise InputError(kernelfile)
[docs] def loadMRSI(self, filename, verbose=True):
""" Load data, error and kernel from mrsi file """
raise Exception('loadMRSI not working after survey class creation.')
print('load file: {}'.format(os.path.abspath(filename)))
idata = loadmat(filename, struct_as_record=False,
squeeze_me=True)['idata']
self.errorWeights = np.sqrt(idata.data.gateL)
self.setGates(idata.data.t + idata.data.effDead,
midpoints=True)
self.setPulses(idata.data.q)
self.setZVector(-np.abs(idata.kernel.z))
self.K = np.conjugate(idata.kernel.K) # convention
self.dcube = idata.data.dcube
self.ecube = idata.data.ecube
print('phase information,', idata.data.phi, np.shape(idata.data.phi))
# load model from matlab file (result of MRSQTInversion) if any
if hasattr(idata, 'inv1Dqt'):
if hasattr(idata.inv1Dqt, 'blockMono'):
sol = idata.inv1Dqt.blockMono.solution
self.model = np.hstack((sol.thk, sol.w, sol.T2))
self.nlay = len(sol.w)
self.setDataAndError() # checks also error
if verbose:
print("loaded file: " + filename)
[docs] def applyBoundsAndTrans(self):
"""
Append the previously given bounaries for the model transformation
to the forward operator.
"""
if self.fop is None:
raise Exception('initialize forward operator before appending '
'proper boundaries for the model transformation')
# set transformation and boundaries
if self.mtype == 'smooth':
# water content
self.fop.region(0).setParameters(self.wcBounds[0],
self.wcBounds[1],
self.wcBounds[2], self.trans[1])
# relaxation times
self.fop.region(1).setParameters(self.t2Bounds[0],
self.t2Bounds[1],
self.t2Bounds[2], self.trans[2])
elif self.mtype == 'block':
# thickness
self.fop.region(0).setParameters(self.thkBounds[0],
self.thkBounds[1],
self.thkBounds[2], self.trans[0])
# water content
self.fop.region(1).setParameters(self.wcBounds[0],
self.wcBounds[1],
self.wcBounds[2], self.trans[1])
# relaxation times
self.fop.region(2).setParameters(self.t2Bounds[0],
self.t2Bounds[1],
self.t2Bounds[2], self.trans[2])
[docs] def setBoundsAndTrans(self,
thkBounds=[10., 1., 30.],
wcBounds=[0.3, 0.0, 0.7],
t2Bounds=[0.2, 0.005, 1.0],
trans=['log', 'cot', 'log']):
"""
Sets the boundarys and transformation for the model domain.
Parameters
----------
thkBounds: list of floats [ [10., 1., 30.] ]
Startvalue, lower and upper boundary for thickness of each layer
in 1D. Ignored for smooth models (or 2D).
wcBounds: list of floats [ [0.3, 0.0, 0.7] ]
Startvalue, lower and upper boundary for water content.
t2Bounds: list of floats [ [0.2, 0.005, 1.0] ]
Startvalue, lower and upper boundary for relaxation times.
trans: list of strings [ ['log', 'cot', 'log'] ]
Defines the type of model transformation. logarithmic ('log') or
cotangens ('cot')
"""
self.thkBounds = thkBounds # thk(start, min, max)
self.wcBounds = wcBounds # wc(start, min, max)
self.t2Bounds = t2Bounds # t2(start, min, max)
self.logpar = True
self.trans = trans
if self.fop is not None:
self.applyBoundsAndTrans()
[docs] def calcModelCovarianceMatrix(self):
"""Compute linear model covariance matrix."""
J = gmat2numpy(self.fop.jacobian()) # (linear) jacobian matrix
D = np.diag(1 / self.error)
DJ = D.dot(J)
JTJ = DJ.T.dot(DJ)
MCM = np.linalg.inv(JTJ) # model covariance matrix
var = np.sqrt(np.diag(MCM)) # standard deviations from main diagonal
di = (1. / var) # variances as column vector
# scaled model covariance (=correlation) matrix
MCMs = di.reshape(len(di), 1) * MCM * di
return var, MCMs
[docs] def calcModelCovarianceMatrixBounds(self):
"""Compute model bounds using covariance matrix diagonals."""
self.mcm = self.calcMCM()[0]
self.modelL = self.model - self.mcm
self.modelU = self.model + self.mcm
[docs] def simulate(self, model, err=250e-9, samplingrate=1000., max_time=1.0,
num_gates=50, verbose=False, debug=False,
**kwargs):
"""
Creates forward operator and calculates a synthetic response to a given
model. Keyword arguments are passed to the function createFOP and to
FOP.response. You can also define the 'Type' to be 'smooth' or 'block'
or let the simulate function analyse the input.
returns datacube, errorcube (both complex) and phaseinformation
(for rotated amplitudes)
Parameters
----------
model: list of lists
Given model of shape [water_content, relaxation_time] if forwarded
to FOP to generate synthetic data set.
"""
# save old fop
initial_fop = self.fop
if err > 1e-5:
log.warning('Very high error detected in simulate: {} Volt'
.format(err))
# initialize FOP solely for simulation
self.fop = None
self.createFOP()
# simulation without frequency offset
self.fid.setFrequencyOffset(0.)
# dont be confused: forward calculation is always on gates
# so in this case we use sampling times as gates and do the real gating
# later
# 04.05.2020: set gates (or times) always without deadtime.
self.fid.setGates(
np.arange(0, max_time + 1./samplingrate,
1./samplingrate), midpoints=False)
# calculate synthetic set of data (complex, raveled)
data = self.fop.forward(model, verbose=verbose, **kwargs)
nmp = len(self.fid.pulses)
nmg = len(self.fid.gates)
# unravel
dcube = data.reshape(nmp, nmg)
log.debug('data: {}'.format(dcube))
ecube = np.ones_like(dcube, dtype=float) * err
log.debug('error: {}'.format(ecube))
phi = getPhiByGridSearch(dcube)
log.info('phases: {}'.format(phi))
# append Noise, random Gaussian (== normal) distribution,
# separately for imag and real
err_i = np.array(
np.random.randn(nmp * nmg)).reshape(nmp, nmg) * err
err_r = np.array(
np.random.randn(nmp * nmg)).reshape(nmp, nmg) * err
Noise = (err_r + 1j * err_i)
if debug:
debug_save = '{}raw_data.npz'.format(
kwargs.pop('raw_debug_praefix', ''))
log.debug('export raw data and noise for debugging: "{}"'
.format(debug_save))
np.savez(debug_save, data=dcube, noise=Noise,
error=ecube, gates=self.fid.gates,
pulses=self.fid.pulses, phases=phi)
dcube += Noise # complex (pulses, gates)
# pre 26.02.2020
# self.fid.setRawDataErrorAndTimes(
# dcube, ecube, self.fid.gates, phases=phi)
# post 26.02.2020: calc on times deadtime-end but raw times start at 0
self.fid.setRawDataErrorAndTimes(
dcube, ecube, self.fid.gates - self.fid.deadtime, phases=phi)
self.fid.gating(num_gates=num_gates)
# set old fop again
self.fop = initial_fop
# return data, error and times
return self.fid.data_gated, self.fid.error_gated, self.fid.phi
[docs] def invert(self, data=None, error=None, phase=None, lam=1000,
runChi1=False, **kwargs):
"""
# TODO!
"""
if data is not None:
self.setDataAndError(data, error, phase)
self.createINV(lam=lam, **kwargs)
log.info('INV: invert(): lam: {}, dtype: {}, mtype: {}'
.format(self.inv.getLambda(), self.dtype, self.mtype))
if runChi1 is True:
self.model = self.inv.runChi1()
elif runChi1 is False:
self.model = self.inv.run()
# TODO! model history!
# if self.modelhistory:
# self.modelhistory.append(self.model)
[docs] def showCube(self, ax=None, vec=None, islog=None, clim=None, clab=None,
cmap='viridis', cbar=True):
"""Plot any data (or response, error, misfit) cube nicely."""
if vec is None:
if not self._hasData:
raise Exception('No data found.')
else:
vec = self.data
if isinstance(vec, pg.Vector):
vec = vec.array()
if len(vec) == len(self.fid.pulses) * len(self.fid.gates) * 2:
# if not real or imag is given, make absolute
vec_2 = int(len(vec)/2)
vec = np.absolute(vec[:vec_2] + 1j*vec[vec_2:]) * \
np.sign(vec[:vec_2])
mul = 1.0
if np.max(vec) < 1e-3: # Volts
mul = 1e9
if ax is None:
fig, ax = plt.subplots(1, 1)
if islog is None:
islog = (np.min(vec) > 0.)
negative = quantile(vec, perc=0.01) < 0
# negative = (np.min(vec) < 0)
if islog:
vec = np.log10(np.abs(vec))
if clim is None:
if negative:
cmax = np.max(np.absolute(vec))
clim = (-cmax, cmax)
else:
cmax = np.max(vec)
if islog:
cmin = cmax - 1.5
else:
cmin = 0.
clim = (cmin, cmax)
xt = range(0, len(self.fid.gates), 10)
xtl = [str(ti) for ti in np.round(self.fid.gates[xt] * 1000.)]
qt = range(0, len(self.fid.pulses), 5)
qtl = [str(qi) for qi in np.round(self.fid.pulses[qt] * 10.) / 10.]
# plot
mat = vec.reshape((len(self.fid.pulses),
len(self.fid.gates))) * mul
im = ax.imshow(mat, interpolation='nearest', aspect='auto',
cmap=plt.get_cmap(cmap))
im.set_clim(clim)
ax.set_xticks(xt)
ax.set_xticklabels(xtl)
ax.set_yticks(qt)
ax.set_yticklabels(qtl)
ax.set_xlabel('time (ms)')
ax.set_ylabel('pulse moment (As)')
if cbar:
cb = plt.colorbar(im, ax=ax, orientation='horizontal')
if clab is not None:
cb.ax.set_title(clab)
return im, cb
return im
[docs]def effectiveNoise(area, noise_lvl=3.6e-3, sample_rate=1000., time=1.0):
""" Calculates the effective noise of a loop for simulation.
noise_lvl = 3.6e-3 nV / m² / sqrt(number_of_samples)
This is a standart noise_lvl from measurements in Schillerslage, Germany.
Output in Volt.
"""
return area * noise_lvl * np.sqrt(sample_rate * time) * 1e-9
[docs]def getPhiByGridSearch(data):
"""
"""
npulses, ntimes = data.shape
power = np.zeros([npulses, 360])
sign = np.zeros([npulses, 360])
testphases = np.linspace(-np.pi, np.pi, 361)[:-1]
for n in range(npulses):
for t, tp in enumerate(testphases):
# mind the minus sign!
# remember we search for the phase, not the angle we have to
# apply in order to rotate it back.
test = data[n, :] * (np.cos(-tp) + 1j * np.sin(-tp))
sign[n, t] = np.sign(np.sum(test.real))
power[n, t] = np.sum(test.imag ** 2)
power[sign <= 0] = np.max(power)
phi = testphases[np.argmin(power, axis=1)]
# import matplotlib.pyplot as plt
# fig, ax = plt.subplots(2, 1, sharex=True)
# ax[0].plot(testphases, power[0, :])
return phi
# so to get the best amplitude we have to rotate the data by -phi!
# return phi
# import matplotlib.pyplot as plt
# for i in range(npulses):
# pp = i
# test = data[pp, :] * (np.cos(phi[pp]) + 1j * np.sin(phi[pp])) * 1e9
# ref = data[pp, :] * 1e9
#
# fig, ax = plt.subplots(2)
# ax[0].plot(ref.real, label='data real: {:.3e}'
# .format(np.sum(ref.real ** 2)))
# ax[0].plot(test.real, label='rot imag: {:.3e}'
# .format(np.sum(test.real ** 2)))
# ax[0].set_title('real data')
# plt.legend()
#
# ax[1].plot(ref.imag, label='data imag: {:.3e}'
# .format(np.sum(ref.imag ** 2)))
# ax[1].plot(test.imag, label='rot imag: {:.3e}'
# .format(np.sum(test.imag ** 2)))
# ax[1].set_title('imag data')
# plt.legend()
# return phi
# The End