mtuadaptive.cpp

Go to the documentation of this file.
00001 #include <fstream>
00002 #include <iostream>
00003 #include <vector>
00004 #include <string>
00005 #include <algorithm>
00006 #include <boost/shared_ptr.hpp>
00007 #include <boost/date_time/posix_time/posix_time.hpp>
00008 #include <boost/program_options.hpp>
00009 #include "LMSCanceller.h"
00010 #include "RLSCanceller.h"
00011 #include "AMRLSCanceller.h"
00012 #include "TimeSeriesData.h"
00013 #include "ApplyFilter.h"
00014 #include "netcdfcpp.h"
00015 #include "NeuralNetwork.h"
00016 
00017 
00018 using namespace std;
00019 using namespace gplib;
00020 
00021 string version = "$Id: mtuadaptive.cpp 1838 2010-03-04 16:19:34Z mmoorkamp $";
00022 
00023 void GetNNSetup(const size_t filterlength, const size_t hiddenlayers, const size_t ntimeseries,
00024     NeuralNetwork::ttypeArray &NNLayers, double &NNmaxinit)
00025   {
00026     NeuralNetwork::ttypeVector typeVector(filterlength,
00027         SigmoidalNeuron::bipolar); // we want filterlength number of bipolar neurons per hidden layer
00028     for (size_t i = 0; i < hiddenlayers; ++i) //intialize the type array for the hidden layers
00029       {
00030         NNLayers.push_back(typeVector); // all layers are the same, so we copy the same vector there
00031       }
00032     typeVector.assign(1, SigmoidalNeuron::identity);
00033     NNLayers.push_back(typeVector); // and then we add it to the type Array
00034   }
00035 
00036 void WriteNetCDF(const string name, const boost::numeric::ublas::matrix<double,
00037     boost::numeric::ublas::column_major> &Matrix)
00038   {
00039     NcFile mtrescdf(name.c_str(), NcFile::Replace);
00040     NcDim* xd = mtrescdf.add_dim("x", Matrix.size1());
00041     NcDim* yd = mtrescdf.add_dim("y", Matrix.size2());
00042     NcVar* x = mtrescdf.add_var("x", ncFloat, xd);
00043     NcVar* y = mtrescdf.add_var("y", ncFloat, yd);
00044     NcVar* z = mtrescdf.add_var("z", ncFloat, xd, yd);
00045     float *xvals = new float[xd->size()];
00046     float *yvals = new float[yd->size()];
00047     float *zvals = new float[xd->size() * yd->size()];
00048     for (size_t i = 0; i < Matrix.size1(); ++i)
00049       {
00050         xvals[i] = i;
00051         for (size_t j = 0; j < Matrix.size2(); ++j)
00052           {
00053             yvals[j] = j;
00054             zvals[(Matrix.size1() - 1 - i) * Matrix.size2() + j] = Matrix(i, j);
00055           }
00056       }
00057     x->put(xvals, xd->size());
00058     y->put(yvals, yd->size());
00059     z->put(zvals, z->edges());
00060   }
00061 
00062 namespace po = boost::program_options;
00063 
00064 int main(int argc, char *argv[])
00065   {
00066     cout
00067         << "This is mtuadaptive: Apply an adaptive filter (LMS or RLS) to cancel noise in MT time series"
00068         << endl << endl;
00069     cout
00070         << " The program will ask for reference and input filename, further settings are read from 'mtuadaptive.conf' "
00071         << endl;
00072     cout
00073         << " Output will be 2 Phoenix format files with ending '.clean' and '.eps' where the chosen channel is overwritten"
00074         << endl << endl;
00075     cout << " This is Version: " << version << endl << endl;
00076 
00077     int filterlength = 0, shift = 0, hiddenlayers;
00078     double mu, lambda, alpha, delta;
00079 
00080     po::options_description desc("General options");
00081     desc.add_options()("help", "produce help message")("filterlength",
00082         po::value<int>(&filterlength)->default_value(10),
00083         "The length of the adaptive filter")("shift",
00084         po::value<int>(&shift)->default_value(0),
00085         "The shift in samples between the time series")("mu",
00086         po::value<double>(&mu)->default_value(1.0),
00087         "Stepsize for LMS adaptive filter")("lambda",
00088         po::value<double>(&lambda)->default_value(1.0), "")("alpha", po::value<
00089         double>(&alpha)->default_value(1.0), "")("delta", po::value<double>(
00090         &delta)->default_value(1.0), "")("hiddenlayers",
00091             po::value<int>(&hiddenlayers)->default_value(1),
00092             "The number of hiddenlayers for the neural network");
00093 
00094     po::variables_map vm;
00095     po::store(po::parse_command_line(argc, argv, desc), vm);
00096     po::notify(vm);
00097 
00098     if (vm.count("help"))
00099       {
00100         std::cout << desc << "\n";
00101         return 1;
00102       }
00103 
00104     TimeSeriesData InputData, ReferenceData;
00105 
00106     int index;
00107     string tsfilename, noisefilename;
00108 
00109     int ntimeseries;
00110     AdaptiveFilter *Filter;
00111     int filterchoice = -1;
00112     cout << "Which type of adaptive filter do you want ?" << endl;
00113     cout << " 1: LMS" << endl;
00114     cout << " 2: RLS" << endl;
00115     cout << " 3: AMRLS" << endl;
00116     cout << " 4: NN" << endl;
00117     cin >> filterchoice;
00118 
00119     cout << "Reference Data: ";
00120     cin >> noisefilename;
00121     ReferenceData.GetData(noisefilename);
00122     cout << "Component Index (Hx,Hy,Ex,Ey): ";
00123     cin >> index;
00124     cout << "Input Time Series Filename: ";
00125     cin >> tsfilename;
00126     InputData.GetData(tsfilename);
00127     cout << "Input Start time: " << InputData.GetData().GetTime().front()
00128         << endl;
00129     cout << "Reference Start time: "
00130         << ReferenceData.GetData().GetTime().front() << endl;
00131     if (InputData.GetData().GetTime().front()
00132         != ReferenceData.GetData().GetTime().front())
00133       {
00134         cerr << "Time series not synchronized !" << endl;
00135         return 100;
00136       }
00137     int lengthdiff = ReferenceData.GetData().Size()
00138         - InputData.GetData().Size();
00139     if (lengthdiff > 0)
00140       {
00141         cout << "Removing " << lengthdiff
00142             << " datapoints from Reference time-series." << endl;
00143         ReferenceData.GetData().erase(ReferenceData.GetData().Size()
00144             - lengthdiff, ReferenceData.GetData().Size());
00145       }
00146     if (lengthdiff < 0)
00147       {
00148         cout << "Removing " << lengthdiff
00149             << " datapoints from Input time-series." << endl;
00150         InputData.GetData().erase(InputData.GetData().Size() + lengthdiff,
00151             InputData.GetData().Size());
00152       }
00153     cout << "Input End time: " << InputData.GetData().GetTime().back() << endl;
00154     cout << "Reference End time: " << ReferenceData.GetData().GetTime().back()
00155         << endl;
00156 
00157     TimeSeriesComponent *RefComp;
00158     switch (index)
00159       {
00160     case 1:
00161       RefComp = &ReferenceData.GetData().GetHx();
00162       break;
00163     case 2:
00164       RefComp = &ReferenceData.GetData().GetHy();
00165       break;
00166     case 3:
00167       RefComp = &ReferenceData.GetData().GetEx();
00168       break;
00169     case 4:
00170       RefComp = &ReferenceData.GetData().GetEy();
00171       break;
00172     default:
00173       cerr << "Component index not valid !";
00174       return 100;
00175       break;
00176       }
00177 
00178     cout << "Number of input channels: ";
00179     cin >> ntimeseries;
00180 
00181     NeuralNetwork::ttypeArray NNLayers;
00182     double NNmaxinit;
00183     switch (filterchoice)
00184       {
00185     case 1:
00186       Filter = new LMSCanceller(filterlength, mu);
00187       break;
00188     case 2:
00189       Filter = new RLSCanceller(filterlength,
00190           delta, lambda);
00191       break;
00192     case 3:
00193       Filter = new AMRLSCanceller(filterlength,
00194           delta, lambda, alpha);
00195       break;
00196     case 4:
00197       GetNNSetup(filterlength, hiddenlayers, ntimeseries, NNLayers, NNmaxinit);
00198       Filter = new NeuralNetwork(filterlength, 1,
00199           mu, NNLayers, NNmaxinit);
00200       break;
00201     default:
00202       cerr << "Invalid number input, don't know what to do and will exit ! ";
00203       return 100;
00204       break;
00205       }
00206 
00207     ApplyFilter Canceller(*Filter, true);
00208     Canceller.AddReferenceChannel(*RefComp);
00209     for (int i = 0; i < ntimeseries; ++i)
00210       {
00211         cout << "Component Index (Hx,Hy,Ex,Ey): ";
00212         cin >> index;
00213         switch (index)
00214           {
00215         case 1:
00216           Canceller.AddInputChannel(InputData.GetData().GetHx());
00217           break;
00218         case 2:
00219           Canceller.AddInputChannel(InputData.GetData().GetHy());
00220           break;
00221         case 3:
00222           Canceller.AddInputChannel(InputData.GetData().GetEx());
00223           break;
00224         case 4:
00225           Canceller.AddInputChannel(InputData.GetData().GetEy());
00226           break;
00227         default:
00228           cerr << "Component index not valid !";
00229           return 100;
00230           break;
00231           }
00232       }
00233 
00234     Canceller.SetWeightSaveIntervall(1000);
00235     Canceller.SetShift(shift);
00236     Canceller.ShowProgress(true);
00237     cout << " First iteration: " << endl << endl;
00238 
00239     Canceller.FilterData();
00240 
00241     cout << endl << endl << " Second iteration: " << endl << endl;
00242 
00243     Canceller.FilterData();
00244 
00245     copy(Canceller.GetOutChannels().front()->GetData().begin(),
00246         Canceller.GetOutChannels().front()->GetData().end(),
00247         RefComp->GetData().begin());
00248     ReferenceData.WriteBack(noisefilename + ".clean");
00249     ofstream cleanfile((noisefilename + ".clean").c_str());
00250     copy(Canceller.GetOutChannels().front()->GetData().begin(),
00251         Canceller.GetOutChannels().front()->GetData().end(), ostream_iterator<
00252             double> (cleanfile, "\n"));
00253     copy(Canceller.GetEpsValues().front().begin(),
00254         Canceller.GetEpsValues().front().end(), RefComp->GetData().begin());
00255 
00256     ReferenceData.WriteBack(noisefilename + ".eps");
00257     ofstream epsfile((noisefilename + ".eps").c_str());
00258     copy(Canceller.GetEpsValues().front().begin(),
00259         Canceller.GetEpsValues().front().end(), ostream_iterator<double> (
00260             epsfile, "\n"));
00261 
00262     WriteNetCDF(noisefilename + ".weights.nc", Canceller.GetWeightHistory());
00263     delete Filter;
00264   }

Generated on Tue May 4 16:52:14 2010 for GPLIB++ by  doxygen 1.5.8