docs: private void trainNN(DataSet dataSet)
      {
          //batchSize
          
          List<List<Vec>> activations = new ArrayList<List<Vec>>(batchSize);
          List<List<Vec>> derivatives = new ArrayList<List<Vec>>(batchSize);
          List<List<Vec>> deltas = new ArrayList<List<Vec>>(batchSize);
          
          List<Matrix> updates = new ArrayList<Matrix>(Ws.size());
          
          List<Vec> cur_x = new ArrayList<Vec>(batchSize);
          List<Vec> prev_x = new ArrayList<Vec>(batchSize);
          
          for(int i = 0; i < batchSize; i++)
          {
              activations.add(new ArrayList<Vec>(Ws.size()));
              derivatives.add(new ArrayList<Vec>(Ws.size()));
              deltas.add(new ArrayList<Vec>(Ws.size()));
              
              for(Matrix w : Ws)
              {
                  int L = w.rows();
                  activations.get(i).add(new DenseVector(L));
                  derivatives.get(i).add(new DenseVector(L));
                  deltas.get(i).add(new DenseVector(L));
                  if(i == 0)
                      updates.add(new DenseMatrix(w.rows(), w.cols()));
              }
          }
          
          IntList iterOrder = new IntList(dataSet.size());
          ListUtils.addRange(iterOrder, 0, dataSet.size(), 1);
          
          final double bSizeInv = 1.0/batchSize;
          
          for(int epoch = 0; epoch < epochs; epoch++)
          {
              Collections.shuffle(iterOrder);
              final double eta = learningRateDecay.rate(epoch, epochs, initialLearningRate);//learningRate;
              double error = 0.0;
              for(int iter = 0; iter < dataSet.size(); iter+=batchSize)
              {
                  if(dataSet.size() - iter < batchSize)
                      continue;//we have run out of enough sampels to do an update
                  
                  cur_x.clear();
                  
                  //Feed batches thought network and get final mistakes
                  for(int bi = 0; bi < batchSize; bi++)
                  {
                      final int idx = iterOrder.get(iter+bi);
                      Vec x = dataSet.getDataPoint(idx).getNumericalValues();
                      cur_x.add(x);
                      feedForward(x, activations.get(bi), derivatives.get(bi));
  
                      
                      //Compution of Deltas
                      Vec delta_out = deltas.get(bi).get(npl.length);
  
                      Vec a_i = activations.get(bi).get(npl.length);
                      Vec d_i = derivatives.get(bi).get(npl.length);
                      
                      error += computeOutputDelta(dataSet, idx, delta_out, a_i, d_i);
                  }
                  
                  //Propigate the collected errors back
                  for(int bi = 0; bi < batchSize; bi++)
                  {
                      for(int i = Ws.size()-2; i >= 0; i--)
                      {
                          Vec delta = deltas.get(bi).get(i);
                          delta.zeroOut();
                          Matrix W = Ws.get(i+1);
                          W.transposeMultiply(1, deltas.get(bi).get(i+1), delta);
                          delta.mutablePairwiseMultiply(derivatives.get(bi).get(i));
                      }
  
                      //Apply weight changes
                      for(int i = 1; i < Ws.size(); i++)
                      {
                          Matrix W = Ws.get(i);
                          Vec b = bs.get(i);
                          W.mutableSubtract(eta*weightDecay, W);
                          
                          if(momentum != 0)
                          {
                              Matrix update = updates.get(i);
                              update.mutableMultiply(momentum);
                              Matrix.OuterProductUpdate(update, deltas.get(bi).get(i), activations.get(bi).get(i-1), -eta*bSizeInv);
                              W.mutableAdd(update);
                          }
                          else//update directly
                          {
                              Matrix.OuterProductUpdate(W, deltas.get(bi).get(i), activations.get(bi).get(i-1), -eta*bSizeInv);
                          }
                          
                          b.mutableAdd(-eta*bSizeInv, deltas.get(bi).get(i));
                      }
  
                      //input layer
                      Matrix W = Ws.get(0);
                      W.mutableSubtract(eta*weightDecay, W);
                      Vec b = bs.get(0);
                      
                      if(momentum != 0)
                      {
                          Matrix update = updates.get(0);
                          update.mutableMultiply(momentum);
  
                          Matrix.OuterProductUpdate(update, deltas.get(bi).get(0), cur_x.get(bi), -eta*bSizeInv);
                          W.mutableAdd(update);
                      }
                      else//update directly
                      {
                          Matrix.OuterProductUpdate(W, deltas.get(bi).get(0), cur_x.get(bi), -eta*bSizeInv);
                      }
                          
                      b.mutableAdd(-eta*bSizeInv, deltas.get(bi).get(0));
                  }
                  
              }
          }
      }