anchor.diff:          let inline r_d_d(a, b) = Add_DM_DM(a, b)
           let inline r_d_c(a, b) = Add_DM_DMCons(a)
           let inline r_c_d(a, b) = Add_DM_DMCons(b)
  +        DM.Op_DM_DM_DM (a, b, ff, fd, df_da, df_db, df_dab, r_d_d, r_d_c, r_c_d)
  +
  +    /// Element-wise addition of `a` and `b`, potentially destructive of the storage of raw matrices in 'b'
  +    static member internal Add_M_M_Inplace (a:DM, b:DM) =
  +        let inline ff(a:number[,], b:number[,]) = Backend.Add_MulM_M_Inplace(N.one, a, b); b
  +        let inline fd(a:DM, b:DM) = a + b
  +        let inline df_da(cp:DM, ap:DM, at:DM) = at
  +        let inline df_db(cp:DM, bp:DM, bt:DM) = bt
  +        let inline df_dab(cp:DM, ap:DM, at:DM, bp:DM, bt:DM) = at + bt
  +        let inline r_d_d(a:DM, b:DM) = Add_DM_DM(a, b)
  +        let inline r_d_c(a:DM, b:DM) = Add_DM_DMCons(a)
  +        let inline r_c_d(a:DM, b:DM) = Add_DM_DMCons(b)
           DM.Op_DM_DM_DM (a, b, ff, fd, df_da, df_db, df_dab, r_d_d, r_d_c, r_c_d)
   
       /// Element-wise subtraction of `a` and `b`
  
positive.code:     
      | Noop
  
  /// Represents the computed adjoints for reverse AD
  and Adjoints() = 
      let dict = Dictionary<int, dobj>()
  
      member internal __.GetD(uniq: int) = dict.[uniq] :?> D
      member internal __.SetD(uniq:int, v:D) = dict.[uniq] <- v
  
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>]
  