From 76fd02b0a3e9fe7fd49a9b1955a56ca25dbed2ac Mon Sep 17 00:00:00 2001
From: breandan <bre@ndan.co>
Date: Fri, 17 Jan 2020 17:19:20 -0500
Subject: [PATCH] initial support for heterogeneous input types #5

---
 .../kotlingrad/experimental/ToyExample.kt     | 232 ++++++++++--------
 .../experimental/ToyMatrixExample.kt          | 110 +++++----
 .../experimental/ToyVectorExample.kt          | 164 ++++++-------
 .../experimental/VariableCapture.kt           |  38 +--
 .../calculus/ExpressionGenerator.kt           |  18 +-
 .../kotlingrad/evaluation/TestSymbolic.kt     |   4 +-
 .../samples/MultilayerPerceptron.kt           |   8 +-
 .../kotlingrad/samples/VisualizeDFG.kt        |   4 +-
 .../samples/physics/DoublePendulum.kt         |  20 +-
 9 files changed, 314 insertions(+), 284 deletions(-)

diff --git a/core/src/main/kotlin/edu/umontreal/kotlingrad/experimental/ToyExample.kt b/core/src/main/kotlin/edu/umontreal/kotlingrad/experimental/ToyExample.kt
index f007c7af..d9d02797 100644
--- a/core/src/main/kotlin/edu/umontreal/kotlingrad/experimental/ToyExample.kt
+++ b/core/src/main/kotlin/edu/umontreal/kotlingrad/experimental/ToyExample.kt
@@ -57,63 +57,82 @@ interface Field<X : Field<X>> : Group<X> {
   fun ln(): X
 }
 
+interface Fun<X: SFun<X>> {
+  val inputs: Inputs<X>
+}
+
+/**
+ * TODO: Figure out how to merge with
+ * @see Bindings
+ */
+
+data class Inputs<X: SFun<X>>(
+  val sVars: Set<Var<X>> = emptySet(),
+  val vVars: Set<VVar<X, *>> = emptySet(),
+  val mVars: Set<MVar<X, *, *>> = emptySet()
+) {
+  constructor(inputs: List<Inputs<X>>): this(
+    inputs.flatMap { it.sVars }.toSet(),
+    inputs.flatMap { it.vVars }.toSet(),
+    inputs.flatMap { it.mVars }.toSet()
+  )
+
+  constructor(vararg funs: Fun<X>): this(funs.map { it.inputs })
+}
+
 /**
  * Scalar function.
  */
 
-sealed class Fun<X : Fun<X>>(open val sVars: Set<Var<X>> = emptySet()) : Field<Fun<X>>, (Bindings<X>) -> Fun<X> {
-  constructor(fn: Fun<X>) : this(fn.sVars)
-  constructor(vararg fns: Fun<X>) : this(fns.flatMap { it.sVars }.toSet())
+sealed class SFun<X: SFun<X>>(override val inputs: Inputs<X>): Fun<X>, Field<SFun<X>>, (Bindings<X>) -> SFun<X> {
+  constructor(vararg funs: Fun<X>) : this(Inputs(*funs))
 
-  override operator fun plus(addend: Fun<X>): Fun<X> = Sum(this, addend)
-  override operator fun times(multiplicand: Fun<X>): Fun<X> = Prod(this, multiplicand)
-  override operator fun div(divisor: Fun<X>): Fun<X> = this * divisor.pow(-One<X>())
+  override operator fun plus(addend: SFun<X>): SFun<X> = Sum(this, addend)
+  override operator fun times(multiplicand: SFun<X>): SFun<X> = Prod(this, multiplicand)
+  override operator fun div(divisor: SFun<X>): SFun<X> = this * divisor.pow(-One<X>())
   open operator fun <E : D1> times(multiplicand: VFun<X, E>): VFun<X, E> = SVProd(this, multiplicand)
   open operator fun <R : D1, C : D1> times(multiplicand: MFun<X, R, C>): MFun<X, R, C> = SMProd(this, multiplicand)
 
-  override operator fun invoke(bnds: Bindings<X>): Fun<X> =
+  override operator fun invoke(bnds: Bindings<X>): SFun<X> =
       Composition(this, bnds).run { if (bnds.isReassignmentFree) evaluate else this }
 
-  open operator fun invoke(): Fun<X> = invoke(Bindings())
+  open operator fun invoke(): SFun<X> = invoke(Bindings())
 
-  open fun d(v1: Var<X>): Fun<X> = Derivative(this, v1)
+  open fun d(v1: Var<X>): SFun<X> = Derivative(this, v1)
   open fun d(v1: Var<X>, v2: Var<X>): Vec<X, D2> = Vec(Derivative(this, v1), Derivative(this, v2))
   open fun d(v1: Var<X>, v2: Var<X>, v3: Var<X>): Vec<X, D3> = Vec(Derivative(this, v1), Derivative(this, v2), Derivative(this, v3))
-  open fun d(vararg vars: Var<X>): Map<Var<X>, Fun<X>> = vars.map { it to Derivative(this, it) }.toMap()
+  open fun d(vararg vars: Var<X>): Map<Var<X>, SFun<X>> = vars.map { it to Derivative(this, it) }.toMap()
 
   open fun <L: D1> d(vVar: VVar<X, L>): VFun<X, L> = Gradient(this, vVar)
   open fun <R: D1, C: D1> d(mVar: MVar<X, R, C>): MFun<X, R, C> = TODO()
 
-  open fun <R : D1, C : D1> d(mv: Mat<X, R, C>) = d(*mv.sVars.toTypedArray()).values.foldIndexed(mutableListOf()) { index, acc: MutableList<MutableList<Fun<X>>>, p ->
+  open fun <R : D1, C : D1> d(mv: Mat<X, R, C>) = d(*mv.inputs.sVars.toTypedArray()).values.foldIndexed(mutableListOf()) { index, acc: MutableList<MutableList<SFun<X>>>, p ->
     if (index % mv.numCols == 0) acc.add(mutableListOf(p)) else acc.last().add(p)
     acc
   }.map { Vec<X, C>(it) }.let { Mat<X, R, C>(it) }
 
-  open fun grad(): Map<Var<X>, Fun<X>> = sVars.map { it to Derivative(this, it) }.toMap()
-
-  override fun ln(): Fun<X> = Log(this)
-  override fun pow(exp: Fun<X>): Fun<X> = Power(this, exp)
-  override fun unaryMinus(): Fun<X> = Negative(this)
-  open fun sqrt(): Fun<X> = this pow (One<X>() / (Two<X>()))
-
-  override fun toString(): String = when {
-    this is Log -> "ln($logarithmand)"
-    this is Negative -> "-($value)"
-    this is Power -> "($base) pow ($exponent)"
-//    this is Prod && right is Sum -> "$left * ($right)"
-//    this is Prod && left is Sum -> "($left) * $right"
-    this is Prod -> "($left) * ($right)"
-    this is Sum && right is Negative -> "$left - ${right.value}"
-    this is Sum -> "$left + $right"
-    this is Var -> name
-    this is Derivative -> "d($fn) / d($vrb)"
-    this is Zero -> "0" //"\uD835\uDFD8" // 𝟘
-    this is One -> "1" //"\uD835\uDFD9"  // 𝟙
-    this is Two -> "2" //"\uD835\uDFDA"  // 𝟚
-    this is E -> "\u2147" // ⅇ
-    this is VMagnitude -> "|$value|"
-    this is DProd -> "($left) dot ($right)"
-    this is Composition -> "($fn) comp $bindings"
+  open fun grad(): Map<Var<X>, SFun<X>> = inputs.sVars.map { it to Derivative(this, it) }.toMap()
+
+  override fun ln(): SFun<X> = Log(this)
+  override fun pow(exp: SFun<X>): SFun<X> = Power(this, exp)
+  override fun unaryMinus(): SFun<X> = Negative(this)
+  open fun sqrt(): SFun<X> = this pow (One<X>() / Two())
+
+  override fun toString(): String = when (this) {
+    is Log -> "ln($logarithmand)"
+    is Negative -> "- ($value)"
+    is Power -> "($base) pow ($exponent)"
+    is Prod -> "($left) * ($right)"
+    is Sum -> if(right is Negative) "$left $right" else "$left + $right"
+    is Var -> name
+    is Derivative -> "d($fn) / d($vrb)"
+    is Zero -> "0" //"\uD835\uDFD8" // 𝟘
+    is One -> "1"  //"\uD835\uDFD9" // 𝟙
+    is Two -> "2"  //"\uD835\uDFDA" // 𝟚
+    is E -> "E()"  //"\u2147"       // ⅇ
+    is VMagnitude -> "|$value|"
+    is DProd -> "($left) dot ($right)"
+    is Composition -> "($fn)($bindings)"
     else -> super.toString()
   }
 
@@ -128,7 +147,7 @@ sealed class Fun<X : Fun<X>>(open val sVars: Set<Var<X>> = emptySet()) : Field<F
   }
 
   fun toGraph(): MutableNode = mutNode(if (this is Var) "$this" else "${hashCode()}").apply {
-    when (this@Fun) {
+    when (this@SFun) {
       is Var -> name
       is Negative -> { value.toGraph() - this; add(Label.of("neg")) }
       is Derivative -> { fn.toGraph() - this; mutNode("$this").apply { add(Label.of(vrb.toString())) } - this; add(Label.of("d")) }
@@ -137,7 +156,7 @@ sealed class Fun<X : Fun<X>>(open val sVars: Set<Var<X>> = emptySet()) : Field<F
       is One -> add(Label.of("one"))
       is Zero -> add(Label.of("zero"))
       is Composition -> { bindings.sMap.entries.map { entry -> mutNode(entry.hashCode().toString()).also { compNode -> entry.key.toGraph() - compNode; entry.value.toGraph() - compNode; compNode.add(Label.of("comp")) } }.map { it - this; add(Label.of("bindings")) } }
-      else -> TODO(this@Fun.javaClass.toString())
+      else -> TODO(this@SFun.javaClass.toString())
     }
   }
 }
@@ -146,16 +165,16 @@ sealed class Fun<X : Fun<X>>(open val sVars: Set<Var<X>> = emptySet()) : Field<F
  * Symbolic operators.
  */
 
-sealed class BiFun<X: Fun<X>>(val left: Fun<X>, val right: Fun<X>): Fun<X>(left, right)
+sealed class BiFun<X: SFun<X>>(val left: SFun<X>, val right: SFun<X>): SFun<X>(left, right)
 
-class Sum<X : Fun<X>>(val addend: Fun<X>, val augend: Fun<X>): BiFun<X>(addend, augend)
+class Sum<X : SFun<X>>(val addend: SFun<X>, val augend: SFun<X>): BiFun<X>(addend, augend)
 
-class Negative<X : Fun<X>>(val value: Fun<X>) : Fun<X>(value)
-class Prod<X : Fun<X>>(val multiplicand: Fun<X>, val multiplicator: Fun<X>): BiFun<X>(multiplicand, multiplicator)
-class Power<X : Fun<X>> internal constructor(val base: Fun<X>, val exponent: Fun<X>) : BiFun<X>(base, exponent)
-class Log<X : Fun<X>> internal constructor(val logarithmand: Fun<X>, val base: Fun<X> = E()) : BiFun<X>(logarithmand, base)
-class Derivative<X : Fun<X>>(val fn: Fun<X>, val vrb: Var<X>) : Fun<X>(fn, vrb) {
-  fun Fun<X>.df(): Fun<X> = when (this) {
+class Negative<X : SFun<X>>(val value: SFun<X>) : SFun<X>(value)
+class Prod<X : SFun<X>>(val multiplicand: SFun<X>, val multiplicator: SFun<X>): BiFun<X>(multiplicand, multiplicator)
+class Power<X : SFun<X>> internal constructor(val base: SFun<X>, val exponent: SFun<X>) : BiFun<X>(base, exponent)
+class Log<X : SFun<X>> internal constructor(val logarithmand: SFun<X>, val base: SFun<X> = E()) : BiFun<X>(logarithmand, base)
+class Derivative<X : SFun<X>>(val fn: SFun<X>, val vrb: Var<X>) : SFun<X>(fn, vrb) {
+  fun SFun<X>.df(): SFun<X> = when (this@df) {
     is Var -> if (this == vrb) One() else Zero()
     is SConst -> Zero()
     is Sum -> left.df() + right.df()
@@ -166,15 +185,15 @@ class Derivative<X : Fun<X>>(val fn: Fun<X>, val vrb: Var<X>) : Fun<X>(fn, vrb)
     is Derivative -> fn.df()
     is DProd -> this().df()
     is VMagnitude -> this().df()
-    is Composition -> bindings.curried().fold(One()) { acc: Fun<X>, binding ->
+    is Composition -> bindings.curried().fold(One()) { acc: SFun<X>, binding ->
       acc * fn.df()(binding) * binding.sMap.entries.first().value.df()
     }
   }
 }
 
-class Composition<X : Fun<X>>(val fn: Fun<X>, val bindings: Bindings<X>) : Fun<X>() {
-  val evaluate by lazy { apply() }
-  override val sVars: Set<Var<X>> by lazy { evaluate.sVars }
+class Composition<X : SFun<X>>(val fn: SFun<X>, val bindings: Bindings<X>) : SFun<X>(fn) {
+  val evaluate by lazy { call() }
+  override val inputs: Inputs<X> by lazy { evaluate.inputs }
 
 //  private fun calculateFixpoint(): Fun<X> {
 //    var result = apply()
@@ -187,64 +206,63 @@ class Composition<X : Fun<X>>(val fn: Fun<X>, val bindings: Bindings<X>) : Fun<X
 //    return result
 //  }
 
-  fun Fun<X>.apply(): Fun<X> =
-    bindings.sMap.getOrElse(this) {
-      when (this) {
-        is Zero -> bindings.zero
-        is One -> bindings.one
-        is Two -> bindings.two
-        is E -> bindings.e
-        is Var -> this
-        is SConst -> this
-        is Prod -> left.apply() * right.apply()
-        is Sum -> left.apply() + right.apply()
-        is Power -> base.apply() pow exponent.apply()
-        is Negative -> -value.apply()
-        is Log -> logarithmand.apply().ln()
-        is Derivative -> df().apply()
-        is DProd -> left(bindings) as Vec<X, D1> dot right(bindings) as Vec<X, D1>
-        is VMagnitude -> value(bindings).magnitude()
-        is Composition -> fn.apply().apply()
-      }
-    }
+  fun SFun<X>.call(): SFun<X> = bindings.sMap.getOrElse(this@call) { bind() }
+
+  fun SFun<X>.bind() = when (this@bind) {
+    is Zero -> bindings.zero
+    is One -> bindings.one
+    is Two -> bindings.two
+    is E -> bindings.e
+    is Var -> this
+    is SConst -> this
+    is Prod -> left.call() * right.call()
+    is Sum -> left.call() + right.call()
+    is Power -> base.call() pow exponent.call()
+    is Negative -> -value.call()
+    is Log -> logarithmand.call().ln()
+    is Derivative -> df().call()
+    is DProd -> left(bindings) as Vec<X, D1> dot right(bindings) as Vec<X, D1>
+    is VMagnitude -> value(bindings).magnitude()
+    is Composition -> fn.call().call()
+  }
 }
 
-data class Bindings<X : Fun<X>>(
-    val sMap: Map<Fun<X>, Fun<X>> = mapOf(),
-    val zero: Fun<X> = Zero(),
-    val one: Fun<X> = One(),
-    val two: Fun<X> = Two(),
-    val e: Fun<X> = E()) {
+data class Bindings<X : SFun<X>>(
+  val sMap: Map<SFun<X>, SFun<X>> = mapOf(),
+  val zero: SFun<X> = Zero(),
+  val one: SFun<X> = One(),
+  val two: SFun<X> = Two(),
+  val e: SFun<X> = E()) {
   //  constructor(sMap: Map<Fun<X>, Fun<X>>,
 //              vMap: Map<VFun<X, *>, VFun<X, *>>,
 //              zero: Fun<X>,
 //              one: Fun<X>,
 //              two: Fun<X>,
 //              E: Fun<X>): this(sMap, zero, one, two, E)
-  val isReassignmentFree = sMap.values.all { it.sVars.isEmpty() }
-  fun determines(fn: Fun<X>) = fn.sVars.all { it in sMap }
+  val isReassignmentFree = sMap.values.all { it.inputs.sVars.isEmpty() }
+  fun fullyDetermines(fn: SFun<X>) = fn.inputs.sVars.all { it in sMap }
   override fun toString() = sMap.toString()
   operator fun contains(v: Var<X>) = v in sMap
   fun curried() = sMap.entries.map { Bindings(mapOf(it.key to it.value), zero, one, two, e) }
 }
 
-class DProd<X: Fun<X>>(val left: VFun<X, *>, val right: VFun<X, *>): Fun<X>(left.sVars + right.sVars)//, left.vVars + right.vVars)
+class DProd<X: SFun<X>>(val left: VFun<X, *>, val right: VFun<X, *>): SFun<X>(left, right)
 
-class VMagnitude<X: Fun<X>>(val value: VFun<X, *>): Fun<X>(value.sVars)//, value.vVars)
+class VMagnitude<X: SFun<X>>(val value: VFun<X, *>): SFun<X>(value)//, value.vVars)
 
 interface Variable { val name: String }
 
-class Var<X : Fun<X>>(override val name: String = "") : Variable, Fun<X>() {
-  override val sVars: Set<Var<X>> = setOf(this)
+class Var<X : SFun<X>>(override val name: String = "") : Variable, SFun<X>() {
+  override val inputs: Inputs<X> = Inputs(setOf(this))
 }
 
-open class SConst<X : Fun<X>> : Fun<X>()
-class Zero<X: Fun<X>> : SConst<X>()
-class One<X: Fun<X>> : SConst<X>()
-class Two<X: Fun<X>> : SConst<X>()
-class E<X: Fun<X>> : SConst<X>()
+open class SConst<X : SFun<X>> : SFun<X>()
+class Zero<X: SFun<X>> : SConst<X>()
+class One<X: SFun<X>> : SConst<X>()
+class Two<X: SFun<X>> : SConst<X>()
+class E<X: SFun<X>> : SConst<X>()
 
-abstract class RealNumber<X : Fun<X>>(open val value: Number) : SConst<X>()
+abstract class RealNumber<X : SFun<X>>(open val value: Number) : SConst<X>()
 
 class DReal(override val value: Double) : RealNumber<DReal>(value) {
   override fun unaryMinus() = DReal(-value)
@@ -255,17 +273,17 @@ class DReal(override val value: Double) : RealNumber<DReal>(value) {
    * Constant propagation.
    */
 
-  override fun plus(addend: Fun<DReal>) = when (addend) {
+  override fun plus(addend: SFun<DReal>) = when (addend) {
     is DReal -> DReal(value + addend.value)
     else -> super.plus(addend)
   }
 
-  override fun times(multiplicand: Fun<DReal>) = when (multiplicand) {
+  override fun times(multiplicand: SFun<DReal>) = when (multiplicand) {
     is DReal -> DReal(value * multiplicand.value)
     else -> super.times(multiplicand)
   }
 
-  override fun pow(exp: Fun<DReal>) = when (exp) {
+  override fun pow(exp: SFun<DReal>) = when (exp) {
     is DReal -> DReal(value.pow(exp.value))
     else -> super.pow(exp)
   }
@@ -290,27 +308,27 @@ class DReal(override val value: Double) : RealNumber<DReal>(value) {
  */
 
 sealed class Protocol<X : RealNumber<X>> {
-  class IndVar<X: Fun<X>> constructor(val fn: Fun<X>)
+  class IndVar<X: SFun<X>> constructor(val fn: SFun<X>)
 
-  class Differential<X: Fun<X>>(private val fx: Fun<X>) {
-    // TODO: make sure this notation works for arbitrary nested functions using the Chain rule
-    infix operator fun div(arg: Differential<X>) = fx.d(arg.fx.sVars.first())
+  class Differential<X: SFun<X>>(private val fx: SFun<X>) {
+    // TODO: ensure correctness for arbitrary nested functions using the Chain rule
+    infix operator fun div(arg: Differential<X>) = fx.d(arg.fx.inputs.sVars.first())
   }
 
-  fun <X: Fun<X>> d(fn: Fun<X>) = Differential(fn)
+  fun <X: SFun<X>> d(fn: SFun<X>) = Differential(fn)
   abstract fun wrap(default: Number): X
 
-  operator fun Number.times(multiplicand: Fun<X>) = multiplicand * wrap(this)
-  operator fun Fun<X>.times(multiplicand: Number) = wrap(multiplicand) * this
+  operator fun Number.times(multiplicand: SFun<X>) = multiplicand * wrap(this)
+  operator fun SFun<X>.times(multiplicand: Number) = wrap(multiplicand) * this
 
-  operator fun Number.plus(addend: Fun<X>) = addend + wrap(this)
-  operator fun Fun<X>.plus(addend: Number) = wrap(addend) + this
+  operator fun Number.plus(addend: SFun<X>) = addend + wrap(this)
+  operator fun SFun<X>.plus(addend: Number) = wrap(addend) + this
 
-  operator fun Number.minus(subtrahend: Fun<X>) = -subtrahend + wrap(this)
-  operator fun Fun<X>.minus(subtrahend: Number) = -wrap(subtrahend) + this
+  operator fun Number.minus(subtrahend: SFun<X>) = -subtrahend + wrap(this)
+  operator fun SFun<X>.minus(subtrahend: Number) = -wrap(subtrahend) + this
 
-  fun Number.pow(exp: Fun<X>) = wrap(this) pow exp
-  infix fun Fun<X>.pow(exp: Number) = this pow wrap(exp)
+  fun Number.pow(exp: SFun<X>) = wrap(this) pow exp
+  infix fun SFun<X>.pow(exp: Number) = this pow wrap(exp)
 }
 
 object DoublePrecision : Protocol<DReal>() {
@@ -323,9 +341,9 @@ object DoublePrecision : Protocol<DReal>() {
 
   fun vrb(name: String) = Var<DReal>(name)
 
-  @JvmName("ValBnd") operator fun Fun<DReal>.invoke(vararg pairs: Pair<Var<DReal>, Number>) =
+  @JvmName("ValBnd") operator fun SFun<DReal>.invoke(vararg pairs: Pair<Var<DReal>, Number>) =
     this(Bindings(pairs.map { (it.first to wrap(it.second)) }.toMap(), zero, one, two, e))
-  @JvmName("FunBnd") operator fun Fun<DReal>.invoke(vararg pairs: Pair<Fun<DReal>, Fun<DReal>>) =
+  @JvmName("FunBnd") operator fun SFun<DReal>.invoke(vararg pairs: Pair<SFun<DReal>, SFun<DReal>>) =
     this(Bindings(pairs.map { (it.first to it.second) }.toMap(), zero, one, two, e))
 
   operator fun <Y : D1> VFun<DReal, Y>.invoke(vararg sPairs: Pair<Var<DReal>, Number>) =
@@ -334,7 +352,7 @@ object DoublePrecision : Protocol<DReal>() {
   operator fun <Rows : D1, Cols: D1> MFun<DReal, Rows, Cols>.invoke(vararg sPairs: Pair<Var<DReal>, Number>) =
     this(Bindings(sPairs.map { (it.first to wrap(it.second)) }.toMap(), zero, one, two, e))
 
-  fun Fun<DReal>.asDouble() = (this as DReal).value
+  fun SFun<DReal>.asDouble() = (this as DReal).value
 
   val x = vrb("X")
   val y = vrb("Y")
diff --git a/core/src/main/kotlin/edu/umontreal/kotlingrad/experimental/ToyMatrixExample.kt b/core/src/main/kotlin/edu/umontreal/kotlingrad/experimental/ToyMatrixExample.kt
index 16acd951..7897af6b 100644
--- a/core/src/main/kotlin/edu/umontreal/kotlingrad/experimental/ToyMatrixExample.kt
+++ b/core/src/main/kotlin/edu/umontreal/kotlingrad/experimental/ToyMatrixExample.kt
@@ -31,7 +31,6 @@ fun main() {
       y * y,
       x * y)
 
-
     val mf2 = Mat1x2(vf2)
 
     val qr = mf2 * Vec(x, y)
@@ -57,11 +56,8 @@ fun main() {
  * Matrix function.
  */
 
-open class MFun<X: Fun<X>, R: D1, C: D1>(
-  open val sVars: Set<Var<X>> = emptySet()
-): (Bindings<X>) -> MFun<X, R, C> {
-  constructor(left: MFun<X, R, *>, right: MFun<X, *, C>): this(left.sVars + right.sVars)
-  constructor(mFun: MFun<X, R, C>): this(mFun.sVars)
+open class MFun<X: SFun<X>, R: D1, C: D1>(override val inputs: Inputs<X>): Fun<X>, (Bindings<X>) -> MFun<X, R, C> {
+  constructor(vararg funs: Fun<X>): this(Inputs(*funs))
 
   open val ᵀ: MFun<X, C, R> by lazy { MTranspose(this) }
 
@@ -76,7 +72,7 @@ open class MFun<X: Fun<X>, R: D1, C: D1>(
       is SMProd -> left(bnds) * right(bnds)
       is MConst -> MZero()
       is Mat -> Mat(rows.map { it(bnds) as Vec<X, C> })
-      else -> throw IllegalArgumentException("Type ${this::class.java.name} unknown")
+      else -> TODO(this::class.java.name)
     }
 
   // Materializes the concrete matrix from the dataflow graph
@@ -84,7 +80,7 @@ open class MFun<X: Fun<X>, R: D1, C: D1>(
 
   open operator fun unaryMinus(): MFun<X, R, C> = MNegative(this)
   open operator fun plus(addend: MFun<X, R, C>): MFun<X, R, C> = MSum(this, addend)
-  open operator fun times(multiplicand: Fun<X>): MFun<X, R, C> = MSProd(this, multiplicand)
+  open operator fun times(multiplicand: SFun<X>): MFun<X, R, C> = MSProd(this, multiplicand)
   open operator fun times(multiplicand: VFun<X, C>): VFun<X, R> = MVProd(this, multiplicand)
 
   // The Hadamard product
@@ -100,24 +96,24 @@ open class MFun<X: Fun<X>, R: D1, C: D1>(
     is HProd -> "$left ʘ $right"
     is MSProd -> "$left * $right"
     is SMProd -> "$left * $right"
-    is MConst -> "TODO()"
+    is MConst -> "${javaClass.name}()"
     is Mat -> "Mat${numRows}x$numCols(${rows.joinToString(", ") { it.contents.joinToString(", ") }})"
     is MDerivative -> "d($mFun) / d($v1)"
-    else -> throw IllegalArgumentException("Type ${this::class.java.name} unknown")
+    else -> TODO(this::class.java.name)
   }
 }
 
-class MNegative<X: Fun<X>, R: D1, C: D1>(val value: MFun<X, R, C>): MFun<X, R, C>(value)
-class MTranspose<X: Fun<X>, R: D1, C: D1>(val value: MFun<X, R, C>): MFun<X, C, R>(value.sVars)
-class MSum<X: Fun<X>, R: D1, C: D1>(val left: MFun<X, R, C>, val right: MFun<X, R, C>): MFun<X, R, C>(left, right)
-class MMProd<X: Fun<X>, R: D1, C1: D1, C2: D1>(val left: MFun<X, R, C1>, val right: MFun<X, C1, C2>): MFun<X, R, C2>(left, right)
-class HProd<X: Fun<X>, R: D1, C: D1>(val left: MFun<X, R, C>, val right: MFun<X, R, C>): MFun<X, R, C>(left, right)
-class MSProd<X: Fun<X>, R: D1, C: D1>(val left: MFun<X, R, C>, val right: Fun<X>): MFun<X, R, C>(left)
-class SMProd<X: Fun<X>, R: D1, C: D1>(val left: Fun<X>, val right: MFun<X, R, C>): MFun<X, R, C>(right)
+class MNegative<X: SFun<X>, R: D1, C: D1>(val value: MFun<X, R, C>): MFun<X, R, C>(value)
+class MTranspose<X: SFun<X>, R: D1, C: D1>(val value: MFun<X, R, C>): MFun<X, C, R>(value)
+class MSum<X: SFun<X>, R: D1, C: D1>(val left: MFun<X, R, C>, val right: MFun<X, R, C>): MFun<X, R, C>(left, right)
+class MMProd<X: SFun<X>, R: D1, C1: D1, C2: D1>(val left: MFun<X, R, C1>, val right: MFun<X, C1, C2>): MFun<X, R, C2>(left, right)
+class HProd<X: SFun<X>, R: D1, C: D1>(val left: MFun<X, R, C>, val right: MFun<X, R, C>): MFun<X, R, C>(left, right)
+class MSProd<X: SFun<X>, R: D1, C: D1>(val left: MFun<X, R, C>, val right: SFun<X>): MFun<X, R, C>(left)
+class SMProd<X: SFun<X>, R: D1, C: D1>(val left: SFun<X>, val right: MFun<X, R, C>): MFun<X, R, C>(right)
 
 // TODO: Generalize tensor derivatives? https://en.wikipedia.org/wiki/Tensor_derivative_(continuum_mechanics)
-class MDerivative<X: Fun<X>, R: D1, C: D1> internal constructor(val mFun: VFun<X, R>, numCols: Nat<C>, val v1: Var<X>): MFun<X, R, C>(mFun.sVars) {
-  fun MFun<X, R, C>.df(): MFun<X, R, C> = when (this) {
+class MDerivative<X: SFun<X>, R: D1, C: D1> internal constructor(val mFun: VFun<X, R>, numCols: Nat<C>, val v1: Var<X>): MFun<X, R, C>(mFun) {
+  fun MFun<X, R, C>.df(): MFun<X, R, C> = when (this@df) {
     is MConst -> MZero()
     is MVar -> MZero()
     is MNegative -> -value.df()
@@ -129,22 +125,38 @@ class MDerivative<X: Fun<X>, R: D1, C: D1> internal constructor(val mFun: VFun<X
     is SMProd -> left.d(v1) * right + left * right.df()
     is HProd -> left.df() ʘ right + left ʘ right.df()
     is Mat -> Mat(rows.map { it.d(v1)() })
-    else -> throw IllegalArgumentException("Unable to differentiate expression of type ${this::class.java.name}")
+    else -> TODO(this@df::class.java.name)
+  }
+}
+
+class MGradient<X : SFun<X>, R: D1, C: D1>(val fn: SFun<X>, val mVar: MVar<X, R, C>): MFun<X, R, C>(fn) {
+  fun df() = fn.df()
+  fun SFun<X>.df(): MFun<X, R, C> = when (this@df) {
+    is MVar<*, *, *> -> if (this == mVar) MOne() else MZero()
+    is Var -> MZero()
+    is SConst -> MZero()
+    is Sum -> left.df() + right.df()
+    is Prod -> left.df() * right + left * right.df()
+    is Power -> this * (exponent * Log(base)).df()
+    is Negative -> -value.df()
+    is Log -> (logarithmand pow -One<X>()) * logarithmand.df()
+//    is Derivative -> fn.df()
+    is DProd -> this().df()
+    is VMagnitude -> this().df()
+    else -> TODO(this@df::class.java.name)
   }
 }
 
-class MVar<X: Fun<X>, R: D1, C: D1>(override val name: String = ""): Variable, MFun<X, R, C>()
-open class MConst<X: Fun<X>, R: D1, C: D1>: MFun<X, R, C>()
+class MVar<X: SFun<X>, R: D1, C: D1>(override val name: String = ""): Variable, MFun<X, R, C>()
+open class MConst<X: SFun<X>, R: D1, C: D1>: MFun<X, R, C>()
 
-class MZero<X: Fun<X>, R: D1, C: D1>: MConst<X, R, C>()
-class MOne<X: Fun<X>, R: D1, C: D1>: MConst<X, R, C>()
+class MZero<X: SFun<X>, R: D1, C: D1>: MConst<X, R, C>()
+class MOne<X: SFun<X>, R: D1, C: D1>: MConst<X, R, C>()
 
-open class Mat<X: Fun<X>, R: D1, C: D1>(override val sVars: Set<Var<X>> = emptySet(),
-                                        val rows: List<Vec<X, C>>): MFun<X, R, C>() {
-  constructor(rows: List<Vec<X, C>>): this(rows.flatMap { it.sVars }.toSet(), rows)
-  constructor(vararg rows: Vec<X, C>): this(rows.flatMap { it.sVars }.toSet(), rows.asList())
+open class Mat<X: SFun<X>, R: D1, C: D1>(val rows: List<Vec<X, C>>): MFun<X, R, C>(*rows.toTypedArray()) {
+  constructor(vararg rows: Vec<X, C>): this(rows.asList())
 
-  val flatContents: List<Fun<X>> by lazy { rows.flatMap { it.contents } }
+  val flatContents: List<SFun<X>> by lazy { rows.flatMap { it.contents } }
 
   val indices = rows.indices
   val cols by lazy { indices.map { i -> Vec<X, R>(rows.map { it[i] }) } }
@@ -169,7 +181,7 @@ open class Mat<X: Fun<X>, R: D1, C: D1>(override val sVars: Set<Var<X>> = emptyS
 
   operator fun get(i: Int): VFun<X, C> = rows[i]
 
-  override operator fun times(multiplicand: Fun<X>): Mat<X, R, C> = Mat(rows.map { it * multiplicand })
+  override operator fun times(multiplicand: SFun<X>): Mat<X, R, C> = Mat(rows.map { it * multiplicand })
 
   override operator fun times(multiplicand: VFun<X, C>): VFun<X, R> =
     when (multiplicand) {
@@ -188,22 +200,22 @@ open class Mat<X: Fun<X>, R: D1, C: D1>(override val sVars: Set<Var<X>> = emptyS
     }
 }
 
-fun <X: Fun<X>> Mat1x1(v0: Vec<X, D1>): Mat<X, D1, D1> = Mat(v0)
-fun <X: Fun<X>> Mat2x1(v0: Vec<X, D1>, v1: Vec<X, D1>): Mat<X, D2, D1> = Mat(v0, v1)
-fun <X: Fun<X>> Mat3x1(v0: Vec<X, D1>, v1: Vec<X, D1>, v2: Vec<X, D1>): Mat<X, D3, D1> = Mat(v0, v1, v2)
-fun <X: Fun<X>> Mat1x2(v0: Vec<X, D2>): Mat<X, D1, D2> = Mat(v0)
-fun <X: Fun<X>> Mat2x2(v0: Vec<X, D2>, v1: Vec<X, D2>): Mat<X, D2, D2> = Mat(v0, v1)
-fun <X: Fun<X>> Mat3x2(v0: Vec<X, D2>, v1: Vec<X, D2>, v2: Vec<X, D2>): Mat<X, D3, D2> = Mat(v0, v1, v2)
-fun <X: Fun<X>> Mat1x3(v0: Vec<X, D3>): Mat<X, D1, D3> = Mat(v0)
-fun <X: Fun<X>> Mat2x3(v0: Vec<X, D3>, v1: Vec<X, D3>): Mat<X, D2, D3> = Mat(v0, v1)
-fun <X: Fun<X>> Mat3x3(v0: Vec<X, D3>, v1: Vec<X, D3>, v2: Vec<X, D3>): Mat<X, D3, D3> = Mat(v0, v1, v2)
-
-fun <X: Fun<X>> Mat1x1(d0: Fun<X>): Mat<X, D1, D1> = Mat(Vec(d0))
-fun <X: Fun<X>> Mat1x2(d0: Fun<X>, d1: Fun<X>): Mat<X, D1, D2> = Mat(Vec(d0, d1))
-fun <X: Fun<X>> Mat1x3(d0: Fun<X>, d1: Fun<X>, d2: Fun<X>): Mat<X, D1, D3> = Mat(Vec(d0, d1, d2))
-fun <X: Fun<X>> Mat2x1(d0: Fun<X>, d1: Fun<X>): Mat<X, D2, D1> = Mat(Vec(d0), Vec(d1))
-fun <X: Fun<X>> Mat2x2(d0: Fun<X>, d1: Fun<X>, d2: Fun<X>, d3: Fun<X>): Mat<X, D2, D2> = Mat(Vec(d0, d1), Vec(d2, d3))
-fun <X: Fun<X>> Mat2x3(d0: Fun<X>, d1: Fun<X>, d2: Fun<X>, d3: Fun<X>, d4: Fun<X>, d5: Fun<X>): Mat<X, D2, D3> = Mat(Vec(d0, d1, d2), Vec(d3, d4, d5))
-fun <X: Fun<X>> Mat3x1(d0: Fun<X>, d1: Fun<X>, d2: Fun<X>): Mat<X, D3, D1> = Mat(Vec(d0), Vec(d1), Vec(d2))
-fun <X: Fun<X>> Mat3x2(d0: Fun<X>, d1: Fun<X>, d2: Fun<X>, d3: Fun<X>, d4: Fun<X>, d5: Fun<X>): Mat<X, D3, D2> = Mat(Vec(d0, d1), Vec(d2, d3), Vec(d4, d5))
-fun <X: Fun<X>> Mat3x3(d0: Fun<X>, d1: Fun<X>, d2: Fun<X>, d3: Fun<X>, d4: Fun<X>, d5: Fun<X>, d6: Fun<X>, d7: Fun<X>, d8: Fun<X>): Mat<X, D3, D3> = Mat(Vec(d0, d1, d2), Vec(d3, d4, d5), Vec(d6, d7, d8))
\ No newline at end of file
+fun <X: SFun<X>> Mat1x1(v0: Vec<X, D1>): Mat<X, D1, D1> = Mat(v0)
+fun <X: SFun<X>> Mat2x1(v0: Vec<X, D1>, v1: Vec<X, D1>): Mat<X, D2, D1> = Mat(v0, v1)
+fun <X: SFun<X>> Mat3x1(v0: Vec<X, D1>, v1: Vec<X, D1>, v2: Vec<X, D1>): Mat<X, D3, D1> = Mat(v0, v1, v2)
+fun <X: SFun<X>> Mat1x2(v0: Vec<X, D2>): Mat<X, D1, D2> = Mat(v0)
+fun <X: SFun<X>> Mat2x2(v0: Vec<X, D2>, v1: Vec<X, D2>): Mat<X, D2, D2> = Mat(v0, v1)
+fun <X: SFun<X>> Mat3x2(v0: Vec<X, D2>, v1: Vec<X, D2>, v2: Vec<X, D2>): Mat<X, D3, D2> = Mat(v0, v1, v2)
+fun <X: SFun<X>> Mat1x3(v0: Vec<X, D3>): Mat<X, D1, D3> = Mat(v0)
+fun <X: SFun<X>> Mat2x3(v0: Vec<X, D3>, v1: Vec<X, D3>): Mat<X, D2, D3> = Mat(v0, v1)
+fun <X: SFun<X>> Mat3x3(v0: Vec<X, D3>, v1: Vec<X, D3>, v2: Vec<X, D3>): Mat<X, D3, D3> = Mat(v0, v1, v2)
+
+fun <X: SFun<X>> Mat1x1(d0: SFun<X>): Mat<X, D1, D1> = Mat(Vec(d0))
+fun <X: SFun<X>> Mat1x2(d0: SFun<X>, d1: SFun<X>): Mat<X, D1, D2> = Mat(Vec(d0, d1))
+fun <X: SFun<X>> Mat1x3(d0: SFun<X>, d1: SFun<X>, d2: SFun<X>): Mat<X, D1, D3> = Mat(Vec(d0, d1, d2))
+fun <X: SFun<X>> Mat2x1(d0: SFun<X>, d1: SFun<X>): Mat<X, D2, D1> = Mat(Vec(d0), Vec(d1))
+fun <X: SFun<X>> Mat2x2(d0: SFun<X>, d1: SFun<X>, d2: SFun<X>, d3: SFun<X>): Mat<X, D2, D2> = Mat(Vec(d0, d1), Vec(d2, d3))
+fun <X: SFun<X>> Mat2x3(d0: SFun<X>, d1: SFun<X>, d2: SFun<X>, d3: SFun<X>, d4: SFun<X>, d5: SFun<X>): Mat<X, D2, D3> = Mat(Vec(d0, d1, d2), Vec(d3, d4, d5))
+fun <X: SFun<X>> Mat3x1(d0: SFun<X>, d1: SFun<X>, d2: SFun<X>): Mat<X, D3, D1> = Mat(Vec(d0), Vec(d1), Vec(d2))
+fun <X: SFun<X>> Mat3x2(d0: SFun<X>, d1: SFun<X>, d2: SFun<X>, d3: SFun<X>, d4: SFun<X>, d5: SFun<X>): Mat<X, D3, D2> = Mat(Vec(d0, d1), Vec(d2, d3), Vec(d4, d5))
+fun <X: SFun<X>> Mat3x3(d0: SFun<X>, d1: SFun<X>, d2: SFun<X>, d3: SFun<X>, d4: SFun<X>, d5: SFun<X>, d6: SFun<X>, d7: SFun<X>, d8: SFun<X>): Mat<X, D3, D3> = Mat(Vec(d0, d1, d2), Vec(d3, d4, d5), Vec(d6, d7, d8))
\ No newline at end of file
diff --git a/core/src/main/kotlin/edu/umontreal/kotlingrad/experimental/ToyVectorExample.kt b/core/src/main/kotlin/edu/umontreal/kotlingrad/experimental/ToyVectorExample.kt
index cf1c78cf..90d42528 100644
--- a/core/src/main/kotlin/edu/umontreal/kotlingrad/experimental/ToyVectorExample.kt
+++ b/core/src/main/kotlin/edu/umontreal/kotlingrad/experimental/ToyVectorExample.kt
@@ -39,9 +39,10 @@ fun main() {
  * Vector function.
  */
 
-sealed class VFun<X: Fun<X>, E: D1>(
-  override val sVars: Set<Var<X>> = emptySet()): MFun<X, E, D1>(sVars) {
-  constructor(vararg vFns: VFun<X, E>): this(vFns.flatMap { it.sVars }.toSet())
+
+
+sealed class VFun<X: SFun<X>, E: D1>(override val inputs: Inputs<X>): Fun<X>, MFun<X, E, D1>(inputs) {
+  constructor(vararg funs: Fun<X>): this(Inputs(*funs))
 
   override operator fun invoke(bnds: Bindings<X>): VFun<X, E> =
     when (this) {
@@ -52,16 +53,17 @@ sealed class VFun<X: Fun<X>, E: D1>(
       is SVProd<X, E> -> left(bnds) * right(bnds)
       is VSProd<X, E> -> left(bnds) * right(bnds)
       is VDerivative -> df()(bnds)
+      is Gradient -> df()(bnds)
       is MVProd<X, *, *> -> left(bnds) as Mat<X, E, E> * (right as Vec<X, E>)(bnds)
       is VMProd<X, *, *> -> (left as Vec<X, E>)(bnds) * (right as Mat<X, E, E>)(bnds)
       is VMap<X, E> -> value(bnds).map(ef)
-      else -> throw IllegalArgumentException("Type ${this::class.java.name} unknown")
+      else -> TODO(this::class.java.name)
     }
 
   // Materializes the concrete vector from the dataflow graph
   operator fun invoke(): Vec<X, E> = invoke(Bindings()) as Vec<X, E>
 
-  open fun map(ef: (Fun<X>) -> Fun<X>): VFun<X, E> = VMap(this, ef)
+  open fun map(ef: (SFun<X>) -> SFun<X>): VFun<X, E> = VMap(this, ef)
 
   fun d(v1: Var<X>) = VDerivative(this, v1)
   fun d(v1: Var<X>, v2: Var<X>) = Jacobian<X, E, D2>(this, v1, v2)
@@ -74,70 +76,71 @@ sealed class VFun<X: Fun<X>, E: D1>(
   fun d(v1: Var<X>, v2: Var<X>, v3: Var<X>, v4: Var<X>, v5: Var<X>, v6: Var<X>, v7: Var<X>, v8: Var<X>, v9: Var<X>) = Jacobian<X, E, D9>(this, v1, v2, v3, v4, v5, v6, v7, v8, v9)
   //...
   fun d(vararg vars: Var<X>): Map<Var<X>, VFun<X, E>> = vars.map { it to VDerivative(this, it) }.toMap()
-  fun grad(): Map<Var<X>, VFun<X, E>> = sVars.map { it to VDerivative(this, it) }.toMap()
+  fun grad(): Map<Var<X>, VFun<X, E>> = inputs.sVars.map { it to VDerivative(this, it) }.toMap()
 
   override operator fun unaryMinus(): VFun<X, E> = VNegative(this)
   open operator fun plus(addend: VFun<X, E>): VFun<X, E> = VSum(this, addend)
   open operator fun minus(subtrahend: VFun<X, E>): VFun<X, E> = VSum(this, -subtrahend)
 
   open infix fun ʘ(multiplicand: VFun<X, E>): VFun<X, E> = VVProd(this, multiplicand)
-  override operator fun times(multiplicand: Fun<X>): VFun<X, E> = VSProd(this, multiplicand)
+  override operator fun times(multiplicand: SFun<X>): VFun<X, E> = VSProd(this, multiplicand)
   open operator fun <Q: D1> times(multiplicand: MFun<X, Q, E>): VFun<X, E> = VMProd(this, multiplicand)//(expand * multiplicand).rows.first()
-  open infix fun dot(multiplicand: VFun<X, E>): Fun<X> = DProd(this, multiplicand)
-
-  open fun magnitude(): Fun<X> = VMagnitude(this)
-
-  override fun toString() =
-    when (this) {
-      is Vec -> contents.joinToString(", ", "[", "]")
-      is VSum -> "$left + $right"
-      is VVProd -> "$left ʘ $right"
-      is SVProd -> "$left * $right"
-      is VSProd -> "$left * $right"
-      is VNegative -> "-($value)"
-      is VDerivative -> "d($vFun) / d($v1)"//d(${v1.joinToString(", ")})"
-      is MVProd<X, E, *> -> "$left * $right"
-      is VMProd<X, *, E> -> "$left * $right"
-      is Gradient -> "($fn).d(${vrbs.joinToString(", ")})"
-      is VMap -> "$value.map { $ef }"
-      is VVar -> "VVar($name)"
-    }
+  open infix fun dot(multiplicand: VFun<X, E>): SFun<X> = DProd(this, multiplicand)
+
+  open fun magnitude(): SFun<X> = VMagnitude(this)
+
+  override fun toString() = when (this) {
+    is Vec -> contents.joinToString(", ", "[", "]")
+    is VSum -> "$left + $right"
+    is VVProd -> "$left ʘ $right"
+    is SVProd -> "$left * $right"
+    is VSProd -> "$left * $right"
+    is VNegative -> "-($value)"
+    is VDerivative -> "d($vFun) / d($v1)"//d(${v1.joinToString(", ")})"
+    is MVProd<X, E, *> -> "$left * $right"
+    is VMProd<X, *, E> -> "$left * $right"
+    is Gradient -> "($fn).d($vVar)"
+    is VMap -> "$value.map { $ef }"
+    is VVar -> "VVar($name)"
+  }
 }
 
-class VNegative<X: Fun<X>, E: D1>(val value: VFun<X, E>): VFun<X, E>(value)
-class VMap<X: Fun<X>, E: D1>(val value: VFun<X, E>, val ef: (Fun<X>) -> Fun<X>): VFun<X, E>(value)
-
-class VSum<X: Fun<X>, E: D1>(val left: VFun<X, E>, val right: VFun<X, E>): VFun<X, E>(left, right)
-
-class VVProd<X: Fun<X>, E: D1>(val left: VFun<X, E>, val right: VFun<X, E>): VFun<X, E>(left, right)
-class SVProd<X: Fun<X>, E: D1>(val left: Fun<X>, val right: VFun<X, E>): VFun<X, E>(right.sVars + right.sVars)
-class VSProd<X: Fun<X>, E: D1>(val left: VFun<X, E>, val right: Fun<X>): VFun<X, E>(left.sVars + right.sVars)
-class MVProd<X: Fun<X>, R: D1, C: D1>(val left: MFun<X, R, C>, val right: VFun<X, C>): VFun<X, R>(left.sVars + right.sVars)
-class VMProd<X: Fun<X>, R: D1, C: D1>(val left: VFun<X, C>, val right: MFun<X, R, C>): VFun<X, C>(left.sVars + right.sVars)
-
-class Gradient<X : Fun<X>, E: D1>: VFun<X, E> {
-
-  val fn: Fun<X>
-  val vrbs: Array<out Var<X>>
-
-  constructor(fn: Fun<X>, vararg vrbs: Var<X>): super(fn.sVars) {
-    this.fn = fn
-    this.vrbs = vrbs
+class VNegative<X: SFun<X>, E: D1>(val value: VFun<X, E>): VFun<X, E>(value)
+class VMap<X: SFun<X>, E: D1>(val value: VFun<X, E>, val ef: (SFun<X>) -> SFun<X>): VFun<X, E>(value)
+
+class VSum<X: SFun<X>, E: D1>(val left: VFun<X, E>, val right: VFun<X, E>): VFun<X, E>(left, right)
+
+class VVProd<X: SFun<X>, E: D1>(val left: VFun<X, E>, val right: VFun<X, E>): VFun<X, E>(left, right)
+class SVProd<X: SFun<X>, E: D1>(val left: SFun<X>, val right: VFun<X, E>): VFun<X, E>(left, right)
+class VSProd<X: SFun<X>, E: D1>(val left: VFun<X, E>, val right: SFun<X>): VFun<X, E>(left, right)
+class MVProd<X: SFun<X>, R: D1, C: D1>(val left: MFun<X, R, C>, val right: VFun<X, C>): VFun<X, R>(left, right)
+class VMProd<X: SFun<X>, R: D1, C: D1>(val left: VFun<X, C>, val right: MFun<X, R, C>): VFun<X, C>(left, right)
+
+class Gradient<X : SFun<X>, E: D1>(val fn: SFun<X>, val vVar: VVar<X, E>): VFun<X, E>(fn) {
+  fun df() = fn.df()
+  fun SFun<X>.df(): VFun<X, E> = when (this@df) {
+    is VVar<*, *> -> if (this == vVar) VOne() else VZero()
+    is Var -> VZero()
+    is SConst -> VZero()
+    is Sum -> left.df() + right.df()
+    is Prod -> left.df() * right + left * right.df()
+    is Power -> this * (exponent * Log(base)).df()
+    is Negative -> -value.df()
+    is Log -> (logarithmand pow -One<X>()) * logarithmand.df()
+//    is Derivative -> fn.df()
+    is DProd -> this().df()
+    is VMagnitude -> this().df()
+    else -> TODO(this@df::class.java.name)
   }
-
-  //  constructor(fn: Fun<X>, vVar: VVar<X, E>): super(fn.sVars, vVar)
-override fun invoke(bnds: Bindings<X>) = Vec<X, E>(vrbs.map { Derivative(fn, it)() })(bnds)
-  //  override fun invoke(bnds: Bindings<X>) =
-  //    if(Vec<X, E>(vrbs.map { Derivative(fn, it)() })(bnds)
 }
 
-class VVar<X: Fun<X>, E: D1>(override val name: String = ""): Variable, VFun<X, E>()
-class Jacobian<X : Fun<X>, R: D1, C: D1>(val vfn: VFun<X, R>, vararg val vrbs: Var<X>): MFun<X, R, C>(vfn.sVars) {
+class VVar<X: SFun<X>, E: D1>(override val name: String = ""): Variable, VFun<X, E>()
+class Jacobian<X : SFun<X>, R: D1, C: D1>(val vfn: VFun<X, R>, vararg val vrbs: Var<X>): MFun<X, R, C>(vfn) {
   override fun invoke(bnds: Bindings<X>) = Mat<X, C, R>(vrbs.map { VDerivative(vfn, it)() }).ᵀ(bnds)
 }
 
-class VDerivative<X : Fun<X>, E: D1> internal constructor(val vFun: VFun<X, E>, val v1: Var<X>) : VFun<X, E>(vFun) {
-  fun VFun<X, E>.df(): VFun<X, E> = when (this) {
+class VDerivative<X : SFun<X>, E: D1> internal constructor(val vFun: VFun<X, E>, val v1: Var<X>) : VFun<X, E>(vFun) {
+  fun VFun<X, E>.df(): VFun<X, E> = when (this@df) {
     is VConst -> VZero()
     is VVar -> VZero()
     is VSum -> left.df() + right.df()
@@ -154,15 +157,12 @@ class VDerivative<X : Fun<X>, E: D1> internal constructor(val vFun: VFun<X, E>,
   }
 }
 
-open class VConst<X: Fun<X>, E: D1>(vararg contents: SConst<X>): Vec<X, E>(emptySet(), contents.asList())
-
-class VZero<X: Fun<X>, E: D1>: VConst<X, E>()
-class VOne<X: Fun<X>, E: D1>: VConst<X, E>()
+open class VConst<X: SFun<X>, E: D1>(vararg contents: SConst<X>): Vec<X, E>(contents.asList())
 
-open class Vec<X: Fun<X>, E: D1>(override val sVars: Set<Var<X>> = emptySet(),
-                                 val contents: List<Fun<X>>): VFun<X, E>() {
-  constructor(contents: List<Fun<X>>): this(contents.flatMap { it.sVars }.toSet(), contents)
+class VZero<X: SFun<X>, E: D1>: VConst<X, E>()
+class VOne<X: SFun<X>, E: D1>: VConst<X, E>()
 
+open class Vec<X: SFun<X>, E: D1>(val contents: List<SFun<X>>): VFun<X, E>(*contents.toTypedArray()) {
   val size = contents.size
   val indices = contents.indices
 
@@ -185,39 +185,39 @@ open class Vec<X: Fun<X>, E: D1>(override val sVars: Set<Var<X>> = emptySet(),
     else -> super.ʘ(multiplicand)
   }
 
-  override fun times(multiplicand: Fun<X>): Vec<X, E> = Vec(contents.map { it * multiplicand })
+  override fun times(multiplicand: SFun<X>): Vec<X, E> = Vec(contents.map { it * multiplicand })
 
   override fun dot(multiplicand: VFun<X, E>) = when(multiplicand) {
     is Vec -> contents.reduceIndexed { index, acc, element -> acc + element * multiplicand[index] }
     else -> super.dot(multiplicand)
   }
 
-  override fun map(ef: (Fun<X>) -> Fun<X>): Vec<X, E> = Vec(contents.map { ef(it) })
+  override fun map(ef: (SFun<X>) -> SFun<X>): Vec<X, E> = Vec(contents.map { ef(it) })
 
   override fun magnitude() = contents.map { it * it }.reduce { acc, p -> acc + p }.sqrt()
 
   override fun unaryMinus(): Vec<X, E> = Vec(contents.map { -it })
 
   companion object {
-    operator fun <T: Fun<T>> invoke(s0: SConst<T>): VConst<T, D1> = VConst(s0)
-    operator fun <T: Fun<T>> invoke(s0: SConst<T>, s1: SConst<T>): VConst<T, D2> = VConst(s0, s1)
-    operator fun <T: Fun<T>> invoke(s0: SConst<T>, s1: SConst<T>, s2: SConst<T>): VConst<T, D3> = VConst(s0, s1, s2)
-    operator fun <T: Fun<T>> invoke(s0: SConst<T>, s1: SConst<T>, s2: SConst<T>, s3: SConst<T>): VConst<T, D4> = VConst(s0, s1, s2, s3)
-    operator fun <T: Fun<T>> invoke(s0: SConst<T>, s1: SConst<T>, s2: SConst<T>, s3: SConst<T>, s4: SConst<T>): VConst<T, D5> = VConst(s0, s1, s2, s3, s4)
-    operator fun <T: Fun<T>> invoke(s0: SConst<T>, s1: SConst<T>, s2: SConst<T>, s3: SConst<T>, s4: SConst<T>, s5: SConst<T>): VConst<T, D6> = VConst(s0, s1, s2, s3, s4, s5)
-    operator fun <T: Fun<T>> invoke(s0: SConst<T>, s1: SConst<T>, s2: SConst<T>, s3: SConst<T>, s4: SConst<T>, s5: SConst<T>, s6: SConst<T>): VConst<T, D7> = VConst(s0, s1, s2, s3, s4, s5, s6)
-    operator fun <T: Fun<T>> invoke(s0: SConst<T>, s1: SConst<T>, s2: SConst<T>, s3: SConst<T>, s4: SConst<T>, s5: SConst<T>, s6: SConst<T>, s7: SConst<T>): VConst<T, D8> = VConst(s0, s1, s2, s3, s4, s5, s6, s7)
-    operator fun <T: Fun<T>> invoke(s0: SConst<T>, s1: SConst<T>, s2: SConst<T>, s3: SConst<T>, s4: SConst<T>, s5: SConst<T>, s6: SConst<T>, s7: SConst<T>, s8: SConst<T>): VConst<T, D9> = VConst(s0, s1, s2, s3, s4, s5, s6, s7, s8)
-
-    operator fun <T: Fun<T>> invoke(t0: Fun<T>): Vec<T, D1> = Vec(arrayListOf(t0))
-    operator fun <T: Fun<T>> invoke(t0: Fun<T>, t1: Fun<T>): Vec<T, D2> = Vec(arrayListOf(t0, t1))
-    operator fun <T: Fun<T>> invoke(t0: Fun<T>, t1: Fun<T>, t2: Fun<T>): Vec<T, D3> = Vec(arrayListOf(t0, t1, t2))
-    operator fun <T: Fun<T>> invoke(t0: Fun<T>, t1: Fun<T>, t2: Fun<T>, t3: Fun<T>): Vec<T, D4> = Vec(arrayListOf(t0, t1, t2, t3))
-    operator fun <T: Fun<T>> invoke(t0: Fun<T>, t1: Fun<T>, t2: Fun<T>, t3: Fun<T>, t4: Fun<T>): Vec<T, D5> = Vec(arrayListOf(t0, t1, t2, t3, t4))
-    operator fun <T: Fun<T>> invoke(t0: Fun<T>, t1: Fun<T>, t2: Fun<T>, t3: Fun<T>, t4: Fun<T>, t5: Fun<T>): Vec<T, D6> = Vec(arrayListOf(t0, t1, t2, t3, t4, t5))
-    operator fun <T: Fun<T>> invoke(t0: Fun<T>, t1: Fun<T>, t2: Fun<T>, t3: Fun<T>, t4: Fun<T>, t5: Fun<T>, t6: Fun<T>): Vec<T, D7> = Vec(arrayListOf(t0, t1, t2, t3, t4, t5, t6))
-    operator fun <T: Fun<T>> invoke(t0: Fun<T>, t1: Fun<T>, t2: Fun<T>, t3: Fun<T>, t4: Fun<T>, t5: Fun<T>, t6: Fun<T>, t7: Fun<T>): Vec<T, D8> = Vec(arrayListOf(t0, t1, t2, t3, t4, t5, t6, t7))
-    operator fun <T: Fun<T>> invoke(t0: Fun<T>, t1: Fun<T>, t2: Fun<T>, t3: Fun<T>, t4: Fun<T>, t5: Fun<T>, t6: Fun<T>, t7: Fun<T>, t8: Fun<T>): Vec<T, D9> = Vec(arrayListOf(t0, t1, t2, t3, t4, t5, t6, t7, t8))
+    operator fun <T: SFun<T>> invoke(s0: SConst<T>): VConst<T, D1> = VConst(s0)
+    operator fun <T: SFun<T>> invoke(s0: SConst<T>, s1: SConst<T>): VConst<T, D2> = VConst(s0, s1)
+    operator fun <T: SFun<T>> invoke(s0: SConst<T>, s1: SConst<T>, s2: SConst<T>): VConst<T, D3> = VConst(s0, s1, s2)
+    operator fun <T: SFun<T>> invoke(s0: SConst<T>, s1: SConst<T>, s2: SConst<T>, s3: SConst<T>): VConst<T, D4> = VConst(s0, s1, s2, s3)
+    operator fun <T: SFun<T>> invoke(s0: SConst<T>, s1: SConst<T>, s2: SConst<T>, s3: SConst<T>, s4: SConst<T>): VConst<T, D5> = VConst(s0, s1, s2, s3, s4)
+    operator fun <T: SFun<T>> invoke(s0: SConst<T>, s1: SConst<T>, s2: SConst<T>, s3: SConst<T>, s4: SConst<T>, s5: SConst<T>): VConst<T, D6> = VConst(s0, s1, s2, s3, s4, s5)
+    operator fun <T: SFun<T>> invoke(s0: SConst<T>, s1: SConst<T>, s2: SConst<T>, s3: SConst<T>, s4: SConst<T>, s5: SConst<T>, s6: SConst<T>): VConst<T, D7> = VConst(s0, s1, s2, s3, s4, s5, s6)
+    operator fun <T: SFun<T>> invoke(s0: SConst<T>, s1: SConst<T>, s2: SConst<T>, s3: SConst<T>, s4: SConst<T>, s5: SConst<T>, s6: SConst<T>, s7: SConst<T>): VConst<T, D8> = VConst(s0, s1, s2, s3, s4, s5, s6, s7)
+    operator fun <T: SFun<T>> invoke(s0: SConst<T>, s1: SConst<T>, s2: SConst<T>, s3: SConst<T>, s4: SConst<T>, s5: SConst<T>, s6: SConst<T>, s7: SConst<T>, s8: SConst<T>): VConst<T, D9> = VConst(s0, s1, s2, s3, s4, s5, s6, s7, s8)
+
+    operator fun <T: SFun<T>> invoke(t0: SFun<T>): Vec<T, D1> = Vec(arrayListOf(t0))
+    operator fun <T: SFun<T>> invoke(t0: SFun<T>, t1: SFun<T>): Vec<T, D2> = Vec(arrayListOf(t0, t1))
+    operator fun <T: SFun<T>> invoke(t0: SFun<T>, t1: SFun<T>, t2: SFun<T>): Vec<T, D3> = Vec(arrayListOf(t0, t1, t2))
+    operator fun <T: SFun<T>> invoke(t0: SFun<T>, t1: SFun<T>, t2: SFun<T>, t3: SFun<T>): Vec<T, D4> = Vec(arrayListOf(t0, t1, t2, t3))
+    operator fun <T: SFun<T>> invoke(t0: SFun<T>, t1: SFun<T>, t2: SFun<T>, t3: SFun<T>, t4: SFun<T>): Vec<T, D5> = Vec(arrayListOf(t0, t1, t2, t3, t4))
+    operator fun <T: SFun<T>> invoke(t0: SFun<T>, t1: SFun<T>, t2: SFun<T>, t3: SFun<T>, t4: SFun<T>, t5: SFun<T>): Vec<T, D6> = Vec(arrayListOf(t0, t1, t2, t3, t4, t5))
+    operator fun <T: SFun<T>> invoke(t0: SFun<T>, t1: SFun<T>, t2: SFun<T>, t3: SFun<T>, t4: SFun<T>, t5: SFun<T>, t6: SFun<T>): Vec<T, D7> = Vec(arrayListOf(t0, t1, t2, t3, t4, t5, t6))
+    operator fun <T: SFun<T>> invoke(t0: SFun<T>, t1: SFun<T>, t2: SFun<T>, t3: SFun<T>, t4: SFun<T>, t5: SFun<T>, t6: SFun<T>, t7: SFun<T>): Vec<T, D8> = Vec(arrayListOf(t0, t1, t2, t3, t4, t5, t6, t7))
+    operator fun <T: SFun<T>> invoke(t0: SFun<T>, t1: SFun<T>, t2: SFun<T>, t3: SFun<T>, t4: SFun<T>, t5: SFun<T>, t6: SFun<T>, t7: SFun<T>, t8: SFun<T>): Vec<T, D9> = Vec(arrayListOf(t0, t1, t2, t3, t4, t5, t6, t7, t8))
   }
 }
 
diff --git a/core/src/main/kotlin/edu/umontreal/kotlingrad/experimental/VariableCapture.kt b/core/src/main/kotlin/edu/umontreal/kotlingrad/experimental/VariableCapture.kt
index 273a0d36..cf033665 100644
--- a/core/src/main/kotlin/edu/umontreal/kotlingrad/experimental/VariableCapture.kt
+++ b/core/src/main/kotlin/edu/umontreal/kotlingrad/experimental/VariableCapture.kt
@@ -77,7 +77,7 @@ open class X<P: Const<P, in Number>>(override val left: BiFn<*>,
     op(left(X), right(X))
 
   private operator fun BiFn<*>.invoke(X: XBnd<P>): P =
-    when (this) {
+    when (this@invoke) {
       is Const<*, *> -> this as P
       is X<*> -> (this as X<P>)(X)
       else -> throw IllegalStateException(toString())
@@ -128,7 +128,7 @@ open class Y<P: Const<P, in Number>>(override val left: BiFn<*>,
     op(left(Y), right(Y))
 
   private operator fun BiFn<*>.invoke(Y: YBnd<P>): P =
-    when (this) {
+    when (this@invoke) {
       is Const<*, *> -> this as P
       is Y<*> -> (this as Y<P>)(Y)
       else -> throw IllegalStateException(toString())
@@ -178,7 +178,7 @@ open class Z<P: Const<P, in Number>>(override val left: BiFn<*>,
   open operator fun invoke(Z: ZBnd<P>): P = op(left(Z), right(Z))
 
   private operator fun BiFn<*>.invoke(Z: ZBnd<P>): P =
-    when (this) {
+    when (this@invoke) {
       is Const<*, *> -> this as P
       is Z<*> -> (this as Z<P>)(Z)
       else -> throw IllegalStateException(toString())
@@ -226,7 +226,7 @@ class XY<P: Const<P, in Number>>(override val left: BiFn<*>,
   operator fun times(that: XYZ<P>): XYZ<P> = XYZ(this, that, mul)
 
   private operator fun BiFn<*>.invoke(X: XBnd<P>, Y: YBnd<P>): P =
-    when (this) {
+    when (this@invoke) {
       is Const<*, *> -> this as P
       is XY<*> -> (this as XY<P>)(X, Y)
       is X<*> -> (this as X<P>)(X)
@@ -237,7 +237,7 @@ class XY<P: Const<P, in Number>>(override val left: BiFn<*>,
   operator fun invoke(X: XBnd<P>, Y: YBnd<P>): P = op(left(X, Y), right(X, Y))
 
   private operator fun BiFn<*>.invoke(X: XBnd<P>): Y<P> =
-    when (this) {
+    when (this@invoke) {
       is Const<*, *> -> Y(this as P)
       is XY<*> -> (this as XY<P>)(X)
       is X<*> -> Y((this as X<P>)(X))
@@ -246,7 +246,7 @@ class XY<P: Const<P, in Number>>(override val left: BiFn<*>,
     }
 
   private operator fun BiFn<*>.invoke(Y: YBnd<P>): X<P> =
-    when (this) {
+    when (this@invoke) {
       is Const<*, *> -> X(this as P)
       is XY<*> -> (this as XY<P>)(Y)
       is Y<*> -> X((this as Y<P>)(Y))
@@ -299,7 +299,7 @@ class XZ<P: Const<P, in Number>>(override val left: BiFn<*>,
   operator fun times(that: XYZ<P>): XYZ<P> = XYZ(this, that, mul)
 
   private operator fun BiFn<*>.invoke(X: XBnd<P>, Z: ZBnd<P>): P =
-    when (this) {
+    when (this@invoke) {
       is Const<*, *> -> this as P
       is XZ<*> -> (this as XZ<P>)(X, Z)
       is X<*> -> (this as X<P>)(X)
@@ -310,7 +310,7 @@ class XZ<P: Const<P, in Number>>(override val left: BiFn<*>,
   operator fun invoke(X: XBnd<P>, Z: ZBnd<P>): P = op(left(X, Z), right(X, Z))
 
   private operator fun BiFn<*>.invoke(X: XBnd<P>): Z<P> =
-    when (this) {
+    when (this@invoke) {
       is Const<*, *> -> Z(this as P)
       is XZ<*> -> (this as XZ<P>)(X)
       is X<*> -> Z((this as X<P>)(X))
@@ -319,7 +319,7 @@ class XZ<P: Const<P, in Number>>(override val left: BiFn<*>,
     }
 
   private operator fun BiFn<*>.invoke(Z: ZBnd<P>): X<P> =
-    when (this) {
+    when (this@invoke) {
       is Const<*, *> -> X(this as P)
       is XZ<*> -> (this as XZ<P>)(Z)
       is Z<*> -> X((this as Z<P>)(Z))
@@ -372,7 +372,7 @@ class YZ<P: Const<P, in Number>>(override val left: BiFn<*>,
   operator fun times(that: XYZ<P>): XYZ<P> = XYZ(this, that, mul)
 
   private operator fun BiFn<*>.invoke(Y: YBnd<P>, Z: ZBnd<P>): P =
-    when (this) {
+    when (this@invoke) {
       is Const<*, *> -> this as P
       is YZ<*> -> (this as YZ<P>)(Y, Z)
       is Y<*> -> (this as Y<P>)(Y)
@@ -384,7 +384,7 @@ class YZ<P: Const<P, in Number>>(override val left: BiFn<*>,
     op(left(Y, Z), right(Y, Z))
 
   private operator fun BiFn<*>.invoke(Y: YBnd<P>): Z<P> =
-    when (this) {
+    when (this@invoke) {
       is Const<*, *> -> Z(this as P)
       is YZ<*> -> (this as YZ<P>)(Y)
       is Y<*> -> Z((this as Y<P>)(Y))
@@ -393,7 +393,7 @@ class YZ<P: Const<P, in Number>>(override val left: BiFn<*>,
     }
 
   private operator fun BiFn<*>.invoke(Z: ZBnd<P>): Y<P> =
-    when (this) {
+    when (this@invoke) {
       is Const<*, *> -> Y(this as P)
       is YZ<*> -> (this as YZ<P>)(Z)
       is Z<*> -> Y((this as Z<P>)(Z))
@@ -444,7 +444,7 @@ class XYZ<P: Const<P, in Number>>(override val left: BiFn<*>,
   operator fun times(that: XYZ<P>): XYZ<P> = XYZ(this, that, mul)
 
   private operator fun BiFn<*>.invoke(X: XBnd<P>, Y: YBnd<P>, Z: ZBnd<P>): P =
-    when (this) {
+    when (this@invoke) {
       is Const<*, *> -> this as P
       is XYZ<*> -> (this as XYZ<P>)(X, Y, Z)
       is XY<*> -> (this as XY<P>)(X, Y)
@@ -460,7 +460,7 @@ class XYZ<P: Const<P, in Number>>(override val left: BiFn<*>,
     op(left(X, Y, Z), right(X, Y, Z))
 
   private operator fun BiFn<*>.invoke(X: XBnd<P>, Z: ZBnd<P>): Y<P> =
-    when (this) {
+    when (this@invoke) {
       is Const<*, *> -> Y(this as P)
       is XYZ<*> -> (this as XYZ<P>)(X, Z)
       is XY<*> -> (this as XY<P>)(X)
@@ -476,7 +476,7 @@ class XYZ<P: Const<P, in Number>>(override val left: BiFn<*>,
     Y(left(X, Z), right(X, Z), op)
 
   private operator fun BiFn<*>.invoke(X: XBnd<P>, Y: YBnd<P>): Z<P> =
-    when (this) {
+    when (this@invoke) {
       is Const<*, *> -> Z(this as P)
       is XYZ<*> -> (this as XYZ<P>)(X, Y)
       is XY<*> -> Z((this as XY<P>)(X, Y))
@@ -492,7 +492,7 @@ class XYZ<P: Const<P, in Number>>(override val left: BiFn<*>,
     Z(left(X, Y), right(X, Y), op)
 
   private operator fun BiFn<*>.invoke(Y: YBnd<P>, Z: ZBnd<P>): X<P> =
-    when (this) {
+    when (this@invoke) {
       is Const<*, *> -> X(this as P)
       is XYZ<*> -> (this as XYZ<P>)(Y, Z)
       is XY<*> -> (this as XY<P>)(Y)
@@ -508,7 +508,7 @@ class XYZ<P: Const<P, in Number>>(override val left: BiFn<*>,
     X(left(Y, Z), right(Y, Z), op)
 
   private operator fun BiFn<*>.invoke(X: XBnd<P>): YZ<P> =
-    when (this) {
+    when (this@invoke) {
       is Const<*, *> -> YZ(this as P)
       is XYZ<*> -> (this as XYZ<P>)(X)
       is XY<*> -> YZ((this as XY<P>)(X))
@@ -523,7 +523,7 @@ class XYZ<P: Const<P, in Number>>(override val left: BiFn<*>,
   operator fun invoke(X: XBnd<P>): YZ<P> = YZ(left(X), right(X), op)
 
   private operator fun BiFn<*>.invoke(Y: YBnd<P>): XZ<P> =
-    when (this) {
+    when (this@invoke) {
       is Const<*, *> -> XZ(this as P)
       is XYZ<*> -> (this as XYZ<P>)(Y)
       is XY<*> -> XZ((this as XY<P>)(Y))
@@ -538,7 +538,7 @@ class XYZ<P: Const<P, in Number>>(override val left: BiFn<*>,
   operator fun invoke(Y: YBnd<P>): XZ<P> = XZ(left(Y), right(Y), op)
 
   private operator fun BiFn<*>.invoke(Z: ZBnd<P>): XY<P> =
-    when (this) {
+    when (this@invoke) {
       is Const<*, *> -> XY(this as P)
       is XYZ<*> -> (this as XYZ<P>)(Z)
       is XY<*> -> this as XY<P>
diff --git a/core/src/test/kotlin/edu/umontreal/kotlingrad/calculus/ExpressionGenerator.kt b/core/src/test/kotlin/edu/umontreal/kotlingrad/calculus/ExpressionGenerator.kt
index 984e025c..5ac62460 100644
--- a/core/src/test/kotlin/edu/umontreal/kotlingrad/calculus/ExpressionGenerator.kt
+++ b/core/src/test/kotlin/edu/umontreal/kotlingrad/calculus/ExpressionGenerator.kt
@@ -7,24 +7,24 @@ import edu.umontreal.kotlingrad.experimental.DoublePrecision.z
 import io.kotlintest.properties.Gen
 import io.kotlintest.properties.shrinking.Shrinker
 
-abstract class ExpressionGenerator<X: Fun<X>>: Gen<Fun<X>> {
+abstract class ExpressionGenerator<X: SFun<X>>: Gen<SFun<X>> {
   companion object: ExpressionGenerator<DReal>() {
     override val variables: List<Var<DReal>> = listOf(x, y, z)
   }
 
-  val sum = { x: Fun<X>, y: Fun<X> -> Sum(x, y) }
-  val mul = { x: Fun<X>, y: Fun<X> -> Prod(x, y) }
+  val sum = { x: SFun<X>, y: SFun<X> -> Sum(x, y) }
+  val mul = { x: SFun<X>, y: SFun<X> -> Prod(x, y) }
 
-  val operators: List<(Fun<X>, Fun<X>) -> Fun<X>> = listOf(sum, mul)
+  val operators: List<(SFun<X>, SFun<X>) -> SFun<X>> = listOf(sum, mul)
   val constants: List<SConst<X>> = listOf(Zero(), One(), Two())
   open val variables: List<Var<X>> = listOf(Var("x"), Var("y"), Var("z"))
 
-  override fun constants(): Iterable<Fun<X>> = constants
+  override fun constants(): Iterable<SFun<X>> = constants
 
-  override fun random(): Sequence<Fun<X>> = generateSequence { randomBiTree() }
+  override fun random(): Sequence<SFun<X>> = generateSequence { randomBiTree() }
 
-  override fun shrinker() = object: Shrinker<Fun<X>> {
-    override fun shrink(failure: Fun<X>): List<Fun<X>> =
+  override fun shrinker() = object: Shrinker<SFun<X>> {
+    override fun shrink(failure: SFun<X>): List<SFun<X>> =
       when(failure) {
         is Sum -> listOf(failure.left, failure.right)
         is Prod -> listOf(failure.left, failure.right)
@@ -32,7 +32,7 @@ abstract class ExpressionGenerator<X: Fun<X>>: Gen<Fun<X>> {
       }
   }
 
-  private fun randomBiTree(level: Int = 1): Fun<X> =
+  private fun randomBiTree(level: Int = 1): SFun<X> =
     if(5 < level)
       if(Math.random() < 0.5) constants.random() else variables.random()
     else operators.random()(randomBiTree(level + 1), randomBiTree(level + 1))
diff --git a/core/src/test/kotlin/edu/umontreal/kotlingrad/evaluation/TestSymbolic.kt b/core/src/test/kotlin/edu/umontreal/kotlingrad/evaluation/TestSymbolic.kt
index 5953b461..917c046e 100644
--- a/core/src/test/kotlin/edu/umontreal/kotlingrad/evaluation/TestSymbolic.kt
+++ b/core/src/test/kotlin/edu/umontreal/kotlingrad/evaluation/TestSymbolic.kt
@@ -14,7 +14,7 @@ import javax.script.SimpleBindings
 class TestSymbolic : StringSpec({
   val engine = ScriptEngineManager().getEngineByExtension("kts")
 
-  fun ktEval(f: Fun<DReal>, vararg kgBnds: Pair<Var<DReal>, Number>) =
+  fun ktEval(f: SFun<DReal>, vararg kgBnds: Pair<Var<DReal>, Number>) =
     engine.run {
       val bindings = kgBnds.map { it.first.name to it.second.toDouble() }.toMap()
       setBindings(SimpleBindings(bindings), ScriptContext.ENGINE_SCOPE)
@@ -23,7 +23,7 @@ class TestSymbolic : StringSpec({
 
   with(DoublePrecision) {
     "test symbolic evaluation" {
-      ExpressionGenerator.assertAll(10) { f: Fun<DReal> ->
+      ExpressionGenerator.assertAll(10) { f: SFun<DReal> ->
         try {
           DoubleGenerator.assertAll(10) { ẋ, ẏ, ż ->
             f(x to ẋ, y to ẏ, z to ż) shouldBeAbout ktEval(f, x to ẋ, y to ẏ, z to ż)
diff --git a/samples/src/main/kotlin/edu/umontreal/kotlingrad/samples/MultilayerPerceptron.kt b/samples/src/main/kotlin/edu/umontreal/kotlingrad/samples/MultilayerPerceptron.kt
index b1ad29c6..f747e402 100644
--- a/samples/src/main/kotlin/edu/umontreal/kotlingrad/samples/MultilayerPerceptron.kt
+++ b/samples/src/main/kotlin/edu/umontreal/kotlingrad/samples/MultilayerPerceptron.kt
@@ -4,7 +4,7 @@ import edu.umontreal.kotlingrad.experimental.*
 import kotlin.random.Random
 
 @Suppress("NonAsciiCharacters")
-class MultilayerPerceptron<T: Fun<T>>(
+class MultilayerPerceptron<T: SFun<T>>(
   val x: Var<T> = Var(),
   val y: Var<T> = Var(),
   val p1v: VVar<T, D3> = VVar(),
@@ -12,12 +12,12 @@ class MultilayerPerceptron<T: Fun<T>>(
   val p3v: VVar<T, D3> = VVar()
 ): (VFun<T, D3>,
     MFun<T, D3, D3>,
-    VFun<T, D3>) -> Fun<T> {
+    VFun<T, D3>) -> SFun<T> {
   override operator fun invoke(
     p1: VFun<T, D3>,
     p2: MFun<T, D3, D3>,
     p3: VFun<T, D3>
-  ): Fun<T> {
+  ): SFun<T> {
     val layer1 = layer(p1 * x)
     val layer2 = layer(p2 * layer1)
     val output = layer2 dot p3
@@ -25,7 +25,7 @@ class MultilayerPerceptron<T: Fun<T>>(
     return lossFun
   }
 
-  private fun sigmoid(x: Fun<T>) = One<T>() / (One<T>() + E<T>().pow(-x))
+  private fun sigmoid(x: SFun<T>) = One<T>() / (One<T>() + E<T>().pow(-x))
 
   private fun layer(x: VFun<T, D3>): VFun<T, D3> = x.map { sigmoid(it) }
 
diff --git a/samples/src/main/kotlin/edu/umontreal/kotlingrad/samples/VisualizeDFG.kt b/samples/src/main/kotlin/edu/umontreal/kotlingrad/samples/VisualizeDFG.kt
index bd8d298c..0620ae3e 100644
--- a/samples/src/main/kotlin/edu/umontreal/kotlingrad/samples/VisualizeDFG.kt
+++ b/samples/src/main/kotlin/edu/umontreal/kotlingrad/samples/VisualizeDFG.kt
@@ -1,7 +1,7 @@
 package edu.umontreal.kotlingrad.samples
 
 import edu.umontreal.kotlingrad.experimental.DoublePrecision
-import edu.umontreal.kotlingrad.experimental.Fun
+import edu.umontreal.kotlingrad.experimental.SFun
 import guru.nidi.graphviz.*
 import guru.nidi.graphviz.attribute.*
 import guru.nidi.graphviz.attribute.Rank.RankDir.LEFT_TO_RIGHT
@@ -20,7 +20,7 @@ fun main() {
 const val DARKMODE = false
 const val THICKNESS = 2
 
-fun Fun<*>.render(filename: String? = null) {
+fun SFun<*>.render(filename: String? = null) {
   val image = graph(directed = true) {
     val color = if (DARKMODE) Color.WHITE else Color.BLACK
 
diff --git a/samples/src/main/kotlin/edu/umontreal/kotlingrad/samples/physics/DoublePendulum.kt b/samples/src/main/kotlin/edu/umontreal/kotlingrad/samples/physics/DoublePendulum.kt
index 077f44f6..199a1fc8 100644
--- a/samples/src/main/kotlin/edu/umontreal/kotlingrad/samples/physics/DoublePendulum.kt
+++ b/samples/src/main/kotlin/edu/umontreal/kotlingrad/samples/physics/DoublePendulum.kt
@@ -19,19 +19,19 @@ import kotlin.math.*
 
 @Suppress("NonAsciiCharacters", "LocalVariableName")
 class DoublePendulum(private val len: Double = 900.0) : Application(), EventHandler<ActionEvent> {
-  var ω1: Fun<DReal> = DoublePrecision.wrap(0.0) // Angular velocities
-  var ω2: Fun<DReal> = DoublePrecision.wrap(0.0)
+  var ω1: SFun<DReal> = DoublePrecision.wrap(0.0) // Angular velocities
+  var ω2: SFun<DReal> = DoublePrecision.wrap(0.0)
   val m1 = 2.0 // Masses
   val m2 = 2.0
-  var G: Fun<DReal> = DoublePrecision.wrap(9.81) // Gravity
-  var µ: Fun<DReal> = DoublePrecision.wrap(0.01) // Friction
+  var G: SFun<DReal> = DoublePrecision.wrap(9.81) // Gravity
+  var µ: SFun<DReal> = DoublePrecision.wrap(0.01) // Friction
   val Gp = 0.01  // Simulate measurement error
   val µp = -0.01
   var r1 = DoublePrecision.Vec(1.0, 0.0) // Polar vector
   var r2 = DoublePrecision.Vec(1.0, 0.0)
   val observationSteps = 30
   var priorVal = 5.0
-  fun step(obs: Fun<DReal>? = null, groundTruth: Pair<VConst<DReal, D2>, VConst<DReal, D2>>? = null) = with(DoublePrecision) {
+  fun step(obs: SFun<DReal>? = null, groundTruth: Pair<VConst<DReal, D2>, VConst<DReal, D2>>? = null) = with(DoublePrecision) {
     val isObserving = false
 //    val priorVal = if(G is Var) G.asDouble()
     if (isObserving) {
@@ -69,10 +69,10 @@ class DoublePendulum(private val len: Double = 900.0) : Application(), EventHand
     }
 
     val r1a = (r1.angle + (ω1 * dt + .5 * α1 * dt * dt)).run {
-      if(G is Var) this(G as Var<DReal> to priorVal).asDouble() else try { asDouble() } catch(e: Exception) {println(this); this(sVars.first() to priorVal).asDouble() }
+      if(G is Var) this(G as Var<DReal> to priorVal).asDouble() else try { asDouble() } catch(e: Exception) {println(this); this(inputs.sVars.first() to priorVal).asDouble() }
     }
     val r2a = (r2.angle + - (ω2 * dt + .5 * α2 * dt * dt)).run {
-      if(G is Var) this(G as Var<DReal> to priorVal).asDouble() else try { asDouble() } catch(e: Exception) {println(this); this(sVars.first() to priorVal).asDouble() }
+      if(G is Var) this(G as Var<DReal> to priorVal).asDouble() else try { asDouble() } catch(e: Exception) {println(this); this(inputs.sVars.first() to priorVal).asDouble() }
     }
 
     if(G is Var) {
@@ -90,11 +90,11 @@ class DoublePendulum(private val len: Double = 900.0) : Application(), EventHand
     }
   }
 
-  fun Fun<DReal>.descend(steps: Int, vinit: Double, gamma: Double, α: Double = 0.1, map: Pair<Var<DReal>, DReal>): Fun<DReal> {
+  fun SFun<DReal>.descend(steps: Int, vinit: Double, gamma: Double, α: Double = 0.1, map: Pair<Var<DReal>, DReal>): SFun<DReal> {
     with(DoublePrecision) {
       val d_dg = this@descend.d(map.first)
-      var G1P: Fun<DReal> = map.second
-      var velocity: Fun<DReal> = wrap(vinit)
+      var G1P: SFun<DReal> = map.second
+      var velocity: SFun<DReal> = wrap(vinit)
       var i = 0
       do {
         velocity = gamma * velocity + d_dg(map.first to G1P) * α