text:
module LogFixPoint16s

    export LogFixPoint16, nan, power2

    include("logfixpoint16.jl")
    include("change_format.jl")

end


text:
@testset verbose=true "wigner_matrices" begin
    include("wigner_matrices/small_d.jl")
    include("wigner_matrices/H.jl")
    include("wigner_matrices/big_D.jl")
    include("wigner_matrices/sYlm.jl")
end

ext: jl


text:
println("Hellgfggfgfggf")

2023+323

θ = 10
θ*5


text:
function create_graph_from_TTC(routes_dict::OrderedDict{Symbol, OrderedDict{Int,Tuple{Float64,Float64}}}
								, n_stops::Int
								, secs_between_stops::Dict{Tuple{Int,Int},Int64}
								, plot_graph::Bool = false
								, write_to_pdf::Bool = false)
	"""
	Description:
		Takes into account real data regarding TTC routes, connections, and distances between stops and create Simple Weighted Directed Graph based on this information.
		Graph info: Vertices --> stops | Edges --> routes | Weights --> seconds between stops
	Returns:
		SimpleWeightedDiGraph with `n_stops` vertices.
	"""
    # basic graph with original routes
    gg = SimpleWeightedDiGraph(n_stops)
    for route in keys(routes_dict)
        for stop in 1:size(collect(keys(routes_dict[route])))[1]-1
            stop_time = secs_between_stops[(collect(keys(routes_dict[route]))[stop], collect(keys(routes_dict[route]))[stop+1])]
            stop_a = collect(keys(routes_dict[route]))[stop]
            stop_b = collect(keys(routes_dict[route]))[stop+1]
            add_edge!(gg, stop_a, stop_b, stop_time)
        end
    end

	# Check if stops on the routes were correctly connected.
	@assert adjacency_matrix(gg)[1,2] > 0.0  "SOMETHING WENT WRONG 1-->2 should have a non-zero weight!"
    @assert adjacency_matrix(gg)[2,1] == 0.0 "SOMETHING WENT WRONG 2-->1 should have zero weight!"
    @assert adjacency_matrix(gg)[616,615] > 0.0 "SOMETHING WENT WRONG 616-->615 should have a non-zero weight!"
    @assert adjacency_matrix(gg)[615,616] == 0.0 "SOMETHING WENT WRONG 615-->616 should have zero weight!"

    ### Adding artificial edges between same location stops ###

    # dictionary mapping node => LLA
    node_id_LLA = Dict{Int,Tuple{Float64,Float64}}()
    for route in collect(values(routes_dict))
        for (stop_id, stop_LLA) in route
           push!(node_id_LLA, stop_id => stop_LLA)
        end
    end

    # identifing same location nodes and assigning transfer times between them
    transfer_time_dict = Dict{Tuple{Int,Int},Int}()
    for (key, value) in node_id_LLA
        obj = [k for k in keys(node_id_LLA) if node_id_LLA[k] == value]
        obj = vec(collect(Iterators.product(obj, obj))) # array
        obj = filter(((xx, yy),)->xx!=yy, obj) # remove same node distances
        obj = Dict(zip(obj,repeat([120], length(obj))))
        merge!(transfer_time_dict, obj)
    end

    # adding fake edges to the graph
    for (key, value) in transfer_time_dict
        add_edge!(gg, key[1], key[2], value)
    end

    # all of the times in 1 object
    merge!(secs_between_stops, transfer_time_dict);

	# PDF or plot graph if needed
    if write_to_pdf == true
        draw(PDF("stops_subway_streetcars_Toronto_graph.pdf", 16cm, 16cm), gplot(gg))
    else
        if plot_graph == true
            gplot(gg)
        end
    end

    println("Graph from TTC data was made")
    return gg
end
println("Graph from TTC functions were read")

max_stars_repo_name: NykPol/EpidemicInUrbanNetworkToronto

max_issues_repo_name: NykPol/EpidemicInUrbanNetworkToronto

max_forks_repo_name: NykPol/EpidemicInUrbanNetworkToronto


text:
abstract type AbstractChargeDensityModel{T <: SSDFloat} end

@inline function ChargeDensityModel(T::DataType, dict::Union{Dict{String, Any}, Dict{Any, Any}}, inputunit_dict::Dict)
    return ChargeDensityModel(T, Val{Symbol(dict["name"])}(), dict, inputunit_dict)
end

"""
    struct LinearChargeDensityModel{T <: SSDFloat} <: AbstractChargeDensityModel{T}

Simple charge density model which assumes a linear gradient in charge density in each dimension.
`offsets::NTuple{3, T}` are the charge densities at 0 and `gradients::NTuple{3, T}` are the linear
slopes in each dimension.
"""
struct LinearChargeDensityModel{T <: SSDFloat} <: AbstractChargeDensityModel{T}
    offsets::NTuple{3, T}
    gradients::NTuple{3, T}
end

function get_charge_density(lcdm::LinearChargeDensityModel{T}, pt::AbstractCoordinatePoint{T})::T where {T <: SSDFloat}
    ρ::T = 0
    for i in eachindex(lcdm.offsets)
        ρ += (lcdm.offsets[i] + pt[i] * lcdm.gradients[i]) #* T(1e16) # * T(1e10) * T(1e6) -> 1/cm^3 -> 1/m^3
    end
    return ρ
end

"""
    struct ZeroChargeDensityModel{T <: SSDFloat} <: AbstractChargeDensityModel{T}

Returns always 0.
"""
struct ZeroChargeDensityModel{T <: SSDFloat} <: AbstractChargeDensityModel{T} end

function get_charge_density(lcdm::ZeroChargeDensityModel{T}, pt::AbstractCoordinatePoint{T})::T where {T <: SSDFloat}
    ρ::T = 0
    return ρ
end

"""
    struct ConstantChargeDensityModel{T <: SSDFloat} <: AbstractChargeDensityModel{T}

Returns always a fixed charge density.
"""
struct ConstantChargeDensityModel{T <: SSDFloat} <: AbstractChargeDensityModel{T}
    ρ::T
end

function get_charge_density(cdm::ConstantChargeDensityModel{T}, pt::AbstractCoordinatePoint{T})::T where {T <: SSDFloat}
    return cdm.ρ
end

function ChargeDensityModel(T::DataType, t::Val{:linear}, dict::Union{Dict{String, Any}, Dict{Any, Any}}, inputunit_dict::Dict)
    unit_factor::T = 1
    gradient_unit_factor::T = 1
    if haskey(inputunit_dict, "length")
        lunit = inputunit_dict["length"]
        unit_factor = inv(ustrip(uconvert( internal_length_unit^3, 1 * lunit^3 )))
        gradient_unit_factor = inv(ustrip(uconvert( internal_length_unit^4, 1 * lunit^4 )))
    end
    return LinearChargeDensityModel{T}( dict, unit_factor, gradient_unit_factor )
end

function LinearChargeDensityModel{T}(dict::Union{Dict{String, Any}, Dict{Any, Any}}, unit_factor::T, gradient_unit_factor::T)::LinearChargeDensityModel{T} where {T <: SSDFloat}
    offsets, gradients = zeros(T,3), zeros(T,3)
    if haskey(dict, "r")     offsets[1] = geom_round(unit_factor * T(dict["r"]["init"]));     gradients[1] = geom_round(gradient_unit_factor * T(dict["r"]["gradient"]))    end
    #if haskey(dict, "phi")   offsets[2] = geom_round(unit_factor * T(dict["phi"]["init"]));   gradients[2] = geom_round(unit_factor * T(dict["phi"]["gradient"]))  end
    if haskey(dict, "z")     offsets[3] = geom_round(unit_factor * T(dict["z"]["init"]));     gradients[3] = geom_round(gradient_unit_factor * T(dict["z"]["gradient"]))    end
    if haskey(dict, "x")     offsets[1] = geom_round(unit_factor * T(dict["x"]["init"]));     gradients[1] = geom_round(gradient_unit_factor * T(dict["x"]["gradient"]))    end
    if haskey(dict, "y")     offsets[2] = geom_round(unit_factor * T(dict["y"]["init"]));     gradients[2] = geom_round(gradient_unit_factor * T(dict["y"]["gradient"]))    end
    LinearChargeDensityModel{T}( NTuple{3, T}(offsets), NTuple{3, T}(gradients) )
end

function ChargeDensityModel(T::DataType, t::Val{:constant}, dict::Union{Dict{String, Any}, Dict{Any, Any}}, inputunit_dict::Dict)
    unit_factor::T = 1
    if haskey(inputunit_dict, "length")
        lunit = inputunit_dict["length"]
        unit_factor = inv(ustrip(uconvert( internal_length_unit^3, 1 * lunit^3 )))
    end
    return ConstantChargeDensityModel{T}( dict, unit_factor )
end
function ConstantChargeDensityModel{T}(dict::Union{Dict{String, Any}, Dict{Any, Any}}, unit_factor::T)::ConstantChargeDensityModel{T} where {T <: SSDFloat}
    ρ::T = if haskey(dict, "charge_density")
        geom_round(unit_factor * T(dict["charge_density"]))
    else
        T(0)
    end
    ConstantChargeDensityModel{T}( ρ )
end


text:
# utility functions that provide Gaussian odometry accumulation

import IncrementalInference: getFactorMean

export getFactorMean
export accumulateDiscreteLocalFrame!, duplicateToStandardFactorVariable, extractDeltaOdo, resetFactor!
export odomKDE

"""
    $SIGNATURES

Return a number increment on symbol.  Example :x2 -> :x3.
"""
nextPose(sym::Symbol, identifier::Union{String, Char}=string(sym)[1]) = Symbol(string(identifier,parse(Int,string(sym)[2:end])+1))

getFactorMean(fct::PriorPose2) = getFactorMean(fct.Z)
getFactorMean(fct::Pose2Pose2) = getFactorMean(fct.z)
getFactorMean(fct::MutablePose2Pose2Gaussian) = getFactorMean(fct.Zij)

"""
    $SIGNATURES

Advance an odometry factor as though integrating an ODE -- i.e. ``X_2 = X_1 ⊕ ΔX``. Accepts continuous domain process noise density `Qc` which is internally integrated to discrete process noise Qd.  ``DX`` is assumed to already be incrementally integrated before this function.  See related `accumulateContinuousLocalFrame!` for fully continuous system propagation.

Notes
- This update stays in the same reference frame but updates the local vector as though accumulating measurement values over time.
- Kalman filter would have used for noise propagation: ``Pk1 = F*Pk*F' + Qdk``
- From Chirikjian, Vol.II, 2012, p.35: Jacobian SE(2), Jr = [cθ sθ 0; -sθ cθ 0; 0 0 1] -- i.e. dSE2/dX' = SE2([0;0;-θ])
- `DX = dX/dt*Dt`
- assumed process noise for `{}^b Qc = {}^b [x;y;yaw] = [fwd; sideways; rotation.rate]`

Dev Notes
- TODO many operations here can be done in-place.

Related

accumulateContinuousLocalFrame!, accumulateDiscreteReferenceFrame!
"""
function accumulateDiscreteLocalFrame!(mpp::MutablePose2Pose2Gaussian,
                                       DX::Vector{Float64},
                                       Qc::Matrix{Float64},
                                       dt::Float64=1.0;
                                       Fk = SE2([0;0;-DX[3]]),
                                       Gk = Matrix{Float64}(LinearAlgebra.I, 3,3),
                                       Phik = SE2(DX) )
  #
  kXk1 = SE2(mpp.Zij.μ)*Phik
  phi, gamma, Qd = cont2disc(Fk, Gk, Qc, dt, Phik)
  Covk1 = Phik*(mpp.Zij.Σ.mat)*(Phik') + Qd
  check = norm(Covk1 - Covk1')
  1e-4 < check ? @warn("Covk1 is getting dangerously non-Hermitian, still forcing symmetric covariance matrix.") : nothing
  @assert check < 1.0
  Covk1 .+= Covk1'
  Covk1 ./= 2
  mpp.Zij = MvNormal(se2vee(kXk1), Covk1)
  nothing
end

"""
    $SIGNATURES

Helper function to duplicate values from a special factor variable into standard factor and variable.  Returns the name of the new factor.

Notes:
- Developed for accumulating odometry in a `MutablePosePose` and then cloning out a standard PosePose and new variable.
- Does not change the original MutablePosePose source factor or variable in any way.
"""
function duplicateToStandardFactorVariable(::Type{Pose2Pose2},
                                           mpp::MutablePose2Pose2Gaussian,
                                           dfg::AbstractDFG,
                                           prevsym::Symbol,
                                           newsym::Symbol;
                                           solvable::Int=1,
                                           autoinit::Bool=true,
                                           cov::Union{Nothing, Matrix{Float64}}=nothing  )::Symbol
  #
  # extract factor values and create PosePose object
  posepose = Pose2Pose2(MvNormal(mpp.Zij.μ, cov==nothing ? mpp.Zij.Σ.mat : cov))

  # modify the factor graph
  addVariable!(dfg, newsym, Pose2, solvable=solvable)
  addFactor!(dfg, [prevsym; newsym], posepose, solvable=solvable, autoinit=autoinit)
  return ls(dfg, newsym)[1]
end

"""
    $SIGNATURES

Reset the transform value stored in a `::MutablePose2Pose2Gaussian` to zero.
"""
function resetFactor!(mpp::MutablePose2Pose2Gaussian)::Nothing
  mpp.Zij = MvNormal(zeros(3), 1e-6*Matrix{Float64}(LinearAlgebra.I, 3,3) )
  nothing
end

"""
    $SIGNATURES

Extract deltas from existing dead reckoning data so that odometry calculations can be repeated later.

Notes
- Useful for reverse engineering data or simulation tools.
"""
function extractDeltaOdo(XX, YY, TH)
  dt = 1.0
  DX = zeros(3,length(XX))
  nXYT__ = zeros(3,size(DX,2))
  nXYT__[:,1] = [XX[1];YY[1];TH[1]]
  for i in 2:length(XX)
    wTbk = SE2([XX[i-1];YY[i-1];TH[i-1]])
    wTbk1 = SE2([XX[i];YY[i];TH[i]])
    bkTbk1 = wTbk\wTbk1
    DX[:,i] = se2vee(bkTbk1)

    # test
    nXYT__[:,i] .= se2vee(SE2(nXYT__[:,i-1])*SE2(DX[:,i]))
    # nXYT__[:,i] .= se2vee(SE2(nXYT__[:,i-1])*bkTbk1)
  end

  return DX
end

## Previous methods

function odomKDE(p1,dx,cov)
  @warn "odomKDE is beig deprecated in its current form, consider using approxConv or predictVariableByFactor instead."
  X = getPoints(p1)
  sig = diag(cov)
  RES = zeros(size(X))
  # increases the number of particles based on the number of modes in the measurement Z
  for i in 1:size(X,2)
      ent = [randn()*sig[1]; randn()*sig[2]; randn()*sig[3]]
      RES[:,i] = addPose2Pose2(X[:,i], dx + ent)
  end
  return manikde!(RES, Pose2)
end


text:
sum=0
for i in 1:1000
    sum+=(i%3==0||i%5==0)*i
end
println("sum is $sum")


text:
# testset = true
using Hijack, Test

@testset "testset" begin
    @test true
    push!(Hijack.RUN, 1)
    include("included_testset.jl")
end

max_issues_repo_name: rfourquet/InlineTest

max_forks_repo_name: rfourquet/InlineTest


text:
module Examples

using ..AlphaZero

include("../games/tictactoe/main.jl")
export Tictactoe

include("../games/toetactic/main.jl")
export Toetactic

include("../games/connect-four/main.jl")
export ConnectFour

include("../games/grid-world/main.jl")
export GridWorld

include("../games/mancala/main.jl")
export Mancala

include("../games/gomoku/main.jl")
export Gomoku

const games = Dict(
    "grid-world" => GridWorld.GameSpec(),
    "tictactoe" => Tictactoe.GameSpec(),
    "toetactic" => Toetactic.GameSpec(),
    "connect-four" => ConnectFour.GameSpec(),
    "mancala" => Mancala.GameSpec(),
    "gomoku" => Gomoku.GameSpec(),
)
# "ospiel_ttt" => OSpielTictactoe.GameSpec()
# ospiel_ttt is added from openspiel_example.jl when OpenSpiel.jl is imported

const experiments = Dict(
    "grid-world" => GridWorld.Training.experiment,
    "tictactoe" => Tictactoe.Training.experiment,
    "toetactic" => Toetactic.Training.experiment,
    "connect-four" => ConnectFour.Training.experiment,
    "mancala" => Mancala.Training.experiment,
    "gomoku" => Gomoku.Training.experiment,
)
# "ospiel_ttt" => OSpielTictactoe.Training.experiment

end


text: >(n,i=0)=n<2 ? i : n%2<1 ? n/2>i+1 : 3n+1>i+1

max_stars_repo_name: pgrepds/Codegolf

max_issues_repo_name: pgrepds/Codegolf

max_forks_repo_name: pgrepds/Codegolf


text:
import Base.Test.@test
using NetCDF
include("intermediate.jl")
include("high.jl")
include("array-like.jl")


text:
# module escape_dimer_one_hit
# export escape_exit_num_dimer
using DifferentialEquations
using Roots
using Printf
# using MATLAB
# using ProgressMeter

include("leap_frog_definitions.jl")

max_hit_q1_back = 1
barrier = 5

H=.25
condition_max_hits_back_one(u, t, integrator) = u[3] == max_hit_q1_back
affect_stop!(integrator) = terminate!(integrator)
function condition_hits_PSS_q1(u, t, integrator) # Event when event_f(u,t) == 0
    u[1]
end
function affect_update_iterator_q1!(integrator)
    q2 = integrator.u[2]
    p2 = integrator.u[5]
    p1 = integrator.u[4]
    P1 = P1_poly(q2, p2, H)
    if !isempty(P1) && abs(p1 - P1) < 1e-5
        integrator.u[3] = integrator.u[3] + 1
    end
end

callback_max_hits_back_one = DiscreteCallback(condition_max_hits_back_one, affect_stop!)
callback_hits_PSS_q1_back_one = ContinuousCallback(condition_hits_PSS_q1, affect_update_iterator_q1!, nothing, rootfind = true)
cb_back_one = CallbackSet(callback_hits_PSS_q1_back_one, callback_max_hits_back_one)
function back_one(Q2, P2, t_end, H)

    P1 = P1_poly(Q2, P2, H)

    if ~isempty(P1)
        q0, p0 = [zeros(3) for i in 1:2]
        q0[1] = 0
        q0[2] = Q2
        q0[3] = 0
        p0[1] = P1[1]
        p0[2] = P2
        p0[3] = 0
       #constructor for ODE
        prob = HamiltonianProblem{true}(Hamiltonian_Dimer_Backwards, q0, p0, (0., t_end));
        # solve ode , save_everystep=false is important to prevent sol to include all points, not just event points
        # sol=solve(prob, Vern9(), maxiters = 1e20, reltol = 1e-13, abstol = 1e-15 ,callback=cb_back_one,save_start=false,save_end=true,save_everystep=false)
        sol = solve(prob, Vern9(), maxiters = 1e20, reltol = 1e-13, abstol = 1e-15, callback = cb_back_one, save_start = false, save_end = true, save_everystep = false)
        # sol = solve(prob, RK4(), maxiters = 1e20, callback = cb_back_one, save_start = false, save_end = true, save_everystep = false)

        # sol=solve(prob, RK4(),maxiters=1e20,reltol=1e-8,abstol=1e-10,callback=cb,save_start=true,save_end=true,save_everystep=false)
        # sol=solve(prob, Vern9(),maxiters=1e20, reltol=1e-13,abstol=1e-15,callback=cb,save_start=true,save_end=false,save_everystep=false)

        uf = zeros(6)
        uf[1] = sol[1,end]
        uf[2] = sol[2,end]
        uf[3] = sol[3,end]
        uf[4] = sol[4,end]
        uf[5] = sol[5,end]
        uf[6] = sol[6,end]
        q1 = uf[1]
        q2 = uf[2]
        p1 = uf[4]
        p2 = uf[5]
        q = [q1, q2]
        p = [p1, p2]
        diff = abs(Hamiltonian_Dimer(q, p, 0) - H)

        if uf[3] == 1 && diff < 1e-5
            return q1, q2, p1, p2

        end
        # if dH<1e-5
        #     if is_it(uf,H,tol_dist)
        #         return sol.u[end][5]
        #     else
        #         return -1
        #     end
        # end
    end
end


text:
#=
---------------------------------------------------------------------
This function advances the solution in time using the LSRK45 Method by Carpenter-Kennedy 1994.

Written by F.X. Giraldo on July 12, 2021
           Department of Applied Mathematics
           Naval Postgraduate School
           Monterey; CA 93943-5216
---------------------------------------------------------------------
=#

include("create_rhs.jl")
include("plot_solution.jl")

function ti_LSRK!(q0,u,coord,M,Dwe,intma,periodicity,time,ntime,dt,space_method,plot_movie,DFloat)

    #Initialize RK coefficients
    RKA = (DFloat(0),
           DFloat(-567301805773)  / DFloat(1357537059087),
           DFloat(-2404267990393) / DFloat(2016746695238),
           DFloat(-3550918686646) / DFloat(2091501179385),
           DFloat(-1275806237668) / DFloat(842570457699 ))

    RKB = (DFloat(1432997174477) / DFloat(9575080441755 ),
           DFloat(5161836677717) / DFloat(13612068292357),
           DFloat(1720146321549) / DFloat(2090206949498 ),
           DFloat(3134564353537) / DFloat(4481467310338 ),
           DFloat(2277821191437) / DFloat(14882151754819))

    RKC = (DFloat(0),
           DFloat(1432997174477) / DFloat(9575080441755),
           DFloat(2526269341429) / DFloat(6820363962896),
           DFloat(2006345519317) / DFloat(3224310063776),
           DFloat(2802321613138) / DFloat(2924317926251))

    Npoin=length(q0)
    dq=zeros(DFloat,Npoin)
    qp=copy(q0)
    stages=length(RKA)

    #Time Integration
    for itime=1:ntime
        time=time + dt
        #RK Stages
        for s = 1:stages
            #Create RHS vector
            #= -----------Students add your CREATE_RHS routine here-----------=#
            #= -----------Students add your CREATE_RHS routine here-----------=#
            R=create_rhs(qp,u,M,Dwe,intma,periodicity,DFloat)
            #= -----------Students add your CREATE_RHS routine here-----------=#
            #= -----------Students add your CREATE_RHS routine here-----------=#

            #Solve System
            @inbounds for I=1:Npoin
                dq[I] = RKA[s]*dq[I] + dt*R[I]
                qp[I] = qp[I] + RKB[s]*dq[I]
            end
            if periodicity[Npoin] == periodicity[1]
                qp[Npoin]=qp[1] #periodicity
            end
        end #s

        #Update Q
        q0 .= qp

        #Plot Solution
        if (plot_movie)
            plot_solution(q0,coord,space_method,time)
        end

         #Print to Screen
         if ( mod(itime,100) == 0 )
            println(" itime = ",itime," time = ",time," qmax = ",maximum(q0)," qmin = ",minimum(q0))
        end

    end #itime

    return(q0)

end


text:
using FdeSolver
using Test
using SpecialFunctions

@testset "FdeSolver.jl" begin

    tSpan = [0, 5]
    y0 = [1, 0.5, 0.3]
    β = [0.5, 0.2, 0.6]

    function F(t, n, β, y)

        F1 = 1 / sqrt(pi) * (((y[n, 2] - 0.5) * (y[n, 3] - 0.3))^(1 / 6) + t[n]^(1 / 2))
        F2 = gamma(2.2) * (y[n, 1] - 1)
        F3 = gamma(2.8) / gamma(2.2) * (y[n, 2] - 0.5)

        return [F1, F2, F3]

    end

    t, Yapp = FDEsolver(F, tSpan, y0, β, nothing)

    @test @isdefined(t)
    @test @isdefined(Yapp)

end


text:
# Metida

using  Test, CSV, DataFrames, StatsModels, StatsBase, LinearAlgebra, CategoricalArrays

path    = dirname(@__FILE__)
include("testdata.jl")

@testset "  Publick API basic tests                                  " begin
    io = IOBuffer();
    transform!(df0, :formulation => categorical, renamecols=false)
    # Basic, no block
    df0.nosubj = ones(size(df0, 1))
    lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0;
    random = Metida.VarEffect(Metida.@covstr(formulation|nosubj), Metida.DIAG),
    )

    # Basic show (before fitting)
    Base.show(io, lmm)
    Metida.fit!(lmm)
    Base.show(io, lmm)
    Base.show(io, lmm.data)
    Base.show(io, lmm.result)
    Base.show(io, lmm.covstr)
    Base.show(io, lmm.covstr.repeated.covtype)
    Base.show(io, Metida.getlog(lmm))
    #
    @test Metida.m2logreml(lmm) ≈ 25.129480634331067 atol=1E-6
    #Verbose
    Metida.fit!(lmm; verbose = 2, io = io)
    @test Metida.m2logreml(lmm) ≈ 25.129480634331067 atol=1E-6
    #Missing
    lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0m;
    random = Metida.VarEffect(Metida.@covstr(formulation|subject), Metida.DIAG),
    )
    Metida.fit!(lmm)
    @test Metida.m2logreml(lmm) ≈ 16.241112644506067 atol=1E-6
    #Basic, Subject block
    lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0;
    random = Metida.VarEffect(Metida.@covstr(formulation|subject), Metida.DIAG),
    )
    Metida.fit!(lmm; aifirst = true)
    @test Metida.m2logreml(lmm) ≈ 16.241112644506067 atol=1E-6

    t3table = Metida.typeiii(lmm;  ddf = :contain) # NOT VALIDATED
    t3table = Metida.typeiii(lmm;  ddf = :residual)
    t3table = Metida.typeiii(lmm)
    Base.show(io, t3table)

    ############################################################################

    ############################################################################
    # API test
    ############################################################################
    l = [0 0 1 0 0 0; 0 0 0 1 0 0; 0 0 0 0 1 0]
    @test Metida.logreml(lmm)   ≈ -8.120556322253035 atol=1E-6
    @test isfitted(lmm) == true
    @test bic(lmm)              ≈ 24.558878811225412 atol=1E-6
    @test aic(lmm)              ≈ 22.241112644506067 atol=1E-6
    @test aicc(lmm)             ≈ 24.241112644506067 atol=1E-6
    @test Metida.caic(lmm)      ≈ 27.558878811225412 atol=1E-6
    @test dof_residual(lmm) == 14

    @test Metida.dof_satter(lmm, 6)   ≈ 5.81896814947982 atol=1E-2
    @test Metida.dof_satter(lmm)[end] ≈ 5.81896814947982 atol=1E-2
    @test Metida.dof_satter(lmm, [0 0 0 0 0 1]) ≈ 5.81896814947982 atol=1E-2
    @test Metida.dof_satter(lmm, l) ≈ 7.575447546211385 atol=1E-2
    @test Metida.fvalue(lmm, l) ≈  0.202727915619993 atol=1E-2
    @test Metida.dof_satter(lmm, Metida.lcontrast(lmm,3)) ≈ 7.575447546211385 atol=1E-2
    @test nobs(lmm) == 20
    @test Metida.thetalength(lmm) == 3
    @test Metida.rankx(lmm) == 6
    @test sum(Metida.gmatrix(lmm, 1)) ≈ 0.3350555603325126 atol=1E-6
    @test sum(Metida.rmatrix(lmm, 1)) ≈ 0.13699999248885292 atol=1E-6
    @test sum(Metida.vmatrix(lmm, 1)) ≈ 1.4772222338189034 atol=1E-6
    @test dof(lmm) == 7
    @test vcov(lmm)[1,1]              ≈ 0.11203611149231425 atol=1E-6
    @test stderror(lmm)[1]            ≈ 0.33471795812641164 atol=1E-6
    @test length(modelmatrix(lmm)) == 120
    @test isa(response(lmm), Vector)
    @test sum(Metida.hessian(lmm))    ≈ 1118.160713481362 atol=1E-2
    @test Metida.nblocks(lmm) == 5
    @test length(coefnames(lmm)) == 6
    @test Metida.confint(lmm)[end][1] ≈ -0.7630380758015894 atol=1E-4
    @test Metida.confint(lmm; ddf = :residual)[end][1] ≈ -0.6740837049617738 atol=1E-4

    Metida.confint(lmm; ddf = :contain)[end][1] #NOT VALIDATED
    @test size(crossmodelmatrix(lmm), 1) == 6
    @test t3table.pval[4]          ≈ 0.7852154468081014 atol=1E-6
    ct = Metida.contrast(lmm, [0 0 1 0 0 0; 0 0 0 1 0 0; 0 0 0 0 1 0])
    @test t3table.pval[3] ≈ ct.pval[1]
    est = Metida.estimate(lmm, [0,0,0,0,0,1]; level = 0.9)
    @test_nowarn Base.show(io, est)
    ############################################################################
    # AI like algo
    Metida.fit!(lmm; aifirst = true, init = Metida.theta(lmm))
    @test Metida.m2logreml(lmm) ≈ 16.241112644506067 atol=1E-6
    # Score
    Metida.fit!(lmm; aifirst = :score)
    @test Metida.m2logreml(lmm) ≈ 16.241112644506067 atol=1E-6
    # AI
    Metida.fit!(lmm; aifirst = :ai)
    @test Metida.m2logreml(lmm) ≈ 16.241112644506067 atol=1E-6
    #Set user coding
    lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0;
    random = Metida.VarEffect(Metida.@covstr(1 + formulation|subject), Metida.CSH; coding = Dict(:formulation => StatsModels.DummyCoding())),
    )
    # Test varlink/rholinkf
    Metida.fit!(lmm; rholinkf = :sqsigm)
    @test Metida.dof_satter(lmm, [0, 0, 0, 0, 0, 1]) ≈ 6.043195705464293 atol=1E-2
    @test Metida.m2logreml(lmm) ≈ 10.314822559210157 atol=1E-6
    @test_nowarn Metida.fit!(lmm; varlinkf = :sq)
    # Repeated effect only
    lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0;
    repeated = Metida.VarEffect(Metida.@covstr(formulation|nosubj)),
    )
    Metida.fit!(lmm)
    @test Metida.m2logreml(lmm) ≈ 25.129480634331063 atol=1E-6

    #BE like
    lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0;
    random = Metida.VarEffect(Metida.@covstr(formulation|subject), Metida.CSH),
    repeated = Metida.VarEffect(Metida.@covstr(formulation|subject), Metida.DIAG),
    )
    Metida.fit!(lmm; aifirst = :score)
    @test Metida.m2logreml(lmm) ≈ 10.065238626765524 atol=1E-6
    #incomplete
    lmm = Metida.LMM(@formula(var~sequence+period+formulation), df1;
    random = Metida.VarEffect(Metida.@covstr(formulation|subject), Metida.CSH),
    repeated = Metida.VarEffect(Metida.@covstr(formulation|subject), Metida.DIAG),
    )
    Metida.fit!(lmm; hes = false)
    @test Metida.m2logreml(lmm) ≈ 14.819463206995163 atol=1E-6
    @test Metida.dof_satter(lmm, 6)   ≈ 3.981102548214154 atol=1E-2

    lmm = Metida.LMM(@formula(var~period*formulation), df0;
    random = Metida.VarEffect(Metida.@covstr(formulation+sequence|nosubj), Metida.SI),
    )
    Metida.fit!(lmm)
    @test Metida.m2logreml(lmm, [0.222283, 0.444566]) ≈ Metida.m2logreml(lmm) atol=1E-6
    # EXPERIMENTAL
    @test Metida.dof_contain(lmm, 1) == 12
    @test Metida.dof_contain(lmm, 5) == 8
end
################################################################################
#                                  df0
################################################################################
@testset "  Model: Only repeated, 0/DIAG                             " begin
    # nowarn
    io = IOBuffer();
    lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0;
    repeated = Metida.VarEffect(Metida.@covstr(formulation|subject), Metida.DIAG),
    )
    Metida.fit!(lmm)
    Base.show(io, lmm)
    @test Metida.m2logreml(lmm) ≈ 25.000777869122338 atol=1E-8
    lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0;
    repeated = Metida.VarEffect(Metida.@covstr(formulation|subject), Metida.DIAG)
    )
    Metida.fit!(lmm)
    @test Metida.m2logreml(lmm) ≈ 25.000777869122338 atol=1E-8
end
@testset "  Model: Only repeated, noblock, 0/CSH (rholinkf = :atan)  " begin
    # nowarn
    lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0;
    repeated = Metida.VarEffect(Metida.@covstr(period|subject), Metida.CSH),
    )
    Metida.fit!(lmm; rholinkf = :atan)
    @test Metida.m2logreml(lmm) ≈ 8.740095378772942 atol=1E-8
end
@testset "  Model: Only random, noblock, SI/SI                       " begin
    # nowarn
    lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0;
    random = Metida.VarEffect(Metida.@covstr(subject|nosubj), Metida.SI),
    )
    Metida.fit!(lmm)
    @test Metida.m2logreml(lmm) ≈ 10.862124583312674 atol=1E-8
end
@testset "  Model: Only random, INT, SI/SI                           " begin
    # nowarn
    lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0;
    random = Metida.VarEffect(Metida.@covstr(1|subject), Metida.SI),
    )
    Metida.fit!(lmm)
    @test Metida.m2logreml(lmm) ≈ 10.862124583312674 atol=1E-8
end
@testset "  Model: Noblock, equal subjects, CSH/CS                   " begin
    lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0;
    random = Metida.VarEffect(Metida.@covstr(formulation|subject), Metida.CSH),
    repeated = Metida.VarEffect(Metida.@covstr(formulation|subject), Metida.CS),
    )
    Metida.fit!(lmm)
    @test Metida.m2logreml(lmm) ≈ 10.3039977509049 atol=1E-6 #need check
end
@testset "  Model: Different subjects, INT, CSH/DIAG                 " begin
    lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0;
    random = Metida.VarEffect(Metida.@covstr(1 + formulation|subject), Metida.CSH; coding = Dict(:formulation => DummyCoding())),
    repeated = Metida.VarEffect(Metida.@covstr(formulation|subject&period), Metida.DIAG),
    )
    Metida.fit!(lmm)
    @test Metida.m2logreml(lmm) ≈ 10.06523870216023 atol=1E-4 #need check
end
@testset "  Model: CSH/DIAG (rholinkf = :psigm)                      " begin
    lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0;
    random = Metida.VarEffect(Metida.@covstr(formulation|subject), Metida.CSH),
    repeated = Metida.VarEffect(Metida.@covstr(formulation|subject), Metida.DIAG),
    )
    Metida.fit!(lmm; rholinkf = :psigm)
    @test Metida.m2logreml(lmm) ≈ 10.065239006121315 atol=1E-6
end
@testset "  Model: Custom covariance type                            " begin
    CCTG = Metida.CovarianceType(Metida.CovmatMethod((q,p) -> (q, 1), (mx, θ, p) -> Metida.gmat_csh!(mx, θ)))
    CCTR = Metida.CovarianceType(Metida.CovmatMethod((q,p) -> (q, 0), (mx, θ, zrv, p) -> Metida.rmatp_diag!(mx, θ, zrv)))
    lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0;
    random = Metida.VarEffect(Metida.@covstr(formulation|subject), CCTG),
    repeated = Metida.VarEffect(Metida.@covstr(formulation|subject), CCTR),
    )

    Metida.fit!(lmm)
    @test Metida.m2logreml(lmm) ≈ 10.065239006121315 atol=1E-6
end
################################################################################
#                                  ftdf / 1fptime.csv
################################################################################
@testset "  Model: Categorical * Continuous predictor, CSH/SI        " begin
    # nowarn
    lmm = Metida.LMM(@formula(response ~1 + factor*time), ftdf;
    random = Metida.VarEffect(Metida.@covstr(1 + time|subject&factor), Metida.CSH),
    )
    Metida.fit!(lmm)
    @test Metida.m2logreml(lmm) ≈ 1300.1807598168923 atol=1E-6
    @test coef(lmm) ≈ [22.13309710783416, 2.000486297455917, 1.1185284725578566, 0.4049714576872601] atol=1E-6
    @test Metida.dof_satter(lmm, [0, 0, 0, 1]) ≈ 37.999999999991786 atol=1E-2
    #Metida.typeiii(lmm)
end

@testset "  Model: Function terms, CSH/SI                            " begin
    ftdf.expresp = exp.(ftdf.response)
    ftdf.exptime = exp.(ftdf.time)
    lmm = Metida.LMM(@formula(log(expresp) ~ 1 + factor*log(exptime)), ftdf;
    random = Metida.VarEffect(Metida.@covstr(1 + log(exptime)|subject&factor), Metida.CSH),
    )
    Metida.fit!(lmm)
    @test Metida.m2logreml(lmm) ≈ 1300.1807598168923 atol=1E-6
end
################################################################################
#                                  ftdf2 / 1freparma.csv
################################################################################
@testset "  Model: Categorical * Continuous predictor, 0/ARMA        " begin
    # nowarn
    # SPSS 715.452856
    io = IOBuffer();
    lmm = Metida.LMM(@formula(response ~ 1 + factor*time), ftdf2;
    repeated = Metida.VarEffect(Metida.@covstr(time|subject&factor), Metida.ARMA),
    )
    Metida.fit!(lmm)
    println(io, lmm.log)
    @test Metida.m2logreml(lmm) ≈ 715.4528559688382 atol = 1E-6
end
@testset "  Model: Categorical * Continuous predictor, DIAG/AR       " begin
    # nowarn
    lmm = Metida.LMM(@formula(response ~ 1 + factor*time), ftdf2;
    random = Metida.VarEffect(Metida.@covstr(factor|subject&factor), Metida.DIAG),
    repeated = Metida.VarEffect(Metida.@covstr(1|subject&factor), Metida.AR),
    )
    Metida.fit!(lmm)
    @test Metida.m2logreml(lmm) ≈ 710.0962305879676 atol=1E-6
end
@testset "  Model: Categorical * Continuous predictor, 0/ARH         " begin
    # nowarn
    lmm = Metida.LMM(@formula(response ~ 1 + factor*time), ftdf2;
    repeated = Metida.VarEffect(Metida.@covstr(factor|subject&factor), Metida.ARH),
    )
    Metida.fit!(lmm)
    @test Metida.m2logreml(lmm) ≈ 731.7794071577566 atol=1E-6
end
################################################################################
#                                  ftdf3 / 2f2rand.csv
################################################################################
@testset "  Model: CS, CS/SI                                         " begin
    lmm = Metida.LMM(@formula(response ~ 1 + factor), ftdf3; contrasts=Dict(:factor => DummyCoding(; base=1.0)),
    random = [Metida.VarEffect(Metida.@covstr(r1|subject), Metida.CS), Metida.VarEffect(Metida.@covstr(r2|subject), Metida.CS)],
    )
    Metida.fit!(lmm)
    @test Metida.m2logreml(lmm)  ≈ 710.4250214813896 atol=1E-8
end
@testset "  Model: AR/SI                                             " begin
    # SPSS 698.879
    # nowarn
    lmm = Metida.LMM(@formula(response ~ 1 + factor), ftdf3; contrasts=Dict(:factor => DummyCoding(; base=1.0)),
    random = Metida.VarEffect(Metida.@covstr(r1|subject), Metida.AR),
    repeated = Metida.VarEffect(Metida.@covstr(p|subject), Metida.DIAG),
    )
    Metida.fit!(lmm)
    @test Metida.m2logreml(lmm)  ≈ 698.8792511057682 atol=1E-8
end

@testset "  Model: ARMA/SI                                           " begin
    io = IOBuffer();
    lmm = Metida.LMM(@formula(response ~ 1 + factor), ftdf3; contrasts=Dict(:factor => DummyCoding(; base=1.0)),
    random = Metida.VarEffect(Metida.@covstr(p|r1&r2), Metida.ARMA),
    )
    Metida.fit!(lmm; verbose = 3, io = io)
    #[1.2964e-5, 0.0299594, 0.0699728, 3.69557]
    println(io, lmm.log)
    @test Metida.m2logreml(lmm)  ≈ 913.9176298311813 atol=1E-8
end

@testset "  Model: ARH/SI (subjects with &)                          " begin
    # SPSS 707.377
    # nowarn
    lmm = Metida.LMM(@formula(response ~ 1 + factor), ftdf3; contrasts=Dict(:factor => DummyCoding(; base=1.0)),
    random = Metida.VarEffect(Metida.@covstr(r1|s2&factor), Metida.ARH),
    )
    Metida.fit!(lmm)
    @test Metida.m2logreml(lmm)  ≈ 707.3765873864152 atol=1E-8
    @test Metida.dof_satter(lmm, [0, 1]) ≈ 23.111983305626193 atol=1E-2
end
@testset "  Model: INT, *, DIAG/SI                                   " begin
    lmm = Metida.LMM(@formula(response ~ 1 + factor), ftdf3;
    random = Metida.VarEffect(Metida.@covstr(1 + r2 * r1|subject), Metida.DIAG; coding=Dict(:r1 => DummyCoding(), :r2 => DummyCoding()))
    )
    Metida.fit!(lmm)
    @test Metida.m2logreml(lmm)  ≈ 713.0655862252027 atol=1E-8
end
@testset "  Model: &, DIAG/SI                                        " begin
    lmm = Metida.LMM(@formula(response ~ 1 + factor), ftdf3;
    random = Metida.VarEffect(Metida.@covstr(r1&r2|subject), Metida.DIAG),
    )
    Metida.fit!(lmm)
    @test Metida.m2logreml(lmm)  ≈ 719.9413776641368 atol=1E-8
end
@testset "  Model: INT, +,  TOEPHP(3)/SI                             " begin
    lmm = Metida.LMM(@formula(response ~ 1 + factor), ftdf3;
    random = Metida.VarEffect(Metida.@covstr(1 + r1 + r2|subject), Metida.TOEPHP(3); coding = Dict(:r1 => DummyCoding(), :r2 => DummyCoding())),
    )
    Metida.fit!(lmm)
    @test Metida.m2logreml(lmm)  ≈ 705.9946274598822 atol=1E-8
end
@testset "  Model: TOEP/SI                                           " begin
    # SPSS 710.200
    # nowarn
    io = IOBuffer();
    lmm = Metida.LMM(@formula(response ~ 1 + factor), ftdf3;
    random = Metida.VarEffect(Metida.@covstr(r1|subject), Metida.TOEP),
    )
    Metida.fit!(lmm)
    Base.show(io, lmm)
    @test Metida.m2logreml(lmm)  ≈ 710.1998669150806 atol=1E-8
end
@testset "  Model: TOEPP(2)/SI                                       " begin
    # nowarn
    io = IOBuffer();
    lmm = Metida.LMM(@formula(response ~ 1 + factor), ftdf3;
    random = Metida.VarEffect(Metida.@covstr(r1|subject), Metida.TOEPP(2)),
    )
    Metida.fit!(lmm)
    Base.show(io, lmm)
    @test Metida.m2logreml(lmm)  ≈ 715.2410264030134 atol=1E-8
end
@testset "  Model: DIAG/TOEPP(3)                                     " begin
    # nowarn
    io = IOBuffer();
    lmm = Metida.LMM(@formula(response ~ 1 + factor), ftdf3;
    random = Metida.VarEffect(Metida.@covstr(r2|subject), Metida.DIAG),
    repeated = Metida.VarEffect(Metida.@covstr(p|subject), Metida.TOEPP(3)),
    )
    Metida.fit!(lmm)
    Base.show(io, lmm)
    @test Metida.m2logreml(lmm)  ≈ 773.9575538254085 atol=1E-8
end
@testset "  Model: TOEPH/SI                                          " begin
    # nowarn
    io = IOBuffer();
    lmm = Metida.LMM(@formula(response ~ 1 + factor), ftdf3;
    random = Metida.VarEffect(Metida.@covstr(r1|subject), Metida.TOEPH),
    )
    Metida.fit!(lmm)
    Base.show(io, lmm)
    @test Metida.m2logreml(lmm)  ≈ 705.7916833009426 atol=1E-8
end
@testset "  Model: SI/TOEPHP(3)                                      " begin
    io = IOBuffer();
    lmm = Metida.LMM(@formula(response ~ 1 + factor), ftdf3;
    random = Metida.VarEffect(Metida.@covstr(r1|subject), Metida.SI),
    repeated = Metida.VarEffect(Metida.@covstr(r1&r2|subject), Metida.TOEPHP(3)),
    )
    Metida.fit!(lmm; optmethod = Metida.LBFGS_OM)
    Base.show(io, lmm)
    @test Metida.m2logreml(lmm)  ≈ 713.5850978377632 atol=1E-8
end
@testset "  Model: BE RDS 1, FDA model                               "  begin
    dfrds        = CSV.File(joinpath(path, "csv", "berds", "rds1.csv"), types = Dict(:PK => Float64, :subject => String, :period => String, :sequence => String, :treatment => String )) |> DataFrame
    dropmissing!(dfrds)
    dfrds.lnpk = log.(dfrds.PK)
    lmm = Metida.LMM(@formula(lnpk~sequence+period+treatment), dfrds;
    random = Metida.VarEffect(Metida.@covstr(treatment|subject), Metida.CSH),
    repeated = Metida.VarEffect(Metida.@covstr(treatment|subject), Metida.DIAG),
    )
    Metida.fit!(lmm)
    @test collect(Metida.confint(lmm)[6]) ≈  [0.05379033790060175, 0.23713821749515449] atol=1E-8
    anovatable = Metida.typeiii(lmm)
    @test anovatable.pval ≈ [3.087934998046721e-63, 0.9176105002577626, 0.6522549061162943, 0.002010933915677479] atol=1E-4

    est = Metida.estimate(lmm, [0,0,0,0,0,1]; level = 0.9)
    @test est.t[1] ≈ 3.12818 atol=1E-4
    @test est.pval[1] ≈ 0.0020 atol=1E-4
    @test est.cil[1] ≈ 0.06863 atol=1E-4
    @test est.ciu[1] ≈ 0.2223 atol=1E-4

    lmm = Metida.LMM(@formula(lnpk~0+sequence+period+treatment), dfrds;
    random = Metida.VarEffect(Metida.@covstr(treatment|subject), Metida.CSH),
    repeated = Metida.VarEffect(Metida.@covstr(treatment|subject), Metida.DIAG),
    )
    Metida.fit!(lmm)
    anovatable = Metida.typeiii(lmm)
    @test anovatable.pval ≈ [0.9176105002855397, 0.6522549061174356, 0.0020109339157131302] atol=1E-4
end
################################################################################
#                                  Errors
################################################################################
@testset "  Errors test                                              " begin
    lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0;
    random = Metida.VarEffect(Metida.@covstr(formulation|nosubj), Metida.DIAG),
    )
    @test_throws ErrorException Metida.fit!(lmm; init = [1.0])
    @test_throws ErrorException Metida.hessian(lmm)
    @test_throws ErrorException Metida.dof_satter(lmm)
    @test_throws ErrorException Metida.confint(lmm)

    @test_throws ErrorException  Metida.LMM(@formula(var~sequence+period+formulation), df0;)

    @test_throws ErrorException  Metida.LMM(@formula(var~sequence+period+formulation), df0;
    random = [Metida.VarEffect(Metida.@covstr(formulation|nosubj), Metida.DIAG), Metida.VarEffect(Metida.@covstr(formulation|nosubj), Metida.RZero())]
    )
    # Error messages
    io = IOBuffer();
    lmm = Metida.LMM(@formula(response ~ 1 + factor*time), ftdf2;
    random = Metida.VarEffect(Metida.@covstr(factor|subject&factor), Metida.DIAG),
    repeated = Metida.VarEffect(Metida.@covstr(1|subject&factor), Metida.ARMA),
    )
    Metida.fit!(lmm)
    println(io, lmm.log)
end
################################################################################
#                                  Sweep test
################################################################################
@testset "  Sweep operator test                                      " begin
    A =
[1.0  2  2  4  1
 2  2  3  3  5
 2  3  3  4  2
 4  3  4  4  5
 1  5  2  5  5]
    iA =  inv(A[1:4, 1:4])
    iAs = Symmetric(-Metida.sweep!(copy(A), 1:4; syrkblas = true)[1][1:4, 1:4])
    B = copy(A)
    for i = 1:4
        Metida.sweep!(B, i; syrkblas = false)
    end
    iAss = Symmetric(-B[1:4, 1:4])
    akk = zeros(5)
    iAb = Symmetric(-Metida.sweepb!(view(akk, 1:5), copy(A), 1:4)[1][1:4, 1:4])
    @test iA  ≈ iAs  atol=1E-6
    @test iA  ≈ iAss atol=1E-6
    @test iAs ≈ iAb  atol=1E-6
end

@testset "  Experimental                                             " begin
    lmm = Metida.LMM(@formula(response ~ 1), ftdf;
    repeated = Metida.VarEffect(Metida.@covstr(response+time|subject), Metida.SPEXP),
    )
    Metida.fit!(lmm)
    @test Metida.m2logreml(lmm) ≈ 1528.7150702624508 atol=1E-6

    @test_nowarn Metida.fit!(lmm; varlinkf = :identity)
end


text:
function ALsearch(userinput; verbose=true)
    try
      checker = replace(userinput[7:end], r"^[\\:\\-\\ ]+"=>"")
      search(checker, replset, verbose = verbose)
      # push!(logicset, replset)
  catch er
      replthrow("Warning! Search Failed: $er")
    end
end


text:
"""
    compound_percentage_occurrence(compound_name::String, impurity_profile_csv, minimum_value::Number; show_not_containing::Bool=true, date_range::Date)
Searches through impurity profile and determines percentage occurrence of compound
Can also search for a combination of compounds if compound_name and minimum_value are given as vector
Give date_range (optional) to only display between given dates format: [yyyy-mm-dd, yyyy-mm-dd]
"""
function compound_percentage_occurrence(compound_name::String, impurity_profile_csv, minimum_value::Number)
    # Read impurity profile and filter on dates
    impurity_profile = CSV.read(impurity_profile_csv, DataFrame)
    date_range = sort(Date.(date_range))

    # Count occurrences
    occurrence_count = count(x -> x >= minimum_value, impurity_profile[!, compound_name])

    # Determine percent occurrence
    total_samples = nrow(impurity_profile)
    percent_occurrence = round(occurrence_count/nrow(impurity_profile) * 100, digits=2)
    println("Percent occurrence: $percent_occurrence, Total occurrence: $occurrence_count, Total number of samples: $total_samples")
    return percent_occurrence, occurrence_count, total_samples
end

# Date range given
function compound_percentage_occurrence(compound_name::String, impurity_profile_csv, minimum_value::Number; date_range::Date)
    # Read impurity profile and filter on dates
    impurity_profile = CSV.read(impurity_profile_csv, DataFrame)
    date_range = sort(Date.(date_range))

    # Filter date
    filter!(row -> Date(row."DateTime (y-m-d)"[1:10]) >= date_range[1] && Date(row."DateTime (y-m-d)"[1:10]) <= date_range[2] , impurity_profile)

    # Count occurrences
    occurrence_count = count(x -> x >= minimum_value, impurity_profile[!, compound_name])

    # Determine percent occurrence
    total_samples = nrow(impurity_profile)
    percent_occurrence = round(occurrence_count/nrow(impurity_profile) * 100, digits=2)
    println("Percent occurrence: $percent_occurrence, Total occurrence: $occurrence_count, Total number of samples: $total_samples")
    return percent_occurrence, occurrence_count, total_samples
end

# Searches for combination of compounds
function compound_percentage_occurrence(compound_name::Vector{String}, impurity_profile_csv, date_range, minimum_value::Vector{Number})
    if length(compound_name) != length(minimum_value)
        @error "compound_name and minimum_value not same length"
    end

    # Read impurity profile and filter on dates
    impurity_profile = CSV.read(impurity_profile_csv, DataFrame)
    date_range = sort(Date.(date_range))
    filter!(row -> Date(row."DateTime (y-m-d)"[1:10]) >= date_range[1] && Date(row."DateTime (y-m-d)"[1:10]) <= date_range[2] , impurity_profile)

    # Count occurrences
    occurrence_count = 0
    for row in eachrow(impurity_profile)
        meets_req = true
        for (compound, min_val) in zip(compound_name, minimum_value)
            if row[compound] < min_val
                meets_req = false
                break
            end
        end
        if meets_req
            occurrence_count += 1
        end
    end

    # Determine percent occurrence
    total_samples = nrow(impurity_profile)
    percent_occurrence = round(occurrence_count/nrow(impurity_profile) * 100, digits=2)
    println("Percent occurrence: $percent_occurrence, Total occurrence: $occurrence_count, Total number of samples: $total_samples")
    return percent_occurrence, occurrence_count, total_samples
end

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

"""
    compound_trends(impurity_profile_csv, compound_group, month_range, minimum_value)
Counts per month how many times each compound is above the minimum value
Outputs this to monthly_data.csv in input folder
show_not_containing defaults to true, displays samples not containing any of the given compounds
"""
function compound_trends(impurity_profile_csv, compound_group, minimum_value::Number; show_not_containing::Bool=true)

    # Read impurity profile and filter out control samples
    impurity_profile = CSV.read(impurity_profile_csv, DataFrame)
    if "item_number" in names(impurity_profile)
        filter!(x -> ismissing(x["item_number"]) || !contains(x["item_number"], "CTRL"), impurity_profile)
    end

    # Remove samples without date
    filter!(x -> !ismissing(x["DateTime (y-m-d)"]), impurity_profile)

    # Group impurity profiles by month
    insertcols!(impurity_profile,1, :yearmonth => yearmonth.(map(x -> DateTime(x[1:end-1]), impurity_profile[!, "DateTime (y-m-d)"])))
    impurity_profiles_grouped = groupby(impurity_profile, :yearmonth, sort=true)

    # Read settings.json and retrieve compound names from group
    json_string = read(joinpath(SETTINGS_JSON_LOCATION, "settings.json"), String)
    settings_json = JSON3.read(json_string)
    compound_names = settings_json[:compound_groups][compound_group]

    # Create DataFrame for output
    monthly_data = DataFrame()
    insertcols!(monthly_data, :year => Int[])
    insertcols!(monthly_data, :month => Int[])
    for compound_name in compound_names
		insertcols!(monthly_data, Symbol(compound_name) => Int32[])
	end
    if show_not_containing
        insertcols!(monthly_data, :no_compounds => Int[])
    end
    insertcols!(monthly_data, :total_samples => Int[])

    # Count times compounds above minimum value for each month
    for group in impurity_profiles_grouped
        total_samples = size(group, 1)
        push!(monthly_data, append!([], group[1, :yearmonth], zeros(size(monthly_data, 2) - 3), total_samples))
        for compound_name in compound_names
            monthly_data[end, compound_name] = count(i -> i >= minimum_value, group[!, compound_name])
        end

        # Count how often no compound is present
        if show_not_containing
            data = select(group, compound_names)
            i = 0
            for row in eachrow(data)
                if any(x -> x > minimum_value, row)
                    i += 1
                end
            end
            monthly_data[end, :no_compounds] = i
        end
    end

    CSV.write(joinpath(dirname(impurity_profile_csv), "monthly_data.csv"), monthly_data)

end

"""
    find_similar_samples(sample, impurity_profile_csv, compound_group, max_deviation_fraction, date_range=0)
Looks for similar samples in impurity profile, samples are determined as similar if for each compound:
{difference values / mean values <= max_deviation_fraction}

sample can be input as DataFrameRow or as index (Int) within the impurity profile
"""
function find_similar_samples(sample::DataFrameRow, impurity_profile_csv, compound_group, max_deviation_fraction, date_range=0)

    # Read impurity profile and filter dates if given
    impurity_profile = CSV.read(impurity_profile_csv, DataFrame)
    if date_range != 0
        date_range = sort(Date.(date_range))
        filter!(row -> Date(row."DateTime (y-m-d)"[1:10]) >= date_range[1] && Date(row."DateTime (y-m-d)"[1:10]) <= date_range[2] , impurity_profile)
    end

    # Read settings.json and retrieve compound names from group
    json_string = read(joinpath(SETTINGS_JSON_LOCATION, "settings.json"), String)
	settings_json = JSON3.read(json_string)
    compound_names = settings_json[:compound_groups][compound_group]

    # filter columns of impurity profile and sample based on compound group
    impurity_profile_filtered = select(impurity_profile, compound_names)
    select!(DataFrame(sample), compound_names)[1, :]

    # Determine similar samples
    similar_samples = DataFrame()
    for (i, row) in enumerate(eachrow(impurity_profile_filtered))
        meet_req = true
        # Compare values of both samples for all compounds
        for compound in compound_names
            samples_mean = mean((row[compound], sample[compound]))
            sample_deviation = abs(samples_mean - sample[compound]) / samples_mean
            if sample_deviation > max_deviation_fraction
                meet_req = false
                break
            end
        end
        if meet_req
            push!(similar_samples, impurity_profile[i, :])
        end
    end

    return similar_samples
end

# Sample input as index within impurity profile
function find_similar_samples(sample::Int, impurity_profile_csv, compound_group, max_deviation_fraction, date_range=0)

    # Read impurity profile and filter dates if given
    impurity_profile = CSV.read(impurity_profile_csv, DataFrame)
    if date_range != 0
        date_range = sort(Date.(date_range))
        filter!(row -> Date(row."DateTime (y-m-d)"[1:10]) >= date_range[1] && Date(row."DateTime (y-m-d)"[1:10]) <= date_range[2] , impurity_profile)
    end

    # Retrieve sample as DataFrameRow
    sample = impurity_profile[sample, :]

    # Read settings.json and retrieve compound names from group
    json_string = read(joinpath(SETTINGS_JSON_LOCATION, "settings.json"), String)
	settings_json = JSON3.read(json_string)
    compound_names = settings_json[:compound_groups][compound_group]

    # filter columns of impurity profile and sample based on compound group
    impurity_profile_filtered = select(impurity_profile, compound_names)
    select!(DataFrame(sample), compound_names)[1, :]

    # Determine similar samples
    similar_samples = DataFrame()
    for (i, row) in enumerate(eachrow(impurity_profile_filtered))
        meet_req = true
        # Compare values of both samples for all compounds
        for compound in compound_names
            samples_mean = mean((row[compound], sample[compound]))
            sample_deviation = abs(samples_mean - sample[compound]) / samples_mean
            if sample_deviation > max_deviation_fraction
                meet_req = false
                break
            end
        end
        if meet_req
            push!(similar_samples, impurity_profile[i, :])
        end
    end

    return similar_samples
end

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


text:
## TODO switch this to a [:,:,:,α] representation?
"""
Computes the *probability* (not charge) current, ∑ fn Im(ψn* ∇ψn)
"""
function compute_current(basis::PlaneWaveBasis, ψ, occupation)
    @assert length(basis.symmetries) == 1  # TODO lift this
    current = [zeros(eltype(basis), basis.fft_size) for α = 1:3]
    for (ik, kpt) in enumerate(basis.kpoints)
        for (n, ψnk) in enumerate(eachcol(ψ[ik]))
            ψnk_real = G_to_r(basis, kpt, ψnk)
            for α = 1:3
                dαψnk = [im * q[α] for q in Gplusk_vectors_cart(basis, kpt)] .* ψnk
                dαψnk_real = G_to_r(basis, kpt, dαψnk)
                current[α] .+= @. basis.kweights[ik] *
                                  occupation[ik][n] *
                                  imag(conj(ψnk_real) * dαψnk_real)
            end
        end
    end
    mpi_sum!.(current, Ref(basis.comm_kpts))
    current
end


text:
using StaticArrays, BenchmarkTools, ForwardDiff
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 10
show_ind_bench = true
println("Benchmarking Jacobian calculation of 3D Discrete systems for different types of e.o.m.")

u = rand(3); un = rand(3); su = SVector{3}(u); mu = MVector{3}(u)
mun = MVector{3}(un)
J = rand(3,3); MJ = MMatrix{3,3}(J); MI = MMatrix{3,3}(J);
mutable struct system1
  J::StaticArrays.MArray{Tuple{3,3},Float64,2,9}
end
sys = system1(MI)

println("\nVersion 1: StaticArray returned: eom_1(u) -> un::SVector")
@inline @inbounds function eom_1(x)
  x1 = x[1]; x2=x[2]; x3 = x[3]
  SVector{3}(
  3.8*x1*(1-x1)-0.05*(x2+0.35)*(1-2*x3),
  0.1*( (x2+0.35)*(1-2*x3)-1 )*(1-1.9*x1),
  3.78*x3*(1-x3)+0.2*x2
  )
end
println("--In-place jacobian")
@inline jac1!(J, x) = ForwardDiff.jacobian!(J, eom_1, x)
bj1 = @benchmark jac1!($sys.J, $su)
show_ind_bench && display(bj1)
println("--Return-value jacobian")
@inline jac1(x) = (j = ForwardDiff.jacobian(eom_1, x); sys.J .= j)
bjr1 = @benchmark jac1($su)
show_ind_bench && display(bjr1)

println("\nVersion 2: Base Array returned eom_2(x) -> xn::Vector ")
@inline @inbounds function eom_2(x)
  x1 = x[1]; x2=x[2]; x3 = x[3];
  xn = zeros(x)
  xn[1] = 3.8*x1*(1-x1)-0.05*(x2+0.35)*(1-2*x3)
  xn[2] = 0.1*( (x2+0.35)*(1-2*x3)-1 )*(1-1.9*x1)
  xn[3] = 3.78*x3*(1-x3)+0.2*x2
  return xn
end
cfg = ForwardDiff.JacobianConfig(eom_2, u)
println("--In-place jacobian")
@inline jac2!(J, x, cfg) = ForwardDiff.jacobian!(J, eom_2, x, cfg)
bj2 = @benchmark jac2!($sys.J, $su, $cfg)
show_ind_bench && display(bj2)
println("--Return-value jacobian")
@inline jac2(x, cfg) = (j = ForwardDiff.jacobian(eom_2, x, cfg); sys.J .= j)
bjr2 = @benchmark jac2($su, $cfg)
show_ind_bench && display(bjr2)

println("\nVersion 3: MVector returned eom_3(x) -> xn::MVector ")
@inline @inbounds function eom_3(x)
  x1 = x[1]; x2=x[2]; x3 = x[3];
  xn = MVector{3}(x)
  xn[1] = 3.8*x1*(1-x1)-0.05*(x2+0.35)*(1-2*x3)
  xn[2] = 0.1*( (x2+0.35)*(1-2*x3)-1 )*(1-1.9*x1)
  xn[3] = 3.78*x3*(1-x3)+0.2*x2
  return xn
end
cfg = ForwardDiff.JacobianConfig(eom_3, mu)
println("--In-place jacobian")
@inline jac3!(J, x, cfg) = ForwardDiff.jacobian!(J, eom_3, x, cfg)
bj3 = @benchmark jac3!($MJ, $mu, $cfg)
show_ind_bench && display(bj3)
println("--Return-value jacobian")
@inline jac3(x, cfg) = (j = ForwardDiff.jacobian(eom_3, x, cfg); sys.J .= j)
bjr3 = @benchmark jac3($mu, $cfg)
show_ind_bench && display(bjr3)

println("\nVersion 4: In place with 2 arguments: eom_towel!(xn, x) (changes xn)")
@inline @inbounds function eom_4!(xn, x)
  xn[1] = 3.8*x[1]*(1-x[1])-0.05*(x[2]+0.35)*(1-2*x[3])
  xn[2] = 0.1*( (x[2]+0.35)*(1-2*x[3])-1 )*(1-1.9*x[1])
  xn[3] = 3.78*x[3]*(1-x[3])+0.2*x[2]
end
cfg = ForwardDiff.JacobianConfig(eom_4!, un, u)
println("--In-place jacobian")
@inline jac4!(J, y, x, cfg) = ForwardDiff.jacobian!(J, eom_4!, y, x, cfg)
bj4 = @benchmark jac4!($J, $un, $u, $cfg)
show_ind_bench && display(bj4)
println("--Return-value jacobian")
@inline jac4(y, x, cfg) = (j = ForwardDiff.jacobian(eom_4!, y, x, cfg); sys.J .= j)
bjr4 = @benchmark jac4($un, $u, $cfg)
show_ind_bench && display(bjr4)

println("\nVersion 5: In place with 2 arguments & using MVector")
@inline @inbounds function eom_5!(xn::MVector, x::MVector)
  xn[1] = 3.8*x[1]*(1-x[1])-0.05*(x[2]+0.35)*(1-2*x[3])
  xn[2] = 0.1*( (x[2]+0.35)*(1-2*x[3])-1 )*(1-1.9*x[1])
  xn[3] = 3.78*x[3]*(1-x[3])+0.2*x[2]
end
cfg5 = ForwardDiff.JacobianConfig(eom_5!, mun, mu)
println("--In-place jacobian")
@inline jac5!(J, y, x, cfg) = ForwardDiff.jacobian!(J, eom_5!, y, x, cfg)
bj5 = @benchmark jac5!($MJ, $mun, $mu, $cfg5)
show_ind_bench && display(bj5)
println("--Return-value jacobian")
@inline jac5(y, x, cfg) = (j = ForwardDiff.jacobian(eom_5!, y, x, cfg); sys.J .= j)
bjr5 = @benchmark jac5($mun, $mu, $cfg5)
show_ind_bench && display(bjr5)

println("comparison of in-place jacobian call:")
for (i, b) in enumerate([bj2,bj3,bj4,bj5])
  println("v$(i+1) versus v1")
  sleep(0.1)
  display(judge(median(b), median(bj1)))
  sleep(0.1)
end

println("Conclusion about ForwardDiff jacobians:")
println("All methods for Jacobian are slower than the first (which uses SVector)")
println("However method 4 comes extremely close to it!")
println("It is actually even faster than the one with MVector")


text:
module Trap
export trap

function trap(
    left_endpt::Float64,
    right_endpt::Float64,
    trap_count::Int64,
    base_len::Float64)

    f(x) = x^2

    estimate = (f(left_endpt) + f(right_endpt)) / 2.0
    for i=1:trap_count - 1
        estimate += f(left_endpt + i * base_len) # type instability
    end
    return estimate * base_len
end
end


text:
using Econometrics, DelimitedFiles, Plots
CO2 = readdlm("../Data/CO2.data")
plot(CO2)
n = size(CO2,1)
trend = 1:n
x = [ones(n) trend]
b, junk, e, junk, junk= ols(CO2,x)
fit = x*b
plot(trend[end-36:end], e[end-36:end], legend=false)
title!("OLS residuals, last 3 years of data")
savefig("CO2Residuals.png")
gui()

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
using ArcadeLearningEnvironment

# For this example you need to obtain the Seaquest ROM file from
# https://atariage.com/system_items.html?SystemID=2600&ItemTypeID=ROM

episodes = 50

ale = ALE_new()
rom = string(ENV["ATARI_DIR"], "SEAQUEST.BIN")
println(rom)
loadROM(ale, rom)

S = Array{Int}(undef, episodes)
TR = Array{Float64}(undef, episodes)
for ei = 1:episodes
    ctr = 0.0

    fc = 0
    while game_over(ale) == false
        actions = getLegalActionSet(ale)
        ctr += act(ale, actions[rand(1:length(actions))])
        fc += 1
    end
    reset_game(ale)
    println("Game $ei ended after $fc frames with total reward $(ctr).")

    S[ei] = fc
    TR[ei] = ctr
end
ALE_del(ale)


text:
using LinearAlgebra
"""
    get_hermite_basis(n = 3, varcoeff = false)

Get Hermite basis function of order `n`.
"""
function get_hermite_basis(n = 3, varcoeff = false)
    global x
    global h
    xa = SymPy.symbols("xa")
    xb = SymPy.symbols("xb")
    if (n < 3) || (n % 2 != 1)
        println("n = ", n)
        error("The degree of Hermite polynomials should be odd and >= 3")
    end
    E = Matrix(1I, n + 1, n + 1)
    symString = ""
    for i = 1:n
        symString *= ("a" * string(i - 1, base=10) * " ")
    end
    symString *= ("a" * string(n, base=10))

    L = symbols(symString)
    basis = Array{SymPy.Sym}(undef, n+1)
    # p is a polynomial of degree n
    p = 0
    for i = 1:n+1
        p = p + L[i]*x^(i-1)
    end

    P = [] # the list containing the derivatives of p
    push!(P, p)
    for i = 1:((n+1)/2 - 1)
        push!(P, diff(p, x, convert(Int64, i)))
    end

    for i = 1:(n + 1)
        Leq = [];
        for j = 1:convert(Int64, (n+1)/2 )
            if (varcoeff)
                push!(Leq, subs(P[j], x, xa))
            else
                push!(Leq, subs(P[j], x, 0))
            end
        end
        for j = 1:convert(Int64, (n+1)/2 )
            if (varcoeff)
                push!(Leq, subs(P[j], x, xb))
            else
                push!(Leq, subs(P[j], x, h))
            end
        end
        Leq = Leq - E[1:(n+1), i]
        LL = solve(convert(Array{SymPy.Sym, 1}, Leq), collect(L[1:(n+1)]))
        LL = convert(Dict{SymPy.Sym, SymPy.Sym}, LL)
        p = 0
        for j = 1:(n+1)
            p = p + LL[L[j]] * x^(j - 1)
        end
        basis[i] = p
    end
    basis
end

"""
    get_hermite_em(p = 3, m = 0, n = 0)

Get Hermite finite elements elementary matrices.

# Arguments
  * `p`: degree of polynomials in the basis.
  * `m`: number of derivatives on the first function.
  * `n`: number of derivatives on the second function.
"""
function get_hermite_em(p = 3, m = 0, n = 0)
    global x
    global h
    M = Array{SymPy.Sym}(undef, p+1, p+1)
    F = get_hermite_basis(p)
    for i = 1:p+1
        for j = 1:p+1
            M[i, j] = simplify(integrate(diff(F[i], x, m) * diff(F[j], x, n), (x, 0, h)))
        end
    end
    M
end

"""
    get_hermite_em_varcoeff(p = 3, m = 0, n = 0, f = 1)
    Get Hermite finite elements elementary matrices for variable coefficients.

# Arguments
  * `p`: degree of polynomials in the basis.
  * `m`: number of derivatives on the first function.
  * `n`: number of derivatives on the second function.
  * `f`: the variable coefficient.
"""
function get_hermite_em_varcoeff(p = 3, m = 0, n = 0, f = 1)
    global x
    global h
    xa = SymPy.symbols("xa")
    xb = SymPy.symbols("xb")
    M = Array{SymPy.Sym}(undef, p+1, p+1)
    F = get_hermite_basis(p, true)
    for i = 1:p+1
        for j = 1:p+1
            M[i, j] = simplify(integrate(diff(F[i], x, m) * diff(F[j], x, n) * f, (x, xa, xb)))
            M[i, j] = simplify(subs(M[i, j], h, xb - xa))
        end
    end
    M
end

"""
    get_square_hermite_em((px, py) = (3, 3), (mx, my) = (0, 0), (nx, ny) = (0, 0))

Get Hermite finite elements elementary matrices for a squared element.

# Arguments
  * `(px, py)` : degree of polynomials in the basis.
  * `(mx, my)` : number of derivatives on the first function wrt x and y
  * `(nx, ny)` : number of derivatives on the second function wrt x and y

"""
function get_square_hermite_em((px, py) = (3, 3), (mx, my) = (0, 0), (nx, ny) = (0, 0))
    Mx = get_hermite_em(px, mx, nx)
    My = get_hermite_em(py, my, ny)

    node_x = [(1, 2, 1, 2, 3, 4, 3, 4, 3, 4, 3, 4, 1, 2, 1, 2)]
    node_y = [(1, 1, 2, 2, 1, 1, 2, 2, 3, 3, 4, 4, 3, 3, 4, 4)]

    kx = div(px, 2)
    ky = div(py, 2)
    p = (px + 1)*(py + 1)
    M = Array{SymPy.Sym}(undef, p, p)

    for i = 1:p
        for j = 1:p
            M[i,j] = Mx[node_x[kx][i], node_x[kx][j]] * My[node_y[ky][i], node_y[ky][j]]
        end
    end
    M
end

"""
    interpolate(fd::Matrix{Float64}, t::Vector{Float64}, ti::Vector{Float64})

Interpolates `fd` from `t` to `ti`.

# Arguments:
  * `fd`: values and derivatives of function to interpolate.
  * `t`: values in which `fd` is known.
  * `ti`: values in which `fd` is interpolated.
"""
function interpolate(fd, t, ti)
    global x
    global h
    p = size(fd, 1) # 2*p -1 is the degree of polynomials
    n = size(fd, 2) # number of samples
    if (n != length(t))
        error("The number of columns in the Matrix fd should be egal to the number of elements in the Vector t.")
    end
    if (norm(diff(diff(t))) > 1e-12 || norm(diff(diff(ti))) > 1e-12 ||
        abs(t[1]-ti[1])>1e-12 || abs(t[end] - ti[end]) > 1e-12)
        error("Vectors t and ti should be uniform partitions of [a,b] interval.")
    end
    F = get_hermite_basis(2*p - 1)
    DT = t[2] - t[1];
    FDT = []
    for i = 1:2*p
        push!(FDT, SymPy.subs(F[i], h, DT))
    end
    ni = length(ti)
    fi = collect(zeros(ni , 1))
    for i = 1:ni-1
        j = floor(ti[i] / DT) + 1
        j = convert(Int64, j)
        xx = ti[i] - t[j]
        S = 0;
        for l = 1:p
            S = S + fd[l, j] * SymPy.subs(FDT[l], x, xx)
            S = S + fd[l, j + 1] * SymPy.subs(FDT[l+p], x, xx)
        end
        fi[i] = S
        if (isnan(S))
            fi[i] = 0
        end

    end
    fi[ni] = fd[1, n]
    fi
end


text:
using Mustache
using Test

@testset " partials " begin

	## The greater-than operator should expand to the named partial.
tpl = """\"{{>text}}\""""

	@test Mustache.render(tpl, Dict{Any,Any}("text"=>"from partial")) == """\"from partial\""""

	## The empty string should be used when the named partial is not found.
tpl = """\"{{>text}}\""""

	@test Mustache.render(tpl, Dict{Any,Any}()) == """\"\""""

	## The greater-than operator should operate within the current context.
tpl = """\"{{>partial}}\""""

	@test Mustache.render(tpl, Dict{Any,Any}("partial"=>"*{{text}}*","text"=>"content")) == """\"*content*\""""

	## The greater-than operator should properly recurse.
tpl = """{{>node}}"""

	@test Mustache.render(tpl, Dict{Any,Any}("nodes"=>Dict{Any,Any}[Dict("nodes"=>Any[],"content"=>"Y")],"content"=>"X","node"=>"{{content}}<{{#nodes}}{{>node}}{{/nodes}}>")) == """X<Y<>>"""

	## The greater-than operator should not alter surrounding whitespace.
tpl = """| {{>partial}} |"""

	@test Mustache.render(tpl, Dict{Any,Any}("partial"=>"\t|\t")) == """| 	|	 |"""

	## Whitespace should be left untouched.
tpl = """  {{data}}  {{> partial}}
"""

	@test Mustache.render(tpl, Dict{Any,Any}("partial"=>">\n>","data"=>"|")) == """  |  >
>
"""

	## "\r\n" should be considered a newline for standalone tags.
tpl = """|
{{>partial}}
|"""

	@test Mustache.render(tpl, Dict{Any,Any}("partial"=>">")) == """|
>|"""

	## Standalone tags should not require a newline to precede them.
tpl = """  {{>partial}}
>"""

	@test_skip Mustache.render(tpl, Dict{Any,Any}("partial"=>">\n>")) == """  >
  >>"""

	## Standalone tags should not require a newline to follow them.
tpl = """>
  {{>partial}}"""

	@test Mustache.render(tpl, Dict{Any,Any}("partial"=>">\n>")) == """>
  >
  >"""

	## Each line of the partial should be indented before rendering.
tpl = """\\
 {{>partial}}
/
"""

	@test Mustache.render(tpl, Dict{Any,Any}("partial"=>"|\n{{{content}}}\n|\n","content"=>"<\n->")) == """\\
 |
 <
->
 |
/
"""

	## Superfluous in-tag whitespace should be ignored.
tpl = """|{{> partial }}|"""

	@test Mustache.render(tpl, Dict{Any,Any}("partial"=>"[]","boolean"=>true)) == """|[]|"""
end


text:
# In this file, we figure out how to link to Python (surprisingly complicated)
# and generate a deps/deps.jl file with the libpython name and other information
# needed for static compilation of PyCall.

# As a result, if you switch to a different version or path of Python, you
# will probably need to re-run Pkg.build("PyCall").

# remove deps.jl if it exists, in case build.jl fails
isfile("deps.jl") && rm("deps.jl")

using Compat

PYTHONIOENCODING = get(ENV, "PYTHONIOENCODING", nothing)
PYTHONHOME = get(ENV, "PYTHONHOME", nothing)

try # save/restore environment vars

# set PYTHONIOENCODING when running python executable, so that
# we get UTF-8 encoded text as output (this is not the default on Windows).
ENV["PYTHONIOENCODING"] = "UTF-8"

#########################################################################

pyconfigvar(python::AbstractString, var::AbstractString) = chomp(readall(`$python -c "import distutils.sysconfig; print(distutils.sysconfig.get_config_var('$var'))"`))
pyconfigvar(python, var, default) = let v = pyconfigvar(python, var)
    v == "None" ? default : v
end

pysys(python::AbstractString, var::AbstractString) = chomp(readall(`$python -c "import sys; print(sys.$var)"`))

#########################################################################

const dlprefix = @windows? "" : "lib"

# return libpython name, libpython pointer
function find_libpython(python::AbstractString)
    # it is ridiculous that it is this hard to find the name of libpython
    v = pyconfigvar(python,"VERSION","")
    libs = [ dlprefix*"python"*v*"."*Libdl.dlext, dlprefix*"python."*Libdl.dlext ]
    lib = pyconfigvar(python, "LIBRARY")
    lib != "None" && unshift!(libs, splitext(lib)[1]*"."*Libdl.dlext)
    lib = pyconfigvar(python, "LDLIBRARY")
    lib != "None" && unshift!(unshift!(libs, basename(lib)), lib)
    libs = unique(libs)

    # it is ridiculous that it is this hard to find the path of libpython
    libpaths = [pyconfigvar(python, "LIBDIR"),
                (@windows ? dirname(pysys(python, "executable")) : joinpath(dirname(dirname(pysys(python, "executable"))), "lib"))]
    @osx_only push!(libpaths, pyconfigvar(python, "PYTHONFRAMEWORKPREFIX"))

    # `prefix` and `exec_prefix` are the path prefixes where python should look for python only and compiled libraries, respectively.
    # These are also changed when run in a virtualenv.
    exec_prefix = pysys(python, "exec_prefix")

    push!(libpaths, exec_prefix)
    push!(libpaths, joinpath(exec_prefix, "lib"))

    if !haskey(ENV, "PYTHONHOME")
        # PYTHONHOME tells python where to look for both pure python
        # and binary modules.  When it is set, it replaces both
        # `prefix` and `exec_prefix` and we thus need to set it to
        # both in case they differ. This is also what the
        # documentation recommends.  However, they are documented
        # to always be the same on Windows, where it causes
        # problems if we try to include both.
        ENV["PYTHONHOME"] = @windows? exec_prefix : pysys(python, "prefix") * ":" * exec_prefix
        # Unfortunately, setting PYTHONHOME screws up Canopy's Python distro?
        try
            run(pipeline(`$python -c "import site"`, stdout=DevNull, stderr=DevNull))
        catch
            pop!(ENV, "PYTHONHOME")
        end
    end

    # TODO: other paths? python-config output? pyconfigvar("LDFLAGS")?

    # find libpython (we hope):
    for lib in libs
        for libpath in libpaths
            libpath_lib = joinpath(libpath, lib)
            if isfile(libpath_lib)
                try
                    return (Libdl.dlopen(libpath_lib,
                                         Libdl.RTLD_LAZY|Libdl.RTLD_DEEPBIND|Libdl.RTLD_GLOBAL),
                            libpath_lib)
                end
            end
        end
    end

    # We do this *last* because the libpython in the system
    # library path might be the wrong one if multiple python
    # versions are installed (we prefer the one in LIBDIR):
    for lib in libs
        lib = splitext(lib)[1]
        try
            return (Libdl.dlopen(lib, Libdl.RTLD_LAZY|Libdl.RTLD_DEEPBIND|Libdl.RTLD_GLOBAL),
                    lib)
        end
    end
    error("Couldn't find libpython; check your PYTHON environment variable")
end

#########################################################################

hassym(lib, sym) = Libdl.dlsym_e(lib, sym) != C_NULL

# call dlsym_e on a sequence of symbols and return the symbol that gives
# the first non-null result
function findsym(lib, syms...)
    for sym in syms
        if hassym(lib, sym)
            return sym
        end
    end
    error("no symbol found from: ", syms)
end

#########################################################################

# need to be able to get the version before Python is initialized
Py_GetVersion(libpy) = bytestring(ccall(Libdl.dlsym(libpy, :Py_GetVersion), Ptr{Uint8}, ()))

const python = get(ENV, "PYTHON", "python")
const (libpython, libpy_name) = find_libpython(python)
const programname = pysys(python, "executable")

# cache the Python version as a Julia VersionNumber
const pyversion = convert(VersionNumber, split(Py_GetVersion(libpython))[1])

info("PyCall is using $python (Python $pyversion) at $programname, libpython = $libpy_name")

if pyversion < v"2.7"
    error("Python 2.7 or later is required for PyCall")
end

# PyUnicode_* may actually be a #define for another symbol, so
# we cache the correct dlsym
const PyUnicode_AsUTF8String =
    findsym(libpython, :PyUnicode_AsUTF8String, :PyUnicodeUCS4_AsUTF8String, :PyUnicodeUCS2_AsUTF8String)
const PyUnicode_DecodeUTF8 =
    findsym(libpython, :PyUnicode_DecodeUTF8, :PyUnicodeUCS4_DecodeUTF8, :PyUnicodeUCS2_DecodeUTF8)

# Python 2/3 compatibility: cache symbols for renamed functions
if hassym(libpython, :PyString_FromString)
    const PyString_FromString = :PyString_FromString
    const PyString_AsString = :PyString_AsString
    const PyString_Size = :PyString_Size
    const PyString_Type = :PyString_Type
else
    const PyString_FromString = :PyBytes_FromString
    const PyString_AsString = :PyBytes_AsString
    const PyString_Size = :PyBytes_Size
    const PyString_Type = :PyBytes_Type
end
if hassym(libpython, :PyInt_Type)
    const PyInt_Type = :PyInt_Type
    const PyInt_FromSize_t = :PyInt_FromSize_t
    const PyInt_FromSsize_t = :PyInt_FromSsize_t
    const PyInt_AsSsize_t = :PyInt_AsSsize_t
else
    const PyInt_Type = :PyLong_Type
    const PyInt_FromSize_t = :PyLong_FromSize_t
    const PyInt_FromSsize_t = :PyLong_FromSsize_t
    const PyInt_AsSsize_t = :PyLong_AsSsize_t
end

# hashes changed from long to intptr_t in Python 3.2
const Py_hash_t = pyversion < v"3.2" ? Clong:Int

# whether to use unicode for strings by default, ala Python 3
const pyunicode_literals = pyversion >= v"3.0"

# some arguments changed from char* to wchar_t* in Python 3
pystring = pyversion.major < 3 ? "bytestring" : "wstring"

open("deps.jl", "w") do f
    print(f, """
          const python = "$(escape_string(python))"
          const libpython = "$(escape_string(libpy_name))"
          const pyprogramname = $pystring("$(escape_string(programname))")
          const pyversion_build = $(repr(pyversion))
          const PYTHONHOME = $pystring("$(escape_string(get(ENV, "PYTHONHOME", "")))")

          const PyUnicode_AsUTF8String = :$PyUnicode_AsUTF8String
          const PyUnicode_DecodeUTF8 = :$PyUnicode_DecodeUTF8

          const PyString_FromString = :$PyString_FromString
          const PyString_AsString = :$PyString_AsString
          const PyString_Size = :$PyString_Size
          const PyString_Type = :$PyString_Type
          const PyInt_Type = :$PyInt_Type
          const PyInt_FromSize_t = :$PyInt_FromSize_t
          const PyInt_FromSsize_t = :$PyInt_FromSsize_t
          const PyInt_AsSsize_t = :$PyInt_AsSsize_t

          const Py_hash_t = $Py_hash_t

          const pyunicode_literals = $pyunicode_literals
          """)
end

#########################################################################

finally # restore env vars

PYTHONIOENCODING != nothing ? (ENV["PYTHONIOENCODING"] = PYTHONIOENCODING) : pop!(ENV, "PYTHONIOENCODING")
PYTHONHOME != nothing ? (ENV["PYTHONHOME"] = PYTHONHOME) : pop!(ENV, "PYTHONHOME", "")

end


text:
function simulate! end

function simulate!(sim::Simulation; pairwise=:cells, kwargs...)
    if pairwise == :cells
        _simulate_cells!(sim; kwargs...)
    else
        _simulate_squared!(sim; kwargs...)
    end

    return nothing
end

function _setup_simulation!(system::System, potential::Potential; kwargs...)
    # Build initial cell lists
    cl = CellListMap.CellList(system.xpos, system.box; parallel=kwargs[:parallel])
    cell_cache = _build_cache(cl)

    # Create a good enough configuration
    packsystem!(system, cell_cache.cell, potential.energy)

    return cell_cache
end

function _simulate_cells!(sim::Simulation; steps=10_000, parallel=false, ishow=10_000)
    @unpack system, ensemble, potential = sim
    # Obtain the energy function for the interaction potential
    cell_cache = _setup_simulation!(system, potential; parallel=parallel)

    # Create the ensemble options
    opts = EnsembleOptions(ensemble)

    for istep in 1:steps
        opts.nattempt += 1
        uener = _mcmove!(system, potential.energy, opts, cell_cache; parallel=parallel)

        if istep % ishow == 0
            @show uener / system.npart
            @show opts.naccept / opts.nattempt
        end

        adjust!(opts)
    end

    return nothing
end

function _simulate_squared!(sim::Simulation; steps=10_000, parallel=false, ishow=10_000)
    @unpack system, ensemble, potential = sim
    # Obtain the energy function for the interaction potential
    _ = _setup_simulation!(system, potential; parallel=parallel)

    # Create the ensemble options
    opts = EnsembleOptions(ensemble)

    for istep in 1:steps
        opts.nattempt += 1
        uener = _mcmove!(system, potential.energy, opts)

        if istep % ishow == 0
            @show uener / system.npart
            @show opts.naccept / opts.nattempt
        end

        adjust!(opts)
    end

    return nothing
end


text:
# # @__NAME__

# PREAMBLE

# PKG_SETUP

# ## Setup

using Pkg # hide
Pkg.activate(@__DIR__) # hide
Pkg.instantiate() # hide
using RigidBodyDynamics
using StaticArrays
using SymPy

# ## Create symbolic parameters
# * Masses: $m_1, m_2$
# * Mass moments of inertia (about center of mass): $I_1, I_2$
# * Link lengths: $l_1, l_2$
# * Center of mass locations (w.r.t. preceding joint axis): $c_1, c_2$
# * Gravitational acceleration: $g$

inertias = @syms m_1 m_2 I_1 I_2 positive = true
lengths = @syms l_1 l_2 c_1 c_2 real = true
gravitational_acceleration = @syms g real = true
params = [inertias..., lengths..., gravitational_acceleration...]
transpose(params)

# ## Create double pendulum `Mechanism`
# A `Mechanism` contains the joint layout and inertia parameters, but no state information.

T = Sym # the 'scalar type' of the Mechanism we'll construct
axis = SVector(zero(T), one(T), zero(T)) # axis of rotation for each of the joints
double_pendulum = Mechanism(RigidBody{T}("world"); gravity = SVector(zero(T), zero(T), g))
world = root_body(double_pendulum) # the fixed 'world' rigid body

# Attach the first (upper) link to the world via a revolute joint named 'shoulder'
inertia1 = SpatialInertia(CartesianFrame3D("upper_link"),
    moment=I_1 * axis * transpose(axis),
    com=SVector(zero(T), zero(T), c_1),
    mass=m_1)
body1 = RigidBody(inertia1)
joint1 = Joint("shoulder", Revolute(axis))
joint1_to_world = one(Transform3D{T}, frame_before(joint1), default_frame(world));
attach!(double_pendulum, world, body1, joint1,
    joint_pose = joint1_to_world);

# Attach the second (lower) link to the world via a revolute joint named 'elbow'
inertia2 = SpatialInertia(CartesianFrame3D("lower_link"),
    moment=I_2 * axis * transpose(axis),
    com=SVector(zero(T), zero(T), c_2),
    mass=m_2)
body2 = RigidBody(inertia2)
joint2 = Joint("elbow", Revolute(axis))
joint2_to_body1 = Transform3D(
    frame_before(joint2), default_frame(body1), SVector(zero(T), zero(T), l_1))
attach!(double_pendulum, body1, body2, joint2,
    joint_pose = joint2_to_body1)

# ## Create `MechanismState` associated with the double pendulum `Mechanism`
# A `MechanismState` stores all state-dependent information associated with a `Mechanism`.

x = MechanismState(double_pendulum);

# Set the joint configuration vector of the MechanismState to a new vector of symbolic variables
q = configuration(x)
for i in eachindex(q)
    q[i] = symbols("q_$i", real = true)
end

# Set the joint velocity vector of the MechanismState to a new vector of symbolic variables
v = velocity(x)
for i in eachindex(v)
    v[i] = symbols("v_$i", real = true)
end

# ## Compute dynamical quantities in symbolic form

# Mass matrix
simplify.(mass_matrix(x))

# Kinetic energy
simplify(kinetic_energy(x))

# Potential energy
simplify(gravitational_potential_energy(x))

max_stars_repo_path: examples/6. Symbolics using SymPy/6. Symbolics using SymPy.jl