GPLIB++
mtinvnet.cpp
Go to the documentation of this file.
1 #include "NeuralNetwork.h"
2 #include <iostream>
3 #include <vector>
4 #include <string>
5 #include <algorithm>
6 #include <iterator>
7 #include <cmath>
8 #include <fstream>
9 #include <boost/progress.hpp>
10 #include "UniformRNG.h"
11 #include "ApplyFilter.h"
12 #include "C1DMTSynthData.h"
13 
14 using namespace std;
15 using namespace gplib;
16 
17 const int frequenzen = 10;
18 const float T[frequenzen] =
19  { 1.06667, 2.15554, 4.33332, 8.66701, 16.000, 32.000, 64.0, 127.992, 256.016,
20  512.033 };
21 
22 int main()
23  {
24  UniformRNG Random;
25  NeuralNetwork::ttypeVector typeVector;
26  NeuralNetwork::ttypeArray typeArray;
27  C1DMTSynthData Synthetic; // create Object for Calculation of Synthetics
28  trealdata frequency;
29 
30  const int inlength = frequenzen;
31  for (int i = 0; i < frequenzen; ++i)
32  frequency.push_back(1. / T[i]);
33  Synthetic.SetFrequencies(frequency);
34 
35  const double maxinit = 0.5;
36  const int maxit = 500;
37  const int predit = 10;
38  const int hlayers = 2;
39  const int hlsize = 50;
40  const int modellayers = 2;
41  const int modelparms = 2;
42  const double mu = 0.1;
43  const double alpha = 0.1;
44  gplib::rvec input(inlength), desired(modelparms), output(modelparms);
45  vector<double> trainingerror;
46  const double minres = 1.0;
47  const double minthick = 5.0;
48  const double maxthick = 15.0;
49  const double noiselevel = 0.02;
50  trealdata resistivities(modellayers), thickness(modellayers);
51 
52  typeVector.assign(hlsize, SigmoidalNeuron::bipolar);
53  for (int i = 0; i < hlayers; ++i)
54  {
55  typeArray.push_back(typeVector);
56  }
57  typeVector.assign(modelparms, SigmoidalNeuron::bipolar);
58  typeArray.push_back(typeVector);
59 
60  NeuralNetwork Network(inlength, modelparms, mu, typeArray, maxinit);
61  Network.SetAlpha(alpha);
62  boost::progress_display progressbar(maxit);
63  for (int iteration = 0; iteration < maxit; ++iteration)
64  {
65  for (int i = 0; i < modellayers; ++i)
66  {
67  resistivities.at(i) = std::pow(10.0, minres + Random.GetNumber(0, 2));
68  thickness.at(i) = minthick; // + (maxthick - minthick) * Random.GetNumber();
69  }
70  Synthetic.SetResistivities(resistivities);
71  Synthetic.SetThicknesses(thickness);
72  Synthetic.GetData();
73  //cout << "Input: " << endl;
74  for (int i = 0; i < inlength; ++i)
75  {
76  double rho = Synthetic.GetMTData().at(i).GetRhoxy();
77  //input(i) = log10(rho + rho*noiselevel * Random.GetNumber(-1,1)) -2.0;
78  input(i) = log10(rho) - 2.0;
79  //cout << input(i) << endl;
80  }
81  //cout << endl;
82 
83  for (int i = 0; i < modellayers; ++i)
84  {
85  desired(i) = log10(resistivities.at(i)) - 2.0;
86 
87  //desired(2*(i+1)) = thickness.at(i)/100.0;
88  }
89  //desired(modellayers -1) = resistivities.at(modellayers-1)/100.0;
90  Network.CalcOutput(input, output);
91  Network.AdaptFilter(input, desired);
92  // cout << "Desired: " << endl;
93  // cout << desired << endl;
94  // cout << endl;
95  // cout << "Output: " << endl;
96  // cout << output << endl;
97  // cout << Network.GetEpsilon() << endl;
98 
99  trainingerror.push_back(inner_prod(Network.GetEpsilon(),
100  Network.GetEpsilon()));
101  ++progressbar;
102  }
103  ofstream terrfile("train.err");
104  copy(trainingerror.begin(), trainingerror.end(), ostream_iterator<double> (
105  terrfile, "\n"));
106  cout << endl << endl;
107  cout << "Finished training !" << endl;
108  for (int iteration = 0; iteration < predit; ++iteration)
109  {
110  for (int i = 0; i < modellayers; ++i)
111  {
112  resistivities.at(i) = std::pow(10.0, minres + Random.GetNumber(0, 2));
113  thickness.at(i) = minthick; // + (maxthick - minthick) * Random.GetNumber();
114  }
115  Synthetic.SetResistivities(resistivities);
116  Synthetic.SetThicknesses(thickness);
117  Synthetic.GetData();
118  for (int i = 0; i < inlength; ++i)
119  {
120  input(i) = log10(Synthetic.GetMTData().at(i).GetRhoxy()) - 2.0;
121  }
122  for (int i = 0; i < modellayers; ++i)
123  {
124  desired(i) = log10(resistivities.at(i)) - 2.0;
125  //desired(2*(i+1)) = thickness.at(i)/100.0;
126  }
127  //desired(modellayers -1) = resistivities.at(modellayers-1)/100.0;
128  Network.CalcOutput(input, output);
129  ++progressbar;
130  cout << "Model: ";
131  copy(desired.begin(), desired.end(), ostream_iterator<double> (cout,
132  " "));
133  cout << endl;
134  cout << "Output: ";
135  copy(output.begin(), output.end(), ostream_iterator<double> (cout, " "));
136  cout << endl;
137  }
138 
139  cout << endl << flush;
140  Synthetic.WriteAsMtt("in"); // Write out Mtt file
141  for (int i = 0; i < modellayers; ++i)
142  {
143  resistivities.at(i) = pow(10.0, output(i) + 2.0);
144  thickness.at(i) = minthick; //output(2*(i+1))*100.0;
145  }
146  //resistivities.at(modellayers-1) = output(modellayers-1)*100.0;
147  Synthetic.SetResistivities(resistivities);
148  Synthetic.SetThicknesses(thickness);
149  Synthetic.GetData();
150  Synthetic.WriteAsMtt("out"); // Write out Mtt file
151  Network.PrintTopology("test.dot");
152 
153  }
const gplib::rvec & GetEpsilon() const
Return the last estimation error.
Generate uniformly distributed random numbers, this is basically a wrapper for the boost random numbe...
Definition: UniformRNG.h:19
int main()
Definition: mtinvnet.cpp:22
void SetAlpha(const double a)
Set the momentum multiplier.
Definition: NeuralNetwork.h:45
std::vector< ttypeVector > ttypeArray
Definition: NeuralNetwork.h:23
virtual void GetData(const std::string filename)
read in data from file, determines format by ending
Definition: MTStation.cpp:597
virtual void AdaptFilter(const gplib::rvec &Input, const gplib::rvec &Desired)
Adapt the Filter with the current input and desired.
std::vector< SigmoidalNeuron::tneurontype > ttypeVector
Definition: NeuralNetwork.h:22
float GetNumber(const float low, const float high)
Return a random float between low and high.
Definition: UniformRNG.cpp:21
void SetThicknesses(const trealdata &thick)
Read only access to the vector of layer thicknesses for the 1D model from top to bottom in km...
void WriteAsMtt(const std::string filename)
Write data in goettingen .mtt format.
Definition: MTStation.cpp:681
const float T[frequenzen]
Definition: mtinvnet.cpp:18
Calculate synthetic MT data for a 1D model using Cagniard's algorithm.
void SetResistivities(const trealdata &res)
Read-write access to the vector of resistivities for the 1D model from top to bottom in Ohmm...
void PrintTopology(std::string filename)
Print the topology and weights of the network for plotting with the dot program.
const std::vector< MTTensor > & GetMTData() const
Get the full vector of Tensor elements read only.
Definition: MTStation.h:119
const int frequenzen
Definition: mtinvnet.cpp:17
void SetFrequencies(const trealdata &freqs)
Set the frequencies of the tensor elements, invalidates the previously stored impedance data...
Definition: MTStation.cpp:144