## Author: Jan
## Created: 2016-11-09

function [w1, w2, w3, errors] = trainNeuralNetworkA (inputFile, eta)
  windowSize = 7; #because of indexing, window size will be 8
  input = dlmread(inputFile);
  records = size(input)(1, 1);
  
  inputLayerSize = 41; #5 x 8 +1 (input types x window size + treshold input)
  hiddenLayerSize1 = 33;
  hiddenLayerSize2 = 9;
  
  intervalM = 0.5;
  intervalA = 0.25;
  
  #initialize weights to random values between -0.25 and 0.25
  w1 = rand(inputLayerSize, hiddenLayerSize1)*intervalM-intervalA;
  w2 = rand(hiddenLayerSize1, hiddenLayerSize2)*intervalM-intervalA;
  w3 = rand(hiddenLayerSize2, 1)*intervalM-intervalA;  
  
  errors = 0;
  for i=1:80
      for j=1:(records-(windowSize+3)) #prediction 3 hours
        inp = [input(j:j+windowSize, 4)', input(j:j+windowSize, 5)', input(j:j+windowSize, 6)', input(j:j+windowSize, 7)', input(j:j+windowSize, 8)', -1];
        #select 8 rows of ScalarB, GSE, GSM, RMS_BZ_GSE, Proton Density
        exp = input(j+windowSize+3, 9); #expected output
        
        HHidden1 = inp*w1;
        VHidden1 = activationFunction(HHidden1);
        VHidden1(1, hiddenLayerSize1) = -1;#treshold input
        
        HHidden2 = VHidden1*w2;
        VHidden2 = activationFunction(HHidden2);
        VHidden2(1, hiddenLayerSize2) = -1;
        
        HOutput = VHidden2*w3;
        VOutput = activationFunction(HOutput);
        
        DOutput = activationFunctionDerivation(HOutput).*(exp-VOutput);
        DHidden2 = activationFunctionDerivation(HHidden2).*(w3*DOutput')';
        DHidden1 = activationFunctionDerivation(HHidden1).*(w2*DHidden2')';
        
        w1 = w1+(eta*inp'*DHidden1);
        w2 = w2+(eta*VHidden1'*DHidden2);        
        w3 = w3+(eta*VHidden2'*DOutput);       
        
        errors(j+windowSize+3) = exp-VOutput;
    endfor
  endfor   
endfunction
