[pktools] 117/375: added pkfilterascii.cc and pkregression_nn.cc

Bas Couwenberg sebastic at xs4all.nl
Wed Dec 3 21:54:05 UTC 2014


This is an automated email from the git hooks/post-receive script.

sebastic-guest pushed a commit to branch upstream-master
in repository pktools.

commit e0883887097346b390b293faf1abd86c533a7cf9
Author: Pieter Kempeneers <kempenep at gmail.com>
Date:   Mon Jun 17 17:45:08 2013 +0200

    added pkfilterascii.cc and pkregression_nn.cc
---
 src/algorithms/Filter.cc          |  20 +++
 src/algorithms/Filter.h           |  53 ++++--
 src/algorithms/myfann_cpp.h       | 136 +++++++++++++--
 src/apps/Makefile.am              |  12 +-
 src/apps/pkclassify_nn.cc         |   1 -
 src/apps/pkfilter.cc              |   4 +-
 src/apps/pkfilterascii.cc         | 233 +++++++++++++++++++++++++
 src/apps/pkfs_nn.cc               |   1 -
 src/apps/pkregression_nn.cc       | 346 ++++++++++++++++++++++++++++++++++++++
 src/base/Vector2d.h               |   8 +
 src/fileclasses/FileReaderAscii.h |  40 +++--
 11 files changed, 806 insertions(+), 48 deletions(-)

diff --git a/src/algorithms/Filter.cc b/src/algorithms/Filter.cc
index 21009de..fc7cb2d 100644
--- a/src/algorithms/Filter.cc
+++ b/src/algorithms/Filter.cc
@@ -38,6 +38,26 @@ void filter::Filter::setTaps(const vector<double> &taps)
   assert(m_taps.size()%2);
 }
 
+// void filter::Filter::dwtForward(std::vector<double>& data, const std::string& wavelet_type, int family){
+//   int nsize=data.size();
+//   assert(nsize);
+//   gsl_wavelet *w;
+//   gsl_wavelet_workspace *work;
+//   w=gsl_wavelet_alloc(getWaveletType(wavelet_type),family);
+//   work=gsl_wavelet_workspace_alloc(nsize);
+//   gsl_wavelet_transform_forward(w,&(data[0]),1,nsize,work);
+// }
+
+// void filter::Filter::dwtInverse(std::vector<double>& data, const std::string& wavelet_type, int family){
+//   int nsize=data.size();
+//   assert(nsize);
+//   gsl_wavelet *w;
+//   gsl_wavelet_workspace *work;
+//   w=gsl_wavelet_alloc(getWaveletType(wavelet_type),family);
+//   work=gsl_wavelet_workspace_alloc(nsize);
+//   gsl_wavelet_transform_inverse(w,&(data[0]),1,nsize,work);
+// }
+
 void filter::Filter::morphology(const ImgReaderGdal& input, ImgWriterGdal& output, const std::string& method, int dim, short down, int offset)
 {
   Vector2d<double> lineInput(input.nrOfBand(),input.nrOfCol());
diff --git a/src/algorithms/Filter.h b/src/algorithms/Filter.h
index 3da2601..ef89a01 100644
--- a/src/algorithms/Filter.h
+++ b/src/algorithms/Filter.h
@@ -22,6 +22,7 @@ along with pktools.  If not, see <http://www.gnu.org/licenses/>.
 
 #include <vector>
 #include <iostream>
+#include <gsl/gsl_wavelet.h>
 #include "StatFactory.h"
 #include "imageclasses/ImgReaderGdal.h"
 #include "imageclasses/ImgWriterGdal.h"
@@ -30,7 +31,7 @@ using namespace std;
 namespace filter
 {
   
-  enum FILTER_TYPE { median=0, var=1 , min=2, max=3, sum=4, mean=5, minmax=6, dilate=7, erode=8, close=9, open=10, homog=11, sobelx=12, sobely=13, sobelxy=14, sobelyx=-14, smooth=15, density=16, majority=17, mixed=18, smoothnodata=19, threshold=20, ismin=21, ismax=22, heterog=23, order=24, stdev=25};
+  enum FILTER_TYPE { median=0, var=1 , min=2, max=3, sum=4, mean=5, minmax=6, dilate=7, erode=8, close=9, open=10, homog=11, sobelx=12, sobely=13, sobelxy=14, sobelyx=-14, smooth=15, density=16, majority=17, mixed=18, smoothnodata=19, threshold=20, ismin=21, ismax=22, heterog=23, order=24, stdev=25, dwtForward=26, dwtInverse=27};
 
 class Filter
 {
@@ -38,9 +39,14 @@ public:
   Filter(void);
   Filter(const vector<double> &taps);
   virtual ~Filter(){};
-  FILTER_TYPE getFilterType(const std::string filterType){
+  static const gsl_wavelet_type* getWaveletType(const std::string waveletType){
+    std::map<std::string, const gsl_wavelet_type* > m_waveletMap;
+    initWaveletMap(m_waveletMap);
+    return m_waveletMap[waveletType];
+  }
+  static FILTER_TYPE getFilterType(const std::string filterType){
     std::map<std::string, FILTER_TYPE> m_filterMap;
-    initMap(m_filterMap);
+    initFilterMap(m_filterMap);
     return m_filterMap[filterType];
   };
   void setTaps(const vector<double> &taps);
@@ -53,19 +59,28 @@ public:
   void doit(const ImgReaderGdal& input, ImgWriterGdal& output, short down=1, int offset=0);
 
   template<class T> double applySrf(const vector<double> &wavelengthIn, const vector<T>& input, const Vector2d<double>& srf, const std::string& interpolationType, T& output, double delta=1.0, bool normalize=false, bool verbose=false);
-  template<class T> double applySrf(const vector<double> &wavelengthIn, const Vector2d<T>& input, const Vector2d<double>& srf, const std::string& interpolationType, vector<T>& output, double delta=1.0, bool normalize=false, int down=1, bool verbose=false);
+  template<class T> double applySrf(const vector<double> &wavelengthIn, const Vector2d<T>& input, const Vector2d<double>& srf, const std::string& interpolationType, vector<T>& output, double delta=1.0, bool normalize=false, int down=1, bool transposeInput=false, bool verbose=false);
 
-  // void applySrf(const vector<double> &wavelengthIn, const ImgReaderGdal& input, const vector< Vector2d<double> > &srf, const std::string& interpolationType, ImgWriterGdal& output, bool verbose=false);
   template<class T> void applyFwhm(const vector<double> &wavelengthIn, const vector<T>& input, const vector<double> &wavelengthOut, const vector<double> &fwhm, const std::string& interpolationType, vector<T>& output, bool verbose=false);
   template<class T> void applyFwhm(const vector<double> &wavelengthIn, const Vector2d<T>& input, const vector<double> &wavelengthOut, const vector<double> &fwhm, const std::string& interpolationType, Vector2d<T>& output, int down=1, bool verbose=false);
-  // void applyFwhm(const vector<double> &wavelengthIn, const ImgReaderGdal& input, const vector<double> &wavelengthOut, const vector<double> &fwhm, const std::string& interpolationType, ImgWriterGdal& output, bool verbose=false);
-// int fir(double* input, int nbandIn, vector<double>& output, int startBand, const string& wavelength, const string& fwhm, bool verbose);
-// int fir(const vector<double>&input, vector<double>& output, int startBand, double fwhm, int ntaps, int down, int offset, bool verbose);
-// int fir(double* input, int nbandIn, vector<double>& output, int startBand, double fwhm, int ntaps, int down, int offset, bool verbose);
+  // void dwtForward(std::vector<double>& data, const std::string& wavelet_type, int family);
+  // void dwtInverse(std::vector<double>& data, const std::string& wavelet_type, int family);
 
 private:
-  static void initMap(std::map<std::string, FILTER_TYPE>& m_filterMap){
-    //initialize selMap
+  static void initWaveletMap(std::map<std::string, const gsl_wavelet_type*> m_waveletMap){
+    //initialize Map
+    m_waveletMap["daubechies"]=gsl_wavelet_daubechies;
+    m_waveletMap["daubechies_centered"]=gsl_wavelet_daubechies_centered;
+    m_waveletMap["haar"]=gsl_wavelet_haar;
+    m_waveletMap["haar_centered"]=gsl_wavelet_haar_centered;
+    m_waveletMap["bspline"]=gsl_wavelet_bspline;
+    m_waveletMap["bspline_centered"]=gsl_wavelet_bspline_centered;
+  }
+
+  static void initFilterMap(std::map<std::string, FILTER_TYPE>& m_filterMap){
+    //initialize Map
+    m_filterMap["dwtForward"]=filter::dwtForward;
+    m_filterMap["dwtInverse"]=filter::dwtInverse;
     m_filterMap["stdev"]=filter::stdev;
     m_filterMap["var"]=filter::var;
     m_filterMap["min"]=filter::min;
@@ -175,13 +190,13 @@ private:
   return(srf[0][maxIndex]);
 }
 
-//input[band][sample], output[sample]
+//input[band][sample], output[sample] (if !transposeInput)
 //returns wavelength for which srf is maximum
-  template<class T> double Filter::applySrf(const vector<double> &wavelengthIn, const Vector2d<T>& input, const Vector2d<double>& srf, const std::string& interpolationType, vector<T>& output, double delta, bool normalize, int down, bool verbose)
+  template<class T> double Filter::applySrf(const vector<double> &wavelengthIn, const Vector2d<T>& input, const Vector2d<double>& srf, const std::string& interpolationType, vector<T>& output, double delta, bool normalize, int down, bool transposeInput, bool verbose)
 {  
   assert(srf.size()==2);//[0]: wavelength, [1]: response function
   int nband=srf[0].size(); 
-  unsigned int nsample=input[0].size();
+  unsigned int nsample=(transposeInput)? input.size():input[0].size();
   output.resize((nsample+down-1)/down);
   double start=floor(wavelengthIn[0]);
   double end=ceil(wavelengthIn.back());
@@ -226,7 +241,10 @@ private:
     if((isample+1+down/2)%down)
       continue;
     vector<T> inputValues;
-    input.selectCol(isample,inputValues);
+    if(transposeInput)
+      inputValues=input[isample];
+    else
+      input.selectCol(isample,inputValues);
     assert(wavelengthIn.size()==inputValues.size());
     vector<double> input_fine;
     vector<double> product(wavelength_fine.size());
@@ -238,7 +256,6 @@ private:
     assert(input_fine.size()==srf_fine.size());
     assert(input_fine.size()==wavelength_fine.size());
     stat.initSpline(splineOut,&(wavelength_fine[0]),&(product[0]),wavelength_fine.size());
-    //hiero
     if(normalize)
       output[isample/down]=gsl_spline_eval_integ(splineOut,start,end,accOut)/norm;
     else
@@ -315,8 +332,8 @@ template<class T> void Filter::applyFwhm(const vector<double> &wavelengthIn, con
   for(double win=floor(wavelengthIn[0]);win<=ceil(wavelengthIn.back());win+=delta)
     wavelength_fine.push_back(win);
   assert(wavelengthOut.size()==fwhm.size());
-  assert(wavelengthIn[0]<wavelengthOut[0]);
-  assert(wavelengthIn.back()>wavelengthOut.back());
+  assert(wavelengthIn[0]<=wavelengthOut[0]);
+  assert(wavelengthIn.back()>=wavelengthOut.back());
   if(verbose){
     for(int index=0;index<wavelength_fine.size();++index)
       std::cout << " " << wavelength_fine[index];
diff --git a/src/algorithms/myfann_cpp.h b/src/algorithms/myfann_cpp.h
index d86958f..1308e22 100644
--- a/src/algorithms/myfann_cpp.h
+++ b/src/algorithms/myfann_cpp.h
@@ -645,7 +645,7 @@ namespace FANN
            Parameters:
              num_data      - The number of training data
              num_input     - The number of inputs per training data
-             num_output    - The number of ouputs per training data
+             num_output    - The number of outputs per training data
              input      - The set of inputs (a pointer to an array of pointers to arrays of floating point data)
              output     - The set of desired outputs (a pointer to an array of pointers to arrays of floating point data)
 
@@ -850,13 +850,13 @@ public:
            Parameters:
              num_data      - The number of training data
              num_input     - The number of inputs per training data
-             num_output    - The number of ouputs per training data
+             num_output    - The number of outputs per training data
              user_function - The user suplied function
 
            Parameters for the user function:
              num        - The number of the training data set
              num_input  - The number of inputs per training data
-             num_output - The number of ouputs per training data
+             num_output - The number of outputs per training data
              input      - The set of inputs
              output     - The set of desired outputs
           
@@ -1485,6 +1485,26 @@ public:
 
       
 
+      void train_on_data(const std::vector< std::vector<fann_type> >& input,
+                         const std::vector< std::vector<fann_type> >& output,
+                         bool initWeights,
+                         unsigned int max_epochs,
+                         unsigned int epochs_between_reports,
+                         float desired_error)
+        {
+          if ((ann != NULL))
+            {
+              training_data data;
+              data.set_train_data(input,output);
+              if(data.train_data != NULL){
+                if(initWeights)
+                  init_weights(data);
+                fann_train_on_data(ann, data.train_data, max_epochs,
+                                   epochs_between_reports, desired_error);
+              }
+            }
+        }
+
         void train_on_data(const std::vector< Vector2d<fann_type> >& input,
                            unsigned int num_data,
                            bool initWeights,
@@ -1505,18 +1525,18 @@ public:
             }
         }
 
+      //cross validation for classification
         float cross_validation(std::vector< Vector2d<fann_type> >& trainingFeatures,
                                unsigned int ntraining,
                                unsigned short cv,
                                unsigned int max_epochs,
-                               unsigned int epochs_between_reports,
                                float desired_error,
-                               std::vector<unsigned short>& input,
-                               std::vector<unsigned short>& output,
+                               std::vector<unsigned short>& referenceVector,
+                               std::vector<unsigned short>& outputVector,
                                short verbose=0)
         {
-          input.clear();
-          output.clear();
+          referenceVector.clear();
+          outputVector.clear();
           assert(cv<ntraining);
           float rmse=0;
           int nclass=trainingFeatures.size();
@@ -1524,6 +1544,8 @@ public:
           int testclass=0;//class to leave out
           int testsample=0;//sample to leave out
           int nrun=(cv>1)? cv : ntraining;
+          if(nrun>ntraining)
+            nrun=ntraining;
           for(int irun=0;irun<nrun;++irun){
             if(verbose>1)
               std::cout << "run " << irun << std::endl;
@@ -1569,13 +1591,13 @@ public:
             if(verbose>1)
               cout << endl << "Set training data" << endl;
             bool initWeights=true;
+            unsigned int epochs_between_reports=0;
             train_on_data(trainingFeatures,ntraining-ntest,initWeights, max_epochs,
                           epochs_between_reports, desired_error);
             //cross validation with testFeatures
             if(verbose>1)
               cout << endl << "Cross validation" << endl;
 
-            //todo: run network and store result in vector
             vector<float> result(nclass);
             int maxClass=-1;
             for(int iclass=0;iclass<testFeatures.size();++iclass){
@@ -1592,8 +1614,8 @@ public:
                   }
                 }
                 assert(maxP>=0);
-                input.push_back(iclass);
-                output.push_back(maxClass);
+                referenceVector.push_back(iclass);
+                outputVector.push_back(maxClass);
               }
             }
 
@@ -1613,6 +1635,98 @@ public:
           return 0;
         }
 
+      //cross validation for regresssion
+        float cross_validation(std::vector< std::vector<fann_type> >& input,
+                               std::vector< std::vector<fann_type> >& output,
+                               unsigned short cv,
+                               unsigned int max_epochs,
+                               float desired_error,
+                               std::vector< std::vector<fann_type> >& referenceVector,
+                               std::vector< std::vector<fann_type> >& outputVector,
+                               short verbose=0)
+        {
+          assert(input.size());
+          assert(output.size()==input.size());
+          unsigned int ntraining=input.size();
+          unsigned int noutput=output[0].size();
+          referenceVector.clear();
+          outputVector.clear();
+          assert(cv<ntraining);
+          float rmse=0;
+          std::vector< std::vector<fann_type> > testInput;
+          std::vector< std::vector<fann_type> > testOutput;
+          int testsample=0;//sample to leave out
+          int nrun=(cv>1)? cv : ntraining;
+          if(nrun>ntraining)
+            nrun=ntraining;
+          for(int irun=0;irun<nrun;++irun){
+            if(verbose>1)
+              std::cout << "run " << irun << std::endl;
+            //reset training sample from last run
+            if(verbose>1)
+              std::cout << "reset training sample from last run" << std::endl;
+            while(testInput.size()){
+              input.push_back(testInput.back());
+              testInput.pop_back();
+            }
+            while(testOutput.size()){
+              output.push_back(testOutput.back());
+              testOutput.pop_back();
+            }
+            assert(testInput.size()==testOutput.size());
+            if(verbose>1){
+              std::cout << "training size: " << input.size() << std::endl;
+              std::cout << "test size: " << testInput.size() << std::endl;
+            }
+            assert(input.size());
+            //create test sample
+            if(verbose>1)
+              std::cout << "create test sample" << std::endl;
+            unsigned int nsample=0;
+            int ntest=(cv>1)? ntraining/cv : 1; //n-fold cross validation or leave-one-out
+            while(nsample<ntest){
+              testInput.push_back(input[0]);
+              testOutput.push_back(output[0]);
+              input.erase(input.begin());
+              output.erase(output.begin());
+              assert(input.size());
+              assert(output.size());
+              assert(input.size()==output.size());
+              ++nsample;
+            }
+            assert(nsample==ntest);
+            assert(testInput.size()==testOutput.size());
+            //training with left out training set
+            if(verbose>1)
+              cout << endl << "Set training data" << endl;
+            bool initWeights=true;
+            unsigned int epochs_between_reports=0;
+            
+            train_on_data(input,output,initWeights, max_epochs,
+                          epochs_between_reports, desired_error);
+            //cross validation with testFeatures
+            if(verbose>1)
+              cout << endl << "Cross validation" << endl;
+
+            vector<fann_type> result(noutput);
+            for(int isample=0;isample<testInput.size();++isample){
+              result=run(testInput[isample]);
+              referenceVector.push_back(testOutput[isample]);
+              outputVector.push_back(result);
+            }
+          }
+          //reset from very last run
+          while(testInput.size()){
+            input.push_back(testInput.back());
+            testInput.pop_back();
+          }
+          while(testOutput.size()){
+            output.push_back(testOutput.back());
+            testOutput.pop_back();
+          }
+          return 0;
+        }
+
         /* Method: train_on_file
            
            Does the same as <train_on_data>, but reads the training data directly from a file.
diff --git a/src/apps/Makefile.am b/src/apps/Makefile.am
index a55b751..6eeebb2 100644
--- a/src/apps/Makefile.am
+++ b/src/apps/Makefile.am
@@ -6,19 +6,22 @@ LDADD = $(GDAL_LDFLAGS) $(top_builddir)/src/algorithms/libalgorithms.la $(top_bu
 ###############################################################################
 
 # the program to build and install (the names of the final binaries)
-bin_PROGRAMS = pkinfo pkcrop pkreclass pkgetmask pksetmask pkcreatect pkdumpimg pkdumpogr pksieve pkstat pkstatogr pkegcs pkextract pkfillnodata pkfilter pkdsm2shadow pkmosaic pkndvi pkpolygonize pkascii2img pkdiff pkclassify_svm pkfs_svm pkascii2ogr pkeditogr
+bin_PROGRAMS = pkinfo pkcrop pkreclass pkgetmask pksetmask pkcreatect pkdumpimg pkdumpogr pksieve pkstat pkstatogr pkegcs pkextract pkfillnodata pkfilter pkfilterascii pkdsm2shadow pkmosaic pkndvi pkpolygonize pkascii2img pkdiff pkclassify_svm pkfs_svm pkascii2ogr pkeditogr
 
 # the program to build but not install (the names of the final binaries)
 #noinst_PROGRAMS =  pkxcorimg pkgeom
 
 if USE_FANN
-bin_PROGRAMS += pkclassify_nn pkfs_nn
+bin_PROGRAMS += pkclassify_nn pkfs_nn pkregression_nn
 pkclassify_nn_SOURCES = $(top_srcdir)/src/algorithms/myfann_cpp.h pkclassify_nn.h pkclassify_nn.cc
 pkclassify_nn_CXXFLAGS = -I$(top_srcdir)/src -I$(top_srcdir)/src/base $(FANN_CFLAGS) -I$(top_srcdir)/src/algorithms $(AM_CXXFLAGS)
 pkclassify_nn_LDADD = $(FANN_LIBS) $(FANN_CFLAGS) $(AM_LDFLAGS)
 pkfs_nn_SOURCES = $(top_srcdir)/src/algorithms/myfann_cpp.h pkfs_nn.cc
 pkfs_nn_CXXFLAGS = -I$(top_srcdir)/src -I$(top_srcdir)/src/base $(FANN_CFLAGS) -I$(top_srcdir)/src/algorithms $(AM_CXXFLAGS)
 pkfs_nn_LDADD = $(GSL_LIBS) $(FANN_LIBS) $(FANN_CFLAGS) $(AM_LDFLAGS)
+pkregression_nn_SOURCES = $(top_srcdir)/src/algorithms/myfann_cpp.h pkregression_nn.cc
+pkregression_nn_CXXFLAGS = -I$(top_srcdir)/src -I$(top_srcdir)/src/base $(FANN_CFLAGS) -I$(top_srcdir)/src/algorithms $(AM_CXXFLAGS)
+pkregression_nn_LDADD = $(GSL_LIBS) $(FANN_LIBS) $(FANN_CFLAGS) $(AM_LDFLAGS)
 endif
 
 if USE_LAS
@@ -43,14 +46,15 @@ pkdumpogr_SOURCES = pkdumpogr.h pkdumpogr.cc
 pksieve_SOURCES = pksieve.cc
 pkstat_SOURCES = $(top_srcdir)/src/algorithms/StatFactory.h pkstat.cc
 pkstat_LDADD = $(GSL_LIBS) $(AM_LDFLAGS)
-#pkstat_LDADD = -llas $(GSL_LIBS) $(AM_LDFLAGS)
 pkstatogr_SOURCES = pkstatogr.cc
 pkegcs_SOURCES = pkegcs.cc
 pkegcs_LDADD = -lgdal $(AM_LDFLAGS) -lgdal
 pkextract_SOURCES = pkextract.cc
 pkfillnodata_SOURCES = pkfillnodata.cc
 pkfilter_SOURCES = pkfilter.cc
-pkfilter_LDADD = $(GSL_LIBS) $(AM_LDFLAGS) -lgsl -lgdal #-lgslwrap
+pkfilter_LDADD = $(GSL_LIBS) $(AM_LDFLAGS) -lgsl -lgdal
+pkfilterascii_SOURCES = pkfilterascii.cc
+pkfilterascii_LDADD = $(GSL_LIBS) $(AM_LDFLAGS) -lgsl -lgdal
 pkdsm2shadow_SOURCES = pkdsm2shadow.cc
 pkmosaic_SOURCES = pkmosaic.cc
 pkndvi_SOURCES = pkndvi.cc
diff --git a/src/apps/pkclassify_nn.cc b/src/apps/pkclassify_nn.cc
index 04a72d2..4ebec2a 100644
--- a/src/apps/pkclassify_nn.cc
+++ b/src/apps/pkclassify_nn.cc
@@ -462,7 +462,6 @@ int main(int argc, char *argv[])
                                             ntraining,
                                             cv_opt[0],
                                             maxit_opt[0],
-                                            0,
                                             desired_error,
                                             referenceVector,
                                             outputVector,
diff --git a/src/apps/pkfilter.cc b/src/apps/pkfilter.cc
index 1a58215..689cb42 100644
--- a/src/apps/pkfilter.cc
+++ b/src/apps/pkfilter.cc
@@ -52,8 +52,8 @@ int main(int argc,char **argv) {
   Optionpk<double> tapz_opt("tapz", "tapz", "taps used for spectral filtering");
   Optionpk<double> fwhm_opt("fwhm", "fwhm", "list of full width half to apply spectral filtering (-fwhm band1 -fwhm band2 ...)");
   Optionpk<std::string> srf_opt("srf", "srf", "list of ASCII files containing spectral response functions (two columns: wavelength response)");
-  Optionpk<double> wavelengthIn_opt("win", "wavelengthIn", "list of wavelengths in input spectrum (-w band1 -w band2 ...)");
-  Optionpk<double> wavelengthOut_opt("wout", "wavelengthOut", "list of wavelengths in output spectrum (-w band1 -w band2 ...)");
+  Optionpk<double> wavelengthIn_opt("win", "wavelengthIn", "list of wavelengths in input spectrum (-win band1 -win band2 ...)");
+  Optionpk<double> wavelengthOut_opt("wout", "wavelengthOut", "list of wavelengths in output spectrum (-wout band1 -wout band2 ...)");
   Optionpk<std::string> interpolationType_opt("interp", "interp", "type of interpolation for spectral filtering (see http://www.gnu.org/software/gsl/manual/html_node/Interpolation-Types.html)","akima");
   Optionpk<std::string>  otype_opt("ot", "otype", "Data type for output image ({Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}). Empty string: inherit type from input image","");
   Optionpk<string>  oformat_opt("of", "oformat", "Output image format (see also gdal_translate). Empty string: inherit from input image");
diff --git a/src/apps/pkfilterascii.cc b/src/apps/pkfilterascii.cc
new file mode 100644
index 0000000..434a5a0
--- /dev/null
+++ b/src/apps/pkfilterascii.cc
@@ -0,0 +1,233 @@
+/**********************************************************************
+pkfilterascii.cc: program to filter data in ASCII file (spectral respons function, dwt)
+Copyright (C) 2008-2013 Pieter Kempeneers
+
+This file is part of pktools
+
+pktools is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+pktools is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with pktools.  If not, see <http://www.gnu.org/licenses/>.
+***********************************************************************/
+#include <assert.h>
+#include <iostream>
+#include <string>
+#include <fstream>
+#include <math.h>
+#include <sys/types.h>
+#include <stdio.h>
+#include "base/Optionpk.h"
+#include "base/Vector2d.h"
+#include "algorithms/Filter.h"
+#include "fileclasses/FileReaderAscii.h"
+
+/*------------------
+  Main procedure
+  ----------------*/
+int main(int argc,char **argv) {
+  Optionpk<std::string> input_opt("i","input","input ASCII file");
+  Optionpk<std::string> output_opt("o", "output", "Output ASCII file");
+  Optionpk<int> inputCols_opt("ic", "inputCols", "input columns (e.g., for three dimensional input data in first three columns use: -ic 0 -ic 1 -ic 2"); 
+  Optionpk<std::string> method_opt("f", "filter", "filter function (to be implemented: dwtForward, dwtInverse)");
+    Optionpk<std::string> wavelet_type_opt("wt", "wavelet", "wavelet type: daubechies,daubechies_centered, haar, haar_centered, bspline, bspline_centered", "daubechies");
+    Optionpk<int> family_opt("wf", "family", "wavelet family (vanishing moment, see also http://www.gnu.org/software/gsl/manual/html_node/DWT-Initialization.html)", 4);
+  Optionpk<int> dimZ_opt("dz", "dz", "filter kernel size in z (band or spectral dimension), must be odd (example: 3).. Set dz>0 if 1-D filter must be used in band domain");
+  Optionpk<double> tapz_opt("tapz", "tapz", "taps used for spectral filtering");
+  Optionpk<double> fwhm_opt("fwhm", "fwhm", "list of full width half to apply spectral filtering (-fwhm band1 -fwhm band2 ...)");
+  Optionpk<std::string> srf_opt("srf", "srf", "list of ASCII files containing spectral response functions (two columns: wavelength response)");
+  Optionpk<int> wavelengthIn_opt("win", "wavelengthIn", "column number of input ASCII file containing wavelengths");
+  Optionpk<double> wavelengthOut_opt("wout", "wavelengthOut", "list of wavelengths in output spectrum (-wout band1 -wout band2 ...)");
+  Optionpk<std::string> interpolationType_opt("interp", "interp", "type of interpolation for spectral filtering (see http://www.gnu.org/software/gsl/manual/html_node/Interpolation-Types.html)","akima");
+  Optionpk<bool> transpose_opt("t", "transpose", "transpose output with samples in rows and wavelengths in cols", false);
+  Optionpk<short> verbose_opt("v", "verbose", "verbose mode if > 0", 0);
+
+  bool doProcess;//stop process when program was invoked with help option (-h --help)
+  try{
+    doProcess=input_opt.retrieveOption(argc,argv);
+    output_opt.retrieveOption(argc,argv);
+    inputCols_opt.retrieveOption(argc,argv);
+    method_opt.retrieveOption(argc,argv);
+    dimZ_opt.retrieveOption(argc,argv);
+    tapz_opt.retrieveOption(argc,argv);
+    fwhm_opt.retrieveOption(argc,argv);
+    srf_opt.retrieveOption(argc,argv);
+    wavelengthIn_opt.retrieveOption(argc,argv);
+    wavelengthOut_opt.retrieveOption(argc,argv);
+    interpolationType_opt.retrieveOption(argc,argv);
+    transpose_opt.retrieveOption(argc,argv);
+    verbose_opt.retrieveOption(argc,argv);
+  }
+  catch(string predefinedString){
+    std::cout << predefinedString << std::endl;
+    exit(0);
+  }
+  if(!doProcess){
+    std::cout << "short option -h shows basic options only, use long option --help to show all options" << std::endl;
+    exit(0);//help was invoked, stop processing
+  }
+
+  Vector2d<double> inputData(inputCols_opt.size());
+  Vector2d<double> filteredData(inputCols_opt.size());
+  vector<double> wavelengthIn;
+  vector<double> wavelengthOut;
+  assert(input_opt.size());
+  FileReaderAscii asciiReader(input_opt[0]);
+  if(wavelengthIn_opt.size())
+    asciiReader.readData(wavelengthIn,wavelengthIn_opt[0]);
+  assert(inputCols_opt.size());
+  asciiReader.readData(inputData,inputCols_opt);
+  if(verbose_opt[0]){
+    std::cout << "wavelengthIn.size(): " << wavelengthIn.size() << std::endl;
+    std::cout << "inputData[0].size(): " << inputData[0].size() << std::endl;
+  }
+  assert(wavelengthIn.size()==inputData[0].size());
+  asciiReader.close();
+  filter::Filter filter1d;
+  if(fwhm_opt.size()){
+    filteredData.resize(inputCols_opt.size(),wavelengthOut_opt.size());
+    assert(wavelengthIn_opt.size());
+    if(verbose_opt[0])
+      std::cout << "spectral filtering to " << fwhm_opt.size() << " bands with provided fwhm " << std::endl;
+    assert(wavelengthOut_opt.size()==fwhm_opt.size());
+    vector<double> fwhmData(wavelengthOut_opt.size());
+    for(int icol=0;icol<inputCols_opt.size();++icol)
+      filter1d.applyFwhm<double>(wavelengthIn,inputData[icol], wavelengthOut_opt,fwhm_opt, interpolationType_opt[0], filteredData[icol],verbose_opt[0]);
+    if(verbose_opt[0])
+      std::cout << "spectra filtered to " << wavelengthOut_opt.size() << " bands" << std::endl;
+    wavelengthOut=wavelengthOut_opt;
+  }
+  else if(srf_opt.size()){
+    wavelengthOut.resize(srf_opt.size());
+    filteredData.resize(inputCols_opt.size(),srf_opt.size());
+    Vector2d<double> srfData(srf_opt.size(),inputCols_opt.size());//transposed output
+    if(verbose_opt[0])
+      std::cout << "spectral filtering to " << srf_opt.size() << " bands with provided SRF " << std::endl;
+    vector< Vector2d<double> > srf(srf_opt.size());//[0] srf_nr, [1]: wavelength, [2]: response
+    ifstream srfFile;
+    for(int isrf=0;isrf<srf_opt.size();++isrf){
+      srf[isrf].resize(2);
+      srfFile.open(srf_opt[isrf].c_str());
+      double v;
+      //add 0 to make sure srf is 0 at boundaries after interpolation step
+      srf[isrf][0].push_back(0);
+      srf[isrf][1].push_back(0);
+      srf[isrf][0].push_back(1);
+      srf[isrf][1].push_back(0);
+      while(srfFile >> v){
+        srf[isrf][0].push_back(v);
+        srfFile >> v;
+        srf[isrf][1].push_back(v);
+      }
+      srfFile.close();
+      //add 0 to make sure srf[isrf] is 0 at boundaries after interpolation step
+      srf[isrf][0].push_back(srf[isrf][0].back()+1);
+      srf[isrf][1].push_back(0);    
+      srf[isrf][0].push_back(srf[isrf][0].back()+1);
+      srf[isrf][1].push_back(0);
+      if(verbose_opt[0])
+        cout << "srf file details: " << srf[isrf][0].size() << " wavelengths defined" << endl;
+      if(verbose_opt[0]>1){
+        for(int iw=0;iw<srf[isrf][0].size();++iw)
+          std::cout << srf[isrf][0][iw] << " " << srf[isrf][1][iw] << std::endl;
+      }
+    }
+    double centreWavelength=0;
+    for(int icol=0;icol<inputCols_opt.size();++icol)
+      filteredData[icol].resize(srf.size());
+
+    for(int isrf=0;isrf<srf.size();++isrf){
+      double delta=1.0;
+      bool normalize=true;
+      centreWavelength=filter1d.applySrf<double>(wavelengthIn,inputData,srf[isrf], interpolationType_opt[0], srfData[isrf], delta, normalize,1,true, verbose_opt[0]);
+      if(verbose_opt[0])
+        std::cout << "centre wavelength srf " << isrf << ": " << centreWavelength << std::endl;
+      wavelengthOut[isrf]=static_cast<int>(centreWavelength+0.5);
+    }
+    srfData.transpose(filteredData);
+    if(verbose_opt[0])
+      std::cout << "spectra filtered to " << srf.size() << " bands" << std::endl;
+  }
+  else{//no filtering
+    if(verbose_opt[0])
+      std::cout << "no filtering selected" << std::endl;
+    wavelengthOut=wavelengthIn;
+    for(int icol=0;icol<inputCols_opt.size();++icol)
+      filteredData[icol]=inputData[icol];
+  }
+  
+  if(method_opt.size()){
+    for(int icol=0;icol<inputCols_opt.size();++icol){
+      switch(filter::Filter::getFilterType(method_opt[0])){
+      // case(filter::dwtForward):
+      //   filter1d.dwtForward(filteredData[icol],wavelet_type_opt[0],family_opt[0]);
+      //   break;
+      // case(filter::dwtInverse):
+      //   filter1d.dwtInverse(filteredData[icol],wavelet_type_opt[0],family_opt[0]);
+      //   break;
+      default:
+        if(verbose_opt[0])
+          std::cout << "method to be implemented" << std::endl;
+        break;
+      }
+    }
+  }
+  ofstream outputStream;
+  if(!output_opt.empty())
+    outputStream.open(output_opt[0].c_str(),ios::out);
+
+  if(transpose_opt[0]){
+    for(int icol=0;icol<inputCols_opt.size();++icol){
+      for(int iband=0;iband<wavelengthOut.size();++iband){
+        if(!output_opt.empty()){
+          outputStream << filteredData[icol][iband];
+          if(iband<wavelengthOut.size()-1)
+            outputStream << " ";
+          else
+            outputStream << std::endl;
+        }
+        else{
+          std::cout << filteredData[icol][iband];
+          if(iband<wavelengthOut.size()-1)
+            std::cout << " ";
+          else
+            std::cout << std::endl;
+        }
+      }
+    }    
+  }
+  else{
+    for(int iband=0;iband<wavelengthOut.size();++iband){
+      if(!output_opt.empty())
+        outputStream << wavelengthOut[iband] << " ";
+      else
+        std::cout << wavelengthOut[iband] << " ";
+      for(int icol=0;icol<inputCols_opt.size();++icol){
+        if(!output_opt.empty()){
+          outputStream << filteredData[icol][iband];
+          if(icol<inputCols_opt.size()-1)
+            outputStream << " ";
+          else
+            outputStream << std::endl;
+        }
+        else{
+          std::cout << filteredData[icol][iband];
+          if(icol<inputCols_opt.size()-1)
+            std::cout << " ";
+          else
+            std::cout << std::endl;
+        }
+      }
+    }    
+  }
+  if(!output_opt.empty())
+    outputStream.close();
+  return 0;
+}
diff --git a/src/apps/pkfs_nn.cc b/src/apps/pkfs_nn.cc
index e87a24c..2fbca98 100644
--- a/src/apps/pkfs_nn.cc
+++ b/src/apps/pkfs_nn.cc
@@ -129,7 +129,6 @@ double getCost(const vector<Vector2d<float> > &trainingFeatures)
                               ntraining,
                               cv_opt[0],
                               maxit_opt[0],
-                              0,
                               desired_error,
                               referenceVector,
                               outputVector,
diff --git a/src/apps/pkregression_nn.cc b/src/apps/pkregression_nn.cc
new file mode 100644
index 0000000..35e0a53
--- /dev/null
+++ b/src/apps/pkregression_nn.cc
@@ -0,0 +1,346 @@
+/**********************************************************************
+pkregression_nn.cc: regression with artificial neural network (multi-layer perceptron)
+Copyright (C) 2008-2013 Pieter Kempeneers
+
+This file is part of pktools
+
+pktools is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+pktools is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with pktools.  If not, see <http://www.gnu.org/licenses/>.
+***********************************************************************/
+#include <vector>
+#include <fstream>
+#include "base/Optionpk.h"
+#include "fileclasses/FileReaderAscii.h"
+#include "floatfann.h"
+#include "myfann_cpp.h"
+
+int main(int argc, char *argv[])
+{
+  //--------------------------- command line options ------------------------------------
+  Optionpk<string> input_opt("i", "input", "input ASCII file"); 
+  Optionpk<string> output_opt("o", "output", "output ASCII file for result"); 
+  Optionpk<int> inputCols_opt("ic", "inputCols", "input columns (e.g., for three dimensional input data in first three columns use: -ic 0 -ic 1 -ic 2"); 
+  Optionpk<int> outputCols_opt("oc", "outputCols", "output columns (e.g., for two dimensional output in columns 3 and 4 (starting from 0) use: -oc 3 -oc 4"); 
+  Optionpk<string> training_opt("t", "training", "training ASCII file (each row represents one sampling unit. Input features should be provided as columns, followed by output)"); 
+  Optionpk<double> from_opt("from", "from", "start from this row in training file (start from 0)",0); 
+  Optionpk<double> to_opt("to", "to", "read until this row in training file (start from 0 or set leave 0 as default to read until end of file)", 0); 
+  Optionpk<short> band_opt("b", "band", "band index (starting from 0, either use band option or use start to end)");
+  Optionpk<double> offset_opt("\0", "offset", "offset value for each spectral band input features: refl[band]=(DN[band]-offset[band])/scale[band]", 0.0);
+  Optionpk<double> scale_opt("\0", "scale", "scale value for each spectral band input features: refl=(DN[band]-offset[band])/scale[band] (use 0 if scale min and max in each band to -1.0 and 1.0)", 0.0);
+  Optionpk<unsigned short> cv_opt("cv", "cv", "n-fold cross validation mode",0);
+  Optionpk<unsigned int> nneuron_opt("\0", "nneuron", "number of neurons in hidden layers in neural network (multiple hidden layers are set by defining multiple number of neurons: -n 15 -n 1, default is one hidden layer with 5 neurons)", 5); 
+  Optionpk<float> connection_opt("\0", "connection", "connection reate (default: 1.0 for a fully connected network)", 1.0); 
+  Optionpk<float> weights_opt("w", "weights", "weights for neural network. Apply to fully connected network only, starting from first input neuron to last output neuron, including the bias neurons (last neuron in each but last layer)", 0.0); 
+  Optionpk<float> learning_opt("l", "learning", "learning rate (default: 0.7)", 0.7); 
+  Optionpk<unsigned int> maxit_opt("\0", "maxit", "number of maximum iterations (epoch) (default: 500)", 500); 
+  Optionpk<short> verbose_opt("v", "verbose", "set to: 0 (results only), 1 (confusion matrix), 2 (debug)",0);
+
+  bool doProcess;//stop process when program was invoked with help option (-h --help)
+  try{
+    doProcess=input_opt.retrieveOption(argc,argv);
+    output_opt.retrieveOption(argc,argv);
+    inputCols_opt.retrieveOption(argc,argv);
+    outputCols_opt.retrieveOption(argc,argv);
+    training_opt.retrieveOption(argc,argv);
+    from_opt.retrieveOption(argc,argv);
+    to_opt.retrieveOption(argc,argv);
+    band_opt.retrieveOption(argc,argv);
+    offset_opt.retrieveOption(argc,argv);
+    scale_opt.retrieveOption(argc,argv);
+    cv_opt.retrieveOption(argc,argv);
+    nneuron_opt.retrieveOption(argc,argv);
+    connection_opt.retrieveOption(argc,argv);
+    weights_opt.retrieveOption(argc,argv);
+    learning_opt.retrieveOption(argc,argv);
+    maxit_opt.retrieveOption(argc,argv);
+    verbose_opt.retrieveOption(argc,argv);
+  }
+  catch(string predefinedString){
+    std::cout << predefinedString << std::endl;
+    exit(0);
+  }
+  if(!doProcess){
+    std::cout << "short option -h shows basic options only, use long option --help to show all options" << std::endl;
+    exit(0);//help was invoked, stop processing
+  }
+
+  unsigned int ninput=inputCols_opt.size();
+  unsigned int noutput=outputCols_opt.size();
+  assert(ninput);
+  assert(noutput);
+  vector< vector<float> > inputUnits;
+  vector< vector<float> > trainingUnits;
+  vector< vector<float> > trainingOutput;
+  FileReaderAscii inputFile;
+  unsigned int inputSize=0;
+  if(input_opt.size()){
+    inputFile.open(input_opt[0]);
+    inputFile.setMinRow(from_opt[0]);
+    inputFile.setMaxRow(to_opt[0]);
+    inputFile.setComment('#');
+    inputFile.readData(inputUnits,inputCols_opt,1,0,true,verbose_opt[0]);
+    inputFile.close();
+    inputSize=inputUnits.size();
+  }
+  FileReaderAscii trainingFile(training_opt[0]);
+  unsigned int sampleSize=0;
+  trainingFile.setMinRow(from_opt[0]);
+  trainingFile.setMaxRow(to_opt[0]);
+  trainingFile.setComment('#');
+  trainingFile.readData(trainingUnits,inputCols_opt,1,0,true,verbose_opt[0]);
+  trainingFile.readData(trainingOutput,outputCols_opt,1,0,true,verbose_opt[0]);
+  trainingFile.close();
+  sampleSize=trainingUnits.size();
+
+  if(verbose_opt[0]>1){
+    std::cout << "sampleSize: " << sampleSize << std::endl;
+    std::cout << "ninput: " << ninput << std::endl;
+    std::cout << "noutput: " << noutput << std::endl;
+    std::cout << "trainingUnits[0].size(): " << trainingUnits[0].size() << std::endl;
+    std::cout << "trainingOutput[0].size(): " << trainingOutput[0].size() << std::endl;
+    std::cout << "trainingUnits.size(): " << trainingUnits.size() << std::endl;
+    std::cout << "trainingOutput.size(): " << trainingOutput.size() << std::endl;
+  }
+
+  assert(ninput==trainingUnits[0].size());
+  assert(noutput==trainingOutput[0].size());
+  assert(trainingUnits.size()==trainingOutput.size());
+
+  //set scale and offset
+  if(offset_opt.size()>1)
+    assert(offset_opt.size()==ninput);
+  if(scale_opt.size()>1)
+    assert(scale_opt.size()==ninput);
+
+  std::vector<float> offset_input(ninput);
+  std::vector<float> scale_input(ninput);
+
+  std::vector<float> offset_output(noutput);
+  std::vector<float> scale_output(noutput);
+
+  for(int iinput=0;iinput<ninput;++iinput){
+    if(verbose_opt[0]>=1)
+      cout << "scaling for input feature" << iinput << endl;
+    offset_input[iinput]=(offset_opt.size()==1)?offset_opt[0]:offset_opt[iinput];
+    scale_input[iinput]=(scale_opt.size()==1)?scale_opt[0]:scale_opt[iinput];
+    //search for min and maximum
+    if(scale_input[iinput]<=0){
+      float theMin=trainingUnits[0][iinput];
+      float theMax=trainingUnits[0][iinput];
+      for(int isample=0;isample<trainingUnits.size();++isample){
+        if(theMin>trainingUnits[isample][iinput])
+          theMin=trainingUnits[isample][iinput];
+        if(theMax<trainingUnits[isample][iinput])
+          theMax=trainingUnits[isample][iinput];
+      }
+      offset_input[iinput]=theMin+(theMax-theMin)/2.0;
+      scale_input[iinput]=(theMax-theMin)/2.0;
+      if(verbose_opt[0]>=1){
+        std::cout << "Extreme image values for input feature " << iinput << ": [" << theMin << "," << theMax << "]" << std::endl;
+        std::cout << "Using offset, scale: " << offset_input[iinput] << ", " << scale_input[iinput] << std::endl;
+        std::cout << "scaled values for input feature " << iinput << ": [" << (theMin-offset_input[iinput])/scale_input[iinput] << "," << (theMax-offset_input[iinput])/scale_input[iinput] << "]" << std::endl;
+      }
+    }
+  }
+
+  for(int ioutput=0;ioutput<noutput;++ioutput){
+    if(verbose_opt[0]>=1)
+      cout << "scaling for output feature" << ioutput << endl;
+    //search for min and maximum
+    float theMin=trainingOutput[0][ioutput];
+    float theMax=trainingOutput[0][ioutput];
+    for(int isample=0;isample<trainingOutput.size();++isample){
+      if(theMin>trainingOutput[isample][ioutput])
+        theMin=trainingOutput[isample][ioutput];
+      if(theMax<trainingOutput[isample][ioutput])
+        theMax=trainingOutput[isample][ioutput];
+    }
+    offset_output[ioutput]=theMin+(theMax-theMin)/2.0;
+    scale_output[ioutput]=(theMax-theMin)/2.0;
+    if(verbose_opt[0]>=1){
+      std::cout << "Extreme image values for output feature " << ioutput << ": [" << theMin << "," << theMax << "]" << std::endl;
+      std::cout << "Using offset, scale: " << offset_output[ioutput] << ", " << scale_output[ioutput] << std::endl;
+      std::cout << "scaled values for output feature " << ioutput << ": [" << (theMin-offset_output[ioutput])/scale_output[ioutput] << "," << (theMax-offset_output[ioutput])/scale_output[ioutput] << "]" << std::endl;
+    }
+  }
+
+
+
+  FANN::neural_net net;//the neural network
+
+  const unsigned int num_layers = nneuron_opt.size()+2;
+  const float desired_error = 0.0003;
+  const unsigned int iterations_between_reports = (verbose_opt[0])? maxit_opt[0]+1:0;
+  if(verbose_opt[0]>=1){
+    cout << "creating artificial neural network with " << nneuron_opt.size() << " hidden layer, having " << endl;
+    for(int ilayer=0;ilayer<nneuron_opt.size();++ilayer)
+      cout << nneuron_opt[ilayer] << " ";
+    cout << "neurons" << endl;
+  }
+
+  switch(num_layers){
+  case(3):
+    net.create_sparse(connection_opt[0],num_layers, ninput, nneuron_opt[0], noutput);
+    break;
+  case(4):
+    net.create_sparse(connection_opt[0],num_layers, ninput, nneuron_opt[0], nneuron_opt[1], noutput);
+    break;
+  default:
+    cerr << "Only 1 or 2 hidden layers are supported!" << endl;
+    exit(1);
+    break;
+  }
+  if(verbose_opt[0]>=1)
+    cout << "network created" << endl;
+  
+  net.set_learning_rate(learning_opt[0]);
+
+  //   net.set_activation_steepness_hidden(1.0);
+  //   net.set_activation_steepness_output(1.0);
+    
+  net.set_activation_function_hidden(FANN::SIGMOID_SYMMETRIC_STEPWISE);
+  net.set_activation_function_output(FANN::SIGMOID_SYMMETRIC_STEPWISE);
+
+  // Set additional properties such as the training algorithm
+  //   net.set_training_algorithm(FANN::TRAIN_QUICKPROP);
+
+  // Output network type and parameters
+  if(verbose_opt[0]>=1){
+    cout << endl << "Network Type                         :  ";
+    switch (net.get_network_type())
+      {
+      case FANN::LAYER:
+        cout << "LAYER" << endl;
+        break;
+      case FANN::SHORTCUT:
+        cout << "SHORTCUT" << endl;
+        break;
+      default:
+        cout << "UNKNOWN" << endl;
+        break;
+      }
+    net.print_parameters();
+  }
+      
+  if(verbose_opt[0]>=1){
+    cout << "Max Epochs " << setw(8) << maxit_opt[0] << ". "
+         << "Desired Error: " << left << desired_error << right << endl;
+  }
+  bool initWeights=true;
+
+  Vector2d<float> trainingFeatures(sampleSize,ninput);
+  for(unsigned int isample=0;isample<sampleSize;++isample){
+    for(unsigned int iinput=0;iinput<ninput;++iinput)
+      trainingFeatures[isample][iinput]=(trainingUnits[isample][iinput]-offset_input[iinput])/scale_input[iinput];
+  }
+
+  Vector2d<float> scaledOutput(sampleSize,noutput);
+  for(unsigned int isample=0;isample<sampleSize;++isample){
+    for(unsigned int ioutput=0;ioutput<noutput;++ioutput)
+      scaledOutput[isample][ioutput]=(trainingOutput[isample][ioutput]-offset_output[ioutput])/scale_output[ioutput];
+  }
+
+  if(cv_opt[0]){
+    if(verbose_opt[0])
+      std::cout << "cross validation" << std::endl;
+    std::vector< std::vector<float> > referenceVector;
+    std::vector< std::vector<float> > outputVector;
+    net.cross_validation(trainingFeatures,
+                         scaledOutput,
+                         cv_opt[0],
+                         maxit_opt[0],
+                         desired_error,
+                         referenceVector,
+                         outputVector);
+    assert(referenceVector.size()==outputVector.size());
+    vector<double> rmse(noutput);
+    for(int isample=0;isample<referenceVector.size();++isample){
+      std::cout << isample << " ";
+      for(int ioutput=0;ioutput<noutput;++ioutput){
+        if(!isample)
+          rmse[ioutput]=0;
+        double ref=scale_output[ioutput]*referenceVector[isample][ioutput]+offset_output[ioutput];
+        double val=scale_output[ioutput]*outputVector[isample][ioutput]+offset_output[ioutput];
+        rmse[ioutput]+=(ref-val)*(ref-val);
+        std::cout << ref << " " << val;
+        if(ioutput<noutput-1)
+          std::cout << " ";
+        else
+          std::cout << std::endl;
+      }
+    }
+    for(int ioutput=0;ioutput<noutput;++ioutput)
+      std::cout << "rmse output variable " << ioutput << ": " << sqrt(rmse[ioutput]/referenceVector.size()) << std::endl;
+  }
+
+
+  net.train_on_data(trainingFeatures,
+                    scaledOutput,
+                    initWeights,
+                    maxit_opt[0],
+                    iterations_between_reports,
+                    desired_error);
+
+
+  if(verbose_opt[0]>=2){
+    net.print_connections();
+    vector<fann_connection> convector;
+    net.get_connection_array(convector);
+    for(int i_connection=0;i_connection<net.get_total_connections();++i_connection)
+      cout << "connection " << i_connection << ": " << convector[i_connection].weight << endl;
+  }
+  //end of training
+
+  ofstream outputStream;
+  if(!output_opt.empty())
+    outputStream.open(output_opt[0].c_str(),ios::out);
+  for(unsigned int isample=0;isample<inputUnits.size();++isample){
+    std::vector<float> inputFeatures(ninput);
+    for(unsigned int iinput=0;iinput<ninput;++iinput)
+      inputFeatures[iinput]=(inputUnits[isample][iinput]-offset_input[iinput])/scale_input[iinput];
+    vector<float> result(noutput);
+    result=net.run(inputFeatures);
+
+    if(!output_opt.empty())
+      outputStream << isample << " ";
+    else
+      std::cout << isample << " ";
+    if(verbose_opt[0]){
+      for(unsigned int iinput=0;iinput<ninput;++iinput){
+        if(output_opt.size())
+          outputStream << inputUnits[isample][iinput] << " ";
+        else
+          std::cout << inputUnits[isample][iinput] << " ";
+      }
+    }
+    for(unsigned int ioutput=0;ioutput<noutput;++ioutput){
+      result[ioutput]=scale_output[ioutput]*result[ioutput]+offset_output[ioutput];
+      if(output_opt.size()){
+        outputStream << result[ioutput];
+        if(ioutput<noutput-1)
+          outputStream << " ";
+        else
+          outputStream << std::endl;
+      }
+      else{
+        std::cout << result[ioutput];
+        if(ioutput<noutput-1)
+          std::cout << " ";
+        else
+          std::cout << std::endl;
+      }
+    }
+  }
+  if(!output_opt.empty())
+    outputStream.close();
+}
diff --git a/src/base/Vector2d.h b/src/base/Vector2d.h
index 81d5909..a8f0ba7 100644
--- a/src/base/Vector2d.h
+++ b/src/base/Vector2d.h
@@ -50,6 +50,14 @@ public:
   void selectCol(int col, T* output) const;
   vector<T> selectCol(int col);
   void selectCols(const list<int> &cols, Vector2d<T> &output) const;
+  void transpose(Vector2d<T> &output) const{
+    output.resize(nCols(),nRows());
+    for(int irow=0;irow<nRows();++irow){
+      for(int icol=0;icol<nCols();++icol){
+        output[icol][irow]=(*this)[irow][icol];
+      }
+    }
+  };
   void selectCols(const list<int> &cols);
   void sort(Vector2d<T>& output);  
   void scale(const vector<double> &scaleVector, const vector<double> &offsetVector, Vector2d<T>& scaledOutput);
diff --git a/src/fileclasses/FileReaderAscii.h b/src/fileclasses/FileReaderAscii.h
index 11edf02..e1e9e83 100644
--- a/src/fileclasses/FileReaderAscii.h
+++ b/src/fileclasses/FileReaderAscii.h
@@ -24,6 +24,7 @@ along with pktools.  If not, see <http://www.gnu.org/licenses/>.
 #include <vector>
 #include <fstream>
 #include "base/Optionpk.h"
+#include <armadillo>
 
 //--------------------------------------------------------------------------
 class FileReaderAscii
@@ -42,8 +43,9 @@ public:
   void setComment(char comment){m_comment=comment;};
   unsigned int nrOfCol(bool checkCols=false, bool verbose=false);
   unsigned int nrOfRow(bool checkCols=false, bool verbose=false);
-  template<class T> unsigned int readData(std::vector<std::vector<T> > &dataVector, const std::vector<int> &cols, double scale=1.0, double offset=0.0, bool verbose=false);
+  template<class T> unsigned int readData(std::vector<std::vector<T> > &dataVector, const std::vector<int> &cols, double scale=1.0, double offset=0.0, bool transpose=false, bool verbose=false);
   template<class T> unsigned int readData(std::vector<T> &dataVector, int col, double scale=1.0, double offset=0, bool verbose=false);
+
   protected:
   std::string m_filename;
   std::ifstream m_ifstream;
@@ -168,10 +170,11 @@ template<class T> unsigned int FileReaderAscii::readData(std::vector<T> &dataVec
   return dataVector.size();
 }
 
-template<class T> unsigned int FileReaderAscii::readData(std::vector<std::vector<T> > &dataVector, const std::vector<int> &cols, double scale, double offset, bool verbose){
+template<class T> unsigned int FileReaderAscii::readData(std::vector<std::vector<T> > &dataVector, const std::vector<int> &cols, double scale, double offset, bool transpose, bool verbose){
   reset();
   dataVector.clear();
-  dataVector.resize(cols.size());
+  if(!transpose)
+    dataVector.resize(cols.size());
   int nrow=0;
   bool withinRange=true;
   if(m_fs>' '&&m_fs<='~'){//field separator is a regular character (minimum ASCII code is space, maximum ASCII code is tilde)
@@ -179,6 +182,7 @@ template<class T> unsigned int FileReaderAscii::readData(std::vector<std::vector
       std::cout << "reading csv file " << m_filename << std::endl;
     std::string csvRecord;
     while(getline(m_ifstream,csvRecord)){//read a line
+      std::vector<T> sampleVector;
       withinRange=true;
       if(nrow<m_minRow)
         withinRange=false;
@@ -207,8 +211,12 @@ template<class T> unsigned int FileReaderAscii::readData(std::vector<std::vector
             if(ncol==cols[icol]){
               T value=scale*string2type<T>(item)+offset;
               // T value=string2type<T>(item);
-              if((value>=m_min&&value<=m_max)||m_max<=m_min)
-                dataVector[icol].push_back(value);
+              if((value>=m_min&&value<=m_max)||m_max<=m_min){
+                if(transpose)
+                  sampleVector.push_back(value);
+                else
+                  dataVector[icol].push_back(value);
+              }
             }
           }
           ++ncol;
@@ -217,9 +225,11 @@ template<class T> unsigned int FileReaderAscii::readData(std::vector<std::vector
         }
         if(verbose)
           std::cout << std::endl;
-        if(dataVector.back().size())
-          assert(ncol>=cols[0]);
+        // if(dataVector.back().size())
+        //   assert(ncol>=cols[0]);
       }
+      if(transpose)
+        dataVector.push_back(sampleVector);
       ++nrow;
     }
     assert(dataVector.size());
@@ -229,6 +239,7 @@ template<class T> unsigned int FileReaderAscii::readData(std::vector<std::vector
       std::cout << "space or tab delimited fields" << std::endl;
     std::string spaceRecord;
     while(!getline(m_ifstream, spaceRecord).eof()){
+      std::vector<T> sampleVector;
       withinRange=true;
       if(nrow<m_minRow)
         withinRange=false;
@@ -260,8 +271,12 @@ template<class T> unsigned int FileReaderAscii::readData(std::vector<std::vector
           // T value=string2type<T>(item);
           for(int icol=0;icol<cols.size();++icol){
             if(ncol==cols[icol]){
-              if((value>=m_min&&value<=m_max)||m_max<=m_min)
-                dataVector[icol].push_back(value);
+              if((value>=m_min&&value<=m_max)||m_max<=m_min){
+                if(transpose)
+                  sampleVector.push_back(value);
+                else
+                  dataVector[icol].push_back(value);
+              }
             }
           }
           ++ncol;
@@ -272,12 +287,15 @@ template<class T> unsigned int FileReaderAscii::readData(std::vector<std::vector
           std::cout << std::endl;
         if(verbose)
           std::cout << "number of columns: " << ncol << std::endl;
-        if(dataVector.back().size())
-          assert(ncol>=cols[0]);
+        // if(dataVector.back().size())
+        //   assert(ncol>=cols[0]);
       }
+      if(transpose)
+        dataVector.push_back(sampleVector);
       ++nrow;
     }
   }
   return dataVector.size();
 }
+
 #endif // _IMGREADERASCII_H_

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-grass/pktools.git



More information about the Pkg-grass-devel mailing list