mtnullspace.cpp

Go to the documentation of this file.

#include "C1DMTObjective.h"
#include "MTRecObjective.h"
#include "SeismicDataComp.h"
#include "GradFunc.h"
#include "CUniformRNG.h"
#include "convert.h"
#include <iostream>
#include <fstream>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/bindings/lapack/lapack.hpp>
#include <boost/numeric/bindings/traits/ublas_matrix.hpp>
#include <boost/numeric/bindings/traits/ublas_vector.hpp>
#include <boost/numeric/ublas/banded.hpp>
#include <algorithm>

using namespace std;
int main()
{
        const int nparam = 6;
        const int nsamples = 10;
        CMTStation MTData;
        SeismicDataComp RecData;
        RecData.GetData("test.rec");
        MTData.GetData("test.mtt"); //
        C1DMTObjective MTObjective(MTData);
        MTObjective.AppendFitParameters(&MTTensor::GetRhoxy,&MTTensor::GetdRhoxy,0.01);
        MTObjective.AppendFitParameters(&MTTensor::GetPhixy,&MTTensor::GetdPhixy,0.01);
        MTRecObjective CombObjective(MTData,RecData,0,0.0,4.0,0.001,0.07);
        
        ublas::vector<double> Model(nparam);
        Model(0) = 2;
        Model(1) = 0.1;
        Model(2) = 4;
        Model(3) = 10;
        Model(4) = 10;
        Model(5) = 10;
        //Model(4) = 3.4;
        //Model(5) = 4.0;
        //Model(6) = 2.8;
        //Model(7) = 3.2;
        double OldRMS = MTObjective.CalcPerformance(Model);
        MTObjective.CalcPerformance(Model);
        ublas::vector<double> OldData(MTObjective.GetSynthData());
        ModelAnalysis Analyser(MTObjective);
        Analyser.SetModel(Model);
        Analyser.SetRelativeSVDThreshold(1.0e-2);
        
        cout << "Singular Values : " << Analyser.GetSingularValues() << endl;
        cout << "Analysis R: " << Analyser.GetModelResolution() << endl;
        cout << "Analysis Null: " << Analyser.GetModelNullspace() << endl;
        
        ublas::vector<double> Projection(nparam), Perturb(nparam), NewModel(nparam);
        CUniformRNG Random;
        for (int sampleno = 0; sampleno < nsamples; ++sampleno)
        {
                for (int i = 0; i < nparam; ++i)
                        Perturb(i) = Random.GetNumber();
                Projection = prod(Analyser.GetModelNullspace(),Perturb);
                //cout << "Projection: " << Projection << endl;
                NewModel = Model+Projection;
                //cout << " Old Model: " << Model << endl;
                //cout << " New Model: " << NewModel << endl;
                double NewRMS = MTObjective.CalcPerformance(NewModel);
                ublas::vector<double> NewData(MTObjective.GetSynthData());
                MTObjective.WritePlot("out"+stringify(sampleno));
        }
        /*cout << "OldRMS: " << OldRMS << " NewRMS: " << NewRMS << endl;
        double difference = 0.0;
        for (int i = 0; i < NewData.size(); ++i)
                difference += pow((NewData(i) - OldData(i))/OldData(i),2);
        //cout << "Old Data: " << OldData << endl;
        //cout << "New Data: " << NewData << endl;
        ofstream diffile("differences");
        ofstream datafile("data");
        for (int i = 0; i < OldData.size(); ++i)
        {
                datafile << i << " " << OldData(i) << " " << NewData(i) << endl;
                diffile << (OldData(i) - NewData(i))/OldData(i) << endl;
        }
        cout << "Difference: " << sqrt(difference/NewData.size()) << endl;*/
}

---