text:
@time (using DamagedShearBand ; const DSB = DamagedShearBand)
#ENV["JULIA_DEBUG"] = DamagedShearBand

# σ₃ = -1e6
# θ = 60.0
# r = DSB.Rheology(D₀=0.2)
# Sc = get_damage_onset(r,σ₃,D)
# σᵢⱼ = DSB.build_principal_stress_tensor(r,5,σ₃,Dc)

# σᵢⱼ_r = DSB.band_coords(σᵢⱼ,θ)

# Dc = DSB.get_KI_mininizer_D(r)
# σᵢⱼ_Dc = DSB.build_principal_stress_tensor(r,100,σ₃,Dc)
# KI_Dc = DSB.compute_KI(r, σᵢⱼ_Dc, Dc)

# σᵢⱼ_Dc_vec = [DSB.build_principal_stress_tensor(r,S,σ₃,Dc) for S in S_vec]
# τ_vec = [DSB.get_τ(dev(σᵢⱼ)) for σᵢⱼ in σᵢⱼ_Dc_vec]
# p_vec = [-1/3*tr(σᵢⱼ) for σᵢⱼ in σᵢⱼ_Dc_vec]
# τ_vec./p_vec
# A,B = DSB.compute_AB(r,Dc)
# A/B
# KI_vec = [DSB.compute_KI(r, σᵢⱼ, Dc) for σᵢⱼ in σᵢⱼ_Dc_vec]

# KI_vec = [DSB.compute_KI(r, σ,τ, Dc) for (σ,τ) in zip(-p_vec,τ_vec)]
# -A.*p_vec .+ B.*τ_vec

# ##
# r = DSB.Rheology(D₀=0.24)
# D_out_i = Dc = DSB.get_KI_mininizer_D(r)
# D_in_i = 0.9
# S_vec = 1:0.1:10
# D_vec = r.D₀:0.01:0.99
# KI_mat = [DSB.compute_KI(r, DSB.build_principal_stress_tensor(r,S,σ₃,D), D) for S in S_vec, D in D_vec]
# DSB.compute_KI(r, DSB.build_principal_stress_tensor(r,10,σ₃,Dc), Dc)
# Sc = get_damage_onset(r,σ₃,D_in_i)

# heatmap(D_vec,S_vec,KI_mat./r.K₁c,clim=(0,0.1))
# xlabel!("D")
# ylabel!("S")

## Test solve
# σ₃ = -1e6
# θ = 60
# r = DSB.Rheology(D₀=0.0)
# # set damage inside the band
# D_in_i = r.D₀
# # find S so that damage starts frowing inside the band
# S = get_damage_onset(r,σ₃,D_in_i)
# # set damage outside so that it is the most unlikely to grow with increase of S
# D_out_i = D_in_i#DSB.get_KI_mininizer_D(r,S,σ₃)
# # build stress tensors and rotate them in the band coordinates system
# σᵢⱼ_in = DSB.build_principal_stress_tensor(r,S,σ₃,D_in_i)
# σᵢⱼ_out = DSB.build_principal_stress_tensor(r,S,σ₃,D_out_i)
# σᵢⱼ_in_ξη = DSB.band_coords(σᵢⱼ_in,θ)
# σᵢⱼ_out_ξη = DSB.band_coords(σᵢⱼ_out,θ)
# # get strains
# ϵᵢⱼ_in_ξη = DSB.compute_ϵij(r,D_in_i,σᵢⱼ_in_ξη)
# ϵᵢⱼ_out_ξη = DSB.compute_ϵij(r,D_out_i,σᵢⱼ_out_ξη)

# # They are not equal on the ξξ componant => two strategies :
# # 1 - initialize with the same initial damage (anything) in and out and only allow the inner damage to grow
# # 2 - initialize with the same initial damage (the lower KI D) and increase D_in just a bit using a modified solve function
# #     so that the deviatoric stress increase will only affect the inside of the band (puts some constraint on ϵ̇)
# Δt = 1
# Ṡ = 1e-6
# σᵢⱼ_in_next = DSB.build_principal_stress_tensor(r,S+Ṡ*Δt,σ₃,D_in_i)
# Δσ = DSB.band_coords(σᵢⱼ_in_next,θ) - σᵢⱼ_in_ξη

# ϵᵢⱼ_in_ξη_next = DSB.compute_ϵij(r,D_in_i,DSB.band_coords(σᵢⱼ_in_next,θ))
# Δϵ = ϵᵢⱼ_in_ξη_next - ϵᵢⱼ_in_ξη

# #C = DSB.compute_damaged_stiffness_tensor(r,(ϵᵢⱼ_in_ξη_next+ϵᵢⱼ_in_ξη)/2,D_in_i)
# #Δσ2 = C ⊡ Δϵ

# #knowns = Δσηη , Δσξη, Δϵξξ
# #unknowns = Δσξξ, Δϵηη, Δϵξη
# knowns = Vec(Δσ[2,2] , Δσ[1,2], Δϵ[1,1])
# unknowns = Vec(Δσ[1,1], Δϵ[2,2], Δϵ[1,2])
# #res = DSB.residual(r, ϵᵢⱼ_in_ξη, D_out_i, knowns, unknowns)

# σᵢⱼ_in_next, ϵᵢⱼ_in_next, knowns, u, D = DSB.solve(r,σᵢⱼ_in_ξη,ϵᵢⱼ_in_ξη,knowns,unknowns,D_out_i,Δt ; abstol=1e-12, maxiter=100)

# KI_vec = [DSB.KI_to_minimize(r,S,σ₃,D) for D in r.D₀:0.01:0.99]
# plot(r.D₀:0.01:0.99,KI_vec)
# Dc = DSB.get_KI_mininizer_D(r,S,σ₃)
σ₂₂
## Second strategy
using DamagedShearBand ; const DSB = DamagedShearBand
using Plots
plotlyjs()

tspan = (0.0,24)
Δt = 10
ϵ̇11 = -1e-4
r = DSB.Rheology(D₀=0.25) # object containing elastic moduli and damage parameters.
σ₃ = -1e6
S = 1
D_i = r.D₀
σᵢⱼ_i = DSB.build_principal_stress_tensor(r,S,σ₃,D_i ; abstol=1e-15) # takes care of the plane strain constraint
ϵᵢⱼ_i = DSB.compute_ϵij(r,D_i,σᵢⱼ_i)
#σᵢⱼ_i2,ϵᵢⱼ_i2,D_i2 = σᵢⱼ_vec_a[end], ϵᵢⱼ_vec_a[end], D_vec_a[end]
#σᵢⱼ_i,ϵᵢⱼ_i,D_i = σᵢⱼ_i2,ϵᵢⱼ_i2,D_i2

#t_vec_a, σᵢⱼ_vec_a, ϵᵢⱼ_vec_a, D_vec_a = DSB.time_integration(r,σᵢⱼ_i,ϵᵢⱼ_i,D_i,ϵ̇11,Δt,tspan ; abstol=1e-12, maxiter=100)
t_vec_a, σᵢⱼ_vec_a, ϵᵢⱼ_vec_a, D_vec_a = DSB.adaptative_time_integration(r,σᵢⱼ_i,ϵᵢⱼ_i,D_i,ϵ̇11,Δt,tspan ; abstol=1e-12, time_maxiter=nothing, newton_maxiter=100, e₀=(D=1e-4, σ=1.0, ϵ=1e-6))

# σ = 1/3 *tr(σᵢⱼ_i)
# I = SymmetricTensor{2,3}(DSB.δ)
# sij = σᵢⱼ_i - σ*I
# -DSB.get_τ(sij)/σ

# DSB.compute_KI(r,σᵢⱼ_i,D_i)
#σᵢⱼnext, ϵᵢⱼnext, D, u = DSB.solve_2(r, σᵢⱼ_i, ϵᵢⱼ_i, D_i, ϵ̇11, Δt ; abstol=1e-12, maxiter=10)

# ϵᵢⱼnext_test = DSB.insert_into(ϵᵢⱼ_i, (ϵᵢⱼ_i[1,1] + ϵ̇11*Δt), (1,1))
# σᵢⱼnext_test = DSB.compute_σij(r,D_i,ϵᵢⱼnext_test)
##
Ce = DSB.get_elastic_stiffness_tensor(r::Rheology)
σᵢⱼnext_e = Ce ⊡ ϵᵢⱼ_vec[end]
σᵢⱼ_vec_elast = Ref(Ce) .⊡ ϵᵢⱼ_vec_a
pente_model = (σᵢⱼ_vec_a[end][1,1]-σᵢⱼ_vec_a[2][1,1]) / (ϵᵢⱼ_vec_a[end][1,1]-ϵᵢⱼ_vec_a[2][1,1])
pente_elast = (σᵢⱼ_vec_elast[end][1,1]-σᵢⱼ_vec_elast[2][1,1]) / (ϵᵢⱼ_vec_a[end][1,1]-ϵᵢⱼ_vec_a[2][1,1])
E = DSB.E_from_Gν(r.G,r.ν)
tan_modulus = E/(1-r.ν^2)
# ϵ̇ij_ana, Ḋ = DSB.compute_ϵ̇ij(r,D_i,σᵢⱼ_i,σᵢⱼnext,Δt)

##
plot(t_vec,D_vec)
  ylims!((D_vec[1],D_vec[end]))
plot(t_vec,[σᵢⱼ[1,1] for σᵢⱼ in σᵢⱼ_vec])

  #plot!(t_vec,[σᵢⱼ[1,1] for σᵢⱼ in σᵢⱼ_vec_elast])
#plot(t_vec,[σᵢⱼ[3,3] for σᵢⱼ in σᵢⱼ_vec])
plot(t_vec_a,D_vec_a)
plot(diff(t_vec_a))
plot(t_vec_a,D_vec_a)
D_vec_a[end]-D_vec_a[1]
plot(t_vec_a,[σᵢⱼ[1,1] for σᵢⱼ in σᵢⱼ_vec_a])
plot(t_vec_a,[ϵᵢⱼ[1,1] for ϵᵢⱼ in ϵᵢⱼ_vec_a])
plot([-ϵᵢⱼ[1,1] for ϵᵢⱼ in ϵᵢⱼ_vec_a],[-σᵢⱼ[1,1] for σᵢⱼ in σᵢⱼ_vec_a])

plot!(t_vec_a[2:end],-diff(t_vec_a).*3e8)

function merge_series!(sp1::Plots.Subplot, sp2::Plots.Subplot)
  append!(sp1.series_list, sp2.series_list)
  Plots.expand_extrema!(sp1[:xaxis], xlims(sp2))
  Plots.expand_extrema!(sp1[:yaxis], ylims(sp2))
  Plots.expand_extrema!(sp1[:zaxis], zlims(sp2))
  return sp1
end

function merge_series!(plt, plts...)
  for (i, sp) in enumerate(plt.subplots)
      for other_plt in plts
          if i in eachindex(other_plt.subplots)
              merge_series!(sp, other_plt[i])
          end
      end
  end
  return plt
end
##
layout = @layout [a ; b]
l2 = @layout [a ; b]

p1 = plot(-ϵ₁₁_vec,-σ₁₁_vec)
#plot(p1,-ϵ₁₁_vec,-σ₁₁_vec)
title!(p1,"D0 = $(r.D₀) & epsdot = $(ϵ̇₁₁) & p=$(-σ₂₂)")
ylabel!("-σ₁₁")

p2 = plot(-ϵ₁₁_vec,D_vec,c=:red)
#plot(p2,-ϵ₁₁_vec,D_vec,c=:red)
ylabel!("D")
xlabel!("-ϵ₁₁")

pl1 = plot(p1,p2,layout=layout)
##
p3 = plot(-ϵ₁₁_vec.+0.0001,-σ₁₁_vec)
#plot(p1,-ϵ₁₁_vec,-σ₁₁_vec)
title!(p1,"D0 = $(r.D₀) & epsdot = $(ϵ̇₁₁) & p=$(-σ₂₂)")
ylabel!("-σ₁₁")

p4 = plot(-ϵ₁₁_vec,D_vec,c=:red)
#plot(p2,-ϵ₁₁_vec,D_vec,c=:red)
ylabel!("D")
xlabel!("-ϵ₁₁")

pl2 = plot(p3,p4,layout=l2)

##
function test_func(args...)
  return sum([args...])
end

test_func(1,2,3)
nt = (a=1,b=2)
@time isdefined(nt,:a)

D_vec = 0.6495327404332681:0.000001:0.6495527404332681
KI_vec = [DSB.compute_KI(r,DSB.build_principal_stress_tensor(r,22.207207207207208,σ₃,D ; abstol=1e-15),D) for D in D_vec]
_,ind = findmin(KI_vec)
D_vec[ind]

op = DSB.OutputParams(save_frequency = n_iter_saved/(tspan[2]-tspan[1]),
                      save_period = 1000,
                      print_frequency = n_iter_printed/(tspan[2]-tspan[1]),
                      print_period = 1000)

sp = DSB.SolverParams(newton_abstol = 1e-12,
                      newton_maxiter = 100,
                      time_maxiter = nothing,
                      e₀ = 1e-4, #(D=1e-2, S=1e-2, σ=100.0, ϵ=1e-2))
                      adaptative_maxrecursions = 50)
p = DSB.Params(op,sp)


text:
module GLVisualizeCI

import GitHub, HttpCommon

dir(paths...) = normpath(joinpath(dirname(@__FILE__), "..", paths...))

function push_status(pr, commit)
    try
        repo = LibGit2.GitRepo(GLVisualizeCI.dir())
        LibGit2.fetch(repo)
        LibGit2.add!(repo, ".")
        LibGit2.commit(repo, "update $pr $commit")
        LibGit2.push(repo, refspecs = ["refs/heads/master"])
    catch e
        warn("couldn't update report for $pr: $e")
    end
end

function report_url(ci, repo, pr)
    "https://github.com/SimiDCI/GLVisualizeCI.jl/blob/master/reports/$ci/$repo/$pr"
end

gitclone!(repo, path) = run(`git clone https://github.com/$(repo).git $(path)`)

function test_pr(package, repo, pr)
    mktempdir() do path
        cd(homedir()) # make sure, we're in a concrete folder
        # init a new julia package repository
        ENV["JULIA_PKGDIR"] = path
        Pkg.init()

        builddir = Pkg.dir(package)
        gitclone!(repo, builddir)
        cd(builddir)
        # Fetch current PR
        try
           run(`git fetch --quiet origin +refs/pull/$(pr)/merge:`)
        catch
           # if there's not a merge commit on the remote (likely due to
           # merge conflicts) then fetch the head commit instead.
           run(`git fetch --quiet origin +refs/pull/$(pr)/head:`)
        end
        run(`git checkout --quiet --force FETCH_HEAD`)
        Pkg.add("GLVisualize") # this checks out the dependencies after a fetch
        Pkg.build("GLVisualize")
        julia_exe = Base.julia_cmd()
        logstd = ENV["CI_REPORT_DIR"] * "/log_sttd.txt"
        logsterr = ENV["CI_REPORT_DIR"] * "/log_sterr.txt"
        testcmd = `$julia_exe -e $("Pkg.test(\"$package\", coverage = true)")`
        coveragecmd = `$julia_exe $(dir("src", "submit_coverage.jl"))`
        try
            run(pipeline(testcmd, stdout = `tee $logstd`, stderr = `tee $logsterr`))
            run(pipeline(coveragecmd, stdout = `tee $logstd`, stderr = `tee $logsterr`))
        catch
            println("Test not succesful!")
            return false
        end
        return true
    end
end

function handle_event(name, event, auth)
    kind, payload, repo = event.kind, event.payload, event.repository

    if kind == "pull_request" && payload["action"] != "closed"
        sha = event.payload["pull_request"]["head"]["sha"]
        pr = string(event.payload["pull_request"]["number"])
        package, jl = splitext(get(repo.name))

        target_url = report_url(name, package, pr)
        @show target_url

        path = dir("reports", name)
        isdir(path) || mkdir(path)
        path = dir("reports", name, package)
        isdir(path) || mkdir(path)
        path = dir("reports", name, package, pr)
        isdir(path) || mkdir(path)

        push_status(pr, sha)

        GitHub.create_status(repo, sha; auth = auth, params = Dict(
            "state" => "pending",
            "context" => name,
            "description" => "Running CI...",
            "target_url" => target_url
        ))
        success = false
        try
            ENV["CI_REPORT_DIR"] = path
            ENV["CI"] = "true"
            success = test_pr(package, get(repo.full_name), pr)
            push_status(pr, sha)
        catch err
            GitHub.create_status(repo, sha; auth = auth, params = Dict(
                "state" => "error",
                "context" => name,
                "description" => "Error!",
                "target_url" => target_url
            ))
            return HttpCommon.Response(500)
        end
        GitHub.create_status(repo, sha; auth = auth, params = Dict(
            "state" => success ? "success" : "failure",
            "context" => name,
            "description" => "Build has finished.",
            "target_url" => target_url
        ))
        return HttpCommon.Response(202, "success")
    else
        return HttpCommon.Response(500)
    end
end

function start(name, func = handle_event;
        host = IPv4(128, 30, 87, 54),
        port = 8000,
        myrepos = [GitHub.Repo("JuliaGL/GLVisualize.jl")],
        myauth = GitHub.authenticate(ENV["GITHUB_AUTH"]),
        mysecret = ENV["GITHUB_SECRET"],
        myevents = ["pull_request"]
    )
    listener = GitHub.EventListener(
        event-> handle_event(name, event, myauth),
        auth = myauth,
        secret = mysecret,
        repos = myrepos,
        events = myevents
    )
    GitHub.run(listener, host = host, port = port)
end

end # module

ext: jl

max_stars_repo_name: SimiDCI/GLVisualizeCI

max_issues_repo_name: SimiDCI/GLVisualizeCI

max_forks_repo_name: SimiDCI/GLVisualizeCI


text:
module Validation

using DifferentialEquations
using Random
using Suppressor

using Model
import Model.Model_meta
using Sim

include("Setting.jl")
using .Setting

include("Sampling.jl")

## SETTING
const S_SPAN = Sampling.S_SPAN
const DEL_CONC=Setting.DEL_CONC
const TRANS_THRESHOLD = Setting.TRANS_THRESHOLD
const tspan = Setting.val_tspan
const SSMETHOD = Setting.SSMETHOD
const DEFAULT_MODEL = Setting.MODEL
##

model_meta = Model_meta.get_model_meta(DEFAULT_MODEL)
const meta = Model_meta.get_model_meta(DEFAULT_MODEL)
const cond, protein_lookup = meta.cond, meta.protein_lookup

export get_output, try_simulate, try_conditions, get_del_u

function try_conditions!(u,p;
                        model=DEFAULT_MODEL,
                        BREAK=false,
                        Cond_Random=true)
        valid = 1
        steady = 1

        num_cond = size(cond)[1]

        if Cond_Random # Shuffling the conditions
            trial_cond_order = shuffle( Int.(range(1,stop=num_cond,length=num_cond ) ) )
        else
            trial_cond_order = Int.(range(1,stop=num_cond,length=num_cond) )
        end

        test_log = Vector{Union{Missing, Int16}}(missing, num_cond)

        # get initial condition
        u[:] = Sim.solve_SSODE(model, u, p) # This step changed the initial value, and reset as the steady ones

        for trial_n in trial_cond_order # Test conditions
            con = cond[trial_n,:]

            if typeof(con.Trans2Nuc) == Missing
                continue # some conditions are not measured
            end

            ud = get_del_u(u, con.rtg1, con.rtg2, con.rtg3, con.mks, con.s, protein_lookup; del_conc=DEL_CONC, s_span=S_SPAN)

            sol = @suppress begin
                    Sim.solve_SSODE(model, ud, p; method=SSMETHOD)
                end

            # Check the steady-state is real
            if sol.retcode==:Failure
                steady = 0
                valid = 0 #unstable system
                break
            end

            trans2nuc = get_output(sol, con.gfp, protein_lookup, threshold=TRANS_THRESHOLD) # 1 means nucleus has significantly higer concentration than cytosol

            if trans2nuc != con.Trans2Nuc
                valid = 0
                test_log[trial_n] = 0
                if BREAK # Early termination when invalid condition is confirmed
                    break
                end
            else
                test_log[trial_n] = 1
            end
        end

        if valid == 1 # Simulation test
            if try_simulate(model,p,u) == 0
                valid = 0
            end
        end

        return valid, test_log, steady
end

"
Compare the cytosolic concentration of either 'rtg1' or 'rtg3' with their nucleus concentrations.
Output:
1 : nucleus concentration is higher than the cytosolic one
0 : otherwise
"
function get_output(sol, gfp, protein_lookup; threshold = TRANS_THRESHOLD)

    if gfp == "rtg1"
        cyt_index = protein_lookup[:Rtg1_c]
        nuc_index = protein_lookup[:Rtg1_n]
    elseif gfp == "rtg3"
        cyt_index = protein_lookup[:Rtg3_c]
        nuc_index = protein_lookup[:Rtg3_n]
    else
        throw(MathodError)
    end

    total_conc_cyt = sum(sol[cyt_index])
    total_conc_nuc = sum(sol[nuc_index])

    if total_conc_nuc > threshold*total_conc_cyt
        output = 1
    else
        output = 0
    end
    return output
end

"
Testing on solving ODE
"
function try_simulate(model, p, u; tspan=tspan)
    st = 1
    try
        println("try simulate")
        oprob = ODEProblem(model, u, tspan, p )
        osol = solve(oprob)
        println("Simulation successed")
    catch y
        println(y) # What to do when an error occurs.
        st = 0
    end
    return st
end

"
Copy vector u, and delete elements with lookup table
"
function get_del_u(u, rtg1, rtg2, rtg3, mks, s, protein_lookup; del_conc=DEL_CONC, s_span=S_SPAN )

    del_info = Dict([:Rtg1, :Rtg2, :Rtg3, :Mks] .=> [rtg1, rtg2, rtg3, mks])
    ud = deepcopy(u)

    for protein in keys(del_info)
        protein_existence = del_info[protein]
        if protein_existence == 0
            ud[protein_lookup[protein]] .= del_conc
        end
    end

    ud[protein_lookup[:s]] = s_span[s+1] #in boolean, s=0 means 1 in span_tuple

    return ud
end

"Qualitative Validation"

"Substraction of nuclear concentration with cytoplasmic concentration"
function get_outputQ(sol, gfp, protein_lookup; threshold = TRANS_THRESHOLD)

    if gfp == "rtg1"
        cyt_index = protein_lookup[:Rtg1_c]
        nuc_index = protein_lookup[:Rtg1_n]
    elseif gfp == "rtg3"
        cyt_index = protein_lookup[:Rtg3_c]
        nuc_index = protein_lookup[:Rtg3_n]
    else
        throw(MathodError)
    end

    total_conc_cyt = sum(sol[cyt_index])
    total_conc_nuc = sum(sol[nuc_index])

    output = total_conc_nuc - threshold*total_conc_cyt

    return output
end

end


text:
using OrdinaryDiffEq
using PyPlot
using LinearAlgebra

pygui(true);

const M = 31 #number of clusters (odd)
N0 = 10000 #total number of atoms

κ = 1
g,Γ,κ,R = (0.002κ, 0.001κ, κ, 0.001κ)

σ = 0.1κ
ξ0 = ones(M) #atomic dephasing
ξ = ξ0*0Γ #zero dephasing

##################################
#Gaussian distribution
dD = 2/(M-1)
Δ0 = [-1:dD:1;] #odd number of clusters with different detunings

#shift = [randn(Integer((M-1)/2)); 0; randn(Integer((M-1)/2))].*dD./10
#Δ0 = Δ0 .+shift

Δ = Δ0.*σ

σ3 = σ
G = 1/(σ3*sqrt(2pi)).*exp.(-1/2 .*Δ.^2 ./σ3^2)
N = G #figure (a)

#N[12]=N[12]+100 #figure (b)

#shift = randn(M).*N0/10M
#N = N .+shift #figure (c)

N=N./sum(N)*N0

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

#ind = []
const ind = Tuple{Int,Int}[] # use specific type and avoid non-constant global
for a=1:M
    for b=a:M
        push!(ind,(a,b))
    end
end

function mj(m, j)
    findfirst(isequal((m, j)), ind)  #findall(isequal((m, j)), x)
end

const ind_spsm = length(ind)
const ind_ata = ind_spsm + M
const ind_asp = ind_ata + 1
const ind_end = ind_asp + M
const N_ = copy(N)

function ff(du,u,p,t)
    Δ,g,Γ,κ,R,ξ = p

        #m - the number of cluster
        #dspsm/dt, m not equal j
        for m=1:M
            for j=m:M
                a = mj(m, j)
                du[a] = -(Γ+R+ξ[m]/2+ξ[j]/2)*u[a] - im*g*u[ind_asp+m] + im*g*u[ind_asp+j]' + 2*im*g*u[ind_asp+m]*u[ind_spsm+j] - 2*im*g*u[ind_asp+j]'*u[ind_spsm+m] -im*Δ[m]*u[a] +im*Δ[j]*u[a]
            end
        end

        #dspsm/dt, m=j
        for m=1:M
            du[ind_spsm+m] = im*g*u[ind_asp+m]' - im*g*u[ind_asp+m] - (Γ+R)*u[ind_spsm+m] + R
        end

        #dat*a/dt
        du[ind_ata+1] = -κ*u[ind_ata+1] + im*g*sum(N_[j]*u[ind_asp+j] for j=1:M) - im*g*sum(N_[j]*u[ind_asp+j]' for j=1:M)

        #da*sp/dt
        for m=1:M
            du[ind_asp+m] = -(κ/2+Γ/2+R/2+ξ[m]/2+im*Δ[m])*u[ind_asp+m] + im*g*u[ind_ata+1] - im*g*u[ind_spsm+m] - 2*im*g*u[ind_ata+1]*u[ind_spsm+m] - im*g*(N_[m]-1)*u[mj(m,m)]
            for j=1:M
                if j>m   du[ind_asp+m] = du[ind_asp+m] - im*g*N_[j]*u[mj(m, j)] end
                if j<m   du[ind_asp+m] = du[ind_asp+m] - im*g*N_[j]*u[mj(j, m)]' end
            end
        end
end

u0 = zeros(ComplexF64, ind_end)
p0 = (Δ,g,Γ,κ,R,ξ)

prob = ODEProblem(ff,u0,(0.0,10000),p0)
sol = solve(prob,RK4())

n = real.(getindex.(sol.u, ind_ata +1))
pe = real.(getindex.(sol.u, ind_spsm+1))
sz = 2*pe .-1

#=

figure(figsize=(9,3))
subplot(121)
title(L"N=5000, \sigma=0.01\kappa, \Delta \in [-\sigma; \sigma], g=0.005\kappa, \Gamma=0.001\kappa, R=0.01\kappa, \kappa=1")
plot(sol.t, n, label = (L"<n>_{clusters}"))
legend()
xlabel(L"\kappa t")
ylabel(L"\langle n \rangle")

subplot(122)
plot(sol.t, sz, label = (L"<S_z>_{clusters}"))
legend()
xlabel(L"\kappa t")
ylabel(L"\langle S_z \rangle")

tight_layout()
=#

###################################################
# Spectrum Laplace transform
pe0=zeros(M)
for m=1:M
    pe0[m] = real.(getindex.(sol.u, ind_spsm+m))[end]
end

u0_spec = zeros(ComplexF64, M+1)
u0_spec[1] = real.(getindex.(sol.u, ind_ata+1))[end]
for m=1:M
    u0_spec[m+1] = getindex.(sol.u, ind_asp+m)[end]
end

a = (0, [-im*g*N[m] for m=1:M]...)
b = (0, [-im*g*(1-2*pe0[m]) for m=1:M]...)
d = [κ/2, [(Γ+R+ξ[m])/2+im*Δ[m] for m=1:M]...]

#A = Matrix{ComplexF64}(I, M+1, M+1).*d
A = zeros(ComplexF64, M+1, M+1)
A[1,:].= a
A[:,1].=b

A = A + Diagonal(d)

w = [-0.11:0.0001:0.11;]

spec2 = zeros(ComplexF64, length(w))

for i=1:length(w)
s = im*w[i]
S = Diagonal(s*ones(M+1))
B = S+A
spec2[i] = (inv(B)*u0_spec)[1]
end

spec2 = real.(spec2)
spec2 = spec2./maximum(spec2)

rc("font", size = 9)
rc("legend", fontsize = 7)
rc("text", usetex = true)
rc("lines", lw = 1.0)

figure(figsize=(5,3))
title(L"N=10^4, M=31, \sigma = 0.1\kappa, \Delta \in [-\sigma:\sigma], g=0.002\kappa, \Gamma=0.001\kappa")
plot(w, spec2, linewidth = 1,label = (L"R=0.001\kappa"), color = "silver")
legend()
xlabel(L"(\omega - \omega_c) / \kappa")
ylabel(L"Intensity")
tight_layout()

max_stars_repo_name: by-anna/Clusters

max_issues_repo_name: by-anna/Clusters

max_forks_repo_name: by-anna/Clusters


text:
#ДАНО: Где-то на неограниченном со всех сторон поле и без внутренних перегородок имеется единственный маркер.
#Робот - в произвольной клетке поля.
#РЕЗУЛЬТАТ: Робот - в клетке с тем маркером.
function poisk_markera(r::Robot)
    check = false
    side = 0
    counter = 1
    c = 0

    while !check
        for _ in 1:counter
            if ismarker(r)
                check = true
                break
            end
            move!(r, HorizonSide(side))
        end
        c+=1
        side = (side+1) % 4
        if c == 2
            counter+=1
            c = 0
        end
    end
end

max_stars_repo_name: listova-nastya/mirea-progs

max_issues_repo_name: listova-nastya/mirea-progs

max_forks_repo_name: listova-nastya/mirea-progs


text:
import Base: error, size

export MagickWand,
    constituteimage,
    exportimagepixels!,
    getblob,
    getimagealphachannel,
    getimagecolorspace,
    getimagedepth,
    getnumberimages,
    #importimagepixels,
    readimage,
    resetiterator,
    setimagecolorspace,
    setimagecompression,
    setimagecompressionquality,
    setimagedelay,
    setimageformat,
    writeimage

const depsfile = joinpath(dirname(@__DIR__), "deps", "deps.jl")
if isfile(depsfile)
    include(depsfile)
else
    error("ImageMagick not properly installed. Please run Pkg.build(\"ImageMagick\") then restart Julia.")
end

magickgenesis() = ccall((:MagickWandGenesis, libwand), Cvoid, ())
magickterminus() = ccall((:MagickWandTerminus, libwand), Cvoid, ())
isinstantiated() = ccall((:IsMagickWandInstantiated, libwand), Cint, ()) == 1

function __init__()
    check_deps()
    magickgenesis()
end

# Constants
# Storage types
const CHARPIXEL    = 1
const DOUBLEPIXEL  = 2
const FLOATPIXEL   = 3
const INTEGERPIXEL = 4
const SHORTPIXEL   = 7
const IMStorageTypes = Union{UInt8,UInt16,UInt32,Float32,Float64}
storagetype(::Type{UInt8})   = CHARPIXEL
storagetype(::Type{UInt16})  = SHORTPIXEL
storagetype(::Type{UInt32})  = INTEGERPIXEL
storagetype(::Type{Float32}) = FLOATPIXEL
storagetype(::Type{Float64}) = DOUBLEPIXEL
storagetype(::Type{T}) where {T<:Normed} = storagetype(FixedPointNumbers.rawtype(T))
storagetype(::Type{CV}) where {CV<:Colorant} = storagetype(eltype(CV))

# Channel types
mutable struct ChannelType
    value::UInt32
end
const UndefinedChannel  = ChannelType(0x00000000)
const RedChannel        = ChannelType(0x00000001)
const GrayChannel       = ChannelType(0x00000001)
const CyanChannel       = ChannelType(0x00000001)
const GreenChannel      = ChannelType(0x00000002)
const MagentaChannel    = ChannelType(0x00000002)
const BlueChannel       = ChannelType(0x00000004)
const YellowChannel     = ChannelType(0x00000004)
const AlphaChannel      = ChannelType(0x00000008)
const MatteChannel      = ChannelType(0x00000008)
const OpacityChannel    = ChannelType(0x00000008)
const BlackChannel      = ChannelType(0x00000020)
const IndexChannel      = ChannelType(0x00000020)
const CompositeChannels = ChannelType(0x0000002F)
const TrueAlphaChannel  = ChannelType(0x00000040)
const RGBChannels       = ChannelType(0x00000080)
const GrayChannels      = ChannelType(0x00000080)
const SyncChannels      = ChannelType(0x00000100)
const AllChannels       = ChannelType(0x7fffffff)
const DefaultChannels   = ChannelType((AllChannels.value | SyncChannels.value) & ~OpacityChannel.value)

# Image type
const IMType = ["BilevelType", "GrayscaleType", "GrayscaleMatteType", "PaletteType", "PaletteMatteType", "TrueColorType", "TrueColorMatteType", "ColorSeparationType", "ColorSeparationMatteType", "OptimizeType", "PaletteBilevelMatteType"]
const IMTypedict = Dict([(IMType[i], i) for i = 1:length(IMType)])

const CStoIMTypedict = Dict("Gray" => "GrayscaleType", "GrayA" => "GrayscaleMatteType", "AGray" => "GrayscaleMatteType", "RGB" => "TrueColorType", "ARGB" => "TrueColorMatteType", "RGBA" => "TrueColorMatteType", "CMYK" => "ColorSeparationType", "I"=>"GrayscaleType", "IA"=>"GrayscaleMatteType", "AI"=>"GrayscaleMatteType", "BGRA"=>"TrueColorMatteType", "ABGR"=>"TrueColorMatteType")

# Colorspace
const IMColorspace = ["RGB", "Gray", "Transparent", "OHTA", "Lab", "XYZ", "YCbCr", "YCC", "YIQ", "YPbPr", "YUV", "CMYK", "sRGB"]
const IMColordict = Dict([(IMColorspace[i], i) for i = 1:length(IMColorspace)])
for AC in vcat(subtypes(AlphaColor), subtypes(ColorAlpha))
    Cstr = ColorTypes.colorant_string(color_type(AC))
    if haskey(IMColordict, Cstr)
        IMColordict[ColorTypes.colorant_string(AC)] = IMColordict[Cstr]
    end
end

flip1(A) = view(A, reverse(axes(A,1)), ntuple(x->Colon(),ndims(A)-1)...)
flip2(A) = view(A, :, reverse(axes(A,2)), ntuple(x->Colon(),ndims(A)-2)...)
flip12(A) = view(A, reverse(axes(A,1)), reverse(axes(A,2)), ntuple(x->Colon(),ndims(A)-2)...)

vertical_major(img::AbstractVector) = img
vertical_major(A) = PermutedDimsArray(A, [2;1;3:ndims(A)])

const orientation_dict = Dict(
    nothing => (A,ph) -> ph ? vertical_major(A) : identity(A),
    "1" => (A,ph) -> ph ? vertical_major(A) : identity(A),
    "2" => (A,ph) -> ph ? vertical_major(flip1(A)) : flip1(A),
    "3" => (A,ph) -> ph ? vertical_major(flip12(A)) : flip12(A),
    "4" => (A,ph) -> ph ? vertical_major(flip2(A)) : flip2(A),
    "5" => (A,ph) -> ph ? identity(A) : vertical_major(A),
    "6" => (A,ph) -> ph ? flip2(A) : vertical_major(flip2(A)),
    "7" => (A,ph) -> ph ? flip12(A) : vertical_major(flip12(A)),
    "8" => (A,ph) -> ph ? flip1(A) : vertical_major(flip1(A)))

function nchannels(imtype::AbstractString, cs::AbstractString, havealpha = false)
    n = 3
    if startswith(imtype, "Grayscale") || startswith(imtype, "Bilevel")
        n = 1
        cs = havealpha ? "GrayA" : "Gray"
    elseif cs == "CMYK"
        n = 4
    else
        cs = havealpha ? "ARGB" : "RGB" # only remaining variants supported by exportimagepixels
    end
    n + havealpha, cs
end

# channelorder = ["Gray" => "I", "GrayA" => "IA", "RGB" => "RGB", "ARGB" => "ARGB", "RGBA" => "RGBA", "CMYK" => "CMYK"]

# Compression
const NoCompression = 1

mutable struct MagickWand
    ptr::Ptr{Cvoid}

    function MagickWand()
        ptr = ccall((:NewMagickWand, libwand), Ptr{Cvoid}, ())
        ptr == C_NULL && throw(OutOfMemoryError())
        obj = new(ptr)
        finalizer(free, obj)
        obj
    end
end

function Base.unsafe_convert(::Type{Ptr{Cvoid}}, wand::MagickWand)
    ptr = wand.ptr
    ptr == C_NULL && throw(UndefRefError())
    ptr
end

function free(wand::MagickWand)
    ptr = wand.ptr
    if ptr != C_NULL
        ccall((:DestroyMagickWand, libwand), Ptr{Cvoid}, (Ptr{Cvoid},), ptr)
    end
    wand.ptr = C_NULL
    nothing
end

mutable struct PixelWand
    ptr::Ptr{Cvoid}

    function PixelWand()
        ptr = ccall((:NewPixelWand, libwand), Ptr{Cvoid}, ())
        ptr == C_NULL && throw(OutOfMemoryError())
        obj = new(ptr)
        finalizer(free, obj)
        obj
    end
end

function Base.unsafe_convert(::Type{Ptr{Cvoid}}, wand::PixelWand)
    ptr = wand.ptr
    ptr == C_NULL && throw(UndefRefError())
    ptr
end

function free(wand::PixelWand)
    ptr = wand.ptr
    if ptr != C_NULL
        ccall((:DestroyPixelWand, libwand), Ptr{Cvoid}, (Ptr{Cvoid},), ptr)
    end
    wand.ptr = C_NULL
    nothing
end

const IMExceptionType = Ref{Cint}()
function error(wand::MagickWand)
    pMsg = ccall((:MagickGetException, libwand), Ptr{UInt8}, (Ptr{Cvoid}, Ptr{Cint}), wand, IMExceptionType)
    msg = unsafe_string(pMsg)
    relinquishmemory(pMsg)
    error(msg)
end
function error(wand::PixelWand)
    pMsg = ccall((:PixelGetException, libwand), Ptr{UInt8}, (Ptr{Cvoid}, Ptr{Cint}), wand, IMExceptionType)
    msg = unsafe_string(pMsg)
    relinquishmemory(pMsg)
    error(msg)
end

function getsize(buffer, channelorder)
    if channelorder == "I"
        return size(buffer, 1), size(buffer, 2), size(buffer, 3)
    else
        return size(buffer, 2), size(buffer, 3), size(buffer, 4)
    end
end
getsize(buffer::AbstractArray{C}, channelorder) where {C<:Colorant} = size(buffer, 1), size(buffer, 2), size(buffer, 3)

colorsize(buffer, channelorder) = channelorder == "I" ? 1 : size(buffer, 1)
colorsize(buffer::AbstractArray{C}, channelorder) where {C<:Colorant} = 1

bitdepth(buffer::AbstractArray{C}) where {C<:Colorant} = 8*sizeof(eltype(C))
bitdepth(buffer::AbstractArray{T}) where {T} = 8*sizeof(T)

# colorspace is included for consistency with constituteimage, but it is not used
function exportimagepixels!(buffer::AbstractArray{T}, wand::MagickWand,  colorspace::String, channelorder::String; x = 0, y = 0) where T<:Unsigned
    cols, rows, nimages = getsize(buffer, channelorder)
    ncolors = colorsize(buffer, channelorder)
    if isa(buffer, Array)
        p = pointer(buffer)
    else
        tmp = similar(buffer)
        p = pointer(tmp)
    end
    for i = 1:nimages
        nextimage(wand)
        status = ccall((:MagickExportImagePixels, libwand), Cint, (Ptr{Cvoid}, Cssize_t, Cssize_t, Csize_t, Csize_t, Ptr{UInt8}, Cint, Ptr{Cvoid}), wand, x, y, cols, rows, channelorder, storagetype(T), p)
        status == 0 && error(wand)
        p += sizeof(T)*cols*rows*ncolors
    end
    isa(buffer, Array) || buffer .= tmp
    buffer
end

# function importimagepixels{T}(buffer::AbstractArray{T}, wand::MagickWand, colorspace::String; x = 0, y = 0)
#     cols, rows = getsize(buffer, colorspace)
#     status = ccall((:MagickImportImagePixels, libwand), Cint, (Ptr{Cvoid}, Cssize_t, Cssize_t, Csize_t, Csize_t, Ptr{UInt8}, Cint, Ptr{Cvoid}), wand, x, y, cols, rows, channelorder[colorspace], storagetype(T), buffer)
#     status == 0 && error(wand)
#     nothing
# end

function constituteimage(buffer::AbstractArray{T}, wand::MagickWand, colorspace::String, channelorder::String; x = 0, y = 0) where T<:Unsigned
    cols, rows, nimages = getsize(buffer, channelorder)
    ncolors = colorsize(buffer, channelorder)
    p = pointer(buffer)
    depth = bitdepth(buffer)
    for i = 1:nimages
        status = ccall((:MagickConstituteImage, libwand), Cint, (Ptr{Cvoid}, Cssize_t, Cssize_t, Ptr{UInt8}, Cint, Ptr{Cvoid}), wand, cols, rows, channelorder, storagetype(T), p)
        status == 0 && error(wand)
        setimagecolorspace(wand, colorspace)
        setimagetype(wand, buffer, channelorder)
        status = ccall((:MagickSetImageDepth, libwand), Cint, (Ptr{Cvoid}, Csize_t), wand, depth)
        status == 0 && error(wand)
        p += sizeof(T)*cols*rows*ncolors
    end
    nothing
end

function getblob(wand::MagickWand, format::AbstractString)
    setimageformat(wand, format)
    len = Ref{Csize_t}(1)
    ptr = ccall((:MagickGetImagesBlob, libwand), Ptr{UInt8}, (Ptr{Cvoid}, Ptr{Csize_t}), wand, len)
    blob = unsafe_wrap(Array, ptr, convert(Int, len[]))
    finalizer(relinquishmemory, blob)
    blob
end

function pingimage(wand::MagickWand, filename::AbstractString)
    status = ccall((:MagickPingImage, libwand), Cint, (Ptr{Cvoid}, Ptr{UInt8}), wand, filename)
    status == 0 && error(wand)
    nothing
end

function readimage(wand::MagickWand, filename::AbstractString)
    status = ccall((:MagickReadImage, libwand), Cint, (Ptr{Cvoid}, Ptr{UInt8}), wand, filename)
    status == 0 && error(wand)
    nothing
end

function readimage(wand::MagickWand, stream::IO)
    status = ccall((:MagickReadImageFile, libwand), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), wand, Libc.FILE(stream).ptr)
    status == 0 && error(wand)
    nothing
end

function readimage(wand::MagickWand, stream::Vector{UInt8})
    status = ccall((:MagickReadImageBlob, libwand), Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Cint), wand, stream, length(stream)*sizeof(eltype(stream)))
    status == 0 && error(wand)
    nothing
end

function writeimage(wand::MagickWand, filename::AbstractString)
    status = ccall((:MagickWriteImages, libwand), Cint, (Ptr{Cvoid}, Ptr{UInt8}, Cint), wand, filename, true)
    status == 0 && error(wand)
    nothing
end

function writeimage(wand::MagickWand, stream::IO)
    status = ccall((:MagickWriteImagesFile, libwand), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), wand, Libc.FILE(stream).ptr)
    status == 0 && error(wand)
    nothing
end

function size(wand::MagickWand)
    height = ccall((:MagickGetImageHeight, libwand), Csize_t, (Ptr{Cvoid},), wand)
    width = ccall((:MagickGetImageWidth, libwand), Csize_t, (Ptr{Cvoid},), wand)
    return convert(Int, width), convert(Int, height)
end

getnumberimages(wand::MagickWand) = convert(Int, ccall((:MagickGetNumberImages, libwand), Csize_t, (Ptr{Cvoid},), wand))

nextimage(wand::MagickWand) = ccall((:MagickNextImage, libwand), Cint, (Ptr{Cvoid},), wand) == 1

resetiterator(wand::MagickWand) = ccall((:MagickResetIterator, libwand), Cvoid, (Ptr{Cvoid},), wand)

newimage(wand::MagickWand, cols::Integer, rows::Integer, pw::PixelWand) = ccall((:MagickNewImage, libwand), Cint, (Ptr{Cvoid}, Csize_t, Csize_t, Ptr{Cvoid}), wand, cols, rows, pw.ptr) == 0 && error(wand)

# test whether image has an alpha channel
getimagealphachannel(wand::MagickWand) = ccall((:MagickGetImageAlphaChannel, libwand), Cint, (Ptr{Cvoid},), wand) == 1

function getimageproperties(wand::MagickWand,patt::AbstractString)
    numbProp = Ref{Csize_t}(0)
    p = ccall((:MagickGetImageProperties, libwand), Ptr{Ptr{UInt8}}, (Ptr{Cvoid}, Ptr{UInt8}, Ptr{Csize_t}), wand, patt, numbProp)
    if p == C_NULL
        error("Pattern not in property names")
    else
        nP = convert(Int, numbProp[])
        ret = Vector{String}(undef, nP)
        for i = 1:nP
            ret[i] = unsafe_string(unsafe_load(p,i))
        end
        ret
    end
end

function getimageproperty(wand::MagickWand, prop::AbstractString, warnuser::Bool=true)
    p = ccall((:MagickGetImageProperty, libwand), Ptr{UInt8}, (Ptr{Cvoid}, Ptr{UInt8}), wand, prop)
    if p == convert(Ptr{UInt8}, C_NULL)
        if warnuser
            possib = getimageproperties(wand,"*")
            @warn("Undefined property, possible names are \"$(join(possib,"\",\""))\"")
        end
        nothing
    else
        unsafe_string(p)
    end
end

# # get number of colors in the image
# magickgetimagecolors(wand::MagickWand) = ccall((:MagickGetImageColors, libwand), Csize_t, (Ptr{Cvoid},), wand)

# get the type
function getimagetype(wand::MagickWand)
    t = ccall((:MagickGetImageType, libwand), Cint, (Ptr{Cvoid},), wand)
    # Apparently the following is necessary, because the type is "potential"
    ccall((:MagickSetImageType, libwand), Cvoid, (Ptr{Cvoid}, Cint), wand, t)
    1 <= t <= length(IMType) || error("Image type ", t, " not recognized")
    IMType[t]
end

# get the colorspace
function getimagecolorspace(wand::MagickWand)
    cs = ccall((:MagickGetImageColorspace, libwand), Cint, (Ptr{Cvoid},), wand)
    1 <= cs <= length(IMColorspace) || error("Colorspace ", cs, " not recognized")
    IMColorspace[cs]
end

function setimagecolorspace(wand::MagickWand, cs::String)
    status = ccall((:MagickSetImageColorspace, libwand), Cint, (Ptr{Cvoid}, Cint), wand, IMColordict[cs])
    status == 0 && error(wand)
    nothing
end

imtype(buffer, cs) = IMTypedict[CStoIMTypedict[cs]]
imtype(buffer::AbstractArray{Bool}, cs) = IMTypedict["BilevelType"]

function setimagetype(wand::MagickWand, buffer, cs::String)
    status = ccall((:MagickSetImageType, libwand), Cint, (Ptr{Cvoid}, Cint), wand, imtype(buffer, cs))
    status == 0 && error(wand)
    nothing
end

# set the compression
function setimagecompression(wand::MagickWand, compression::Integer)
    status = ccall((:MagickSetImageCompression, libwand), Cint, (Ptr{Cvoid}, Cint), wand, Int32(compression))
    status == 0 && error(wand)
    nothing
end

function setimagecompressionquality(wand::MagickWand, quality::Integer)
    0 < quality <= 100 || error("quality setting must be in the (inclusive) range 1-100.\nSee http://www.imagemagick.org/script/command-line-options.php#quality for details")
    status = ccall((:MagickSetImageCompressionQuality, libwand), Cint, (Ptr{Cvoid}, Cint), wand, quality)
    status == 0 && error(wand)
    nothing
end

# set fps (for GIF-type images)
function getimagetickspersecond(wand::MagickWand)
    ccall((:MagickGetImageTicksPerSecond, libwand), Cint, (Ptr{Cvoid},), wand)
end

function getimagedelay(wand::MagickWand)
    ccall((:MagickGetImageDelay, libwand), Cint, (Ptr{Cvoid},), wand)
end

function setimagedelay(wand::MagickWand, fps)
    tps = getimagetickspersecond(wand)
    delay = round(Int, tps/fps)
    for i = 1:getnumberimages(wand)+1  # not clear why +1
        status = ccall((:MagickSetImageDelay, libwand), Cint, (Ptr{Cvoid}, Csize_t), wand, delay)
        status == 0 && error(wand)
        nextimage(wand)
    end
    resetiterator(wand)
    nothing
end

# set the image format
function setimageformat(wand::MagickWand, format::String)
    status = ccall((:MagickSetImageFormat, libwand), Cint, (Ptr{Cvoid}, Ptr{UInt8}), wand, format)
    status == 0 && error(wand)
    nothing
end

# get the pixel depth
getimagedepth(wand::MagickWand) = convert(Int, ccall((:MagickGetImageDepth, libwand), Csize_t, (Ptr{Cvoid},), wand))

# pixel depth for given channel type
getimagechanneldepth(wand::MagickWand, channelType::ChannelType) = convert(Int, ccall((:MagickGetImageChannelDepth, libwand), Csize_t, (Ptr{Cvoid}, UInt32), wand, channelType.value))

pixelsetcolor(wand::PixelWand, colorstr::String) = ccall((:PixelSetColor, libwand), Csize_t, (Ptr{Cvoid}, Ptr{UInt8}), wand, colorstr) == 0 && error(wand)

relinquishmemory(p) = ccall((:MagickRelinquishMemory, libwand), Ptr{UInt8}, (Ptr{UInt8},), p)

# get library information
# If you pass in "*", you get the full list of options
function queryoptions(pattern::AbstractString)
    nops = Ref{Cint}()
    pops = ccall((:MagickQueryConfigureOptions, libwand), Ptr{Ptr{UInt8}}, (Ptr{UInt8}, Ptr{Cint}), pattern, nops)
    ret = Vector{String}(nops[])
    for i = 1:nops[]
        ret[i] = unsafe_string(unsafe_load(pops, i))
    end
    ret
end

# queries the value of a particular option
function queryoption(option::AbstractString)
    p = ccall((:MagickQueryConfigureOption, libwand), Ptr{UInt8}, (Ptr{UInt8},), option)
    unsafe_string(p)
end


text:
using Statistics, Plots
# two stage, 500 epochs.
time_500 = [0.990310 , 0.973435 ,1.085917 ,1.083980 , 1.064565 ]
print(mean(time_500))
# two stage, training time 500,1000,...,5000 epochs
x = [500,1000,2500,5000]
y = [1.0396414, 2.041506, 4.725699, 8.181069]
scatter(x,y,label = "Training time [s] for number of epochs")
savefig("training_time.pdf")

# difference: 9 to 10 epochs -> one training epoch.
data_sizes = [30, 150, 300,1000]
two_stage_time = [time_30_2, time_150_2, time_300_2,time_1000_2]
l2_time =  [time_30, time_150, time_300,time_1000]

time_30 = 11.8 - 11.3
time_150 = 51.089 - 49.663
time_300 = 03.083- 00.793
time_1000 = 17.73 - 11.828
time_30_2 = 40.305 - 40.303
time_150_2 = 19.282- 19.271
time_300_2 = 04.662- 04.638
time_1000_2 = 28.201 - 28.047

scatter(data_sizes, two_stage_time,label = "data size with time per epoch: two stage")
scatter!(data_sizes, l2_time,label="data size with time per epoch:L2")
savefig("data_size.pdf")


text:
"""
    _include_my_codes(base_path::String, code_file_array::Array{String,1})

Function to include Julia files in a directory
"""
function _include_my_codes(base_path::String, code_file_array::Array{String,1})
    for code_file in code_file_array
        path_to_code = joinpath(base_path,code_file)
        include(path_to_code)
    end
end

# define constants here -
const alphavantage_api_url_string = "https://www.alphavantage.co/query"
const path_to_package = dirname(pathof(@__MODULE__))
const _PATH_TO_BASE = joinpath(path_to_package,"base")
const _PATH_TO_STS = joinpath(path_to_package, "sts")
const _PATH_TO_TI = joinpath(path_to_package, "ti")
const _PATH_TO_CRYPTO = joinpath(path_to_package, "crypto")
const _PATH_TO_FUNDAMENTALS = joinpath(path_to_package, "fundamentals")

# load offical packages here -
using DataFrames
using CSV
using HTTP
using JSON
using Dates
using Logging
using Reexport
@reexport using PooksoftBase

# need to update this syntax -
# load my base codes -
my_base_codes = [
    "Types.jl", "Network.jl", "User.jl", "Checks.jl",
    "Log.jl", "Handlers.jl", "Filesystem.jl", "Datastore.jl", "General.jl"
];
_include_my_codes(_PATH_TO_BASE, my_base_codes)

# stock time series -
my_sts_codes = [
    "STSDaily.jl", "STSWeekly.jl", "STSMonthly.jl", "Quote.jl", "Search.jl", "STSIntraday.jl"
]
_include_my_codes(_PATH_TO_STS, my_sts_codes)

# technical indicators -
my_ti_codes = [
    "SMA.jl", "EMA.jl", "RSI.jl"
]
_include_my_codes(_PATH_TO_TI, my_ti_codes)

# fundementals -
my_fundamental_codes = [
    "Overview.jl", "Earnings.jl", "Income.jl"
]
_include_my_codes(_PATH_TO_FUNDAMENTALS, my_fundamental_codes)


text:
using Utilities,AReg,StaticNets,JLD,MLBase,StatsBase
using PyCall; pygui(:qt); using PyPlot

## Load data
#half period
loadFileName = "/home/Domenico/Dropbox/Dynamic_Networks/data/emid_data/juliaFiles/Weekly_eMid_Data_from_2012_03_12_to_2015_02_27.jld"
#whole period
#loadFileName = "/home/Domenico/Dropbox/Dynamic_Networks/data/emid_data/juliaFiles/Weekly_eMid_Data_from_2009_06_22_to_2015_02_27.jld"
@load(loadFileName, AeMidWeekly_T,banksIDs,inactiveBanks, YeMidWeekly_T,weekInd,datesONeMid ,degsIO_T,strIO_T)

## Methods definitions

## Application to Data
N2,T = size(degsIO_T);N = round(Int,N2/2)
#estimate different versions
snapModData = SnapSeqNetDirBin1(degsIO_T)
singSnapEstNoIdent_T = estimate(snapModData)
singSnapEstIdent_T = estimate(snapModData; identPost = true,identIter=false)
singSnapEstIdentIter_T = estimate(snapModData; identPost = false,identIter=true)
# plot(singSnapEstNoIdent_T[1:30,:]')
# plot(singSnapEstIdent_T[1:3,:]')
# plot(singSnapEstIdentIter_T[1:3,:]')

## check that equal degrees are associated to equal fitnesses.
    # singSnapEst_T = singSnapEstIdent_T
    # fit,degs = [singSnapEst_T[1:N,1] singSnapEst_T[1+N:2N,1]],[degsIO_T[1:N,1] degsIO_T[1+N:2N,1]]
    # [unique(degs,1) unique(fit,1)]
    # # actually is a pair (InDeg,OutDeg) that is associated with a pair (InFit,OutFit). That is because
    # # there cannot be single loops. Hence the value of a fitness has to take into account
    # # that links cannot come from the node itself. So all things equal a node with
    # # (1,1) should have an higher  InFit than a (1,0) because it needs to attract more
    # # from the other ones
    # ##
    # D = maximum(degsIO_T)
    # storeFitDistrPerDeg =  Matrix{Vector{Float64}}(D,T)
    # meansFit = Matrix{Float64}(D,T)
    # stdFit = Matrix{Float64}(D,T)
    # testSingSnapEst_T = singSnapEst_T
    # t=1;d=1
    # for t=1:T
    #     degsIOt = degsIO_T[:,t]
    #     degsUniIOt = sort(unique(degsIOt))
    #     degsUniIOt = degsUniIOt[2:end]
    #     println(length(degsUniIOt))
    #     for d=1:length(degsUniIOt)
    #         deg = degsUniIOt[d]
    #         indsDeg = degsIOt[1:N].== deg
    #         storeFitDistrPerDeg[d,t] = testSingSnapEst_T[1:N,t][indsDeg]
    #         meansFit[d,t] = mean(testSingSnapEst_T[1:N,t][indsDeg])
    #         stdFit[d,t] = std(testSingSnapEst_T[1:N,t][indsDeg])
    #     end
    #
    # end
    #
    # deg = 1:4
    # plot(meansFit[deg,:]',"-");plot(meansFit[deg,:]' + stdFit[deg,:]',"--");plot(meansFit[deg,:]' - stdFit[deg,:]',"--")
    # tmp = 1./sumSq(degsIO_T,1)
    # cor(meansFit[2,:],tmp)
    # mean(meansFit[4,:])
    # #plot(tmp)
    # # the mean value of the fitnesses associated to a given degree is highly correlated
    # # with the inverse of total number of links

## lets estimate some AR(1)
forecastEvalAR1Net(singSnapEstIdent_T,AeMidWeekly_T)
#Compute rocCurve for each method
rocCurve(singSnapEstIdent_T,AeMidWeekly_T)
rocCurve(singSnapEstNoIdent_T,AeMidWeekly_T)
rocCurve(singSnapEstIdentIter_T,AeMidWeekly_T)

## Test that the naive model is Biased in the estimate of B (as Piero mentioned)
#estimate one AR1 for each fitness
using AReg,RCall
testSSFitSeq = singSnapEstIdent_T
(estAR,isCon) = estManyAR(testSSFitSeq)
# check results with R estimates
testSSFitSeqNnc =  testSSFitSeq[.!isCon,:]
# run estimates in R
@rput testSSFitSeqNnc
R"estAR_Rnnc = matrix(,nrow = dim(testSSFitSeqNnc)[1],ncol = 2)"
R"estARvar_Rnnc = matrix(,nrow = dim(testSSFitSeqNnc)[1],ncol = 1)"
R"for(n in 1:dim(testSSFitSeqNnc)[1]){ x = ar(testSSFitSeqNnc[n,],order.max = 1,demean = TRUE); estAR_Rnnc[n,] = c(x$x.mean,x$ar); estARvar_Rnnc[n] = c(x$var.pred) }"
@rget estAR_Rnnc estARvar_Rnnc
estAR_Rnnc[:,1] = estAR_Rnnc[:,1].*(1-estAR_Rnnc[:,2]) # the R function used computes the mean i use w
estAR_R = zeros(estAR); estAR_R[.!isCon,:] = estAR_Rnnc
estARvar_R = zeros(N2); estARvar_R[.!isCon,:] = estARvar_Rnnc

diff = (estAR_R .-estAR)[.!isCon,:]
plot(diff[:,1])
#plt[:hist](diff)[:,2] # Histogram
plot(estAR[:,2])
indRprob = estAR_R[:,2].<0
plot(degsIO_T[indRprob,:])
sum(degsIO_T[indRprob,:].>0,2)
estAR_R[indRprob,2] = 0
estAR[indRprob,2] = 0
estAR[indRprob,1] = -25
indJprob = estAR[:,2].<0
(sum(indRprob), sum(indJprob))
estAR[indRprob & (.!indJprob),:]
rocCurve(testSSFitSeq;estAndInd = (estAR_R,isCon))
rocCurve(testSSFitSeq;estAndInd = (estAR,isCon))
# My AR estimates (OLS) are in line with the Yule Walker ones from R. There are a few cases that
# give problems to both. due to almost constant degrees and, as consequence, fit
# some persistency parameters are estimated negative. I choose to put such cases
# to B = 0  for the moment, because this helps the out of sample performances of
# the naive model (keep calling it naive but we'll see who wins...)

## Simulate and estimate the AR1 on single snapshots to test for bias in the estimates
testSSFitSeq = singSnapEstIdent_T
D,T = size(testSSFitSeq)
i=1
using AReg
#estimate simulated data
testSSFitReEst =  estimate(snapModData; identPost = true,identIter=false)

#
function simulFit(Model::SnapSeqNetDirBin1,dgpARparAndVar::Array{<:Real,2},T::Int)
    N2 = length(dgpARparAndVar[:,1])
    simFitT = zeros(N2,T)
    for i=1:N2 simFitT[i,:]= AReg.simulateAR1(dgpARparAndVar[i,1:2],dgpARparAndVar[i,3];T=T) end
    return simFitT
 end

reEstDgpARparAndVar = [estAR estARvar_R]; reEstDgpARparAndVar[:,3] = 1
Ntest = 10
storeReEstAR = zeros(N2,3,Ntest)
for j=1:Ntest
    #use as AR parameters those estimated on the data
    #simulate AR paths for fit
    fitSimulT = simulFit(fooSnapSeqNetDirBin1,reEstDgpARparAndVar,T)
    #plot(fitSimulT')
    #given fit_T simulate one network
    sampledA_T = sampl(fooSnapSeqNetDirBin1,fitSimulT)
    samplDegIO_T = [sumSq(sampledA_T,2); sumSq(sampledA_T,1)]
    singSnapReEstIdent_T = estimate(fooSnapSeqNetDirBin1;degsIO_T = samplDegIO_T,  identPost = true,identIter=true)

    #estimate(fooErgmDirBin1;degIO = samplDegIO_T[:,47],targetErr = 1e-4)#,identIter=identIter)
    #estimate the B
    (reEstAR,isCon) = estManyAR(singSnapReEstIdent_T)

    storeReEstAR[:,:,j] = [reEstAR isCon]
end
##

plt[:hist](diff) # Histogram
plot(diff)
diffMat = storeReEstAR[:,2,:] .-reEstDgpARparAndVar[:,2]

plot(reEstDgpARparAndVar[:,2],mean(diffMat,2),".")

# check Bias in my AR! estimates. There is not!!
Nsim = 10000
    arp = AReg.ARp(ones(1+1),.1,[])
    tmp = zeros(Nsim,2)
    parSim = [-1.5,0.85]
    for t=1:Nsim
        tmp[t,:] = AReg.fitARp(AReg.simulateAR1(parSim,20.0;T=T),arp).estimates .- parSim
    end
    plt[:hist](tmp[:,2],50)

#

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
# override some methods to use Plotlyjs/Blink

import Plotlyjs

function _create_plot(pkg::PlotlyPackage; kw...)
    d = Dict(kw)
    # TODO: create the window/canvas/context that is the plot within the backend (call it `o`)
    # TODO: initialize the plot... title, xlabel, bgcolor, etc
    o = Plotlyjs.Plot()

    Plot(o, pkg, 0, d, Dict[])
end

function _add_series(::PlotlyPackage, plt::Plot; kw...)
    d = Dict(kw)

    # add to the data array
    pdict = plotly_series(d)
    gt = Plotlyjs.GenericTrace(pdict[:type], pdict)
    push!(plt.o.data, gt)
    if !isnull(plt.o.window)
        Plotlyjs.addtraces!(plt.o, gt)
    end

    push!(plt.seriesargs, d)
    plt
end

# TODO: override this to update plot items (title, xlabel, etc) after creation
function _update_plot(plt::Plot{PlotlyPackage}, d::Dict)
    pdict = plotly_layout(d)
    plt.o.layout = Plotlyjs.Layout(pdict)
    if !isnull(plt.o.window)
        Plotlyjs.relayout!(plt.o, pdict...)
    end
end

function Base.display(::PlotsDisplay, plt::Plot{PlotlyPackage})
    dump(plt.o)
    show(plt.o)
end

function Base.display(::PlotsDisplay, plt::Subplot{PlotlyPackage})
    error()
end


text:
# slt_configload.jl
#
# Jeremy Rüffer
# Thünen Institut
# Institut für Agrarklimaschutz
# Junior Research Group NITROSPHERE
# Julia 0.7
# 09.12.2016
# Last Edit: 18.04.2019

"""# slt_configload

Load a list of .CFG files

`slt_configload(files::Array{T,1}) where T <: String`\n
* **files**::Array{String,1} = Array of CFG file names"""
function slt_configload(files::Array{T,1})::Array{Dict,1} where T <: String
	# Dict Prep
	Time = DateTime(9999)			# DateTime
	FileName = ""					# String
	Sonic = ""						# String
	Analyzer = ""					# String
	Sonic_Alignment = 0.0			# Float64
	Frequency = ""					# String
	Average_Time = 0.0				# Float64
	Station_Height = ""				# String
	Measurement_Height = ""			# String
	Vegetation_Height = ""			# String
	Inductance_CO2 = ""				# String
	Inductance_H2O = ""				# String
	Coordinate_Rotation = ""		# String
	Webb_Correction = false			# Bool
	Linear_Detrend = false			# Bool
	Analog_Inputs = zeros(Int,6)	# Array{Int,1}
	Analog_Names = Array{String}(undef,6)		# Array{String,1}
	Analog_Lower = zeros(Int,6)					# Array{Int,1}
	Analog_Upper = zeros(Int,6)					# Array{Int,1}
	Analog_Units = Array{String}(undef,6)		# Array{String,1}
	Analog_Delay = zeros(Int,6)					# Array{Int,1}
	Slope = zeros(Float64,6)			# Array{Float64,1}
	Analog_Count = 0					# Int

	configs = Array{Dict}(undef,length(files))

	# Fill Fields with Single Values
	for i=1:1:length(files)
		# Initialize empty arrays
		Analog_Inputs = zeros(Int,6)
		Analog_Names = Array{String}(undef,0)
		Analog_Lower = zeros(Int,6)
		Analog_Upper = zeros(Int,6)
		Analog_Units = Array{String}(undef,0)
		Analog_Delay = zeros(Int,6)
		Slope = zeros(Float64,6)

		FileName = files[i]
		Time = DateTime(Meta.parse(files[i][end-14:end-11])) + Dates.Day(Meta.parse(files[i][end-10:end-8])) - Dates.Day(1) + Dates.Hour(Meta.parse(files[i][end-7:end-6])) + Dates.Minute(Meta.parse(files[i][end-5:end-4]))

		# Load Config File
		analog = Array{String}(undef,6)
		fid = open(files[i],"r")
		try
			readline(fid)
			Sonic = readline(fid)[8:end]
			Analyzer = readline(fid)[11:end]
			Sonic_Alignment = Meta.parse(readline(fid)[18:end-4])
			Frequency = readline(fid)[19:end-3]
			Average_Time = Meta.parse(readline(fid)[15:end-4])
			analog[1] = readline(fid)
			analog[2] = readline(fid)
			analog[3] = readline(fid)
			analog[4] = readline(fid)
			analog[5] = readline(fid)
			analog[6] = readline(fid)
			Station_Height = readline(fid)[17:end-2]
			readline(fid)
			Measurement_Height = readline(fid)[21:end-2]
			Vegetation_Height = readline(fid)[20:end-2]
			Inductance_CO2 = readline(fid)[21:end-1]
			Inductance_H2O = readline(fid)[21:end-1]
			Coordinate_Rotation = readline(fid)[22:end]
			Webb_Correction = readline(fid)[18:end] == "yes"
			Linear_Detrend = readline(fid)[20:end] == "yes"
		catch e
			println("Problem loading configuration: " * files[1])
			println(e)
		end
		close(fid)

		# Parse Analog Inputs
		for j=1:1:length(analog)
			if analog[j][end] == 'E'
				f = readdlm(IOBuffer(analog[j]),',')

				Analog_Inputs[j] = j
				Analog_Names = [Analog_Names;String(strip(f[1]))]
				Analog_Lower[j] = f[2]
				Analog_Units = [Analog_Units;String(strip(f[4]))]
				Analog_Upper[j] = f[3]
				Analog_Delay[j] = f[5]
				Slope[j] = (f[3] - f[2])/5000
			end
		end
		Analog_Count = length(Analog_Names)

		configs[i] = Dict("Time"=>Time,"FileName"=>FileName,"Sonic"=>Sonic,"Analyzer"=>Analyzer,"Sonic_Alignment"=>Sonic_Alignment,"Frequency"=>Frequency,"Average_Time"=>Average_Time,"Station_Height"=>Station_Height,"Measurement_Height"=>Measurement_Height,"Vegetation_Height"=>Vegetation_Height,"Inductance_CO2"=>Inductance_CO2,"Inductance_H2O"=>Inductance_H2O,"Coordinate_Rotation"=>Coordinate_Rotation,"Webb_Correction"=>Webb_Correction,"Linear_Detrend"=>Linear_Detrend,"Analog_Inputs"=>Analog_Inputs,"Analog_Names"=>Analog_Names,"Analog_Lower"=>Analog_Lower,"Analog_Upper"=>Analog_Upper,"Analog_Units"=>Analog_Units,"Analog_Delay"=>Analog_Delay,"Slope"=>Slope,"Analog_Count"=>Analog_Count)
	end

	# Sort Config Files
	I = sortperm(get.(configs,"Time",DateTime(9999)))
	configs = configs[I]
	return configs
end # slt_configload(files::Array{T,1}) where T <: String

max_stars_repo_name: JeremyRueffer/ClimateDataIO

max_issues_repo_name: JeremyRueffer/ClimateDataIO


text:
module Market

import ..Clarus
function caplets(;params...)
  return Clarus.api_request("Market","Caplets";params...)
end

function curveusage(;params...)
  return Clarus.api_request("Market","CurveUsage";params...)
end

function df(;params...)
  return Clarus.api_request("Market","DF";params...)
end

function fixings(;params...)
  return Clarus.api_request("Market","Fixings";params...)
end

function futures(;params...)
  return Clarus.api_request("Market","Futures";params...)
end

function fxrates(;params...)
  return Clarus.api_request("Market","FXRates";params...)
end

function pardv01(;params...)
  return Clarus.api_request("Market","ParDV01";params...)
end

function parrates(;params...)
  return Clarus.api_request("Market","ParRates";params...)
end

function zerorates(;params...)
  return Clarus.api_request("Market","ZeroRates";params...)
end

end #MODULE-END


text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni
#

"""
    (renderer::OnOffRenderer)(ray::Ray) :: RGB{Float32}

Return the `renderer` color (default WHITE) if the ray hit an object,
else return the `renderer` background color (default BLACK).

See also: [`OnOffRenderer`](@ref), [`Ray`](@ref)
"""
function (renderer::OnOffRenderer)(ray::Ray)
    ray_intersection(renderer.world, ray) ≠ nothing ?
        renderer.color :
        renderer.background_color
end

"""
    (renderer::FlatRenderer)(ray::Ray) :: RGB{Float32}

Return the `renderer` background color (default BLACK) if ray
doesn't hit any objects, else evaluate the color depending on
the material and the self luminosity.

See also: [`FlatRenderer`](@ref), [`Ray`](@ref)
"""
function (renderer::FlatRenderer)(ray::Ray)
    hit = ray_intersection(renderer.world, ray)
    !(isnothing(hit)) || (return renderer.background_color)

    mat = hit.shape.Material
    col1 = get_color(mat.brdf.pigment, hit.surface_point)
    col2 = get_color(mat.emitted_radiance, hit.surface_point)

    return (col1 + col2)
end

"""
    (renderer::PathTracer)(ray::Ray) :: RGB{Float32}

A simple path-tracing renderer.

Given the `renderer` parameters (`world`, `background_color`, `pcg`,
`num_of_rays`, `max_depth` and `russian_roulette_limit`) evaluate the
radiance contribute for each secondary `Ray` according to them.

See also: [`PathTracer`](@ref), [`Ray`](@ref)
"""
function (renderer::PathTracer)(ray::Ray)
    !(ray.depth > renderer.max_depth) || (return BLACK)

    hit_record = ray_intersection(renderer.world, ray)
    !isnothing(hit_record) || (return renderer.background_color)

    hit_material = hit_record.shape.Material
    hit_surface_point = hit_record.surface_point
    hit_color = get_color(hit_material.brdf.pigment, hit_surface_point)
    emitted_radiance = get_color(hit_material.emitted_radiance, hit_surface_point)
    hit_color_lum = max(hit_color.r, hit_color.g, hit_color.b)

    # Russian Roulette
    if ray.depth >= renderer.russian_roulette_limit
        if random(renderer.pcg) > hit_color_lum
            # Keep the recursion going, but compensate
            # for other potentially discarded rays
            hit_color *= 1.0 / (1.0 - hit_color_lum)
        else
            # Terminate prematurely
            return emitted_radiance
        end
    end

    # Monte Carlo integration

    cum_radiance =  RGB{Float32}(0.0, 0.0, 0.0)
    num_of_rays = renderer.num_of_rays
    # Only do costly recursions if it's worth it
    if hit_color_lum > 0.0
        for ray_index in 1:num_of_rays
            new_ray = scatter_ray(
                            typeof(hit_material.brdf),
                            renderer.pcg,
                            hit_record.ray.dir,
                            hit_record.world_point,
                            hit_record.normal,
                            ray.depth + 1,
                        )

            # Recursive call
            new_radiance = renderer(new_ray)

            cum_radiance += hit_color * new_radiance
        end
    end

    return emitted_radiance + cum_radiance * (1.0 / num_of_rays)
end

function (renderer::PointLightRenderer)(ray::Ray)
    hit_record = ray_intersection(renderer.world, ray)
    !isnothing(hit_record) || (return renderer.background_color)

    hit_material = hit_record.shape.Material
    flag_bg = hit_record.shape.flag_background
    flag_pl = hit_record.shape.flag_pointlight
    result_color = renderer.ambient_color

    for cur_light in renderer.world.point_lights
        if is_point_visible(renderer.world, cur_light.position, hit_record.world_point) && flag_bg==false && flag_pl==false
            distance_vec = hit_record.world_point - cur_light.position
            distance = norm(distance_vec)
            in_dir = distance_vec * (1.0 / distance)
            #cos_theta = max(0.0, -normalize(ray.dir)⋅hit_record.normal)
            cos_theta = max(0.0, -normalize(distance_vec)⋅hit_record.normal)

            distance_factor =
                (cur_light.linear_radius > 0) ?
                (cur_light.linear_radius / distance)^2 :
                1.0

            emitted_color = get_color(hit_material.emitted_radiance, hit_record.surface_point)
            brdf_color = evaluate(
                            hit_material.brdf,
                            hit_record.normal,
                            in_dir,
                            -ray.dir,
                            hit_record.surface_point,
                        )
            result_color += (emitted_color + brdf_color) * cur_light.color * cos_theta * distance_factor

        elseif flag_pl==true && flag_bg==false
            distance_vec = hit_record.world_point - cur_light.position
            distance = norm(distance_vec)
            in_dir = distance_vec * (1.0 / distance)
            cos_theta = max(0.0, -normalize(ray.dir)⋅hit_record.normal)
            #cos_theta = max(0.0, -normalize(distance_vec)⋅hit_record.normal)

            distance_factor =
                (cur_light.linear_radius > 0) ?
                (cur_light.linear_radius / distance)^2 :
                1.0

            emitted_color = get_color(hit_material.emitted_radiance, hit_record.surface_point)
            brdf_color = evaluate(
                            hit_material.brdf,
                            hit_record.normal,
                            in_dir,
                            -ray.dir,
                            hit_record.surface_point,
                        )
            result_color += (emitted_color + brdf_color) * cur_light.color * cos_theta * distance_factor

        elseif flag_bg==true
            result_color += get_color(hit_material.brdf.pigment, hit_record.surface_point)

        else
            result_color += renderer.dark_parameter*get_color(hit_material.brdf.pigment, hit_record.surface_point)
        end
    end

    return result_color
end


text:
using ModelingToolkit
using LightGraphs
using DiffEqBase
using Test
using UnPack

# Define some variables
@parameters t L g
@variables x(t) y(t) w(t) z(t) T(t) xˍt(t) yˍt(t) xˍˍt(t) yˍˍt(t)
D = Differential(t)

eqs2 = [D(D(x)) ~ T*x,
        D(D(y)) ~ T*y - g,
        0 ~ x^2 + y^2 - L^2]
pendulum2 = ODESystem(eqs2, t, [x, y, T], [L, g], name=:pendulum)
lowered_sys = ModelingToolkit.ode_order_lowering(pendulum2)

lowered_eqs = [D(xˍt) ~ T*x,
               D(yˍt) ~ T*y - g,
               D(x) ~ xˍt,
               D(y) ~ yˍt,
               0 ~ x^2 + y^2 - L^2,]
@test ODESystem(lowered_eqs, t, [xˍt, yˍt, x, y, T], [L, g], name=:pendulum) == lowered_sys
@test isequal(equations(lowered_sys), lowered_eqs)

# Simple pendulum in cartesian coordinates
eqs = [D(x) ~ w,
       D(y) ~ z,
       D(w) ~ T*x,
       D(z) ~ T*y - g,
       0 ~ x^2 + y^2 - L^2]
pendulum = ODESystem(eqs, t, [x, y, w, z, T], [L, g], name=:pendulum)

pendulum = initialize_system_structure(pendulum)
sss = structure(pendulum)
@unpack graph, fullvars, varassoc = sss
@test StructuralTransformations.matching(sss, varassoc .== 0) == map(x -> x == 0 ? StructuralTransformations.UNASSIGNED : x, [1, 2, 3, 4, 0, 0, 0, 0, 0])

sys, assign, eqassoc = StructuralTransformations.pantelides!(pendulum)
sss = structure(sys)
@unpack graph, fullvars, varassoc = sss
scc = StructuralTransformations.find_scc(graph, assign)
@test sort(sort.(scc)) == [
                           [1],
                           [2],
                           [3, 4, 7, 8, 9],
                           [5],
                           [6],
                          ]

@test graph.fadjlist == [[1, 7], [2, 8], [3, 5, 9], [4, 6, 9], [5, 6], [1, 2, 5, 6], [1, 3, 7, 10], [2, 4, 8, 11], [1, 2, 5, 6, 10, 11]]
@test varassoc == [10, 11, 0, 0, 1, 2, 3, 4, 0, 0, 0]
#1: D(x) ~ w
#2: D(y) ~ z
#3: D(w) ~ T*x
#4: D(z) ~ T*y - g
#5: 0 ~ x^2 + y^2 - L^2
# ----
#6: D(eq:5) -> 0 ~ 2xx'+ 2yy'
#7: D(eq:1) -> D(D(x)) ~ D(w) -> D(xˍt) ~ D(w) -> D(xˍt) ~ T*x
#8: D(eq:2) -> D(D(y)) ~ D(z) -> D(y_t) ~ T*y - g
#9: D(eq:6) -> 0 ~ 2xx'' + 2x'x' + 2yy'' + 2y'y'
#                 [1, 2, 3, 4, 5, 6, 7, 8, 9]
@test eqassoc == [7, 8, 0, 0, 6, 9, 0, 0, 0]

using ModelingToolkit
@parameters t L g
@variables x(t) y(t) w(t) z(t) T(t) xˍt(t) yˍt(t)
D = Differential(t)
idx1_pendulum = [D(x) ~ w,
                 D(y) ~ z,
                 #0 ~ x^2 + y^2 - L^2,
                 D(w) ~ T*x,
                 D(z) ~ T*y - g,
                 # intermediate 1: 0 ~ 2x*D(x) + 2y*D(y) - 0,
                 # intermediate 2(a): 0 ~ 2x*w + 2y*z - 0, (substitute D(x) and D(y))
                 #0 ~ 2x*w + 2y*z,
                 # D(D(x)) ~ D(w) and substitute the rhs
                 D(xˍt) ~ T*x,
                 # D(D(y)) ~ D(z) and substitute the rhs
                 D(yˍt) ~ T*y - g,
                 # 2x*D(D(x)) + 2*D(x)*D(x) + 2y*D(D(y)) + 2*D(y)*D(y) and
                 # substitute the rhs
                 0 ~ 2x*(T*x) + 2*xˍt*xˍt + 2y*(T*y - g) + 2*yˍt*yˍt]
@named idx1_pendulum = ODESystem(idx1_pendulum, t, [x, y, w, z, xˍt, yˍt, T], [L, g])
first_order_idx1_pendulum = ode_order_lowering(idx1_pendulum)

using OrdinaryDiffEq
using LinearAlgebra
prob = ODEProblem(ODEFunction(first_order_idx1_pendulum),
        #  [x, y, w, z, xˍt, yˍt, T]
           [1, 0, 0, 0,   0,   0, 0.0],# 0, 0, 0, 0],
           (0, 10.0),
           [1, 9.8],
           mass_matrix=calculate_massmatrix(first_order_idx1_pendulum))
sol = solve(prob, Rodas5());
#plot(sol, vars=(1, 2))

new_sys = dae_index_lowering(ModelingToolkit.ode_order_lowering(pendulum2))

prob_auto = ODEProblem(new_sys,
                       [D(x)=>0,
                        D(y)=>0,
                        x=>1,
                        y=>0,
                        T=>0.0],
                       (0, 100.0),
                       [1, 9.8])
sol = solve(prob_auto, Rodas5());
#plot(sol, vars=(x, y))

# Define some variables
@parameters t L g
@variables x(t) y(t) T(t)
D = Differential(t)

eqs2 = [D(D(x)) ~ T*x,
        D(D(y)) ~ T*y - g,
        0 ~ x^2 + y^2 - L^2]
pendulum2 = ODESystem(eqs2, t, [x, y, T], [L, g], name=:pendulum)

# Turn into a first order differential equation system
first_order_sys = ModelingToolkit.ode_order_lowering(pendulum2)

# Perform index reduction to get an Index 1 DAE
new_sys = dae_index_lowering(first_order_sys)

u0 = [
  D(x)    => 0.0,
  D(y)    => 0.0,
  x       => 1.0,
  y       => 0.0,
  T       => 0.0
]

p = [
    L => 1.0,
    g => 9.8
]

prob_auto = ODEProblem(new_sys,u0,(0.0,10.0),p)
sol = solve(prob_auto, Rodas5());
#plot(sol, vars=(D(x), y))

###
### More BLT/SCC tests
###

# Test Tarjan (1972) Fig. 3
g = [
     [2],
     [3,8],
     [4,7],
     [5],
     [3,6],
     Int[],
     [4,6],
     [1,7],
    ]
graph = StructuralTransformations.BipartiteGraph(8, 8)
for (eq, vars) in enumerate(g), var in vars
    add_edge!(graph, eq, var)
end
scc = StructuralTransformations.find_scc(graph)
@test scc == [
              [6],
              [3, 4, 5, 7],
              [1, 2, 8],
             ]


text:
module MITprof

using NetCDF, Dates

"""
    read(f::String="MITprof/MITprof_mar2016_argo9506.nc")

Standard Depth Argo Data Example.

Here we read the `MITprof` standard depth data set from `https://doi.org/10.7910/DVN/EE3C40`
For more information, please refer to Forget, et al 2015 (`http://dx.doi.org/10.5194/gmd-8-3071-2015`)
The produced figure shows the number of profiles as function of time for a chosen file
    and maps out the locations of Argo profiles collected for a chosen year.

```
using OceanRobots, Plots

fi="MITprof/MITprof_mar2016_argo9506.nc"
(lo,la,ye)=MITprof.read(fi)

h = histogram(ye,bins=20,label=fi[end-10:end],title="Argo profiles")

ye0=2004; ye1=ye0.+1
kk=findall((ye.>ye0) .* (ye.<ye1))
scatter(lo[kk],la[kk],label=fi[end-10:end],title="Argo profiles count")
```
"""
function read(f::String="MITprof/MITprof_mar2016_argo9506.nc")
    #i = ncinfo(f)
    lo = ncread(f, "prof_lon")
    la = ncread(f, "prof_lat")
    x = ncread(f, "prof_date")
    t = julian2datetime.( datetime2julian(DateTime(0)) .+ x )
    ye = year.(t) + dayofyear.(t) ./ 365.0 #neglecting leap years ...
    return Float64.(lo),Float64.(la),Float64.(ye)
end

"""
    loop(pth::String="profiles/")

Standard Depth Argo Data Collection -- see `?MITprof.read` for detail.
"""
function loop(pth::String="profiles/")
    λ=("MITprof_mar2016_argo9506.nc","MITprof_mar2016_argo0708.nc",
    "MITprof_mar2016_argo0910.nc","MITprof_mar2016_argo1112.nc",
    "MITprof_mar2016_argo1314.nc","MITprof_mar2016_argo1515.nc")

    lo=[]; la=[]; ye=[];
    for f=1:length(λ)
        (tmplo,tmpla,tmpye)=read(pth*λ[f])
        lo=cat(lo,tmplo,dims=1)
        la=cat(la,tmpla,dims=1)
        ye=cat(ye,tmpye,dims=1)
    end

    return Float64.(lo),Float64.(la),Float64.(ye)
end

end


text:
module Sparse23push
using SparseArrays
export PreLapl2,PreLapl3
function PreLapl2(size::Int)
    ind(i::Int,j::Int)=(j-1)*size+i
    I=Int[]
    J=Int[]
    V=Float64[]
    h=1.0/(size-1)
    size2=size*size
    h2=h*h
    cd=-4.0/h2
    hd=1.0/h2
    nc=0
    for i=1:size
        for j=1:size
            ii=ind(i,j)
            push!(I,ii)
            push!(J,ii)
            push!(V,cd)
            nc+=1
            for i1 in [-1,1]
                if i+i1>0 && i+i1<=size
                    ii1=ind(i+i1,j)
                    push!(I,ii)
                    push!(J,ii1)
                    push!(V,hd)
                    nc+=1
                end
                if j+i1>0 && j+i1<=size
                    ii1=ind(i,j+i1)
                    push!(I,ii)
                    push!(J,ii1)
                    push!(V,hd)
                    nc+=1
                end
            end

        end
    end

    M=sparse(I,J,V)
    return M,size*size,nc
end
function PreLapl3(size)
    size2=size*size
    ind(i::Int,j::Int,k::Int)=(k-1)*size2+(j-1)*size+i
    I=Int[]
    J=Int[]
    V=Float64[]
    h=1.0/(size-1)
    size2=size*size
    h2=h*h
    cd=-6.0/h2
    hd=1.0/h2
    nc=0
    for i=1:size
        for j=1:size
            for k=1:size
                ii=ind(i,j,k)
                push!(I,ii)
                push!(J,ii)
                push!(V,cd)
                nc+=1
                for i1 in [-1,1]
                    if i+i1>0 && i+i1<=size
                        ii1=ind(i+i1,j,k)
                        push!(I,ii)
                        push!(J,ii1)
                        push!(V,hd)
                        nc+=1
                    end
                   if j+i1>0 && j+i1<=size
                       ii1=ind(i,j+i1,k)
                       push!(I,ii)
                       push!(J,ii1)
                       push!(V,hd)
                       nc+=1
                   end
                    if k+i1>0 && k+i1<=size
                        ii1=ind(i,j,k+i1)
                        push!(I,ii)
                        push!(J,ii1)
                        push!(V,hd)
                        nc+=1
                    end
                end
            end
        end

    end
    M=sparse(I,J,V)
    return M,size^3,nc
end
end

max_stars_repo_name: danielcort/benchmarks-python-julia-c

max_issues_repo_name: danielcort/benchmarks-python-julia-c

max_forks_repo_name: danielcort/benchmarks-python-julia-c


text:
using Test
import MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities

mock = MOIU.MockOptimizer(MOIU.Model{Float64}())
config = MOIT.TestConfig()

@testset "Linear" begin
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1.0, 0.0, 2.0],
                          (MOI.VectorAffineFunction{Float64}, MOI.Zeros)        => [[-3, -1]])
    MOIT.lin1vtest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1.0, 0.0, 2.0],
                          (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) => [[0, 2, 0]],
                          (MOI.VectorAffineFunction{Float64}, MOI.Zeros)        => [[-3, -1]])
    MOIT.lin1ftest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [-4, -3, 16, 0],
                          (MOI.VectorAffineFunction{Float64}, MOI.Zeros)        => [[7, 2, -4]])
    MOIT.lin2vtest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [-4, -3, 16, 0],
                          (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) => [[0]],
                          (MOI.VectorAffineFunction{Float64}, MOI.Nonpositives) => [[0]],
                          (MOI.VectorAffineFunction{Float64}, MOI.Zeros)        => [[7, 2, -4], [7]])
    MOIT.lin2ftest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, MOI.INFEASIBLE, MOI.INFEASIBLE_POINT, MOI.INFEASIBILITY_CERTIFICATE)
    MOIT.lin3test(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, MOI.INFEASIBLE)
    MOIT.lin3test(mock, MOIT.TestConfig(infeas_certificates=false))
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, MOI.INFEASIBLE_OR_UNBOUNDED)
    MOIT.lin3test(mock, MOIT.TestConfig(infeas_certificates=false))
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, MOI.INFEASIBLE, MOI.INFEASIBLE_POINT, MOI.INFEASIBILITY_CERTIFICATE)
    MOIT.lin4test(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, MOI.INFEASIBLE)
    MOIT.lin4test(mock, MOIT.TestConfig(infeas_certificates=false))
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, MOI.INFEASIBLE_OR_UNBOUNDED)
    MOIT.lin4test(mock, MOIT.TestConfig(infeas_certificates=false))
end
@testset "NormInf" begin
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock,
        [1, 0.5, 1],
        (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[-1], [-1]],
        (MOI.VectorOfVariables, MOI.NormInfinityCone) => [[1, 0, -1]])
    MOIT.norminf1vtest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock,
        [1, 0.5, 1],
        (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[-1], [-1]],
        (MOI.VectorAffineFunction{Float64}, MOI.NormInfinityCone) => [[1, 0, -1]])
    MOIT.norminf1ftest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, MOI.INFEASIBLE, MOI.INFEASIBLE_POINT, MOI.INFEASIBILITY_CERTIFICATE)
    MOIT.norminf2test(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock,
        [2, -1, -1, -1],
        (MOI.VectorAffineFunction{Float64}, MOI.NormInfinityCone) => [vcat(1, fill(-inv(3), 3))],
        (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) => [fill(inv(3), 3)])
    MOIT.norminf3test(mock, config)
end
@testset "NormOne" begin
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock,
        [1, 0.5, 0.5],
        (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[-1], [0]],
        (MOI.VectorOfVariables, MOI.NormOneCone) => [[1, -1, -1]])
    MOIT.normone1vtest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock,
        [1, 0.5, 0.5],
        (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[-1], [0]],
        (MOI.VectorAffineFunction{Float64}, MOI.NormOneCone) => [[1, -1, -1]])
    MOIT.normone1ftest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, MOI.INFEASIBLE, MOI.INFEASIBLE_POINT, MOI.INFEASIBILITY_CERTIFICATE)
    MOIT.normone2test(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock,
        [4, -1, -1, -1],
        (MOI.VectorAffineFunction{Float64}, MOI.NormOneCone) => [vcat(1, fill(-1, 3))],
        (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) => [ones(3)])
    MOIT.normone3test(mock, config)
end
@testset "SOC" begin
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1.0, 1/√2, 1/√2],
                          (MOI.VectorAffineFunction{Float64}, MOI.Zeros)           => [[-√2]])
    MOIT.soc1vtest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1.0, 1/√2, 1/√2],
                          (MOI.VectorAffineFunction{Float64}, MOI.SecondOrderCone) => [[√2, -1, -1]],
                          (MOI.VectorAffineFunction{Float64}, MOI.Zeros)           => [[-√2]])
    MOIT.soc1ftest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [-1/√2, 1/√2, 1.],
                          (MOI.VectorAffineFunction{Float64}, MOI.SecondOrderCone) => [[√2, 1, -1]],
                          (MOI.VectorAffineFunction{Float64}, MOI.Zeros)           => [[√2]],
                          (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives)    => [[1.0]])
    MOIT.soc2ntest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [-1/√2, 1/√2, 1.],
                          (MOI.VectorAffineFunction{Float64}, MOI.SecondOrderCone) => [[√2, 1, -1]],
                          (MOI.VectorAffineFunction{Float64}, MOI.Zeros)           => [[√2]],
                          (MOI.VectorAffineFunction{Float64}, MOI.Nonpositives)    => [[-1.0]])
    MOIT.soc2ptest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, MOI.INFEASIBLE, MOI.INFEASIBLE_POINT, MOI.INFEASIBILITY_CERTIFICATE)
    MOIT.soc3test(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1.0, 2/√5, 1/√5, 2/√5, 1/√5],
                          (MOI.VectorAffineFunction{Float64}, MOI.Zeros)           => [[-√5, -2.0, -1.0]])
    MOIT.soc4test(mock, config)
end
@testset "RSOC" begin
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [0.5, 1.0, 1/√2, 1/√2],
                          (MOI.SingleVariable,                MOI.EqualTo{Float64})       => [-√2, -1/√2],
                          (MOI.VectorOfVariables,             MOI.RotatedSecondOrderCone) => [[√2, 1/√2, -1.0, -1.0]])
    # double variable bounds on a and b variables
    mock.eval_variable_constraint_dual = false
    MOIT.rotatedsoc1vtest(mock, config)
    mock.eval_variable_constraint_dual = true
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1/√2, 1/√2],
                          (MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone) => [[√2, 1/√2, -1.0, -1.0]])
    MOIT.rotatedsoc1ftest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, MOI.INFEASIBLE, tuple(),
                          (MOI.SingleVariable,                MOI.LessThan{Float64})      => [-1],
                          (MOI.SingleVariable,                MOI.EqualTo{Float64})       => [-1],
                          (MOI.SingleVariable,                MOI.GreaterThan{Float64})   => [1],
                          (MOI.VectorOfVariables            , MOI.RotatedSecondOrderCone) => [[1, 1, -1]])
    # double variable bounds on x, y, z variables
    mock.eval_variable_constraint_dual = false
    MOIT.rotatedsoc2test(mock, config)
    mock.eval_variable_constraint_dual = true
    n = 2
    ub = 3.0
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1.0; zeros(n-1); ub; √ub; ones(2)],
                          (MOI.SingleVariable,                MOI.EqualTo{Float64})       => [-√ub/4, -√ub/4],
                          (MOI.VectorOfVariables,             MOI.Nonnegatives)           => [zeros(n)],
                          (MOI.SingleVariable,                MOI.GreaterThan{Float64})   => [0.0],
                          (MOI.SingleVariable,                MOI.LessThan{Float64})      => [-1/(2*√ub)],
                          (MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone) => [[√ub/(2*√2); √ub/(2*√2); -√ub/2; zeros(n-1)], [√ub/√2, 1/√(2*ub), -1.0]])
    # double variable bounds on u
    mock.eval_variable_constraint_dual = false
    MOIT.rotatedsoc3test(mock, config)
    mock.eval_variable_constraint_dual = true

    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, ones(4),
        (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64})       => [-1.0],
        (MOI.VectorOfVariables, MOI.RotatedSecondOrderCone) => [[1.0, 1.0, -1.0, -1.0]])
    MOIT.rotatedsoc4test(mock, config)
end
@testset "GeoMean" begin
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, ones(4),
        (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [-inv(3)])
    MOIT.geomean1vtest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, ones(4),
        (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [-inv(3)],
        (MOI.VectorAffineFunction{Float64}, MOI.GeometricMeanCone) => [vcat(-1.0, fill(inv(3), 3))])
    MOIT.geomean1ftest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, ones(10),
        (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) => fill(-inv(9), 9))
    MOIT.geomean2vtest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, ones(10),
        (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) => fill(-inv(9), 9),
        (MOI.VectorAffineFunction{Float64}, MOI.GeometricMeanCone) => [vcat(-1.0, fill(inv(9), 9))])
    MOIT.geomean2ftest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [2.0, 2.0],
        (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [-2.0])
    MOIT.geomean3vtest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [2.0, 2.0],
        (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [-2.0],
        (MOI.VectorAffineFunction{Float64}, MOI.GeometricMeanCone) => [[-2.0, 2.0]])
    MOIT.geomean3ftest(mock, config)
end
@testset "Exponential" begin
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1., 2., 2exp(1/2)],
                          (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})  => [1 + exp(1/2), 1 + exp(1/2)/2])
    MOIT.exp1vtest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1., 2., 2exp(1/2)],
                          (MOI.VectorAffineFunction{Float64}, MOI.ExponentialCone)   => [[-exp(1/2), -exp(1/2)/2, 1.]],
                          (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})  => [1 + exp(1/2), 1 + exp(1/2)/2])
    MOIT.exp1ftest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [0., -0.3, 0., exp(-0.3), exp(-0.3), exp(-0.3), 0., 1.0, 0.],
                          (MOI.VectorAffineFunction{Float64}, MOI.ExponentialCone)   => [[-exp(-0.3)/2, -1.3exp(-0.3)/2, 0.5], [-exp(-0.3)/2, -1.3exp(-0.3)/2, 0.5]],
                          (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})  => [-1.0, exp(-0.3)*0.3],
                          (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [-exp(-0.3)*0.3],
                          (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives)      => [[0.0, exp(-0.3), exp(-0.3)/2], [0.0, 0.0, exp(-0.3)/2]])
    MOIT.exp2test(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [log(5), 5.],
                          (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [0.],
                          (MOI.VectorAffineFunction{Float64}, MOI.ExponentialCone)   => [[-1., log(5)-1, 1/5]])
    MOIT.exp3test(mock, config)
end
@testset "Dual Exponential" begin
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [-exp(1 / 2), -exp(1 / 2) / 2, 1.0, 1 + exp(1 / 2), 1 + exp(1 / 2) / 2],
                          (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})  => [-1.0, -2.0, -2exp(1 / 2)])
    MOIT.dualexp1vtest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [-exp(1 / 2), -exp(1 / 2) / 2, 1.0, 1 + exp(1 / 2), 1 + exp(1 / 2) / 2],
                          (MOI.VectorAffineFunction{Float64}, MOI.DualExponentialCone)   => [[1., 2.0, 2exp(1 / 2)]],
                          (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})  => [-1.0, -2.0, -2exp(1 / 2)])
    MOIT.dualexp1ftest(mock, config)
end
@testset "Dual Power" begin
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [0.9, 0.1, 1.0, -0.9, -0.1],
                          (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})  => [-1.0, -1.0, 1.0])
    MOIT.dualpow1vtest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [0.9, 0.1, 1.0, -0.9, -0.1],
                          (MOI.VectorAffineFunction{Float64}, MOI.DualPowerCone{Float64})   => [[1.0, 1.0, -1.0]],
                          (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})  => [-1.0, -1.0, 1.0])
    MOIT.dualpow1ftest(mock, config)
end
@testset "PSD" begin
    # PSD0
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, ones(3),
                          (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})                 => [2])
    MOIT.psdt0vtest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, ones(3),
                          (MOI.VectorAffineFunction{Float64}, MOI.PositiveSemidefiniteConeTriangle) => [[1.0, -1.0, 1.0]],
                          (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})                 => [2])
    MOIT.psdt0ftest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, ones(4),
                          (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})                 => [2])
    MOIT.psds0vtest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, ones(4),
                          (MOI.VectorAffineFunction{Float64}, MOI.PositiveSemidefiniteConeSquare)   => [[1.0, -2.0, 0.0, 1.0]],
                          (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})                 => [2])
    MOIT.psds0ftest(mock, config)
    # PSD1
    δ = √(1 + (3*√2+2)*√(-116*√2+166) / 14) / 2
    ε = √((1 - 2*(√2-1)*δ^2) / (2-√2))
    y2 = 1 - ε*δ
    y1 = 1 - √2*y2
    obj = y1 + y2/2
    k = -2*δ/ε
    x2 = ((3-2obj)*(2+k^2)-4) / (4*(2+k^2)-4*√2)
    α = √(3-2obj-4x2)/2
    β = k*α
    Xv = [α^2, α*β, β^2, α^2, α*β, α^2]
    xv = [√2*x2, x2, x2]
    cX0 = 1+(√2-1)*y2
    cX1 = 1-y2
    cX2 = -y2
    cXv = [cX0, cX1, cX0, cX2, cX1, cX0]
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [Xv; xv],
                          (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})                 => [y1, y2])
    MOIT.psdt1vtest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [Xv; xv],
                          (MOI.VectorAffineFunction{Float64}, MOI.PositiveSemidefiniteConeTriangle) => [cXv],
                          (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})                 => [y1, y2])
    MOIT.psdt1ftest(mock, config)
    Xv = [α^2, α*β, α^2, α*β, β^2, α*β, α^2, α*β, α^2]
    cXv = [cX0, cX1, cX2, cX1, cX0, cX1, cX2, cX1, cX0]
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [Xv; xv],
                          (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})                 => [y1, y2])
    MOIT.psds1vtest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [Xv; xv],
                          (MOI.VectorAffineFunction{Float64}, MOI.PositiveSemidefiniteConeSquare)   => [cXv],
                          (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})                 => [y1, y2])
    MOIT.psds1ftest(mock, config)
    # PSD2
    η = 10.0
    α = 0.8
    δ = 0.9
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [2η/3, 0, η/3, 0, 0, 0, η*δ*(1 - 1/√3)/2],
                          (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64})             => [δ*(1-1/√3)/2],
                          (MOI.VectorAffineFunction{Float64}, MOI.Nonpositives)                     => [[0, -α/√3+δ/(2*√6)*(2*√2-1), 0, -3δ*(1-1/√3)/8, -3δ*(1-1/√3)/8, -δ*(3 - 2*√3 + 1/√3)/8]],
                          (MOI.VectorAffineFunction{Float64}, MOI.PositiveSemidefiniteConeTriangle) => [[(1-1/√3)/2, 1/√6, (1+1/√3)/2]],
                          (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})                 => [0])
    MOIT.psdt2test(mock, config)
    # PSD3
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, ones(1),
        (MOI.VectorAffineFunction{Float64}, MOI.PositiveSemidefiniteConeTriangle) => [[2, -1, 2, -1, -1, 2] / 6])
    MOIT.psdt3test(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, ones(1),
        (MOI.VectorAffineFunction{Float64}, MOI.PositiveSemidefiniteConeSquare) => [[1, 0, 0, -1, 1, 0, -1, -1, 1] / 3])
    MOIT.psds3test(mock, config)
end
@testset "LogDet" begin
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [0, 1, 0, 1, 1],
        (MOI.SingleVariable, MOI.EqualTo{Float64}) => [2],
        (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) => [[1, 1]],
        (MOI.VectorOfVariables, MOI.LogDetConeTriangle) => [[-1, -2, 1, 0, 1]])
    mock.eval_variable_constraint_dual = false
    MOIT.logdett1vtest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [0, 1, 0, 1, 1],
        (MOI.SingleVariable, MOI.EqualTo{Float64}) => [2],
        (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) => [[1, 1]],
        (MOI.VectorAffineFunction{Float64}, MOI.LogDetConeTriangle) => [[-1, -2, 1, 0, 1]])
    MOIT.logdett1ftest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [log(5)],
        (MOI.VectorAffineFunction{Float64}, MOI.LogDetConeTriangle) => [[-1, log(5) - 3, 1, -1, 1.6, 0, -0.2, 0.4]])
    MOIT.logdett2test(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [0, 1, 0, 0, 1, 1],
        (MOI.SingleVariable, MOI.EqualTo{Float64}) => [2],
        (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) => [[1, 1]],
        (MOI.VectorOfVariables, MOI.LogDetConeSquare) => [[-1, -2, 1, 0, 0, 1]])
    mock.eval_variable_constraint_dual = false
    MOIT.logdets1vtest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [0, 1, 0, 0, 1, 1],
        (MOI.SingleVariable, MOI.EqualTo{Float64}) => [2],
        (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) => [[1, 1]],
        (MOI.VectorAffineFunction{Float64}, MOI.LogDetConeSquare) => [[-1, -2, 1, 0, 0, 1]])
    MOIT.logdets1ftest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [log(5)],
        (MOI.VectorAffineFunction{Float64}, MOI.LogDetConeSquare) => [[-1, log(5) - 3, 1, -1, 0, -1, 1.6, -0.2, 0, -0.2, 0.4]])
    MOIT.logdets2test(mock, config)
end
@testset "RootDet" begin
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1, 1, 0, 1],
        (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) => [[0.5, 0.5]],
        (MOI.VectorOfVariables, MOI.RootDetConeTriangle) => [[-1.0, 0.5, 0.0, 0.5]])
    MOIT.rootdett1vtest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1, 1, 0, 1],
        (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) => [[0.5, 0.5]],
        (MOI.VectorAffineFunction{Float64}, MOI.RootDetConeTriangle) => [[-1.0, 0.5, 0.0, 0.5]])
    MOIT.rootdett1ftest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [5 ^ inv(3)],
        (MOI.VectorAffineFunction{Float64}, MOI.RootDetConeTriangle) => [vcat(-1, [1, -1, 1.6, 0, -0.2, 0.4] / 3 * (5 ^ inv(3)))])
    MOIT.rootdett2test(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1, 1, 0, 0, 1],
        (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) => [[0.5, 0.5]],
        (MOI.VectorOfVariables, MOI.RootDetConeSquare) => [[-1.0, 0.5, 0.0, 0.0, 0.5]])
    MOIT.rootdets1vtest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1, 1, 0, 0, 1],
        (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) => [[0.5, 0.5]],
        (MOI.VectorAffineFunction{Float64}, MOI.RootDetConeSquare) => [[-1.0, 0.5, 0.0, 0.0, 0.5]])
    MOIT.rootdets1ftest(mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [5 ^ inv(3)],
        (MOI.VectorAffineFunction{Float64}, MOI.RootDetConeSquare) => [vcat(-1, [1, -1, 0, -1, 1.6, -0.2, 0, -0.2, 0.4] / 3 * (5 ^ inv(3)))])
    MOIT.rootdets2test(mock, config)
end


text:
#imports
using NearestNeighbors, Statistics, LinearAlgebra, Random, Flux#master
using Flux: onehotbatch, onecold, onehot, crossentropy, throttle, NNlib, @functor
using Zygote: @nograd, Params
using Statistics: mean
using Base.Iterators: partition

#args
include(joinpath(@__DIR__,"config.jl"))

# possible datset {modelnet,3dmnist}
if "$(ARGS[1])" in ["modelnet", "3dmnist"]
    dataset = "$(ARGS[1])"
    root = "./$(dataset)/"
elseif "$(ARGS[1])" == ""
    print("using dataset from the config $(dataset)")
else
   print("$(ARGS[1]) dataset not supported, currently supported dataset are modelnet and 3dmnist")
   exit()
end

if dataset == "3dmnist"
    num_classes=10
end

#data and model
include(joinpath(@__DIR__,"data.jl"))
include(joinpath(@__DIR__,"model.jl"))

# input: width*height*channel*minibatch

# Fetching the train and validation data and getting them into proper shape
print("Preparing Dataset\n")
train, (valX, valY) = get_data(num_classes)
print("Dataset loaded\n")

# Defining the loss and accuracy functions

m = DGCNN_cls(num_classes, K, npoints)

loss(x, y) = crossentropy(m(x), y)

accuracy(x, y) = mean(onecold(cpu(m(x)), 1:num_classes) .== onecold(cpu(y), 1:num_classes))

# Defining the callback and the optimizer

opt = ADAM(0.003)

# Starting to train models

function custom_train!(loss, ps, data, opt, epochs)
    ps = Zygote.Params(ps)
    for epoch in 1:epochs
        running_loss = 0
        for d in data
        gs = gradient(ps) do
            training_loss = loss(d...)
            running_loss += training_loss
            return training_loss
        end
        Flux.update!(opt, ps, gs)
        end
        print("Epoch: $(epoch), epoch_loss: $(running_loss), accuracy: $(accuracy(valX, valY))\n")
    end
  end

print("Training Start\n")
custom_train!(loss, params(m), train, opt, epochs)
print("Training Finished\n")

# Accuracy of valset
@show(accuracy(valX, valY))

max_stars_repo_name: nirmal-suthar/DGCNN

max_issues_repo_name: nirmal-suthar/DGCNN

max_forks_repo_name: nirmal-suthar/DGCNN


text:
##### Generic methods #####

"""
    length(d::MultivariateDistribution) -> Int

Return the sample dimension of distribution `d`.
"""
length(d::MultivariateDistribution)

"""
    size(d::MultivariateDistribution)

Return the sample size of distribution `d`, *i.e* `(length(d),)`.
"""
size(d::MultivariateDistribution)

## sampling

# multiple multivariate, must allocate matrix
# TODO: inconsistency with other `ArrayLikeVariate`s and `rand(s, (n,))` - maybe remove?
rand(rng::AbstractRNG, s::Sampleable{Multivariate}, n::Int) =
    @inbounds rand!(rng, sampler(s), Matrix{eltype(s)}(undef, length(s), n))
rand(rng::AbstractRNG, s::Sampleable{Multivariate,Continuous}, n::Int) =
    @inbounds rand!(rng, sampler(s), Matrix{float(eltype(s))}(undef, length(s), n))

## domain

"""
    insupport(d::MultivariateDistribution, x::AbstractArray)

If ``x`` is a vector, it returns whether x is within the support of ``d``.
If ``x`` is a matrix, it returns whether every column in ``x`` is within the support of ``d``.
"""
insupport{D<:MultivariateDistribution}(d::Union{D, Type{D}}, x::AbstractArray)

function insupport!(r::AbstractArray, d::Union{D,Type{D}}, X::AbstractMatrix) where D<:MultivariateDistribution
    n = length(r)
    size(X) == (length(d),n) ||
        throw(DimensionMismatch("Inconsistent array dimensions."))
    for i in 1:n
        @inbounds r[i] = insupport(d, view(X, :, i))
    end
    return r
end

insupport(d::Union{D,Type{D}}, X::AbstractMatrix) where {D<:MultivariateDistribution} =
    insupport!(BitArray(undef, size(X,2)), d, X)

## statistics

"""
    mean(d::MultivariateDistribution)

Compute the mean vector of distribution `d`.
"""
mean(d::MultivariateDistribution)

"""
    var(d::MultivariateDistribution)

Compute the vector of element-wise variances for distribution `d`.
"""
var(d::MultivariateDistribution)

"""
    entropy(d::MultivariateDistribution)

Compute the entropy value of distribution `d`.
"""
entropy(d::MultivariateDistribution)

"""
    entropy(d::MultivariateDistribution, b::Real)

Compute the entropy value of distribution ``d``, w.r.t. a given base.
"""
entropy(d::MultivariateDistribution, b::Real) = entropy(d) / log(b)

"""
    cov(d::MultivariateDistribution)

Compute the covariance matrix for distribution `d`. (`cor` is provided based on `cov`).
"""
cov(d::MultivariateDistribution)

"""
    cor(d::MultivariateDistribution)

Computes the correlation matrix for distribution `d`.
"""
function cor(d::MultivariateDistribution)
    C = cov(d)
    n = size(C, 1)
    @assert size(C, 2) == n
    R = Matrix{eltype(C)}(undef, n, n)

    for j = 1:n
        for i = 1:j-1
            @inbounds R[i, j] = R[j, i]
        end
        R[j, j] = 1.0
        for i = j+1:n
            @inbounds R[i, j] = C[i, j] / sqrt(C[i, i] * C[j, j])
        end
    end

    return R
end

##### Specific distributions #####

for fname in ["dirichlet.jl",
              "multinomial.jl",
              "dirichletmultinomial.jl",
              "mvnormal.jl",
              "mvnormalcanon.jl",
              "mvlognormal.jl",
              "mvtdist.jl",
              "product.jl",
              "vonmisesfisher.jl"]
    include(joinpath("multivariate", fname))
end


text:
import Term.Colors:
    is_named_color,
    is_rgb_color,
    is_hex_color,
    is_color,
    is_background,
    get_color,
    NamedColor,
    RGBColor,
    BitColor

@testset "\e[34mANSI detect color" begin
    for name in ["red", "blue", "black", "grey42", "deep_sky_blue3", "dark_blue"]
        @test is_named_color(name)
        @test is_color(name)
        @test !is_rgb_color(name)
        @test !is_hex_color(name)
    end

    for name in ["255", "1", "21", "6", "128"]
        @test is_named_color(name)
        @test is_color(name)
        @test !is_rgb_color(name)
        @test !is_hex_color(name)
    end

    rgbs = ["(1, .8, .0)", "(.2,.1,.33)", "(.5,0.1, .99)", "(255, 128,128)", "(128,58,200)"]
    for rgb in rgbs
        @test is_rgb_color(rgb)
        @test is_color(rgb)
        @test !is_hex_color(rgb)
        @test !is_named_color(rgb)
    end

    hexes = ["#ffffff", "#000000", "#dadada", "#123123"]
    for hex in hexes
        @test is_hex_color(hex)
        @test is_color(hex)
        @test !is_named_color(hex)
        @test !is_rgb_color(hex)
    end
end

@testset "\e[34mANSI get color" begin
    for name in ["red", "blue", "black"]
        @test get_color(name) isa NamedColor
    end

    for name in ["255", "1", "21", "6", "128", "grey42", "deep_sky_blue3", "dark_blue"]
        @test get_color(name) isa BitColor
    end

    for rgb in
        ["(1, .8, .0)", "(.2,.1,.33)", "(.5,0.1, .99)", "(255, 128,128)", "(128,58,200)"]
        @test get_color(rgb) isa RGBColor
    end

    hexes = ["#ffffff", "#000000", "#dada1a", "#123123"]
    for hex in hexes
        @test get_color(hex) isa RGBColor
    end
end

@testset "\e[34mANSI background color" begin
    for name in ["on_red", "on_blue", "on_black"]
        @test is_background(name)
        @test !is_color(name)
        @test get_color(name; bg = true) isa NamedColor
    end

    for name in ["on_255", "on_1", "on_21", "on_6", "on_128"]
        @test is_background(name)
        @test !is_color(name)
        @test get_color(name; bg = true) isa BitColor
    end

    for rgb in [
        "on_(1, .8, .0)",
        "on_(.2,.1,.33)",
        "on_(.5,0.1, .99)",
        "on_(255, 128,128)",
        "on_(128,58,200)",
    ]
        @test is_background(rgb)
        @test !is_color(rgb)
        @test get_color(rgb; bg = true) isa RGBColor
    end

    hexes = ["on_#ffffff", "on_#000000", "on_#dada2a", "on_#123123"]
    for hex in hexes
        @test is_background(hex)
        @test get_color(hex; bg = true) isa RGBColor
    end
end