testiterdecon.cpp

Go to the documentation of this file.

#include <iostream>
#include <algorithm>
#include <fstream>
#include "miscfunc.h"
#include "IterDecon.h"
#include "convert.h"
#include <boost/numeric/ublas/vector.hpp>


namespace ublas=boost::numeric::ublas;
int main()
{
        int length = 500; 
    double dt = 0;
    double b = 0;
    
        std::ifstream radial("smp.r"), vertical("smp.z");
        std::ofstream predictfile,diffile,filterfile("filter");
        radial >> length >> dt >> b;
        //std::cout << "Length: " << length << " dt: " << dt << " b: " << b << std::endl;
        vertical >> length >> dt >> b;
        //std::cout << "Length: " << length << " dt: " << dt << " b: " << b << std::endl;
        gplib::rvec TimeSeries1(2*length);
        gplib::rvec TimeSeries2(2*length);
        gplib::rvec Output(2*length);
        for (int i = length/2; i < 1.5*length; ++i)
        {
                radial >> TimeSeries2(i);
                vertical >> TimeSeries1(i);
        }
        for (int i= 0; i < length/2; ++i)
        {
                TimeSeries1(i) = TimeSeries1(int(floor(i+1.5*length))) = 0;
                TimeSeries2(i) = TimeSeries2(int(floor(i+1.5*length))) = 0;
        }
        //std::cout << "lenght1: " << TimeSeries1.size() << " length2: " << TimeSeries2.size() << std::endl; 
        //double mean1=std::accumulate(TimeSeries1.begin(),TimeSeries1.end(),0.)/length;
        //double mean2=std::accumulate(TimeSeries2.begin(),TimeSeries2.end(),0.)/length;
        //std::cout << "Mean1: " << mean1 << " Mean2: " << mean2 << std::endl;
        //std::transform(TimeSeries2.begin(),TimeSeries2.end(),TimeSeries2.begin(),bind2nd(std::minus<double>(),mean2));
        //std::transform(TimeSeries1.begin(),TimeSeries1.end(),TimeSeries1.begin(),bind2nd(std::minus<double>(),mean1));
        gplib::rvec Filter(TimeSeries1.size());
    double variance = 2;
    int shift = 0;
    double normfactor = 1./(variance * sqrt(2.* acos(-1.)));
    for (unsigned int i = 0; i <= Filter.size()/2; ++i)
    {
        Filter(i) = normfactor * exp( -pow((i-shift),2)/(2. * variance * variance));
        Filter(Filter.size() -i-1) = normfactor * exp( -pow((-i-1-shift),2)/(2. * variance * variance));
    }
    std::copy(Filter.begin(),Filter.end(),std::ostream_iterator<double>(filterfile,"\n"));
        Convolve(TimeSeries1,Filter,TimeSeries1);
        Convolve(TimeSeries2,Filter,TimeSeries2);
        std::ofstream ts1filt("smpf.z"),ts2filt("smpf.r");
        std::copy(TimeSeries1.begin(),TimeSeries1.end(),std::ostream_iterator<double>(ts1filt,"\n"));
        std::copy(TimeSeries2.begin(),TimeSeries2.end(),std::ostream_iterator<double>(ts2filt,"\n"));
        gplib::rvec Current(TimeSeries2);
        IterDecon Decon(TimeSeries1.size());
        const double power = ublas::prec_inner_prod(TimeSeries2,TimeSeries2);
        double error = power;
        double previouserror = 1.;
        std::string iter;
        std::string name;
        for (int i = 0; i < 10; ++i)
        {
                std::cout << std::endl << "Iteration: " << i << std::endl;
         iter = stringify(i);
         Decon.AdaptFilter(TimeSeries1,Current);
         Decon.CalcOutput(TimeSeries1,Output);
      
         gplib::rvec Corr(TimeSeries2);
         Correl(TimeSeries1,Current,Corr);
            name = "corr"+iter;
         predictfile.open(name.c_str());
         std::copy(Corr.begin(),Corr.end(),std::ostream_iterator<double>(predictfile,"\n"));
        predictfile.close();    
      
         Current = TimeSeries2-Output;
         error = ublas::prec_inner_prod(Current,Current) /power;
         std::cout << "Error: " << error *100 << " Improv: " << (previouserror - error)*100 << std::endl;
         name = "pred"+iter;
         predictfile.open(name.c_str());
         name = "diff"+iter;
         diffile.open(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);
        }
}

---