mtinvnet.cpp

Go to the documentation of this file.

#include "NeuralNetwork.h"
#include <iostream>
#include <vector>
#include <string>
#include <algorithm>
#include <iterator>
#include <cmath>
#include <fstream>
#include <boost/progress.hpp>
#include "UniformRNG.h"
#include "ApplyFilter.h"
#include "C1DMTSynthData.h"

using namespace std;
using namespace gplib;

const int frequenzen = 10;
const float T[frequenzen] =
  { 1.06667, 2.15554, 4.33332, 8.66701, 16.000, 32.000, 64.0, 127.992, 256.016,
      512.033 };

int main()
  {
    UniformRNG Random;
    NeuralNetwork::ttypeVector typeVector;
    NeuralNetwork::ttypeArray typeArray;
    C1DMTSynthData Synthetic; // create Object for Calculation of Synthetics
    trealdata frequency;

    const int inlength = frequenzen;
    for (int i = 0; i < frequenzen; ++i)
      frequency.push_back(1. / T[i]);
    Synthetic.SetFrequencies(frequency);

    const double maxinit = 0.5;
    const int maxit = 500;
    const int predit = 10;
    const int hlayers = 2;
    const int hlsize = 50;
    const int modellayers = 2;
    const int modelparms = 2;
    const double mu = 0.1;
    const double alpha = 0.1;
    gplib::rvec input(inlength), desired(modelparms), output(modelparms);
    vector<double> trainingerror;
    const double minres = 1.0;
    const double minthick = 5.0;
    const double maxthick = 15.0;
    const double noiselevel = 0.02;
    trealdata resistivities(modellayers), thickness(modellayers);

    typeVector.assign(hlsize, SigmoidalNeuron::bipolar);
    for (int i = 0; i < hlayers; ++i)
      {
        typeArray.push_back(typeVector);
      }
    typeVector.assign(modelparms, SigmoidalNeuron::bipolar);
    typeArray.push_back(typeVector);

    NeuralNetwork Network(inlength, modelparms, mu, typeArray, maxinit);
    Network.SetAlpha(alpha);
    boost::progress_display progressbar(maxit);
    for (int iteration = 0; iteration < maxit; ++iteration)
      {
        for (int i = 0; i < modellayers; ++i)
          {
            resistivities.at(i) = std::pow(10.0, minres + Random.GetNumber(0, 2));
            thickness.at(i) = minthick; // + (maxthick - minthick) * Random.GetNumber();
          }
        Synthetic.SetResistivities(resistivities);
        Synthetic.SetThicknesses(thickness);
        Synthetic.GetData();
        //cout << "Input: " << endl;
        for (int i = 0; i < inlength; ++i)
          {
            double rho = Synthetic.GetMTData().at(i).GetRhoxy();
            //input(i) = log10(rho + rho*noiselevel * Random.GetNumber(-1,1)) -2.0;
            input(i) = log10(rho) - 2.0;
            //cout << input(i) << endl;
          }
        //cout  << endl;

        for (int i = 0; i < modellayers; ++i)
          {
            desired(i) = log10(resistivities.at(i)) - 2.0;

            //desired(2*(i+1)) = thickness.at(i)/100.0;
          }
        //desired(modellayers -1) = resistivities.at(modellayers-1)/100.0;
        Network.CalcOutput(input, output);
        Network.AdaptFilter(input, desired);
        //              cout << "Desired: " << endl;
        //              cout << desired << endl;
        //              cout  << endl;
        //              cout << "Output: " << endl;
        //              cout << output << endl;
        //              cout << Network.GetEpsilon() << endl;

        trainingerror.push_back(inner_prod(Network.GetEpsilon(),
            Network.GetEpsilon()));
        ++progressbar;
      }
    ofstream terrfile("train.err");
    copy(trainingerror.begin(), trainingerror.end(), ostream_iterator<double> (
        terrfile, "\n"));
    cout << endl << endl;
    cout << "Finished training !" << endl;
    for (int iteration = 0; iteration < predit; ++iteration)
      {
        for (int i = 0; i < modellayers; ++i)
          {
            resistivities.at(i) = std::pow(10.0, minres + Random.GetNumber(0, 2));
            thickness.at(i) = minthick; // + (maxthick - minthick) * Random.GetNumber();
          }
        Synthetic.SetResistivities(resistivities);
        Synthetic.SetThicknesses(thickness);
        Synthetic.GetData();
        for (int i = 0; i < inlength; ++i)
          {
            input(i) = log10(Synthetic.GetMTData().at(i).GetRhoxy()) - 2.0;
          }
        for (int i = 0; i < modellayers; ++i)
          {
            desired(i) = log10(resistivities.at(i)) - 2.0;
            //desired(2*(i+1)) = thickness.at(i)/100.0;
          }
        //desired(modellayers -1) = resistivities.at(modellayers-1)/100.0;
        Network.CalcOutput(input, output);
        ++progressbar;
        cout << "Model: ";
        copy(desired.begin(), desired.end(), ostream_iterator<double> (cout,
            " "));
        cout << endl;
        cout << "Output: ";
        copy(output.begin(), output.end(), ostream_iterator<double> (cout, " "));
        cout << endl;
      }

    cout << endl << flush;
    Synthetic.WriteAsMtt("in"); // Write out Mtt file
    for (int i = 0; i < modellayers; ++i)
      {
        resistivities.at(i) = pow(10.0, output(i) + 2.0);
        thickness.at(i) = minthick; //output(2*(i+1))*100.0;
      }
    //resistivities.at(modellayers-1) = output(modellayers-1)*100.0;
    Synthetic.SetResistivities(resistivities);
    Synthetic.SetThicknesses(thickness);
    Synthetic.GetData();
    Synthetic.WriteAsMtt("out"); // Write out Mtt file
    Network.PrintTopology("test.dot");

  }