text:
using Raylib

function main()

    screenWidth = 800
    screenHeight = 450

    Raylib.InitWindow(screenWidth, screenHeight, "raylib [core] example - gamepad input")

    Raylib.SetTargetFPS(60)

    while !Raylib.WindowShouldClose()

        Raylib.BeginDrawing()

        Raylib.ClearBackground(Raylib.RAYWHITE)

        if Raylib.IsGamepadAvailable(0)

            gamepad_name = unsafe_string(Raylib.GetGamepadName(0))
            Raylib.DrawText("GP1: $(gamepad_name)", 10, 10, 10, Raylib.BLACK)

            Raylib.DrawText("- GENERIC GAMEPAD -", 280, 180, 20, Raylib.GRAY)

            Raylib.DrawText("DETECTED AXIS [$(Raylib.GetGamepadAxisCount(0))]:", 10, 50, 10, Raylib.MAROON)

            for i in 0:Raylib.GetGamepadAxisCount(0)-1
                Raylib.DrawText("AXIS $i: $(Raylib.GetGamepadAxisMovement(0, i))", 20, 70 + 20*i, 10, Raylib.DARKGRAY)
            end

            if Raylib.GetGamepadButtonPressed() != -1
                Raylib.DrawText("DETECTED BUTTON: $(Raylib.GetGamepadButtonPressed())", 10, 430, 10, Raylib.RED)
            else
                Raylib.DrawText("DETECTED BUTTON: NONE", 10, 430, 10, Raylib.GRAY)
            end

        else
            Raylib.DrawText("GP1: NOT DETECTED", 10, 10, 10, Raylib.GRAY)
        end

        Raylib.EndDrawing()

    end

    Raylib.CloseWindow()

end


text:
using PyPlot
using Animations

function easingplots()

    efuncs = [noease, stepease, linear, sineio, saccadic, expin, expout, polyin, polyout]
    params = [[nothing], [0.25, 0.5, 0.75], [nothing], [nothing], [1, 2, 4], [0.2, 2, 8], [0.2, 2, 8], [2, 3, 6], [2, 3, 6]]

    n = length(efuncs)
    cols = 3
    rows = Int(ceil(n / cols))
    fig, axes = subplots(rows, cols, figsize=(5, rows/3*4), constrained_layout=true)

    xx = 0:0.005:1

    for (i, (ef, ps)) in enumerate(zip(efuncs, params))

        a = axes[(divrem(i - 1, cols) .+ (1, 1))...]
        a.axis("off")

        a.plot([-0.25, 0], [0, 0], color=(0.5, 0.5, 0.5), linestyle="dotted")
        a.plot([1, 1.25], [1, 1], color=(0.5, 0.5, 0.5), linestyle="dotted")
        a.set_title("$(Symbol(ef))", fontweight=600)
        for p in ps
            easing = isnothing(p) ? ef() : ef(p)
            anim = Animation([0, 1], [0.0, 1.0], easing)
            yy = at.(anim, xx)
            if isnothing(p)
                a.plot(xx, yy)
            else
                a.plot(xx, yy, label="$p")
            end
        end
        if ps != [nothing]
            a.legend(loc="center", bbox_to_anchor=(1, 0.5), frameon=false, borderpad=0, borderaxespad=0, handlelength=0.5)
        end
    end

    fig.savefig(joinpath(dirname(pathof(Animations)), "../misc/easingplots.svg"))

end
easingplots()

ext: jl


text:
"""
    function kfold(f, df, k, r = 1, shuffle=true; kwargs...)

    Provide a simple `k` fold(s) cross-validation, repeated `r` time(s).
    `kwargs` arguments are passed to the `regress` function call.
    This feature overlap in part with the PRESS statistics.
"""
function kfold(f, df, k, r = 1, shuffle=true; kwargs...)
    totalrows = nrow(df)
    gindexes = Vector(undef, k)
    resvec = Vector(undef, k*r)
    sdf = @view df[!, :]

    for cr in 1:r
        if shuffle
            sdf = @view df[shuffle(axes(df, 1)), :]
        end
        for ck in 1:k
            gindexes[ck] = ck:k:totalrows
        end
        for ck in 1:k
            training_range = reduce(vcat, map(x->gindexes[x], filter(x->x!=ck, 1:k)))
            training = @view sdf[training_range, :]
            testing = @view sdf[gindexes[ck], :]

            lr = regress(f, training; kwargs...)
            cres = predict_in_sample(lr, testing, req_stats=[:residuals])
            t_mse= mean(cres.residuals .^2)
            resvec[ (k * (cr -1)) + ck] = (R2 = lr.R2, ADJR2 = lr.ADJR2, TRAIN_MSE = lr.MSE,
                        TRAIN_RMSE = lr.RMSE,
                        TEST_MSE= t_mse , TEST_RMSE= √t_mse )
        end
    end
    return DataFrame(resvec)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "lm_Sensors"
version = v"3.5.0"

# Collection of sources required to build lm_sensors
sources = [
    "https://github.com/lm-sensors/lm-sensors/archive/V$(version.major)-$(version.minor)-$(version.patch).tar.gz" =>
    "f671c1d63a4cd8581b3a4a775fd7864a740b15ad046fe92038bcff5c5134d7e0",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/lm-sensors-*/
make -j${nproc} PREFIX=${prefix}
make install PREFIX=${prefix}
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [p for p in supported_platforms() if p isa Linux]

# The products that we will ensure are always built
products = [
    LibraryProduct("libsensors", :libsensors),
    ExecutableProduct("sensors", :sensors),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"8")


text:
# definitions of natural parameter messages from factors to nodes

@defnaturals LogNormalFactor x Normal begin
    Eμ, Eτ = E(μ), E(τ)
    (Eμ * Eτ, -Eτ/2)
end

@defnaturals LogNormalFactor μ Normal begin
    Ex, Eτ = E(x), E(τ)
    (Ex * Eτ, -Eτ/2)
end

@defnaturals LogNormalFactor τ Gamma begin
    v = V(x) + V(μ) + (E(x) - E(μ))^2
    (1/2, -v/2)
end

@defnaturals LogGammaFactor x Gamma begin
    (E(α) - 1, -1/E(θ))
end

@defnaturals LogGammaFactor x Gamma begin
    (E(α) - 1, -1/E(θ))
end

@defnaturals LogGammaCanonFactor x Gamma begin
    (E(α) - 1, -E(β))
end

nats_mvn{T <: Number}(μ::Vector{T}, Λ::Matrix{T}, x::MvNormalCanon) =
    (Λ * μ, -Λ/2)
nats_mvn{T <: Number}(μ::Number, Λ::Matrix{T}, x::MvNormalCanon) =
    (μ * sum(Λ, 2), -Λ/2)
nats_mvn{T <: Number}(μ::Vector{T}, τ::Vector{T}, x::MvNormalCanon) =
    (τ .* μ, -diagm(τ)/2)
nats_mvn{T <: Number}(μ::Number, τ::Vector{T}, x::MvNormalCanon) =
    (τ * μ, -diagm(τ)/2)
nats_mvn{T <: Number}(μ::Vector{T}, τ::Number, x::MvNormalCanon) =
    (τ * μ, -τ * eye(length(μ))/2)
nats_mvn(μ::Number, τ::Number, x::MvNormalCanon) =
    (μ * τ * ones(length(x)), -τ * eye(length(x))/2)
nats_mvn{T <: Number}(μ::Vector{T}, Λ::Matrix{T}, x::Normal) =
    (sum(Λ * μ), -sum(Λ)/2)
nats_mvn{T <: Number}(μ::Number, Λ::Matrix{T}, x::Normal) =
    (μ * sum(Λ), -sum(Λ)/2)
nats_mvn{T <: Number, D <: Normal}(μ::Vector{T}, τ::Vector{T}, x::Vector{D}) =
    [nats_mvn(m, t) for (m, t) in zip(μ, τ)]
nats_mvn{T <: Number, D <: Normal}(μ::Vector{T}, τ::Number, x::Vector{D}) =
    [nats_mvn(m, τ) for m in μ]
nats_mvn{T <: Number, D <: Normal}(μ::Number, τ::Vector{T}, x::Vector{D}) =
    [nats_mvn(μ, t) for t in τ]
nats_mvn(μ, τ, x::Normal) =
    reduce(add_nats, nats_mvn(μ, τ, [x]))
nats_mvn(μ::Number, τ::Number) = (μ * τ, -τ/2)
nats_mvn{D <: Gamma}(v::Vector, x::Vector{D}) =
    Tuple{Float64, Float64}[(1/2, vv/2) for vv in v]
nats_mvn(v::Vector, x::Gamma) =
    reduce(add_nats, nats_mvn(v, [x]))

@defnaturals LogMvNormalCanonFactor x MvNormalCanon begin
    nats_mvn(E(μ), E(Λ), x)
end

@defnaturals LogMvNormalCanonFactor μ MvNormalCanon begin
    nats_mvn(E(x), E(Λ), μ)
end

@defnaturals LogMvNormalCanonFactor x Normal begin
    nats_mvn(E(μ), E(Λ), x)
end

@defnaturals LogMvNormalCanonFactor μ Normal begin
    nats_mvn(E(x), E(Λ), μ)
end

@defnaturals LogMvNormalCanonFactor Λ Wishart begin
    δ = E(x) - E(μ)
    v = C(x) .+ C(μ) .+ δ * δ'
    (1/2, -v/2)
end

@defnaturals LogMvNormalDiagCanonFactor x Normal begin
    nats_mvn(E(μ), E(τ), x)
end

@defnaturals LogMvNormalDiagCanonFactor x MvNormalCanon begin
    nats_mvn(E(μ), E(τ), x)
end

@defnaturals LogMvNormalDiagCanonFactor μ Normal begin
    nats_mvn(E(x), E(τ), μ)
end

@defnaturals LogMvNormalDiagCanonFactor μ MvNormalCanon begin
    nats_mvn(E(x), E(τ), μ)
end

@defnaturals LogMvNormalDiagCanonFactor τ Gamma begin
    δ = E(x) - E(μ)
    v = V(x) + V(μ) + δ.^2
    nats_mvn(v, τ)
end

@defnaturals LogWishartFactor X Wishart begin
    n = E(ν)
    Vinv = Einv(V)
    p = size(Vinv, 1)
    ((n - p - 1)/2, -Vinv/2)
end

@defnaturals LogMarkovChainFactor π0 Dirichlet begin
    (E(z)[:, 1], )
end

@defnaturals LogMarkovChainFactor A MarkovMatrix begin
    (slice(sum(C(z), 3), :, :), )
end

@defnaturals LogMarkovChainFactor z HMM begin
    (zero_like(z.ψ), Elog(π0), Elog(A))
end

@defnaturals LogDirichletFactor x Dirichlet begin
    (E(α) - 1,)
end

@defnaturals LogMarkovMatrixFactor x MarkovMatrix begin
    (E(A) - 1,)
end


text:
# The bulk of LLVM's object model consists of values, which comprise a very rich type
# hierarchy.

export Value

# subtypes are expected to have a 'ref::API.LLVMValueRef' field
abstract type Value end

Base.unsafe_convert(::Type{API.LLVMValueRef}, val::Value) = val.ref

const value_kinds = Vector{Type}(fill(Nothing, typemax(API.LLVMValueKind)+1))
function identify(::Type{Value}, ref::API.LLVMValueRef)
    kind = API.LLVMGetValueKind(ref)
    typ = @inbounds value_kinds[kind+1]
    typ === Nothing && error("Unknown value kind $kind")
    return typ
end
function register(T::Type{<:Value}, kind::API.LLVMValueKind)
    value_kinds[kind+1] = T
end

function refcheck(::Type{T}, ref::API.LLVMValueRef) where T<:Value
    ref==C_NULL && throw(UndefRefError())
    if Base.JLOptions().debug_level >= 2
        T′ = identify(Value, ref)
        if T != T′
            error("invalid conversion of $T′ value reference to $T")
        end
    end
end

# Construct a concretely typed value object from an abstract value ref
function Value(ref::API.LLVMValueRef)
    ref == C_NULL && throw(UndefRefError())
    T = identify(Value, ref)
    return T(ref)
end

## general APIs

export llvmtype, name, name!, replace_uses!, isconstant, isundef, context

llvmtype(val::Value) = LLVMType(API.LLVMTypeOf(val))

name(val::Value) = unsafe_string(API.LLVMGetValueName(val))
name!(val::Value, name::String) = API.LLVMSetValueName(val, name)

function Base.show(io::IO, val::Value)
    output = unsafe_message(API.LLVMPrintValueToString(val))
    print(io, output)
end

replace_uses!(old::Value, new::Value) = API.LLVMReplaceAllUsesWith(old, new)

isconstant(val::Value) = convert(Core.Bool, API.LLVMIsConstant(val))

isundef(val::Value) = convert(Core.Bool, API.LLVMIsUndef(val))

context(val::Value) = Context(API.LLVMGetValueContext(val))

## user values

include("value/user.jl")

## constants

include("value/constant.jl")

## usage

export Use, user, value

@checked struct Use
    ref::API.LLVMUseRef
end

Base.unsafe_convert(::Type{API.LLVMUseRef}, use::Use) = use.ref

user(use::Use) =  Value(API.LLVMGetUser(     use))
value(use::Use) = Value(API.LLVMGetUsedValue(use))

# use iteration

export uses

struct ValueUseSet
    val::Value
end

uses(val::Value) = ValueUseSet(val)

Base.eltype(::ValueUseSet) = Use

function Base.iterate(iter::ValueUseSet, state=API.LLVMGetFirstUse(iter.val))
    state == C_NULL ? nothing : (Use(state), API.LLVMGetNextUse(state))
end

Base.IteratorSize(::ValueUseSet) = Base.SizeUnknown()


text:
#=
  Steiner triplets in ConstraintSolver.jl

  http://www.probp.com/examples/clpset/steiner.pl
  """
  The ternary Steiner problem of order n is to find n(n-1)/6 sets of elements
  in {1,2,...,n} such that each set contains three elements and any two
  sets have at most one element in common.

  For example, the following shows a solution for size n=7:

      {1,2,3}, {1,4,5}, {1,6,7}, {2,4,6}, {2,5,7}, {3,4,7}, {3,5,6}

  Problem taken from:
  C. Gervet: Interval Propagation to Reason about Sets: Definition and
             Implementation of a PracticalLanguage,
             Constraints, An International Journal, vol.1, pp.191-246, 1997.
  """

  Note: This model uses triplets of integers (in comparison to
        binary representation in steiner.jl) and is faster!

  This model manage to solve n=15 (timeout 16s)
n:3 nb:1
[1 2 1]
  0.159783 seconds (428.18 k allocations: 25.673 MiB, 99.24% compilation time)

n:7 nb:7
[1 4 6; 1 3 5; 3 4 7; 1 2 7; 2 5 6; 4 5 7; 4 2 2]
  0.067136 seconds (503.62 k allocations: 26.409 MiB)

n:9 nb:12
[1 4 6; 4 5 7; 5 6 8; 3 6 8; 4 8 9; 1 6 9; 3 4 9; 6 7 9; 1 7 8; 3 7 9; 1 3 5; 2 2 2]
  2.441683 seconds (11.63 M allocations: 495.798 MiB, 44.58% gc time)

n:13 nb:26
[1 3 11; 8 10 12; 1 2 5; 2 10 11; 5 7 10; 2 7 9; 3 6 13; 1 8 9; 3 5 12; 2 4 13; 7 10 12; 8 11 13; 9 10 13; 3 9 10; 7 11 13; 3 7 8; 1 6 7; 5 6 8; 4 11 12; 6 10 11; 6 9 12; 5 9 11; 1 12 13; 2 8 12; 2 3 6; 4 1 4]
  4.297591 seconds (17.44 M allocations: 869.541 MiB, 33.34% gc time)

n:15 nb:35
[1 2 8; 7 10 15; 5 8 15; 3 6 9; 2 3 5; 3 4 8; 1 6 7; 1 12 15; 6 7 14; 1 13 14; 6 11 15; 4 9 15; 3 13 15; 2 14 15; 8 11 12; 8 10 14; 1 4 5; 1 3 10; 2 12 13; 7 8 13; 2 10 11; 6 10 12; 5 11 14; 5 10 13; 2 7 9; 5 7 12; 1 11 12; 5 6 8; 4 12 14; 4 11 13; 4 7 10; 3 12 14; 3 7 11; 2 4 6; 9 1 9]
  6.824711 seconds (29.03 M allocations: 1.556 GiB, 21.72% gc time)

n:19 nb:57
status:TIME_LIMIT
 16.978271 seconds (53.21 M allocations: 3.464 GiB, 25.02% gc time)

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#
using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function steiner(n=7,print_solutions=true,all_solutions=true,timeout=6)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            # "traverse_strategy"=>:BFS,
                                                            "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            "branch_split"=>:Biggest,
                                                            # "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            # "branch_strategy" => :IMPS, # default
                                                            "branch_strategy" => :ABS, # Activity Based Search
                                                            "activity.decay" => 0.999, # default 0.999
                                                            "activity.max_probes" => 1, # default, 10
                                                            "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>timeout,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))
    if !(n % 6 == 1 || n % 6 == 3)
        println("n must be (1|3) modulo 6")
        return
    end

    nb = round(Int,(n * (n-1)) / 6) # number of sets
    println("n:$n nb:$nb")

    @variable(model, 1 <= x[1:nb,1:3] <= n, Int)
    @constraint(model, x[1,1] == 1) # symmetry breaking

    # atmost 1 element in common
    for i in 1:nb, j in i+1:nb
        increasing_strict(model, x[i,:])
        b = @variable(model, [1:3,1:3], Bin)
        for p in 1:3, q in 1:3
            @constraint(model, b[p,q] := {x[i,p] == x[j,q]} )
        end
        @constraint(model, sum(b[:]) <= 1)
    end

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    num_sols = 0
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        # println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                println(x_val)
            end
        end
    else
        println("status:$status")
    end

    return status, num_sols
end

# @time steiner(7,true,false)
for n in [i for i in 3:20 if (i % 6 == 1 || i % 6 == 3)]
    @time status, num_sols = steiner(n,true,false,16)
    println()
    if status == MOI.TIME_LIMIT
        break
    end
end

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
module DistributedBoosting
using HDF5

type Unibeta
	unibeta::Array{Float64,1}
	n::Array{Int,1}
	labels::Array{ASCIIString,1}
end

type Covarmat
	covarmat::Array{Float64,2}
	labels::Array{ASCIIString,1}
end

abstract Abstractdiscohort

type Discohort <: Abstractdiscohort
	x::Array{Float64,2}
	y::Array{Float64,1}
	labels::Array{ASCIIString,1}
end

type Remotecohort <: Abstractdiscohort #hier nochmal schauen
	filename::ASCIIString
	unibeta::Unibeta
	covarmat::Covarmat
end

type Boostscratch
	actualstepno::Int
	actualbeta::Array{Float64,1}
	actualscore::Array{Float64,1}
	pooledunibeta::Unibeta
	pooledcovarmat::Covarmat
	actualupdates::Array{Float64,1}
	selections::Array{ASCIIString,1}
	wantedlabels::Array{ASCIIString,1}
	nu::Float64
	stepno::Int
	actualsel::Int
	actualnom::Array{Float64,1}
	a::Int
	x::Int
end

include("boost.jl")
include("boosting.jl")
include("calc.jl")
include("fetch_pool.jl")
include("getselections.jl")
include("pooling.jl")
include("readwrite.jl")

end

max_stars_repo_name: danielazoeller/DistributedBoosting

max_issues_repo_name: danielazoeller/DistributedBoosting

max_forks_repo_name: danielazoeller/DistributedBoosting


text:
# [MLJ.jl]: https://github.com/alan-turing-institute/MLJ.jl
# [RDatasets.jl]: https://github.com/JuliaStats/RDatasets.jl
# [NearestNeighbors.jl]: https://github.com/KristofferC/NearestNeighbors.jl
#
# ## Tuning a single hyperparameter
#
# In MLJ, tuning is implemented as a model wrapper.
# After wrapping a model in a _tuning strategy_ (e.g. cross-validation) and binding the wrapped model to data in a _machine_, fitting the machine initiates a search for optimal model hyperparameters.
#
# Let's use a decision tree classifier and tune the maximum depth of the tree.
# As usual, start by loading data and the model

using MLJ
MLJ.color_off() #src
X, y = @load_iris
@load DecisionTreeClassifier

# ### Specifying a range of value
#
# To specify a range of value, you can use the `range` function:

dtc = DecisionTreeClassifier()
r   = range(dtc, :max_depth, lower=1, upper=5)

# As you can see, the range function takes a model (`dtc`), a symbol for the hyperparameter of interest (`:max_depth`) and indication of how to samples values.
# You can either manually specify an iterator over values by using the `values=` keyword or, where appropriate, use the `lower` and `upper` keywords which help define such an iterator between bounds.
#
# Once a range of value has been defined, you can then wrap the model in a `TunedModel` specifying the tuning strategy:

tm = TunedModel(model=dtc, ranges=[r, ], measure=cross_entropy)

# ### Fitting and inspecting a tuned model
#
# To fit a tuned model, you can use the usual syntax:

m = machine(tm, X, y)
fit!(m)

# In order to inspect the best model, you can use the function `fitted_params` on the machine and inspect the `best_model` field:

fitted_params(m).best_model.max_depth

# ## Tuning nested hyperparameters

# Let's generate simple dummy regression data

X = (x1=rand(100), x2=rand(100), x3=rand(100))
y = 2X.x1 - X.x2 + 0.05 * randn(100);

# Let's then build a simple ensemble model with decision tree regressors:

dtr = @load DecisionTreeRegressor
forest = EnsembleModel(atom=dtr)

# Such a model has *nested* hyperparameters in that the ensemble has hyperparameters (e.g. the `:bagging_fraction`) and the atom has hyperparameters (e.g. `:n_subfeatures` or `:max_depth`).
#
# Range for nested hyperparameters are specified using dot syntax, the rest is done in much the same way as before:

r1 = range(forest, :(atom.n_subfeatures), lower=1, upper=3)
r2 = range(forest, :bagging_fraction, lower=0.4, upper=1.0)
tm = TunedModel(model=forest, tuning=Grid(resolution=12),
                resampling=CV(nfolds=6), ranges=[r1, r2],
                measure=rms)
m = machine(tm, X, y)
fit!(m);

# A useful function to inspect a model after fitting it is the `report` function which collects information on the model and the tuning, for instance you can use it to recover the best measurement:

r = report(m)
r.best_measurement

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials


text:
struct BecamePinger
    depth::Int
end

struct BecamePonger
    depth::Int
end

mutable struct TestChecker
    results
    TestChecker() = new([])
end

struct BecamePing end
struct BecamePong end

struct TestFn
    fn::Function
    results_to::Addr
end

@ctx function (me::TestChecker)(msg)
    push!(me.results, msg)
end

@ctx function (me::BecamePinger)(msg::BecamePong)
    depth = me.depth
    if depth > 0
        become(BecamePonger(depth))
        send(self(), BecamePing())
    end
    return nothing
end

@ctx function (me::BecamePonger)(msg::BecamePing)
    depth = me.depth
    if depth > 0
        become(BecamePinger(depth - 1))
        send(self(), BecamePong())
    end
    return nothing
end

@ctx function (me::BecamePinger)(msg::TestFn)
    send(msg.results_to, msg.fn(me))
end

@ctx function (me::BecamePonger)(msg::TestFn)
    send(msg.results_to, msg.fn(me))
end

const PING_COUNT = 50_000

function run_becomepingpong_test(lib::ActorLib)
    return @testset "ActorInterfaceTests/became_pingpong: changing type with become" begin
        @test ex_actorcount(lib) == 0
        becamer = BecamePinger(PING_COUNT)
        becameraddr = ex_spawn!(lib, becamer)
        checker = TestChecker()
        checkeraddr = ex_spawn!(lib, checker)
        @test ex_actorcount(lib) == 2
        ex_send!(lib, becameraddr, TestFn(me -> me.depth == 0, checkeraddr))
        ex_send!(lib, becameraddr, BecamePong())
        println("Becoming $(PING_COUNT * 2) times to a different type and delivering the same amount of messages")
        ex_runtofinish(lib)
        ex_send!(lib, becameraddr, TestFn(me -> me.depth == 0, checkeraddr))
        ex_runtofinish(lib)
        @test checker.results == [false, true]
    end
end


text:
using Luxor, IntervalLinearAlgebra, LazySets

A = [2..4 -2..1;-1..2 2..4]
b = [-2..2, -2..2]

polytopes = solve(A, b, LinearOettliPrager())

# convert h polytopes to Luxor objects
function luxify(P::HPolytope)

    # convert to vertex representation
    V = convert(VPolygon, P)
    out = vertices_list(P)

    out_luxor = [Luxor.Point(Tuple(p)) for p in out]

    return out_luxor
end

luxor_polytopes = luxify.(polytopes)

# draw logo
function logo(polytopes, fname)
    Drawing(500, 500, fname)
    origin()
    # transparent sethue("grey30")
    # squircle(O, 248, 248, :fill, rt=0.1)
    cols = [Luxor.julia_green, Luxor.julia_blue, Luxor.julia_purple, Luxor.julia_red]
    sizes = [1, 2.5, 4]
    Luxor.scale(1, -1)
    for i in 1:4
        sethue("black")
        poly(polytopes[i]*60, :stroke, close=true)

        sethue(cols[i])
        poly(polytopes[i]*60, :fill, close=true)
        #setline(rescale(sizes[i], 1, 4, 18, 18))
    end

    finish()
    preview()
end

function lockup(logo, fname)
    Drawing(1000, 300, fname)
    # background("grey90")     # otherwise transparent
    origin()
    logosvg = readsvg(logo) # this requires Luxor#master
    panes = Table([300], [300, 700])
    # place SVG
    @layer begin
        Luxor.translate(panes[1])
        Luxor.scale(0.5)
        placeimage(logosvg, O, centered=true)
    end
    # text
    @layer begin
        sethue("black")#sethue("rebeccapurple")
        setline(0.75)
        Luxor.translate(boxmiddleleft(BoundingBox(box(panes, 2))))
        fontsize(50) #80
        fontface("JuliaMono Bold")
        titlex = -40
        Luxor.text("IntervalLinearAlgebra.jl", Point(titlex, 0))
        @layer begin
            sethue("grey60")
            textoutlines("IntervalLinearAlgebra.jl", Point(titlex, 0), :stroke)
        end
        fontsize(40) #45
        fontface("JuliaMono")
        subx = -38
        Luxor.text("Linear algebra done rigorously", Point(subx, 45))
        @layer begin
            sethue("grey60")
            textoutlines("Linear algebra done rigorously", Point(subx, 45), :stroke)
        end
    end
    finish()
    preview()
end

logo(luxor_polytopes, "logo.svg")
lockup("logo.svg","logo-text.svg")


text:
module DA_NlayerTests

using Test

@testset "DA_Nlayer_functionTests" begin
    include("DA_Nlayer_functionTests.jl")
end

@testset "DA_Nlayer_kernelTests" begin
    include("DA_Nlayer_kernelTests.jl")
end

end # module

max_issues_repo_name: heltonmc/LightPropagation


text:
using Metatheory
using Metatheory.EGraphs
using Test

open_term = @theory begin
    # if-true
    cond(true, then, alt) => then
    cond(false, then, alt) => alt
    # if-elim
    cond(var(x) == e, then, alt) |>
        if addexpr!(_egraph, :(llet($x,$e,$then))) ==
            addexpr!(_egraph, :(llet($x,$e,$alt)))
            alt
        else _lhs_expr end
    a + b => b + a
    a + (b + c) => (a + b) + c
    (a == b) => (b == a)
end

subst_intro = @theory begin
    fix(v, e) => llet(v, fix(v,e), e)
    # beta reduction
    app(λ(v, body), e) => llet(v, e, body)
end

subst_prop = @theory begin
    # let-app
    llet(v, e, app(a, b)) => app(llet(v,e,a), llet(v,e,b))
    # let-add
    llet(v, e, a + b) => llet(v,e,a) + llet(v,e,b)
    # let-eq
    llet(v, e, a == b) => llet(v,e,a) == llet(v,e,b)
    # let-cond (let-if)
    llet(v, e, cond(guard, then, alt)) =>
        cond(llet(v,e,guard), llet(v,e,then), llet(v,e,alt))
end

subst_elim = @theory begin
    # let-const
    llet(v, e, c::Any) => c
    # let-var-same
    llet(v1, e, var(v1)) => e
    # TODO fancy let-var-diff
    llet(v1, e, var(v2)) |>
        if find(_egraph, v1) != find(_egraph, v2)
            :(var($v2))
        else _lhs_expr end
    # let-lam-same
    llet(v1, e, λ(v1, body)) => λ(v1, body)
    # let-lam-diff #TODO captureavoid
    llet(v1, e, λ(v2, body)) |>
        if v2.id ∈ getdata(e, FreeVarAnalysis, Set()) # is free
            :(λ($fresh, llet($v1, $e, llet($v2, var($fresh), $body))))
        else
            :(λ($v2, llet($v1, $e, $body)))
        end
end

λT = open_term ∪ subst_intro ∪ subst_prop ∪ subst_elim

ex = :(λ(x, 4 + app(λ(y, var(y)), 4)))
g = EGraph(ex)
# analyze!(g, FreeVarAnalysis)
saturate!(g, λT)
display(g.classes); println()
extract!(g, astsize)

@test @areequal λT 2 app(λ(x, var(x)), 2)

abstract type FreeVarAnalysis <: AbstractAnalysis end

function EGraphs.make(an::Type{FreeVarAnalysis}, g::EGraph, n::ENode)
    free = Set{Int64}()
    if n.head == :var
        push!(free, n.args[1])
    elseif n.head == :llet
        v,a,b = n.args[1:3]
        adata = getdata(g[a], an, Set{Int64}())
        bdata = getdata(g[a], an, Set{Int64}())
        union!(free, adata)
        delete!(free, v)
        union!(free, bdata)
    elseif n.head == :λ
        v,b = n.args[1:2]
        bdata = getdata(g[b], an, Set{Int64}())
        union!(free, bdata)
        delete!(free, v)
    end

    return free
end

function EGraphs.join(an::Type{FreeVarAnalysis}, from, to)
    union(from, to)
end

islazy(an::Type{FreeVarAnalysis}) = false


text:
function KrewerasEigenfunction(domain::SphericalTriangle)
    if domain.angles[1] != fmpq(2 // 3) ||
       domain.angles[2] != fmpq(2 // 3) ||
       domain.angles[3] != fmpq(2 // 3)
        throw(
            ArgumentError(
                "domain must be the Kreweras triangle with angles (2π/3, 2π/3, 2π/3)",
            ),
        )
    end
    KrewerasEigenfunction(domain, arb[])
end

function Base.show(io::IO, u::KrewerasEigenfunction)
    println(io, "Kreweras eigenfunction")
    if !haskey(io, :compact) || !io[:compact]
        println(io, "domain: $(u.domain)")
        print(io, "number of set coefficients: $(length(u.coefficients))")
    end
end

function active_boundaries(domain::SphericalTriangle, u::KrewerasEigenfunction)
    if domain === u.domain
        return 1:1
    else
        return 1:3
    end
end

function (u::KrewerasEigenfunction)(
    xyz::AbstractVector{T},
    λ::arb,
    k::Integer;
    boundary = nothing,
    notransform::Bool = false,
) where {T<:Union{arb,arb_series}}
    ν::arb = -0.5 + sqrt(0.25 + λ)
    if k % 2 == 1
        μ::arb = u.domain.parent(-k * inv(u.domain.angles[1]))

        # Create a linear map for rotating so that the next vertex ends up
        # at the north pole
        α = -acos(vertex(u.domain, 2)[3])
        β = angle(u.domain, 2) - u.domain.parent(π)
        L = LinearMap(RotZY(β, α))

        res = u.domain.parent(0)

        if isnothing(boundary) || boundary == 1
            ϕ = atan(xyz[2], xyz[1])
            res += legendre_p_safe(ν, μ, xyz[3]) * sin(μ * ϕ)
        end

        xyz = L(xyz)
        if isnothing(boundary) || boundary == 2
            ϕ = atan(xyz[2], xyz[1])
            res += legendre_p_safe(ν, μ, xyz[3]) * sin(μ * ϕ)
        end

        if isnothing(boundary) || boundary == 3
            xyz = L(xyz)
            ϕ = atan(xyz[2], xyz[1])
            res += legendre_p_safe(ν, μ, xyz[3]) * sin(μ * ϕ)
        end

        return res
    else
        μ = u.domain.parent(3div(k - 2, 2))

        (θ, φ) = spherical(center(u.domain))
        L = LinearMap(RotYZ(-θ, -φ))
        xyz = L(xyz)
        ϕ = atan(xyz[2], xyz[1])
        return legendre_p_safe(ν, μ, xyz[3]) * cos(μ * ϕ)
    end
end

function (u::KrewerasEigenfunction)(
    xyz::AbstractVector{T},
    λ::arb;
    boundary = nothing,
    notransform::Bool = false,
) where {T<:Union{arb,arb_series}}
    res = u.domain.parent(0)

    for k = 1:length(u.coefficients)
        res += u.coefficients[k] * u(xyz, λ, k, boundary = boundary)
        if (T == arb && !isfinite(res)) || (T == arb_series && !isfinite(res[end]))
            return res
        end
    end

    res
end


text:
println("Loading packages. This may take a few minutes ...")
flush(stdout)

# NOTE: These libs have to be imported in main process and worker processes.
import Pkg; Pkg.activate(joinpath(@__DIR__, "../../../../"))  # CytofRepFAM
using CytofRepFAM
using Random, Distributions, BSON, Dates, DataFrames, CSV

println("Finished loading packages. This may take a few minutes ...")
flush(stdout)


text:
using Distributed: nprocs, @everywhere
using Folds
using Test
using Transducers: DistributedEx  # TODO: stop special-casing
using Transducers: Eduction, SequentialEx, extract_transducer

const Comprehension = Union{
    Iterators.Generator,
    Iterators.Filter,
    Iterators.Flatten,
}

function down(xs::Union{Eduction,Comprehension})
    xf, foldable = extract_transducer(xs)
    return xf(down(foldable))
end

down(xs) = collect(xs)
down(xs::AbstractDict) = Dict(xs)
down(xs::AbstractSet) = Set(xs)

mapbottom(f, itr) = f(itr)

function mapbottom(f, itr::Base.Generator)
    ys = mapbottom(f, itr.iter)
    ys === nothing && return nothing
    return Base.Generator(itr.f, ys)
end

function mapbottom(f, itr::Iterators.Filter)
    ys = mapbottom(f, itr.itr)
    ys === nothing && return nothing
    return Iterators.Filter(itr.flt, ys)
end

function mapbottom(f, itr::Iterators.Flatten)
    ys = mapbottom(f, itr.it)
    ys === nothing && return nothing
    return Iterators.flatten(ys)
end

function mapbottom(f, itr::Tuple)
    ys = map(x -> mapbottom(f, x), itr)::Tuple
    any(isnothing, ys) && return nothing
    return ys
end

function getlabel((i, example),)
    if example isa NamedTuple
        return example.label
    else
        return "$i"
    end
end

function getdata((_, example),)
    if example isa NamedTuple
        return example.data
    else
        return example
    end
end

function getequality((_, example),)
    if example isa NamedTuple && haskey(example, :eq)
        return example.eq
    else
        return isequal
    end
end

macro data_str(s)
    data = include_string(__module__, s)
    return (label = s, data = data)
end

TESTCASES_WITH_SEQUENTIAL_RAWDATA = """
all(isodd(x) for x in 1:10)
all(isodd(x) for x in 1:2:10)
any(isodd(x) for x in 1:10)
any(isodd(x) for x in 2:2:10)
copy(Set, 1:10)                                                   # nodist,FIXME
copy(Set, Dict(zip('a':'j', 1:10)))                               # nodist,FIXME
copy(Dict, (x => x^2 for x in 1:10))
copy(Dict, (x => x^2 for x in Set(1:10)))                         # nodist,FIXME
copy(Vector, 1:10)
count(isodd(x) for x in 1:10)
dict(x => x^2 for x in 1:10)
extrema((x - 5)^2 for x in 1:10)
findall(isodd, 1:10)                                              # nodist,FIXME
findfirst(x -> x > 3, 1:10)                                       # nodist,FIXME
findlast(x -> x < 3, 1:10)                                        # nodist,FIXME
map(x -> x^2, 1:10)
map(+, 1:10, 11:20)                                               # nodist,FIXME
map(+, 1:10, 11:20, 21:30)                                        # nodist,FIXME
mapreduce(identity, +, 1:0)
mapreduce(x -> x^2, +, 1:10)
mapreduce(*, +, 1:10, 11:20)                                      # nodist,FIXME
mapreduce(*, +, 1:10, 11:20, 21:30)                               # nodist,FIXME
maximum(0:9)
maximum(9:-1:0)
maximum([2, 3, 0, 3, 4, 0, 5, 7, 4, 2])
maximum([1:10; [missing]])
minimum(0:9)
minimum(9:-1:0)
minimum([2, 3, 0, 3, 4, 0, 5, 7, 4, 2])
minimum([1:10; [missing]])
prod(1:2:10)
prod([1:10; [missing]])
prod(([x -x; -x x] for x in 1:2:19))
prod(y for x in 1:11 if isodd(x) for y in 1:x:x^2; init = 1)
reduce(+, 1:0)
set(1:10)                                                         # nodist,FIXME
set(x^2 for x in 1:10)
sum(1:10)
sum([1:10; [missing]])
sum(x^2 for x in 1:11)
sum(x^2 for x in 1:11 if isodd(x); init = 0)
unique(gcd(x, 42) for x in 1:30)
unique(x for x in 1:100 if x < 3)
issorted([1:5; 5:-1:0])                                           # nodist,FIXME
issorted(1:10)                                                    # nodist,FIXME
"""
#=
issorted(x^2 for x in 1:10)
issorted(x^2 for x in 1:10 if isodd(x))
=#

args_and_kwargs(args...; kwargs...) =
    (preargs = args[1:end-1], data = args[end], kwargs = (; kwargs...))

function parse_tests(str, _module)
    return map(split(str, "\n", keepempty = false)) do x
        @debug "Parsing: $x"
        fstr, rest = split(x, "(", limit = 2)
        ex = Meta.parse("DUMMY($rest")
        ex.args[1] = args_and_kwargs
        testcase = Base.eval(_module, ex)
        f = getproperty(Folds, Symbol(fstr))
        if (m = match(r"^(.*?) *# *(.*?) *$", x)) !== nothing
            label = m[1]
            tags = map(Symbol, split(m[2], ","))
        else
            label = x
            tags = Symbol[]
        end
        return (; label = label, tags = tags, f = f, testcase...)
    end
end

function _eval_collection_data(str, _module)
    return map(split(str, "\n", keepempty = false)) do x
        if (m = match(r"^(.*?) *# *(.*?) *$", x)) !== nothing
            label = m[1]
            tags = map(Symbol, split(m[2], ","))
        else
            label = x
            tags = Symbol[]
        end
        (label = label, data = Base.include_string(_module, x), tags = tags)
    end
end

eval_call_data(str, preamble::Union{Expr,Nothing} = nothing) =
    _eval_data(parse_tests, str, preamble)

eval_collection_data(str, preamble::Union{Expr,Nothing} = nothing) =
    _eval_data(_eval_collection_data, str, preamble)

function _eval_data(f, str, preamble)
    load_me_everywhere()
    @gensym folds_testing
    preamble = QuoteNode(preamble)
    _module = QuoteNode(folds_testing)
    @everywhere $_remote_eval_data($f, $str, $preamble, $_module)
    return getfield(getfield(Main, folds_testing), :DATA)
end

function _remote_eval_data(f, str, preamble, _module)
    @eval Main module $_module
        $preamble
        const DATA = $f($str, $_module)
    end
end

function _reeval()
    global TESTCASES_WITH_SEQUENTIAL_DEFAULT =
        parse_tests(TESTCASES_WITH_SEQUENTIAL_RAWDATA, @__MODULE__)
end

_reeval()

always(_) = true

function testcases_with_sequential(tests; upload = identity, include_test = always)
    tests = deepcopy(tests)
    tests = filter(include_test, tests)
    tests = map(tests) do ex
        (; ex..., data = mapbottom(upload, ex.data))
    end
    tests = filter(ex -> ex.data !== nothing, tests)
    return tests
end

test_with_sequential(executors; kwargs...) =
    test_with_sequential(TESTCASES_WITH_SEQUENTIAL_DEFAULT, executors; kwargs...)

function test_with_sequential(tests, executors; kwargs...)
    tests = testcases_with_sequential(tests; kwargs...)
    @testset "$(getlabel(x))" for x in enumerate(tests)
        @debug "test_with_sequential $(getlabel(x))"
        i, testcase = x
        f(args...) = testcase.f(testcase.preargs..., args...; testcase.kwargs...)
        ==′ = getequality(x)
        @testset "$ex" for ex in executors
            @debug "test_with_sequential $(getlabel(x)) $ex"
            :nodist in testcase.tags && ex isa DistributedEx && continue
            xs0 = down(getdata(x))
            xs1 = getdata(x)
            @test f(xs1, ex) ==′ f(xs0, SequentialEx())
        end
    end
end

function spinfor(d)
    limit = time_ns() + d
    while limit > time_ns()
        ccall(:jl_cpu_pause, Cvoid, ())
    end
end

function test_uses_threads(ex)
    @testset "test_uses_threads" begin
        if Threads.nthreads() == 1
            error("require `Threads.nthreads() > 1`")
        else
            ids = Folds.map(1:2*Threads.nthreads(), ex) do _
                spinfor(1000_000)
                Threads.threadid()
            end
            @test length(Set(ids)) > 1
        end
    end
end

is_ci() = get(ENV, "CI", "false") == "true"

function test_uses_processes(ex; autoskip = !is_ci())
    @testset "test_uses_processes" begin
        if nprocs() == 1
            if autoskip
                @warn """
                Skipping `test_uses_processes` since there is no Distributed workers.
                Run the test with environment variable `CI` set to `true` to run this test.
                """
            else
                error("require `nprocs() > 1`")
            end
        else
            load_me_everywhere()
            ids = Folds.map(1:2*nprocs(), ex) do _
                getpid()
            end
            @test length(Set(ids)) > 1
        end
    end
end

function load_me_everywhere()
    pkgid = Base.PkgId(@__MODULE__)
    @everywhere Base.require($pkgid)
end


text:
using KmeansAlgorithmTask
using Test

@testset "KmeansAlgorithmTask.jl" begin
    # Write your own tests here.
end


text:
using Plots
using Statistics, LinearAlgebra, Random
Random.seed!(1)
using Robotlib, StaticArrays, Test
import Robotlib: Rt2T, T2R, T2t, I4
# include(joinpath(dirname(@__FILE__),"..","src","calibNAXP.jl"))

const TM = SMatrix{4,4,Float64,16}

mutable struct Plane{T1,T2,T3,T4,T5}
    N::T1
    N1::T2
    d::T3
    Pn::T4
    Qn::T5
end

function generateRandomPlane(no)
    if no == 1
        N_RB = [1, 0, 0]
    elseif no == 2
        N_RB = [0, 1, 0]
    elseif no == 3
        N_RB = [0, 0, 1]
    else
        N_RB = normalize(randn(3)) # Normal of calibration plane, given in RB
    end
    d_RB = 1          # abs(0.001*randn);
    Pn = N_RB * N_RB' # Projection matrix onto normal
    Qn = I - Pn       # Projection matrix onto plane

    Plane(N_RB, [N_RB; 1], d_RB, Pn, Qn)
end

function generateRandomPoses(N_poses, σₙ, T_TF_S, plane)
    LN_S      = [0, 0, 1] # Normal of laser plane in sensor frame
    N_RB      = plane.N
    d_RB      = plane.d
    T_RB_TF   = Array{Float64}(undef, 4, 4, N_poses)
    lines_RB  = zeros(3, N_poses)
    lines_S   = zeros(3, N_poses)
    points_RB = zeros(3, N_poses)
    points_S  = zeros(4, N_poses)
    i = 1
    while i <= N_poses
        T_RB_TF[:, :, i] = Rt2T(
            rpy2R(1.5 * randn(), 1.5 * randn(), 1.5 * randn()),
            0.6 * randn(3) + 2 * N_RB,
        )
        T_RB_S               = T_RB_TF[:, :, i] * T_TF_S |> TM
        LN_RB                = T2R(T_RB_S) * LN_S
        intersection_line_RB = normalize(cross(N_RB, LN_RB))
        intersection_line_S  = inv(T2R(T_RB_S)) * intersection_line_RB
        lines_RB[:, i]       = intersection_line_RB
        lines_S[:, i]        = intersection_line_S
        laser_orig_RB        = T_RB_S[1:3, 4]
        laser_y_RB           = T_RB_S[1:3, 2]
        λ                    = (d_RB - N_RB'laser_orig_RB) / (N_RB'laser_y_RB)
        # λ is the distance from laser origin to plane, should be positive and small to be realistic, also, the angle with the plane normal should not be too large. In paper, the following (slightly unrealistic) condition was used
        # if abs(λ) > 2
        if abs(λ) > 0.2 || λ < 0 || (180 / pi .* acos(intersection_line_S[2])) > 120
            continue
        end
        points_RB[:, i] = laser_orig_RB + λ * laser_y_RB
        points_S[:, i] = inv(T_RB_S) * [points_RB[:, i]; 1]
        i += 1
    end
    points_S = points_S[1:3, :]
    points_S .+= σₙ * randn(size(points_S))

    # Verify poses
    ei          = eigen(cov(points_RB'))
    nhat_points = ei.vectors[:, 1]
    V           = eigen(cov(lines_RB')).vectors
    nhat_lines  = V[:, 1]
    if 1 - abs(nhat_points'nhat_lines) > 0.001 || ei.values[1, 1] > 1e-5
        @warn("Something is wrong")
        @show 1 - abs(nhat_points'nhat_lines), ei.values[1, 1]
    end

    points = zeros(4, N_poses)
    for i = 1:N_poses
        T_RB_S = T_RB_TF[:, :, i] * T_TF_S |> TM
        #         points(:,end+1) = T_RB_S*[(points_S[:,i]);1];
        points[:, i] = T_RB_S * [points_S[:, i]; 1]
        #         points(:,end+1) = inv(T_RB_S)*[(points_S[:,i] - 0.01*lines_S[:,i]);1];
    end
    if any(abs.(points[1:3, :] - points_RB) .> 1e-2)
        @warn("Something is wrong 2")
    end
    points_S, lines_S, T_RB_TF
end

MC = 10   # In paper 100
N_planes = 3    # In paper 3
iters = 80   # In paper 30
N_poses = 20   # In paper 10
σₙ = 5e-4 # In paper: 5e-4
function run_calib(verbose = false)
    SSEStart  = zeros(MC)
    normStart = zeros(MC)
    distStart = zeros(MC)
    rotStart  = zeros(MC)
    distEnd   = zeros(MC)
    rotEnd    = zeros(MC)
    SSEMC     = zeros(MC, iters)
    normMC    = zeros(MC, iters)
    distEnd2  = zeros(MC)
    rotEnd2   = zeros(MC)
    SSEMC2    = zeros(MC, iters)
    normMC2   = zeros(MC, iters)
    mc = 1
    for mc = 1:MC
        @info "MC iterations $(mc)/$(MC)"
        # Simulate measurement data
        T_TF_S = Rt2T(rpy2R(randn(3)), 0.5 * randn(3)) |> TM #TCP to sensor
        points_Sv = Matrix{Float64}[]
        lines_Sv = Matrix{Float64}[]
        T_RB_TFv = Array{Float64,3}[]
        plane = Plane[]
        SSE = zeros(N_planes)
        verbose && println("Genereating planes and poses")
        for j = 1:N_planes
            push!(plane, generateRandomPlane(j))
            points_St, lines_St, T_RB_TFt =
                generateRandomPoses(N_poses, σₙ, T_TF_S, plane[j])
            push!(points_Sv, points_St)
            push!(lines_Sv, lines_St)
            push!(T_RB_TFv, T_RB_TFt)
            SSE[j] = Robotlib.pointDiff(T_TF_S, T_RB_TFt, points_St)[1]
        end
        planes = repeat((1:N_planes)', N_poses)[:]
        points_S = cat(points_Sv..., dims = 2)
        lines_S = cat(lines_Sv..., dims = 2)
        T_RB_TF = cat(T_RB_TFv..., dims = 3)

        any(abs.(SSE .> 1e-3)) && @error("Points does not seem to lie on a plane")
        dist = [
            norm(T_RB_TF[1:4, 4, i] - T_TF_S * [points_S[:, i]; 1])
            for i = 1:size(T_RB_TF, 3)
        ]
        #     hist(dist,15)

        ## Calibration
        # display('----------------------------------------')
        # Generate nominal transformation
        T_TF_S_real = copy(T_TF_S)
        if true
            T_TF_S0 =
                Rt2T(rpy2R(60π / 180 * (rand(3) .- 0.5)), 0.4 * (rand(3) .- 0.5)) *
                T_TF_S_real
            # T_TF_S0 = Rt2T(rpy2R(10π/180*(rand(3).-0.5)),0.1*(rand(3).-0.5))*T_TF_S_real
            #         T_TF_S0 = rt2tr(rpy2R(2*(rand(1,3).-0.5),'deg'),0.01*(rand(3,1).-0.5))*T_TF_S_real;
        else
            T_TF_S0 = copy(T_TF_S_real)
        end

        for j = 1:N_planes
            ind = findall(planes .== j)
            SSE[j] = sqrt(Robotlib.pointDiff(T_TF_S0, T_RB_TF[:, :, ind], points_S[1:3, ind])[1])
        end

        SSEStart[mc] = mean(SSE)
        normStart[mc] = norm(vec(T_TF_S0 - T_TF_S_real))
        distStart[mc] = norm(T_TF_S0[1:3, 4] - T_TF_S_real[1:3, 4])
        rotStart[mc] = norm(180 / π * R2rpy(T_TF_S0 \ T_TF_S_real))
        # display(T_TF_S_real)
        T_TF_S, RMScalibs, ALLcalibs, norms = calibNAXP(
            points_S,
            lines_S,
            T_RB_TF,
            T_TF_S0,
            planes,
            iters,
            trueT_TF_S = T_TF_S_real,
            variance = (1, 1),
            doplot = true,
        )
        T_TF_S2, RMScalibs2, ALLcalibs2, norms2 = T_TF_S, RMScalibs, ALLcalibs, norms#Robotlib.Calibration.calibNAXP_bootstrap(points_S, lines_S, T_RB_TF, T_TF_S0, planes,  iters, trueT_TF_S=T_TF_S_real, nbootstrap=500)

        SSEMC[mc, :] = RMScalibs
        normMC[mc, :] = norms
        if RMScalibs[end] > 1e-2
            println("Bad result")
        end
        distEnd[mc] = norm(T_TF_S[1:3, 4] - T_TF_S_real[1:3, 4])
        rotEnd[mc] = norm(180 / π * R2rpy(T_TF_S \ T_TF_S_real))

        SSEMC2[mc, :] = RMScalibs2
        normMC2[mc, :] = norms2
        if RMScalibs2[end] > 1e-2
            println("Bad result")
        end
        distEnd2[mc] = norm(T_TF_S2[1:3, 4] - T_TF_S_real[1:3, 4])
        rotEnd2[mc] = norm(180 / π * R2rpy(T_TF_S2 \ T_TF_S_real))
    end
    SSEStart,
    normStart,
    distStart,
    rotStart,
    distEnd,
    rotEnd,
    SSEMC,
    normMC,
    distEnd2,
    rotEnd2,
    SSEMC2,
    normMC2
end

# @testset "CalibNAXP" begin

@test isapprox(Robotlib.pointDiff(I4, cat(fill(I4, 5)..., dims = 3), zeros(3, 5))[], 0, atol = 1e-10)

SSEStart, normStart, distStart, rotStart, distEnd, rotEnd, SSEMC, normMC, distEnd2, rotEnd2, SSEMC2, normMC2 = run_calib()

# iters = size(SSEMC,2)
# using StatPlots
# gr(legend=false)
# plot(0:iters,copy([normStart normMC]'),yscale=:log10,c=:black, xlabel="Number of iterations", layout=2, subplot=1)
# # plot!(0:iters,copy([normStart normMC2]'),yscale=:log10,c=:green, subplot=1)
# hline!([σₙ, σₙ],l=:dash, c=:red, subplot=1)
#
# plot!(0:iters,[SSEStart SSEMC]',yscale=:log10,c=:black, xlabel="Number of iterations",title="RMS distance from points to plane [m]", subplot=2)
# # plot!(0:iters,[SSEStart SSEMC2]',c=:green, subplot=2)
# hline!([σₙ σₙ],l=:dash,c=:red, subplot=2)
#
#
# boxplot(["Before" "After" "After WTLS"],([distStart distEnd distEnd2]), title="Distance error [m]", yscale=:log10, layout=2, subplot=1)
# hline!([σₙ σₙ],c=:red)
# boxplot!(["Before" "After" "After WTLS"],([rotStart rotEnd rotEnd2]),title="Rotation error [degree]", yscale=:log10, subplot=2)

# end

max_forks_repo_name: baggepinnen/Robotlib


text:
module PCREO

export lsdir, reldiff, to_point_vector, dict_push!,
    constrain_sphere!, spherical_riesz_gradient!,
    spherical_riesz_gradient, spherical_riesz_hessian,
    run!, refine,
    convex_hull_facets, incidence_degrees, adjacency_structure,
    nearest_neighbor_graph, write_graph6, canonical_graph6,
    spherical_circumcenter, covering_radius,
    PCREORecord

using DZOptimization: LBFGSOptimizer, step!, normalize_columns!
using DZOptimization.ExampleFunctions:
    riesz_energy, riesz_gradient!, constrain_riesz_gradient_sphere!,
    riesz_hessian!, constrain_riesz_hessian_sphere!
using Graphs: AbstractGraph, SimpleGraph, SimpleGraphFromIterator, nv
using Graphs.SimpleGraphs: SimpleEdge
using MultiFloats: Float64x2, Float64x3
using Printf: @printf
using StaticArrays: SVector, dot, cross, norm

############################################################## GENERAL UTILITIES

lsdir(path...; join=false) = filter(
    !startswith('.') ∘ basename,
    readdir(joinpath(path...); join=join, sort=false)
)

function reldiff(old, new)
    tmp = new - old
    return (tmp + tmp) / (abs(old) + abs(new))
end

function to_point_vector(::Val{N}, points::AbstractMatrix{T}) where {T,N}
    dimension, num_points = size(points)
    @assert dimension == N
    return [SVector{N,T}(view(points, :, i)) for i = 1 : num_points]
end

function dict_incr!(d::Dict{K,Int}, k::K) where {K}
    if haskey(d, k)
        d[k] += 1
    else
        d[k] = 1
    end
    return d[k]
end

function dict_push!(d::Dict{K,Vector{T}}, k::K, v::T) where {K,T}
    if haskey(d, k)
        push!(d[k], v)
    else
        d[k] = [v]
    end
    return d[k]
end

######################################################## RIESZ ENERGY ON SPHERES

function constrain_sphere!(points)
    normalize_columns!(points)
    return true
end

function spherical_riesz_gradient!(grad, points)
    riesz_gradient!(grad, points)
    constrain_riesz_gradient_sphere!(grad, points)
    return grad
end

spherical_riesz_gradient(points) =
    spherical_riesz_gradient!(similar(points), points)

function spherical_riesz_hessian(points::AbstractMatrix{T}) where {T}
    unconstrained_grad = riesz_gradient!(similar(points), points)
    hess = Array{T,4}(undef, size(points)..., size(points)...)
    riesz_hessian!(hess, points)
    constrain_riesz_hessian_sphere!(hess, points, unconstrained_grad)
    return reshape(hess, length(points), length(points))
end

################################################################### OPTIMIZATION

function run!(opt; quiet::Bool=true, framerate=10)
    if quiet
        while !opt.has_converged[]
            step!(opt)
        end
        return opt
    else
        last_print_time = time_ns()
        frame_time = round(Int, 1_000_000_000 / framerate)
        while !opt.has_converged[]
            step!(opt)
            if time_ns() >= last_print_time + frame_time
                println(opt.iteration_count[], '\t',
                        opt.current_objective_value[])
                last_print_time += frame_time
            end
        end
        println(opt.iteration_count[], '\t',
                opt.current_objective_value[])
        return opt
    end
end

function refine(points::Matrix{T}, initial_step_size::T,
                history_length::Int) where {T}
    constrain_sphere!(points)
    energy = riesz_energy(points)
    @assert isfinite(energy)
    while true
        opt = LBFGSOptimizer(
            riesz_energy, spherical_riesz_gradient!, constrain_sphere!,
            points, initial_step_size, history_length
        )
        run!(opt; quiet=!(stdout isa Base.TTY))
        next_energy = opt.current_objective_value[]
        next_points = opt.current_point
        if next_energy < energy
            energy = next_energy
            points = next_points
        else
            break
        end
    end
    return (points, energy)
end

###################################################################### ADJACENCY

function convex_hull_facets(points::AbstractVector{SVector{N,T}}) where {T,N}
    @assert length(points) >= 4
    num_retries = 0
    while true
        buffer = IOBuffer()
        process = open(`qconvex i`, buffer, write=true)
        println(process, N)
        println(process, length(points))
        for point in points
            for coord in point
                print(process, ' ', Float64(coord))
            end
            println(process)
        end
        close(process)
        wait(process)
        first = true
        num_facets = 0
        result = Vector{Int}[]
        seek(buffer, 0)
        for line in eachline(buffer)
            if first
                num_facets = parse(Int, line)
                first = false
            else
                push!(result, [parse(Int, s) + 1 for s in split(line)])
            end
        end
        if (num_facets > 0) && (num_facets == length(result))
            return result
        else
            # Sometimes, qconvex doesn't work, and we need to try again.
            num_retries += 1
            if num_retries >= 10
                error("qconvex failed to return a result after ten tries.")
            end
        end
    end
end

function incidence_degrees(facets::AbstractVector{Vector{Int}})
    degrees = Dict{Int,Int}()
    for facet in facets
        for vertex in facet
            dict_incr!(degrees, vertex)
        end
    end
    return degrees
end

function adjacency_structure(facets::AbstractVector{Vector{Int}})
    pair_dict = Dict{Tuple{Int,Int},Vector{Int}}()
    for (k, facet) in enumerate(facets)
        n = length(facet)
        for i = 1 : n-1
            for j = i+1 : n
                @inbounds dict_push!(pair_dict, minmax(facet[i], facet[j]), k)
            end
        end
    end
    adjacent_vertices = Vector{Tuple{Int,Int}}()
    adjacent_facets = Vector{Tuple{Int,Int}}()
    for (vertex_pair, incident_facets) in pair_dict
        if length(incident_facets) >= 2
            @assert length(incident_facets) == 2
            push!(adjacent_vertices, vertex_pair)
            push!(adjacent_facets, minmax(incident_facets...))
        end
    end
    return (adjacent_vertices, adjacent_facets)
end

###################################################################### GRAPH I/O

function nearest_neighbor_graph(filepath::AbstractString)
    contents = read(filepath, String)
    entries = split(contents, "\n\n")
    @assert length(entries) == 5
    pair_dict = Dict{Tuple{Int,Int},Int}()
    for line in split(entries[4], '\n')
        facet = [parse(Int, entry) for entry in split(line, ',')]
        n = length(facet)
        for i = 1 : n-1
            for j = i+1 : n
                @inbounds dict_incr!(pair_dict, minmax(facet[i], facet[j]))
            end
        end
    end
    return SimpleGraphFromIterator(
        SimpleEdge(vertex_pair)
        for (vertex_pair, num_incident_facets) in pair_dict
        if num_incident_facets >= 2
    )
end

@inline write_graph6_header(io::IO) = write(io, ">>graph6<<")

@inline write_graph6_byte(io::IO, b::UInt8) = write(io, b + UInt8(63))

function write_graph6_num_vertices(io::IO, n::Int)
    @assert 0 < n < Int64(2)^36
    if n < 63
        write_graph6_byte(io, UInt8(n))
    elseif n < 258048
        write(io, UInt8(126))
        write_graph6_byte(io, UInt8((0x000000000003f000 & n) >> 12))
        write_graph6_byte(io, UInt8((0x0000000000000fc0 & n) >>  6))
        write_graph6_byte(io, UInt8((0x000000000000003f & n) >>  0))
    else
        write(io, UInt8(126))
        write(io, UInt8(126))
        write_graph6_byte(io, UInt8((0x0000000fc0000000 & n) >> 30))
        write_graph6_byte(io, UInt8((0x000000003f000000 & n) >> 24))
        write_graph6_byte(io, UInt8((0x0000000000fc0000 & n) >> 18))
        write_graph6_byte(io, UInt8((0x000000000003f000 & n) >> 12))
        write_graph6_byte(io, UInt8((0x0000000000000fc0 & n) >>  6))
        write_graph6_byte(io, UInt8((0x000000000000003f & n) >>  0))
    end
    return nothing
end

function write_graph6_adjacency_matrix(io::IO, g::SimpleGraph{T}) where {T}
    n = nv(g)
    bits = UInt8(0)
    k = 0
    @inbounds for j = 2 : n
        next = 1
        for i in g.fadjlist[j]
            if i < j
                while next < i
                    bits <<= 1
                    k += 1
                    if k >= 6
                        write_graph6_byte(io, bits)
                        bits = UInt8(0)
                        k = 0
                    end
                    next += 1
                end
                bits = (bits << 1) | UInt8(1)
                k += 1
                if k >= 6
                    write_graph6_byte(io, bits)
                    bits = UInt8(0)
                    k = 0
                end
                next += 1
            else
                break
            end
        end
        while next < j
            bits <<= 1
            k += 1
            if k >= 6
                write_graph6_byte(io, bits)
                bits = UInt8(0)
                k = 0
            end
            next += 1
        end
    end
    if k > 0
        write(io, (bits << (6 - k)) + UInt8(63))
    end
    return nothing
end

function write_graph6(io::IO, g::SimpleGraph{T}) where {T}
    write_graph6_header(io)
    write_graph6_num_vertices(io, nv(g))
    write_graph6_adjacency_matrix(io, g)
    return nothing
end

############################################################# GRAPH CANONIZATION

function canonical_graph6(g::AbstractGraph)
    num_retries = 0
    while true
        inputbuf = IOBuffer()
        write_graph6(inputbuf, g)
        println(inputbuf)
        seek(inputbuf, 0)
        outputbuf = IOBuffer()
        run(pipeline(`labelg`; stdin=inputbuf, stdout=outputbuf,
                               stderr=devnull); wait=true)
        result = String(take!(outputbuf))
        if startswith(result, ">>graph6<<")
            return chomp(result)
        else
            num_retries += 1
            if num_retries >= 10
                error("labelg failed to return a result after ten tries.")
            end
        end
    end
end

canonical_graph6(filepath::AbstractString) =
    canonical_graph6(nearest_neighbor_graph(filepath))

####################################################################### GEOMETRY

function spherical_circumcenter(a::SVector{3,T}, b::SVector{3,T},
                                c::SVector{3,T}) where {T}
    p = cross(b - a, c - a)
    p /= norm(p)
    if signbit(dot(p, a))
        @assert signbit(dot(p, b))
        @assert signbit(dot(p, c))
        return -p
    else
        @assert !signbit(dot(p, b))
        @assert !signbit(dot(p, c))
        return +p
    end
end

function spherical_circumcenter(points::AbstractVector{SVector{3,T}},
                                facet::AbstractVector{Int}) where {T}
    n = length(facet)
    result = zero(SVector{3,T})
    for i = 1 : n-2
        for j = i+1 : n-1
            for k = j+1 : n
                @inbounds a = points[facet[i]]
                @inbounds b = points[facet[j]]
                @inbounds c = points[facet[k]]
                result += spherical_circumcenter(a, b, c)
            end
        end
    end
    return result / norm(result)
end

function covering_radius(points::AbstractVector{SVector{3,T}},
                         facets::AbstractVector{Vector{Int}}) where {T}
    result = zero(T)
    for facet in facets
        center = spherical_circumcenter(points, facet)
        for i in facet
            @inbounds result = max(result, norm(center - points[i]))
        end
    end
    return result
end

####################################################################### FILE I/O

struct PCREORecord
    dimension::Int
    num_points::Int
    symmetry_group::String
    energy::Float64x3
    first_eigenvalue::Float64
    packing_radius::Float64x2
    covering_radius::Float64x2
    num_gradient_bits::Float64
    num_hessian_bits::Float64
    facet_distance::Float64x2
    points::Matrix{Float64x2}
    facets::Vector{Vector{Int}}
    initial_points::Matrix{Float64}
end

function Base.show(io::IO, rec::PCREORecord)
    @assert size(rec.points) == (rec.dimension, rec.num_points)
    @assert size(rec.initial_points, 1) == rec.dimension
    println(io, rec.dimension)
    println(io, rec.num_points)
    println(io, '"', rec.symmetry_group, '"')
    println(io, rec.energy)
    println(io, rec.first_eigenvalue)
    @printf(io, "%.25f\n", rec.packing_radius)
    @printf(io, "%.25f\n", rec.covering_radius)
    println(io)
    @printf(io, "%.3f\n", rec.num_gradient_bits)
    @printf(io, "%.3f\n", rec.num_hessian_bits)
    @printf(io, "%.25f\n", rec.facet_distance)
    println(io)
    for j = 1 : rec.num_points
        for i = 1 : rec.dimension
            if i > 1
                print(io, ", ")
            end
            @printf(io, "%+.25f", rec.points[i,j])
        end
        println(io)
    end
    println(io)
    num_digits = length(string(rec.num_points))
    for facet in rec.facets
        println(io, join([lpad(index, num_digits, ' ')
                          for index in facet], ", "))
    end
    println(io)
    for j = 1 : size(rec.initial_points, 2)
        for i = 1 : rec.dimension
            if i > 1
                print(io, ", ")
            end
            @printf(io, "%+.17f", rec.initial_points[i,j])
        end
        println(io)
    end
    return nothing
end

function PCREORecord(points::AbstractMatrix{T},
                     initial_points::AbstractMatrix{U}) where {T,U}
    dimension, num_points = shape(points)

    return PCREORecord(
        dimension,
        num_points,
        "",
        Float64x3(0.0),
        0.0,
        Float64x2(0.0),
        Float64x2(0.0),
        0,0,
        0.0,
        Float64x2.(points),
        Vector{Int}[],
        Float64.(initial_points)
    )
end

function PCREORecord(filepath::AbstractString)

    # A PCREO Record is a plain text file containing five sections, delimited
    # by blank lines, that describes an energy-minimizing point configuration
    # on the unit sphere in d-dimensional Euclidean space.
    contents = read(filepath, String)
    entries = split(contents, "\n\n")
    @assert length(entries) == 5

    # The first section contains exactly seven lines, each containing a single
    # number or string describing a geometric property of the point
    # configuration.
    geometric_properties = split(entries[1], '\n')
    @assert length(geometric_properties) == 7

    # The first two lines contain the dimension of the point configuration
    # (i.e., the number of coordinates in each point) and the total number
    # of points in the configuration. We use Cartesian coordinates, so the
    # intrinsic topological dimension of the unit sphere being considered is
    # one lower than the dimension stored here (e.g., S^2 embedded in R^3).
    dimension = parse(Int, geometric_properties[1])
    @assert dimension > 0
    num_points = parse(Int, geometric_properties[2])
    @assert num_points > 0

    # The third line contains the symmetry group of the point configuration
    # in Schoenflies notation. The string must begin and end with double
    # quotes (e.g., "C_2v").
    symmetry_group = geometric_properties[3]
    @assert length(symmetry_group) > 2
    @assert symmetry_group[1] == '"'
    @assert symmetry_group[end] == '"'
    symmetry_group = symmetry_group[2:end-1]

    # The fourth line contains the Riesz 1-energy of the point configuration.
    energy = Float64x3(geometric_properties[4])

    # The fifth line contains the first (expected) nonzero eigenvalue of the
    # Hessian matrix of the Riesz energy function. This measures the steepness
    # of the "energy bowl" this point configuration lies in. Note that many
    # eigenvalues of the Hessian matrix are expected to be zero due to the
    # constraint that each point lie on the unit sphere, together with the
    # rotational symmetry of the sphere.
    first_eigenvalue = parse(Float64, geometric_properties[5])

    # The sixth and seventh lines contain the packing and covering radii for
    # the point configuration, respectively.
    packing_radius = Float64x2(geometric_properties[6])
    covering_radius = Float64x2(geometric_properties[7])

    # The second section contains exactly three lines, each containing a single
    # number describing the precision of the point configuration. PCREO
    # produces energy-minimizing point configurations using iterative nonlinear
    # optimization algorithms, so these numbers serve as convergence measures.
    convergence_measures = split(entries[2], '\n')
    @assert length(convergence_measures) == 3

    num_gradient_bits = parse(Float64, convergence_measures[1])
    num_hessian_bits = parse(Float64, convergence_measures[2])
    facet_distance = Float64x2(convergence_measures[3])

    # The third section contains the Cartesian coordinates of the points in the
    # energy-minimizing configuration. Each line contains the Cartesian
    # coordinates of one point, separated by commas.
    points = reshape(
        [
            Float64x2(strip(entry))
            for line in split(entries[3], '\n')
            for entry in split(line, ',')
            if !isempty(strip(line))
        ],
        (dimension, num_points)
    )

    # The fourth section contains the facets of the convex hull of the
    # energy-minimizing point configuration. Each facet is described by a
    # comma-separated list of its vertices. Instead of listing the coordinates
    # of the vertices, we describe each vertex by its index in the point
    # configuration. We use 1-based indexing, so there is no point 0.
    facets = [
        [
            parse(Int, strip(entry))
            for entry in split(line, ',')
        ]
        for line in split(entries[4], '\n')
        if !isempty(strip(line))
    ]

    # The fifth section contains the initial (uniformly random) points that
    # this energy-minimizing point configuration was produced from.
    initial_points = reshape(
        [
            parse(Float64, strip(entry))
            for line in split(entries[5], '\n')
            for entry in split(line, ',')
            if !isempty(strip(line))
        ],
        # Point configurations produced with the symmetry-constrained optimizer
        # have a smaller number of initial points than total points, so we
        # leave the second dimension free to vary.
        (dimension, :)
    )

    return PCREORecord(
        dimension,
        num_points,
        symmetry_group,
        energy,
        first_eigenvalue,
        packing_radius,
        covering_radius,
        num_gradient_bits,
        num_hessian_bits,
        facet_distance,
        points,
        facets,
        initial_points
    )
end

end # module PCREO

# ########################################################### FILE NAMES AND PATHS

# # const PCREO_DIRNAME_REGEX = Regex(
# #     "^PCREO-([0-9]{2})-([0-9]{4})-([0-9a-fA-F]{8}-[0-9a-fA-F]{4}-" *
# #     "[0-9a-fA-F]{4}-[0-9a-fA-F]{4}-[0-9a-fA-F]{12})\$")

# # const PCREO_FILENAME_REGEX = Regex(
# #     "^PCREO-([0-9]{2})-([0-9]{4})-([0-9a-fA-F]{8}-[0-9a-fA-F]{4}-" *
# #     "[0-9a-fA-F]{4}-[0-9a-fA-F]{4}-[0-9a-fA-F]{12})\\.csv\$")

# ########################################################### GEOMETRIC PROPERTIES

# function packing_radius(points::Vector{SVector{N,T}},
#                         facets::Vector{Vector{Int}}) where {T,N}
#     adjacent_vertices, _ = adjacency_structure(facets)
#     result = typemax(T)
#     for (i, j) in adjacent_vertices
#         result = min(result, norm(points[i] - points[j]))
#     end
#     return half(T) * result
# end

# middle(x::AbstractVector) = x[(length(x) + 1) >> 1]

# ############################################################ CONVERGENCE TESTING

# function symmetrize!(mat::Matrix{T}) where {T}
#     m, n = size(mat)
#     @assert m == n
#     @inbounds for i = 1 : n-1
#         @simd ivdep for j = i+1 : n
#             sym = half(T) * (mat[i, j] + mat[j, i])
#             mat[i, j] = mat[j, i] = sym
#         end
#     end
#     return mat
# end

# function facet_normal_vector(points::Vector{SVector{3,T}}) where {T}
#     n = length(points)
#     result = zero(SVector{3,T})
#     for i = 1 : n-2
#         for j = i+1 : n-1
#             for k = j+1 : n
#                 normal = cross(points[j] - points[i],
#                                points[k] - points[i])
#                 normal /= norm(normal)
#                 positive = all(!signbit, normal' * p for p in points)
#                 negative = all(signbit, normal' * p for p in points)
#                 @assert xor(positive, negative)
#                 if positive
#                     result += normal
#                 else
#                     result -= normal
#                 end
#             end
#         end
#     end
#     return result / norm(result)
# end

# function parallel_facet_distance(points::Vector{SVector{3,T}},
#                                  facets::Vector{Vector{Int}}) where {T}
#     _, adjacent_facets = adjacency_structure(facets)
#     normals = [facet_normal_vector(points[facet]) for facet in facets]
#     return minimum(one(T) - normals[i]' * normals[j]
#                    for (i, j) in adjacent_facets)
# end

# ############################################################ TOPOLOGICAL DEFECTS

# function defect_graph(facets::Vector{Vector{Int}})
#     adjacent_vertices, _ = adjacency_structure(facets)
#     adjacency_lists = Dict{Int,Vector{Int}}()
#     for (v, w) in adjacent_vertices
#         dict_push!(adjacency_lists, v, w)
#         dict_push!(adjacency_lists, w, v)
#     end
#     degrees = Dict(v => length(l) for (v, l) in adjacency_lists)
#     @assert degrees == incidence_degrees(facets)
#     hexagonal_vertices = [v for (v, d) in degrees if d == 6]
#     for k in hexagonal_vertices
#         delete!(adjacency_lists, k)
#         delete!(degrees, k)
#     end
#     for (v, l) in adjacency_lists
#         deleteat!(adjacency_lists[v],
#             [i for (i, w) in enumerate(l)
#              if w in hexagonal_vertices])
#     end
#     return (adjacency_lists, degrees)
# end

# function connected_components(adjacency_lists::Dict{V,Vector{V}}) where {V}
#     visited = Dict{V,Bool}()
#     for (v, l) in adjacency_lists
#         visited[v] = false
#     end
#     components = Vector{V}[]
#     for (v, l) in adjacency_lists
#         if !visited[v]
#             visited[v] = true
#             current_component = [v]
#             to_visit = copy(l)
#             while !isempty(to_visit)
#                 w = pop!(to_visit)
#                 if !visited[w]
#                     visited[w] = true
#                     push!(current_component, w)
#                     append!(to_visit, adjacency_lists[w])
#                 end
#             end
#             push!(components, current_component)
#         end
#     end
#     @assert allunique(vcat(components...))
#     return components
# end

# function defect_classes(facets::Vector{Vector{Int}})
#     adjacency_lists, defect_degrees = defect_graph(facets)
#     defect_components = connected_components(adjacency_lists)
#     defect_counts = Dict{Vector{Tuple{Int,Tuple{Int,Int}}},Int}()
#     for component in defect_components
#         shape_counts = Dict{Tuple{Int,Int},Int}()
#         for v in component
#             dict_incr!(shape_counts,
#                 (length(adjacency_lists[v]), defect_degrees[v]))
#         end
#         shape_table = [(num, shape) for (shape, num) in shape_counts]
#         sort!(shape_table; rev=true)
#         dict_incr!(defect_counts, shape_table)
#     end
#     defect_table = [(num, defect) for (defect, num) in defect_counts]
#     sort!(defect_table; rev=true)
#     return defect_table
# end

# function shape_code(n::Int)
#     if n == 3
#         return 'T'
#     elseif n == 4
#         return 'S'
#     elseif n == 5
#         return 'P'
#     elseif n == 6
#         @assert n != 6
#     elseif n == 7
#         return 'H'
#     elseif n == 8
#         return 'O'
#     elseif n == 9
#         return 'N'
#     elseif n == 10
#         return 'D'
#     elseif n == 11
#         return 'U'
#     else
#         @assert false
#     end
# end

# unicode_subscript_string(n::Int) = foldl(replace,
#     map(Pair, Char.(0x30:0x39), Char.(0x2080:0x2089));
#     init=string(n))

# html_subscript_string(n::Int) = "<sub>" * string(n) * "</sub>"

# unicode_defect_string(shape_table::Vector{Tuple{Int,Tuple{Int,Int}}}) =
#     join([
#         (num == 1 ? "" : string(num)) *
#             shape_code(shape) * unicode_subscript_string(degree)
#         for (num, (degree, shape)) in shape_table])

# html_defect_string(shape_table::Vector{Tuple{Int,Tuple{Int,Int}}}) =
#     join([
#         (num == 1 ? "" : string(num)) *
#             shape_code(shape) * html_subscript_string(degree)
#         for (num, (degree, shape)) in shape_table])

# ####################################################### SYMMETRIZED RIESZ ENERGY

# # Benchmarked in Julia 1.5.3 for zero allocations or exceptions.

# # @benchmark symmetrized_riesz_energy(points, group, external_points) setup=(
# #     points=randn(3, 10); group=chiral_tetrahedral_group(Float64);
# #     external_points=SVector{3,Float64}.(eachcol(randn(3, 5))))

# # view_asm(symmetrized_riesz_energy,
# #     Matrix{Float64},
# #     Vector{SArray{Tuple{3,3},Float64,2,9}},
# #     Vector{SVector{3,Float64}})

# function symmetrized_riesz_energy(
#         points::AbstractMatrix{T},
#         group::Vector{SArray{Tuple{N,N},T,2,M}},
#         external_points::Vector{SArray{Tuple{N},T,1,N}}) where {T,N,M}
#     dim, num_points = size(points)
#     group_size = length(group)
#     num_external_points = length(external_points)
#     energy = zero(T)
#     for i = 1 : num_points
#         @inbounds p = SVector{N,T}(view(points, 1:N, i))
#         for j = 2 : group_size
#             @inbounds g = group[j]
#             energy += 0.5 * inv(norm(g*p - p))
#         end
#     end
#     for i = 2 : num_points
#         @inbounds p = SVector{N,T}(view(points, 1:N, i))
#         for g in group
#             gp = g * p
#             for j = 1 : i-1
#                 @inbounds q = SVector{N,T}(view(points, 1:N, j))
#                 energy += inv(norm(gp - q))
#             end
#         end
#     end
#     energy *= group_size
#     for i = 1 : num_points
#         @inbounds p = SVector{N,T}(view(points, 1:N, i))
#         for g in group
#             gp = g * p
#             for j = 1 : num_external_points
#                 @inbounds q = external_points[j]
#                 energy += inv(norm(gp - q))
#             end
#         end
#     end
#     return energy
# end

# # Benchmarked in Julia 1.5.3 for zero allocations or exceptions.

# # @benchmark symmetrized_riesz_gradient!(
# #     grad, points, group, external_points) setup=(
# #     points=randn(3, 10); grad=similar(points);
# #     group=chiral_tetrahedral_group(Float64);
# #     external_points=SVector{3,Float64}.(eachcol(randn(3, 5))))

# # view_asm(symmetrized_riesz_gradient!,
# #     Matrix{Float64}, Matrix{Float64},
# #     Vector{SArray{Tuple{3,3},Float64,2,9}},
# #     Vector{SVector{3,Float64}})

# function symmetrized_riesz_gradient!(
#         grad::AbstractMatrix{T},
#         points::AbstractMatrix{T},
#         group::Vector{SArray{Tuple{N,N},T,2,M}},
#         external_points::Vector{SArray{Tuple{N},T,1,N}}) where {T,N,M}
#     dim, num_points = size(points)
#     group_size = length(group)
#     num_external_points = length(external_points)
#     for i = 1 : num_points
#         @inbounds p = SVector{N,T}(view(points, 1:N, i))
#         force = zero(SVector{N,T})
#         for j = 2 : group_size
#             @inbounds r = group[j] * p - p
#             force += r / norm(r)^3
#         end
#         for j = 1 : num_points
#             if i != j
#                 @inbounds q = SVector{N,T}(view(points, 1:N, j))
#                 for g in group
#                     r = g * q - p
#                     force += r / norm(r)^3
#                 end
#             end
#         end
#         force *= group_size
#         for j = 1 : num_external_points
#             @inbounds q = external_points[j]
#             for g in group
#                 r = g * q - p
#                 force += r / norm(r)^3
#             end
#         end
#         @simd ivdep for j = 1 : N
#             @inbounds grad[j,i] = force[j]
#         end
#     end
#     return grad
# end

# struct SymmetrizedRieszEnergyFunctor{T}
#     group::Vector{SArray{Tuple{3,3},T,2,9}}
#     external_points::Vector{SArray{Tuple{3},T,1,3}}
#     external_energy::T
# end

# struct SymmetrizedRieszGradientFunctor{T}
#     group::Vector{SArray{Tuple{3,3},T,2,9}}
#     external_points::Vector{SArray{Tuple{3},T,1,3}}
# end

# function (sref::SymmetrizedRieszEnergyFunctor{T})(
#           points::AbstractMatrix{T}) where {T}
#     return sref.external_energy + symmetrized_riesz_energy(
#         points, sref.group, sref.external_points)
# end

# function (srgf::SymmetrizedRieszGradientFunctor{T})(
#           grad::AbstractMatrix{T}, points::AbstractMatrix{T}) where {T}
#     symmetrized_riesz_gradient!(grad, points, srgf.group, srgf.external_points)
#     constrain_riesz_gradient_sphere!(grad, points)
#     return grad
# end

# function symmetrized_riesz_functors(
#         ::Type{T}, group_function::Function,
#         orbit_functions::Vector{Function}) where {T}
#     group = group_function(T)::Vector{SArray{Tuple{3,3},T,2,9}}
#     external_points = vcat([orbit_function(T)::Vector{SArray{Tuple{3},T,1,3}}
#                             for orbit_function in orbit_functions]...)
#     external_points_matrix = Matrix{T}(undef, 3, length(external_points))
#     for (i, point) in enumerate(external_points)
#         @simd ivdep for j = 1 : 3
#             @inbounds external_points_matrix[j,i] = point[j]
#         end
#     end
#     external_energy = riesz_energy(external_points_matrix)
#     return (SymmetrizedRieszEnergyFunctor{T}(group, external_points,
#                                              external_energy),
#             SymmetrizedRieszGradientFunctor{T}(group, external_points))
# end


text:
__precompile__()
module SeisIO
using Blosc, Dates, DSP, FFTW, LightXML, LinearAlgebra, Printf, Sockets
using DelimitedFiles: readdlm
using Glob: glob
using HTTP: request, Messages.statustext
using Statistics: mean
Blosc.set_compressor("lz4")
Blosc.set_num_threads(Sys.CPU_THREADS)
path = Base.source_dir()

# DO NOT CHANGE IMPORT ORDER
include("imports.jl")
include("constants.jl")

# =========================================================
# CoreUtils: SeisIO needs these for core functions
# DO NOT CHANGE ORDER OF INCLUSIONS
include("CoreUtils/ls.jl")
include("CoreUtils/time.jl")
include("CoreUtils/namestrip.jl")
include("CoreUtils/typ2code.jl")
include("CoreUtils/file_io.jl")
include("CoreUtils/poly.jl")
include("CoreUtils/calculus.jl")
include("CoreUtils/mkchans.jl")
include("CoreUtils/svn.jl")

# =========================================================
# Types and methods
# DO NOT CHANGE ORDER OF INCLUSIONS

# Back-end types
include("Types/KWDefs.jl")

# Prereqs for custom types
include("Types/InstPosition.jl")
include("Types/InstResp.jl")

# IO buffer including location and instrument response buffers
include("Types/SeisIOBuf.jl")

# Abstract types
include("Types/GphysData.jl")
include("Types/GphysChannel.jl")

# Custom types
include("Types/SeisData.jl")
include("Types/SeisChannel.jl")

for i in readdir(path*"/Types/Methods")
  if endswith(i, ".jl")
    include(joinpath("Types/Methods",i))
  end
end

# =========================================================
# Utilities that may require SeisIO types to work
for i in ls(path*"/Utils/")
  if endswith(i, ".jl")
    include(i)
  end
end
include("Utils/Parsing/streams.jl")
include("Utils/Parsing/strings.jl")
include("Utils/Parsing/buffers.jl")

# =========================================================
# Data processing operations
for i in ls(path*"/Processing/*")
  if endswith(i, ".jl")
    include(i)
  end
end

# =========================================================
# Data formats
for i in ls(path*"/Formats/*")
  if endswith(i, ".jl")
    include(i)
  end
end

# =========================================================
# Web clients
for i in ls(path*"/Web/")
  if endswith(i, ".jl")
    include(i)
  end
end

# =========================================================
# Submodules
include("Submodules/FormatGuide.jl")
import .Formats: formats

include("Submodules/SEED.jl")
using .SEED: mseed_support, parsemseed!, parserec!, read_dataless, read_mseed_file!, read_seed_resp!, read_seed_resp, RESP_wont_read, seed_cleanup!, seed_support
export mseed_support, read_dataless, read_seed_resp!, read_seed_resp, RESP_wont_read, seed_support

include("Submodules/Quake.jl")
using .Quake:read_qml
import .Quake:convert
export read_qml, write_qml
include("Submodules/RandSeis.jl")

include("Submodules/SeisHDF.jl")
using .SeisHDF: read_hdf5, read_hdf5!, scan_hdf5, write_hdf5
export read_hdf5, read_hdf5!, scan_hdf5, write_hdf5

include("Submodules/SUDS.jl")
include("Submodules/UW.jl")
using .UW: formats

# =========================================================
# Wrappers
for i in ls(path*"/Wrappers/")
  if endswith(i, ".jl")
    include(i)
  end
end

# We need these types for the native file format
using .Quake: EQLoc, EQMag, EventChannel, EventTraceData, PhaseCat, SeisEvent, SeisHdr, SeisPha, SeisSrc, SourceTime
formats["list"] = collect(keys(formats))

# Last steps
include("Last/splat.jl")
include("Last/native_file_io.jl")
include("Last/read_legacy.jl")
include("Last/set_file_ver.jl")

# Module ends
end


text:
import PowerModelsACDC;
const _PMACDC = PowerModelsACDC;
import PowerModels;
const _PM = PowerModels;
import InfrastructureModels;
const _IM = InfrastructureModels;
import JuMP
import Gurobi
using MAT
using XLSX
using JLD2
using Statistics
include("basencont_nw.jl")

Total_sample = 1  # sample per year
total_yr = 1# the years in horizon, data coming from excels

Prot_system = "Permanentloss"
# Prot_system_coll = ["FS_HDCCB", "NS_CB", "Permanentloss"]
curtailed_gen = [1] #geneartor numbers # change also constraint max(), generating power,file name
syncarea = 2
max_curt = 0
# for proti = 1:3
#     Prot_system = Prot_system_coll[proti]
file = "./test/data/4bus_OPF_PLdim.m"
data_sp = _PM.parse_file(file)
_PMACDC.process_additional_data!(data_sp)

include("bkgrnd_cal.jl")
gurobi = JuMP.optimizer_with_attributes(Gurobi.Optimizer, "Presolve" => -1)

    no_nw = 2 * Total_sample*total_yr
    data_mp = multi_network(file, no_nw)
    # year = ["NAT_2050_Generation"]
    yr = 1

    data_cont, Cont_list, base_list = mp_contignecy_nocl_sensitivity(deepcopy(data_mp), Total_sample*total_yr, 1)
    year_base = year_base_networks(Total_sample, total_yr, base_list)

    conv_rate = Int(data_cont["nw"]["1"]["convdc"]["1"]["Pacmax"]*100)
    Clearingtime = [0.15 0.3 0.45 0.6]
    column = ["A" "B" "C" "D" "E" "F" "G" "H"]
    filepath = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\plots\\OPF\\4bus\\sensitivity\\",conv_rate,"_",Prot_system,".xlsx")
    XLSX.openxlsx(filepath, mode="w") do xf
        for i = 1:length(Clearingtime)
            for (n,nw) in data_cont["nw"]
                for (r,reserves) in nw["reserves"]
                     reserves["Tcl"] = Clearingtime[i]
                end
            end
    if Prot_system == "FS_HDCCB" || Prot_system == "FS_MDCCB"
    s = Dict("output" => Dict("branch_flows" => true), "conv_losses_mp" => true, "process_data_internally" => false, "FSprotection" => true, "NSprotection" => false,
    "Permanentloss" => false, "Cont_list" => Cont_list,"base_list" => base_list,"Total_sample" =>Total_sample, "curtailed_gen" => curtailed_gen, "max_curt" => max_curt, "syncarea" => syncarea, "year_base" => year_base, "total_yr" => total_yr)
    elseif Prot_system == "NS_CB" || Prot_system == "NS_FB"
    s = Dict("output" => Dict("branch_flows" => true), "conv_losses_mp" => true, "process_data_internally" => false, "FSprotection" => false, "NSprotection" => true,
    "Permanentloss" => false, "Cont_list" => Cont_list,"base_list" => base_list,"Total_sample" =>Total_sample, "curtailed_gen" => curtailed_gen, "max_curt" => max_curt, "syncarea" => syncarea,  "year_base" => year_base, "total_yr" => total_yr)
    else
    s = Dict("output" => Dict("branch_flows" => true), "conv_losses_mp" => true, "process_data_internally" => false, "FSprotection" => false, "NSprotection" => false,
    "Permanentloss" => true, "Cont_list" => Cont_list,"base_list" => base_list,"Total_sample" =>Total_sample, "curtailed_gen" => curtailed_gen, "max_curt" => max_curt, "syncarea" => syncarea, "year_base" => year_base, "total_yr" => total_yr)
    end
    @assert ( s["FSprotection"] == true && (Prot_system == "FS_HDCCB" || Prot_system == "FS_MDCCB") ) || (s["NSprotection"] == true && (Prot_system == "NS_CB" || Prot_system == "NS_FB"))|| (Prot_system == "Permanentloss" && s["NSprotection"] == false && s["FSprotection"] == false)
     resultDC1 = _PMACDC.run_acdcscopf_sensitivity(data_cont, _PM.DCPPowerModel, gurobi, multinetwork=true; setting = s)
     # curtail, maxFFR, maxFCR, meanFFR, meanFCR = curtailment(data_cont, base_list, resultDC1, curtailed_gen)
     display_keyindictionary_OPF_PLdim(resultDC1, "isbuilt", "Pgg")
     sheet = xf[1]
     cell_no = string(column[i],i)
     sheet["$(string(column[i],1))"] = data_cont["nw"]["1"]["reserves"]["2"]["Tcl"]
     sheet["$(string(column[i],2))"] = resultDC1["solution"]["nw"]["1"]["FFR_Reserves"]
     sheet["$(string(column[i],3))"] = resultDC1["solution"]["nw"]["1"]["FCR_Reserves"]
     sheet["$(string(column[i],4))"] = resultDC1["solution"]["nw"]["1"]["Gen_cost"]
     # sheet["$(string("A",5))"] = resultDC1["solution"]["nw"]["1"]["Cont"]
     sheet["$(string(column[i],5))"] = mean(resultDC1["solution"]["nw"]["1"]["Curt"])
     sheet["$(string(column[i],6))"] = resultDC1["objective"]
     sheet["$(string(column[i],7))"] = resultDC1["solution"]["nw"]["2"]["reserves"]["2"]["z1"]
     sheet["$(string(column[i],8))"] = resultDC1["solution"]["nw"]["2"]["reserves"]["2"]["z2"]
     sheet["$(string(column[i],9))"] = resultDC1["solution"]["nw"]["2"]["reserves"]["2"]["z3"]
     sheet["$(string(column[i],10))"] = resultDC1["solution"]["nw"]["2"]["reserves"]["2"]["z4"]
     sheet["$(string(column[i],11))"] =  resultDC1["solution"]["nw"]["1"]["branchdc"]["1"]["pf"]
     sheet["$(string(column[i],12))"] =  resultDC1["solution"]["nw"]["2"]["branchdc"]["1"]["pf"]

     sheet["$(string(column[i],13))"] = resultDC1["solution"]["nw"]["2"]["reserves"]["2"]["k11_dup"]
     sheet["$(string(column[i],14))"] = resultDC1["solution"]["nw"]["2"]["reserves"]["2"]["k12_dup"]
     sheet["$(string(column[i],15))"] = resultDC1["solution"]["nw"]["2"]["reserves"]["2"]["k21_dup"]
     sheet["$(string(column[i],16))"] = resultDC1["solution"]["nw"]["2"]["reserves"]["2"]["k22_dup"]
     sheet["$(string(column[i],17))"] = resultDC1["solution"]["nw"]["2"]["reserves"]["2"]["k23_dup"]
     sheet["$(string(column[i],18))"] = resultDC1["solution"]["nw"]["2"]["reserves"]["2"]["k31"]
     sheet["$(string(column[i],19))"] = resultDC1["solution"]["nw"]["2"]["reserves"]["2"]["k32"]
     sheet["$(string(column[i],20))"] = resultDC1["solution"]["nw"]["2"]["reserves"]["2"]["k33"]
     sheet["$(string(column[i],21))"] = resultDC1["solution"]["nw"]["2"]["reserves"]["2"]["k31_dup"]
     sheet["$(string(column[i],22))"] = resultDC1["solution"]["nw"]["2"]["reserves"]["2"]["k32_dup"]
     sheet["$(string(column[i],23))"] = resultDC1["solution"]["nw"]["2"]["reserves"]["2"]["k43_dup"]

     if round(Int64,resultDC1["solution"]["nw"]["2"]["reserves"]["2"]["z3"]) == 1
         t = nadir_time_t3(resultDC1, data_cont)
         display("t3")
     elseif round(Int64,resultDC1["solution"]["nw"]["2"]["reserves"]["2"]["z4"]) == 1
         t = nadir_time_t4(resultDC1, data_cont)
         display("t4")
     else
         t = 500
     end
     sheet["$(string(column[i],24))"] =  t
    end
end
# end

function nadir_time_t3(resultDC1, data_cont)
    Phvdcoaux = resultDC1["solution"]["nw"]["2"]["reserves"]["2"]["Phvdcoaux"]
    Phvdccaux = resultDC1["solution"]["nw"]["2"]["reserves"]["2"]["Phvdccaux"]
    display("Phvdcoaux:$Phvdcoaux")
    display("Phvdccaux:$Phvdccaux")
    Pg = resultDC1["solution"]["nw"]["2"]["reserves"]["2"]["Pgg"]/data_cont["nw"]["1"]["load"]["1"]["pd"]
    Pf = resultDC1["solution"]["nw"]["2"]["reserves"]["2"]["Pff"]/data_cont["nw"]["1"]["load"]["1"]["pd"]
    Tcl = data_cont["nw"]["1"]["reserves"]["2"]["Tcl"]
    Td = data_cont["nw"]["1"]["reserves"]["2"]["Td"]
    Tg = data_cont["nw"]["1"]["reserves"]["2"]["Tg"]
    Tf = data_cont["nw"]["1"]["reserves"]["2"]["Tf"]
    display("Pg:$Pg")
    display("Pf:$Pf")
    display(data_cont["nw"]["1"]["load"]["1"]["pd"])
    t3 = (Phvdcoaux  - Pf +  Phvdccaux*Tcl/Td)/((Pg/Tg) + (Phvdccaux/Td))
    display("t3:$t3")
    return t3
end

function nadir_time_t4(resultDC1, data_cont)
    Phvdcoaux = resultDC1["solution"]["nw"]["2"]["reserves"]["2"]["Phvdcoaux"]
    Phvdccaux = resultDC1["solution"]["nw"]["2"]["reserves"]["2"]["Phvdccaux"]
    display("Phvdcoaux:$Phvdcoaux")
    display("Phvdccaux:$Phvdccaux")
    Pg = resultDC1["solution"]["nw"]["2"]["reserves"]["2"]["Pgg"]/data_cont["nw"]["1"]["load"]["1"]["pd"]
    Pf = resultDC1["solution"]["nw"]["2"]["reserves"]["2"]["Pff"]/data_cont["nw"]["1"]["load"]["1"]["pd"]
    Tcl = data_cont["nw"]["1"]["reserves"]["2"]["Tcl"]
    Td = data_cont["nw"]["1"]["reserves"]["2"]["Td"]
    Tg = data_cont["nw"]["1"]["reserves"]["2"]["Tg"]
    Tf = data_cont["nw"]["1"]["reserves"]["2"]["Tf"]
    display("Pg:$Pg")
    display("Pf:$Pf")
    display(data_cont["nw"]["1"]["load"]["1"]["pd"])
    t4 = (Phvdcoaux -  Phvdccaux - Pf)*Tg/Pg
    display("t4:$t4")
    return t4
end


text:
using Distributed

# to run in parallel set environment variable JULIA_NUM_THREADS to > 1
if haskey(ENV, "JULIA_NUM_THREADS")
    nw = parse(Int, ENV["JULIA_NUM_THREADS"])
    addprocs(nw; exeflags="--project")
end

using VortexCollisions

const nt = 10000

const M = 64
const N = 64

const ℳcut = 0
const 𝒩cut = 0

const a = 1
const b = 2

function u_init_exp_cos(x,y)
    exp(a*cos(x-π) + b*cos(y-π))
end

function u_init_sinx4(x,y)
    sin(x)^4
end

function u_init_sinx4siny4(x,y)
    sin(x)^4 * sin(y)^4
end

function u_init_sinx4siny4etc(x,y)
    sin(x)^4 * sin(y)^4 + 0.2 * sin(2x)^4 * sin(0.5 * y)^4
end

gr = Grid2d(M,N)
ft = FourierTransform(gr, ℳcut=ℳcut, 𝒩cut=𝒩cut)
op = FokkerPlanckOperator(gr, ft)

# test
run_simulation(op, 10, 1E-6, u_init_exp_cos,       "FokkerPlanck_exp_cos_dt1E-6_nt10.h5")
# run_simulation(op, 10, 1E-3, u_init_sinx4,         "FokkerPlanck_sinx4_dt1E-3_nt10.h5")
# run_simulation(op, 10, 1E-3, u_init_sinx4siny4,    "FokkerPlanck_sinx4siny4_dt1E-3_nt10.h5")
# run_simulation(op, 10, 1E-3, u_init_sinx4siny4etc, "FokkerPlanck_sinx4siny4etc_dt1E-3_nt10.h5")

# production
# run_simulation(op, nt, 1E-6, u_init_exp_cos,       "FokkerPlanck_exp_cos_dt1E-6_nt"      * string(nt) * ".h5")
# run_simulation(op, nt, 1e-3, u_init_sinx4,         "FokkerPlanck_sinx4_dt1E-3_nt"        * string(nt) * ".h5")
# run_simulation(op, nt, 1e-3, u_init_sinx4siny4,    "FokkerPlanck_sinx4siny4_dt1E-3_nt"   * string(nt) * ".h5")
# run_simulation(op, nt, 1e-3, u_init_sinx4siny4etc, "FokkerPlanck_sinx4sinyetc_dt1E-3_nt" * string(nt) * ".h5")


text:
include("../../module/CustomStencil.jl")
include("../../misc/misc.jl")
include("../../misc/cpu_stencils.jl")
## Star Stencil definition with radius = 2
# first deriv  4 acccuracy
coefs = [0;2/3;-1/12]
fwd = @def_stencil_expression c[0]D[x,y,z] +
            @sum(i, 1, 2,  c[i]*(D[x+i,y,z] + D[x,y+i,z] + D[x,y,z+i])) +
            @sum(i, 1, 2, -c[i]*(D[x-i,y,z] + D[x,y-i,z] + D[x,y,z-i]))
st_def = CreateStencilDefinition(fwd, coefs)
st_inst = NewStencilInstance(st_def, m_step=false, prev_time_coeff=1)  #don't combine time steps

## Input Data size Definition
radius = 2
nx = 5
ny = 5
nz = 5
dx = 1<<(nx)
dy = 1<<(ny)
dz = 1<<(nz)

data = CreateData(dx,dy,dz)

t_steps = 20
## GPU Implementation
gpu_out = ApplyStencil(st_inst, data, t_steps)

## CPU result
cpu_out = cpu_star_stencil(data, radius, coefs, t_steps)

## Check
mean_er = sum(abs.(cpu_out .- gpu_out))./(reduce(*, size(gpu_out)))
@assert(mean_er < 0.001)

## Visualize Error
## Visualize Error
using Plots

shift = 0
Int(round(dz/2))+shift
heatmap(gpu_out[:,:,Int(round(dz/2))+shift]-cpu_out[:,:,Int(round(dz/2))+shift])
heatmap(cpu_out[:,:,Int(round(dz/2))+shift])
heatmap(gpu_out[:,:,Int(round(dz/2))+shift])


text:
# Labels
include("Tags/Tags.jl")

# Ket
include("Ket.jl")

# Sites
include("Core.jl")
include("Sites.jl")
include("SubSites.jl")

# Space
include("Space/Space.jl")


text:
using QuantumClifford
using Yao
using LinearAlgebra
using Plots

N = 12
M = 1000
sts_clifford = [random_stabilizer(N) for _ = 1:M]
sts_haar = [rand_state(N) for _ = 1:M]
kernel_mat_clifford = ones(M, M)
for i = 1:M, j = i+1:M
    kernel_mat_clifford[i,j] = sts_clifford[i]⋅sts_clifford[j]
    @show i,j
end
for i = 1:M, j = i+1:M
    kernel_mat_clifford[j,i] = kernel_mat_clifford[i,j]
end

kernel_mat_haar = ones(M, M)
for i = 1:M, j = i+1:M
    kernel_mat_haar[i,j] = fidelity(sts_haar[i], sts_haar[j])
    @show i,j
end
for i = 1:M, j = i+1:M
    kernel_mat_haar[j,i] = kernel_mat_haar[i,j]
end

begin
    histogram(eigvals(kernel_mat_haar); label = "haar")
    histogram!(eigvals(kernel_mat_clifford); label = "clifford")
end


text:
#Write a guessing game where the user has to guess a secret number. After every guess the program tells the user whether their number was too large or too small. At the end the number of tries needed should be printed. It counts only as one try if they input the same number multiple times consecutively.

rand_num = rand(1:100)
println(rand_num)

function guess_number()
    println("Please guess a number between 1-100")
    num = parse(Int64,readline())
    return num
end

function give_hint(num, rand_num = rand_num)
    if num > rand_num
        println("Number is less than your guess")
    else num < rand_num
        println("Number is greater than your guess")
    end
end

function is_correct_num(num,rand_num=rand_num)
    corr = 0
    if num == rand_num
        println("Correct!!!")
        corr += 1
        return corr
    else
        println("Please try again!")
        corr += 0
        return corr
    end
end

while true
    num = guess_number()
    give_hint(num)
    is_correct_num(num) == 1 && break
end

max_stars_repo_name: atott/Simple-Programming-Problems-Julia

max_issues_repo_name: atott/Simple-Programming-Problems-Julia

max_forks_repo_name: atott/Simple-Programming-Problems-Julia