anchor.diff:      /// Propagates the adjoint `v` backwards through the evaluation trace of `d`. The adjoints in the trace are reset before the push.
       let rec reverseProp (adjoints: Adjoints) (v:dobj) (d:dobj) =
           let fanouts = reverseReset adjoints d
  -        let inline bx v d = bxd v, bxd d
  +        let inline bxd (x : dobj) = x
  +        let inline bxdelta (x : delta) = x
  +
  +        let inline bd (v: Delta) (d:D) = bxdelta v, bxd d
  +        let inline bdv (v: DeltaV) (d:DV) = bxdelta v, bxd d
  +        let inline bdm (v: DeltaM) (d:DM) = bxdelta v, bxd d
  +
  +        let inline bx (v: D) (d:D) = bd (X v) d
  +        let inline bxv (v: DV) (d:DV) = bdv (XV v) d
  +        let inline bxm (v: DM) (d:DM) = bdm (XM v) d
  +
           // Note, this uses an explicit worklist over (D*D|DV*DV|DM*DM) to make it tail-recursive
  -        let rec pushRec (ds:(dobj*dobj) list) =
  +        let rec pushRec (ds:(delta*dobj) list) =
               match ds with
               | [] -> ()
               | (v, d) :: t ->
  -                match d, v with
  -                | (:? D as d), (:? D as v) ->
  +                match v, d with
  +                | (:? Delta as delta),  (:? D as d) ->
                       match d with
                       | DR(_, o, _, uniq) ->
  -                        let dA = adjoints.GetD(uniq) + v
  +                        let dA = adjoints.ApplyDelta(uniq, delta) 
                           let fanout = fanouts.[uniq] - 1u
  -                        adjoints.SetD(uniq,dA)
                           fanouts.[uniq] <- fanout
                           // If all incoming parts of the adjoint have been received, then proceed to the parent
                           if fanout = 0u then
  
positive.code:     override __.ToString() = sprintf "(%d computed adjoints)" dict.Count
  
  
  /// Functional-oriented operations on vectors. Implementing functionality similar to FSharp.Collections.Array.
  [<RequireQualifiedAccess>]
  [<CompilationRepresentation (CompilationRepresentationFlags.ModuleSuffix)>]
  module DV =
  
      // Note: map operations are not implemented on purpose. To benefit from the performance of BLAS ops, supplied element-wise operations are used. For example: "exp v" instead of "DV.map exp v"
      /// Creates a vector from array `a`
      let inline ofArray a = DV.OfArray(a)
  
      /// Converts vector `v` into an array
      let inline toArray (v:DV) = v.ToArray()
  
      /// Converts vector `v` into a row matrix
      let inline toRowDM (v:DV) = v.ToRowDM()
  
      /// Converts vector `v` into a column matrix
      let inline toColDM (v:DV) = v.ToColDM()
  
      /// Creates a vector with `n` elements, each with value `v`
      let inline create n (v:'a) = 
          let at = typeof<'a>
          if at.Equals(typeof<D>) then DV.OfArray(Array.create n (unbox<D>(box v)))
          elif at.Equals(typeof<number>) then DV (Array.create n (unbox<number>(box v)))
          elif at.Equals(typeof<int>) then DV (Array.create n (unbox<int>(box v) |> N.toNumber))
          else N.failWithInvalidTypeMessage ()
  
      /// Creates a vector with `n` zero elements
      let inline zeroCreate n = DV.ZeroN n
  
      /// Empty vector
      let empty = DV.Zero
  
      /// Creates a vector of `n` elements, where each element is defined by function `f`
      let inline init n (f:int->'a) = 
          let at = typeof<'a>
          if at.Equals(typeof<D>) then DV.OfArray(Array.init n (unbox<int->D>(box f)))
          elif at.Equals(typeof<number>) then DV (Array.init n (unbox<int->number>(box f)))
          elif at.Equals(typeof<int>) then DV ((Array.init n (unbox<int->int>(box f))) |> Array.map N.toNumber)
          else N.failWithInvalidTypeMessage ()
  
      /// Returns true if vector `v` is empty, otherwise returns false
      let isEmpty (v:DV) = v.Length = 0
  
      /// Iterates function `f` over the elements of vector `v`
      let inline iter (f:D->unit) (v:DV) = v |> toArray |> Array.iter f
  
      /// Iterates function `f` over the elements of vector `v`. An element index is also supplied to `f`.
      let inline iteri (f:int->D->unit) (v:DV) = v |> toArray |> Array.iteri f
  
      /// Iterates function `f` over the elements of vectors `v1` and `v2`
      let inline iter2 (f:D->D->unit) (v1:DV) (v2:DV) = Array.iter2 f (v1 |> toArray) (v2 |> toArray)
  
      /// Iterates function `f` over the elements of vectors `v1` and `v2`. An element index is also supplied to `f`.
      let inline iteri2 (f:int->D->D->unit) (v1:DV) (v2:DV) = Array.iteri2 f (v1 |> toArray) (v2 |> toArray)
  
      /// Length of vector `v`
      let inline length (v:DV) = v.Length
  
      /// L1 norm of vector `v`
      let inline l1norm (v:DV) = DV.L1Norm(v)
  
      /// L2 norm of vector `v`
      let inline l2norm (v:DV) = DV.L2Norm(v)
  
      /// Squared L2 norm of vector `v`
      let inline l2normSq (v:DV) = DV.L2NormSq(v)
  
      /// Maximum of the elements of vector `v`
      let inline max (v:DV) = DV.Max(v)
  
      /// Index of the maximum element of vector `v`
      let inline maxIndex (v:DV) = DV.MaxIndex(v)
  
      /// Minimum of the elements of vector `v`
      let inline min (v:DV) = DV.Min(v)
  
      /// Index of the minimum element of vector `v`
      let inline minIndex (v:DV) = DV.MinIndex(v)
  
      /// Mean of vector `v`
      let inline mean (v:DV) = DV.Mean(v)
  
      /// Average of vector `v`. Same with mean.
      let average = mean
  
      /// Standard deviation of vector `v`
      let inline standardDev (v:DV) = DV.StandardDev(v)
  
      /// Variance of vector `v`
      let inline variance (v:DV) = DV.Variance(v)
  
      /// Shift and scale the elements of vector `v` to have zero mean and unit variance
      let inline standardize (v:DV) = DV.Standardize(v)
  
      /// Shift and scale the elements of vector `v` to be in the range [0, 1]
      let inline normalize (v:DV) = DV.Normalize(v)
  
      /// L2 norm of vector `v`. Same with DV.l2norm.
      let inline norm (v:DV) = DV.L2Norm(v)
  
      /// Squared L2 norm of vector `v`. Same with DV.l2normSq.
      let inline normSq(v:DV) = DV.L2NormSq(v)
  
      // TODO: implement supNorm (infinity norm, with BLAS IDAMAX)
      /// Creates a vector where elements of `v1` are followed by elements of `v2`
      let inline append (v1:DV) (v2:DV) = DV.Append(v1, v2)
  
      /// Creates a vector where elements of `v2` are followed by elements of `v1`
      let inline prepend (v1:DV) (v2:DV) = DV.Append(v2, v1)
  
      /// Concatenates the given sequence of vectors `v` into one vector
      let inline concat (v:seq<DV>) = Seq.fold append DV.Zero v
  
      /// Splits vector `v` into a sequence of subvectors whose lengths are given in sequence `n`
      let inline split (n:seq<int>) (v:DV) = DV.Split(v, n)
  
      /// Splits vector `v` into `n` subvectors of equal length. The length of vector `v` must be an integer multiple of `n`.
      let inline splitEqual (n:int) (v:DV) = DV.Split(v, Array.create n (v.Length / n))
  
      /// Sums the elements of vector `v`
      let inline sum (v:DV) = DV.Sum(v)
  
      /// Creates a vector with `n` elements where the `i`-th element is one and the rest of the elements are zero
      let inline standardBasis (n:int) (i:int) = DV(standardBasis n i)
  
      /// Creates a vector with `n` elements where the `i`-th element has value `v` and the rest of the elements are zero
      let inline standardBasisVal (n:int) (i:int) (v:number) = DV(standardBasisVal n i v)
  
      /// Gets the unit vector codirectional with vector `v`
      let inline unitDV (v:DV) = v / DV.L2Norm(v)
  
      /// Converts matrix `m` into a vector by stacking its rows
      let inline ofDM (m:DM) = DM.ReshapeToDV(m)
  
      /// Creates a matrix with `m` rows from vector `v`
      let inline toDM (m:int) (v:DV) = DV.ReshapeToDM(m, v)
  
      // Experimental
      let inline toString (v:DV) = v.ToString()
      let inline visualize (v:DV) = v.Visualize()
      let inline visualizeAsDM (m:int) (v:DV) = DV.ReshapeToDM(m, v).Visualize()
  
  
  /// Functional-oriented operations on matrices. Implementing functionality similar to FSharp.Collections.Array2D.
  [<RequireQualifiedAccess>]
  [<CompilationRepresentation (CompilationRepresentationFlags.ModuleSuffix)>]
  module DM =
  
      /// Creates a matrix from 2D array `a`
      let inline ofArray2D a = DM.OfArray2D(a)
  
      /// Converts matrix `m` into a 2D array
      let inline toArray2D (m:DM) = m.GetRows() |> Seq.map DV.toArray |> array2D
  
      /// Creates a matrix with `m` rows from array `a`
      let inline ofArray m a = DM.OfArray(m, a)
  
      /// Converts matrix `m` into an array by stacking its rows
      let inline toArray (m:DM) = DM.ReshapeToDV(m) |> DV.toArray
  
      /// Transpose of matrix `m`
      let inline transpose (m:DM) = DM.Transpose(m)
  
      /// Creates a matrix from a sequence of row vectors `s`
      let inline ofRows s = DM.OfRows(s)
  
      /// Creates a matrix from a sequence of column vectors `s`
      let inline ofCols (s:seq<DV>) = s |> ofRows |> transpose
  
      /// Gets the sequence of row vectors in matrix `m`
      let inline toRows (m:DM) = m.GetRows()
  
      /// Gets the sequence of column vectors in matrix `m`
      let inline toCols (m:DM) = m.GetCols()
  
      /// Converts matrix `m` into a vector by stacking its rows
      let inline toDV (m:DM) = DM.ReshapeToDV(m)
  
      /// Creates a matrix with `m` rows from vector `v`
      let inline ofDV (m:int) (v:DV) = DV.ReshapeToDM(m, v)
  
      /// Gets the column with index `j` of matrix `m`
      let inline col (j:int) (m:DM) = m.[*, j]
  
      /// Gets the row with index `i` of matrix `m`
      let inline row (i:int) (m:DM) = m.[i, *]
  
      /// Number of columns in matrix `m`
      let inline cols (m:DM) = m.Cols
  
      /// Number of rows in matrix `m`
      let inline rows (m:DM) = m.Rows
  
      /// Creates a matrix with `m` rows and `n` columns, where all entries have value `v`
      let inline create m n (v:'a) = 
          let at = typeof<'a>
          if at.Equals(typeof<D>) then DM.OfArray2D(Array2D.create m n (unbox<D>(box v)))
          elif at.Equals(typeof<number>) then DM (Array2D.create m n (unbox<number>(box v)))
          elif at.Equals(typeof<int>) then DM (Array2D.create m n (unbox<int>(box v) |> N.toNumber))
          else N.failWithInvalidTypeMessage ()
  
      /// Creates a matrix with `m` rows, where all rows are equal to `v`
      let inline createRows (m:int) (v:DV) = DM.OfRows(m, v)
  
      /// Creates a matrix with `n` columns, where all columns are equal to `v`
      let inline createCols (n:int) (v:DV) = DM.OfCols(n, v)
  
      /// Creates a matrix with `m` rows and `n` columns, where all entries are zero
      let inline zeroCreate m n = DM.ZeroMN m n
  
      /// Gets the diagonal of matrix `m`
      let inline diagonal (m:DM) = DM.Diagonal(m)
  
      /// Zero matrix
      let empty = DM.Zero
  
      /// Returns true if matrix `m` is empty, otherwise returns false
      let isEmpty (m:DM) = m.Length = 0
  
      /// Creates a matrix with `m` rows and `n` columns, where each element is given by function `f`
      let inline init m n (f:int->int->'a) = 
          let at = typeof<'a>
          if at.Equals(typeof<D>) then DM.OfArray2D(Array2D.init m n (unbox<int->int->D>(box f)))
          elif at.Equals(typeof<number>) then DM (Array2D.init m n (unbox<int->int->number>(box f)))
          elif at.Equals(typeof<int>) then DM ((Array2D.init m n (unbox<int->int->int>(box f))) |> Array2D.map N.toNumber)
          else N.failWithInvalidTypeMessage ()
  
      /// Creates a matrix with `m` rows, where each row is given by `f` as a vector
      let inline initRows (m:int) (f:int->DV) = Seq.init m f |> ofRows
  
      /// Creates a matrix with `n` columns, where each column is given by `f` as a vector
      let inline initCols (n:int) (f:int->DV) = Seq.init n f |> ofCols
  
      /// Inverse of matrix `m`
      let inline inverse (m:DM) = DM.Inverse(m)
  
      /// Iterates function `f` over the entries of matrix `m`
      let inline iter (f:D->unit) (m:DM) = m |> toDV |> DV.iter f
  
      /// Iterates function `f` over the entries of matrices `m1` and `m2`
      let inline iter2 (f:D->D->unit) (m1:DM) (m2:DM) = DV.iter2 f (m1 |> toDV) (m2 |> toDV)
  
      /// Iterates function `f` over the entries of matrix `m`. Indices are also supplied to `f`.
      let inline iteri (f:int->int->D->unit) (m:DM) = m |> toArray2D |> Array2D.iteri f
  
      /// Iterates function `f` over the columns of matrix `m`
      let inline iterCols (f:DV->unit) (m:DM) = m |> toCols |> Seq.iter f
  
      /// Iterates function `f` over the rows of matrix `m`
      let inline iterRows (f:DV->unit) (m:DM) = m |> toRows |> Seq.iter f
  
      /// Iterates function `f` over the columns of matrix `m`. Column indices are also supplied to `f`.
      let inline iteriCols (f:int->DV->unit) (m:DM) = m |> toCols |> Seq.iteri f
  
      /// Iterates function `f` over the rows of matrix `m`. Row indices are also supplied to `f`.
      let inline iteriRows (f:int->DV->unit) (m:DM) = m |> toRows |> Seq.iteri f
  
      /// Iterates function `f` over the columns of matrices `m1` and `m2`
      let inline iter2Cols (f:DV->DV->unit) (m1:DM) (m2:DM) = Seq.iter2 f (m1 |> toCols) (m2 |> toCols)
  
      /// Iterates function `f` over the rows of matrices `m1` and `m2
      let inline iter2Rows (f:DV->DV->unit) (m1:DM) (m2:DM) = Seq.iter2 f (m1 |> toRows) (m2 |> toRows)
  
      /// Iterates function `f` over the columns of matrices `m1` and `m2`. Column indices are also supplied to `f`.
      let inline iteri2Cols (f:int->DV->DV->unit) (m1:DM) (m2:DM) = Seq.iteri2 f (m1 |> toCols) (m2 |> toCols)
  
      /// Iterates function `f` over the rows of matrices `m1` and `m2`. Row indices are also supplied to `f`.
      let inline iteri2Rows (f:int->DV->DV->unit) (m1:DM) (m2:DM) = Seq.iteri2 f (m1 |> toRows) (m2 |> toRows)
  
      /// Total number of elements in matrix `m`
      let inline length (m:DM) = m.Length
  
      /// Number of rows in matrix `m`. Same with DM.rows.
      let inline length1 (m:DM) = m.Rows
  
      /// Number of columns in matrix `m`. Same with DM.cols.
      let inline length2 (m:DM) = m.Cols
  
      /// Determinant of matrix `m`
      let inline det (m:DM) = DM.Det(m)
  
      /// Maps function `f` to the columns of matrix `m`
      let inline mapCols (f:DV->DV) (m:DM) = m |> toCols |> Seq.map f |> ofCols
  
      /// Maps function `f` to the rows of matrix `m`
      let inline mapRows (f:DV->DV) (m:DM) = m |> toRows |> Seq.map f |> ofRows
  
      /// Maps function `f` to the columns of matrix `m`. Column indices are also supplied to `f`.
      let inline mapiCols (f:int->DV->DV) (m:DM) = m |> toCols |> Seq.mapi f |> ofCols
  
      /// Maps function `f` to the rows of matrix `m`. Row indices are also supplied to `f`.
      let inline mapiRows (f:int->DV->DV) (m:DM) = m |> toRows |> Seq.mapi f |> ofRows
  
      /// Maps function `f` to the columns of matrices `m1` and `m2`
      let inline map2Cols (f:DV->DV->DV) (m1:DM) (m2:DM) = Seq.map2 f (m1 |> toCols) (m2 |> toCols) |> ofCols
  
      /// Maps function `f` to the rows of matrices `m1` and `m2`
      let inline map2Rows (f:DV->DV->DV) (m1:DM) (m2:DM) = Seq.map2 f (m1 |> toRows) (m2 |> toRows) |> ofRows
  
      /// Maps function `f` to the columns of matrices `m1` and `m2`. Column indices are also supplied to `f`.
      let inline mapi2Cols (f:int->DV->DV->DV) (m1:DM) (m2:DM) = Seq.mapi2 f (m1 |> toCols) (m2 |> toCols) |> ofCols
  
      /// Maps function `f` to the rows of matrices `m1` and `m2`. Row indices are also supplied to `f`.
      let inline mapi2Rows (f:int->DV->DV->DV) (m1:DM) (m2:DM) = Seq.mapi2 f (m1 |> toRows) (m2 |> toRows) |> ofRows
  
      /// Maximum of the entries of matrix `m`
      let inline max (m:DM) = DM.Max(m)
  
      /// Index of the maximum entry of matrix `m`
      let inline maxIndex (m:DM) = DM.MaxIndex(m)
  
      /// Minimum of the entries of matrix `m`
      let inline min (m:DM) = DM.Min(m)
  
      /// Index of the minimum entry of matrix `m`
      let inline minIndex (m:DM) = DM.MinIndex(m)
  
      /// Mean of matrix `m`
      let inline mean (m:DM) = DM.Mean(m)
  
      /// Average of matrix `m`. Same with mean.
      let average = mean
  
      /// Standard deviation of matrix `m`
      let inline standardDev (m:DM) = DM.StandardDev(m)
  
      /// Variance of matrix `m`
      let inline variance (m:DM) = DM.Variance(m)
  
      /// Shift and scale the elements of matrix `m` to have zero mean and unit variance
      let inline standardize (m:DM) = DM.Standardize(m)
  
      /// Shift and scale the elements of matrix `m` to be in the range [0, 1]
      let inline normalize (m:DM) = DM.Normalize(m)
  
      /// Solve a system of linear equations Ax = b, where the coefficient matrix `m` has general form
      let inline solve (m:DM) (v:DV) = DM.Solve(m, v)
  
      /// Solve a system of linear equations Ax = b, where the coefficient matrix `m` is symmetric
      let inline solveSymmetric (m:DM) (v:DV) = DM.SolveSymmetric(m, v)
  
      /// Sums the elements of matrix `m`
      let inline sum (m:DM) = DM.Sum(m)
  
      /// Trace of matrix `m`
      let inline trace (m:DM) = DM.Trace(m)
  
      /// Append row `v` to matrix `m`
      let inline appendRow (v:DV) (m:DM) = let rows = m |> toRows in Seq.append rows (seq [v]) |> ofRows
  
      /// Prepend row `v` to matrix `m`
      let inline prependRow (v:DV) (m:DM) = let rows = m |> toRows in Seq.append (seq [v]) rows |> ofRows
  
      /// Append column `v` to matrix `m`
      let inline appendCol (v:DV) (m:DM) = let cols = m |> toCols in Seq.append cols (seq [v]) |> ofCols
  
      /// Prepend column `v` to matrix `m`
      let inline prependCol (v:DV) (m:DM) = let cols = m |> toCols in Seq.append (seq [v]) cols |> ofCols
  
      /// Experimental
      let inline toString (m:DM) = m.ToString()
      let inline visualize (m:DM) = m.Visualize()
      let inline visualizeAsDV (m:DM) = DM.ReshapeToDV(m).Visualize()
  
  
  /// D, DV, DM operations (automatically opened)
  [<AutoOpen>]
  module DOps =
  
negative.code: and DV =
      | DV of number[] // Primal
      | DVF of DV * DV * uint32 // Primal, tangent, layer tag
      | DVR of DV * TraceOp * uint32 * int32 // Primal, parent operation, layer tag, unique
  
      interface dobj
  
  