Skip to content
Snippets Groups Projects
Commit da54245a authored by Boris Martin's avatar Boris Martin
Browse files

routine for computing the fitting coefficients and testing it on simple cases

parent 1d2fee32
No related branches found
No related tags found
1 merge request!10Draft: Source estimation
......@@ -51,6 +51,12 @@ public:
gmshFem.userDefinedParameter(remesh, "remesh");
_remesh = ((bool) remesh);
};
ConfigurationInterface() : _name("dummy"), _parametrization(nullptr), _mc(1,0.)
{
_remesh = true;
};
virtual ~ConfigurationInterface() {delete _parametrization;};
virtual void mesh() const;
......
......@@ -182,4 +182,39 @@ Data<T_Physic> arg(const Data<T_Physic>& d)
return out;
};
/**
* @brief Returns a vector ([f][s]) that tells for each shot what coefficient minimizes the L² misfit for that shot.
* If d0 is the computed value for a unit shot and d the reference data, we look for a complexe number "alpha" such that
* alpha * d0 is the best approximation to d.
* It is given by alpha = <d0, d> / <d0, d0> with < , > the complex inner product
*/
template<Physic T_Physic>
std::vector<std::vector<std::complex<double>>> leastSquareShotIntensity(const Data<T_Physic>& d0, const Data<T_Physic>& d)
{
areCompatible(d0, d);
std::vector<std::vector<std::complex<double>>> alpha;
alpha.resize(d.nf());
for (auto &shots: alpha)
shots.resize(d.ns());
for (unsigned int f = 0; f < d.nf(); f++)
{
for (unsigned int s = 0; s < d.ns(); s++)
{
std::complex<double> num = 0.;
std::complex<double> den = 0.;
for (unsigned int r = 0; r < d.nr(s); r++)
{
num += conj<T_Physic>(d0.value(f,s,r)) * d.value(f,s,r);
den += conj<T_Physic>(d0.value(f,s,r)) * d0.value(f,s,r);
}
alpha[f][s] = num / den;
}
}
return alpha;
};
#endif // H_COMMON_DATA_ELEMENT
......@@ -7,6 +7,10 @@
#include <gmshfem/Formulation.h>
#include <gmshfem/GmshFem.h>
#include "../common/wave/element.h"
#include "../common/data/element.h"
#include "../specific/data/element.h"
#include "dummyAcquisition.h"
#define INIT_GMSHFEM() \
namespace geo = gmsh::model::geo; \
......@@ -102,4 +106,26 @@ BOOST_AUTO_TEST_CASE(dummy_test)
BOOST_CHECK(1 + 1 == 2);
}
BOOST_AUTO_TEST_CASE(LeastSquareShotIntensityFitAcoustic)
{
std::complex<double> im(0, 1);
// A Data needs a configuration with some shots
DummyConfigurationOneSource config(2, {1, 2}); // Check perfect fit if one rec, least square fit for 2
std::vector<double> freqs{0};
Data<Physic::acoustic> target(freqs, &config);
Data<Physic::acoustic> unitShot(freqs, &config);
target.value(0, 0, 0, im);
unitShot.value(0, 0, 0, 1.); // Expected coef is "i"
target.value(0, 1, 0, 1. + im);
target.value(0, 1, 1, 1. - im);
unitShot.value(0, 1, 0, -2);
unitShot.value(0, 1, 1, -2); // Expected coef is -0.5
auto alphas = leastSquareShotIntensity(unitShot, target);
BOOST_CHECK_SMALL(std::norm(alphas[0][0] - im), 1e-4);
BOOST_CHECK_SMALL(std::norm(alphas[0][1] - (-0.5)), 1e-4);
}
BOOST_AUTO_TEST_SUITE_END()
\ No newline at end of file
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment