Index: trunk/Ohana/src/opihi/cmd.data/nnet.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/nnet.c	(revision 40317)
+++ trunk/Ohana/src/opihi/cmd.data/nnet.c	(revision 40317)
@@ -0,0 +1,52 @@
+# include "data.h"
+
+int nnet_list (int argc, char **argv);
+int nnet_create (int argc, char **argv);
+int nnet_delete (int argc, char **argv);
+int nnet_train (int argc, char **argv);
+int nnet_apply (int argc, char **argv);
+int nnet_show (int argc, char **argv);
+
+static Command nnet_commands[] = {
+  {1, "list",     nnet_list,     "list nnets"},
+  {1, "init",     nnet_init,     "initialize a nnet"},
+  {1, "delete",   nnet_delete,   "delete a nnet"},
+  {1, "show",     nnet_show,     "display nnet values"},
+  {1, "create",   nnet_create,   "create a nnet"},
+  {1, "train",    nnet_train,    "train nnet on a set of data"},
+  {1, "apply",    nnet_apply,    "apply nnet to a set of data"},
+};
+
+int nnet_command (int argc, char **argv) {
+
+  int i, N, status;
+
+  if (argc < 2) {
+    gprint (GP_ERR, "USAGE: nnet (command)\n");
+    gprint (GP_ERR, "    nnet list          : list nnets\n");
+    gprint (GP_ERR, "    nnet init   (nnet) : removes all pages from nnet\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 create (nnet) (Ninput) [Nnodes] [Nnodes] ... (Noutput) : create a nnet\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");
+    return (FALSE);
+  }
+
+  N = sizeof (nnet_commands) / sizeof (Command);
+
+  /* find the nnet sub-command which matches */
+  for (i = 0; i < N; i++) {
+    if (!strcmp (nnet_commands[i].name, argv[1])) {
+      status = (*nnet_commands[i].func) (argc - 1, argv + 1);
+      return (status);
+    }
+  }
+
+  gprint (GP_ERR, "unknown nnet command %s\n", argv[1]);
+  return (FALSE);
+}
+
+/* nnet is called with the command "nnet".  
+   the command line word "nnet" is meant to be followed the one of several 
+   possible options list above */
Index: trunk/Ohana/src/opihi/cmd.data/nnet_apply.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/nnet_apply.c	(revision 40317)
+++ trunk/Ohana/src/opihi/cmd.data/nnet_apply.c	(revision 40317)
@@ -0,0 +1,84 @@
+# include "data.h"
+
+int nnet_apply (int argc, char **argv) {
+
+  if (argc < 4) {
+    gprint (GP_ERR, "USAGE: nnet apply (nnet) [input] [input] ... [output] [output] ...\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;
+  }
+
+  // save the number of input and output nodes; these need to match the supplied vectors
+  int Nlayer   = nnet[0].Nlayer;
+  int Ninput   = nnet[0].Nnodes[0];
+  int Noutput  = nnet[0].Nnodes[Nlayer - 1];  // number of input + output nodes
+
+  // we need to have the right number of input and output vectors
+  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);
+    return FALSE;
+  }
+
+  // grab the input and output vectors from the argument list
+  ALLOCATE_PTR (inVec, Vector *, Ninput);
+  ALLOCATE_PTR (outVec, Vector *, Noutput);
+
+  // check that the vectors exist and that lengths match
+  int Ntrial = 0;
+  for (int i = 0; i < Ninput; i++) {
+    if ((inVec[i] = SelectVector (argv[i + 2], OLDVECTOR, FALSE)) == NULL) {
+      gprint (GP_ERR, "unknown input vector %s\n", argv[i+2]);
+      gprint (GP_ERR, "USAGE: nnet train (nnet) [input] [input] ... [output] [output] ...\n");
+      free (inVec);
+      free (outVec);
+      return (FALSE);    
+    }
+    if (Ntrial && (inVec[i][0].Nelements != Ntrial)) {
+      gprint (GP_ERR, "input vectors have inconsistent lengths: %d vs %d for %s\n", Ntrial, inVec[i][0].Nelements, inVec[i][0].name);
+      free (inVec);
+      free (outVec);
+      return (FALSE);    
+    }
+    if (!Ntrial) Ntrial = inVec[i][0].Nelements;
+    if (!Ntrial) {
+      gprint (GP_ERR, "trial vectors must be non-zero length\n");
+      free (inVec);
+      free (outVec);
+      return (FALSE);    
+    }
+  }    
+  for (int i = 0; i < Noutput; i++) {
+    if ((outVec[i] = SelectVector (argv[i + 2 + Ninput], ANYVECTOR, FALSE)) == NULL) {
+      gprint (GP_ERR, "invalid output vector %s\n", argv[i+2+Ninput]);
+      gprint (GP_ERR, "USAGE: nnet train (nnet) [input] [input] ... [output] [output] ...\n");
+      free (inVec);
+      free (outVec);
+      return (FALSE);    
+    }
+    ResetVector (outVec[i], OPIHI_FLT, Ntrial);
+  }    
+
+  for (int i = 0; i < Ntrial; i++) {
+
+    // store the input values for this row (trial) in the input value vector
+    for (int j = 0; j < Ninput; j++) {
+      nnet[0].svalue[0][j] = inVec[j][0].elements.Flt[i];
+    }
+    
+    nnet_feedforward (nnet);
+
+    // store the output value vector values in output vectors for this row
+    for (int j = 0; j < Noutput; j++) {
+      outVec[j][0].elements.Flt[i] = nnet[0].svalue[Nlayer-1][j];
+    }
+  }    
+
+  free (inVec);
+  free (outVec);
+  return TRUE;
+}
Index: trunk/Ohana/src/opihi/cmd.data/nnet_commands.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/nnet_commands.c	(revision 40317)
+++ trunk/Ohana/src/opihi/cmd.data/nnet_commands.c	(revision 40317)
@@ -0,0 +1,98 @@
+# include "data.h"
+
+int nnet_list (int argc, char **argv) {
+  OHANA_UNUSED_PARAM(argv);
+  if (argc != 1) {
+    gprint (GP_ERR, "USAGE: nnet list\n");
+    return FALSE;
+  }
+
+  ListNnets();
+  return TRUE;
+}
+
+/** Fix These ...
+int nnet_init (int argc, char **argv) {
+
+  int status;
+  Nnet *nnet;
+
+  if (argc != 2) {
+    gprint (GP_ERR, "USAGE: nnet init (nnet)\n");
+    return FALSE;
+  }
+
+  nnet = FindNnet (argv[1]);
+  if (nnet != NULL) {
+      status = DeleteNnet (nnet);
+      if (!status) abort ();
+  }
+
+  CreateNnet (argv[1]);
+  return TRUE;
+}
+
+int nnet_delete (int argc, char **argv) {
+
+  int status;
+  Nnet *nnet;
+
+  if (argc != 2) {
+    gprint (GP_ERR, "USAGE: nnet delete (nnet)\n");
+    return FALSE;
+  }
+
+  nnet = FindNnet (argv[1]);
+  if (nnet == NULL) {
+    gprint (GP_ERR, "nnet %s not found\n", argv[1]);
+    return FALSE;
+  }
+
+  status = DeleteNnet (nnet);
+  if (!status) abort ();
+  return TRUE;
+}
+
+int nnet_show (int argc, char **argv) {
+
+  Nnet *nnet;
+
+  if (argc != 2) {
+    gprint (GP_ERR, "USAGE: nnet listnnet (nnet)\n");
+    return FALSE;
+  }
+
+  nnet = FindNnet (argv[1]);
+  if (nnet == NULL) {
+    gprint (GP_ERR, "nnet %s not found\n", argv[1]);
+    return FALSE;
+  }
+
+  ListPages (nnet);
+  return TRUE;
+}
+**/
+
+int nnet_create (int argc, char **argv) {
+
+  if (argc < 5) {
+    gprint (GP_ERR, "USAGE: nnet create (nnet) (Ninput) [Nnodes] [Nnodes] ... (Noutput)\n");
+    return FALSE;
+  }
+
+  int Nlayer = argc - 2;
+
+  // nnet is guaranteed to be initialized with only name and the containers (but not vectors) allocated
+  Nnet *nnet = CreateNnet (argv[1], Nlayer);
+
+  // define the size of each layer
+  for (int i = 0; i < Nlayer; i++) {
+    nnet[0].Nnodes[i] = atoi (argv[i + 2]); // input (0), hidden layer nodes, output (Nlayer)
+  }
+
+  // allocate the vector and matrix data arrays (and init)
+  CreateNnetData (nnet);
+
+  return TRUE;
+}
+
Index: trunk/Ohana/src/opihi/cmd.data/nnet_train.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/nnet_train.c	(revision 40317)
+++ trunk/Ohana/src/opihi/cmd.data/nnet_train.c	(revision 40317)
@@ -0,0 +1,293 @@
+# include "data.h"
+void sortseq (float *X, int *Y, int N);
+
+int nnet_train (int argc, char **argv) {
+
+  int Nepoch = 10;
+  if ((N = get_argument (argc, argv, "-Nepoch"))) {
+    remove_argument (N, &argc, argv);
+    Nepoch = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  
+  int Nmini = 100;
+  if ((N = get_argument (argc, argv, "-Nmini"))) {
+    remove_argument (N, &argc, argv);
+    Nmini = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  
+  if (argc < 4) {
+    gprint (GP_ERR, "USAGE: nnet train (nnet) [input] [input] ... [output] [output] ...\n");
+    gprint (GP_ERR, "OPTIONS: -Nepoch [N] -Nmini [N]\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;
+  }
+
+  // save the number of input and output nodes; these need to match the supplied vectors
+  int Nlayer   = nnet[0].Nlayer;
+  int Ninput   = nnet[0].Nnodes[0];
+  int Noutput  = nnet[0].Nnodes[Nlayer - 1];  // number of input + output nodes
+
+  // we need to have the right number of input and output vectors
+  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);
+    return FALSE;
+  }
+
+  // grab the input and output vectors from the argument list
+  ALLOCATE_PTR (inVec, Vector *, Ninput);
+  ALLOCATE_PTR (outVec, Vector *, Noutput);
+
+  // check that the vectors exist and that lengths match
+  int Ntrial = 0;
+  for (int i = 0; i < Ninput; i++) {
+    if ((inVec[i] = SelectVector (argv[i + 2], OLDVECTOR, FALSE)) == NULL) {
+      gprint (GP_ERR, "unknown input vector %s\n", argv[i+2]);
+      gprint (GP_ERR, "USAGE: nnet train (nnet) [input] [input] ... [output] [output] ...\n");
+      free (inVec);
+      free (outVec);
+      return (FALSE);    
+    }
+    if (Ntrial && (inVec[i][0].Nelements != Ntrial)) {
+      gprint (GP_ERR, "input vectors have inconsistent lengths: %d vs %d for %s\n", Ntrial, inVec[i][0].Nelements, inVec[i][0].name);
+      free (inVec);
+      free (outVec);
+      return (FALSE);    
+    }
+    if (!Ntrial) Ntrial = inVec[i][0].Nelements;
+    if (!Ntrial) {
+      gprint (GP_ERR, "trial vectors must be non-zero length\n");
+      free (inVec);
+      free (outVec);
+      return (FALSE);    
+    }
+  }    
+  for (int i = 0; i < Noutput; i++) {
+    if ((outVec[i] = SelectVector (argv[i + 2 + Ninput], OLDVECTOR, FALSE)) == NULL) {
+      gprint (GP_ERR, "unknown output vector %s\n", argv[i+2+Ninput]);
+      gprint (GP_ERR, "USAGE: nnet train (nnet) [input] [input] ... [output] [output] ...\n");
+      free (inVec);
+      free (outVec);
+      return (FALSE);    
+    }
+    if (outVec[i][0].Nelements != Ntrial) {
+      gprint (GP_ERR, "output and input vectors have inconsistent lengths: %d vs %d for %s\n", Ntrial, outVec[i][0].Nelements, outVec[i][0].name);
+      free (inVec);
+      free (outVec);
+      return (FALSE);    
+    }
+  }    
+
+  // core of the training process:
+
+  // index to hold a random sequence
+  ALLOCATE_PTR (seq, int,   Ntrial);
+  ALLOCATE_PTR (rnd, float, Ntrial);
+
+  // train for Nepochs
+  // this recreates 'SGD' from http://neuralnetworksanddeeplearning.com/chap1.html
+  for (int epoch = 0; epoch < Nepoch; epoch ++) {
+    
+    // generate a random sequence : used to select random mini batches
+    for (int i = 0; i < Ntrial; i++) { seq[i] = i; rnd[i] = drand48(); }
+    sortseq (rnd, seq, Ntrial);
+
+    int Npass = Ntrial / Nmini;
+    for (int pass = 0; pass < Npass; pass ++) {
+      // for a given pass select seq elements pass*Nmini to pass*Nmini + Nmini - 1
+      // update the weights and biases using the mini batch subset
+      nnet_descent_step (nnet, inVec, outVec, seq, pass, Nmini, eta);
+      gprint (GP_ERR, "epoch %d of %d, pass %d of %d\n", epoch, Nepoch, pass, Npass);
+    }
+  }  
+
+  free (seq);
+  free (rnd);
+  free (inVec);
+  free (outVec);
+  return TRUE;
+}
+
+// this recreates 'update_mini_batch' from http://neuralnetworksanddeeplearning.com/chap1.html
+void nnet_descent_step (NNet *nnet, Vector **inVec, Vector **outVec, int *seq, int pass, int Nmini, float eta) {
+
+  int Ntrial = inVec[0][0].Nelements;
+
+  nnet_reset_Nabla (nnet);
+
+  for (i = 0; (i < Nmini) && (pass*Nmini + i < Ntrial); i++) {
+
+    // N is the element of the mini batch on which we are currently operating
+    N = seq[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);
+  }
+
+  nnet_apply_Nabla (nnet, Nmini, eta);
+}
+
+void nnet_backprop (NNet *nnet, Vector **inVec, Vector **outVec, int N) {
+
+  // start with the input values
+
+  // store the input values for this row (trial) in the input value vector
+  int Nlayer   = nnet[0].Nlayer;
+  int Ninput   = nnet[0].Nnodes[0];
+  int Noutput  = nnet[0].Nnodes[Nlayer - 1];  // number of input + output nodes
+
+  for (int j = 0; j < Ninput; j++) {
+    nnet[0].svalue[0][j] = inVec[j][0].elements.Flt[N];
+  }
+
+  // z = w * activation + bias 
+  // save z [one per non-input layer]
+  // activation = sigmoid (z)
+
+  nnet_feedforward (nnet);
+
+  // backward pass (note for now these are vector operations:
+  for (L = Nlayer - 1; L > 0; L--) {
+
+    // sp is array of same size as zvalue for each layer
+    // sp = sigmoid_prime (nnet[0].zvalue[L]); sprime is precalculated in feedforward
+
+    if (L == Nlayer - 1) {
+      // starting point uses cost_derivative to compare last svalue set with truth output
+      // delta = cost_derivative(nnet[0].svalue[L], outVec, N) * sprime[L];
+      for (int j = 0; j < nnet[0].Nnodes[L]; j++) {
+	nnet[0].delta[L][j] = cost_derivative(nnet[0].svalue[L][j], outVec[j][0].elements.Flt[N]) * nnet[0].sprime[L][j];
+      }
+    } else {
+      // delta = (delta DOT transpose(weights[L+1])) * sp;
+      for (int j = 0; j < nnet[0].Nnodes[L]; j++) {
+	tmpdelta[L][j] = 0;
+	for (int i = 0; i < nnet[0].Nnodes[L+1]; i++) {
+	  int k = j + i*nnet[0].Nnodes[L];
+	  tmpdelta[L][j] += nnet[0].weight[L][k] * nnet[0].sprime[L][i]; // XXX check on the index values
+	}
+	tmpdelta[L][j] * nnet[0].sprime[L][j];
+      }
+    }					       
+
+    // UPDATE
+    Nabla_b[L] = delta;
+    
+    // Nabla_w[L] = delta DOT transpose(nnet[0].svalue[L-1]);
+    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[n-1];
+	nnet[0]. Nabla_w[n][k] = nnet[0].svalue[L-1][i] * nnet[0].delta[L][j];
+      }
+    }
+
+  }
+}
+
+// support functions to loop over the Nabla entries
+void nnet_apply_Nabla (NNet *nnet, int Nmini, float eta) {
+  for (int n = 1; n < nnet[0].Nlayer; n++) {
+
+    for (int j = 0; j < nnet[0].Nnodes[n]; j++) {
+      nnet[0].biases[n][j] -= (eta / Nmini) * nnet[0].Nabla_b[n][j];
+
+      for (int i = 0; i < nnet[0].Nnodes[n-1]; i++) {
+	int k = i + j*nnet[0].Nnodes[n-1];
+	nnet[0].weights[n][k] -= (eta / Nmini) * nnet[0].Nabla_w[n][k];
+      }
+    }
+  }
+}
+void nnet_update_Nabla (NNet *nnet) {
+  for (int n = 1; n < nnet[0].Nlayer; n++) {
+
+    for (int j = 0; j < nnet[0].Nnodes[n]; j++) {
+      nnet[0]. Nabla_b[n][j] += nnet[0].dNabla_b[n][j];
+
+      for (int i = 0; i < nnet[0].Nnodes[n-1]; i++) {
+	int k = i + j*nnet[0].Nnodes[n-1];
+	nnet[0]. Nabla_w[n][k] += nnet[0].dNabla_w[n][k];
+      }
+    }
+  }
+}
+void nnet_reset_Nabla (NNet *nnet) {
+  for (int n = 1; n < nnet[0].Nlayer; n++) {
+
+    for (int j = 0; j < nnet[0].Nnodes[n]; j++) {
+
+      nnet[0]. Nabla_b[n][j] = 0;
+      nnet[0].dNabla_b[n][j] = 0;
+
+      for (int i = 0; i < nnet[0].Nnodes[n-1]; i++) {
+	int k = i + j*nnet[0].Nnodes[n-1];
+	nnet[0]. Nabla_w[n][k] = 0;
+	nnet[0].dNabla_w[n][k] = 0;
+      }
+    }
+  }
+}
+
+// the input values must already be copied to the input layer svalue[]
+void nnet_feedforward (NNet *nnet) {
+
+  for (int i = 1; i < nnet[0].Nlayer; i++) {
+    nnet_onelayer (nnet, i);
+  }
+  return;
+}
+
+// calcularte z, sigmoid(z) for each layer (z = w*value + bias)
+int nnet_onelayer (NNet *nnet, int n) {
+
+  if (n < 1) return FALSE; // abort here?
+  if (n >= nnet[0].Nlayer) return FALSE; // abort here?
+
+  // evaluating a single layer [n], n > 0, n < Nlayer:
+  int Ninput   = nnet[0].Nnodes[n - 1];
+  int Noutput  = nnet[0].Nnodes[n];
+
+  // input layer is [n-1], output layer is [n]
+  for (int j = 0; j < Noutput; j ++) {
+    float sum = 0;
+    for (int i = 0; i < Ninput; i++) {
+      // weight matrix order is (0, 1, ... Ninput-1, Ninput, Ninput + 1, ... Ninput * Noutput - 1)
+      int k = j * Ninput + i;
+      sum += nnet[0].weight[n][k]*nnet[0].svalue[n-1][i];
+    }
+    sum += nnet[0].biases[n][j];
+    nnet[0].zvalue[n][j] = sum;
+    nnet[0].svalue[n][j] = nnet_sigmoid(sum);
+    nnet[0].sprime[n][j] = nnet[0].svalue[n][j] * (1 - nnet[0].svalue[n][j]);
+    // note that d sigmoid / dz = sigmoid * (1 - sigmoid)
+  }
+  return TRUE;
+}
+
+float nnet_sigmoid (float value) {
+  return 1.0 / (1.0 + exp(-sum));
+}
+
+void sortseq (float *X, int *Y, int N) {
+
+# define SWAPFUNC(A,B){ float ftmp; int itmp; \
+    ftmp = X[A]; X[A] = X[B]; X[B] = ftmp;      \
+    itmp = Y[A]; Y[A] = Y[B]; Y[B] = itmp;       \
+}
+# define COMPARE(A,B)(X[A] < X[B])
+
+  OHANA_SORT (N, COMPARE, SWAPFUNC);
+
+# undef SWAPFUNC
+# undef COMPARE
+
+}
+
Index: trunk/Ohana/src/opihi/lib.data/nnet.c
===================================================================
--- trunk/Ohana/src/opihi/lib.data/nnet.c	(revision 40317)
+++ trunk/Ohana/src/opihi/lib.data/nnet.c	(revision 40317)
@@ -0,0 +1,224 @@
+# include "data.h"
+
+typedef struct {
+  char *name;
+  int Nlayers; // Nlayers = input layer + output layer + hidden layers
+  int *Nnodes; // number of nodes per layer
+  float **weights; // a matrix between each layer
+  float **biases; // a vector for each layer
+
+  float **zvalue; // a vector of z values for each layer (= w*input + b)
+  float **svalue; // a vector of s values for each layer (= sigmoid(z))
+  float **sprime; // 
+  float **delta;  // 
+
+  float ** Nabla_b; // a vector of Nabla_b values for each layer
+  float **dNabla_b; // a vector of Nabla_b values for each layer
+
+  float ** Nabla_w; // a matrix of Nabla_w values for each layer
+  float **dNabla_w; // a matrix of Nabla_w values for each layer
+} NNet;
+
+static Nnet **nnets  = NULL; /* nnet to store the list of all nnets */
+static int    Nnnets = 0;    /* number of currently defined nnets */
+static int    NNNETS = 0;    /* number of currently allocated nnets */
+
+void InitNnets () {
+  Nnnets = 0;
+  NNNETS = 16;
+  ALLOCATE (nnets, Nnet *, NNNETS); 
+}
+
+void FreeNnets () {
+
+  int i;
+
+  for (i = 0; i < Nnnets; i++) {
+    FreeNnet (nnets[i]);
+  }
+  free (nnets);
+}
+
+// InitNnetData requires nnet to unassigned or initialized
+// if you have an existing Nnet, call FreeNnetData first
+void InitNnetData (Nnet *nnet, char *name) {
+
+  nnet[0].name = strcreate (name);
+
+  nnet[0].Nlayers = 0;
+  nnet[0].Nnodes  = NULL;
+  nnet[0].weights = NULL;
+  nnet[0].biases  = NULL;
+  nnet[0].avalue  = NULL;
+  nnet[0].zvalue  = NULL;
+
+  nnet[0]. Nabla_b  = NULL;
+  nnet[0].dNabla_b  = NULL;
+  nnet[0]. Nabla_w  = NULL;
+  nnet[0].dNabla_w  = NULL;
+}
+
+void FreeNnetData (Nnet *nnet) {
+
+    int i;
+
+    free (nnet[0].name);
+    for (i = 0; i < nnet[0].Nlayers; i++) {
+      free (nnet[0].weights[i]);
+      free (nnet[0].biases[i]);
+      free (nnet[0].svalue[i]);
+      free (nnet[0].zvalue[i]);
+      free (nnet[0].sprime[i]);
+      free (nnet[0].delta [i]);
+
+      free (nnet[0]. Nabla_b[i]);
+      free (nnet[0].dNabla_b[i]);
+      free (nnet[0]. Nabla_w[i]);
+      free (nnet[0].dNabla_w[i]);
+    }
+    free (nnet[0].Nnodes);
+    free (nnet[0].weights);
+    free (nnet[0].biases);
+    free (nnet[0].svalue);
+    free (nnet[0].zvalue);
+    free (nnet[0].sprime);
+    free (nnet[0].delta );
+
+    free (nnet[0]. Nabla_b);
+    free (nnet[0].dNabla_b);
+    free (nnet[0]. Nabla_w);
+    free (nnet[0].dNabla_w);
+}
+
+/* return the given nnet */
+Nnet *GetNnet (int where) {
+
+  if (where < 0) where += Nnnets;
+  if (where < 0) return NULL;
+  if (where >= Nnnets) return NULL;
+  return (nnets[where]);
+}
+
+/* return the given nnet */
+Nnet *FindNnet (char *name) {
+
+  int i;
+
+  for (i = 0; i < Nnnets; i++) {
+    if (!strcmp (nnets[i][0].name, name)) {
+      return (nnets[i]);
+    }
+  }
+  return (NULL);
+}
+
+/* make a new named nnet with Nlayers */
+Nnet *CreateNnet (char *name, int Nlayer) {
+
+  int N;
+  Nnet *nnet;
+
+  nnet = FindNnet (name);
+  if (nnet != NULL) {
+    FreeNnetData (nnet);
+    InitNnetData (nnet);
+    return (nnet);
+  }
+
+  N = Nnnets;
+  Nnnets ++;
+  CHECK_REALLOCATE (nnets, Nnet *, NNNETS, Nnnets, 16);
+  ALLOCATE (nnet, Nnet, 1);
+  InitNnetData (nnet, name);
+
+  nnet[0].Nlayer = Nlayer;
+  ALLOCATE (nnet[0].Nnodes, int, Nlayer);
+  ALLOCATE (nnet[0].weights, float *, Nlayer);
+  ALLOCATE (nnet[0].biases, float *, Nlayer);
+  ALLOCATE (nnet[0].svalue, float *, Nlayer);
+  ALLOCATE (nnet[0].zvalue, float *, Nlayer);
+  ALLOCATE (nnet[0].sprime, float *, Nlayer);
+  ALLOCATE (nnet[0].delta , float *, Nlayer);
+
+  ALLOCATE (nnet[0]. Nabla_b, float *, Nlayer);
+  ALLOCATE (nnet[0].dNabla_b, float *, Nlayer);
+  ALLOCATE (nnet[0]. Nabla_w, float *, Nlayer);
+  ALLOCATE (nnet[0].dNabla_w, float *, Nlayer);
+
+  // save this nnet in the static array of nnets
+  nnets[N] = nnet;
+  return (nnet);
+}
+
+void CreateNnetData (Nnet *nnet) {
+
+  ohana_gaussdev_init ();
+
+  for (int i = 1; i < Nlayer; i++) {
+    ALLOCATE (nnet[0].biases[i], float, nnet[0].Nnodes[i]);  // biases for each node in the hidden and output layers only
+    for (int j = 0; j < nnet[0].Nnodes[i]; j++) {
+      nnet[0].biases[i][j] = ohana_gaussdev_rnd (0.0, 1.0);
+    }
+
+    ALLOCATE (nnet[0].weights[i], float, nnet[0].Nnodes[i-1]*nnet[0].Nnodes[i]);  // weights connected each node in the previous layer to the current layer (excludes input layer)
+    for (int j = 0; j < nnet[0].Nnodes[i-1]*nnet[0].Nnodes[i]; j++) {
+      nnet[0].weights[i][j] = ohana_gaussdev_rnd (0.0, 1.0);
+    }
+
+    ALLOCATE (nnet[0].svalue[i], float, nnet[0].Nnodes[i]);  // vectors for holding results / values for each node in the input, hidden, output layers
+    ALLOCATE (nnet[0].zvalue[i], float, nnet[0].Nnodes[i]);  // vectors for holding results / values for each node in the input, hidden, output layers
+    ALLOCATE (nnet[0].sprime[i], float, nnet[0].Nnodes[i]);  // vectors for holding results / values for each node in the input, hidden, output layers
+    ALLOCATE (nnet[0].delta [i], float, nnet[0].Nnodes[i]);  // vectors for holding results / values for each node in the input, hidden, output layers
+
+    ALLOCATE (nnet[0]. Nabla_b, float, nnet[0].Nnodes[i]);
+    ALLOCATE (nnet[0].dNabla_b, float, nnet[0].Nnodes[i]);
+
+    ALLOCATE (nnet[0]. Nabla_w, float, nnet[0].Nnodes[i-1]*nnet[0].Nnodes[i]);
+    ALLOCATE (nnet[0].dNabla_w, float, nnet[0].Nnodes[i-1]*nnet[0].Nnodes[i]);
+  }
+}
+
+/* delete a nnet */
+int DeleteNnet (Nnet *nnet) {
+
+  int i, N, NNNETS_2;
+
+  /* find nnet in nnet list */
+  N = -1;
+  for (i = 0; i < Nnnets; i++) {
+    if (nnets[i] == nnet) {
+      N = i;
+      break;
+    }
+  }
+  if (N == -1) return (FALSE);
+
+  for (i = N; i < Nnnets - 1; i++) {
+    nnets[i] = nnets[i + 1];
+  }
+  Nnnets --;
+  NNNETS_2 = MAX (16, NNNETS / 2);
+  if (Nnnets < NNNETS_2) {
+    NNNETS = NNNETS_2;
+    REALLOCATE (nnets, Nnet *, NNNETS);
+  }
+
+  FreeNnet (nnet);
+  return (TRUE);
+}
+
+/* list known nnets */
+void ListNnets () {
+
+  int i, j;
+
+  for (i = 0; i < Nnnets; i++) {
+    gprint (GP_ERR, "%-15s :", nnets[i][0].name);
+    for (j = 0; j < nnets[i][0].Nlayers; i++) {
+      gprint (GP_ERR, " %3d", nnets[i][0].Nnodes[j]);
+    }
+    gprint (GP_ERR, "\n");
+  }  
+  return;
+}
+