pktools 2.6.7
Processing Kernel for geospatial data
pkregann.cc
1/**********************************************************************
2pkregann.cc: regression with artificial neural network (multi-layer perceptron)
3Copyright (C) 2008-2014 Pieter Kempeneers
4
5This file is part of pktools
6
7pktools is free software: you can redistribute it and/or modify
8it under the terms of the GNU General Public License as published by
9the Free Software Foundation, either version 3 of the License, or
10(at your option) any later version.
11
12pktools is distributed in the hope that it will be useful,
13but WITHOUT ANY WARRANTY; without even the implied warranty of
14MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15GNU General Public License for more details.
16
17You should have received a copy of the GNU General Public License
18along with pktools. If not, see <http://www.gnu.org/licenses/>.
19***********************************************************************/
20#include <vector>
21#include <fstream>
22#include "base/Optionpk.h"
23#include "fileclasses/FileReaderAscii.h"
24#include "floatfann.h"
25#include "algorithms/myfann_cpp.h"
26/******************************************************************************/
70using namespace std;
71
72int main(int argc, char *argv[])
73{
74 //--------------------------- command line options ------------------------------------
75 Optionpk<string> input_opt("i", "input", "input ASCII file");
76 Optionpk<string> output_opt("o", "output", "output ASCII file for result");
77 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");
78 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");
79 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)");
80 Optionpk<double> from_opt("from", "from", "start from this row in training file (start from 0)",0);
81 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);
82 Optionpk<double> offset_opt("\0", "offset", "offset value for each spectral band input features: refl[band]=(DN[band]-offset[band])/scale[band]", 0.0);
83 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);
84 Optionpk<unsigned short> cv_opt("cv", "cv", "n-fold cross validation mode",0);
85 Optionpk<unsigned int> nneuron_opt("nn", "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);
86 Optionpk<float> connection_opt("\0", "connection", "connection reate (default: 1.0 for a fully connected network)", 1.0);
87 // 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);
88 Optionpk<float> learning_opt("l", "learning", "learning rate (default: 0.7)", 0.7);
89 Optionpk<unsigned int> maxit_opt("\0", "maxit", "number of maximum iterations (epoch) (default: 500)", 500);
90 Optionpk<short> verbose_opt("v", "verbose", "set to: 0 (results only), 1 (confusion matrix), 2 (debug)",0,2);
91
92 offset_opt.setHide(1);
93 scale_opt.setHide(1);
94 connection_opt.setHide(1);
95 learning_opt.setHide(1);
96 maxit_opt.setHide(1);
97
98 bool doProcess;//stop process when program was invoked with help option (-h --help)
99 try{
100 doProcess=input_opt.retrieveOption(argc,argv);
101 training_opt.retrieveOption(argc,argv);
102 inputCols_opt.retrieveOption(argc,argv);
103 outputCols_opt.retrieveOption(argc,argv);
104 output_opt.retrieveOption(argc,argv);
105 from_opt.retrieveOption(argc,argv);
106 to_opt.retrieveOption(argc,argv);
107 cv_opt.retrieveOption(argc,argv);
108 nneuron_opt.retrieveOption(argc,argv);
109 offset_opt.retrieveOption(argc,argv);
110 scale_opt.retrieveOption(argc,argv);
111 connection_opt.retrieveOption(argc,argv);
112 // weights_opt.retrieveOption(argc,argv);
113 learning_opt.retrieveOption(argc,argv);
114 maxit_opt.retrieveOption(argc,argv);
115 verbose_opt.retrieveOption(argc,argv);
116 }
117 catch(string predefinedString){
118 std::cout << predefinedString << std::endl;
119 exit(0);
120 }
121 if(!doProcess){
122 cout << endl;
123 cout << "Usage: pkregann -i input -t training [-ic col]* [-oc col]* -o output" << endl;
124 cout << endl;
125 std::cout << "short option -h shows basic options only, use long option --help to show all options" << std::endl;
126 exit(0);//help was invoked, stop processing
127 }
128
129 unsigned int ninput=inputCols_opt.size();
130 unsigned int noutput=outputCols_opt.size();
131 assert(ninput);
132 assert(noutput);
133 vector< vector<float> > inputUnits;
134 vector< vector<float> > trainingUnits;
135 vector< vector<float> > trainingOutput;
136 FileReaderAscii inputFile;
137 unsigned int inputSize=0;
138 if(input_opt.size()){
139 inputFile.open(input_opt[0]);
140 inputFile.setMinRow(from_opt[0]);
141 inputFile.setMaxRow(to_opt[0]);
142 inputFile.setComment('#');
143 inputFile.readData(inputUnits,inputCols_opt,1,0,true,verbose_opt[0]);
144 inputFile.close();
145 inputSize=inputUnits.size();
146 }
147 FileReaderAscii trainingFile(training_opt[0]);
148 unsigned int sampleSize=0;
149 trainingFile.setMinRow(from_opt[0]);
150 trainingFile.setMaxRow(to_opt[0]);
151 trainingFile.setComment('#');
152 trainingFile.readData(trainingUnits,inputCols_opt,1,0,true,verbose_opt[0]);
153 trainingFile.readData(trainingOutput,outputCols_opt,1,0,true,verbose_opt[0]);
154 trainingFile.close();
155 sampleSize=trainingUnits.size();
156
157 if(verbose_opt[0]>1){
158 std::cout << "sampleSize: " << sampleSize << std::endl;
159 std::cout << "ninput: " << ninput << std::endl;
160 std::cout << "noutput: " << noutput << std::endl;
161 std::cout << "trainingUnits[0].size(): " << trainingUnits[0].size() << std::endl;
162 std::cout << "trainingOutput[0].size(): " << trainingOutput[0].size() << std::endl;
163 std::cout << "trainingUnits.size(): " << trainingUnits.size() << std::endl;
164 std::cout << "trainingOutput.size(): " << trainingOutput.size() << std::endl;
165 }
166
167 assert(ninput==trainingUnits[0].size());
168 assert(noutput==trainingOutput[0].size());
169 assert(trainingUnits.size()==trainingOutput.size());
170
171 //set scale and offset
172 if(offset_opt.size()>1)
173 assert(offset_opt.size()==ninput);
174 if(scale_opt.size()>1)
175 assert(scale_opt.size()==ninput);
176
177 std::vector<float> offset_input(ninput);
178 std::vector<float> scale_input(ninput);
179
180 std::vector<float> offset_output(noutput);
181 std::vector<float> scale_output(noutput);
182
183 for(int iinput=0;iinput<ninput;++iinput){
184 if(verbose_opt[0]>=1)
185 cout << "scaling for input feature" << iinput << endl;
186 offset_input[iinput]=(offset_opt.size()==1)?offset_opt[0]:offset_opt[iinput];
187 scale_input[iinput]=(scale_opt.size()==1)?scale_opt[0]:scale_opt[iinput];
188 //search for min and maximum
189 if(scale_input[iinput]<=0){
190 float theMin=trainingUnits[0][iinput];
191 float theMax=trainingUnits[0][iinput];
192 for(int isample=0;isample<trainingUnits.size();++isample){
193 if(theMin>trainingUnits[isample][iinput])
194 theMin=trainingUnits[isample][iinput];
195 if(theMax<trainingUnits[isample][iinput])
196 theMax=trainingUnits[isample][iinput];
197 }
198 offset_input[iinput]=theMin+(theMax-theMin)/2.0;
199 scale_input[iinput]=(theMax-theMin)/2.0;
200 if(verbose_opt[0]>=1){
201 std::cout << "Extreme image values for input feature " << iinput << ": [" << theMin << "," << theMax << "]" << std::endl;
202 std::cout << "Using offset, scale: " << offset_input[iinput] << ", " << scale_input[iinput] << std::endl;
203 std::cout << "scaled values for input feature " << iinput << ": [" << (theMin-offset_input[iinput])/scale_input[iinput] << "," << (theMax-offset_input[iinput])/scale_input[iinput] << "]" << std::endl;
204 }
205 }
206 }
207
208 for(int ioutput=0;ioutput<noutput;++ioutput){
209 if(verbose_opt[0]>=1)
210 cout << "scaling for output feature" << ioutput << endl;
211 //search for min and maximum
212 float theMin=trainingOutput[0][ioutput];
213 float theMax=trainingOutput[0][ioutput];
214 for(int isample=0;isample<trainingOutput.size();++isample){
215 if(theMin>trainingOutput[isample][ioutput])
216 theMin=trainingOutput[isample][ioutput];
217 if(theMax<trainingOutput[isample][ioutput])
218 theMax=trainingOutput[isample][ioutput];
219 }
220 offset_output[ioutput]=theMin+(theMax-theMin)/2.0;
221 scale_output[ioutput]=(theMax-theMin)/2.0;
222 if(verbose_opt[0]>=1){
223 std::cout << "Extreme image values for output feature " << ioutput << ": [" << theMin << "," << theMax << "]" << std::endl;
224 std::cout << "Using offset, scale: " << offset_output[ioutput] << ", " << scale_output[ioutput] << std::endl;
225 std::cout << "scaled values for output feature " << ioutput << ": [" << (theMin-offset_output[ioutput])/scale_output[ioutput] << "," << (theMax-offset_output[ioutput])/scale_output[ioutput] << "]" << std::endl;
226 }
227 }
228
229
230
231 FANN::neural_net net;//the neural network
232
233 const unsigned int num_layers = nneuron_opt.size()+2;
234 const float desired_error = 0.0003;
235 const unsigned int iterations_between_reports = (verbose_opt[0])? maxit_opt[0]+1:0;
236 if(verbose_opt[0]>=1){
237 cout << "creating artificial neural network with " << nneuron_opt.size() << " hidden layer, having " << endl;
238 for(int ilayer=0;ilayer<nneuron_opt.size();++ilayer)
239 cout << nneuron_opt[ilayer] << " ";
240 cout << "neurons" << endl;
241 }
242
243 switch(num_layers){
244 case(3):{
245 unsigned int layers[3];
246 layers[0]=ninput;
247 layers[1]=nneuron_opt[0];
248 layers[2]=noutput;
249 net.create_sparse_array(connection_opt[0],num_layers,layers);
250 // net.create_sparse(connection_opt[0],num_layers, ninput, nneuron_opt[0], noutput);
251 break;
252 }
253 case(4):{
254 unsigned int layers[3];
255 layers[0]=ninput;
256 layers[1]=nneuron_opt[0];
257 layers[2]=nneuron_opt[1];
258 layers[3]=noutput;
259 net.create_sparse_array(connection_opt[0],num_layers,layers);
260 // net.create_sparse(connection_opt[0],num_layers, ninput, nneuron_opt[0], nneuron_opt[1], noutput);
261 break;
262 }
263 default:
264 cerr << "Only 1 or 2 hidden layers are supported!" << endl;
265 exit(1);
266 break;
267 }
268 if(verbose_opt[0]>=1)
269 cout << "network created" << endl;
270
271 net.set_learning_rate(learning_opt[0]);
272
273 // net.set_activation_steepness_hidden(1.0);
274 // net.set_activation_steepness_output(1.0);
275
276 net.set_activation_function_hidden(FANN::SIGMOID_SYMMETRIC_STEPWISE);
277 net.set_activation_function_output(FANN::SIGMOID_SYMMETRIC_STEPWISE);
278
279 // Set additional properties such as the training algorithm
280 // net.set_training_algorithm(FANN::TRAIN_QUICKPROP);
281
282 // Output network type and parameters
283 if(verbose_opt[0]>=1){
284 cout << endl << "Network Type : ";
285 switch (net.get_network_type())
286 {
287 case FANN::LAYER:
288 cout << "LAYER" << endl;
289 break;
290 case FANN::SHORTCUT:
291 cout << "SHORTCUT" << endl;
292 break;
293 default:
294 cout << "UNKNOWN" << endl;
295 break;
296 }
297 net.print_parameters();
298 }
299
300 if(verbose_opt[0]>=1){
301 cout << "Max Epochs " << setw(8) << maxit_opt[0] << ". "
302 << "Desired Error: " << left << desired_error << right << endl;
303 }
304 bool initWeights=true;
305
306 Vector2d<float> trainingFeatures(sampleSize,ninput);
307 for(unsigned int isample=0;isample<sampleSize;++isample){
308 for(unsigned int iinput=0;iinput<ninput;++iinput)
309 trainingFeatures[isample][iinput]=(trainingUnits[isample][iinput]-offset_input[iinput])/scale_input[iinput];
310 }
311
312 Vector2d<float> scaledOutput(sampleSize,noutput);
313 for(unsigned int isample=0;isample<sampleSize;++isample){
314 for(unsigned int ioutput=0;ioutput<noutput;++ioutput)
315 scaledOutput[isample][ioutput]=(trainingOutput[isample][ioutput]-offset_output[ioutput])/scale_output[ioutput];
316 }
317
318 if(cv_opt[0]){
319 if(verbose_opt[0])
320 std::cout << "cross validation" << std::endl;
321 std::vector< std::vector<float> > referenceVector;
322 std::vector< std::vector<float> > outputVector;
323 net.cross_validation(trainingFeatures,
324 scaledOutput,
325 cv_opt[0],
326 maxit_opt[0],
327 desired_error,
328 referenceVector,
329 outputVector);
330 assert(referenceVector.size()==outputVector.size());
331 vector<double> rmse(noutput);
332 for(int isample=0;isample<referenceVector.size();++isample){
333 std::cout << isample << " ";
334 for(int ioutput=0;ioutput<noutput;++ioutput){
335 if(!isample)
336 rmse[ioutput]=0;
337 double ref=scale_output[ioutput]*referenceVector[isample][ioutput]+offset_output[ioutput];
338 double val=scale_output[ioutput]*outputVector[isample][ioutput]+offset_output[ioutput];
339 rmse[ioutput]+=(ref-val)*(ref-val);
340 std::cout << ref << " " << val;
341 if(ioutput<noutput-1)
342 std::cout << " ";
343 else
344 std::cout << std::endl;
345 }
346 }
347 for(int ioutput=0;ioutput<noutput;++ioutput)
348 std::cout << "rmse output variable " << ioutput << ": " << sqrt(rmse[ioutput]/referenceVector.size()) << std::endl;
349 }
350
351
352 net.train_on_data(trainingFeatures,
353 scaledOutput,
354 initWeights,
355 maxit_opt[0],
356 iterations_between_reports,
357 desired_error);
358
359
360 if(verbose_opt[0]>=2){
361 net.print_connections();
362 vector<fann_connection> convector;
363 net.get_connection_array(convector);
364 for(int i_connection=0;i_connection<net.get_total_connections();++i_connection)
365 cout << "connection " << i_connection << ": " << convector[i_connection].weight << endl;
366 }
367 //end of training
368
369 ofstream outputStream;
370 if(!output_opt.empty())
371 outputStream.open(output_opt[0].c_str(),ios::out);
372 for(unsigned int isample=0;isample<inputUnits.size();++isample){
373 std::vector<float> inputFeatures(ninput);
374 for(unsigned int iinput=0;iinput<ninput;++iinput)
375 inputFeatures[iinput]=(inputUnits[isample][iinput]-offset_input[iinput])/scale_input[iinput];
376 vector<float> result(noutput);
377 result=net.run(inputFeatures);
378
379 if(!output_opt.empty())
380 outputStream << isample << " ";
381 else
382 std::cout << isample << " ";
383 if(verbose_opt[0]){
384 for(unsigned int iinput=0;iinput<ninput;++iinput){
385 if(output_opt.size())
386 outputStream << inputUnits[isample][iinput] << " ";
387 else
388 std::cout << inputUnits[isample][iinput] << " ";
389 }
390 }
391 for(unsigned int ioutput=0;ioutput<noutput;++ioutput){
392 result[ioutput]=scale_output[ioutput]*result[ioutput]+offset_output[ioutput];
393 if(output_opt.size()){
394 outputStream << result[ioutput];
395 if(ioutput<noutput-1)
396 outputStream << " ";
397 else
398 outputStream << std::endl;
399 }
400 else{
401 std::cout << result[ioutput];
402 if(ioutput<noutput-1)
403 std::cout << " ";
404 else
405 std::cout << std::endl;
406 }
407 }
408 }
409 if(!output_opt.empty())
410 outputStream.close();
411}