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

function [w1, w2, w3, errors] = trainNeuralNetworkB (inputFile, eta)
  windowSize = 7; #because of indexing, window size will be 8
  input = dlmread(inputFile);
  records = size(input)(1, 1);
  storms = records/146;
  
  inputLayerSize = 49; #5 x 8 + 8 (input types x window size)
  #special input from generated values = 8
  hiddenLayerSize1 = 33;
  hiddenLayerSize2 = 15;
  
  intervalM = 1;
  intervalA = 0.5;
  
  #initialize weights to random values between -0.5 and 0.5
  w1 = rand(inputLayerSize, hiddenLayerSize1)*intervalM-intervalA;
  w2 = rand(hiddenLayerSize1, hiddenLayerSize2)*intervalM-intervalA;
  w3 = rand(hiddenLayerSize2, 1)*intervalM-intervalA;
  
  errors = 0;
  outputs = zeros(1, records);
  outputs(1, 1:10) = input(1:10, 9);
  #I consider this step neccessary, because network is working with previously generated
  #values, but first 10 values (8 window + 3 prediction) are zeros, and this can make traoubles
  for i=1:100
      for j=1:storms;#(records-(windowSize+3)) #prediction 3 hours               
        for k=1:146-(windowSize+3)
          c = (j-1)*146+k;
          inp = [input(c:c+windowSize, 4)', input(c:c+windowSize, 5)', input(c:c+windowSize, 6)', input(c:c+windowSize, 7)', input(c:c+windowSize, 8)', outputs(1, c:c+windowSize), -1];
          #select 8 rows of ScalarB, GSE, GSM, RMS_BZ_GSE, Proton Density, Dst(predicted)
          exp = input(c+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(c+windowSize+3) = exp-VOutput;
          outputs(c+windowSize+3) = VOutput;#for next iteration input
      endfor
    endfor
  endfor   
endfunction
