dataiterdecon.cpp

Go to the documentation of this file.

#include <iostream>
#include <algorithm>
#include <fstream>
#include "miscfunc.h"
#include "IterDecon.h"
#include "CSeismicDataComp.h"
#include "convert.h"
#include <numeric>
#include "VecMat.h"
using namespace std;
namespace ublas=boost::numeric::ublas;
int main()
{
         SeismicDataComp Radial,Vertical;
         std::string radfilename, verfilename;
         cout << "Radial Component: ";
         cin >> radfilename;
         cout << "Vertical Component: ";
         cin >> verfilename;
         
        Radial.GetData(radfilename);
        Vertical.GetData(verfilename);
        //Radial.GetData().insert(Radial.GetData().end(),Radial.GetData().size(),0.0);
        //Vertical.GetData().insert(Vertical.GetData().end(),Vertical.GetData().size(),0.0);
        std::vector<double> Filter(Radial.GetData().size());
    double variance = 10;
    int shift = 0;
    double normfactor = 1./(variance * sqrt(2.* acos(-1.)));
    for (unsigned int i = 0; i <= Filter.size()/2; ++i)
    {
        Filter.at(i) = normfactor * exp( -double((i-shift)*(i-shift))/(2. * variance * variance));
        Filter.at(Filter.size() -(i+1)) = normfactor * exp( -double((-i-1-shift)*(-i-1-shift))/(2. * variance * variance));
    }
    //std::ofstream filterfile("filter");
    //std::copy(Filter.begin(),Filter.end(),std::ostream_iterator<double>(filterfile,"\n"));
        Convolve(Vertical.GetData(),Filter,Vertical.GetData());
        Convolve(Radial.GetData(),Filter,Radial.GetData());
        //std::ofstream ts1filt("smpf.z"),ts2filt("smpf.r");
        //std::copy(Vertical.GetData().begin(),Vertical.GetData().end(),std::ostream_iterator<double>(ts1filt,"\n"));
        //std::copy(Radial.GetData().begin(),Radial.GetData().end(),std::ostream_iterator<double>(ts2filt,"\n"));
        
        IterDecon Decon(Vertical.GetData().size());
        const double power = inner_product(Radial.GetData().begin(),Radial.GetData().end(),Radial.GetData().begin(),0.0);
        double error = power;
        double previouserror = 1.;
        std::string iter;
        std::string name;
        gplib::rvec Output(Vertical.GetData().size());
        gplib::rvec Current(Radial.GetData().size());
        gplib::rvec Ver(Vertical.GetData().size());
        gplib::rvec Rad(Radial.GetData().size());
        copy(Vertical.GetData().begin(),Vertical.GetData().end(),Ver.begin());
        copy(Radial.GetData().begin(),Radial.GetData().end(),Rad.begin());
        copy(Radial.GetData().begin(),Radial.GetData().end(),Current.begin());
        //copy(Radial.GetData().begin(),Radial.GetData().end(),Output.begin());
        for (int i = 0; i < 50; ++i)
        {
                //std::cout << std::endl << "Iteration: " << i << std::endl;
         iter = stringify(i);
         Decon.AdaptFilter(Ver,Current);
         Decon.CalcOutput(Ver,Output);
      
         ublas::vector<double> Corr(Radial.GetData().size());
         Correl(Ver,Current,Corr);
          name = "corr"+iter;
         //ofstream corrfile(name.c_str());
         //std::copy(Corr.begin(),Corr.end(),std::ostream_iterator<double>(corrfile,"\n"));     
      
         Current = Rad-Output;
         error = inner_product(Current.begin(),Current.end(),Current.begin(),0.0) /power;
        
         /*name = "pred"+iter;
         ofstream predictfile(name.c_str());
         name = "diff"+iter;
         ofstream diffile(name.c_str());
         std::copy(Current.begin(),Current.end(),std::ostream_iterator<double>(diffile,"\n"));
                 std::copy(Output.begin(),Output.end(),std::ostream_iterator<double>(predictfile,"\n"));
                 diffile.close();
                 predictfile.close();*/ 
                 previouserror = error;
                 //std::ofstream outfile(("out"+iter).c_str());
                 //Decon.PrintWeights(outfile);
        }
        string outfilename = radfilename + ".decon.rec";
        std::cout << "Name: " << outfilename << std::endl;
        std::cout << "Error: " << error *100 << " Improv: " << (previouserror - error)*100 << std::endl;
        SeismicDataComp Receiver;
        Receiver.GetData().assign(Decon.GetWeights().size(),0.0);
        copy(Decon.GetWeights().begin(),Decon.GetWeights().end(),Receiver.GetData().begin());
        Convolve(Receiver.GetData(),Filter,Receiver.GetData());
        Receiver.GetData().erase(Receiver.GetData().begin()+Receiver.GetData().size()/2,Receiver.GetData().end());
        Receiver.SetDt(Radial.GetDt());
        Receiver.SetB(0.0); 
        
        Receiver.WriteAsSac(outfilename);
}