Modifications for inversion

examples/CMakeLists.txt

@@ -10,6 +10,9 @@ target_link_libraries(  KV0-3loops  "lemmacore" "fdem1d" "merlin")
 add_executable( ModelAligner ModelAligner.cpp ) 
 target_link_libraries(  ModelAligner  "lemmacore" "fdem1d" "merlin")
 
+add_executable( KernelAligner KernelAligner.cpp ) 
+target_link_libraries(  KernelAligner  "lemmacore" "fdem1d" "merlin")
+
 add_executable( ForwardFID ForwardFID.cpp  )
 target_link_libraries(  ForwardFID  "lemmacore" "fdem1d" "merlin")
 

examples/ForwardFID.cpp

@@ -36,6 +36,7 @@ int main(int argc, char** argv) {
         //Forward->SetWindowEdges( VectorXr::LinSpaced(10,0,1) );
         Forward->SetLogSpacedWindows(10,1000,30);
         Forward->SetKernel(Kernel);
+        Forward->SetNoiseFloor(50); // nV
         auto FID = Forward->ForwardModel(Model);
     std::cout << *FID << std::endl;
 

examples/KernelAligner.cpp

@@ -0,0 +1,32 @@
+/* This file is part of Lemma, a geophysical modelling and inversion API.
+ * More information is available at http://lemmasoftware.org
+ */
+
+/* This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
+ */
+
+/**
+ * @file
+ * @date      12/01/2017 10:29:26 PM
+ * @version   $Id$
+ * @author    Trevor Irons (ti)
+ * @email     tirons@egi.utah.edu
+ * @copyright Copyright (c) 2017, University of Utah
+ * @copyright Copyright (c) 2017, Trevor Irons & Lemma Software, LLC
+ */
+
+#include "Merlin"
+
+using namespace Lemma;
+
+int main(int argc, char** argv) {
+    if (argc<2) {
+        std::cout << "useage\n"
+                  << "./KernelAligner  <AkvoData> \n";
+        exit(EXIT_SUCCESS);
+    }
+    auto K0 = KernelV0::NewSP();
+    std::cout << *K0 << std::endl;
+}

examples/ModelAligner.cpp

@@ -34,7 +34,7 @@ int main(int argc, char** argv) {
 
     auto Model = LayeredEarthMR::NewSP();
         Model->AlignWithKernel(Kernel);
-        Model->SetT2StarBins(10, 500, 20);
+        Model->SetT2StarBins(atoi(argv[2]), atoi(argv[3]), atoi(argv[4]));
 
     std::cout << *Model << std::endl;
 }

include/DataFID.h

@@ -129,6 +129,14 @@ namespace Lemma {
         /** Holds the actual data */
         MatrixXcr   FIDData;
 
+        /** Estimate of noise for each datum */
+        MatrixXr   NoiseEstimate;
+
+        /** Forward model matrix */
+        MatrixXcr   Qt;
+
+        VectorXr    TrueModel;
+
         VectorXr    WindowCentres;
 
         VectorXr    WindowEdges;

include/ForwardFID.h

@@ -21,6 +21,7 @@
 #define  FORWARDFID_INC
 
 #pragma once
+#include <random>
 #include "LayeredEarthEM.h"
 #include "PolygonalWireAntenna.h"
 #include "EMEarth1D.h"
@@ -132,6 +133,12 @@ namespace Lemma {
          */
         void SetKernel( std::shared_ptr< KernelV0 > K0 );
 
+        /**
+         *  @param[in] floor is the standard deviation of the noise (zero mean),
+         *             defaults to zero
+         */
+        void SetNoiseFloor( const Real& floor );
+
         // ====================  INQUIRY       =======================
         /**
          *  Returns the name of the underlying class, similiar to Python's type
@@ -148,7 +155,7 @@ namespace Lemma {
         /** Copy is disabled */
         ForwardFID( const ForwardFID& ) = delete;
 
-        // ====================  DATA MEMBERS  =========================
+
         private:
 
         /**
@@ -157,12 +164,16 @@ namespace Lemma {
          */
         void CalcQTMatrix( VectorXr T2StarBins );
 
+        // ====================  DATA MEMBERS  =========================
         /** ASCII string representation of the class name */
         static constexpr auto CName = "ForwardFID";
 
         /** Imaging kernel used in calculation */
         std::shared_ptr< KernelV0 >   Kernel = nullptr;
 
+        /** Noise floor for additive Gaussian noise */
+        Real NoiseFloor = 0.;
+
         /** Time gate windows */
         VectorXr WindowEdges;
 

src/DataFID.cpp

@@ -25,7 +25,7 @@ namespace Lemma {
     // ====================  FRIEND METHODS  =====================
 
     std::ostream &operator << (std::ostream &stream, const DataFID &ob) {
-        stream << ob.Serialize()  << "\n---\n"; // End of doc ---
+        stream <<"%YAML 1.2\n---\n" << ob.Serialize()  << "\n"; // End of doc ---
         return stream;
     }
 
@@ -77,6 +77,12 @@ namespace Lemma {
         // FILL IN CLASS SPECIFICS HERE
         node["FID_REAL"] = FIDData.real().eval();
         node["FID_IMAG"] = FIDData.imag().eval();
+        //node["Qt_REAL"] = Qt.real().eval();
+        //node["Qt_IMAG"] = Qt.imag().eval();
+        MatrixXr QTM = Qt.array().abs();
+        //node["Qt_MOD"] = QTM;
+        //node["TrueModel"] = TrueModel;
+        node["NoiseEstimate"] = NoiseEstimate;
         node["WindowCentres"] = WindowCentres;
         node["WindowEdges"] =  WindowEdges;
         node["PulseMoment"] = PulseMoment;

src/ForwardFID.cpp

@@ -107,6 +107,17 @@ namespace Lemma {
 
     //--------------------------------------------------------------------------------------
     //       Class:  ForwardFID
+    //      Method:  SetNoiseFloor
+    //--------------------------------------------------------------------------------------
+    void ForwardFID::SetNoiseFloor ( const Real& floor ) {
+        assert (floor > 0.);
+        NoiseFloor = floor;
+        return ;
+    }		// -----  end of method ForwardFID::SetNoiseFloor  -----
+
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  ForwardFID
     //      Method:  ForwardModel
     //--------------------------------------------------------------------------------------
     std::shared_ptr<DataFID> ForwardFID::ForwardModel ( std::shared_ptr<LayeredEarthMR> Mod ) {
@@ -120,14 +131,33 @@ namespace Lemma {
         // Forward calculation is just a matrix vector product
         VectorXcr data = QT*Mod->GetModelVector();
 
+        // rearrange solution back into a matrix
+        MatrixXcr B(Eigen::Map<MatrixXcr>(data.data(), nt, nq)); // Complex Data
+        MatrixXr N = MatrixXr::Zero(nt, nq);  // Noise
+        B.real() = B.array().abs();
+        B.imag() *= 0.;
+
         // TODO add noise
+        if (NoiseFloor > 1e-5) {
+            std::random_device rd;
+            std::mt19937 gen(rd());
+            std::normal_distribution<> d(0,1e-9*NoiseFloor);
+            for (int iq=0; iq<nq; ++iq) {
+                for (int it=0; it<nt; ++it) {
+                    B(it, iq) += Complex(d(gen), d(gen)) /
+                        std::sqrt( WindowEdges(it+1)-WindowEdges(it) );
+                    N(it,iq) = (1e-9*NoiseFloor) / std::sqrt(WindowEdges(it+1)-WindowEdges(it));
+                }
+            }
+        }
 
-        // rearrange solution back into a matrix
-        MatrixXcr B(Eigen::Map<MatrixXcr>(data.data(), nt, nq));
         //std::cout << B.imag().transpose() <<std::endl;
         auto FID = DataFID::NewSP();
 
         FID->FIDData = B;
+        FID->NoiseEstimate = N;
+        FID->Qt = QT;
+        FID->TrueModel = Mod->GetModelVector();
         FID->WindowEdges = WindowEdges;
         FID->WindowCentres = WindowCentres;
         FID->PulseMoment = Kernel->GetPulseCurrent()*Kernel->GetPulseDuration();

src/LayeredEarthMR.cpp

@@ -24,7 +24,7 @@ namespace Lemma {
     // ====================  FRIEND METHODS  =====================
 
     std::ostream &operator << (std::ostream &stream, const LayeredEarthMR &ob) {
-        stream << ob.Serialize()  << "\n---\n"; // End of doc ---
+        stream << ob.Serialize()  << "\n"; // End of doc ---
         return stream;
     }
 
@@ -170,7 +170,10 @@ namespace Lemma {
     //      Method:  GetModelVector
     //--------------------------------------------------------------------------------------
     VectorXr LayeredEarthMR::GetModelVector (  ) {
-        VectorXr B(Eigen::Map<VectorXr>(ModelMat.data(), ModelMat.cols()*ModelMat.rows()));
+        //VectorXr B(Eigen::Map<VectorXr>(ModelMat.data(), ModelMat.cols()*ModelMat.rows()));
+        MatrixXr MT = ModelMat.transpose();
+        VectorXr B(Eigen::Map<VectorXr>(MT.data(), MT.cols()*MT.rows()));
+        //VectorXr B(Eigen::Map<VectorXr>(ModelMat.transpose().data(), ModelMat.cols()*ModelMat.rows()));
         return B;
     }		// -----  end of method LayeredEarthMR::GetModelVector  -----