text:
const testcube_datadep = DataDeps.DataDep(
"VAMPIRES_data_cube",
"""
A single temporal data cube from the VAMPIRES subinstrument on Subaru/SCExAO. This data was taken at high frequency (~40 Hz) and is used to test lucky imaging techniques.
""",
"https://zenodo.org/record/5805552/files/VAMPIRES_cube.fits",
"586f9b24f6d6b2bf886766814bdf4f44cf6c119236b8c4a106963edca37dee12"
)

"""
testcube()

Return the filepath of the test cube artifact. This needs to be loaded, using [FITSIO.jl](https://github.com/JuliaAstro/FITSIO.jl), for example. This data is a sequence of frames captured on the VAMPIRES instrument on Subaru/SCExAO.[^2]

[^2]: [VAMPIRES](https://www.naoj.org/Projects/SCEXAO/scexaoWEB/030openuse.web/040vampires.web/indexm.html)

# Examples

```julia
julia> using FITSIO

julia> cube = read(FITS(testcube())[1])
```
"""
testcube() = joinpath(DataDeps.datadep"VAMPIRES_data_cube", "VAMPIRES_cube.fits")


text:
@defcomp emissions begin
	SIG0	= Variable()				#Carbon intensity 2010-2015 (kgCO2 per output 2010 USD)
	GSIG	= Variable(index=[time])	#Change in sigma (cumulative improvement of energy efficiency)
	SIGMA	= Variable(index=[time])	#CO2-equivalent-emissions output ratio
    EIND    = Variable(index=[time])    #Industrial emissions (GtCO2 per year)
	ETREE	= Variable(index=[time])	#Emissions from deforestation
    E       = Variable(index=[time])    #Total CO2 emissions (GtCO2 per year)
	CUMETREE= Variable(index=[time])	#Cumulative from land
    CCA     = Variable(index=[time])    #Cumulative industrial emissions
	CCATOT	= Variable(index=[time])	#Cumulative total carbon emissions

	gsigma1 = Parameter()				#Initial growth of sigma (per year)
	dsig	= Parameter()				#Decline rate of decarbonization (per period)
	eland0	= Parameter()				#Carbon emissions from land 2015 (GtCO2 per year)
	deland	= Parameter()				#Decline rate of land emissions (per period)
	e0		= Parameter()				#Industrial emissions 2015 (GtCO2 per year)
    MIU     = Parameter(index=[time])   #Emission control rate GHGs
    YGROSS  = Parameter(index=[time])   #Gross world product GROSS of abatement and damages (trillions 2010 USD per year)
    cca0    = Parameter()               #Initial cumulative industrial emissions
	cumetree0=Parameter()				#Initial emissions from deforestation (see GAMS code)

    function run_timestep(p, v, d, t)
		#Define SIG0
			v.SIG0 = p.e0/(p.YGROSS[TimestepIndex(1)] * (1 - p.MIU[TimestepIndex(1)]))

		#Define function for GSIG
		if is_first(t)
			v.GSIG[t] = p.gsigma1
		else
			v.GSIG[t] = v.GSIG[t-1] * ((1 + p.dsig)^5)
		end

		#Define function for SIGMA
		if is_first(t)
			v.SIGMA[t] = v.SIG0
		else
			v.SIGMA[t] = v.SIGMA[t-1] * exp(v.GSIG[t-1] * 5)
		end

        #Define function for EIND
		v.EIND[t] = v.SIGMA[t] * p.YGROSS[t] * (1 - p.MIU[t])

		#Define function for ETREE
		if is_first(t)
			v.ETREE[t] = p.eland0
		else
			v.ETREE[t] = v.ETREE[t - 1] * (1 - p.deland)
		end

        #Define function for E
        v.E[t] = v.EIND[t] + v.ETREE[t]

		#Define function for CUMETREE
		if is_first(t)
			v.CUMETREE[t] = p.cumetree0
		else
			v.CUMETREE[t] = v.CUMETREE[t-1] + v.ETREE[t-1]
		end

        #Define function for CCA
        if is_first(t)
            v.CCA[t] = p.cca0
        else
            v.CCA[t] = v.CCA[t-1] + v.EIND[t-1] * 5/3.666
        end

		#Define function for CCATOT
			v.CCATOT[t] = v.CCA[t] + v.CUMETREE[t]

    end
end

ext: jl


text:
# use Julia1 for this
using PyPlot
#  using StatPlots
using Statistics
using CSV

type_list = [Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64]
#  table = CSV.read("logs/test_gconvert.csv")
#  table = CSV.read("logs/net_transition_discount_vary_reference_solver_training.csv",types=type_list,weakrefstrings=false,rows=10)
#  table = CSV.read("logs/n_vary_reference_solver_training.csv")
table = CSV.read("logs/mturk_trusted_condition_1.csv")
#  table = CSV.read("logs/mturk_candidate_condition_4.csv")
#  for c in names(table)
    #  println(typeof(table[c]))
#  end
tbl_ary = convert(Array{Float64},table[:,1:22])
println(typeof(tbl_ary))

# normalize all data
#  for i = 1:size(tbl_ary,2)
    #  tbl_ary[:,i] = tbl_ary[:,i]./maximum(tbl_ary[:,i])
#  end

col_list = []

# detect nearly identical columns and leave them out
#  for i = 1:size(tbl_ary,2)
    #  if sum(isequal.(diff(tbl_ary[:,i]),0.0))/length(tbl_ary[:,i]) > 0.40
        #  # this is an (nearly) identical column
        #  continue
    #  else
        #  push!(col_list,i)
    #  end
#  end

println(col_list)
println(size(tbl_ary))

fig,ax = PyPlot.subplots(2,1)
fig.set_size_inches(8,8)
fsize = 15
#
ax[1].scatter(tbl_ary[:,3],tbl_ary[:,21],marker=".",c=:black,label="Monte-Carlo Results")
#  ax[1].scatter(tbl_ary[:,3],tbl_ary[:,21],c=tbl_ary[:,12])
ax[1].hlines(y=0,xmin=0,xmax=1,colors=:gray,label="rwd=0.0")
ax[1].set_title("transition probability vs. reward",size=fsize)
ax[1].set_xlabel("transition probability",size=fsize)
ax[1].set_ylabel("Reward",size=fsize)
ax[1].legend()

mean_rwd = []
for u in unique(tbl_ary[:,12])
    uniqueN = tbl_ary[:,12].==u

    append!(mean_rwd,mean(tbl_ary[:,21][uniqueN]))
end

mean_sort = sortperm(unique(tbl_ary[:,12]))

ax[2].scatter(tbl_ary[:,12],tbl_ary[:,21],marker=".",c=:black,label="Monte-Carlo Results")
ax[2].plot(unique(tbl_ary[:,12])[mean_sort],mean_rwd[mean_sort],c=:blue,label="Mean Reward")
ax[2].hlines(y=0,xmin=16,xmax=43,colors=:gray,label="rwd=0.0")
ax[2].set_title("N vs. mean reward",size=fsize)
ax[2].set_xlabel("N",size=fsize)
ax[2].set_ylabel("Reward",size=fsize)
ax[2].legend()

fig.tight_layout()

savefig("raw_data.pdf",transparent=true,dpi=300)

#  fig2,ax2 = PyPlot.subplots(1,1)
#  fig.set_size_inches(5,5)
#  fsize = 15
#
#  ax2.scatter(tbl_ary[:,3],tbl_ary[:,12],c=tbl_ary[:,21])
#  ax2.set_title("ptrans vs N",size=fsize)
#  ax2.set_xlabel("ptrans",size=fsize)
#  ax2.set_ylabel("N",size=fsize)
#
#  fig3,ax3 = PyPlot.subplots(1,1)
#  fig3.set_size_inches(5,5)
#  fsize = 15
#
#  mean_rwd = []
#  for u in unique(tbl_ary[:,12])
    #  uniqueN = tbl_ary[:,12].==u
#
    #  append!(mean_rwd,mean(tbl_ary[:,21][uniqueN]))
#  end
#
#  ax3.scatter(mean_rwd,unique(tbl_ary[:,12]))
#  ax3.set_title("mean reward vs N",size=fsize)
#  ax3.set_xlabel("mean reward",size=fsize)
#  ax3.set_ylabel("N",size=fsize)
#
#  fig4,ax4 = PyPlot.subplots(1,1)
#  fig4.set_size_inches(5,5)
#  fsize = 15

max_stars_repo_name: COHRINT/FaMSeC

max_issues_repo_name: COHRINT/FaMSeC

max_forks_repo_name: COHRINT/FaMSeC


text:
function cpp_emu_lasertag(r::Int; kwargs...)
    if r == 4
        obstacles = Set{Coord}(Coord(c) for c in ([5,7], [2,6], [4,3], [3,2], [4,2], [8,2], [10,2], [4,1]))
    else
        warn("r = $r not recognized.")
    end
    return gen_lasertag(7, 11; obstacles=obstacles, kwargs...)
end

#=
function gen_lasertag(n_rows::Int,
                      n_cols::Int,
                      n_obstacles::Int,
                      reading_std::Float64=2.5;
                      discrete=false,
                      rng=Base.GLOBAL_RNG,
                      kwargs...)

    return gen_lasertag(n_rows, n_cols;
                        obstacles=gen_obstacles(n_rows, n_cols, n_obstacles, rng),
                        discrete=discrete,
                        reading_std=reading_std,
                        rng=rng,
                        kwargs...
                       )
end
=#

function gen_lasertag(n_rows::Int,
                      n_cols::Int,
                      n_obstacles::Int,
                      obs_model::ObsModel=DESPOTEmu(Floor(n_rows, n_cols), 2.5);
                      rng=Random.GLOBAL_RNG,
                      kwargs...
                     )

    return gen_lasertag(n_rows,
                        n_cols,
                        obstacles=gen_obstacles(n_rows, n_cols, n_obstacles, rng),
                        obs_model=obs_model,
                        rng=rng;
                        kwargs...
                       )
end

function gen_obstacles(n_rows::Int, n_cols::Int, n_obstacles::Int, rng::AbstractRNG=Random.GLOBAL_RNG)
    f = Floor(n_rows, n_cols)
    obs_inds = randperm(rng, n_pos(f))[1:n_obstacles] # XXX inefficient
    obs_subs = CartesianIndices((n_cols, n_rows))[obs_inds]
    obstacles = Set{Coord}(Coord(p.I) for p in obs_subs)
    for c in obstacles
        if !inside(f, c)
            @show c
            @show obs_inds
            @show obs_subs
            error("not inside")
        end
    end
    return obstacles
end

function gen_lasertag(n_rows::Int=7,
                      n_cols::Int=11;
                      rng=Random.GLOBAL_RNG,
                      obstacles=gen_obstacles(n_rows, n_cols, 8, rng),
                      obs_model::ObsModel=DESPOTEmu(Floor(n_rows, n_cols), 2.5),
                      robot_position_known::Bool=false,
                      kwargs...)

    f = Floor(n_rows, n_cols)
    if robot_position_known
        r = Coord(rand(rng, 1:f.n_cols), rand(rng, 1:f.n_rows))
    else
        r = nothing
    end
    M = typeof(obs_model)
    O = obs_type(M)
    return LaserTagPOMDP{M, O}(;floor=f,
                                obstacles=obstacles,
                                robot_init=r,
                                obs_model=obs_model,
                                kwargs...)
end


text:
@testset "FBM" begin
    rng = MersenneTwister(42)
    h = 0.3
    k = FBMKernel(; h=h)
    v1 = rand(rng, 3)
    v2 = rand(rng, 3)
    @test k(v1, v2) ≈
          (
        sqeuclidean(v1, zero(v1))^h + sqeuclidean(v2, zero(v2))^h -
        sqeuclidean(v1 - v2, zero(v1 - v2))^h
    ) / 2 atol = 1e-5
    @test repr(k) == "Fractional Brownian Motion Kernel (h = $(h))"

    test_interface(k)
    @test repr(k) == "Fractional Brownian Motion Kernel (h = $(h))"
    test_ADs(FBMKernel; ADs=[:ReverseDiff])

    # Tests failing for ForwardDiff and Zygote@0.6.
    # Related to: https://github.com/FluxML/Zygote.jl/issues/1036
    f(x, y) = x^y
    @test_broken !isinf(
        Zygote.gradient((x, y) -> sum(f.(x, y)), zeros(1), fill(0.9, 1))[1][1]
    )

    test_params(k, ([h],))
end


text:
module SlicedNormals

using CovarianceEstimation
using Distributions
using DynamicPolynomials
using IntervalArithmetic
using LinearAlgebra
using Memoize
using MinimumVolumeEllipsoids
using NLopt
using PDMats
using TransitionalMCMC

import Base: rand

export SlicedNormal, Z, rand, pdf, fit_baseline, fit_scaling

struct SlicedNormal
    d::Integer
    μ::AbstractVector
    P::AbstractPDMat
    Δ::IntervalBox
    δ::AbstractMatrix
end

function pdf(sn::SlicedNormal, δ::AbstractVector, normalize::Bool=true)
    if δ ∈ sn.Δ
        f = exp(-_ϕ(δ, sn.μ, sn.P, sn.d))

        if normalize
            return f / c(sn.μ, sn.P, sn.d, sn.δ)
        else
            return f
        end
    else
        return 0
    end
end

function Distributions.pdf(sn::SlicedNormal, δ::AbstractMatrix)
    n, m = size(δ)
    if n == 1 || m == 1
        return pdf(sn, vec(δ))
    end
    if n < m
        return [pdf(sn, c) for c in eachcol(δ)]
    end
    return return [pdf(sn, c) for c in eachrow(δ)]
end

function Z(δ::AbstractVector, d::Integer)
    x = @polyvar x[1:length(δ)]
    z = monomials(x..., 1:d)

    # double reverse to achieve graded lexographic order
    return reverse(map(p -> p(reverse(δ)), z))
end

@memoize function c(
    μ::AbstractVector, P::AbstractPDMat, d::Integer, x::AbstractMatrix, b::Integer=10000
)
    ϵ = minimum_volume_ellipsoid(x')

    V = volume(ϵ)
    u = rand(ϵ, b)

    return V / b * sum([exp(-_ϕ(δ, μ, P, d)) for δ in eachcol(u)])
end

function bounds(Δ::IntervalBox)
    lb = getfield.(Δ, :lo)
    ub = getfield.(Δ, :hi)

    return lb, ub
end

function mean_and_covariance(x::AbstractMatrix, d::Integer)
    z = mapreduce(r -> transpose(Z(r, d)), vcat, eachrow(x))

    μ = vec(mean(z; dims=1))
    P = PDMat(inv(cov(LinearShrinkage(ConstantCorrelation()), z)))

    return μ, P
end

function fit_baseline(x::AbstractMatrix, d::Integer)
    lb = vec(minimum(x; dims=1))
    ub = vec(maximum(x; dims=1))

    Δ = IntervalBox(interval.(lb, ub)...)

    return fit_baseline(x, d, Δ)
end

function fit_baseline(x::AbstractMatrix, d::Integer, Δ::IntervalBox)
    μ, P = mean_and_covariance(x, d)
    D = sum([_ϕ(δ, μ, P, d) for δ in eachrow(x)])

    m = size(x, 1)

    cΔ = c(μ, P, d, x)
    lh = m * log(1 / cΔ) - D

    return SlicedNormal(d, μ, P, Δ, x), lh
end

function fit_scaling(x::AbstractMatrix, d::Integer)
    lb = vec(minimum(x; dims=1))
    ub = vec(maximum(x; dims=1))

    Δ = IntervalBox(interval.(lb, ub)...)

    return fit_scaling(x, d, Δ)
end

function fit_scaling(x::AbstractMatrix, d::Integer, Δ::IntervalBox, b::Integer=10000)
    μ, P = mean_and_covariance(x, d)

    D = sum([_ϕ(δ, μ, P, d) for δ in eachrow(x)])

    m = size(x, 1)

    ϵ = minimum_volume_ellipsoid(x')

    U = rand(ϵ, b)
    V = volume(ϵ)

    opt = Opt(:LN_NELDERMEAD, 1)
    opt.lower_bounds = [0.0]
    opt.xtol_rel = 1e-5
    opt.min_objective =
        (s, _) -> begin
            cΔ = V / b * sum([exp(-_ϕ(δ, μ, s[1] * P, d)) for δ in eachcol(U)])
            lh = m * log(1 / cΔ) - s[1] * D
            return -1 * lh
        end

    (lh, factor, _) = optimize(opt, [1.0])
    return SlicedNormal(d, μ, factor[1] * P, Δ, x), -1 * lh
end

function rand(sn::SlicedNormal, n::Integer)
    lb, ub = bounds(sn.Δ)

    prior = Uniform.(lb, ub)

    logprior(x) = sum(logpdf.(prior, x))
    sampler(n) = mapreduce(u -> rand(u, n), hcat, prior)
    loglikelihood(x) = log(SlicedNormals.pdf(sn, x, false))

    samples, _ = tmcmc(loglikelihood, logprior, sampler, n)

    return samples
end

function _ϕ(δ, μ, P, d)
    z = Z(δ, d)
    return ((z - μ)' * P * (z - μ)) / 2
end

end # module


text:
using Optim, HTTP, GLM, LinearAlgebra, Random, Statistics, DataFrames, DataFramesMeta, CSV

#Worked with Ahmed (our class), Aleeze, and Waleed.
function ps6()
include("create_grids.jl")

#Question 1
url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS5-ddc/busdata.csv"
df = CSV.read(HTTP.get(url).body,DataFrame)

# create bus id variable
df = @transform(df, bus_id = 1:size(df,1))

#---------------------------------------------------
# reshape from wide to long (must do this twice be-
# cause DataFrames.stack() requires doing it one
# variable at a time)
#---------------------------------------------------
# first reshape the decision variable
dfy = @select(df, :bus_id,:Y1,:Y2,:Y3,:Y4,:Y5,:Y6,:Y7,:Y8,:Y9,:Y10,:Y11,:Y12,:Y13,:Y14,:Y15,:Y16,:Y17,:Y18,:Y19,:Y20,:RouteUsage,:Branded)
dfy_long = DataFrames.stack(dfy, Not([:bus_id,:RouteUsage,:Branded]))
rename!(dfy_long, :value => :Y)
dfy_long = @transform(dfy_long, time = kron(collect([1:20]...),ones(size(df,1))))
select!(dfy_long, Not(:variable))

# next reshape the odometer variable
dfx = @select(df, :bus_id,:Odo1,:Odo2,:Odo3,:Odo4,:Odo5,:Odo6,:Odo7,:Odo8,:Odo9,:Odo10,:Odo11,:Odo12,:Odo13,:Odo14,:Odo15,:Odo16,:Odo17,:Odo18,:Odo19,:Odo20)
dfx_long = DataFrames.stack(dfx, Not([:bus_id]))
rename!(dfx_long, :value => :Odometer)
dfx_long = @transform(dfx_long, time = kron(collect([1:20]...),ones(size(df,1))))
select!(dfx_long, Not(:variable))

# join reshaped df's back together
df_long = leftjoin(dfy_long, dfx_long, on = [:bus_id,:time])
sort!(df_long,[:bus_id,:time])

#2
flx = glm(@formula(Y ~ Odometer*(Odometer)^2*RouteUsage*(RouteUsage)^2*Branded*time*(time)^2), df_long, Binomial(), LogitLink())
println(flx)

#3A
url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS5-ddc/busdata.csv"
df2 = CSV.read(HTTP.get(url).body, DataFrame)
Y = Matrix(df[:,[:Y1,:Y2,:Y3,:Y4,:Y5,:Y6,:Y7,:Y8,:Y9,:Y10,:Y11,:Y12,:Y13,:Y14,:Y15,:Y16,:Y17,:Y18,:Y19,:Y20]])
X = Matrix(df[:,[:Odo1,:Odo2,:Odo3,:Odo4,:Odo5,:Odo6,:Odo7,:Odo8,:Odo9,:Odo10,:Odo11,:Odo12,:Odo13,:Odo14,:Odo15,:Odo16,:Odo17,:Odo18,:Odo19,:Odo20]])
Z = Vector(df[:,:RouteUsage])
B = Vector(df[:,:Branded])
N = size(Y,1)
T = size(Y,2)
Xstate = Matrix(df[:,[:Xst1,:Xst2,:Xst3,:Xst4,:Xst5,:Xst6,:Xst7,:Xst8,:Xst9,:Xst10,:Xst11,:Xst12,:Xst13,:Xst14,:Xst15,:Xst16,:Xst17,:Xst18,:Xst19,:Xst20]])
Zstate = Vector(df[:,:Zst])

zval,zbin,xval,xbin,xtran = create_grids()

D = DataFrame()
                D.Odometer = kron(ones(zbin), xval)
                D.RouteUsage = kron(zval, ones(xbin))
                D.time= zeros(size(xtran,1))
                D.Branded = zeros(size(xtran,1))

function values(din,Zstate,Xstate,xtran,zbin,xbin,xval)

    FVARRY=zeros(size(xtran,1),2,T+1)
    for t=2:T
        for b=0:1
            din.time[t]=t
            din.Branded[t]=b
            p_0 = predict(θ, din)

            FVARRY[:, b+1, t] = - .9 .*log.(p_0)
        end
    end
    FTV1 = zeros(size(df,1), T)
    for i=1:size(df,1)
        row0 = (Zstate[i]-1)*xbin+1
        for t=1:T
            row1  = Xstate[i,t] + (Zstate[i]-1)*xbin
            FTV1[i,t] = (xtran[row1,:].-xtran[row0,:])'*FV[row0:row0+xbin-1,df.Branded[i]+1,t+1]
        end
    end

    return FVT1'[:]
end

ftv1 = futval(flexlog, Zstate,Xstate,xtran,zbin,xbin,E)
df_long = @transform(df_long, fv=ftv1)

theta_hat_ccp_glm = glm(@formula(Y ~ Odometer + Branded),
                    df_long, Binomial(), LogitLink(),
                    offset=df_long.fv)
println(theta_ccp)
end

@time ps6()


text:
#!/usr/bin/env julia
using Gtk, Gtk.ShortNames

mutable struct TextModule
    name::AbstractString
    fileplace::AbstractString

    TextModule() = new("", "")
end

text = TextModule()

ui = Builder(filename=(@__DIR__) * "/ui.glade")
showall(ui["win"])

writetext!(textview::GtkTextView, s::AbstractString) = get_gtk_property(textview, :buffer, GtkTextBuffer) |> x -> set_gtk_property!(x, :text, s)
gettext(textview::GtkTextView) = get_gtk_property(textview, :buffer, GtkTextBuffer) |> x -> get_gtk_property(x, :text, AbstractString)

function file_new()
    textcontent = gettext(ui["textview"])
    if !isempty(textcontent)
        info_dialog("中身消しちゃったけど許してね!")
    end
    writetext!(ui["textview"], "")
    set_gtk_property!(ui["win"], :title, "")
    text.name, text.fileplace = "", ""

    return nothing
end

function file_open()
    text.fileplace = open_dialog("Open file", ui["win"], ("*.txt", "*",))
    if !isempty(text.fileplace)
        textcontent = open(io->read(io, String), text.fileplace)
        writetext!(ui["textview"], textcontent)
        text.name = split(text.fileplace, "/")[end]
        set_gtk_property!(ui["win"], :title, text.name)
    end

    return nothing
end

function file_save()
    if isempty(text.name)
        file_save_as()
    else
        textcontent = gettext(ui["textview"])
        write(text.fileplace, textcontent)
    end

    return nothing
end

function file_save_as()
    fileplace = save_dialog("Save file")
    if !isempty(fileplace)
        text.fileplace = fileplace
        textcontent = gettext(ui["textview"])
        write(text.fileplace, textcontent)
        text.name = split(text.fileplace, "/")[end]
        set_gtk_property!(ui["win"], :title, text.name)
    end

    return nothing
end

signal_connect((x,y)->file_new(), ui["file_new"], :activate, Nothing, (), false)
signal_connect((x,y)->file_open(), ui["file_open"], :activate, Nothing, (), false)
signal_connect((x,y)->file_save(), ui["file_save"], :activate, Nothing, (), false)
signal_connect((x,y)->file_save_as(), ui["file_save_as"], :activate, Nothing, (), false)
signal_connect((x,y)->exit(), ui["file_quit"], :activate, Nothing, (), false)

if !isinteractive()
    c = Condition()
    signal_connect(ui["win"], :destroy) do widget
        notify(c)
    end
    wait(c)
end

max_stars_repo_name: goropikari/JuliaGUIAppExample

max_issues_repo_name: goropikari/JuliaGUIAppExample

max_forks_repo_name: goropikari/JuliaGUIAppExample


text:
# generateur d'Instances pour la RODD projet 4

P = 1
m = 15
n = 15
w1 = 1
w2 = 5
L = 3
g = 1.26157

# Ouvrir le fichier "output.txt" dans lequel on pourra écrire
fout = open("output30.txt", "w")

println(fout, "P = " * string(P))
println(fout, "m = " * string(m))
println(fout, "n = " * string(n))
println(fout, "w1 = " * string(w1))
println(fout, "w2 = " * string(w2))
println(fout, "L = " * string(L))
println(fout, "g = " * string(g))
println(fout, "title")

for i in 1:n
  myLine = string(i)
  for j in 1:(m)
    myLine = myLine * " " * string(rand(70:100))
  end
  println(fout, myLine)
end

close(fout)

max_stars_repo_name: boumpteryx/RODD

max_issues_repo_name: boumpteryx/RODD

max_forks_repo_name: boumpteryx/RODD


text:
using Random
"""
    digraph(E)

Create `DiGraph` from edge-list.
"""
function digraph(E)
    d = maximum(flatten(E))
    g = DiGraph(d)
    for (i, j) in E
        add_edge!(g, i, j)
    end
    g
end

"""
    graph(E)

Create `Graph` from edge-list.
"""
function graph(E)
    d = maximum(flatten(E))
    g = Graph(d)
    for (i, j) in E
        add_edge!(g, i, j)
    end
    g
end

"""
    vpairs(g)

Return the edge-list as `Pair`s.
"""
vpairs(g) = map(Pair, collect(edges(g)))

"""
    skel_oracle(g)

Compute the `skeleton` using the `dseporacle` for the DAG `g`.
"""
skel_oracle(g) = skeleton(nv(g), dseporacle, g)

"""
    pc_oracle(g)

Compute CPDAG using the PC algorithm using the `dseporacle` on the DAG `g`.
"""
pc_oracle(g) = pcalg(nv(g), dseporacle, g)

"""
    randdag(n, alpha = 0.1)

Create random DAG from randomly permuted random triangular matrix with
edge probability `alpha`.
"""
function randdag(n, alpha = 0.1)
    g = DiGraph(n)
    p = randperm(n)
    for i in 1:n
        for j in i+1:n
            rand() < alpha && add_edge!(g, p[i], p[j])
        end
    end
    g
end


text:
function process(realname::String, simname::String)
    output = open("$(basename(realname)).csv", "w")
    simreader = FASTA.Reader(open(simname, "r"))
    try
        println(output, "Seq1, Seq2, dN, dS, Real")
        sequence_records = open(realname, "r") do file
            collect(FASTA.Reader(file))
        end
        println("Processing real dataset...")
        process_seq_records(output, sequence_records, true)
        println("Processing simulated dataset...")
        while !eof(simreader)
            @inbounds for i in 1:endof(sequence_records)
                read!(simreader, sequence_records[i])
            end
            process_seq_records(output, sequence_records, false)
        end
        println("Finished.")
    finally
        close(output)
        close(simreader)
    end
end

@inline function process_seq_records(output, records, isreal)
    write_rep_to_file(output,
                      FASTA.identifier.(records),
                      pairwise_dNdS(NG86,
                      FASTA.sequence.(records)),
                      isreal)
end

@inline function write_rep_to_file(io, names, results, real)
    @inbounds for i ∈ 1:endof(names), j ∈ (i + 1):endof(names)
        dN, dS = results[i, j]
        println(io, names[i], ", ", names[j], ", ", dN, ", ", dS, ", ", real)
    end
end


text:
using TimeSeries, MarketData
FactCheck.setstyle(:compact)
FactCheck.onlystats(true)

facts("readwrite parses IO correctly") do

io  = "DateTime,Open
       2010/1/4|9:01:00,7300
       2010/1/4|9:02:00,7316
       2010/1/4|9:03:00,7316
       2010/1/4|9:04:00,7316
       2010/1/4|9:05:00,7316"

tm0 = readtimearray(IOBuffer(io), format="yyyy/mm/dd|HH:MM:SS")

    context("Input as IOBuffer and specifying DateTime string format for reading") do
        @fact length(tm0)      --> 5
        @fact size(tm0.values) --> (5,1)
        @fact tm0.timestamp[4] --> DateTime(2010,1,4,9,4)
    end
end

facts("readwrite parses csv correctly") do

tm1 = readtimearray(joinpath(dirname(@__FILE__), "data/datetime3.csv"), format="yyyy/mm/dd|HH:MM:SS")
tm2 = readtimearray(joinpath(dirname(@__FILE__), "data/datetime3.csv"), format="yyyy/mm/dd|HH:MM:SS", meta="foo")
tm3 = readtimearray(joinpath(dirname(@__FILE__), "data/read_example_delim.csv"), format="dd/mm/yyyy HH:MM", delim=';')

    context("Specifying DateTime string format for reading") do
        @fact length(tm1)      --> 5
        @fact size(tm1.values) --> (5,1)
        @fact tm1.timestamp[4]  --> DateTime(2010,1,4,9,4)
        @fact_throws readtimearray(joinpath(dirname(@__FILE__), "data/datetime3.csv"))
    end

    context("readtimearray accepts meta field") do
        @fact tm2.meta --> "foo"
    end

    context("readtimearray works with arbitrary delimiters") do
        @fact length(tm3)      --> 2
        @fact size(tm3.values) --> (2,2)
        @fact tm3.timestamp[2] --> DateTime(2015,1,1,1,0)
        @fact tm3.values[1,1]  --> 10.42
    end
end

facts("readwrite parses MarketData objects correctly") do

    context("1d values array works") do
        @fact typeof(cl.values) --> Array{Float64,1}
    end

    context("2d values array works") do
        @fact typeof(ohlc.values) --> Array{Float64,2}
    end

    context("timestamp parses to correct type") do
        @fact typeof(cl.timestamp)        --> Vector{Date}
        @fact typeof(datetime1.timestamp) --> Vector{DateTime}
    end
end

facts("writetimearray method works") do
    filename = "$(randstring()).csv"
    uohlc    = uniformspace(ohlc)
    writetimearray(uohlc, filename)
    readback = readtimearray(filename)

    context("writetimearray output can be round-tripped") do
        @fact uohlc.colnames                         --> readback.colnames
        @fact uohlc.timestamp                        --> readback.timestamp
        @fact isequal(uohlc.values, readback.values) --> true
        rm(filename)
    end
end


text:
module BranchTests

using Test

## Runtime

struct Tag
    name :: String
    idx :: Int
end
Base.show(io :: IO, v :: Tag) = print(io, "Tag[$(v.name)($(v.idx))]")

Base.@kwdef mutable struct Node
    parent   :: Union{Nothing, Node} = nothing
    children :: Vector{Node}         = []
    tag      :: Tag
end
function Base.show(io :: IO, v :: Node)
    print(io, "Node($(v.tag))[")
    if !isnothing(v.parent)
        print(io, "p=$(v.parent)")
    end
    if !isempty(v.children)
        print(io, "c={")
        for c in v.children
            print(io, c.tag)
            print(io, ",")
        end
        print(io, "}")
    end
    print(io, "]")
end

tag(n :: Node) = n.tag
is_branch(n :: Node) = !isempty(n.children)
is_leaf(n :: Node) = !is_branch(n)

function collect_children!(res, n :: Node)
    append!(res, n.children)
    for c in n.children
        collect_children!(res, c)
    end
end

# includes `n`
function descendants(n :: Node)
    res = []
    push!(res, n)
    collect_children!(res, n)
    res
end

function collect_parents!(res, n :: Node)
    if isnothing(n.parent)
        return
    end
    push!(res, n.parent)
    collect_parents!(res, n.parent)
end

# includes `n`
function ancestors(n :: Node)
    res = []
    push!(res, n)
    collect_parents!(res, n)
    res
end

function collect_leaves!(res, n :: Node)
    if isempty(n.children)
        push!(res, n)
    else
        for c in n.children
            collect_leaves!(res, c)
        end
    end
end

struct Tree
    root :: Node
end

function leaves(t :: Tree)
    res = []
    collect_leaves!(res, t.root)
    res
end

function node(t :: Tree, tag :: Tag)
    res = descendants(t.root)
    res[findfirst(n -> n.tag == tag, res)]
end

Base.@kwdef mutable struct Run
    tree     :: Tree
    failed   :: Set{Tag} = Set()
    finished :: Set{Tag} = Set()
end
Run(tree :: Tree) = Run(tree = tree)
Base.length(r :: Run) = length(leaves(r.tree))

function should_run(r :: Run, t :: Tag)
    n = node(r.tree, t)
    if isempty(intersect(map(tag, ancestors(n)), r.failed)) # not failed
        (is_branch(n) && !isempty(setdiff(map(tag, descendants(n)), r.finished))) || !(t in r.finished)
    else
        false
    end
end
should_run(r :: Run) = !isempty(setdiff(map(tag, leaves(r.tree)), union(r.failed, r.finished)))
failed!(r :: Run, t :: Tag) = union!(r.failed, map(tag, descendants(node(r.tree, t))))
finished!(r :: Run, t :: Tag) = push!(r.finished, t)

struct BranchTestSet{T <: Test.AbstractTestSet} <: Test.AbstractTestSet
    node    :: Tag
    testset :: T
end

BranchTestSet(desc :: AbstractString; node :: Expr, testsettype = Test.DefaultTestSet, kwargs...) =
    # TODO: what's the proper way to pass `node` as Tag instead of Expr to avoid `eval`?
    BranchTestSet{testsettype}(eval(node), testsettype(desc; kwargs...))

Test.record(ts :: BranchTestSet, t :: Test.AbstractTestSet) = Test.record(ts.testset, t)

function Test.record(ts :: BranchTestSet, t :: Union{Test.Broken, Test.Pass, Test.Fail, Test.Error})
    Test.record(ts.testset, t)
    if isa(t, Test.Fail) || isa(t, Test.Error)
        failed!(task_local_storage(:__TESTBRANCHRUN__), ts.node)
    end
    t
end

function Test.finish(ts :: BranchTestSet)
    finished!(task_local_storage(:__TESTBRANCHRUN__), ts.node)
    Test.finish(ts.testset)
end

function run(f, r :: Run)
    testsets = []
    try
        task_local_storage(:__TESTBRANCHRUN__, r)
        n = 0
        while should_run(r) && n <= length(r)
            n += 1
            push!(testsets, f(r))
        end
    finally
        task_local_storage(:__TESTBRANCHRUN__, nothing)
    end
    return testsets
end

## Macros

function parse_testbranch_args(args)
    ln = nothing
    name = nothing
    testsetargs = []
    for arg in args
        if isa(arg, Symbol) || (isa(arg, Expr) && arg.head === Symbol(".")) # (qualified) testset name
            push!(testsetargs, Expr(Symbol("="), :testsettype, arg))
        elseif isa(arg, LineNumberNode)
            ln = arg
        elseif isa(arg, AbstractString) #|| (isa(arg, Expr) && arg.head === :string) && name = esc(arg)
            # TODO: name should be a Union{String, Expr} or Expr
            name = arg
        elseif isa(arg, Expr) && arg.head === :(=)
            push!(testsetargs, arg)
        else
            error("Unexpected argument $arg to @testset")
        end
    end
    (ln, name, testsetargs)
end

function testbranch_expr(input_expr)
    node_idx = 1
    parent = nothing
    node_ = nothing
    root = nothing
    @gensym state_ done_

    process = expr -> begin
        if expr.head == :macrocall && expr.args[1] == Symbol("@testbranch")
            if length(expr.args) < 2
                error("Too few arguments to @testbranch: $(expr.args[1:end])")
            end

            (ln, name, testsetargs) = parse_testbranch_args(expr.args[2:end-1])
            name = isnothing(name) ? "branch test set" : name
            body = expr.args[end]
            tag = Tag(name, node_idx)
            node_idx += 1

            node_ = Node(parent = parent, tag = tag)
            if isnothing(root)
                root = node_
            end
            if !isnothing(parent)
                push!(parent.children, node_)
            end
            parent = node_

            new_body = process(body)

            node_ = parent
            if !isnothing(parent)
                parent = parent.parent
            end

            inner = quote
                Test.@testset BranchTestSet node=:(Tag($$(tag.name), $$(tag.idx))) $(testsetargs...) $name begin
                    try
                        $new_body
                    finally
                        $(if is_leaf(node_) quote $done_ = true end end)
                    end
                end
            end

            Expr(:block, ln, quote
                $(if isnothing(parent) quote local $done_ = false end end)
                if !$done_ && BranchTests.should_run($state_, $tag)
                    $(if isnothing(parent) Expr(:return, inner) else inner end)
                end
            end)
        else
            new_expr = copy(expr)
            for (idx, arg) in enumerate(expr.args)
                if isa(arg, Expr)
                    new_expr.args[idx] = process(arg)
                end
            end
            new_expr
        end
    end
    expr = process(input_expr)
    tree = BranchTests.Tree(root)
    testset_name = "Branches - $(root.tag.name)($(length(leaves(tree))))"
    quote
        $state_ = $tree
        Test.@testset $testset_name begin
            BranchTests.run(BranchTests.Run($state_)) do $state_
                $expr
            end
        end
    end
end

"""
@testbranch [CustomTestSet] [option=val  ...] ["description"] begin ... end

Starts a new test branch.

If no custom testset type is given it defaults to creating a
`Test.DefaultTestSet`. See `Test.@testset` documentation.

Returns an array of passing test sets.
"""
macro testbranch(args...)
    if args[end].head == :Expr
        error("Expected a begin..end block")
    end
    esc(testbranch_expr(:(@testbranch $(args...))))
end

export BranchTestSet, @testbranch

end


text:
@testset "datafreetree" begin
    data = rand(2,100)
    data2 = rand(2,100)
    data3 = rand(3,100)
    t = DataFreeTree(KDTree, data)
    @test_throws ArgumentError injectdata(t, data2)
    @test_throws DimensionMismatch injectdata(t, data3)
    for typ in [KDTree, BallTree]
        dfilename = tempname()
        rfilename = tempname()
        d = 2
        n = 100
        data = Mmap.mmap(dfilename, Matrix{Float32}, (d, n))
        data[:] = rand(Float32, d, n)
        reorderbuffer = Mmap.mmap(rfilename, Matrix{Float32}, (d, n))
        t = injectdata(DataFreeTree(typ, data), data)
        tr = injectdata(DataFreeTree(typ, data, reorderbuffer = reorderbuffer), reorderbuffer)
        for i = 1:n
            @test knn(t, data[:,i], 3) == knn(tr, data[:,i], 3)
        end
        rm(dfilename)
        rm(rfilename)
    end

    data = rand(2,1000)
    buf = zeros(data)
    for typ in [KDTree, BallTree]
        t = injectdata(DataFreeTree(typ, data, indicesfor = :data), data)
        t2 = injectdata(DataFreeTree(typ, data, reorderbuffer = buf, indicesfor = :reordered), buf)
        @test data[:,knn(t, data[:,1], 3)[1]] == buf[:,knn(t2, data[:,1], 3)[1]]
    end
end


text:
using AHPQ
using Test

@testset "builder" begin
    using Statistics: norm

    n_dp = 1000
    n_dim = 4
    data = rand(n_dim, n_dp)
    data = data ./ mapslices(norm, data, dims=1)

    stopcond = 1e-1

    @info("\nTesting preclustering")
    ahpq = builder(deepcopy(data); T=0.2, a=10, increment_steps = 0, verbose=true, stopcond=stopcond);
    ahpq = builder(deepcopy(data); T=0, a=10, increment_steps = 0, verbose=true, stopcond=stopcond);
    @test true

    @info("\nTesting Using different loss functions for preclustering")
    ahpq = builder(deepcopy(data); T=0.2, a=10, increment_steps = 0, verbose=true, stopcond=stopcond, T_preclustering=0.3); # Anisotropic, T=0.3
    ahpq = builder(deepcopy(data); T=0.2, a=10, increment_steps = 0, verbose=true, stopcond=stopcond, T_preclustering=-1);  # Anisotropic, T=T
    ahpq = builder(deepcopy(data); T=0.2, a=10, increment_steps = 0, verbose=true, stopcond=stopcond, T_preclustering=0);   # L2
    @test true

    @info("\nTesting L2 initialisation and incremental initialisation")
    ahpq = builder(deepcopy(data); T=0.2, a=0, increment_steps = 2, training_points=0, verbose=true, stopcond=stopcond);
    ahpq = builder(deepcopy(data); T=0.2, a=0, increment_steps = 0, initialise_with_l2_loss=true, verbose=true, stopcond=stopcond);
    @test true

    @info("\nTesting Inverted Index ")
    ahpq = builder(deepcopy(data); T=0, a=0, increment_steps = 0, verbose=true, stopcond=stopcond, inverted_index=true);
    ahpq = builder(deepcopy(data); T=0.2, a=0, increment_steps = 0, verbose=true, stopcond=stopcond, inverted_index=true);
    @test true

    @info("\nTesting different optimisation methods")
    ahpq = builder(deepcopy(data); T=0, a=0, increment_steps = 0, verbose=true, stopcond=stopcond, optimisation="exact");
    ahpq = builder(deepcopy(data); T=0, a=0, increment_steps = 0, verbose=true, stopcond=stopcond, optimisation="nesterov");
    ahpq = builder(deepcopy(data); T=0.2, a=0, increment_steps = 0, verbose=true, stopcond=stopcond, optimisation="exact");
    ahpq = builder(deepcopy(data); T=0.2, a=0, increment_steps = 0, verbose=true, stopcond=stopcond, optimisation="nesterov");
    @test true

    @info("\nTesting multithreading")
    ahpq = builder(deepcopy(data); T=0, a=0, increment_steps = 0, verbose=true, stopcond=stopcond, optimisation="exact", multithreading=true);
    ahpq = builder(deepcopy(data); T=0, a=0, increment_steps = 0, verbose=true, stopcond=stopcond, optimisation="nesterov", multithreading=true);
    ahpq = builder(deepcopy(data); T=0.2, a=0, increment_steps = 0, verbose=true, stopcond=stopcond, optimisation="exact", multithreading=true);
    ahpq = builder(deepcopy(data); T=0.2, a=0, increment_steps = 0, verbose=true, stopcond=stopcond, optimisation="nesterov", multithreading=true);
    @test true

    @info("\nTesting reorder function")
    ahpq = builder(deepcopy(data); T=0.2, a=0, increment_steps = 0, verbose=true, stopcond=stopcond, optimisation="nesterov", reorder=250);
    @test true
end


text:
using DataFrames

A = DataFrame(Age = [27, 18, 28, 18, 28], Name = ["Jonah", "Alan", "Glory", "Popeye", "Alan"])
B = DataFrame(Name = ["Jonah", "Jonah", "Alan", "Alan", "Glory"],
    Nemesis = ["Whales", "Spiders", "Ghosts", "Zombies", "Buffy"])
AB = join(A, B, on = :Name)

@show A B AB

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
const buf8bytes = zeros(UInt8, 8)
const buf5int32 = zeros(Int32, 5)

function readvariablerecord!(io::IO, nms::Vector{Symbol}, typs::Vector{DataType},
                            widths::Vector{Int32}, labels::Vector{String},
                            msngvals::Vector{Vector{Float64}}, readwrite::Vector{Int32})
    v = read!(io, buf5int32)
    width = v[1]
    width < 0 && return
    push!(typs, iszero(width) ? Float64 : String)
    push!(widths, iszero(width) ? 8 : width)
    read!(io, buf8bytes)
    if 0x00 in buf8bytes
        error("got name $buf8bytes at position $(length(typs))")
    else
        push!(nms, Symbol(strip(String(buf8bytes))))
    end
    label = ""
    if Bool(v[2])
        label_len = read(io, Int32)
        label = strip(String(read(io, label_len)))
        if (rem = mod(label_len, 4)) ≠ 0
            skip(io, 4 - rem)
        end
    end
    push!(labels, label)
    push!(msngvals, read!(io, Vector{Float64}(v[3])))
    push!(readwrite, v[4])
    push!(readwrite, v[5])
end

struct VariableDictionary
    nms::Vector{Symbol}
    invnms::Dict{Symbol,Int}  # dictionary mapping names to positions
    typs::Vector{DataType}
    widths::Vector{Int32}
    labels::Vector{String}
    valuelabels::Dict{Symbol, Union{Dict{Float64,String},Dict{UInt8,String}}}
    msngvals::Vector{Vector{Float64}}

    readwritefmt::Matrix{Int32}
end

function readvariabledictionary(io::IO, casesz)
    nms = sizehint!(Symbol[], casesz)
    invnms = sizehint!(Dict{Symbol,Int}(), casesz)
    typs = sizehint!(DataType[], casesz)
    wdths = sizehint!(Int32[], casesz)
    labels = sizehint!(String[], casesz)
    msngvals = sizehint!(Vector{Float64}[], casesz)
    rdwrt = sizehint!(Int32[], 2 *casesz)
    rectyp = Int32(2)
    while rectyp == 2
        readvariablerecord!(io, nms, typs, wdths, labels, msngvals, rdwrt)
        rectyp = read(io, Int32)
    end
    if all(isempty, msngvals)
        msngvals = Vector{Float64}[]
    end
    for (i,v) in enumerate(nms)
        invnms[v] = i
    end
    vdicts = [Dict{T,String}() for T in typs]
    rdwrt = reshape(rdwrt, (2, length(typs)))
    valuelabels = Dict{Symbol,Union{Dict{Float64,String},Dict{UInt8,String}}}()
    rectyp, VariableDictionary(nms, invnms, typs, wdths, labels, valuelabels, msngvals, rdwrt)
end


text:
export @creator, @extractor, @updater

"""
    @creator

Put `@creator` inside `@aml` to add a custom code to the creator function (DOM creation when the struct is instanciated).

The location of the macro specifies its location in the function. For example, putting it at the begining, adds the code to begining of creator function.

This macro only affects creation (not extraction/updating), but can be used in combination with other macros.

In the following example `IQ` and `average` are calculated automatically.

# Example
```julia
@aml mutable struct Student "student"

    # add custom code to the begining of creator function
    # the following automatically fills IQ based on the name
    @creator begin
        if occursin("smart", name)
            name = replace(name, "-smart" => "") # remove smart from name
            IQ = "smart"
        elseif occursin("genius", name)
            name = replace(name, "-genius" => "") # remove smart from name
            IQ = "genius"
        else
            error("Give a smart student.")
        end
    end

    name::String, "~"
    GPA::Float64, "~"
    IQ::UN{String} = nothing, att"~" # default to nothing, but filled automatically by first @cretor macro
end

@aml mutable struct MathClass "math-class"
    @creator begin
        GPAsum = 0.0
        for student in students
            GPAsum = GPAsum + student.GPA
        end
        average = GPAsum / length(students) # fill average field
    end

    students::Vector{Student}, "student"
    average::UN{Float64} = nothing, "average" # calculated automatically
end

smarts = [Student(name = "Jack-smart", GPA = 2.0), Student(name = "Sara-genius", GPA = 5.0)]
mathclass = MathClass(students = smarts)
```
"""
macro creator(exp)
    return :creator, exp
end

"""
    @extractor

Put `@extractor` inside `@aml` to add a custom code to the extractor function (DOM parsing when a html/xml text is used for instanciation of a struct).

The location of the macro specifies its location in the function. For example, putting it at the begining, adds the code to begining of extractor function.

This macro only affects creation (not creation/updating), but can be used in combination with other macros.

Be careful that setting struct fields using `@extractor` only changes the struct field and not the xml code.

In the following example `log` is filled automatically (which doesn't have an associated xml element).

# Example
```julia
@aml mutable struct Student "student"

    name::String, "~"
    GPA::Float64, "~"

    log::UN{String} = nothing, "~"

    # add custom code to the end of extractor function
    @extractor begin
        if GPA > 4.0
            log = "A genius with a GPA of \$GPA is found" # setting fields using @extractor only changes the field and not the xml code
        end
    end
end

@aml mutable struct MathClass "math-class"
    students::Vector{Student}, "student"
end

smarts = [Student(name = "Jack-smart", GPA = 2.0), Student(name = "Sara-genius", GPA = 5.0)]
mathclass = MathClass(students = smarts)
```
"""
macro extractor(exp)
    return :extractor, exp
end

"""
    @updater

Put `@updater` inside `@aml` to add a custom code to the updater function (DOM updating after instanciation of a struct).

The location of the macro specifies its location in the function. For example, putting it at the begining, adds the code to begining of updater function.

This macro only affects updating (not creation/extraction), but can be used in combination with other macros.

See [`@creator`](@ref) and [`@extractor`](@ref) for some examples.
"""
macro updater(exp)
    return :updater, exp
end


text:
module VentisController
# app/resources/books/BooksController.jl
# app/resources/books/BooksController.jl
using Genie.Renderer.Html
# using Genie.@params
# using Genie.serve.io
using Images,TestImages
using Genie.Renderer
using ImageShow,ImageView
using ImageMorphology
using Images
using FileIO

img1=testimage("lighthouse")
save("img.png",img1)
struct Book
  title::String
  author::String
end

const BillGatesBooks = Book[
  Book("The Best We Could Do", "Thi Bui"),
  Book("Evicted: Poverty and Profit in the American City", "Matthew Desmond"),
  Book("Believe Me: A Memoir of Love, Death, and Jazz Chickens", "Eddie Izzard"),
  Book("The Sympathizer", "Viet Thanh Nguyen"),
  Book("Energy and Civilization, A History", "Vaclav Smil")
]

# function billgatesbooks()
#   "
#
Bill Gates' list of recommended books

#

#     $(["
$(book.title) by $(book.author)
" for book in BillGatesBooks]...)
#

#   "
# end

function upload_img()

end

function billgatesbooks()
  for (name,file) in @params(:FILES)
    write("./public/img/genie/$(file.name)", IOBuffer(file.data))
  end

  println("After writing")
  img_final = img1
  for (name,file) in @params(:FILES)
    img_final = load("./public/img/genie/$(file.name)")
    println("hello")
  end

  println(typeof(img_final))
  println("hello")
  # What thresholding algorithm the user has chosen. WORKS
  # algorithm = @params(:algorithm)
  # println(algorithm)
  # # The image that the user has uploaded. ERROR IS HERE
  #img = @params(:image)
  res = img1
  println(typeof(img1))
  img_erode = Gray.(Gray.(img1) .< 0.8); # keeps white objects white
  res = erode(img_erode)
  save("image.png",img1)
  imshow(res)
  for (name,file) in @params(:FILES)
    write("./public/img/genie/$(file.name)", IOBuffer(file.data))
  end
  # img_erosion2 = erode(erode(img_erode))
  # img_erosion3 = erode(erode(erode(img_erode)))
  #mosaicview(img_erode, img_erosion1, img_erosion2, img_erosion3; nrow=1)

  result =1;
  algorithm=1;
  # img=testimage("lighthouse")
  # save("img.png",img)
  # save("appe.png",img)
  println("Test Pass 2")
  # data ="appe.png"
  println(pwd())
  html(:ventis, :billgatesbooks, books = BillGatesBooks,image = img1,layout=:new1)
end

# function threshold()
#   println("threshold")
#   # What thresholding algorithm the user has chosen. WORKS
#   algorithm = @params(:algorithm)

#   # # The image that the user has uploaded. ERROR IS HERE
#   #img = @params(:image)
#   result =1;
#   algorithm=1;
#   img=testimage("lighthouse")
#   # save("appe.png",img)
#   println("Test Pass 2")
#   data ="appe.png"
#   println(pwd())
#   html(:ventis, :billgatesbooks, res = result,image =img, alg=algorithm,layout=:new1)
# end
# # BooksController.jl
# function billgatesbooks()
#   html(:books, "billgatesbooks.jl.md", books = BillGatesBooks)
# end

function showimg()
   img= testimage("lighthouse")
end

# function billgatesbooks1()

#   [
#     Html.h1() do
#       "Image Restoration"
#     end
#     Html.ul() do
#       @foreach(BillGatesBooks) do book
#         Html.li() do
#           book.title * " by 1 " * book.author
#         end
#       end
#     end
#   ]
# end

# function billgatesbooks()
#   "
#
Bill Gates' list of recommended books

#

#     $(["
$(book.title) by $(book.author)
" for book in BillGatesBooks]...)
#

#   "
# end

end


text:
#!/usr/bin/env bash
# -*- mode: julia -*-
#=
exec julia --startup-file=no -q --project=.
=#

using PNML, EzXML
#using SnoopCompileCore
#using Profile, ProfileView

str = """<?xml version="1.0"?>
<pnml xmlns="http://www.pnml.org/version-2009/grammar/pnml">
<net id="net0" type="nonstandard">
<page id="page0">
    <place id="rabbits"> <initialMarking> <text>100.0</text> </initialMarking> </place>
    <place id="wolves">  <initialMarking> <text>10.0</text> </initialMarking> </place>
    <transition id ="birth">     <rate> <text>0.3</text> </rate> </transition>
    <transition id ="death">     <rate> <text>0.7</text> </rate> </transition>
    <transition id ="predation"> <rate> <text>0.015</text> </rate> </transition>
    <arc id="a1" source="rabbits"   target="birth">     <inscription><text>1</text> </inscription> </arc>
    <arc id="a2" source="birth"     target="rabbits">   <inscription><text>2</text> </inscription> </arc>
    <arc id="a3" source="wolves"    target="predation"> <inscription><text>1</text> </inscription> </arc>
    <arc id="a4" source="rabbits"   target="predation"> <inscription><text>1</text> </inscription> </arc>
    <arc id="a5" source="predation" target="wolves">    <inscription><text>2</text> </inscription> </arc>
    <arc id="a6" source="wolves"    target="death">     <inscription><text>1</text> </inscription> </arc>
</page>
</net>
</pnml>
"""
xmlroot(s::String) = root(EzXML.parsexml(s))
const pnmlroot = xmlroot(str)
#const reg = PNML.IDRegistry()

#tinf = @snoopi_deep begin end
function test1(proot, n=10)
    for i in 1:n
        model = parse_pnml(proot; reg=PNML.IDRegistry())
    end
end

#test1(pnmlroot)
#PNML.reset_registry!(reg)
#ProfileView.
#@profview test1(pnmlroot)
#@profview PNML.IDRegistry()


text:
using LinearAlgebra, ProgressMeter
using KitBase

begin
    cd(@__DIR__)
    cells, points = KitBase.read_mesh("../assets/mesh/linesource.su2")
    cellid = KitBase.extract_cell(cells)
    facePoints, faceCells, cellNeighbors = KitBase.mesh_connectivity_2D(cellid)
    cellType = KitBase.mesh_cell_type(cellNeighbors)
    cellArea = KitBase.mesh_area_2D(points, cellid)
    cellCenter = KitBase.mesh_center_2D(points, cellid)
    faceCenter = KitBase.mesh_face_center(points, facePoints)
    cellFaces = KitBase.mesh_cell_face(cellid, faceCells)
    faceType = mesh_face_type(faceCells, cellType)
    ps = KitBase.UnstructPSpace(
        cells,
        points,
        cellid,
        cellType,
        cellNeighbors,
        cellFaces,
        cellCenter,
        cellArea,
        facePoints,
        faceCells,
        faceCenter,
        faceType,
    )
end
ps = KitBase.UnstructPSpace("../assets/mesh/linesource.su2")

begin
    # quadrature
    quadratureorder = 10
    points, weights = KitBase.legendre_quadrature(quadratureorder)
    #points, triangulation = KitBase.octa_quadrature(quadratureorder)
    #weights = KitBase.quadrature_weights(points, triangulation)
    nq = size(points, 1)
    vs = KitBase.UnstructVSpace(-1.0, 1.0, nq, points, weights)

    # IC
    s2 = 0.03^2
    init_field(x, y) = max(1e-4, 1.0 / (4.0 * π * s2) * exp(-(x^2 + y^2) / 4.0 / s2))

    # particle
    SigmaS = ones(size(cellid, 1))
    SigmaA = zeros(size(cellid, 1))
    SigmaT = SigmaS + SigmaA

    # time
    tspan = (0.0, 0.3)
    cfl = 0.7
end

ctr = Array{KitBase.ControlVolumeUS1F}(undef, size(ps.cellid, 1))
for i in eachindex(ctr)
    n = Vector{Float64}[]
    for j = 1:3
        push!(
            n,
            KitBase.unit_normal(
                ps.points[facePoints[cellFaces[i, j], 1], :],
                ps.points[facePoints[cellFaces[i, j], 2], :],
            ),
        )

        if dot(ps.faceCenter[ps.cellFaces[i, j], :] .- ps.cellCenter[i, :], n[j]) < 0
            n[j] .= -n[j]
        end
    end

    phi = zeros(nq)
    phi .= init_field(ps.cellCenter[i, 1], ps.cellCenter[i, 2])
    #phi = ones(nq) .* 1e-4
    #if -0.01 <= ps.cellCenter[i, 1] <= 0.01 && -0.01 <= ps.cellCenter[i, 2] <= 0.01
    #    phi .= 1.0
    #end

    w = sum(weights .* phi)
    dx = [
        KitBase.point_distance(
            cellCenter[i, :],
            ps.points[ps.cellid[i, 1], :],
            ps.points[ps.cellid[i, 2], :],
        ),
        KitBase.point_distance(
            cellCenter[i, :],
            ps.points[ps.cellid[i, 2], :],
            ps.points[ps.cellid[i, 3], :],
        ),
        KitBase.point_distance(
            cellCenter[i, :],
            ps.points[ps.cellid[i, 3], :],
            ps.points[ps.cellid[i, 1], :],
        ),
    ]

    ctr[i] = KitBase.ControlVolumeUS1F(n, cellCenter[i, :], dx, w, w, phi)
end

face = Array{KitBase.Interface2D1F}(undef, size(ps.facePoints, 1))
for i in eachindex(face)
    len = norm(ps.points[facePoints[i, 1], :] .- ps.points[facePoints[i, 2], :])
    n = KitBase.unit_normal(ps.points[facePoints[i, 1], :], ps.points[facePoints[i, 2], :])

    if !(-1 in ps.faceCells[i, :])
        n0 = ps.cellCenter[ps.faceCells[i, 2], :] .- ps.cellCenter[ps.faceCells[i, 1], :]
    else
        n0 = zero(n)
    end
    if dot(n, n0) < 0
        n .= -n
    end

    fw = 0.0
    ff = zeros(nq)

    face[i] = KitBase.Interface2D1F(len, n[1], n[2], fw, ff)
end

dt = 1.2 / 150 * cfl
nt = tspan[2] ÷ dt |> Int
@showprogress for iter = 1:nt
    @inbounds Threads.@threads for i in eachindex(face)
        velo = vs.u[:, 1] .* face[i].n[1] + vs.u[:, 2] .* face[i].n[2]
        if !(-1 in ps.faceCells[i, :])
            KitBase.flux_kfvs!(
                face[i].ff,
                ctr[ps.faceCells[i, 1]].f,
                ctr[ps.faceCells[i, 2]].f,
                velo,
                dt,
            )
        end
    end

    @inbounds Threads.@threads for i in eachindex(ctr)
        if ps.cellType[i] == 0
            for j = 1:3
                dirc = sign(dot(ctr[i].n[j], face[ps.cellFaces[i, j]].n))
                @. ctr[i].f -=
                    dirc * face[ps.cellFaces[i, j]].ff * face[ps.cellFaces[i, j]].len /
                    ps.cellArea[i]
            end

            integral = KitBase.discrete_moments(ctr[i].f, vs.weights)
            integral /= 4.0 * π
            @. ctr[i].f += (integral - ctr[i].f) * dt
            ctr[i].w = sum(ctr[i].f .* vs.weights)
        end
    end
end

#KitBase.write_vtk(ps, ctr)

cdata = zeros(length(ctr), 1)
for i in eachindex(ctr)
    cdata[i, 1] = ctr[i].w
end
write_vtk(ps.points, ps.cellid, cdata)


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

abstract type AbstractPenalization end
abstract type AbstractPenalizer end
# Penalization API #
# ------------------------------------------------------------
penalty(ph::AbstractProgressiveHedging) = penalty(ph, ph.penalization)
initialize_penalty!(ph::AbstractProgressiveHedging) = initialize_penalty!(ph, ph.penalization)
update_penalty!(ph::AbstractProgressiveHedging) = update_penalty!(ph, ph.penalization)
# ------------------------------------------------------------
# Attributes #
# ------------------------------------------------------------
"""
    RawPenalizationParameter

An optimizer attribute used for raw parameters of the penalizer. Defers to `RawParameter`.
"""
struct RawPenalizationParameter <: PenalizationParameter
    name::Any
end
include("common.jl")
include("fixed.jl")
include("adaptive.jl")


text:
# the code in here is intend to help make the move lookups table

function rowvec2uint16(col::Vector{Int})
    row = zero(UInt16)
    row |= UInt8(col[1])
    for i in 2:4
        row <<= 4
        row |= UInt8(col[i])
    end
    row
end

# u
function make_row_lookup()
    left_lookup = Vector{UInt16}(undef, 2^16)
    left_reward_lookup = Vector{Int32}(undef, 2^16)
    right_lookup = Vector{UInt16}(undef, 2^16)
    right_reward_lookup = Vector{Int32}(undef, 2^16)

    vecs = vec([[v...] for v in Iterators.product(0:15, 0:15, 0:15, 0:15)])
    left_lookup_arr = deepcopy(vecs)
    right_lookup_arr = deepcopy(vecs)

    for (i, v) in enumerate(vecs)
        idx = rowvec2uint16(v)
        _, reward, _ = move_up!(left_lookup_arr[i])
        left_reward_lookup[idx+1] = reward
        left_lookup[idx+1] = rowvec2uint16(left_lookup_arr[i])

        rla = right_lookup_arr[i]
        _, reward, _ = move_up!(@view rla[4:-1:1])
        right_reward_lookup[idx+1] = reward
        right_lookup[idx+1] = rowvec2uint16(right_lookup_arr[i])
    end

    left_lookup, left_reward_lookup, right_lookup, right_reward_lookup
end


text:
using ACTRSimulators, Test, ACTRModels, Random, DataFrames
import ACTRSimulators: start!, press_key!
Random.seed!(8985)
include("task.jl")

scheduler = ACTRScheduler(;model_trace=true, store=true)
task = SimpleTask(;scheduler)
procedural = Procedural()
T = vo_to_chunk() |> typeof
visual_location = VisualLocation(buffer=T[])
visual = Visual(buffer=T[])
visicon = VisualObject[]
motor = Motor()
memory = [Chunk(;animal=:dog), Chunk(;animal=:cat)]
declarative = Declarative(;memory)
actr = ACTR(;scheduler, procedural, visual_location, visual, motor, declarative, visicon)

function can_retrieve(actr)
    c1(actr) = !actr.declarative.state.busy
    c2(actr) = actr.declarative.state.empty
    return c1, c2
end

function can_stop(actr)
    c1(actr) = !actr.declarative.state.empty
    return (c1,)
end

function retrieve_chunk(actr, task)
    retrieving!(actr; animal=:dog)
end

function stop(actr, task)
    stop!(actr.scheduler)
end

rule1 = Rule(;conditions=can_retrieve, action=retrieve_chunk, actr, task, name="Retrieve")
push!(procedural.rules, rule1)

rule2 = Rule(;conditions=can_stop, action=stop, actr, task, name="Stop")
push!(procedural.rules, rule2)
run!(actr, task)
chunk = actr.declarative.buffer[1]
@test chunk.slots == (animal=:dog,)

observed = map(x->x.description, scheduler.complete_events)
expected = [
    "Starting",
    "Selected Retrieve",
    "Retrieve",
    "Selected Stop"
]
@test expected == observed


text:
if VERSION < v"0.7"
    testdir = joinpath(Pkg.dir("MathProgBase"), "test")
else
    import MathProgBase
    testdir = joinpath(dirname(pathof(MathProgBase)), "..", "test")
end
include(joinpath(testdir, "linprog.jl"))
linprogtest(ClpSolver())

include(joinpath(testdir, "linproginterface.jl"))
linprogsolvertest(ClpSolver())

include(joinpath(testdir, "conicinterface.jl"))
coniclineartest(ClpSolver())


text:
module ContingencyTest

using Test
using ReactiveMP
using Distributions
using Random

@testset "Contingency" begin

    @testset "common" begin
        @test Contingency <: Distribution
        @test Contingency <: ContinuousDistribution
        @test Contingency <: MatrixDistribution

        @test value_support(Contingency) === Continuous
        @test variate_form(Contingency)  === Matrixvariate
    end

    @testset "contingency_matrix" begin
        @test ReactiveMP.contingency_matrix(Contingency(ones(3, 3))) == ones(3, 3)
        @test ReactiveMP.contingency_matrix(Contingency(ones(4, 4))) == ones(4, 4)
    end

    @testset "vague" begin
        @test_throws MethodError vague(Contingency)

        d1 = vague(Contingency, 3)

        @test typeof(d1) <: Contingency
        @test ReactiveMP.contingency_matrix(d1) ≈ ones(3, 3) ./ 9

        d2 = vague(Contingency, 4)

        @test typeof(d2) <: Contingency
        @test ReactiveMP.contingency_matrix(d2) ≈ ones(4, 4) ./ 16
    end

    @testset "entropy" begin
        @test entropy(Contingency([ 0.1 0.9; 0.9 0.1 ]))     ≈ 0.6501659467828964
        @test entropy(Contingency([ 0.2 0.8; 0.8 0.2 ]))     ≈ 1.0008048470763757
        @test entropy(Contingency([ 0.45 0.75; 0.55 0.25 ])) ≈ 1.2504739583323967
    end

end

end


text:
using Test
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using CutCellDG
include("../useful_routines.jl")

function exact_displacement(x, e11, e22)
    return [x[1] * e11, x[2] * e22]
end

function test_merged_simple_tension()
    polyorder = 1
    numqp = 2

    widths = [2.0, 1.0]
    nelmts = [5, 5]

    elasticitybasis = LagrangeTensorProductBasis(2, polyorder)
    levelsetbasis = HermiteTensorProductBasis(2)

    dim, nf = size(interpolation_points(levelsetbasis))
    refpoints = interpolation_points(elasticitybasis)

    cgmesh = CutCellDG.CGMesh([0.0, 0.0], widths, nelmts, nf)
    mesh = CutCellDG.DGMesh([0.0, 0.0], widths, nelmts, refpoints)

    interfaceangle = 40.0
    normal = [cosd(interfaceangle), sind(interfaceangle)]
    x0 = [1.8, 0.0]
    levelset = CutCellDG.LevelSet(
        x -> plane_distance_function(x, normal, x0),
        cgmesh,
        levelsetbasis,
    )

    cutmesh = CutCellDG.CutMesh(mesh, levelset)
    cellquads = CutCellDG.CellQuadratures(cutmesh, levelset, numqp)
    interfacequads = CutCellDG.InterfaceQuadratures(cutmesh, levelset, numqp)
    facequads = CutCellDG.FaceQuadratures(cutmesh, levelset, numqp)

    mergedmesh =
        CutCellDG.MergedMesh!(cutmesh, cellquads, facequads, interfacequads)
    @assert CutCellDG.has_merged_cells(mergedmesh)

    lambda, mu = (1.0, 2.0)
    penalty = 1.0
    eta = 1
    dx = 0.1
    e11 = dx / 2.0
    e22 = -lambda / (lambda + 2mu) * e11
    dy = e22
    s11 = (lambda + 2mu) * e11 + lambda * e22

    stiffness = CutCellDG.HookeStiffness(lambda, mu, lambda, mu)

    sysmatrix = CutCellDG.SystemMatrix()
    sysrhs = CutCellDG.SystemRHS()

    CutCellDG.assemble_displacement_bilinear_forms!(
        sysmatrix,
        elasticitybasis,
        cellquads,
        stiffness,
        mergedmesh,
    )
    CutCellDG.assemble_interelement_condition!(
        sysmatrix,
        elasticitybasis,
        facequads,
        stiffness,
        mergedmesh,
        penalty,
        eta,
    )
    CutCellDG.assemble_coherent_interface_condition!(
        sysmatrix,
        elasticitybasis,
        interfacequads,
        stiffness,
        mergedmesh,
        penalty,
        eta,
    )
    CutCellDG.assemble_penalty_displacement_component_bc!(
        sysmatrix,
        sysrhs,
        x -> 0.0,
        elasticitybasis,
        facequads,
        stiffness,
        mergedmesh,
        x -> x[1] ≈ 0.0,
        [1.0, 0.0],
        penalty,
    )
    CutCellDG.assemble_penalty_displacement_component_bc!(
        sysmatrix,
        sysrhs,
        x -> 0.0,
        elasticitybasis,
        facequads,
        stiffness,
        mergedmesh,
        x -> x[2] ≈ 0.0,
        [0.0, 1.0],
        penalty,
    )
    CutCellDG.assemble_penalty_displacement_component_bc!(
        sysmatrix,
        sysrhs,
        x -> dx,
        elasticitybasis,
        facequads,
        stiffness,
        mergedmesh,
        x -> x[1] ≈ widths[1],
        [1.0, 0.0],
        penalty,
    )

    matrix = CutCellDG.sparse_displacement_operator(sysmatrix, mergedmesh)
    rhs = CutCellDG.displacement_rhs_vector(sysrhs, mergedmesh)

    solution = matrix \ rhs
    displacement = reshape(solution, 2, :)

    err = mesh_L2_error(
        displacement,
        x -> exact_displacement(x, e11, e22),
        elasticitybasis,
        cellquads,
        mergedmesh,
    )

    @test allapprox(err, zeros(2), 1e2eps())
end

test_merged_simple_tension()


text:
function normout(A::SparseMatrixCSC{T,Int64}) where T
  # NORMOUT Normalize the outdegrees of the matrix A.
  #
  # P = normout(A)
  #
  #   P has the same non-zero structure as A, but is normalized such that the
  #   sum of each row is 1, assuming that A has non-negative entries.
  #

  # compute the row-sums/degrees
  #   d = sum(A,2)
  #   d = squeeze(d',1)
  #   id = 1./d
  #   P = spdiagm(id)*A
  #   return P

    S = vec(sum(A,dims=2))
    ei,ej,ev = findnz(A)
    m,n = size(A)
    P = sparse(ei,ej,ev./S[ei],m,n)
    return P

end


text:
# # plotting functions ------------------------------------------------------------
function plot_omegahat(chain::Chains, s::NamedTuple; prior_only = false, legend = :bottomleft)
    w, w1, w2 = summarize_omega(chain)
    p1 = plot(s.timesteps, w.μ', ribbon = [w.lb', w.ub'],
                fillalpha = 0.3, linewidth = 2, color = :blue,
                title = "Total velocity",
                label = "Estimated",
                xlabel = "Time (s)",
                ylabel = "Angular velocity (deg/s)")

    plot!(p1, s.timesteps, s.ω,
            linestyle = :dash,
            linewidth = 2,
            color = :cornflowerblue,
            label = "True")
    if !prior_only
        scatter!(p1, s.timesteps, s.y,
            label = "Observations",
            markershape = :circle,
            markersize = 2,
            markerstrokealpha = 0.8,
            markeralpha = 0.2,
            markercolor = :black)
    end

    p2 = plot(s.timesteps, w1.μ', ribbon = [w1.lb', w1.ub'],
            fillalpha = 0.3, linewidth = 2, color = :darkorange,
            title = "Active",
            label = "Estimated")
    plot!(p2, s.timesteps, s.ωᵤ,
            linestyle = :dash,
            linewidth = 2,
            color = :orangered,
            label = "True")

    p3 = plot(s.timesteps, w2.μ', ribbon = [w2.lb', w2.ub'],
            fillalpha = 0.3, linewidth = 2, color = :palegreen,
            title = "Passive",
            label = "Estimated",
            xlabel = "Time (s)")
    plot!(p3, s.timesteps, s.ωₑ,
            linestyle = :dash,
            linewidth = 2,
            color = :darkgreen,
            label = "True")

    l = @layout [a [b; c]]
    p = plot(p1, p2, p3, layout = l,
            link=:all,
            titleloc = :right,
            legend = legend,
            foreground_color_legend = nothing )
    return p
end

function plot_omegahat_simple(chain::Chains, s::NamedTuple;
                                prior_only = false,
                                spaghetti = false,
                                legend = :bottomleft)
    if spaghetti === true
        out = Array(group(chain, :ω))
        nrows, ncols = size(out)
        idx = sample(1:nrows, 100)
        samples = out[idx, :]

        p1 = plot(s.timesteps, samples',
                        linewidth = 2,
                    # color = :blue,
                    alpha = 0.1, color = :steelblue,
                    title = "Total velocity",
                    legend = legend,
                    label = "Estimated",
                    xlabel = "Time (s)",
                    ylabel = "Angular velocity (deg/s)")
    else
        w_array = Array(group(chain, :ω))
        w = compute_stats(w_array)
        p1 = plot(s.timesteps, w.μ', ribbon = [w.lb', w.ub'],
                fillalpha = 0.3,
                linewidth = 2,
                color = :blue,
                title = "Total velocity",
                legend = legend,
                label = "Estimated",
                xlabel = "Time (s)",
                ylabel = "Angular velocity (deg/s)")
    end

    plot!(p1, s.timesteps, s.ω,
            linestyle = :dash,
            linewidth = 2,
            color = :cornflowerblue,
            label = "True")
    if !prior_only
        scatter!(p1, s.timesteps, s.y,
            label = "Observations",
            markershape = :circle,
            markersize = 2,
            markerstrokealpha = 0.8,
            markeralpha = 0.2,
            markercolor = :black)
    end
    return p1
end

function plot_omegatrue(s::NamedTuple;
                        displacement = false,
                        legend = :bottomleft)
    plot(s.timesteps, [s.ωᵤ, s.ωₑ], label = ["self" "external"],
        linewidth = 3,
        linestyle = :dot,
        linecolor = ["red" "steelblue"],
        legend = legend,
        foreground_color_legend = nothing)

    plot!(s.timesteps, s.ω,
        label = "total velocity",
        lw = 3,
        color = "purple")

    if displacement===true
        plot!(s.timesteps, s.θ,
            label = "total displacement",
            lw = 3,
            color = "orange")
    end

    plot!(s.timesteps, s.ω - s.c,
        label = "lag (cupula dynamics)",
        lw = 3,
        color = :grey)

    scatter!(s.timesteps, s.y,
        xlabel = "Time (s)",
        ylabel = "Angular velocity (deg/s)",
        label = "observations",
        markershape = :circle,
        markersize = 3,
        markeralpha = 0.5,
        markercolor = :black)
end

function plot_omegahat_simulation(chain::Chains, s::NamedTuple; prior_only = false, legend = :topright)
    w_array = Array(group(chain, :ω))
    w = compute_stats(w_array)
    p1 = plot(s.timesteps, w.μ', ribbon = [w.lb', w.ub'],
                fillalpha = 0.3,
                linewidth = 2,
                color = :blue,
                title = "Total velocity",
                legend = legend,
                label = "Estimated",
                xlabel = "Time (s)",
                ylabel = "Angular velocity (deg/s)")

    if !prior_only
        scatter!(p1, s.timesteps, s.y,
            label = "Observations",
            markershape = :circle,
            markersize = 2,
            markerstrokealpha = 0.8,
            markeralpha = 0.2,
            markercolor = :black)
    end
    return p1
end

function plot_nigmatullina_trial(motion)
    p = scatter(motion.timesteps, motion.y,
    xlabel = "Time (s)",
    ylabel = "Angular velocity (deg/s)",
    legend = :none,
    markershape = :circle,
    markersize = 2,
    markeralpha = 0.1,
    markercolor = :black)
    plot!(motion.timesteps, motion.α,
        linewidth = 4,
        linecolor = :steelblue)
    vline!([motion.onset], linewidth = 2, linestyle = :dot, linecolor = "grey")
    plot!(p, motion.timesteps, motion.ω,
        linewidth = 4,
        linecolor = :orange)
end


text:
push!(LOAD_PATH, "/home/elsuizo/Robotics/src/")
ENV["JULIA_EDITOR"] = "/home/elsuizo/neovim/build/bin/./nvim"


text:
using Documenter
using Graphics

DocMeta.setdocmeta!(Graphics, :DocTestSetup, :(using Graphics); recursive=true)

makedocs(
    sitename = "Graphics",
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    modules = [Graphics],
    pages = ["index.md", "reference.md"],
    clean = false,
    checkdocs = :exports
)

deploydocs(
    repo = "github.com/JuliaGraphics/Graphics.jl.git",
    push_preview = true
)


text:
using Zygote, Tracker

### Zygote

Zygote.@nograd rsimilar
Zygote.@nograd _rand_nograd

### Tracker

Base.eltype(x::Tracker.TrackedArray) = eltype(Tracker.data(x))
rsimilar(rng, f!, x::Tracker.TrackedArray, n) = rsimilar(rng, f!, Tracker.data(x), n)

function u2gumbelback(u, g)
    return 1 ./ exp.(-g) ./ (u .+ eps(u))
end

u2gumbel(u::Tracker.TrackedArray) = Tracker.track(u2gumbel, u)

Tracker.@grad function u2gumbel(u)
    g = u2gumbel(Tracker.data(u))
    return g, Δ -> (Δ .* u2gumbelback(u, g),)
end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module TestVariableRSOCtoPSD

using Test

using MathOptInterface
const MOI = MathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

include("../utilities.jl")

# Can come from a SOC of dimension 2 which makes more sense
# FIXME should it be moved to contconic or is RSOC of dimension 2 too exotic ?
function test_RSOC_of_dimension_2()
    mock = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
    bridged_mock = MOI.Bridges.Variable.RSOCtoPSD{Float64}(mock)
    xy, cxy = MOI.add_constrained_variables(
        bridged_mock,
        MOI.RotatedSecondOrderCone(2),
    )
    x, y = xy
    c = MOI.add_constraint(bridged_mock, 1.0x + 1.0y, MOI.LessThan(1.0))
    obj = 1.0y
    MOI.set(bridged_mock, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(bridged_mock, MOI.ObjectiveFunction{typeof(obj)}(), obj)
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [0.0, 2.0],
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) =>
                [-1.0],
        )
    MOI.optimize!(bridged_mock)
    @test MOI.get(bridged_mock, MOI.ObjectiveValue()) == 1.0
    @test MOI.get(bridged_mock, MOI.DualObjectiveValue()) == 1.0
    @test MOI.get(bridged_mock, MOI.VariablePrimal(), xy) == [0.0, 1.0]
    @test MOI.get(bridged_mock, MOI.ConstraintPrimal(), cxy) == [0.0, 1.0]
    @test MOI.get(bridged_mock, MOI.ConstraintDual(), cxy) == [1.0, 0.0]
    @test MOI.get(bridged_mock, MOI.ConstraintPrimal(), c) == 1.0
    @test MOI.get(bridged_mock, MOI.ConstraintDual(), c) == -1.0
    MOI.set(mock, MOI.ConstraintName(), c, "c")
    var_names = ["a", "b"]
    MOI.set(
        mock,
        MOI.VariableName(),
        MOI.get(mock, MOI.ListOfVariableIndices()),
        var_names,
    )
    nonneg = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{MOI.VectorOfVariables,MOI.Nonnegatives}(),
    )
    @test length(nonneg) == 1
    MOI.set(mock, MOI.ConstraintName(), nonneg[1], "cab")
    s = """
    variables: a, b
    cab: [a, b] in MathOptInterface.Nonnegatives(2)
    c: a + 0.5b <= 1.0
    maxobjective: 0.5b
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(mock, model, var_names, ["cab", "c"])
    var_names = ["x", "y"]
    MOI.set(
        bridged_mock,
        MOI.VariableName(),
        MOI.get(bridged_mock, MOI.ListOfVariableIndices()),
        var_names,
    )
    MOI.set(bridged_mock, MOI.ConstraintName(), cxy, "cxy")
    s = """
    variables: x, y
    cxy: [x, y] in MathOptInterface.RotatedSecondOrderCone(2)
    c: x + y <= 1.0
    maxobjective: 1.0y
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        bridged_mock,
        model,
        var_names,
        ["cxy", "c"],
    )
    _test_delete_bridged_variables(
        bridged_mock,
        xy,
        MOI.RotatedSecondOrderCone,
        2,
        (
            (MOI.VectorOfVariables, MOI.Nonnegatives, 0),
            (MOI.VectorOfVariables, MOI.PositiveSemidefiniteConeTriangle, 0),
            (MOI.VariableIndex, MOI.EqualTo{Float64}, 0),
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}, 0),
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}, 1),
        ),
    )
    return
end

function test_RSOC4()
    mock = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
    bridged_mock = MOI.Bridges.Variable.RSOCtoPSD{Float64}(mock)
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [1.0, 1.0, 2.0, 1.0, 0.0, 2.0],
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) =>
                [0.25],
            (MOI.VariableIndex, MOI.EqualTo{Float64}) => [-0.5],
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) =>
                [-1.0],
            (MOI.VectorOfVariables, MOI.PositiveSemidefiniteConeTriangle) =>
                [[1.0, -0.5, 0.25, -0.5, 0.25, 0.25]],
        )
    mock.eval_variable_constraint_dual = false
    MOI.Test.test_conic_RotatedSecondOrderCone_out_of_order(
        bridged_mock,
        MOI.Test.Config(),
    )
    mock.eval_variable_constraint_dual = true
    var_names = ["Q$i$j" for j in 1:3 for i in 1:j]
    MOI.set(
        mock,
        MOI.VariableName(),
        MOI.get(mock, MOI.ListOfVariableIndices()),
        var_names,
    )
    psd = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorOfVariables,
            MOI.PositiveSemidefiniteConeTriangle,
        }(),
    )
    @test length(psd) == 1
    MOI.set(mock, MOI.ConstraintName(), psd[1], "psd")
    off_diag = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.EqualTo{Float64}}(),
    )
    @test length(off_diag) == 1
    diag = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{
            MOI.ScalarAffineFunction{Float64},
            MOI.EqualTo{Float64},
        }(),
    )
    @test length(diag) == 1
    MOI.set(mock, MOI.ConstraintName(), diag[1], "diag33")
    c = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{
            MOI.ScalarAffineFunction{Float64},
            MOI.LessThan{Float64},
        }(),
    )
    @test length(c) == 1
    MOI.set(mock, MOI.ConstraintName(), c[1], "c")

    s = """
    variables: Q11, Q12, Q13, Q22, Q23, Q33
    psd: [Q11, Q12, Q22, Q13, Q23, Q33] in MathOptInterface.PositiveSemidefiniteConeTriangle(3)
    Q23 == 0.0
    diag33: Q22 + -1.0Q33 == 0.0
    c: Q11 + 0.5Q22 <= 2.0
    maxobjective: Q12 + Q13
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        mock,
        model,
        var_names,
        ["psd", "diag33", "c"],
        [("Q23", MOI.EqualTo{Float64}(0.0))],
    )

    var_names = ["t", "u", "x", "y"]
    MOI.set(
        bridged_mock,
        MOI.VariableName(),
        MOI.get(bridged_mock, MOI.ListOfVariableIndices()),
        var_names,
    )
    rsoc = MOI.get(
        bridged_mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorOfVariables,
            MOI.RotatedSecondOrderCone,
        }(),
    )
    @test length(rsoc) == 1
    MOI.set(bridged_mock, MOI.ConstraintName(), rsoc[1], "rsoc")
    c = MOI.get(
        bridged_mock,
        MOI.ListOfConstraintIndices{
            MOI.ScalarAffineFunction{Float64},
            MOI.LessThan{Float64},
        }(),
    )
    @test length(c) == 1
    MOI.set(bridged_mock, MOI.ConstraintName(), c[1], "c")

    s = """
    variables: t, u, x, y
    rsoc: [t, u, x, y] in MathOptInterface.RotatedSecondOrderCone(4)
    c: t + u <= 2.0
    maxobjective: x + y
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        bridged_mock,
        model,
        var_names,
        ["rsoc", "c"],
    )

    v = MOI.get(bridged_mock, MOI.ListOfVariableIndices())
    @test length(v) == 4

    message = string(
        "Cannot delete variable as it is constrained with other",
        " variables in a `MOI.VectorOfVariables`.",
    )
    for i in 1:4
        err = MOI.DeleteNotAllowed(v[i], message)
        @test_throws err MOI.delete(bridged_mock, v[i])
    end

    _test_delete_bridged_variables(
        bridged_mock,
        v,
        MOI.RotatedSecondOrderCone,
        4,
        (
            (MOI.VectorOfVariables, MOI.Nonnegatives, 0),
            (MOI.VectorOfVariables, MOI.PositiveSemidefiniteConeTriangle, 0),
            (MOI.VariableIndex, MOI.EqualTo{Float64}, 0),
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}, 0),
        ),
    )
    return
end

end  # module

TestVariableRSOCtoPSD.runtests()

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
const Z = UInt8('0')

ndigits_decimal(x::Signed) =
    ndigits_decimal(unsigned(abs(x)))

function ndigits_decimal(x::Unsigned)
    n = 0
    while true
        if x < 10
            n += 1
            break
        elseif x < 100
            n += 2
            break
        elseif x < 1000
            n += 3
            break
        else
            n += 4
            x < 10000 && break
            x = div(x, oftype(x, 10000))
        end
    end
    return n
end

ndigits_decimal(x::BigInt) = ndigits(x, base = 10)

# inlined copy of static data
macro copy(dst, p, src::String)
    block = Expr(:block)
    n = ncodeunits(src)
    for i in 1:n
        push!(block.args, :($dst[$p+$i-1] = $(codeunit(src, i))))
    end
    push!(block.args, :(p + $n))
    return esc(block)
end

# write sign mark
function sign(data, p, x, sgn)
    signed = false
    if sgn == SIGN_NONE
        # pass
    elseif signbit(x)
        data[p] = UInt8('-')
        signed = true
    elseif sgn == SIGN_PLUS
        data[p] = UInt8('+')
        signed = true
    elseif sgn == SIGN_SPACE
        data[p] = UInt8(' ')
        signed = true
    end
    return p + signed, signed
end

# calculate the number of leading zeros for padding
function number_of_leading_zeros(m, k, minwidth)
    n = div(m - 1, k)  # number of separators
    m + n < minwidth || return 0
    width = minwidth + (rem(minwidth, k + 1) == 0)
    ndigits = width - div(width - 1, k + 1)
    return ndigits - m
end

function insert_zeros(data, start, p, n)
    n > 0 || return p
    copyto!(data, start + n, data, start, p - start)
    for i in 1:n
        data[start+i-1] = Z
    end
    return p + n
end

# post-hoc alignment
function aligncontent(data, p, start, width::Int, fill::Char, align::Alignment, minwidth::Int)
    pw = max(minwidth - width, 0)  # padding width
    ps = pw * ncodeunits(fill)     # padding size
    size = p - start
    if align == ALIGN_RIGHT
        copyto!(data, start + ps, data, start, size)
        pad(data, start, fill, pw)
        p += ps
    elseif align == ALIGN_CENTER
        offset = ncodeunits(fill) * (pw ÷ 2)
        copyto!(data, start + offset, data, start, size)
        pad(data, start, fill, pw ÷ 2)
        pad(data, start + offset + size, fill, pw - pw ÷ 2)
        p += ps
    else
        p = pad(data, p, fill, pw)
    end
    return p
end

@inline function paddingwidth(s::Spec, width::Int)
    @assert s.width isa Int || s.width isa Nothing
    return s.width isa Int ? max(s.width - width, 0) : 0
end

paddingsize(s::Spec, width::Int) =
    s.fill === nothing ? 0 : paddingwidth(s, width) * ncodeunits(s.fill)

@inline function padleft(data::Vector{UInt8}, p::Int, fill::Char, align::Alignment, pw::Int)
    @assert isspecified(align)
    if align == ALIGN_RIGHT
        p = pad(data, p, fill, pw)
    elseif align == ALIGN_CENTER
        p = pad(data, p, fill, pw ÷ 2)
    end
    return p
end

@inline function padright(data::Vector{UInt8}, p::Int,fill::Char, align::Alignment, pw::Int)
    @assert isspecified(align)
    if align == ALIGN_LEFT
        p = pad(data, p, fill, pw)
    elseif align == ALIGN_CENTER
        p = pad(data, p, fill, pw - pw ÷ 2)
    end
    return p
end

function pad(data::Vector{UInt8}, p::Int, fill::Char, w::Int)
    for _ in 1:w
        p = char(data, p, fill)
    end
    return p
end

function char(data::Vector{UInt8}, p::Int, char::Char)
    m = ncodeunits(char)
    x = reinterpret(UInt32, char) >> 8(4 - m)
    if m == 1
        data[p]   =  x        % UInt8
    elseif m == 2
        data[p+1] =  x        % UInt8
        data[p]   = (x >> 8)  % UInt8
    elseif m == 3
        data[p+2] =  x        % UInt8
        data[p+1] = (x >>  8) % UInt8
        data[p]   = (x >> 16) % UInt8
    else
        @assert m == 4
        data[p+3] =  x        % UInt8
        data[p+2] = (x >>  8) % UInt8
        data[p+1] = (x >> 16) % UInt8
        data[p]   = (x >> 24) % UInt8
    end
    return p + m
end


text:
using DelimitedFiles, Unicode

"""

# read_samples

Read sample output files created by cmdstan.

### Method
```julia
read_samples(m::Stanmodel)
```

### Required arguments
```julia
* `m::Stanmodel`    : Stanmodel object
```

"""
function read_samples(m::Stanmodel, diagnostics=false, warmup_samples=false)

  local monitors, index, idx, indvec, ftype, noofsamples

  # a3d will contain the samples such that a3d[s, i, c] where

  #   s:  the sample index (corrected warmup samples and thinning by cmdstan)
  #   i:  variable name (either from monitors or read from .csv file produiced by cmdstan)
  #   c:  chain number (from m.chains)

  if isa(m.method, Sample)
    ftype = diagnostics ? "diagnostics" : "samples"

    # The sample index s runs from 1:noofsamples

    if m.method.save_warmup
      noofsamples = floor(Int, (m.method.num_samples+m.method.num_warmup)/m.method.thin)
    else
      noofsamples = floor(Int, m.method.num_samples/m.method.thin)
    end
  else
    ftype = lowercase(string(typeof(m.method)))
    if (ftype[1:7] == "cmdstan")
      ftype = ftype[9:end]
    end

    # The noofsamples is obtained from method.output_samples

    noofsamples = m.method.output_samples
  end

  # Read .csv files created by each chain

  for i in 1:m.nchains
    file_path = joinpath(m.tmpdir, "$(m.name)_$(ftype)_$(i).csv")
    if isfile(file_path)

      instream = open(file_path)

      # Skip initial set of commented lines, e.g. containing
      # cmdstan version info, etc.

      skipchars(isspace, instream, linecomment='#')

      # First non-comment line contains names of variables

      line = Unicode.normalize(readline(instream), newline2lf=true)
      idx = split(strip(line), ",")
      index = [idx[k] for k in 1:length(idx)]

      # Only continue with monitored variables

      if length(m.monitors) == 0
        indvec = 1:length(index)
      else
        indvec = findall((in)(m.monitors), index)
      end

      # Preserve cnames and create a3d

      if i == 1
        global cnames = convert.(String, idx[indvec])
        global a3d = fill(0.0, noofsamples, length(indvec), m.nchains)
      end

      # Read in the samples for all chains

      skipchars(isspace, instream, linecomment='#')
      for j in 1:noofsamples
        skipchars(isspace, instream, linecomment='#')
        line = Unicode.normalize(readline(instream), newline2lf=true)
        if eof(instream) && length(line) < 2
          return
        else
          flds = parse.(Float64, split(strip(line), ","))
          flds = reshape(flds[indvec], 1, length(indvec))
          a3d[j,:,i] = flds
        end
      end   # read in samples
    end   # select next file if it exists
  end   # loop over all chains

  (a3d, cnames)

end   # end of read_samples


text:
using NeuralNetwork
import NeuralNetwork: accuracy

using RDatasets: dataset

iris = dataset("datasets", "iris")

X = Array(iris[:, 1:4])
y = zeros(Int, size(iris, 1))
for (i, label) in enumerate(["versicolor", "setosa", "virginica"])
    y[iris[:, 5] .== label] = i
end

indices = shuffle(1:size(iris, 1))

X = X[indices, :]'
y = y[indices]

N = Int(round(size(X, 2)*0.8))
X_train, X_test = X[:, 1:N], X[:, N+1:end]
y_train, y_test = y[1:N], y[N+1:end]

layer = LogisticRegression(4, 3, 5e-3, 1000)
fit!(layer, X_train, onehot(y_train))
Ŷ = predict(layer, X_test)

ŷ = maxindices(Ŷ)
println("Accuracy: $(accuracy(y_test, ŷ))")

max_stars_repo_name: IshitaTakeshi/NeuralNetwork

max_issues_repo_name: IshitaTakeshi/NeuralNetwork

max_forks_repo_name: IshitaTakeshi/NeuralNetwork


text:
#Multi-threaded:
#julia -t <number-of-threads> --project=@. main.jl
#Single thread:
#julia --project=@. main.jl
using NEOs
using Dates
using TaylorIntegration
using JLD
using PlanetaryEphemeris
@show Threads.nthreads()

#script parameters (TODO: use ArgParse.jl instead)
varorder = 5 # varorder is the order corresponding to the jet transport perturbation
nv = 8 #number of TaylorN variables
objname = "Apophis"
maxsteps = 10000
nyears_bwd = -18.0 # years in backward integration
nyears_fwd = 9.0 # years in forward integration
dense = false #true
quadmath = false # use quadruple precision
lyap = false #true # compute Lyapunov exponents
neosjlpath = pkgdir(NEOs)
dynamics = RNp1BP_pN_A_J23E_J2S_ng_eph_threads!

# integration parameters
order = 25
abstol = 1.0E-20

### observation data files (ra/dec, del/dop)
optfile_bwd = joinpath(neosjlpath, "data", "99942_2004_2020.dat")
optfile_fwd = joinpath(neosjlpath, "data", "99942_2020_2021.dat")
radarfile_bwd = joinpath(neosjlpath, "data", "99942_RADAR_2005_2013.dat")
radarfile_fwd = joinpath(neosjlpath, "data", "99942_RADAR_2021.dat")

### path to local Solar System ephemeris file
ss_eph_file = "./sseph343ast016_p31y_et.jld"
ss16asteph_et = JLD.load(ss_eph_file, "ss16ast_eph")

### TaylorN variables setup
if lyap
    ### setup TaylorN variables with order=1, numvars=nv
    TNvars = set_variables("δx", order=1, numvars=nv)
    # dq corresponding to solution OR7 for Apophis
    dq_OR7 = [-8.053250543083672e-7, -6.4934453239292154e-9, -3.552581604396334e-8, 2.382431039885935e-9, -1.3384789262277344e-8, -4.6746457798167725e-9, -2.892614243659006e-14, 0.0]
    dq = dq_OR7
else
    #### dq: perturbation to nominal initial condition (Taylor1 jet transport)
    #dq = Taylor1.(zeros(7), varorder)
    #dq[end][1] = 1e-14
    # dq: perturbation to nominal initial condition (TaylorN jet transport)
    dq = set_variables("δx", order=varorder, numvars=nv)
    for i in 1:6
        dq[i][1][i] = 1e-8
    end
    if get_numvars() == 8
        dq[7][1][7] = 1e-14
        dq[8][1][8] = 1e-13
    end
end

### initial conditions from Apophis JPL solution #197 at the Dec-17-2020.0 (TDB) epoch
jd0 = datetime2julian(DateTime(2020,12,17)) # JDTDB = 2459200.5
q00 = [-0.18034747703273316, 0.9406910666200128, 0.3457360259054398, -0.016265942170279046, 4.392889725556651e-5, -0.00039519931615139716] ### JPL solution #197 at 2020Dec17.0 (TDB)
q0 = vcat(q00, 0.0, 0.0) .+ dq

####integrator warmup
propagate(objname, dynamics, 1, jd0, nyears_fwd, ss16asteph_et, output=false, dense=dense, q0=q0, quadmath=quadmath, lyap=lyap, order=order, abstol=abstol)
println("*** Finished warmup")

######Full jet transport integration
propagate(objname*"_bwd", dynamics, maxsteps, jd0, nyears_bwd, ss16asteph_et, dense=dense, q0=q0, quadmath=quadmath, lyap=lyap, order=order, abstol=abstol, radarobsfile=radarfile_bwd, opticalobsfile=optfile_bwd, tord=10, niter=5)
propagate(objname*"_fwd", dynamics, maxsteps, jd0, nyears_fwd, ss16asteph_et, dense=dense, q0=q0, quadmath=quadmath, lyap=lyap, order=order, abstol=abstol, radarobsfile=radarfile_fwd, opticalobsfile=optfile_fwd, tord=10, niter=5)
println("*** Finished asteroid ephemeris integration")


text:
using ModelingToolkitStandardLibrary.Mechanical.Rotational, ModelingToolkit, OrdinaryDiffEq, Test
import ModelingToolkitStandardLibrary.Blocks
# using Plots

@parameters t
D = Differential(t)

@testset "two inertias" begin
    @named fixed = Fixed()
    @named inertia1 = Inertia(J=2) # this one is fixed
    @named spring = Spring(c=1e4)
    @named damper = Damper(d=10)
    @named inertia2 = Inertia(J=2, phi_start=pi/2)

    connections = [
        connect(fixed.flange, inertia1.flange_b)
        connect(inertia1.flange_b, spring.flange_a, damper.flange_a)
        connect(spring.flange_b, damper.flange_b, inertia2.flange_a)
    ]

    @named model = ODESystem(connections, t, systems=[fixed, inertia1, inertia2, spring, damper])
    sys = structural_simplify(model)
    prob = DAEProblem(sys, D.(states(sys)) .=> 0.0, [], (0, 10.0))
    sol = solve(prob, DFBDF())

    # Plots.plot(sol; vars=[inertia1.w, inertia2.w])

    @test sol.retcode == :Success
    @test_skip all(sol[inertia1.w] .== 0)
    @test sol[inertia2.w][end] ≈ 0 atol=1e-3 # all energy has dissipated
end

@testset "two inertias with driving torque" begin
    amplitude = 10 # Amplitude of driving torque
    frequency = 5 # Frequency of driving torque
    J_motor = 0.1 # Motor inertia

    @named fixed = Fixed()
    @named torque = Torque(use_support=true)
    @named inertia1 = Inertia(J=2, phi_start=pi/2)
    @named spring = Spring(c=1e4)
    @named damper = Damper(d=10)
    @named inertia2 = Inertia(J=4)
    @named sine = Blocks.Sine(amplitude=amplitude, frequency=frequency)

    connections = [
        connect(sine.output, torque.tau)
        connect(torque.support, fixed.flange)
        connect(torque.flange, inertia1.flange_a)
        connect(inertia1.flange_b, spring.flange_a, damper.flange_a)
        connect(spring.flange_b, damper.flange_b, inertia2.flange_a)
    ]

    @named model = ODESystem(connections, t, systems=[fixed, torque, inertia1, inertia2, spring, damper, sine])
    sys = structural_simplify(model)
    #prob = ODAEProblem(sys, Pair[], (0, 1.0))
    prob = DAEProblem(sys, D.(states(sys)) .=> 0.0, [D(D(inertia2.phi)) => 1.0], (0, 10.0))
    sol = solve(prob, DFBDF())

    # Plots.plot(sol; vars=[inertia1.w, -inertia2.w*2])

    @test_skip begin
        @test sol.retcode == :Success
        @test all(isapprox.(sol[inertia1.w], -sol[inertia2.w]*2, atol=1)) # exact opposite oscillation with smaller amplitude J2 = 2*J1
        @test_broken all(sol[torque.flange.tau] .== -sol[sine.output.u]) # torque source is equal to negative sine
    end
end

# see: https://doc.modelica.org/Modelica%204.0.0/Resources/helpWSM/Modelica/Modelica.Mechanics.Rotational.Examples.First.html
@testset "first example" begin
    amplitude = 10 # Amplitude of driving torque
    frequency = 5 # Frequency of driving torque
    J_motor = 0.1 # Motor inertia
    J_load = 2 # Load inertia
    ratio = 10 # Gear ratio
    damping = 10 # Damping in bearing of gear

    @named fixed = Fixed()
    @named torque = Torque(use_support=true)
    @named inertia1 = Inertia(J=J_motor)
    @named idealGear = IdealGear(ratio=ratio, use_support=true)
    @named inertia2 = Inertia(J=2)
    @named spring = Spring(c=1e4)
    @named inertia3 = Inertia(J=J_load)
    @named damper = Damper(d=damping)
    @named sine = Blocks.Sine(amplitude=amplitude, frequency=frequency)

    connections = [
        connect(inertia1.flange_b, idealGear.flange_a)
        connect(idealGear.flange_b, inertia2.flange_a)
        connect(inertia2.flange_b, spring.flange_a)
        connect(spring.flange_b, inertia3.flange_a)
        connect(damper.flange_a, inertia2.flange_b)
        connect(damper.flange_b, fixed.flange)
        connect(sine.output, torque.tau)
        connect(torque.support, fixed.flange)
        connect(idealGear.support, fixed.flange)
        connect(torque.flange, inertia1.flange_a)
    ]

    @named model = ODESystem(connections, t, systems=[fixed, torque, inertia1, idealGear, inertia2, spring, inertia3, damper, sine])
    sys = structural_simplify(model)
    @test_broken prob = ODAEProblem(sys, Pair[], (0, 1.0))
    # sol = solve(prob, Rodas4())
    # @test sol.retcode == :Success
    # Plots.plot(sol; vars=[inertia2.w, inertia3.w])
end