anchor.diff:      
       | Noop
   
  -/// Represents the computed adjoints for reverse AD
  -and Adjoints() = 
  +
  +/// A constraint used to ensure the evaluation stack is only over D, DV or DM
  +and dobj = interface end
  +
  +
  +/// 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()
  +
  +
  +// Scripts for adjusting the adjoint
  +type Delta = 
  +    | X of D
  +    //| XNeg of DeltaV
  +    interface delta
  +and DeltaV = 
  +    | XV of DV
  +    //| XNegV of DeltaV
  +    interface delta
  +and DeltaM = 
  +    | XM of DM
  +    | XNegM of DeltaM
  +    interface delta
  +and delta = interface end
  +
  +/// Represents the computed adjoints for reverse AD.  This is a table indexed by node ID.
  +/// The table is destructively updated as the adjoints are accumulated.
  +type Adjoints() = 
       let dict = Dictionary<int, dobj>()
  +
  +    let rec eval d = 
  +        match d with 
  +        | X d -> d
  +    and evalV d =
  +        match d with 
  +        | XV d -> d
  +    and evalM d = 
  +        match d with 
  +        | XM v -> v
  +        | XNegM v -> -(evalM v)
   
       member internal __.GetD(uniq: int) = dict.[uniq] :?> D
       member internal __.SetD(uniq:int, v:D) = dict.[uniq] <- v
  
positive.code:     let inline log10Val<'T>              = log10ValFloat32
  
  
  /// A constraint used to ensure the evaluation stack is only over D, DV or DM
  type dobj = interface end
  let bxd (x : dobj) = x
  
  /// Scalar numeric type keeping dual numbers for forward mode and adjoints and tapes for reverse mode AD, 
  /// with nesting capability, using tags to avoid perturbation confusion
  
negative.code:                         else pushRec t
                      | _ -> pushRec t
  
                  | (:? DeltaV as delta), (:? DV as d) ->
                      match d with
                      | DVR(_, o, _, uniq) ->
                          let dA = adjoints.ApplyDeltaV(uniq,delta)
                          let fanout = fanouts.[uniq] - 1u
                          fanouts.[uniq] <- fanout
                          // If all incoming parts of the adjoint have been received, then proceed to the parent
                          if fanout = 0u then
                              match o with
                              | Add_DV_DV(a, b) -> pushRec ((bxv dA a) :: (bxv dA b) :: t)
                              | Add_DV_DVCons(a) -> pushRec ((bxv dA a) :: t)
                              | Add_DV_D(a, b) -> pushRec ((bxv dA a) :: (bx (DV.Sum(dA)) b) :: t)
                              | Add_DV_DCons(a) -> pushRec ((bxv dA a) :: t)
                              | Add_DVCons_D(b) -> pushRec ((bx (DV.Sum(dA)) b) :: t)
                              | Sub_DV_DV(a, b) -> pushRec ((bxv dA a) :: (bxv -dA b) :: t)
                              | Sub_DV_DVCons(a) -> pushRec ((bxv dA a) :: t)
                              | Sub_DVCons_DV(a) -> pushRec ((bxv -dA a) :: t)
                              | Sub_DV_D(a, b) -> pushRec ((bxv dA a) :: (bx -(DV.Sum(dA)) b) :: t)
                              | Sub_DV_DCons(a) -> pushRec ((bxv dA a) :: t)
                              | Sub_DVCons_D(b) -> pushRec ((bx -(DV.Sum(dA)) b) :: t)
                              | Sub_D_DV(a, b) -> pushRec ((bx (DV.Sum(dA)) a) :: (bxv -dA b) :: t)
                              | Sub_D_DVCons(a) -> pushRec ((bx (DV.Sum(dA)) a) :: t)
                              | Sub_DCons_DV(b) -> pushRec ((bxv -dA b) :: t)
                              | Mul_Had_DV_DV(a, b) -> pushRec ((bxv (dA .* b.P) a) :: (bxv (dA .* a.P) b) :: t)
                              | Mul_Had_DV_DVCons(a, cons) -> pushRec ((bxv (dA .* cons) a) :: t)
                              | Mul_DV_D(a, b) -> pushRec ((bxv (dA * b.P) a) :: (bx (dA * a.P) b) :: t)
                              | Mul_DV_DCons(a, cons) -> pushRec ((bxv (dA * cons) a) :: t)
                              | Mul_DVCons_D(cons, b) -> pushRec ((bx (dA * cons) b) :: t)
                              | Mul_DM_DV(a, b) -> pushRec ((bxm (dA &* b.P) a) :: (bxv (DM.Transpose(a.P) * dA) b) :: t)
                              | Mul_DM_DVCons(a, cons) -> pushRec ((bxm (dA &* cons) a) :: t)
                              | Mul_DMCons_DV(cons, b) -> pushRec ((bxv (DM.Transpose(cons) * dA) b) :: t)
                              | Mul_DV_DM(a, b) -> pushRec ((bxv (dA * DM.Transpose(b.P)) a) :: (bxm (a.P &* dA) b) :: t)
                              | Mul_DV_DMCons(a, cons) -> pushRec ((bxv (dA * DM.Transpose(cons)) a) :: t)
                              | Mul_DVCons_DM(cons, b) -> pushRec ((bxm (cons &* dA) b) :: t)
                              | Div_Had_DV_DV(a, b) -> pushRec ((bxv (dA ./ b.P) a) :: (bxv (dA .* (-a.P ./ (b.P .* b.P))) b) :: t)
                              | Div_Had_DV_DVCons(a, cons) -> pushRec ((bxv (dA ./ cons) a) :: t)
                              | Div_Had_DVCons_DV(cons, b) -> pushRec ((bxv (dA .* (-cons ./ (b.P .* b.P))) b) :: t)
                              | Div_DV_D(a, b) -> pushRec ((bxv (dA / b.P) a) :: (bx (dA * (-a.P / (b.P * b.P))) b) :: t)
                              | Div_DV_DCons(a, cons) -> pushRec ((bxv (dA / cons) a) :: t)
                              | Div_DVCons_D(cons, b) -> pushRec ((bx (dA * (-cons / (b.P * b.P))) b) :: t)
                              | Div_D_DV(a, b) -> pushRec ((bx (DV.Sum(dA ./ b.P)) a) :: (bxv (dA .* (-a.P / (b.P .* b.P))) b) :: t)
                              | Div_D_DVCons(a, cons) -> pushRec ((bx (DV.Sum(dA ./ cons)) a) :: t)
                              | Div_DCons_DV(cons, b) -> pushRec ((bxv (dA .* (-cons / (b.P .* b.P))) b) :: t)
                              | Pow_DV_DV(a, b) -> pushRec ((bxv (dA .* (a.P ** (b.P - D.One)) .* b.P) a) :: (bxv (dA .* (a.P ** b.P) .* log a.P) b) :: t)
                              | Pow_DV_DVCons(a, cons) -> pushRec ((bxv (dA .* (a.P ** (cons - D.One)) .* cons) a) :: t)
                              | Pow_DVCons_DV(cons, b) -> pushRec ((bxv (dA .* (cons ** b.P) .* log cons) b) :: t)
                              | Atan2_DV_DV(a, b) -> let denom = (a.P .* a.P) + (b.P .* b.P) in pushRec ((bxv (dA .* b.P ./ denom) a) :: (bxv (dA .* (-a.P) ./ denom) b) :: t)
                              | Atan2_DV_DVCons(a, cons) -> pushRec ((bxv (dA .* cons ./ ((a.P .* a.P) + (cons .* cons))) a) :: t)
                              | Atan2_DVCons_DV(cons, b) -> pushRec ((bxv (dA .* (-cons) ./ ((cons .* cons) + (b.P .* b.P))) b) :: t)
                              | Pow_DV_D(a, b) -> pushRec ((bxv (dA .* (a.P ** (b.P - D.One)) * b.P) a) :: (bx (DV.Sum(dA .* (a.P ** b.P) .* log a.P)) b) :: t)
                              | Pow_DV_DCons(a, cons) -> pushRec ((bxv (dA .* (a.P ** (cons - D.One)) * cons) a) :: t)
                              | Pow_DVCons_D(cons, b) -> pushRec ((bx (DV.Sum(dA .* (cons ** b.P) .* log cons)) b) :: t)
                              | Pow_D_DV(a, b) -> pushRec ((bx (DV.Sum(dA .* (DV.Pow(a.P, b.P - D.One)) .* b.P)) a) :: (bxv (dA .* (DV.Pow(a.P, b.P)) * log a.P) b) :: t)
                              | Pow_D_DVCons(a, cons) -> pushRec ((bx (DV.Sum(dA .* (DV.Pow(a.P, cons - D.One)) .* cons)) a) :: t)
                              | Pow_DCons_DV(cons, b) -> pushRec ((bxv (dA .* (DV.Pow(cons, b.P)) * log cons) b) :: t)
                              | Atan2_DV_D(a, b) -> let denom = (a.P .* a.P) + (b.P * b.P) in pushRec ((bxv (dA * b.P ./ denom) a) :: (bx (DV.Sum(dA .* (-a.P) ./ denom)) b) :: t)
                              | Atan2_DV_DCons(a, cons) -> pushRec ((bxv (dA * cons ./ ((a.P .* a.P) + (cons * cons))) a) :: t)
                              | Atan2_DVCons_D(cons, b) -> pushRec ((bx (DV.Sum(dA .* (-cons) ./ ((cons .* cons) + (b.P * b.P)))) b) :: t)
                              | Atan2_D_DV(a, b) -> let denom = (a.P * a.P) + (b.P .* b.P) in pushRec ((bx (DV.Sum(dA .* b.P ./ denom)) a) :: (bxv (dA * (-a.P) ./ denom) b) :: t)
                              | Atan2_D_DVCons(a, cons) -> pushRec ((bx (DV.Sum(dA .* cons ./ ((a.P * a.P) + (cons .* cons)))) a) :: t)
                              | Atan2_DCons_DV(cons, b) -> pushRec ((bxv (dA * (-cons) ./ ((cons * cons) + (b.P .* b.P))) b) :: t)
                              | Log_DV(a) -> pushRec ((bxv (dA ./ a.P) a) :: t)
                              | Log10_DV(a) -> pushRec ((bxv (dA ./ (a.P * N.log10Val)) a) :: t)
                              | Exp_DV(a) -> pushRec ((bxv (dA .* d.P) a) :: t) // d.P = exp a.P
                              | Sin_DV(a) -> pushRec ((bxv (dA .* cos a.P) a) :: t)
                              | Cos_DV(a) -> pushRec ((bxv (-dA .* sin a.P) a) :: t)
                              | Tan_DV(a) -> let seca = D.One / cos a.P in pushRec ((bxv (dA .* seca .* seca) a) :: t)
                              | Neg_DV(a) -> pushRec ((bxv -dA a) :: t)
                              | Sqrt_DV(a) -> pushRec ((bxv (dA ./ (N.two * d.P)) a) :: t) // d.P = sqrt a.P
                              | Sinh_DV(a) -> pushRec ((bxv (dA .* cosh a.P) a) :: t)
                              | Cosh_DV(a) -> pushRec ((bxv (dA .* sinh a.P) a) :: t)
                              | Tanh_DV(a) -> let secha = D.One / cosh a.P in pushRec ((bxv (dA .* secha .* secha) a) :: t)
                              | Asin_DV(a) -> pushRec ((bxv (dA ./ sqrt (D.One - (a.P .* a.P))) a) :: t)
                              | Acos_DV(a) -> pushRec ((bxv (-dA ./ sqrt (D.One - (a.P .* a.P))) a) :: t)
                              | Atan_DV(a) -> pushRec ((bxv (dA ./ (D.One + (a.P .* a.P))) a) :: t)
                              | Abs_DV(a) -> pushRec ((bxv (dA .* DV.Sign a.P) a) :: t)
                              | Sign_DV(a) -> pushRec ((bxv DV.Zero a) :: t)
                              | Floor_DV(a) -> pushRec ((bxv DV.Zero a) :: t)
                              | Ceil_DV(a) -> pushRec ((bxv DV.Zero a) :: t)
                              | Round_DV(a) -> pushRec ((bxv DV.Zero a) :: t)
                              | Make_DV_ofDs(a) -> pushRec (t |> List.append (a |> Array.mapi (fun i v -> (bx dA.[i] v)) |> List.ofArray))
                              | SliceRow_DM(a, i, j) ->
                                  adjoints.[a] <- DM.AddSubMatrix(adjoints.[a], i, j, dA.ToRowDM())
                                  pushRec ((bxm DM.Zero a) :: t)
                              | SliceCol_DM(a, i, j) ->
                                  adjoints.[a] <- DM.AddSubMatrix(adjoints.[a], i, j, dA.ToColDM())
                                  pushRec ((bxm DM.Zero a) :: t)
                              | Solve_DM_DV(a, b) -> let ba = DM.Solve(DM.Transpose(a), dA) in pushRec ((bxm (-ba &* dA) a) :: (bxv (ba) b) :: t)
                              | Solve_DM_DVCons(a, cons) -> let ba = DM.Solve(DM.Transpose(a), dA) in pushRec ((bxm (-ba &* dA) a) :: t)
                              | Solve_DMCons_DV(cons, b) -> let ba = DM.Solve(DM.Transpose(cons), dA) in pushRec ((bxv ba b) :: t)
                              | Append_DV_DV(a, b) ->
                                  adjoints.[a] <- adjoints.[a] + dA.[..(a.Length - 1)]
                                  adjoints.[b] <- adjoints.[b] + dA.[a.Length..]
                                  pushRec ((bxv DV.Zero a) :: (bxv DV.Zero b) :: t)
                              | Append_DV_DVCons(a) ->
                                  adjoints.[a] <- adjoints.[a] + dA.[..(a.Length - 1)]
                                  pushRec ((bxv DV.Zero a) :: t)
                              | Append_DVCons_DV(b) ->
                                  adjoints.[b] <- adjoints.[b] + dA.[(d.Length - b.Length)..]
                                  pushRec ((bxv DV.Zero b) :: t)
                              | Split_DV(a, i) ->
                                  adjoints.[a] <- DV.AddSubVector(adjoints.[a], i, dA)
                                  pushRec ((bxv DV.Zero a) :: t)
                              | AddItem_DV_D(a, i, b) -> pushRec ((bxv dA a) :: (bx (dA.[i]) b) :: t)
                              | AddItem_DV_DCons(a) -> pushRec ((bxv dA a) :: t)
                              | AddItem_DVCons_D(i, b) -> pushRec ((bx dA.[i] b) :: t)
                              | AddSubVector_DV_DV(a, i, b) -> pushRec ((bxv dA a) :: (bxv (dA.[i..(i + b.Length - 1)]) b) :: t)
                              | AddSubVector_DV_DVCons(a) -> pushRec ((bxv dA a) :: t)
                              | AddSubVector_DVCons_DV(i, b) -> pushRec ((bxv (dA.[i..(i + b.Length - 1)]) b) :: t)
                              | ReshapeCopy_DM_DV(a) -> pushRec ((bxm (DV.ReshapeToDM(a.Rows, dA)) a) :: t)
                              | Slice_DV(a, i) ->
                                  adjoints.[a] <- DV.AddSubVector(adjoints.[a], i, dA)
                                  pushRec ((bxv DV.Zero a) :: t)
                              | Diagonal_DM(a) -> 
                                  adjoints.[a] <- DM.AddDiagonal(adjoints.[a], dA)
                                  pushRec ((bxm DM.Zero a) :: t)
                              | ReLU_DV(a) -> pushRec ((bxv (dA .* ((DV.Sign(a.P) + N.one) / N.two)) a) :: t)
                              | Sigmoid_DV(a) -> pushRec ((bxv (dA .* d.P .* (N.one - d.P)) a) :: t) // d.P = DV.Sigmoid(a.P)
                              | _ -> pushRec t
                          else pushRec t
                      | _ -> pushRec t
  
                  | (:? DeltaM as delta), (:? DM as d) ->
                      match d with
                      | DMR(_, o, _, uniq) ->
                          let dA = adjoints.ApplyDeltaM(uniq,delta)
                          let fanout = fanouts.[uniq] - 1u
                          fanouts.[uniq] <- fanout
                          // If all incoming parts of the adjoint have been received, then proceed to the parent
                          if fanout = 0u then
                              match o with
                              | Add_DM_DM(a, b) -> pushRec ((bxm dA a) :: (bxm dA b) :: t)
                              | Add_DM_DMCons(a) -> pushRec ((bxm dA a) :: t)
  
                              // When pushing "-dA" as adjoint increment for b, the operation 
                              //    "b.Adjoint <- -1.0 * dA + b.Adjoint" 
                              // can be performed directly in-place. Instead of pushing a D|DV|DM we should a 
                              // structured expression about how to compute the D|DV|DM which can be interpreted 
                              // to do an in-place update
                              | Sub_DM_DM(a, b) -> pushRec ((bxm dA a) :: (bdm (XNegM (XM dA)) b) :: t)
  
                              // TODO: also avoid the inplace operations in most of the below. 
                              | Sub_DM_DMCons(a) -> pushRec ((bxm dA a) :: t)
                              | Sub_DMCons_DM(a) -> pushRec ((bxm -dA a) :: t)
                              | Mul_DM_DM(a, b) -> pushRec ((bxm (dA * DM.Transpose(b.P)) a) :: (bxm (DM.Transpose(a.P) * dA) b) :: t)
                              | Mul_DM_DMCons(a, cons) -> pushRec ((bxm (dA * DM.Transpose(cons)) a) :: t)
                              | Mul_DMCons_DM(cons, b) -> pushRec ((bxm (DM.Transpose(cons) * dA) b) :: t)
                              | Mul_Had_DM_DM(a, b) -> pushRec ((bxm (dA .* b.P) a) :: (bxm (dA .* a.P) b) :: t)
                              | Mul_Had_DM_DMCons(a, cons) -> pushRec ((bxm (dA .* cons) a) :: t)
                              | Mul_DM_D(a, b) -> pushRec ((bxm (dA * b.P) a) :: (bx (DM.Sum(dA .* a.P)) b) :: t)
                              | Mul_DM_DCons(a, cons) -> pushRec ((bxm (dA * cons) a) :: t)
                              | Mul_DMCons_D(cons, b) -> pushRec ((bx (DM.Sum(dA .* cons)) b) :: t)
                              | Mul_Out_DV_DV(a, b) -> pushRec ((bxv (dA * b.P) a) :: (bxv (DM.Transpose(dA) * a.P) b) :: t)
                              | Mul_Out_DV_DVCons(a, cons) -> pushRec ((bxv (dA * cons) a) :: t)
                              | Mul_Out_DVCons_DV(cons, b) -> pushRec ((bxv (DM.Transpose(dA) * cons) b) :: t)
                              | Div_Had_DM_DM(a, b) -> pushRec ((bxm (dA ./ b.P) a) :: (bxm (dA .* (-a.P ./ (b.P .* b.P))) b) :: t)
                              | Div_Had_DM_DMCons(a, cons) -> pushRec ((bxm (dA ./ cons) a) :: t)
                              | Div_Had_DMCons_DM(cons, b) -> pushRec ((bxm (dA .* (-cons ./ (b.P .* b.P))) b) :: t)
                              | Pow_DM_DM(a, b) -> pushRec ((bxm (dA .* (a.P ** (b.P - D.One)) .* b.P) a) :: (bxm (dA .* (a.P ** b.P) .* log a.P) b) :: t)
                              | Pow_DM_DMCons(a, cons) -> pushRec ((bxm (dA .* (a.P ** (cons - D.One)) .* cons) a) :: t)
                              | Pow_DMCons_DM(cons, b) -> pushRec ((bxm (dA .* (cons ** b.P) .* log cons) b) :: t)
                              | Atan2_DM_DM(a, b) -> let denom = (a.P .* a.P) + (b.P .* b.P) in pushRec ((bxm (dA .* b.P ./ denom) a) :: (bxm (dA .* (-a.P) ./ denom) b) :: t)
                              | Atan2_DM_DMCons(a, cons) -> pushRec ((bxm (dA .* cons ./ ((a.P .* a.P) + (cons .* cons))) a) :: t)
                              | Atan2_DMCons_DM(cons, b) -> pushRec ((bxm (dA .* (-cons) ./ ((cons .* cons) + (b.P .* b.P))) b) :: t)
                              | Add_DM_D(a, b) -> pushRec ((bxm dA a) :: (bx (DM.Sum(dA)) b) :: t)
                              | Add_DM_DCons(a) -> pushRec ((bxm dA a) :: t)
                              | Add_DMCons_D(b) -> pushRec ((bx (DM.Sum(dA)) b) :: t)
                              | Add_DMCols_DV(a, b) ->
                                  dA.GetCols() |> Seq.iter (fun v -> adjoints.[b] <- adjoints.[b] + v)
                                  pushRec ((bxm dA a) :: (bxv DV.Zero b) :: t)
                              | Add_DMCols_DVCons(a) ->
                                  pushRec ((bxm dA a) :: t)
                              | Add_DMColsCons_DV(b) ->
                                  dA.GetCols() |> Seq.iter (fun v -> adjoints.[b] <- adjoints.[b] + v)
                                  pushRec ((bxv DV.Zero b) :: t)
                              | Sub_DM_D(a, b) -> pushRec ((bxm dA a) :: (bx -(DM.Sum(dA)) b) :: t)
                              | Sub_DM_DCons(a) -> pushRec ((bxm dA a) :: t)
                              | Sub_DMCons_D(b) -> pushRec ((bx -(DM.Sum(dA)) b) :: t)
                              | Sub_D_DM(a, b) -> pushRec ((bx (DM.Sum(dA)) a) :: (bxm -dA b) :: t)
                              | Sub_D_DMCons(a) -> pushRec ((bx (DM.Sum(dA)) a) :: t)
                              | Sub_DCons_DM(b) -> pushRec ((bxm -dA b) :: t)
                              | Div_DM_D(a, b) -> pushRec ((bxm (dA / b.P) a) :: (bx (DM.Sum (dA .* (-a.P / b.P * b.P))) b) :: t)
                              | Div_DM_DCons(a, cons) -> pushRec ((bxm (dA / cons) a) :: t)
                              | Div_DMCons_D(cons, b) -> pushRec ((bx (DM.Sum (dA .* (-cons / (b.P * b.P)))) b) :: t)
                              | Div_D_DM(a, b) -> pushRec ((bx (DM.Sum(dA ./ b.P)) a) :: (bxm (dA .* (-a.P / (b.P .* b.P))) b) :: t)
                              | Div_D_DMCons(a, cons) -> pushRec ((bx (DM.Sum(dA ./ cons)) a) :: t)
                              | Div_DCons_DM(cons, b) -> pushRec ((bxm (dA .* (-cons / (b.P .* b.P))) b) :: t)
                              | Pow_DM_D(a, b) -> pushRec ((bxm (dA .* (a.P ** (b.P - D.One)) * b.P) a) :: (bx (DM.Sum(dA .* (a.P ** b.P) .* log a.P)) b) :: t)
                              | Pow_DM_DCons(a, cons) -> pushRec ((bxm (dA .* (a.P ** (cons - D.One)) * cons) a) :: t)
                              | Pow_DMCons_D(cons, b) -> pushRec ((bx (DM.Sum(dA .* (cons ** b.P) .* log cons)) b) :: t)
                              | Pow_D_DM(a, b) -> pushRec ((bx (DM.Sum(dA .* (DM.Pow(a.P, b.P - D.One)) .* b.P)) a) :: (bxm (dA .* (DM.Pow(a.P, b.P)) * log a.P) b) :: t)
                              | Pow_D_DMCons(a, cons) -> pushRec ((bx (DM.Sum(dA .* (DM.Pow(a.P, cons - D.One)) .* cons)) a) :: t)
                              | Pow_DCons_DM(cons, b) -> pushRec ((bxm (dA .* (DM.Pow(cons, b.P)) * log cons) b) :: t)
                              | Atan2_DM_D(a, b) -> let denom = (a.P .* a.P) + (b.P * b.P) in pushRec ((bxm (dA * b.P ./ denom) a) :: (bx (DM.Sum(dA .* (-a.P) ./ denom)) b) :: t)
                              | Atan2_DM_DCons(a, cons) -> pushRec ((bxm (dA * cons ./ ((a.P .* a.P) + (cons * cons))) a) :: t)
                              | Atan2_DMCons_D(cons, b) ->pushRec ((bx (DM.Sum(dA .* (-cons) ./ ((cons .* cons) + (b.P * b.P)))) b) :: t)
                              | Atan2_D_DM(a, b) -> let denom = (a.P * a.P) + (b.P .* b.P) in pushRec ((bx (DM.Sum(dA .* b.P ./ denom)) a) :: (bxm (dA * (-a.P) ./ denom) b) :: t)
                              | Atan2_D_DMCons(a, cons) -> pushRec ((bx (DM.Sum(dA .* cons ./ ((a.P * a.P) + (cons .* cons)))) a) :: t)
                              | Atan2_DCons_DM(cons, b) -> pushRec ((bxm (dA * (-cons) ./ ((cons * cons) + (b.P .* b.P))) b) :: t)
                              | Log_DM(a) -> pushRec ((bxm (dA ./ a.P) a) :: t)
                              | Log10_DM(a) -> pushRec ((bxm (dA ./ (a.P * N.log10Val)) a) :: t)
                              | Exp_DM(a) -> pushRec ((bxm (dA .* d.P) a) :: t) // d.P = exp a.P
                              | Sin_DM(a) -> pushRec ((bxm (dA .* cos a.P) a) :: t)
                              | Cos_DM(a) -> pushRec ((bxm (-dA .* sin a.P) a) :: t)
                              | Tan_DM(a) -> let seca = D.One / cos a.P in pushRec ((bxm (dA .* seca .* seca) a) :: t)
                              | Neg_DM(a) -> pushRec ((bxm -dA a) :: t)
                              | Sqrt_DM(a) -> pushRec ((bxm (dA ./ (N.two * d.P)) a) :: t) // d.P = sqrt a.P
                              | Sinh_DM(a) -> pushRec ((bxm (dA .* cosh a.P) a) :: t)
                              | Cosh_DM(a) -> pushRec ((bxm (dA .* sinh a.P) a) :: t)
                              | Tanh_DM(a) -> let secha = D.One / cosh a.P in pushRec ((bxm (dA .* secha .* secha) a) :: t)
                              | Asin_DM(a) -> pushRec ((bxm (dA ./ sqrt (D.One - (a.P .* a.P))) a) :: t)
                              | Acos_DM(a) -> pushRec ((bxm (-dA ./ sqrt (D.One - (a.P .* a.P))) a) :: t)
                              | Atan_DM(a) -> pushRec ((bxm (dA ./ (D.One + (a.P .* a.P))) a) :: t)
                              | Abs_DM(a) -> pushRec ((bxm (dA .* DM.Sign a.P) a) :: t)
                              | Sign_DM(a) -> pushRec ((bxm DM.Zero a) :: t)
                              | Floor_DM(a) -> pushRec ((bxm DM.Zero a) :: t)
                              | Ceil_DM(a) -> pushRec ((bxm DM.Zero a) :: t)
                              | Round_DM(a) -> pushRec ((bxm DM.Zero a) :: t)
                              | Transpose_DM(a) -> pushRec ((bxm (DM.Transpose(dA)) a) :: t)
                              | Make_DM_ofDs(a) -> pushRec (t |> List.append (List.map2 (fun v dd -> (bx v dd)) (dA |> DM.toDV |> DV.toArray |> Array.toList) (a |> Array2D.toArray |> List.ofArray)))
                              | Make_DMRows_ofDV(a) ->
                                  dA.GetRows() |> Seq.iter (fun v -> adjoints.[a] <- adjoints.[a] + v)
                                  pushRec ((bxv DV.Zero a) :: t)
                              | Make_DMCols_ofDV(a) ->
                                  dA.GetCols() |> Seq.iter (fun v -> adjoints.[a] <- adjoints.[a] + v)
                                  pushRec ((bxv DV.Zero a) :: t)
                              | Make_DMRows_ofDVs(a) -> pushRec (t |> List.append (a |> List.ofArray |> List.mapi (fun i v -> (bxv dA.[i, *] v))))
                              | AddItem_DM_D(a, i, j, b) -> pushRec ((bxm dA a) :: (bx (dA.[i, j]) b) :: t)
                              | AddItem_DM_DCons(a) -> pushRec ((bxm dA a) :: t)
                              | AddItem_DMCons_D(i, j, b) -> pushRec ((bx dA.[i, j] b) :: t)
                              | AddSubMatrix_DM_DM(a, i, j, b) -> pushRec ((bxm dA a) :: (bxm (dA.[i..(i + b.Rows - 1), j..(j + b.Cols - 1)]) b) :: t)
                              | AddSubMatrix_DM_DMCons(a) -> pushRec ((bxm dA a) :: t)
                              | AddSubMatrix_DMCons_DM(i, j, b) -> pushRec ((bxm (dA.[i..(i + b.Rows - 1), j..(j + b.Cols - 1)]) b) :: t)
                              | Slice_DM(a, i, j) ->
                                  adjoints.[a] <- DM.AddSubMatrix(adjoints.[a], i, j, dA)
                                  pushRec ((bxm DM.Zero a) :: t)
                              | RowMatrix_DV(a) -> pushRec ((bxv (dA.[0, *]) a) :: t)
                              | AddDiagonal_DM_DV(a, b) -> pushRec ((bxm dA a) :: (bxv (DM.Diagonal(dA)) b) :: t)
                              | AddDiagonal_DM_DVCons(a) -> pushRec ((bxm dA a) :: t)
                              | AddDiagonal_DMCons_DV(b) -> pushRec ((bxv (DM.Diagonal(dA)) b) :: t)
                              | ReshapeCopy_DV_DM(a) -> pushRec ((bxv (DM.ReshapeToDV(dA)) a) :: t)
                              | Inverse_DM(a) -> let dpt = DM.Transpose(d.P) in pushRec ((bxm (-dpt * dA * dpt) a) :: t) // d.P = DM.Inverse(a.P)
                              | ReLU_DM(a) -> pushRec ((bxm (dA .* ((DM.Sign(a.P) + N.one) / N.two)) a) :: t)
                              | Sigmoid_DM(a) -> pushRec ((bxm (dA .* d.P .* (N.one - d.P)) a) :: t) // d.P = DM.Sigmoid(a.P)
                              | _ -> pushRec t
                          else pushRec t
                      | _ -> pushRec t
                  | _ -> pushRec t
          let initialv = 
              match v with 
              | :? D as v -> bxdelta (X v)
              | :? DV as v -> bxdelta (XV v)
              | :? DM as v -> bxdelta (XM v)
              | _ -> failwith "invalid dobj"
          pushRec [(initialv, d)]
  
  /// Forward and reverse differentiation operations module (automatically opened)
  [<AutoOpen>]
  