Index: trunk/Ohana/src/opihi/cmd.data/nnet.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/nnet.c	(revision 40336)
+++ trunk/Ohana/src/opihi/cmd.data/nnet.c	(revision 40371)
@@ -4,8 +4,10 @@
   {1, "list",     nnet_list,     "list nnets"},
   {1, "delete",   nnet_delete,   "delete a nnet"},
-  // {1, "show",     nnet_show,     "display nnet values"},
+  {1, "print",    nnet_print,    "display nnet values"},
   {1, "create",   nnet_create,   "create a nnet"},
   {1, "set",      nnet_set,      "set nnet node values"},
   {1, "get",      nnet_get,      "get nnet node values"},
+  {1, "read",     nnet_read,     "read nnet values from a file"},
+// {1, "write",    nnet_write,    "write nnet values to a file"},
   {1, "train",    nnet_train,    "train nnet on a set of data"},
   {1, "apply",    nnet_apply,    "apply nnet to a set of data"},
@@ -20,8 +22,9 @@
     gprint (GP_ERR, "    nnet list          : list nnets\n");
     gprint (GP_ERR, "    nnet delete (nnet) : delete a nnet\n");
-    gprint (GP_ERR, "    nnet show   (nnet) : show values for a nnet\n");
+    gprint (GP_ERR, "    nnet print  (nnet) : print values for a nnet\n");
     gprint (GP_ERR, "    nnet create (nnet) (Ninput) [Nnodes] [Nnodes] ... (Noutput) : create a nnet\n");
     gprint (GP_ERR, "    nnet set    (nnet) [weights] [biases] ... [weights] [biases] : set nnet weights (images) and biases (vectors)\n");
     gprint (GP_ERR, "    nnet get    (nnet) [weights] [biases] ... [weights] [biases] : get nnet weights (images) and biases (vectors)\n");
+    gprint (GP_ERR, "    nnet read   (nnet) (filename) : set nnet weights and biases using a data file\n");
     gprint (GP_ERR, "    nnet train  (nnet) [input] [input] ... [output] [output] ... : train nnet on data from a set of vectors\n");
     gprint (GP_ERR, "    nnet apply  (nnet) [input] [input] ... [output] [output] ... : apply nnet to input data and generate output\n");
Index: trunk/Ohana/src/opihi/cmd.data/nnet_commands.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/nnet_commands.c	(revision 40336)
+++ trunk/Ohana/src/opihi/cmd.data/nnet_commands.c	(revision 40371)
@@ -62,4 +62,21 @@
   CreateNnetData (nnet, LargeWeightInit);
 
+  return TRUE;
+}
+
+int nnet_print (int argc, char **argv) {
+
+  if (argc != 2) {
+    gprint (GP_ERR, "USAGE: nnet print (nnet)\n");
+    return FALSE;
+  }
+
+  Nnet *nnet = FindNnet (argv[1]);
+  if (nnet == NULL) {
+    gprint (GP_ERR, "nnet %s not found, create it first\n", argv[1]);
+    return FALSE;
+  }
+
+  PrintNnet (nnet);
   return TRUE;
 }
@@ -146,4 +163,99 @@
 }
 
+int nnet_read (int argc, char **argv) {
+
+  # define D_LINE 0x10000
+  char word[128];
+  char line[D_LINE];
+
+  if (argc < 3) {
+    gprint (GP_ERR, "USAGE: nnet read (nnet) (filename) : set nnet weights and biases based on a file\n");
+    gprint (GP_ERR, "the first line of the file specifies the number of layers, each layer in the network is written as a matrix of numbers (weights) and a vector\n");
+    return FALSE;
+  }
+
+  // open the file and read the number of layers
+  FILE *f = fopen (argv[2], "r");
+  if (f == NULL) {
+    gprint (GP_ERR, "USAGE: nnet read (nnet) (filename) : set nnet weights and biases based on a file\n");
+    gprint (GP_ERR, "file %s could not be opened\n", argv[2]);
+    return FALSE;
+  }
+
+  // read the number of layers
+  // NLAYER 3
+  int Nlayer;
+  scan_line_maxlen (f, line, D_LINE);
+  sscanf (line, "%127s %d", word, &Nlayer);
+  if (strcmp(word, "NLAYER")) {
+    gprint (GP_ERR, "warning: NLAYER keyword not found\n");
+  }
+
+  Nnet *nnet = CreateNnet (argv[1], Nlayer);
+
+  // read the number of nodes
+  // LAYERS 2 4 2
+  scan_line_maxlen (f, line, D_LINE);
+  char *tmpword = getword (line);
+  if (strcmp (tmpword, "LAYERS")) {
+    gprint (GP_ERR, "warning: LAYERS keyword not found\n");
+  }
+  FREE (tmpword);
+  for (int i = 0; i < Nlayer; i++) {
+    int Nnode;
+    int status = iparse (&Nnode, i + 2, line); // numbering is fields 1 2 3
+    if (!status) {
+      gprint (GP_ERR, "error: failed to find all Nnode values\n");
+      fclose (f);
+      DeleteNnet (nnet);
+      return FALSE;
+    }
+    nnet[0].Nnodes[i] = Nnode;
+  }
+
+  // this creates the data for each node and inits with gaussian weights
+  CreateNnetData (nnet, FALSE);
+
+  // read the weights and biases from the data file
+  for (int L = 1; L < nnet[0].Nlayer; L++) {
+    
+    char word1[128], word2[128], word3[128];
+    int Nx, Ny, layer;
+    scan_line_maxlen (f, line, D_LINE);
+    sscanf (line, "%127s %d %127s %d %127s %d", word1, &layer, word2, &Nx, word3, &Ny);
+    gprint (GP_ERR, "LAYER: %d, Nx: %d, Ny: %d\n", layer, Nx, Ny);
+
+    if (layer != L - 1) {
+      gprint (GP_ERR, "warning: expect layer = %d, got %d\n", L - 1, layer);
+    }
+    if (Nx != nnet[0].Nnodes[L - 1]) {
+      gprint (GP_ERR, "warning: expect Nx = %d, got %d\n", nnet[0].Nnodes[L - 1], Nx);
+    }
+    if (Ny != nnet[0].Nnodes[L]) {
+      gprint (GP_ERR, "warning: expect Ny = %d, got %d\n", nnet[0].Nnodes[L], Ny);
+    }
+
+    double value;
+    for (int j = 0; j < Ny; j++) {
+      scan_line_maxlen (f, line, D_LINE);
+      for (int i = 0; i < Nx; i++) {
+	int k = i + j*nnet[0].Nnodes[L-1];
+	int status = dparse (&value, i + 1, line);
+	if (!status) {
+	  gprint (GP_ERR, "error: failed to read entry from a line (layer: %d, i: %d, j: %d)\n", layer, i, j);
+	  fclose (f);
+	  DeleteNnet (nnet);
+	  return FALSE;
+	}
+	nnet[0].weight[L][k] = value;
+      }
+      dparse (&value, Nx + 1, line);
+      nnet[0].biases[L][j] = value;
+    }
+  }
+
+  return TRUE;
+}
+
 int nnet_get (int argc, char **argv) {
 
Index: trunk/Ohana/src/opihi/cmd.data/nnet_train.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/nnet_train.c	(revision 40336)
+++ trunk/Ohana/src/opihi/cmd.data/nnet_train.c	(revision 40371)
@@ -10,4 +10,6 @@
 void  nnet_backprop (Nnet *nnet, Vector **inVec, Vector **outVec, int N);
 void  nnet_descent_step (Nnet *nnet, Vector **inVec, Vector **outVec, int *seq, int pass, int Nmini, float eta, float lambda);
+void  nnet_print_Nabla (Nnet *nnet);
+void  nnet_write_Nabla (char *filename, Nnet *nnet);
 
 static int QUADRATIC_COST = 0;
@@ -76,5 +78,5 @@
     gprint (GP_ERR, "USAGE: nnet train (nnet) [input] [input] ... [output] [output] ...\n");
     gprint (GP_ERR, "OPTIONS: -Nepoch [N] -Nmini [N]\n");
-    FREE (resid);
+    // FREE (resid);
     return FALSE;
   }
@@ -83,5 +85,5 @@
   if (nnet == NULL) {
     gprint (GP_ERR, "nnet %s not found, create it first\n", argv[1]);
-    FREE (resid);
+    // FREE (resid);
     return FALSE;
   }
@@ -95,5 +97,5 @@
   if (argc != Ninput + Noutput + 2) {
     gprint (GP_ERR, "need %d input and %d output vectors, but we have %d total\n", nnet[0].Nnodes[0], nnet[0].Nnodes[Nlayer - 1], argc - 2);
-    FREE (resid);
+    // FREE (resid);
     return FALSE;
   }
@@ -111,5 +113,5 @@
       free (inVec);
       free (outVec);
-      FREE (resid);
+      // FREE (resid);
       return FALSE;    
     }
@@ -118,5 +120,5 @@
       free (inVec);
       free (outVec);
-      FREE (resid);
+      // FREE (resid);
       return FALSE;    
     }
@@ -126,5 +128,5 @@
       free (inVec);
       free (outVec);
-      FREE (resid);
+      // FREE (resid);
       return FALSE;    
     }
@@ -136,5 +138,5 @@
       free (inVec);
       free (outVec);
-      FREE (resid);
+      // FREE (resid);
       return FALSE;    
     }
@@ -143,5 +145,5 @@
       free (inVec);
       free (outVec);
-      FREE (resid);
+      // FREE (resid);
       return FALSE;    
     }
@@ -181,4 +183,6 @@
   }
 
+  // PrintNnet (nnet);
+
   // train for Nepochs
   // this recreates 'SGD' from http://neuralnetworksanddeeplearning.com/chap1.html
@@ -194,5 +198,7 @@
       // update the weights and biases using the mini batch subset
       nnet_descent_step (nnet, inVec, outVec, seq, pass, Nmini, eta, lambda);
-    }
+      // return TRUE; // XXX short-circuit at one step
+    }
+    // PrintNnet (nnet);
 
     if (resid) {
@@ -218,10 +224,13 @@
       float mean = s1 / Npts;
       float sigma = sqrt(s2 / Npts - mean*mean);
-      if (epoch % 10 == 0) gprint (GP_ERR, "epoch %d of %d, %f +/- %f\n", epoch, Nepoch, mean, sigma);
+      // if (epoch % 10 == 0) gprint (GP_ERR, "epoch %d of %d, %f +/- %f\n", epoch, Nepoch, mean, sigma);
+      gprint (GP_ERR, "epoch %d of %d, %f +/- %f\n", epoch, Nepoch, mean, sigma);
       resid[0].elements.Flt[epoch] = sigma;
     } else {
-      if (epoch % 10 == 0) gprint (GP_ERR, "epoch %d of %d\n", epoch, Nepoch);
+      // if (epoch % 10 == 0) gprint (GP_ERR, "epoch %d of %d\n", epoch, Nepoch);
+      gprint (GP_ERR, "epoch %d of %d\n", epoch, Nepoch);
     }
   }  
+  // PrintNnet (nnet);
 
   if (result) {
@@ -265,13 +274,25 @@
 
     // N is the element of the mini batch on which we are currently operating
-    int N = seq[pass*Nmini + i];
+    // int N = seq[pass*Nmini + i]; // XXX uncomment to turn on random shuffle
+    int N = pass*Nmini + i;
 
     // backprop generates a dNabla_b, dNabla_w pair for the element N of the input and output vectors
     nnet_backprop (nnet, inVec, outVec, N);
-
     nnet_update_Nabla (nnet);
-  }
-
+    // gprint (GP_ERR, ". ");
+
+    // nnet_print_Nabla (nnet); // XXX print nabla for each epoch
+    // XXX uncomment to dump nablas after one step, one element
+    // nnet_write_Nabla ("test.nabla.op.dat", nnet); // XXX print nabla for each epoch
+    // return; 
+  }
+  // gprint (GP_ERR, " done mini batch\n");
+
+  // nnet_print_Nabla (nnet);
   nnet_apply_Nabla (nnet, Nmini, eta, lambda, Ntrial);
+
+  // XXX uncomment to dump nablas after one mini batch
+  // nnet_write_Nabla ("test.nabla.op.dat", nnet); // XXX print nabla for each epoch
+  // PrintNnet (nnet);
 }
 
@@ -317,4 +338,14 @@
       }
     } else {
+      // XXX TEST PRINTS to catch code errors compared to python implementation
+      // gprint (GP_ERR, "z: ");
+      // for (int j = 0; j < nnet[0].Nnodes[L]; j++) {
+      // 	gprint (GP_ERR, "%f ", 	nnet[0].zvalue[L][j]);
+      // } gprint (GP_ERR, "\n");
+      // gprint (GP_ERR, "sp: ");
+      // for (int j = 0; j < nnet[0].Nnodes[L]; j++) {
+      // 	gprint (GP_ERR, "%f ", 	nnet[0].sprime[L][j]);
+      // } gprint (GP_ERR, "\n");
+
       // delta = DOT(delta, transpose(weight[L+1])) * sprime;
       for (int j = 0; j < nnet[0].Nnodes[L]; j++) {
@@ -322,5 +353,7 @@
 	for (int i = 0; i < nnet[0].Nnodes[L+1]; i++) {
 	  int k = j + i*nnet[0].Nnodes[L]; // note order of (i,j) : j is [L+1] direction 
-	  tmpdelta += nnet[0].weight[L][k] * nnet[0].delta[L+1][i];
+	  myAssert (k < nnet[0].Nnodes[L]*nnet[0].Nnodes[L+1], "overflow");
+	  tmpdelta += nnet[0].weight[L+1][k] * nnet[0].delta[L+1][i];
+	  // gprint (GP_ERR, "%e %e\n", nnet[0].weight[L+1][k], nnet[0].delta[L+1][i]);
 	}
 	nnet[0].delta[L][j] = tmpdelta * nnet[0].sprime[L][j];
@@ -335,5 +368,5 @@
     // Nabla_b[L] = delta;
     for (int j = 0; j < nnet[0].Nnodes[L]; j++) {
-      nnet[0]. Nabla_b[L][j] = nnet[0].delta[L][j];
+      nnet[0]. dNabla_b[L][j] = nnet[0].delta[L][j];
     }
     
@@ -342,5 +375,6 @@
       for (int i = 0; i < nnet[0].Nnodes[L-1]; i++) {
 	int k = i + j*nnet[0].Nnodes[L-1];
-	nnet[0]. Nabla_w[L][k] = nnet[0].svalue[L-1][i] * nnet[0].delta[L][j];
+	myAssert (k < nnet[0].Nnodes[L-1]*nnet[0].Nnodes[L], "overflow");
+	nnet[0]. dNabla_w[L][k] = nnet[0].svalue[L-1][i] * nnet[0].delta[L][j];
       }
     }
@@ -349,4 +383,35 @@
 
 // support functions to loop over the Nabla entries
+void nnet_reset_Nabla (Nnet *nnet) {
+  for (int L = 1; L < nnet[0].Nlayer; L++) {
+
+    for (int j = 0; j < nnet[0].Nnodes[L]; j++) {
+
+      nnet[0]. Nabla_b[L][j] = 0;
+      nnet[0].dNabla_b[L][j] = 0;
+
+      for (int i = 0; i < nnet[0].Nnodes[L-1]; i++) {
+	int k = i + j*nnet[0].Nnodes[L-1];
+	myAssert (k < nnet[0].Nnodes[L-1]*nnet[0].Nnodes[L], "overflow");
+	nnet[0]. Nabla_w[L][k] = 0;
+	nnet[0].dNabla_w[L][k] = 0;
+      }
+    }
+  }
+}
+void nnet_update_Nabla (Nnet *nnet) {
+  for (int L = 1; L < nnet[0].Nlayer; L++) {
+
+    for (int j = 0; j < nnet[0].Nnodes[L]; j++) {
+      nnet[0]. Nabla_b[L][j] += nnet[0].dNabla_b[L][j];
+
+      for (int i = 0; i < nnet[0].Nnodes[L-1]; i++) {
+	int k = i + j*nnet[0].Nnodes[L-1];
+	myAssert (k < nnet[0].Nnodes[L-1]*nnet[0].Nnodes[L], "overflow");
+	nnet[0]. Nabla_w[L][k] += nnet[0].dNabla_w[L][k];
+      }
+    }
+  }
+}
 void nnet_apply_Nabla (Nnet *nnet, int Nmini, float eta, float lambda, int Ntrial) {
   for (int L = 1; L < nnet[0].Nlayer; L++) {
@@ -357,4 +422,5 @@
       for (int i = 0; i < nnet[0].Nnodes[L-1]; i++) {
 	int k = i + j*nnet[0].Nnodes[L-1];
+	myAssert (k < nnet[0].Nnodes[L-1]*nnet[0].Nnodes[L], "overflow");
 	// nnet[0].weight[L][k] -= (eta / Nmini) * nnet[0].Nabla_w[L][k];
 	// with lambda > 0.0, we have L2 regularization.  if lambda = 0.0, we recover the default implementation
@@ -364,32 +430,46 @@
   }
 }
-void nnet_update_Nabla (Nnet *nnet) {
+
+void nnet_print_Nabla (Nnet *nnet) {
+
   for (int L = 1; L < nnet[0].Nlayer; L++) {
-
-    for (int j = 0; j < nnet[0].Nnodes[L]; j++) {
-      nnet[0]. Nabla_b[L][j] += nnet[0].dNabla_b[L][j];
-
+    gprint (GP_ERR, " ----- Nabla %d -----\n", L);
+    for (int j = 0; j < nnet[0].Nnodes[L]; j++) {
       for (int i = 0; i < nnet[0].Nnodes[L-1]; i++) {
-	int k = i + j*nnet[0].Nnodes[L-1];
-	nnet[0]. Nabla_w[L][k] += nnet[0].dNabla_w[L][k];
-      }
-    }
-  }
-}
-void nnet_reset_Nabla (Nnet *nnet) {
+	int k = j * nnet[0].Nnodes[L-1] + i;
+	myAssert (k < nnet[0].Nnodes[L-1]*nnet[0].Nnodes[L], "overflow");
+	gprint (GP_ERR, "%10.3e ", nnet[0].Nabla_w[L][k]);
+      }
+      gprint (GP_ERR, " : %10.3e\n", nnet[0].Nabla_b[L][j]);
+    }
+  }
+  return;
+}
+
+void nnet_write_Nabla (char *filename, Nnet *nnet) {
+
+  FILE *f = fopen (filename, "w");
+
+  fprintf (f, "NLAYER %d\n", nnet[0].Nlayer);
+  fprintf (f, "LAYERS ");
+  for (int L = 0; L < nnet[0].Nlayer; L++) {
+    fprintf (f, "%d ", nnet[0].Nnodes[L]);
+  }
+  fprintf (f, "\n");
+
   for (int L = 1; L < nnet[0].Nlayer; L++) {
-
-    for (int j = 0; j < nnet[0].Nnodes[L]; j++) {
-
-      nnet[0]. Nabla_b[L][j] = 0;
-      nnet[0].dNabla_b[L][j] = 0;
-
+    fprintf (f, "LAYER %d NX %d NY %d\n", L - 1, nnet[0].Nnodes[L-1], nnet[0].Nnodes[L]);
+    for (int j = 0; j < nnet[0].Nnodes[L]; j++) {
       for (int i = 0; i < nnet[0].Nnodes[L-1]; i++) {
-	int k = i + j*nnet[0].Nnodes[L-1];
-	nnet[0]. Nabla_w[L][k] = 0;
-	nnet[0].dNabla_w[L][k] = 0;
-      }
-    }
-  }
+	int k = j * nnet[0].Nnodes[L-1] + i;
+	myAssert (k < nnet[0].Nnodes[L-1]*nnet[0].Nnodes[L], "overflow");
+	fprintf (f, "%.9f ", nnet[0].Nabla_w[L][k]);
+      }
+      fprintf (f, "%.9f\n", nnet[0].Nabla_b[L][j]);
+    }
+  }
+
+  fclose (f);
+  return;
 }
 
@@ -410,5 +490,5 @@
 
   // evaluating a single layer [L], n > 0, n < Nlayer:
-  int Ninput   = nnet[0].Nnodes[L - 1];
+  int Ninput   = nnet[0].Nnodes[L-1];
   int Noutput  = nnet[0].Nnodes[L];
 
@@ -419,4 +499,5 @@
       // weight matrix order is (0, 1, ... Ninput-1, Ninput, Ninput + 1, ... Ninput * Noutput - 1)
       int k = j * Ninput + i;
+      myAssert (k < Ninput*Noutput, "overflow");
       sum += nnet[0].weight[L][k]*nnet[0].svalue[L-1][i];
     }
Index: trunk/Ohana/src/opihi/cmd.data/test/nnet.sh
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/test/nnet.sh	(revision 40336)
+++ trunk/Ohana/src/opihi/cmd.data/test/nnet.sh	(revision 40371)
@@ -1,133 +1,15 @@
 
-macro test
-
-  nnet create t0 5 10 15 3
-  nnet list
-
-  nnet create t1 6 4 5
-  nnet list
-
-  nnet create t1 5 8
-  nnet list
-
-  # vlist v0 1 2 3 4 5
-  vlist v1 1 2 3 4 5 6 7 8
-
-  mcreate m1 5 8
-  
-  # generate a weight vector of length v1[]
-  vlist my 9 8 7 6 5 4 3 2
-
-  for i 0 5
-    mset m1 my -y $i
-    set my = my + $i
-  end
-  
-  # set weights and biases
-  nnet set t1 m1 v1
-
-  # get weights and biases
-  nnet get t1 M1 V1
-
-  # compare
-  set dv = v1 - V1
-  vstat dv
-
-  set dm = m1 - M1
-  stat dm
-  
-end
-
-macro test2
-
-  nnet create t1 2 2
-  nnet list
-
-  # output = sigmoid (sum(weight[i][j] * input[i]) + bias[j])
-
-  # biases of 
-  vlist v1 0.8 0.2
-
-  mcreate m1 2 2
-  m1[0][0] =  0.5
-  m1[1][0] =  0.1
-  m1[0][1] = -0.2
-  m1[1][1] =  0.8
-
-  # set weights and biases
-  nnet set t1 m1 v1
-
-  vlist vin0 1 2
-  vlist vin1 2 1
-
-  # get weights and biases
-  nnet apply t1 vin0 vin1 vout0 vout1
-  
-  vlist z0 {m1[0][0]*vin0[0] + m1[1][0]*vin1[0] + v1[0]} {m1[0][0]*vin0[1] + m1[1][0]*vin1[1] + v1[0]} 
-  vlist z1 {m1[0][1]*vin0[0] + m1[1][1]*vin1[0] + v1[1]} {m1[0][1]*vin0[1] + m1[1][1]*vin1[1] + v1[1]}
-  set t0 = 1 / (1 + exp(-1*z0))
-  set t1 = 1 / (1 + exp(-1*z1))
-
-  # compare
-  set dt0 = t0 - vout0
-  vstat dt0
-
-  set dt1 = t1 - vout1
-  vstat dt1
-end
-
-macro test3
-
-  memory all
-
-  nnet create t0 5 10 15 3
-  nnet list
-
-  nnet create t1 5 8
-  nnet list
-
-  # vlist v0 1 2 3 4 5
-  vlist v1 1 2 3 4 5 6 7 8
-
-  mcreate m1 5 8
-  
-  # generate a weight vector of length v1[]
-  vlist my 9 8 7 6 5 4 3 2
-
-  for i 0 5
-    mset m1 my -y $i
-    set my = my + $i
-  end
-  
-  # set weights and biases
-  nnet set t1 m1 v1
-
-  # get weights and biases
-  nnet get t1 M1 V1
-
-  # compare
-  set dv = v1 - V1
-  vstat dv
-
-  set dm = m1 - M1
-  stat dm
-  
-  nnet delete t0
-  nnet delete t1
-
-  delete m1 M1 dm
-  delete v1 V1 dv my
-
-  vectors
-  buffers
-
-  memory all
-end
-
-macro test.simple
+macro test.python
   if ($0 != 4)
-    echo "USAGE: test.simple (Nmini) (Nepoch) (eta)"
+    echo "USAGE: test.python (Nmini) (Nepoch) (eta)"
     break
   end
+
+  # test I/O
+
+  create x 0 100
+  set y = dsin(x/10)
+  write test.dat x y
+  break
 
   local Nmini Nepoch eta
@@ -187,4 +69,191 @@
 end
 
+macro test
+
+  nnet create t0 5 10 15 3
+  nnet list
+
+  nnet create t1 6 4 5
+  nnet list
+
+  nnet create t1 5 8
+  nnet list
+
+  # vlist v0 1 2 3 4 5
+  vlist v1 1 2 3 4 5 6 7 8
+
+  mcreate m1 5 8
+  
+  # generate a weight vector of length v1[]
+  vlist my 9 8 7 6 5 4 3 2
+
+  for i 0 5
+    mset m1 my -y $i
+    set my = my + $i
+  end
+  
+  # set weights and biases
+  nnet set t1 m1 v1
+
+  # get weights and biases
+  nnet get t1 M1 V1
+
+  # compare
+  set dv = v1 - V1
+  vstat dv
+
+  set dm = m1 - M1
+  stat dm
+  
+end
+
+macro test2
+
+  nnet create t1 2 2
+  nnet list
+
+  # output = sigmoid (sum(weight[i][j] * input[i]) + bias[j])
+
+  # biases of 
+  vlist v1 0.8 0.2
+
+  mcreate m1 2 2
+  m1[0][0] =  0.5
+  m1[1][0] =  0.1
+  m1[0][1] = -0.2
+  m1[1][1] =  0.8
+
+  # set weights and biases
+  nnet set t1 m1 v1
+
+  vlist vin0 1 2
+  vlist vin1 2 1
+
+  # get weights and biases
+  nnet apply t1 vin0 vin1 vout0 vout1
+  
+  vlist z0 {m1[0][0]*vin0[0] + m1[1][0]*vin1[0] + v1[0]} {m1[0][0]*vin0[1] + m1[1][0]*vin1[1] + v1[0]} 
+  vlist z1 {m1[0][1]*vin0[0] + m1[1][1]*vin1[0] + v1[1]} {m1[0][1]*vin0[1] + m1[1][1]*vin1[1] + v1[1]}
+  set t0 = 1 / (1 + exp(-1*z0))
+  set t1 = 1 / (1 + exp(-1*z1))
+
+  # compare
+  set dt0 = t0 - vout0
+  vstat dt0
+
+  set dt1 = t1 - vout1
+  vstat dt1
+end
+
+macro test3
+
+  memory all
+
+  nnet create t0 5 10 15 3
+  nnet list
+
+  nnet create t1 5 8
+  nnet list
+
+  # vlist v0 1 2 3 4 5
+  vlist v1 1 2 3 4 5 6 7 8
+
+  mcreate m1 5 8
+  
+  # generate a weight vector of length v1[]
+  vlist my 9 8 7 6 5 4 3 2
+
+  for i 0 5
+    mset m1 my -y $i
+    set my = my + $i
+  end
+  
+  # set weights and biases
+  nnet set t1 m1 v1
+
+  # get weights and biases
+  nnet get t1 M1 V1
+
+  # compare
+  set dv = v1 - V1
+  vstat dv
+
+  set dm = m1 - M1
+  stat dm
+  
+  nnet delete t0
+  nnet delete t1
+
+  delete m1 M1 dm
+  delete v1 V1 dv my
+
+  vectors
+  buffers
+
+  memory all
+end
+
+macro test.simple
+  if ($0 != 4)
+    echo "USAGE: test.simple (Nmini) (Nepoch) (eta)"
+    break
+  end
+
+  local Nmini Nepoch eta
+  $Nmini  = $1
+  $Nepoch = $2
+  $eta    = $3
+
+  # create a nnet with just 2 inputs and 2 outputs
+  nnet create t1 2 2
+
+  # biases:
+  vlist v1 0.5 0.0
+
+  # weights
+  mcreate m1 2 2
+  m1[0][0] =  0.5
+  m1[1][0] = -0.2
+  m1[0][1] = -0.5
+  m1[1][1] =  0.2
+
+  # set weights and biases
+  nnet set t1 m1 v1
+
+  # generate an input set spanning the range -5 to +5
+
+  create n 0 1000
+  set x0 = 10*rnd(n) - 5.0
+  set x1 = 10*rnd(n) - 5.0
+
+  # get output vectors the these input vectors
+  nnet apply t1 x0 x1 y0 y1
+
+# -large-weight-initializer : revert to original implementation
+# nnet create tn 2 2 
+  nnet create tn 2 2 -large-weight-initializer
+
+# -quadratic-cost : revert to original implementation
+  nnet train tn x0 x1 y0 y1 -Nmini $Nmini -Nepoch $Nepoch -eta $eta -lambda 0.0 -quadratic-cost -resid resid -result result
+# nnet train tn x0 x1 y0 y1 -Nmini $Nmini -Nepoch $Nepoch -eta $eta -lambda 0.0
+# nnet train tn x0 x1 y0 y1 -Nmini $Nmini -Nepoch $Nepoch -eta $eta -lambda 2.0
+
+  nnet apply tn x0 x1 Y0 Y1
+
+  set dy0 = y0 - Y0
+  set dy1 = y1 - Y1
+
+  vstat dy0
+  $yp = $MAX
+  $ym = $MIN
+
+  vstat dy1
+  $yp = max ($MAX , $yp)
+  $ym = min ($MIN , $ym)
+  $yp = 0.02; $ym = -0.02
+
+  style -pt circle -sz 2; lim x0 $ym $yp; clear; box; plot -c blue x0 dy0; plot -c red x0 dy1
+end
+
 macro test.bilevel
   if ($0 != 4)
@@ -284,2 +353,285 @@
 end
 
+macro test.bilevel.small
+  if ($0 != 4)
+    echo "USAGE: test.bilevel (Nmini) (Nepoch) (eta)"
+    break
+  end
+
+  local Nmini Nepoch eta
+  $Nmini  = $1
+  $Nepoch = $2
+  $eta    = $3
+
+  # create a nnet with just 2 inputs and 2 outputs
+  nnet create t1 2 3 2
+
+  # biases:
+  vlist v1 -0.2 0.0 0.2
+  vlist v2 -0.3 0.3
+
+  # weights
+  mcreate m1 2 3
+  m1[0][0] =  0.2
+  m1[1][0] = -0.2
+  m1[0][1] = -0.6
+  m1[1][1] =  0.5
+  m1[0][2] = -0.2
+  m1[1][2] =  0.2
+
+  mcreate m2 3 2
+  m2[0][0] = -0.6
+  m2[1][0] = -0.4
+  m2[2][0] = -0.5
+
+  m2[0][1] = -0.2
+  m2[1][1] = -0.5
+  m2[2][1] =  0.1
+
+  # set weights and biases
+  nnet set t1 m1 v1 m2 v2
+
+  # generate an input set spanning the range -5 to +5
+
+  create n 0 1000
+  set x0 = 10*rnd(n) - 5.0
+  set x1 = 10*rnd(n) - 5.0
+
+  # get output vectors the these input vectors
+  nnet apply t1 x0 x1 y0 y1
+
+# nnet create tn 3 4 3 -large-weight-initializer
+  nnet create tn 2 3 2
+
+# make the starting point close to the solution
+  m2[0][0] = -0.58
+# v2[1] = 0.25
+  nnet set tn m1 v1 m2 v2
+
+# nnet train tn x0 x1 x2 y0 y1 y2 -Nmini $Nmini -Nepoch $Nepoch -eta $eta -lambda 0.1
+# nnet train tn x0 x1 y0 y1 -Nmini $Nmini -Nepoch $Nepoch -eta $eta -lambda 0.0 -quadratic-cost -resid dS -result result
+  nnet train tn x0 x1 y0 y1 -Nmini $Nmini -Nepoch $Nepoch -eta $eta -lambda 0.0 -resid dS -result result
+
+  nnet apply tn x0 x1 Y0 Y1
+
+  set dy0 = y0 - Y0
+  set dy1 = y1 - Y1
+
+  vstat dy0
+  $yp = $MAX
+  $ym = $MIN
+
+  vstat dy1
+  $yp = max ($MAX , $yp)
+  $ym = min ($MIN , $ym)
+
+  dev -n 0 
+  style -pt circle -sz 2; 
+  lim x0 $ym $yp; clear; box; 
+  plot -c blue  x0 dy0; 
+  plot -c red   x0 dy1
+
+  set n = ramp(dS)
+  lim -n 1 n dS; clear; box; plot n dS
+end
+
+macro test.bilevel.grid.wt
+  if ($0 != 4)
+    echo "USAGE: test.bilevel.grid.wt (level) (x) (y)"
+    break
+  end
+
+  local Level ix iy
+  $Level  = $1
+  $ix     = $2
+  $iy     = $3
+
+  # create a nnet with just 2 inputs and 2 outputs
+  nnet create t1 2 3 2
+
+  # biases:
+  vlist v1 -0.2 0.0 0.2
+  vlist v2 -0.3 0.3
+
+  # weights
+  mcreate m1 2 3
+  m1[0][0] =  0.2
+  m1[1][0] = -0.2
+  m1[0][1] = -0.6
+  m1[1][1] =  0.5
+  m1[0][2] = -0.2
+  m1[1][2] =  0.2
+
+  mcreate m2 3 2
+  m2[0][0] = -0.6
+  m2[1][0] = -0.4
+  m2[2][0] = -0.5
+
+  m2[0][1] = -0.2
+  m2[1][1] = -0.5
+  m2[2][1] =  0.1
+
+  # set weights and biases
+  nnet set t1 m1 v1 m2 v2
+
+  # generate an input set spanning the range -5 to +5
+
+  create n 0 1000
+  set x0 = 10*rnd(n) - 5.0
+  set x1 = 10*rnd(n) - 5.0
+
+  # get output vectors from these input vectors
+  nnet apply t1 x0 x1 y0 y1
+
+  # create a test nnet
+  nnet create tn 2 3 2
+
+  # 1D chi-square grid about truth
+  
+  set M1 = m1
+  set M2 = m2
+  set V1 = v1
+  set V2 = v2
+
+  # disturb one element to see cross terms
+  M2[0][0] = -0.58
+
+  $To = M$Level[$ix][$iy]
+
+  delete -q value svec0 svec1 sigvec
+  for frac 0.90 1.10 0.005
+    M$Level[$ix][$iy] = $frac * $To
+
+    nnet set tn M1 V1 M2 V2
+    nnet apply tn x0 x1 Y0 Y1
+
+    set dy0 = y0 - Y0
+    set dy1 = y1 - Y1
+
+    vstat -q dy0
+    $S0 = $SIGMA  
+
+    vstat -q dy1
+    $S1 = $SIGMA
+
+    concat {$frac * $To} value
+    concat $S0 svec0
+    concat $S1 svec1
+
+    concat {sqrt($S0^2 + $S1^2)} sigvec
+  end
+
+  vstat -q sigvec
+
+  lim -n 1 value -0.0001 $MAX; clear; box; 
+  plot -c black value sigvec; 
+  plot -c blue value svec0 -pt ocir -sz 2 ; 
+  plot -c red value svec1 -pt ocir -sz 2
+end
+
+macro test.bilevel.grid.wt.2d
+  if ($0 != 7)
+    echo "USAGE: test.bilevel.grid.wt.2d (level) (x) (y) (level) (x) (y)"
+    break
+  end
+
+  local LevelA ixA iyA LevelB ixB iyB
+  $LevelA  = $1
+  $ixA     = $2
+  $iyA     = $3
+  $LevelB  = $4
+  $ixB     = $5
+  $iyB     = $6
+
+  # create a nnet with just 2 inputs and 2 outputs
+  nnet create t1 2 3 2
+
+  # biases:
+  vlist v1 -0.2 0.0 0.2
+  vlist v2 -0.3 0.3
+
+  # weights
+  mcreate m1 2 3
+  m1[0][0] =  0.2
+  m1[1][0] = -0.2
+  m1[0][1] = -0.6
+  m1[1][1] =  0.5
+  m1[0][2] = -0.2
+  m1[1][2] =  0.2
+
+  mcreate m2 3 2
+  m2[0][0] = -0.6
+  m2[1][0] = -0.4
+  m2[2][0] = -0.5
+
+  m2[0][1] = -0.2
+  m2[1][1] = -0.5
+  m2[2][1] =  0.1
+
+  # set weights and biases
+  nnet set t1 m1 v1 m2 v2
+
+  # generate an input set spanning the range -5 to +5
+
+  create n 0 1000
+  set x0 = 10*rnd(n) - 5.0
+  set x1 = 10*rnd(n) - 5.0
+
+  # get output vectors from these input vectors
+  nnet apply t1 x0 x1 y0 y1
+
+  # create a test nnet
+  nnet create tn 2 3 2
+
+  # 1D chi-square grid about truth
+  
+  set M1 = m1
+  set M2 = m2
+  set V1 = v1
+  set V2 = v2
+
+  # disturb one element to see cross terms
+  M2[0][0] = -0.58
+
+  $ToA = M$LevelA[$ixA][$iyA]
+  $ToB = M$LevelB[$ixB][$iyB]
+
+  delete -q valueA valueB sigvec
+
+  mcreate sigbuf 41 41
+
+  $ix = 0
+  for fracA 0.90 1.10 0.005
+    M$LevelA[$ixA][$iyA] = $fracA * $ToA
+
+    $iy = 0   
+    for fracB 0.90 1.10 0.005
+   
+      M$LevelB[$ixB][$iyB] = $fracB * $ToB
+
+      nnet set tn M1 V1 M2 V2
+      nnet apply tn x0 x1 Y0 Y1
+      
+      set dy0 = y0 - Y0
+      set dy1 = y1 - Y1
+      
+      vstat -q dy0
+      $S0 = $SIGMA  
+      
+      vstat -q dy1
+      $S1 = $SIGMA
+      
+      $sigval = sqrt($S0^2 + $S1^2)
+      concat {$fracA * $ToA} valueA
+      concat {$fracB * $ToB} valueB
+      concat $sigval sigvec
+
+      sigbuf[$ix][$iy] = $sigval
+      $iy ++
+    end
+    $ix ++
+  end
+
+  stat -q sigbuf
+  tv -n tv sigbuf $MIN {$MAX - $MIN}
+end