text:
## GAPM - UC
## 1. Main GAPM Loop
# Available methods and stats
# GAPM: Adaptative dual-mean comparison
# GAPM_Vec: Full vector comparison
# AffProp: Affinity Propagation
Methods = ["GAPM" "GAPM_Vec" "AffProp"]
# Similarity distances for Affinity Propagation
dists = getdists()
Dict_AffProp = Dict()
ALL_DUALS = false

Meth, mem = "GAPM_Vec", 0
    Meth == "AffProp" ? nN = length(dists) : nN=1
for mem in ["0"]#["0", "1", "Inf"]
    for n in 1:nN
        global Dict_AffProp, RES
        d = dists[n]
        n == 1 ? println("Here we go..") : nothing
        # try
        stat = "L1"
        L = 3

        stepprint = true#(Meth != "AffProp")
        # Parameters for GAPM
        α, γ = 1/2, 1/2
        RES = SolarData[:,1:100]
            N = size(RES,2)
            S = 1:N
            α_0 = α
            μ = mean(RES, dims=2) # Vector with mean RES values
            k, K = 1, 20
            NP = 1
            Ω = 1:NP
            ρ = [1]
            Ind = ones(1,N)'
            R = mean(RES, dims=2)

            Gap = 100
            UB = Inf
            LB = -Inf
            OldGAP = Inf
            GapVec = Array{Float64}[]
            TimeVec = Array{Float64}[]
            Bounds = Array{Float64,2}[]

        # Meth == "AffProp" ? println("Starting UC - "*Meth*": "*stat) : println("Starting UC-"*Meth)
        Meth == "AffProp" ? nothing : println("Starting UC-"*Meth)

            Meth == "GAPM" ? println("α: $(α) - γ: $(γ)") : nothing

        start = now()
            timeT = @elapsed while (k ≤ K)&&(Gap≥0.1)
            global Ind_old, NP_old, Q, duals#, k, ρ, Ω, NP, Ind, Gap, LB, UB, α
            global inds, m#, OldGAP, GapVec, TimeVec, Bounds

            OldGAP = Gap

                if k > 1
                    Ind_old = Ind
                    duals_old = duals
                    if Meth=="GAPM"
                        α = α*γ
                        Ind = GAPM(Ind_old, duals, α_val = α)
                    elseif Meth=="GAPM_Vec"
                        Ind = GAPM_Vec(Ind_old, duals)
                        # Ind2 = Ind
                    elseif Meth=="AffProp"
                        if mem == "Inf"
                            Ind = memoryPart(Ind_old, duals, d, memory = mem)
                        else
                            Ind = memoryPart(duals_old, duals, d, memory = mem)
                            # println()
                            println("$Ind")
                        end

                    end
                    NP = Integer(maximum(Ind))
                    Ω = 1:NP # Scenarios' set
                    R = zeros(length(T),NP)
                    ρ = zeros(NP) # Partition probability
                    R_extreme = zeros(length(T), 2, NP)# Extreme points of partitions
                    for n in 1:NP
                        Ns = sum(Ind .== n)
                        inds = S[Ind[:].==n]
                        R[:,n] = mean(RES[:,inds], dims=2) # mean RES generation in partition
                        ρ[n] = Ns/N
                    end
                    # println("k: $k")
                end
                # println("R: $R")
                V, Pexp, Rup, Rdw, Obj1, LB1, m = UC(Ω, ρ, R)
                LB = max(LB, LB1)

                duals, Q = UCEvaluation(V, Pexp, Rup, Rdw, RES)

                    UB = min(UB, (Obj1 + Q))
                    Gap = (100*(UB-LB)/LB)

                    # Compare the whole vector or its mean
                    Meth == "GAPM" ? duals = mean(duals,dims = 2) : nothing

                t = now()
                    spentTime = 3600*(Dates.hour(t)-Dates.hour(start))+60*(Dates.minute(t)-Dates.minute(start))+(Dates.second(t)-Dates.second(start))+(Dates.millisecond(t)-Dates.millisecond(start))/1000

                    if stepprint == true
                        println("$k & $NP & $(round(LB,digits=3)) & $(round(UB,digits=3)) & $(round(Gap,digits=3)) & $(spentTime) \\\\")
                    end

                GapVec = cat(dims=1, GapVec, [Gap])
                Bounds = cat(dims=1,Bounds, [LB UB])
                TimeVec = cat(dims=1, TimeVec, [spentTime])

            # Controlling Gap stop
            if k>1
                OldGAP ≤ Gap ? break : k+=1
                else
                    Gap = Inf
                    k+=1
                end
            end

            if stepprint == true
                Meth == "GAPM" ? println("α_end: $(α*γ)") : nothing
                Meth == "AffProp" ? println(Meth*": "*d) : println(Meth)
                Meth == "GAPM" ? println("α: $(α_0) - γ: $(γ)") : nothing
            end
        if Meth == "AffProp"
            println()
            println("\\multicolumn{2}{c|}{$(replace(d, "_"=>" "))} & $NP & $(round(Gap,digits=3)) & $(k) & $(round(timeT,digits=2)) \\\\ %$n")
            elseif Meth == "GAPM"
                println()
                println("$α_0 & $γ & $NP & $(round(Gap,digits=3)) & $(k) & $(round(timeT,digits=2)) \\\\")
            elseif Meth == "GAPM_Vec"
                println()
                println("\\multicolumn{6}{c}{GAPM--FullVec} & $NP & $(round(Gap,digits=3)) & $(k) & $(round(timeT,digits=2)) \\\\")
            end

        # Saving in Dictionary
        if Meth == "AffProp"
            Dict_AffProp[d*"_mem$(mem)_Gap"], Dict_AffProp[d*"_mem$(mem)_NP"], Dict_AffProp[d*"_mem$(mem)_TimeVec"] = GapVec, NP, TimeVec
            elseif Meth == "GAPM_Vec"
                Dict_AffProp["GAPMVec_Gap"], Dict_AffProp["GAPMVec_NP"], Dict_AffProp["GAPMVec_TimeVec"] = GapVec, NP, TimeVec
            elseif Meth == "GAPM"
                AlphaGamma = "α$(round(α_0, digits=3))_γ$(round(γ, digits=3))"
                AlphaGamma = replace.((AlphaGamma),"."=>"")
                Dict_AffProp["GAPM_"*AlphaGamma*"_Gap"], Dict_AffProp["GAPM_"*AlphaGamma*"_NP"], Dict_AffProp["GAPM_"*AlphaGamma*"_TimeVec"] = GapVec, NP, TimeVec
            end

        # catch; println("\\multicolumn{2}{c|}{$(replace(d, "_"=>" "))} & -- & -- & -- & -- \\\\ % $d did not work.")
        # end
    end
end

# Plotting Results Gap vs Time
DictAffProp = replace.(keys(Dict_AffProp),"_NP"=>"")
    DictAffProp = replace.((DictAffProp),"_Gap"=>"")
    DictAffProp = replace.((DictAffProp),"_TimeVec"=>"")
    DictAffProp = replace.((DictAffProp),"_mem0"=>"")
    DictAffProp = replace.((DictAffProp),"_mem1"=>"")
    DictAffProp = replace.((DictAffProp),"_memInf"=>"")
    DictAffProp = replace.((DictAffProp),"mem0"=>"")
    unique!(DictAffProp)
    sort!(DictAffProp)
    # show(DictAffProp)

    # collect(DictAffProp[i] for i in 1:length(DictAffProp) if DictAffProp[i][1:4] == "GAPM")
    plotAffProp(Dict_AffProp, chosen = DictAffProp, ymin= 0.0005)

plotAffProp(Dict_AffProp, chosen = DictAffProp, ymin= 0.0005)

DictAffProp[1][1:4]

# Full UC
N = 10
    @elapsed V, Pexp, Rup, Rdw, Obj1, LB = UC(1:N, ones(N)./N, SolarData[:,1:N],printing = 1)
LB

V, Pexp, Rup, Rdw, Obj1, LB = UC(1:1, [1], mean(SolarData,dims=2)[:,:,1], printing=0)
    println("LB: $LB")
    println("LB: $LB")
    # plot!(Pexp')

    plot(sum(Pd,dims=1)', label = "Load", alpha= 0.5)
    plot!(1:24,sum(Pexp,dims=1)'+mean(SolarData,dims=2), alpha = 0.5, labels = "Gen")
heatmap(V')
plot(Pexp')
plot(Rup', label="")
plot(Rdw', label="")

UCEvaluation(V, Pexp, Rup, Rdw, RES)[2]
Obj1

plot!(maximum(SolarData,dims=2), label = "RES")
    hline!([sum(GenData.Pmax)], label = "GenCap")

max_stars_repo_name: gonzalez-alvaro/MultidimensionalGAPM

max_issues_repo_name: gonzalez-alvaro/MultidimensionalGAPM

max_forks_repo_name: gonzalez-alvaro/MultidimensionalGAPM


text:
using Ito
using Ito.Integration

#All the following use the default convergence limit of 1e-9, and max 50 iterations

epsilon = 1e-9

r=integrate((x)->1, AdaptiveSimpsonsIntegration(), 0.0,1)
@assert abs(1 - r ) < epsilon

r=integrate((x)->x, AdaptiveSimpsonsIntegration(), 0.0,1)
@assert abs(1/2 - r ) < epsilon

r=integrate((x)->x*x, AdaptiveSimpsonsIntegration(), 0.0,1)
@assert abs(1/3 - r ) < epsilon

r=integrate(sin, AdaptiveSimpsonsIntegration(), 0,pi)
@assert abs(2-r) < epsilon

r=integrate(cos, AdaptiveSimpsonsIntegration(), 0,pi)
@assert abs(0-r) < epsilon

r=integrate((x)->(sin(x))^2+(sin(x))^2, AdaptiveSimpsonsIntegration(epsilon, 30), 0,2pi)
@assert abs(2pi-r) < epsilon

r=integrate((x)->Distributions.pdf(Distributions.Normal(),x), AdaptiveSimpsonsIntegration(), -10,10)
@assert abs(1-r) < epsilon

epsilon = 1e-6
r=integrate((x)->1/x, AdaptiveSimpsonsIntegration(), 0.01,1)
@assert abs(4.60517 - r ) < epsilon

r=integrate((x)->(cos(x))^8, AdaptiveSimpsonsIntegration(epsilon, 30), 0,2pi)
@assert abs(35pi/64 - r ) < epsilon

r=integrate((x)->sin(x) - sin(x^2) + sin(x^3), AdaptiveSimpsonsIntegration(), pi,2pi)
@assert abs(-1.830467 - r) < epsilon

ext: jl


text:
#=

Implement a job scheduler which takes in a function f and an integer n, and calls f after n milliseconds.

=#

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

@test round(Int, @elapsed job_scheduler(job, 10000)) == 10

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
abstract type «name» end
abstract type «name» <: «supertype» end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function test_raw_blas(T)
  println("-- Testing RawBLAS{$T} Utilities")

  eps = 1e-10
  M, N, K = (108, 6, 36)
  A = rand(T, M, K)
  B = rand(T, N, K)

  C = rand(T, M, N)
  C2 = rand(T, M, N)
  RawBLAS.gemm!('N', 'T', M, N, K, convert(T,1.0), A, M, B, N, convert(T,0.0), C, M)
  BLAS.gemm!('N', 'T', convert(T,1.0), A, B, convert(T,0.0), C2)

  @test all(abs(C - C2) .< eps)

  A2 = rand(T, M, 2K)
  B2 = rand(T, N, 3K)
  C = rand(T, M, N)
  C2 = A2[:, 1:K] * B2[:, 1:K]'
  RawBLAS.gemm!('N', 'T', M, N, K, convert(T,1.0), A2, M, B2, N, convert(T,0.0), C, M)
  @test all(abs(C - C2) .< eps)

  C2 = A2[:, K+1:2K] * B2[:, 2K+1:3K]'
  C = rand(T, M, 2N)
  RawBLAS.gemm!('N', 'T', M, N, K, convert(T, 1.0),
      pointer(A2) + M*K*sizeof(T), M,
      pointer(B2) + 2N*K*sizeof(T), N,
      convert(T,0.0), pointer(C) + M*N*sizeof(T), M)
  @test all(abs(C2 - C[1:M,N+1:2N]) .< eps)
end

test_raw_blas() = begin
  test_raw_blas(Float32)
  test_raw_blas(Float64)
end

test_raw_blas()


text:
function plot(plon,plat,data; cmap = "RdYlBu_r", coastlinecolor="magenta")
    pcolormesh(plon,plat,data,cmap=cmap)
    lon,lat,lsmask = GeoDatasets.landseamask(;resolution='l',grid=5)
    contour(lon,lat,lsmask',[0.5],linewidths=[1.],colors=coastlinecolor);
    xlim(extrema(plon))
    ylim(extrema(plat))
    return nothing
end

"""
    makeplots(wavname,satellite_name)

Decodes the APT signal in `wavname` as recorded by gqrx using
wide FM-mono demodulation.
The file name `wavname` should  have the followng structure:
`string_date_time_frequency.wav` like `gqrx_20190811_075102_137620000.wav`.
Date and time of the file name are in UTC (not local time). `satellite_name` is
the name of the satellite (generally `"NOAA 15"`, `"NOAA 18"`, `"NOAA 19"`).

# Example:

```julia
import APTDecoder

wavname = "gqrx_20190825_182745_137620000.wav"
APTDecoder.makeplots(wavname,"NOAA 15")
```

"""
function makeplots(wavname,satellite_name;
                   starttime = nothing,
                   eop = nothing,
                   prefix = replace(wavname,r".wav$" => ""),
                   qrange = (0.01,0.99),
                   coastlinecolor = "magenta",
                   cmap = "RdYlBu_r",
                   tles = get_tle(:weather),
                   dpi = 150)

    if starttime == nothing
        starttime = APTDecoder.starttimename(basename(wavname))
    end

    y,Fs,nbits,opt = FileIO.load(wavname)

    datatime,(channelA,channelB),data = APTDecoder.decode(y,Fs)

    vmin,vmax = quantile(view(data,:),[qrange[1],qrange[2]])
    data[data .> vmax] .= vmax;
    data[data .< vmin] .= vmin;

    # save raw image
    imagenames = (
        rawname = prefix * "_raw.png",
        channel_a = prefix * "_channel_a.png",
        channel_b = prefix * "_channel_b.png")

    FileIO.save(imagenames.rawname, colorview(Gray, data[:,1:3:end]./maximum(data)))
    grid_color = [0,0.7,0.6]

    fig = figure()
    plt.style.use("dark_background")
    Alon,Alat,Adata = APTDecoder.georeference(
        channelA,satellite_name,datatime,starttime, eop = eop, tles = tles)
    APTDecoder.plot(Alon,Alat,Adata; coastlinecolor=coastlinecolor, cmap=cmap)
    plt.grid(linestyle = "--",color=grid_color)
    savefig(imagenames.channel_a,dpi=dpi,pad_inches=0, bbox_inches="tight", transparent=false)
    fig.clf()

    Blon,Blat,Bdata = APTDecoder.georeference(
        channelB,satellite_name,datatime,starttime, eop = eop, tles = tles)
    APTDecoder.plot(Blon,Blat,Bdata; coastlinecolor=coastlinecolor, cmap=cmap)
    plt.grid(linestyle = "--",color=grid_color)
    savefig(imagenames.channel_b,dpi=dpi,pad_inches=0, bbox_inches="tight", transparent=false)
    close(fig)

    return imagenames
end

max_stars_repo_name: Alexander-Barth/APTDecoder

max_issues_repo_name: Alexander-Barth/APTDecoder

max_forks_repo_name: Alexander-Barth/APTDecoder


text:
@testset "formula" begin

    using StatsModels: hasresponse, hasintercept, omitsintercept

    y, x1, x2, x3, a, b, c, onet = term.((:y, :x1, :x2, :x3, :a, :b, :c, 1))

    ## totally empty
    t = @formula(0 ~ 0)
    @test !hasresponse(t)
    @test !hasintercept(t)
    @test omitsintercept(t)
    @test t.rhs == ConstantTerm(0)
    @test issetequal(terms(t), [ConstantTerm(0)])

    ## empty lhs, intercept on rhs
    t = @formula(0 ~ 1)
    @test !hasresponse(t)
    @test hasintercept(t)
    @test !omitsintercept(t)

    ## empty RHS
    t = @formula(y ~ 0)
    @test hasintercept(t) == false
    @test omitsintercept(t) == true
    @test hasresponse(t)
    @test t.rhs == ConstantTerm(0)
    @test issetequal(terms(t), term.((:y, 0)))

    t = @formula(y ~ -1)
    @test hasintercept(t) == false
    @test omitsintercept(t) == true

    ## intercept-only
    t = @formula(y ~ 1)
    @test hasresponse(t) == true
    @test hasintercept(t) == true
    @test t.rhs == onet
    @test issetequal(terms(t), (onet, y))

    ## terms add
    t = @formula(y ~ 1 + x1 + x2)
    @test hasintercept(t) == true
    @test t.rhs == (onet, x1, x2)
    @test issetequal(terms(t), [y, onet, x1, x2])

    ## implicit intercept behavior: NO intercept after @formula
    t = @formula(y ~ x1 + x2)
    @test hasintercept(t) == false
    @test omitsintercept(t) == false
    @test t.rhs == (x1, x2)
    @test issetequal(terms(t), [y, x1, x2])

    ## no intercept
    t = @formula(y ~ 0 + x1 + x2)
    @test hasintercept(t) == false
    @test omitsintercept(t) == true
    @test t.rhs == term.((0, :x1, :x2))

    t = @formula(y ~ -1 + x1 + x2)
    @test hasintercept(t) == false
    @test omitsintercept(t) == true
    @test t.rhs == term.((-1, :x1, :x2))

    t = @formula(y ~ x1 & x2)
    @test hasintercept(t) == false
    @test omitsintercept(t) == false
    @test t.rhs == x1&x2
    @test issetequal(terms(t), [y, x1, x2])

    ## `*` expansion
    t = @formula(y ~ x1 * x2)
    @test hasintercept(t) == false
    @test omitsintercept(t) == false
    @test t.rhs == (x1, x2, x1&x2)
    @test issetequal(terms(t), [y, x1, x2])

    ## associative rule:
    ## +
    t = @formula(y ~ x1 + x2 + x3)
    @test t.rhs == (x1, x2, x3)

    ## &
    t = @formula(y ~ x1 & x2 & x3)
    @test t.rhs == x1&x2&x3
    @test issetequal(terms(t), [y, x1, x2, x3])

    ## distributive property of + and &
    t = @formula(y ~ x1 & (x2 + x3))
    @test t.rhs == (x1&x2, x1&x3)
    @test issetequal(terms(t), [y, x1, x2, x3])

    ## ordering of interaction terms is preserved across distributive
    t = @formula(y ~ (x2 + x3) & x1)
    @test t.rhs == x2&x1 + x3&x1

    ## distributive with *
    t = @formula(y ~ (a + b) * c)
    @test t.rhs == (a, b, c, a&c, b&c)

    ## three-way *
    t = @formula(y ~ a * b * c)
    @test t.rhs == (a, b, c, a&b, a&c, b&c, a&b&c)
    @test issetequal(terms(t), (y, a, b, c))

    ## Interactions with `1` reduce to main effect.
    t = @formula(y ~ 1 & x1)
    @test t.rhs == x1

    t = @formula(y ~ (1 + x1) & x2)
    @test t.rhs == (x2, x1&x2)

    ## PR #54 breaks formula-level equality because original (un-lowered)
    ## expression is kept on Formula struct.  but functional (RHS) equality
    ## should be maintained
    using StatsModels: drop_term

    @test drop_term(@formula(foo ~ 1 + bar + baz), term(:bar)) ==
        @formula(foo ~ 1 + baz)
    @test drop_term(@formula(foo ~ 1 + bar + baz), term(1)) ==
        @formula(foo ~ bar + baz)

    # drop_term no longer checks for whether term is found...
    @test_broken drop_term(@formula(foo ~ bar + baz), term(0))
    @test_broken drop_term(@formula(foo ~ bar + baz), term(:boz))

    form = @formula(foo ~ 1 + bar + baz)
    @test form == @formula(foo ~ 1 + bar + baz)
    @test drop_term(form, term(:bar)) == @formula(foo ~ 1 + baz)
    # drop_term creates a new formula:
    @test form != @formula(foo ~ 1 + baz)

    # Incorrect formula separator
    @test_throws LoadError @eval @formula(y => x + 1)

    # copying formulas
    f = @formula(foo ~ 1 + bar)
    @test f == deepcopy(f)

    f = @formula(foo ~ bar)
    @test f == deepcopy(f)

end


text:
# Types -------

abstract Authorization

immutable BasicAuth <: Authorization
    user::String
    password::String
end

function show(io::IO, a::BasicAuth)
    pw_str = repeat("*", 8)
    print(io, "GitHub Authorization ($(a.user), $pw_str))")
end

immutable OAuth2 <: Authorization
    token::String
end

function show(io::IO, a::OAuth2)
    token_str = a.token[1:6] * repeat("*", length(a.token) - 6)
    print(io, "GitHub Authorization ($token_str)")
end

immutable AnonymousAuth <: Authorization
end

# Interface -------

function authenticate(user::String, password::String)
    auth = BasicAuth(user, password)
end

function authenticate(token::String)
    auth = OAuth2(token)

    r = get(API_ENDPOINT; query = { "access_token" => auth.token })
    if r.status < 200 || r.status >= 300
        data = JSON.parse(r.data)
        throw(AuthError(r.status, get(data, "message", ""), get(data, "documentation_url", "")))
    end

    auth
end

# Utility -------

function authenticate_headers(headers, auth::OAuth2)
    headers["Authorization"] = "token $(auth.token)"
end

function authenticate_headers(headers, auth::BasicAuth)
    error("authentication with BasicAuth is not fully supported")
end

function authenticate_headers(headers, auth::AnonymousAuth)
    # nothing to be done
end


text:
function lazySolve(mdp::MDP, mcts::MCTS)
  statespace, actionspace = getspaces(mdp, mcts)
  mcts.stategrid = RectangleGrid(statespace...)
  mcts.actiongrid = RectangleGrid(actionspace...)
  return internalsolve(mdp, mcts)
end

function getspaces(mdp::MDP, mcts::MCTS)
  statedim = length(mdp.statemap)
  actiondim = length(mdp.actionmap)

  stateargs = mdp.reward.argnames[1:statedim]
  actionargs = mdp.reward.argnames[1 + statedim:end]

  statespace = getspace(statedim, stateargs, mdp.statemap, mcts.statemap)
  actionspace = getspace(actiondim, actionargs, mdp.actionmap, mcts.actionmap)

  return statespace, actionspace
end

function getspace(
    dim::Int64,
    args::Vector{ASCIIString},
    lazymap::Dict{AbstractString, LazyVar},
    discmap::Dict{AbstractString, LazyDiscrete})

  space = Array(Vector{Float64}, dim)
  for ivar in 1:length(args)
    var = args[ivar]
    lazy = lazymap[var]

    if isa(lazy, RangeVar)
      space[ivar] = collect(lazy.minval : discmap[var].step : lazy.maxval)
    elseif isa(lazy, ValuesVar)
      space[ivar] = map(Float64, collect(1:length(lazy.values)))  # map to indices
    else
      error(string(
        "unknown state/action variable definition type for ", var))
    end
  end

  return space

end

include("mcts_serial.jl")


text:
"""
    MCGP(args...; kwargs...)

Monte-Carlo Gaussian Process

## Arguments
- `X::AbstractArray` : Input features, if `X` is a matrix the choice of colwise/rowwise is given by the `obsdim` keyword
- `y::AbstractVector` : Output labels
- `kernel::Kernel` : Covariance function, can be any kernel from KernelFunctions.jl
- `likelihood` : Likelihood of the model. For compatibilities, see [`Likelihood Types`](@ref likelihood_user)
- `inference` : Inference for the model, at the moment only [`GibbsSampling`](@ref) is available (see the [`Compatibility Table`](@ref compat_table))

## Keyword arguments
- `verbose::Int` : How much does the model print (0:nothing, 1:very basic, 2:medium, 3:everything)
- `optimiser` : Optimiser used for the kernel parameters. Should be an Optimiser object from the [Flux.jl](https://github.com/FluxML/Flux.jl) library, see list here [Optimisers](https://fluxml.ai/Flux.jl/stable/training/optimisers/) and on [this list](https://github.com/theogf/AugmentedGaussianProcesses.jl/tree/master/src/inference/optimisers.jl). Default is `ADAM(0.001)`
- `atfrequency::Int=1` : Choose how many variational parameters iterations are between hyperparameters optimization
- `mean=ZeroMean()` : PriorMean object, check the documentation on it [`MeanPrior`](@ref meanprior)
- `obsdim::Int=1` : Dimension of the data. 1 : X ∈ DxN, 2: X ∈ NxD
"""
mutable struct MCGP{
    T<:Real,
    TLikelihood<:AbstractLikelihood{T},
    TInference<:AbstractInference{T},
    TData<:AbstractDataContainer,
    N,
} <: AbstractGP{T,TLikelihood,TInference,N}
    data::TData
    f::NTuple{N,SampledLatent{T}} # Vector of latent GPs
    likelihood::TLikelihood
    inference::TInference
    verbose::Int64 #Level of printing information
    atfrequency::Int64
    trained::Bool
end

function MCGP(
    X::AbstractArray,
    y::AbstractVector,
    kernel::Kernel,
    likelihood::Union{AbstractLikelihood,Distribution},
    inference::AbstractInference;
    verbose::Int=0,
    optimiser=ADAM(0.01),
    atfrequency::Integer=1,
    mean::Union{<:Real,AbstractVector{<:Real},PriorMean}=ZeroMean(),
    obsdim::Int=1,
)
    X, T = wrap_X(X, obsdim)
    y, nLatent, likelihood = check_data!(y, likelihood)

    inference isa SamplingInference || error(
        "The inference object should be of type `SamplingInference` : either `GibbsSampling` or `HMCSampling`",
    )
    !isa(likelihood, GaussianLikelihood) || error(
        "For a Gaussian Likelihood you should directly use the `GP` model or the `SVGP` model for large datasets",
    )
    implemented(likelihood, inference) ||
        error("The $likelihood is not compatible or implemented with the $inference")
    !isa(likelihood, Distribution) ||
        error("Using Distributions.jl distributions is unfortunately not yet implemented")
    data = wrap_data(X, y)
    nFeatures = nSamples(data)

    if isa(optimiser, Bool)
        optimiser = optimiser ? ADAM(0.01) : nothing
    end

    if typeof(mean) <: Real
        mean = ConstantMean(mean)
    elseif typeof(mean) <: AbstractVector{<:Real}
        mean = EmpiricalMean(mean)
    end

    latentf = ntuple(_ -> SampledLatent(T, nFeatures, kernel, mean), nLatent)

    likelihood = init_likelihood(likelihood, inference, nLatent, nSamples(data))
    xview = view_x(data, 1:nSamples(data))
    yview = view_y(likelihood, data, 1:nSamples(data))
    inference = tuple_inference(
        inference, nLatent, nSamples(data), nSamples(data), nSamples(data), xview, yview
    )
    return MCGP{T,typeof(likelihood),typeof(inference),typeof(data),nLatent}(
        data, latentf, likelihood, inference, verbose, atfrequency, false
    )
end

function Base.show(io::IO, model::MCGP)
    return print(
        io,
        "Monte Carlo Gaussian Process with a $(likelihood(model)) sampled via $(inference(model)) ",
    )
end

Zviews(model::MCGP) = [input(model)]
objective(::MCGP{T}) where {T} = NaN

@traitimpl IsFull{MCGP}


text:
using Test
using TreeTools

print("empty node: \n")
n = TreeNode()
empty_t = node2tree(n);
print_tree_ascii(IO, empty_t)

print("single node: \n")
t1 = parse_newick_string("A:1")
print_tree_ascii(IO, t1)

print("nodes with branch length: \n")
t1 = parse_newick_string("(A:1, B:2)R:0")
print_tree_ascii(IO, t1)

print("unnamed nodes and no branch length: \n")
t2 = parse_newick_string("(A:1,)")
print_tree_ascii(IO, t2)

print("unnamed node with branch length: \n")
t2 = parse_newick_string("(A:1,:2):0")
print_tree_ascii(IO, t2)

print("node with no branch length: \n")
t2 = parse_newick_string("(A:1, B)R")
print_tree_ascii(IO, t2)

print("tree with branch lengths from nwk file")
t3 = node2tree(TreeTools.read_newick("$(dirname(pathof(TreeTools)))/../test/iterators/tree1.nwk"))
print_tree_ascii(IO, t3)


text:
# - - parameters - - - - - - - - - - -
system["saveU_dir"] = ""
system["verboselevel"] = 2
system["L"] = (4, 4, 4, 4)
system["Nwing"] = 1
system["Nsteps"] = 100
system["quench"] = true
system["logfile"] = "Heatbath_L04040404_beta5.7_quenched.txt"
system["initial"] = "cold"
system["Dirac_operator"] = nothing
system["log_dir"] = "./logs"
system["Nthermalization"] = 10
system["update_method"] = "HB"
system["randomseed"] = 111
system["NC"] = 3
system["BoundaryCondition"] = [1, 1, 1, -1]
system["saveU_format"] = nothing
system["β"] = 5.7
actions["use_autogeneratedstaples"] = false
actions["couplingcoeff"] = Any[]
actions["couplinglist"] = Any[]
md["Δτ"] = 0.05
md["N_SextonWeingargten"] = 2
md["SextonWeingargten"] = false
md["MDsteps"] = 20
measurement["measurement_methods"] = Dict[Dict{Any,Any}("fermiontype" => nothing,"measure_every" => 1,"methodname" => "Polyakov_loop"), Dict{Any,Any}("fermiontype" => nothing,"measure_every" => 1,"methodname" => "Plaquette")]
measurement["measurement_dir"] = "Heatbath_L04040404_beta5.7_quenched"
measurement["measurement_basedir"] = "./measurements"

actions["use_autogeneratedstaples"] = true
actions["couplingcoeff"] = [5.7,1.0,0.5]
#actions["couplinglist"] = ["plaq","rect","polyx"]

function makeplaq()
    loopset = []
    for μ=1:4
        for ν=μ:4
            if ν == μ
                continue
            end
            push!(loopset,[(μ,1),(ν,1),(μ,-1),(ν,-1)])
        end
    end
    return loopset
end

function makechair()
    loopset = []
    set1 = (1,2,3)
    set2 = (1,2,4)
    set3 = (2,3,4)
    set4 = (1,3,4)

    for set in (set1,set2,set3,set4)
        mu,nu,rho = set
        loop = [(mu,1),(nu,1),(rho,1),(mu,-1),(rho,-1),(nu,-1)]
        push!(loopset,loop)

        mu,rho,nu = set
        loop = [(mu,1),(nu,1),(rho,1),(mu,-1),(rho,-1),(nu,-1)]
        push!(loopset,loop)

        nu,rho,mu = set
        loop = [(mu,1),(nu,1),(rho,1),(mu,-1),(rho,-1),(nu,-1)]
        push!(loopset,loop)

        nu,mu,rho = set
        loop = [(mu,1),(nu,1),(rho,1),(mu,-1),(rho,-1),(nu,-1)]
        push!(loopset,loop)

        rho,mu,nu = set
        loop = [(mu,1),(nu,1),(rho,1),(mu,-1),(rho,-1),(nu,-1)]
        push!(loopset,loop)

        rho,nu,mu = set
        loop = [(mu,1),(nu,1),(rho,1),(mu,-1),(rho,-1),(nu,-1)]
        push!(loopset,loop)

    end
    return loopset
end

plaqloops = makeplaq()
rectloops = make_rectloops(2)
chairloops = makechair()
L = system["L"]
poly1loops = make_polyakovloops(1,L[1])
poly2loops = make_polyakovloops(2,L[2])
poly3loops = make_polyakovloops(3,L[3])
poly4loops = make_polyakovloops(4,L[4])

actions["coupling_loops"] = [plaqloops,rectloops,chairloops]

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


text:
"""
This struct contains the financial data of projects.
    investment_cost: Yearly capex.
    effective_investment_cost: Used to fix yearly capex.
    lag_time: Time (in years) for the project to get built.
    life_time: Total life_time of the project.
    capex_years: Number of years for capital cost recovery calculation.
    fixed_OM_cost: Annual investment cost.
    queue_cost: Vector of yearly queue costs.
    scenario_profit: Array of profits for each product through the project's life_time for each scenario.
    realized_profit: AxisArray of realized profits.
    discount_rate: Discount rate (%) for investment in the project.
    scenario_npv: Net Present Value of the project for each scenario.
    expected_npv: Expected value of NPV.
    scenario_utility: Utility of the project for each scenario.
    expected_utility: Expected value of the utility.
    annual_cashflow: Annual realized cashflow.
    ownedby: Name of investor.
"""
mutable struct Finance
    investment_cost::Vector{Float64}
    effective_investment_cost::Float64
    lag_time::Int64
    life_time::Int64
    capex_years::Int64
    fixed_OM_cost::Float64
    queue_cost::Vector{Float64}
    scenario_profit::Dict{String, Vector{AxisArrays.AxisArray{Float64, 2}}}
    realized_profit::AxisArrays.AxisArray{Float64, 2}
    discount_rate::Float64
    scenario_npv::Dict{String, Vector{Float64}}
    expected_npv::Vector{Float64}
    scenario_utility::Dict{String, Vector{Float64}}
    expected_utility::Vector{Float64}
    annual_cashflow::Vector{Float64}
    ownedby::String
end

get_investment_cost(finance_data::Finance) = finance_data.investment_cost
get_effective_investment_cost(finance_data::Finance) = finance_data.effective_investment_cost
get_lag_time(finance_data::Finance) = finance_data.lag_time
get_life_time(finance_data::Finance) = finance_data.life_time
get_capex_years(finance_data::Finance) = finance_data.capex_years
get_fixed_OM_cost(finance_data::Finance) = finance_data.fixed_OM_cost
get_queue_cost(finance_data::Finance) = finance_data.queue_cost
get_scenario_profit(finance_data::Finance) = finance_data.scenario_profit
get_realized_profit(finance_data::Finance) = finance_data.realized_profit
get_discount_rate(finance_data::Finance) = finance_data.discount_rate
get_scenario_npv(finance_data::Finance) = finance_data.scenario_npv
get_expected_npv(finance_data::Finance) = finance_data.expected_npv
get_scenario_utility(finance_data::Finance) = finance_data.scenario_utility
get_expected_utility(finance_data::Finance) = finance_data.expected_utility
get_annual_cashflow(finance_data::Finance) = finance_data.annual_cashflow
get_ownedby(finance_data::Finance) = finance_data.ownedby

function set_scenario_npv!(finance_data::Finance, scenario_name::String, iteration_year::Int64, npv::Float64)
    finance_data.scenario_npv[scenario_name][iteration_year] = npv
    return
end

function set_scenario_npv!(finance_data::Finance, scenario_name::String, npv::Vector{Float64})
    finance_data.scenario_npv[scenario_name] = npv
    return
end

function set_expected_npv!(finance_data::Finance, year::Int64, npv::Float64)
    finance_data.expected_npv[year] = npv
    return
end

function set_expected_npv!(finance_data::Finance, npv_vector::Vector{Float64})
    finance_data.expected_npv = npv_vector
end

function set_scenario_utility!(finance_data::Finance, scenario_name::String, iteration_year::Int64, utility::Float64)
    finance_data.scenario_utility[scenario_name][iteration_year] = utility
    return
end

function set_scenario_utility!(finance_data::Finance, scenario_name::String, utility::Vector{Float64})
    finance_data.scenario_utility[scenario_name] = utility
    return
end

function set_expected_utility!(finance_data::Finance, year::Int64, utility::Float64)
    finance_data.expected_utility[year] = utility
    return
end

function set_expected_utility!(finance_data::Finance, utility_vector::Vector{Float64})
    finance_data.expected_utility = utility_vector
end

function set_scenario_profit!(finance_data::Finance, scenario_name::String, profit::Vector{AxisArrays.Array{Float64, 2}})
    finance_data.scenario_profit[scenario_name] = profit
    return
end

function set_scenario_profit!(finance_data::Finance, scenario_name::String, iteration_year::Int64,  profit::AxisArrays.AxisArray{Float64, 2})
    finance_data.scenario_profit[scenario_name][iteration_year] = profit
    return
end

function set_scenario_profit!(finance_data::Finance,
                    product_name::Symbol,
                    iteration_year::Int64,
                    year_index::Int64,
                    profit::Nothing)
    return
end

function set_scenario_profit!(finance_data::Finance,
                    scenario_name::String,
                    product_name::Symbol,
                    iteration_year::Int64,
                    year_index::Int64,
                    profit::Float64)

    finance_data.scenario_profit[scenario_name][iteration_year][product_name, year_index] = profit
    return
end

function set_scenario_profit!(finance_data::Finance,
                    scenario_name::String,
                    product_name::Symbol,
                    iteration_year::Int64,
                    start_year::Int64,
                    end_year::Int64,
                    profit::Vector{Float64})

    finance_data.scenario_profit[scenario_name][iteration_year][product_name, start_year:end_year] = profit
    return
end

function set_realized_profit!(finance_data::Finance, profit::AxisArrays.AxisArray{Float64, 2})
    finance_data.realized_profit = profit
    return
end

function set_realized_profit!(finance_data::Finance,
                    product_name::Symbol,
                    year_index::Int64,
                    profit::Nothing)
    return
end

function set_realized_profit!(finance_data::Finance,
                    product_name::Symbol,
                    year_index::Int64,
                    profit::Float64)
    finance_data.realized_profit[product_name, year_index] = profit
    return
end

function set_realized_profit!(finance_data::Finance,
                    product_name::Symbol,
                    start_year::Int64,
                    end_year::Int64,
                    profit::Vector{Float64})
    finance_data.realized_profit[product_name, start_year:end_year] = profit
    return
end

function set_annual_cashflow!(finance_data::Finance, year::Int64, cashflow::Float64)
    finance_data.annual_cashflow[year] = cashflow
end


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

__precompile__(true)

module RadiationSpectra

    using DelimitedFiles
    using Distributions
    using LinearAlgebra
    using Statistics

    using DensityInterface
    using EmpiricalDistributions
    using IntervalSets
    using NamedTupleTools
    using Optim
    using Random123
    using RecipesBase
    using Requires
    using SpecialFunctions
    using StaticArrays
    using StaticUnivariatePolynomials
    using StatsBase
    using ValueShapes

    import StatsBase: fit

    """
        get_example_spectrum()::Histogram

    Returns an uncalibrated radiation spectrum for testing and demonstrating purpose.
    """
    function get_example_spectrum()::Histogram
        weights::Array{Int, 1} = readdlm(joinpath(@__DIR__, "..", "examples", "data", "hpge-spectrum.txt"))[:]
        h = Histogram(0:100:100 * length(weights), :left)
        h.weights = weights
        return h
    end

    include("general_tools.jl")

    include("ModelDensities.jl")

    include("fit.jl")
    include("UnivariateFit.jl")

    include("UvModelDensities/NormalPeakUvD.jl")

    include("PeakFinder/PeakFinder.jl")
    include("AutoCalibration/AutoCalibration.jl")

    function __init__()
        @require BAT="c0cd4b16-88b7-57fa-983b-ab80aecada7e" include("BAT.jl")
    end

    include("plot_recipes.jl")

end # module


text:
"""
    TypeD()

Lorem ipsum dolor sit amet.
"""
struct TypeD

end


text:
const SocialCompetence = Normed{UInt16,16}
const CategoryIdx = UInt8
const InCategoryIdx = PersonIdx
const CategoryProb = Float32

const FriendshipCategorySelector = AliasSampler{CategoryIdx, CategoryProb}
const FriendshipCategorySampler = AliasSampler{InCategoryIdx, CategoryProb}

struct FriendshipSampler
    categories_selectors::Vector{FriendshipCategorySelector}
    category_samplers::Vector{FriendshipCategorySampler}
    categories::Vector{Vector{PersonIdx}}
end

const max_age = 120
const num_friendship_categories = 2*max_age+2

@inline friendship_to_idx(age::Integer, gender::Bool) = gender ? UInt(age+max_age+2) : UInt(age+1)
@inline friendship_to_age_gender(idx::Integer) = (idx <= max_age+1) ? (UInt(idx-1), false) : (UInt(idx-max_age-2), true)
@inline friendship_phi(age::Integer, alpha::Real) = age <= 20 ? float(age) : float(20+(age-20.0)^alpha)

@inline function friendship_g(age1::Integer, age2::Integer, H::AbstractVector{T} where T<:Real, alpha::Real, beta::Real)
    nom = float(H[age1+1]) * float(H[age2+1]) * exp( -0.08 * (age1 + age2) )
    denom = 1.0 + 0.2 * abs(friendship_phi(age1, alpha) - friendship_phi(age2, alpha)) ^ beta
    return nom/denom
end

function FriendshipSampler(
  ages::AbstractVector{T} where T<:Integer,
  genders::AbstractVector{Bool},
  social_competences::AbstractVector{T} where T <: Real,
  alpha::Real = 0.75,
  beta::Real = 1.6)

  num_individuals = length(ages)
  @assert num_individuals == length(genders) == length(social_competences)

  categories = [Vector{PersonIdx}() for _ in 1:num_friendship_categories]
  for ii in 1:num_individuals
      push!(categories[friendship_to_idx(ages[ii], genders[ii])], ii)
  end

  H = Vector{Rational{PersonIdx}}(undef, max_age+1)
  for age in 0:max_age
    category_f = categories[friendship_to_idx(age, false)]
    category_t = categories[friendship_to_idx(age, true)]
    H[age+1] = (length(category_f) + length(category_t)) // num_individuals #should fit despite being Rational{Int}
  end

  categories_selectors = Vector{FriendshipCategorySelector}(undef, num_friendship_categories)

  P = Vector{CategoryProb}(undef, num_friendship_categories)
  for idx in 1:num_friendship_categories
      age, gender = friendship_to_age_gender(idx)
      for idx2 in 1:num_friendship_categories
          age2, gender2 = friendship_to_age_gender(idx2)
          P[idx2] = friendship_g(age, age2, H, alpha, beta) * (gender == gender2 ? 1.2 : 0.8)
      end
      categories_selectors[idx] = AliasSampler(CategoryIdx, P)#FriendshipCategorySelector(P)
  end

  category_samplers = Vector{FriendshipCategorySampler}(undef, num_friendship_categories)
  for cat_idx in 1:num_friendship_categories
      category = categories[cat_idx]
      alias_sampler_par = Vector{CategoryProb}(undef, length(category))
      for pers_idx in 1:length(category)
          person_id = category[pers_idx]
          alias_sampler_par[pers_idx] = social_competences[person_id]
      end
      category_samplers[cat_idx] = AliasSampler(InCategoryIdx, alias_sampler_par)
  end

  return FriendshipSampler(categories_selectors, category_samplers, categories)
end

FriendshipSampler(population::DataFrame, alpha::Real = 0.75, beta::Real = 1.6) =
    FriendshipSampler(population.age, population.gender, population.social_competence, alpha, beta)

function friend_sample(fs::FriendshipSampler, age::Integer, gender::Bool, rng=Random.GLOBAL_RNG)
    category = a_sample(fs.categories_selectors[friendship_to_idx(age, gender)], rng)
    return fs.categories[category][a_sample(fs.category_samplers[category], rng)]
end

struct FriendshipKernelParams
  kernel_constant::Real
  social_competences::Vector{SocialCompetence}
  sampler::FriendshipSampler
end

function FriendshipKernelParams(
  kernel_constant::Real,
  ages::AbstractVector{Age},
  genders::AbstractVector{Bool},
  social_competences::Vector{SocialCompetence}) #concrete type for direct storage

  @assert all(0 .<= ages .<= max_age)

  FriendshipKernelParams(
    kernel_constant,
    social_competences,
    FriendshipSampler(
      ages,
      genders,
      social_competences))
end
socialcompetence(p::FriendshipKernelParams, person_id::Integer) = p.social_competences[person_id]


text:
function show(io::IO, object::testOutput)
    if object.method != nothing
        @printf("%s\n", object.method)
    end
    if object.h0 != nothing
        @printf(" Alternative hypothesis:   % s\n", object.h0)
    end
    if object.n != nothing
        @printf(" Sample size:         ")
        for i = 1:length(object.n)
            @printf("% d  ", object.n[i])
        end
        @printf("\n")
    end
    if object.df != nothing
            @printf(" Degrees of freedom:  ")
        for i = 1:length(object.df)
            if iseltype(object.df[i], Integer)
                @printf(" %d  ", object.df[i])
            elseif iseltype(object.df, FloatingPoint)
                @printf("% .2f  ", object.df[i])
            end
        end
       @printf("\n")
    end
    if object.estimate != nothing
        @printf(" Estimate:            % .6f\n", object.estimate)
    end
    if object.se != nothing
        @printf(" Standard error:      % .6f\n", object.se)
    end
    if object.statistic != nothing
        if ndims(object.statistic)==0
            @printf(" Statistic:           % .6f\n", object.statistic)
        else ndims(object.statistic)==2
            @printf(" Statistic:            % 7s  % 7s  % 9s\n", "Group 1", "Group 2", "Statistic")
            for i = 1:size(object.statistic, 1)
                @printf("                       % 7d  % 7d  %  7.6f\n",
                        object.statistic[i,1], object.statistic[i,2], object.statistic[i,3])
            end
        end
    end
    if object.crit != nothing
        @printf(" Critical value:      % .6f\n", object.crit)
    end
    if object.ci != nothing
        if iseltype(object.ci, Integer)
            @printf(" Confidence interval: % d      % d\n", object.ci[1], object.ci[2])
        else
            @printf(" Confidence interval: % .6f      % .6f\n", object.ci[1], object.ci[2])
        end
    end
    if object.p != nothing
        if object.p < 10e-16
            @printf(" p value:             < 10e-16\n")
        elseif object.p < 1
            @printf(" p value:             % .6f\n", object.p)
        else
            @printf(" p value:              1.0")
        end
    end
end

type binomciOutput
    phat::FloatingPoint
    confint::Array
    n::Int
end

function show(io::IO, object::outOutput)
    if object.method!=nothing
        println(io, "$(object.method)")
    end
    @printf("Outlier ID:         ")
    if length(object.outid) == 0
        @printf("\n")
    else
        for i = 1:length(object.outid)
            if i < length(object.outid)
                @printf("%d, ", object.outid[i])
            else
                @printf("%d\n", object.outid[i])
            end
        end
    end
    @printf("Outlier value:      ")
    if length(object.outid) == 0
        @printf("\n")
    else
        if iseltype(object.outval, Integer)
            for i = 1:length(object.outval)
                if i < length(object.outval)
                    @printf("%d, ", object.outval[i])
                else
                    @printf("%d\n", object.outval[i])
                end
            end
        else
            for i=1:length(object.outval)
                if i < length(object.outval)
                    @printf("%.5f, ", object.outval[i])
                else
                    @printf("%.5f\n", object.outval[i])
                end
            end
        end
    end
    @printf("Number of outliers: %d\n", object.nout)
    @printf("Non-outlier ID:     ")
    for i=1:length(object.keepid)
        if i < length(object.keepid)
            @printf("%d, ", object.keepid[i])
        else
            @printf("%d\n", object.keepid[i])
        end
    end
end

function show(io::IO, object::indtOutput)
    if object.method!=nothing
        println(io, "$(object.method)")
    end
    if object.flag!=2
            @printf("Kolmogorov-Smirnov test statistic:  % .6f\n", object.dstat)
        if object.pval_d >=1
            @printf("Kolmogorov-Smirnov p value:          1.0\n")
        elseif object.pval_d < 10e-16
            @printf("Kolmogorov-Smirnov p value:         < 10e-16\n")
        else
            @printf("Kolmogorov-Smirnov p value:         % .6f\n", object.pval_d)
        end
    end
    if object.flag!=1
            @printf("Cramer-von Mises test statistic:    % .6f\n", object.wstat)
        if object.pval_w >=1
            @printf("Cramer-von Mises p value:            1.0\n")
        elseif object.pval_w < 10e-16
            @printf("Cramer-von Mises p value:           < 10e-16\n")
        else
            @printf("Cramer-von Mises p value:           % .6f\n", object.pval_w)
        end
    end
end

function show(io::IO, object::idealfOutput)
    if object.method!=nothing
        println(io, "$(object.method)")
    end
    @printf("Lower quartile:  %.6f\n", object.lower_quartile)
    @printf("Upper quartile:  %.6f\n", object.upper_quartile)
end

function show(io::IO, object::binomciOutput)
   @printf(" phat:                %.4f\n", object.phat)
   @printf(" confidence interval: %.4f   %.4f\n", object.confint[1], object.confint[2])
   @printf(" Sample size          %d\n", object.n)
end

function show(io::IO, object::indirectTestOutput)
    @printf("\nTESTS OF INDIRECT EFFECT\n\n")
    @printf("Sample size:  %d\n", object.n)
    @printf("Number of bootstrap samples: %d\n\n", object.nboot)
    @printf("DIRECT AND TOTAL EFFECTS\n")
    @printf("          % 8s  %9s    %8s   %8s\n", "Estimate", "Std.Error", "t value", "Pr(>|t|)")
    @printf("b(YX):   % 8.6f   %8.6f   % 8.6f   %8.6f\n",
            object.regfit[1,1], object.regfit[1,2], object.regfit[1,3], object.regfit[1,4])
    @printf("b(MX):   % 8.6f   %8.6f   % 8.6f   %8.6f\n",
            object.regfit[2,1], object.regfit[2,2], object.regfit[2,3], object.regfit[2,4])
    @printf("b(YM.X): % 8.6f   %8.6f   % 8.6f   %8.6f\n",
            object.regfit[3,1], object.regfit[3,2], object.regfit[3,3], object.regfit[3,4])
    @printf("b(YX.M): % 8.6f   %8.6f   % 8.6f   %8.6f\n\n",
            object.regfit[4,1], object.regfit[4,2], object.regfit[4,3], object.regfit[4,4])

    @printf("INDIRECT EFFECT AND SIGNIFICANCE USING NORMAL DISTRIBUTION\n")
    @printf("          % 8s  %9s    %8s   %8s   %8s  %8s \n",
            "Estimate", "Std.Error", "z score", "CI lo", "CI hi", "Pr(>|z|)")
    @printf("Sobel:   % 8.6f   %8.6f   % 8.6f  % 8.6f  % 8.6f  %8.6f\n\n",
            object.sobel[1,1], object.sobel[1,2], object.sobel[1,3],
            object.sobel[1,4], object.sobel[1,5], object.sobel[1,6])

    @printf("BOOTSTRAP RESULTS OF INDIRECT EFFECT\n")
    @printf("          % 8s  %9s    %8s   %8s   %8s \n",
            "Estimate", "Std.Error", "CI lo", "CI hi", "P value")
    @printf("Effect:  % 8.6f   %8.6f   % 8.6f  % 8.6f   %8.6f\n",
            object.bootest, object.bootse, object.bootci[1], object.bootci[2], object.p)
end


text:
for N in [Float64, Rational{Int}, Float32]
    # random ball
    rand(Ball1)

    # 1D Ball1
    b = Ball1(N[0], N(1))
    # dimension
    @test dim(b) == 1
    # support vector
    d = N[1]
    @test σ(d, b) == N[1]
    d = N[-1]
    @test σ(d, b) == N[-1]

    # 2D Ball1
    b = Ball1(N[0, 0], N(1))
    # dimension
    @test dim(b) == 2
    # support vector
    d = N[1, 0]
    @test σ(d, b) == N[1, 0]
    d = N[-1, 0]
    @test σ(d, b) == N[-1, 0]
    d = N[0, 1]
    @test σ(d, b) == N[0, 1]
    d = N[0, -1]
    @test σ(d, b) == N[0, -1]
    d = N[0, 0]
    @test σ(d, b) == N[0, 0]

    # 2D Ball1 not 0-centered
    b = Ball1(N[1, 2], N(1))
    # dimension
    @test dim(b) == 2
    # support vector
    d = N[1, 0]
    @test σ(d, b) == N[2, 2]
    d = N[-1, 0]
    @test σ(d, b) == N[0, 2]
    d = N[0, 1]
    @test σ(d, b) == N[1, 3]
    d = N[0, -1]
    @test σ(d, b) == N[1, 1]

    # 2D Ball1 radius != 1
    b = Ball1(N[0, 0], N(2))
    # dimension
    @test dim(b) == 2
    # support vector
    d = N[1, 0]
    @test σ(d, b) == N[2, 0]
    d = N[-1, 0]
    @test σ(d, b) == N[-2, 0]
    d = N[0, 1]
    @test σ(d, b) == N[0, 2]
    d = N[0, -1]
    @test σ(d, b) == N[0, -2]

    # center
    @test center(b) == N[0, 0]

    # boundedness
    @test isbounded(b)

    # isempty
    @test !isempty(b)

    # isuniversal
    answer, w = isuniversal(b, true)
    @test !isuniversal(b) && !answer && w ∉ b

    # an_element & membership function
    e = an_element(b)
    answer = e ∈ b
    @test answer == ∈(e, b, true) == ∈(e, b, false) == true
    e += N[2 * b.radius, 1]
    answer = e ∈ b
    @test answer == ∈(e, b, true) == ∈(e, b, false) == false

    # translation
    @test translate(b, N[1, 2]) == Ball1(N[1, 2], N(2))

    # vertices_list
    vl = vertices_list(b)
    @test ispermutation(vl, [N[2, 0], N[0, 2], N[-2, 0], N[0, -2]])

    # check that vertices_list for zero radius doesn't repeat vertices
    b = Ball1(N[1, 2], N(0))
    vl = vertices_list(b)
    @test vl == [center(b)]

    # list of constraints
    b = Ball1(N[0, 0], N(1))
    clist = constraints_list(b)
    @test ispermutation(clist, [HalfSpace(N[1, 1], N(1)),  # x + y <= 1
                                HalfSpace(N[-1, 1], N(1)),  # -x + y <= 1
                                HalfSpace(N[1, -1], N(1)), # x - y <= 1
                                HalfSpace(N[-1, -1], N(1))])  # x + y >= -1
end


text:
# Minimal sanity checks for MakieLayout
@testset "Layoutables constructors" begin
    fig = Figure()
    ax = Axis(fig[1, 1])
    cb = Colorbar(fig[1, 2])
    gl2 = fig[2, :] = MakieLayout.GridLayout()
    bu = gl2[1, 1] = Button(fig)
    sl = gl2[1, 2] = Slider(fig)

    scat = scatter!(ax, rand(10))
    le = gl2[1, 3] = Legend(fig, [scat], ["scatter"])

    to = gl2[1, 4] = Toggle(fig)
    te = fig[0, :] = Label(fig, "A super title")
    me = fig[end + 1, :] = Menu(fig, options=["one", "two", "three"])
    tb = fig[end + 1, :] = Textbox(fig)
    is = fig[end + 1, :] = IntervalSlider(fig)
    @test true
end

@testset "deleting from axis" begin
    f = Figure()
    ax = Axis(f[1, 1])
    sc = scatter!(ax, randn(100, 2))
    li = lines!(ax, randn(100, 2))
    hm = heatmap!(ax, randn(20, 20))
    # axis contains 3 + 1 plots, one for the zoomrectangle
    @test length(ax.scene.plots) == 4
    delete!(ax, sc)
    @test length(ax.scene.plots) == 3
    @test sc ∉ ax.scene.plots
    empty!(ax)
    @test isempty(ax.scene.plots)
end

@testset "zero heatmap" begin
    xs = LinRange(0, 20, 10)
    ys = LinRange(0, 15, 10)
    zs = zeros(length(xs), length(ys))

    fig = Figure()
    _, hm = heatmap(fig[1, 1], xs, ys, zs)
    cb = Colorbar(fig[1, 2], hm)

    @test hm.attributes[:colorrange][] == (-.5, .5)
    @test cb.limits[] == (-.5, .5)

    hm.attributes[:colorrange][] = Float32.((-1, 1))
    @test cb.limits[] == (-1, 1)

    cb.limits[] = Float32.((-2, 2))
    @test hm.attributes[:colorrange][] == (-2, 2)
end

@testset "Axis limits basics" begin
    f = Figure()
    ax = Axis(f[1, 1], limits = (nothing, nothing))
    ax.targetlimits[] = BBox(0, 10, 0, 20)
    @test ax.finallimits[] == BBox(0, 10, 0, 20)
    @test ax.limits[] == (nothing, nothing)
    xlims!(ax, -10, 10)
    @test ax.limits[] == ((-10, 10), nothing)
    @test ax.finallimits[] == BBox(-10, 10, 0, 20)
    ylims!(ax, -20, 30)
    @test ax.limits[] == ((-10, 10), (-20, 30))
    @test ax.finallimits[] == BBox(-10, 10, -20, 30)
    limits!(ax, -5, 5, -10, 10)
    @test ax.finallimits[] == BBox(-5, 5, -10, 10)
    @test ax.limits[] == ((-5, 5), (-10, 10))
    ax.limits[] = (nothing, nothing)
    ax.xautolimitmargin = (0, 0)
    ax.yautolimitmargin = (0, 0)
    scatter!(Point2f[(0, 0), (1, 2)])
    @test ax.limits[] == (nothing, nothing)
    @test ax.targetlimits[] == BBox(0, 1, 0, 2)
    @test ax.finallimits[] == BBox(0, 1, 0, 2)
    scatter!(Point2f(3, 4))
    @test ax.limits[] == (nothing, nothing)
    @test ax.targetlimits[] == BBox(0, 3, 0, 4)
    @test ax.finallimits[] == BBox(0, 3, 0, 4)
    limits!(ax, -1, 1, 0, 2)
    @test ax.limits[] == ((-1, 1), (0, 2))
    @test ax.targetlimits[] == BBox(-1, 1, 0, 2)
    @test ax.finallimits[] == BBox(-1, 1, 0, 2)
    scatter!(Point2f(5, 6))
    @test ax.limits[] == ((-1, 1), (0, 2))
    @test ax.targetlimits[] == BBox(-1, 1, 0, 2)
    @test ax.finallimits[] == BBox(-1, 1, 0, 2)
    autolimits!(ax)
    @test ax.limits[] == (nothing, nothing)
    @test ax.targetlimits[] == BBox(0, 5, 0, 6)
    @test ax.finallimits[] == BBox(0, 5, 0, 6)
    xlims!(ax, [-10, 10])
    @test ax.limits[] == ([-10, 10], nothing)
    @test ax.targetlimits[] == BBox(-10, 10, 0, 6)
    @test ax.finallimits[] == BBox(-10, 10, 0, 6)
    scatter!(Point2f(11, 12))
    @test ax.limits[] == ([-10, 10], nothing)
    @test ax.targetlimits[] == BBox(-10, 10, 0, 12)
    @test ax.finallimits[] == BBox(-10, 10, 0, 12)
    autolimits!(ax)
    ylims!(ax, [5, 7])
    @test ax.limits[] == (nothing, [5, 7])
    @test ax.targetlimits[] == BBox(0, 11, 5, 7)
    @test ax.finallimits[] == BBox(0, 11, 5, 7)
    scatter!(Point2f(-5, -7))
    @test ax.limits[] == (nothing, [5, 7])
    @test ax.targetlimits[] == BBox(-5, 11, 5, 7)
    @test ax.finallimits[] == BBox(-5, 11, 5, 7)
end

@testset "Colorbar plot object kwarg clash" begin
    for attr in (:colormap, :limits)
        f, ax, p = scatter(1:10, 1:10, color = 1:10, colorrange = (1, 10))
        Colorbar(f[2, 1], p)
        @test_throws ErrorException Colorbar(f[2, 1], p; Dict(attr => nothing)...)
    end

    for attr in (:colormap, :limits, :highclip, :lowclip)
        for F in (heatmap, contourf)
            f, ax, p = F(1:10, 1:10, randn(10, 10))
            Colorbar(f[1, 2], p)
            @test_throws ErrorException Colorbar(f[1, 3], p; Dict(attr => nothing)...)
        end
    end
end

@testset "Tick functions" begin
    automatic = Makie.automatic
    Automatic = Makie.Automatic

    get_ticks = MakieLayout.get_ticks
    get_tickvalues = MakieLayout.get_tickvalues
    get_ticklabels = MakieLayout.get_ticklabels

    for func in [identity, log, log2, log10, Makie.logit]
        tup = ([1, 2, 3], ["a", "b", "c"])
        @test get_ticks(tup, func, automatic, 0, 5) == tup

        rng = 1:5
        @test get_ticks(rng, func, automatic, 0, 5) == ([1, 2, 3, 4, 5], ["1", "2", "3", "4", "5"])

        numbers = [1.0, 1.5, 2.0]
        @test get_ticks(numbers, func, automatic, 0, 5) == (numbers, ["1.0", "1.5", "2.0"])

        @test get_ticks(numbers, func, xs -> string.(xs) .* "kg", 0, 5) == (numbers, ["1.0kg", "1.5kg", "2.0kg"])

        @test get_ticks(WilkinsonTicks(5), identity, automatic, 1, 5) == ([1, 2, 3, 4, 5], ["1", "2", "3", "4", "5"])
    end
end

@testset "Colorbars" begin
    fig = Figure()
    hmap = heatmap!(Axis(fig[1, 1]), rand(4, 4))
    cb1 = Colorbar(fig[1,2], hmap; height = Relative(0.65))
    @test cb1.height[] == Relative(0.65)
end


text:
include("cross_terms.jl")
include("level_bools.jl")

function create_matrices(level, h, h_lag, head_line)
    ## Load full Data Set
    data = readtable("../../../../data/source/PCEPI_Detail.csv")::DataFrame
    PCE = data[4:end, :]::DataFrame
    agg = convert(Array{Int}, data[2, 2:end])
    term = convert(Array{Int}, data[3, 2:end])
    data_agg = PCE[:, find([1 level_bools(level, agg, term)])]

    ## Compute h-step ahead headline inflation as in Gamber and Smith (2016)
    ## Define left-hand-side lagged variables, divided by (h_lag / 12) to annualize
    PCE_hl = convert(Array{Float64}, PCE[:, [2]])
    π_hl = ((PCE_hl[h_lag + h + 1:end] ./ PCE_hl[h_lag + 1:end - h]) - 1) * (100 / (h / 12))
    ## Define right-hand-side lagged variables
    π_hl_lag = ((PCE_hl[h_lag + 1:end] ./ PCE_hl[1:end - h_lag]) - 1) * (100 / (h_lag / 12))

    ## Compute h-step ahead headline inflation as in Gamber and Smith (2016)
    ## Define left hand-side-variables, divided by (h_lag / 12) to annualize
    PCE_core = convert(Array{Float64}, PCE[:, [3]])
    π_core = ((PCE_core[h_lag + h + 1:end] ./ PCE_core[h_lag + 1:end - h]) - 1) * 100 / (h / 12)
    ## Define right hand side lagged variables
    π_core_lag = ((PCE_core[h_lag + 1:end] ./ PCE_core[1:end - h_lag]) - 1) * (100 / (h_lag / 12))

    ## Define the right-hand-side disaggregate variables
    PCE_agg = convert(Array{Float64}, data_agg[:, 2:end])
    π_agg_lag = ((PCE_agg[h_lag + 1:end, :] ./ PCE_agg[1:end - h_lag, :]) - 1) * (100 / (h_lag / 12))

    Y = head_line .* π_hl + (1 - head_line) .* π_core
    Y_lag = head_line .* π_hl_lag + (1 - head_line) .* π_core_lag
    X_lag = π_agg_lag
    X_lag2 = second_order_cross(X_lag)

    return Y, Y_lag, X_lag, X_lag2
end

max_stars_repo_name: jungangl/forecast-machine-learning

max_issues_repo_name: jungangl/forecast-machine-learning

max_forks_repo_name: jungangl/forecast-machine-learning


text:
module DistMesh2D

using GeometricalPredicates, Deldir, DataFrames

include("scaler.jl")
include("util.jl")
include("distancefunctions.jl")
include("triangle.jl")
include("edges.jl")
include("pointstoforces.jl")
include("finalpoints.jl")
include("distmesh.jl")

export distmesh2d, drectangle, dcircle

end

max_stars_repo_name: jstarczewski/julia-num

max_issues_repo_name: jstarczewski/julia-num

max_forks_repo_name: jstarczewski/julia-num


text:
using ProperOrthogonalDecomposition
using Test
using StableRNGs
using Statistics
using DelimitedFiles

# Define test matrix X with dimensions n×m, where n is number of data poitns and
# m is the number of snapshots
# W is the weight matrix representing the cell volume

rng = StableRNGs.StableRNG(1)
X = rand(rng,1000,10)
X .+= 10
rng = StableRNGs.StableRNG(1)
W = rand(rng,1:0.2:10,1000)

@testset "Standard POD" begin

    PODbase, Σ = POD(X)
    PODbaseEig, Σeig = PODeigen(X)
    PODbaseSvd, Σsvd = PODsvd(X)

    meanPODbaseEig, meanΣeig = PODeigen(copy(X), subtractmean = true)
    meanPODbaseSvd, meanΣsvd = PODsvd(copy(X), subtractmean = true)

    Σ₁ = 1050.0370061587332
    Σ₂ = 8.413650551681984
    meanΣ₁ = 9.964071096830972
    meanΣ₂ = 8.42358075601378

    @testset "POD using eigenvalue decomposition" begin
        @test Σeig[1] ≈ Σ₁
        @test Σeig[end] ≈ Σ₂
        @test meanΣeig[1] ≈ meanΣ₁
        @test meanΣeig[end-1] ≈ meanΣ₂

    end

    @testset "POD using SVD" begin
        @test Σsvd[1] ≈ Σ₁
        @test Σsvd[end] ≈ Σ₂
        @test meanΣsvd[1] ≈ meanΣ₁
        @test meanΣsvd[end-1] ≈ meanΣ₂
    end

    @testset "Default method" begin
        @test maximum(abs.(abs.(PODbaseSvd.coefficients) .- abs.(PODbase.coefficients))) ≈ 0 atol=1e-8
        @test maximum(abs.(abs.(PODbaseSvd.modes) .- abs.(PODbase.modes))) ≈ 0 atol=1e-8
        @test maximum(abs.(Σsvd.-Σ)) ≈ 0 atol=1e-8
    end

    @testset "Method equality of SVD and Eig" begin
        @test maximum(abs.(abs.(PODbaseSvd.coefficients) .- abs.(PODbaseEig.coefficients))) ≈ 0 atol=1e-8
        @test maximum(abs.(abs.(PODbaseSvd.modes) .- abs.(PODbaseEig.modes))) ≈ 0 atol=1e-8
        @test maximum(abs.(abs.(meanPODbaseSvd.coefficients) .- abs.(meanPODbaseEig.coefficients))) ≈ 0 atol=1e-7
        @test maximum(abs.(abs.(meanPODbaseSvd.modes[:,1:end-1]) .- abs.(meanPODbaseEig.modes[:,1:end-1]))) ≈ 0 atol=1e-7
        @test maximum(abs.(Σeig.-Σsvd)) ≈ 0 atol=1e-8
    end

    @testset "Rebuild solution" begin
        @test PODbaseEig.modes*PODbaseEig.coefficients ≈ X
        @test PODbaseSvd.modes*PODbaseSvd.coefficients ≈ X
        @test meanPODbaseEig.modes*meanPODbaseEig.coefficients ≈ X .- mean(X,dims=2)
        @test meanPODbaseSvd.modes*meanPODbaseSvd.coefficients ≈ X .- mean(X,dims=2)
    end
end

@testset "Weighted POD" begin

    PODbase, Σ = POD(X, W)
    PODbaseEig, Σeig = PODeigen(X,W)
    PODbaseSvd, Σsvd = PODsvd(X,W)

    meanPODbaseEig, meanΣeig = PODeigen(copy(X), W, subtractmean = true)
    meanPODbaseSvd, meanΣsvd = PODsvd(copy(X), W, subtractmean = true)

    Σ₁ = 2462.447359829031
    Σ₂ = 19.570713247455668
    meanΣ₁ = 23.470360231176947
    meanΣ₂ = 19.595332053746322

    @testset "POD using eigenvalue decomposition" begin
        @test Σeig[1] ≈ Σ₁
        @test Σeig[end] ≈ Σ₂
        @test meanΣeig[1] ≈ meanΣ₁
        @test meanΣeig[end-1] ≈ meanΣ₂

    end

    @testset "POD using SVD" begin
        @test Σsvd[1] ≈ Σ₁
        @test Σsvd[end] ≈ Σ₂
        @test meanΣsvd[1] ≈ meanΣ₁
        @test meanΣsvd[end-1] ≈ meanΣ₂
    end

    @testset "Default method with weights" begin
        @test maximum(abs.(abs.(PODbaseSvd.coefficients) .- abs.(PODbase.coefficients))) ≈ 0 atol=1e-8
        @test maximum(abs.(abs.(PODbaseSvd.modes) .- abs.(PODbase.modes))) ≈ 0 atol=1e-8
        @test maximum(abs.(Σsvd.-Σ)) ≈ 0 atol=1e-8
    end

    @testset "Method equality of SVD and Eig with weights" begin
        @test maximum(abs.(abs.(PODbaseSvd.coefficients) .- abs.(PODbaseEig.coefficients))) ≈ 0 atol=1e-8
        @test maximum(abs.(abs.(PODbaseSvd.modes) .- abs.(PODbaseEig.modes))) ≈ 0 atol=1e-8
        @test maximum(abs.(Σeig.-Σsvd)) ≈ 0 atol=1e-8
    end

    @testset "Rebuild solution with weights" begin
        @test PODbaseEig.modes*PODbaseEig.coefficients ≈ X
        @test PODbaseSvd.modes*PODbaseSvd.coefficients ≈ X
    end
end

@testset "Mode convergence" begin
    A = [0.010164609073873544 0.006921747306424847 0.004321454495431182 0.0;
         0.7274414547632297 0.6102769856096316 0.2273767387120779 0.0;
         1.1710785171906266 1.1955870815537555 0.12812520121296264 0.0]
    W₂ = ones(1000)

    testdataPath = joinpath(dirname(pathof(ProperOrthogonalDecomposition)),"..","test","testdata.csv")

    @testset "Convergence of number of included snapshots" begin

        @test modeConvergence(X,PODeigen,[1:4,1:6,1:8,1:10],3) ≈ A
        @test modeConvergence(X,PODsvd,[1:4,1:6,1:8,1:10],3) ≈ A
        @test modeConvergence!(()->readdlm(testdataPath, ','),PODeigen!,[1:4,1:6,1:8,1:10],3) ≈ A
        @test modeConvergence!(()->readdlm(testdataPath, ','),PODsvd!,[1:4,1:6,1:8,1:10],3) ≈ A

    end

    @testset "Convergence of number of included snapshots with one weights" begin

        @test modeConvergence(X,x->PODeigen(x,W₂),[1:4,1:6,1:8,1:10],3) ≈ A
        @test modeConvergence(X,x->PODsvd(x,W₂),[1:4,1:6,1:8,1:10],3) ≈ A
        @test modeConvergence!(()->readdlm(testdataPath, ','),x->PODeigen!(x,W₂),[1:4,1:6,1:8,1:10],3) ≈ A
        @test modeConvergence!(()->readdlm(testdataPath, ','),x->PODsvd!(x,W₂),[1:4,1:6,1:8,1:10],3) ≈ A

    end

end


text:
module MenuItems

Base.Experimental.@optlevel 2

using Redux
using CImGui

## actions
"""
    AbstractGenericMenuItemAction <: AbstractSyncAction
Abstract supertype for all kinds of menu item actions. Literally, its subtype could be anything.
"""
abstract type AbstractGenericMenuItemAction <: AbstractSyncAction end

"""
    AbstractMenuItemAction <: AbstractGenericMenuItemAction
Abstract supertype for all menu item actions. Any subtypes of this type should has a
`label` field or at least `x.label` should work.
"""
abstract type AbstractMenuItemAction <: AbstractGenericMenuItemAction end

# for determining whether we should apply the reducer function to certain action types.
struct Valid end
struct Invalid end

is_valid(a::AbstractSyncAction) = Invalid
is_valid(a::AbstractGenericMenuItemAction) = Invalid
is_valid(a::AbstractMenuItemAction) = Valid

_get_label(a::AbstractSyncAction, ::Type{Invalid}) = "___REDUXIMGUI_RESERVED_DUMMY_ACTION_LABEL"
_get_label(a::AbstractSyncAction, ::Type{Valid}) = a.label
get_label(a::AbstractSyncAction) = _get_label(a, is_valid(a))

"""
    Rename(label, new_label)
Change widget's label to `new_label`.
Note that renaming may also change the identifier, please refer to `help?> CImGui.PushID`.
"""
struct Rename <: AbstractMenuItemAction
    label::String
    new_label::String
end

"""
    SetTriggeredTo(label, is_triggered)
Update widget's `is_triggered` state.
"""
struct SetTriggeredTo <: AbstractMenuItemAction
    label::String
    is_triggered::Bool
end

"""
    ChangeShortcut(label, shortcut)
Change shortcut string to `shortcut`.
"""
struct ChangeShortcut <: AbstractMenuItemAction
    label::String
    shortcut::String
end

"""
    Enable(label)
Enable the menu item.
"""
struct Enable <: AbstractMenuItemAction
    label::String
end

"""
    Disable(label)
Disable the menu item.
"""
struct Disable <: AbstractMenuItemAction
    label::String
end

"""
    AddMenuItem(label::AbstractString, shortcut = "", is_enabled = true)
Add a menu item.
"""
struct AddMenuItem <: AbstractMenuItemAction
    label::String
    shortcut::String
    is_enabled::Bool
end
AddMenuItem(label::AbstractString, shortcut = "") = AddMenuItem(label, shortcut, true)

"""
    DeleteMenuItem(label::AbstractString)
Delete the menu item.
"""
struct DeleteMenuItem <: AbstractMenuItemAction
    label::String
end

## state
"""
    AbstractGenericMenuItem <: AbstractImmutableState
Abstract supertype for all kinds of menu items. Literally, its subtype could be anything.
"""
abstract type AbstractGenericMenuItem <: AbstractImmutableState end

"""
    AbstractMenuItem <: AbstractGenericMenuItem
Abstract supertype for all menu items. Any subtypes of this type should has a
`label` field or at least `x.label` should work.
"""
abstract type AbstractMenuItem <: AbstractGenericMenuItem end

is_valid(s::AbstractImmutableState) = Invalid
is_valid(s::AbstractGenericMenuItem) = Invalid
is_valid(s::AbstractMenuItem) = Valid

_get_label(s::AbstractImmutableState, ::Type{Invalid}) = "___REDUXIMGUI_RESERVED_DUMMY_STATE_LABEL"
_get_label(s::AbstractImmutableState, ::Type{Valid}) = s.label
get_label(s::AbstractImmutableState) = _get_label(s, is_valid(s))

"""
    MenuItem(label::AbstractString, shortcut = "", is_enabled = true)
A menu item state which contains a label a.k.a the identifier, a shortcut string,
a flag value `is_enabled` and a flag value `is_triggered` that records the state of the latest poll events.
"""
struct MenuItem <: AbstractMenuItem
    label::String
    shortcut::String
    is_enabled::Bool
    is_triggered::Bool
end
MenuItem(label::AbstractString, shortcut = "", is_enabled = true) =
    MenuItem(label, shortcut, is_enabled, false)

get_shortcut(s::MenuItem) = s.shortcut
is_triggered(s::MenuItem) = s.is_triggered
is_enabled(s::MenuItem) = s.is_enabled

# reducers
reducer(state::AbstractState, action::AbstractAction) = state
reducer(state::Vector{<:AbstractState}, action::AbstractAction) = state
reducer(state::Dict{String,<:AbstractState}, action::AbstractAction) = state

reducer(s::MenuItem, a::Rename) = MenuItem(a.new_label, s.shortcut, s.is_enabled, s.is_triggered)
reducer(s::MenuItem, a::ChangeShortcut) = MenuItem(s.label, a.shortcut, s.is_enabled, s.is_triggered)
reducer(s::MenuItem, a::Enable) = MenuItem(s.label, s.shortcut, true, s.is_triggered)
reducer(s::MenuItem, a::Disable) = MenuItem(s.label, s.shortcut, false, s.is_triggered)
reducer(s::MenuItem, a::SetTriggeredTo) = MenuItem(s.label, s.shortcut, s.is_enabled, a.is_triggered)

reducer(s::Dict{String,<:AbstractImmutableState}, a::AbstractMenuItemAction) =
    Dict(k => (get_label(v) == get_label(a) ? reducer(v, a) : v) for (k, v) in s)

reducer(s::Vector{<:AbstractImmutableState}, a::AbstractMenuItemAction) = map(s) do s
    get_label(s) == get_label(a) ? reducer(s, a) : s
end
reducer(s::Vector{<:AbstractImmutableState}, a::AddMenuItem) =
    [s..., MenuItem(a.label, a.shortcut, a.is_enabled, false)]
reducer(s::Vector{<:AbstractImmutableState}, a::DeleteMenuItem) =
    filter(s -> get_label(s) != get_label(a), s)

# UI
"""
    (::MenuItem)(store::AbstractStore, get_state=Redux.get_state, chain_action=identity) -> Bool
Return `true` when triggered/activated.
`get_state` is a router function that tells how to find the target state from `store`.
`chain_action` is for chaining upstream actions.
"""
function (::MenuItem)(store::AbstractStore, get_state=Redux.get_state, chain_action=identity)
    s = get_state(store)
    is_activated = CImGui.MenuItem(get_label(s), s.shortcut, false, s.is_enabled)
    dispatch!(
        store,
        SetTriggeredTo(get_label(s), is_activated) |> chain_action,
    )
    return is_activated
end

end # module


text:
using Distributions, DataFrames, DataFramesMeta, Statistics, BusinessDays, Dates, CSV, Arrow

##
BusinessDays.initcache(:USNYSE)

##

dates = listbdays(:USNYSE, Date(1980), Date(2020))

mkt_mat = [
  0.000116248  -1.03659e-5   1.1843e-5   -1.04678e-5;
 -1.03659e-5    3.49562e-5  -2.37018e-6   1.9103e-6;
  1.1843e-5    -2.37018e-6   3.77871e-5  -1.05367e-5;
 -1.04678e-5    1.9103e-6   -1.05367e-5   5.98879e-5;
]

df_mkt_benchmark = DataFrame(
        rand(MvNormal(mkt_mat), length(dates))' |> Matrix, [:mkt, :smb, :hml, :umd])
df_mkt_benchmark[!, :date] = dates
select!(df_mkt_benchmark, [:date, :mkt, :smb, :hml, :umd])

##
count = 10000

coef_means = [
        0.0003416652538979692,
        0.7557301603393168,
        0.6035365140401592,
        0.13841355976862427,
        -0.0917780093857371,
    ]

coef_cov = [
    2.52326e-6  -4.71181e-5   3.31182e-5  -7.05013e-5    1.32415e-5;
    -4.71181e-5   0.27552      0.173033     0.0742338    -0.0448335;
     3.31182e-5   0.173033     0.389159     0.0460528    -0.0280572;
    -7.05013e-5   0.0742338    0.0460528    0.40226      -0.000262376;
     1.32415e-5  -0.0448335   -0.0280572   -0.000262376   0.183439;
    ]
d = Gamma(0.5792772418365709, 0.00448646291687775)

df_firm_resp = DataFrame(hcat(
    1:count |> collect,
    rand(MvNormal(coef_means, coef_cov), count)' |> Matrix,
    rand(d, count)
), [:firm_id, :int, :mkt, :smb, :hml, :umd, :var]
)
@transform!(df_firm_resp, :firm_id = Int.(:firm_id))
##

main_mkt_mat = df_mkt_benchmark[:, [:mkt, :smb, :hml, :umd]] |> Matrix
firm_mats = collect.(Tuple.(Tables.rowtable(df_firm_resp[:, [:mkt, :smb, :hml, :umd]])))
firm_errors = rand.(Normal.(0, df_firm_resp.var), nrow(df_mkt_benchmark))

##
df_rand = flatten(DataFrame(
    firm_id = 1:count,
    ret = Ref(main_mkt_mat) .* firm_mats .+ firm_errors
), :ret)

df_rand = @chain df_rand begin
    groupby(:firm_id)
    @transform(:date = df_mkt_benchmark.date)
    select([:firm_id, :date, :ret])
end
##

CSV.write(joinpath("data", "mkt_ret.csv"), df_mkt_benchmark)
CSV.write(joinpath("data", "firm_ret.csv"), df_rand)

##

df_events = DataFrame(
    firm_id = rand(1:count, 1_000_000),
    event_date = rand(Date(1982):Day(1):Date(2019, 12), 1_000_000)
)
df_events = @chain df_events begin
    @rtransform(
        :event_window_start = advancebdays("USNYSE", :event_date, -10),
        :event_window_end = advancebdays("USNYSE", :event_date, 10),
        :est_window_start = :event_date - Year(1),
        :est_window_end = advancebdays("USNYSE", :event_date, -11),
    )
end

CSV.write(joinpath("data", "event_dates.csv"), df_events)


text:
import HypothesisTests
using StatsBase: autocov
using Distributions: Normal

struct JointDistributionTest <: HypothesisTest
    statistic::Float64
end

function JointDistributionTest(
    g, θ::AbstractArray, y::AbstractArray, θ̃::AbstractArray, ỹ::AbstractArray
)
    M₁ = length(θ)
    @assert length(y) == M₁
    M₂ = length(θ̃)
    @assert length(ỹ) == M₂

    g₁ = g.(θ, y)
    g₂ = g.(θ̃, ỹ)
    ḡ₁ = mean(g)
    ḡ₂ = mean(g̃)
    σ̂² = var(g₁; mean = ḡ₁)
    τ̂² = mcvar(g₂)

    statistic = (ḡ₁ - ḡ₂) / (σ̂² / M₁ + τ̂² / M₂)
    return JointDistributionTest(statistic)
end

pvalue(test::JointDistributionTest) = pvalue(Normal(), test.statistic)

function mcvar(x::Vector{T}) where {T<:Real}
    # https://github.com/TuringLang/MCMCChains.jl/blob/9ab389648728064c85c8d331c1efbabc38649107/src/mcse.jl
    n = length(x)
    m = div(n - 2, 2)
    x_ = map(Float64, x)
    ghat = autocov(x_, [0, 1])
    Ghat = sum(ghat)
    value = -ghat[1] + 2 * Ghat
    for i in 1:m
        Ghat = min(Ghat, sum(autocov(x_, [2 * i, 2 * i + 1])))
        Ghat > 0 || break
        value += 2 * Ghat
    end
    return value / n
end


text:
module SDCO

using SparseArrays, DataStructures, LinearAlgebra, Printf

include("types.jl")

include("algebra_E.jl")
include("evaluate.jl")

include("input.jl")
include("io.jl")

include("NT.jl")
include("predcorr.jl")

include("bigM_types.jl")
include("bigM_solver.jl")

include("HSD_solver.jl")

include("test_cases.jl")

end


text:
using CompScienceMeshes
using BEAST

tetrs = CompScienceMeshes.tetmeshsphere(1.0, 0.10)
# tetrs = CompScienceMeshes.tetmeshsphere(1.0, 0.35)
@show numcells(tetrs)

bndry = boundary(tetrs)
edges = skeleton(tetrs, 1)

bndry_edges = [sort(c) for c in cells(skeleton(bndry, 1))]
function is_interior(edge)
    !(sort(edge) in bndry_edges)
end

interior_edges = submesh(is_interior, edges)
@assert numcells(interior_edges) + numcells(skeleton(bndry,1)) == numcells(edges)

X = BEAST.nedelecc3d(tetrs, interior_edges)
@assert numfunctions(X) == numcells(interior_edges)

Id = BEAST.Identity()
Y = curl(X)
A1 = assemble(Id, Y, Y)
A2 = assemble(Id, X, X)
A = A1 - A2

using LinearAlgebra
f = BEAST.ScalarTrace(x -> point(1,0,0) * exp(-norm(x)^2/4))
b = assemble(f, X)

u = A \ b

using Plots

pts = [point(0,t,0) for t in range(-2,2,length=200)];
vals = BEAST.grideval(pts, u, X)
vals2 = BEAST.grideval(pts, u, Y)
vals2 = [point(0,1,0) × val for val in vals2]
plot(getindex.(real(vals),1), m=2)
plot(norm.(vals2), m=2)

# Inspect the value of a single basis function
p = 120
tet = chart(tetrs, cells(tetrs)[p])
nbd = neighborhood(tet, [0.5, 0.5, 0.0])
edg = simplex(tet.vertices[1], tet.vertices[2])
tgt = normalize(edg.vertices[2] - edg.vertices[1])

vals = refspace(X)(nbd)
@assert dot(tgt,vals[1].value) ≈ 1/volume(edg) atol=1e-8

# check also the value of the raviart-thomas 3d in the center
# of their defining face
Z = BEAST.nedelecd3d(tetrs)
@assert numfunctions(Z) == numcells(skeleton(tetrs,2))

fce = simplex(tet.vertices[2], tet.vertices[3], tet.vertices[4])
nbd_rt = neighborhood(tet, [0.0, 1/3, 1/3])
vals_rt = refspace(Z)(nbd_rt)

# carttobary(fce, cartesian(nbd_rt))

#
dot(vals_rt[1].value, normal(fce)) * volume(fce)
o, x, y, z = euclidianbasis(3)
@assert volume(simplex(x,y,z,o)) * 6 ≈ 1

G = boundary(tetrs)
Z = BEAST.ttrace(curl(X), G)

fcr, geo = facecurrents(u, Z)
import Plotly
Plotly.plot(patch(geo, norm.(fcr)))

const CSM = CompScienceMeshes
hemi = submesh(tet -> cartesian(CSM.center(chart(tetrs,tet)))[3] < 0, tetrs)
# error()
bnd_hemi = boundary(hemi)

Xhemi = BEAST.restrict(X, hemi)
tXhemi = BEAST.ttrace(Xhemi, bnd_hemi)

fcr, geo = facecurrents(u, tXhemi)
Plotly.plot(patch(geo, norm.(fcr)))


text:
function extract_params!(pm)
    if haskey(pm, "user_defined_params")
        params = Dict{Symbol,Dict{String,Any}}()
        for key in keys(pm["user_defined_params"])
            params[Symbol(key)] = pm["user_defined_params"][key]
        end
        delete!(pm, "user_defined_params")
    end
    return params
end

function check_powermodels_data!(pm)
    if pm["correct_pm_network_data"]
        _PM.correct_network_data!(pm)
    end
    if haskey(pm, "simplify_net")
        if pm["simplify_net"]
            _PM.simplify_network!(pm)
            _PM.deactivate_isolated_components!(pm)
            _PM.propagate_topology_status!(pm)
        end
    end
    return pm
end

function active_powermodels_silence!(pm)
    if pm["silence"]
        _PM.silence()
    end
end

function check_current_limit!(pm)
    cl = 0
    for (i, branch) in pm["branch"]
        if "c_rating_a" in keys(branch)
            cl += 1
        end
    end
    return cl
end

function set_pq_values_from_timeseries(pm)
    # This function iterates over multinetwork entries and sets p, q values
    # of loads and "sgens" (which are loads with negative P and Q values)
    steps = pm["time_series"]["to_time_step"]-pm["time_series"]["from_time_step"]
    mn = _PM.replicate(pm, steps)

    for (step, network) in mn["nw"]
        step_1=string(parse(Int64,step) - 1)
        load_ts = pm["time_series"]["load"]
        network = delete!(network, "user_defined_params")
        for (idx, load) in network["load"]
            if haskey(load_ts, idx)
                load["pd"] = load_ts[idx]["p_mw"][step_1]
                if haskey(load_ts[idx], "q_mvar")
                    load["qd"] = load_ts[idx]["q_mvar"][step_1]
                end
            end
        end

        gen_ts = pm["time_series"]["gen"]
        for (idx, gen) in network["gen"]
            if haskey(gen_ts, idx)
                gen["pg"] = gen_ts[idx]["p_mw"][step_1]
                if haskey(gen_ts[idx], "max_p_mw")
                    gen["pmax"] = gen_ts[idx]["max_p_mw"][step_1]
                else
                    gen["pmax"] = gen_ts[idx]["p_mw"][step_1]
                end
                if haskey(gen_ts[idx], "min_p_mw")
                    gen["pmin"] = gen_ts[idx]["min_p_mw"][step_1]
                else
                    gen["pmin"] = gen_ts[idx]["p_mw"][step_1]
                end
                if haskey(gen_ts[idx], "max_q_mvar")
                    gen["qmax"] = gen_ts[idx]["max_q_mvar"][step_1]
                end
                if haskey(gen_ts[idx], "min_q_mvar")
                    gen["qmin"] = gen_ts[idx]["min_q_mvar"][step_1]
                end
                if haskey(gen_ts[idx], "q_mvar")
                    gen["qg"] = gen_ts[idx]["q_mvar"][step_1]
                end
            end
        end
    end
    return mn
end


text:
import GPLinearODEMaker.powers_of_negative_one

"""
    cos_kernel(hyperparameters, δ, dorder; shift_ind=0)

Created by kernel_coder(). Requires 1 hyperparameters.
Likely created using cos_kernel_base() as an input.
Use with include("src/kernels/cos_kernel.jl").

# Arguments
- `hyperparameters::Vector`: The hyperparameter values. For this kernel, they should be `["λ"]`
- `δ::Real`: The difference between the inputs (e.g. `t1 - t2`)
- `dorder::Vector{<:Integer}`: How many times to differentiate with respect to the inputs and the `hyperparameters` (e.g. `dorder=[0, 1, 0, 2]` would correspond to differentiating once w.r.t the second input and twice w.r.t `hyperparameters[2]`)
- `shift_ind::Integer=0`: If changed, the index of which hyperparameter is the `δ` shifting one
"""
function cos_kernel(
    hyperparameters::Vector{<:Real},
    δ::Real,
    dorder::Vector{<:Integer};
    shift_ind::Integer=0)

    include_shift = shift_ind!=0
    @assert length(hyperparameters)==1+Int(include_shift) "hyperparameters is the wrong length"
    dorder_len = 3 + Int(include_shift)
    @assert length(dorder)==dorder_len "dorder is the wrong length"
    @assert maximum(dorder) < 3 "No more than two derivatives for each time or hyperparameter can be calculated"

    dorder2 = dorder[2]
    dorder[2] += dorder[1]

    if include_shift
        dorder[2] += dorder[shift_ind+2]
        δ += hyperparameters[shift_ind]
        hyperparameters_view = view(hyperparameters, 1:2 .!= shift_ind)
        dorder_view = view(dorder, 2:dorder_len)
        dorder_view = view(dorder_view, 1:(dorder_len-1) .!= (shift_ind+1))
    else
        hyperparameters_view = hyperparameters
        dorder_view = view(dorder, 2:dorder_len)
    end

    λ = hyperparameters_view[1]

    if dorder_view==[6, 2]
        func = -2688*cos(2*δ*π/λ)*π^6/λ^8 + 1792*δ*sin(2*δ*π/λ)*π^7/λ^9 + 256*δ^2*cos(2*δ*π/λ)*π^8/λ^10
    end

    if dorder_view==[5, 2]
        func = -960*sin(2*δ*π/λ)*π^5/λ^7 - 768*δ*cos(2*δ*π/λ)*π^6/λ^8 + 128*δ^2*sin(2*δ*π/λ)*π^7/λ^9
    end

    if dorder_view==[4, 2]
        func = 320*cos(2*δ*π/λ)*π^4/λ^6 - 320*δ*sin(2*δ*π/λ)*π^5/λ^7 - 64*δ^2*cos(2*δ*π/λ)*π^6/λ^8
    end

    if dorder_view==[3, 2]
        func = 96*sin(2*δ*π/λ)*π^3/λ^5 + 128*δ*cos(2*δ*π/λ)*π^4/λ^6 - 32*δ^2*sin(2*δ*π/λ)*π^5/λ^7
    end

    if dorder_view==[2, 2]
        func = -24*cos(2*δ*π/λ)*π^2/λ^4 + 48*δ*sin(2*δ*π/λ)*π^3/λ^5 + 16*δ^2*cos(2*δ*π/λ)*π^4/λ^6
    end

    if dorder_view==[1, 2]
        func = -4*sin(2*δ*π/λ)*π/λ^3 - 16*δ*cos(2*δ*π/λ)*π^2/λ^4 + 8*δ^2*sin(2*δ*π/λ)*π^3/λ^5
    end

    if dorder_view==[0, 2]
        func = -4*δ*sin(2*δ*π/λ)*π/λ^3 - 4*δ^2*cos(2*δ*π/λ)*π^2/λ^4
    end

    if dorder_view==[6, 1]
        func = 384*cos(2*δ*π/λ)*π^6/λ^7 - 128*δ*sin(2*δ*π/λ)*π^7/λ^8
    end

    if dorder_view==[5, 1]
        func = 160*sin(2*δ*π/λ)*π^5/λ^6 + 64*δ*cos(2*δ*π/λ)*π^6/λ^7
    end

    if dorder_view==[4, 1]
        func = -64*cos(2*δ*π/λ)*π^4/λ^5 + 32*δ*sin(2*δ*π/λ)*π^5/λ^6
    end

    if dorder_view==[3, 1]
        func = -24*sin(2*δ*π/λ)*π^3/λ^4 - 16*δ*cos(2*δ*π/λ)*π^4/λ^5
    end

    if dorder_view==[2, 1]
        func = 8*cos(2*δ*π/λ)*π^2/λ^3 - 8*δ*sin(2*δ*π/λ)*π^3/λ^4
    end

    if dorder_view==[1, 1]
        func = 2*sin(2*δ*π/λ)*π/λ^2 + 4*δ*cos(2*δ*π/λ)*π^2/λ^3
    end

    if dorder_view==[0, 1]
        func = 2*δ*sin(2*δ*π/λ)*π/λ^2
    end

    if dorder_view==[6, 0]
        func = -64*cos(2*δ*π/λ)*π^6/λ^6
    end

    if dorder_view==[5, 0]
        func = -32*sin(2*δ*π/λ)*π^5/λ^5
    end

    if dorder_view==[4, 0]
        func = 16*cos(2*δ*π/λ)*π^4/λ^4
    end

    if dorder_view==[3, 0]
        func = 8*sin(2*δ*π/λ)*π^3/λ^3
    end

    if dorder_view==[2, 0]
        func = -4*cos(2*δ*π/λ)*π^2/λ^2
    end

    if dorder_view==[1, 0]
        func = -2*sin(2*δ*π/λ)*π/λ
    end

    if dorder_view==[0, 0]
        func = cos(2*δ*π/λ)
    end

    dorder[2] = dorder2  # resetting dorder[2]
    return powers_of_negative_one(dorder2) * float(func)  # correcting for amount of t2 derivatives

end

return cos_kernel, 1  # the function handle and the number of kernel hyperparameters

max_stars_repo_name: christiangil/GPLinearODEMaker

max_issues_repo_name: christiangil/GPLinearODEMaker

max_forks_repo_name: christiangil/GPLinearODEMaker


text:
using JuMPTutorials
using Test
using Statistics

cd("introduction")
@testset "Introduction" begin
    include("introduction/an_introduction_to_julia.jl")
    include("introduction/getting_started_with_JuMP.jl")
    @test objective_value(model) ≈ 205
    include("introduction/variables_constraints_objective.jl")
    @test objective_value(vector_model) ≈ 4.9230769230769225
    include("introduction/solvers_and_solutions.jl")
end
cd("..")

cd("using_JuMP")
@testset "Using JuMP" begin
    include("using_JuMP/working_with_data_files.jl")
    @test objective_value(model) == 23
    @test countryindex == [1, 5, 9, 10, 38, 39, 55, 63, 64, 75, 78, 81, 89, 104, 107, 130, 138, 158, 162, 167, 182, 188, 190]
    include("using_JuMP/problem_modification.jl")
    include("using_JuMP/nonlinear_modelling.jl")
end
cd("..")

cd("optimization_concepts")
@testset "Optimization Concepts" begin
    include("optimization_concepts/benders_decomposition.jl")
    @test objective_value(master_problem_model) == -4.0
    @test value.(x) == [0.0, 1.0]
    include("optimization_concepts/integer_programming.jl")
    include("optimization_concepts/conic_programming.jl")
end
cd("..")

cd("modelling")
@testset "Modelling Examples" begin
    include("modelling/n-queens.jl")
    @test sum([sum(value.(x)[i, :]) == 1 for i=1:N]) == N
    @test sum([sum(value.(x)[:, i]) == 1 for i=1:N]) == N
    @test sum([sum(diag(value.(x),i)) <= 1 for i in -(N-1):(N-1)]) == length(-(N-1):(N-1))
    @test sum([sum(diag(reverse(value.(x),dims=1), i)) <=1 for i in -(N-1):(N-1)]) == length(-(N-1):(N-1))
    @test JuMP.termination_status(model) == MOI.OPTIMAL

    include("modelling/sudoku.jl")
    @test sol == [
        5  3  4  6  7  8  9  1  2;
        6  7  2  1  9  5  3  4  8;
        1  9  8  3  4  2  5  6  7;
        8  5  9  7  6  1  4  2  3;
        4  2  6  8  5  3  7  9  1;
        7  1  3  9  2  4  8  5  6;
        9  6  1  5  3  7  2  8  4;
        2  8  7  4  1  9  6  3  5;
        3  4  5  2  8  6  1  7  9]

    include("modelling/problems_on_graphs.jl")
    @test value.(y) == [0.0, 1.0, 1.0, 1.0, 0.0, 0.0]
    @test value.(x) == [0.0, 1.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0]
    @test value.(m) == [
        0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0;
        1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0
    ]
    @test value(sum(z)) ≈ 3
    @test sum(c) ≈ nv(g)

    include("modelling/network_flows.jl")
    @test objective_value(shortest_path) == 55.0
    @test value.(x) == [
        0.0 0.0 1.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 1.0 0.0;
        0.0 1.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0]
    @test objective_value(assignment) ≈ 20.0
    @test value.(y) == [0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0; 1.0 0.0 0.0 0.0; 0.0 0.0 0.0 1.0]
    @test objective_value(max_flow) ≈ 6.0

    include("modelling/finance.jl")
    @test objective_value(financing) ≈ 92.49694915254233
    @test objective_value(auction) ≈ 21.0
    @test value.(y) == [1.0, 1.0, 1.0, 0.0, 0.0, 0.0]
    @test objective_value(portfolio) ≈ 22634.41784988414
    @test value.(x) ≈ [497.0455298498642, 0.0, 502.95448015948074]
    include("modelling/power_systems.jl")
    @test g_opt == [1000.0, 300.0]
    @test w_opt == 200.0
    @test w_f - w_opt == 0
    @test obj == 90000.0
    include("modelling/geometric_problems.jl")
    @test value.(p) ≈ [
        0.44790964261631827    0.0468981793661497;
        -0.03193526635198919   -0.6706136210384356;
        0.404335805799056     -0.45130815913688466;
        -0.39429534726904925   -0.13282535401213758;
        0.02532703978422118    0.4124207687120701;
        -0.0016520566419420052 -0.43954821308159137;
        1.0         1.0;
        1.0        -1.0;
        -1.0        -1.0;
        -1.0         1.0;
        1.0        -0.5;
        -1.0        -0.2;
        -0.2        -1.0;
        0.1         1.0]
    include("modelling/rocket_control.jl")
    @test objective_value(rocket) ≈ 1.0128340648308058
end
cd("..")


text:
using Documenter, CellularAutomata

makedocs(sitename = "CellularAutomata Module Documentation", format = Documenter.HTML(prettyurls = false))


text:
#
# Denoise an SPD Example with Cyclic Proximal Point applied to the
#
# L2-TV functional with anisotropic TV
#
# where the example is the same data as for the corresponding CP algorithm
#
using Manopt, Manifolds
using Images, CSV, DataFrames, LinearAlgebra, JLD2
#
# Settings
ExportResult = true
ExportOrig = true
ExportTable = true
asy_render_detail = 2
#
# Manifold and Data
f = artificial_SPD_image2(32)
pixelM = SymmetricPositiveDefinite(3)
resultsFolder = "examples/Total_Variation/SPD_TV/"
experimentName = "ImageCPPA"
if !isdir(resultsFolder)
    mkdir(resultsFolder)
end
if ExportOrig
    asymptote_export_SPD(
        resultsFolder * experimentName * "-orig.asy";
        data = f,
        scaleAxes = (7.5, 7.5, 7.5),
    )
    render_asymptote(resultsFolder * experimentName * "-orig.asy", render=asy_render_detail)
end
#
# Parameters
α = 6.0
maxIterations = 4000
#
# Build Problem for L2-TV
M = PowerManifold(pixelM, NestedPowerRepresentation(), size(f)...)
d = length(size(f))
rep(d) = (d > 1) ? [ones(Int, d)..., d] : d
fidelity(x) = 1 / 2 * distance(M, x, f)^2
Λ(x) = forward_logs(M, x) # on T_xN
prior(x) = norm(norm.(Ref(pixelM), repeat(x, rep(d)...), Λ(x)), 1)
#
# Setup and Optimize
cost(x) = fidelity(x) + α * prior(x)
proxes = [(λ, x) -> prox_distance(M, λ, f, x, 2), (λ, x) -> prox_TV(M, α * λ, x, 1)]
x0 = f
@time o = cyclic_proximal_point(
    M,
    cost,
    proxes,
    x0;
    debug = [
        :Iteration,
        " | ",
        DebugProximalParameter(),
        " | ",
        :Change,
        " | ",
        :Cost,
        "\n",
        100,
        :Stop,
    ],
    record = [:Iteration, :Iterate, :Cost],
    stopping_criterion = StopAfterIteration(maxIterations),
    return_options = true,
)
y = get_solver_result(o)
yRec = get_record(o)
#
# Results
if ExportResult
    asymptote_export_SPD(
        resultsFolder *
        experimentName *
        "-result-$(maxIterations)-α$(replace(string(α), "." => "-")).asy";
        data = y,
        scaleAxes = (7.5, 7.5, 7.5),
    )
    render_asymptote(
        resultsFolder *
        experimentName *
        "-result-$(maxIterations)-α$(replace(string(α), "." => "-")).asy";
        render=asy_render_detail,
    )
end
if ExportTable
    A = cat([y[1] for y in yRec], [y[3] for y in yRec]; dims = 2)
    CSV.write(
        string(resultsFolder * experimentName * "ResultCost.csv"),
        DataFrame(A),
        writeheader = false,
    )
    save(
        resultsFolder * experimentName * "-CostValue.jld2",
        Dict("compareCostFunctionValue" => last(yRec)[3]),
    )
end


text:
export result_beam_force,result_beam_displacement

"""
    result_beam_force(assembly,lc_id,beam_id)
读取已求解的assembly实例的单个beam单元内力
# Arguments
- `assembly`: Assembly类型实例
- `lc_id`: 工况id
- `beam_id`: 单元id
# 返回
- 包含梁单元坐标系下的12个节点力向量
"""
function result_beam_force(assembly,lc_id,beam_id)
    path=assembly.working_path
    beam_id=string(beam_id)
    if !(id in keys(assembly.structure.beams))
        throw("beam id "*string(beam_id)*" doesn't exists!")
    end
    lc_id=string(lc_id)
    if !(lc_id in keys(assembly.lcset.cases))
        throw("loadcase with id "*lc_id*" doesn't exist!")
    end
    elm=assembly.structure.beams[beam_id]
    i=elm.node1.hid
    j=elm.node2.hid
    iDOFs=6*i-5:6*i
    jDOFs=6*j-5:6*j
    u=read_vector(path*"/.analysis",lc_id*"_u.v")
    T=elm.T
    uᵉ=T*[d[iDOFs];u[jDOFs]]

    rDOF=findall(x->x==true,elm.release)
    Kᵉ=integrateK(elm)
    K̃ᵉ,P̃ᵉ=static_condensation(K̃ᵉ,zeros(12),rDOF)
    return K̃ᵉ*uᵉ
end

"""
    result_beam_displacement(assembly,lc_id,beam_id)
读取已求解的assembly实例的单个beam单元两端节点位移
# 参数
- `assembly`: Assembly类型实例
- `lc_id`: 工况id
- `beam_id`: 单元id
# 返回
- 包含梁单元坐标系下的12个节点位移向量
"""
function result_beam_displacement(assembly,lc_id,beam_id)
    path=assembly.working_path
    beam_id=string(beam_id)
    if !(id in keys(assembl.structure.beams))
        throw("beam id "*string(beam_id)*" doesn't exists!")
    end
    lc_id=string(lc_id)
    if !(lc_id in keys(assembly.lcset.cases))
        throw("loadcase with id "*lc_id*" doesn't exist!")
    end
    loadcase=assembly.lcset.cases[lc_id]
    u=read_vector(path*"/.analysis",lc_id*"_u.v")
    elm=assembly.structure.beams[beam_id]
    i=elm.node1.hid
    j=elm.node2.hid
    iDOFs=6*i-5:6*i
    jDOFs=6*j-5:6*j
    if beam_id in keys(loadcase.beam_forces)
        Pᵉ=zeros(12)
    else
        force=loadcase.beam_forces[beam_id]
        Pᵉ=FEBeam.integrateP!(elm,force)
    end
    T=elm.T
    uᵉ=T*[u[iDOFs];u[jDOFs]]
    K̄ᵉ=integrateK(elm)
    i=.!elm.release
    j=elm.release
    Kᵢᵢ=K̄ᵉ[i,i]
    Kᵢⱼ=K̄ᵉ[i,j]
    Kⱼᵢ=K̄ᵉ[j,i]
    Kⱼⱼ=K̄ᵉ[j,j]
    Pⱼ=Pᵉ[j]
    ūᵉ=zeros(12)
    ūᵉ[i]=uᵉ[i]
    ūᵉ[j]=inv(Kⱼⱼ)*(Pⱼ-Kⱼᵢ*uᵉ[i])
    return ūᵉ
end

"""
    result_beam_displacement(assembly,lc_id)
读取已求解的assembly实例的所有beam单元两端节点位移
# 参数
- `assembly`: Assembly类型实例
- `lc_id`: 工况id
# 返回
- 所有梁在单元坐标系下的12个节点位移向量
"""
function result_beam_displacement(assembly,lc_id)
    path=assembly.working_path
    lc_id=string(lc_id)
    if !(lc_id in keys(assembly.lcset.cases))
        throw("loadcase with id "*lc_id*" doesn't exist!")
    end
    loadcase=assembly.lcset.cases[lc_id]
    u=read_vector(path*"/.analysis",lc_id*"_u.v")
    res=Dict{String,Array{Float64}}()
    for elm in values(assembly.structure.beams)
        i=elm.node1.hid
        j=elm.node2.hid
        iDOFs=6*i-5:6*i
        jDOFs=6*j-5:6*j
        if !(elm.id in keys(loadcase.beam_forces))
            Pᵉ=zeros(12)
        else
            force=loadcase.beam_forces[elm.id]
            Pᵉ=BeamModule.integrateP!(elm,force)
        end
        T=elm.T
        uᵉ=T*[u[iDOFs];u[jDOFs]]
        K̄ᵉ=integrateK(elm)
        i=.!elm.release
        j=elm.release
        Kᵢᵢ=K̄ᵉ[i,i]
        Kᵢⱼ=K̄ᵉ[i,j]
        Kⱼᵢ=K̄ᵉ[j,i]
        Kⱼⱼ=K̄ᵉ[j,j]
        Pⱼ=Pᵉ[j]
        ūᵉ=zeros(12)
        ūᵉ[i]=uᵉ[i]
        ūᵉ[j]=inv(Kⱼⱼ)*(Pⱼ-Kⱼᵢ*uᵉ[i])
        res[elm.id]=ūᵉ
    end
    return res
end

function result_element()
    Δdᵉ=(elm.T)'*[Δd[6*hid1-5:6*hid1];Δd[6*hid2-5:6*hid2]]
    Δϵᵉ=Bᵉ*Δdᵉ
    Δσᵉ=∫D*Δϵᵉ
    σᵉ+=Δσᵉ
end


text:
StandaloneInteriorEigenfunction(
    domain::AbstractPlanarDomain,
    orientation = fmpq(0);
    stride::Integer = 1,
    offset::Integer = 0,
    even = false,
    odd = false,
) = StandaloneInteriorEigenfunction(
    center(domain),
    orientation,
    domain.parent;
    stride,
    offset,
    even,
    odd,
)

function Base.show(io::IO, u::StandaloneInteriorEigenfunction)
    println(
        io,
        "Standalone interior eigenfunction" *
        ifelse(u.even, " - even", "") *
        ifelse(u.odd, " - odd", ""),
    )
    if !haskey(io, :compact) || !io[:compact]
        println(io, "vertex: $(u.vertex)")
        println(io, "orientation: $(u.orientation)")
        print(io, "number of set coefficients: $(length(u.coefficients))")
    end
end

function set_eigenfunction!(u::StandaloneInteriorEigenfunction, coefficients::Vector)
    resize!(u.coefficients, length(coefficients))
    copy!(u.coefficients, u.parent.(coefficients))
end

"""
    coordinate_transformation(u::StandaloneInteriorEigenfunction, xy::AbstractVector)

Takes a 2-element vector `xy` representing a point in the plane in
Cartesian coordinates and makes a (affine) change of coordinates so
that `u.vertex` is put at the origin.
"""
function coordinate_transformation(
    u::StandaloneInteriorEigenfunction{T},
    xy::AbstractVector,
) where {T}
    iszero(u.orientation) && return xy .- u.vertex

    if T == fmpq
        s, c = sincospi(u.orientation, u.parent)
    elseif T == arb
        s, c = sincos(u.orientation)
    end
    M = SMatrix{2,2}(c, s, -s, c)
    return M * (xy .- u.vertex)
end

function (u::StandaloneInteriorEigenfunction)(
    xy::AbstractVector{T},
    λ::arb,
    ks::UnitRange{Int};
    boundary = nothing,
    notransform::Bool = false,
) where {T<:Union{arb,arb_series}}
    if !notransform
        xy = coordinate_transformation(u, xy)
    end

    r, θ = polar_from_cartesian(xy)
    if u.even
        νs = u.stride * (ks.start-1:ks.stop-1) .+ u.offset
    elseif u.odd
        νs = u.stride * (ks.start-1:ks.stop-1) .+ 1 .+ u.offset
    else
        νs = u.stride * (div(ks.start, 2):div(ks.stop, 2)) .+ u.offset
    end
    bessel_js = let sqrtλr = sqrt(λ) * r
        Dict(ν => bessel_j(u.parent(ν), sqrtλr) for ν in νs)
    end

    res = similar(ks, T)
    for i in eachindex(ks)
        k = ks[i]
        if u.even
            ν = u.stride * (k - 1) + u.offset

            res[i] = bessel_js[ν] * cos(ν * θ)
        elseif u.odd
            ν = u.stride * (k - 1) + 1 + u.offset

            res[i] = bessel_js[ν] * sin(ν * θ)
        else
            ν = u.stride * div(k, 2) + u.offset

            if isodd(k)
                res[i] = bessel_js[ν] * cos(ν * θ)
            else
                res[i] = bessel_js[ν] * sin(ν * θ)
            end
        end
    end

    return res
end

# TODO: Figure out how to handle this
set_domain!(u::StandaloneInteriorEigenfunction, ::AbstractDomain) = u

# TODO: Figure out how to handle this
recompute!(u::StandaloneInteriorEigenfunction) = u


text:
module ChirpsValidator

using SearchLight, Validation

function not_empty(field::Symbol, m::T, args::Vararg{Any})::ValidationResult where {T<:AbstractModel}
  isempty(getfield(m, field)) && return ValidationResult(invalid, :not_empty, "should not be empty")

  ValidationResult(valid)
end

function minimum_length(field::Symbol, m::T, args::Vararg{Any})::ValidationResult where {T<:AbstractModel}
  length(getfield(m, field)) < 20 && return ValidationResult(invalid, :minimum_length, "should be at least 20 letters long")

  ValidationResult(valid)
end

end

max_stars_repo_name: essenciary/chirper

max_issues_repo_name: essenciary/chirper

max_forks_repo_name: essenciary/chirper


text:
#Copyright (c) 2015 Rex Computing and Isaac Yonemoto

#see LICENSE.txt

#this work was supported in part by DARPA Contract D15PC00135

#hilbert.jl
#tests inverting a hilbert matrix.  All of the entries *should* be integers.

vtypes = [Float16, Float32, Float64, BigFloat]

#generate a hilbert matrix
function hilbert(T::Type, n)
  [1 / (convert(T, x + y - 1)) for x=1:n, y=1:n]
end

for T in vtypes
  for i in 1:10
    h = hilbert(T, i)

    invh = h^(-1)

    println("hilbert matrix $(i) upper corner in $(T) should be $(i ^ 2) but is: $(invh[1,1])")
  end
end

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl