Commit 6d6a575a authored by thomas.forbriger's avatar thomas.forbriger
Browse files

[MERGE] (master|vendor_Seitosh): vendor import of 23.10.2015

Merge branch 'vendor_Seitosh'

Conflicts:
	contrib/aff/README.export
	contrib/fourier/README.export
	contrib/stfinv/README.export

Commands executed on branch master (as a guideline for future imports):

  git merge -X theirs vendor_Seitosh
  git rm contrib/*/README.export
  git commit
parents 02cfa44c 25190a67
......@@ -3,7 +3,6 @@
*
* ----------------------------------------------------------------------------
*
* $Id: stfinvfixedstf.h 3985 2011-05-29 15:43:24Z tforb $
* \author Thomas Forbriger
* \date 06/05/2011
*
......@@ -38,9 +37,7 @@
#ifndef STFINV_STFINVFIXEDSTF_H_VERSION
#define STFINV_STFINVFIXEDSTF_H_VERSION \
"STFINV_STFINVFIXEDSTF_H V1.0 "
#define STFINV_STFINVFIXEDSTF_H_CVSID \
"$Id: stfinvfixedstf.h 3985 2011-05-29 15:43:24Z tforb $"
"STFINV_STFINVFIXEDSTF_H V1.0"
#include<stfinv/stfinvfourier.h>
......
......@@ -3,7 +3,6 @@
*
* ----------------------------------------------------------------------------
*
* $Id: stfinvfourier.cc 4968 2013-02-01 13:58:05Z lrehor $
* \author Thomas Forbriger
* \date 08/05/2011
*
......@@ -32,16 +31,17 @@
* - 08/05/2011 V1.0 Thomas Forbriger
* - 30/09/2011 V1.1 implemented handling of additional time series pairs
* - 04/10/2011 V1.2 correction in debug message
* - 14/10/2015 V1.3 new end-user usage functions
*
* ============================================================================
*/
#define STFINV_STFINVFOURIER_CC_VERSION \
"STFINV_STFINVFOURIER_CC V1.2"
#define STFINV_STFINVFOURIER_CC_CVSID \
"$Id: stfinvfourier.cc 4968 2013-02-01 13:58:05Z lrehor $"
"STFINV_STFINVFOURIER_CC V1.3"
#include <sstream>
#include <stfinv/stfinvfourier.h>
#include <stfinv/stfinvfourier_summary_usage.h>
#include <stfinv/stfinvfourier_description_usage.h>
#include <stfinv/debug.h>
#include <aff/subarray.h>
#include <aff/slice.h>
......@@ -103,6 +103,13 @@ namespace stfinv {
STFFourierDomainEngine::classhelp(os);
} // void STFFourierDomainEngine::help(std::ostream& os) const
/*----------------------------------------------------------------------*/
void STFFourierDomainEngine::usage(std::ostream& os) const
{
STFFourierDomainEngine::classusage(os);
} // void STFFourierDomainEngine::usage(std::ostream& os) const
/*----------------------------------------------------------------------*/
......@@ -110,6 +117,7 @@ namespace stfinv {
{
return("STFFourierDomainEngine");
} // const char const* STFFourierDomainEngine::name() const
/*----------------------------------------------------------------------*/
/*! \brief online help text giving information on options
......@@ -119,23 +127,25 @@ namespace stfinv {
*/
void STFFourierDomainEngine::classhelp(std::ostream& os)
{
os << "Options and parameters in common for Fourier engines:\n"
<< "fpad=f padding factor (default: 1.5)\n"
<< "fpow2 use power of two\n"
<< "fdiv=d use integer multiple of d\n"
<< "tshift=d delay source correction filter wavelet by \"d\"\n"
<< " seconds in order to expose acausal parts\n"
<< "These options define the number of samples N used for the FFT.\n"
<< "This should be larger than the number of samples M in the\n"
<< "original series to avoid wrap-around. N=M*f at least. If fpow2\n"
<< "is set, N will be the next power of 2 larger than M*f. Else if\n"
<< "fdiv is set, N will be the next integer multiple of d larger\n"
<< "than M*f." << std::endl;
Tbase::classhelp(os);
os << stfinvfourier_summary_usage;
} // void STFFourierDomainEngine::classhelp(std::ostream& os)
/*----------------------------------------------------------------------*/
/*! \brief online help text giving information on options
*
* This must be kept synchronous with the options used by
* STFFourierDomainEngine::initialize()
*/
void STFFourierDomainEngine::classusage(std::ostream& os)
{
os << stfinvfourier_description_usage;
os << std::endl;
Tbase::classusage(os);
} // void STFFourierDomainEngine::classusage(std::ostream& os)
/*----------------------------------------------------------------------*/
/*! \brief Create FFT engines
*
* \par Number of samples
......
......@@ -3,7 +3,6 @@
*
* ----------------------------------------------------------------------------
*
* $Id: stfinvfourier.h 4968 2013-02-01 13:58:05Z lrehor $
* \author Thomas Forbriger
* \date 08/05/2011
*
......@@ -30,6 +29,7 @@
* REVISIONS and CHANGES
* - 08/05/2011 V1.0 Thomas Forbriger
* - 30/09/2011 V1.1 implemented handling of additional time series pairs
* - 14/10/2015 V1.2 new end-user usage functions
*
* ============================================================================
*/
......@@ -38,9 +38,7 @@
#ifndef STFINV_STFINVFOURIER_H_VERSION
#define STFINV_STFINVFOURIER_H_VERSION \
"STFINV_STFINVFOURIER_H V1.1"
#define STFINV_STFINVFOURIER_H_CVSID \
"$Id: stfinvfourier.h 4968 2013-02-01 13:58:05Z lrehor $"
"STFINV_STFINVFOURIER_H V1.2"
#include <stfinv/stfinvbase.h>
#include <aff/array.h>
......@@ -114,6 +112,10 @@ namespace stfinv {
virtual void help(std::ostream& os=std::cout) const;
//! \brief print online help
static void classhelp(std::ostream& os=std::cout);
//! \brief print detailed description
virtual void usage(std::ostream& os=std::cout) const;
//! \brief print detailed description
static void classusage(std::ostream& os=std::cout);
//! \brief return name of engine
virtual const char* name() const;
protected:
......
......@@ -3,7 +3,6 @@
*
* ----------------------------------------------------------------------------
*
* $Id: stfinvidentity.cc 3986 2011-05-29 16:01:09Z tforb $
* \author Thomas Forbriger
* \date 07/05/2011
*
......@@ -29,23 +28,31 @@
*
*
* REVISIONS and CHANGES
* - 07/05/2011 V1.0 Thomas Forbriger
* - 07/05/2011 V1.0 Thomas Forbriger (thof)
* - 21/02/2014 V1.1 implemented scaling to average weighted energy
* - 14/10/2015 V1.2 new end-user usage functions
* - 19/10/2015 V1.3 weight factor is defined in terms amplitude; apply
* square of weight to signal energy
*
* ============================================================================
*/
#define STFINV_STFINVIDENTITY_CC_VERSION \
"STFINV_STFINVIDENTITY_CC V1.0 "
#define STFINV_STFINVIDENTITY_CC_CVSID \
"$Id: stfinvidentity.cc 3986 2011-05-29 16:01:09Z tforb $"
"STFINV_STFINVIDENTITY_CC V1.3"
#include <stfinv/stfinvidentity.h>
#include <stfinv/stfinvidentity_summary_usage.h>
#include <stfinv/stfinvidentity_description_usage.h>
#include <aff/functions/sqrsum.h>
#include <aff/seriesoperators.h>
#include <cmath>
#include <stfinv/tools.h>
namespace stfinv {
const char* const STFEngineIdentity::ID="ident";
const char* const STFEngineIdentity::description
="identity: just apply a scalar factor";
="scale with amplitude factor";
/*----------------------------------------------------------------------*/
......@@ -56,6 +63,13 @@ namespace stfinv {
/*----------------------------------------------------------------------*/
void STFEngineIdentity::usage(std::ostream& os) const
{
STFEngineIdentity::classusage(os);
} // void STFEngineIdentity::usage(std::ostream& os) const
/*----------------------------------------------------------------------*/
const char* STFEngineIdentity::name() const
{
return("STFEngineIdentity");
......@@ -73,16 +87,42 @@ namespace stfinv {
void STFEngineIdentity::exec()
{
STFINV_assert(!Mscaleenergy,
"ERROR: energy scaling not yet implemented");
// effective amplitude factor
double fac=1.;
// scale to reproduce average energy if requested
if (Mscaleenergy)
{
double recording_sqrsum=0.;;
double synthetic_sqrsum=0.;;
for (unsigned int i=0; i<this->nreceivers(); ++i)
{
synthetic_sqrsum
+= aff::func::sqrsum(this->synthetic(i))
* this->weight(i) * this->weight(i);
recording_sqrsum
+= aff::func::sqrsum(this->recording(i))
* this->weight(i) * this->weight(i);
fac = std::sqrt(recording_sqrsum/synthetic_sqrsum);
}
} // if (Mscaleenergy)
Tseries stf=this->stf();
stf=0.;
stf(0)=1.;
stf(0)=fac/this->dt();
for (unsigned int i=0; i<this->nreceivers(); ++i)
{
Tseries::Tcoc synthetic=this->synthetic(i);
Tseries convolvedsynthetic=this->convolvedsynthetic(i);
convolvedsynthetic.copyin(synthetic);
if (Mscaleenergy) { convolvedsynthetic *= fac; }
}
for (unsigned int i=0; i<this->npairs(); ++i)
{
Tseries::Tcoc series=this->series(i);
Tseries convolvedseries=this->convolvedseries(i);
convolvedseries.copyin(series);
if (Mscaleenergy) { convolvedseries *= fac; }
}
} // void STFEngineIdentity::exec()
......@@ -90,19 +130,22 @@ namespace stfinv {
void STFEngineIdentity::classhelp(std::ostream& os)
{
os << "class STFEngineIdentity ("
<< STFEngineIdentity::ID << ")\n";
os << STFEngineIdentity::description << "\n" << std::endl;
os << "This engine convolves the synthetic data with a discrete delta\n"
<< "pulse so to speak. Optionally the delta-peak ist scale such that\n"
<< "the convolved synthetics will be of equal scaled energy as the\n"
<< "recordings.\n";
os << "Options and parameters in common for Fourier engines:\n"
<< "scaleenergy if flag is set: scale energy"
<< std::endl;
Tbase::classhelp(os);
os << stfinvidentity_summary_usage;
os << std::endl;
stfinv::tools::report_engine_ID<STFEngineIdentity>(os);
} // void STFEngineIdentity::classhelp(std::ostream& os)
/*----------------------------------------------------------------------*/
void STFEngineIdentity::classusage(std::ostream& os)
{
os << stfinvidentity_description_usage;
os << std::endl;
Tbase::classusage(os);
os << std::endl;
stfinv::tools::report_engine_ID<STFEngineIdentity>(os);
} // void STFEngineIdentity::classusage(std::ostream& os)
} // namespace stfinv
/* ----- END OF stfinvidentity.cc ----- */
......@@ -3,7 +3,6 @@
*
* ----------------------------------------------------------------------------
*
* $Id: stfinvidentity.h 3986 2011-05-29 16:01:09Z tforb $
* \author Thomas Forbriger
* \date 07/05/2011
*
......@@ -30,6 +29,7 @@
*
* REVISIONS and CHANGES
* - 07/05/2011 V1.0 Thomas Forbriger
* - 14/10/2015 V1.1 new end-user usage functions
*
* ============================================================================
*/
......@@ -38,9 +38,7 @@
#ifndef STFINV_STFINVIDENTITY_H_VERSION
#define STFINV_STFINVIDENTITY_H_VERSION \
"STFINV_STFINVIDENTITY_H V1.0 "
#define STFINV_STFINVIDENTITY_H_CVSID \
"$Id: stfinvidentity.h 3986 2011-05-29 16:01:09Z tforb $"
"STFINV_STFINVIDENTITY_H V1.1"
#include<stfinv/stfinvbase.h>
......@@ -71,6 +69,15 @@ namespace stfinv {
: Tbase(triples, stf, parameters),
Mscaleenergy(false)
{ this->initialize(); }
/*! \brief Constructor.
*/
STFEngineIdentity(const stfinv::Tvectoroftriples& triples,
const stfinv::Waveform& stf,
const stfinv::Tvectorofpairs& pairs,
const std::string& parameters)
: Tbase(triples, stf, pairs, parameters),
Mscaleenergy(false)
{ this->initialize(); }
//! \brief abstract base requires virtual destructor
virtual ~STFEngineIdentity() { }
//! \brief Start engine
......@@ -79,6 +86,10 @@ namespace stfinv {
virtual void help(std::ostream& os=std::cout) const;
//! \brief print online help
static void classhelp(std::ostream& os=std::cout);
//! \brief print detailed description
virtual void usage(std::ostream& os=std::cout) const;
//! \brief print detailed description
static void classusage(std::ostream& os=std::cout);
//! \brief return name of engine
virtual const char* name() const;
private:
......
......@@ -3,7 +3,6 @@
*
* ----------------------------------------------------------------------------
*
* $Id: stfinvnormalize.cc 4164 2011-10-04 09:00:57Z tforb $
* \author Thomas Forbriger
* \date 08/05/2011
*
......@@ -35,10 +34,11 @@
*/
#define STFINV_STFINVNORMALIZE_CC_VERSION \
"STFINV_STFINVNORMALIZE_CC V1.0"
#define STFINV_STFINVNORMALIZE_CC_CVSID \
"$Id: stfinvnormalize.cc 4164 2011-10-04 09:00:57Z tforb $"
#include <stfinv/stfinvnormalize.h>
#include <aff/functions/sqrsum.h>
#include <stfinv/debug.h>
#include <cmath>
namespace stfinv {
......@@ -65,14 +65,70 @@ namespace stfinv {
void STFEngineNormalize::initialize()
{
STFINV_illegal;
// scale energy
std::istringstream is (this->parameter("waterlevel","1.e-3"));
is >> Mwaterlevel;
STFINV_debug(Mdebug&1, "STFEngineNormalize::initialize()",
"Mwaterlevel=" << Mwaterlevel);
STFINV_assert(Mwaterlevel > 0,
"ERROR: parameter for option \"waterlevel\" not larger than 0");
} // void STFEngineNormalize::initialize()
/*----------------------------------------------------------------------*/
void STFEngineNormalize::exec()
{
STFINV_illegal;
// read signals and calculate FFT
this->fftinput();
// cycle through frequencies
// calculate weighted energy of all synthetic traces
Tseries syntheticenergy(0,this->nfreq()-1);
syntheticenergy=0.;
double totalsyntheticenergy=0.;
for (unsigned int i=0; i<this->nfreq(); ++i)
{
TAspectrum syntheticref=this->syntheticcoeff(i);
for (unsigned int j=0; j<this->nreceivers(); ++j)
{
syntheticenergy(i) += this->weight(j)*this->weight(j)
*std::norm(syntheticref(j));
}
totalsyntheticenergy += syntheticenergy(i);
}
// calculate waterlevel
double waterlevel=Mwaterlevel*totalsyntheticenergy/this->nfreq();
double weightsum=aff::func::sqrsum(this->weights());
for (unsigned int i=0; i<this->nfreq(); ++i)
{
double recordingenergy=0.;
TAspectrum recordingref=this->recordingcoeff(i);
TAspectrum syntheticref=this->syntheticcoeff(i);
for (unsigned int j=0; j<this->nreceivers(); ++j)
{
recordingenergy += this->weight(j)*this->weight(j)
*std::norm(recordingref(j));
}
double modulus
=std::sqrt(recordingenergy/(syntheticenergy(i)+waterlevel));
double phase = 0.;
for (unsigned int j=0; j<this->nreceivers(); ++j)
{
phase += this->weight(j)*
(std::arg(recordingref(j))-std::arg(syntheticref(j)));
}
phase /= weightsum;
this->stfcoeff(i)=std::polar(modulus, phase);
} // for (unsigned int i=0; i<this->nfreq(); ++i)
// provide results to user
this->fftoutput();
} // void STFEngineNormalize::exec()
/*----------------------------------------------------------------------*/
......@@ -82,8 +138,13 @@ namespace stfinv {
os << "class STFEngineNormalize ("
<< STFEngineNormalize::ID << ")\n";
os << STFEngineNormalize::description << "\n" << std::endl;
os << "Options and parameters in common for Fourier engines:\n"
<< "NOT YET IMPLEMENTED"
STFINV_illegal;
os << "DESCRIBE HERE\n"
<< "A waterlevel as a fraction of the signal energy of the\n"
<< "input synthetics is applied. If per receiver scaling is\n"
<< "selected, the receivers will be weighted in the deconvolution.\n";
os << "Options and parameters:\n"
<< "waterlevel=l waterlevel to be applied for regularization."
<< std::endl;
Tbase::classhelp(os);
} // void STFEngineNormalize::classhelp(std::ostream& os)
......
......@@ -3,7 +3,6 @@
*
* ----------------------------------------------------------------------------
*
* $Id: stfinvnormalize.h 4968 2013-02-01 13:58:05Z lrehor $
* \author Thomas Forbriger
* \date 08/05/2011
*
......@@ -40,8 +39,6 @@
#define STFINV_STFINVNORMALIZE_H_VERSION \
"STFINV_STFINVNORMALIZE_H V1.1"
#define STFINV_STFINVNORMALIZE_H_CVSID \
"$Id: stfinvnormalize.h 4968 2013-02-01 13:58:05Z lrehor $"
#include<stfinv/stfinvfourier.h>
......@@ -70,6 +67,7 @@ namespace stfinv {
* \f$f_l\f$ and receiver \f$k\f$ at offset \f$r_k\f$ and
* - \f$s_{lk}\f$ is the Fourier coefficient of the corresponding
* synthetics.
*
* Then we seek a source correction filter with Fourier coefficients
* \f$q_l\f$ such that
* \f[
......@@ -85,6 +83,45 @@ namespace stfinv {
* \sum\limits_{k}\Im(\ln(s_{lk}\,q_l))
* \quad\forall\, l.
* \f]
*
* Actually we will apply weights \f$w_k\f$ to data from different offsets
* \f$r_k\f$ such that
* \f[
* \sum\limits_{k}\left|w_k\,d_{lk}\right|^2
* =
* \sum\limits_{k}\left|w_k\,s_{lk}\,q_l\right|^2
* \quad\forall\, l
* \f]
* and
* \f[
* \sum\limits_{k}w_k\,\Im(\ln(d_{lk}))
* =
* \sum\limits_{k}w_k\,\Im(\ln(s_{lk}\,q_l))
* \quad\forall\, l.
* \f]
* Additionally we will apply a waterlevel \f$\epsilon^2\f$, such that
* Fourier coefficients for which the average weighted energy of the
* synthetics is smaller than a fraction \f$\epsilon^2\f$ of the total
* average energy of the synthetics, will be damped.
*
* These conditions are satisfied by
* \f[
* q_l=A_l\,e^{i\Phi_{l}}
* \f]
* with
* \f[
* A_l=\sqrt{\frac{\sum\limits_{k}\left|w_k\,d_{lk}\right|^2}{
* \sum\limits_{k}\left|w_k\,s_{lk}\right|^{2}
* +
* \frac{\epsilon^2}{N_{f}}\,\sum\limits_{k}\sum\limits_{l=0}^{N_{f}-1}
* \left|w_k\,s_{lk}\right|^2}}
* \f]
* and
* \f[
* \Phi_{l}=\frac{\sum\limits_{k}w_{k}\left(
* \Im\left(\ln(d_{lk})\right)-\Im\left(\ln(s_{lk})\right)
* \right)}{\sum\limits_k w_{k}}
* \f]
*/
class STFEngineNormalize: public stfinv::STFFourierDomainEngine {
public:
......@@ -117,6 +154,8 @@ namespace stfinv {
// member data
private:
//! \brief waterlevel
double Mwaterlevel;
}; // class STFEngineNormalize
} // namespace stfinv
......
......@@ -67,9 +67,5 @@ Information concerning copyright and license:
along with this program; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
P.S.: Please have a look at trunk/src/ts/wf/testcases/Makefile.soutifu to find
P.S.: Please have a look at src/ts/wf/testcases/Makefile.soutifu to find
some rules to create test data
# Makefile for Test_libstfinv
#
# Copyright (c) 2011 by Thomas Bohlen (KIT Karlsruhe) and Lisa Rehor (KIT Karlsruhe)
#
# This file was copied from the DENISE code and adjusted for its use
# in this test programme.
#
# ----
# DENISE 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 2 of the License, or
# (at your option) any later version.
#
# DENISE 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, write to the Free Software
# Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
# ----
#--------------------------------------------------------
# edit here:
......
/*------------------------------------------------------------------------
* Test_libstfinv.c
*
* Copyright (c) 2011 by Lisa Rehor (KIT Karlsruhe)
* Copyright (c) 2011 by Lisa Rehor (KIT Karlsruhe) and Martin Schaefer (KIT Karlsruhe)
*
* ----
* 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 2 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, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
* ---------------------------------------------------------------------*/
......@@ -171,7 +186,7 @@ for(tracl1=0;tracl1<ntr;tracl1++){
/* storing the traces itself in a vector */
for(j=0;j<ns;j++){
dump=tr_add[1].data[j];
dump=tr_add[0].data[j];
add[j]=dump;
}
......
/*------------------------------------------------------------------------
* fd.h - include file for FD programs
*
* Copyright (c) 2011 by Thomas Bohlen (KIT Karlsruhe) and Lisa Rehor (KIT Karlsruhe)
*
* This file was copied from the DENISE code and adjusted for its use
* in this test programme.
* For copyright information see the DENISE code.
*
* ----
* DENISE 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 2 of the License, or
* (at your option) any later version.
*
* DENISE 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, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*
* ---------------------------------------------------------------------*/
/* files to include */
......