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
00145
00146
00147 #include <iostream>
00148 #include <fstream>
00149 #include <algorithm>
00150 #include <numeric>
00151 #include <sstream>
00152 #include <string>
00153 #include "GA.h"
00154 #include "Iso1DMTObjective.h"
00155 #include "MTStation.h"
00156 #include "C1dInvGaConf.h"
00157 #include "AbsVelRecObjective.h"
00158 #include "Multi1DRecObjective.h"
00159 #include "SeismicDataComp.h"
00160 #include "CombinedRoughness.h"
00161 #include "SeismicModelDiff.h"
00162 #include "ResPkModel.h"
00163 #include "SurfaceWaveObjective.h"
00164 #include "SurfaceWaveData.h"
00165 #include "Adaptors.h"
00166 #include "MTFitSetup.h"
00167 #include <boost/bind.hpp>
00168 #include <boost/shared_ptr.hpp>
00169 #include <boost/filesystem.hpp>
00170 #include <boost/date_time/posix_time/posix_time.hpp>
00171 #include "Util.h"
00172 #include "MTRoughness.h"
00173 #include "SurfaceWaveData.h"
00174 #include "FileAbort.h"
00175 #include "SeisTools.h"
00176 enum tgatype
00177 { pareto,anneal};
00178
00179 typedef boost::shared_ptr<GeneralObjective> pCGeneralObjective;
00180 using namespace std;
00181
00182 void SetupAnnealingGA(boost::shared_ptr<GeneralGA> &GA,
00183 const C1dInvGaConf Configuration)
00184 {
00185 double InitTemperature = Configuration.inittemp;
00186 if (Configuration.inittemp < 0)
00187 {
00188 double avgcombfit = 0;
00189 GA->DoIteration(0, false);
00190
00191
00192 avgcombfit = accumulate(GA->GetAvgFit().begin(), GA->GetAvgFit().end(), 0.0);
00193 InitTemperature = avgcombfit * abs(Configuration.inittemp);
00194 }
00195 AnnealingGA *AnnealGA = dynamic_cast<AnnealingGA*>(GA.get());
00196 AnnealGA->SetParams(InitTemperature, Configuration.annealinggeneration,
00197 Configuration.coolingratio);
00198 }
00199
00200 void SetupRecObjective(const string recinfofile,
00201 Multi1DRecObjective &RecObjective)
00202 {
00203 ifstream recinfo(recinfofile.c_str());
00204 double slowness, omega, sigma, cc, recerror, recweight, absweight;
00205 int shift;
00206 string recinputdata, absinputdata;
00207 char wantabsvel;
00208 while (recinfo.good())
00209 {
00210 wantabsvel = 'r';
00211 recinfo >> slowness >> omega >> sigma >> cc >> recerror >> shift
00212 >> wantabsvel >> recinputdata;
00213 if (wantabsvel == 'a')
00214 {
00215 recinfo >> recweight >> absweight >> absinputdata;
00216 }
00217 if (recinfo.good())
00218 {
00219 boost::shared_ptr<const SeismicDataComp>
00220 RecData(new SeismicDataComp(recinputdata, SeismicDataComp::sac));
00221 Normalize(const_cast<SeismicDataComp*>(RecData.get())->GetData());
00222 if (wantabsvel == 'a')
00223 {
00224 SurfaceWaveData AbsVelData;
00225 AbsVelData.ReadAscii(absinputdata);
00226 RecObjective.AddAbsVelFunction(RecData, AbsVelData, shift,
00227 omega, sigma, cc, slowness, RecCalc::iterdecon, recerror,
00228 absweight, recweight);
00229 }
00230 else
00231 {
00232 RecObjective.AddRecFunction(RecData, shift, omega, sigma, cc,
00233 slowness, RecCalc::iterdecon, recerror);
00234 }
00235 }
00236 }
00237 }
00238
00239 int main(int argc, char *argv[])
00240 {
00241 boost::posix_time::ptime starttime =
00242 boost::posix_time::microsec_clock::local_time();
00243 string version = "$Id: 1dinvga.cpp 1696 2008-05-20 12:08:38Z mmoorkamp $";
00244 cout << endl << endl;
00245 cout << "Program " << version << endl;
00246 cout << "Genetic algorithm to jointly invert seismic and MT data " << endl;
00247 cout << "look at 1dinvga.conf for configuration and setup " << endl;
00248 cout << "The root name of log-files can be overwritten by using " << endl;
00249 cout << "1dinvga newrootname " << endl;
00250
00251 C1dInvGaConf Configuration;
00252 MTStation MTData;
00253
00254 SurfaceWaveData RFAbsVel;
00255 SurfaceWaveData DispData;
00256 UniformRNG Random;
00257 Configuration.GetData("1dinvga.conf");
00258
00259 ttranscribed transcribed;
00260 tgatype GAtype = pareto;
00261
00262 string outputbase;
00263
00264 if (argc > 1)
00265 outputbase = argv[1];
00266 else
00267 outputbase = Configuration.outputbase;
00268
00269 ofstream misfitfile;
00270 ofstream valuefile;
00271 ofstream probfile;
00272 ofstream distancefile;
00273 ofstream fitstat((outputbase+".ftst").c_str());
00274 ofstream uniquepop((outputbase+".unpop").c_str());
00275 ofstream rankfile((outputbase+".rank").c_str());
00276 ofstream frontfile((outputbase+".front").c_str());
00277 ofstream bestfile((outputbase+".best").c_str());
00278
00279 if (Configuration.verbose)
00280 {
00281 misfitfile.open((outputbase+".msft").c_str());
00282 valuefile.open((outputbase+".vals").c_str());
00283 probfile.open((outputbase+".probs").c_str());
00284 distancefile.open((outputbase+".dist").c_str());
00285 }
00286
00287 if ( (Configuration.thickbase.size() != Configuration.resbase.size())
00288 || (Configuration.thickbase.size() != Configuration.svelbase.size()))
00289 {
00290 cerr << "Configuration of genes is not consistent ! ";
00291 return 100;
00292 }
00293 const int nparam = 3;
00294 const int genes = (Configuration.thickbase.size()
00295 + Configuration.resbase.size() + Configuration.svelbase.size());
00296 const int popsize = Configuration.popsize;
00297 const int maxgen = Configuration.generations;
00298 const int paramlength = genes/nparam;
00299 ttranscribed basevalues(genes);
00300 ttranscribed stepsizes(genes);
00301 tsizev genesizes(genes);
00302
00303 misfitfile << popsize << " " << maxgen << " " << endl;
00304 valuefile << genes << " "<< popsize << " " << maxgen << endl;
00305
00306 copy(Configuration.resbase.begin(), Configuration.resbase.end(),
00307 basevalues.begin());
00308 copy(Configuration.resstep.begin(), Configuration.resstep.end(),
00309 stepsizes.begin());
00310 copy(Configuration.ressizes.begin(), Configuration.ressizes.end(),
00311 genesizes.begin());
00312
00313 copy(Configuration.thickbase.begin(), Configuration.thickbase.end(),
00314 basevalues.begin()+paramlength);
00315 copy(Configuration.thickstep.begin(), Configuration.thickstep.end(),
00316 stepsizes.begin()+paramlength);
00317 copy(Configuration.thicksizes.begin(), Configuration.thicksizes.end(),
00318 genesizes.begin()+paramlength);
00319
00320 copy(Configuration.svelbase.begin(), Configuration.svelbase.end(),
00321 basevalues.begin()+2*paramlength);
00322 copy(Configuration.svelstep.begin(), Configuration.svelstep.end(),
00323 stepsizes.begin()+2*paramlength);
00324 copy(Configuration.svelsizes.begin(), Configuration.svelsizes.end(),
00325 genesizes.begin()+2*paramlength);
00326
00327 const int totalgenesize = accumulate(genesizes.begin(), genesizes.end(), 0);
00328 BinaryPopulation Population(popsize, totalgenesize, Random, true);
00329 BinaryTournamentSelect Select(Random, boost::bind(
00330 &GeneralPopulation::GetProbabilities, boost::ref(Population)),
00331 boost::bind(&GeneralPopulation::GetCrowdingDistances,
00332 boost::ref(Population)));
00333 StandardPropagation Propagator(&Select, &Population, &Random);
00334 double indmutationprob = 1.0 - pow(1.0 - Configuration.mutationprob, 1.
00335 /totalgenesize);
00336 cout << "Individual Mutation set to: " << indmutationprob << endl;
00337 Propagator.SetParams(indmutationprob, Configuration.crossoverprob);
00338 GrayTranscribe Transcribe(basevalues, stepsizes, genesizes);
00339
00340
00341 try
00342 {
00343 boost::shared_ptr<Multi1DRecObjective> RecObjective(new Multi1DRecObjective());
00344
00345 if (Configuration.weights.at(0)> 0)
00346 {
00347 MTData.GetData(Configuration.mtinputdata);
00348 }
00349 if (Configuration.weights.at(1)> 0)
00350 {
00351 SetupRecObjective(Configuration.recinfofile,*RecObjective.get());
00352 }
00353 if (Configuration.weights.at(2)> 0)
00354 {
00355 DispData.ReadFile(Configuration.dispdata);
00356 }
00357
00358 boost::shared_ptr<Iso1DMTObjective> MTObjective(new Iso1DMTObjective(MTData));
00359 SetupMTFitParameters(Configuration,*MTObjective);
00360 MTObjective->SetFitExponent(Configuration.mtfitexponent);
00361
00362
00363 RecCalc::trfmethod rfmethod = RecCalc::specdiv;
00364 if (Configuration.recmethod == "iterdecon")
00365 {
00366 rfmethod = RecCalc::iterdecon;
00367 }
00368
00369 RecObjective->SetPoisson(Configuration.poisson);
00370 RecObjective->SetTimeWindow(Configuration.starttime,Configuration.endtime);
00371
00372
00373 boost::shared_ptr<SurfaceWaveObjective> SurfObjective(new SurfaceWaveObjective(DispData));
00374 SurfObjective->SetFitExponent(Configuration.surffitexponent);
00375 SurfObjective->SetPoisson(Configuration.poisson);
00376 SurfObjective->SetErrorLevel(Configuration.surferror);
00377
00378
00379 boost::shared_ptr<GeneralObjective> RecRegularization;
00380 if (!Configuration.usevrefmodel)
00381 {
00382 RecRegularization = boost::shared_ptr<GeneralObjective>(new CombinedRoughness);
00383 }
00384 else
00385 {
00386 if (!boost::filesystem::exists(Configuration.vrefmodel))
00387 {
00388 cerr << "Cannot find reference model for regularization ! Exiting !" << endl;
00389 return 100;
00390 }
00391 ResPkModel VRefModel;
00392 VRefModel.GetData(Configuration.vrefmodel);
00393 RecRegularization = boost::shared_ptr<GeneralObjective>(new SeismicModelDiff(VRefModel));
00394 }
00395 boost::shared_ptr<GeneralObjective> MTRegularization (new MTRoughness());
00396 vector<pCGeneralObjective> ObjVector;
00397 ObjVector.push_back(MTObjective);
00398 ObjVector.push_back(RecObjective);
00399 ObjVector.push_back(SurfObjective);
00400 ObjVector.push_back(RecRegularization);
00401 ObjVector.push_back(MTRegularization);
00402
00403 boost::shared_ptr<GeneralGA> GA;
00404 if (Configuration.gatype == "pareto")
00405 {
00406 GAtype = pareto;
00407 boost::shared_ptr<ParetoGA> Pareto(new ParetoGA(&Propagator,&Population,&Transcribe,ObjVector));
00408 GA = Pareto;
00409 }
00410 else if (Configuration.gatype == "anneal")
00411 {
00412 GAtype = anneal;
00413 boost::shared_ptr<AnnealingGA> Annealing(new AnnealingGA(&Propagator,&Population,&Transcribe,ObjVector));
00414 GA = Annealing;
00415 }
00416 else
00417 {
00418 cerr << "Invalid GA type in configuration file." << endl;
00419 cerr << " Has to be either 'pareto' or 'anneal' !" << endl;
00420 return 100;
00421 }
00422 GA->SetWeights(Configuration.weights);
00423
00424 GeneralGA::tparamindv ParamIndices;
00425 const int nmtparam = 2 * Configuration.resbase.size();
00426 std::vector<int> LocalIndices;
00427
00428 generate_n(back_inserter(LocalIndices),nmtparam,IntSequence(0));
00429 ParamIndices.push_back(LocalIndices);
00430 const int nseisparam = 2 * Configuration.thickbase.size();
00431 LocalIndices.clear();
00432
00433 generate_n(back_inserter(LocalIndices),nseisparam,IntSequence(Configuration.thickbase.size()));
00434 ParamIndices.push_back(LocalIndices);
00435 ParamIndices.push_back(LocalIndices);
00436 const int nmodelparam = 3 * Configuration.thickbase.size();
00437 LocalIndices.clear();
00438
00439 generate_n(back_inserter(LocalIndices),nmodelparam,IntSequence(0));
00440 ParamIndices.push_back(LocalIndices);
00441 ParamIndices.push_back(LocalIndices);
00442 GA->SetParameterIndices(ParamIndices);
00443 GA->SetThreads(Configuration.threads);
00444 GA->SetElitist(Configuration.elitist);
00445
00446 if (GAtype == anneal)
00447 {
00448 SetupAnnealingGA(GA,Configuration);
00449 Population.InitPop();
00450 }
00451
00452 for (int i = 0; i < maxgen; ++i)
00453 {
00454 GA->DoIteration(i, i == (maxgen-1));
00455 if (Configuration.verbose)
00456
00457 {
00458 GA->PrintMisfit(misfitfile);
00459 GA->PrintTranscribed(valuefile);
00460 Population.PrintProbabilities(probfile);
00461 Population.PrintDistances(distancefile);
00462 }
00463 GA->PrintBestMisfit(frontfile);
00464
00465 fitstat << i << " ";
00466 GA->PrintFitStat(fitstat);
00467 if (WantAbort())
00468
00469 {
00470 RemoveAbort();
00471 i = maxgen;
00472 }
00473 }
00474
00475 int bestcount = 0;
00476 UniquePop FinalUnique;
00477 tfitvec bestfit(ObjVector.size());
00478 tpopmember member(Population.GetGenesize());
00479
00480 vector<int> BestIndices(GA->GetBestModelIndices());
00481 const int noutmodels = GA->GetNBestmodels();
00482 const int nobjective = ObjVector.size();
00483 GeneralGA::tparamvector LocalParameters(nobjective);
00484 for (int i = 0; i < nobjective; ++i)
00485 {
00486 LocalParameters.at(i).resize(ParamIndices.at(i).size(),false);
00487 }
00488 ofstream finalmodels((outputbase+".final").c_str());
00489 for (int h= 0; h < noutmodels; ++h)
00490 {
00491 member = row(Population.GetPopulation(),BestIndices.at(h));
00492 cout << "Best: " << h << " Index : " << BestIndices.at(h) << " ";
00493
00494 transcribed = Transcribe.GetValues(member);
00495 GA->SetupParams(transcribed,LocalParameters);
00496 for (int f = 0; f < nobjective; ++f)
00497 {
00498 if (Configuration.weights.at(f)> 0)
00499 {
00500 bestfit(f) = ObjVector.at(f)->CalcPerformance(LocalParameters.at(f));
00501 }
00502 else
00503 {
00504 bestfit(f) = 0;
00505 }
00506 cout << bestfit(f) << " ";
00507 }
00508
00509 if (FinalUnique.Insert(bestfit,transcribed))
00510 {
00511 copy(bestfit.begin(),bestfit.end(),ostream_iterator<double>(finalmodels," "));
00512 cout << " " << h;
00513 finalmodels << " " << bestcount << endl;
00514 ostringstream filename;
00515 filename << "best_" << bestcount;
00516 cout << " Saved as : " << filename.str();
00517 if (Configuration.weights.at(0)> 0)
00518 {
00519 MTObjective->WriteModel(filename.str()+"_mt.mod");
00520 MTObjective->WritePlot(filename.str()+"_mt.plot");
00521 MTObjective->WriteData(filename.str());
00522 }
00523 if (Configuration.weights.at(1)> 0)
00524 {
00525 RecObjective->WriteData(filename.str()+".rec");
00526 RecObjective->WriteModel(filename.str()+"_rec.mod");
00527 RecObjective->WritePlot(filename.str()+"_rec");
00528 }
00529 if (Configuration.weights.at(2)> 0)
00530 {
00531 SurfObjective->WriteData(filename.str()+".disp");
00532 SurfObjective->WriteModel(filename.str()+"_disp.mod");
00533 SurfObjective->WritePlot(filename.str()+"_disp.plot");
00534 }
00535 filename.str("");
00536 bestcount++;
00537 }
00538 cout << endl;
00539 }
00540 FinalUnique.PrintAll(bestfile);
00541 GA->PrintUniquePop(uniquepop);
00542 boost::posix_time::ptime endtime = boost::posix_time::microsec_clock::local_time();
00543 std::cout << "Total Runtime: " << endtime - starttime << std::endl;
00544 }
00545 catch(CFatalException &e)
00546 {
00547 cerr << e.what() << endl;
00548 return -1;
00549 }
00550 }