text:
# Local Projection (LP) stabilization

import SummationByParts: UnaryFunctor, Add, Subtract

#=
"""
  Computes the LPS matrix P = I - L*L^T*H
"""
function getLPSMatrix(sbp::AbstractOperator)

  #TODO: do this correctly
  P = zeros(sbp.numnodes, sbp.numnodes)
  for i=1:length(P)
    P[i] = 0.01*i
  end

  return P
end
=#
function getLPSMatrix(sbp::AbstractOperator)

  diss = zeros(sbp.numnodes, sbp.numnodes)

  x = calcnodes(sbp)
  V = SummationByParts.calcvandermondproriol(x.', sbp.degree)

  # Minv = diagm(1./diag(V.'*diagm(sbp.w)*V))
  # Minv is the idenity, since the V's are orthogonal in L2
  diss[:,:] = eye(sbp.numnodes) - (V*V.')*diagm(sbp.w)

  return diss
end

"""
  Main function for applying LPS stabilization

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
"""
function applyLPStab(mesh, sbp, eqn, opts)

  lps_data = eqn.params.lps_data
  entropy_vars = lps_data.entropy_vars

  _applyLPStab(mesh, sbp, eqn, opts, lps_data, entropy_vars)

  return nothing
end

@noinline function _applyLPStab(mesh, sbp, eqn, opts, data::LPSData{Tsol},
                       entropy_vars::AbstractVariables) where {Tsol}

  w_i = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerElement)
  for i=1:mesh.numEl
    q_i = ro_sview(eqn.q, :, :, i)
    dxidx = ro_sview(mesh.dxidx, :, :, :, i)
    res_i = sview(eqn.res, :, :, i)

    # convert to entropy variables
    @simd for j=1:mesh.numNodesPerElement
      q_j = ro_sview(q_i, :, j)
      w_j = sview(w_i, :, j)
      convertToEntropy(entropy_vars, eqn.params, q_j, w_j)
    end

    applyLPSKernel(eqn.params, data, sbp, entropy_vars, q_i, w_i, dxidx,res_i, i)
  end

  return nothing
end

"""
  Compute res += op( P^T * H * A * P * u_in)  for each component of the system.

  **Inputs**

   * params: ParamType
   * data: [`LPSData`](@ref)
   * sbp: SBP operator
   * q: conservative variables for the element, `numDofPerNode` x
        `numNodesPerElement`
   * u_in: values to multiply against, same size as `q`
   * op: UnaryFunctor specifying whether to add or subtract from `res`.  Default
         subtract.

  **Inputs/Outputs**

   * res: array to update with result, same size as `q`

  **Aliasing**

  Note that `q` and `u_in` can alias
"""
function applyLPSKernel(params::ParamType{Tdim}, data::LPSData, sbp::AbstractOperator,
                        entropy_vars::AbstractVariables,
                        q::AbstractMatrix{Tsol}, u_in::AbstractMatrix,
                        dxidx::Abstract3DArray{Tmsh},
                        res::AbstractMatrix{Tres}, elnum, op::UnaryFunctor=Subtract()
                        ) where {Tsol, Tres, Tmsh, Tdim}
  # For scalar equations, the operator is applied -epsilon * P^T M P * u
  # For vector equations, P needs to be applied to all equations as once:
  # utmp^T = P*u^T
  # Instead, work with utmp = (P*u^T)^T = u*P^T

  @unpack data t1 t2 A0 dir
  numDofPerNode, numNodesPerElement = size(q)

  # apply P
  smallmatmatT!(u_in, data.P, t1)

  # apply mass matrix and scaling
  @simd for i=1:numNodesPerElement
    q_i = ro_sview(q, :, i)
    getA0(entropy_vars, params, q_i, A0)
    t1_i = ro_sview(t1, :, i)
    t2_i = sview(t2, :, i)

    # scaling by wave speed
    lambda = zero(Tres)
    @simd for d1=1:Tdim
      @simd for d2=1:Tdim
        dir[d2] = dxidx[d1, d2, i]
      end
      lambda += getLambdaMax(params, q_i, dir)
    end
    fac = lambda*sbp.w[i]/Tdim

    # multiply by A0 and fac at the same time
    smallmatvec_kernel!(A0, t1_i, t2_i, fac, 0)
  end

  # apply P^T
  smallmatmat!(t2, data.P, t1)
  @simd for i=1:numNodesPerElement
    @simd for j=1:numDofPerNode
      res[j, i] += op(data.alpha*t1[j, i])
    end
  end

  return nothing
end

#------------------------------------------------------------------------------
# Differentiated version

"""
  Differentiated version of `applyLPStab`
"""
function applyLPStab_diff(mesh, sbp, eqn, opts, assem::AssembleElementData)

  lps_data = eqn.params.lps_data
  entropy_vars = lps_data.entropy_vars

  _applyLPStab_diff(mesh, sbp, eqn, opts, lps_data, entropy_vars, assem)

  return nothing
end

@noinline function _applyLPStab_diff(mesh::AbstractMesh, sbp::AbstractOperator,
                                 eqn::EulerData, opts,
                                 data::LPSData{Tsol, Tres},
                                 entropy_vars::AbstractVariables,
                                 assem::AssembleElementData) where {Tsol, Tres}

  res_jac = eqn.params.calc_volume_integrals_data.res_jac

  w_in = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerElement)
  w_jac = zeros(Tsol, mesh.numDofPerNode, mesh.numDofPerNode,
                      mesh.numNodesPerElement)
  for i=1:mesh.numEl
    q_i = ro_sview(eqn.q, :, :, i)
    dxidx = ro_sview(mesh.dxidx, :, :, :, i)
    # convert to entropy variables
    @simd for j=1:mesh.numNodesPerElement
      q_j = ro_sview(q_i, :, j)
      w_j = sview(w_in, :, j)
      w_jac_j = sview(w_jac, :, :, j)
      convertToEntropy(entropy_vars, eqn.params, q_j, w_j)
      getA0inv(entropy_vars, eqn.params, q_j, w_jac_j)
    end

    fill!(res_jac, 0)
    applyLPSKernel_diff(eqn.params, data, sbp, entropy_vars, q_i, w_in, w_jac, dxidx,
                        res_jac, i)

    assembleElement(assem, mesh, i, res_jac)
  end

  return nothing
end

"""
  Differentiated version of `applyEntropyKernel`

  **Inputs**

   * params
   * data: [`LPSData`](@ref)
   * sbp
   * q
   * u_in_jac: 3D jacobian of `u_in` wrt `q`, `numDofPerNode` x `numDofPerNode`
               x `numNodesPerElement`
   * op

  **Inputs/Outputs**

   * res_jac: 4D jacobian to be updated with the result.
"""
function applyLPSKernel_diff(params::ParamType{Tdim}, data::LPSData,
                        sbp::AbstractOperator, entropy_vars::AbstractVariables,
                        q_el::AbstractMatrix{Tsol}, u_in::AbstractMatrix,
                        u_in_jac::Abstract3DArray,
                        dxidx::Abstract3DArray{Tmsh},
                        res_jac::Abstract4DArray{Tres}, elnum, op::UnaryFunctor=Subtract()
                        ) where {Tsol, Tres, Tmsh, Tdim}

  numDofPerNode, numNodesPerElement = size(q_el)

  @unpack data t1 t1_jac t2 t2_jac t3_jac A0 A0_diff q_dot lambdas lambdas_jac
  @unpack data dir lambda_dot_p

  # compute required primal quantities (t1 and t2)
  smallmatmatT!(u_in, data.P, t1)

  @simd for i=1:numNodesPerElement
    q_i = ro_sview(q_el, :, i)
    getA0(entropy_vars, params, q_i, A0)
    t1_i = ro_sview(t1, :, i)
    t2_i = sview(t2, :, i)
    fac = 1

    # multiply by A0 and fac at the same time
    smallmatvec_kernel!(A0, t1_i, t2_i, fac, 0)
  end

  # compute lambda and its derivative at all volume nodes
  fill!(lambdas, 0); fill!(lambdas_jac, 0)
  @simd for p=1:numNodesPerElement
    q_p = ro_sview(q_el, :, p)
    @simd for d1=1:Tdim
      @simd for d2=1:Tdim
        dir[d2] = dxidx[d1, d2, p]
      end
      lambdas[p] += getLambdaMax_diff(params, q_p, dir, lambda_dot_p)
      @simd for i=1:numDofPerNode
        lambdas_jac[i, p] += lambda_dot_p[i]
      end
    end
    lambdas[p] /= Tdim
    @simd for i=1:numDofPerNode
      lambdas_jac[i, p] /= Tdim
    end
  end

  #--------------------------
  # compute derivative

  @simd for q=1:numNodesPerElement

    # t1
    @simd for p=1:numNodesPerElement
      @simd for j=1:numDofPerNode
        @simd for i=1:numDofPerNode
          t1_jac[i, j, p] = data.P[p, q]*u_in_jac[i, j, q]
        end
      end
    end

    # A0 * t1_dot
    @simd for p=1:numNodesPerElement
      q_p = ro_sview(q_el, :, p)
      getA0(entropy_vars, params, q_p, A0)
      @simd for j=1:numDofPerNode
        @simd for i=1:numDofPerNode
          t2_jac[i, j, p] = 0
          @simd for k=1:numDofPerNode
            # this should be A0[i, k], but A0 is symmetric
            t2_jac[i, j, p] += A0[k, i]*t1_jac[k, j, p]
          end
        end
      end
    end  # end p

    # contribution of d A0/dq * t1
    q_q = ro_sview(q_el, :, q)
    getA0_diff(entropy_vars, params, q_q, q_dot, A0, A0_diff)
    @simd for j=1:numDofPerNode
      @simd for i=1:numDofPerNode
        @simd for k=1:numDofPerNode
          # A0 is symmetric, so d A0[i, k]/dq = d A0[k, i]/dq
          t2_jac[i, j, q] += A0_diff[k, i, j]*t1[k, q]
        end
      end
    end

    # multiply by lambda and H
    @simd for p=1:numNodesPerElement
      fac = lambdas[p]*sbp.w[p]
      @simd for j=1:numDofPerNode
        @simd for i=1:numDofPerNode
          t3_jac[i, j, p] = fac*t2_jac[i, j, p]
        end
      end
    end  # end p

    @simd for i=1:numDofPerNode
      @simd for j=1:numDofPerNode
        t3_jac[i, j, q] += sbp.w[q]*lambdas_jac[j, q]*t2[i, q]
      end
    end

    # multiply by P^T
    @simd for p=1:numNodesPerElement
      @simd for k=1:numNodesPerElement
        @simd for j=1:numDofPerNode
          @simd for i=1:numDofPerNode
            res_jac[i, j, p, q] += op(data.alpha*data.P[k, p]*t3_jac[i, j, k])
          end
        end
      end
    end

  end  # end q

  return nothing
end

#------------------------------------------------------------------------------
# revq

"""
  Main function for applying LPS stabilization revq

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
"""
function applyLPStab_revq(mesh, sbp, eqn, opts)

  lps_data = eqn.params.lps_data
  entropy_vars = lps_data.entropy_vars

  _applyLPStab_revq(mesh, sbp, eqn, opts, lps_data, entropy_vars)

  return nothing
end

@noinline function _applyLPStab_revq(mesh, sbp, eqn, opts, data::LPSData{Tsol, Tres},
                       entropy_vars::AbstractVariables) where {Tsol, Tres}

  w_i = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerElement)
  w_i_bar = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerElement)
  A0inv = zeros(Tsol, mesh.numDofPerNode, mesh.numDofPerNode)
  for i=1:mesh.numEl
    q_i = ro_sview(eqn.q, :, :, i); q_i_bar = sview(eqn.q_bar, :, :, i)
    dxidx = ro_sview(mesh.dxidx, :, :, :, i)
    res_i = sview(eqn.res, :, :, i); res_i_bar = sview(eqn.res_bar, :, :, i)

    # convert to entropy variables
    @simd for j=1:mesh.numNodesPerElement
      q_j = ro_sview(q_i, :, j)
      w_j = sview(w_i, :, j)
      convertToEntropy(entropy_vars, eqn.params, q_j, w_j)
    end

    fill!(w_i_bar, 0)
    applyLPSKernel_revq(eqn.params, data, sbp, entropy_vars, q_i, q_i_bar,
                        w_i, w_i_bar, dxidx, res_i, res_i_bar, i)

    for j=1:mesh.numNodesPerElement
      q_j = ro_sview(q_i, :, j)
      getA0inv(entropy_vars, eqn.params, q_j, A0inv)
      smallmatvec_kernel!(A0inv, sview(w_i_bar, :, j), sview(q_i_bar, :, j), 1, 1)
    end
  end

  return nothing
end

function applyLPSKernel_revq(params::ParamType{Tdim}, data::LPSData,
                        sbp::AbstractOperator,
                        entropy_vars::AbstractVariables,
                        q::AbstractMatrix{Tsol}, q_bar::AbstractMatrix,
                        u_in::AbstractMatrix, u_in_bar::AbstractMatrix,
                        dxidx::Abstract3DArray{Tmsh},
                        res::AbstractMatrix{Tres}, res_bar::AbstractMatrix,
                        elnum, op::UnaryFunctor=Subtract()
                        ) where {Tsol, Tres, Tmsh, Tdim}
  # For scalar equations, the operator is applied -epsilon * P^T M P * u
  # For vector equations, P needs to be applied to all equations as once:
  # utmp^T = P*u^T
  # Instead, work with utmp = (P*u^T)^T = u*P^T

  @unpack data t1 t2 A0 dir
  numDofPerNode, numNodesPerElement = size(q)

  # apply P
  smallmatmatT!(u_in, data.P, t1)
#=
  # apply mass matrix and scaling
  @simd for i=1:numNodesPerElement
    q_i = ro_sview(q, :, i)
    getA0(entropy_vars, params, q_i, A0)
    t1_i = ro_sview(t1, :, i)
    t2_i = sview(t2, :, i)

    # scaling by wave speed
    lambda = zero(Tres)
    @simd for d1=1:Tdim
      @simd for d2=1:Tdim
        dir[d2] = dxidx[d1, d2, i]
      end
      lambda += getLambdaMax(params, q_i, dir)
    end
    fac = lambda*sbp.w[i]/Tdim

    # multiply by A0 and fac at the same time
    smallmatvec_kernel!(A0, t1_i, t2_i, fac, 0)
  end
=#
#=
  # apply P^T
  smallmatmat!(t2, data.P, t1)
  @simd for i=1:numNodesPerElement
    @simd for j=1:numDofPerNode
      res[j, i] += op(data.alpha*t1[j, i])
    end
  end
=#
  #--------------------------
  # reverse sweep
  @unpack data t1_bar t2_bar tmp tmp_bar A0_bar

  for i=1:numNodesPerElement
    for j=1:numDofPerNode
      t1_bar[j, i] = op(data.alpha*res_bar[j, i])
    end
  end

  smallmatmatT!(t1_bar, data.P, t2_bar)

  fill!(t1_bar, 0)

  for i=1:numNodesPerElement

    # need fac, tmp, A0
    q_i = ro_sview(q, :, i); q_i_bar = sview(q_bar, :, i)
    t1_i = ro_sview(t1, :, i)

    getA0(entropy_vars, params, q_i, A0)

    lambda = zero(Tres)
    @simd for d1=1:Tdim
      @simd for d2=1:Tdim
        dir[d2] = dxidx[d1, d2, i]
      end
      lambda += getLambdaMax(params, q_i, dir)
    end
    fac = lambda*sbp.w[i]/Tdim

    smallmatvec_kernel!(A0, t1_i, tmp, 1, 0)

    # reverse mode
    fill!(tmp_bar, 0)
    fac_bar = zero(Tres)
    for j=1:numDofPerNode
      tmp_bar[j] = fac*t2_bar[j, i]
      fac_bar   += tmp[j]*t2_bar[j, i]
    end

    t1_i_bar = sview(t1_bar, :, i)
    smallmatvec!(A0, tmp_bar, t1_i_bar)
    for j=1:numDofPerNode
      for k=1:numDofPerNode
        A0_bar[k, j] = tmp_bar[k]*t1[j, i]
      end
    end

    lambda_bar = sbp.w[i]*fac_bar/Tdim
    @simd for d1=1:Tdim
      @simd for d2=1:Tdim
        dir[d2] = dxidx[d1, d2, i]
      end
      getLambdaMax_revq(params, q_i, q_i_bar, dir, lambda_bar)
    end

    getA0_revq(entropy_vars, params, q_i, q_i_bar, A0, A0_bar)
  end

  smallmatmat!(t1_bar, data.P, u_in_bar)

  return nothing
end

#------------------------------------------------------------------------------
# revm

"""
  Main function for applying LPS stabilization revm

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
"""
function applyLPStab_revm(mesh, sbp, eqn, opts)

  lps_data = eqn.params.lps_data
  entropy_vars = lps_data.entropy_vars

  _applyLPStab_revm(mesh, sbp, eqn, opts, lps_data, entropy_vars)

  return nothing
end

@noinline function _applyLPStab_revm(mesh, sbp, eqn, opts, data::LPSData{Tsol},
                       entropy_vars::AbstractVariables) where {Tsol}

  w_i = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerElement)
  for i=1:mesh.numEl
    q_i = ro_sview(eqn.q, :, :, i)
    dxidx = ro_sview(mesh.dxidx, :, :, :, i);
    dxidx_bar = sview(mesh.dxidx_bar, :, :, :, i)
    res_i = sview(eqn.res, :, :, i); res_i_bar = sview(eqn.res_bar, :, :, i)

    # convert to entropy variables
    @simd for j=1:mesh.numNodesPerElement
      q_j = ro_sview(q_i, :, j)
      w_j = sview(w_i, :, j)
      convertToEntropy(entropy_vars, eqn.params, q_j, w_j)
    end

    applyLPSKernel_revm(eqn.params, data, sbp, entropy_vars, q_i, w_i,
                        dxidx, dxidx_bar, res_i, res_i_bar, i)
  end

  return nothing
end

function applyLPSKernel_revm(params::ParamType{Tdim}, data::LPSData,
                      sbp::AbstractOperator,
                      entropy_vars::AbstractVariables,
                      q::AbstractMatrix{Tsol}, u_in::AbstractMatrix,
                      dxidx::Abstract3DArray{Tmsh}, dxidx_bar::Abstract3DArray,
                      res::AbstractMatrix{Tres}, res_bar::AbstractMatrix,
                      elnum, op::UnaryFunctor=Subtract()
                        ) where {Tsol, Tres, Tmsh, Tdim}
  # For scalar equations, the operator is applied -epsilon * P^T M P * u
  # For vector equations, P needs to be applied to all equations as once:
  # utmp^T = P*u^T
  # Instead, work with utmp = (P*u^T)^T = u*P^T

  @unpack data t1 t2 A0 dir
  numDofPerNode, numNodesPerElement = size(q)

  # apply P
  smallmatmatT!(u_in, data.P, t1)
#=
  # apply mass matrix and scaling
  @simd for i=1:numNodesPerElement
    q_i = ro_sview(q, :, i)
    getA0(entropy_vars, params, q_i, A0)
    t1_i = ro_sview(t1, :, i)
    t2_i = sview(t2, :, i)

    # scaling by wave speed
    lambda = zero(Tres)
    @simd for d1=1:Tdim
      @simd for d2=1:Tdim
        dir[d2] = dxidx[d1, d2, i]
      end
      lambda += getLambdaMax(params, q_i, dir)
    end
    fac = lambda*sbp.w[i]/Tdim

    # multiply by A0 and fac at the same time
    smallmatvec_kernel!(A0, t1_i, t2_i, fac, 0)
  end
=#
#=
  # apply P^T
  smallmatmat!(t2, data.P, t1)
  @simd for i=1:numNodesPerElement
    @simd for j=1:numDofPerNode
      res[j, i] += op(data.alpha*t1[j, i])
    end
  end
=#
  #--------------------------
  # reverse sweep
  @unpack data t1_bar t2_bar tmp dir_bar

  for i=1:numNodesPerElement
    for j=1:numDofPerNode
      t1_bar[j, i] = op(data.alpha*res_bar[j, i])
    end
  end

  smallmatmatT!(t1_bar, data.P, t2_bar)

  for i=1:numNodesPerElement

    # need fac, tmp, A0
    q_i = ro_sview(q, :, i)
    t1_i = ro_sview(t1, :, i)

    getA0(entropy_vars, params, q_i, A0)
    smallmatvec_kernel!(A0, t1_i, tmp, 1, 0)

    # reverse mode
    fac_bar = zero(Tres)
    for j=1:numDofPerNode
      fac_bar   += tmp[j]*t2_bar[j, i]
    end

    lambda_bar = sbp.w[i]*fac_bar/Tdim
    @simd for d1=1:Tdim
      @simd for d2=1:Tdim
        dir[d2] = dxidx[d1, d2, i]
        dir_bar[d2] = 0
      end
      getLambdaMax_revm(params, q_i, dir, dir_bar, lambda_bar)

      @simd for d2=1:Tdim
        dxidx_bar[d1, d2, i] += dir_bar[d2]
      end
    end
  end

  return nothing
end


text:
# Julia wrapper for header: IviFgen.h
# Automatically generated using Clang.jl

function IviFgen_init(logicalName, IDQuery, resetDevice, vi)
    ccall((:IviFgen_init, libnifgen), ViStatus, (ViRsrc, ViBoolean, ViBoolean, Ptr{ViSession}), logicalName, IDQuery, resetDevice, vi)
end

function IviFgen_close(vi)
    ccall((:IviFgen_close, libnifgen), ViStatus, (ViSession,), vi)
end

function IviFgen_reset(vi)
    ccall((:IviFgen_reset, libnifgen), ViStatus, (ViSession,), vi)
end

function IviFgen_self_test(vi, selfTestResult, selfTestMessage)
    ccall((:IviFgen_self_test, libnifgen), ViStatus, (ViSession, Ptr{ViInt16}, Ptr{ViChar}), vi, selfTestResult, selfTestMessage)
end

function IviFgen_error_query(vi, errorCode, errorMessage)
    ccall((:IviFgen_error_query, libnifgen), ViStatus, (ViSession, Ptr{ViInt32}, Ptr{ViChar}), vi, errorCode, errorMessage)
end

function IviFgen_error_message(vi, statusCode, message)
    ccall((:IviFgen_error_message, libnifgen), ViStatus, (ViSession, ViStatus, Ptr{ViChar}), vi, statusCode, message)
end

function IviFgen_revision_query(vi, driverRev, instrRev)
    ccall((:IviFgen_revision_query, libnifgen), ViStatus, (ViSession, Ptr{ViChar}, Ptr{ViChar}), vi, driverRev, instrRev)
end

function IviFgen_InvalidateAllAttributes(vi)
    ccall((:IviFgen_InvalidateAllAttributes, libnifgen), ViStatus, (ViSession,), vi)
end

function IviFgen_ResetWithDefaults(vi)
    ccall((:IviFgen_ResetWithDefaults, libnifgen), ViStatus, (ViSession,), vi)
end

function IviFgen_Disable(vi)
    ccall((:IviFgen_Disable, libnifgen), ViStatus, (ViSession,), vi)
end

function IviFgen_InitWithOptions(logicalName, IDQuery, resetDevice, optionString, vi)
    ccall((:IviFgen_InitWithOptions, libnifgen), ViStatus, (ViRsrc, ViBoolean, ViBoolean, ViConstString, Ptr{ViSession}), logicalName, IDQuery, resetDevice, optionString, vi)
end

function IviFgen_GetAttributeViInt32(vi, channelName, attributeId, value)
    ccall((:IviFgen_GetAttributeViInt32, libnifgen), ViStatus, (ViSession, ViConstString, ViAttr, Ptr{ViInt32}), vi, channelName, attributeId, value)
end

function IviFgen_SetAttributeViInt32(vi, channelName, attributeId, value)
    ccall((:IviFgen_SetAttributeViInt32, libnifgen), ViStatus, (ViSession, ViConstString, ViAttr, ViInt32), vi, channelName, attributeId, value)
end

function IviFgen_CheckAttributeViInt32(vi, channelName, attributeId, value)
    ccall((:IviFgen_CheckAttributeViInt32, libnifgen), ViStatus, (ViSession, ViConstString, ViAttr, ViInt32), vi, channelName, attributeId, value)
end

function IviFgen_GetAttributeViInt64(vi, channelName, attributeId, value)
    ccall((:IviFgen_GetAttributeViInt64, libnifgen), ViStatus, (ViSession, ViConstString, ViAttr, Ptr{ViInt64}), vi, channelName, attributeId, value)
end

function IviFgen_SetAttributeViInt64(vi, channelName, attributeId, value)
    ccall((:IviFgen_SetAttributeViInt64, libnifgen), ViStatus, (ViSession, ViConstString, ViAttr, ViInt64), vi, channelName, attributeId, value)
end

function IviFgen_CheckAttributeViInt64(vi, channelName, attributeId, value)
    ccall((:IviFgen_CheckAttributeViInt64, libnifgen), ViStatus, (ViSession, ViConstString, ViAttr, ViInt64), vi, channelName, attributeId, value)
end

function IviFgen_GetAttributeViReal64(vi, channelName, attributeId, value)
    ccall((:IviFgen_GetAttributeViReal64, libnifgen), ViStatus, (ViSession, ViConstString, ViAttr, Ptr{ViReal64}), vi, channelName, attributeId, value)
end

function IviFgen_SetAttributeViReal64(vi, channelName, attributeId, value)
    ccall((:IviFgen_SetAttributeViReal64, libnifgen), ViStatus, (ViSession, ViConstString, ViAttr, ViReal64), vi, channelName, attributeId, value)
end

function IviFgen_CheckAttributeViReal64(vi, channelName, attributeId, value)
    ccall((:IviFgen_CheckAttributeViReal64, libnifgen), ViStatus, (ViSession, ViConstString, ViAttr, ViReal64), vi, channelName, attributeId, value)
end

function IviFgen_GetAttributeViString(vi, channelName, attributeId, bufferSize, value)
    ccall((:IviFgen_GetAttributeViString, libnifgen), ViStatus, (ViSession, ViConstString, ViAttr, ViInt32, Ptr{ViChar}), vi, channelName, attributeId, bufferSize, value)
end

function IviFgen_SetAttributeViString(vi, channelName, attributeId, value)
    ccall((:IviFgen_SetAttributeViString, libnifgen), ViStatus, (ViSession, ViConstString, ViAttr, ViConstString), vi, channelName, attributeId, value)
end

function IviFgen_CheckAttributeViString(vi, channelName, attributeId, value)
    ccall((:IviFgen_CheckAttributeViString, libnifgen), ViStatus, (ViSession, ViConstString, ViAttr, ViConstString), vi, channelName, attributeId, value)
end

function IviFgen_GetAttributeViBoolean(vi, channelName, attributeId, value)
    ccall((:IviFgen_GetAttributeViBoolean, libnifgen), ViStatus, (ViSession, ViConstString, ViAttr, Ptr{ViBoolean}), vi, channelName, attributeId, value)
end

function IviFgen_SetAttributeViBoolean(vi, channelName, attributeId, value)
    ccall((:IviFgen_SetAttributeViBoolean, libnifgen), ViStatus, (ViSession, ViConstString, ViAttr, ViBoolean), vi, channelName, attributeId, value)
end

function IviFgen_CheckAttributeViBoolean(vi, channelName, attributeId, value)
    ccall((:IviFgen_CheckAttributeViBoolean, libnifgen), ViStatus, (ViSession, ViConstString, ViAttr, ViBoolean), vi, channelName, attributeId, value)
end

function IviFgen_GetAttributeViSession(vi, channelName, attributeId, value)
    ccall((:IviFgen_GetAttributeViSession, libnifgen), ViStatus, (ViSession, ViConstString, ViAttr, Ptr{ViSession}), vi, channelName, attributeId, value)
end

function IviFgen_SetAttributeViSession(vi, channelName, attributeId, value)
    ccall((:IviFgen_SetAttributeViSession, libnifgen), ViStatus, (ViSession, ViConstString, ViAttr, ViSession), vi, channelName, attributeId, value)
end

function IviFgen_CheckAttributeViSession(vi, channelName, attributeId, value)
    ccall((:IviFgen_CheckAttributeViSession, libnifgen), ViStatus, (ViSession, ViConstString, ViAttr, ViSession), vi, channelName, attributeId, value)
end

function IviFgen_LockSession(vi, callerHasLock)
    ccall((:IviFgen_LockSession, libnifgen), ViStatus, (ViSession, Ptr{ViBoolean}), vi, callerHasLock)
end

function IviFgen_UnlockSession(vi, callerHasLock)
    ccall((:IviFgen_UnlockSession, libnifgen), ViStatus, (ViSession, Ptr{ViBoolean}), vi, callerHasLock)
end

function IviFgen_GetError(vi, errorCode, bufferSize, description)
    ccall((:IviFgen_GetError, libnifgen), ViStatus, (ViSession, Ptr{ViStatus}, ViInt32, Ptr{ViChar}), vi, errorCode, bufferSize, description)
end

function IviFgen_ClearError(vi)
    ccall((:IviFgen_ClearError, libnifgen), ViStatus, (ViSession,), vi)
end

function IviFgen_GetNextInterchangeWarning(vi, bufferSize, warning)
    ccall((:IviFgen_GetNextInterchangeWarning, libnifgen), ViStatus, (ViSession, ViInt32, Ptr{ViChar}), vi, bufferSize, warning)
end

function IviFgen_ClearInterchangeWarnings(vi)
    ccall((:IviFgen_ClearInterchangeWarnings, libnifgen), ViStatus, (ViSession,), vi)
end

function IviFgen_ResetInterchangeCheck(vi)
    ccall((:IviFgen_ResetInterchangeCheck, libnifgen), ViStatus, (ViSession,), vi)
end

function IviFgen_GetNextCoercionRecord(vi, bufferSize, record)
    ccall((:IviFgen_GetNextCoercionRecord, libnifgen), ViStatus, (ViSession, ViInt32, Ptr{ViChar}), vi, bufferSize, record)
end

function IviFgen_GetSpecificDriverCHandle(vi, specificDriverCHandle)
    ccall((:IviFgen_GetSpecificDriverCHandle, libnifgen), ViStatus, (ViSession, Ptr{ViSession}), vi, specificDriverCHandle)
end

function IviFgen_GetSpecificDriverIUnknownPtr(vi, specificDriverIUnknownPtr)
    ccall((:IviFgen_GetSpecificDriverIUnknownPtr, libnifgen), ViStatus, (ViSession, Ptr{Cvoid}), vi, specificDriverIUnknownPtr)
end

function IviFgen_ConfigureOutputMode(vi, outputMode)
    ccall((:IviFgen_ConfigureOutputMode, libnifgen), ViStatus, (ViSession, ViInt32), vi, outputMode)
end

function IviFgen_ConfigureOperationMode(vi, channelName, operationMode)
    ccall((:IviFgen_ConfigureOperationMode, libnifgen), ViStatus, (ViSession, ViConstString, ViInt32), vi, channelName, operationMode)
end

function IviFgen_ConfigureRefClockSource(vi, refClockSource)
    ccall((:IviFgen_ConfigureRefClockSource, libnifgen), ViStatus, (ViSession, ViInt32), vi, refClockSource)
end

function IviFgen_ConfigureOutputImpedance(vi, channelName, outputImpedance)
    ccall((:IviFgen_ConfigureOutputImpedance, libnifgen), ViStatus, (ViSession, ViConstString, ViReal64), vi, channelName, outputImpedance)
end

function IviFgen_ConfigureOutputEnabled(vi, channelName, enabled)
    ccall((:IviFgen_ConfigureOutputEnabled, libnifgen), ViStatus, (ViSession, ViConstString, ViBoolean), vi, channelName, enabled)
end

function IviFgen_GetChannelName(vi, index, bufferSize, name)
    ccall((:IviFgen_GetChannelName, libnifgen), ViStatus, (ViSession, ViInt32, ViInt32, Ptr{ViChar}), vi, index, bufferSize, name)
end

function IviFgen_InitiateGeneration(vi)
    ccall((:IviFgen_InitiateGeneration, libnifgen), ViStatus, (ViSession,), vi)
end

function IviFgen_AbortGeneration(vi)
    ccall((:IviFgen_AbortGeneration, libnifgen), ViStatus, (ViSession,), vi)
end

function IviFgen_ConfigureStandardWaveform(vi, channelName, waveform, amplitude, dcOffset, frequency, startPhase)
    ccall((:IviFgen_ConfigureStandardWaveform, libnifgen), ViStatus, (ViSession, ViConstString, ViInt32, ViReal64, ViReal64, ViReal64, ViReal64), vi, channelName, waveform, amplitude, dcOffset, frequency, startPhase)
end

function IviFgen_QueryArbWfmCapabilities(vi, maxNumWfms, wfmQuantum, minWfmSize, maxWfmSize)
    ccall((:IviFgen_QueryArbWfmCapabilities, libnifgen), ViStatus, (ViSession, Ptr{ViInt32}, Ptr{ViInt32}, Ptr{ViInt32}, Ptr{ViInt32}), vi, maxNumWfms, wfmQuantum, minWfmSize, maxWfmSize)
end

function IviFgen_CreateArbWaveform(vi, wfmSize, wfmData, wfmHandle)
    ccall((:IviFgen_CreateArbWaveform, libnifgen), ViStatus, (ViSession, ViInt32, Ptr{ViReal64}, Ptr{ViInt32}), vi, wfmSize, wfmData, wfmHandle)
end

function IviFgen_ConfigureSampleRate(vi, sampleRate)
    ccall((:IviFgen_ConfigureSampleRate, libnifgen), ViStatus, (ViSession, ViReal64), vi, sampleRate)
end

function IviFgen_ConfigureArbWaveform(vi, channelName, wfmHandle, arbGain, arbOffset)
    ccall((:IviFgen_ConfigureArbWaveform, libnifgen), ViStatus, (ViSession, ViConstString, ViInt32, ViReal64, ViReal64), vi, channelName, wfmHandle, arbGain, arbOffset)
end

function IviFgen_ClearArbWaveform(vi, wfmHandle)
    ccall((:IviFgen_ClearArbWaveform, libnifgen), ViStatus, (ViSession, ViInt32), vi, wfmHandle)
end

function IviFgen_ConfigureArbFrequency(vi, channelName, frequency)
    ccall((:IviFgen_ConfigureArbFrequency, libnifgen), ViStatus, (ViSession, ViConstString, ViReal64), vi, channelName, frequency)
end

function IviFgen_QueryArbSeqCapabilities(vi, maxNumSeqs, minSeqLength, maxSeqLength, maxLoopCount)
    ccall((:IviFgen_QueryArbSeqCapabilities, libnifgen), ViStatus, (ViSession, Ptr{ViInt32}, Ptr{ViInt32}, Ptr{ViInt32}, Ptr{ViInt32}), vi, maxNumSeqs, minSeqLength, maxSeqLength, maxLoopCount)
end

function IviFgen_CreateArbSequence(vi, seqLength, wfmHandle, wfmLoopCount, seqHandle)
    ccall((:IviFgen_CreateArbSequence, libnifgen), ViStatus, (ViSession, ViInt32, Ptr{ViInt32}, Ptr{ViInt32}, Ptr{ViInt32}), vi, seqLength, wfmHandle, wfmLoopCount, seqHandle)
end

function IviFgen_ConfigureArbSequence(vi, channelName, seqHandle, arbGain, arbOffset)
    ccall((:IviFgen_ConfigureArbSequence, libnifgen), ViStatus, (ViSession, ViConstString, ViInt32, ViReal64, ViReal64), vi, channelName, seqHandle, arbGain, arbOffset)
end

function IviFgen_ClearArbSequence(vi, seqHandle)
    ccall((:IviFgen_ClearArbSequence, libnifgen), ViStatus, (ViSession, ViInt32), vi, seqHandle)
end

function IviFgen_ClearArbMemory(vi)
    ccall((:IviFgen_ClearArbMemory, libnifgen), ViStatus, (ViSession,), vi)
end

function IviFgen_ConfigureTriggerSource(vi, channelName, trigSource)
    ccall((:IviFgen_ConfigureTriggerSource, libnifgen), ViStatus, (ViSession, ViConstString, ViInt32), vi, channelName, trigSource)
end

function IviFgen_ConfigureStartTrigger(vi, channelName, source, slope)
    ccall((:IviFgen_ConfigureStartTrigger, libnifgen), ViStatus, (ViSession, ViConstString, ViConstString, ViInt32), vi, channelName, source, slope)
end

function IviFgen_ConfigureStopTrigger(vi, channelName, source, slope)
    ccall((:IviFgen_ConfigureStopTrigger, libnifgen), ViStatus, (ViSession, ViConstString, ViConstString, ViInt32), vi, channelName, source, slope)
end

function IviFgen_SendSoftwareStopTrigger(vi)
    ccall((:IviFgen_SendSoftwareStopTrigger, libnifgen), ViStatus, (ViSession,), vi)
end

function IviFgen_ConfigureHoldTrigger(vi, channelName, source, slope)
    ccall((:IviFgen_ConfigureHoldTrigger, libnifgen), ViStatus, (ViSession, ViConstString, ViConstString, ViInt32), vi, channelName, source, slope)
end

function IviFgen_SendSoftwareHoldTrigger(vi)
    ccall((:IviFgen_SendSoftwareHoldTrigger, libnifgen), ViStatus, (ViSession,), vi)
end

function IviFgen_ConfigureResumeTrigger(vi, channelName, source, slope)
    ccall((:IviFgen_ConfigureResumeTrigger, libnifgen), ViStatus, (ViSession, ViConstString, ViConstString, ViInt32), vi, channelName, source, slope)
end

function IviFgen_SendSoftwareResumeTrigger(vi)
    ccall((:IviFgen_SendSoftwareResumeTrigger, libnifgen), ViStatus, (ViSession,), vi)
end

function IviFgen_ConfigureAdvanceTrigger(vi, channelName, source, slope)
    ccall((:IviFgen_ConfigureAdvanceTrigger, libnifgen), ViStatus, (ViSession, ViConstString, ViConstString, ViInt32), vi, channelName, source, slope)
end

function IviFgen_SendSoftwareAdvanceTrigger(vi)
    ccall((:IviFgen_SendSoftwareAdvanceTrigger, libnifgen), ViStatus, (ViSession,), vi)
end

function IviFgen_ConfigureInternalTriggerRate(vi, rate)
    ccall((:IviFgen_ConfigureInternalTriggerRate, libnifgen), ViStatus, (ViSession, ViReal64), vi, rate)
end

function IviFgen_SendSoftwareTrigger(vi)
    ccall((:IviFgen_SendSoftwareTrigger, libnifgen), ViStatus, (ViSession,), vi)
end

function IviFgen_ConfigureBurstCount(vi, channelName, count)
    ccall((:IviFgen_ConfigureBurstCount, libnifgen), ViStatus, (ViSession, ViConstString, ViInt32), vi, channelName, count)
end

function IviFgen_ConfigureAMEnabled(vi, channelName, enabled)
    ccall((:IviFgen_ConfigureAMEnabled, libnifgen), ViStatus, (ViSession, ViConstString, ViBoolean), vi, channelName, enabled)
end

function IviFgen_ConfigureAMSource(vi, channelName, source)
    ccall((:IviFgen_ConfigureAMSource, libnifgen), ViStatus, (ViSession, ViConstString, ViInt32), vi, channelName, source)
end

function IviFgen_ConfigureAMInternal(vi, amdepth, amWaveform, amFrequency)
    ccall((:IviFgen_ConfigureAMInternal, libnifgen), ViStatus, (ViSession, ViReal64, ViInt32, ViReal64), vi, amdepth, amWaveform, amFrequency)
end

function IviFgen_ConfigureFMEnabled(vi, channelName, enabled)
    ccall((:IviFgen_ConfigureFMEnabled, libnifgen), ViStatus, (ViSession, ViConstString, ViBoolean), vi, channelName, enabled)
end

function IviFgen_ConfigureFMSource(vi, channelName, source)
    ccall((:IviFgen_ConfigureFMSource, libnifgen), ViStatus, (ViSession, ViConstString, ViInt32), vi, channelName, source)
end

function IviFgen_ConfigureFMInternal(vi, fmdeviation, fmWaveform, fmFrequency)
    ccall((:IviFgen_ConfigureFMInternal, libnifgen), ViStatus, (ViSession, ViReal64, ViInt32, ViReal64), vi, fmdeviation, fmWaveform, fmFrequency)
end

function IviFgen_ConfigureSampleClock(vi, source)
    ccall((:IviFgen_ConfigureSampleClock, libnifgen), ViStatus, (ViSession, ViInt32), vi, source)
end

function IviFgen_ConfigureSampleClockOutputEnabled(vi, enabled)
    ccall((:IviFgen_ConfigureSampleClockOutputEnabled, libnifgen), ViStatus, (ViSession, ViBoolean), vi, enabled)
end

function IviFgen_QueryArbWfmCapabilities64(vi, maxNumWfms, wfmQuantum, minWfmSize, maxWfmSize)
    ccall((:IviFgen_QueryArbWfmCapabilities64, libnifgen), ViStatus, (ViSession, Ptr{ViInt32}, Ptr{ViInt32}, Ptr{ViInt64}, Ptr{ViInt64}), vi, maxNumWfms, wfmQuantum, minWfmSize, maxWfmSize)
end

function IviFgen_CreateChannelArbWaveform(vi, channelName, wfmSize, wfmData, wfmHandle)
    ccall((:IviFgen_CreateChannelArbWaveform, libnifgen), ViStatus, (ViSession, ViConstString, ViInt64, Ptr{ViReal64}, Ptr{ViInt32}), vi, channelName, wfmSize, wfmData, wfmHandle)
end

function IviFgen_CreateChannelArbWaveform16(vi, channelName, wfmSize, wfmData, wfmHandle)
    ccall((:IviFgen_CreateChannelArbWaveform16, libnifgen), ViStatus, (ViSession, ViConstString, ViInt64, Ptr{ViInt16}, Ptr{ViInt32}), vi, channelName, wfmSize, wfmData, wfmHandle)
end

function IviFgen_CreateChannelArbWaveform32(vi, channelName, wfmSize, wfmData, wfmHandle)
    ccall((:IviFgen_CreateChannelArbWaveform32, libnifgen), ViStatus, (ViSession, ViConstString, ViInt64, Ptr{ViInt32}, Ptr{ViInt32}), vi, channelName, wfmSize, wfmData, wfmHandle)
end

function IviFgen_ConfigureDataMarker(vi, name, sourceChannel, bitPosition, destination)
    ccall((:IviFgen_ConfigureDataMarker, libnifgen), ViStatus, (ViSession, ViConstString, ViConstString, ViInt32, ViConstString), vi, name, sourceChannel, bitPosition, destination)
end

function IviFgen_GetDataMarkerName(vi, index, nameBufferSize, name)
    ccall((:IviFgen_GetDataMarkerName, libnifgen), ViStatus, (ViSession, ViInt32, ViInt32, Ptr{ViChar}), vi, index, nameBufferSize, name)
end

function IviFgen_DisableAllDataMarkers(vi)
    ccall((:IviFgen_DisableAllDataMarkers, libnifgen), ViStatus, (ViSession,), vi)
end

function IviFgen_ConfigureSparseMarker(vi, name, wfmHandle, numIndexes, indexes, destination)
    ccall((:IviFgen_ConfigureSparseMarker, libnifgen), ViStatus, (ViSession, ViConstString, ViInt32, ViInt64, Ptr{ViInt64}, ViConstString), vi, name, wfmHandle, numIndexes, indexes, destination)
end

function IviFgen_GetSparseMarkerName(vi, index, nameBufferSize, name)
    ccall((:IviFgen_GetSparseMarkerName, libnifgen), ViStatus, (ViSession, ViInt32, ViInt32, Ptr{ViChar}), vi, index, nameBufferSize, name)
end

function IviFgen_GetSparseMarkerIndexes(vi, name, indexesArraySize, indexes, indexesActualSize)
    ccall((:IviFgen_GetSparseMarkerIndexes, libnifgen), ViStatus, (ViSession, ViConstString, ViInt64, Ptr{ViInt64}, Ptr{ViInt64}), vi, name, indexesArraySize, indexes, indexesActualSize)
end

function IviFgen_SetSparseMarkerIndexes(vi, name, indexesArraySize, indexes)
    ccall((:IviFgen_SetSparseMarkerIndexes, libnifgen), ViStatus, (ViSession, ViConstString, ViInt64, Ptr{ViInt64}), vi, name, indexesArraySize, indexes)
end

function IviFgen_DisableAllSparseMarkers(vi)
    ccall((:IviFgen_DisableAllSparseMarkers, libnifgen), ViStatus, (ViSession,), vi)
end

ext: jl


text:
@inline function DRA(CellData,s,f)
    """
    Discrete Realiastion Algorithm
    # Add License
    # Add Ins and Outs
        # Cell Data
        # Frequency Vector
        # Discretisation Locations
        # Electrode Definition
    """

    DRA_Debug = 0

    # Loop Call Transfer Functions ---------------------------------
    tfft = @. (1/CellData.RA.Fs)*f
    OrgT = @. CellData.RA.SamplingT*(0:floor(tfft[end]/CellData.RA.SamplingT))
    i=Int64(1)

    #Initialise Loop Variables
    puls = Array{Float64}(undef,0,size(OrgT,1)-1)
    D = Array{Float64}(undef,0,1)
    Dtt = Array{String}(undef,0,1)
    C_Aug = Array{Float64}(undef,0,1)
    #tf = Vector{Array}
    # stpsum__ = Array{Float64}(undef,0,length(s))
    #tf__ = Array{Float64}(undef,0,length(s))
    # jk__ = Array{Float64}(undef,length(s),0)
    # testifft__ = Array{Float64}(undef,length(s),0)
    #D_term = Array{Float64}(undef,0,1)
    #res0 = Array{Float64}(undef,0,1)

    for run in CellData.Transfer.tfs[:,1]
        if CellData.Transfer.tfs[i,2] == "Pos"
            tf, Di, res0, Dti = run(CellData,s,CellData.Transfer.tfs[i,3],"Pos") #high compute line
        elseif CellData.Transfer.tfs[i,2] == "Neg"
            tf, Di, res0, Dti = run(CellData,s,CellData.Transfer.tfs[i,3],"Neg") #high compute line
        else
            tf, Di, res0, Dti = run(CellData,s,CellData.Transfer.tfs[i,3]) #high compute line
        end

        #tf = reverse!(tf, dims=2)
        #c = plan_ifft(tf,2)
        #jk = CellData.RA.Fs*real((c*tf)') # inverse fourier transform tranfser function response
        jk = CellData.RA.Fs*real(ifft(tf,2)') # inverse fourier transform tranfser function response (Large Compute)
        stpsum = (cumsum(jk, dims=1).*(1/CellData.RA.Fs))' # cumulative sum of tf response * sample time
        samplingtf = Array{Float64}(undef,size(stpsum,1),length(OrgT))
        # Interpolate H(s) to obtain h_s(s) to obtain discrete-time impulse response
            for Output in 1:size(stpsum,1)
                spl1 = Spline1D(tfft,stpsum[Output,:]; k=3) #Large compute - First to improve
                samplingtf[Output,:]= evaluate(spl1,OrgT)
            end
            puls = [puls; diff(samplingtf, dims=2)]
            D = [D; Di]
            Dtt = [Dtt; Dti]
            C_Aug = [C_Aug; res0]
            i += 1

        if Debug == 1
            println("jk:",jk, "\n")
            println("stpsum:",size(stpsum))
            println("tf:",tf, "\n")
            println("D:",D)
            println("nR:",nR)
        end

        if DRA_Debug == 1
            # stpsum__ = [stpsum__; stpsum]
            # tf__ = [tf__; tf]
            # jk__ = [jk__ jk]
            # testifft__ = [testifft__ testifft]
        end
    end

    if DRA_Debug == 1
        println("tf__:",size(tf__))
        display("text/plain", tf__)

        println("testifft__:",size(testifft__))
        display("text/plain", testifft__)

        println("jk__:",size(jk__))
        display("text/plain", jk__)

        println("stpsum__:")
        display("text/plain", stpsum__)
    end

    #Scale Transfer Functions in Pulse Response
    SFactor = sqrt.(sum(puls.^2,dims=2))
    puls .= puls./SFactor[:,ones(Int64,size(puls,2))]

    #Hankel Formation, perform svd to determine the highest order singular values
    Puls_L = size(puls,1)
    Hank1 = Array{Float64}(undef,length(CellData.RA.H1)*Puls_L,length(CellData.RA.H2))
    Hank2 = Array{Float64}(undef,length(CellData.RA.H1)*Puls_L,length(CellData.RA.H2))

    for lp1 in 1:length(CellData.RA.H2)
        for lp2 in 1:length(CellData.RA.H1)
             Hank1[Puls_L*(lp2-1)+1:Puls_L*lp2,lp1] = @view puls[:,CellData.RA.H2[lp1]+CellData.RA.H1[lp2]+1] #High compute line
             Hank2[Puls_L*(lp2-1)+1:Puls_L*lp2,lp1] = @view puls[:,CellData.RA.H2[lp1]+CellData.RA.H1[lp2]+2] #High compute line
        end
    end

    #Truncated SVD of Hank1 Matrix
    T = svds(Hank1; nsv=CellData.RA.M)[1]
    #U,S,V = tsvd(Hank1, k=CellData.RA.M)

    # Create Observibility and Control Matrices -> Create A, B, and C
    S_ = sqrt(diagm(T.S[1:CellData.RA.M]))
    Observibility = (@view T.U[:,1:CellData.RA.M])*S_
    Control = S_*(@view T.V[:,1:CellData.RA.M])'

    A = Matrix{Float64}(I,CellData.RA.M+1,CellData.RA.M+1)
    A[2:end,2:end] = (Observibility\Hank2)/Control #High compute line (Second)

    #Error check
    if any(i -> i>1., real(eigvals(A)))
        println("Oscilating System: A has indices of values > 1")
    end

    if any(i -> i<0., real(eigvals(A)))
        println("Unstable System: A has indices of negative values")
    end

    B = [CellData.RA.SamplingT; Control[:,1:CellData.RA.N]]
    C = [C_Aug SFactor[:,ones(Int64,CellData.RA.M)].*Observibility[1:size(puls,1),:]]

    #  Final State-Space Form
    d, S = eigen(A)
    A = Diagonal(d)
    B = S\B
    C = C*S
    # sys = diagonalize(ss(A,B,C,D,CellData.RA.SamplingT))
    # ss_B = sys.B
    # ss_C = sys.C

    #Scale C and tansform B to improve real-time linearisation performance
    # ss_C = ss_C.*ss_B'
    # ss_B = ones(size(ss_B))

    C = C.*B'
    B = ones(size(B))

return real(A), real(B), real(C), D
end

max_stars_repo_name: BradyPlanden/LIBRA

max_issues_repo_name: BradyPlanden/LIBRA

max_forks_repo_name: BradyPlanden/LIBRA


text:
module ShearFunctions

using StaticArrays
using ..PhysicalVectors

export b1g_shear, b2g_shear, bxg_shears!

"""
    b1g_shear( eval_r::Vector2D, bob::Vector2D )

Calculate the ``B_{1g}`` from an edge dislocation for a given Burgers vector `bob` at site `eval_r`.

# Additional Information
* The edge dislocation is aligned along the ``z`` axis at the origin.
* This is in units of ``1/(1-σ)`` for a Poisson ratio ``σ``.

# Examples
```jldocstest
julia> b1g_shear( Vector2D(1.,1.), Vector2D(1.,0.) )
-0.07957747154594767
```
"""
function b1g_shear( eval_r::Vector2D, bob::Vector2D )
    return -1/π * ( eval_r.vec[1] * eval_r.vec[2] / magnitude2(eval_r) ) * (bob ⋅ eval_r) / magnitude2(eval_r)
end

"""
    b2g_shear( eval_r::Vector2D, bob::Vector2D )

Calculate the ``B_{2g}`` from an edge dislocation for a given Burgers vector `bob` at site `eval_r`.

# Additional Information
* The edge dislocation is aligned along the ``z`` axis at the origin.
* This is in units of ``1/(1-σ)`` for a Poisson ratio ``σ``.

# Examples
```jldocstest
julia> b2g_shear( Vector2D(1.,1.), Vector2D(1.,0.) )
0.15915494309189535
```
"""
function b2g_shear( bob::Vector2D, eval_r::Vector2D )
    return 1/(2*π) * ( ( eval_r.vec[1]^2 - eval_r.vec[2]^2 ) / magnitude2(eval_r) ) * (bob ⋅ eval_r) / magnitude2(eval_r)
end

"""
    bxg_shears!(eval_r::Vector2D, bob::Vector2D; diff::Function, source_r::Vector2D = Vector2D(0.,0.)) -> Tuple{T, T} where {T <: Real}

Calculate the two shears `(b1g, b2g)` from the ``B_{1g}`` and ``B_{2g}`` channels from and edge dislocation situated at the `source_r`
location and aligned along the ``z`` axis with a Burger's vector `bob`.

# Additional Information
* This function mutates both and `eval_r`. It should be used in cases where `eval_r` is a temporary vector.
* The keyword argument `diff <: Function` should be used in cases where `Vector2D` subtraction has a different definition than the normal one expected, for example with periodic boundary conditions.

# Examples
```jldocstest
julia> using StaticArrays

julia> A = MVector{2}(0., 0.)
2-element MVector{2, Float64} with indices SOneTo(2):
 0.0
 0.0

julia> bxg_shears!( A, Vector2D(1.,1.), Vector2D(1.,0.) )

julia> A
2-element MVector{2, Float64} with indices SOneTo(2):
 -0.07957747154594767
  0.15915494309189535
```
"""
function bxg_shears!( eval_r::Vector2D, bob::Vector2D; diff::Function = subtract!, source_r::Vector2D = Vector2D(0.,0.) )
    diff(eval_r, source_r)
    b1g = b1g_shear(eval_r, bob)
    b2g = b2g_shear(eval_r, bob)
    return (b1g, b2g)
end

end # module ShearFunctions


text:
type VectorF{F,SizeType} <: Function
  f::F
  sizeu::SizeType
end
function (p::VectorF)(t,uprev,du)
  p.f(t,reshape(uprev,p.sizeu...),reshape(du,p.sizeu...))
  uprev = vec(uprev)
  du = vec(du)
end

type VectorFReturn{F,SizeType} <: Function
  f::F
  sizeu::SizeType
end
function (p::VectorFReturn)(t,uprev,du)
  p.f(t,reshape(uprev,p.sizeu...),reshape(du,p.sizeu...))
  vec(du)
end

type TimeGradientWrapper{VFType,uType,rateType} <: Function
  vf::VFType
  uprev::uType
  du2::rateType
end
(p::TimeGradientWrapper)(t) = p.vf(t,p.uprev,p.du2)

type UJacobianWrapper{VFRType,tType} <: Function
  vfr::VFRType
  t::tType
end

(p::UJacobianWrapper)(du1,uprev) = p.vfr(p.t,uprev,du1)

type TimeDerivativeWrapper{F,uType} <: Function
  f::F
  u::uType
end
(p::TimeDerivativeWrapper)(t) = p.f(t,p.u)

type UDerivativeWrapper{F,tType} <: Function
  f::F
  t::tType
end
(p::UDerivativeWrapper)(u) = p.f(p.t,u)


text:
using ModelStatisticsLH, Random, Test

ms = ModelStatisticsLH;

reduceFct(x :: Vector{T}) where T <: AbstractFloat =
    T(sum(x) / length(x));
reduceFct(x :: Vector{T}) where T <: Integer =
    round(T, sum(x) / length(x));

function reduce_one_test()
    rng = MersenneTwister(123);
    @testset "Reduce one field" begin
        n = 5;
        statV = [ms.make_test_model_stats(:x, 4, rng)  for j = 1 : n];
        varNameV = var_names(statV[1]);
        for varName in varNameV
            xMean = reduce_one_field(statV, varName, reduceFct);
            if !isnothing(xMean)
                @test typeof(xMean) == typeof(get_values(statV[1], varName));
            end
        end
	end
end

function reduce_stats_test()
    rng = MersenneTwister(43);
    @testset "reduce stats" begin
        n = 3;
        scV = [ms.make_test_scollection(:x, rng)  for j = 1 : n];
        sc = reduce_stats(scV, reduceFct);
        @test isequal(get_groups(sc), get_groups(scV[1]));

        for grp in get_groups(sc)
            mStats = get_mstats(sc, grp);
            varNameV = var_names(mStats);
            mStatsV = [get_mstats(scV[j], grp)  for j = 1 : n];
            for varName in varNameV
                # valueV = [get_stats(mStatsV[j], varName)  for j = 1 : n];
                xr = reduce_one_field(mStatsV, varName, reduceFct);
                xr2 = get_stats(sc, grp, varName);
                if isnothing(xr)
                    @test isnothing(xr2)
                elseif xr isa Array
                    for j in eachindex(xr)
                        @test isapprox(xr[j], xr2[j])
                    end
                else
                    @test isapprox(xr, xr2)
                end
            end
        end
	end
end

@testset "Reductions" begin
    reduce_one_test();
    reduce_stats_test();
end

# -----------


text:
using ODEFilters
using Test
using LinearAlgebra

using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_lotkavoltera

h = rand()
σ = rand()

@testset "Test vanilla (ie. non-preconditioned) IBM (d=2,q=2)" begin
    # Tests are broken, because there is no option right now to get A and Q without
    # the preconditioning
    d, q = 2, 2

    A!, Q! = ODEFilters.vanilla_ibm(d, q)
    Ah = diagm(0 => ones(d*(q+1)))
    Qh = zeros(d*(q+1), d*(q+1))
    A!(Ah, h)
    Q!(Qh, h, σ^2)

    AH_22_IBM = [1 0 h 0 h^2/2 0;
                 0 1 0 h 0 h^2/2;
                 0 0 1 0 h 0;
                 0 0 0 1 0 h;
                 0 0 0 0 1 0;
                 0 0 0 0 0 1]
    @test AH_22_IBM ≈ Ah

    QH_22_IBM = σ^2 .* [h^5/20  0       h^4/8  0      h^3/6  0;
                        0       h^5/20  0      h^4/8  0      h^3/6;
                        h^4/8   0       h^3/3  0      h^2/2  0;
                        0       h^4/8   0      h^3/3  0      h^2/2;
                        h^3/6   0       h^2/2  0      h      0;
                        0       h^3/6   0      h^2/2  0      h]
    @test QH_22_IBM ≈ Qh
end

@testset "Test IBM with preconditioning (d=1,q=2)" begin
    d, q = 1, 2

    A, Q = ODEFilters.ibm(d, q)
    Qh = Q * σ^2

    AH_21_PRE = [1  1  0.5
                 0  1  1
                 0  0  1]

    QH_21_PRE = σ^2 * [1/20 1/8 1/6
                       1/8  1/3 1/2
                       1/6  1/2 1]

    @test AH_21_PRE ≈ A
    @test QH_21_PRE ≈ Qh
end

@testset "Verify correct prior dim" begin
    prob = prob_ode_lotkavoltera
    d = length(prob.u0)
    for q in 1:5
        integ = init(prob, EK0(order=q), initialize_derivatives=false)
        @test length(integ.cache.x.μ) == d*(q+1)
        sol = solve!(integ)
        @test length(integ.cache.x.μ) == d*(q+1)
        @test length(sol.x[end].μ) == d*(q+1)
    end
end


text:
@testset "UniformMatroid" begin
    @testset "Uniform matroid" begin
        r = Float64[5, 4, 3]

        @testset "Interface" begin
            @testset "Global interface" begin
                i = UniformMatroidInstance(r, 1)
                @test objective(i) == Maximise()
            end

            @testset "Invalid parameter" begin
                @test_throws ErrorException UniformMatroidInstance(r, -1)
                @test_throws ErrorException UniformMatroidInstance(r, 0)
            end

            @testset "Copy" begin
                i = UniformMatroidInstance(Float64[5, 4, 3], 1)
                i2 = copy(i)
                @test i.m == i2.m
                @test i.rewards == i2.rewards
                @test i.objective == i2.objective
            end

            @testset "Approximation" begin
                i = UniformMatroidInstance(r, 1)

                @test approximation_ratio(i, GreedyAlgorithm()) == 1.0
                @test approximation_term(i, GreedyAlgorithm()) == 0.0

                @test approximation_ratio(i, DynamicProgramming()) == 1.0
                @test approximation_term(i, DynamicProgramming()) == 0.0

                @test approximation_ratio(i, DefaultLinearFormulation(Cbc.Optimizer)) == 1.0
                @test approximation_term(i, DefaultLinearFormulation(Cbc.Optimizer)) == 0.0
            end

            @testset "Make solution" begin
                i = UniformMatroidInstance(r, 1)
                d = Dict(1 => 0.4, 2 => 0.8, 3 => 0.7)
                s = make_solution(i, d)

                @test s.instance == i
                @test length(s.variables) == 2
                @test 1 ∉ s.variables
                @test 2 ∈ s.variables
                @test 3 ∈ s.variables
            end

            @testset "Solution helpers" begin
                i = UniformMatroidInstance(r, 1)

                # Infeasible solution.
                s = UniformMatroidSolution(i, Int[])
                @test value(s) == -Inf

                # Valid solution.
                d = Dict(1 => 0.4, 2 => 0.8, 3 => 0.7)
                s = make_solution(i, d)

                @test value(s) == 7
            end
        end

        @testset "Basic" begin
            m = 2
            i = UniformMatroidInstance(r, m)
            g = solve(i, GreedyAlgorithm())
            d = solve(i, DynamicProgramming())
            l = solve(i, DefaultLinearFormulation(Cbc.Optimizer))

            @test i.m == m

            @test value(g) == value(d)
            @test length(g.variables) <= m
            @test length(d.variables) <= m

            @test g.instance == i
            @test d.instance == i

            if !is_travis
                @test value(g) == value(l)
                @test length(l.variables) <= m
                @test l.instance == i
            end
        end
    end

    @testset "Budgeted uniform matroid" begin
        @testset "Interface" begin
            r = Float64[5, 4, 3]
            w = Int[1, 1, 1]

            @testset "Global interface" begin
                i = MinimumBudget(UniformMatroidInstance(r, 1), w, 1)
                @test objective(i) == Maximise()
            end

            @testset "Errors when building the UniformMatroidInstance should bubble up to MinimumBudget" begin
                @test_throws ErrorException MinimumBudget(
                    UniformMatroidInstance(r, -1),
                    w,
                )
                @test_throws ErrorException MinimumBudget(
                    UniformMatroidInstance(r, 0),
                    w,
                )
            end

            @testset "Negative weight or budget" begin
                @test_throws ErrorException MinimumBudget(UniformMatroidInstance(r, 2), w, -1)
                @test_throws ErrorException MinimumBudget(UniformMatroidInstance(r, 2), Int[1, -1, 1])
            end

            @testset "Different number of items between the matroid and the new constraint" begin
                @test_throws ErrorException MinimumBudget(
                    UniformMatroidInstance(r, 0),
                    w[1:2],
                )
                @test_throws ErrorException MinimumBudget(
                    UniformMatroidInstance(r[1:2], 0),
                    w,
                )
            end

            @testset "Copy" begin
                i = MinimumBudget(UniformMatroidInstance(r, 1), w, 1)
                i2 = copy(i)

                @test i.weights == i2.weights
                @test i.min_budget == i2.min_budget
                @test i.compute_all_values == i2.compute_all_values
                @test i.instance.m == i2.instance.m
                @test i.instance.rewards == i2.instance.rewards
                @test i.instance.objective == i2.instance.objective
            end

            @testset "Approximation" begin
                i = MinimumBudget(UniformMatroidInstance(r, 1), w, 1)

                @test approximation_ratio(i, DynamicProgramming()) == 1.0
                @test approximation_term(i, DynamicProgramming()) == 0.0

                @test approximation_ratio(i, DefaultLinearFormulation(Cbc.Optimizer)) == 1.0
                @test approximation_term(i, DefaultLinearFormulation(Cbc.Optimizer)) == 0.0
            end

            @testset "Make solution" begin
                i = MinimumBudget(UniformMatroidInstance(r, 1), w, 1)
                d = Dict(1 => 0.4, 2 => 0.8, 3 => 0.7)
                s = make_solution(i, d)

                @test s.instance == i
                @test length(s.variables) == 2
                @test 1 ∉ s.variables
                @test 2 ∈ s.variables
                @test 3 ∈ s.variables
            end

            @testset "Solution helpers" begin
                i = MinimumBudget(UniformMatroidInstance(r, 1), w, 1)

                # A first solution, with no per-budget information.
                d = Dict(1 => 0.4, 2 => 0.8, 3 => 0.7)
                s = make_solution(i, d)
                @test value(s) == 7
                @test_throws ErrorException items(s, 0)
                @test_throws ErrorException items_all_budgets(s, 0)
                @test_throws ErrorException value(s, 0)

                # An infeasible solution.
                s = MinBudgetedUniformMatroidSolution(
                    i,
                    Int[],
                    Dict{Int, Vector{Int}}()
                )
                @test value(s) == -Inf
                @test_throws ErrorException items(s, 0)
                @test_throws ErrorException items_all_budgets(s, 0)
                @test_throws ErrorException value(s, 0)

                # A more complete solution.
                d = Dict(0 => [1, 2], 1 => [1, 3], 2 => [2, 3], 3 => Int[])
                s = MinBudgetedUniformMatroidSolution(
                    i,
                    [2, 3],
                    d
                )
                @test value(s) == 7 # For the solution [2, 3].
                @test items(s, 0) == [1, 2]
                @test value(s, 0) == 9
                @test items(s, 1) == [1, 3]
                @test value(s, 1) == 8
                @test items(s, 2) == [2, 3]
                @test value(s, 2) == 7
                @test items(s, 3) == Int[]
                @test value(s, 3) == -Inf

                @test items_all_budgets(s, 0) == Dict(0 => [1, 2])
                @test items_all_budgets(s, 1) == Dict(0 => [1, 2], 1 => [1, 3])
                @test items_all_budgets(s, 2) == Dict(0 => [1, 2], 1 => [1, 3], 2 => [2, 3])
                @test items_all_budgets(s, 3) == Dict(0 => [1, 2], 1 => [1, 3], 2 => [2, 3], 3 => Int[])
            end
        end

        function test_solution_at(
            s::MinBudgetedUniformMatroidSolution{Float64, Int},
            kv::Dict{Int, Float64},
        )
            for (k, v) in kv
                @test value(s, k) ≈ v
            end
        end

        function test_items_at(
            s::MinBudgetedUniformMatroidSolution{Float64, Int},
            kv::Dict{Int, Vector{Int}},
        )
            for (k, v) in kv
                @test items(s, k) == v
            end
        end

        @testset "Basic (all budget values)" begin
            m = 2
            i = MinimumBudget(
                UniformMatroidInstance(Float64[5, 4, 3], m),
                Int[1, 1, 1],
                3,
                compute_all_values=true,
            )
            d = solve(i, DynamicProgramming())
            l = solve(i, DefaultLinearFormulation(Cbc.Optimizer))

            expected_items = Dict{Int, Vector{Int}}(
                0 => [1, 2],
                1 => [1, 2],
                2 => [1, 2],
                3 => [-1],
            )
            expected =
                Dict{Int, Float64}(0 => 9.0, 1 => 9.0, 2 => 9.0, 3 => -Inf)

            @test i.instance.m == m

            @testset "Dynamic programming" begin
                @test d.instance == i

                @test d.solutions[0] == [1, 2]
                @test d.solutions[1] == [1, 2]
                @test d.solutions[2] == [1, 2]
                @test d.solutions[3] == [-1]

                test_items_at(d, expected_items)
                test_solution_at(d, expected)
            end

            @testset "Linear programming" begin
                @test l.instance == i

                for i in 0:3
                    @test l.solutions[i] == d.solutions[i]
                end

                test_items_at(l, expected_items)
                test_solution_at(l, expected)
            end
        end

        @testset "Basic (one budget value)" begin
            m = 2
            i = MinimumBudget(
                UniformMatroidInstance(Float64[5, 4, 3], m),
                Int[1, 1, 1],
                3,
                compute_all_values=false,
            )
            d = solve(i, DynamicProgramming())
            l = solve(i, DefaultLinearFormulation(Cbc.Optimizer))

            expected_items = Dict(3 => [-1])
            expected = Dict{Int, Float64}(3 => -Inf)

            @test i.instance.m == m

            @testset "Dynamic programming" begin
                @test d.instance == i
                @test d.solutions[3] == [-1]

                test_items_at(d, expected_items)
                test_solution_at(d, expected)
            end

            @testset "Linear programming" begin
                @test l.instance == i
                @test d.solutions[3] == [-1]

                test_items_at(l, expected_items)
                test_solution_at(l, expected)
            end
        end

        # More advanced tests to ensure the algorithm works as expected.

        @testset "Conformity 1" begin
            a = 3.612916190062782
            b = 7.225832380125564
            v = Float64[a, a, b, b, b, b, b, b, b, b]
            w = Int[32, 32, 32, 32, 0, 32, 32, 32, 0, 32]
            m = 3

            i = MinimumBudget(
                UniformMatroidInstance(v, m),
                w,
                320,
                compute_all_values=true,
            )
            d = solve(i, DynamicProgramming())
            l = solve(i, DefaultLinearFormulation(Cbc.Optimizer))

            expected = Dict{Int, Float64}(
                0 => 3 * b,
                96 => 3 * b,
                96 + 1 => -Inf,
                320 => -Inf,
            ) # No more solutions after 96.

            @testset "Dynamic programming" begin
                test_solution_at(d, expected)
            end

            @testset "Linear programming" begin
                test_solution_at(l, expected)
            end
        end

        @testset "Conformity 2" begin
            v = [
                7.840854066284411,
                3.9204270331422055,
                7.840854066284411,
                3.9204270331422055,
                7.840854066284411,
                7.840854066284411,
                7.840854066284411,
                7.840854066284411,
                15.681708132568822,
                5.227236044189607,
            ]
            w = [16, 32, 16, 24, 16, 16, 0, 16, 0, 32]
            m = 3

            i = MinimumBudget(
                UniformMatroidInstance(v, m),
                w,
                320,
                compute_all_values=true,
            )
            d = solve(i, DynamicProgramming())
            l = solve(i, DefaultLinearFormulation(Cbc.Optimizer))

            a = 31.363416265137644
            b = 28.74979824304284
            c = 24.829371209900632
            e = 16.988517143616225
            expected = Dict{Int, Float64}(
                0 => a,
                4 => a,
                20 => a,
                24 => a,
                25 => a,
                32 => a,
                33 => b,
                40 => b,
                41 => b,
                48 => b,
                49 => c,
                72 => e,
                73 => e,
                96 => -Inf,
                97 => -Inf,
                280 => -Inf,
                319 => -Inf,
                320 => -Inf,
            )

            @testset "Dynamic programming" begin
                test_solution_at(d, expected)
            end

            @testset "Linear programming" begin
                test_solution_at(l, expected)
            end
        end
    end
end


text:
using Bridge, Distributions, StaticArrays
using Test
#import Bridge: b, σ, a, transitionprob
n = 1000
TT = 5
tt = 0.0:TT/n:TT
m = 1000
S = SVector{2,Float64}
M = SArray{Tuple{2,2},Float64,2,4}
B = SMatrix{2,2,Float64}([-1 0.1; -0.2 -1])
mu = 0.1*S([0.2, 0.3])
sigma = SMatrix{2,2,Float64}(2*[-0.212887  0.0687025
  0.193157  0.388997 ])
a = sigma*sigma'

P = LinPro(B, mu, sigma)

u = S([1.0, 0.0])
v = S([.5, 0.0])

@test (norm(Matrix(-P.lambda*B' - B*P.lambda - a))) < eps()

t = 0.5
T = 2.0

@test Bridge.bderiv(t, v, P) == B
@test Bridge.σderiv(t, v, P) == zero(sigma)

dt = 1e-6
@test norm( (Bridge.V(t+dt, T, v, P)-Bridge.V(t, T, v, P))/dt - Bridge.dotV(t, T, v, P) ) < 10*dt

n2 = 150
tt = range(t, stop=T, length=n2)
K = SamplePath(tt, zeros(M, length(tt)))
Bridge.gpHinv!(K, P)
V = SamplePath(tt, zeros(S, length(tt)))
Bridge.gpHinv!(K, P) # warm up
Bridge.gpV!(V, P, v)
Mu = SamplePath(tt, zeros(S, length(tt)))
Mu2 = SamplePath(tt, zeros(S, length(tt)))

solve!(Bridge.R3(), Bridge.b, Mu, u, P)
solve!(BS3(), Bridge.b, Mu2, u, P)
@test (@allocated Bridge.gpHinv!(K, P)) <= 32
@test (@allocated Bridge.gpV!(V, P, v)) == 0
@test norm(K.yy[1]*Bridge.H(t, T, P) - I) < 10/n2^3
@test norm(V.yy[1] - Bridge.V(t, T, v, P)) < 10/n2^3

@test norm(Mu.yy[end] - Bridge.mu(t, u, T, P)) < 10/n2^3
@test norm(Mu2.yy[end] - Bridge.mu(t, u, T, P)) < 10/n2^3

@test norm(Bridge.K(t, T, P) - solve(Bridge.R3(), Bridge._dK, tt, zero(M), P)) < 10/n2^3

# Normal(mu, lambda) is the stationary distribution. check by starting in stationary distribution and evolve 20 time units
X = Bridge.mat(S[solve(EulerMaruyama(), mu + cholupper(P.lambda)*randn(S), sample(tt, Wiener{S}()),P).yy[end] - mu for i in 1:m])

@test supnorm(cov(X, dims=2) - Matrix(P.lambda)) < .1

n = 1000
TT = 0.5
tt = 0.0:TT/n:TT

X = solve(EulerMaruyama(), mu, sample(tt, Wiener{S}()),P)
@test supnorm(quvar(X)-TT*a) < 0.06

theta = 0.7
@test_skip  begin X = Bridge.mat(S[thetamethod(mu + chol(P.lambda)*randn(S), sample(tt, Wiener{S}()), P, theta).yy[end] - mu for i in 1:m])
  supnorm(cov(X,2) - Matrix(P.lambda)) < 1.0
end

Pt = Bridge.Ptilde(Bridge.CSpline(tt[1], tt[end], 1.0, 0.0, 0.0, 1.0), sigma)

@test Bridge.gamma(t, v, Pt) ≈ inv(sigma*sigma')


text:
# BLAS operations
import Base: dot, A_mul_Bt, At_mul_B, At_mul_Bt, A_mul_Bc,
    Ac_mul_B, Ac_mul_Bc, transpose, ctranspose, transpose!, ctranspose!

dot{T,S}(lhs::AFAbstractArray{T}, rhs::AFAbstractArray{S}) =
        AFArray{af_promote(T,S)}(af_dot(lhs, rhs))

# Matmul
*{T,S}(a::AFAbstractArray{T}, b::AFAbstractArray{S}) =
    AFArray{af_promote(T,S)}(af_matmul(a,b))
*{T,S,V}(a::AFAbstractArray{T}, b::AFAbstractArray{S}, c::AFAbstractArray{V}) =
    AFArray{af_promote(af_promote(T,S), V)}(af_matmul3(a,b,c))
*{T,S,V,W}(a::AFAbstractArray{T}, b::AFAbstractArray{S}, c::AFAbstractArray{V}, d::AFAbstractArray{W}) =
    AFArray{af_promte(af_promote(af_promote(T,S), V), W)}(af_matmul4(a,b,c,d))

function _matmul(a::AFAbstractArray, b::AFAbstractArray;
    lhsProp = AF_MAT_NONE, rhsProp = AF_MAT_NONE)
    out = af_matmul_flags(a, b, lhsProp, rhsProp)
    AFArray{backend_eltype(out)}(out)
end

# with transpose
function A_mul_Bt(a::AFAbstractArray, b::AFAbstractArray)
    out = af_matmulNT(a,b)
    AFArray{backend_eltype(out)}(out)
end

function At_mul_B(a::AFAbstractArray, b::AFAbstractArray)
    out = af_matmulTN(a,b)
    AFArray{backend_eltype(out)}(out)
end

function At_mul_Bt(a::AFAbstractArray, b::AFAbstractArray)
    out = af_matmulTT(a,b)
    AFArray{backend_eltype(out)}(out)
end

# with complex conjugate
A_mul_Bc(a::AFAbstractArray, b::AFAbstractArray) =
    _matmul(a,b,rhsProp=AF_MAT_CTRANS)
Ac_mul_B(a::AFAbstractArray, b::AFAbstractArray) =
    _matmul(a,b,lhsProp=AF_MAT_CTRANS)
Ac_mul_Bc(a::AFAbstractArray, b::AFAbstractArray) =
    _matmul(a,b,lhsProp=AF_MAT_CTRANS,rhsProp=AF_MAT_CTRANS)

# transpose
transpose{T}(x::AFAbstractArray{T}) = AFArray{T}(af_transpose(x))
ctranspose{T}(x::AFAbstractArray{T}) = AFArray{T}(af_ctranspose(x))

transpose!{T}(x::AFAbstractArray{T}) = af_transposeInPlace(x)
ctranspose!{T}(x::AFAbstractArray{T}) = af_ctransposeInPlace(x)

# solve

# TODO : The documentation says only AF_MAT_LOWER/AF_MAT_UPPER are supported
# once AF_MAT_(C)TRANS is supported this could be useful for A_rdiv, etc
# TODO : Think about integrating solveLU in `\` so it becomes a poly algorithm like base.

import Base: \

\{S,T}(a::AFAbstractArray{S}, b::AFAbstractArray{T}) = AFArray{af_promote(T,S)}(af_solve(a, b);)

# Factorizations

import Base.LinAlg: chol, chol!, PosDefException

#Cholesky
function _chol{T}(a::AFMatrix{T}, is_upper::Bool)
    out = AFArray()
    info = af_cholesky(out, a, is_upper)
    info > 0 && throw(PosDefException(info))
    out = is_upper ? (AFArray{T}(af_upper(out))) : (AFArray{T}(af_lower(out)))
end

function _chol!{T}(a::Matrix{T}, is_upper::Bool)
    info = af_choleskyInPlace(a, is_upper)
    info > 0 && throw(PosDefException(info))
    b = is_upper ? (AFArray{T}(af_upper(a))) : (AFArray{T}(af_lower(a)))
    return b
end

chol(a::AFMatrix, ::Type{Val{:U}}) = _chol(a, true)
chol(a::AFMatrix, ::Type{Val{:L}}) = _chol(a, false)
chol(a::AFMatrix) = chol(a,Val{:U})

chol!(a::AFMatrix, ::Type{Val{:U}}) = _chol!(a, true)
chol!(a::AFMatrix, ::Type{Val{:L}}) = _chol!(a, false)
chol!(a::AFMatrix) = chol!(a,Val{:U})

#LU
function lu(a::AFMatrix)
    l = AFArray()
    u = AFArray()
    p = AFArray()
    af_lu(l, u, p, a)
    AFArray{backend_eltype(l)}(l), AFArray{backend_eltype(u)}(u), (AFArray{backend_eltype(p)}(p) + 1)
end

#QR
function qr(a::AFMatrix)
    q = AFArray()
    r = AFArray()
    tau = AFArray()
    af_qr(q, r, tau, a)
    AFArray{backend_eltype(q)}(q), AFArray{backend_eltype(r)}(r), AFArray{backend_eltype(tau)}(tau)
end

#SVD
function svd(a::AFMatrix)
    u = AFArray()
    s = AFArray()
    vt = AFArray()
    af_svd(u, s, vt, a)
    AFArray{backend_eltype(u)}(u), AFArray{backend_eltype(s)}(s), AFArray{backend_eltype(vt)}(vt)
end

#Mat ops

import Base: det, inv, norm

function det{T}(a::AFMatrix{T})
    if ndims(a) != 2
        throw(DimensionMismatch("Input isn't a matrix"))
    else
        return af_det(a)
    end
end

function inv(a::AFMatrix)
    if ndims(a) != 2
        throw(DimensionMismatch("Input isn't a matrix"))
    else
        return AFArray{Float32}(af_inverse(a))
    end
end

norm(a::AFAbstractArray) = af_norm(a)
rank(a::AFMatrix) = af_rank(a)


text:
# Plasmo Interface to Dsp
module PlasmoDspInterface

include("DspCInterface.jl")
using MPI
using SparseArrays
import .DspCInterface
import Plasmo.PlasmoModelGraph
import JuMP

export dsp_solve

#Build up Dsp Model using structure from Plasmo
# This function is hooked by JuMP (see block.jl)
function dsp_solve(graph::PlasmoModelGraph.ModelGraph,master_node::PlasmoModelGraph.ModelNode,children_nodes::Vector{PlasmoModelGraph.ModelNode};
                        probabilities = Dict(zip(1:length(children_nodes),fill(1/length(children_nodes),length(children_nodes)))),suppress_warnings = false, options...)
    master = PlasmoModelGraph.getmodel(master_node)
    submodels = [PlasmoModelGraph.getmodel(child) for child in children_nodes]
    linkconstraints = PlasmoModelGraph.getlinkconstraints(graph) #get all of the link constraints in the graph

    dspmodel = DspCInterface.DspModel()
    DspCInterface.freeModel(dspmodel)
    # parse options
    for (optname, optval) in options
        if optname == :param
            DspCInterface.readParamFile(dspmodel, optval)
        elseif optname == :solve_type
            if optval in [:Dual, :Benders, :Extensive]
                dspmodel.solve_type = optval
            else
                @warn("solve_type $optval is not available.")
            end
        else
            @warn("Options $optname is not available.")
        end
    end

    nblocks = length(submodels)
    DspCInterface.setBlockIds(dspmodel, nblocks)
    #setBlockIds(dspmodel, nblocks)
    #println(dspmodel)
    # load the plasmo models to Dsp
    loadProblem(dspmodel,graph, master, submodels,probabilities)
    # println(DspCInterface.getTotalNumCols(dspmodel))
    # println(DspCInterface.getTotalNumRows(dspmodel))
    # println(DspCInterface.getNumCouplingRows(dspmodel))
    # solve
    DspCInterface.solve(dspmodel)

    # solution status
    statcode = DspCInterface.getSolutionStatus(dspmodel)
    stat = parseStatusCode(statcode)

    # Extract solution from the solver
    dspmodel.numRows = DspCInterface.getNumRows(dspmodel, 0) + DspCInterface.getNumRows(dspmodel, 1) * DspCInterface.getNumScenarios(dspmodel)
    dspmodel.numCols = DspCInterface.getTotalNumCols(dspmodel)

    #this should be the graph objective and column values?
    #TODO
    master.objVal = NaN
    master.colVal = fill(NaN, dspmodel.numCols)

    if stat != :Optimal
        suppress_warnings || @warn("Not solved to optimality, status: $stat")
    end

    if !(stat == :Infeasible || stat == :Unbounded)
        getDspSolution(dspmodel,master,submodels)  #add solution data to master and children models
    end
    PlasmoModelGraph._setobjectivevalue(graph,dspmodel.primVal)
    # Return the solve status
    return stat
    #return dspmodel
end

function loadProblem(dsp::DspCInterface.DspModel,graph::PlasmoModelGraph.ModelGraph, master::JuMP.Model, subproblems::Vector{JuMP.Model}, dedicatedMaster::Bool,probabilities)
    DspCInterface.check_problem(dsp)

    #if haskey(model.ext, :DspBlocks) #if there are children models given.....
    if length(subproblems) > 0
        loadStochasticProblem(dsp, graph, master,subproblems, dedicatedMaster,probabilities)
        #loadStochasticProblem(dsp, graph, master,subproblems, dedicatedMaster,probabilities = Dict(zip(1:length(subproblems),fill(1.0,length(subproblems)))))

    else #Do I even need to support this?
        @warn("No blocks were defined.")
        loadDeterministicProblem(dsp, master)
    end
end
loadProblem(dsp::DspCInterface.DspModel,graph::PlasmoModelGraph.ModelGraph,master::JuMP.Model, subproblems::Vector{JuMP.Model},probabilites::Dict) = loadProblem(dsp,graph, master, subproblems, true,probabilites);

function loadStochasticProblem(dsp::DspCInterface.DspModel, graph::PlasmoModelGraph.ModelGraph, master::JuMP.Model, subproblems::Vector{JuMP.Model}, dedicatedMaster::Bool, probabilities::Dict)
    # model was a Dsp JuMP model
    # get DspBlocks
    #blocks = model.ext[:DspBlocks]    #this is a blockstructure with ids and weights
    #blocks = subproblems #maybe this will work?

    nscen  = dsp.nblocks
    ncols1 = master.numCols
    nrows1 = length(master.linconstr)
    ncols2 = 0
    nrows2 = 0
    #for s in values(blocks.children)  #these are models
    #Do I need to figure out which variables are actually first stage?
    for s in subproblems
        node = PlasmoModelGraph.getnode(s)
        link_cons = PlasmoModelGraph.getlinkconstraints(node)[graph]  #link constraints specific to this subproblem
        ncols2 = s.numCols  #this is normally the number of second stage variables, but I'm lifting so....
        nrows2 = length(s.linconstr) + length(link_cons)
        break #this only runs the loop once?
    end

    # set scenario indices for each MPI processor
    if dsp.comm_size > 1
        ncols2 = MPI.allreduce([ncols2], MPI.MAX, dsp.comm)[1]
        nrows2 = MPI.allreduce([nrows2], MPI.MAX, dsp.comm)[1]
    end

    #Need to include linkconstraints on subproblems in dimensions
    # println()
    # @show nscen
    # @show ncols1
    # @show ncols2
    # @show nrows1
    # @show nrows2

    DspCInterface.@dsp_ccall("setNumberOfScenarios", Nothing, (Ptr{Nothing}, Cint), dsp.p, convert(Cint, nscen))
    DspCInterface.@dsp_ccall("setDimensions", Nothing,
        (Ptr{Nothing}, Cint, Cint, Cint, Cint),
        dsp.p, convert(Cint, ncols1), convert(Cint, nrows1), convert(Cint, ncols2), convert(Cint, nrows2))

    # get problem data
    # this is the master model
    start, index, value, clbd, cubd, ctype, obj, rlbd, rubd = getDataFormat(master)

    # println()
    # @show start
    # @show index
    # @show value
    # @show clbd
    # @show cubd
    # @show ctype
    # @show obj
    # @show rlbd
    # @show rubd

    DspCInterface.@dsp_ccall("loadFirstStage", Nothing,
        (Ptr{Nothing}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble},
            Ptr{Cdouble}, Ptr{UInt8}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
            dsp.p, start, index, value, clbd, cubd, ctype, obj, rlbd, rubd)

    for id in dsp.block_ids
        # model and probability
        #blk = blocks.children[id] #the sub model
        blk = subproblems[id]  #need to be careful with block ids here
        node = PlasmoModelGraph.getnode(blk)
        #probability = blocks.weight[id] #probability of this scenario
        probability = probabilities[id]
        linkcons = PlasmoModelGraph.getlinkconstraints(node)[graph]
        # get model data
        start, index, value, clbd, cubd, ctype, obj, rlbd, rubd = getDataFormat(master,blk,linkcons)

        # println()
        # @show start
        # @show index
        # @show value
        # @show clbd
        # @show cubd
        # @show ctype
        # @show obj
        # @show typeof(obj)
        # @show rlbd
        # @show rubd
        # @show typeof(rlbd)

        #@show probability

        DspCInterface.@dsp_ccall("loadSecondStage", Nothing,
            (Ptr{Nothing}, Cint, Cdouble, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble},
                Ptr{Cdouble}, Ptr{UInt8}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
            dsp.p, id-1, probability, start, index, value, clbd, cubd, ctype, obj, rlbd, rubd)
    end
    #println(dsp)
end

#TODO
function loadDeterministicProblem(dsp::DspCInterface.DspModel, model::JuMP.Model)
    ncols = convert(Cint, model.numCols)
    nrows = convert(Cint, length(model.linconstr))
    start, index, value, clbd, cubd, ctype, obj, rlbd, rubd = getDataFormat(model)
    numels = length(index)
    DspCInterface.@dsp_ccall("loadDeterministic", Nothing,
        (Ptr{Nothing}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Cint, Cint, Cint,
            Ptr{Cdouble}, Ptr{Cdouble}, Ptr{UInt8}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
            dsp.p, start, index, value, numels, ncols, nrows, clbd, cubd, ctype, obj, rlbd, rubd)
end

#This is used in getDSPSolution
function getNumBlockCols(dsp::DspCInterface.DspModel, blocks::Dict{Int,JuMP.Model})
    DspCInterface.check_problem(dsp)
    #blocks = m.ext[:DspBlocks].children
    #Might need to be smarter about his
    #blocks = collect(zip(1:length(subproblems),subproblems))
    # get number of block columns
    numBlockCols = Dict{Int,Int}()  #{block:num_cols}
    if dsp.comm_size > 1
        num_proc_blocks = convert(Vector{Cint}, MPI.Allgather(length(blocks), dsp.comm))
        #@show num_proc_blocks
        #@show collect(keys(blocks))
        block_ids = MPI.Allgatherv(collect(keys(blocks)), num_proc_blocks, dsp.comm)
        #@show block_ids
        ncols_to_send = Int[blocks[i].numCols for i in keys(blocks)]
        #@show ncols_to_send
        ncols = MPI.Allgatherv(ncols_to_send, num_proc_blocks, dsp.comm)
        #@show ncols
        for i in 1:dsp.nblocks
            setindex!(numBlockCols, ncols[i], block_ids[i])
        end
    else
        for b in blocks
            setindex!(numBlockCols, b.second.numCols, b.first)
        end
    end
    return numBlockCols
end

#this is used for a master model
function getDataFormat(model::JuMP.Model)

    # Column wise sparse matrix
    mat = JuMP.prepConstrMatrix(model)
    # Tranpose; now I have row-wise sparse matrix
    #mat = sparse(Matrix(mat'))
    mat = permutedims(mat)

    # sparse description
    start = convert(Vector{Cint}, mat.colptr .- 1) #column
    index = convert(Vector{Cint}, mat.rowval .- 1) #row index
    value = mat.nzval

    # column type
    ctype = ""
    for i = 1:length(model.colCat)
        if model.colCat[i] == :Int
            ctype = ctype * "I";
        elseif model.colCat[i] == :Bin
            ctype = ctype * "B";
        else
            ctype = ctype * "C";
        end
    end
    ctype = unsafe_wrap(Vector{UInt8}, ctype)

    # objective coefficients
    obj = JuMP.prepAffObjective(model)
    rlbd, rubd = JuMP.prepConstrBounds(model)

    # set objective sense
    if model.objSense == :Max
        obj *= -1
    end

    return start, index, value, model.colLower, model.colUpper, ctype, obj, rlbd, rubd
end

#Get the data format for a child node.  We also need to pass the child's linkconstraints
function getDataFormat(master::JuMP.Model,child::JuMP.Model,linkcons::Vector{PlasmoModelGraph.LinkConstraint})
    # Column wise sparse matrix
    mat = prepChildConstrMatrix(master,child,linkcons)
    # Tranpose; now I have row-wise sparse matrix
    mat = permutedims(mat)

    # sparse description
    start = convert(Vector{Cint}, mat.colptr .- 1) #column
    index = convert(Vector{Cint}, mat.rowval .- 1) #row index
    value = mat.nzval

    # column type
    ctype = ""
    for i = 1:length(child.colCat)
        if child.colCat[i] == :Int
            ctype = ctype * "I";
        elseif child.colCat[i] == :Bin
            ctype = ctype * "B";
        else
            ctype = ctype * "C";
        end
    end
    ctype = unsafe_wrap(Vector{UInt8}, ctype)

    # objective coefficients
    obj = JuMP.prepAffObjective(child)

    linconstr = child.linconstr
    allconstr = [linconstr;linkcons]

    child_linear_lb = Array{Float64}(undef,0)
    child_linear_ub = Array{Float64}(undef,0)

    #rlbd, rubd = JuMP.prepConstrBounds(model)
    for con in allconstr
        push!(child_linear_lb,con.lb)
        push!(child_linear_ub,con.ub)
    end
    rlbd = child_linear_lb
    rubd = child_linear_ub

    # set objective sense
    if child.objSense == :Max
        obj *= -1
    end

    return start, index, value, child.colLower, child.colUpper, ctype, obj, rlbd, rubd
end

function prepChildConstrMatrix(master::JuMP.Model,child::JuMP.Model,linkconstraints::Vector{PlasmoModelGraph.LinkConstraint})
    rind = Int[]
    cind = Int[]
    value = Float64[]
    linconstr = child.linconstr              #these should be local model constraints
    #linkconstr = getlinkconstraints(node)[graph]#the constraints that connect this model to the parent
    allconstr = [linconstr;linkconstraints]
    for (nrow,con) in enumerate(allconstr)
        aff = con.terms
        for (var,id) in zip(reverse(aff.vars), length(aff.vars):-1:1)
            push!(rind, nrow)  #each variable has a row index
            if allconstr[nrow].terms.vars[id].m == master   #if the variable belongs to the master, use the master column index
                push!(cind, var.col)
            elseif allconstr[nrow].terms.vars[id].m == child          #if it's a local variable to the subproblem
                push!(cind, master.numCols + var.col)          #push the variable passed the master column
            end
            push!(value, aff.coeffs[id])
            # splice!(aff.vars, id)   #remove the variable and coefficient from aff
            # splice!(aff.coeffs, id)
        end
    end
    return sparse(rind, cind, value, length(allconstr), master.numCols + child.numCols)
end

function parseStatusCode(statcode::Integer)
    stat = :NotSolved
    if statcode == 3000
        stat = :Optimal
    elseif statcode == 3001
        stat = :Infeasible
    elseif statcode == 3002
        stat = :Unbounded
    elseif statcode in [3004,3007,3010]
        stat = :IterOrTimeLimit
    elseif statcode == 3005
        stat = :GapTolerance
    elseif statcode == 3006
        stat = :NodeLimit
    elseif statcode in [3008,3009,3014,3015,3016]
        stat = :UserLimit
    elseif statcode in [3011,3012,3013,3999]
        stat = :Error
    else
        stat = :Unknown
        @warn("Unknown status: $statcode")
    end

    stat
end

function getDspSolution(dspmodel,master,subproblems)
    dspmodel.primVal = DspCInterface.getPrimalBound(dspmodel)
    dspmodel.dualVal = DspCInterface.getDualBound(dspmodel)
    #println(dspmodel.primVal)
    if dspmodel.solve_type == :Dual
        dspmodel.rowVal = DspCInterface.getDualSolution(dspmodel)
        #Need to get primal solution
    else
        dspmodel.colVal = DspCInterface.getSolution(dspmodel)
    #    println(dspmodel.colVal)
        if master != nothing
            # parse solution to each block
            #start with the master
            n_start = 1
            n_end = master.numCols
            master.colVal = dspmodel.colVal[n_start:n_end]
            n_start += master.numCols
            #if haskey(m.ext, :DspBlocks) == true
            #blocks = m.ext[:DspBlocks].children
            blocks = Dict(zip(1:length(subproblems),subproblems))
            numBlockCols = getNumBlockCols(dspmodel,blocks)
            for i in 1:dspmodel.nblocks
                n_end += numBlockCols[i]
                if haskey(blocks, i)
                    # @show b
                    # @show n_start
                    # @show n_end
                    blocks[i].colVal = dspmodel.colVal[n_start:n_end]
                end
                n_start += numBlockCols[i]
            end
        end
    end
    #set the graph objective
    if master != nothing
        master.objVal = dspmodel.primVal
        # maximization?
        if master.objSense == :Max
            master.objVal *= -1
            dspmodel.primVal *= -1
            dspmodel.dualVal *= -1
        end
    end
end

end #end module


text:
module MinCut

import TreeLas.Cxx: lib

struct WeightedArc{I, F}
    head::I
    tail::I
    weight::F
end

WeightedArc{I, F}(t::Tuple{I2,I2,F2}) where {I, F, I2, F2} =
    WeightedArc{I,F}(I(t[1]), I(t[2]), F(t[3]))

function cxx_min_cut!(
    cut::Vector{Bool},
    n::Integer,
    warcs::Vector{WeightedArc{Cint, Cdouble}},
    source::Integer,
    target::Integer,
)
    resize!(cut, n)
    m = length(warcs)
    return ccall(
        (:min_cut_f64_i32, lib),
        Cdouble,
        (
            Ref{Bool},
            Csize_t,
            Csize_t,
            Ref{WeightedArc{Cint, Cdouble}},
            Cint,
            Cint,
        ),
        cut,
        n,
        m,
        warcs,
        source,
        target,
    )
end

end

module TestMinCut

import ..MinCut
using Test

@testset "MinCut" begin
    @testset "wiki1" begin
        warcs = MinCut.WeightedArc{Cint, Cdouble}.([
            (0, 1, 5.0),
            (0, 3, 10),
            (1, 2, 10),
            (3, 1, 15),
            (3, 2, 5),
        ])
        source = 0
        target = 2
        n = 4
        cut = Vector{Bool}()
        MinCut.cxx_min_cut!(cut, n, warcs, source, target)
        @test Set(i for i=1:n if cut[i]) == Set(1)
    end
    @testset "chengw1005" begin
        warcs = MinCut.WeightedArc{Cint, Cdouble}.([
            (0, 2, 13),
            (0, 1, 16),
            (1, 2, 10),
            (1, 3, 12),
            (2, 1, 4),
            (2, 4, 14),
            (3, 2, 9),
            (3, 5, 20),
            (4, 3, 7),
            (4, 5, 4),
        ])
        source, target = 0, 5
        n = 6
        cut = Vector{Bool}()
        MinCut.cxx_min_cut!(cut, n, warcs, source, target)
        @test Set(i for i=1:n if cut[i]) == Set(1 .+ [0, 1, 2, 4])
    end
end

end

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
module hill_decay_diffeq
using FinEtools
using FinEtoolsHeatDiff
using FinEtoolsHeatDiff.AlgoHeatDiffModule
import LinearAlgebra: cholesky
using UnicodePlots
using DifferentialEquations

function hill_decay_t3()
	thermal_conductivity =  [i==j ? 0.2 : zero(FFlt) for i=1:2, j=1:2]; # conductivity matrix
	Width = 60.0
	Height = 40.0
	N = 50
	specific_heat = 1.0
	dt = 2.0; # time step
	tend = 20*dt; # length of the time interval
	T0(x, y) = 500.0/(x^2+y^2+5);# Initial distribution of temperature

	tolerance =Width/N/100;

    fens,fes =T3block(Width, Height, N, N)
    for i in 1:count(fens)
    	fens.xyz[i, 1] -= Width / 2
    	fens.xyz[i, 2] -= Height / 2
    end
    l1  = selectnode(fens; box=[Width/2 Width/2 Height/2 Height/2],  inflate = tolerance)

    geom = NodalField(fens.xyz)
    Temp = NodalField(zeros(size(fens.xyz,1),1))

    for i in 1:count(fens)
    	Temp.values[i] = T0(fens.xyz[i, :]...)
    end

    applyebc!(Temp)
    numberdofs!(Temp)

    cornerdof = Temp.dofnums[l1]

    material = MatHeatDiff(thermal_conductivity, specific_heat)
    femm = FEMMHeatDiff(IntegDomain(fes, TriRule(3)), material)

    # The ODE to solve is C*dT/dt + K*T = 0
    K = conductivity(femm, geom, Temp)
    C = capacity(femm, geom, Temp)
    CF = cholesky(C)

    Tn = gathersysvec(Temp)

    tstart = time()
	f(dT, T, p, t) = begin dT .= -(CF \ (K * T)); end
	tspan = (0.0, tend)
	prob = ODEProblem(f, Tn, tspan)
	sol = solve(prob, KenCarp4(autodiff = false), abstol=1.0e-4, reltol=1.0e-4)
	# @show sol
	println("Time = $(time() - tstart)")

	# This is postprocessing  to extract the data for the plot.
	ts = Float64[]
	Corner_T = Float64[];
	for t in sol.t
		T = sol(t)
		push!(ts, t)
		push!(Corner_T, T[cornerdof][1])
	end

	@show ts
	@show Corner_T
	@show minimum(Corner_T), maximum(Corner_T)
	plt = lineplot(vec(ts), vec(Corner_T), canvas = DotCanvas, title = "Transient temperature at the corner", name = "T", xlabel = "Time", ylabel = "T")
	display(plt)

    true
end

end # module hill_decay_diffeq

using .hill_decay_diffeq
hill_decay_diffeq.hill_decay_t3()

max_stars_repo_name: PetrKryslUCSD/FinEtoolsHeatDiff

max_issues_repo_name: PetrKryslUCSD/FinEtoolsHeatDiff

max_forks_repo_name: PetrKryslUCSD/FinEtoolsHeatDiff


text:
# color vision deficiency (CVD) -  (color blind) - friendly
# https://zenodo.org/record/3381072
# Creative Commons Attribution 4.0 International
# 8, 12, and 15 color palettes for deuteranopia color blindness
# Martin Krzywinski martink@bcgsc.ca
# Methods and details: http://mkweb.bcgsc.ca/colorblind

#From Paul Tol: https://personal.sron.nl/~pault/
loadcolorscheme(:tol_bright, [
    colorant"#EE6677",
    colorant"#228833",
    colorant"#4477AA",
    colorant"#CCBB44",
    colorant"#66CCEE",
    colorant"#AA3377",
    colorant"#BBBBBB"
], "cvd", "cvd friendly, color blind friendly")

loadcolorscheme(:tol_muted, [
    colorant"#88CCEE",
    colorant"#44AA99",
    colorant"#117733",
    colorant"#332288",
    colorant"#DDCC77",
    colorant"#999933",
    colorant"#CC6677",
    colorant"#882255",
    colorant"#AA4499",
    colorant"#DDDDDD"
], "cvd", "cvd friendly, color blind friendly")

loadcolorscheme(:tol_light, [
    colorant"#BBCC33",
    colorant"#AAAA00",
    colorant"#77AADD",
    colorant"#EE8866",
    colorant"#EEDD88",
    colorant"#FFAABB",
    colorant"#99DDFF",
    colorant"#44BB99",
    colorant"#DDDDDD"
], "cvd", "cvd friendly, color blind friendly")

loadcolorscheme(:tol_sunset, [
    colorant"#364B9A",
    colorant"#4A7BB7",
    colorant"#6EA6CD",
    colorant"#98CAE1",
    colorant"#C2E4EF",
    colorant"#EAECCC",
    colorant"#FEDA8B",
    colorant"#FDB366",
    colorant"#F67E4B",
    colorant"#DD3D2D",
    colorant"#A50026",
], "cvd", "cvd friendly, color blind friendly")

loadcolorscheme(:tol_bu_rd, [
    colorant"#2166AC",
    colorant"#4393C3",
    colorant"#92C5DE",
    colorant"#D1E5F0",
    colorant"#F7F7F7",
    colorant"#FDDBC7",
    colorant"#F4A582",
    colorant"#D6604D",
    colorant"#B2182B",
], "cvd", "cvd friendly, color blind friendly")

loadcolorscheme(:tol_prgn, [
    colorant"#762A83",
    colorant"#9970AB",
    colorant"#C2A5CF",
    colorant"#E7D4E8",
    colorant"#F7F7F7",
    colorant"#D9F0D3",
    colorant"#ACD39E",
    colorant"#5AAE61",
    colorant"#1B7837",
], "cvd", "cvd friendly, color blind friendly")

loadcolorscheme(:tol_ylorbr, [
    colorant"#FFFFE5",
    colorant"#FFF7BC",
    colorant"#FEE391",
    colorant"#FEC44F",
    colorant"#FB9A29",
    colorant"#EC7014",
    colorant"#CC4C02",
    colorant"#993404",
    colorant"#662506",
], "cvd", "cvd friendly, color blind friendly")

loadcolorscheme(:ground_cover, [
    colorant"#5566AA",
    colorant"#117733",
    colorant"#44AA66",
    colorant"#55AA22",
    colorant"#668822",
    colorant"#99BB55",
    colorant"#558877",
    colorant"#88BBAA",
    colorant"#AADDCC",
    colorant"#44AA88",
    colorant"#DDCC66",
    colorant"#FFDD44",
    colorant"#FFEE88",
    colorant"#BB0011",
], "cvd", "cvd friendly, color blind friendly")

loadcolorscheme(:tol_YlOrBr, [
    colorant"#FFFFE5",
    colorant"#FFF7BC",
    colorant"#FEE391",
    colorant"#FEC44F",
    colorant"#FB9A29",
    colorant"#EC7014",
    colorant"#CC4C02",
    colorant"#993404",
    colorant"#662506",
], "cvd", "cvd friendly, color blind friendly")

loadcolorscheme(:Iridescent, [
    colorant"#FEFBE9",
    colorant"#FCF705",
    colorant"#F5F3C1",
    colorant"#EAF0B5",
    colorant"#DDECBF",
    colorant"#D0E7CA",
    colorant"#C2E3D2",
    colorant"#B5DDD8",
    colorant"#A8D8DC",
    colorant"#9BD2E1",
    colorant"#8DCBE4",
    colorant"#81C4E7",
    colorant"#7BBCE7",
    colorant"#7EB2E4",
    colorant"#88A5DD",
    colorant"#9398D2",
    colorant"#9B8AC4",
    colorant"#9D7DB2",
    colorant"#9A709E",
    colorant"#906388",
    colorant"#805770",
    colorant"#684957",
    colorant"#46353A",
], "cvd", "cvd friendly, color blind friendly")

loadcolorscheme(:tol_rainbow, [
    colorant"#E8ECFB",
    colorant"#D9CCE3",
    colorant"#D1BBD7",
    colorant"#CAACCB",
    colorant"#BA8BD4",
    colorant"#AE76A3",
    colorant"#AA6F9E",
    colorant"#994F88",
    colorant"#882E72",
    colorant"#1965B0",
    colorant"#437DBF",
    colorant"#5289C7",
    colorant"#6195CF",
    colorant"#7BAFDE",
    colorant"#4EB265",
    colorant"#90C987",
    colorant"#CAE0AB",
    colorant"#F7F056",
    colorant"#F7CB45",
    colorant"#F6C141",
    colorant"#F4A736",
    colorant"#F1932D",
    colorant"#EE8026",
    colorant"#E8601C",
    colorant"#E65518",
    colorant"#DC050C",
    colorant"#A5170E",
    colorant"#72190E",
    colorant"#42150A",
], "cvd", "cvd friendly, color blind friendly")

#From Color Universal Design (CUD): https://jfly.uni-koeln.de/color/
loadcolorscheme(:okabe_ito, [
    colorant"#E69F00",
    colorant"#56B4E9",
    colorant"#009E73",
    colorant"#F0E442",
    colorant"#0072B2",
    colorant"#D55E00",
    colorant"#CC79A7",
    colorant"#000000"
], "cvd", "cvd friendly, color blind friendly")

# 8, 12, and 15 color palettes for deuteranopia color blindness
# Martin Krzywinski martink@bcgsc.ca
# Methods and details: http://mkweb.bcgsc.ca/colorblind

loadcolorscheme(:mk_8, [
    colorant"#000000",
    colorant"#2271B2",
    colorant"#3DB7E9",
    colorant"#F748A5",
    colorant"#359B73",
    colorant"#d55e00",
    colorant"#e69f00",
    colorant"#f0e442",
], "cvd", "cvd friendly, color blind friendly")

loadcolorscheme(:mk_12, [
    colorant"#9F0162",
    colorant"#009F81",
    colorant"#FF5AAF",
    colorant"#00FCCF",
    colorant"#8400CD",
    colorant"#008DF9",
    colorant"#00C2F9",
    colorant"#FFB2FD",
    colorant"#A40122",
    colorant"#E20134",
    colorant"#FF6E3A",
    colorant"#FFC33B",
], "cvd", "cvd friendly, color blind friendly")

loadcolorscheme(:mk_15, [
    colorant"#68023F",
    colorant"#008169",
    colorant"#EF0096",
    colorant"#00DCB5",
    colorant"#FFCFE2",
    colorant"#003C86",
    colorant"#9400E6",
    colorant"#009FFA",
    colorant"#FF71FD",
    colorant"#7CFFFA",
    colorant"#6A0213",
    colorant"#008607",
    colorant"#F60239",
    colorant"#00E307",
    colorant"#FFDC3D",
], "cvd", "cvd friendly, color blind friendly")