diff --git a/src/experimental/symbolic/gauge.nim b/src/experimental/symbolic/gauge.nim
new file mode 100644
index 0000000..a933a32
--- /dev/null
+++ b/src/experimental/symbolic/gauge.nim
@@ -0,0 +1,246 @@
+import experimental/symbolic/graph
+import ../../gauge, physics/qcdTypes
+
+#
+# gauge support
+#
+
+# TODO: needs more thought
+type Gauge = seq[DLatticeColorMatrixV]
+
+type SymNodeValueGauge = ref object of SymNodeValueConcrete
+  gaugeValue: Gauge
+
+method getGauge*(v: SymNodeValue): Gauge {.base.} =
+  raiseValueError("Custom method required for value: " & $v)
+
+method getGauge*(v:SymNodeValueGauge): Gauge = v.gaugeValue
+
+method `$`*(v: SymNodeValueGauge): string = "gaugeValue"
+
+proc assign*(z: SymNode, v: Gauge) =
+  z.assign SymNodeValueGauge(gaugeValue: v)
+
+method copySymNodeValue*(v: SymNodeValueGauge): SymNodeValueGauge =
+  # TODO: we don't need this, if we don't take gradient after assign
+  let z = v.gaugeValue[0].l.newGauge
+  threads:
+    for mu in 0..<z.len:
+      z[mu] := v.gaugeValue[mu]
+  SymNodeValueGauge(gaugeValue: z)
+
+method identSymNodeValue*(z: SymNodeValueGauge, x: SymNodeValueGauge) =
+  threads:
+    for mu in 0..<z.gaugeValue.len:
+      z.gaugeValue[mu] := x.gaugeValue[mu]
+
+method identAllocateSymNodeValue*(z: SymNode, x: SymNodeValueGauge) =
+  # TODO: leave it uninitialized?
+  z.value = SymNodeValueGauge(gaugeValue: x.gaugeValue[0].l.newGauge)
+
+method iaddSymNodeValue*(z: SymNodeValueGauge, x: SymNodeValueGauge, y: SymNodeValueGauge) =
+  threads:
+    for mu in 0..<z.gaugeValue.len:
+      z.gaugeValue[mu] := x.gaugeValue[mu] + y.gaugeValue[mu]
+
+method iaddAllocateSymNodeValue*(z: SymNode, x: SymNodeValueGauge, y: SymNodeValueGauge) =
+  z.value = SymNodeValueGauge(gaugeValue: x.gaugeValue[0].l.newGauge)
+
+method norm2SymNodeValue*(z: SymNodeValue, x: SymNodeValue) {.base.} =
+  raiseErrorBaseMethod("args:\n  " & z.repr & "\n  " & x.repr)
+
+method norm2AllocateSymNodeValue*(z: SymNode, x: SymNodeValue) {.base.} =
+  raiseErrorBaseMethod("args:\n  " & z.nodeRepr & "\n  " & x.repr)
+
+method imulSymNodeValue*(z: SymNodeValue, x: SymNodeValue, y: SymNodeValue) {.base.} =
+  raiseErrorBaseMethod("args:\n  " & z.repr & "\n  " & x.repr & "\n  " & y.repr)
+
+method imulAllocateSymNodeValue*(z: SymNode, x: SymNodeValue, y: SymNodeValue) {.base.} =
+  raiseErrorBaseMethod("args:\n  " & z.nodeRepr & "\n  " & x.repr & "\n  " & y.repr)
+
+method norm2SymNodeValue*(z: SymNodeValueFloat, x: SymNodeValueGauge) =
+  threads:
+    var t = 0.0
+    for mu in 0..<x.gaugeValue.len:
+      t += x.gaugeValue[mu].norm2
+    threadMaster: z.floatValue = t
+
+method norm2AllocateSymNodeValue*(z: SymNode, x: SymNodeValueGauge) =
+  z.value = SymNodeValueFloat()
+
+method imulSymNodeValue*(z: SymNodeValueFloat, x: SymNodeValueFloat, y: SymNodeValueFloat) =
+  z.floatValue = x.floatValue * y.floatValue
+
+method imulAllocateSymNodeValue*(z: SymNode, x: SymNodeValueFloat, y: SymNodeValueFloat) =
+  z.value = SymNodeValueFloat()
+
+method imulSymNodeValue*(z: SymNodeValueGauge, x: SymNodeValueFloat, y: SymNodeValueGauge) =
+  threads:
+    for mu in 0..<z.gaugeValue.len:
+      z.gaugeValue[mu] := x.floatValue * y.gaugeValue[mu]
+
+method imulAllocateSymNodeValue*(z: SymNode, x: SymNodeValueFloat, y: SymNodeValueGauge) =
+  z.value = SymNodeValueGauge(gaugeValue: y.gaugeValue[0].l.newGauge)
+
+method imulSymNodeValue*(z: SymNodeValueGauge, x: SymNodeValueGauge, y: SymNodeValueFloat) =
+  threads:
+    for mu in 0..<z.gaugeValue.len:
+      z.gaugeValue[mu] := x.gaugeValue[mu] * y.floatValue
+
+method imulAllocateSymNodeValue*(z: SymNode, x: SymNodeValueGauge, y: SymNodeValueFloat) =
+  z.value = SymNodeValueGauge(gaugeValue: x.gaugeValue[0].l.newGauge)
+
+method imulSymNodeValue*(z: SymNodeValueGauge, x: SymNodeValueGauge, y: SymNodeValueGauge) =
+  threads:
+    for mu in 0..<z.gaugeValue.len:
+      z.gaugeValue[mu] := x.gaugeValue[mu] * y.gaugeValue[mu]
+
+method imulAllocateSymNodeValue*(z: SymNode, x: SymNodeValueGauge, y: SymNodeValueGauge) =
+  z.value = SymNodeValueGauge(gaugeValue: x.gaugeValue[0].l.newGauge)
+
+#
+# more functions
+#
+
+proc mul*(x: SymNode, y: SymNode): SymNode
+
+proc mulForward(z: SymNode) =
+  imulSymNodeValue(z.value, z.inputs[0].value, z.inputs[1].value)
+
+proc mulAllocate(z: SymNode) =
+  imulAllocateSymNodeValue(z, z.inputs[0].value, z.inputs[1].value)
+
+proc mulBackward(z: SymNode, i: int, dep: SymNode): SymNode =
+  if i != 0 and i != 1:
+    raiseError("mul has 2 operands, got i = " & $i)
+  let g = z.gradientDependentOrNil dep
+  if g == nil:
+    if i == 0:
+      return z.inputs[1]
+    else:
+      return z.inputs[0]
+  else:
+    # assume noncommuting operands here
+    if i == 0:
+      return g.mul z.inputs[1]
+    else:
+      return z.inputs[0].mul g
+
+proc mul*(x: SymNode, y: SymNode): SymNode =
+  newSymNode(
+    inputs = @[x, y],
+    forward = mulForward,
+    allocateValue = mulAllocate,
+    backward = mulBackward,
+    name = "mul")
+
+proc norm2Forward(z: SymNode) =
+  norm2SymNodeValue(z.value, z.inputs[0].value)
+
+proc norm2Allocate(z: SymNode) =
+  norm2AllocateSymNodeValue(z, z.inputs[0].value)
+
+proc norm2Backward(z: SymNode, i: int, dep: SymNode): SymNode =
+  if i != 0:
+    raiseError("norm2 has 1 operand, got i = " & $i)
+  let g = z.gradientDependentOrNil dep
+  let two = newSymNode(value = SymNodeValueFloat(floatValue: 2.0), name = "Two[norm2]")
+  if g == nil:
+    return two.mul z.inputs[0]
+  else:
+    let m = two.mul g
+    m.tagUpdate
+    return m.mul z.inputs[0]
+
+proc norm2*(x: SymNode): SymNode =
+  newSymNode(
+    inputs = @[x],
+    forward = norm2Forward,
+    allocateValue = norm2Allocate,
+    backward = norm2Backward,
+    name = "norm2")
+
+#
+# SymNode and float
+#
+
+# TODO perhaps make a generic converter: anyValue -> SymNode
+
+proc `*`(x: float, y: SymNode): SymNode =
+  return newSymNode(value = SymNodeValueFloat(floatValue: x)).mul y
+
+proc `*`(x: SymNode, y: float): SymNode =
+  return x.mul newSymNode(value = SymNodeValueFloat(floatValue: y))
+
+when isMainModule:
+  let p = newSym("p")
+  let k = 0.5 * p.norm2
+  let dkdp = k.gradient p
+
+  import qex
+  import os, strutils
+
+  qexInit()
+  type GaugeActType = enum ActWilson, ActAdjoint, ActRect, ActSymanzik, ActIwasaki, ActDBW2
+  converter toGaugeActType(s:string):GaugeActType = parseEnum[GaugeActType](s)
+  letParam:
+    gaugefile = ""
+    savefile = "config"
+    savefreq = 10
+    lat =
+      if fileExists(gaugefile):
+        getFileLattice gaugefile
+      else:
+        if gaugefile.len > 0:
+          qexWarn "Nonexistent gauge file: ", gaugefile
+        @[4,4,4,8]
+    gact:GaugeActType = "ActWilson"
+    beta = 6.0
+    adjFac = -0.25
+    rectFac = -1.4088
+    seed:uint64 = 4321
+  echoParams()
+  echo "rank ", myRank, "/", nRanks
+  threads: echo "thread ", threadNum, "/", numThreads
+
+  let
+    gc = case gact
+      of ActWilson: GaugeActionCoeffs(plaq: beta)
+      of ActAdjoint: GaugeActionCoeffs(plaq: beta, adjplaq: beta*adjFac)
+      of ActRect: gaugeActRect(beta, rectFac)
+      of ActSymanzik: Symanzik(beta)
+      of ActIwasaki: Iwasaki(beta)
+      of ActDBW2: DBW2(beta)
+    lo = lat.newLayout
+    vol = lo.physVol
+
+  echo gc
+  var r = lo.newRNGField(MRG32k3a, seed)
+  var R:MRG32k3a  # global RNG
+  R.seed(seed, 987654321)
+
+  var mom = lo.newgauge
+  threads:
+    mom.randomTAH r
+
+  p.assign mom
+  optimize(output = @[k,dkdp], variables = @[p])
+
+  k.allocate
+  dkdp.allocate
+  echo k.treerepr
+  echo dkdp.treerepr
+
+  k.eval
+  dkdp.eval
+
+  echo k.value
+  let dp = dkdp.value.getGauge
+  var d2 = 0.0
+  threads:
+    var d2t = 0.0
+    for mu in 0..<mom.len:
+      mom[mu] -= dp[mu]
+      d2t += mom[mu].norm2
+    threadMaster: d2 = d2t
+  echo d2
diff --git a/src/experimental/symbolic/graph.nim b/src/experimental/symbolic/graph.nim
index dd63d5b..244f695 100644
--- a/src/experimental/symbolic/graph.nim
+++ b/src/experimental/symbolic/graph.nim
@@ -50,16 +50,16 @@ type
     sntVisited, sntAssigned, sntNeedUpdate, sntNeedGradient, sntFixedGradient
     # sntReusable, ...
   SymNodeTags = set[SymNodeTag]
-  SymNodeValue = ref object of RootObj  ## Represent unallocated symbolic value
-  SymNodeValueConcrete = ref object of SymNodeValue  ## For any concrete values
+  SymNodeValue* = ref object of RootObj  ## Represent unallocated symbolic value
+  SymNodeValueConcrete* = ref object of SymNodeValue  ## For any concrete values
   SymNodeGradient = object
     ## for a particular variable v
     dependent: SymNode  ## a variable v that depends on this node, x
     gradient: SymNode  ## dv/dx
-  SymNode = ref object
+  SymNode* = ref object
     # This object can be cyclic, because gradients refer to ancestor nodes
-    value: SymNodeValue
-    inputs: seq[SymNode]
+    value*: SymNodeValue
+    inputs*: seq[SymNode]
     forward: proc(z: SymNode)  ## runs the actual compute
     arg: SymNodeValue  ## extra argument forward/backward uses
     runCount: int
@@ -72,34 +72,34 @@ type
   SymNodeError = object of Defect
   SymNodeValueError = object of Defect
 
-template raiseError(msg: string) =
+template raiseError*(msg: string) =
   raise newException(SymNodeError, msg)
 
-template raiseValueError(msg: string) =
+template raiseValueError*(msg: string) =
   raise newException(SymNodeValueError, msg)
 
-template raiseErrorBaseMethod(msg: string) =
+template raiseErrorBaseMethod*(msg: string) =
   raise newException(
     SymNodeError,
     "Base method invoked: " & msg &
       "\nMake sure to pass `--multimethods:on` and check there is a custom method for each derived type.")
 
-method `$`(v: SymNodeValue): string {.base.} = "SymNodeValue"
+method `$`*(v: SymNodeValue): string {.base.} = "SymNodeValue"
 
-func `$`(z: SymNode): string =
+func `$`*(z: SymNode): string =
   z.name & "#" & $z.id
 
-func nodeRepr(z: SymNode): string
-func treeRepr(z: SymNode): string
+func nodeRepr*(z: SymNode): string
+func treeRepr*(z: SymNode): string
 
-method copySymNodeValue(v: SymNodeValue): SymNodeValue {.base.} =
+method copySymNodeValue*(v: SymNodeValue): SymNodeValue {.base.} =
   ## nothing to copy
   v
 
-method copySymNodeValue(v: SymNodeValueConcrete): SymNodeValueConcrete =
+method copySymNodeValue*(v: SymNodeValueConcrete): SymNodeValueConcrete =
   raiseValueError("Custom method required for concrete value: " & $v)
 
-proc newSymNode(
+proc newSymNode*(
     value = SymNodeValue(),
     inputs: seq[SymNode] = @[],
     forward: proc(z: SymNode) = nil,
@@ -116,13 +116,13 @@ proc newSymNode(
     allocateValue: allocateValue, backward: backward, gradients: gradients, name: name, tag: tag, id: id)
   id.inc
 
-proc copySymNode(z: SymNode): SymNode =
+proc copySymNode*(z: SymNode): SymNode =
   newSymNode(value = z.value.copySymNodeValue, inputs = z.inputs, forward = z.forward,
     arg = if z.arg != nil: z.arg.copySymNodeValue else: nil, runCount = z.runCount,
     allocateValue = z.allocateValue, backward = z.backward, gradients = z.gradients,
     name = z.name, tag = z.tag)
 
-proc assignSymNode(z: SymNode, x: SymNode) =
+proc assignSymNode*(z: SymNode, x: SymNode) =
   z.value = x.value.copySymNodeValue
   z.inputs = x.inputs
   z.forward = x.forward
@@ -135,7 +135,7 @@ proc assignSymNode(z: SymNode, x: SymNode) =
   z.name = x.name
   z.tag = x.tag
 
-proc gradientDependentOrNil(z: SymNode, dep: SymNode): SymNode =
+proc gradientDependentOrNil*(z: SymNode, dep: SymNode): SymNode =
   ## May return nil if not exists.
   for g in z.gradients:
     if dep == g.dependent:
@@ -151,7 +151,7 @@ proc gradientDependentAssign(z: SymNode, dep: SymNode, grad: SymNode) =
       return
   z.gradients.add SymNodeGradient(dependent: dep, gradient: grad)
 
-proc assign(z: SymNode, v: SymNodeValueConcrete) =
+proc assign*(z: SymNode, v: SymNodeValueConcrete) =
   z.value = v
   z.tag.incl sntAssigned
 
@@ -159,54 +159,54 @@ proc assign(z: SymNode, v: SymNodeValueConcrete) =
 # generic symbol support
 #
 
-proc newSym(s: string): SymNode =
+proc newSym*(s: string): SymNode =
   newSymNode(name = s)
 
-method identSymNodeValue(z: SymNodeValue, x: SymNodeValue) {.base.} =
+method identSymNodeValue*(z: SymNodeValue, x: SymNodeValue) {.base.} =
   raiseErrorBaseMethod("args:\n  " & z.repr & "\n  " & x.repr)
 
-method identAllocateSymNodeValue(z: SymNode, x: SymNodeValue) {.base.} =
+method identAllocateSymNodeValue*(z: SymNode, x: SymNodeValue) {.base.} =
   raiseErrorBaseMethod("args:\n  " & z.nodeRepr & "\n  " & x.repr)
 
-method iaddSymNodeValue(z: SymNodeValue, x: SymNodeValue, y: SymNodeValue) {.base.} =
+method iaddSymNodeValue*(z: SymNodeValue, x: SymNodeValue, y: SymNodeValue) {.base.} =
   raiseErrorBaseMethod("args:\n  " & z.repr & "\n  " & x.repr & "\n  " & y.repr)
 
-method iaddAllocateSymNodeValue(z: SymNode, x: SymNodeValue, y: SymNodeValue) {.base.} =
+method iaddAllocateSymNodeValue*(z: SymNode, x: SymNodeValue, y: SymNodeValue) {.base.} =
   raiseErrorBaseMethod("args:\n  " & z.nodeRepr & "\n  " & x.repr & "\n  " & y.repr)
 
 #
 # float support
 #
 
-type SymNodeValueFloat = ref object of SymNodeValueConcrete
-  floatValue: float
+type SymNodeValueFloat* = ref object of SymNodeValueConcrete
+  floatValue*: float
 
-method `$`(v: SymNodeValueFloat): string = $v.floatValue
+method `$`*(v: SymNodeValueFloat): string = $v.floatValue
 
-proc assign(z: SymNode, v: float) =
+proc assign*(z: SymNode, v: float) =
   z.assign SymNodeValueFloat(floatValue: v)
 
-method copySymNodeValue(v: SymNodeValueFloat): SymNodeValueFloat =
+method copySymNodeValue*(v: SymNodeValueFloat): SymNodeValueFloat =
   SymNodeValueFloat(floatValue: v.floatValue)
 
-method identSymNodeValue(z: SymNodeValueFloat, x: SymNodeValueFloat) =
+method identSymNodeValue*(z: SymNodeValueFloat, x: SymNodeValueFloat) =
   z.floatValue = x.floatValue
 
-method identAllocateSymNodeValue(z: SymNode, x: SymNodeValueFloat) =
+method identAllocateSymNodeValue*(z: SymNode, x: SymNodeValueFloat) =
   z.value = SymNodeValueFloat()
 
-method iaddSymNodeValue(z: SymNodeValueFloat, x: SymNodeValueFloat, y: SymNodeValueFloat) =
+method iaddSymNodeValue*(z: SymNodeValueFloat, x: SymNodeValueFloat, y: SymNodeValueFloat) =
   z.floatValue = x.floatValue + y.floatValue
 
-method iaddAllocateSymNodeValue(z: SymNode, x: SymNodeValueFloat, y: SymNodeValueFloat) =
+method iaddAllocateSymNodeValue*(z: SymNode, x: SymNodeValueFloat, y: SymNodeValueFloat) =
   z.value = SymNodeValueFloat()
 
 #
 # minimum algebra for the nodes
 #
 
-proc ident(x:SymNode): SymNode
-proc add(x: SymNode, y: SymNode): SymNode
+proc ident*(x:SymNode): SymNode
+proc add*(x: SymNode, y: SymNode): SymNode
 
 proc identForward(z: SymNode) =
   identSymNodeValue(z.value, z.inputs[0].value)
@@ -215,13 +215,15 @@ proc identAllocate(z: SymNode) =
   identAllocateSymNodeValue(z, z.inputs[0].value)
 
 proc identBackward(z: SymNode, i: int, dep: SymNode): SymNode =
+  if i != 0:
+    raiseError("ident has 1 operand, got i = " & $i)
   let g = z.gradientDependentOrNil dep
   if g == nil:
     return newSymNode(value = SymNodeValueFloat(floatValue: 1.0), name = "One[ident]")
   else:
     return g.ident
 
-proc ident(x:SymNode): SymNode =
+proc ident*(x:SymNode): SymNode =
   newSymNode(
     inputs = @[x],
     forward = identForward,
@@ -236,13 +238,15 @@ proc addAllocate(z: SymNode) =
   iaddAllocateSymNodeValue(z, z.inputs[0].value, z.inputs[1].value)
 
 proc addBackward(z: SymNode, i: int, dep: SymNode): SymNode =
+  if i != 0 and i != 1:
+    raiseError("add has 2 operands, got i = " & $i)
   let g = z.gradientDependentOrNil dep
   if g == nil:
     return newSymNode(value = SymNodeValueFloat(floatValue: 1.0), name = "One[add]")
   else:
     return g.ident
 
-proc add(x: SymNode, y: SymNode): SymNode =
+proc add*(x: SymNode, y: SymNode): SymNode =
   newSymNode(
     inputs = @[x, y],
     forward = addForward,
@@ -271,10 +275,14 @@ proc allocateRec(z: SymNode) =
         raiseError("undefined allocateValue for node: " & z.nodeRepr)
       z.allocateValue z
 
-proc allocate(z: SymNode) =
+proc allocate*(z: SymNode) =
   z.allocateRec
   z.tagClearRec sntVisited
 
+proc tagUpdate*(z: SymNode) =
+  ## call this for newly created inner nodes in *Backward
+  z.tag.incl sntNeedUpdate
+
 proc tagUpdateRec(z: SymNode) =
   if sntVisited in z.tag:
     return
@@ -305,7 +313,7 @@ proc evalRec(z: SymNode) =
       raiseError("inputs.len: " & $z.inputs.len & ", but no forward function defined for:\n" & z.nodeRepr)
     z.tag.excl sntNeedUpdate
 
-proc eval(z: SymNode) =
+proc eval*(z: SymNode) =
   z.tagUpdateRec
   z.tagClearRec sntVisited
   z.evalRec
@@ -361,7 +369,7 @@ proc gradientRec(z: SymNode, dep: SymNode) =
           input.tag.incl sntFixedGradient
         input.gradientRec dep
 
-proc gradient(dep: SymNode, x: SymNode): SymNode =
+proc gradient*(dep: SymNode, x: SymNode): SymNode =
   if sntNeedGradient notin x.tag:
     x.tag.incl sntNeedGradient
   dep.tagUpdateNeedGradientRec
@@ -372,6 +380,10 @@ proc gradient(dep: SymNode, x: SymNode): SymNode =
   dep.tagClearRec sntFixedGradient
   x.gradientDependentOrNil dep
 
+proc optimize*(output: seq[SymNode], variables: seq[SymNode]) =
+  # TODO
+  discard
+
 proc sharedNodesRec(shared: var seq[SymNode], z: SymNode) =
   if sntVisited in z.tag:
     var found = false
@@ -390,7 +402,7 @@ proc sharedNodesRec(shared: var seq[SymNode], z: SymNode) =
 # to string
 #
 
-func nodeRepr(z: SymNode): string =
+func nodeRepr*(z: SymNode): string =
   result = $z & $z.tag & ": " & $z.value & ", run: " & $z.runCount
   if z.arg != nil:
     result &= ", arg: " & $z.arg
@@ -424,7 +436,7 @@ func toStringRec(z: SymNode, pre: string, shared: seq[SymNode]): string =
     for i in z.inputs:
       result &= "\n" & i.toStringRec(pre & "  ", shared)
 
-func treeRepr(z: SymNode): string =
+func treeRepr*(z: SymNode): string =
   var shared = newseq[SymNode]()
   shared.sharedNodesRec z
   z.tagClearRec sntVisited