text:
using KernelAbstractions: @index, @kernel
using KernelAbstractions.Extras.LoopInfo: @unroll
using Oceananigans.Architectures: device_event, architecture
using Oceananigans.Utils: launch!
using Oceananigans.Grids
using Oceananigans.Operators: Δzᶜᶜᶜ

@kernel function _compute_light_growth!(light, grid, h, P, light_function, light_growth_function, simulation_time,
                                         average_code, shading, chl2c, Kc)
    i, j = @index(Global, NTuple)

    h_ij = @inbounds h[i, j]

    chlinteg = 0

    @unroll for k in grid.Nz : -1 : 1
        if shading
            chlinteg = chlinteg + @inbounds P[i, j, k] * chl2c * Δzᶜᶜᶜ(i, j, k, grid)
        end
        z_center = znode(Center(), Center(), Center(), i, j, k, grid)

        if average_code==1
            @inbounds light[i, j, k] = light_growth_function(light_function(simulation_time, z_center)/exp(chlinteg*Kc))
        elseif (average_code==0)|(average_code==2)
            @inbounds light[i, j, k] = light_function(simulation_time, z_center)/exp(chlinteg*Kc)
        end

    end

    if ~(average_code==0)
        light_sum = 0
        dz_sum = 0

        @unroll for k in grid.Nz : -1 : 1 # scroll from surface to bottom
            z_center = znode(Center(), Center(), Center(), i, j, k, grid)

            if z_center > - h_ij

                Δz_ijk = Δzᶜᶜᶜ(i, j, k, grid)

                light_sum = light_sum + @inbounds light[i, j, k] * Δz_ijk
                dz_sum = dz_sum + Δz_ijk

            end
        end

        @unroll for k in grid.Nz : -1 : 1 # scroll to point just above the bottom
            z_center = znode(Center(), Center(), Center(), i, j, k, grid)

            if z_center > - h_ij
                @inbounds light[i, j, k] = light_sum/dz_sum
            end
        end
    end

    if average_code!=:1
        @unroll for k in grid.Nz : -1 : 1
            @inbounds light[i, j, k] = light_growth_function(light[i, j, k])
        end
    end

end

"""
    LightGrowth!(light, h, P, light_function, light_growth_function, simulation_time, average=nothing, shading=false, chl2c=0.02, Kc = 0.041)

Compute and update the light growth for Phytoplankton.
"""
function LightGrowth!(light, h, P, light_function, light_growth_function, simulation_time,
                      average=nothing, shading=false, chl2c=0.02, Kc = 0.041)
    grid = h.grid
    arch = architecture(grid)

    if average==nothing
        average_code = 0
    elseif average==:growth
        average_code = 1
    elseif average==:light
        average_code = 2
    else
        error("Unsupported average ($average).")
    end

    event = launch!(arch, grid, :xy,
                    _compute_light_growth!, light, grid, h, P, light_function, light_growth_function,
                    simulation_time, average_code, shading, chl2c, Kc,
                    dependencies = device_event(arch))

    wait(device_event(arch), event)

    return nothing
end

max_stars_repo_name: iuryt/biocenanigans

max_issues_repo_name: iuryt/biocenanigans


text:
using LinearAlgebra, MATLAB, ForwardDiff, SparseArrays, Printf

include("/Users/kevintracy/.julia/dev/CPEG/src/qp_solver.jl")

ext: jl


text:
type Node
  value
  parent::Nullable
  children::Array
end

# Outer Constructor Method for convenience
Node(v) = Node(v, Nullable{Node}(), Array(Node, 0))

# build tree

function add_children!(parent::Node, children::Array)
  for child in children
    child.parent = Nullable(parent)
    push!(parent.children, child)
  end
end

# get Node

function get_parent(node::Node)
  if !isnull(node.parent)
    get(node.parent)
  end
end

function get_root(node::Node)
  if is_root(node)
    return node
  else
    get_root(get_parent(node))
  end
end

# helpers

function is_leaf(node::Node)
  length(node.children) <= 0
end

function is_root(node::Node)
  isnull(node.parent)
end

# check if 1st node (possible_ancestor) is ancestor (older) of 2nd node (possible_descendant)
function is_ancestor(possible_ancestor::Node, possible_descendant::Node)
  if isnull(possible_descendant.parent)
    return false
  elseif get(possible_descendant.parent) == possible_ancestor
    return true
  else
    is_ancestor(possible_ancestor, get(possible_descendant.parent))
  end
end

# check if 1st node (possible_descendant) is descendant (younger) of 2nd node (possible_ancestor)
function is_descendant(possible_descendant::Node, possible_ancestor::Node)
  is_ancestor(possible_ancestor, possible_descendant)
end

# check if 1st node (possible_sibling) is a sibling (same parents) of 2nd node (check_node)
function is_sibling(possible_sibling::Node, check_node::Node)
  if isnull(possible_sibling.parent)
    return false
  else
    parent = get(possible_sibling.parent)
    if length(parent.children) <= 0
      return false
    else
      return in(check_node, parent.children)
    end
  end
end

# get distance from node to root (top)
function depth(node::Node, count=0)
  if isnull(node.parent)
    return count
  else
    depth(get_parent(node), count+1)
  end
end

# get longest distance from node to longest child/grandchild (bottom)
function height(node::Node)
  if length(node.children) <= 0
    return 0
  else
    return 1 + maximum( map(height, node.children) )
  end
end


text:
function State(c::AbstractArray, t::Timestep, iₚ::Integer)
    size = length(c)
    State{size}(SVector{size}(c), t, iₚ)
end

function State(m::PMDP, c::AbstractArray, t::Timestep, product::AbstractArray)
    size = length(c)
    product==empty_product(m) ? iₚ = empty_product_id(m) : iₚ = findfirst(x-> x.res==product, products(m))
    State{size}(SVector{size}(c), t, iₚ)
end

function product(m::PMDP, s::State)
    s.iₚ==empty_product_id(m) ? empty_product(m) : products(m)[s.iₚ]
end

function show(io::IO, s::State)
    print(io, "c:$(s.c)_t:$(s.t)_p:$(s.iₚ)")
end


text:
module RulesMvNormalMeanCovarianceMeanTest

using Test
using ReactiveMP
using Random

import ReactiveMP: @test_rules

@testset "rules:MvNormalMeanCovariance:mean" begin

    @testset "Belief Propagation: (m_out::PointMass, m_Σ::PointMass)" begin

        @test_rules [ with_float_conversions = true ] MvNormalMeanCovariance(:μ, Marginalisation) [
            (input = (m_out = PointMass([ -1.0 ]), m_Σ = PointMass([ 2.0 ])), output = MvNormalMeanCovariance([ -1.0 ], [ 2.0 ])),
            (input = (m_out = PointMass([  1.0 ]), m_Σ = PointMass([ 2.0 ])), output = MvNormalMeanCovariance([  1.0 ], [ 2.0 ])),
            (input = (m_out = PointMass([ 2.0 ]),  m_Σ = PointMass([ 1.0 ])), output = MvNormalMeanCovariance([  2.0 ], [ 1.0 ])),
            (input = (m_out = PointMass([ 1.0, 3.0 ]),  m_Σ = PointMass([ 3.0 2.0; 2.0 4.0 ])),   output = MvNormalMeanCovariance([ 1.0, 3.0 ], [ 3.0 2.0; 2.0 4.0 ])),
            (input = (m_out = PointMass([ -1.0, 2.0 ]), m_Σ = PointMass([ 7.0 -1.0; -1.0 9.0 ])), output = MvNormalMeanCovariance([ -1.0, 2.0 ], [ 7.0 -1.0; -1.0 9.0 ])),
            (input = (m_out = PointMass([ 0.0, 0.0 ]),  m_Σ = PointMass([ 1.0 0.0; 0.0 1.0 ])),   output = MvNormalMeanCovariance([ 0.0, 0.0 ], [ 1.0 0.0; 0.0 1.0 ]))
        ]

    end

    @testset "Belief Propagation: (m_out::MultivariateNormalDistributionsFamily, m_Σ::PointMass)" begin

        @test_rules [ with_float_conversions = true ] MvNormalMeanCovariance(:μ, Marginalisation) [
            (input = (m_out = MvNormalMeanPrecision([ 2.0, 1.0 ], [ 3.0 2.0; 2.0 4.0 ]),   m_Σ = PointMass([ 6.0 4.0; 4.0 8.0 ])),    output = MvNormalMeanCovariance([ 2.0, 1.0 ], [ 13/2 15/4; 15/4 67/8 ])),
            (input = (m_out = MvNormalMeanPrecision([ 0.0, 0.0 ], [ 7.0 -1.0; -1.0 9.0 ]), m_Σ = PointMass([ 12.0 -2.0; -2.0 7.0 ])), output = MvNormalMeanCovariance([ 0.0, 0.0 ], [ 753/62 -123/62; -123/62 441/62 ])),
            (input = (m_out = MvNormalMeanPrecision([ 3.0, -1.0 ], [ 1.0 0.0; 0.0 1.0 ]),  m_Σ = PointMass([ 1.0 0.0; 0.0 1.0 ])),    output = MvNormalMeanCovariance([ 3.0, -1.0 ], [ 2.0 0.0; 0.0 2.0 ])),

        ]

        @test_rules [ with_float_conversions = true ] MvNormalMeanCovariance(:μ, Marginalisation) [
            (input = (m_out = MvNormalMeanCovariance([ 2.0, 1.0 ], [ 3.0 2.0; 2.0 4.0 ]),  m_Σ = PointMass([ 6.0 4.0; 4.0 8.0 ])),     output = MvNormalMeanCovariance([ 2.0, 1.0 ], [ 9.0 6.0; 6.0 12.0 ])),
            (input = (m_out = MvNormalMeanCovariance([ 0.0, 0.0 ], [ 7.0 -1.0; -1.0 9.0 ]), m_Σ = PointMass([ 12.0 -2.0; -2.0 7.0 ])), output = MvNormalMeanCovariance([ 0.0, 0.0 ], [ 19.0 -3.0; -3.0 16.0 ])),
            (input = (m_out = MvNormalMeanCovariance([ 3.0, -1.0 ], [ 1.0 0.0; 0.0 1.0 ]),  m_Σ = PointMass([ 1.0 0.0; 0.0 1.0 ])),    output = MvNormalMeanCovariance([ 3.0, -1.0 ], [ 2.0 0.0; 0.0 2.0 ]))
        ]

        @test_rules [ with_float_conversions = true ] MvNormalMeanCovariance(:μ, Marginalisation) [
            (input = (m_out = MvNormalWeightedMeanPrecision([ 2.0, 1.0 ], [ 3.0 2.0; 2.0 4.0 ]),   m_Σ = PointMass([ 6.0 4.0; 4.0 8.0 ])),    output = MvNormalMeanCovariance([ 3/4, -1/8 ], [ 13/2 15/4; 15/4 67/8 ])),
            (input = (m_out = MvNormalWeightedMeanPrecision([ 0.0, 0.0 ], [ 7.0 -1.0; -1.0 9.0 ]), m_Σ = PointMass([ 12.0 -2.0; -2.0 7.0 ])), output = MvNormalMeanCovariance([ 0.0, 0.0 ], [ 753/62 -123/62; -123/62 441/62 ])),
            (input = (m_out = MvNormalWeightedMeanPrecision([ 3.0, -1.0 ], [ 1.0 0.0; 0.0 1.0 ]),  m_Σ = PointMass([ 1.0 0.0; 0.0 1.0 ])),    output = MvNormalMeanCovariance([ 3.0, -1.0 ], [ 2.0 0.0; 0.0 2.0 ])),
        ]

    end

    @testset "Variational: (q_out::Any, q_Σ::Any)" begin

        @test_rules [ with_float_conversions = true ] MvNormalMeanCovariance(:μ, Marginalisation) [
            (input = (q_out = PointMass([ -1.0 ]), q_Σ = PointMass([ 2.0 ])), output = MvNormalMeanCovariance([ -1.0 ], [ 2.0 ])),
            (input = (q_out = PointMass([  1.0 ]), q_Σ = PointMass([ 2.0 ])), output = MvNormalMeanCovariance([  1.0 ], [ 2.0 ])),
            (input = (q_out = PointMass([ 2.0 ]),  q_Σ = PointMass([ 1.0 ])), output = MvNormalMeanCovariance([  2.0 ], [ 1.0 ])),
            (input = (q_out = PointMass([ 1.0, 3.0 ]),  q_Σ = PointMass([ 3.0 2.0; 2.0 4.0 ])),   output = MvNormalMeanCovariance([ 1.0, 3.0 ], [ 3.0 2.0; 2.0 4.0 ])),
            (input = (q_out = PointMass([ -1.0, 2.0 ]), q_Σ = PointMass([ 7.0 -1.0; -1.0 9.0 ])), output = MvNormalMeanCovariance([ -1.0, 2.0 ], [ 7.0 -1.0; -1.0 9.0 ])),
            (input = (q_out = PointMass([ 0.0, 0.0 ]),  q_Σ = PointMass([ 1.0 0.0; 0.0 1.0 ])),   output = MvNormalMeanCovariance([ 0.0, 0.0 ], [ 1.0 0.0; 0.0 1.0 ]))
        ]

        @test_rules [ with_float_conversions = true ] MvNormalMeanCovariance(:μ, Marginalisation) [
            (input = (q_out = MvNormalMeanCovariance([ 2.0, 1.0 ], [ 3.0 2.0; 2.0 4.0 ]), q_Σ = PointMass([ 6.0 4.0; 4.0 8.0 ])), output = MvNormalMeanCovariance([ 2.0, 1.0 ], [ 6.0 4.0; 4.0 8.0 ])),
            (input = (q_out = MvNormalMeanPrecision([ 2.0, 1.0 ], [ 3.0 2.0; 2.0 4.0 ]), q_Σ = PointMass([ 6.0 4.0; 4.0 8.0 ])), output = MvNormalMeanCovariance([ 2.0, 1.0 ], [ 6.0 4.0; 4.0 8.0 ])),
            (input = (q_out = MvNormalWeightedMeanPrecision([ 2.0, 1.0 ], [ 3.0 2.0; 2.0 4.0 ]), q_Σ = PointMass([ 6.0 4.0; 4.0 8.0 ])), output = MvNormalMeanCovariance([ 3/4, -1/8 ], [ 6.0 4.0; 4.0 8.0 ])),
        ]

    end

    @testset "Structured variational: (m_out::PointMass, q_Σ::Any)" begin

        @test_rules [ with_float_conversions = true ] MvNormalMeanCovariance(:μ, Marginalisation) [
            (input = (m_out = PointMass([ 1.0, 3.0 ]),  q_Σ = PointMass([ 3.0 2.0; 2.0 4.0 ])),   output = MvNormalMeanCovariance([ 1.0, 3.0 ], [ 3.0 2.0; 2.0 4.0 ])),
            (input = (m_out = PointMass([ -1.0, 2.0 ]), q_Σ = PointMass([ 7.0 -1.0; -1.0 9.0 ])), output = MvNormalMeanCovariance([ -1.0, 2.0 ], [ 7.0 -1.0; -1.0 9.0 ])),
            (input = (m_out = PointMass([ 0.0, 0.0 ]),  q_Σ = PointMass([ 1.0 0.0; 0.0 1.0 ])),   output = MvNormalMeanCovariance([ 0.0, 0.0 ], [ 1.0 0.0; 0.0 1.0 ])),
        ]

    end

    @testset "Structured variational: (m_out::MultivariateNormalDistributionsFamily, q_Σ::Any)" begin

        @test_rules [ with_float_conversions = true ] MvNormalMeanCovariance(:μ, Marginalisation) [
            (input = (m_out = MvNormalMeanPrecision([ 2.0, 1.0 ], [ 3.0 2.0; 2.0 4.0 ]),   q_Σ = PointMass([ 6.0 4.0; 4.0 8.0 ])),    output = MvNormalMeanCovariance([ 2.0, 1.0 ], [ 13/2 15/4; 15/4 67/8 ])),
            (input = (m_out = MvNormalMeanPrecision([ 0.0, 0.0 ], [ 7.0 -1.0; -1.0 9.0 ]), q_Σ = PointMass([ 12.0 -2.0; -2.0 7.0 ])), output = MvNormalMeanCovariance([ 0.0, 0.0 ], [ 753/62 -123/62; -123/62 441/62 ])),
            (input = (m_out = MvNormalMeanPrecision([ 3.0, -1.0 ], [ 1.0 0.0; 0.0 1.0 ]),  q_Σ = PointMass([ 1.0 0.0; 0.0 1.0 ])),    output = MvNormalMeanCovariance([ 3.0, -1.0 ], [ 2.0 0.0; 0.0 2.0 ])),
        ]

        @test_rules [ with_float_conversions = true ] MvNormalMeanCovariance(:μ, Marginalisation) [
            (input = (m_out = MvNormalMeanCovariance([ 2.0, 1.0 ], [ 3.0 2.0; 2.0 4.0 ]),  q_Σ = PointMass([ 6.0 4.0; 4.0 8.0 ])),     output = MvNormalMeanCovariance([ 2.0, 1.0 ], [ 9.0 6.0; 6.0 12.0 ])),
            (input = (m_out = MvNormalMeanCovariance([ 0.0, 0.0 ], [ 7.0 -1.0; -1.0 9.0 ]), q_Σ = PointMass([ 12.0 -2.0; -2.0 7.0 ])), output = MvNormalMeanCovariance([ 0.0, 0.0 ], [ 19.0 -3.0; -3.0 16.0 ])),
            (input = (m_out = MvNormalMeanCovariance([ 3.0, -1.0 ], [ 1.0 0.0; 0.0 1.0 ]),  q_Σ = PointMass([ 1.0 0.0; 0.0 1.0 ])),    output = MvNormalMeanCovariance([ 3.0, -1.0 ], [ 2.0 0.0; 0.0 2.0 ]))
        ]

        @test_rules [ with_float_conversions = true ] MvNormalMeanCovariance(:μ, Marginalisation) [
            (input = (m_out = MvNormalWeightedMeanPrecision([ 2.0, 1.0 ], [ 3.0 2.0; 2.0 4.0 ]),   q_Σ = PointMass([ 6.0 4.0; 4.0 8.0 ])),    output = MvNormalMeanCovariance([ 3/4, -1/8 ], [ 13/2 15/4; 15/4 67/8 ])),
            (input = (m_out = MvNormalWeightedMeanPrecision([ 0.0, 0.0 ], [ 7.0 -1.0; -1.0 9.0 ]), q_Σ = PointMass([ 12.0 -2.0; -2.0 7.0 ])), output = MvNormalMeanCovariance([ 0.0, 0.0 ], [ 753/62 -123/62; -123/62 441/62 ])),
            (input = (m_out = MvNormalWeightedMeanPrecision([ 3.0, -1.0 ], [ 1.0 0.0; 0.0 1.0 ]),  q_Σ = PointMass([ 1.0 0.0; 0.0 1.0 ])),    output = MvNormalMeanCovariance([ 3.0, -1.0 ], [ 2.0 0.0; 0.0 2.0 ])),
        ]

    end

end

end


text:
Base.@deprecate_binding StringVector Base.StringVector false

# PR #17302
# Provide a non-deprecated version of `@vectorize_(1|2)arg` macro which defines
# deprecated version of the function so that the depwarns can be fixed without
# breaking users.
# Packages are expected to use this to maintain the old API until all users
# of the deprecated vectorized function have migrated.
# These macros should be dropped when the support for `0.6` is dropped from `Compat`.
# Modified based on the version copied from 0.6 Base.
if VERSION < v"0.7.0-DEV.1211"
    macro dep_vectorize_1arg(S, f)
        S = esc(S)
        f = esc(f)
        T = esc(:T)
        x = esc(:x)
        AbsArr = esc(:AbstractArray)
        Base.depwarn("Implicit vectorized function is deprecated in favor of compact broadcast syntax.",
            Symbol("@dep_vectorize_1arg"))
        :(@deprecate $f{$T<:$S}($x::$AbsArr{$T}) @compat($f.($x)))
    end

    macro dep_vectorize_2arg(S, f)
        S = esc(S)
        f = esc(f)
        T1 = esc(:T1)
        T2 = esc(:T2)
        x = esc(:x)
        y = esc(:y)
        AbsArr = esc(:AbstractArray)
        Base.depwarn("Implicit vectorized function is deprecated in favor of compact broadcast syntax.",
            Symbol("@dep_vectorize_2arg"))
        quote
            @deprecate $f{$T1<:$S}($x::$S, $y::$AbsArr{$T1}) @compat($f.($x,$y))
            @deprecate $f{$T1<:$S}($x::$AbsArr{$T1}, $y::$S) @compat($f.($x,$y))
            @deprecate $f{$T1<:$S,$T2<:$S}($x::$AbsArr{$T1}, $y::$AbsArr{$T2}) @compat($f.($x,$y))
        end
    end
else
    macro dep_vectorize_1arg(S, f)
        AbstractArray = GlobalRef(Base, :AbstractArray)
        Base.depwarn("Implicit vectorized function is deprecated in favor of compact broadcast syntax.",
            Symbol("@dep_vectorize_1arg"))
        return esc(:(@deprecate $f(x::$AbstractArray{T}) where {T<:$S} $f.(x)))
    end

    macro dep_vectorize_2arg(S, f)
        AbstractArray = GlobalRef(Base, :AbstractArray)
        Base.depwarn("Implicit vectorized function is deprecated in favor of compact broadcast syntax.",
            Symbol("@dep_vectorize_2arg"))
        return esc(quote
            @deprecate $f(x::$S, y::$AbstractArray{T1}) where {T1<:$S} $f.(x, y)
            @deprecate $f(x::$AbstractArray{T1}, y::$S) where {T1<:$S} $f.(x, y)
            @deprecate $f(x::$AbstractArray{T1}, y::$AbstractArray{T2}) where {T1<:$S, T2<:$S} $f.(x, y)
        end)
    end
end

# compatibility with https://github.com/JuliaLang/julia/pull/26156
Base.@deprecate trunc(x, digits; base = 10) Compat.trunc(x, digits = digits, base = base) false
Base.@deprecate floor(x, digits; base = 10) Compat.floor(x, digits = digits, base = base) false
Base.@deprecate ceil(x, digits; base = 10) Compat.ceil(x, digits = digits, base = base) false
Base.@deprecate round(x, digits; base = 10) Compat.round(x, digits = digits, base = base) false
Base.@deprecate signif(x, digits; base = 10) Compat.round(x, sigdigits = digits, base = base) false

if VERSION >= v"1.1.0-DEV.506"
    # deprecation of range(start, stop) for earlier versions is done in Compat.jl
    # This method is restricted to Number, since we don't
    # want to overwrite the (::Any, ::Any) method in Base.
    function range(start::Number, stop::Number; kwargs...)
        rangedepwarn(;kwargs...)
        range(start; stop=stop, kwargs...)
    end
end

# this was defined for use with Julia versions prior to 0.5
# (see https://github.com/JuliaLang/Compat.jl/pull/316)
macro dotcompat(x)
    Base.depwarn("`@dotcompat x` is deprecated, use `@compat @. x` instead.", Symbol("@dotcompat"))
    esc(:(Compat.@compat @. $x))
end
export @dotcompat

# Compat.Random.uuid1, uuid4, uuid_version are deprecated in Compat.jl


text:
# sort results

using DelimitedFiles
using DistributedFactorGraphs

@show ARGS

data = readdlm(ARGS[1], ',')
Base.cp(ARGS[1], ARGS[1]*"_backup")

usPoses = Symbol.(data[:,1])

poses = usPoses |> sortDFG

sortperm = zeros(Int, length(poses))
for i in 1:length(poses)
  sortperm[i] = findfirst(x->x==poses[i], usPoses)
end

sortperm

sdata = data[sortperm, :]

writedlm(ARGS[1], sdata, ',')

#


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

include("dgsem_tree/dg.jl")
include("dgsem_structured/dg.jl")
include("dgsem_unstructured/dg.jl")
include("dgsem_p4est/dg.jl")
include("dg_common.jl")
include("fdsbp_tree/fdsbp_2d.jl")

end # @muladd


text:
using LCPSim
using RigidBodyDynamics
using JuMP
using Base.Test

srand(1)
@testset "linearized mechanism" begin
    @testset "truly linear mechanism" begin
        mechanism = rand_chain_mechanism(Float64, [Prismatic{Float64} for i = 1:5]...)
        x0 = MechanismState{Float64}(mechanism)
        x = MechanismState{Float64}(mechanism)
        for i in 1:10
            set_configuration!(x0, randn(num_positions(x0)))
            set_velocity!(x0, randn(num_velocities(x0)))
            set_configuration!(x, randn(num_positions(x0)))
            set_velocity!(x, randn(num_velocities(x0)))
            x_linear = LinearizedState(x0, x)
            @test center_of_mass(x) ≈ linearized(center_of_mass, x_linear)
            @test !(center_of_mass(x) ≈ center_of_mass(x0))

            for body in bodies(mechanism)
                frame = default_frame(body)
                @test transform_to_root(x, frame) ≈ linearized(x -> transform_to_root(x, frame), x_linear)

                p = Point3D(frame, randn(3))
                @test transform_to_root(x, frame) * p ≈ linearized(x -> transform_to_root(x, frame) * p, x_linear)
            end
        end
    end

    @testset "nonlinear mechanism" begin
        mechanism = rand_chain_mechanism(Float64, [QuaternionFloating{Float64}; [Revolute{Float64} for i = 1:5]]...)
        x0 = MechanismState{Float64}(mechanism)
        x = MechanismState{Float64}(mechanism)
        for i in 1:10
            set_configuration!(x0, randn(num_positions(x0)))
            set_velocity!(x0, randn(num_velocities(x0)))

            set_configuration!(x, configuration(x0) + 1e-5 * (rand(num_positions(x0)) .- 0.5))
            set_velocity!(x, velocity(x0) + 1e-5 * (rand(num_velocities(x0)) .- 0.5))

            x_linear =  LCPSim.Linear.LinearizedState(x0, x)

            @test center_of_mass(x) ≈ linearized(center_of_mass, x_linear)
            @test !(center_of_mass(x0) ≈ center_of_mass(x))  # make sure we've perturbed enough for linearization to help

            for (i, body) in enumerate(bodies(mechanism))
                frame = default_frame(body)
                @test transform_to_root(x, frame) ≈ linearized(x -> transform_to_root(x, frame), x_linear)

                if i > 1
                    @test !(transform_to_root(x0, frame) ≈ transform_to_root(x, frame))
                end

                p = Point3D(frame, randn(3))
                @test transform_to_root(x, frame) * p ≈ linearized(x -> transform_to_root(x, frame) * p, x_linear)
            end
        end
    end

    @testset "jump state" begin
        mechanism = rand_chain_mechanism(Float64, [Prismatic{Float64} for i = 1:5]...)
        m = Model()
        @variable m q[1:num_positions(mechanism)]
        @variable m v[1:num_velocities(mechanism)]
        setvalue(q, randn(length(q)))
        setvalue(v, randn(length(v)))
        x_var = MechanismState{Variable, Float64, AffExpr}(mechanism, q, v, eltype(q)[])
        x_linear = LinearizedState{Float64}(x_var)
        @test configuration(x_linear.linearization_state) == getvalue(q)
        @test velocity(x_linear.linearization_state) == getvalue(v)
    end
end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    pathgraph, cyclegraph,
    completebipartite, completegraph,
    laddergraph, circularladder, moebiusladder,
    squaregrid,
    generalizedpetersen, petersengraph, dodecahedralgraph

"""
    pathgraph(n::Int; mutable=false) -> PlainGraph

Generate path graph ``P_n``.
"""
function pathgraph(n::Int; mutable=false)
    n >= 2 || throw(DomainError(n, "n should be 2 or more"))
    f = mutable ? plaingraph : immutableplaingraph
    return f(n, (i, i + 1) for i in 1:(n-1))
end

"""
    cyclegraph(length::Int; mutable=false) -> PlainGraph

Generate cycle graph ``C_n``.
"""
function cyclegraph(n::Int; mutable=false)
    n >= 3 || throw(DomainError(n, "n should be 3 or more"))
    edges = [(i, i + 1) for i in 1:(n-1)]
    push!(edges, (n, 1))
    f = mutable ? plaingraph : immutableplaingraph
    return f(n, edges)
end

"""
    completebipartite(m::Int,n::Int; mutable=false) -> PlainGraph

Generate bipartite graph ``K_{m,n}``.
"""
function completebipartite(m::Int, n::Int; mutable=false)
    m >= 1 || throw(DomainError(m, "m should not be 1 or more"))
    n >= 1 || throw(DomainError(n, "n should not be 1 or more"))
    edges = Tuple{Int,Int}[]
    for i in 1:n
        for j in n+1:n+m
            push!(edges, (i, j))
        end
    end
    f = mutable ? plaingraph : immutableplaingraph
    return f(n + m, edges)
end

"""
    completegraph(length::Int; mutable=false) -> PlainGraph

Generate complete graph ``K_n``.
"""
function completegraph(n::Int; mutable=false)
    n >= 0 || throw(DomainError(n, "n should not be negative"))
    f = mutable ? plaingraph : immutableplaingraph
    n == 0 && return f()
    n == 1 && return f(1, [])
    return f(n, ((u, v) for (u, v) in combinations(n)))
end

"""
    laddergraph(n::Int; mutable=false) -> PlainGraph

Generate ladder graph ``L_n``.
"""
function laddergraph(n::Int; mutable=false)
    n >= 1 || throw(DomainError(n, "n should be 1 or more"))
    edges = [(1, 2)]
    for i in 1:n-1
        push!(edges, (2i+1, 2i+2))
        push!(edges, (2i-1, 2i+1))
        push!(edges, (2i, 2i+2))
    end
    f = mutable ? plaingraph : immutableplaingraph
    return f(n * 2, edges)
end

"""
    circularladder(n::Int; mutable=false) -> PlainGraph

Generate circular ladder graph ``CL_n``.
"""
function circularladder(n::Int; mutable=false)
    n >= 3 || throw(DomainError(n, "n should be 3 or more"))
    graph = laddergraph(n, mutable=true)
    addedge!(graph, 1, 2n-1)
    addedge!(graph, 2, 2n)
    if !mutable
        graph = immutableplaingraph(graph)
    end
    return graph
end

"""
    moebiusladder(n::Int; mutable=false) -> PlainGraph

Generate Möbius ladder graph ``ML_n``.
"""
function moebiusladder(n::Int; mutable=false)
    n >= 3 || throw(DomainError(n, "n should be 3 or more"))
    graph = laddergraph(n, mutable=true)
    addedge!(graph, 1, 2n)
    addedge!(graph, 2, 2n-1)
    if !mutable
        graph = immutableplaingraph(graph)
    end
    return graph
end

"""
    squaregrid(m::Int,n::Int; mutable=false) -> PlainGraph

Generate ``m \\times n`` square grid graph.

Use `cartesianproduct` for higher dimensional grid graphs.
"""
function squaregrid(m::Int, n::Int; mutable=false)
    m >= 1 || throw(DomainError(m, "m should not be 1 or more"))
    n >= 1 || throw(DomainError(n, "n should not be 1 or more"))
    edges = Tuple{Int,Int}[]
    for i in 1:(m*n-1)
        mod(i, n) == 0 && continue
        push!(edges, (i, i + 1))
    end
    for i in 1:((m-1)*n)
        push!(edges, (i, i + n))
    end
    f = mutable ? plaingraph : immutableplaingraph
    return f(n * m, edges)
end

"""
    generalizedpetersen(n::Int, k::Int; mutable=false) -> PlainGraph

Generate generalized petersen graph ``G(n,k)``.
"""
function generalizedpetersen(n::Int, k::Int; mutable=false)
    n >= 3 || throw(DomainError(n, "n should be 3 or more"))
    k < n / 2 || throw(DomainError(n, "k should be less than n/2"))
    edges = Tuple{Int,Int}[]
    for i in 0:(n-1)
        push!(edges, (i + 1, mod(i + 1, n) + 1))
        push!(edges, (i + 1, i + n + 1))
        push!(edges, (i + n + 1, mod(i + n + k, n) + n + 1))
    end
    f = mutable ? plaingraph : immutableplaingraph
    return f(2n, edges)
end

petersengraph(; kwargs...) = generalizedpetersen(5, 2; kwargs...)
dodecahedralgraph(; kwargs...) = generalizedpetersen(10, 2; kwargs...)


text:
using RegularExpressions
import Documenter: makedocs, deploydocs

makedocs(
    sitename = "RegularExpressions.jl",
    strict = true,
    modules = [RegularExpressions]
)

deploydocs(
    repo = "github.com/bramtayl/RegularExpressions.jl.git"
)


text:
using Test
import CryptoGroups: PGroup, Hash, specialize, value, ROPRG, MODP

tobig(x) = parse(BigInt, bytes2hex(reverse(x)), base=16)
interpret(::Type{BigInt}, x::Vector{UInt8}) = tobig(reverse(x))

interpret(::Type{Vector{UInt8}}, x::Integer) = reverse(reinterpret(UInt8, [x]))

function leaf(x::Vector{UInt8})

    N = UInt32(length(x))

    LEAF = UInt8(1)

    Nbin = interpret(Vector{UInt8}, N)
    bin = UInt8[LEAF, Nbin..., x...]

    return bin
end

leaf(x::String) = leaf(Vector{UInt8}(x))

### Let's make the setup complete. From repo I ahve a following public parameters:
nr = 100
rohash = Hash("sha256")
prghash = Hash("sha256")

ρ = hex2bytes("15e6c97600bbe30125cbc08598dcde01a769c15c8afe08fe5b7f5542533159e9")

# Need to unmarchal this into numbers

# group_spec = "00000000020100000020636f6d2e766572696669636174756d2e61726974686d2e4d6f645047726f757000000000040100000041009a91c3b704e382e0c772fa7cf0e5d6363edc53d156e841555702c5b6f906574204bf49a551b695bed292e0218337c0861ee649d2fe4039174514fe2c23c10f6701000000404d48e1db8271c17063b97d3e7872eb1b1f6e29e8ab7420aaab8162db7c832ba1025fa4d2a8db4adf69497010c19be0430f7324e97f201c8ba28a7f1611e087b3010000004100300763b0150525252e4989f51e33c4e6462091152ef2291e45699374a3aa8acea714ff30260338bddbb48fc7446b273aaada90e3ee8326f388b582ea8a073502010000000400000001"

p = 8095455969267383450536091939011431888343052744233774808515030596297626946131583007491317242326621464546135030140184007406503191036604066436734360427237223
q = 4047727984633691725268045969505715944171526372116887404257515298148813473065791503745658621163310732273067515070092003703251595518302033218367180213618611

h_str = "(1da949a3dfbeb316e9b225bc7d75b78d0ddd5e44fc382e74f3de95ad10eac798c4cc7be7e57d3afb259964c90fe7eb7e28a7673228d6b35a789dabd0d8351675,2937c1c4771b70c8b3dc935681aa8cef45ed24cf0c74cd2f9599ea5876850936168c5a092270d6396e634b9e46ab59836f509d0a68a65edc1426a87fd157798f,3767a5e108ecc067c31cd8bf544bcebfa3e3002f2af499ad3568e6fdb2775f4147bcf9485aeea3180b81aad365e3f4375c280941fdae10812bf1ba445030ab00,5c83d2d4f9cf3b3c2267788dd95a14740964e89f177690d06c43adfd137a101698c6bdc7ea9565fd0b18b3cfa89a2d5ea860de3912af4ffbae2c1f722912f6b9,2c453b082a201f3410d0b2906e7f2998a2c50d84f8eb0af017bf76124424bd385c15446eaad77b3ef1831798093f688331475c754f9e737d1485d80f0ef8289b,534897b93e335b0de44f3426aa6597170857f610ce0e4ff8af8efc5010ab5dc205ec53da9759e0c5ad0d2fdbf3d116e1ad7a94629c2ee331e8ceb29417983c50,7d09ab7d0704e68376e933ccfbbaa0a2e8b10a270ce72a5add597ae66e692bd930b019cc9be30b25e859dd6ba5e3bd74a872588b1e14885bacaeca3b254cc86a,33ff08f3aae9132accec73805ec2cd98ce871779369d69d135a7e1fb75b653d415480be64116feda3176dc83b2d54cf8240f8c6f1d85168081efb4f013106cdb,921ce5123ccb3b891f9e1f314d84b8ca827c44dc2027dd4fd8272189ba076eba7c2cfe7e8f8a8915b263fc72232e5d5f193f15ff6aece5b9c55f4ca9c1887fcf,410494b95d3e453e46d0b8915765f2316fcf4db7d7eff0349fe674a6498d3aaa2c3bae2b4e3c7ec1a5d5fd90d88a1e9886914f4b8d039cd470d69642b0203b73)"

𝐡 = interpret.(BigInt, hex2bytes.(split(h_str[2:end-1], ",")))

d = [ρ..., leaf("generators")...]
roprg = ROPRG(d, rohash, prghash)
prg = roprg(UInt8[]) # d is a better argument than x

sp = MODP(p; q)

𝐡′ = rand(prg, sp, 10; nr)

@test 𝐡 == 𝐡′


text:
using FundamentalsNumericalComputation
using Test
using LinearAlgebra

@testset "Chapter 1" begin
	@test FNC.horner([-1,3,-3,1],1.6) ≈ 0.6^3
end

@testset "Chapter 2" begin
	A = [ 1 2 3 0; -1 1 2 -1; 3 1 2 4; 1 1 1 1 ]
	L,U = FNC.lufact(A)
	@test norm(L*U - A) < 100eps()
	@test norm(U - triu(U)) < 100eps()
	@test norm(L - tril(L)) < 100eps()
	b = [1,10,0,-1] / 5;
	@test norm(L\b - FNC.forwardsub(L,b)) < 100eps()
	@test norm(U\b - FNC.backsub(U,b)) < 100eps()
	L,U,p = FNC.plufact(A)
	@test norm(L*U - A[p,:]) < 100eps()
	@test norm(U - triu(U)) < 100eps()
	@test norm(L - tril(L)) < 100eps()
end

@testset "Chapter 3" begin
	A = [3 4 5;-1 0 1;4 2 0; 1 1 2; 3 -4 1]
	b = 5:-1:1
	@test FNC.lsnormal(A,b) ≈ A\b
	@test FNC.lsqrfact(A,b) ≈ A\b
	Q,R = qr(A)
	QQ,RR = FNC.qrfact(A)
	@test Q ≈ QQ
	@test R ≈ RR[1:3,:]
end

@testset "Chapter 4" begin

	for c in [2,4,7.5,11]
		f = x -> exp(x) - x - c;
		dfdx = x -> exp(x) - 1;
		x = FNC.newton(f,dfdx,1.0);  r = x[end];
		@test abs(f(r)) < 100eps()
	end

	for c in [2,4,7.5,11]
		f = x -> exp(x) - x - c;
		dfdx = x -> exp(x) - 1;
		x = FNC.secant(f,3,0.5);  r = x[end];
		@test abs(f(r)) < 100eps()
	end

	function nlfun(x)
		f = zeros(3)
		f[1] = exp(x[2]-x[1]) - 2;
		f[2] = x[1]*x[2] + x[3];
		f[3] = x[2]*x[3] + x[1]^2 - x[2];
		return f
	end
	function nljac(x)
		J = zeros(3,3)
		J[1,:] = [-exp(x[2]-x[1]),exp(x[2]-x[1]), 0]
		J[2,:] = [x[2], x[1], 1]
		J[3,:] = [2*x[1], x[3]-1, x[2]]
		return J
	end

	x = FNC.newtonsys(nlfun,nljac,[0,0,0]);
	@test norm(nlfun(x[end])) < 100eps()
	x = FNC.newtonsys(nlfun,nljac,[1,2,3]);
	@test norm(nlfun(x[end])) < 100eps()

	x = FNC.levenberg(nlfun,[10,-4,-3])
	@test norm(nlfun(x[end])) < 1e-12

end

@testset "Chapter 5" begin
	f = t->cos(5t)
	Q,t = FNC.intadapt(f,-1,3,1e-8)
	@test Q ≈ (sin(15)+sin(5))/5 rtol = 1e-5
	T,_ = FNC.trapezoid(f,-1,3,820)
	@test T ≈ (sin(15)+sin(5))/5 rtol = 1e-4

	t = [-2,-0.5,0,1,1.5,3.5,4]/10
	w = FNC.fdweights(t.-0.12,2)
	f = x->cos(3x)
	@test dot(w,f.(t)) ≈ -9cos(0.36) rtol = 1e-3

	t = [-2,-0.5,0,1,1.5,3.5,4]/10
	H = FNC.hatfun(t,5)
	@test H(0.22) ≈ (0.22-t[5])/(t[6]-t[5])
	@test H(.38) ≈ (t[7]-.38)/(t[7]-t[6])
	@test H(0.06)==0
	@test H(t[6])==1
	@test H(t[7])==0
	@test H(t[1])==0
	p = FNC.plinterp(t,f.(t))
	@test p(0.22) ≈ f(t[5]) + (f(t[6])-f(t[5]))*(0.22-t[5])/(t[6]-t[5])
	S = FNC.spinterp(t,exp.(t))
	x = [-.17,-0.01,0.33,.38]
	@test S.(x) ≈ exp.(x) rtol = 1e-5
	@test S.(t) ≈ exp.(t) rtol = 1e-11
end

@testset "Chapter 6" begin
	f = (u,p,t) -> u + p*t^2
	û = exp(1.5) - 2*(-2 + 2*exp(1.5) - 2*1.5 - 1.5^2)
	ivp = ODEProblem(f,1,(0,1.5),-2)
	t,u = FNC.euler(ivp,4000)
	@test û ≈ u[end] rtol = 0.005
	t,u = FNC.am2(ivp,4000)
	@test û ≈ u[end] rtol = 0.005

	g = (u,p,t) -> [t+p-sin(u[2]),u[1]]
	ivp = ODEProblem(g,[-1.,4],(1.,2.),-6)
	sol = solve(ivp,Tsit5())
	t,u = FNC.euler(ivp,4000)
	@test u[end] ≈ sol.u[end] rtol=0.004
	t,u = FNC.ie2(ivp,4000)
	@test u[end] ≈ sol.u[end] rtol=0.0005
	t,u = FNC.rk4(ivp,800)
	@test u[end] ≈ sol.u[end] rtol=0.0005
	t,u = FNC.ab4(ivp,800)
	@test u[end] ≈ sol.u[end] rtol=0.0005
	t,u = FNC.rk23(ivp,1e-4)
	@test u[end] ≈ sol.u[end] rtol=0.0005
	t,u = FNC.am2(ivp,2000)
	@test u[end] ≈ sol.u[end] rtol=0.0005
end

@testset "Chapter 8" begin
	V = randn(4,4)
	D = diagm([-2,0.4,-0.1,0.01])
	A = V*D/V;

	γ,x = FNC.poweriter(A,30)
	@test γ[end] ≈ -2 rtol=1e-10
	@test abs( dot(x,V[:,1])/(norm(V[:,1])*norm(x)) ) ≈ 1 rtol=1e-10

	γ,x = FNC.inviter(A,0.37,15)
	@test γ[end] ≈ 0.4 rtol=1e-10
	@test abs( dot(x,V[:,2])/(norm(V[:,2])*norm(x)) ) ≈ 1 rtol=1e-10

	Q,H = FNC.arnoldi(A,ones(4),4)
	@test A*Q[:,1:4] ≈ Q*H

	x,res = FNC.gmres(A,ones(4),3)
	@test norm(ones(4) - A*x) ≈ res[end]
	x,res = FNC.gmres(A,ones(4),4)
	@test A*x ≈ ones(4)
end

@testset "Chapter 9" begin
	f = x -> exp(sin(x)+x^2)
	t = [-cos(k*π/40) for k in 0:40 ]
	p = FNC.polyinterp(t,f.(t))
	@test p(-0.12345) ≈ f(-0.12345)

	f = x -> exp(sin(π*x))
	n = 30
	t = [ 2k/(2n+1) for k in -n:n ]
	p = FNC.triginterp(t,f.(t))
	@test p(-0.12345) ≈ f(-0.12345)
	t = [ k/n for k in -n:n-1 ]
	p = FNC.triginterp(t,f.(t))
	@test p(-0.12345) ≈ f(-0.12345)

	F = x -> tan(x/2-0.2)
	f = x -> 0.5*sec(x/2-0.2)^2
	@test FNC.ccint(f,40)[1] ≈ F(1)-F(-1)
	@test FNC.glint(f,40)[1] ≈ F(1)-F(-1)

	f = x -> 1/(32+2x^4)
	@test FNC.intinf(f,1e-9)[1] ≈ sqrt(2)*π/32 rtol=1e-5

	f = x -> (1-x)/( sin(x)^0.5 )
	@test FNC.intsing(f,1e-8)[1] ≈ 1.34312 rtol=1e-5
end

@testset "Chapter 10" begin
	λ = 0.6
	phi = (r,w,dwdr) -> λ/w^2 - dwdr/r;
	a = eps();  b = 1;
	g₁ = (u,du)->du    # Neumann at left
	g₂ = (u,du)->u-1   # Dirichlet at right

	r,w,dwdx = FNC.shoot(phi,(a,b),g₁,g₂,[0.8,0])
	@test w[1] ≈ 0.78776 rtol=1e-4

	f = x -> exp(x^2-3x)
	df = x -> (2x-3)*f(x)
	ddf = x -> ((2x-3)^2+2)*f(x)

	t,D,DD = FNC.diffmat2(400,(-0.5,2))
	@test df.(t) ≈ D*f.(t) rtol=1e-3
	@test ddf.(t) ≈ DD*f.(t) rtol=1e-3
	t,D,DD = FNC.diffcheb(80,(-0.5,2))
	@test df.(t) ≈ D*f.(t) rtol=1e-7
	@test ddf.(t) ≈ DD*f.(t) rtol=1e-7

	exact = x -> exp(sin(x));
	p = x -> -cos(x);
	q = sin;
	r = x -> 0;
	x,u = FNC.bvplin(p,q,r,[0,pi/2],1,exp(1),300);
	@test u ≈ exact.(x) rtol=1e-3

	ϕ = (t,θ,ω) -> -0.05*ω - sin(θ);
	g1(u,du) = u-2.5
	g2(u,du) = u+2
	init = collect(range(2.5,-2,length=101));

	t,θ = FNC.bvp(ϕ,[0,5],g1,g2,init)
	@test θ[7] ≈ 2.421850016880724 rtol=1e-10

	c = x -> x^2;
	 q = x -> 4;
	f = x -> sin(π*x);
	x,u = FNC.fem(c,q,f,0,1,100)
	@test u[33] ≈ 0.1641366907307196 rtol=1e-10
end

@testset "Chapter 11" begin
	s = x -> sin(π*(x-0.2))
	c = x -> cos(π*(x-0.2))
	f = x -> 1 + s(x)^2
	df = x -> 2π*s(x)*c(x)
	ddf = x -> 2π^2*(c(x)^2 - s(x)^2)

	t,D,DD = FNC.diffper(400,(0,2))
	@test df.(t) ≈ D*f.(t) rtol=1e-3
	@test ddf.(t) ≈ DD*f.(t) rtol=1e-3

	phi = (t,x,u,uₓ,uₓₓ) -> uₓₓ + t*u
	g1(u,ux) = ux;
	g2(u,ux) = u-1;
	init(x) = x^2;
	x,u = FNC.parabolic(phi,(0,1),40,g1,g2,(0,2),init);
	@test u(0.5)[21] ≈ 0.845404 rtol=1e-4
	@test u(1)[end] ≈ 1 rtol=1e-4
	@test u(2)[1] ≈ 2.45692 rtol=1e-4
end

@testset "Chapter 13" begin

	f = (x,y) -> -sin(3*x.*y-4*y)*(9*y^2+(3*x-4)^2);
	g = (x,y) -> sin(3*x*y-4*y);
	xspan = [0,1];  yspan = [0,2];
	x,y,U = FNC.poissonfd(f,g,60,xspan,60,yspan);
	X = [x for x in x,y in y];
	Y = [y for x in x,y in y];
	@test g.(X,Y) ≈ U rtol=1e-3

	λ = 1.5
	function pde(X,Y,U,Ux,Uxx,Uy,Uyy)
		return @. Uxx + Uyy - λ/(U+1)^2   # residual
	end
	g = (x,y) -> x+y     # boundary condition
	u = FNC.elliptic(pde,g,30,[0,2.5],24,[0,1]);
	@test u(1.25,0.5) ≈ 1.7236921361 rtol = 1e-6
	@test u(1,0) ≈ 1 rtol = 1e-6
end


text:
push!(LOAD_PATH, joinpath(@__DIR__, "..", ".."))

using BenchmarkTools, Plots

proj = Base.active_project()
julia_cmd = `$(Base.julia_cmd()) --startup-file=no --project=$(proj)`
run_script = joinpath(@__DIR__, "run_ref_thread.jl")

thread_counts = [1, 2, 4, 8, 16]
runtimes = Float64[]
result = IOBuffer()
for n in thread_counts
    cmd = deepcopy(julia_cmd)
    append!(cmd.exec, ["--threads=$(n)", run_script])
    Base.run(pipeline(cmd, stdout = result); wait = true)
    push!(runtimes, parse(Float64, String(take!(result))))
end

using Plots
ENV["GKSwstype"] = "nul"

plt = plot(
    xlabel = "#threads",
    ylabel = "time (secs)",
    title = "Volume thread scaling",
)
plot!(plt, thread_counts, runtimes, label = "Runtimes w/threads")
png(plt, joinpath(@__DIR__, "vol_thread_scaling.png"))

function linkfig(figpath, alt = "")
    # buildkite-agent upload figpath
    # link figure in logs if we are running on CI
    if get(ENV, "BUILDKITE", "") == "true"
        artifact_url = "artifact://$figpath"
        print("\033]1338;url='$(artifact_url)';alt='$(alt)'\a\n")
    end
end

linkfig("benchmarks/bickleyjet/vol_thread_scaling.png", "Volume thread scaling")


text:
export @nview, apy2, nrm2, swap!, scal!, axpy!, xpby!, axpby!

function mach_acc() # @code_warntype ✓
    sfmin = nextfloat(0.0)
    small = inv(prevfloat(Inf))
    if small ≥ sfmin
        # Use "small" plus 1 bit, to avoid the possibility of rounding
        # causing overflow when computing 1/sfmin.
        sfmin = small * (1.0 + eps())
    end
    return sfmin
end

const SAFMIN = mach_acc() / eps()

macro nview(arr::Symbol, N::Int)
    a = Vector{Expr}(undef, N)
    @inbounds for i in 1:N
        a[i] = Expr(:(=), Symbol("k", i), :(view($arr, :, $i)))
    end
    return Expr(:escape, Expr(:block, a...))
end

macro nview(arr::Symbol, vars::Symbol...)
    N = length(vars)
    a = Vector{Expr}(undef, N)
    @inbounds for i in 1:N
        a[i] = Expr(:(=), vars[i], :(view($arr, :, $i)))
    end
    return Expr(:escape, Expr(:block, a...))
end

function apy2(x::Real, y::Real) # @code_warntype ✓
    isnan(x) && return x
    isnan(y) && return y
    # general case
    xabs = abs(x)
    yabs = abs(y)
    w = max(xabs, yabs)
    z = min(xabs, yabs)
    iszero(z) && return w
    return w * sqrt(1.0 + abs2(z / w))
end

function swap!(v::VecI, i::Int, j::Int) # @code_warntype ✓
    @inbounds temp = v[i]
    @inbounds v[i] = v[j]
    @inbounds v[j] = temp
    return nothing
end

function swap!(v::MatI, i1::Int, j1::Int, i2::Int, j2::Int) # @code_warntype ✓
    @inbounds temp     = v[i1,j1]
    @inbounds v[i1,j1] = v[i2,j2]
    @inbounds v[i2,j2] = temp
    return nothing
end

function nrm2(x::VecI{Tx}, y::VecI{Ty}, b::VecB{Tb}) where {Tx<:Real,Ty<:Real,Tb<:Real} # @code_warntype ✓
    @simd for i in eachindex(b)
        @inbounds b[i] = abs2(x[i] - y[i])
    end
    return sqrt(sum(b))
end

function scal!(a::Real, x::VecO{Tx}) where Tx<:Real # @code_warntype ✓
    isone(a) && return nothing
    if iszero(a)
        zeroTx = zero(Tx)
        @simd for i in eachindex(x)
            @inbounds x[i] = zeroTx
        end
        return nothing
    end
    @simd for i in eachindex(x)
        @inbounds x[i] *= a
    end
end

# a * x + y → y, @code_warntype ✓
function axpy!(a::Real, x::VecI{Tx}, y::VecO{Ty}) where {Tx<:Real,Ty<:Real}
    iszero(a) && return nothing
    if isone(a)
        @simd for i in eachindex(y)
            @inbounds y[i] += x[i]
        end
        return nothing
    end
    @simd for i in eachindex(y)
        @inbounds y[i] += a * x[i]
    end
end

# x + b * y → y, @code_warntype ✓
function xpby!(x::VecI{Tx}, b::Real, y::VecO{Ty}) where {Tx<:Real,Ty<:Real}
    if iszero(b)
        @simd for i in eachindex(y)
            @inbounds y[i] = x[i]
        end
        return nothing
    end
    if isone(b)
        @simd for i in eachindex(y)
            @inbounds y[i] += x[i]
        end
        return nothing
    end
    @simd for i in eachindex(y)
        @inbounds y[i] = x[i] + b * y[i]
    end
end

# a * x + b * y → y, @code_warntype ✓
function axpby!(a::Real, x::VecI{Tx}, b::Real, y::VecO{Ty}) where {Tx<:Real,Ty<:Real}
    isone(b) && return axpy!(a, x, y)
    if iszero(b)
        iszero(a) && return scal!(0, y)
        if isone(a)
            @simd for i in eachindex(y)
                @inbounds y[i] = x[i]
            end
        else
            @simd for i in eachindex(y)
                @inbounds y[i] = a * x[i]
            end
        end
        return nothing
    end
    iszero(a) && return scal!(b, y)
    if isone(a)
        @simd for i in eachindex(y)
            @inbounds y[i] = x[i] + b * y[i]
        end
    else
        @simd for i in eachindex(y)
            @inbounds y[i] = a * x[i] + b * y[i]
        end
    end
end


text:
## SingFun stieltjes

export logabs
logabs(z) = log(abs2(z))/2
sqrtabs(z) = sqrt(abs(z))

# sqrtx2 is analytic continuation of sqrt(z^2-1)
# with the oriented branch cut [-1,1]
# we need to reverseorientation, as sqrt has branch cut [0,∞)
sqrtx2(z::Directed) = sqrt(reverseorientation(z)-1)*sqrt(z.x+1)
sqrtx2(z::Number) = sqrt(z-1)*sqrt(z+1)
sqrtx2(x::Real) = sign(x)*sqrt(x^2-1)

# these mimc logabs.  They are useful because they are continuous functions
# and are used in logkernel for evaluating on the branch cut
# sqrtx2im is omitted as it's discontinuous on the branch cut
sqrtx2abs(z) = sqrtabs(z-1)*sqrtabs(z+1)
sqrtx2real(z) = sqrtx2abs(z)*cos((angle(z-1)+angle(z+1))/2)
#  real(z̄*sqrt(z-1)*sqrt(z+1))
# sqrtx2abs(z)*abs(z)*cos(angle(z-1)+angle(z+1))/2-angle(z))

x̄sqrtx2real(z) = sqrtx2abs(z)*abs(z)*cos((angle(z-1)+angle(z+1))/2-angle(z))

function sqrtx2(f::Fun)
    B = Evaluation(leftendpoint(domain(f)))
    A = Derivative()-f*differentiate(f)/(f^2-1)
    \([B,A],[sqrtx2(first(f)),0];tolerance=ncoefficients(f)*10E-15)
end

# these are two inverses of the joukowsky map
# the first maps the slit plane to the inner circle, the second to the outer circle
#
# it is more accurate near infinity to do 1/J_- than z - sqrtx2(z) as it avoids round off
joukowskyinverse(::Type{Val{true}},z) = 1/joukowskyinverse(Val{false},z)
joukowskyinverse(::Type{Val{false}},z) = undirected(z)+sqrtx2(z)

joukowskyinverseabs(::Type{Val{true}},z) = 1/joukowskyinverseabs(Val{false},z)
joukowskyinverseabs(::Type{Val{false}},z) = sqrt(abs2(z)+2x̄sqrtx2real(z)+sqrtx2abs(z)^2)

joukowskyinversereal(::Type{Val{true}},z) =
    joukowskyinversereal(Val{false},z)/joukowskyinverseabs(Val{false},z)^2
joukowskyinversereal(::Type{Val{false}},z) = real(z)+sqrtx2real(z)

function hornersum(cfs::AbstractVector{S},y::V) where {S<:Number,V<:Number}
    N,P = length(cfs),Base.promote_op(*,S,V)
    ret = N > 0 ? convert(P,cfs[N]) : zero(P)
    for k=N-1:-1:1
        ret = muladd(y,ret,cfs[k])
    end
    y*ret
end

function divkhornersum(cfs::AbstractVector{S},ζ::T,ζs::U,s::V) where {S<:Number,T<:Number,U<:Number,V<:Number}
    N,P = length(cfs),promote_type(S,T,U,V)
    ret = N > 0 ? convert(P,cfs[N]/(N+s)) : zero(P)
    for k=N-1:-1:1
        ret=muladd(ζ,ret,cfs[k]/(k+s))
    end
    ζ*ζs*ret
end

realdivkhornersum(cfs::AbstractVector{S},ζ,ζs,s) where {S<:Real} = real(divkhornersum(cfs,ζ,ζs,s))
realdivkhornersum(cfs::AbstractVector{S},ζ,ζs,s) where {S<:Complex} =
    complex(real(divkhornersum(real(cfs),ζ,ζs,s)),
            real(divkhornersum(imag(cfs),ζ,ζs,s)))

function stieltjes(sp::JacobiWeight{<:PolynomialSpace,<:IntervalOrSegment},u,z)
    d=domain(sp)

    if sp.α == sp.β == .5
        cfs = coefficients(u,sp.space,Ultraspherical(1,d))
        π*hornersum(cfs,joukowskyinverse(Val{true},mobius(sp,z)))
    elseif sp.α == sp.β == -.5
        cfs = coefficients(u,sp.space,ChebyshevDirichlet{1,1}(d))
        z=mobius(sp,z)

        sx2z=sqrtx2(z)
        sx2zi=1/sx2z
        Jm=1/(undirected(z)+sx2z)  # joukowskyinverse(true,z)

        if length(cfs) ≥1
            ret = π*cfs[1]*sx2zi

            if length(cfs) ≥2
                ret += cfs[2]*π*(undirected(z).*sx2zi-1)
            end

            ret - 2π*hornersum(cfs[3:end],Jm)
        else
            zero(undirected(z))
        end
    elseif isapproxinteger(sp.α) && isapproxinteger(sp.β)
        stieltjes(sp.space,coefficients(u,sp,sp.space),z)
    else
        if d==ChebyshevInterval()
            S2=JacobiWeight(sp.β,sp.α,Jacobi(sp.β,sp.α))  # convert and then use recurrence
            stieltjesintervalrecurrence(S2,coefficients(u,sp,S2),z)
        else
            # project to interval
            stieltjes(setdomain(sp,ChebyshevInterval()),u,mobius(sp,z))
        end
    end
end

##
#  hilbert on JacobiWeight space
#  hilbert is always equal to im*(C^+ + C^-)
#  hilbert(f,z)=im*(cauchy(true,f,z)+cauchy(false,f,z))
##

function hilbert(sp::JacobiWeight{Chebyshev{DD},DD},u) where {DD<:IntervalOrSegment}
    d=domain(u)

    if sp.α == sp.β == .5
        # Corollary 5.7 of Olver&Trogdon
        cfs=coefficients(u,sp.space,Ultraspherical(1))
        Fun(d,[0.;-cfs])
    elseif sp.α == sp.β == -.5
        # Corollary 5.11 of Olver&Trogdon
        cfs= coefficients(u,sp.space,ChebyshevDirichlet{1,1})
        if length(cfs)≥2
            Fun(d,[cfs[2];2cfs[3:end]])
        else
            Fun(d,zeros(eltype(cfs),1))
        end
    else
        error("hilbert only implemented for Chebyshev weights")
    end
end

integratejin(c,cfs,ζ) =
    0.5*(-(cfs[1] == 0 ? cfs[1] : cfs[1]*(log(ζ)+log(c))) +
            divkhornersum(cfs,ζ,ζ,1) -
            divkhornersum(view(cfs,2:length(cfs)),ζ,zero(ζ)+1,0))
realintegratejin(c,cfs,ζ) =
    0.5*(-(cfs[1] == 0 ? cfs[1] : cfs[1]*(logabs(ζ)+logabs(c))) +
        realdivkhornersum(cfs,ζ,ζ,1) -
        realdivkhornersum(view(cfs,2:length(cfs)),ζ,zero(ζ)+1,0))

#########
# stieltjesintegral is an indefinite integral of stieltjes
# normalized so that there is no constant term
# logkernel is the real part of stieljes normalized by π.
#####

function logkernel(sp::JacobiWeight{<:PolynomialSpace,<:IntervalOrSegment},u,z)
    d=domain(sp)
    a,b= endpoints(d)

    if sp.α == sp.β == .5
        cfs=coefficients(u,sp.space,Ultraspherical(1,d))
        z=mobius(sp,z)
        x,r = joukowskyinversereal(Val{true},z),joukowskyinverseabs(Val{true},z)
        y =  sqrt(abs(r^2 - x^2))  # dummy variable, choice of ± in arg doesn't matter
        ζ = x + im*y
        arclength(d)*realintegratejin(4/(b-a),cfs,ζ)/2
    elseif  sp.α == sp.β == -.5
        cfs = coefficients(u,sp.space,ChebyshevDirichlet{1,1}(d))
        z=mobius(sp,z)

        x,r = joukowskyinversereal(Val{true},z),joukowskyinverseabs(Val{true},z)
        y = sqrt(abs(r^2 - x^2) )
        ζ = x + im*y

        if length(cfs) ≥1
            ret = -cfs[1]*arclength(d)*(log(r)+logabs(4/(b-a)))/2

            if length(cfs) ≥2
                ret += -arclength(d)*cfs[2]*x/2
            end

            if length(cfs) ≥3
                ret - arclength(d)*realintegratejin(4/abs(b-a),view(cfs,3:length(cfs)),ζ)
            else
                ret
            end
        else
            zero(z)
        end
    elseif domain(sp) == ChebyshevInterval()
        DS=WeightedJacobi(sp.β+1,sp.α+1)
        D=Derivative(DS)[2:end,:]

        f=Fun(Fun(sp,u),WeightedJacobi(sp.β,sp.α))  # convert to Legendre expansion
        uu=D\(f|(2:∞))   # find integral, dropping first coefficient of f

        (f.coefficients[1]*logabsjacobimoment(sp.α,sp.β,z) + real(stieltjes(uu,z)))/π
    else # map to canonical
        c = 1/abs(fromcanonicalD(d))
        g = setcanonicaldomain(Fun(sp,u))
        (logkernel(g,mobius(d,z))/c - logabs(c)sum(g)/(c*π))
    end
end

function stieltjesintegral(sp::JacobiWeight{<:PolynomialSpace,<:IntervalOrSegment},u,z)
    d=domain(sp)
    a,b=endpoints(d)

    if sp.α == sp.β == .5
        cfs=coefficients(u,sp.space,Ultraspherical(1,d))
        y=joukowskyinverse(Val{true},mobius(sp,z))
        π*complexlength(d)*integratejin(4/(b-a),cfs,y)/2
    elseif  sp.α == sp.β == -0.5
        cfs = coefficients(u,sp.space,ChebyshevDirichlet{1,1}(d))
        z=mobius(sp,z)
        ζ=joukowskyinverse(Val{true},z)

        if length(cfs) ≥1
            ret = -cfs[1]*π*complexlength(d)*(log(ζ)+log(4/abs(b-a)))/2

            if length(cfs) ≥2
                ret += -π*complexlength(d)*cfs[2]*ζ/2
            end

            if length(cfs) ≥3
                ret - π*complexlength(d)*integratejin(4/abs(b-a),view(cfs,3:length(cfs)),ζ)
            else
                ret
            end
        else
            zero(z)
        end
    else
        error("stieltjes integral not implemented for parameters β="*string(sp.β)*", α="*string(sp.α))
    end
end


text:
srand(0)

functionalitytests = [
    "sintest1"
  ]

println("============================")
println("Running functionality tests:")
println("============================")

for t in functionalitytests
  tfile = t*".jl"
  println("  * $(tfile) *")
  include(string("functionality/",tfile))
end


text:
using Oxygen

println("Testing Operations")

@testset "MolOps" begin
   ethyl_acetate = smilestomol("CC(=O)OCC")
   @test Oxygen.getbonds(ethyl_acetate) == [
      (1, 2, 1.0)
      (2, 3, 2.0)
      (2, 4, 1.0)
      (4, 5, 1.0)
      (5, 6, 1.0)
   ]
   @test Oxygen.getadjmatrix(ethyl_acetate) == [
      0.0 1.0 0.0 0.0 0.0 0.0
      1.0 0.0 2.0 1.0 0.0 0.0
      0.0 2.0 0.0 0.0 0.0 0.0
      0.0 1.0 0.0 0.0 1.0 0.0
      0.0 0.0 0.0 1.0 0.0 1.0
      0.0 0.0 0.0 0.0 1.0 0.0
   ]
   @test Oxygen.valencies(ethyl_acetate) == [
      1  0  0  0  0  0
      0  0  1  0  0  0
      1  0  0  0  0  0
      0  1  0  0  0  0
      0  1  0  0  0  0
      1  0  0  0  0  0
   ]
   @test Oxygen.atomtypes(ethyl_acetate) == [
      1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
      1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
      0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
      0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
      1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
      1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
   ]
   @test molfeatures(ethyl_acetate) == [
      1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
      1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0
      0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0
      0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0
      1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0
      1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
   ]

   benzene = smilestomol("c1ccccc1")
   @test Oxygen.getbonds(benzene) == [
      (1, 2, 1.5)
      (1, 6, 1.5)
      (2, 3, 1.5)
      (3, 4, 1.5)
      (4, 5, 1.5)
      (5, 6, 1.5)
   ]
   @test Oxygen.getadjmatrix(benzene) == [
      0.0  1.5  0.0  0.0  0.0  1.5
      1.5  0.0  1.5  0.0  0.0  0.0
      0.0  1.5  0.0  1.5  0.0  0.0
      0.0  0.0  1.5  0.0  1.5  0.0
      0.0  0.0  0.0  1.5  0.0  1.5
      1.5  0.0  0.0  0.0  1.5  0.0
   ]
   @test Oxygen.valencies(benzene) == [
      0  1  0  0  0  0
      0  1  0  0  0  0
      0  1  0  0  0  0
      0  1  0  0  0  0
      0  1  0  0  0  0
      0  1  0  0  0  0
   ]
   @test Oxygen.atomtypes(benzene) == [
      1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
      1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
      1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
      1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
      1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
      1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
   ]
   @test molfeatures(benzene) == [
      1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0
      1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0
      1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0
      1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0
      1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0
      1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0
   ]
end

@testset "AtomOps" begin
   ethyl_acetate = smilestomol("CC(=O)OCC")
   @test map(Oxygen.valence_electrons, ethyl_acetate.atoms) == [4, 4, 6, 6, 4, 4]
   @test Oxygen.neighbors(ethyl_acetate, 2) == [OxygenAtom("C"), OxygenAtom("O"), OxygenAtom("O")]

   benzene = smilestomol("c1ccccc1")
   @test map(Oxygen.valence_electrons, benzene.atoms) == [4, 4, 4, 4, 4, 4]
   @test Oxygen.neighbors(benzene, 2) == [OxygenAtom("C"), OxygenAtom("C")]
end

@testset "GraphUtils" begin
   benzene = smilestomol("c1ccccc1")
   @test Oxygen.dfs(benzene) == [(1, 12), (2, 11), (3, 10), (4, 9), (5, 8), (6, 7)]

   ethyl_acetate = smilestomol("CC(=O)OCC")
   @test Oxygen.dfs(ethyl_acetate) == [(1, 12), (2, 11), (3, 4), (5, 10), (6, 9), (7, 8)]
end


text:
using Documenter, UiPathOrchestratorJobSchedulingPlanCreate1

makedocs(;
    modules=[UiPathOrchestratorJobSchedulingPlanCreate1],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "Getting Started" => "getting_started.md",
        "Function" => "function.md",
        "Excel" => "excel.md",
    ],
    repo="https://github.com/wakakusa/UiPathOrchestratorJobSchedulingPlanCreate1.jl/blob/{commit}{path}#L{line}",
    sitename="UiPathOrchestratorJobSchedulingPlanCreate1.jl",
    authors="wakakusa",
    assets=String[],
)

deploydocs(;
    repo="github.com/wakakusa/UiPathOrchestratorJobSchedulingPlanCreate1.jl",
)


text:
using ADCME
using PyPlot
using Random
using PyCall
np = pyimport("numpy")

# ! need kx, ky and nn_type, and ky_nn
# ! E is always 1
# ! loss is based only on y[:,1]
function hidden_function(model_type, x, x_, y_)
    if model_type == "Plasticity"
        # x  = [ε_n+1]
        # x_ = [ε_n]
        # y_ = [σ_n, α_n] or [σ_n, α_n, q_n]
        # The yield condition is
        #      f = |sigma - q| - (σY + alpha * K)
        # here K is the plastic modulus, , σY is the flow stress
        # D_eps_p = gamma df/dsigma
        # D_alpha = |D_eps_p|
        # D_q     = B * D_eps_p
        # f  = 0    or f  < 0
        #trial stress

        σY,  E,  K, B =  0.006, 1.0, 0.1, 0.01
        if length(y_) == 3
            σ0, α0, q0 = y_[1], y_[2], y_[3]
        elseif length(y_) == 2
            σ0, α0, q0 = y_[1], y_[2], 0.0
            B = 0.0
        end
        ε, ε0 = x[1], x_[1]
        Δγ = 0.0
        σ = σ0 + E*(ε - ε0)
        α = α0 + abs(Δγ)
        q = q0 + B*Δγ
        ξ = σ - q

        r2 = abs(ξ) - (σY + K*α)
        if r2 <= 0
            σ = σ0 + E*(ε - ε0)
            dΔσdΔε = E

        else

            Δγ = r2/(B + E + K)
            q += B * Δγ * sign(ξ)
            α += Δγ
            σ -= Δγ * E * sign(ξ)
            dΔσdΔε = E*(B + K)/(B + E + K)
        end
        if length(y_) == 3
            y = [σ; α + Δγ; q]
        elseif length(y_) == 2
            y = [σ; α + Δγ]
        end

    elseif model_type == "PlasticityLawBased"
        # x  = [ε_n+1]
        # x_ = [ε_n]
        # y_ = [σ_n]
        # The yield condition is
        #      f = σ_vm - (α + β * ε_vm^pow)
        # D_eps_p = gamma df/dsigma
        # f  = 0    or f  < 0

        E, α, β, pow =  1.0, 0.006, 0.018, 0.4
        σ0 = y_[1]
        ε, ε0 = x[1], x_[1]
        Δγ = 0.0

        #trial stress
        σ = σ0 + E*(ε - ε0)
        ε_vm = α + β*2.0/3.0 * abs(ε)
        r2 = abs(σ) - ε_vm
        if r2 <= 0
            σ = σ0 + E*(ε - ε0)
            dΔσdΔε = E

        else

            Δγ = r2/E
            σ -= Δγ * E * sign(σ)
        end

        y = [σ]

    else
        error("model_type ", model_type, " have not implemented yet ")
    end

    return y
end

function generate_sequence(model_type, xs, y0)
    n = size(xs,1)
    ys = zeros(n,ky)
    ys[1,:] = y0
    for i = 2:n
        ys[i,:] = hidden_function(model_type, xs[i,:], xs[i-1,:], ys[i-1,:])
        @show "generate data ", i, ys[i,:]
    end
    ys
end

@doc """
    generate data
    m set of data, each is a serial sequency of n time
"""->
function generate_data(model_type, m = 2, n = 100)
    xs_set, ys_set = [], []

    # generate xs_set
    if model_type == "Plasticity" || model_type == "PlasticityLawBased"
        T = 0.1
        #
        t = np.linspace(0.0, T, n)
        A = 0.02
        xs = A * reshape(sin.(π*t/(3.0*T)), :, kx)
        push!(xs_set, xs)

        #
        if m >= 2
            t = np.linspace(0.0, T, n)
            A = 0.02
            xs = A * reshape(sin.(π*t/(T)), :, kx)
            push!(xs_set, xs)
        end

        if m >= 3
            t = np.linspace(0.0, T, n)
            A = 0.01
            xs = A * reshape(sin.(π*t/(1.2*T)), :, kx)
            push!(xs_set, xs)
        end

        if m >= 4
            t = np.linspace(0.0, T, n)
            A = 0.01
            xs = A * reshape(sin.(π*t/(1.0*T)), :, kx)
            push!(xs_set, xs)
        end

    else
        error("model_type ", model_type, " have not implemented yet ")
    end

    # generate ys_set array
    for i = 1:m
        y0 = zeros(ky)
        ys = generate_sequence(model_type, xs_set[i], y0)
        push!(ys_set, ys)
    end

    return xs_set, ys_set
end

function sequence_test(xs_set, sess)
    m = length(xs_set)
    ys_pred_set = []
    for i_set = 1:m
        xs = xs_set[i_set]
        n = size(xs,1)
        ys_pred = zeros((n, ky_nn))

        plx = placeholder(zeros(kx))
        plx_ = placeholder(zeros(kx))
        ply_ = placeholder(zeros(ky_nn))
        res = nn(plx, plx_, ply_)

        for i = 2:n
            ys_pred[i,:] = run(sess,res, feed_dict = Dict(
            plx=>xs[i,:],
            plx_=>xs[i-1,:],
            ply_=>ys_pred[i-1,:]))
            @show i, ys_pred[i,:]
        end
        push!(ys_pred_set, ys_pred)
    end

    return ys_pred_set

end

function point2point_test(xs_set, ys_set, sess)
    m = length(xs_set)
    ys_pred_set = []
    for i_set = 1:m
        xs, ys = xs_set[i_set], ys_set[i_set]
        n = size(xs,1)
        ys_pred = zeros((n, ky))

        ys_pred[2:end,1:ky_nn] = run(sess, nn_all(xs[2:end,:], xs[1:end-1,:], ys[1:end-1, 1:ky_nn]))
        push!(ys_pred_set, ys_pred)
    end

    return ys_pred_set

end

######################################################################
function compute_sequence_loss(xs_set, ys_set, nn)
    loss = constant(0.0)
    m = length(xs_set)

    for i_set = 1:m
        xs, ys = xs_set[i_set], ys_set[i_set][:, 1:ky_nn]
        n = size(xs,1)
        y = constant(ys[1,:])
        for i = 2:n
            y = nn(constant(xs[i,:]), constant(xs[i-1,:]), y)
            #todo change the loss function
            #loss += (ys[i,:] - y)^2[1]

            @show i, size(y)

            loss += 100000.0* ((ys[i,:] - y)^2)[1]

            #loss += (ys[i,1] - y[1])^2

        end
    end
    return loss
end

######################################################################
function compute_point2point_loss(xs_set, ys_set, nn)
    loss = constant(0.0)
    m = length(xs_set)

    for i_set = 1:m
        xs, ys = xs_set[i_set], ys_set[i_set]

        y = nn_all(xs[2:end,:], xs[1:end-1,:], ys[1:end-1, 1:ky_nn])

    loss += sum((y - ys[2:end, 1:ky_nn])^2)
    end
    return loss
end

config=[20,ky_nn]

function nn(x, x_, y_)

    #ipt = reshape([x;x_;y_], 1, :)
    if nn_type == "ae"
        ipt = reshape([x;x_;y_], :, 2*kx + ky_nn)

        out = ae(ipt, config)
        squeeze(out, dims=1)

    elseif nn_type=="piecewise"

        E = 1.0
        threshold = 1.e-5
        ipt = reshape([x;x_;y_], :, 2*kx + ky_nn)

        σnn = ae(ipt, config)
        σH = E * (x - x_)
        z = sum(x^2)

        i = sigmoid(1e6*(z-threshold))

        out = σnn * i + σH * (1-i)  + y_

        squeeze(out, dims=1)

    elseif nn_type=="piecewise2"

        @assert(ky_nn == 1)
        E = 1.0
        threshold = 1.e-5
        ipt = reshape([x;x_;y_], :, 2*kx + ky_nn)

        σnn = reshape(ae(ipt, config)* (x - x_), 1, 1)

        σH = E * (x - x_)
        z = sum(y_^2)

        i = sigmoid(1e6*(z-threshold))

        out = σnn * i + σH * (1-i)  + y_

        squeeze(out, dims=1)
    else
        error("$nn_type does not exist")
    end
end

function sigmoid_(z)

    return 1.0 ./ (1.0 .+ exp.(-z))

end

function nn_all(x, x_, y_)
    if nn_type == "ae"

        ipt = [x x_ y_]

        #out = ae(ipt, [20,ky_nn])
        out = ae(ipt, config)
    elseif nn_type=="piecewise"

        E = 1.0
        threshold = 1.e-5
        ipt = [x x_ y_]

        σnn = ae(ipt, config)
        σH = E * (x - x_)
        z = x.^2

        i = sigmoid_(1e6*(z.-threshold))

        out = σnn .* i + σH .* (1 .- i)  .+ y_

    elseif nn_type=="piecewise2"

        E = 1.0
        threshold = 1.e-5
        ipt = [x x_ y_]

        σnn = ae(ipt, config) .* (x - x_)

        σH = E * (x - x_)
        z = y_.^2

        i = sigmoid_(1e6*(z.-threshold))

        out = σnn .* i + σH .* (1 .- i)  .+ y_

    else
        error("$nn_type does not exist")
    end
end


text:
mutable struct ErgodicPolicySE2 <: FEBOL.Policy
    execution_horizon::Int      # how long should we execute
    horizon_index::Int          #

    em::ErgodicManagerSE2
    tm::TrajectoryManager

    cache::Array{Float64,4}
    visualize::Bool

    xd::VVF
    ud::VVF
end

function ErgodicPolicySE2(f::DF;
                       K::Int = 5,
                       planning_horizon::Int = 50,
                       execution_horizon::Int = 5,
                       visualize::Bool=true
                      )

    # ergodic manager
    N = planning_horizon
    d = Domain([0,0,0], [1,1,2*pi], [f.n,f.n,36])
    em = ErgodicManagerSE2(d, K)

    # trajectory manager
    ci = ConstantInitializer([0.0,0.0])
    h = 1.
    x0 = [0.5,0.5]
    tm = TrajectoryManager(x0,h,N,ci)
    tm.barrier_cost = 1.0

    println("tm.R = ", tm.R)
    tm.R = 0.2*eye(3)
    tm.Q = eye(3)

    cache = make_cache2(f)

    return ErgodicPolicySE2(execution_horizon,
                         9999999,       # index should just start very large
                         em,
                         tm,
                         cache,
                         visualize,
                         [[0.0]],
                         [[0.0]]
                        )
end

function FEBOL.action(m::SearchDomain, x::Vehicle, o, f::DF, p::ErgodicPolicySE2)

    L = m.length
    if p.horizon_index >= p.execution_horizon

        # generate EID from filter
        tic()
        copy!( p.em.phi, mutual_information(f, p.cache) )
        t = toq()
        println("time to eid = ", t)

        # normlize and decompose
        tic()
        normalize!(p.em)
        t = toq()
        println("time to normalize = ", t)
        tic()
        decompose!(p.em)
        t = toq()
        println("time to decompose = ", t)

        # right now, let's just make a new trajectory each time
        p.tm.x0 = [x.x/L, x.y/L, deg2rad(x.heading)]
        tic()
        p.xd, p.ud = pto_trajectory(p.em, p.tm)
        t = toq()
        println("solve time = ", t)

        p.horizon_index = 0
    end

    p.horizon_index += 1

    if p.visualize
        xd2 = deepcopy(p.xd)
        for xi in xd2
            xi[1] *= L
            xi[2] *= L
        end
        figure("Simulation")
        plot_trajectory(xd2)

        figure("MI")
        ErgodicControlPlots.cla()
        plot(p.em, p.xd)
    end

    a = (p.ud[p.horizon_index][1] * L, p.ud[p.horizon_index][2] * L, 0.0)
    println("a = ", a)
    return a
end


text:
# # Introduction
#
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__notebooks/01_Introduction.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__notebooks/01_Introduction.ipynb)
#
# Let's start with a classic:

message = "Hello World!"

#-

println(message)

# Julia is a general purpose programming language designed with scientific
# computing and mathematics in mind. For this reason, some Julia expressions
# look like mathematics.
#
# For example, let's write a Julia function to calculate
# a circle’s circumference $C$ as a function of the radius $r$:
#
# $C = 2 \pi r$

r = 10

C = 2π*r

# `π` is one of the
# [mathematical constants](https://docs.julialang.org/en/v1.1/base/numbers/#General-Number-Functions-and-Constants-1)
# defined in Julia. You can write `π` by typing `\pi` and pressing `<TAB>` in
# the REPL:

π

#-

pi

# The multiplication operator `*` is not needed when a literal number
# (e.g. `2`) is placed just before a variable or constant (e.g. `π`),
# see [*Numeric Literal Coefficients*](https://docs.julialang.org/en/v1.1/manual/integers-and-floating-point-numbers/#man-numeric-literal-coefficients-1).
# That makes polynomial expressions much cleaner:

f(x) = 1.56 + 2.24x + 3.47x^2

#-

f(0.04)

# This one-line function definition, i.e. `function_name(arguments) = body`, is
# the simplest we can use in Julia.

# #### Exercise 1
#
# Define an one-line function to calculate the length of the  hypotenuse of a
# right-angled triangle with catheti (legs) $a$ and $b$ as
# $\sqrt{a^{2} + b^{2}}$ using the `sqrt()` function:

## hypotenuse(...

# If your solution is correct, the following test should pass without errors.

using Test

# The `Test` module exports the `@test` macro, useful for writing unit tests of
# the code. We are going to use `≈`, `\approx<TAB>`, (`a ≈ b` or
# `isapprox(a, b)`) test if two values can be equal taking into account
# floating point errors. :

@test hypotenuse(2, 5) ≈ hypot(2, 5)

# In real life, you can use the function `hypot` to avoid
# [underflow and overflow problems](https://www.johndcook.com/blog/2010/06/02/whats-so-hard-about-finding-a-hypotenuse/).

# ## Control Flow
#
# ### `if`

distance = hypot(1.0, 3.0, 2.5)

#-

if distance <= 2.08  # Å
    "Disulphide"
elseif distance <= 6  # Å
    "Contact"
else
    "Not interacting"
end

# Each condition should be a boolean variable (`true` or `false`) or an
# expression that evaluates to a boolean value.

# ### Short-circuit evaluation

condition = true
condition && println("It's true!")

# #### Exercise 2
#
# Can you write the equivalent expression using `if` ?

## if ...

# ### Ternary operator

condition = true
result = condition ? "👍" : "👎"

# This is equivalent to:

result = if condition
        "👍"
    else
        "👎"
end

#-

result

# ## `for` loops

dna = "ATGCAT"  # dna is a string

for base in dna
    println(base)  # base is a character
end

# `for base = dna` and `for base ∈ dna` are also a possible notations, i.e.
# using `=` or `∈` instead of `in`:

for base ∈ dna
    println(base)
end

# Julia has a nice unicode support and unicode characters can be used in
# variables names, operators, etc. You can write `∈` in the REPL by typing `\in`
# and pressing the `<TAB>` key.

# ## Functions
#
# ### One-line syntax

"Return the number of codons in the sequence."
codon_number(seq) = div(length(seq), 3) # This is a Julia comment

# A strings before the function declaration is used as docstring. You can access
# the function documentation by typing `?` in the REPL and the name of the
# function, e.g.: `?codon_number`

## Try to look at the documentation of condon_number

#-

codon_number(dna)

# ### Multi-line syntax

# Functions with longer bodies are defined using the `function` keyword:

"Returns the GC number in a given sequence."
function count_gc(string)
    count = 0
    for char in string
        if char == 'C' || char == 'G'  # || is the short-circuiting or.
            ## ' ' (instead of " ") is used to define a character.
            count += 1  # i.e. count = count + 1
        end
    end
    count  # i.e. return count
end

#-

count_gc(dna)

# By default, a function returns the result of the last evaluated expression,
# `count` in the previous example. Otherwise, the keyword `return` should be
# used:

"""
`is_dna` iterates a string and returns:
- `false` if a character different from **A, C, T or G** is found
- `true` otherwise
"""
function is_dna(string)
    for char in string
        if char != 'A' && char != 'C' && char != 'T' && char != 'G'
            return false
        end
    end
    true
end

#-

is_dna("ACHL")

# Here, the `return` keyword is used inside the loop to return `false` as soon
# as an incorrect character is found.
#
# `"""` defines multiline strings. In this example, the multiline string is the
# documentation string of the function. Markdown syntax can be used in
# docstrings.

#md ?is_dna
#nb ?is_dna

# #### Exercise 3
#
# Write a function that returns `true` if the given string (argument) has at
# least one `N` using the keywords `function`, `for`, `if` and `return`.

## has_n(...

# `@testset` is useful to aggregate tests for a particular functionality:

using Test
@testset "Exercise 3" begin
    @test has_n("ACTGN")
    @test !has_n("ACTG")
end

# #### Material
#
# You can find more complete information about
# [functions](https://docs.julialang.org/en/v1.1/manual/functions/)
# and [control flow](https://docs.julialang.org/en/v1.1/manual/control-flow/)
# statements in the Julia manual.

max_stars_repo_name: diegozea/JuliaForBioinformatics

max_issues_repo_name: diegozea/JuliaForBioinformatics

max_forks_repo_name: diegozea/JuliaForBioinformatics


text:
using TidalFluxConfigurations, ADCPDataProcessing, PIEMetData, TidalFluxPlots, TidalFluxExampleData, Plots
using Base.Test

TidalFluxConfigurations.config[:_ADCPDATA_DIR] = Pkg.dir("TidalFluxExampleData","data","adcp")
TidalFluxConfigurations.config[:_METDATA_DIR] = Pkg.dir("TidalFluxExampleData","data","met")

creek = Creek{:sweeney}()
deps = parse_deps(creek)
ad1 = load_data(deps[1])

plot(ad1)

cs = parse_cs(creek)
csd = load_data(cs)

# Test the cross-section plot
plot(csd)

h1,Q1 = computedischarge(ad1,csd)

# Test a few Quantities
plot(h1)
plot(Q1)

cals = parse_cals(creek)
cald = load_data(cals[1])

plot(cald)
calscatter([cald])


text:
using JuAFEM
using UnicodePlots

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  grid = generate_grid(Quadrilateral, (20,20))
  addnodeset!(grid, "boundary", x -> abs(x[1]) ≈ 1 ||  abs(x[2]) ≈ 1);

  dim = 2
  ip = Lagrange{dim, RefCube, 1}()
  qr = QuadratureRule{dim, RefCube}(2)
  cellvalues = CellScalarValues(qr, ip);

  dh = DofHandler(grid)
  push!(dh, :T, 1) # Add a temperature field
  close!(dh)

  dbc = DirichletBoundaryConditions(dh)
  add!(dbc, :T, getnodeset(grid, "boundary"), (x,t) -> 0.0)
  close!(dbc)
  update!(dbc, 0.0)

  K = create_sparsity_pattern(dh);
  fill!(K.nzval, 1.0);
  spy(K)

  function doassemble(cellvalues::CellScalarValues{dim}, K::SparseMatrixCSC, dh::DofHandler) where dim
      b = 1.0
      f = zeros(ndofs(dh))
      assembler = start_assemble(K, f)

      n_basefuncs = getnbasefunctions(cellvalues)
      global_dofs = zeros(Int, ndofs_per_cell(dh))

      fe = zeros(n_basefuncs) # Local force vector
      Ke = zeros(n_basefuncs, n_basefuncs) # Local stiffness mastrix

      @inbounds for (cellcount, cell) in enumerate(CellIterator(dh))
          fill!(Ke, 0)
          fill!(fe, 0)

          reinit!(cellvalues, cell)
          for q_point in 1:getnquadpoints(cellvalues)
              dΩ = getdetJdV(cellvalues, q_point)
              for i in 1:n_basefuncs
                  δT = shape_value(cellvalues, q_point, i)
                  ∇δT = shape_gradient(cellvalues, q_point, i)
                  fe[i] += (δT * b) * dΩ
                  for j in 1:n_basefuncs
                      ∇T = shape_gradient(cellvalues, q_point, j)
                      Ke[i, j] += (∇δT ⋅ ∇T) * dΩ
                  end
              end
          end

          celldofs!(global_dofs, cell)
          assemble!(assembler, global_dofs, fe, Ke)
      end
      return K, f
  end;

  K, f = doassemble(cellvalues, K, dh);

  apply!(K, f, dbc)
  T = K \ f;

  vtkfile = vtk_grid("heat", dh, T)
  vtk_save(vtkfile);

  Base.Test.@test maximum(T) ≈ 0.29526786377073544
  println("Heat square successful")

end


text:
using ModelingToolkit: operation, istree, arguments
# Method of lines discretization scheme
struct MOLFiniteDifference{T,T2} <: DiffEqBase.AbstractDiscretization
    dxs::T
    time::T2
    upwind_order::Int
    centered_order::Int
end

# Constructors. If no order is specified, both upwind and centered differences will be 2nd order
MOLFiniteDifference(dxs, time; upwind_order = 1, centered_order = 2) =
    MOLFiniteDifference(dxs, time, upwind_order, centered_order)

function SciMLBase.symbolic_discretize(pdesys::ModelingToolkit.PDESystem,discretization::DiffEqOperators.MOLFiniteDifference)

    pdeeqs = pdesys.eqs isa Vector ? pdesys.eqs : [pdesys.eqs]
    t = discretization.time
    nottime = filter(x->~isequal(x.val, t.val),pdesys.indvars)
    nspace = length(nottime)
    depvars = pdesys.depvars

    # Discretize space
    space = map(nottime) do x
        xdomain = pdesys.domain[findfirst(d->isequal(x.val, d.variables),pdesys.domain)]
        @assert xdomain.domain isa IntervalDomain
        dx = discretization.dxs[findfirst(dxs->isequal(x.val, dxs[1].val),discretization.dxs)][2]
        dx isa Number ? (xdomain.domain.lower:dx:xdomain.domain.upper) : dx
    end
    # Get tspan
    tdomain = pdesys.domain[findfirst(d->isequal(t.val, d.variables),pdesys.domain)]
    @assert tdomain.domain isa IntervalDomain
    tspan = (tdomain.domain.lower,tdomain.domain.upper)

    # Build symbolic variables
    indices = CartesianIndices(((axes(s)[1] for s in space)...,))
    depvarsdisc = map(depvars) do u
        [Num(Variable{Symbolics.FnType{Tuple{Any}, Real}}(Base.nameof(ModelingToolkit.operation(u.val)),II.I...))(t) for II in indices]
    end
    spacevals = map(y->[Pair(nottime[i],space[i][y.I[i]]) for i in 1:nspace],indices)

    ### INITIAL AND BOUNDARY CONDITIONS ###
    # Build symbolic maps for boundaries
    edges = reduce(vcat,[[vcat([Colon() for j in 1:i-1],1,[Colon() for j in i+1:nspace]),
      vcat([Colon() for j in 1:i-1],size(depvarsdisc[1],i),[Colon() for j in i+1:nspace])] for i in 1:nspace])

    #edgeindices = [indices[e...] for e in edges]
    edgevals = reduce(vcat,[[nottime[i]=>first(space[i]),nottime[i]=>last(space[i])] for i in 1:length(space)])
    edgevars = [[d[e...] for e in edges] for d in depvarsdisc]

    bclocs = map(e->substitute.(pdesys.indvars,e),edgevals) # location of the boundary conditions e.g. (t,0.0,y)
    edgemaps = Dict(bclocs .=> [spacevals[e...] for e in edges])
    initmaps = substitute.(depvars,[t=>tspan[1]])

    # Generate map from variable (e.g. u(t,0)) to discretized variable (e.g. u₁(t))
    subvar(depvar) = substitute.((depvar,),edgevals)
    # depvarbcmaps will dictate what to replace the variable terms with in the bcs
    depvarbcmaps = reduce(vcat,[subvar(depvar) .=> edgevar for (depvar, edgevar) in zip(depvars, edgevars)])
    # depvarderivbcmaps will dictate what to replace the Differential terms with in the bcs
    if nspace == 1
        # 1D system
        left_weights(j) = DiffEqOperators.calculate_weights(discretization.upwind_order, space[j][1], space[j][1:2])
        right_weights(j) = DiffEqOperators.calculate_weights(discretization.upwind_order, space[j][end], space[j][end-1:end])
        central_neighbor_idxs(i,j) = [i+CartesianIndex([ifelse(l==j,-1,0) for l in 1:nspace]...),i,i+CartesianIndex([ifelse(l==j,1,0) for l in 1:nspace]...)]
        left_idxs = central_neighbor_idxs(CartesianIndex(2),1)[1:2]
        right_idxs(j) = central_neighbor_idxs(CartesianIndex(length(space[j])-1),1)[end-1:end]
        # Constructs symbolic spatially discretized terms of the form e.g. au₂ - bu₁
        derivars = [[dot(left_weights(1),depvar[left_idxs]), dot(right_weights(1),depvar[right_idxs(1)])]
                    for depvar in depvarsdisc]
        # Create list of all the symbolic Differential terms evaluated at boundary e.g. Differential(x)(u(t,0))
        subderivar(depvar,s) = substitute.((Differential(s)(depvar),),edgevals)
        # Create map of symbolic Differential terms with symbolic spatially discretized terms
        depvarderivbcmaps = reduce(vcat,[subderivar(depvar, s) .=> derivars[i]
                                       for (i, depvar) in enumerate(depvars) for s in nottime])
   else
        # Higher dimension
        # TODO: Fix Neumann and Robin on higher dimension
        depvarderivbcmaps = []
   end

    # Generate initial conditions and bc equations
    u0 = []
    bceqs = []
    for bc in pdesys.bcs
        if ModelingToolkit.operation(bc.lhs) isa Sym && ~any(x -> isequal(x, t.val), ModelingToolkit.arguments(bc.lhs))
            # initial condition
            # Assume in the form `u(...) ~ ...` for now
            push!(u0,vec(depvarsdisc[findfirst(isequal(bc.lhs),initmaps)] .=> substitute.((bc.rhs,),spacevals)))
        else
            # Algebraic equations for BCs
            i = findfirst(x->occursin(x.val,bc.lhs),first.(depvarbcmaps))
            bcargs = first(depvarbcmaps[i]).val.arguments
            # Replace Differential terms in the bc lhs with the symbolic spatially discretized terms
            # TODO: Fix Neumann and Robin on higher dimension
            lhs = nspace == 1 ? substitute(bc.lhs,depvarderivbcmaps[i]) : bc.lhs

            # Replace symbol in the bc lhs with the spatial discretized term
            lhs = substitute(lhs,depvarbcmaps[i])
            rhs = substitute.((bc.rhs,),edgemaps[bcargs])
            lhs = lhs isa Vector ? lhs : [lhs] # handle 1D
            push!(bceqs,lhs .~ rhs)
        end
    end

    u0 = reduce(vcat,u0)
    bceqs = reduce(vcat,bceqs)

    ### PDE EQUATIONS ###
    interior = indices[[2:length(s)-1 for s in space]...]
    eqs = vec(map(Base.product(interior,pdeeqs)) do p
        II,eq = p

        # Create a stencil in the required dimension centered around 0
        # e.g. (-1,0,1) for 2nd order, (-2,-1,0,1,2) for 4th order, etc
        order = discretization.centered_order
        stencil(j) = CartesianIndices(Tuple(map(x -> -x:x, (1:nspace.==j) * (order÷2))))

        # TODO: Generalize central difference handling to allow higher even order derivatives
        # The central neighbour indices should add the stencil to II, unless II is too close
        # to an edge in which case we need to shift away from the edge
        # Calculate buffers
        I1 = oneunit(first(indices))
        Imin = first(indices) + I1 * (order÷2)
        Imax = last(indices) - I1 * (order÷2)
        # Use max and min to apply buffers
        central_neighbor_idxs(II,j) = stencil(j) .+ max(Imin,min(II,Imax))
        central_neighbor_space(II,j) = vec(space[j][map(i->i[j],central_neighbor_idxs(II,j))])
        central_weights(II,j) = DiffEqOperators.calculate_weights(2, space[j][II[j]], central_neighbor_space(II,j))
        central_deriv(II,j,k) = dot(central_weights(II,j),depvarsdisc[k][central_neighbor_idxs(II,j)])

        central_deriv_rules = [(Differential(s)^2)(u) => central_deriv(II,j,k) for (j,s) in enumerate(nottime), (k,u) in enumerate(depvars)]
        valrules = vcat([depvars[k] => depvarsdisc[k][II] for k in 1:length(depvars)],
                        [nottime[j] => space[j][II[j]] for j in 1:nspace])

        # TODO: Use rule matching for nonlinear Laplacian

        # TODO: upwind rules needs interpolation into `@rule`
        #forward_weights(i,j) = DiffEqOperators.calculate_weights(discretization.upwind_order, 0.0, [space[j][i[j]],space[j][i[j]+1]])
        #reverse_weights(i,j) = DiffEqOperators.calculate_weights(discretization.upwind_order, 0.0, [space[j][i[j]-1],space[j][i[j]]])
        #upwinding_rules = [@rule(*(~~a,(Differential(nottime[j]))(u),~~b) => IfElse.ifelse(*(~~a..., ~~b...,)>0,
        #                        *(~~a..., ~~b..., dot(reverse_weights(i,j),depvarsdisc[k][central_neighbor_idxs(i,j)[1:2]])),
        #                        *(~~a..., ~~b..., dot(forward_weights(i,j),depvarsdisc[k][central_neighbor_idxs(i,j)[2:3]]))))
        #                        for j in 1:nspace, k in 1:length(depvars)]

        rules = vcat(vec(central_deriv_rules),valrules)
        substitute(eq.lhs,rules) ~ substitute(eq.rhs,rules)
    end)

    # Finalize
    defaults = pdesys.ps === nothing || pdesys.ps === SciMLBase.NullParameters() ? u0 : vcat(u0,pdesys.ps)
    ps = pdesys.ps === nothing || pdesys.ps === SciMLBase.NullParameters() ? Num[] : first.(pdesys.ps)
    # Combine PDE equations and BC equations
    sys = ODESystem(vcat(eqs,unique(bceqs)),t,vec(reduce(vcat,vec(depvarsdisc))),ps,defaults=Dict(defaults))
    sys, tspan
end

function SciMLBase.discretize(pdesys::ModelingToolkit.PDESystem,discretization::DiffEqOperators.MOLFiniteDifference)
    sys, tspan = SciMLBase.symbolic_discretize(pdesys,discretization)
    simpsys = structural_simplify(sys)
    prob = ODEProblem(simpsys,Pair[],tspan)
end


text:
include("solver_preamble.jl")
import SDPAFamily
factory = with_optimizer(SDPAFamily.Optimizer{Float64}, params=(gammaStar=0.8, maxIteration=200),presolve=true)
config = MOI.Test.TestConfig(atol=1e-5, rtol=1e-5, query=false)
@testset "Linear" begin
    # With `dsos_horn`, the termination status is
    # `INFEASIBLE_OR_UNBOUNDED` instead of `INFEASIBLE`.
    Tests.linear_test(factory, config, ["dsos_horn"])
end
@testset "SDP" begin
    # With `maxcut` and `quartic_infeasible_lyapunov_switched_system`, the dual status is `UNKNOWN_RESULT_STATUS` instead of `MOI.INFEASIBILITY_CERTIFICATE`
    Tests.sd_test(factory, config, ["maxcut", "quartic_infeasible_lyapunov_switched_system"])

end