flag=0
proc err () { xpanel("**** ERROR ****") xlabel($s1) xpanel() }
strdef tstr
if (unix_mac_pc() == 3) {
   chdir(neuronhome())
   if (chdir("fctb")==0) flag=1
   if (chdir("..\fctb")==0) flag=1
   if (flag==0) {
     sprint(tstr,"ERROR: FCTB directory not found; should be in %s or %s/..\n",neuronhome(),neuronhome()) err(tstr)
   } else if (default_dll_loaded_== 0) {
     flag = nrn_load_dll("nrnmech.dll")
     if (flag==1) { default_dll_loaded_ = 1 } else {
       sprint(tstr,"ERROR loading %s/nrnmech.dll; remake with mknrndll.\n",getcwd()) err(tstr) }
   }
}

objref helpbox, helpdeck
proc help_list() {
  helpbox=new HBox(2)
  helpbox.intercept(1)
  xpanel("")
  xbutton("General","helpdeck.flip_to(0)")
  xbutton("Button actions","helpdeck.flip_to(1)")
  xbutton("Exercises","helpdeck.flip_to(2)")
  xpanel()

 helpdeck=new Deck()
helpdeck.intercept(1)
  xpanel("")
    xlabel("Model for Fig. 9-7. uses back-prop to calculate weights need to produce ")
    xlabel("VOR mapping (eyes/ears -> eye movement).  System will learn mapping")
    xlabel("between inputs and outputs by adjusting weights.")
    xlabel("")
  xpanel()
  xpanel("")
    xlabel("Stop processing: checkbox to interrupt \"learn\" loop")
    xlabel("Max iter: number of times to present patterns before stopping")
    xlabel("Sum err: average error across all patterns during learning")
    xlabel("Targ err: target for \"Sum err\" -- learning will stop when this")
    xlabel("  either this value is reached or iterations reach \"Max iter\"")
    xlabel("Forward prop: present input pattern (# given in value field) and")
    xlabel("  propagate through hidden to produce output compared to target")
    xlabel("Back prop: propagate error backwards and change weights")
    xlabel("Weights (pull-down menu):")
    xlabel("  Draw one: click on pre- then on postsynaptic unit to show weight")
    xlabel("    and then monitor in \"Weight monitor\" below ")
    xlabel("  Draw all: not very useful since puts lines everywhere; lines coming")
    xlabel("    from side are from bias units")
    xlabel("  Show in->hid matrix: pops up the input to hidden unit weight matrix")
    xlabel("    for editing; note that first row is not used")
    xlabel("  Show hid->out matrix: pops up the hidden to output unit weight")
    xlabel("    matrix for editing")
    xlabel("  Randomize weights: reset all weights to random values")
    xlabel("  Save/Read weight matrix: save/read file")
    xlabel("Patterns&Units")
    xlabel("  Clear all patterns")
    xlabel("  Load patterns: read from a file")
    xlabel("  Save patterns: save to a file")
    xlabel("  Show values: write the values on the units")
    xlabel("  Set inputs: Allows you to set the values for the input units")
    xlabel("    and feed them forward (use \"edit pattern\" to change the inputs")
    xlabel("    presented during learning).")
    xlabel("Edit pattern#: change one of the input/target pairs")
    xlabel("Weight monitor: allows you to monitor a single weight selected with")
    xlabel("  ")
    xlabel("")
  xpanel()
  xpanel("")
    xlabel("Working with large networks can be confusing.  You need to keep in")
    xlabel("mind the distinctions between input patterns and inputs and between")
    xlabel("outputs and targets.  It's also valuable to get a feel for changes")
    xlabel("in the weight matrices and their effects on the feedforward process.")
    xlabel("In addition, a complicated simulation requires a complicated GUI,")
    xlabel("which will take some learning to use.  Among other things, the")
    xlabel("exercises below try to lead you through button usage.")
    xlabel("")
    xlabel("As you go through the exercises that use the network to learn, you")
    xlabel("may wish to save the weight matrices (pull down \"Weights\"->\"Save")
    xlabel("weights\").  Use different file names (eg Exercise7.mat) and keep")
    xlabel("notes so that you know what is what later.")
    xlabel("")
    xlabel("1. Randomize the weights (pull down \"Weights\"->\"Randomize\").  Look")
    xlabel("at the Input->Hidden weight matrix (pull down \"Weights\"->\"Show")
    xlabel("in->hid matrix\").  Now draw the weights (pull down \"Weights\"->\"Draw")
    xlabel("all\").  Convince yourself that red and blue on the diagram")
    xlabel("correspond to positive and negative in the matrix.  Notice that the")
    xlabel("first row is all zero.  This is a confusing implementation detail:")
    xlabel("the hidden unit that the zeroth row projects to is a virtual unit")
    xlabel("that is used as the bias for the output units.  This virtual unit is")
    xlabel("always set to one and is not updated.")
    xlabel("")
    xlabel("2. Clear the graph (pull down \"Weights\"->\"Clear graph\").  Look")
    xlabel("through the patterns using \"Patterns&Units\"->\"Show patterns\".  The")
    xlabel("pattern number is given below next to \"Edit pattern #\".  After each")
    xlabel("of the 4 patterns is displayed, use \"Patterns&Units\"->\"Show values\"")
    xlabel("to write the numerical state values on top of the color-coded")
    xlabel("circles.  Put the labels back in (pull down \"Patterns&Units\"->\"Label")
    xlabel("layers\").  Notice that the top layer (targets) and bottom layer of")
    xlabel("inputs are the only ones updated.  As you go through the patterns")
    xlabel("notice that each pattern is also printed on the console (the hoc")
    xlabel("interpreter with the \"oc>\" prompt).")
    xlabel("")
    xlabel("3. Use \"Edit pattern #\" and change pattern #2's inputs to (1 0 0 0).")
    xlabel("Leave the outputs alone.  Unfortunately, the display does not update")
    xlabel("to show the change.  Cycle through all of the patterns")
    xlabel("(\"Patterns&Units\"->\"Show patterns\") to make sure that the pattern")
    xlabel("you saved was in fact saved.  Now go to \"Forward prop,\" set the value")
    xlabel("to '2' and hit the button.  Compare the row of the weight matrix")
    xlabel("(pull down \"Weights\"->\"Show in->hid matrix\") to the values in the")
    xlabel("units (\"Patterns&Units\"->\"Show values\").  Demonstrate that the")
    xlabel("values in the hidden units reflect the weights connecting the single")
    xlabel("set unit to each -- this is the divergence or projective field for")
    xlabel("this input unit.  Display individual weights on the graph")
    xlabel("(\"Weights\"->\"Draw one (click on units)\" -- click first on the active")
    xlabel("input unit and then click on one of the hidden units to show the")
    xlabel("weight in between).  Confirm that these are the same values as are")
    xlabel("shown in the weight matrix.  Notice that the weight you have chosen")
    xlabel("also shows up at the bottom of the panel next to \"Weight monitor\".")
    xlabel("")
    xlabel("4. Note the output values associated with the single set input unit")
    xlabel("from Exercise 3.  This is the projective field of that input unit at")
    xlabel("the output layer.  Demonstrate the projective field for each input")
    xlabel("unit by directly altering the units (use \"Patterns&Units\"->\"Set")
    xlabel("inputs\"; this is easier than editing an input pattern using \"Edit")
    xlabel("pattern #\" as done in the previous Exercise).  Trace the feedforward")
    xlabel("activity so as to relate the weight matrix patterns to the")
    xlabel("projective fields.  How does knowledge of the projective fields help")
    xlabel("you predict the receptive fields of the output units?  (Definition:")
    xlabel("the receptive field of a unit is the pattern of inputs that best")
    xlabel("activates that unit.)")
    xlabel("")
    xlabel("5. Enter the following patterns into the system: (1 0 0 0)->(1 0);")
    xlabel("(0 1 0 0)->(1 0); (0 0 1 0)->(0 1); (0 0 0 1)->(0 1).  You can do")
    xlabel("this first by using \"Edit pattern #\" for each.  For subsequent")
    xlabel("exercises however, it will be easier if you edit a 'ptn' file on")
    xlabel("your computer.  To do this, first save the patterns to a '.ptn' file")
    xlabel("(\"Patterns&Units\"->\"Save patterns\"), noting the path used to save")
    xlabel("them.  Then find the file using an editing program.  (If you use a")
    xlabel("word processor instead of an editor, make sure that you save in")
    xlabel("straight ascii without any additional formating.)  Implementation")
    xlabel("detail: the values for input and hidden units, Ni and Nh, as as")
    xlabel("given in the file, are increased by 1 to account for the bias units.")
    xlabel("Lines in the file list each input pattern (4 values) and output")
    xlabel("pattern (2 values) in order.  Edit this.  Clear all of the patterns")
    xlabel("(\"Patterns&Units\"->\"\"Clear all patterns\") and then read in your file")
    xlabel("(\"Patterns&Units\"->\"Load patterns\").  Use \"Patterns&Units\"->\"Show")
    xlabel("patterns\" to demonstrate that it has been loaded properly.")
    xlabel("Alternatively, use the command-line function \"showpatts()\" at the")
    xlabel("hoc prompt (\"oc>\").")
    xlabel("")
    xlabel("6. Use the patterns from Exercise 5.  Monitor one of the weights")
    xlabel("between hidden and output.  To monitor a weight, first draw it using")
    xlabel("\"Weights\"->\"Draw one (click on units)\" -- click first on a hidden")
    xlabel("unit and then click on one of the output units to show the weight in")
    xlabel("between.  The value of this weight will appear in the \"Weight")
    xlabel("monitor\" window and on the display.  Using the first pattern (#0) do")
    xlabel("single forward propagations (\"Forward prop\" button with 0 in the")
    xlabel("value box) followed by single back propagations (\"Back prop\"")
    xlabel("button).  Understand why the weight is being pushed in the direction")
    xlabel("that it is changing.  After this try monitoring a weight between")
    xlabel("input and hidden.  It may help to look at the intermediate")
    xlabel("calculation values being spit out on the console at the hoc prompt.")
    xlabel("")
    xlabel("7. Use the patterns from Exercise 5.  Learn these patterns (\"Learn\"")
    xlabel("button near top of panel).  After learning, document how the")
    xlabel("projective fields have changed for each input unit and how the")
    xlabel("receptive fields have changed for each output unit.  What are the")
    xlabel("projective and receptive fields like for each of the hidden units?")
    xlabel("This can be inferred from feedforward activity and then checked by")
    xlabel("observing the weights (\"Weights\"->\"Show ...\").  Write down the")
    xlabel("results of your hidden unit analysis.  Save the weights (\"Weights\"->\"Save")
    xlabel("weights\") so that you can recheck your findings later.")
    xlabel("")
    xlabel("8. Use the patterns from Exercise 5 (clear and read from file).")
    xlabel("Rerandomize the weights (\"Weights\"->\"Randomize weights\") and")
    xlabel("relearn.  Notice that the fields for inputs (projective field) and")
    xlabel("outputs (receptive field) are the same as before.  Look at the")
    xlabel("fields for the hidden units.  Are they the same as before?  Why or")
    xlabel("why not?  Restore the prior weights (\"Weights\"->\"Read weights\") to")
    xlabel("perform more comparisons.")
    xlabel("")
    xlabel("9. Go back to the default patterns (\"Reset\" button).  Randomize the")
    xlabel("weights (\"Weights\"->\"Randomize weights\").  Learn.  Think of a couple")
    xlabel("of examples where you would want the network should generalize to")
    xlabel("produce a sensible answer (eg what if you turn your head to the")
    xlabel("right too slowly while following a target to the right?).  Test the")
    xlabel("network (set the units and then feedforward or produce a file of")
    xlabel("testing patterns, read them in and test each in turn).  Does the")
    xlabel("network behave sensibly?  Can you find examples where the network")
    xlabel("does not behave sensibly?")
    xlabel("")
    xlabel("10. Randomize the weights.  Set \"Max iter\" to 10.  Learn.  How")
    xlabel("adequately has the system learned?  Does it generalize at all?  Push")
    xlabel("\"Learn\" again.  Each time you push it you will run another 10")
    xlabel("back-prop iterations and learn the patterns a little better.  See")
    xlabel("how well the system does after each learning set.  Look at the")
    xlabel("console output to see the numerical results.  Try again from scratch")
    xlabel("(rerandomize).  Is there any consistency to the pattern of learning?")
    xlabel("In this example one would not expect to find any pattern since all")
    xlabel("of the pattern inputs are equivalent in power.")
    xlabel("")
    xlabel("11. Give a behavioral-level description (eg head, target and eye")
    xlabel("movements) for the following patterns: (0 0 1 1)->(0.5 0.5); (1 1 0")
    xlabel("0)->0.5 0.5; (0 0.5 0.5 1)->(1 0); (1 0.5 0.5 0)->(0 1).  Enter")
    xlabel("these patterns.  Learn them.  Does the system learn adequately?")
    xlabel("Does it generalize?  Relearn from scratch 10 iterations at a time.")
    xlabel("Are some patterns learned before other patterns?  Why?")
    xlabel("")
    xlabel("12. Go back to the default patterns (\"Reset\" button).  Learn.")
    xlabel("Analyze the hidden unit receptive fields in terms of vestibular and")
    xlabel("target movement input.  Find a hidden unit in which these are")
    xlabel("paradoxical (eg responds to target to the right and also responds to")
    xlabel("head to the right).  Analyze the hidden unit projective fields in")
    xlabel("terms of eye movement.  Find a paradoxical hidden unit (eg, one that")
    xlabel("pushes strongly right and also pushes strongly left).  Are the")
    xlabel("receptive field paradoxical units and the projective field")
    xlabel("paradoxical units one and the same?")
  xpanel()
helpdeck.intercept(0)
  helpdeck.map()
  helpdeck.flip_to(0)
  helpbox.intercept(0)
  helpbox.map()
}
// Created 07/24/09 17:50:30 by "/usr/site/scripts/loadfiles -q -a SUP_VOR.hoc"
//================================================================
// INSERTED SUP_VOR.hoc
// =Id=  SUP_VOR.hoc,v 1.14 2003/01/27 19:10:50 billl Exp 

// load_file("SUP_VOR.hoc")
//================================================================
// INSERTED string2.hoc
// =Id=  string2.hoc,v 1.1 2002/09/07 12:55:06 billl Exp 

//*String2 template
begintemplate String2
public s,t
strdef s,t
proc init() {
    if (numarg() == 1) { s=$s1 }
    if (numarg() == 2) { s=$s1 t=$s2 }
  }
endtemplate String2
// END string2.hoc
//================================================================
//================================================================
// INSERTED bkutils.hoc
// =Id=  bkutils.hoc,v 1.73 2004/04/30 15:35:32 billl Exp 

// load_file("stdgui.hoc")
// load_file("setup.hoc")
load_file("stdrun.hoc")
if (name_declared("install_matrix")) execute("install_matrix()") 
if (name_declared("install_vecst")) execute("install_vecst()")

//* objects
objref  g[4], stim[2], nc, gp, sns, ind, vec0, vec[2], rdm[2], XO, YO, cvode
objref st,wt,scr,smat,box[4],insr,nil,sref,tmpfile,tmpobj,tmpvec,tmplist,sfunc,deck
strdef mesg,tstr,tstr2,temp_string_,filename,section
tmpfile = new File()
tmplist = new List()
fchooser_flag = 0
graph_flag = 1
double x[4],y[4]
rdm = new Random()
cvode = new CVode()
sfunc = new StringFunctions()
{ ind=new Vector() vec[0]=ind.c vec[1]=ind.c vec0=ind.c}

//* iterators and templates
//** string iterator
iterator case() { local i
    i1 = 0
    for i = 2, numarg() {
        $&1 = $i
        iterator_statement
        i1+=1
      }
  }

iterator scase() { local i
    i1 = 0
    for i = 1, numarg() {
        temp_string_ = $si
        iterator_statement
        i1+=1
      }
  }

// eg for scase2("a","b","c","d","e","f") print temp_string_,temp_string2_
iterator scase2() { local i
    i1 = 0
    if (numarg()%2==1) {print "ERROR: scase2 needs even number of args" return }
    for i = 1, numarg() {
        tmpobj=new String2()
        tmpobj.s=$si i+=1  tmpobj.t=$si
        iterator_statement
        i1+=1
      }
  }

// like perl chop -- removes the last character
proc chop () { sfunc.left($s1,sfunc.len($s1)-1) }

//** list iterator ltr
// usage 'for ltr(XO, tmplist) { print XO }'
iterator ltr() { local i
    if (numarg()==3) {$&3=0} else {i1 = 0}
    for i = 0, $o2.count() - 1 {
        $o1 = $o2.object(i)
        iterator_statement
        if (numarg()==3) { $&3+=1 } else { i1+=1 }
      }
    $o1 = nil
  }

//** list pairwise iterator ltrp
// usage 'for ltrp(XO, YO, list) { print XO,YO }' takes them pairwise
iterator ltrp() { local i
    if (numarg()==4) {$&4=0} else {i1 = 0}
    for (i=0;i<$o3.count()-1;i+=2) {
        $o1 = $o3.object(i) $o2 = $o3.object(i+1)
        iterator_statement
        if (numarg()==4) { $&4+=1 } else { i1+=1 }
      }
    $o1=nil $o2=nil
  }

//** vector iterator vtr
// usage 'for vtr(&x, vec) { print x }'
iterator vtr() { local i
    if (numarg()==3) {$&3=0} else {i1 = 0}
    for i = 0, $o2.size() - 1 {
        $&1 = $o2.x[i]
        iterator_statement
        if (numarg()==3) { $&3+=1 } else { i1+=1 }
      }
  }

//* vlk(vec) -- display vector slice -- more flexible than vec.printf
// vlk(vec,max)
// vlk(vec,min,max)
// prints out a segment of a vector
vlk_width=20
proc vlk () { local i,j,min,max,dual,wdh,nonl,nl
    j=dual=0 nl=1 wdh=vlk_width
    if (numarg()==1) { min=0 max=$o1.size-1 }
    if (numarg()==2) if ($2==0) { nl=min=0 max=$o1.size-1 // vlk(vec,0) flag to suppress new lines
          } else if ($2>0) { min=0 max=$2-1 } else { min=$o1.size+$2 max=$o1.size-1 }
    if (numarg()==3) if ($3>-1) { min=$2 max=$3 } else { min=0 max=$o1.size-1 dual=1 }
    if (numarg()==4) { min=$3 max=$4 dual=1 }
    if (min<0) min=0
    if (max>$o1.size-1) max=$o1.size-1
    if (dual) if (max>$o2.size-1) max=$o2.size-1
    for i=min,max { 
        if (dual) printf("%g:%g ",$o1.x[i],$o2.x[i]) else printf("%g ",$o1.x[i]) 
        if ((j=j+1)%vlk_width==0 && nl) { print "" }
      }
    if (nl) print ""
  }

//** savevec([list,]vec1[,vec2,...]) add vector onto veclist or other list if given as 1st arg
objref veclist
veclist = new List()
proc savevec () { local i,flag,beg
    if (isobj($o1,"List")) beg=2 else beg=1
    for i=beg, numarg() { 
        tmpvec = new Vector($oi.size)
        tmpvec.copy($oi)
        if (beg==2) $o1.append(tmpvec) else veclist.append(tmpvec)
        tmpvec = nil
      }
  }

//* graphics; symbols and colors
{symnum = 7 colnum = 9}
objectvar cl[colnum], sym[symnum]
for ii=0,colnum-1 cl[ii]=new String()
{ cl[0].s ="white" cl[1].s ="black"  cl[2].s ="red"  cl[3].s ="blue"  cl[4].s ="green"  
  cl[5].s ="orange" cl[6].s ="brown"  cl[7].s ="violet"  cl[8].s ="yellow"  }
for ii=0,symnum-1 sym[ii]=new String()
{ sym[0].s = "o" sym[1].s = "t" sym[2].s = "s"  
    sym[3].s = "O" sym[4].s = "T" sym[5].s = "S"  sym[6].s = "+"}

gcl=0  // graph color
func color () {
    gcl = (gcl+1)%colnum
    if (gcl==0) gcl=1 // throw out white
    if (numarg()==1) $o1.color(gcl)
    return gcl
  }

//** qhp() basic buttons
proc qhp () {
    if (numarg()==1) xpanel("CONTROLS",1)
    xbutton("Quit","quit()")
    xbutton("Help","help_list()")
    xbutton("Print panel","pwman_place(50,50)")  
    xbutton("Reset","reset()")
    if (numarg()==1) xvarlabel(mesg)
    if (numarg()==1) xpanel(0,0)
    mesg="" 
  }

//** allgraphs()
proc remgrs () { local ii
    for ii=0,3 graphList[ii].remove_all()
    objref graphItem
  }
proc allgraphs () { for ii=0,3 for ltr(XO,graphList[ii]) XO.exec_menu($s1) }
proc vp() { allgraphs("View = plot") }
proc allgrop () { 
    for ii=0,3 for ltr(XO,graphList[ii]) {
        sprint(tstr,"%s.%s",XO,$s1)
        execute(tstr)
      }
  }

//** grrtsize() use view=plot and then pad a little
proc grrtsize () { local h,w,frac
    if (numarg()>=1) tmpobj=$o1 else tmpobj=graphItem
    if (numarg()>=2) frac = $2 else frac=.05
    tmpobj.exec_menu("View = plot")
    tmpobj.size(&x)
    w=frac*(x[1]-x[0])  h=frac*(x[3]-x[2])
    x[0]-=2*w x[1]+=w x[2]-=4*h x[3]+=h // need extra padding on bottom
    tmpobj.size(x[0],x[1],x[2],x[3])
  }

//** isobj(o1,s2) checks whether object $o1 is of type $s2
func isobj () {
    sprint(temp_string_,"%s",$o1)
    if (sfunc.substr(temp_string_,$s2)==0) {
        return 1
      } else {
        return 0
      }
  }

// cbin(min,max): binary colormap
proc cbin () {
    $o1.colormap(2)
    $o1.colormap(0, 255, 0, 0)
    $o1.colormap(1, 255, 255, 0)
    $o1.scale($2, $3)
  }

// my version of newPlotV() leaves off the label
proc newplotv () {
    newplot()
    graphItem.addvar("","v(.5)",1,3)
  }

proc redraw0 () { plrdrf=0 } // stub with flag
proc plcback () { } // callback
objref oxv,oyv // coordinates: 0,1 line start; 2,3 last loc
oxv=new Vector(2)
oyv=new Vector(2)
plrdrf=1 // flag will be unset if using redraw stub
plhorz=0 // 1 for horizontal, 2 for verticle lines

proc pl () {
    if ($1==2) { // mouse down -- new line
        redraw0() // unset plrdrf if redraw0() not overwritten
        oxv.x[0]=oxv.x[1]=$2 oyv.x[0]=oyv.x[1]=$3
      } else {
        if (plrdrf) graphItem.erase_all
        oyv.line(graphItem,oxv,0,2)
        oxv.x[1]=$2 oyv.x[1]=$3
        if (plhorz==1) oyv.x[1]=oyv.x[0]
        if (plhorz==2) oxv.x[1]=oxv.x[0]
        oyv.line(graphItem,oxv,2,2)
        redraw0() // will redraw if defined
        if ($1==3) plcback(oxv,oyv)
      }
  }

// newplot([width,height,pointer]) OR newplot([pointer])
proc newplot () { local wd,ht
    graphItem = new Graph(0)
    if (numarg()==1) $o1=graphItem // make a pointer
    if (numarg()==3) $o3=graphItem
    if (numarg()>1) { wd=$1 ht=$2 } else { wd=500 ht=300 }
    graphItem.save_name("graphList[0].")
    graphList[0].append(graphItem)
    graphItem.view(0,-90,tstop,150,600,200,wd,ht)
  }

//* cbw(min,max): b/w binary colormap
proc cbw () {
    $o1.colormap(2)
    $o1.colormap(0, 0, 0, 0)
    $o1.colormap(1, 255, 255, 255)
    $o1.scale($2, $3)
  }

//* I/O
//** svrd(0/1,name,file_ext[,fflag]) save-0 and read-1
// eg svrd(1,"Image","img") will look for file.img files to read
// fflag -- return filename only -- don't call parser
proc svrd () { local ii,num,cnt,fflag
    rdflag=$1
    if (numarg()==4) fflag=$4 else fflag=0
    sprint(tstr2,"*.%s",$s3)
    if (rdflag==0) { // write
        sprint(tstr,"Write %s",$s2)
        tmpfile.chooser("w",tstr,tstr2,"WRITE","","")
        if (tmpfile.chooser()==1) {
            tmpfile.getname(filename)
            if (!fflag) parsr(rdflag)
          }
      } else if (rdflag==1) { // read
        sprint(tstr,"Read %s",$s2)
        tmpfile.chooser("r",tstr,tstr2,"READ","","")
        if (tmpfile.chooser()==1) {
            tmpfile.getname(filename)
            if (!fflag) parsr(rdflag)
          }
      }
    mesg=filename
    tmpfile.close
  }

//* xgetargs, d2b
//** range(val,min,max) -- return val only if in proper range
// eg stim.noise=range(stim.noise,0,1)
func range () {
    if ($1<$2) return $2 
    if ($1>$3) return $3
    return $1
  }

//** err()
proc err () {
    if (numarg()==1) mesg=$s1
    print mesg
  }

//** xvarstr(): display a list of strings in an xpanel using xvarlabel
proc xvarstr () { local ii,cnt
    // ivoc_style("*font", "fixed")
    if (unix_mac_pc() == 1) {
        ivoc_style("*font", "*helvetica-bold-r-normal*--14*")
      }
    xpanel($s1)
    for ltr(XO,$o2) xvarlabel(XO.s)
    xpanel()
  }

//** xgetargs(panel_name,command,arg1[,arg2,...],defaults)
// xgetargs("Random session","newrand","# of patts","patt size  ","overlap   ","5,33,7")
objref argv
argv = new Vector()
proc xgetargs () { local i,args
    args=numarg()-3  i=numarg()
    argv.resize(0)
    sprint(temp_string_,"argv.append(%s)",$si)
    execute(temp_string_)
    if (argv.size!=args) argv.resize(args)
    xpanel($s1)
    mesg=$s1
    xvarlabel(mesg)
    for i=3,numarg()-1 {
        sprint(temp_string_,"argv.x[%d]",i-3)
        xvalue($si,temp_string_)
      }
    sprint(temp_string_,"xgetexec(\"%s\",%d)",$s2,args)
    xbutton("Execute",temp_string_)
    xpanel()
  }

proc xgetexec () { local i,args
    args = $2
    if (argv.size!=args) { mesg="Error-close & relaunch panel" return }
    sprint(temp_string_,"%s(",$s1)
    for i=0,args-2 sprint(temp_string_,"%s%g,",temp_string_,argv.x[i])
    sprint(temp_string_,"%s%g)",temp_string_,argv.x[i])
    print temp_string_
    execute(temp_string_)
  }
  

//** d2b(): decimal to binary, gives a binary appearing number for showing
// eg for ii=0,100 printf("%010d\n",d2b(ii))
func d2b () { local num,ii,rem,res
    num=$1 ii=0 res=0
    while (num>0) {
        rem = int(2*(num/2-int(num/2))) // right digit
        num=int(num/2) // remove right digit
        res+=rem*10^ii
        ii+=1
      }
    return res
  }

//* matrix stuff

//** cvec(),rvec() put out col bzw row vecs in latex format
proc cvec () { local i
    printf("\\bpm ")
    for vtr(&x,$o1) if (i1<$o1.size-1) printf("%g \\\\ ",x)
    printf("%g \\epm\n",x)
  }

proc rvec () { local i
    printf("\\bpm ")
    for vtr(&x,$o1) if (i1<$o1.size-1) printf("%g & ",x)
    printf("%g \\epm\n",x)
  }

//** matp(mat,rows,cols) prints out matrix in latex format
proc matp () { local i,j,rows,cols
    rows=$2 cols=$3
    if ($o1.size!=rows*cols) {print "ERR: matp - mat.size != rows*cols" return}
    printf("\\bpm\n")
    for (i=0;i<rows;i+=1) {
        for (j=0;j<cols-1;j+=1) {
            printf("%4g & ",$o1.x[i*cols+j])
          }
        printf("%4g \\\\ \n",$o1.x[i*cols+j])
      }
    printf("\\epm\n")
  }

//** same as matp but not latex format
proc mpf () { local i,j,rows,cols
    rows=$2 cols=$3
    if ($o1.size!=rows*cols) {print "ERR: - mat.size != rows*cols" return}
    for (i=0;i<rows;i+=1) {
        for (j=0;j<cols-1;j+=1) {
            printf("%12.3g  ",$o1.x[i*cols+j])
          }
        printf("%12.04g \n",$o1.x[i*cols+j])
      }
  }

//** print out a matrix*vec multiplication
proc mpfv () { local i,j,rows,cols
    cols=$o2.size rows=$o1.size/cols
    for (i=0;i<rows;i+=1) {
        printf("|  ")
        for (j=0;j<cols;j+=1) {
            printf("%-12.03g  ",$o1.x[i*cols+j])
          }
        if (i<$o2.size) printf("  |   %12.04g \n",$o2.x[i]) else printf("  |\n")
      }
    if (cols>rows) for i=rows,cols-1 printf("%70.04g \n",$o2.x[i])
  }

//** mmlt(mdest,m1,m2,r2) matrix mult between 2 mats
// r2==c1 gives shared dimension
proc mmlt () { local i,j,r1,c1,r2,c2
    c1 = r2 = $4 // shared dimension (cols of 1st/rows of second)
    r1=$o2.size/c1 c2=$o3.size/r2
    printf("%dx%d * %dx%d = %dx%d\n",r1,c1,r2,c2,r1,c2)
    $o1.resize(r1*c2)
    vec[0].resize(c1) // length of a row of mat A
    vec[1].resize(r2) // length of a col of mat B
    for i=0,r1-1 for j=0,c2-1 {
        vec[0].mrow($o2,i,c1) vec[1].mcol($o3,j,c2)
        $o1.x[i*c2+j]=vec[0].dot(vec[1])
      }
    $o1.mprintf(r1,c2)
    dealloc(a)
  }

//* membrane and pp's
//** runbutton
running_ = 1
proc runbutton () {
    if (running_ == 0) { stoprun = 1  return }
    if (t>0 && t<tstop-0.1) continuerun(tstop) else run()
    stoprun=1 running_=0
  }

//** proc print_pp_location(PP), from doc/html/refman/nocmodl.html
proc print_pp_location() { local x //arg1 must be a point process
     x = $o1.get_loc()
     sectionname(section)
     printf("%s located at %s(%g)\n", $o1, section, x)
     pop_section()
  }

func get_pp () { local x
     x = $o1.get_loc()
     sectionname(section)
     section sref = new SectionRef()
     pop_section()
     return x
  }
// END bkutils.hoc
//================================================================
//================================================================
// INSERTED ncirc.hoc
// =Id=  ncirc.hoc,v 1.36 2002/10/02 13:41:02 billl Exp 

// load_file("ncirc.hoc")
// eg
// objref nc
// nc = new NCIRC(10,10,90,10,10,90,90,90)
// vec.resize(4)
// for ii=0,20 { vec.fill(1) vec.x[ii%4]=0 nc.disp(vec) system("sleep .3") nc.g.flush doEvents()}
// nc.gcmap(-1,1)
// rdm.uniform(-1,1)
// for ii=0,20 { vec.setrand(rdm) vec.printf nc.dispcol(vec) system("sleep .2") nc.g.flush doEvents()}
ncview=1
begintemplate NCIRC
public vals,g,disp,sz,coord,clrsym,symsz,setsym,clrcol,setcol,colr,labl,line,clrall,clkset
public low,high,gcmap,dispcol,lscale,which,mapg
external ncview

objref vals[3], xloc[3], yloc[3], g, mapg
strdef clrsym,setsym,symb,tstr
double low[1],high[1]

proc init () { local i
    g=new Graph(0)
    g.menu_tool("Set units", "clkset")
    g.size(0,100,0,100)
    g.xaxis(3)
    for ii=0,2 { vals[ii]=new Vector(0) xloc[ii]=vals.c yloc[ii]=vals.c }
  
    clrcol=1 setcol=2 colr=1 symsz=40 lscale=2
    clrsym="o" setsym="O"
    sz = numarg() i=1
    vals.resize(sz/2) xloc.resize(0) yloc.resize(0)
    while (i<sz) { 
        xloc.append($i)
        i+=1
        yloc.append($i)
        i+=1
      }
    sz /= 2
    if (ncview) g.view(0,0,100,100,600,200,400,300) 
  }

proc coord () { for ii=0,sz-1 print xloc.x[ii],yloc.x[ii] }

proc disp () {
    if (numarg()==1) {
        if (vals.size!=$o1.size) { printf("NCIRC ERR: Vec arg must be size %d.\n",vals.size) return }
        vals.copy($o1,0)
      }
    vals[1].indvwhere(vals,"==",BVBASE)
    vals[2].indvwhere(vals,"!=",BVBASE)
    xloc[1].index(xloc,vals[1])  // clear bits
    yloc[1].index(yloc,vals[1])
    xloc[2].index(xloc,vals[2])  // set bits
    yloc[2].index(yloc,vals[2])
    g.erase_all()
    yloc[1].mark(g,xloc[1],clrsym,symsz,clrcol,3)
    yloc[2].mark(g,xloc[2],setsym,symsz,setcol,3)
  }

// dispcol() for display of analog vlaues
// uses low and high and divides between colors
// 0 (white) and 1 (black) are reserved
proc cmap () {
    $o1.color(10,"#0000cf")  $o1.color(11,"#0000ff")  $o1.color(12,"#002fff")  $o1.color(13,"#005fff")
    $o1.color(14,"#008fff")  $o1.color(15,"#00bfff")  $o1.color(16,"#00efff")  $o1.color(17,"#1fffef")
    $o1.color(18,"#4fffbf")  $o1.color(19,"#7fff8f")  $o1.color(20,"#afff5f")  $o1.color(21,"#dfff2f")
    $o1.color(22,"#ffff00")  $o1.color(23,"#ffcf00")  $o1.color(24,"#ff9f00")  $o1.color(25,"#ff6f00")
    $o1.color(26,"#ff3f00")  $o1.color(27,"#ff0f00")  $o1.color(28,"#df0000")
  }

proc gcmap () { local step,a,ii
    maxcol=19
    low=$1 high=$2 step=($2-$1)/(maxcol-1)
    mapg = new Graph(0)
    cmap(mapg)
    for ii=10,28 {
        mapg.mark(1,low+(ii-10)*step,"S",20,ii,2)
        mapg.mark(1,low+(ii-10)*step,"s",20,1,2)
      } 
    mapg.view(-1,low,3,high-low,600,200,10,200)
  }

proc dispcol () {
    if (numarg()==1) {
        if (vals.size!=$o1.size) { printf("NCIRC ERR2: Vec arg must be size %d.\n",vals.size) return }
        vals.copy($o1,0)
      }
    g.erase_all()
    for ii=0,vals.size-1 {
        g.mark(xloc.x[ii],yloc.x[ii],"O",symsz,pkcol(vals.x[ii]),1)
        g.mark(xloc.x[ii],yloc.x[ii],"o",symsz,1,3) // black border
      }
  }

// pkcol() -- pick the right color
func pkcol () { local y
    y=$1
    if (y<=low) return 10  // 10 is min since don't want to overwrite built-in colors
    if (y>=high) return maxcol+10-1
    return int(maxcol*(y-low)/(high-low))+10
  }

// put labels in each circle
proc labl () { local i,j
    if (numarg()==0) { // label all units with values in val
        for ii=0,sz-1 {
            sprint(tstr,"%0.1f",vals.x[ii])
            g.label(xloc.x[ii],yloc.x[ii],tstr,1,lscale,.5,.5,colr)      
          }
      } else if (numarg()==2 && xloc.size!=2) { // label one unit
        g.label(xloc.x[$1],yloc.x[$1],$s2,1,lscale,.5,.5,colr)
      } else {
        for i=1,numarg() {  // label all the units with consecutive arg strings
            j=i-1
            g.label(xloc.x[j],yloc.x[j],$si,1,lscale,.5,.5,colr)
          }
      }
  }

// draw a lines btwn units with a spot on the end
// line(from,to[,color,symbol,symbol_size,x_offset,y_offset)]
func min () { if ($1<$2) return $1 else return $2 }
proc line () { local i,b,e,ssz,offx,offy
    b=$1 e=$2
    if (numarg()>2) colr=$3 // else as set
    if (numarg()>3) symb=$s4 else symb="O"
    if (numarg()>4) ssz=$5 else ssz=15
    if (numarg()>5) {offx=$6 offy=$7} else {offx=0 offy=0}
    if (numarg()>7) wth=$8 else wth=0
    g.color(colr)
    g.beginline()
    g.line(xloc.x[b]+offx,yloc.x[b]+offy)
    g.line(xloc.x[e]+offx,yloc.x[e]+offy)
    g.color(1)
    g.mark(xloc.x[e]+offx,yloc.x[e]+offy, symb , ssz, colr, 1)
    if (wth) { // weight label
        xa=min(xloc.x[b],xloc.x[e])
        ya=min(yloc.x[b],yloc.x[e])
        sprint(tstr,"%1.1f",wth)
        g.label(xa+offx+abs(xloc.x[e]-xloc.x[b])/2,ya+offy+abs(yloc.x[e]-yloc.x[b])/2,tstr,1,1,.5,.5,colr)
      }
  }

//** clkset() binary change unit state which click 
proc clkset () { local min,dist,whc,x,y
    if ($1==2) { // mouse down only  
        min=1e10 x=$2 y=$3
        for ii=0,sz-1 {
            dist=abs((xloc.x[ii]-x)*(xloc.x[ii]-x) + (yloc.x[ii]-y)*(yloc.x[ii]-y))
            if (dist<min) {min=dist whc=ii}
          }
        if (vals.x[whc]==BVBASE) vals.x[whc]=1 else vals.x[whc]=BVBASE
        disp()
      }
  }

//** which() return id of a unit by coordinates
func which () { local min,dist,whc,x,y
    x=$1 y=$2 min=1e10
    for ii=0,sz-1 {
        dist=abs((xloc.x[ii]-x)*(xloc.x[ii]-x) + (yloc.x[ii]-y)*(yloc.x[ii]-y))
        if (dist<min) {min=dist whc=ii}
      }
    return whc
  }

proc clrall () { vals.fill(BVBASE) }

endtemplate NCIRC
// END ncirc.hoc
//================================================================
//================================================================
// INSERTED xtabl.hoc
// =Id=  xtabl.hoc,v 1.29 2002/09/20 14:49:46 billl Exp 

// load_file("xtabl.hoc")
// new XTABL(rows,cols,vec[,key,"NAME",arrows])
// vec.size must == rows*cols
// key identifies which XTABL in case of callbacks
// arrows = 1 puts in checkbox and scroll arrows

// notify(vec,index,row,col,key) -- new value is vec.x[index]
// index is redudant since = row*cols+col
proc notify () {} // stub for callbacks
// template XTABL(rows,cols)
begintemplate XTABL
public rows,cols,key,label,mapdv,arrow
external notify,XTABL
double rows[1],cols[1],key[1],arrow[1]

objref boxl,bxh,bxv,this
objref mapdv // map an external vector instead of maintaining double
strdef tstr,tstr2,label

proc init () { local ii,jj
    rows=$1 cols=$2  
    if ($o3.size!=rows*cols) { 
        printf("XTAB ERR0: sz %d rows %d cols %d\n",$o3.size,rows,cols) return 
      }
    mapdv=$o3
    if (numarg()>3) key=$4
    if (numarg()>4) label=$s5 else label="TABLE"
    if (numarg()>5) arrow=$6 else arrow=0
    boxl=new List()
    boxl.remove_all
    bxv=new VBox()
    boxl.append(bxv)
    bxv.intercept(1)
    for ii=0,rows-1 {
        bxh=new HBox()
        boxl.append(bxh)
        bxh.intercept(1)
        for jj=0,cols-1 {
            bxv=new VBox()
            boxl.append(bxv)
            bxv.intercept(1)
            xpanel("")
            if (rows==1) sprint(tstr,"%d",jj) else if (cols==1) sprint(tstr,"%d",ii) else {
                sprint(tstr,"%d,%d",ii,jj) }
            sprint(tstr2,"notify(%s,%d,%d,%d,%d)",this.mapdv,ii*cols+jj,ii,jj,key)
            xpvalue(tstr,&mapdv.x(ii*cols+jj),arrow,tstr2)
            xpanel()
            bxv.intercept(0)
            bxv.map()
          }
        bxh.intercept(0)
        bxh.map()
      }
    boxl.object(0).intercept(0)
    boxl.object(0).map(label)
  }
endtemplate XTABL

// proc notify () { $o1.x[$2] *= 2 } // eg -- multiply any new input by 2
// proc notify () { printf("Table#%d:%s %d,%d CHANGED TO %g BY USER\n",$5,$o1,$3,$4,$o1.x[$2]) }
// objref a
// vec.resize(15)
// a=new XTABL(5,3,vec,0,"hello")
// vec.indgen(3)
// END xtabl.hoc
//================================================================
//================================================================
// INSERTED delta.hoc
// =Id=  delta.hoc,v 1.104 2003/01/27 18:25:41 billl Exp  

// load_file("delta.hoc")
// pattl.remove_all
// newpatts("VOR.in") 
// bplearn()
// rdpatts("VORtest.in")
// bprecg()
// load_file("string2.hoc")
// load_file("bkutils.hoc")
//* parameters and objects
objref sti,sth,sto,pattl,wt,wtih,wtho,wtdel,vo,vh,del,wtptr,wtio
objref odelih,odelho // old weight deltas for momentum

// i=input, h=hidden, o=output
{Ni=10  Nh=10 No=10}
{vo  = new Vector(No)     vh  = new Vector(Nh)     del = new Vector(No)}
{sti = new Vector(Ni)     sth = new Vector(Nh)     sto = new Vector(No)}
{wtih= new Vector(Ni*Nh)  wtho= new Vector(Nh*No) wtdel = new Vector(Ni*No) }
{odelih= new Vector(Ni*Nh)  odelho= new Vector(Nh*No)}
wt = wtho // use this for delta
wtptr=wtih // use this for monitoring a weight
rdm.uniform(-.5,.5)
pattl=new List()
wpind=0
BVBASE=0
alpha = 1
eta = 0
thresh= .5
iters=20000
npatts=pattnum=0

objref mso[7]
for ii=0,6 mso[ii]=new Vector(0)
{ID=0 EE=1 GG=2 HH=3 TT=4 xo=5 yo=6}
objref targ
targ=new Vector(No)

//* svpatts(filename) writes pattern file
proc svpatts () {  local ii,i,j
    tmpfile.wopen($s1)
    tmpfile.printf("Ni=%d\nNh=%d\nNo=%d\nnpatts=%d\nPATTERNS:\n",Ni,Nh,No,npatts)
    for ltrp(XO,YO,pattl) { XO.printf(tmpfile,"%g ") YO.printf(tmpfile,"%g ") }
    tmpfile.close
  }

//* newpatts() reads input files and chooses algo -- delta or bp
// eg newpatts("VOR/XOR.in")
first=0
proc newpatts () {  local ii,i,j,min
    if (Nh!=0 && first==0) { Ni+=1 Nh+=1 first=1} // bias units; units for learning threshold
    if (numarg()==1) {
        tmpfile.ropen($s1)
        temp_string_="//"  Ni=Nh=No=npatts=0
        while (strcmp(temp_string_,"PATTERNS:")!=0) {
            if (sfunc.substr(temp_string_,"//")!=0) execute(temp_string_)
            tmpfile.gets(temp_string_)
            chop(temp_string_)
          }
        if (npatts==0||Ni==0||No==0) {
            sprint(mesg,"File must define Ni,No,npatts\n",$s1)
            err(mesg)
            tmpfile.close return 
          } else {
            pattl.remove_all
            for ii=1,npatts {
                mso[ID].scanf(tmpfile,Ni) savevec(pattl,mso[ID])
                mso[ID].scanf(tmpfile,No) savevec(pattl,mso[ID])
              }
          }
      } else { // just read numbers for values
        i=1 j=1
        pattl.remove_all
        if ((Nh>0 && numarg()!=npatts*(Ni-1+No)) || (Nh==0 && numarg()!=npatts*(Ni+No))) { 
            sprint(mesg,"DELTA ERR0: Must set npatts first: %d",numarg()) err(mesg)  return 
          }
        if (Nh>0) min=1 else min=0
        for ii=1,npatts {
             mso[ID].resize(0)
            for j=min,Ni-1 { mso[ID].append($i) i+=1 } // need err checking for i<=numargs()
            savevec(pattl,mso[ID])
            mso[ID].resize(0)
            for j=1,No { mso[ID].append($i) i+=1 } // need err checking for i<=numargs()
            savevec(pattl,mso[ID])
          }
      } 
    chkpattsz()
    if (Nh==0) { // delta rule
        vo.resize(No) sti.resize(Ni) sto.resize(No)
        wt.resize(Ni*No) wtdel.resize(Ni*No)
        dinit()
      } else {
        vo.resize(No) vh.resize(Nh) del.resize(No)
        sti.resize(Ni) sth.resize(Nh) sto.resize(No)
        sti.x[0]=1 sth.x[0]=1 // bias unit for threshold learning
        wtih.resize(Ni*Nh) wtho.resize(Nh*No) wtdel.resize(Ni*No)
        odelih.resize(Ni*Nh) odelho.resize(Nh*No)
        bpinit()
      }
  }

//** rdpatts() reads in new patterns to test network on
// rdpatts(file[,flag]) -- flag==1 replace patts else append
// rdpatts(npatts,flag,patt_vals...) -- flag==1 replace patts else append
proc rdpatts () {  local ii, npsav, flag, i,j, nh, ni
    nh=Nh
    if (Nh!=0 && first==0) { Ni+=1 Nh+=1 first=1} // bias units; units for learning threshold
    if (Nh>0) ni=Ni-1  // don't count in the bias unit
    if (numarg()>1) flag=$2 else flag=0
    if (flag) pattl.remove_all  // flag: remove previous patterns
    temp_string_="//"  
    npsav = (pattl.count/2)
    if (numarg()<=2) {
        if (tmpfile.ropen($s1)==0) { sprint(mesg,"ERR: Can't open ",$s1) err(mesg)  return }
        while (strcmp(temp_string_,"PATTERNS:")!=0) {
            if (sfunc.substr(temp_string_,"//")==0) {} // ignore comments
            if (sfunc.substr(temp_string_,"Ni")>=0) first=0
            print temp_string_
            execute1(temp_string_) // set parameters
            tmpfile.gets(temp_string_)
            chop(temp_string_)
          }
        if (!first) { ni=Ni-1 first=1} // Ni has been reset
        if (npatts==0) {
            sprint(mesg,"File must define npatts")
            err(mesg)
            tmpfile.close return 
          } else {
            for ii=1,npatts {
                mso[ID].scanf(tmpfile,ni) savevec(pattl,mso[ID])
                mso[ID].scanf(tmpfile,No) savevec(pattl,mso[ID])
              }
          }
      } else {
        npatts=$1
        i=3 j=1
        for ii=1,npatts {
            mso[ID].resize(0)
            for j=1,ni { mso[ID].append($i) i+=1 } // need err checking for i<=numargs()
            savevec(pattl,mso[ID])
            mso[ID].resize(0)
            for j=1,No { mso[ID].append($i) i+=1 } // need err checking for i<=numargs()
            savevec(pattl,mso[ID])
          }
      }
    if (! flag) npatts += npsav
    chkpattsz(ni)
    showpatts(0)
  }

//* delta rule learning
proc dinit () {
    if (numarg()==2) rdm.uniform($1,$2) else rdm.uniform(-.5,.5)
    wt.setrand(rdm)
    mesg=""
    rdm.uniform(0,npatts) // for picking pattern to present
  }

proc dlearn () { local ii,res
    if (Nh!=0) { printf("Nh=%d: need bplearn()\n",Nh) return }
    res=0
    for(ii=0;res<npatts && ii<iters;ii+=1) res=dlearn1()
    print ii
  }

//** dlearn1()
func dlearn1 () { local cnt,only
    cnt=0
    if (numarg()==1) only = $1 else only = -1
    for ltrp(XO,YO,pattl) {
        if (i1==only || only==-1) {
            sti.copy(XO) sto.copy(YO)
            vo.mmult(wt,sti)  
            vo.thresh(thresh)
            vo.sub(sto) // subtract target
            if (vo.max==vo.min && vo.min==0) { cnt+=1 // all correct
              } else { // wtdel = alpha*(sto - W*sti) X sti ; [where X is outer product]
                vo.mul(-alpha)
                wtdel.outprod(vo,sti)
                wt.add(wtdel)
                // printf("%d: ",i1) wt.mprintf(No,Ni)
              }
            // vlk(wt)
          }
      }
    return cnt
  }

//** drecall()
proc drecall () { local cnt,only,np
    if (Nh!=0) { printf("Nh=%d: need bplearn()\n",Nh) return }
    if (numarg()==1) only = $1 else only = -1
    cnt=0
    for ltrp(XO,YO,pattl) {
        if (i1==only || only==-1) {
            sti.copy(XO) sto.copy(YO)
            vo.mmult(wt,sti)  
            vo.thresh(thresh)
            vlk(sti,0) printf(" -> ") vlk(vo,0) 
            if (sto.eq(vo)) cnt+=1 else {  printf("( ") vlk(sto,0) printf(")\n") }
          }
      }
    if (only>-1) { np=1 if (cnt==1) mesg="Correct" else mesg="Wrong" } else np=npatts
    printf("Correct: %d/%d \nW matrix:\n",cnt,np) wt.mprintf(No,Ni)
  }

//* backprop rule learning
//**  bpinit(), chkpattsz()
//*** chkpattsz()
func chkpattsz () { local err, ni
    err=-1
    if (numarg()==1) ni=$1 else {if (Nh>0) ni=Ni-1 else ni=Ni}
    for ltrp(XO,YO,pattl) if (ni!=XO.size || No!=YO.size) err=i1
    if (err!=-1) { 
        printf("ERROR: all patts must be same size: %d,%d: \n",ni,No) 
        showpatts(err) 
        return 0 } else return 1
  }

//*** bpinit()
proc bpinit () { local ii
    mesg=""
    if (numarg()==2) rdm.uniform($1,$2) else rdm.uniform(-.5,.5)
    wtih.setrand(rdm)
    for ii=0,Ni-1 wtih.mset(0,ii,Ni,0) // clear the row projecting to the bias unit
    wtho.setrand(rdm)
    showpatts(pattnum=0)
    rdm.uniform(0,npatts) // for picking pattern to present
  }

//** bpforward
// implement learning of the theta term as an additional input unit that is always 1 and an 
// additional hidden unit = 1
proc bpforward () { local cnt
    sth.mmult(wtih,sti)  
    sth.apply("sigmoid")
    sth.x[0]=1 // also sti.x[0]=1 // "theta" weights
    sto.mmult(wtho,sth)  
    sto.apply("sigmoid")
  }

//** bplearn()
{chkpt=500 errt=1e-7 bpstop=0}
func bplearn () { local cnt,ii,jj,err,savgr
    bpstop=graph_flag=0
    if (Nh==0) { printf("Nh=0: need dlearn()\n") return }
    if (numarg()==1) if ($1==0) bpinit() else iters=$1
    for(jj=0;jj<iters && bpstop==0;jj+=1) {
        // ii=int(rdm.repick)
        ii = jj%npatts
        XO=pattl.object(2*ii) YO=pattl.object(2*ii+1)
        sti.copy(XO,1)
        bpforward()
        bpbp(YO) // backpropagate error
        chwts()
        if (jj%chkpt==chkpt-1) { 
            if (bprecall(-1)<errt) bpstop=1 
            sprint(mesg,"ITER# %d",jj)
            doNotify()
          }
      }
    printf("ITERS: %d\nERR: %g\n",jj,bprecall())
    graph_flag=1
    return(jj)
  }

// expand sigmoid out so it's easy to reach 0 and 1 but now f'(net) != o*(1-o)
{sigmul=1.4 sigsub=(sigmul-1)/2}
func sigmoid () { return sigmul/(1+exp(-$1))-sigsub }

//*** bpbp() -- back-propagation step
// reverse the expansion in back prop so that get correct 1st deriv
proc bpbp () { local tmp
    if (numarg()==1) targ.copy($o1)
    del.resize(0) del.copy(targ) del.sub(sto) // targ-outp
    if (graph_flag) {printf("Output vs. target deltas: ") vlk(del)}
    // vo=sto*(1-sto)*del
    vo.bpeval(sto,del,sigsub,sigmul)
    wtdel.resize(wtho.size)
    wtdel.transpose(wtho,No,Nh) // N_o rows and N_h columns
    vh.mmult(wtdel,vo) // backprop error
    del.resize(0) del.copy(vh) vh.bpeval(sth,del,sigsub,sigmul)
    vh.x[0]=0
    if (graph_flag) {printf("Hidden unit deltas: ") vlk(vh)}
  }

//*** chwts()
proc chwts () {
    if (eta>0) { wtdel.copy(odelho) wtdel.mul(eta) }
    odelho.outprod(vo,sth)
    odelho.mul(alpha)
    wtho.add(odelho)
    if (eta>0) { wtho.add(wtdel)
        wtdel.copy(odelih) wtdel.mul(eta) }
    odelih.outprod(vh,sti)
    odelih.mul(alpha)
    wtih.add(odelih)
    if (eta>0) wtih.add(wtdel)
  }  

//*** bpdiffr(res,outp,del) res=outp*(1-outp)*del
// replaced with compiled function: res.bpeval(outp,del)
proc bpdiffr () { local tv
    $o1.copy($o2) $o1.mul(-1) $o1.add(1) // (1-outp)
    $o1.mul($o3) // (1-outp)*del
    $o1.mul($o2) // outp*(1-outp)*del
  }

//** bprecall()
err0=.01
// bprecall(arg) arg=-1 means quiet else arg means just run that pattern
func bprecall () { local cnt,error,er,sh,verbose,tpatt
    verbose=1
    if (numarg()==1) {sh=$1 if (sh==-1) verbose=0} else sh=-1
    if (Nh==0) { printf("Nh=0: need dlearn()\n") return }
    cnt=0 error=0 
    tpatt = pattl.count
    for ltrp(XO,YO,pattl) {
        if (sh==-1 || i1==sh) {
            sti.copy(XO,1) // sti.x[0]=1 reserved for thresh weight (bias)
            bpforward()    // forward pro
            del.copy(YO) del.sub(sto) // targ-outp
            er=0.5*del.dot(del)  error+=er // so that 1st deriv will get rid of 1/2
            if (error<err0) cnt+=1
            if (verbose) { vlk(XO,0) printf(" -> ") vlk(sto,0) 
                printf("( ") vlk(YO,0) printf("\tERR: %g)\n",er) }
          }
        bpcallback()
      }
    if (verbose) {
        if (sh>-1) {printf("HIDDEN: ") vlk(sth)}
        printf("Correct: %d/%d \n",cnt,npatts)
      }
    return sumerr=error/tpatt
  }
proc bpcallback () {}

proc bppresent () {
    sti.copy($o1,1)
    bpforward()
    vlk(sto)
  }

proc bpsstep () { local cnt,error,er,sh,ii
    if (Nh==0) { printf("Nh=0: need dlearn()\n") return }
    sh=$1
    XO=pattl.object(2*sh) YO=pattl.object(2*sh+1)
    sti.copy(XO,1) // sti.x[0]=1 reserved for thresh weight
    bpforward()    // forward pro
    printf("INPUT: ")
    for vtr(&x,sti) if (i1>0) printf("%12.04g",x)
    printf("\n________________________________________________________________\nWT_ih*S_i\n")
    mpfv2(wtih,sti,0)
    printf("________________________________________________________________\nHIDDEN: ")
    for vtr(&x,sth) if (i1>0) printf("%12.04g",x)
    printf("\n________________________________________________________________\nWT_ho*S_h\n")
    mpfv2(wtho,sth)
    printf("________________________________________________________________\nOUTPUT: ")
    for vtr(&x,sto) printf("%12.04g",x)
    printf("\n________________________________________________________________\nTARGET")
    for vtr(&x,YO) printf("%12.04g",x)
    del.copy(YO) del.sub(sto) // targ-outp
    printf("\n________________________________________________________________\nDELTA")
    for vtr(&x,del) printf("%12.04g",x)
    printf("\n________________________________________________________________\n")
    er=0.5*del.dot(del)
    printf("ERR: %12.04g\n",er)
  }

//*** bprecg() graph recall results
// eye movement (circle), visual target movement (triangle),
// head movement (square), goal movement (filled circle) -- global YO
objref legl
legl=new List()
proc bprecg () { local ii,jj,flag,max
    if (numarg()==2) flag=$2 else flag=0
    legl.remove_all
    for scase("actual","desired","head","duck") legl.append(new String2(temp_string_,sym[i1+2].s))
    $o1.erase_all()
    if (npatts==0) return
    $o1.size(0,npatts+3,-2,2)
    $o1.xaxis(3)
    $o1.xaxis(0,npatts+3,0,1,0,0,0)
    $o1.yaxis(-2,2)
    grrtsize($o1)
    if (flag!=2) for ii=ID,TT mso[ii].resize(0)
    mso[ID].indgen(1,npatts,1)
    sprint(temp_string_,"proc bpcallback() { %s.append(sto.x[1]-sto.x[0]) ",mso[EE])
      sprint(temp_string_,"                 %s %s.append(sti.x[3]-sti.x[2])",temp_string_,mso[TT])
      sprint(temp_string_,"                 %s %s.append(sti.x[4]-sti.x[1])",temp_string_,mso[HH])
      sprint(temp_string_,"                 %s %s.append( YO.x[1]- YO.x[0])}",temp_string_,mso[GG])
    execute1(temp_string_)
    if (flag!=2) bprecall(-1) // else just fake it with values placed in mso[] vectors
    if (flag==1) max=GG else max=TT
    for (ii=max;ii>=EE;ii-=1) {
        jj=ii-EE
        mso[ii].mark($o1,mso[ID],sym[jj+2].s,20,jj+1,3) // y-axis is VOR; x is pursuit
      }
    legend($o1,legl)
    execute1("proc bpcallback() { }")
  }

//** legend()
proc legend () { local xv,yv,yv2,xw,yw,sz1,sz2,sz3,sz4,wd,ht,xexpand,yexpand
    sz1=$o1.size(1) sz2=$o1.size(2) sz3=$o1.size(3) sz4=$o1.size(4)
    xexpand=1.2 yexpand=1.02
    sz1*=xexpand sz2*=xexpand sz3*=yexpand sz4*=yexpand
    wd=sz2-sz1 ht=sz4-sz3 // need to translate between world and view coordinates
    $o1.align(.5,1)
    for ltr(XO,$o2) { // string2 list of strings and marks
        xv=.8 yv=1-.08*i1
        xv=.8 yv2=.98-.07*i1
        $o1.label(xv,yv2-.02,XO.s)
        xw=sz1+xv*wd  yw=sz3+yv*ht
        $o1.mark(xw,yw,XO.t,20,i1+1,3)
      }
    $o1.label(0.05,0.85,"Right")
    $o1.label(0.05,0.35,"Left")
  }

//** print out a matrix*vec multiplication
proc mpfv2 () { local i,j,rows,cols,beg
    cols=$o2.size rows=$o1.size/cols
    if (numarg()==3) beg=1 else beg=0
    for (i=beg;i<rows;i+=1) {
        for (j=0;j<cols;j+=1) {
            if (j==1) printf("|  ")
            printf("%-12.03g  ",$o1.x[i*cols+j])
          }
        if (i<$o2.size) printf("  |   %12.04g \n",$o2.x[i]) else printf("  |\n")
      }
    if (cols>rows) for i=rows,cols-1 printf("%70.04g \n",$o2.x[i])
  }

proc bpshowhid () {
    for ltrp(XO,YO,pattl) {
        sti.copy(XO,1) // sti.x[0]=1 reserved for thresh weight
        bpforward()    // forward pro
        for vtr(&x,sti) if (i1>0) printf("%g ",x)
        printf("  ->")
        for vtr(&x,sth) if (i1>0) printf("%12.04g",x)
        printf("  ->")
        for vtr(&x,sto) printf("%12.04g",x)
        print ""
      }
  }

//** showwts(), showpatts()
proc showwts () {
    if (Nh==0) {
        printf("WT I-O: \n") wt.mprintf(No,Ni)
      } else {
        printf("WT I-H: \n") wtih.mprintf(Nh,Ni)
        printf("WT H-O: \n") wtho.mprintf(No,Nh)
      }
  }

proc showpatts () { local sh,flag
    flag=0
    if (numarg()==1) sh=pattnum=$1 else flag=1
    if (sh>pattl.count/2-1) sh=pattnum=pattl.count/2-1
    if (sh<0) sh=pattnum=0
    if (pattl.count()==0) { err("No patterns defined") return }
    for ltrp(XO,YO,pattl) { 
        if (flag || i1==sh) { 
            if (Nh>0) {
                sti.copy(XO,1) targ.copy(YO)
                ind.resize(Ni+Nh+2*No)
                ind.fill(0)
                ind.copy(XO,1)
                ind.copy(YO,Ni+Nh+No)
                nc.dispcol(ind)
                doEvents()
                vlk(XO,0) printf(" -> ") vlk(YO) 
              } else {
                sti.copy(XO,0) targ.copy(YO)
                ind.resize(Ni+2*No)
                ind.fill(0)
                ind.copy(XO,0)
                ind.copy(YO,Ni+No)
                nc.disp(ind)
              }
          }
      }
  }

//*** bpgrwts() graph weights
proc bpgrwts () {
    mso[ID].resize(0) mso[xo].resize(0) mso[yo].resize(0)
    if (numarg()==3) { mso[xo].copy($o2) mso[yo].copy($o3) } else {
        bpgetwts(mso[ID],mso[xo],mso[yo]) }
    mso[yo].mark($o1,mso[xo],"O",10,gcl,1)
    mso[xo].resize(0) mso[yo].resize(0)
    mso[xo].append(1,-1) mso[yo].append(-1,1)
    mso[yo].line($o1,mso[xo])
  
    graphItem.label(0.85,0.55,"duck")
    graphItem.label(0.4,0.05,"head")
    graphItem.label(0.95,0.55,"R")
    graphItem.label(0.05,0.55,"L")
    graphItem.label(0.55,0.95,"R")
    graphItem.label(0.55,0.05,"L")
  
    grrtsize()
  }

//*** bpgetwts(scr,VOR,pursuit) returns vectors with wt differences
proc bpgetwts () { local ii
    $o1.resize(Ni)
    for ii=1,Nh-1 {
        $o1.mrow(wtih,ii,Ni)
        $o2.append($o1.x[3]-$o1.x[2]) // pursuit influence
        $o3.append($o1.x[4]-$o1.x[1]) // VOR influence
      }
  }

//* experiments
proc svinitwts () {
    rdm.ACG(3804485) // seed
    veclist.remove_all
    for ii=0,9 { bpinit() savevec(wtih,wtho) }
  }

proc bpmulti () { local ii
    vec0.resize(0)
    rdm.ACG(3804485) // seed
    for ii=0,9 { 
        print ii
        wtih.copy(veclist.object(2*ii)) // restore the weights
        wtho.copy(veclist.object(2*ii+1)) // restore the weights  
        vec0.append(bplearn())
      }
  }    

//* adapted from MPNfdfwd1.hoc33
objref box[3],nc,tb[2]

lmv=.7 // movement for drawing connection lines
mvmax=5
colr=2

//* mkitbp()
ncview=0
proc layers() {} // for labeling layers
proc extra () {} // for adding extra buttons
proc mkitbp () { 
    box[2]=new VBox() // main
    box[0]=new HBox() // main
    box[2].intercept(1)
    qhp(1)
    box[0].intercept(1)
    drawnet()
    nc.g.view(0,0,100,100,600,200,500,300) 
    nc.gcmap(BVBASE,1)  // graph of color map
    xpanel("")
    xvarlabel(mesg)
    // xcheckbox("Show units",&graph_flag)
    xcheckbox("Stop processing",&bpstop)
    xbutton("Learn","bplearn()")
    xvalue("Max iter","iters")
    xvalue("Sum Err","sumerr")
    xvalue("Targ err","errt",1)
    xvalue("Forward prop","pattnum",1,"showpatts(pattnum) sprint(mesg,\"ERR: %g\",bprecall(pattnum)) mapsts()")
    xbutton("Back prop","bpbp() chwts()")
    extra()
    xmenu("Weights")
    nc.g.menu_tool("Draw weight","drawwt")
    xbutton("Save weights","svrd(0,\"weights\",\"mat\")")
    xbutton("Read weights","svrd(1,\"weights\",\"mat\")")
    xbutton("Draw one (click on units)","nc.g.exec_menu(\"Draw weight\")")
    xbutton("Draw all","shmat(3)")
    xbutton("Show in->hid matrix","shmat(0)")
    xbutton("Show hid->out matrix","shmat(1)")
    xbutton("Randomize weights","bpinit()")
    xbutton("Clear graph","shmat(4)")
    xmenu()
    xmenu("Patterns&Units")
    xbutton("Clear all patterns","pattl.remove_all() npatts=0")
    xbutton("Load patterns","svrd(1,\"Patterns\",\"ptn\",1) rdpatts(filename)")
    xbutton("Save patterns","svrd(0,\"Patterns\",\"ptn\",1) svpatts(filename)")
    xbutton("Show patterns","showpatts(pn=(pn+1)%npatts)")
    xbutton("Show values","label()")
    xbutton("Label layers","layers()")
    xbutton("Set inputs","shmat(2)")
    xbutton("Clear graph","shmat(4)")
  
    xmenu()
    xvalue("Edit pattern #","pn",1,"pnum()")
    xvalue("Weight monitor","wtptr.x[wpind]")
    xpanel()
    box[0].intercept(0)
    box[0].map("Back-prop network")
    box[2].intercept(0)
    box[2].map("Back-prop network")
    // initvals()
  }

//** mkitdelta()
proc mkitdelta () { 
    box[2]=new VBox() // main
    box[0]=new HBox() // main
    box[2].intercept(1)
    qhp(1)
    box[0].intercept(1)
    extra0()
    xpanel("")
    xvarlabel(mesg)
    // xcheckbox("Show units",&graph_flag)
    xvalue("Recall","pattnum",1,"delrec()")
    xvalue("Delta","pattnum",1,"delshow()")
    extra()
    xmenu("Weights")
    xbutton("Save weights","svrd(0,\"weights\",\"mat\")")
    xbutton("Read weights","svrd(1,\"weights\",\"mat\")")
    xbutton("Draw weights","shmat(3)")
    xbutton("Show weight matrix","shmat(0)")
    xbutton("Randomize weights","dinit()")
    xbutton("Clear graph","shmat(4)")
    xmenu()
    xmenu("Patterns&Units")
    xbutton("Clear all patterns","pattl.remove_all() npatts=0")
    xbutton("Load patterns","svrd(1,\"Patterns\",\"ptn\",1) rdpatts(filename)")
    xbutton("Save patterns","svrd(0,\"Patterns\",\"ptn\",1) svpatts(filename)")
    xbutton("Show values","label()")
    xbutton("Set inputs","shmat(2)")
    xmenu()
    xvalue("Edit pattern #","pn",1,"pnum()")
    xpanel()
    box[0].intercept(0)
    box[0].map("Delta-rule learning: AND ")
    box[2].intercept(0)
    box[2].map("Delta-rule learning: AND ")
    // initvals()
  }

proc delrec () {
    if (pattnum<0) pattnum=0
    if (pattnum>npatts-1) pattnum=npatts-1
    showpatts(pattnum)
    drecall(pattnum) 
    mapstd()
  }

proc delshow () {
    if (pattnum<0) pattnum=0
    if (pattnum>npatts-1) pattnum=npatts-1
    showpatts(pattnum) 
    dlearn1(pattnum) 
    drecall(pattnum) 
    mapstd()
  }

//** drawwt()
wtclk=1e10
proc drawwt () { local ii,wh,wht
    if ($1==2) { 
        wh=nc.which($2,$3) 
        if (wtclk==1e10) { wtclk=wh print "Click another" } else {
            ii=wtclk jj=wh wtclk=1e10
            if (ii<Ni && jj>=Ni && jj<Ni+Nh) {
                wht=wtih.mget(jj-Ni,ii,Ni)
                wtptr=wtih  wpind=(jj-Ni)*Ni+ii
                printf("Monitoring input %d to hid %d: %g\n",ii,jj-Ni,wht)
              } else if (ii>=Ni && ii<Ni+Nh && jj>=Ni+Nh && jj<Ni+Nh+No) {
                wht=wtho.mget(jj-Nh-Ni,ii-Ni,Nh) 
                wtptr=wtho wpind=(jj-Nh-Ni)*Nh+ii-Ni
                printf("Monitoring hid %d to output %d: %g\n",ii-Ni,jj-Nh-Ni,wht)
              } else { print "Not connected" return }
            nc.line(ii,jj,1,"O",8,0,0,wht) 
          }
      }
  }

//** pnum() -- change patterns
pn=1
proc pnum () { local n
    if (numarg()==1) n=$1 else n=pn
    pattnum=pn
    if (n<0) n=0  if(n>npatts)n=npatts
    if (Nh>0) {
        tb[0]=new XTABL(Ni-1,1,pattl.object(2*n),2,"INPUT PATTERN",1)
      } else {
        tb[0]=new XTABL(Ni,1,pattl.object(2*n),2,"INPUT PATTERN",1)
      }
    tb[1]=new XTABL(No,1,pattl.object(2*n+1),3,"OUTPUT PATTERN",1)
  }

proc reset () { 
    initialize() 
    if (isobj(box[1],"VBox")) box[1].unmap
  }

proc notify () {
    if ($5==0 && $3==0) {
        err("Unused row (projection to bias unit)")
        $o1.x[$2]=0
      }
  }

//* shmat() show weight matrix
proc shmat () { local num,ii,wht,min
    num=$1
    if (num==0) {
        if (Nh==0) {
            tb[0]=new XTABL(No,Ni,wt,1,"Weight Matrix",1)
          } else {
            tb[0]=new XTABL(Nh,Ni,wtih,0,"INPUT TO HIDDEN (TOP ROW INACTIVE)",0)
            mesg="in->hid matrix: TOP ROW INACTIVE"
          }
      } else if (num==1) {
        tb[1]=new XTABL(No,Nh,wtho,1,"HIDDEN TO OUTPUT WEIGHT MATRIX",1)
      } else if (num==2) {
        xpanel("Set input units")
        xbutton("Forward prop","bpforward() mapsts()")
        if (Nh>0) min=2 else min=1
        for ii=min,Ni { // ignore first unit which is bias unit
            sprint(tstr,"Input #%d",ii-1)
            sprint(tstr2,"sti.x[%d]",ii-1)
            if (Nh>0) temp_string_="mapsts()" else temp_string_="mapstd()"
            xvalue(tstr,tstr2,1,temp_string_)
          }
        xpanel()
      } else if (num==3) { // show weights
        if (Nh==0) {
            for ii=0,Ni-1 for jj=0,No-1 {
                wht=wt.mget(jj,ii,Ni)
                if (wht<0) colr=3 else colr=2
                nc.line(ii,jj+Ni,colr,"O",8,-lmv*rdm.repick,mvmax*lmv-2*lmv*rdm.repick,wht) 
              }
          } else {
            for ii=0,Ni-1 for jj=0,Nh-1 {
                wht=wtih.mget(jj,ii,Ni)
                if (wht<0) colr=3 else colr=2
                nc.line(ii,jj+Ni,colr,"O",8,-lmv*rdm.repick,mvmax*lmv-2*lmv*rdm.repick,wht) 
              }
            for ii=0,Nh-1 for jj=0,No-1 { 
                wht=wtho.mget(jj,ii,Nh)
                if (wht<0) colr=3 else colr=2
                nc.line(ii+Ni,jj+Ni+Nh,colr,"O",8,-lmv*rdm.repick,mvmax*lmv-2*lmv*rdm.repick,wht) 
              } 
          }
      } else if (num==4) {
        nc.g.erase_all
        nc.vals.fill(0)
        if (Nh>0) nc.vals.x[0]=nc.vals.x[Ni]=1 // bias units
        nc.disp()
      }
  }

// just draw the connectivity without the bias units
proc drawconn () {
    for ii=1,Ni-1 for jj=1,Nh-1 {
        wht=wtih.mget(jj,ii,Ni)
        if (wht<0) colr=3 else colr=2
        nc.line(ii,jj+Ni,colr,"O",8,-lmv*rdm.repick,mvmax*lmv-2*lmv*rdm.repick)
      }
    for ii=1,Nh-1 for jj=0,No-1 { 
        wht=wtho.mget(jj,ii,Nh)
        if (wht<0) colr=3 else colr=2
        nc.line(ii+Ni,jj+Ni+Nh,colr,"O",8,-lmv*rdm.repick,mvmax*lmv-2*lmv*rdm.repick)
      } 
  }


proc initvals () {
    mkmat()
    nc.vals.fill(0)
    // for ii=0,1 nc.vals.x[ii]=1
    nc.disp
  }

stp=0
//* dispwait () -- display values and wait
proc dispwait () { local ii
    mesg=$s1
    if (numarg()==2) for ii=0,$o2.size-1 sprint(mesg,"%s %g ",mesg,$o2.x[ii])
    print mesg
    if (stp==0) box[1]=new VBox()
    box[1].unmap()
    box[1].intercept(1)
    xpanel("")
    if (stp==0) xbutton("Continue","cont()")
    xlabel(mesg)
    xpanel()
    box[1].intercept(0)
    box[1].map("Message",0,0,-1,-1)
  }

//** parsr()
proc parsr () { local flag
    flag=$1
    if (numarg()==2) { filename=$s2 
        if (flag==0) tmpfile.wopen(filename) else tmpfile.ropen(filename) }
    if (flag==0) { 
        vec.resize(Ni)
        for ii=0,Nh-1 {
            vec.mrow(wtih,ii,Ni) 
            vec.printf(tmpfile,"%f ")
          }
        vec.resize(Nh)
        for ii=0,No-1 {
            vec.mrow(wtho,ii,Nh) 
            vec.printf(tmpfile,"%f ")
          }
      } else { 
        wtih.scanf(tmpfile,Ni*Nh) 
        wtho.scanf(tmpfile,Nh*No) 
      }
    tmpfile.close
  }

//* mapsts()  -- display state values on the picture
proc mapsts () {
    ind.copy(sti,0)
    ind.copy(sth,Ni)  
    ind.copy(sto,Ni+Nh)
    nc.dispcol(ind)
  }

//** for delta
proc mapstd () {
    ind.resize(Ni+2*No)
    ind.copy(sti,0)
    ind.copy(vo,Ni)
    nc.disp(ind)
    shmat(3)
  }
// END delta.hoc
//================================================================

{Ni=4 Nh=9 No=2}

proc drawnet () {
    nc=new NCIRC(-20,90, \
                 15,10, 35,10, 65,10, 85,10,\
                 -20,110,\
                 0,50, 12,50, 24,50, 36,50, 48,50, 60,50, 72,50, 84,50, 96,50,\
                 30,90, 70,90,\
                 30,100, 70,100) // targets
    nc.low=0 nc.high=1
    nc.symsz=30
  }

//* layers() -- label the layers
proc layers () {
    nc.g.label(50,0,"INPUTS",   1,2,.5,.5,1)
    nc.g.label( 7,12,"VOR  L",   1,2,.5,.5,1)
    nc.g.label(93,12,"R  VOR",   1,2,.5,.5,1)
    nc.g.label(50,12,"L  pursuit  R",   1,2,.5,.5,1)
    nc.g.label(50,40,"HIDDEN",  1,2,.5,.5,1)
    nc.g.label(30,92,"L",  1,2,.5,.5,1)
    nc.g.label(70,92,"R",  1,2,.5,.5,1)
    nc.g.label(50,90,"OUTPUT",  1,2,.5,.5,1)
    nc.g.label(50,82,"(eye movement)",  1,2,.5,.5,1)
    nc.g.label(50,100,"TARGETS",1,2,.5,.5,1)
  }

proc label () { nc.lscale=nc.colr=1 nc.labl() }
proc reset () { np() bpinit() layers() }

proc mkmat () {
    wtih.resize(Nh*Ni)
    wtho.resize(No*Nh)
  }

mkitbp()
// newpatts("VOR.in") 
{alpha = 1 Ni=4 Nh=9 No=2 npatts=4 sigmul=3 sigsub=(sigmul-1)/2}
proc np () { newpatts(0,.5,.5,1,1,0,1,.5,.5,0,0,1,.5,1,0,.5,1,0,.5,0,1,.5,0,1) }
reset()
// END SUP_VOR.hoc
//================================================================