GPLIB++
mtuadaptive.cpp
Go to the documentation of this file.
1 #include <fstream>
2 #include <iostream>
3 #include <vector>
4 #include <string>
5 #include <algorithm>
6 #include <boost/shared_ptr.hpp>
7 #include <boost/date_time/posix_time/posix_time.hpp>
8 #include <boost/program_options.hpp>
9 #include "LMSCanceller.h"
10 #include "RLSCanceller.h"
11 #include "AMRLSCanceller.h"
12 #include "TimeSeriesData.h"
13 #include "ApplyFilter.h"
14 #include "netcdfcpp.h"
15 #include "NeuralNetwork.h"
16 
17 
18 using namespace std;
19 using namespace gplib;
20 
21 string version = "$Id: mtuadaptive.cpp 1838 2010-03-04 16:19:34Z mmoorkamp $";
22 
23 void GetNNSetup(const size_t filterlength, const size_t hiddenlayers, const size_t ntimeseries,
24  NeuralNetwork::ttypeArray &NNLayers, double &NNmaxinit)
25  {
26  NeuralNetwork::ttypeVector typeVector(filterlength,
27  SigmoidalNeuron::bipolar); // we want filterlength number of bipolar neurons per hidden layer
28  for (size_t i = 0; i < hiddenlayers; ++i) //intialize the type array for the hidden layers
29  {
30  NNLayers.push_back(typeVector); // all layers are the same, so we copy the same vector there
31  }
32  typeVector.assign(1, SigmoidalNeuron::identity);
33  NNLayers.push_back(typeVector); // and then we add it to the type Array
34  }
35 
36 void WriteNetCDF(const string name, const boost::numeric::ublas::matrix<double,
37  boost::numeric::ublas::column_major> &Matrix)
38  {
39  NcFile mtrescdf(name.c_str(), NcFile::Replace);
40  NcDim* xd = mtrescdf.add_dim("x", Matrix.size1());
41  NcDim* yd = mtrescdf.add_dim("y", Matrix.size2());
42  NcVar* x = mtrescdf.add_var("x", ncFloat, xd);
43  NcVar* y = mtrescdf.add_var("y", ncFloat, yd);
44  NcVar* z = mtrescdf.add_var("z", ncFloat, xd, yd);
45  float *xvals = new float[xd->size()];
46  float *yvals = new float[yd->size()];
47  float *zvals = new float[xd->size() * yd->size()];
48  for (size_t i = 0; i < Matrix.size1(); ++i)
49  {
50  xvals[i] = i;
51  for (size_t j = 0; j < Matrix.size2(); ++j)
52  {
53  yvals[j] = j;
54  zvals[(Matrix.size1() - 1 - i) * Matrix.size2() + j] = Matrix(i, j);
55  }
56  }
57  x->put(xvals, xd->size());
58  y->put(yvals, yd->size());
59  z->put(zvals, z->edges());
60  }
61 
62 namespace po = boost::program_options;
63 
64 int main(int argc, char *argv[])
65  {
66  cout
67  << "This is mtuadaptive: Apply an adaptive filter (LMS or RLS) to cancel noise in MT time series"
68  << endl << endl;
69  cout
70  << " The program will ask for reference and input filename, further settings are read from 'mtuadaptive.conf' "
71  << endl;
72  cout
73  << " Output will be 2 Phoenix format files with ending '.clean' and '.eps' where the chosen channel is overwritten"
74  << endl << endl;
75  cout << " This is Version: " << version << endl << endl;
76 
77  int filterlength = 0, shift = 0, hiddenlayers;
78  double mu, lambda, alpha, delta;
79 
80  po::options_description desc("General options");
81  desc.add_options()("help", "produce help message")("filterlength",
82  po::value<int>(&filterlength)->default_value(10),
83  "The length of the adaptive filter")("shift",
84  po::value<int>(&shift)->default_value(0),
85  "The shift in samples between the time series")("mu",
86  po::value<double>(&mu)->default_value(1.0),
87  "Stepsize for LMS adaptive filter")("lambda",
88  po::value<double>(&lambda)->default_value(1.0), "")("alpha", po::value<
89  double>(&alpha)->default_value(1.0), "")("delta", po::value<double>(
90  &delta)->default_value(1.0), "")("hiddenlayers",
91  po::value<int>(&hiddenlayers)->default_value(1),
92  "The number of hiddenlayers for the neural network");
93 
94  po::variables_map vm;
95  po::store(po::parse_command_line(argc, argv, desc), vm);
96  po::notify(vm);
97 
98  if (vm.count("help"))
99  {
100  std::cout << desc << "\n";
101  return 1;
102  }
103 
104  TimeSeriesData InputData, ReferenceData;
105 
106  int index;
107  string tsfilename, noisefilename;
108 
109  int ntimeseries;
110  AdaptiveFilter *Filter;
111  int filterchoice = -1;
112  cout << "Which type of adaptive filter do you want ?" << endl;
113  cout << " 1: LMS" << endl;
114  cout << " 2: RLS" << endl;
115  cout << " 3: AMRLS" << endl;
116  cout << " 4: NN" << endl;
117  cin >> filterchoice;
118 
119  cout << "Reference Data: ";
120  cin >> noisefilename;
121  ReferenceData.GetData(noisefilename);
122  cout << "Component Index (Hx,Hy,Ex,Ey): ";
123  cin >> index;
124  cout << "Input Time Series Filename: ";
125  cin >> tsfilename;
126  InputData.GetData(tsfilename);
127  cout << "Input Start time: " << InputData.GetData().GetTime().front()
128  << endl;
129  cout << "Reference Start time: "
130  << ReferenceData.GetData().GetTime().front() << endl;
131  if (InputData.GetData().GetTime().front()
132  != ReferenceData.GetData().GetTime().front())
133  {
134  cerr << "Time series not synchronized !" << endl;
135  return 100;
136  }
137  int lengthdiff = ReferenceData.GetData().Size()
138  - InputData.GetData().Size();
139  if (lengthdiff > 0)
140  {
141  cout << "Removing " << lengthdiff
142  << " datapoints from Reference time-series." << endl;
143  ReferenceData.GetData().erase(ReferenceData.GetData().Size()
144  - lengthdiff, ReferenceData.GetData().Size());
145  }
146  if (lengthdiff < 0)
147  {
148  cout << "Removing " << lengthdiff
149  << " datapoints from Input time-series." << endl;
150  InputData.GetData().erase(InputData.GetData().Size() + lengthdiff,
151  InputData.GetData().Size());
152  }
153  cout << "Input End time: " << InputData.GetData().GetTime().back() << endl;
154  cout << "Reference End time: " << ReferenceData.GetData().GetTime().back()
155  << endl;
156 
157  TimeSeriesComponent *RefComp;
158  switch (index)
159  {
160  case 1:
161  RefComp = &ReferenceData.GetData().GetHx();
162  break;
163  case 2:
164  RefComp = &ReferenceData.GetData().GetHy();
165  break;
166  case 3:
167  RefComp = &ReferenceData.GetData().GetEx();
168  break;
169  case 4:
170  RefComp = &ReferenceData.GetData().GetEy();
171  break;
172  default:
173  cerr << "Component index not valid !";
174  return 100;
175  break;
176  }
177 
178  cout << "Number of input channels: ";
179  cin >> ntimeseries;
180 
181  NeuralNetwork::ttypeArray NNLayers;
182  double NNmaxinit;
183  switch (filterchoice)
184  {
185  case 1:
186  Filter = new LMSCanceller(filterlength, mu);
187  break;
188  case 2:
189  Filter = new RLSCanceller(filterlength,
190  delta, lambda);
191  break;
192  case 3:
193  Filter = new AMRLSCanceller(filterlength,
194  delta, lambda, alpha);
195  break;
196  case 4:
197  GetNNSetup(filterlength, hiddenlayers, ntimeseries, NNLayers, NNmaxinit);
198  Filter = new NeuralNetwork(filterlength, 1,
199  mu, NNLayers, NNmaxinit);
200  break;
201  default:
202  cerr << "Invalid number input, don't know what to do and will exit ! ";
203  return 100;
204  break;
205  }
206 
207  ApplyFilter Canceller(*Filter, true);
208  Canceller.AddReferenceChannel(*RefComp);
209  for (int i = 0; i < ntimeseries; ++i)
210  {
211  cout << "Component Index (Hx,Hy,Ex,Ey): ";
212  cin >> index;
213  switch (index)
214  {
215  case 1:
216  Canceller.AddInputChannel(InputData.GetData().GetHx());
217  break;
218  case 2:
219  Canceller.AddInputChannel(InputData.GetData().GetHy());
220  break;
221  case 3:
222  Canceller.AddInputChannel(InputData.GetData().GetEx());
223  break;
224  case 4:
225  Canceller.AddInputChannel(InputData.GetData().GetEy());
226  break;
227  default:
228  cerr << "Component index not valid !";
229  return 100;
230  break;
231  }
232  }
233 
234  Canceller.SetWeightSaveIntervall(1000);
235  Canceller.SetShift(shift);
236  Canceller.ShowProgress(true);
237  cout << " First iteration: " << endl << endl;
238 
239  Canceller.FilterData();
240 
241  cout << endl << endl << " Second iteration: " << endl << endl;
242 
243  Canceller.FilterData();
244 
245  copy(Canceller.GetOutChannels().front()->GetData().begin(),
246  Canceller.GetOutChannels().front()->GetData().end(),
247  RefComp->GetData().begin());
248  ReferenceData.WriteBack(noisefilename + ".clean");
249  ofstream cleanfile((noisefilename + ".clean").c_str());
250  copy(Canceller.GetOutChannels().front()->GetData().begin(),
251  Canceller.GetOutChannels().front()->GetData().end(), ostream_iterator<
252  double> (cleanfile, "\n"));
253  copy(Canceller.GetEpsValues().front().begin(),
254  Canceller.GetEpsValues().front().end(), RefComp->GetData().begin());
255 
256  ReferenceData.WriteBack(noisefilename + ".eps");
257  ofstream epsfile((noisefilename + ".eps").c_str());
258  copy(Canceller.GetEpsValues().front().begin(),
259  Canceller.GetEpsValues().front().end(), ostream_iterator<double> (
260  epsfile, "\n"));
261 
262  WriteNetCDF(noisefilename + ".weights.nc", Canceller.GetWeightHistory());
263  delete Filter;
264  }
const std::vector< std::vector< double > > & GetEpsValues()
Return the vector of channel approximation errors.
Definition: ApplyFilter.h:47
std::vector< double > & GetData()
Access for data vector, for ease of use and efficiency we return a reference.
void AddReferenceChannel(TimeSeriesComponent &Channel)
Add a reference channel to the filter, some AdaptiveFilter objects require only one reference...
Definition: ApplyFilter.cpp:45
Apply an adaptive filter to a time-series.
Definition: ApplyFilter.h:15
TimeSeries & GetData()
return a reference to the actual object stored in the pointer
TimeSeriesComponent & GetEx()
Definition: TimeSeries.h:47
gplib::rmat GetWeightHistory()
Return a matrix with the values of the weights at iterations specified by weightsaveintervall.
Definition: ApplyFilter.cpp:23
void ShowProgress(const bool what)
Do we want visual feedback of the progess on the screen, if yes draw a simple progress indicator in t...
Definition: ApplyFilter.h:30
void GetNNSetup(const size_t filterlength, const size_t hiddenlayers, const size_t ntimeseries, NeuralNetwork::ttypeArray &NNLayers, double &NNmaxinit)
Definition: mtuadaptive.cpp:23
void erase(const int startindex, const int endindex)
Erase data between startindex and endindex.
Definition: TimeSeries.cpp:59
void AddInputChannel(TimeSeriesComponent &Channel)
Add an input channel to the filter.
Definition: ApplyFilter.cpp:36
std::vector< ttypeVector > ttypeArray
Definition: NeuralNetwork.h:23
An implementation of the Recursive Least Squares filter with adptive memory as described in Hakin...
TimeSeriesComponent & GetHy()
Definition: TimeSeries.h:39
Implements a recursive least-squares adaptive filter, as described in Haykin, p. 443.
Definition: RLSCanceller.h:16
void SetWeightSaveIntervall(const int intervall)
Set the distance between iterations at which the weights are saved.
Definition: ApplyFilter.h:37
A generic base class for all types of adaptive filters.
TimeSeriesComponent is the base storage class for all types of time series data.
std::vector< SigmoidalNeuron::tneurontype > ttypeVector
Definition: NeuralNetwork.h:22
void FilterData()
Filter the input channels with the current settings.
Definition: ApplyFilter.cpp:57
int main(int argc, char *argv[])
Definition: mtuadaptive.cpp:64
void SetShift(const int theshift)
Set the shift between the input time series and the reference time series.
Definition: ApplyFilter.h:52
void WriteBack(std::string filename_base)
Write in the format it was originally read in.
Implements a LMS adaptive filter.
Definition: LMSCanceller.h:19
TimeSeriesData stores a pointer to the different components of magnetotelluric data and provides func...
ttimedata & GetTime()
Definition: TimeSeries.h:55
size_t Size()
Return the size of the time series, throws if one of the components has a different size...
Definition: TimeSeries.cpp:74
const std::vector< boost::shared_ptr< TimeSeriesComponent > > & GetOutChannels()
Return the vector of output channels.
Definition: ApplyFilter.h:42
string version
Definition: mtuadaptive.cpp:21
TimeSeriesComponent & GetEy()
Definition: TimeSeries.h:51
void WriteNetCDF(const string name, const boost::numeric::ublas::matrix< double, boost::numeric::ublas::column_major > &Matrix)
Definition: mtuadaptive.cpp:36
TimeSeriesComponent & GetHx()
Access function for Hx, returns reference for efficiency.
Definition: TimeSeries.h:35