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 <kempe...@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
_______________________________________________
Pkg-grass-devel mailing list
Pkg-grass-devel@lists.alioth.debian.org
http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-grass-devel
