text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations
gamma = 5/3
equations = CompressibleEulerEquations2D(gamma)

# Initial condition adopted from
# - Yong Liu, Jianfang Lu, and Chi-Wang Shu
#   An oscillation free discontinuous Galerkin method for hyperbolic systems
#   https://tinyurl.com/c76fjtx4
# Mach = 2000 jet
function initial_condition_astro_jet(x, t, equations::CompressibleEulerEquations2D)
  @unpack gamma = equations
  rho = 0.5
  v1 = 0
  v2 = 0
  p =  0.4127
  # add inflow for t>0 at x=-0.5
  # domain size is [-0.5,+0.5]^2
  if (t > 0) && (x[1] ≈ -0.5) && (abs(x[2]) < 0.05)
    rho = 5
    v1 = 800 # about Mach number Ma = 2000
    v2 = 0
    p = 0.4127
  end
  return prim2cons(SVector(rho, v1, v2, p), equations)
end
initial_condition = initial_condition_astro_jet

boundary_conditions = (
                       x_neg=BoundaryConditionDirichlet(initial_condition_astro_jet),
                       x_pos=BoundaryConditionDirichlet(initial_condition_astro_jet),
                       y_neg=boundary_condition_periodic,
                       y_pos=boundary_condition_periodic,
                      )

surface_flux = flux_lax_friedrichs # HLLC needs more shock capturing (alpha_max)
volume_flux  = flux_ranocha # works with Chandrashekar flux as well
polydeg = 3
basis = LobattoLegendreBasis(polydeg)

# shock capturing necessary for this tough example
indicator_sc = IndicatorHennemannGassner(equations, basis,
                                         alpha_max=0.3,
                                         alpha_min=0.0001,
                                         alpha_smooth=true,
                                         variable=density_pressure)
volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
                                                 volume_flux_dg=volume_flux,
                                                 volume_flux_fv=surface_flux)

solver = DGSEM(basis, surface_flux, volume_integral)

coordinates_min = (-0.5, -0.5)
coordinates_max = ( 0.5,  0.5)

mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=6,
                periodicity=(false,true),
                n_cells_max=100_000)
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 0.001)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 5000
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=5000,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

amr_indicator = IndicatorHennemannGassner(semi,
                                          alpha_max=1.0,
                                          alpha_min=0.0001,
                                          alpha_smooth=false,
                                          variable=Trixi.density)

amr_controller = ControllerThreeLevelCombined(semi, amr_indicator, indicator_sc,
                                              base_level=2,
                                              med_level =0, med_threshold=0.0003, # med_level = current level
                                              max_level =8, max_threshold=0.003,
                                              max_threshold_secondary=indicator_sc.alpha_max)

amr_callback = AMRCallback(semi, amr_controller,
                           interval=1,
                           adapt_initial_condition=true,
                           adapt_initial_condition_only_refine=true)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        amr_callback, save_solution)

# positivity limiter necessary for this tough example
stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds=(5.0e-6, 5.0e-6),
                                                     variables=(Trixi.density, pressure))

###############################################################################
# run the simulation
# use adaptive time stepping based on error estimates, time step roughly dt = 1e-7
sol = solve(ode, SSPRK43(stage_limiter!),
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
module Inspire

using Codex.Questionnaires: Items, get_scores
using DataFrames

export mike2scores

# The scoring is determined by guessing by me.
# I expect to have made a few mistakes, but overall the scoring was reasonably after knowing
# the number of items per predictor; thanks to
# https://repository.tudelft.nl/view/tno/uuid:98d97d5c-ff8e-4f26-9fba-f993bc585d39.
# For coping, the last 4 items seem to be reversed and added later by someone.
# Those are ignored.

coping_flexibility = Items(1:36, [])

emotional_stability = Items(
    [41, 44, 49, 50, 52, 54, 56],
    [42, 43, 45, 46, 47, 48, 51, 53, 55]
)

optimism = Items(
    [61, 63, 69],
    [58, 59, 67]
)

social_competence = Items(
    [57, 60, 62, 64, 65, 66, 68, 70],
    []
)

self_efficacy = Items(71:83, [])

self_reflection = Items(
    [88, 90, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103],
    [84, 85, 86, 87, 91, 92, 89]
)

function mike2scores(df::DataFrame)
    new_df = DataFrame(df)
    if !("completed_at" in names(new_df))
        # Cleaning up the 2018 data.
        new_df[:, :completed_at] = repeat([""], nrow(new_df))
        select!(new_df, :id, :completed_at, Not([:id, :completed_at]))
        new_df = filter(:id => x -> !contains(x, "unknown"), new_df)
    end

    new_df[:, :coping_flexibility] = get_scores(new_df, coping_flexibility)
    new_df[:, :emotional_stability] = get_scores(new_df, emotional_stability)
    new_df[:, :optimism] = get_scores(new_df, optimism)
    new_df[:, :social_competence] = get_scores(new_df, social_competence)
    new_df[:, :self_efficacy] = get_scores(new_df, self_efficacy)
    new_df[:, :self_reflection] = get_scores(new_df, self_reflection)

    scored = select(new_df, :id, :completed_at,
        :coping_flexibility, :emotional_stability, :optimism,
        :social_competence, :self_efficacy, :self_reflection
    )
    # This could have been done above, but I mistakenly used `disallowmissing`
    # instead and went from there.
    dropmissing!(scored)
    scored
end

end # module

ext: jl


text:
@testset "Test Basis" begin
    bt = [Spline(), Cheb(), Lin()]
    n = [9, 7, 11]
    a = [-3, 1e-4, -4]
    b = [3, 10.0, 2]
    params = [SplineParams(linspace(-3, 3, 9), 0, 1), ChebParams(7, 1e-4, 10.0), LinParams(linspace(-4, 2, 11), 0)]

    # use univariate constructors
    b1, b2, b3 = map(Basis, bt, n, a, b, params)

    # directly construct multivariate basis using (default) inner constructor
    b_all = Basis{3,BasisFamily,BasisParams}(bt, n, a, b, params)

    # should be equivalent to Basis(b1, b1)
    b_spline3d = Basis(Spline(), [n[1], n[1]], [a[1], a[1]], [b[1], b[1]], [1, 1])

    # should be equivalent to Basis(b2, b2)
    b_cheb3d = Basis(Cheb(), [n[2], n[2]], [a[2], a[2]], [b[2], b[2]])

    # should be equivalent to Basis(b3, b3)
    b_lin3d = Basis(Lin(), [n[3], n[3]], [a[3], a[3]], [b[3], b[3]])

    @testset "constructors" begin
        @test map(x->isa(x, Basis{1}), (b1, b2, b3)) == (true, true, true)
        @test isa(b_all, Basis{3}) == true

        # use convenience outer constructor to not supply type parameter when
        # constructing multi-varaite basis
        @test Basis(bt, n, a, b, params) == b_all

        # use Basis(params) constructor
        @test Basis(params[1]) == b1
        @test Basis(params[2]) == b2
        @test Basis(params[3]) == b3

        # test Basis(b1, b2, b3) format
        @test Basis(b1, b2, b3) == b_all

        # test Basis(params1, params2, params3)
        @test Basis(params...) == b_all

        # test fundefn type methods
        @test b_spline3d == Basis(b1, b1)
        @test b_cheb3d == Basis(b2, b2)
        @test b_lin3d == Basis(b3, b3)

        # test that basis of different dimensions are not equal
        @test (  ==(b_spline3d, b1) ) == false
        @test (  ==(b_cheb3d, b2) ) == false
        @test (  ==(b_lin3d, b3) ) == false

        # test concrete types come out when possible
        @test isa(Basis(b1, b1), Basis{2,Spline,SplineParams}) == true
        @test isa(Basis(b2, b2), Basis{2,Cheb,ChebParams}) == true
        @test isa(Basis(b3, b3), Basis{2,Lin,LinParams}) == true
        @test isa(Basis(b2, b2, b1), Basis{3,BasisFamily,BasisParams}) == true
        @test isa(Basis(b3, b2, b1), Basis{3,BasisFamily,BasisParams}) == true
    end

    @testset "getindex and combining preserves basis" begin
        @test b_all[1] == b1
        @test b_all[2] == b2
        @test b_all[3] == b3
        @test Basis(b_all[1], b_all[2], b_all[3]) == b_all
        for b in [b1, b2, b3]
            @test b[1] == b
            @test_throws ErrorException b[2]
        end
    end

    @testset "ndims, length, and size" begin
        for b in [b1, b2, b3]
            @test ndims(b) == 1
            @test length(b) == b.n[1]
            @test size(b) == (b.n[1],)
            @test size(b, 1) == b.n[1]
            @test_throws ErrorException size(b, 2)
        end

        @test ndims(b_all) == 3
        @test length(b_all) == prod(n)
        @test size(b_all) == tuple(n...)
        @test size(b_all, 1) == n[1]

    end

    @testset "test nodes" begin
        # extract nodes from individual basis
        n1 = nodes(b1)[1]
        n2 = nodes(b2)[1]
        n3 = nodes(b3)[1]

        # extract product nodes as well as nodes along both dimensions using
        # the 3d basis
        n_all, (n1_all, n2_all, n3_all) = nodes(b_all)

        # test the nodes from basis 1 are what we expect, are the same as the
        # corresponding nodes from the 3d basis and have the correct length
        @test n1 == collect(linspace(-3, 3, 9))
        @test n1 == n1_all
        @test length(n1) == b1.n[1]

        # test that the nodes from basis 2 are the same as corresponding nodes
        # from 3d basis and have correct length
        @test n2 == n2_all
        @test length(n2) == b2.n[1]

        # test that the nodes from basis 3 are the same as corresponding nodes
        # from 3d basis and have correct length
        @test n3 == collect(linspace(-4, 2, 11))
        @test length(n3) == b3.n[1]

        # verify that the nodes from combined 3d basis is correct size
        @test size(n_all) == (length(n1)*length(n2)*length(n3), 3)
    end

    # call show (which calls writemime) so I can get 100% coverage :)
    @testset "Printing" begin
        iob = IOBuffer()
        show(iob, ChebParams(10, -1, 1))
        show(iob, SplineParams(10, -1, 1))
        show(iob, LinParams(10, -1, 1))
        show(iob, b_all)
    end
end


text:
const _max_keys = Ref(10_000)

"""
	updatemaxkeys!(n::Int)

	limits the maximum number of keys in statistics of nodes in JSON. Default value is 10_000.
"""
updatemaxkeys!(n::Int) = _max_keys[] = n
max_keys() = _max_keys[]

const _max_len = Ref(10_000)

"""
	updatemaxlen!(n::Int)

	limits the maximum size of string values in statistics of nodes in JSON. Default value is 10_000.
	Longer strings will be trimmed and their length and hash will be appended to retain the uniqueness.
	This is due to some strings being very long and causing the schema to be even order of magnitute larger than needed.
"""
updatemaxlen!(n::Int) = _max_len[] = n
max_len() = _max_len[]


text:
using sb, LinearAlgebra, StaticArrays

# This file implements mean field many body second moment CEID

"""

Structure to hold data for CEID dynamics

"""
mutable struct CEIDOps
    No :: Int64
    ρ  :: SArray{Tuple{2,2},Complex{Float64},2,4}
    μ  :: Vector{SArray{Tuple{2,2},Complex{Float64},2,4}}
    λ  :: Vector{SArray{Tuple{2,2},Complex{Float64},2,4}}
    CRR :: Array{Float64,2}
    CPR :: Array{Float64,2}
    CPP :: Array{Float64,2}
    p  :: Vector{Float64}
    q  :: Vector{Float64}
    ke :: Float64
    function CEIDOps(ρ0, sbm :: SBModel)
        μ = Vector{SArray{Tuple{2,2},Complex{Float64},2,4}}(undef,sbm.No)
        λ = Vector{SArray{Tuple{2,2},Complex{Float64},2,4}}(undef,sbm.No)
        CRR = zeros(Float64,sbm.No,sbm.No)
        CPR = zeros(Float64,sbm.No,sbm.No)
        CPP = zeros(Float64,sbm.No,sbm.No)
        p = zeros(Float64,sbm.No)
        q = zeros(Float64,sbm.No)
        ke = 0.0
         for ν in 1:sbm.No
             μ[ν] = SA[0.0im 0.0im; 0.0im 0.0im]
             λ[ν] = SA[0.0im 0.0im; 0.0im 0.0im]
             if sbm.β == Inf # zero temperature
                 CPP[ν,ν] = 0.5 * HBAR * sbm.ωs[ν]
                 CRR[ν,ν] = 0.5 * HBAR / sbm.ωs[ν]
             else
                 hav = HBAR*sbm.ωs[ν]*(0.5 + 1.0/(exp(HBAR*sbm.ωs[ν]*sbm.β)-1.0))
                 CPP[ν,ν] = hav
                 CRR[ν,ν] = hav / sbm.ωs[ν]^2
             end
         end
        new(sbm.No,ρ0,μ,λ,CRR,CPR,CPP,p,q,ke)
    end
end

"""

Leapfrog bootstrap of data (CEID)

"""
function CEIDbootstrap!(ops :: CEIDOps, oldops :: CEIDOps,
                      dotops :: CEIDOps, dt :: Float64)
    oldops.ρ = ops.ρ - dt * dotops.ρ
    for i in 1:ops.No
        oldops.μ[i] = ops.μ[i] - dt * dotops.μ[i]
        oldops.λ[i] = ops.λ[i] - dt * dotops.λ[i]
    end
    oldops.CRR .= ops.CRR .- dt .* dotops.CRR
    oldops.CPR .= ops.CPR .- dt .* dotops.CPR
    oldops.CPP .= ops.CPP .- dt .* dotops.CPP
    oldops.p .= ops.p .- dt .* dotops.p
    oldops.q .= ops.q .- dt .* dotops.q
    oldops.ke = ops.ke - dt * dotops.ke

    return nothing
end

"""

Leapfrog step forward (CEID)

"""
function CEIDforward!(ops :: CEIDOps, oldops :: CEIDOps,
                    dotops :: CEIDOps, dt :: Float64)
    oldops.ρ = oldops.ρ + 2 * dt * dotops.ρ
    for i in 1:ops.No
        oldops.μ[i] = oldops.μ[i] + 2 * dt * dotops.μ[i]
        oldops.λ[i]= oldops.λ[i] + 2 * dt * dotops.λ[i]
    end
    oldops.CRR .= oldops.CRR .+ (2 * dt) .* dotops.CRR
    oldops.CPR .= oldops.CPR .+ (2 * dt) .* dotops.CPR
    oldops.CPP .= oldops.CPP .+ (2 * dt) .* dotops.CPP
    oldops.p .= oldops.p .+ (2 * dt) .* dotops.p
    oldops.q .= oldops.q .+ (2 * dt) .* dotops.q
    oldops.ke = oldops.ke + 2 * dt * dotops.ke

    return nothing
end

"""

Calculate CEID EOM RHS for K with no coupling between oscillators

"""
function CEIDcalcdots!(dotops :: CEIDOps, ops :: CEIDOps,
                    dt :: Float64, sbm :: SBModel; thermostat = false)

    # Time derivatives of generalized coordinates and momenta ------------------
    fbar = zeros(Float64,sbm.No)
    for ν in 1:ops.No
        fbar[ν] = real(tr(F(sbm,ops.q,ν)*ops.ρ)) -
                  real(tr(K(sbm,ops.q,ν,ν)*ops.μ[ν]))
        dotops.p[ν] = fbar[ν]
        dotops.q[ν] = ops.p[ν]
    end
    #---------------------------------------------------------------------------

    # Time derivative of ρ -----------------------------------------------------
    dotops.ρ = IOHBAR * comm(H(sbm,ops.q,ops.p),ops.ρ)
    for ν in 1:ops.No
        dotops.ρ += - IOHBAR * comm(F(sbm,ops.q,ν),ops.μ[ν]) +
                    0.5 * IOHBAR * ops.CRR[ν,ν] * comm(K(sbm,ops.q,ν,ν),ops.ρ)
    end
    #---------------------------------------------------------------------------

    # Time derivative of μ -----------------------------------------------------
    Threads.@threads for ν in 1:ops.No
        dotops.μ[ν] = ops.λ[ν] + IOHBAR * comm(H(sbm,ops.q,ops.p),ops.μ[ν])
        for νp in 1:ops.No
            dotops.μ[ν] += - IOHBAR * ops.CRR[ν,νp] * comm(F(sbm,ops.q,νp),ops.ρ)
        end
    # Time derivative of λ -----------------------------------------------------
        dotops.λ[ν] = IOHBAR * comm(H(sbm,ops.q,ops.p),ops.λ[ν]) +
                      0.5 * acomm(F(sbm,ops.q,ν),ops.ρ) - tr(F(sbm,ops.q,ν)*ops.ρ)*ops.ρ -
                      0.5 * acomm(K(sbm,ops.q,ν,ν),ops.μ[ν])
        for νp in 1:ops.No
            dotops.λ[ν] += - IOHBAR * ops.CPR[ν,νp] * comm(F(sbm,ops.q,νp),ops.ρ)
        end
    # Time derivatives of mean field second order correlators ------------------
        if !thermostat
            for νp in 1:ops.No
                dotops.CRR[ν,νp] = ops.CPR[ν,νp] + ops.CPR[νp,ν]
                dotops.CPR[ν,νp] = ops.CPP[ν,νp] + real(tr(F(sbm,ops.q,ν)*ops.μ[νp])) -
                                   ops.CRR[ν,νp] * real(tr(K(sbm,ops.q,ν,ν)*ops.ρ))
                dotops.CPP[ν,νp] = real(tr(F(sbm,ops.q,ν)*ops.λ[νp])) +
                                   real(tr(F(sbm,ops.q,νp)*ops.λ[ν])) -
                                   ops.CPR[ν,νp] * real(tr(K(sbm,ops.q,νp,νp)*ops.ρ)) -
                                   ops.CPR[νp,ν] * real(tr(K(sbm,ops.q,ν,ν)*ops.ρ))
            end
        else
            dotops.CRR .= 0.0
            dotops.CPR .= 0.0
            dotops.CPP .= 0.0
        end
    end
    #---------------------------------------------------------------------------

    # Time derivative of quantum kinetic energy --------------------------------
    dotops.ke = 0.0
    for ν in 1:ops.No
        for νp in 1:ops.No
            dotops.ke += real(tr(K(sbm,ops.q,ν,νp)*ops.ρ)) * ops.CPR[ν,νp] +
                            real(0.5 * IOHBAR * ops.CRR[νp,ν] *
                            tr(K(sbm,ops.q,ν,νp)*comm(H(sbm,ops.q,ops.p),ops.ρ)))
        end
    end
    #---------------------------------------------------------------------------

    return nothing
end

"""

Calculate CEID EOM RHS

"""
function CEIDcalcdots_full!(dotops :: CEIDOps, ops :: CEIDOps,
                    dt :: Float64, sbm :: SBModel; thermostat = false)

    # Time derivatives of generalized coordinates and momenta ------------------
    fbar = zeros(Float64,sbm.No)
    for ν in 1:ops.No
        fbar[ν] = real(tr(F(sbm,ops.q,ν)*ops.ρ))
        for νp in 1:ops.No
            fbar[ν] += - real(tr(K(sbm,ops.q,ν,νp)*ops.μ[νp]))
        end
        dotops.p[ν] = fbar[ν]
        dotops.q[ν] = ops.p[ν]
    end
    #---------------------------------------------------------------------------

    # Time derivative of ρ -----------------------------------------------------
    dotops.ρ = IOHBAR * comm(H(sbm,ops.q,ops.p),ops.ρ)
    for ν in 1:ops.No
        dotops.ρ += - IOHBAR * comm(F(sbm,ops.q,ν),ops.μ[ν])
    end
    for ν in 1:ops.No
        for νp in 1:ops.No
            dotops.ρ += 0.5 * IOHBAR * ops.CRR[ν,νp] * comm(K(sbm,ops.q,ν,νp),ops.ρ)
        end
    end
    #---------------------------------------------------------------------------

    # Time derivative of μ -----------------------------------------------------
    Threads.@threads for ν in 1:ops.No
        dotops.μ[ν] = ops.λ[ν] + IOHBAR * comm(H(sbm,ops.q,ops.p),ops.μ[ν])
        for νp in 1:ops.No
            dotops.μ[ν] += - IOHBAR * ops.CRR[ν,νp] * comm(F(sbm,ops.q,νp),ops.ρ)
        end
    # Time derivative of λ -----------------------------------------------------
        dotops.λ[ν] = IOHBAR * comm(H(sbm,ops.q,ops.p),ops.λ[ν])
        dotops.λ[ν] += 0.5 * acomm(F(sbm,ops.q,ν),ops.ρ) - tr(F(sbm,ops.q,ν)*ops.ρ)*ops.ρ
        for νp in 1:ops.No
            dotops.λ[ν] += - IOHBAR * ops.CPR[ν,νp] * comm(F(sbm,ops.q,νp),ops.ρ)
            dotops.λ[ν] += - 0.5 * acomm(K(sbm,ops.q,ν,νp),ops.μ[νp])
        end
    # Time derivatives of mean field second order correlators ------------------
        if !thermostat
            for νp in 1:ops.No
                dotops.CRR[ν,νp] = ops.CPR[ν,νp] + ops.CPR[νp,ν]
                dotops.CPR[ν,νp] = ops.CPP[ν,νp] + real(tr(F(sbm,ops.q,ν)*ops.μ[νp]))
                dotops.CPP[ν,νp] = real(tr(F(sbm,ops.q,ν)*ops.λ[νp])) +
                                   real(tr(F(sbm,ops.q,νp)*ops.λ[ν]))
                for νpp in 1:ops.No
                    dotops.CPR[ν,νp] += - ops.CRR[νpp,νp] * real(tr(K(sbm,ops.q,ν,νpp)*ops.ρ))
                    dotops.CPP[ν,νp] += - ops.CPR[ν,νpp] * real(tr(K(sbm,ops.q,νpp,νp)*ops.ρ)) -
                                          ops.CPR[νp,νpp] * real(tr(K(sbm,ops.q,νpp,ν)*ops.ρ))
                end
            end
        else
            dotops.CRR .= 0.0
            dotops.CPR .= 0.0
            dotops.CPP .= 0.0
        end
    end
    #---------------------------------------------------------------------------

    # Time derivative of quantum kinetic energy --------------------------------
    dotops.ke = 0.0
    for ν in 1:ops.No
        for νp in 1:ops.No
            dotops.ke += real(tr(K(sbm,ops.q,ν,νp)*ops.ρ)) * ops.CPR[ν,νp] +
                            real(0.5 * IOHBAR * ops.CRR[νp,ν] *
                            tr(K(sbm,ops.q,ν,νp)*comm(H(sbm,ops.q,ops.p),ops.ρ)))
        end
    end
    #---------------------------------------------------------------------------

    return nothing
end

"""

CEID total energy

"""
function CEIDenergy(ops :: CEIDOps, sbm :: SBModel)
    energy = real(tr(H(sbm,ops.q,ops.p)*ops.ρ)) + ops.ke[1]
    for ν in 1:ops.No
        energy += 0.5 * (ops.p[ν]^2 + ops.CPP[ν,ν])
        energy += - real(tr(F(sbm,ops.q,ν)*ops.μ[ν]))
    end
    return energy
end

"""

This function integrates the CEID EOMs starting from the
initial condition provided in ops for nsteps of dt.
The store! callback function is called at every timestep.

"""
function RunCEID!(sbm,ops,nsteps,dt,store!,storage; thermostat = false, knondiag = false)

    dotops = CEIDOps(SA[0.0im 0.0im; 0.0im 0.0im],sbm)
    oldops = CEIDOps(SA[0.0im 0.0im; 0.0im 0.0im],sbm)

    # Bootstrap the integrator
    if !knondiag
        sb.CEIDcalcdots!(dotops, ops, dt, sbm; thermostat = thermostat)
    else
        sb.CEIDcalcdots_full!(dotops, ops, dt, sbm; thermostat = thermostat)
    end
    sb.CEIDbootstrap!(ops,oldops,dotops,dt)

    # call results storage function at t = 0
    t = 0.0
    store!(storage,t,ops,dotops,sbm)

    # integrate
    for i in 1:nsteps-1
        t = i*dt
        if !knondiag
            sb.CEIDcalcdots!(dotops, ops, dt, sbm; thermostat = thermostat)
        else
            sb.CEIDcalcdots_full!(dotops, ops, dt, sbm; thermostat = thermostat)
        end
        sb.CEIDforward!(ops,oldops,dotops,dt)
        (ops,oldops) = (oldops,ops)
        store!(storage,t,ops,dotops,sbm)
    end

    return nothing
end

#


text:
export SEdge,SElement,SVertex,subdMesh

mutable struct SElement
    RingNodes::Vector{Int64}
    Edges::Vector{Int64}
end

mutable struct SEdge
    EndPoints::Vector{Int64}
    RelatedPoints::Vector{Int64}
    TwoElems::Vector{Int64}
    orientation::Vector{Int64}
end

mutable struct SVertex
    Elements::Vector{Int64}
    Valence::Int64
    Edges::Vector{Int64}
    Coords::SVector{3,Float64}
end

mutable struct subdMesh
    edges::Vector{SEdge}
    vertices::Vector{SVertex}
    elements::Vector{SElement}
    mesh::Mesh
end
"""
    GSubdMesh(mesh::Mesh)

    Given a linear polygon mesh, construct a data structure for subdivision surfaces
"""
function GSubdMesh(mesh::Mesh)
    D1 = 3
    edges = skeleton(mesh,1)
    faces = skeleton(mesh,2)
    connectMat = connectivity(edges,faces)
    NV, NE, NF = numvertices(mesh), numcells(edges), numcells(faces)
    # nv = NV + NE
    Svertices = Vector{SVertex}(undef,NV)
    Sedges = Vector{SEdge}(undef,NE)
    Selements = Vector{SElement}(undef,NF)
    for V in 1 : numvertices(mesh)
        ElementIndices = Int[]
        EdgeIndices = Int[]
        for (F, Face) in enumerate(cells(faces))
            for iv = 1 : 3
                # if V is a vertex of the face add all the vertices to the mask
                if Face[iv] == V
                    append!(ElementIndices, F)
                end
            end
        end
        for (E, Edge) in enumerate(cells(edges))
            for ie = 1 : 2
                if Edge[ie] == V
                    append!(EdgeIndices,E)
                end
            end
        end
        coords = mesh.vertices[V]
        # delete duplicate vertices
        # RingverticesIndices = unique(RingverticesIndices)
        # calculate valence for vertex V
        valence = length(ElementIndices)
        SV=SVertex(ElementIndices,valence,EdgeIndices,coords)
        # append!(Svertices, SV)
        Svertices[V] = SV
    end
    for (E, Edge) in enumerate(cells(edges))
        EndVerticesIndices = Edge
        FacesIndex = []
        RelatedVerticesIndices = []
        orientation = []
        for (F, Face) in enumerate(cells(faces))
            if connectMat[F,E] != 0
                    # Find out two faces in the mask
                append!(FacesIndex, F)
                append!(orientation,connectMat[F,E])
                    ## Loop over vertices in the faces find out related vertices
                for iv = 1 : length(Face)
                    if Face[iv] != Edge[1] && Face[iv] != Edge[2]
                            append!(RelatedVerticesIndices, Face[iv])
                    end
                end
            end
        end
        SE=SEdge(EndVerticesIndices,RelatedVerticesIndices,FacesIndex,orientation)
        # append!(Sedges,SE)
        Sedges[E] = SE
    end
    for (F, Face) in enumerate(cells(faces))
        VertexIndices = Face
        irreg_count = 0
        irreg_index = 0
        irreg_valence = 0
        EdgesIndices,orientation = find_edges(F,Face,edges,connectMat)
        for iv = 1 : length(Face)
            ivalence = Svertices[VertexIndices[iv]].Valence
            if ivalence != 6
                irreg_count += 1
                irreg_index = iv
                irreg_valence = ivalence
            end
        end
        if irreg_count > 1
            error("An element has more than 1 extraordinary point, please refine mesh")
        elseif irreg_count == 0
            neighborFaces,Nn = find_neighbor(faces,edges,F,EdgesIndices,Sedges)
            neighborFaces2,Mn = find_neighbor2(neighborFaces,faces,edges,VertexIndices,EdgesIndices,Sedges,connectMat)
            # reorder nodes
            # Ringnodes=[Mn[1],Mn[6],Nn[1],VertexIndices[1],Nn[3],Mn[2],VertexIndices[2],VertexIndices[3],Mn[5],Mn[3],Nn[2],Mn[4]]
            Ringnodes=[VertexIndices[1],VertexIndices[2],Nn[1],Mn[1],Mn[6],Nn[3],VertexIndices[3],Nn[2],Mn[3],Mn[2],Mn[4],Mn[5]]
            SF=SElement(Ringnodes,EdgesIndices)
        else
            # print("element $(F) has 1 extraordinary point!")
            Ringnodes = Vector{Int64}(undef,irreg_valence+6)
            Ringnodes[1] = VertexIndices[irreg_index]
            next = [2,3,1]
            next2 = [3,1,2]
            Ringnodes[2] = VertexIndices[next[irreg_index]]
            Ringnodes[irreg_valence+1] = VertexIndices[next2[irreg_index]]
            # if irreg_index == 1
                # Ringnodes[2] = VertexIndices[irreg_index+1]
                # Ringnodes[irreg_valence+1] = VertexIndices[irreg_index+2]
            # elseif irreg_index == 2
            #     Ringnodes[2] = VertexIndices[irreg_index+1]
            #     Ringnodes[irreg_valence+1] = VertexIndices[irreg_index-1]
            # elseif irreg_index == 3
            #     Ringnodes[2] = VertexIndices[irreg_index-2]
            #     Ringnodes[irreg_valence+1] = VertexIndices[irreg_index-1]
            # end
            #start a loop to find the ring of nodes 1:valence-1
            firstEdge = Sedges[EdgesIndices[irreg_index]]
            firstF = F
            for ir = 1:irreg_valence - 2
                twoelement =firstEdge.TwoElems
                tworelatedpoints = firstEdge.RelatedPoints
                neighborF = 0
                neighborV = 0
                if twoelement[1] == firstF
                    neighborF = twoelement[2]
                    neighborV = tworelatedpoints[2]
                elseif twoelement[2] == firstF
                    neighborF = twoelement[1]
                    neighborV = tworelatedpoints[1]
                end
                Ringnodes[ir+2] = neighborV
                neighborface = cells(faces)[neighborF]
                firstF = neighborF
                edgesIndices_n,orientation_n = find_edges(neighborF,neighborface,edges,connectMat)
                for E in edgesIndices_n
                    Edge = cells(edges)[E]
                    if (Edge[1] == Ringnodes[1] && Edge[2] == neighborV) || (Edge[2] == Ringnodes[1] && Edge[1] == neighborV)
                        firstEdge = Sedges[E]
                    end
                end
            end
            secondEdge = Sedges[EdgesIndices[next[irreg_index]]]
            firstF = F
            for ir = 1:3
                twoelement =secondEdge.TwoElems
                tworelatedpoints = secondEdge.RelatedPoints
                neighborF = 0
                neighborV = 0
                if twoelement[1] == firstF
                    neighborF = twoelement[2]
                    neighborV = tworelatedpoints[2]
                elseif twoelement[2] == firstF
                    neighborF = twoelement[1]
                    neighborV = tworelatedpoints[1]
                end
                Ringnodes[irreg_valence + 1 + ir] = neighborV
                neighborface = cells(faces)[neighborF]
                firstF = neighborF
                edgesIndices_n,orientation_n = find_edges(neighborF,neighborface,edges,connectMat)
                for E in edgesIndices_n
                    Edge = cells(edges)[E]
                    if (Edge[1] == Ringnodes[2] && Edge[2] == neighborV)||(Edge[2] == Ringnodes[2] && Edge[1] == neighborV)
                        secondEdge = Sedges[E]
                    end
                end
            end

            thirdEdge = Sedges[EdgesIndices[next2[irreg_index]]]
            firstF = F
            for ir = 1:3
                twoelement =thirdEdge.TwoElems
                tworelatedpoints = thirdEdge.RelatedPoints
                neighborF = 0
                neighborV = 0
                if twoelement[1] == firstF
                    neighborF = twoelement[2]
                    neighborV = tworelatedpoints[2]
                elseif twoelement[2] == firstF
                    neighborF = twoelement[1]
                    neighborV = tworelatedpoints[1]
                end
                if ir == 2
                    Ringnodes[irreg_valence + 6] = neighborV
                elseif ir ==3
                    Ringnodes[irreg_valence + 5] = neighborV
                end

                neighborface = cells(faces)[neighborF]
                firstF = neighborF
                edgesIndices_n,orientation_n = find_edges(neighborF,neighborface,edges,connectMat)
                for E in edgesIndices_n
                    Edge = cells(edges)[E]
                    if (Edge[1] == Ringnodes[irreg_valence+1] && Edge[2] == neighborV)||(Edge[2] == Ringnodes[irreg_valence+1] && Edge[1] == neighborV)
                        thirdEdge = Sedges[E]
                    end
                end
            end
            SF=SElement(Ringnodes,EdgesIndices)
        end
        # append!(Selements,SF)
        Selements[F] = SF
    end
    subdmesh=subdMesh(Sedges,Svertices,Selements,mesh)
end
"""
    find_neighbor(faces,edges,F,EdgesIndices,orientation,Sedges)

    Given a face find out the neighbor elements (share edge) and vertices.
"""
function find_neighbor(faces,edges,F,EdgesIndices,Sedges)
    neighborFaces = Vector{Int64}(undef,3)
    neighborVertices = Vector{Int64}(undef,3)
    for iE = 1:3
        E = EdgesIndices[iE]
        twoelement = Sedges[E].TwoElems
        tworelatedpoints = Sedges[E].RelatedPoints
        if twoelement[1] == F
            neighborFaces[iE] = twoelement[2]
            neighborVertices[iE] = tworelatedpoints[2]
        elseif twoelement[2] == F
            neighborFaces[iE] = twoelement[1]
            neighborVertices[iE] = tworelatedpoints[1]
        end
    end
    neighborFaces,neighborVertices
end

function find_neighbor2(neighborFaces,faces,edges,VertexIndices,EdgesIndices,Sedges,connectMat)
    neighborFaces2 = Vector{Int64}(undef,6)
    neighborVertices2 = Vector{Int64}(undef,6)
    C=[1,2,2,3,3,1]
    for i = 1 : 3
        F = neighborFaces[i]
        Face = cells(faces)[F]
        edgesIndices_n,orientation = find_edges(F,Face,edges,connectMat)
        for E in edgesIndices_n
            Edge = cells(edges)[E]
            if E != EdgesIndices[i] && (Edge[1] == VertexIndices[C[2*i-1]] || Edge[2] == VertexIndices[C[2*i-1]])
                twoelement = Sedges[E].TwoElems
                tworelatedpoints = Sedges[E].RelatedPoints
                if twoelement[1] == F
                    neighborFaces2[2*i-1] = twoelement[2]
                    neighborVertices2[2*i-1] = tworelatedpoints[2]
                elseif twoelement[2] == F
                    neighborFaces2[2*i-1] = twoelement[1]
                    neighborVertices2[2*i-1] = tworelatedpoints[1]
                end
            elseif E != EdgesIndices[i] && (Edge[1] == VertexIndices[C[2*i]] || Edge[2] == VertexIndices[C[2*i]])
                twoelement = Sedges[E].TwoElems
                tworelatedpoints = Sedges[E].RelatedPoints
                if twoelement[1] == F
                    neighborFaces2[2*i] = twoelement[2]
                    neighborVertices2[2*i] = tworelatedpoints[2]
                elseif twoelement[2] == F
                    neighborFaces2[2*i] = twoelement[1]
                    neighborVertices2[2*i] = tworelatedpoints[1]
                end
            end
        end
    end
    neighborFaces2,neighborVertices2
end

function find_edges(F,Face,edges,connectMat)
    EdgesIndices = Vector{Int64}(undef,3)
    orientation = Vector{Int64}(undef,3)
    for (E, Edge) in enumerate(cells(edges))
        if connectMat[F,E] == 1
            if Edge[1] == Face[1]||Edge[2] == Face[2]
                EdgesIndices[1] = E
                orientation[1] = 1
            elseif Edge[1] == Face[2]||Edge[2] == Face[3]
                EdgesIndices[2] = E
                orientation[1] = 1
            elseif Edge[1] == Face[3]||Edge[2] == Face[1]
                EdgesIndices[3] = E
                orientation[1] = 1
            end
        elseif connectMat[F,E] == -1
            if Edge[1] == Face[2]||Edge[2] == Face[1]
                EdgesIndices[1] = E
                orientation[1] = -1
            elseif Edge[1] == Face[3]||Edge[2] == Face[2]
                EdgesIndices[2] = E
                orientation[1] = -1
            elseif Edge[1] == Face[1]||Edge[2] == Face[3]
                EdgesIndices[3] = E
                orientation[1] = -1
            end
        end
    end
    EdgesIndices,orientation
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -3.1544947018849654 <= q <= 4.066926790342683)

                     add_NL_constraint(m, :(swish(-0.7544717811096784 + -0.6289249660430758*swish(-0.401848603561894 + 0.9879850116297502*swish(0.7471741672319268 + 0.1562161389212755*$(x[1]) + -0.06104000939937393*$(x[2]) + -0.24396039580062512*$(x[3]) + -0.35002310972719863*$(x[4]) + -0.47890036372503086*$(x[5])) + -0.6833595346541679*swish(0.07659783311726898 + 0.1400377608602934*$(x[1]) + -0.5920144548543678*$(x[2]) + -0.3715997206073691*$(x[3]) + -0.08423829417739581*$(x[4]) + -0.48744433011047184*$(x[5]))) + -0.656030629474813*swish(-0.9665555719654946 + 0.014645656858022438*swish(0.7471741672319268 + 0.1562161389212755*$(x[1]) + -0.06104000939937393*$(x[2]) + -0.24396039580062512*$(x[3]) + -0.35002310972719863*$(x[4]) + -0.47890036372503086*$(x[5])) + -0.5583055948877655*swish(0.07659783311726898 + 0.1400377608602934*$(x[1]) + -0.5920144548543678*$(x[2]) + -0.3715997206073691*$(x[3]) + -0.08423829417739581*$(x[4]) + -0.48744433011047184*$(x[5])))) + swish(-0.3137145052115069 + -0.6627216812930383*swish(-0.401848603561894 + 0.9879850116297502*swish(0.7471741672319268 + 0.1562161389212755*$(x[1]) + -0.06104000939937393*$(x[2]) + -0.24396039580062512*$(x[3]) + -0.35002310972719863*$(x[4]) + -0.47890036372503086*$(x[5])) + -0.6833595346541679*swish(0.07659783311726898 + 0.1400377608602934*$(x[1]) + -0.5920144548543678*$(x[2]) + -0.3715997206073691*$(x[3]) + -0.08423829417739581*$(x[4]) + -0.48744433011047184*$(x[5]))) + -0.39516873731093893*swish(-0.9665555719654946 + 0.014645656858022438*swish(0.7471741672319268 + 0.1562161389212755*$(x[1]) + -0.06104000939937393*$(x[2]) + -0.24396039580062512*$(x[3]) + -0.35002310972719863*$(x[4]) + -0.47890036372503086*$(x[5])) + -0.5583055948877655*swish(0.07659783311726898 + 0.1400377608602934*$(x[1]) + -0.5920144548543678*$(x[2]) + -0.3715997206073691*$(x[3]) + -0.08423829417739581*$(x[4]) + -0.48744433011047184*$(x[5])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
function TimeDataFrame(f::Function, tdf::TimeDataFrame, args...)
    df = dataframe(tdf)
    f1(x) = f(x, args...)
    fp = firstperiod(tdf)
    return TimeDataFrame(mapcols(f1, df), fp)
end

function lag(x::AbstractVector{T}, k::Int64) where T
    n = length(x)
    y = Vector{Union{T, Missing}}(undef, n)
    view(y, 1:k) .= missing
    view(y, k+1:n) .= view(x, 1:n-k)
    return y
end

lag(tdf::TimeDataFrame, k::Int64) = TimeDataFrame(lag, tdf, k)
lag(x) = lag(x, 1)

function lead(x::AbstractVector{T}, k::Int64) where T
    n = length(x)
    y = Vector{Union{T, Missing}}(undef, n)
    view(y, 1:n-k) .= view(x, k+1:n)
    view(y, n-k+1:n) .= missing
    return y
end

lead(tdf::TimeDataFrame, k::Int64) = TimeDataFrame(lead, tdf, k)
lead(x) = lead(x, 1)

function addmissingfirst!(tdf, k::DatePeriod)
    kv = k.value
    tmp = Vector{Union{Float64, Missing}}(undef, size(tdf,1) + kv)
    for c in eachcol(tdf)
        tmp .= append!(c, repeat([missing], kv))
        c .= circshift(tmp, kv)
    end
end

function addmissingend!(tdf, k::DatePeriod)
    kv = k.value
    for c in eachcol(tdf)
        append!(c, repeat([missing], kv))
    end
end

function addperiodsfirst!(periods, k::DatePeriod)
    kv = k.value
    append!(periods,
            collect(range(periods[1] - typeof(k)(1), length = kv, step=-typeof(k)(1))))
    tmp = circshift(periods, kv)
    periods .= tmp
end

addperiodsend!(periods, k::DatePeriod) =
    append!(periods,
            collect(range(periods[end] + typeof(k)(1), length = k.value, step=typeof(k)(1))))

function iscontinuous(tdfs...)
    for tdf in tdfs
        iscontinuous(tdf) || return false
    end
    return true
end

function iscontinuous(tdf, tdfs...)
    iscontinuous(tdf) || return false
    for t in tdfs
        iscontinuous(t) || return false
    end
    return true
end

function havesameperiods(tdfs::TimeDataFrame...)
    havesameperiods(Val(iscontinuous(tdfs)), tdfs)
end

function havesameperiods(::Val{true}, tdf::TimeDataFrame, tdfs::TimeDataFrame...)
    p1 = firstperiod(tdf)
    p2 = lastperiod(tdf)
    for t in tdfs
        firstperiod(t) == p1 || return false
        lastperiod(t) == p2 || return false
    end
    return true
end

function havesameperiods(::Val{false}, tdf::TimeDataFrame, tdfs::TimeDataFrame...)
    p = periods(tdf)
    for t in tdfs
        periods(t) == p || return false
    end
    return true
end

function align_!(tdf, tdfs...)
    (iscontinuous(tdf, tdfs...)) ? align_!(Val(true), tdf, tdfs...) : align_!(Val(false), tdf, tdfs...)
end

function align_!(::Val{true}, tdf, tdfs...)
    p = periods(tdf)
    p1 = p[1]
    p2 = p[end]
    tp = typeof(p[1])
    for t in tdfs
        p = periods(t)
        if typeof(p[1]) != tp
            error("TimeDataFrames must have the same frequency")
        end
        p1 = min(p1, p[1])
        p2 = max(p2, p[end])
    end
    zeroperiod = p[1] - p[1]
    for t in push!([tdfs...], tdf)
        p = periods(t)
        m = p[1] - p1
        if m > zeroperiod
            allowmissing!(dataframe(t))
            addmissingfirst!(t, m)
            addperiodsfirst!(p, m)
        end
        m = p2 - p[end]
        if m > zeroperiod
            allowmissing!(dataframe(t))
            addmissingend!(t, m)
            addperiodsend!(p, m)
        end
    end
end

function align_!(::Val{false}, tdf, tdfs...)
    error("Not implemented")
end

function align(tdf, tdfs...)
    tdf1 = copy(tdf)
    tdfs1 = [copy(t) for t in tdfs]
    align_!(tdf1, tdfs1...)
    return(tdf1, tdfs1)
end

function align(tdf1::TimeDataFrame, tdf2::TimeDataFrame)
    tdf1a = copy(tdf1)
    tdf2a = copy(tdf2)
    align_!(tdf1a, tdf2a)
    return (tdf1a, tdf2a)
end

function align!(tdf, tdfs...)
    align_!(tdf, tdfs...)
    return(tdf, tdfs)
end

import Base: +, -

function +(x1::TimeDataFrame, y1::TimeDataFrame)
    (TimeDataFrames.ncol(x1) != TimeDataFrames.ncol(y1)) && error("TimeDataFrames must have the same number of columns")
    (x, y) = align(x1, y1)
    if TimeDataFrames.ncol(x) == 1
        return TimeDataFrame( Matrix(x) + Matrix(y), ["$(names(x)[1])_+_$(names(y)[1])"], firstperiod(x))
    else
        names(x) != names(y) && error("Multicolumn TimeDataFrames must have the same column names")
        S = Matrix{Union{Float64, Missing}}(undef, size(x))
        dfx = dataframe(x)
        dfy = dataframe(y)
        for (i, n) in enumerate(names(x))
            view(S, :, i) = getproperty(dfx, Symbol(n)) .+ getproperty(dfy, Symbol(n))
        end
        return TimeDataFrame(S, names(x), firstperiod(x))
    end
end

function -(x1::TimeDataFrame, y1::TimeDataFrame)
    (TimeDataFrames.ncol(x1) != TimeDataFrames.ncol(y1)) && error("TimeDataFrames must have the same number of columns")
    (x, y) = align(x1, y1)
    if TimeDataFrames.ncol(x) == 1
        return TimeDataFrame( Matrix(x) - Matrix(y), ["$(names(x)[1])_-_$(names(y)[1])"], firstperiod(x))
    else
        names(x) != names(y) && error("Multicolumn TimeDataFrames must have the same column names")
        S = Matrix{Union{Float64, Missing}}(undef, size(x))
        dfx = dataframe(x)
        dfy = dataframe(y)
        for (i, n) in enumerate(names(x))
            view(S, :, i) = getproperty(dfx, Symbol(n)) .- getproperty(dfy, Symbol(n))
        end
        return TimeDataFrame(S, names(x), firstperiod(x))
    end
end

function index(tdf::TimeDataFrame, base::ExtendedDates.SimpleDate)
    baseindex = findall(base .== periods(tdf))
    return TimeDataFrame(x -> 100 .* x./x[baseindex], tdf)
end


text:
export graph_monomial, graph_monomial_degopt

"""
     (graph,crefs)=graph_monomial(a; input=:A, polyname=:P)

Generates the graph for the polynomial using the monomial basis coefficients.
More precisely, it corresponds to the evaluation of the polynomial

    p(s) = a[1] + a[2]*s + ... + a[n]*s^(n-1),

where s^k is naively evaluated as s^k=s*s^(k-1) for k=2,3,...,n-1.
The kwarg `polyname` specifies the name of intermediate variables.
"""
function graph_monomial(a; input = :A, polyname = :P)

    # Initial setup
    n = length(a)
    cref = Vector{Tuple{Symbol,Int}}(undef, n)
    graph = Compgraph(eltype(a))
    d = n - 1 # Polynomial degree

    # Degenerate case d = 0
    if d == 0
        outkey = Symbol("$(polyname)2")
        add_lincomb!(graph, outkey, a[1], :I, 0.0, :I)
        cref[1] = (outkey, 1)
        add_output!(graph, outkey)
        return (graph, cref)
    end

    outkey = Symbol("$(polyname)$n")
    nodelist = Vector{Symbol}(undef, n)
    nodelist[1] = :I
    nodelist[2] = Symbol("$input")

    # Create monomial basis
    key = input
    for i = 2:d
        prevkey = key
        key = Symbol("$input$i")
        add_mult!(graph, key, input, prevkey)
        nodelist[i+1] = key
    end

    # Sum up the polynomial
    cref[:] .= add_sum!(graph, outkey, a, nodelist, polyname)

    add_output!(graph, outkey)
    return (graph, cref)
end

"""
     (graph,crefs)=graph_monomial_degopt(a; input=:A)

Generates the same polynomial as [`graph_monomial`](@ref), in the monomial
basis. However, it does so by wrapping a call to [`graph_degopt`](@ref),
resulting in more degrees of freedom in `crefs`.
"""
function graph_monomial_degopt(a; input = :A)
    n = length(a)
    d = n - 1
    T = eltype(a)

    if d == 0
        error("Does not implement degree-zero polynomial.")
    end

    x = Vector{Tuple{Vector{T},Vector{T}}}(undef, d - 1)
    for i = 1:d-1
        x[i] =
            (vcat(zeros(T, i), one(T)), vcat(zero(T), one(T), zeros(T, i - 1)))
    end

    return graph_degopt(x, a, input = input)
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

using Random, StatsBase, Distributions, DensityInterface, ValueShapes
using HypothesisTests

import NestedSamplers

@testset "test_ellipsoidal_nested_sampling" begin
    p = Uniform(-80,80)
    prior = BAT.NamedTupleDist(a=[p,p,p])

    dists = [
        MixtureModel(Normal, [(-50.0,2.5),(50.0,2.5)],[0.2,0.8]),
        MixtureModel(Normal, [(-20.0,2.5),(20.0,2.5)],[0.7,0.3]),
        MixtureModel(Normal, [(-40.0,5.0),(10.0,1.0)],[0.5,0.5])
    ]

    dist = NamedTupleDist(a = product_distribution(dists))

    likelihood = logfuncdensity(logdensityof(dist))

    posterior = PosteriorDensity(likelihood, prior)
    algorithm = EllipsoidalNestedSampling(max_ncalls = 10^5)
    r = BAT.sample_and_verify(posterior, algorithm, dist)
    @test r.verified

    smpls = r.result
    @test logdensityof(posterior).(smpls.v) ≈ smpls.logd

    logz_expected = -log(prod(maximum.(prior.a.v) .- minimum.(prior.a.v)))
    @test isapprox(r.logintegral.val, logz_expected, atol = 100 * r.logintegral.err)

    @test r.ess > 50
end


text:
const treated = [85, 88, 75, 66, 25, 29, 83, 39, 97]
const control = [68, 41, 10, 49, 16, 65, 32, 92, 28, 98]

(better, worse) = permutation_test(treated, control)

tot = better + worse

println("Permutation test using the following data:")
println("Treated:  ", treated)
println("Control:  ", control)
println("\nThere are $tot different permuted groups of these data.")
@printf("%8d, %5.2f%% showed better than actual results.\n", better, 100 * better / tot)
print(@sprintf("%8d, %5.2f%% showed equalivalent or worse results.", worse, 100 * worse / tot))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module Subscripts

export super, sub, recover

super_chars = ['⁰', '¹', '²', '³', '⁴', '⁵', '⁶', '⁷', '⁸', '⁹']
  sub_chars = ['₀', '₁', '₂', '₃', '₄', '₅', '₆', '₇', '₈', '₉']

function super(s)
    res = ""
    for c in s
        if c >= '0' && c <= '9'
            res *= super_chars[c - '0' + 1]
        end
    end
    return res
end

function sub(s)
    res = ""
    for c in s
        if c >= '0' && c <= '9'
            res *= sub_chars[c - '0' + 1]
        end
    end
    return res
end

function recover(s::String)
    res = ""
    for c in s
        for i in 1:length(super_chars)
            if c == super_chars[i] || c == sub_chars[i]
                res *= '0' + i - 1
            end
        end
    end
    return res
end

end # module


text:
using Optim
using Devectorize

"""
Linear Regression Cost function with Regularisation

Parameters:
	params      = vector containing rolled up X and Θ
		X       = array of features for each movie
		Θ       = array of weights for each user
	Y           = array of ratings for each movie
	R           = array of the existance of ratings for each movie
	numUsers    = number of users
	numMovies   = number of movies
	numFeatures = number of features
	λ           = regularisation constant

Returns:
	J = linear regression cost of the given weights with the given regularisation constant
"""
function cost(params, Y::Array{Int}, numUsers, numMovies, numFeatures, λ)
	# Unfold params -> X, Θ
	X, Θ = unroll(params, numUsers, numMovies, numFeatures)

	# Calculate cost
	J = 0
 	@simd for i = 1:size(Y)[2]
		@inbounds J += (dot(X[Y[2, i], :][:], Θ[Y[1, i], :][:]) - Y[3, i])^2
	end
	J = 1/2 * J
	# add regularisation
	J += λ/2 * (sum(Θ.^2) + sum(X.^2))

	return J
end

"""
Linear Regression Gradient with Regularisation

Parameters:
	params      = vector containing rolled up X and Θ
		X       = array of features for each movie
		Θ       = array of weights for each user
	gradient    = vector containing rolled up Xgrad and ΘGrad from previous iteration
	Y           = array of ratings for each movie
	R           = array of the existance of ratings for each movie
	numUsers    = number of users
	numMovies   = number of movies
	numFeatures = number of features
	λ           = regularisation constant

Modifies:
	gradient = XGrad and ΘGrad are updated based on the given weights and regularisation costant
"""
function grad!(params, gradient, Y, numUsers, numMovies, numFeatures, λ)
	@printf "hi"
	# Unfold params -> X, Θ
	X, Θ = unroll(params, numUsers, numMovies, numFeatures)

	XGrad = zeros(size(X))
	ΘGrad = zeros(size(Θ))

	for n = 1:size(Y)[2]
		i = Y[2, n]
		j = Y[1, n]
		@devec x = X[i, :]
		@devec þ = Θ[j, :]
		temp = dot(x[:], þ[:]) - Y[3, n]
		@devec XGrad[i, :] += temp .* þ + λ.*x
		@devec ΘGrad[j, :] += temp .* x + λ.*þ
	end

	gradient[:] = rollup(XGrad, ΘGrad)
end

function costGrad!(params, gradient, Y, numUsers, numMovies, numFeatures, λ)
	# Unfold params -> X, Θ
	X, Θ = unroll(params, numUsers, numMovies, numFeatures)

	XGrad = zeros(size(X))
	ΘGrad = zeros(size(Θ))

	# Calculate cost
	J = 0
	start = rand(1:size(Y)[2] - 5_000_000)
	for n = start:start+5_000_000
		i = Y[2, n]
		j = Y[1, n]
		@devec x = X[i, :]
		@devec þ = Θ[j, :]
		temp = dot(x[:], þ[:]) - Y[3, n]
		J += temp^2
		@devec XGrad[i, :] += temp .* þ + λ.*x
		@devec ΘGrad[j, :] += temp .* x + λ.*þ
	end
	J = 1/2 * J
	# add regularisation
	J += λ/2 * (sum(Θ.^2) + sum(X.^2))

	gradient[:] = rollup(XGrad, ΘGrad)

	return J
end

"""
Unroll parameter vector into X and Θ

Parameters:
	params = vector containing rolled up X and Θ
	numUsers    = number of users
	numMovies   = number of movies
	numFeatures = number of features

Returns:
	X = array of features for each movie
	Θ = array of weights for each user
"""
function unroll(params, numUsers, numMovies, numFeatures)
	X = reshape(params[1 : numMovies * numFeatures], numMovies, numFeatures)
	Θ = reshape(params[numMovies * numFeatures + 1 : end], numUsers, numFeatures)
	return X, Θ
end

"""
Roll up X and Θ into single vector

Parameters:
	X = array of features for each movie
	Θ = array of weights for each user

Returns:
	params = vector containing rolled up X and Θ
"""
function rollup(X, Θ)
	return params = vcat(X..., Θ...)
end

"""
Entry point

Params:
	dataPath = path to() CSVs
"""
function main()
	dataPath = ""
	if ARGS[1] == "small"
		dataPath = "../ml-latest-small/fixed/"
	elseif ARGS[1] == "full"
		dataPath = "../ml-latest/fixed/"
	end

	# Load data from CSV - using wc to speed up loading by preallocating
	# N.B. requires cat and wc to be installed - a shell should not be nessecary though
	lines = parse(readall(pipeline(`cat $[dataPath]ratings.csv`, `wc -l`)))
	@time ratings = readdlm(dataPath * "ratings.csv", ',', header=true, dims=(lines, 4))
	ratingsHeader = ratings[2]
	ratings = ratings[1]

	lines = parse(readall(pipeline(`cat $[dataPath]movies.csv`, `wc -l`)))
	@time movies = readdlm(dataPath * "movies.csv", ',', header=true, dims=(lines, 3))
	moviesHeader = movies[2]
	movies = movies[1]

	lines = parse(readall(pipeline(`cat X-results.csv`, `wc -l`)))
	@time X = readdlm("X-results.csv", ',', header=false, dims=(lines, 18))

	lines = parse(readall(pipeline(`cat Theta-results.csv`, `wc -l`)))
	@time Θ = readdlm("Theta-results.csv", ',', header=false, dims=(lines, 18))

	lines = parse(readall(pipeline(`cat $[dataPath]train.csv`, `wc -l`)))
	@time train = readdlm(dataPath * "train.csv", ',', header=true, dims=(lines, 4))[1]

	lines = parse(readall(pipeline(`cat $[dataPath]test.csv`, `wc -l`)))
	@time test = readdlm(dataPath * "test.csv", ',', header=true, dims=(lines, 4))[1]

	@printf "Finished loading, begin processing\n"

	train[:, 3] *= 2
	train = Array{Int}(train[:, 1:3]')
	test[:, 3] *= 2
	test = Array{Int}(test[:, 1:3]')

	# Count users, movies and features
	numUsers = length(union(ratings[:,1]))
	numMovies = size(movies)[1]

	numFeatures = 18 # Leave this for now, 18 initiall chosen b/c 18 listed genres
	λ = 0.05

	@printf "Finished processing, begin learning\n"

	function costTemp(x)
		cost(x, train, numUsers, numMovies, numFeatures, λ)
	end

	function gradTemp!(x, y)
		grad!(x, y, train, numUsers, numMovies, numFeatures, λ)
	end

	i = 0
	function costGradTemp!(x, y)
		i += 1
		if i % 500 == 0
			@printf "%i/%i \t %f \t %s\n" i 20 * 1000 costTemp(x) now()
		end
		costGrad!(x, y, train, numUsers, numMovies, numFeatures, λ)
	end

	@time @show results = optimize(DifferentiableFunction(costTemp, gradTemp!, costGradTemp!), rollup(X,Θ), iterations=1000)

	@printf "Finished learning, begin saving\n"

	writecsv("X-results.csv", X)
	writecsv("Theta-results.csv", Θ)

	@show score = cost(results.minimum, test, numUsers, numMovies, numFeatures, 0) / size(test)[2]

	X, Θ = unroll(results.minimum, numUsers, numMovies, numFeatures)

	f = open("score.txt", "w")
	write(f, "$(score)", "\n")
	close(f)

	@printf "Done :D\n"
end

main()


text:
@testset "mean_functions" begin
    @testset "CustomMean" begin
        rng, N, D = MersenneTwister(123456), 11, 2
        x = randn(rng, N)
        foo_mean = x->sum(abs2, x)
        f = CustomMean(foo_mean)

        @test map(f, x) == map(foo_mean, x)
        # differentiable_mean_function_tests(f, randn(rng, N), x)
    end
    @testset "ZeroMean" begin
        rng, P, Q, D = MersenneTwister(123456), 3, 2, 4
        # X = ColVecs(randn(rng, D, P))
        x = randn(rng, P)
        f = ZeroMean{Float64}()

        for x in [x]
            @test map(f, x) == zeros(size(x))
            # differentiable_mean_function_tests(f, randn(rng, P), x)
        end
    end
    @testset "ConstMean" begin
        rng, D, N = MersenneTwister(123456), 5, 3
        # X = ColVecs(randn(rng, D, N))
        x = randn(rng, N)
        c = randn(rng)
        m = ConstMean(c)

        for x in [x]
            @test map(m, x) == fill(c, N)
            # differentiable_mean_function_tests(m, randn(rng, N), x)
        end
    end
end


text:
module CompositionsBase

# Use README as the docstring of the module:
@doc let path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    replace(read(path, String), r"^```julia"m => "```jldoctest README")
end CompositionsBase

if !isdefined(Base, Symbol("@var_str"))
    macro var_str(x)
        return Symbol(x)
    end
end

const compose = ∘

# Since `⨟` may not be rendered correctly in all environments, let's
# use ASCII version as the definition and then `⨟` as an alias.  This
# is not symmetric with how `compose` is defined but appropriately
# "opposite."
opcompose(fs...) = ∘(reverse(fs)...)
const var"⨟" = opcompose

"""
    g ⨟ f
    opcompose(g, f)

The opposite composition operator defined as

    g ⨟ f = f ∘ g
    ⨟(f) = f
    ⨟(fs...) = ∘(reverse(fs)...)
"""
(var"⨟", opcompose)

export ∘, compose, opcompose
@eval export $(Symbol("⨟"))  # for Julia 1.0

if isdefined(Base, :ComposedFunction)
    export decompose, deopcompose
    """
        decompose(comp)

    Destructure a composition into its pieces:
    ```jldoctest
    julia> using CompositionsBase

    julia> decompose(sin ∘ cos)
    (sin, cos)

    julia> decompose(sin∘cos∘identity)
    (sin, cos, identity)

    julia> decompose(sin)
    (sin,)
    ```

    Requires at least Julia version 1.6. See also [deopcompose](@ref).
    """
    decompose(f) = (f,)
    decompose(o::Base.ComposedFunction) = (decompose(o.outer)..., decompose(o.inner)...)

    """
        deopcompose(comp)

    Destructure a composition into its pieces:
    ```jldoctest
    julia> using CompositionsBase

    julia> deopcompose(sin ⨟cos⨟tan)
    (sin, cos, tan)

    julia> decompose(sin⨟cos⨟tan)
    (tan, cos, sin)
    ```

    Requires at least Julia version 1.6. See also [decompose](@ref).
    """
    deopcompose(f) = (f,)
    deopcompose(o::Base.ComposedFunction) = (deopcompose(o.inner)..., deopcompose(o.outer)...)
end

end # module


text:
using Clang.Generators
using Clang.JLLEnvs
using CImGuiPack_jll

cd(@__DIR__)

include_dir = joinpath(CImGuiPack_jll.artifact_dir, "include")

cimgui_h = joinpath(include_dir, "cimgui.h") |> normpath
cimplot_h = joinpath(include_dir, "cimplot.h") |> normpath
cimnodes_h = joinpath(include_dir, "cimnodes.h") |> normpath

for target in JLLEnvs.JLL_ENV_TRIPLES
    @info "processing $target"

    options = load_options(joinpath(@__DIR__, "generator.toml"))
    options["general"]["output_file_path"] = joinpath(@__DIR__, "..", "lib", "$target.jl")

    args = get_default_args(target)
    push!(args, "-I$include_dir", "-DCIMGUI_DEFINE_ENUMS_AND_STRUCTS")

    ctx = create_context([cimgui_h, cimplot_h, cimnodes_h], args, options)

    build!(ctx)
end


text:
# Juliaを使った株価の取得と通知
using JSON, DataFrames, Dates, CSV

# 株式名or番号，取得期間，取得間隔を指定してURLを生成
function yahoofinurl(stock, period, interval)
    prefix = "https://query1.finance.yahoo.com/v7/finance/chart/"
    url = prefix * stock * "?range=" * period * "&interval=" * interval * "&indicators=quote&includeTimestamps=true"
    println(url)
    return url
end

# Yahoo!FinanceのJSONからDataFrameに変換
function stockjson2df(json; dropnan=true)
    # 時刻の変換
    unixtime = json["chart"]["result"][1]["timestamp"]
    dt = Dates.unix2datetime.(unixtime) + Dates.Hour(9) # 時刻は日本時間としている．米国時間の場合"+Dates.Hour(9)"を削る

    # DataFrameに代入
    df = DataFrame()
    df."timestamp" = dt
    df."open" = json["chart"]["result"][1]["indicators"]["quote"][1]["open"]
    df."low" = json["chart"]["result"][1]["indicators"]["quote"][1]["low"]
    df."high" = json["chart"]["result"][1]["indicators"]["quote"][1]["high"]
    df."close" = json["chart"]["result"][1]["indicators"]["quote"][1]["close"]
    df."volume" = json["chart"]["result"][1]["indicators"]["quote"][1]["volume"]

    # nothingを削除あるいはNaNに変換
    idx = (df."open" .== nothing) # nothingのある行
    if dropnan
        # 1557.TなどNothingが含まれる物がある．DataFrameからnothingを含む行を削除する．
        df = df[.!idx,:]
    else
        # nothingのある行をNaNにする
        df[idx,2:end] = ones(sum(idx), 5) .* NaN
    end

    # 数値をFloat64にする
    df2 = hcat(df."timestamp", convert.(Float64, df[:,2:end]))
    rename!(df2, "x1" => "timestamp")

    return df2
end

# 株価をYahoo!Financeから取得する．
# stock: 株価ティッカー
# period: 株価取得期間．"1d"なら1日，"1y"なら1年
# interval: 株価間隔．"1m"なら1分
function getstockprice(stock, period="1d", interval="1m")
    # URLの作成
    url = yahoofinurl(stock, period, interval)

    # データ取得・JSON->DataFrameへ変換
    str = read(download(url), String) # urlからjson文字列の取得
    json = JSON.parse(str) # json文字列をパースしてDict型に変換
    df = stockjson2df(json) # JSONからDataFrameに変換
end

# 過去の1日単位のS&P500関連の株価を最大期間(30年分)取得する
function getstockpricepast()
    # 取得対象の銘柄の指定．東証の株価は".T"をつける．
    stocks = ["1557.T" "SPY" "SPXL" "SPXS"]
    periods = ["10y" "30y" "30y" "30y"]    # 取得期間 1y, 2y, 1d, 2dという指定方法．1557は10yまでしか無い模様．
    savestockprice.(stocks, periods)
end

# 1日足データ，過去最大期間のデータを取得し，CSVで保存する．
function savestockprice(stock, period)
    interval = "1d" # 取得データの間隔
    df = getstockprice(stock, period, interval)

    # ファイル名用の時刻文字列
    now = Dates.format(df."timestamp"[end], "_yyyymmdd_HHMMSS")

    # CSV書き込み前に時刻のフォーマットを変更しておく
    df."timestamp" = Dates.format.(df."timestamp", "yyyy/mm/dd HH:MM:SS")

    # CSVへの書き込み
    println("write stock price: " * stock * now * ".csv")
    CSV.write(stock * now * ".csv", df)
end

# 現在のS&P500関連の株価を取得する
function getstockpricenow(stock)
    price = getstockprice(stock, "10m", "1m")."close"[end]
end


text:
export IDAPBCProblem, PDELossKinetic, PDELossPotential
export kineticpde, potentialpde
export trainable, unstack, paramstack
export hamiltonian, hamiltoniand, interconnection, controller
export gradient

const UFA = Union{T,C} where {T<:Function, C<:Chain}
const J2Chain = NTuple{N,<:Chain} where {N}
const J2FChain = NTuple{N,<:FastChain} where {N}

struct IDAPBCProblem{TJ2,M,MD,V,VD,GT,GP}
    N::Int
    M⁻¹::M
    Md⁻¹::MD
    V::V
    Vd::VD
    J2::TJ2
    G::GT
    G⊥::GP
    function IDAPBCProblem(N,M⁻¹,Md⁻¹,V,Vd,J2,G,G⊥)
        if isa(M⁻¹,Matrix) && (isa(Md⁻¹,Matrix) || isa(Md⁻¹,PSDMatrix)) && !isnothing(J2)
            error("Constant M and Md requires J2=nothing.")
        end
        if isa(M⁻¹,Function) && isa(Md⁻¹,Matrix)
            error("M is a function, so Md cannot be constant.")
        end
        for f in (isnothing(J2) ? (Md⁻¹, Vd) : (Md⁻¹, Vd, J2...))
            if isa(f, Chain)
                valid = (
                    isa(f.layers, DenseLayers),
                    isa(front(f.layers), DenseLayers),
                    isa(front(tail(f.layers)), DenseLayers)
                )
                if !any(valid)
                    error("Unsupported type $(f)")
                end
            elseif isa(f, FastChain)
                valid = (
                    isa(f.layers, FastDenseLayers),
                    isa(front(f.layers), FastDenseLayers),
                    isa(front(tail(f.layers)), FastDenseLayers)
                )
                if !any(valid)
                    error("Unsupported type $(f)")
                end
            end
        end
        new{typeof(J2),typeof(M⁻¹),typeof(Md⁻¹),typeof(V),typeof(Vd),typeof(G),typeof(G⊥)}(
            N,M⁻¹,Md⁻¹,V,Vd,J2,G,G⊥
        )
    end
end

function IDAPBCProblem(N,M⁻¹,Md⁻¹,V,Vd,G,G⊥)
    IDAPBCProblem(N,M⁻¹,Md⁻¹,V,Vd,nothing,G,G⊥)
end

function Base.show(io::IO, p::IDAPBCProblem)
    println(io, "IDAPBCProblem [q ∈ ℜ^$(p.N)]");
    println(io, "M⁻¹  => $(typeof(p.M⁻¹).name.name)");
    println(io, "Md⁻¹ => $(typeof(p.Md⁻¹).name.name)");
    println(io, "Vd   => $(typeof(p.Vd).name.name)");
    println(io, "J2   => $(typeof(p.J2).name.name)");
end
function Base.getindex(P::IDAPBCProblem, sym::Symbol)
    if sym === :Md
        if isimplicit(P)
            return (q,::Any=missing)->inv(_Md⁻¹(P,q))
        else
            return (q,ps)->inv(_Md⁻¹(P,q,ps))
        end
    elseif sym === :Mdinv
        if isimplicit(P)
            return (q,::Any=missing)->_Md⁻¹(P,q)
        else
            return (q,ps)->_Md⁻¹(P,q,ps)
        end
    elseif sym === :Vd
        if isimplicit(P)
            return (q,::Any=missing)->_Vd(P,q)
        else
            return (q,ps)->_Vd(P,q,ps)
        end
    elseif sym === :Hd
        if isimplicit(P)
            function (x,::Any=missing)
                q, qdot = x[1:P.N], x[P.N+1:end]
                momentum = _M⁻¹(P,q) \ qdot
                hamiltoniand(P, [q; momentum])
            end
        else
            function (x,ps)
                q, qdot = x[1:P.N], x[P.N+1:end]
                momentum = _M⁻¹(P,q) \ qdot
                hamiltoniand(P, [q; momentum], ps)
            end
        end
    elseif sym === :J2
        isnothing(P.J2) && return nothing
        if isimplicit(P)
            return (q,::Any=missing)->interconnection(P,q)
        else
            return (q,ps)->interconnection(P,q,ps)
        end
    end
end

hasfreevariables(::IDAPBCProblem{J2}) where {J2} = J2===Nothing
isimplicit(P::IDAPBCProblem) = any(f->isa(f,Chain), trainable(P))

function kineticpde(M⁻¹, Md⁻¹, ∇M⁻¹, ∇Md⁻¹, G⊥, J2=0)
    return G⊥ * (∇M⁻¹' - (Md⁻¹\M⁻¹)*∇Md⁻¹' + J2*Md⁻¹)
end

function potentialpde(M⁻¹, Md⁻¹, ∇V, ∇Vd, G⊥)
    return first(G⊥ * (∇V - (Md⁻¹\M⁻¹)*∇Vd))
end

trainable(p::IDAPBCProblem{Nothing,M,MD,V,VD}) where
    {M,V,MD<:UFA,VD<:UFA} = (p.Md⁻¹, p.Vd)
trainable(p::IDAPBCProblem{Nothing,M,MD,V,VD}) where
    {M,V,MD<:Matrix,VD<:UFA} = (p.Vd,)
trainable(p::IDAPBCProblem{J,M,MD,V,VD}) where
    {M,V,MD<:UFA,VD<:UFA,N,J<:NTuple{N,<:UFA}} = (p.Md⁻¹, p.Vd, p.J2...)
function unstack(ts::NTuple{N,Function}, ps) where N
    n = DiffEqFlux.paramlength(first(ts))
    res = ps[1:n]
    return (res, unstack(tail(ts), ps[n+1:end])...)
end
unstack(::Tuple{}, ::Any) = ()
unstack(p::IDAPBCProblem, ps) = unstack(trainable(p), ps)
unstack(::IDAPBCProblem, ::Nothing) = nothing
paramstack(p::IDAPBCProblem) = vcat(DiffEqFlux.initial_params.(trainable(p))...)

_M⁻¹(P::IDAPBCProblem{J2,M}, ::Any) where {J2,M<:Matrix} = P.M⁻¹
_M⁻¹(P::IDAPBCProblem{J2,M}, q) where {J2,M<:Function} = P.M⁻¹(q)
_Md⁻¹(P::IDAPBCProblem{J2,M,MD}, q, ::Any=nothing) where {J2,M,MD<:Matrix} = P.Md⁻¹
_Md⁻¹(P::IDAPBCProblem{J2,M,MD}, q) where {J2,M,MD<:Chain} = P.Md⁻¹(q)
function _Md⁻¹(P::IDAPBCProblem{J2,M,MD}, q) where {J2,M,MD<:Function}
    P.Md⁻¹(q)
end
function _Md⁻¹(P::IDAPBCProblem{J2,M,MD}, q, ps) where {J2,M,MD<:Function}
    θMd = first(unstack(P, ps))
    P.Md⁻¹(q, θMd)
end
_∇M⁻¹(P::IDAPBCProblem{J2,M}, q) where {J2,M<:Matrix} = collect(0P.M⁻¹ for _=1:P.N)
_∇M⁻¹(P::IDAPBCProblem{J2,M}, q) where {J2,M<:Function} = jacobian(P.M⁻¹,q)
_∇Md⁻¹(P::IDAPBCProblem{J2,M,MD}, q, ::Any=nothing) where {J2,M,MD<:Matrix} = collect(0P.M⁻¹ for _=1:P.N)
_∇Md⁻¹(P::IDAPBCProblem{J2,M,MD}, q) where {J2,M,MD<:Chain} = jacobian(P.Md⁻¹,q)
function _∇Md⁻¹(P::IDAPBCProblem{J2,M,MD}, q, ps) where {J2,M,MD<:Function}
    θMd = first(unstack(P, ps))
    jacobian(P.Md⁻¹, q, θMd)
end
function _∇Md⁻¹(P::IDAPBCProblem{J2,M,MD}, q) where {J2,M,MD<:Function}
    jacobian(P.Md⁻¹, q)
end
_Vd(P::IDAPBCProblem{J2,M,MD,V,VD}, q) where {J2,M,MD,V,VD<:Chain} = P.Vd(q)
function _Vd(P::IDAPBCProblem{J2,M,MD,V,VD}, q, ps) where {J2,M,MD<:Matrix,V,VD<:Function}
    P.Vd(q, ps)
end
function _Vd(P::IDAPBCProblem{J2,M,MD,V,VD}, q, ps) where {J2,M,MD<:Function,V,VD<:Function}
    _, θVd = unstack(P, ps)
    P.Vd(q, θVd)
end
function _∇Vd(P::IDAPBCProblem{J2,M,MD,V,VD}, q) where {J2,M,MD,V,VD<:Chain}
    if isa(last(P.Vd.layers), Flux.Dense) && isequal(size(last(P.Vd.layers).weight,1), P.N)
        return P.Vd(q)
    else
        return jacobian(P.Vd, q)[:]
    end
end
function _∇Vd(P::IDAPBCProblem{J2,M,MD,V,VD}, q, ps) where {J2,M,V,MD<:Matrix,VD<:FastChain}
    if isa(last(P.Vd.layers), DiffEqFlux.FastDense) && isequal(last(P.Vd.layers).out, P.N)
        return P.Vd(q, ps)
    else
        return jacobian(P.Vd, q, ps)[:]
    end
end
function _∇Vd(P::IDAPBCProblem{J2,M,MD,V,VD}, q, ps) where {J2,M,V,MD<:Function,VD<:FastChain}
    _, θVd = unstack(P,ps)
    if isa(last(P.Vd.layers), DiffEqFlux.FastDense) && isequal(last(P.Vd.layers).out, P.N)
        return P.Vd(q, θVd)
    else
        return jacobian(P.Vd, q, θVd)[:]
    end
end
function _∇Vd(P::IDAPBCProblem{J2,M,MD,V,VD}, q, ps) where {J2,M,V,MD<:Matrix,VD<:Function}
    return jacobian(P.Vd, q, ps)[:]
end
function _∇Vd(P::IDAPBCProblem{J2,M,MD,V,VD}, q, ps) where {J2,M,V,MD<:Function,VD<:Function}
    _, θVd = unstack(P,ps)
    return jacobian(P.Vd, q, θVd)[:]
end

interconnection(P::IDAPBCProblem{Nothing}, x, ::Any=nothing) = 0
function interconnection(P::IDAPBCProblem{<:J2Chain}, x)
    N = P.N
    q, p = x[1:N], x[N+1:2N]
    mapreduce(+, P.J2, p) do Uk, pk
        Uk(q)*pk/2
    end
end
function interconnection(P::IDAPBCProblem{<:J2FChain}, x, ps)
    N = P.N
    q, p = x[1:N], x[N+1:2N]
    θ = last(unstack(P, ps), 2)
    mapreduce(+, P.J2, p, θ) do Uk, pk, θk
        Uk(q,θk)*pk/2
    end
end

function hamiltonian(P::IDAPBCProblem, qp)
    N = P.N
    q, p = qp[1:N], qp[N+1:2N]
    dot(p, _M⁻¹(P, q)*p)/2 + first(P.V(q))
end
function ∇H(P::IDAPBCProblem, qp)
    N = P.N
    q, p = qp[1:N], qp[N+1:2N]
    ∇M⁻¹ = _∇M⁻¹(P, q)
    (mapreduce(*,+,∇M⁻¹,p)'*p)/2 + jacobian(P.V, q)[:]
end

function hamiltoniand(P::IDAPBCProblem, qp)
    N = P.N
    q, p = qp[1:N], qp[N+1:2N]
    dot(p, _Md⁻¹(P, q)*p)/2 + first(_Vd(P,q))
end
function hamiltoniand(P::IDAPBCProblem, qp, ps)
    N = P.N
    q, p = qp[1:N], qp[N+1:2N]
    dot(p, _Md⁻¹(P, q, ps)*p)/2 + first(_Vd(P,q,ps))
end
function ∇Hd(P::IDAPBCProblem, qp)
    N = P.N
    q, p = qp[1:N], qp[N+1:2N]
    ∇Md⁻¹ = _∇Md⁻¹(P, q)
    (mapreduce(*,+,∇Md⁻¹,p)'*p)/2 + _∇Vd(P, q)
end
function ∇Hd(P::IDAPBCProblem, qp, ps)
    N = P.N
    q, p = qp[1:N], qp[N+1:2N]
    ∇Md⁻¹ = _∇Md⁻¹(P, q, ps)
    (mapreduce(*,+,∇Md⁻¹,p)'*p)/2 + _∇Vd(P, q, ps)
end

function controller(P::IDAPBCProblem, x; kv=1)
    N = P.N
    q, qdot = x[1:N], x[N+1:2N]
    G⁻¹ = (P.G'*P.G)\(P.G')
    M⁻¹ = _M⁻¹(P,q)
    Md⁻¹ = _Md⁻¹(P,q)
    MdM⁻¹ = Md⁻¹ \ M⁻¹
    momentum = M⁻¹ \ qdot
    qp = [q; momentum]
    u_es = G⁻¹ * ( ∇H(P,qp) - MdM⁻¹*∇Hd(P,qp) + interconnection(P,qp)*Md⁻¹*momentum )
    u_di = -kv * dot(P.G, Md⁻¹*momentum)
    return u_es + u_di
end
function controller(P::IDAPBCProblem, x, ps; kv=1)
    N = P.N
    q, qdot = x[1:N], x[N+1:2N]
    G⁻¹ = (P.G'*P.G)\(P.G')
    M⁻¹ = _M⁻¹(P,q)
    Md⁻¹ = _Md⁻¹(P,q,ps)
    MdM⁻¹ = Md⁻¹ \ M⁻¹
    momentum = M⁻¹ \ qdot
    qp = [q; momentum]
    u_es = G⁻¹ * ( ∇H(P,qp) - MdM⁻¹*∇Hd(P,qp,ps) + interconnection(P,qp,ps)*Md⁻¹*momentum )
    u_di = -kv * dot(P.G, Md⁻¹*momentum)
    return u_es + u_di
end

abstract type IDAPBCLoss{P} end

function gradient(l::IDAPBCLoss{P}, q, ps::AbstractVector) where
    {J2,M,MD,V,VD<:Function,P<:IDAPBCProblem{J2,M,MD,V,VD}}
    ReverseDiff.gradient(_2->l(q,_2), ps)
end
function gradient(l::IDAPBCLoss{P}, q, ps::Flux.Params) where
    {J2,M,MD,V,VD<:Chain,P<:IDAPBCProblem{J2,M,MD,V,VD}}
    Flux.gradient(()->l(q), ps)
end

struct PDELossKinetic{hasJ2,P} <: IDAPBCLoss{P}
    prob::P
    function PDELossKinetic(prob::P) where {J2,P<:IDAPBCProblem{J2}}
        new{J2!==Nothing,P}(prob)
    end
end
function (L::PDELossKinetic{false})(q)
    p = L.prob
    M⁻¹ = _M⁻¹(p, q)
    ∇M⁻¹ = _∇M⁻¹(p, q)
    Md⁻¹ = _Md⁻¹(p, q)
    ∇Md⁻¹ = _∇Md⁻¹(p, q)
    r = map(∇M⁻¹, ∇Md⁻¹) do _1, _2
        sum(abs, kineticpde(M⁻¹, Md⁻¹, _1, _2, p.G⊥))
    end
    sum(r)
end
function (L::PDELossKinetic{true})(q)
    p = L.prob
    M⁻¹ = _M⁻¹(p, q)
    ∇M⁻¹ = _∇M⁻¹(p, q)
    Md⁻¹ = _Md⁻¹(p, q)
    ∇Md⁻¹ = _∇Md⁻¹(p, q)
    r = map(∇M⁻¹, ∇Md⁻¹, p.J2) do _1, _2, Uk
        sum(abs, kineticpde(M⁻¹, Md⁻¹, _1, _2, p.G⊥, Uk(q)))
    end
    sum(r)
end
function (L::PDELossKinetic{false})(q,ps)
    p = L.prob
    M⁻¹ = _M⁻¹(p, q)
    ∇M⁻¹ = _∇M⁻¹(p, q)
    Md⁻¹ = _Md⁻¹(p, q, ps)
    ∇Md⁻¹ = _∇Md⁻¹(p, q, ps)
    r = map(∇M⁻¹, ∇Md⁻¹) do _1, _2
        sum(abs, kineticpde(M⁻¹, Md⁻¹, _1, _2, p.G⊥))
    end
    sum(r)
end
function (L::PDELossKinetic{true})(q,ps)
    p = L.prob
    M⁻¹ = _M⁻¹(p, q)
    ∇M⁻¹ = _∇M⁻¹(p, q)
    Md⁻¹ = _Md⁻¹(p, q, ps)
    ∇Md⁻¹ = _∇Md⁻¹(p, q, ps)
    J2ps = last(unstack(p, ps), 2)
    r = map(∇M⁻¹, ∇Md⁻¹, p.J2, J2ps) do _1, _2, Uk, Ukps
        sum(abs, kineticpde(M⁻¹, Md⁻¹, _1, _2, p.G⊥, Uk(q,Ukps)))
    end
    sum(r)
end

struct PDELossPotential{P} <: IDAPBCLoss{P}
    prob::P
end
function (L::PDELossPotential)(q)
    p = L.prob
    ∇V = jacobian(p.V, q)[:]
    ∇Vd = _∇Vd(p, q)
    M⁻¹ = _M⁻¹(p, q)
    Md⁻¹ = _Md⁻¹(p, q)
    abs(potentialpde(M⁻¹, Md⁻¹, ∇V, ∇Vd, p.G⊥))
end
function (L::PDELossPotential)(q,ps)
    p = L.prob
    ∇V = jacobian(p.V, q)[:]
    ∇Vd = _∇Vd(p, q, ps)
    M⁻¹ = _M⁻¹(p, q)
    Md⁻¹ = _Md⁻¹(p, q, ps)
    abs(potentialpde(M⁻¹, Md⁻¹, ∇V, ∇Vd, p.G⊥))
end

struct PotentialHessianSymLoss{P} <: IDAPBCLoss{P}
    prob::P
    function PotentialHessianSymLoss(p::IDAPBCProblem{J,M,MD,V,VD}) where
        {J,M,MD,V,VD<:Chain}
        if isequal(length(last(p.Vd.layers).bias), p.N)
            new{typeof(p)}(p)
        else
            error("Not applicable for this type of IDAPBCProblem")
        end
    end
    function PotentialHessianSymLoss(p::IDAPBCProblem{J,M,MD,V,VD}) where
        {J,M,MD,V,VD<:FastChain}
        if isequal(last(p.Vd.layers).out, p.N)
            new{typeof(p)}(p)
        else
            error("Not applicable for this type of IDAPBCProblem")
        end
    end
end
function (L::PotentialHessianSymLoss)(q)
    J = jacobian(L.prob.Vd, q)
    mapreduce(abs, +, J - J')
end
function (L::PotentialHessianSymLoss{P})(q, ps) where {J2,M,MD<:Matrix,P<:IDAPBCProblem{J2,M,MD}}
    J = jacobian(L.prob.Vd, q, ps)
    mapreduce(abs, +, J - J')
end
function (L::PotentialHessianSymLoss{P})(q, ps) where {J2,M,MD<:Function,P<:IDAPBCProblem{J2,M,MD}}
    _, θVd = unstack(P,ps)
    J = jacobian(L.prob.Vd, q, θVd)
    mapreduce(abs, +, J - J')
end


text:
# TODO
# Only for 1dof joints
mutable struct LQR{T,N,NK} <: Controller
    K::Vector{Vector{SMatrix{1,NK,T,NK}}} # for each time step and each eqc

    xd::Vector{SVector{3,Float64}}
    vd::Vector{SVector{3,Float64}}
    qd::Vector{UnitQuaternion{T}}
    ωd::Vector{SVector{3,Float64}}

    eqcids::Vector{Integer}
    Fτd::Vector{SVector{1,Float64}}

    control!::Function

    function LQR(A, Bu, Bλ, G, Q, R, horizon, eqcids, xd, vd, qd, ωd, Fτd, Δt, ::Type{T}; controlfunction::Function = control_lqr!) where {T}
        Q = cat(Q...,dims=(1,2))*Δt
        R = cat(R...,dims=(1,2))*Δt

        N = horizon/Δt
        if N<Inf
            N = Integer(ceil(horizon/Δt))
            Ntemp = N
        else
            Ntemp = Integer(ceil(10/Δt)) # 10 second time horizon as maximal horizon for convergence for Inf
        end

        # calculate K
        if size(G)[1] == 0
            @assert size(Bλ)[2] ==0
            if N==Inf
                K = dlqr(A, Bu, Q, R, N)
                Ku = [[K[i:i,:] for i=1:size(K)[1]]]
            else
                Ku = dlqr(A,Bu,zeros(size(A)[1],0),zeros(0,size(A)[1]),Q,R,N)
            end
        else
            Ku = dlqr(A, Bu, Bλ, G, Q, R, Ntemp)
            if N == Inf
                Ku[1] != Ku[2] && @info "Riccati recursion did not converge."
                Ku = [Ku[1]]
            end
        end

        new{T, N, size(Ku[1][1])[2]}(Ku, xd, vd, qd, ωd, eqcids, Fτd, controlfunction)
    end

    function LQR(mechanism::Mechanism{T,Nn,Nb}, bodyids::AbstractVector{<:Integer}, eqcids::AbstractVector{<:Integer},
            Q::Vector{<:AbstractMatrix{T}}, R::Vector{<:AbstractMatrix{T}}, horizon;
            xd::Vector{<:AbstractVector{T}} = [SA{T}[0; 0; 0] for i=1:Nb],
            vd::Vector{<:AbstractVector{T}} = [SA{T}[0; 0; 0] for i=1:Nb],
            qd::Vector{UnitQuaternion{T}} = [one(UnitQuaternion{T}) for i=1:Nb],
            ωd::Vector{<:AbstractVector{T}} = [SA{T}[0; 0; 0] for i=1:Nb],
            Fτd::Vector{<:AbstractVector{T}} = [SA{T}[0] for i=1:length(eqcids)],
            controlfunction::Function = control_lqr!
        ) where {T, Nn, Nb}

        @assert length(bodyids) == length(Q) == length(xd) == length(vd) == length(qd) == length(ωd) == Nb "Missmatched length for bodies"
        @assert length(eqcids) == length(R) == length(Fτd) "Missmatched length for constraints"

        # linearize
        A, Bu, Bλ, G = linearsystem(mechanism, xd, vd, qd, ωd, Fτd, bodyids, eqcids)

        LQR(A, Bu, Bλ, G, Q, R, horizon, eqcids, xd, vd, qd, ωd, Fτd, mechanism.Δt, T, controlfunction = controlfunction)
    end

    function LQR(mechanism::Mechanism{T,Nn,Nb}, controlledids::AbstractVector{<:Integer}, controlids::AbstractVector{<:Integer},
        Q::AbstractVector{T}, R::AbstractVector{T}, horizon;
        xθd::AbstractVector{T} = szeros(T,length(controlledids)),
        vωd::AbstractVector{T} = szeros(T,length(controlledids)),
        Fτd::AbstractVector{T} = szeros(T,length(controlids)),
        controlfunction::Function = control_lqr!
    ) where {T, Nn, Nb}

        @assert length(controlledids) == length(Q) == length(xθd) == length(vωd) == Nb "Missmatched length for bodies"
        @assert length(controlids) == length(R) == length(Fτd) "Missmatched length for constraints"

        # linearize
        A, Bu, Bλ, G, xd, vd, qd, ωd = linearsystem(mechanism, xθd, vωd, Fτd, controlledids, controlids)

        Q = [diagm(ones(12))*Q[i] for i=1:length(Q)]
        R = [diagm(ones(1))*R[i] for i=1:length(R)]

        LQR(A, Bu, Bλ, G, Q, R, horizon, controlids, xd, vd, qd, ωd, [[Fτd[i]] for i=1:length(Fτd)], mechanism.Δt, T, controlfunction = controlfunction)
    end
end

function control_lqr!(mechanism::Mechanism{T,Nn,Nb}, lqr::LQR{T,N}, k) where {T,Nn,Nb,N}
    Δz = zeros(T,Nb*12)
    qvm = QuatVecMap()
    for (id,body) in pairs(mechanism.bodies)
        colx = (id-1)*12+1:(id-1)*12+3
        colv = (id-1)*12+4:(id-1)*12+6
        colq = (id-1)*12+7:(id-1)*12+9
        colω = (id-1)*12+10:(id-1)*12+12

        state = body.state
        Δz[colx] = state.xsol[2]-lqr.xd[id]
        Δz[colv] = state.vsol[2]-lqr.vd[id]
        Δz[colq] = rotation_error(state.qsol[2],lqr.qd[id],qvm)
        Δz[colω] = state.ωsol[2]-lqr.ωd[id]
    end

    if k<N
        for (i,id) in enumerate(lqr.eqcids)
            u = lqr.Fτd[i] - lqr.K[k][i]*Δz
            setForce!(mechanism, geteqconstraint(mechanism, id), u)
        end
    end

    return
end

function control_lqr!(mechanism::Mechanism{T,Nn,Nb}, lqr::LQR{T,Inf}, k) where {T,Nn,Nb}
    Δz = zeros(T,Nb*12)
    qvm = QuatVecMap()
    for (id,body) in pairs(mechanism.bodies)
        colx = (id-1)*12+1:(id-1)*12+3
        colv = (id-1)*12+4:(id-1)*12+6
        colq = (id-1)*12+7:(id-1)*12+9
        colω = (id-1)*12+10:(id-1)*12+12

        state = body.state
        Δz[colx] = state.xsol[2]-lqr.xd[id]
        Δz[colv] = state.vsol[2]-lqr.vd[id]
        Δz[colq] = rotation_error(state.qsol[2],lqr.qd[id],qvm)
        Δz[colω] = state.ωsol[2]-lqr.ωd[id]
    end

    for (i,id) in enumerate(lqr.eqcids)
        u = lqr.Fτd[i] - lqr.K[1][i]*Δz
        setForce!(mechanism, geteqconstraint(mechanism, id), u)
    end

    return
end

function dlqr(A,Bu,Bλ,G,Q,R,N)
    mx = size(A)[2]
    mu = size(Bu)[2]
    mλ = size(Bλ)[2]
    Ku = [[zeros(1,size(Q)[1]) for j=1:mu] for i=1:N-1]
    Kλ = [[zeros(1,size(Q)[1]) for j=1:mλ] for i=1:N-1]
    Pk = Q

    k = 0
    for outer k=N-1:-1:1
        D = Bu - Bλ/(G*Bλ)*G*Bu
        M11 = R + D'*Pk*Bu
        M12 = D'*Pk*Bλ
        M21 = G*Bu
        M22 = G*Bλ

        M = [M11 M12;M21 M22]
        b = [D'*Pk;G]*A

        Kk = M\b

        for i=1:mu
            Ku[k][i] = Kk[i:i,:]
        end

        Kuk = Kk[1:mu,:]
        Kλk = Kk[mu+1:mu+mλ,:]

        Abar = A-Bu*Kuk-Bλ*Kλk
        Pkp1 = Q + Kuk'*R*Kuk + Abar'*Pk*Abar

        if norm(Pk-Pkp1) < 1e-5
            break
        end

        Pk = Pkp1
    end

    for k2=k-1:-1:1
        Ku[k2] = Ku[k2+1]
    end

    return Ku
end


text:
using Documenter
using QuranTree

makedocs(
    sitename = "QuranTree.jl",
    format = Documenter.HTML(
        assets = ["assets/logo.ico"]
    ),
    modules = [QuranTree],
    authors = "Al-Ahmadgaid B. Asaad",
    pages = [
        "Home" => "index.md",
        "Getting Started" => "man/getting_started.md",
        "Indexing" => "man/indexing.md",
        "Transliteration" => "man/transliteration.md",
        "Morphological Features" => "man/morphological_features.md",
        "Data Processing" => "man/data_processing.md",
        "Natural Language Processing" => [
            "Introduction" => "man/nlp/nlp.md",
            "Text Summarization" => "man/nlp/text_summarization.md",
            "Topic Modeling" => "man/nlp/topic_modeling.md",
        ],
        "CAMeL Tools" => [
            "Getting Started" => "man/camel/external.md",
            "Morphological Analysis" => "man/camel/analysis.md",
            "POS Tagging" => "man/camel/tagger.md",
            "Disambiguation" => "man/camel/disambig.md"
        ],
        "API" => "man/api.md"
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/alstat/QuranTree.jl.git"
)


text:
# Autogenerated wrapper script for openPMD_api_jll for armv7l-linux-gnueabihf-cxx03-julia_version+1.6.0
export libopenPMD, libopenPMD_jl

using CompilerSupportLibraries_jll
using MPICH_jll
using libcxxwrap_julia_jll
JLLWrappers.@generate_wrapper_header("openPMD_api")
JLLWrappers.@declare_library_product(libopenPMD, "libopenPMD.so")
JLLWrappers.@declare_library_product(libopenPMD_jl, "libopenPMD.jl.so")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll, MPICH_jll, libcxxwrap_julia_jll)
    JLLWrappers.@init_library_product(
        libopenPMD,
        "lib/libopenPMD.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopenPMD_jl,
        "lib/libopenPMD.jl.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Autogenerated wrapper script for SARSOP_jll for x86_64-linux-gnu-cxx11
export polgraph, pomdpsim, pomdpeval, pomdpconvert, pomdpsol

## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"

# Relative path to `polgraph`
const polgraph_splitpath = ["bin", "polgraph"]

# This will be filled out by __init__() for all products, as it must be done at runtime
polgraph_path = ""

# polgraph-specific global declaration
function polgraph(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(polgraph_path)
    end
end

# Relative path to `pomdpsim`
const pomdpsim_splitpath = ["bin", "pomdpsim"]

# This will be filled out by __init__() for all products, as it must be done at runtime
pomdpsim_path = ""

# pomdpsim-specific global declaration
function pomdpsim(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(pomdpsim_path)
    end
end

# Relative path to `pomdpeval`
const pomdpeval_splitpath = ["bin", "pomdpeval"]

# This will be filled out by __init__() for all products, as it must be done at runtime
pomdpeval_path = ""

# pomdpeval-specific global declaration
function pomdpeval(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(pomdpeval_path)
    end
end

# Relative path to `pomdpconvert`
const pomdpconvert_splitpath = ["bin", "pomdpconvert"]

# This will be filled out by __init__() for all products, as it must be done at runtime
pomdpconvert_path = ""

# pomdpconvert-specific global declaration
function pomdpconvert(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(pomdpconvert_path)
    end
end

# Relative path to `pomdpsol`
const pomdpsol_splitpath = ["bin", "pomdpsol"]

# This will be filled out by __init__() for all products, as it must be done at runtime
pomdpsol_path = ""

# pomdpsol-specific global declaration
function pomdpsol(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(pomdpsol_path)
    end
end

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"SARSOP")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    # We first need to add to LIBPATH_list the libraries provided by Julia
    append!(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)])
    global polgraph_path = normpath(joinpath(artifact_dir, polgraph_splitpath...))

    push!(PATH_list, dirname(polgraph_path))
    global pomdpsim_path = normpath(joinpath(artifact_dir, pomdpsim_splitpath...))

    push!(PATH_list, dirname(pomdpsim_path))
    global pomdpeval_path = normpath(joinpath(artifact_dir, pomdpeval_splitpath...))

    push!(PATH_list, dirname(pomdpeval_path))
    global pomdpconvert_path = normpath(joinpath(artifact_dir, pomdpconvert_splitpath...))

    push!(PATH_list, dirname(pomdpconvert_path))
    global pomdpsol_path = normpath(joinpath(artifact_dir, pomdpsol_splitpath...))

    push!(PATH_list, dirname(pomdpsol_path))
    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(LIBPATH_list, ':')

    # Add each element of LIBPATH to our DL_LOAD_PATH (necessary on platforms
    # that don't honor our "already opened" trick)
    #for lp in LIBPATH_list
    #    push!(DL_LOAD_PATH, lp)
    #end
end  # __init__()


text:
# Some of the common imports that all demonstrations are nice to have.

using LinearAlgebra, SparseArrays
using Plots
using CSV, ProgressMeter


text:
#=
Given a list of integers S and a target number k, write a function that returns a subset of S that adds up to k. If such a subset cannot be made, then return null.

Integers can appear more than once in the list. You may assume all numbers in the list are positive.

For example, given S = [12, 1, 61, 5, 9, 2] and k = 24, return [12, 9, 2, 1] since it sums up to 24.
=#

using Test
include("Solutions/problem42_subset_adds_to_k.jl")

@test subarraywithsum([12, 1, 61, 5, 9, 2], 24) == [12, 9, 2, 1]
@test subarraywithsum([12, 1, 61, 5, 9, 2], 19) == [12, 5, 2]
@test isnothing(subarraywithsum([12, 1, 61, 5, 9, 2], 40))

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem