#include <DinverCore.h>#include <DinverDCCore.h>#include <QGpCoreTools.h>#include <QGpCoreWave.h>#include "gpdcmisfitVersion.h"#include "gpdcmisfitInstallPath.h"Enumerations | |
| enum | AppMode { Report, SurfaceModel, RefractionModel } |
Functions | |
| ApplicationHelp * | help () |
| int | main (int argc, char **argv) |
| PACKAGE_INFO (gpdcmisfit, GPDCMISFIT) | |
| int | refractionModelMode () |
| int | reportMode () |
| int | surfaceModelMode () |
Variables | |
| int | akaikeDof = 1 |
| int | bestOutputCount = 0 |
| QStringList | inputFileList |
| double | maxInputMisfit = 1e99 |
| int | maxOutputCount = 0 |
| double | maxOutputMisfit = 1e99 |
| AppMode | mode = Report |
| QList< DCModelInfo * > | models |
| ReportWriter::Action | reportAction = ReportWriter::Ask |
| QString | reportFile |
| TargetList * | targets = 0 |
| enum AppMode |
| ApplicationHelp* help | ( | ) |
| int main | ( | int | argc, |
| char ** | argv | ||
| ) |
References QGpCoreTools::Akaike, akaikeDof, QGpCoreTools::AkaikeFewSamples, DinverDCCore::TargetList::autocorrTarget(), bestOutputCount, DinverDCCore::TargetList::dispersionTarget(), DinverDCCore::TargetList::ellipticityCurveTarget(), DinverDCCore::TargetList::ellipticityPeakTarget(), QGpCoreTools::endl(), help(), inputFileList, QGpCoreTools::L2_Normalized, DinverDCCore::TargetList::magnetoTelluricTarget(), maxInputMisfit, maxOutputCount, maxOutputMisfit, mode, RefractionModel, refractionModelMode(), DinverDCCore::TargetList::refractionVpTarget(), DinverDCCore::TargetList::refractionVsTarget(), Report, reportAction, reportFile, reportMode(), DinverDCCore::Target::selected(), DinverDCCore::Target::setMinimumMisfit(), DinverDCCore::TargetList::setMisfitType(), DinverDCCore::Target::setSelected(), SurfaceModel, surfaceModelMode(), targets, QGpCoreTools::tr(), DinverDCCore::TargetList::validateTargets(), and QGpCoreTools::XMLHeader::xml_restoreFile().
{
CoreApplication a(argc, argv, help);
targets=new TargetList;
// Options
QString targetFile; // Target file to consider in misfit computation
int i, j=1;
// Check arguments for Target file only
for(i=1; i<argc; i++) {
QByteArray arg=argv[i];
if(arg[0]=='-') {
if(arg=="-target") {
CoreApplication::checkOptionArg(i, argc, argv);
targetFile=argv[i];
} else {
argv[j++]=argv[i];
}
} else {
argv[j++]=argv[i];
}
}
if(j < argc) {
argv[j]=0;
argc=j;
}
if(targetFile.isEmpty()) {
App::stream() << tr("gpdcmisfit: option -target is mandatory") << endl;
delete targets;
return 2;
}
XMLVirtualPlugin plugin(targets, "DispersionCurve");
XMLDinverHeader hdr(&plugin);
if(hdr.xml_restoreFile(targetFile)!=XMLClass::NoError) {
App::stream() << tr("gpdcmisfit: error parsing target file") << endl;
delete targets;
return 2;
}
// By default no misfit component is selected
targets->dispersionTarget().setSelected(false);
targets->autocorrTarget().setSelected(false);
targets->ellipticityCurveTarget().setSelected(false);
targets->ellipticityPeakTarget().setSelected(false);
targets->refractionVpTarget().setSelected(false);
targets->refractionVsTarget().setSelected(false);
MisfitType misfitType=L2_Normalized;
bool forceMisfitType=false;
// Check arguments
j=1;
for(i=1; i<argc; i++) {
QByteArray arg=argv[i];
if(arg[0]=='-') {
if(arg=="-report") {
mode=Report;
} else if(arg=="-surf-model") {
mode=SurfaceModel;
} else if(arg=="-refra-model") {
mode=RefractionModel;
} else if(arg=="-akaike") {
CoreApplication::checkOptionArg(i, argc, argv);
misfitType=Akaike;
akaikeDof=atoi(argv[i] );
forceMisfitType=true;
} else if(arg=="-akaike-few") {
CoreApplication::checkOptionArg(i, argc, argv);
misfitType=AkaikeFewSamples;
akaikeDof=atoi(argv[i] );
forceMisfitType=true;
} else if(arg=="-o") {
CoreApplication::checkOptionArg(i, argc, argv);
reportFile=argv[i];
} else if(arg=="-m-in") {
CoreApplication::checkOptionArg(i, argc, argv);
maxInputMisfit=atof(argv[i])*(1.0+1e-15);
} else if(arg=="-m-out") {
CoreApplication::checkOptionArg(i, argc, argv);
maxOutputMisfit=atof(argv[i])*(1.0+1e-15);
} else if(arg=="-best-out") {
CoreApplication::checkOptionArg(i, argc, argv);
bestOutputCount=atoi(argv[i] );
} else if(arg=="-n-out") {
CoreApplication::checkOptionArg(i, argc, argv);
maxOutputCount=atoi(argv[i]);
} else if(arg=="-all") {
CoreApplication::checkOptionArg(i, argc, argv);
targets->dispersionTarget().setSelected(true);
targets->autocorrTarget().setSelected(true);
targets->ellipticityCurveTarget().setSelected(true);
targets->ellipticityPeakTarget().setSelected(true);
targets->refractionVpTarget().setSelected(true);
targets->refractionVsTarget().setSelected(true);
double m=atof(argv[i]);
targets->dispersionTarget().setMinimumMisfit(m);
targets->autocorrTarget().setMinimumMisfit(m);
targets->ellipticityCurveTarget().setMinimumMisfit(m);
targets->ellipticityPeakTarget().setMinimumMisfit(m);
targets->refractionVpTarget().setMinimumMisfit(m);
targets->refractionVsTarget().setMinimumMisfit(m);
} else if(arg=="-disp") {
CoreApplication::checkOptionArg(i, argc, argv);
targets->dispersionTarget().setSelected(true);
targets->dispersionTarget().setMinimumMisfit(atof(argv[i]));
} else if(arg=="-autocorr") {
CoreApplication::checkOptionArg(i, argc, argv);
targets->autocorrTarget().setSelected(true);
targets->autocorrTarget().setMinimumMisfit(atof(argv[i]));
} else if(arg=="-ellcurve") {
CoreApplication::checkOptionArg(i, argc, argv);
targets->ellipticityCurveTarget().setSelected(true);
targets->ellipticityCurveTarget().setMinimumMisfit(atof(argv[i]));
} else if(arg=="-ellpeak") {
CoreApplication::checkOptionArg(i, argc, argv);
targets->ellipticityPeakTarget().setSelected(true);
targets->ellipticityPeakTarget().setMinimumMisfit(atof(argv[i]));
} else if(arg=="-refravp") {
CoreApplication::checkOptionArg(i, argc, argv);
targets->refractionVpTarget().setSelected(true);
targets->refractionVpTarget().setMinimumMisfit(atof(argv[i]));
} else if(arg=="-refravs") {
CoreApplication::checkOptionArg(i, argc, argv);
targets->refractionVsTarget().setSelected(true);
targets->refractionVsTarget().setMinimumMisfit(atof(argv[i]));
} else if(arg=="-mt") {
CoreApplication::checkOptionArg(i, argc, argv);
targets->magnetoTelluricTarget().setSelected(true);
targets->magnetoTelluricTarget().setMinimumMisfit(atof(argv[i]));
} else if(arg=="-f") {
reportAction=ReportWriter::Overwrite;
} else if(arg=="-resume") {
reportAction=ReportWriter::Append;
} else {
App::stream() << tr("gpdcmisfit: bad option %1, see --help").arg(argv[i]) << endl;
delete targets;
return 2;
}
} else {
argv[j++]=argv[i];
}
}
if(j < argc) {
argv[j]=0;
argc=j;
}
if(!targets->dispersionTarget().selected() &&
!targets->autocorrTarget().selected() &&
!targets->ellipticityCurveTarget().selected() &&
!targets->ellipticityPeakTarget().selected() &&
!targets->refractionVpTarget().selected() &&
!targets->refractionVsTarget().selected() &&
!targets->magnetoTelluricTarget().selected()) {
App::stream() << tr("gpdcmisfit: no sub-target selected, see --help") << endl;
delete targets;
return 2;
}
for(int i =1;i<argc;i++) {
inputFileList.append(argv[i]);
}
if(forceMisfitType) {
targets->setMisfitType(misfitType);
}
targets->validateTargets();
int exitCode=0;
switch(mode) {
case Report:
exitCode=reportMode();
break;
case SurfaceModel:
exitCode=surfaceModelMode();
break;
case RefractionModel:
exitCode=refractionModelMode();
break;
}
delete targets;
return exitCode;
}
| PACKAGE_INFO | ( | gpdcmisfit | , |
| GPDCMISFIT | |||
| ) |
| int refractionModelMode | ( | ) |
References akaikeDof, QGpCoreTools::endl(), QGpCoreWave::RefractionDippingModel::fromStream(), inputFileList, QGpCoreWave::RefractionDippingModel::layerCount(), DinverDCCore::TargetList::refractionMisfit(), QGpCoreWave::RefractionDippingModel::toStream(), and QGpCoreTools::tr().
Referenced by main().
{
RefractionDippingModel m;
QFile * f;
QStringList::iterator itFile=inputFileList.begin();
if(itFile==inputFileList.end()) {
f=new QFile;
f->open(stdin,QIODevice::ReadOnly);
CoreApplication::instance()->debugUserInterrupts(false);
} else {
f=0;
}
while(true) {
if(f && f->atEnd()) {
delete f;
f=0;
// Even if stdin was not in use, it does not hurt.
CoreApplication::instance()->debugUserInterrupts(true);
}
if(!f) {
if(itFile!=inputFileList.end()) {
f=new QFile(*itFile);
if(!f->open(QIODevice::ReadOnly)) {
App::stream() << tr("Cannot open file %1 for reading.").arg(*itFile) << endl;
delete f;
return 2;
}
itFile++;
} else break;
}
QString comments;
QTextStream s(f);
if(!m.fromStream(s, &comments)) {
break;
}
if(m.layerCount()>0) {
double totalMisfit=0.0;
double totalWeight=0.0;
bool ok=true;
targets->refractionMisfit(totalMisfit, totalWeight, &m, akaikeDof,ok);
if(ok && totalWeight>0) {
s << totalMisfit/totalWeight << endl;
m.toStream(s);
}
} else {
break;
}
}
delete f;
return 0;
}
| int reportMode | ( | ) |
References DinverCore::ReportWriter::addModel(), QGpCoreTools::SharedObject::addReference(), akaikeDof, bestOutputCount, QGpCoreWave::Profile::depths(), QGpCoreTools::endl(), iModel, DinverDCCore::DCModelInfo::indexInReport(), inputFileList, DinverDCCore::TargetList::magnetoTelluricMisfit(), maxInputMisfit, maxOutputCount, maxOutputMisfit, misfit(), DinverCore::ReportReader::misfit(), DinverDCCore::DCModelInfo::misfit(), models, DinverCore::ReportReader::open(), DinverCore::ReportWriter::open(), outputReport, DinverDCCore::DCReportBlock::pitch(), DinverDCCore::DCReportBlock::readElectricModel(), DinverDCCore::DCReportBlock::readProfiles(), QGpCoreWave::Profile::readReport(), DinverDCCore::DCReportBlock::readSeismicModel(), DinverDCCore::TargetList::refractionMisfit(), DinverDCCore::DCModelInfo::report(), reportAction, reportFile, reports, QGpCoreWave::Profile::resample(), DinverDCCore::DCModelInfo::setIndexInReport(), DinverDCCore::DCModelInfo::setMisfit(), DinverCore::Parameter::setRealValue(), DinverDCCore::DCModelInfo::setReport(), sOut(), DinverCore::ReportReader::stream(), DinverDCCore::TargetList::surfaceMisfit(), DinverCore::ReportReader::synchronize(), QGpCoreTools::tr(), QGpCoreWave::Profile::values(), and QGpCoreWave::Seismic1DModel::vpProfile().
Referenced by main().
{
ReportWriter * outputReport;
if(!reportFile.isEmpty()) {
if(!ReportWriter::initReport(reportFile, tr("Open report for writing"), reportAction)) {
return 2;
}
outputReport=new ReportWriter(reportFile);
if(!outputReport->open()) {
App::stream() << tr("gpdcmisfit: cannot open report file for writing") << endl;
delete outputReport;
return 2;
}
} else {
outputReport=0;
}
QTextStream sOut(stdout);
int iModel=0;
//report->addModel(totalMisfit/totalWeight, checksum, nDimensions, _varParamList);
sOut << "# Index in report | New misfit | Old misfit\n";
QList<ReportReader *> reports;
for(QStringList::iterator it=inputFileList.begin(); it!=inputFileList.end(); it++) {
ReportReader * report=new ReportReader(*it);
if(!report->open()) {
qDeleteAll(models);
foreach(ReportReader * report, reports) {
ReportReader::removeReference(report);
}
delete report;
delete outputReport;
return 2;
}
report->addReference();
reports << report;
report->synchronize();
int nReportModels=report->nModels();
sOut << "# Report=" << *it << "\n";
sOut << "# N Models=" << nReportModels << "\n";
QDataStream& s=report->stream();
for(int iReportModel=0; iReportModel<nReportModels; iReportModel++) {
double oldMisfit=report->misfit(iReportModel);
if(oldMisfit<maxInputMisfit) {
int nParams;
uint cs;
double val;
qint64 blockOffset=s.device()->pos();
s >> nParams;
s >> cs;
Parameter * params[nParams];
for(int id=0 ; id<nParams; id++) {
s >> val;
params[id]=new Parameter;
params[id]->setRealValue(val);
}
// Read layered model
s.device()->seek(blockOffset);
int reportDCVersion=report->userBlockVersion("DISP");
if(reportDCVersion>=0) {
DCReportBlock dcBlock(s);
dcBlock.readProfiles(reportDCVersion);
Seismic1DModel * sm=dcBlock.readSeismicModel();
Resistivity1DModel * em=dcBlock.readElectricModel();
// Compute new misfit
double totalMisfit=0.0;
double totalWeight=0.0;
bool ok=true;
if(sm) {
targets->surfaceMisfit(totalMisfit, totalWeight, sm, akaikeDof, ok);
if(dcBlock.pitch()) {
Profile pitch;
pitch.readReport(s);
pitch.resample(sm->vpProfile().depths());
targets->refractionMisfit(totalMisfit, totalWeight, sm, pitch.values(), akaikeDof, ok);
}
}
if(em) {
// No static shift: TODO provide a way to identify parameters values
targets->magnetoTelluricMisfit(totalMisfit, totalWeight, em, 1.0, akaikeDof, ok);
}
if(!sm && !em) {
// Really nothing to retrieve skip this model
for(int id=0; id<nParams; id++) delete params[id];
continue;
}
if(ok && totalWeight>0) {
double misfit=totalMisfit/totalWeight;
if(misfit<maxOutputMisfit) {
if(bestOutputCount>0) {
DCModelInfo * info=new DCModelInfo;
info->setReport(report);
info->setIndexInReport(iReportModel);
info->setMisfit(misfit);
info->addReference();
models << info;
} else {
sOut << QString("%1 %2 %3").arg(iReportModel+1).arg(oldMisfit).arg(misfit) << endl;
if(outputReport) {
// Save new misfit and the rest of the model entry (param, disp,...)
outputReport->addModel(misfit, cs, nParams, params);
s.device()->seek(blockOffset);
DCReportBlock::write(outputReport, report);
}
iModel++;
if(iModel==maxOutputCount) {
for(int id=0 ; id<nParams; id++) {
delete params[id];
}
delete sm;
delete em;
break;
}
}
}
delete sm;
delete em;
}
for(int id=0; id<nParams; id++) {
delete params[id];
}
}
}
}
}
if(bestOutputCount>0) {
qSort(models.begin(), models.end(), DCModelInfo::misfitLessThan);
int iStop=models.count() - maxOutputCount;
if(iStop<0) iStop=0;
for(int i=models.count()-1; i>=iStop; i-- ) {
DCModelInfo * info=models.at(i);
ReportReader * report=info->report();
double oldMisfit=report->misfit(info->indexInReport());
sOut << QString("%1 %2 %3\n").arg(info->indexInReport()+1).arg(oldMisfit).arg(info->misfit());
if(outputReport) {
QDataStream& s=report->stream();
int nParams;
uint cs;
double val;
qint64 blockOffset=s.device()->pos();
s >> nParams;
s >> cs;
Parameter * params[nParams];
for(int id=0 ; id<nParams; id++ ) {
s >> val;
params[id]=new Parameter;
params[id]->setRealValue(val);
}
// Save new misfit and the rest of the model entry (param, disp,...)
outputReport->addModel(info->misfit(), cs, nParams, params);
for(int id=0 ; id<nParams; id++ ) delete params[id];
s.device()->seek(blockOffset);
DCReportBlock::write(outputReport, report);
}
iModel++;
}
}
qDeleteAll(models);
foreach(ReportReader * report, reports) ReportReader::removeReference(report);
if(outputReport) {
sOut << QString("# Saved %1 misfits in %2").arg(iModel).arg(reportFile) << endl;
delete outputReport;
}
return 0;
}
| int surfaceModelMode | ( | ) |
References akaikeDof, QGpCoreTools::endl(), QGpCoreWave::Seismic1DModel::fromStream(), inputFileList, QGpCoreWave::Seismic1DModel::layerCount(), DinverDCCore::TargetList::surfaceMisfit(), QGpCoreWave::Seismic1DModel::toStream(), and QGpCoreTools::tr().
Referenced by main().
{
Seismic1DModel m;
QFile * f;
QStringList::iterator itFile=inputFileList.begin();
if(itFile==inputFileList.end()) {
f=new QFile;
f->open(stdin,QIODevice::ReadOnly);
CoreApplication::instance()->debugUserInterrupts(false);
} else {
f=0;
}
while(true) {
if(f && f->atEnd()) {
delete f;
f=0;
// Even if stdin was not in use, it does not hurt.
CoreApplication::instance()->debugUserInterrupts(true);
}
if(!f) {
if(itFile!=inputFileList.end()) {
f=new QFile(*itFile);
if(!f->open(QIODevice::ReadOnly)) {
App::stream() << tr("Cannot open file %1 for reading.").arg(*itFile) << endl;
delete f;
return 2;
}
itFile++;
} else break;
}
QString comments;
QTextStream s(f);
if(!m.fromStream(s, &comments)) {
break;
}
if(m.layerCount()>0) {
double totalMisfit=0.0;
double totalWeight=0.0;
bool ok=true;
targets->surfaceMisfit(totalMisfit, totalWeight, &m, akaikeDof, ok);
if(ok && totalWeight>0) {
s << totalMisfit/totalWeight << endl;
m.toStream(s);
}
} else {
break;
}
}
delete f;
return 0;
}
| int akaikeDof = 1 |
Referenced by main(), refractionModelMode(), reportMode(), and surfaceModelMode().
| int bestOutputCount = 0 |
Referenced by main(), and reportMode().
| QStringList inputFileList |
Referenced by main(), refractionModelMode(), reportMode(), and surfaceModelMode().
| double maxInputMisfit = 1e99 |
Referenced by main(), and reportMode().
| int maxOutputCount = 0 |
Referenced by main(), and reportMode().
| double maxOutputMisfit = 1e99 |
Referenced by main(), and reportMode().
Referenced by QGpCompatibility::CompatMultiModalData::allocatesData(), QGpCompatibility::CompatMultiModalCurves::allocatesValues(), QGpCompatibility::CompatDispersion::checkSlopes(), QGpCompatibility::CompatMultiModalData::checkStdDev(), QGpCompatibility::CompatDispersionData::closestModeMisfit(), QGpCompatibility::CompatDispersionData::convertStddev(), QGpCompatibility::CompatMultiModalData::dataToReport(), QGpCompatibility::CompatMultiModalData::deleteData(), QGpCompatibility::CompatMultiModalCurves::deleteValues(), dispersion_curve_love_(), dispersion_curve_rayleigh_(), QGpCompatibility::CompatMultiModalData::isSameData(), main(), QGpCompatibility::CompatDispersionData::maxDataFrequency(), QGpCompatibility::CompatMultiModalData::measurement(), QGpCompatibility::CompatDispersionData::minDataFrequency(), QGpCompatibility::CompatDispersionData::misfit(), QGpCompatibility::CompatMultiModalData::reportToData(), QGpCompatibility::CompatMultiModalData::reportToDataWeight(), QGpCompatibility::CompatMultiModalCurves::reportToValues(), QGpCompatibility::CompatEllipticity::resetValues(), QGpCompatibility::CompatDispersion::resetValues(), QGpCompatibility::CompatAutocorrCurves::resetValues(), QGpCompatibility::CompatMultiModalData::setMean(), QGpCompatibility::CompatMultiModalData::setStddev(), QGpCompatibility::CompatMultiModalCurves::setValue(), QGpCompatibility::CompatMultiModalData::stddev(), QGpCompatibility::CompatMultiModalCurves::toStream(), QGpCompatibility::CompatMultiModalCurves::value(), QGpCompatibility::CompatMultiModalData::valuesToData(), and QGpCompatibility::CompatMultiModalCurves::valuesToReport().
| QList<DCModelInfo *> models |
| ReportWriter::Action reportAction = ReportWriter::Ask |
Referenced by main(), modeImportanceSampling(), modeNeighborhoodOptimization(), modeSnoopOptimization(), and reportMode().
| QString reportFile |
Referenced by main(), and reportMode().
| TargetList* targets = 0 |
Referenced by DinverCore::XMLDinverContext::all(), and main().
