Main Lemma Repository
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

inversesolvertem1d.cpp 11KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419
  1. /* This file is part of Lemma, a geophysical modelling and inversion API */
  2. /* This Source Code Form is subject to the terms of the Mozilla Public
  3. * License, v. 2.0. If a copy of the MPL was not distributed with this
  4. * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
  5. /**
  6. @file
  7. @author M. Andy Kass
  8. @date 03/23/2011
  9. @version $Id: inversesolvertem1d.cpp 193 2014-11-10 23:51:41Z tirons $
  10. **/
  11. #include "inversesolvertem1d.h"
  12. namespace Lemma {
  13. std::ostream &operator<<(std::ostream &stream,
  14. const InverseSolverTEM1D &ob) {
  15. stream << *(LemmaObject*)(&ob);
  16. stream << "Wire antenna source address: " << ob.Antenna << std::endl;
  17. stream << "Dipole source address: " << ob.Dipole << std::endl;
  18. stream << *ob.Receiver << std::endl;
  19. stream << *ob.StartingMod << std::endl;
  20. stream << *ob.RecoveredMod << std::endl;
  21. stream << *ob.InstrumentFor << std::endl;
  22. return stream;
  23. }
  24. // ==================== LIFECYCLE =======================
  25. InverseSolverTEM1D::InverseSolverTEM1D(const std::string &name) :
  26. InverseSolver(name),Antenna(NULL), Dipole(NULL),
  27. Receiver(NULL), StartingMod(NULL), RecoveredMod(NULL),
  28. InstrumentFor(NULL) {
  29. this->InstrumentFor = InstrumentTem::New();
  30. //this->SCondIndex.resize(42);
  31. //this->SThickIndex.resize(42);
  32. //this->SCondIndex(0) = 42;
  33. //this->SThickIndex(0) = 42;
  34. }
  35. InverseSolverTEM1D::~InverseSolverTEM1D() {
  36. if (NumberOfReferences !=0) throw DeleteObjectWithReferences(this);
  37. if (this->Antenna != NULL) this->Antenna->DetachFrom(this);
  38. if (this->Dipole != NULL) this->Dipole->DetachFrom(this);
  39. if (this->Receiver != NULL) this->Receiver->DetachFrom(this);
  40. if (this->StartingMod != NULL) this->StartingMod->DetachFrom(this);
  41. if (this->RecoveredMod != NULL) this->RecoveredMod->DetachFrom(this);
  42. this->InstrumentFor->Delete();
  43. }
  44. InverseSolverTEM1D* InverseSolverTEM1D::New() {
  45. InverseSolverTEM1D* Obj =
  46. new InverseSolverTEM1D( "InverseSolverTEM1D" );
  47. Obj->AttachTo(Obj);
  48. return Obj;
  49. }
  50. void InverseSolverTEM1D::Delete() {
  51. this->DetachFrom(this);
  52. }
  53. void InverseSolverTEM1D::Release() {
  54. delete this;
  55. }
  56. // ==================== OPERATIONS =======================
  57. void InverseSolverTEM1D::Calculate() {
  58. MatrixXr jac;
  59. VectorXr newthickness;
  60. VectorXcr newconductivity;
  61. VectorXr oldthickness;
  62. VectorXcr oldconductivity;
  63. MatrixXr olddata;
  64. MatrixXr newdata;
  65. VectorXr deltathick;
  66. VectorXcr deltacond;
  67. MatrixXr H;
  68. VectorXr g;
  69. VectorXr p;
  70. VectorXr s;
  71. VectorXr y;
  72. //int nrec=42;
  73. int nlayers=42;
  74. this->InstrumentFor->EMEarthModel(this->RecoveredMod);
  75. this->InstrumentFor->SetTransmitLoop(this->Antenna);
  76. this->InstrumentFor->SetReceiver(this->Receiver);
  77. this->InstrumentFor->SetTimes(this->Times);
  78. nlayers = this->RecoveredMod->GetNumberOfLayers();
  79. //nrec = this->Receiver->GetNumberOfReceivers();
  80. //check for parameters set
  81. if (this->SCondIndex.size() == 0 /*&& SCondIndex(0) == 42*/) {
  82. //std::cout << "Cond Index is null" << std::endl;
  83. this->SCondIndex.resize(nlayers-1);
  84. for (int ii=1;ii<nlayers;ii++) {
  85. this->SCondIndex(ii-1)=ii;
  86. }
  87. }
  88. if (this->SThickIndex.size()==0 /*&& SThickIndex(0)==42*/) {
  89. //std::cout << "Thickness Index is null" << std::endl;
  90. this->SThickIndex.resize(nlayers-2);
  91. for (int ii=0;ii<nlayers-2;ii++) {
  92. this->SThickIndex(ii)=ii+1;
  93. }
  94. }
  95. std::cout << "Thickness indices to solve for: " << this->SThickIndex
  96. << std::endl;
  97. //std::cout << "Conductivity indices to solve for: " << this->SCondIndex
  98. // << std::endl;
  99. //Real alpha;
  100. // while \nabla g > \epsilon
  101. // Forward model from starting
  102. this->InstrumentFor->MakeLaggedCalculation(FHTKEY201);
  103. olddata = this->InstrumentFor->GetMeasurements();
  104. // Get A Jacobian Matrix
  105. jac = this->CalcJac(this->RecoveredMod);
  106. // Output jacobian matrix
  107. std::ofstream outfile2;
  108. outfile2.open("jacobian.out");
  109. for (int ii=0; ii<jac.rows();ii++) {
  110. for (int jj=0;jj<jac.cols();jj++) {
  111. outfile2 << jac(ii,jj) << " ";
  112. }
  113. outfile2 << std::endl;
  114. }
  115. outfile2.close();
  116. H = jac.transpose()*jac;
  117. std::cout << "H: " << std::endl;
  118. std::cout << H << std::endl;
  119. //Temporary solution for only one linear step
  120. //delta m = (JTJ)^-1JTd
  121. //VectorXr x0;
  122. //VectorXr md;
  123. //x0.resize(jac.cols());
  124. //x0.setZero();
  125. //VectorXr jtd;
  126. //md = this->MeasuredData.col(0);
  127. //std::cout << md << std::endl;
  128. //jtd = jac.transpose()*md;
  129. //std::cout << jtd.size() << std::endl;
  130. //Real tol;
  131. //int iter;
  132. //tol = 1.e-8;
  133. //iter = 1000;
  134. //this->RecModStep = CG(H,x0,jtd,iter,tol);
  135. //std::cout << this->RecModStep << std::endl;
  136. // Gauss-Newton Solution
  137. //
  138. // Quasi newton solution
  139. // while ||\nabla g|| > \epsilon
  140. // g = jac.transpose()*(olddata(:,1)-this->MeasuredData(:,1));
  141. // p = H*g*(-1);
  142. // newmod = oldmod + alpha * p
  143. // s = alpha * p
  144. // y = g(k+1) - g(k)
  145. // H(k+1) = H - (H*s*s.transpose()*H)/(s.T*H*s) + (y*y.T)/(y.T*s)
  146. // Solve for model perturbation
  147. // end while
  148. }
  149. MatrixXr InverseSolverTEM1D::CalcJac(LayeredEarthEM* m0) {
  150. MatrixXr jac;
  151. int nrec;
  152. int nlayers;
  153. int nparams;
  154. int nfsigs;
  155. //int nflays;
  156. Real step = 0.05;
  157. VectorXr currentsigma;
  158. VectorXr currentthick;
  159. VectorXr newsigma;
  160. VectorXr newthick;
  161. Complex ctemp;
  162. VectorXcr vctemp;
  163. VectorXr params;
  164. MatrixXr currentmeasurements;
  165. MatrixXr perturbedmeasurements;
  166. MatrixXr diffs;
  167. nrec = this->Receiver->GetNumberOfReceivers();
  168. nlayers = this->StartingMod->GetNumberOfLayers();
  169. nfsigs = this->SCondIndex.size();
  170. //nflays = this->SThickIndex.size();
  171. nparams = this->SCondIndex.size() + this->SThickIndex.size();
  172. jac.resize(nrec*this->Times.size(),nparams); /// TODO memory leak
  173. params.resize(nparams);
  174. std::cout << this->SThickIndex << std::endl;
  175. //Need to set currentsigma and currentthick
  176. //
  177. currentsigma.resize(nlayers);
  178. currentthick.resize(nlayers);
  179. currentsigma = m0->GetLayerConductivity().array().real();
  180. for (int ii=0;ii<nlayers;ii++) {
  181. currentthick(ii) = m0->GetLayerThickness(ii);
  182. }
  183. newsigma = currentsigma;
  184. newthick = currentthick;
  185. currentmeasurements = this->InstrumentFor->GetMeasurements();
  186. // Concatenate thicknesses and conductivities into one vector (params)
  187. // Here is where I write over the original parameters
  188. for (int ii=0;ii<SCondIndex.size();ii++) {
  189. params(ii)=currentsigma(this->SCondIndex(ii));
  190. }
  191. int jj = 0;
  192. for (int ii=SCondIndex.size();ii<nparams;ii++) {
  193. params(ii) = currentthick(this->SThickIndex(jj)); //this was +1
  194. jj++;
  195. }
  196. std::cout << "Size Params: " << params.size() << std::endl;
  197. std::cout << "Size newsigma: " << newsigma.size() << std::endl;
  198. std::cout << "Size newthick: " << newthick.size() << std::endl;
  199. std::cout << "Size currentsigma: " << currentsigma.size() << std::endl;
  200. std::cout << "Size currentthick: " << currentthick.size() << std::endl;
  201. std::cout << "Thickness vector : " << currentthick << std::endl;
  202. std::cout << "Sigma vector: " << currentsigma << std::endl;
  203. //std::cout << "Jacobian size: " << jac.rows() << " x " <<
  204. // jac.cols() << std::endl;
  205. Real tparam;
  206. for (int jj=0;jj<nparams;jj++) {
  207. tparam = params(jj)*step;
  208. params(jj) = params(jj) + tparam;
  209. newsigma = currentsigma;
  210. newthick = currentthick;
  211. std::cout << "New sigma, in loop: " << newsigma << std::endl;
  212. std::cout << "New thick, in loop: " << newthick << std::endl;
  213. //plop params back into newsigma and newthick
  214. for (int kk=0;kk<nfsigs;kk++) {
  215. std::cout << "Sigma index: " << this->SCondIndex(kk) << std::endl;
  216. newsigma(this->SCondIndex(kk)) = params(kk);
  217. }
  218. int ii = 0;
  219. std::cout << "Params, in loop: " << params << std::endl;
  220. for (int kk=nfsigs;kk<nparams;kk++) {
  221. std::cout << "Thick index: " << this->SThickIndex(ii) << std::endl;
  222. newthick(this->SThickIndex(ii)) = params(kk); //here is the problem
  223. //params is loaded incorrectly
  224. ii++;
  225. }
  226. std::cout << "New thickness: " << newthick << std::endl;
  227. std::cout << "New cond: " << newsigma << std::endl;
  228. m0->SetLayerThickness(newthick.segment(1,newthick.size()-2));
  229. m0->SetLayerConductivity(newsigma.cast<Complex>());
  230. std::cout << *m0 << std::endl;
  231. this->InstrumentFor->MakeLaggedCalculation(FHTKEY201);
  232. perturbedmeasurements = this->InstrumentFor->GetMeasurements();
  233. diffs = currentmeasurements - perturbedmeasurements;
  234. diffs = diffs / tparam;
  235. //Put the values back
  236. params(jj) = params(jj) - tparam;
  237. for (int ii=0;ii<jac.rows();ii++) {
  238. jac(ii,jj)=diffs(ii,1);
  239. }
  240. std::cout << "Current Jacobian: " << std::endl;
  241. std::cout << jac << std::endl;
  242. }
  243. return jac;
  244. }
  245. void InverseSolverTEM1D::AttachStartMod(LayeredEarthEM* inpstartmod) {
  246. if (this->StartingMod != NULL) {
  247. this->StartingMod->DetachFrom(this);
  248. }
  249. inpstartmod->AttachTo(this);
  250. this->StartingMod = inpstartmod;
  251. }
  252. void InverseSolverTEM1D::AttachRecMod(LayeredEarthEM* inprecmod) {
  253. if (this->RecoveredMod != NULL) {
  254. this->RecoveredMod->DetachFrom(this);
  255. }
  256. inprecmod->AttachTo(this);
  257. this->RecoveredMod = inprecmod;
  258. }
  259. void InverseSolverTEM1D::AttachWireAntenna(WireAntenna* inpwireantenna) {
  260. if (this->Antenna != NULL) {
  261. this->Antenna->DetachFrom(this);
  262. }
  263. inpwireantenna->AttachTo(this);
  264. this->Antenna = inpwireantenna;
  265. }
  266. void InverseSolverTEM1D::AttachDipoleSource(DipoleSource* inpdipsource) {
  267. if (this->Dipole != NULL) {
  268. this->Dipole->DetachFrom(this);
  269. }
  270. inpdipsource->AttachTo(this);
  271. this->Dipole = inpdipsource;
  272. }
  273. void InverseSolverTEM1D::AttachReceiver(ReceiverPoints* inpreceiver) {
  274. if (this->Receiver != NULL) {
  275. this->Receiver->DetachFrom(this);
  276. }
  277. inpreceiver->AttachTo(this);
  278. this->Receiver = inpreceiver;
  279. }
  280. void InverseSolverTEM1D::SetTimes(VectorXr &inptimes) {
  281. this->Times = inptimes;
  282. }
  283. void InverseSolverTEM1D::ShowSoln() {
  284. int ntimes;
  285. ntimes = this->InstrumentFor->GetNumberOfTimes();
  286. for (int ii=0;ii<ntimes;ii++) {
  287. std::cout<<" "<<this->InstrumentFor->GetMeasurements()(ii,0)
  288. <<" "<<this->InstrumentFor->GetMeasurements()(ii,1)<<std::endl;
  289. }
  290. }
  291. void InverseSolverTEM1D::WriteSoln(const std::string &filename) {
  292. int nlayers;
  293. VectorXcr conductivities;
  294. VectorXr thicknesses;
  295. nlayers = this->RecoveredMod->GetNumberOfLayers();
  296. conductivities.resize(nlayers);
  297. thicknesses.resize(nlayers);
  298. conductivities = this->RecoveredMod->GetLayerConductivity();
  299. for (int ii=0;ii<nlayers;ii++) {
  300. thicknesses(ii) = this->RecoveredMod->GetLayerThickness(ii);
  301. }
  302. std::ofstream outfile1;
  303. outfile1.open(filename.c_str());
  304. // Model file has format:
  305. // thickness (m) conductivity (S/m)
  306. for (int ii=0;ii<nlayers;ii++) {
  307. outfile1 << thicknesses(ii) << " " << conductivities(ii) <<
  308. " " << std::endl;
  309. }
  310. outfile1.close();
  311. }
  312. void InverseSolverTEM1D::SetFreeParams(const VectorXi &layers,
  313. const VectorXi &cond) {
  314. this->SCondIndex = cond;
  315. this->SThickIndex = layers;
  316. }
  317. void InverseSolverTEM1D::AttachMeasuredData(const MatrixXr &data) {
  318. this->MeasuredData = data;
  319. }
  320. int InverseSolverTEM1D::NumberOfIterations() {
  321. return 42;
  322. }
  323. bool InverseSolverTEM1D::Success() {
  324. return false;
  325. }
  326. VectorXr InverseSolverTEM1D::GetPhiMVector() {
  327. VectorXr temp;
  328. temp.resize(1);
  329. temp(0)=42;
  330. return temp;
  331. }
  332. VectorXr InverseSolverTEM1D::GetPhiDVector() {
  333. VectorXr temp;
  334. temp.resize(1);
  335. temp(0) = 42;
  336. return temp;
  337. }
  338. void InverseSolverTEM1D::FillInG(const Vector3r& pos,
  339. const Vector3r& step) {
  340. }
  341. }