Forward modelling of FID data is possible, but final format of files not nailed down.

examples/CMakeLists.txt

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

examples/ForwardFID.cpp

@@ -22,24 +22,22 @@ using namespace Lemma;
 
 int main(int argc, char** argv) {
 
-    if (argc<5) {
-        std::cout << "./ForwardFID Kernel.yaml TxString RxString  vtkoutput<true/false> \n";
-    //    return(EXIT_FAILURE);
+    if (argc<3) {
+        std::cout << "./ForwardFID Kernel.yaml Model.yaml  \n";
+        return(EXIT_FAILURE);
     }
     auto Kernel = KernelV0::DeSerialize(YAML::LoadFile(argv[1]));
     //    std::cout << *Kernel;
 
-    auto Model = LayeredEarthMR::NewSP();
-        Model->AlignWithKernel(Kernel);
-        Model->SetT2StarBins(10, 500, 20);
-        std::cout << *Model << std::endl;
+    auto Model = LayeredEarthMR::DeSerialize(YAML::LoadFile(argv[2]));
+    //    std::cout << *Model << std::endl;
 
-/*
     auto Forward = ForwardFID::NewSP();
         //Forward->SetWindowEdges( VectorXr::LinSpaced(10,0,1) );
         Forward->SetLogSpacedWindows(10,1000,30);
-        //auto FID =  Forward->ForwardModel(Model);
-    //std::cout << *Forward << std::endl;
-*/
+        Forward->SetKernel(Kernel);
+        auto FID = Forward->ForwardModel(Model);
+    std::cout << *FID << std::endl;
+
     return EXIT_SUCCESS;
 }

examples/ModelAligner.cpp

@@ -0,0 +1,42 @@
+/* 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      08/30/2017 04:08:53 AM
+ * @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<5) {
+        std::cout << "ModelAligner aligns a dummy model with a pre-calculated"
+                  << "imaging kernel.\n\n"
+                  << "./ModelAligner Kernel.yaml T2Low T2High nT2  \n";
+
+        return(EXIT_FAILURE);
+    }
+
+    auto Kernel = KernelV0::DeSerialize(YAML::LoadFile(argv[1]));
+
+    auto Model = LayeredEarthMR::NewSP();
+        Model->AlignWithKernel(Kernel);
+        Model->SetT2StarBins(10, 500, 20);
+
+    std::cout << *Model << std::endl;
+}
+
+

include/DataFID.h

@@ -36,6 +36,8 @@ namespace Lemma {
 
         friend std::ostream &operator<<(std::ostream &stream, const DataFID &ob);
 
+        friend class ForwardFID;
+
         protected:
         /*
          *  This key is used to lock the constructor. It is protected so that inhereted
@@ -124,6 +126,15 @@ namespace Lemma {
         /** ASCII string representation of the class name */
         static constexpr auto CName = "DataFID";
 
+        /** Holds the actual data */
+        MatrixXcr   FIDData;
+
+        VectorXr    WindowCentres;
+
+        VectorXr    WindowEdges;
+
+        VectorXr    PulseMoment;
+
     }; // -----  end of class  DataFID  -----
 }  // -----  end of namespace Lemma ----
 

include/ForwardFID.h

@@ -24,9 +24,11 @@
 #include "LayeredEarthEM.h"
 #include "PolygonalWireAntenna.h"
 #include "EMEarth1D.h"
+
 #include "MerlinObject.h"
 #include "DataFID.h"
 #include "KernelV0.h"
+#include "LayeredEarthMR.h"
 
 namespace Lemma {
 
@@ -107,7 +109,7 @@ namespace Lemma {
          *  Performs forward model calculation based on input parameters
          *  @return Merlin class representing data
          */
-        std::shared_ptr<DataFID> ForwardModel();
+        std::shared_ptr<DataFID> ForwardModel( std::shared_ptr<LayeredEarthMR> );
 
         // ====================  ACCESS        =======================
         /**
@@ -125,6 +127,11 @@ namespace Lemma {
          */
         void SetLogSpacedWindows(const Real& first, const Real& last, const int& n);
 
+        /**
+         *   @param[in] K0 is the initial amplitude imaging kernel
+         */
+        void SetKernel( std::shared_ptr< KernelV0 > K0 );
+
         // ====================  INQUIRY       =======================
         /**
          *  Returns the name of the underlying class, similiar to Python's type
@@ -144,7 +151,11 @@ namespace Lemma {
         // ====================  DATA MEMBERS  =========================
         private:
 
-        void CalcQTMatrix(  );
+        /**
+         * Calculates the QT matrix
+         * @param[in] T2StarBins are the T2* bins to use
+         */
+        void CalcQTMatrix( VectorXr T2StarBins );
 
         /** ASCII string representation of the class name */
         static constexpr auto CName = "ForwardFID";
@@ -158,6 +169,9 @@ namespace Lemma {
         /** Time gate centres */
         VectorXr WindowCentres;
 
+        /** QT matrix */
+        MatrixXcr QT;
+
         /** Include RDP effects? */
         bool RDP = false;
 

include/LayeredEarthMR.h

@@ -102,7 +102,20 @@ namespace Lemma {
 
         // ====================  ACCESS        =======================
 
+        /**
+         * @return a VectorXr of the T2* bins
+         */
+        VectorXr GetT2StarBins();
 
+        /**
+         * @return a flattened version of the model matrix
+         */
+        VectorXr GetModelVector();
+
+        /**
+         * @return the model matrix
+         */
+        MatrixXr GetModelMatrix();
 
         /**
          *  Sets the T2StarBins to solve for, these are log spaced

src/DataFID.cpp

@@ -75,6 +75,11 @@ namespace Lemma {
         YAML::Node node = MerlinObject::Serialize();
         node.SetTag( GetName() );
         // FILL IN CLASS SPECIFICS HERE
+        node["FID_REAL"] = FIDData.real().eval();
+        node["FID_IMAG"] = FIDData.imag().eval();
+        node["WindowCentres"] = WindowCentres;
+        node["WindowEdges"] =  WindowEdges;
+        node["PulseMoment"] = PulseMoment;
         return node;
     }		// -----  end of method DataFID::Serialize  -----
 

src/ForwardFID.cpp

@@ -109,8 +109,28 @@ namespace Lemma {
     //       Class:  ForwardFID
     //      Method:  ForwardModel
     //--------------------------------------------------------------------------------------
-    std::shared_ptr<DataFID> ForwardFID::ForwardModel (  ) {
+    std::shared_ptr<DataFID> ForwardFID::ForwardModel ( std::shared_ptr<LayeredEarthMR> Mod ) {
+
+        MatrixXcr K0 = Kernel->GetKernel();
+        int nq = K0.cols();
+        int nt = WindowCentres.size();
+
+        CalcQTMatrix( Mod->GetT2StarBins() );
+
+        // Forward calculation is just a matrix vector product
+        VectorXcr data = QT*Mod->GetModelVector();
+
+        // TODO add noise
+
+        // 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->WindowEdges = WindowEdges;
+        FID->WindowCentres = WindowCentres;
+        FID->PulseMoment = Kernel->GetPulseCurrent()*Kernel->GetPulseDuration();
         return FID;
     }		// -----  end of method ForwardFID::ForwardModel  -----
 
@@ -131,11 +151,56 @@ namespace Lemma {
         return;
     }		// -----  end of method ForwardFID::LogSpaced  -----
 
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  ForwardFID
+    //      Method:  SetKernel
+    //--------------------------------------------------------------------------------------
+    void ForwardFID::SetKernel ( std::shared_ptr< KernelV0 > K0 ) {
+        Kernel = K0;
+        return ;
+    }		// -----  end of method ForwardFID::SetKernel  -----
+
+
     //--------------------------------------------------------------------------------------
     //       Class:  ForwardFID
     //      Method:  CalcQTMatrix
     //--------------------------------------------------------------------------------------
-    void ForwardFID::CalcQTMatrix (  ) {
+    void ForwardFID::CalcQTMatrix ( VectorXr T2Bins ) {
+        MatrixXcr K0 = Kernel->GetKernel();
+        VectorXcr K0r(Eigen::Map<VectorXcr>(K0.data(), K0.cols()*K0.rows()));
+
+        // K0 = nq     \times nlay
+        // Qt = nq*nt  \times nlay*nT2
+
+        int nLay = K0.rows();
+        int nq = K0.cols();
+        int nt = WindowCentres.size();
+        int nT2 = T2Bins.size();
+        //std::cout << "# nLay " << nLay << std::endl;
+        //std::cout << "# nq " << nq << std::endl;
+        //std::cout << "# nt " << nt << std::endl;
+        //std::cout << "# nT2 " << nT2 << std::endl;
+
+        QT = MatrixXcr::Zero( nq*nt, nLay*nT2 );
+//        std::cout << "K0 " << K0.rows() << "\t" << K0.cols() << std::endl;
+//        std::cout << "QT " << QT.rows() << "\t" << QT.cols() << std::endl;
+
+        // Ugly!
+        int ir=0;
+        for (int iq=0; iq<nq; ++iq) {
+            for (int it=0; it<nt; ++it) {
+                int ic=0;
+                for (int ilay=0; ilay<nLay; ++ilay) {
+                    for (int it2=0; it2<nT2; ++it2) {
+                        QT(ir, ic) = K0(ilay, iq) * std::exp( -WindowCentres[it]/T2Bins[it2] );
+                        ++ic;
+                    }
+                }
+                ++ir;
+            }
+        }
+//        std::cout << QT.imag() << std::endl;
         return ;
     }		// -----  end of method ForwardFID::CalcQTMatrix  -----
 

src/KernelV0.cpp

@@ -85,6 +85,12 @@ namespace Lemma {
             }
         }
 
+        if (node["K0"]) {
+            Kern = MatrixXcr::Zero( Interfaces.size()-1, PulseI.size()  ).array() + 1.;
+            for ( int ilay=0; ilay<Interfaces.size()-1; ++ilay ) {
+                Kern.row(ilay) = node["K0"]["layer-" + to_string(ilay) ].as<VectorXcr>();
+            }
+        }
     }  // -----  end of method KernelV0::KernelV0  (constructor)  -----
 
     //--------------------------------------------------------------------------------------
@@ -137,6 +143,7 @@ namespace Lemma {
         node["IntegrationSize"] = Size;
         node["IntegrationOrigin"] = Origin;
 
+        // TODO, use better matrix encapulation
         if (Kern.array().abs().any() > 1e-16) {
             for ( int ilay=0; ilay<Interfaces.size()-1; ++ilay ) {
                 node["K0"]["layer-" + to_string(ilay) ] = static_cast<VectorXcr>(Kern.row(ilay));

src/LayeredEarthMR.cpp

@@ -146,6 +146,35 @@ namespace Lemma {
         return ;
     }		// -----  end of method LayeredEarthMR::InitModelMat  -----
 
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  LayeredEarthMR
+    //      Method:  GetT2StarBins
+    //--------------------------------------------------------------------------------------
+    VectorXr LayeredEarthMR::GetT2StarBins (  ) {
+        return T2StarBins ;
+    }		// -----  end of method LayeredEarthMR::GetT2StarBins  -----
+
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  LayeredEarthMR
+    //      Method:  GetModelMatrix
+    //--------------------------------------------------------------------------------------
+    MatrixXr LayeredEarthMR::GetModelMatrix (  ) {
+        return ModelMat ;
+    }		// -----  end of method LayeredEarthMR::GetModelMatrix  -----
+
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  LayeredEarthMR
+    //      Method:  GetModelVector
+    //--------------------------------------------------------------------------------------
+    VectorXr LayeredEarthMR::GetModelVector (  ) {
+        VectorXr B(Eigen::Map<VectorXr>(ModelMat.data(), ModelMat.cols()*ModelMat.rows()));
+        return B;
+    }		// -----  end of method LayeredEarthMR::GetModelVector  -----
+
+
 } // ----  end of namespace Lemma  ----
 
 /* vim: set tabstop=4 expandtab: */