00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144 #include <iostream>
00145 #include <iomanip>
00146 #include <fstream>
00147 #include <algorithm>
00148 #include <numeric>
00149 #include <sstream>
00150 #include <string>
00151 #include "GA.h"
00152 #include "Iso1DMTObjective.h"
00153 #include "MTStation.h"
00154 #include "C1dInvGaConf.h"
00155 #include "AbsVelRecObjective.h"
00156 #include "Multi1DRecObjective.h"
00157 #include "SeismicDataComp.h"
00158 #include "CombinedRoughness.h"
00159 #include "SeismicModelDiff.h"
00160 #include "ResPkModel.h"
00161 #include "SurfaceWaveObjective.h"
00162 #include "SurfaceWaveData.h"
00163 #include "Adaptors.h"
00164 #include "MTFitSetup.h"
00165 #include <boost/bind.hpp>
00166 #include <boost/shared_ptr.hpp>
00167 #include <boost/filesystem.hpp>
00168 #include <boost/date_time/posix_time/posix_time.hpp>
00169 #include "Util.h"
00170 #include "MTRoughness.h"
00171 #include "SurfaceWaveData.h"
00172 #include "FileAbort.h"
00173 #include "SeisTools.h"
00174
00175 using namespace std;
00176 using namespace gplib;
00177
00178 enum tgatype
00179 {
00180 pareto, anneal
00181 };
00182
00183 typedef boost::shared_ptr<GeneralObjective> pCGeneralObjective;
00184
00185
00186
00187 void SetupAnnealingGA(boost::shared_ptr<GeneralGA> &GA,
00188 const C1dInvGaConf Configuration)
00189 {
00190
00191 double InitTemperature = Configuration.inittemp;
00192
00193 if (Configuration.inittemp < 0)
00194 {
00195 double avgcombfit = 0;
00196 GA->DoIteration(0, false);
00197
00198 avgcombfit = accumulate(GA->GetAvgFit().begin(), GA->GetAvgFit().end(),
00199 0.0);
00200
00201 InitTemperature = avgcombfit * abs(Configuration.inittemp);
00202 }
00203 AnnealingGA *AnnealGA = dynamic_cast<AnnealingGA*> (GA.get());
00204
00205 AnnealGA->SetParams(InitTemperature, Configuration.annealinggeneration,
00206 Configuration.coolingratio);
00207 }
00208
00209
00210 void SetupRecObjective(const string recinfofile,
00211 Multi1DRecObjective &RecObjective, const bool normrec,
00212 RecCalc::trfmethod rfmethod)
00213 {
00214 ifstream recinfo(recinfofile.c_str());
00215 double slowness, sigma, cc, recerror, recweight, absweight;
00216 int shift;
00217 string recinputdata, absinputdata;
00218 char wantabsvel;
00219
00220 while (recinfo.good())
00221 {
00222
00223 wantabsvel = 'r';
00224
00225 recinfo >> slowness >> sigma >> cc >> recerror >> shift >> wantabsvel
00226 >> recinputdata;
00227
00228 if (wantabsvel == 'a')
00229 {
00230 recinfo >> recweight >> absweight >> absinputdata;
00231 }
00232
00233 if (recinfo.good())
00234 {
00235
00236 boost::shared_ptr<const SeismicDataComp> RecData(
00237 new SeismicDataComp(recinputdata, SeismicDataComp::sac));
00238
00239 if (wantabsvel == 'a')
00240 {
00241
00242 SurfaceWaveData AbsVelData;
00243 AbsVelData.ReadAscii(absinputdata);
00244
00245 RecObjective.AddAbsVelFunction(RecData, AbsVelData, shift,
00246 sigma, cc, slowness, rfmethod, recerror, normrec,
00247 absweight, recweight);
00248 }
00249
00250 else
00251 {
00252
00253 RecObjective.AddRecFunction(RecData, shift, sigma, cc,
00254 slowness, rfmethod, recerror, normrec);
00255 }
00256 }
00257 }
00258 }
00259
00260 int main(int argc, char *argv[])
00261 {
00262
00263 boost::posix_time::ptime starttime =
00264 boost::posix_time::microsec_clock::local_time();
00265
00266 string version = "$Id: 1dinvga.cpp 1816 2009-09-07 11:28:35Z mmoorkamp $";
00267 cout << endl << endl;
00268 cout << "Program " << version << endl;
00269 cout << "Genetic algorithm to jointly invert seismic and MT data " << endl;
00270 cout << "look at 1dinvga.conf for configuration and setup " << endl;
00271 cout << "The root name of log-files can be overwritten by using " << endl;
00272 cout << "1dinvga newrootname " << endl;
00273
00274 C1dInvGaConf Configuration;
00275 MTStation MTData;
00276 SurfaceWaveData RFAbsVel;
00277 SurfaceWaveData DispData;
00278 UniformRNG Random;
00279 Configuration.GetData("1dinvga.conf");
00280
00281 ttranscribed transcribed;
00282 tgatype GAtype = pareto;
00283
00284 string outputbase;
00285
00286 if (argc > 1)
00287 outputbase = argv[1];
00288 else
00289 outputbase = Configuration.outputbase;
00290
00291
00292 ofstream fitstat((outputbase + ".ftst").c_str());
00293 ofstream rankfile((outputbase + ".rank").c_str());
00294 ofstream frontfile((outputbase + ".front").c_str());
00295
00296 ofstream misfitfile;
00297 ofstream valuefile;
00298 ofstream probfile;
00299 ofstream distancefile;
00300
00301
00302 if (Configuration.verbose)
00303 {
00304 misfitfile.open((outputbase + ".msft").c_str());
00305 valuefile.open((outputbase + ".vals").c_str());
00306 probfile.open((outputbase + ".probs").c_str());
00307 distancefile.open((outputbase + ".dist").c_str());
00308 }
00309
00310 if ((Configuration.thickbase.size() != Configuration.resbase.size())
00311 || (Configuration.thickbase.size() != Configuration.svelbase.size()))
00312 {
00313 cerr << "Configuration of genes is not consistent ! ";
00314 return 100;
00315 }
00316
00317 const int nparam = 3;
00318
00319 const int genes = (Configuration.thickbase.size()
00320 + Configuration.resbase.size() + Configuration.svelbase.size());
00321 const int popsize = Configuration.popsize;
00322 const int maxgen = Configuration.generations;
00323 const int paramlength = genes / nparam;
00324 ttranscribed basevalues(genes);
00325 ttranscribed stepsizes(genes);
00326 tsizev genesizes(genes);
00327
00328 misfitfile << popsize << " " << maxgen << " " << endl;
00329 valuefile << genes << " " << popsize << " " << maxgen << endl;
00330
00331 copy(Configuration.resbase.begin(), Configuration.resbase.end(),
00332 basevalues.begin());
00333 copy(Configuration.resstep.begin(), Configuration.resstep.end(),
00334 stepsizes.begin());
00335 copy(Configuration.ressizes.begin(), Configuration.ressizes.end(),
00336 genesizes.begin());
00337
00338 copy(Configuration.thickbase.begin(), Configuration.thickbase.end(),
00339 basevalues.begin() + paramlength);
00340 copy(Configuration.thickstep.begin(), Configuration.thickstep.end(),
00341 stepsizes.begin() + paramlength);
00342 copy(Configuration.thicksizes.begin(), Configuration.thicksizes.end(),
00343 genesizes.begin() + paramlength);
00344
00345 copy(Configuration.svelbase.begin(), Configuration.svelbase.end(),
00346 basevalues.begin() + 2 * paramlength);
00347 copy(Configuration.svelstep.begin(), Configuration.svelstep.end(),
00348 stepsizes.begin() + 2 * paramlength);
00349 copy(Configuration.svelsizes.begin(), Configuration.svelsizes.end(),
00350 genesizes.begin() + 2 * paramlength);
00351
00352 const int totalgenesize = accumulate(genesizes.begin(), genesizes.end(), 0);
00353
00354 BinaryPopulation Population(popsize, totalgenesize, Random, true);
00355 BinaryTournamentSelect Select(Random, boost::bind(
00356 &GeneralPopulation::GetProbabilities, boost::ref(Population)),
00357 boost::bind(&GeneralPopulation::GetCrowdingDistances, boost::ref(
00358 Population)));
00359 StandardPropagation Propagator(&Select, &Population, &Random);
00360 double indmutationprob = 1.0 - pow(1.0 - Configuration.mutationprob, 1.
00361 / totalgenesize);
00362 cout << "Individual Mutation set to: " << indmutationprob << endl;
00363 Propagator.SetParams(indmutationprob, Configuration.crossoverprob);
00364 GrayTranscribe Transcribe(basevalues, stepsizes, genesizes);
00365
00366
00367 try
00368 {
00369 boost::shared_ptr<Multi1DRecObjective> RecObjective(
00370 new Multi1DRecObjective());
00371
00372 if (Configuration.weights.at(0) > 0)
00373 {
00374 MTData.GetData(Configuration.mtinputdata);
00375 }
00376 if (Configuration.weights.at(1) > 0)
00377 {
00378 RecCalc::trfmethod rfmethod = RecCalc::specdiv;
00379 if (Configuration.recmethod == "iterdecon")
00380 {
00381 rfmethod = RecCalc::iterdecon;
00382 }
00383
00384 SetupRecObjective(Configuration.recinfofile, *RecObjective.get(),
00385 Configuration.normrec, rfmethod);
00386 }
00387 if (Configuration.weights.at(2) > 0)
00388 {
00389 DispData.ReadFile(Configuration.dispdata);
00390 }
00391
00392 boost::shared_ptr<Iso1DMTObjective> MTObjective(new Iso1DMTObjective(
00393 MTData));
00394 SetupMTFitParameters(Configuration, *MTObjective);
00395 MTObjective->SetFitExponent(Configuration.mtfitexponent);
00396
00397
00398 RecObjective->SetPoisson(Configuration.poisson);
00399 RecObjective->SetTimeWindow(Configuration.starttime,
00400 Configuration.endtime);
00401
00402
00403 boost::shared_ptr<SurfaceWaveObjective> SurfObjective(
00404 new SurfaceWaveObjective(DispData));
00405 SurfObjective->SetFitExponent(Configuration.surffitexponent);
00406 SurfObjective->SetPoisson(Configuration.poisson);
00407 SurfObjective->SetErrorLevel(Configuration.surferror);
00408
00409
00410 boost::shared_ptr<GeneralObjective> RecRegularization;
00411
00412
00413 if (!Configuration.usevrefmodel)
00414 {
00415 RecRegularization = boost::shared_ptr<GeneralObjective>(
00416 new CombinedRoughness);
00417 }
00418
00419 else
00420 {
00421
00422 if (!boost::filesystem::exists(Configuration.vrefmodel))
00423 {
00424 cerr
00425 << "Cannot find reference model for regularization ! Exiting !"
00426 << endl;
00427 return 100;
00428 }
00429
00430 ResPkModel VRefModel;
00431 VRefModel.ReadModel(Configuration.vrefmodel);
00432
00433 RecRegularization = boost::shared_ptr<GeneralObjective>(
00434 new SeismicModelDiff(VRefModel));
00435 }
00436
00437
00438 boost::shared_ptr<GeneralObjective> MTRegularization(new MTRoughness());
00439
00440 vector<pCGeneralObjective> ObjVector;
00441 ObjVector.push_back(MTObjective);
00442 ObjVector.push_back(RecObjective);
00443 ObjVector.push_back(SurfObjective);
00444 ObjVector.push_back(RecRegularization);
00445 ObjVector.push_back(MTRegularization);
00446
00447 boost::shared_ptr<GeneralGA> GA;
00448 if (Configuration.gatype == "pareto")
00449 {
00450 GAtype = pareto;
00451 boost::shared_ptr<ParetoGA> Pareto(new ParetoGA(&Propagator,
00452 &Population, &Transcribe, ObjVector, Configuration.threads));
00453 GA = Pareto;
00454 }
00455 else if (Configuration.gatype == "anneal")
00456 {
00457 GAtype = anneal;
00458 boost::shared_ptr<AnnealingGA> Annealing(new AnnealingGA(
00459 &Propagator, &Population, &Transcribe, ObjVector,
00460 Configuration.threads));
00461 GA = Annealing;
00462 }
00463 else
00464 {
00465 cerr << "Invalid GA type in configuration file." << endl;
00466 cerr << " Has to be either 'pareto' or 'anneal' !" << endl;
00467 return 100;
00468 }
00469 GA->SetWeights(Configuration.weights);
00470
00471 GeneralGA::tparamindv ParamIndices;
00472 const int nmtparam = 2 * Configuration.resbase.size();
00473 std::vector<int> LocalIndices;
00474
00475 generate_n(back_inserter(LocalIndices), nmtparam, IntSequence(0));
00476 ParamIndices.push_back(LocalIndices);
00477 const int nseisparam = 2 * Configuration.thickbase.size();
00478 LocalIndices.clear();
00479
00480 generate_n(back_inserter(LocalIndices), nseisparam, IntSequence(
00481 Configuration.thickbase.size()));
00482 ParamIndices.push_back(LocalIndices);
00483 ParamIndices.push_back(LocalIndices);
00484 const int nmodelparam = 3 * Configuration.thickbase.size();
00485 LocalIndices.clear();
00486
00487 generate_n(back_inserter(LocalIndices), nmodelparam, IntSequence(0));
00488 ParamIndices.push_back(LocalIndices);
00489 ParamIndices.push_back(LocalIndices);
00490 GA->SetParameterIndices(ParamIndices);
00491 GA->SetElitist(Configuration.elitist);
00492
00493 if (GAtype == anneal)
00494 {
00495 SetupAnnealingGA(GA, Configuration);
00496 Population.InitPop();
00497 }
00498
00499 for (int i = 0; i < maxgen; ++i)
00500 {
00501 GA->DoIteration(i, i == (maxgen - 1));
00502 if (Configuration.verbose)
00503 {
00504 GA->PrintMisfit(misfitfile);
00505 GA->PrintTranscribed(valuefile);
00506 Population.PrintProbabilities(probfile);
00507 Population.PrintDistances(distancefile);
00508 }
00509 GA->PrintBestMisfit(frontfile);
00510
00511 fitstat << i << " ";
00512 GA->PrintFitStat(fitstat);
00513 if (WantAbort())
00514 {
00515 RemoveAbort();
00516
00517
00518
00519 i = maxgen - 2;
00520 }
00521 }
00522
00523 int bestcount = 0;
00524
00525
00526 UniquePop FinalUnique;
00527 tfitvec bestfit(ObjVector.size());
00528 tpopmember member(Population.GetGenesize());
00529
00530 vector<int> BestIndices(GA->GetBestModelIndices());
00531 const int noutmodels = GA->GetNBestmodels();
00532 const int nobjective = ObjVector.size();
00533
00534
00535 GeneralGA::tparamvector LocalParameters(nobjective);
00536
00537 for (int i = 0; i < nobjective; ++i)
00538 {
00539 LocalParameters.at(i).resize(ParamIndices.at(i).size(), false);
00540 }
00541
00542 ofstream finalmodels((outputbase + ".final").c_str());
00543
00544 for (int h = 0; h < noutmodels; ++h)
00545 {
00546
00547 member = row(Population.GetPopulation(), BestIndices.at(h));
00548
00549 transcribed = Transcribe.GetValues(member);
00550
00551 GA->SetupParams(transcribed, LocalParameters);
00552
00553 if (!FinalUnique.Find(transcribed, bestfit))
00554 {
00555
00556 cout << "Best: " << h << " Index : " << BestIndices.at(h)
00557 << " ";
00558
00559
00560 for (int f = 0; f < nobjective; ++f)
00561 {
00562 if (Configuration.weights.at(f) > 0)
00563 {
00564 bestfit(f) = ObjVector.at(f)->CalcPerformance(
00565 LocalParameters.at(f));
00566 }
00567 else
00568 {
00569 bestfit(f) = 0;
00570 }
00571 cout << setw(6) << setfill(' ') << setprecision(4)
00572 << bestfit(f) << " ";
00573 }
00574
00575
00576 copy(bestfit.begin(), bestfit.end(), ostream_iterator<double> (
00577 finalmodels, " "));
00578 cout << " " << h;
00579 finalmodels << " " << bestcount << endl;
00580 ostringstream filename;
00581 filename << "best_" << bestcount;
00582
00583
00584
00585 cout << " Saved as : " << filename.str();
00586 if (Configuration.weights.at(0) > 0)
00587 {
00588 MTObjective->WriteModel(filename.str() + "_mt.mod");
00589 MTObjective->WritePlot(filename.str() + "_mt.plot");
00590 MTObjective->WriteData(filename.str());
00591 }
00592 if (Configuration.weights.at(1) > 0)
00593 {
00594 RecObjective->WriteData(filename.str() + ".rec");
00595 RecObjective->WriteModel(filename.str() + "_rec.mod");
00596 RecObjective->WritePlot(filename.str() + "_rec");
00597 }
00598 if (Configuration.weights.at(2) > 0)
00599 {
00600 SurfObjective->WriteData(filename.str() + ".disp");
00601 SurfObjective->WriteModel(filename.str() + "_disp.mod");
00602 SurfObjective->WritePlot(filename.str() + "_disp.plot");
00603 }
00604 filename.str("");
00605
00606 FinalUnique.Insert(bestfit, transcribed);
00607 bestcount++;
00608 cout << endl;
00609 }
00610 }
00611
00612 if (Configuration.verbose)
00613 {
00614
00615 ofstream bestfile((outputbase + ".best").c_str());
00616 FinalUnique.PrintAll(bestfile);
00617
00618 ofstream uniquepop((outputbase + ".unpop").c_str());
00619 GA->PrintUniquePop(uniquepop);
00620 }
00621
00622 boost::posix_time::ptime endtime =
00623 boost::posix_time::microsec_clock::local_time();
00624 std::cout << "Total Runtime: " << endtime - starttime << std::endl;
00625 } catch (FatalException &e)
00626 {
00627 cerr << e.what() << endl;
00628 return -1;
00629 }
00630 }