text:
export
    set_x0!,
    is_constrained,
    update_problem,
    continuous

"""$(TYPEDEF) Trajectory Optimization Problem.
Contains the full definition of a trajectory optimization problem, including:
* dynamics model (`Model`). Can be either continuous or discrete.
* objective (`Objective`)
* constraints (`Constraints`)
* initial and final states
* Primal variables (state and control trajectories)
* Discretization information: knot points (`N`), time step (`dt`), and total time (`tf`)

# Constructors:
```julia
Problem(model, obj X0, U0; integration, constraints, x0, xf, dt, tf, N)
Problem(model, obj, U0; integration, constraints, x0, xf, dt, tf, N)
Problem(model, obj; integration, constraints, x0, xf, dt, tf, N)
```
# Arguments
* `model`: Dynamics model. Can be either `Discrete` or `Continuous`
* `obj`: Objective
* `X0`: Initial state trajectory. If omitted it will be initialized with NaNs, to be later overwritten by the solver.
* `U0`: Initial control trajectory. If omitted it will be initialized with zeros.
* `x0`: Initial state. Defaults to zeros.
* `xf`: Final state. Defaults to zeros.
* `dt`: Time step
* `tf`: Final time. Set to zero to specify a time penalized problem.
* `N`: Number of knot points. Defaults to 51, unless specified by `dt` and `tf`.
* `integration`: One of the defined integration types to discretize the continuous dynamics model. Defaults to `:none`, which will pass in the continuous dynamics (eg. for DIRCOL)
Both `X0` and `U0` can be either a `Matrix` or a `Vector{Vector}`, but must be the same.
At least 2 of `dt`, `tf`, and `N` need to be specified (or just 1 of `dt` and `tf`).
"""
mutable struct Problem{T<:AbstractFloat,D<:DynamicsType}
    model::AbstractModel
    obj::AbstractObjective
    constraints::Constraints
    x0::Vector{T}
    xf::Vector{T}
    X::VectorTrajectory{T}
    U::VectorTrajectory{T}
    N::Int
    dt::T
    tf::T

    function Problem(model::Model{M,D}, obj::AbstractObjective, constraints::Constraints,
        x0::Vector{T},xf::Vector{T}, X::VectorTrajectory, U::VectorTrajectory, N::Int, dt::Real, tf::Real) where {M,T,D}

        n,m = model.n, model.m
        # TODO these checks break for infeasible, minimum time -> do a post check
        # @assert length(X[1]) == n
        # @assert length(U[1]) == m
        # @assert length(x0) == n
        @assert length(obj) == N
        @assert length(constraints) == N
        @assert length(X) == N

        # if length(U) == N
        #     @info "Length of U should be N-1, not N. Trimming last entry"
        #     U = U[1:end-1]
        # end
        # @assert length(U) == N-1
        if !(dt > 0)
            throw(ArgumentError("dt must be strictly positive"))
        end

        new{T,D}(model,obj,constraints,x0,xf,X,U,N,dt,tf)
    end
end

function Problem(model::Model{M,Continuous}, obj::AbstractObjective, X0::VectorTrajectory{T}, U0::VectorTrajectory{T};
        N::Int=length(obj), constraints::Constraints=Constraints(N),
        x0::Vector{T}=zeros(model.n), xf::Vector{T}=zeros(model.n),
        dt=NaN, tf=NaN, integration=:none) where {M,T}
    N, tf, dt = _validate_time(N, tf, dt)
    if integration == :none
        model = model
    elseif isdefined(TrajectoryOptimization,integration)
        discretizer = eval(integration)
        model = discretizer(model)
    else
        throw(ArgumentError("$integration is not a defined integration scheme"))
    end
    Problem(model, obj, constraints, x0, xf, deepcopy(X0), deepcopy(U0), N, dt, tf)
end

function Problem(model::Model{M,Discrete}, obj::AbstractObjective, X0::VectorTrajectory{T}, U0::VectorTrajectory{T};
        N::Int=length(obj), constraints::Constraints=Constraints(N), x0::Vector{T}=zeros(model.n),
        xf::Vector{T}=zeros(model.n), dt=NaN, tf=NaN) where {M,T}
    N, tf, dt = _validate_time(N, tf, dt)
    Problem(model, obj, constraints, x0, xf, deepcopy(X0), deepcopy(U0), N, dt, tf)
end

Problem(model::Model, obj::Objective, X0::Matrix{T}, U0::Matrix{T}; kwargs...) where T =
    Problem(model, obj, to_dvecs(X0), to_dvecs(U0); kwargs...)

function Problem(model::Model, obj::AbstractObjective, U0::VectorTrajectory{T}; kwargs...) where T
    N = length(obj)
    X0 = empty_state(model.n, N)
    Problem(model, obj, X0, U0; kwargs...)
end
Problem(model::Model, obj::AbstractObjective, U0::Matrix{T}; kwargs...) where T =
    Problem(model, obj, to_dvecs(U0); kwargs...)

function Problem(model::Model, obj::AbstractObjective; kwargs...)
    N = length(obj)
    U0 = [zeros(model.m) for k = 1:N-1]
    X0 = empty_state(model.n, N)
    Problem(model, obj, X0, U0; kwargs...)
end

# QUESTION: Remove these constructors?
"$(SIGNATURES) Pass in a cost instead of an objective"
function Problem(model::Model{Nominal,Discrete}, cost::CostFunction, U0::VectorTrajectory{T}; kwargs...) where T
    N = length(U0) + 1
    obj = Objective(cost, N)
    Problem(model, obj, U0; kwargs...)
end

Problem(model::Model{Nominal,Discrete}, cost::CostFunction, U0::Matrix{T}; kwargs...) where T =
    Problem(model, cost, to_dvecs(U0); kwargs...)

"""```julia
update_problem(prob; kwargs...)
```
Create a new problem from another, specifing all fields as keyword argumen# if length(U) == N
        #     @info "Length of U should be N-1, not N. Trimming last entry"
        #     U = U[1:end-1]
        # end
        # @assert length(U) == N-1ts
The `newProb` argument can be set to true if a the primal variables are to be copied, otherwise they will be passed to the modified problem.
"""
function update_problem(p::Problem;
    model=p.model, obj=p.obj, constraints=p.constraints, x0=p.x0, xf=p.xf, X=p.X, U=p.U,
    N=p.N, dt=p.dt, tf=p.tf, newProb=true)

    if newProb
        Problem(model,obj,constraints,x0,xf,deepcopy(X),deepcopy(U),N,dt,tf)
    else
        Problem(model,obj,constraints,x0,xf,X,U,N,dt,tf)
    end
end

"$(SIGNATURES) Set the initial control trajectory for a problem. U0 can be either a `Matrix` or a `Vector{Vector}`"
initial_controls!(prob::Problem{T}, U0::AbstractVectorTrajectory{T}) where T = copyto!(prob.U, U0[1:prob.N-1])
initial_controls!(prob::Problem{T}, U0::Matrix{T}) where T = initial_controls!(prob, to_dvecs(U0))

"$(SIGNATURES) Set the initial state trajectory for a problem. X0 can be either a `Matrix` or a `Vector{Vector}`"
initial_states!(prob::Problem{T}, X0::AbstractVectorTrajectory{T}) where T = copyto!(prob.X, X0)
initial_states!(prob::Problem{T}, X0::Matrix{T}) where T = initial_state!(prob, to_dvecs(X0))

"$(SIGNATURES) Set the initial state"
set_x0!(prob::Problem{T}, x0::Vector{T}) where T = copyto!(prob.x0, x0)

# TODO: think about how to do this properly now that objective and constraints depend on N
# function change_N(prob::Problem, N::Int)
#     tf = final_time(prob)
#     dt = tf/(N-1)
#     X, U = interp_traj(N, tf, prob.X, prob.U)
#     @show length(X)
#     Problem(prob.model, prob.obj, prob.constraints, prob.x0, X, U, N, dt, tf)
# end

function _validate_time(N,tf,dt)
    if N == -1 && isnan(dt) && isnan(tf)
        throw(ArgumentError("Must specify at least 2: N, dt, or tf"))
    end
    # Check for minimum time
    if (tf == 0) || (tf == :min)
        tf = 0.0
        if N==-1
            throw(ArgumentError("N must be specified for a minimum-time problem"))
        end
    elseif tf > 0
        # Handle combination of N and dt
        if isnan(dt) && N > 0
            dt = tf / (N-1)
        elseif ~isnan(dt) && N == -1
            N, dt = calc_N(tf, dt)
        elseif isnan(dt) && N==-1
            @warn "Neither dt or N were specified. Setting N = 51"
            N = 51
            dt = tf / (N-1)
        elseif ~isnan(dt) && N > 0
            if dt !== tf/(N-1)
                throw(ArgumentError("Specified time step, number of knot points, and final time do not agree ($dt ≢ $(obj.tf)/$(N-1))"))
            end
        end
        if dt == 0
            throw(ArgumentError("dt must be non-zero for non-minimum time problems"))
        end
    elseif isnan(tf)
        if dt > 0
            if N == -1
                @warn "Neither tf or N were specified. Setting N = 51"
                N = 51
            end
            tf = dt*(N-1)
        else
            throw(ArgumentError("dt must be positive for a non-minimum-time problem"))
        end
    else
        throw(ArgumentError("Invalid input for tf"))
    end

    # Check N, dt for valid entries
    if N < 0
        err = ArgumentError("$N is not a valid entry for N. Number of knot points must be a positive integer.")
        throw(err)
    elseif dt < 0
        err = ArgumentError("$dt is not a valid entry for dt. Time step must be positive.")
        throw(err)
    end
    return N,tf,dt
end

"""```
n,m,N = size(p::Problem)
```
Return the number of states (n), number of controls (m), and the number of knot points (N)"""
Base.size(p::Problem)::NTuple{3,Int} = (p.model.n,p.model.m,p.N)

"$(TYPEDSIGNATURES) Copy a problem"
Base.copy(p::Problem) = Problem(p.model, p.obj, copy(p.constraints), copy(p.x0), copy(p.xf),
    deepcopy(p.X), deepcopy(p.U), p.N, p.dt, p.tf)

empty_state(n::Int,N::Int) = [ones(n)*NaN32 for k = 1:N]

"$(SIGNATURES) Checks if a problem has any constraints"
is_constrained(prob::Problem{T}) where T = !all(isempty.(prob.constraints.C))

TrajectoryOptimization.num_constraints(prob::Problem) = num_constraints(prob.constraints)

"$(TYPEDSIGNATURES) Evaluate the current cost for the problem"
cost(prob::Problem{T}) where T = cost(prob.obj, prob.X, prob.U, get_dt_traj(prob))::T

function max_violation(prob::Problem{T})::T where T
    if is_constrained(prob)
        N = prob.N
        c_max = 0.0
        for k = 1:N-1
            if num_stage_constraints(prob.constraints[k]) > 0
                stage_con = stage(prob.constraints[k])
                c = PartedVector(stage_con)
                evaluate!(c,stage_con,prob.X[k],prob.U[k])
                max_E = norm(c.equality,Inf)
                max_I = maximum(pos.(c))
                c_max = max(c_max,max(max_E,max_I))
            end
        end
        if num_terminal_constraints(prob.constraints[N]) > 0
            c = PartedVector(prob.constraints[N],:terminal)
            evaluate!(c,prob.constraints[N],prob.X[N])
            max_E = norm(c.equality,Inf)
            max_I = maximum(pos.(c))
            c_max = max(c_max,max(max_E,max_I))
        end
        return c_max
    else
        return 0.
    end
end

function Expansion(prob::Problem{T}) where T
    n = prob.model.n; m = prob.model.m
    Expansion(zeros(T,n),zeros(T,m),zeros(T,n,n),zeros(T,m,m),zeros(T,m,n))
end

function Expansion(prob::Problem{T},exp::Symbol) where T
    n = prob.model.n; m = prob.model.m
    if exp == :x
        return Expansion(zeros(T,n),zeros(T,0),zeros(T,n,n),zeros(T,0,0),zeros(T,0,0))
    elseif exp == :u
        return Expansion(zeros(T,0),zeros(T,m),zeros(T,0,0),zeros(T,m,m),zeros(T,0,0))
    else
        error("Invalid expansion components requested")
    end
end

midpoint(prob::Problem{T,Continuous}) where T = update_problem(prob, model=midpoint(prob.model))
rk3(prob::Problem{T,Continuous}) where T = update_problem(prob, model=rk3(prob.model))
rk4(prob::Problem{T,Continuous}) where T = update_problem(prob, model=rk4(prob.model))

continuous(prob::Problem{T,Discrete}) where T<:AbstractFloat = update_problem(prob, model=continuous(prob.model))

"Return timestep"
function get_dt(prob::Problem,k::Int)
    if prob.tf == 0.
        return prob.U[k][end]^2
    else
        return prob.dt
    end
end

function get_dt(prob::Problem,U::AbstractArray)
    if prob.tf == 0.
        return U[end]^2
    else
        return prob.dt
    end
end

function get_dt_traj(prob::Problem)
    [get_dt(prob,k) for k = 1:prob.N-1]
end

function get_dt_traj(prob::Problem,U::Trajectory)
    [get_dt(prob,U[k]) for k = 1:prob.N-1]
end

"$(SIGNATURES) Get the total time for tje trajectory (applicable for time-penalized problems)"
final_time(prob::Problem) = sum(get_dt_traj(prob))

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
using FractalFlow, StaticArrays, PyCall, Polynomials
using FractalFlow: BoundingBox, Segment, BoxCounter

@pyimport matplotlib.pyplot as plt
@pyimport matplotlib.collections as plt_col
@pyimport matplotlib.patches as plt_patch

"""Plot list of lines as line collection"""
function plot{T}(segments::Vector{Segment{2, T}}, ax)
    segments = [[seg.A, seg.B] for seg in segments]
    lines = plt_col.LineCollection(segments, colors="k")
    py"$(ax).add_collection($(lines))"
    ax
end

"""Plot outline of bounding boxes"""
function plot{T}(boxes::Vector{BoundingBox{2, T}}, ax)
    boxes = [plt_patch.Rectangle(box.low, (box.high - box.low)...)
             for box in boxes]
    boxes = plt_col.PatchCollection(boxes, edgecolor="b", facecolor="none")
    py"$(ax).add_collection($(boxes))"
    ax
end

volume{D, T}(p::SVector{D, T}) = (prod(p))^(1./D)

function fd_estimation{T}(c::Curve{2, T}, nlevel::Int)
    counts = []
    sizes = []

    bc = BoxCounter(cc)
    state = start(bc)
    @show state

    fig = plt.figure()
    ax = py"$(fig).gca()"

    plot(bc.segments, ax)

    boxes = map(last, bc.seg_box)
    # Extract dim
    min_ = first(boxes).low - 0.1*ones(2)
    max_ = first(boxes).high + 0.1*ones(2)

    xlim = map(first, (min_, max_))
    ylim = map(last, (min_, max_))

    plot(boxes, ax)
    plt.xlim(xlim)
    plt.ylim(ylim)
    plt.axis("off")

    for level in 1:7
        state = next(bc, state)

        fig = plt.figure()
        ax = py"$(fig).gca()"

        plot(bc.segments, ax)

        boxes = map(last, bc.seg_box)
        plot(boxes, ax)

        plt.xlim(xlim)
        plt.ylim(ylim)
        plt.axis("off")

        @show state
        push!(counts, first(state))
        push!(sizes, volume(last(state)))
    end

    fig = plt.figure()
    plt.loglog(sizes, counts, "bx-")
    plt.loglog(sizes, sizes.^(-1.)*(last(counts)/(last(sizes)^(-1.))), "r")

    (counts, sizes)
end

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

# segments = [(SVector(0., 0.), SVector(1., 0)),
#             (SVector(1., 0.), SVector(1., -1)),
#             (SVector(1., -1.), SVector(0., -1)),
#             (SVector(0., -1.), SVector(0., 0)),
#             (SVector(0., 0.), SVector(-1., 0)),
#             (SVector(-1., 0.), SVector(-1., 1)),
#             (SVector(-1., 1.), SVector(0., 1)),
#             (SVector(0., 1.), SVector(0., 0)),
#             (SVector(-1., 0), SVector(0., 1)),
#             (SVector(0., -1.), SVector(1., 0))]

# cc = Curve(segments)

cc = koch_q1_curve(6)

(counts, sizes) = fd_estimation(cc, 8)

# Use the last 3 for eastimation
fit = Polynomials.polyfit(log.(counts[end-3:end]), log.(sizes[end-3:end]), 1)
dim = -last(fit.a)
println("Estimated fractal dim $(dim)")

# FIXME: Would be fun to spin this in parallel
dims = [first(FractalFlow.fractal_dim(cc, 10)) for _ in 1:10]
@show dims

dim = mean(dims)
dim_std = std(dims)
@show (dim, dim_std)

ext: jl


text:
@testset "MPI" begin
    db = setup_db("test/data/TEST_PowerFlowNetworks_SQLite.sqlite"; delete_if_exists=true)
    load_instance_in_db!(db, "case6ww", 0, "test/data/case6ww.m", "MATPOWERM", Dates.now())
    load_instance_in_db!(db, "C2FEN02312", 3, "test/data/C2FEN02312/scenario_3/", "RAWGO", Dates.now())
    load_instance_in_db!(db, "C2S6N02045", 1, "test/data/C2S6N02045/scenario_1/", "RAWGO", Dates.now())

    execute_process_mpi(db, "save_basic_features_instances", "test/log/")
    execute_process_mpi(db, "save_features_instances", "test/log/")
    execute_process_mpi(db, "serialize_instances", "test/log/"; serialize_path="test/data/networks_serialized", graphs_path="test/data/graphs")
    execute_process_mpi(db, "generate_decompositions", "test/log";
                        cliques_path="test/data/cliques",
                        cliquetrees_path="test/data/cliquetrees",
                        graphs_path="test/data/graphs",
                        extension_alg="minimum_degree",
                        preprocess_path="test/data/configs/preprocess_add_random_50.json")
    execute_process_mpi(db, "merge_decompositions", "test/log";
                        heuristic=["molzahn"],
                        heuristic_switch=[0],
                        treshold_name="clique_nv_up",
                        merge_kwargs=Dict("treshold_percent" => 0.5))
    execute_process_mpi(db, "combine_decompositions", "test/log";
                        how="clique_intersect", extension_alg="cholesky")
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
using Flux, Flux.Data.MNIST, Statistics, BSON, Random, StatsPlots; pyplot()
using Flux: onehotbatch, onecold, crossentropy, @epochs

epochs = 30
eta = 1e-3
batchSize = 200
trainRange, validateRange = 1:1000, 1001:5000

function minibatch(x, y, idxs)
    xBatch = Array{Float32}(undef, size(x[1])..., 1, length(idxs))
    for i in 1:length(idxs)
        xBatch[:, :, :, i] = Float32.(x[idxs[i]])
    end
    return (xBatch, onehotbatch(y[idxs], 0:9))
end

trainLabels = MNIST.labels()[trainRange]
trainImgs = MNIST.images()[trainRange]
mbIdxs = Iterators.partition(1:length(trainImgs), batchSize)
trainSet = [minibatch(trainImgs, trainLabels, i) for i in mbIdxs]

validateLabels = MNIST.labels()[validateRange]
validateImgs = MNIST.images()[validateRange]
validateSet = minibatch(validateImgs, validateLabels, 1:length(validateImgs))

accuracy(x, y, model) = mean(onecold(model(x)) .== onecold(y))
loss(x, y, model) = crossentropy(model(x), y)

function evalAccuracy(dropP)
    model= Chain(Conv((5, 5), 1=>8, relu), MaxPool((2,2)),
                    Conv((3, 3), 8=>16, relu), MaxPool((2,2)),
                    flatten,
                    Dense(400, 200,relu), Dropout(dropP),
                    Dense(200, 200,relu), Dropout(dropP),
                    Dense(200, 200,relu), Dropout(dropP),
                    Dense(200, 10,relu), Dropout(dropP),
                    softmax)
    opt = ADAM(eta);
    @epochs epochs Flux.train!((x,y)->loss(x,y,model),params(model),trainSet,opt)
    accuracy(validateSet..., model)
end

pToTest = [0.0, 0.25, 0.5, 0.75]
n = 10
results = [[evalAccuracy(p) for _ in 1:n] for p in pToTest]
bestAcc, bestI = findmax(median.(results))
println("The best dropout probability is $(pToTest[bestI]).")
println("It achieves $(bestAcc) accuracy on average.")

boxplot(results,label="",
       xticks=([1:1:4;],string.(pToTest)),
       xlabel="Dropout Probability", ylabel = "Accuracy", legend = false)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
"""    ABODDetector(n_neighbors = 5,
                      method = "fast")
$(make_docs_link("abod"))
"""
@pymodel mutable struct ABODDetector <: UnsupervisedDetector
    n_neighbors::Integer = 5::(_ > 0)
    method::String = "fast"::(_ in ("default", "fast"))
end

"""    CBLOFDetector(n_clusters = 8,
                       alpha = 0.9,
                       beta = 5,
                       use_weights = false,
                       random_state = nothing,
                       n_jobs = 1)
$(make_docs_link("cblof"))
"""
@pymodel mutable struct CBLOFDetector <: UnsupervisedDetector
    n_clusters::Integer = 8::(_ > 1)
    alpha::Real = 0.9::(0.5 < _ < 1)
    beta::Real = 5::(_ > 1)
    use_weights::Bool = false
    random_state::Union{Nothing, Integer} = nothing
    n_jobs::Integer = 1::(_ >= -1)
end

"""    COFDetector(n_neighbors = 5)
$(make_docs_link("cof"))
"""
@pymodel mutable struct COFDetector <: UnsupervisedDetector
    n_neighbors::Integer = 5::(_ > 0)
end

"""    COPODDetector()
$(make_docs_link("copod"))
"""
@pymodel mutable struct COPODDetector <: UnsupervisedDetector end

"""    HBOSDetector(n_bins = 10,
              alpha = 0.1,
              tol = 0.5)
$(make_docs_link("hbos"))
"""
@pymodel mutable struct HBOSDetector <: UnsupervisedDetector
    n_bins::Integer = 10::(_ > 1)
    alpha::Real = 0.1::(0 < _ < 1)
    tol::Real = 0.5::(0 < _ < 1)
end

"""    IForestDetector(n_estimators = 100,
                         max_samples = "auto",
                         max_features = 1.0
                         bootstrap = false,
                         behaviour = "new",
                         random_state = nothing,
                         verbose = 0,
                         n_jobs = 1)
$(make_docs_link("iforest"))
"""
@pymodel mutable struct IForestDetector <: UnsupervisedDetector
    n_estimators::Integer = 100::(_ > 0)
    max_samples::Union{String, Real} = "auto"
    max_features::Real = 1.0
    bootstrap::Bool = false
    behaviour::String = "new"
    random_state::Union{Nothing, Integer} = nothing
    verbose::Integer = 0::(0 <= _ <= 2)
    n_jobs::Integer = 1::(_ >= -1)
end

"""    KNNDetector(n_neighbors = 5,
             method = "largest",
             radius = 1.0,
             algorithm = "auto",
             leaf_size = 30,
             metric = "minkowski",
             p = 2,
             metric_params = nothing,
             n_jobs = 1)
$(make_docs_link("knn"))
"""
@pymodel mutable struct KNNDetector <: UnsupervisedDetector
    n_neighbors::Integer = 5::(_ > 0)
    method::String = "largest"::(_ in ("largest", "mean", "median"))
    radius::Real = 1.0
    algorithm::String = "auto"::(_ in ("auto", "ball_tree", "kd_tree", "brute"))
    leaf_size::Integer = 30::(_ > 0)
    metric::String = "minkowski"::(_ in ("cityblock", "cosine", "euclidean", "l1", "l2", "manhatten", "braycurtis", "canberra", "chebyshev", "correlation", "dice", "hamming", "jaccard", "kulsinski", "mahalanobis", "matching", "minkowski", "rogerstanimoto", "russellrao", "seuclidean", "sokalmichener", "sokalsneath", "sqeuclidean", "yule"))
    p::Union{Nothing, Integer} = 2
    metric_params::Union{Nothing, Any} = nothing
    n_jobs::Integer = 1::(_ >= -1)
end

"""    LMDDDetector(n_iter = 50,
                      dis_measure = "aad",
                      random_state = nothing)
$(make_docs_link("lmdd"))
"""
@pymodel mutable struct LMDDDetector <: UnsupervisedDetector
    n_iter::Integer = 50::(_ > 0)
    dis_measure::String = "aad"::(_ in ("aad", "var", "iqr"))
    random_state::Union{Nothing, Integer} = nothing
end

"""    LODADetector(n_bins = 10,
                      n_random_cuts = 100)
$(make_docs_link("loda"))
"""
@pymodel mutable struct LODADetector <: UnsupervisedDetector
    n_bins::Integer = 10::(_ > 1)
    n_random_cuts::Integer = 100::(_ > 0)
end

"""    LOFDetector(n_neighbors = 5,
                     method = "largest",
                     algorithm = "auto",
                     leaf_size = 30,
                     metric = "minkowski",
                     p = 2,
                     metric_params = nothing,
                     n_jobs = 1)
$(make_docs_link("lof"))
"""
@pymodel mutable struct LOFDetector <: UnsupervisedDetector
    n_neighbors::Integer = 5::(_ > 0)
    algorithm::String = "auto"::(_ in ("auto", "ball_tree", "kd_tree", "brute"))
    leaf_size::Integer = 30::(_ > 0)
    metric::String = "minkowski"::(_ in ("cityblock", "cosine", "euclidean", "l1", "l2", "manhatten", "braycurtis", "canberra", "chebyshev", "correlation", "dice", "hamming", "jaccard", "kulsinski", "mahalanobis", "matching", "minkowski", "rogerstanimoto", "russellrao", "seuclidean", "sokalmichener", "sokalsneath", "sqeuclidean", "yule"))
    p::Union{Nothing, Integer} = 2
    metric_params::Union{Nothing, Any} = nothing
    n_jobs::Integer = 1::(_ >= -1)
end

"""    LOCIDetector(alpha = 0.5,
                      k = 3)
$(make_docs_link("loci"))
"""
@pymodel mutable struct LOCIDetector <: UnsupervisedDetector
    alpha::Real = 0.5::(0 < _ < 1)
    k::Real = 3::(_ > 0)
end

"""    MCDDetector(store_precision = true,
                     assume_centered = false,
                     support_fraction = nothing,
                     random_state = nothing)
$(make_docs_link("mcd"))
"""
@pymodel mutable struct MCDDetector <: UnsupervisedDetector
    store_precision::Bool = true
    assume_centered::Bool = false
    support_fraction::Union{Nothing, Real} = nothing
    random_state::Union{Nothing, Integer} = nothing
end

"""    OCSVMDetector(kernel = "rbf",
                       degree = 3,
                       gamma = "auto",
                       coef0 = 0.0,
                       tol = 0.001,
                       nu = 0.5,
                       shrinking = true,
                       cache_size = 200,
                       verbose = false,
                       max_iter = -1)
$(make_docs_link("ocsvm"))
"""
@pymodel mutable struct OCSVMDetector <: UnsupervisedDetector
    kernel::String = "rbf"::(_ in ("linear", "poly", "rbf", "sigmoid", "precomputed"))
    degree::Integer = 3::(_ > 1)
    gamma::Union{String, Real} = "auto"
    coef0::Real = 0.0
    tol::Real = 0.001
    nu::Real = 0.5::(0 < _ <= 1)
    shrinking::Bool = true
    cache_size::Integer = 200::(_ > 0)
    verbose::Bool = false
    max_iter::Integer = -1
end

"""    PCADetector(n_components = nothing,
                     n_selected_components = nothing,
                     copy = true,
                     whiten = false,
                     svd_solver = "auto",
                     tol = 0.0
                     iterated_power = "auto",
                     standardization = true,
                     weighted = true,
                     random_state = nothing)
$(make_docs_link("pca"))
"""
@pymodel mutable struct PCADetector <: UnsupervisedDetector
    n_components::Union{Nothing, Real} = nothing
    n_selected_components::Union{Nothing, Integer} = nothing
    copy::Bool = true
    whiten::Bool = false
    svd_solver::String = "auto"::(_ in ("auto", "full", "arpack", "randomized"))
    tol::Real = 0.0
    iterated_power::Union{String, Integer} = "auto"
    standardization::Bool = true
    weighted::Bool = true
    random_state::Union{Nothing, Integer} = nothing
end

"""    RODDetector(parallel_execution = false)
$(make_docs_link("rod"))
"""
@pymodel mutable struct RODDetector <: UnsupervisedDetector
    parallel_execution::Bool = false
end

"""    SODDetector(n_neighbors = 5,
                     ref_set = 10,
                     alpha = 0.8)
$(make_docs_link("sod"))
"""
@pymodel mutable struct SODDetector <: UnsupervisedDetector
    n_neighbors::Integer = 20::(_ > 0)
    ref_set::Integer = 10::(_ > 0)
    alpha::Real = 0.8::(0 < _ < 1)
end

"""    SOSDetector(perplexity = 4.5,
             metric = "minkowski",
             eps = 1e-5)
$(make_docs_link("sos"))
"""
@pymodel mutable struct SOSDetector <: UnsupervisedDetector
    perplexity::Real = 4.5::(_ > 0)
    metric::String = "minkowski"::(_ in ("cityblock", "cosine", "euclidean", "l1", "l2", "manhatten", "braycurtis", "canberra", "chebyshev", "correlation", "dice", "hamming", "jaccard", "kulsinski", "mahalanobis", "matching", "minkowski", "rogerstanimoto", "russellrao", "seuclidean", "sokalmichener", "sokalsneath", "sqeuclidean", "yule"))
    eps::Real = 1e-5::(_ > 0)
end


text:
# https://adventofcode.com/2020/day/1

# find the two entries that sum to 2020 and then multiply those two numbers together

function part1(list_of_numbers, total=2020)
    for (i, x) in enumerate(list_of_numbers)
        y = total - x
        if y in list_of_numbers[(i + 1):end]
            return x * y
        end
    end
end

example = [1721, 979, 366, 299, 675, 1456]
@assert part1(example) == 514579

# Using the above example again, the three entries that sum to 2020 are
# 979, 366, and 675. Multiplying them together produces the answer, 241861950.

function part2(list_of_numbers)
    for (i, x) in enumerate(list_of_numbers)
        partial = 2020 - x
        tmp = filter(y -> y <= partial, list_of_numbers[(i + 1):end])
        y = part1(tmp, partial)
        if !isnothing(y)
            return x * y
        end
    end
end

@assert part2(example) == 241861950

function read_array(string::AbstractString)::Vector{Int}
    return parse.(Int, split(string, '\n', keepempty=false))
end

test = read_array(read("data/day-01.txt", String))
println("Part 1: $(result1 = part1(test))")
println("Part 2: $(result2 = part2(test))")

@assert result1 == 326211
@assert result2 == 131347190

max_stars_repo_name: twolodzko/advent-of-code

max_issues_repo_name: twolodzko/advent-of-code

max_forks_repo_name: twolodzko/advent-of-code


text:
include("helium_hf_fem_scf.jl")
using Printf
using .Helium_HF_FEM_SCF

const RESULT_FILENAME = "result.csv"

function main()
    res = Helium_HF_FEM_SCF.scfloop()
    if res != nothing
        hfem_param, hfem_val, energy = res
        Helium_HF_FEM_SCF.save_result(hfem_param, hfem_val, RESULT_FILENAME)
        @printf "SCF計算が収束しました: energy = %.14f (Hartree)、計算結果を%sに書き込みました\n" energy RESULT_FILENAME
    else
        @printf "SCF計算が収束しませんでした\n"
    end
end

@time main()


text:
import Pkg
Pkg.activate(temp=true)
Pkg.add("GLMakie")

using Printf
using GLMakie

# Slider

fig = Figure()
angle_selector = Slider(fig[2, 1], range = 0:0.01:2pi, startvalue = pi/2)
angles = lift(angle -> range(0, angle, length=100), angle_selector.value)
points = @lift([Point2(cos(angle), sin(angle))  for angle in $angles])
lines(fig[1,1], points, axis=(aspect=DataAspect(), limits=(-1.5, 1.5, -1.5, 1.5)))

# Slider with label

fig = Figure()
_angle_label(angle) = @sprintf "%1.2fπ" angle/pi
angle_selector = labelslider!(fig, "angle", 0:0.01:2pi, format=_angle_label)
set_close_to!(angle_selector.slider, pi/2)
fig[2, 1] = angle_selector.layout
angles = lift(angle -> range(0, angle, length=100), angle_selector.slider.value)
points = @lift([Point2(cos(angle), sin(angle))  for angle in $angles])
lines(fig[1,1], points, axis=(aspect=DataAspect(), limits=(-1.5, 1.5, -1.5, 1.5)))

# Mouse interaction

fig = Figure()
axs = Axis(fig[1,1], aspect=DataAspect(), limits=(-1.5, 1.5, -1.5, 1.5))
angle = Observable(pi/2)
points = @lift([Point2(cos(a), sin(a)) for a in range(0, $angle, length=100)])
lines!(points)

fig_events = events(axs.scene)
on(fig_events.mouseposition) do event
    x, y = mouseposition(axs.scene)
    println("$x, $y")
    angle[] = y > 0 ? atan(y, x) : 2pi + atan(y, x)
end


text:
struct StandardCubicPolynomial <: AbstractParameters end

Base.@kwdef struct GeneralizedCubicPolynomial{T} <: AbstractParameters
    CS0 :: T = 0.0
    CS1 :: T = 1.0
end

shape(d, p::StandardCubicPolynomial) = d * (1-d)
shape(d, p::GeneralizedCubicPolynomial) = d * (1-d) * ( p.CS0 + p.CS1*(1-d) )


text:
module MangalTestNodes
    using Mangal
    using Test
    include("setup.jl")

    @test typeof(nodes()) <: Vector{MangalNode}

    n1 = node(2158)
    @test typeof(n1) <: MangalNode
    @test typeof(n1.taxon) <: MangalReferenceTaxon

    b1 = backbone(1)
    nb1 = nodes(b1)
    @test first(nb1).taxon == b1

end


text:
"""
Smooth jump function ``\\Phi_{\\omega, T}(t)``

    Returns a planned trajectory.

    smoothjump(t :: Real, T:: Real; ω = 2.0 :: Real)

`t` : actual time

`T` : final time

`ω` : steepness factor of trajectory
"""
function smoothjump(t :: Real, T:: Real; ω = 2.0 :: Real)

    if t <= 0
        return  0.0
    elseif t >= T
        return 1.0
    else

        numer, err1 = quadgk(x-> pulsebell(x, T, ω=ω), 0, t)
        denom, err2 = quadgk(x-> pulsebell(x, T, ω=ω), 0, T)

        return numer/denom
    end
end

"""
Bell-shaped pulse function ``\\Omega_{\\omega, T}(t)``

`t` : actual time

`T` : final time

`ω` : steepness factor of trajectory
"""
function pulsebell(t :: Real, T:: Real; ω = 2.0 :: Real)

    if t < 0 || t > T
        return 0
    else
        exp( -1 / ( (1 - t/T)*(t/T))^ω )
    end

end


text:
include("graph.jl")
include("opinion.jl")

function row_sum_max(M::Matrix{WEIGHT})
    (n, m) = size(M)
    acc = fill(0.0, n)
    for i = 1:n
        for j = 1:n
            acc[i] += abs(M[i, j])
        end
    end
    return maximum(acc)
end

function main()
    dataset = "cloister"
    filepath = "./data/raw/$(dataset).edgelist"
    @show filepath
    (W, U, n) = read_edgelist(filepath)
    m = 2
    normalize_edgelist!(W)
    normalize_edgelist!(U)
    y = generate_opinions(n, m)
    α, β, γ = generate_factors(n)
    M = (β .* W) - (γ .* U)
    Q = get_fundamental_matrix(β, W, γ, U)
    z = get_equlibrium(α, Q, y)
    z_next = get_next_opinion(α,
                              β,
                              γ,
                              W, U, y,z)
    M = Matrix(β .* W - γ .* U)

    y_trans = y'
    z_trans = z'

    p = scatter(y_trans[1,:], y_trans[2,:], fmt = :png)
    savefig("before_equlilbrium_$(dataset).png")
    p = scatter(z_trans[1,:], z_trans[2,:], fmt = :png)
    savefig("after_equlilbrium_$(dataset).png")

    return y, Q, z
end

y, Q, z = main()

df = DataFrame(X = Float64.(y[:,1]), Y = Float64.(y[:, 2]))

df1 = DataFrame(X = Float64.(z[:,1]), Y = Float64.(z[:, 2]))

max_stars_repo_name: InzamamRahaman/SignedDegroot

max_issues_repo_name: InzamamRahaman/SignedDegroot

max_forks_repo_name: InzamamRahaman/SignedDegroot


text:
function parse_map(fn::String)::Vector{Asteroid}
    map = open(fn) do f
        readlines(f)
    end

    asteroids = []

    for (y, i) in enumerate(map)
        for (x, j) in enumerate(i)
            if map[y][x] == '#'
                push!(asteroids, Asteroid([x,y]))
            end
        end
    end

    return asteroids
end

struct Station
    coordinates::Vector{Int}
    asteroids::Int
end

mutable struct Asteroid
    coordinates::Vector{Int}
    distance::Union{Float64, Nothing}
    angle::Union{Float64, Nothing}
end
Asteroid(c) = Asteroid(c, nothing, nothing)

Base.show(io::IO, s::Station) = print(io, "Station at (",s.coordinates[1],",",s.coordinates[2],") with ",s.asteroids," asteroids visible")
Base.show(io::IO, s::Asteroid) = print(io, "Asteroid at (",s.coordinates[1],",",s.coordinates[2],"), angle ",round(rad2deg(s.angle), digits=2),", distance ",round(s.distance, digits=2))

function find_best_station(asteroids::Vector{Asteroid})::Station
    candidate = Station([0,0], 0)

    for a in asteroids
        angles = Set()
        others = filter(x -> x != a, asteroids)

        for o in others
            v = o.coordinates - a.coordinates
            angle = atan(v[1],v[2]) - atan(1,0)
            len = sqrt(v[1]^2+v[2]^2)
            # println(Asteroid(o, len, angle))
            push!(angles, angle)
        end

        if length(angles) > candidate.asteroids
            candidate = Station(a.coordinates, length(angles))
        end
        # println(sort(collect(Float64, angles)))
        # println(candidate)
    end

    return candidate
end

function nuke(asteroids::Vector{Asteroid}, station::Station)
    laser_direction = pi/2 # -> 0deg, ccw
    others = filter(x -> x.coordinates != station.coordinates, asteroids)
    for o in others
        v = o.coordinates - station.coordinates
        # angle = acos(dot(x,v)/(norm(x)*norm(v)))
        o.angle = atan(v[1],v[2]) - atan(1,0)
        if o.angle < 0
            o.angle = o.angle + 2*pi
        end
        o.distance = sqrt(v[1]^2+v[2]^2)
        # a = Asteroid(o.coordinates, len, angle)
        # println(o)
    end
    counter = 0
    # laser_direction = 0
    println("Targets: ", length(others))
    while length(others) > 0 # && counter < 100
        # println("laser direction ", round(laser_direction, digits=2))
        targets = sort(filter(x -> x.angle == laser_direction, others), by = x -> x.distance)
        # println("Targets: ", targets)
        if length(targets) > 0
            target = targets[1]
            filter!(x -> x != target, others)
            counter = counter + 1
            println("Target ",counter,": ", target)
            # println(length(others)," targets left")
            if length(others) == 0
                break
            end
        else
            println("No target?!")
        end
        angles = sort(filter(x -> x < laser_direction, map(x -> x.angle, others)), rev = true)
        if length(angles) == 0
            laser_direction = laser_direction + 2*pi
            angles = sort(filter(x -> x < laser_direction, map(x -> x.angle, others)), rev = true)
        end
        # println(angles)
        laser_direction = angles[1]
        # break
    end
end

asteroids = parse_map("asteroids.txt")
station = find_best_station(asteroids)
println(station)
nuke(asteroids, station)


text:
```
MIT License

Copyright (c) 2021 Qi Bao

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
```

include("Graph.jl")
include("Tools.jl")

using SparseArrays
using Laplacians
using Printf
using Statistics

function Exact(G, s) # Returns ci, d, p, idc
	T = time()
	L = getSparseL(G)
	W = getW(G)
	avgs = sum(s)/G.n
	s2 = zeros(G.n)
	for i = 1 : G.n
		s2[i] = s[i] - avgs
	end
	z = W * s
	z2 = W * s2
	# calculate C_I(G)
	lz = L * z
	ci = lz' * lz
	# calculate D(G)
	d = z2' * L * z2
	# calculate P(G)
	p = z2' * z2
	# calculate C(G)
	c = z' * z
	# calculate I_dc(G)
	idc = d + c
	# END CALCULATION
	T = time() - T
	return T, ci, d, p, idc
end

function Approx(G, s; eps = 1e-6) # Returns aci, ad, ap, aidc
	T = time()
	IpL = getSparseIpL(G)
	sL = getSparseL(G)
	f = approxchol_sddm(IpL, tol=0.1*eps)
	avgs = sum(s)/G.n
	s2 = zeros(G.n)
	for i = 1 : G.n
		s2[i] = s[i] - avgs
	end
	z = f(s)
	z2 = f(s2)

	# calculate aC_I(G)
	alz = sL * z
	aci = alz' * alz
	# calculate aD(G)
	ad = z2' * sL * z2
	# calculate aP(G)
	ap = z2' * z2
	# calculate aC(G)
	ac = z' * z
	# calculate aI_dc(G)
	aidc = ad + ac
	# END CALCULATION
	T = time() - T
	return T, aci, ad, ap, aidc
end

# experiment for smaller graphs
function doExp(G, s, lg)
	T, ci, d, p, idc = Exact(G, s)
	T2, aci, ad, ap, aidc = Approx(G, s)
	println(lg, "Exact  Time : ", T)
	println(lg, "Approx Time : ", T2)

	println(lg, "ci : ", ci)
	println(lg, "d : ", d)
	println(lg, "p : ", p)
	println(lg, "idc : ", idc)

	println(lg, "ERROR of ci : ", abs(aci-ci)/ci)
	println(lg, "ERROR of d : ", abs(ad-d)/d)
	println(lg, "ERROR of p  : ", abs(ap-p)/p)
	println(lg, "ERROR of idc : ", abs(aidc-idc)/idc)
	println(lg)
end

# experiment for the large graphs
function doLarge(G, s, lg)
	T2, aci, ad, ap, aidc = Approx(G, s)
	println(lg, "Approx Time : ", T2)

	println(lg, "aci : ", aci)
	println(lg, "ad : ", ad)
	println(lg, "ap : ", ap)
	println(lg, "aidc : ", aidc)

	println(lg)
end

# New Methods
function Approx_marketing(G, s_flag, p, epsilon, t=10)
	z_list, aci_list, ad_list, ap_list, aidc_list = [], [], [], [], [], []
	aci, ad, ap, aidc = 0, 0, 0, 0

	if s_flag == 0
		s = [0 for _ in 1:G.n]
		s = epsilon * p .+ s
		s[s .> 1] .= 1 # set elements of s1 > 1 back to 1

		_, aci, ad, ap, aidc = Approx(G, s)
		push!(z_list, sum(s))
		push!(aci_list, aci)
		push!(ad_list, ad)
		push!(ap_list, ap)
		push!(aidc_list, aidc)
	else
		for _ = 1:t
			if s_flag == 1
				s = Uniform(G.n)
			elseif s_flag == 2
				s = Exponential(G.n)
			elseif s_flag == 3
				s = powerLaw(G.n)
			end

			z_init = sum(s)
			_, aci_init, ad_init, ap_init, aidc_init = Approx(G, s)

			s = epsilon * p .+ s
			s[s .> 1] .= 1 # set elements of s1 > 1 back to 1

			_, aci, ad, ap, aidc = Approx(G, s)
			push!(z_list, (sum(s) - z_init) / z_init)
			push!(aci_list, (aci - aci_init) / aci_init)
			push!(ad_list, (ad - ad_init) / ad_init)
			push!(ap_list, (ap - ap_init) / ap_init)
			push!(aidc_list, (aidc - aidc_init) / aidc_init)
		end
	end
	z_sum = mean(z_list)
	z_std = std(z_list)
	aci = mean(aci_list)
	aci_std = std(aci_list)
	ad = mean(ad_list)
	ad_std = std(ad_list)
	ap = mean(ap_list)
	ap_std = std(ap_list)
	aidc = mean(aidc_list)
	aidc_std = std(aidc_list)

	return z_sum, aci, ad, ap, aidc, z_std, aci_std, ad_std, ap_std, aidc_std
end

function Approx_backfire(G, s_flag, p, epsilon, beta, t=10)
	z_list, aci_list, ad_list, ap_list, aidc_list = [], [], [], [], [], []
	aci, ad, ap, aidc = 0, 0, 0, 0

	for _ = 1:t
		if s_flag == 0
			s = [0.5 for _ in 1:G.n]
		elseif s_flag == 1
			s = Uniform(G.n)
		elseif s_flag == 2
			s = Exponential(G.n)
		elseif s_flag == 3
			s = powerLaw(G.n)
		end

		z_init = sum(s)
		_, aci_init, ad_init, ap_init, aidc_init = Approx(G, s)

		s2 = deepcopy(s)
		for sn in 1:G.n
			if s[sn] >= beta
				s2[sn] += epsilon * p[sn]
			else
				s2[sn] -= epsilon * p[sn]
			end
		end

		s2[s2 .> 1] .= 1 # set elements of s1 > 1 back to 1
		s2[s2 .< 0] .= 0

		_, aci, ad, ap, aidc = Approx(G, s2)
		push!(z_list, (sum(s2) - z_init) / z_init)
		push!(aci_list, (aci - aci_init) / aci_init)
		push!(ad_list, (ad - ad_init) / ad_init)
		push!(ap_list, (ap - ap_init) / ap_init)
		push!(aidc_list, (aidc - aidc_init) / aidc_init)

		if s_flag == 0
			break
		end
	end

	z_sum = mean(z_list)
	z_std = std(z_list)
	aci = mean(aci_list)
	aci_std = std(aci_list)
	ad = mean(ad_list)
	ad_std = std(ad_list)
	ap = mean(ap_list)
	ap_std = std(ap_list)
	aidc = mean(aidc_list)
	aidc_std = std(aidc_list)

	return z_sum, aci, ad, ap, aidc, z_std, aci_std, ad_std, ap_std, aidc_std
end


text:
mutable struct MultigridSolver{T}

    p :: Vector{LinearProblem{T}}  # linear problems
    r :: Vector{AffineMap{T}}  # restrictors
    i :: Vector{AffineMap{T}}   # interpolators
    n :: Int64

end
function MultigridSolver{T}(L::LinearProblem{T}, max_layers::Int64) where {T}
    p = Vector{LinearProblem{T}}(undef, max_layers)
    r = Vector{AffineMap{T}}(undef, max_layers-1)
    i = Vector{AffineMap{T}}(undef, max_layers-1)

    p[1] = L
    return MultigridSolver{T}(p, r, i, 1)
end
MultigridSolver(args...) = MultigridSolver{Float64}(args...)

function add_layer!(M::MultigridSolver{T}, Restr::AffineMap{T}, Interp::AffineMap{T}) where {T}

    Lc = coarsen_problem!(M.p[M.n], Restr, Interp)
    M.p[M.n+1] = Lc
    M.r[M.n] = Restr
    M.i[M.n] = Interp
    M.n += 1
end

function vcycle!(M::MultigridSolver{T}, nj::Vector{Int64}) where {T}

    # initial smoothing
    jsmooth!(M.p[1], nj[1])

    for i=2:M.n
        M.p[i].x[:] = M.r[i-1] * M.p[i-1].x
        compute_residual!(M.p[i])
        jsmooth!(M.p[i], nj[i])
    end

    for i=M.n-1:1
        M.p[i].x[:] = M.i[i] * M.p[i+1].x
        compute_residual!(M.p[i])
        jsmooth!(M.p[i], nj[i])
    end
end


text:
module Model

using Random

# model solution
c_star(w::Real, e::Real, θ::NamedTuple) = θ.γ*(1 - θ.τ)*w + θ.γ*e

l_star(w::Real, e::Real, θ::NamedTuple) = 1 - θ.γ + (1 - θ.γ)/(1 - θ.τ)*(e/w)

# objective function
function obj_fun(θ::NamedTuple,
                 w::AbstractVector,
                 mom_data::Tuple{Float64,Float64,Float64},
                 moments_fun::Function,
                 rng::AbstractRNG=MersenneTwister(893245) # Allow possing in an RNG to use when drawing random numbers
                 )

    S = 100 # number of simulation draws for each observed individual
    n = length(w)    # number of individuals
    w = repeat(w, S) # stack observed wages S times
    e = θ.σ*randn(rng, n*S)

    con = c_star.(w, e, Ref(θ))
    lab = l_star.(w, e, Ref(θ))

    # Calculate moments based on simulated data for the value of θ
    mom_sim = moments_fun(w, con, lab)

    # Calculate objective function as squared difference
    Q = sum(abs2(data - sim) for (data, sim) in zip(mom_data, mom_sim))

    return Q
end

end


text:
#InspectDR: Top-level interface for Cairo layer
#-------------------------------------------------------------------------------

#Render base annotation (titles, axis labels, ...)
#-------------------------------------------------------------------------------
function render_baseannotation(canvas::PCanvas2D, graphinfo::Graph2DInfo,
		xticklabels::Bool, plot::Plot2D, istrip::Int)
	strip = plot.strips[istrip] #WANTCONST
	if plot.layout.values.enable_legend
		legend_render(canvas, plot, istrip)
	end

	#Draw data area & grid lines
Cairo.save(canvas.ctx)
	drawrectangle(canvas.ctx, canvas.graphbb, plot.layout.values.frame_data)
Cairo.restore(canvas.ctx)
	render_grid(canvas, plot.layout.values, graphinfo.grid)

	render_axes(canvas, plot.layout.values, graphinfo.grid, plot.xscale, strip.yscale, xticklabels)
	return
end
function render_baseannotation(ctx::CairoContext, bb::BoundingBox, plotinfo::Plot2DInfo,
		databb::BoundingBox, graphbblist::Vector{BoundingBox}, plot::Plot2D)

Cairo.save(ctx)
	drawrectangle(ctx, bb, plot.layout.values.frame_canvas)
Cairo.restore(ctx)

	#Render titles, axis labels, ...
	#TODO: use Plot2DInfo instead of databb/graphbblist??:
	render(ctx, plot.annotation, bb, databb, graphbblist, plot.layout.values)

	nstrips = length(plot.strips)
	for i in 1:nstrips
		graphinfo = Graph2DInfo(plotinfo, i)
		canvas = PCanvas2D(ctx, bb, graphinfo)
		render_xticklabels = (i == nstrips) #Only bottom-most graph
		render_baseannotation(canvas, graphinfo, render_xticklabels, plot, i)
	end
end

#Render actual plot data
#-------------------------------------------------------------------------------
function render_data(canvas::PCanvas2D, plot::Plot2D, istrip::Int)
Cairo.save(canvas.ctx)
	setclip(canvas.ctx, canvas.graphbb)

	render(canvas, plot.display_data, istrip)
	if plot.displayNaN
		rendernans(canvas, plot.data, istrip)
	end
Cairo.restore(canvas.ctx)
end
function render_data(ctx::CairoContext, bb::BoundingBox, plotinfo::Plot2DInfo, plot::Plot2D)
	nstrips = length(plot.strips)
	for i in 1:nstrips
		graphinfo = Graph2DInfo(plotinfo, i)
		canvas = PCanvas2D(ctx, bb, graphinfo)
		render_data(canvas, plot, i)
	end
end

#Render secondary plot annotation & redraw frame:
#-------------------------------------------------------------------------------
function render_userannotation(canvas::PCanvas2D, graphinfo::Graph2DInfo, plot::Plot2D, istrip::Int)
Cairo.save(canvas.ctx)
	setclip(canvas.ctx, canvas.graphbb)
	render(canvas, plot.userannot, graphinfo, istrip)
	render(canvas, plot.parentannot, graphinfo, istrip)
Cairo.restore(canvas.ctx)

	render_graphframe(canvas, plot.layout.values.frame_data) #Redraw frame
	#TODO: should we be more careful about setclip() - instead of redrawing frame?
end
function render_userannotation(ctx::CairoContext, bb::BoundingBox, plotinfo::Plot2DInfo, plot::Plot2D)
	nstrips = length(plot.strips)
	for i in 1:nstrips
		graphinfo = Graph2DInfo(plotinfo, i)
		canvas = PCanvas2D(ctx, bb, graphinfo)
		render_userannotation(canvas, graphinfo, plot, i)
	end
end

#Render entire plot within provided bounding box:
#-------------------------------------------------------------------------------
#TODO: Should base API be written such that user provides databb, and have
#      InspectDR calculate outwards instead?
function render(ctx::CairoContext, plot::Plot2D, bb::BoundingBox)
	lyt = plot.layout.values #WANTCONST
	databb = databounds(bb, lyt, grid1(plot)) #WANTCONST
	nstrips = length(plot.strips) #WANTCONST

	graphbblist = graphbounds_list(databb, lyt, nstrips)
	refreshed = update_ddata(plot) #Also computes new extents
	plotinfo = Plot2DInfo(plot, bb) #Compute new extents before calling.

	_reset(ctx) #Ensure known state.
	render_baseannotation(ctx, bb, plotinfo, databb, graphbblist, plot)
	render_data(ctx, bb, plotinfo, plot)
	render_userannotation(ctx, bb, plotinfo, plot)

	return plotinfo #So parent/caller can position objects/format data
end

#Render entire plot plot using CairoBufferedPlot:
#-------------------------------------------------------------------------------
#TODO: split in parts so that GUI can be refreshed before refreshing data
#TODO: Or maybe send in final context so it can be done here.
function render(bplot::CairoBufferedPlot, plot::Plot2D, bb::BoundingBox, refreshdata::Bool)
	lyt = plot.layout.values #WANTCONST
	databb = databounds(bb, lyt, grid1(plot)) #WANTCONST
	nstrips = length(plot.strips) #WANTCONST
	refreshed = false

	graphbblist = graphbounds_list(databb, lyt, nstrips)
	if refreshdata
		refreshed = update_ddata(plot) #Also computes new extents
		#TODO: don't refresh data if !refresh && bounding box did not change??
	end
	plotinfo = Plot2DInfo(plot, bb) #Compute new extents before calling.

	ctx = CairoContext(bplot.surf)
	_reset(ctx); clear(ctx)
	render_baseannotation(ctx, bb, plotinfo, databb, graphbblist, plot)

	if refreshdata #Update data surface:
		ctxdata = CairoContext(bplot.data)
		_reset(ctxdata); clear(ctxdata)
		render_data(ctxdata, bb, plotinfo, plot)
		Cairo.destroy(ctxdata)
	end

	#Now render data portion onto main context:
	Cairo.set_source_surface(ctx, bplot.data, 0, 0)
	Cairo.paint(ctx) #Applies contents of bplot.data

	render_userannotation(ctx, bb, plotinfo, plot)
	Cairo.destroy(ctx)

	return plotinfo #So parent/caller can position objects/format data
end

#Last line


text:
abstract AbstractCell
abstract AbstractCells{T}

function styles{T <: AbstractStyleElement}(::Type{T}, tbl::Table)
  n = length(tbl.styles)
  styleels = Array(AbstractStyleElement, n)
  visited = fill(false, n)

  for rng in tbl.layout
    datael_id = rng.datael_id
    style_id  = rng.style_id
    if !visited[style_id]
      styleels[style_id] = style(T, tbl.data[datael_id], tbl.styles[style_id])
      visited[style_id]  = true
    end
  end
  return styleels
end

tcell(::Type{AbstractCell},
      styleel::AbstractStyleElement,
      datael::AbstractDataElement,
      i, j) = error("cell type must be concrete")

function writerange!{T <: AbstractCell}(cells::Matrix{T}, tbl::Table, rng::CellRange)
    box    = rng.box
    styleel  = tbl.styles[rng.style_id]
    datael = tbl.data[rng.datael_id]

    offset_x = xmin(box)
    offset_y = ymin(box)
    n, m = size(datael)

    for i=1:n
        for j=1:m
            cells[offset_y + i,offset_x + j] = tcell(T, styleel, datael, i, j)
        end
    end
    return nothing
end

function convert{T <: AbstractCell}(::Type{Matrix{T}}, tbl::Table)
    cells = Array(T, size(tbl))
    for rng in tbl.layout
        writerange!(cells, tbl, rng)
    end
    return cells
end

renderformatspecific{T <: AbstractCell}(
  io::IO, ::Type{T}, cells::AbstractCells) = error("eltype(cells) must be concrete")

function render{T <: AbstractCell, S <: AbstractStyleElement}(
    io::IO, ::Type{T}, ::Type{S}, tbl::Table)

  styleels = styles(S, tbl)
  newtbl = Table(tbl.data, tbl.layout, styleels)

  cells = convert(Matrix{T}, newtbl)
  renderformatspecific(io, T, cells)
end


text:
module PolarFact

using LinearAlgebra
using Printf

import LinearAlgebra: transpose!

export polarfact

include("common.jl")     # common functions and types
include("newton.jl")     # using Newton's method
include("svd.jl")        # using the SVD method
include("halley.jl")     # using Halley's method and the QDWH method
include("hybrid.jl")     # using hybrid method

function polarfact(A::Matrix{T};
                   alg::Symbol=:newton,
                   maxiter::Integer=:100,
                   tol::Real = cbrt(eps(T)),
                   verbose::Bool=false) where {T}

    # choose algorithm
    algorithm =
       alg == :newton ? NewtonAlg{T}(; maxiter=maxiter, tol=tol, verbose=verbose) :
       alg == :qdwh ?  QDWHAlg{T}(; maxiter=maxiter, tol=tol, verbose=verbose) :
       alg == :halley ? HalleyAlg{T}(; maxiter=maxiter, tol=tol, verbose=verbose) :
       alg == :svd ? SVDAlg() :
       alg == :schulz ? NewtonSchulzAlg{T}(; maxiter=maxiter, tol=tol, verbose=verbose) :
       alg == :hybrid ? NewtonHybridAlg{T}(; maxiter=maxiter, tol=tol, verbose=verbose) :
       error("Invalid algorithm.")

    # Initialization: if m > n, do QR factorization
    m, n = size(A)
    if m > n
        if alg in [:newton, :qdwh, :halley, :schulz, :hybrid]
            m = n
            A = qrfact(A)[:R]
        end
    elseif m < n
        error("The row dimension of the input matrix must be
              greater or equal to column dimension.")
    end
    U = Array{T}(undef, m, n)
    H = Array{T}(undef, n, n)
    # solve for polar factors
    solve!(algorithm, A, U, H)
end

end # module


text:
# Sparse Bitmapped Tries
# ======================

abstract SparseBitmappedTrie{T} <: BitmappedTrie{T}

immutable SparseNode{T} <: SparseBitmappedTrie{T}
    arr::Vector{SparseBitmappedTrie{T}}
    shift::Int
    length::Int
    maxlength::Int
    bitmap::Int
end
SparseNode(T::Type) = SparseNode{T}(SparseBitmappedTrie{T}[], shiftby*7, 0, trielen^7, 0)

immutable SparseLeaf{T} <: SparseBitmappedTrie{T}
    arr::Vector{T}
    bitmap::Int

    SparseLeaf(arr::Vector, bitmap::Int) = new(arr, bitmap)
    SparseLeaf() = new(T[], 0)
end

arrayof(    n::SparseNode) = n.arr
shift(      n::SparseNode) = n.shift
maxlength(  n::SparseNode) = n.maxlength
Base.length(n::SparseNode) = n.length

arrayof(    l::SparseLeaf) = l.arr
shift(       ::SparseLeaf) = 0
maxlength(  l::SparseLeaf) = trielen
Base.length(l::SparseLeaf) = length(arrayof(l))

function demoted{T}(n::SparseNode{T})
    shift(n) == shiftby ?
    SparseLeaf{T}(T[], 0) :
    SparseNode{T}(SparseBitmappedTrie{T}[],
                  shift(n) - shiftby,
                  0,
                  round(Int, maxlength(n) / trielen), 0)
end

bitpos(  t::SparseBitmappedTrie, i::Int) = 1 << (mask(t, i) - 1)
hasindex(t::SparseBitmappedTrie, i::Int) = t.bitmap & bitpos(t, i) != 0
index(   t::SparseBitmappedTrie, i::Int) =
    1 + count_ones(t.bitmap & (bitpos(t, i) - 1))

function update{T}(l::SparseLeaf{T}, i::Int, el::T)
    hasi = hasindex(l, i)
    bitmap = bitpos(l, i) | l.bitmap
    idx = index(l, i)
    if hasi
        newarr = arrayof(l)[:]
        newarr[idx] = el
    else
        newarr = vcat(arrayof(l)[1:idx-1], [el], arrayof(l)[idx:end])
    end
    (SparseLeaf{T}(newarr, bitmap), !hasi)
end
function update{T}(n::SparseNode{T}, i::Int, el::T)
    bitmap = bitpos(n, i) | n.bitmap
    idx = index(n, i)
    if hasindex(n, i)
        newarr = arrayof(n)[:]
        updated, inc = update(newarr[idx], i, el)
        newarr[idx] = updated
    else
        child, inc = update(demoted(n), i, el)
        newarr = vcat(arrayof(n)[1:idx-1], [child], arrayof(n)[idx:end])
    end
    (SparseNode{T}(newarr,
                   n.shift,
                   inc ? n.length + 1 : n.length,
                   n.maxlength, bitmap),
     inc)
end

Base.get(n::SparseLeaf, i::Int, default) =
    hasindex(n, i) ? arrayof(n)[index(n, i)] : default
Base.get(n::SparseNode, i::Int, default) =
    hasindex(n, i) ? get(arrayof(n)[index(n, i)], i, default) : default

function Base.start(t::SparseBitmappedTrie)
    t.length == 0 && return true
    ones(Int, 1 + round(Int, t.shift / shiftby))
end

function directindex(t::SparseBitmappedTrie, v::Vector{Int})
    isempty(v) && return arrayof(t)
    local node = arrayof(t)
    for i=v
        node = node[i]
        node = isa(node, SparseBitmappedTrie) ? arrayof(node) : node
    end
    node
end

Base.done(t::SparseBitmappedTrie, state) = state == true

function Base.next(t::SparseBitmappedTrie, state::Vector{Int})
    item = directindex(t, state)
    while true
        index = pop!(state)
        node = directindex(t, state)
        if length(node) > index
            push!(state, index + 1)
            return item, vcat(state, ones(Int, 1 + round(Int, t.shift / shiftby) -
                                               length(state)))
        elseif is(node, arrayof(t))
            return item, true
        end
    end
end

# Persistent Hash Maps
# ====================

immutable PersistentHashMap{K, V} <: PersistentMap{K, V}
    trie::SparseBitmappedTrie{PersistentArrayMap{K, V}}
    length::Int

    PersistentHashMap(trie, length) = new(trie, length)
    PersistentHashMap() = new(SparseNode(PersistentArrayMap{K, V}), 0)
end

function PersistentHashMap(itr)
    if length(itr) == 0
        return PersistentHashMap()
    end
    if VERSION >= v"0.4.0-dev"
        K, V = typejoin(map(typeof, itr)...).types
    else
        K, V = typejoin(map(typeof, itr)...)
    end
    m = PersistentHashMap{K, V}()
    for (k, v) in itr
        m = assoc(m, k, v)
    end
    m
end

function PersistentHashMap(kvs::(@compat Tuple{Any, Any})...)
    PersistentHashMap([kvs...])
end

function PersistentHashMap(; kwargs...)
    isempty(kwargs) ?
    PersistentHashMap{Any, Any}() :
    PersistentHashMap(kwargs...)
end

Base.length(m::PersistentHashMap) = m.length
Base.isempty(m::PersistentHashMap) = length(m) == 0

# avoid anon function JuliaLang/julia#1864
tup_isequal(x) = isequal(x[1], x[2])
Base.isequal(m1::PersistentHashMap, m2::PersistentHashMap) =
    length(m1) == length(m2) && all(tup_isequal, zip(m1, m2))

tup_eq(x) = x[1] == x[2]
==(m1::PersistentHashMap, m2::PersistentHashMap) =
    length(m1) == length(m2) && all(tup_eq, zip(m1, m2))

function _update{K, V}(f::Function, m::PersistentHashMap{K, V}, key)
    keyhash = reinterpret(Int, hash(key))
    arraymap = get(m.trie, keyhash, PersistentArrayMap{K, V}())
    newmap = f(arraymap)
    newtrie, _ = update(m.trie, keyhash, newmap)
    PersistentHashMap{K, V}(newtrie,
                            m.length + (length(newmap) < length(arraymap) ? -1 :
                                        length(newmap) > length(arraymap) ? 1 :
                                        0))
end

function assoc{K, V}(m::PersistentHashMap{K, V}, key, value)
    _update(m, key) do arraymap
        assoc(arraymap, key, value)
    end
end

function dissoc(m::PersistentHashMap, key)
    _update(m, key) do arraymap
        dissoc(arraymap, key)
    end
end

function Base.getindex(m::PersistentHashMap, key)
    val = get(m.trie, reinterpret(Int, hash(key)), NotFound())
    is(val, NotFound()) && error("key not found")
    val[key]
end

Base.get(m::PersistentHashMap, key) = m[key]
function Base.get(m::PersistentHashMap, key, default)
    val = get(m.trie, reinterpret(Int, hash(key)), NotFound())
    is(val, NotFound()) && return default
    val[key]
end

function Base.haskey(m::PersistentHashMap, key)
    get(m.trie, reinterpret(Int, hash(key)), NotFound()) != NotFound()
end

function Base.start(m::PersistentHashMap)
    state = start(m.trie)
    done(m.trie, state) && return (Any[], state)
    arrmap, triestate = next(m.trie, state)
    (arrmap.kvs, triestate)
end
Base.done(m::PersistentHashMap, state) =
    isempty(state[1]) && done(m.trie, state[2])

function Base.next(m::PersistentHashMap, state)
    kvs, triestate = state
    if isempty(kvs)
        arrmap, triestate = next(m.trie, triestate)
        next(m, (arrmap.kvs, triestate))
    else
        (convert(Tuple, kvs[1]), (kvs[2:end], triestate))
    end
end

function Base.map(f::(@compat Union{Function, DataType}), m::PersistentHashMap)
    PersistentHashMap([f(kv) for kv in m]...)
end

function Base.filter{K, V}(f::Function, m::PersistentHashMap{K, V})
    arr = Array((@compat Tuple{K, V}), 0)
    for el in m
        f(el) && push!(arr, el)
    end
    isempty(arr) ? PersistentHashMap{K, V}() : PersistentHashMap(arr...)
end

# Suppress ambiguity warning while allowing merging with array
function _merge(d::PersistentHashMap, others...)
    acc = d
    for other in others
        for (k, v) in other
            acc = assoc(acc, k, v)
        end
    end
    acc
end

# This definition suppresses ambiguity warning
Base.merge(d::PersistentHashMap, others::Associative...) =
    _merge(d, others...)
Base.merge(d::PersistentHashMap, others...) =
    _merge(d, others...)

@compat function Base.show{K, V}(io::IO, ::MIME"text/plain", m::PersistentHashMap{K, V})
    print(io, "Persistent{$K, $V}[")
    print(io, join(["$k => $v" for (k, v) in m], ", "))
    print(io, "]")
end


text:
using Statistics

export bachelierFormula

"Price of a vanilla option under the Bachelier/Normal model"
function bachelierFormula(isCall::Bool, strike::T, forward::T, variance::T, discountDf::T) where {T}
    sign = 1.0
    if !isCall
        sign = -1.0
    end
    if variance == 0
        return Max(sign * (forward - strike), 0) * discountDf
    end
    sqrtvar = sqrt(variance)
    d = sign * (forward - strike) / sqrtvar
    if forward == strike
        return sqrtvar / sqrt(2 * pi)
    end
    normal = Normal()
    Nd = cdf(normal, d)
    nd = pdf(normal, d)
    return discountDf * (sign * (forward - strike) * Nd + sqrtvar * nd)
end


text:
import HTTP2
using HTTP2.Frame
using Test
using Sockets
using MbedTLS

function sslconnect(dest, port, certhostname)
    println("Connecting over SSL ...")
    buffer = connect(dest, port)
    sslconfig = MbedTLS.SSLConfig(false)
    sslbuffer = MbedTLS.SSLContext()
    MbedTLS.setup!(sslbuffer, sslconfig)
    MbedTLS.set_bio!(sslbuffer, buffer)
    MbedTLS.hostname!(sslbuffer, certhostname)
    MbedTLS.handshake!(sslbuffer)
    return sslbuffer
end

function show_response(headers, body)
    for header in headers
        @info("Result header: " * String(header[1]) * ": " * String(header[2]))
    end
    @info("Result body: " * String(body))
end

function test_request(dest, port, url, certhostname=nothing)
    for conn_id in 1:3
        @info("Opening connection", conn_id)
        buffer = (certhostname === nothing) ? connect(dest, port) : sslconnect(dest, port, certhostname)
        connection = HTTP2.Session.new_connection(buffer; isclient=true)
        headers = [(":method", "GET"),
                   (":path", url),
                   (":scheme", "http"),
                   (":authority", "127.0.0.1:9000"),
                   ("accept", "*/*"),
                   ("accept-encoding", "gzip, deflate"),
                   ("user-agent", "HTTP2.jl")]

        for req_id in 1:3
            @info("Sending request", req_id)
            HTTP2.Session.put_act!(connection, HTTP2.Session.ActSendHeaders(HTTP2.Session.next_free_stream_identifier(connection), headers, true))

            resp_headers = HTTP2.Session.take_evt!(connection).headers
            @show resp_headers
            resp_data = HTTP2.Session.take_evt!(connection).data
            @show resp_data
            show_response(resp_headers, resp_data)
        end

        @info("Closing connection", conn_id)
        close(connection)
        sleep(2) # wait for close message to percolate
        close(buffer)
    end
end

if length(ARGS) == 1
    test_request(ip"127.0.0.1", 8000, "/", ARGS[1])
else
    test_request(ip"127.0.0.1", 8000, "/")
end


text:
# findall cherche toutes les rexexp dans une chaine
# (devrait apparaitre avec cette API dans julia-1.3
#
# Exemple d'utilisation
#   findall("ing", "Spinning laughing dancing")
#   =>
#   3-element Array{UnitRange{Int64},1}:
#    6:8
#    15:17
#    23:25
#
#   txt "Spinning laughing dancing"
#   rpat = r"\w+"
#   ranges = [txt[r] for r in findall(rpat, txt)]
#   map(r->s[r], ranges)
#   =>
#   3-element Array{String,1}:
#    "Spinning"
#    "laughing"
#    "dancing"
#
#   txt "Spinning laughing dancing"
#   rpat = r"\w+"
#   words = [txt[r] for r in findall(rpat, txt)]
#   =>
#   3-element Array{String,1}:
#    "Spinning"
#    "laughing"
#    "dancing"
#
if !hasmethod(Base.findall, Tuple{Regex, String})
    # Julia pas assez récent (< julia-1.3)
    @info "On définit : Base.findall pour Tuple{Regex, String}"
    function Base.findall(t::Union{AbstractString,Regex}, s::AbstractString;
                     overlap::Bool=false)
        found = UnitRange{Int}[]
        i, e = firstindex(s), lastindex(s)
        while true
            r = findnext(t, s, i)
            isnothing(r) && return found
            push!(found, r)
            j = overlap || isempty(r) ? first(r) : last(r)
            j > e && return found
            @inbounds i = nextind(s, j)
        end
    end
else
    @info "Déjà définie : Base.findall pour Tuple{Regex, String}"
end

#./


text:
const PROMOTE_OPTIONS = ["Rook", "Knight", "Bishop", "Queen"]
const COLORS = Dict(1 => "White", -1 => "Black")

function log(msg)
  buf = IOBuffer()
  ctx = IOContext(buf, :color => true)
  printstyled(ctx, "pacoplay ", color = 245)
  println(String(take!(buf)) * msg)
end

function log(match :: Int, msg)
  buf = IOBuffer()
  ctx = IOContext(buf, :color => true)
  printstyled(ctx, "pacoplay:", color = 245)
  printstyled(ctx, "$match", color = 245)
  printstyled(ctx, " ", color = 245)
  println(String(take!(buf)) * msg)
end

function logerr(msg)
  buf = IOBuffer()
  ctx = IOContext(buf, :color => true)
  printstyled(ctx, "pacoplay ", color = 245)
  printstyled(ctx, "Error: $msg\n", color = red)
  println(stderr, String(take!(buf)) * msg)
end

gamehash(game) = "#" * string(Game.hash(game))[1:5]

# transform an absolute field index on the chess board (always from white's
# perspective) to the perspective of the current players
function adapt_perspective(index, current_player :: Int)
  if current_player == 1
    index
  else
    x = index % 8
    y = div(index, 8)
    x + 8 * (7 - y)
  end
end

# Convert pacoplay action json to action index
# Examples would be {"Lift" : 5}, {"Place" : 7} or {"Promote" : "Bishop"}
function parse_action(json, current_player) :: Int
  if "Lift" in keys(json)
    1 + adapt_perspective(json["Lift"], current_player)
  elseif "Place" in keys(json)
    1 + 64 + adapt_perspective(json["Place"], current_player)
  elseif "Promote" in keys(json)
    piece_type = json["Promote"]
    offset = findfirst(isequal(piece_type), PROMOTE_OPTIONS)
    if !isnothing(offset)
      2 * 64 + offset
    else
      msg = "Piece type '$piece_type' in JSON is not understood (promotion)."
      throw(ArgumentError(msg))
    end
  else
    throw(ArgumentError("Action in JSON is not understood."))
  end
end

function parse_game(state_json)
  if !("actions" in keys(state_json))
    throw(ArgumentError("Provided JSON has no field 'actions'"))
  end
  game = PacoSako()
  for pacoplay_action in state_json["actions"]
    action = parse_action(pacoplay_action, Game.current_player(game))
    Game.apply_action!(game, action)
  end
  victory_string = nothing
  if Game.is_over(game)
    if Game.status(game) in [-1, 1]
      victory_string = "Winner: $(COLORS[Game.status(game)])"
    else
      victory_string = "Draw"
    end
  elseif "TimeoutVictory" in keys(state_json["victory_state"])
    winner = state_json["victory_state"]["TimeoutVictory"]
    victory_string = "Winner: $winner (by timeout)"
  end
  game, victory_string
end

# Convert action index to pacoplay action json string
function action_string(id :: Int, current_player :: Int) :: String
  if 0 < id <= 64
    body = "{\"Lift\":$(adapt_perspective(id - 1, current_player))}"
  elseif 64 < id <= 128
    body = "{\"Place\":$(adapt_perspective(id - 1 - 64, current_player))}"
  elseif 128 < id <=132
    body = "{\"Promote\":\"$(PROMOTE_OPTIONS[id - 128])\"}"
  else
    throw(ArgumentError("Action index $id is not supported"))
  end
end

function do_action_string(match, id, current_player :: Int) :: String
  action = action_string(id, current_player)
  "{\"DoAction\":{\"key\":\"$match\",\"action\":$action}}"
end

function subscription_string(match :: Int) :: String
  "{\"SubscribeToMatch\":{\"key\":\"$match\"}}"
end

function websocket_url(server_url :: String)
  parts = splitpath(server_url)
  if parts[1] in ["ws:", "http:"]
    "ws://" * joinpath(parts[2:end]..., "websocket")
  elseif parts[1] in ["wss:", "https:"]
    "wss://" * joinpath(parts[2:end]..., "websocket")
  else
    "wss://" * joinpath(parts..., "websocket")
  end
end

function subscribe(ws, match :: Int)
  write(ws, subscription_string(match))
  response = readavailable(ws)
  json = nothing
  try
    json = LazyJSON.parse(String(response))
  catch _
    logerr("Could not parse server response")
    return nothing
  end
  if "MatchConnectionSuccess" in keys(json)
    return json["MatchConnectionSuccess"]
  else
    logerr("Could not understand server response")
    return nothing
  end
end

function check_updates(games :: Channel)
  new_victory = nothing
  new_state = nothing
  while isready(games)
    new_state, new_victory = take!(games)
  end
  new_state, new_victory
end

function wait_updates(games :: Channel)
  wait(games)
  check_updates(games)
end

function play_match( server_url :: String
                   , match :: Int
                   , color :: Int = 1
                   ; delay :: Float64 = 0.1
                   , player = Player.MCTSPlayer(power = 1000))

  @assert color in [-1, 1] "AI must either play as White (1) or Black (-1)"

  ws_url = websocket_url(server_url)

  log("Connecting...")
  HTTP.WebSockets.open(ws_url) do ws

    log("Subscribing...")
    # Try to connect to the websocket and hope that it responds as anticipated
    json = subscribe(ws, match)
    if isnothing(json) return end

    # We connected successfully
    log("Assuming color $(COLORS[color]) in match $match on $server_url")

    game, victory = parse_game(json["state"])
    games = Channel(100)
    log(match, "Received first state $(gamehash(game))")

    # Now we let the player play if it has the right color
    @sync begin

      exiting = false

      # listen to pacoplay changes of the game state and fill the games channel
      @async while true
        try
          msg = readavailable(ws)
          json = LazyJSON.parse(String(msg))
          state = parse_game(json["CurrentMatchState"])
          put!(games, state)
          log(match, "Received state $(gamehash(state[1]))")
        catch err
          if err isa EOFError && !exiting
            logerr("Websocket has closed")
          elseif err isa LazyJSON.ParseError
            logerr("Could not parse server message")
          elseif err isa KeyError
            logerr("Could not understand server message")
          elseif err isa ArgumentError
            logerr(err.msg)
          elseif err isa InvalidStateException && !exiting
            # the games channel has been closed
          end
          exiting = true
          close(ws)
          close(games)
          break
        end
      end

      # think about next action(s) to take and submit them
      @async while true

        if !isnothing(victory)
          log(match, "Match is over. " * victory)
          exiting = true
          close(ws)
          close(games)
          break
        end

        # the following could raise acceptions when waiting on games
        # or when trying to write to ws after it was closed
        try
          if color != Game.current_player(game)
            log(match, "Waiting for opponent...")
            game, victory = wait_updates(games)
            continue
          end

          log(match, "Thinking...")
          actions = Player.decide_chain(player, game)

          # since finding the decision might have taken some time, we
          # check if the game state has changed by human intervention
          # before submitting
          sleep(delay)
          current_game, vs = check_updates(games)

          # the game has changed. don't submit
          if !isnothing(vs) || (!isnothing(current_game) && (current_game != game))
            log(match, "Info: Game state was changed while player was thinking")
            game = current_game
            victory = vs
            continue
          end

          # the game has not changed, and we try to submit all actions
          log(match, "Submitting actions $(join(actions, ", "))")
          for action in actions
            # do action
            write(ws, do_action_string(match, action, Game.current_player(game)))
            Game.apply_action!(game, action)
            # wait for feedback from the server
            sleep(delay)
            current_game, vs = wait_updates(games)
            if current_game != game
              # manual intervention changed the game, not our submitted action
              log(match, "Info: Game state was changed while player was acting")
              break
            end
            game = current_game
            victory = vs
          end
        catch err
          if err isa InvalidStateException && !exiting
            logerr("Games channel was closed unexpectedly")
          elseif !exiting
            logerr(err)
          end
          exiting = true
          close(ws)
          close(games)
          break
        end
      end
    end
  end
  log("Exiting...")
end

max_stars_repo_name: kreibaum/pacosako

max_issues_repo_name: kreibaum/pacosako

max_forks_repo_name: kreibaum/pacosako


text:
#! /usr/bin/julia

# Rosetta Code, Execute a system command

run(`ls`)

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
#---------------------------#
# RESTRICTED 4-BODY SYSTEMS #
#---------------------------#

default_reference_frame(::Abstract_R4BPModel) = SynodicFrame()

Base.show(io::IO, x::Abstract_R4BPModel) = print(io, "$(nameof(typeof(x)))$(x.props)")

#-------------------#
# SYSTEM PROPERTIES #
#-------------------#
struct R4BPSystemProperties{N<:Number}
    b1 :: Symbol  # Identifier of body 1 (central body)
    b2 :: Symbol  # Identifier of body 2 (smaller body)
    b3 :: Symbol  # Identifier of body 3 (large distance body)
    m3 :: N       # Normalised mass of body 3 (relative to mass of R3BP)
    a3 :: N       # Normalised distance between body 3 and R3BP barycenter
    ω3 :: N       # Angular velocity of body 3 in synodic coordinates
    μ  :: N       # Mass ratio of the r3bp system
    μ2 :: N       # Mass ratio of the whole system
    α0 :: N       # Initial phase angle of body 3
end

# @memoize function R4BPSystemProperties(a::Symbol, b::Symbol, c::Symbol)
function R4BPSystemProperties(a::Symbol, b::Symbol, c::Symbol, α0::Number=0.; kwargs...)
    # NOTE: c is the 3rd body (e.g. the Sun), a and b are the primaries
    # Get the 3-body problem properties
    r3bp = R3BPSystemProperties(a, b)

    # Run-time compute remaining values
    GM = @. SpiceUtils.get_GM((a, b, c))u"km^3/s^2"
    m3 = GM[3] / (GM[1] + GM[2])
    state_3 = Array(SpiceUtils.get_state(0., a, c))
    elements_3 = oscltx(state_3, 0., ustrip(u"km^3/s^2", GM[3]))
    a3 = get(kwargs, :a3, elements_3[10]u"km" / r3bp.L)
    ω3 = get(kwargs, :ω3, √((1 + m3) / a3^3) - 1)  # [DeiTos2017 eq.18]
    μ2 = get(kwargs, :μ2, sum(GM) / (GM[1] + GM[2]))
    μ = get(kwargs, :μ, r3bp.μ)

    R4BPSystemProperties(a, b, c, m3, a3, ω3, μ, μ2, α0)
end

Base.show(io::IO, x::R4BPSystemProperties) = print(io, parameters(x))#(a=x.b1, b=x.b2, c=x.b3, α0=x.α0))
ModelingToolkit.parameters(props::R4BPSystemProperties) = NamedTuple([i => getfield(props, i) for i in propertynames(props)])
ModelingToolkit.parameters(model::Abstract_R4BPModel) = [getfield(model.props, i.name) for i in ModelingToolkit.parameters(model.ode.ode_system)]

#---------#
# METHODS #
#---------#

primary_body(m::Abstract_R4BPModel) = m.props.b1
secondary_body(m::Abstract_R4BPModel) = m.props.b2

#----------#
# INCLUDES #
#----------#
include("bc4bp.jl")


text:
using Base.Test

function test_output(L,iters,m)

    const mint = 2
    const maxt = 20
    run(`./exe/annealer -l $L -t $mint -T $maxt -i $iters -m $m 42 debugfile`);
    io = open("debugfile")

    # Firstly, there should be some fields with information.
    header = map(Char, read(io,Int8,8)) |> join
    @test header == "ANNEALER"

    is_big_endian = read(io,Int64)
    @test is_big_endian ∈ [0,1]

    rL = read(io,Int64)
    @test rL == L

    ntws = read(io,Int64)
    tws = read(io,Int64,ntws)
    @test minimum(tws) == 1
    @test maximum(tws) == iters
    @test all(1≤tw≤iters for tw in tws)

    nts = read(io,Int64)
    ts = read(io,Int64,nts)
    @test minimum(ts) == mint
    @test maximum(ts) == maxt
    @test all(mint≤t≤maxt for t in ts)

    total_meas = read(io, Int64)
    @test total_meas ≤ ntws*nts

    T = read(io, Float64)
    @test T ≈ 42

    Jup = read(io,UInt64,L*L*L)
    Jright = read(io,UInt64,L*L*L)
    Jfront = read(io,UInt64,L*L*L)

    for i in 1:L*L*L # Check J∈0,1
        @test all(J∈('0','1') for J in oct(Jup[i]))
        @test all(J∈('0','1') for J in oct(Jright[i]))
        @test all(J∈('0','1') for J in oct(Jfront[i]))
    end

    mcs = Int64[]
    for m in 1:total_meas
        mc = read(io,Int64)
        push!(mcs,mc)
        spins = read(io,UInt64,L*L*L)
        for i in 1:L*L*L # Check σ∈0,1
            @test all(σ∈('0','1') for σ in oct(spins[i]))
        end
    end

    @test minimum(mcs) == 1
    @test maximum(mcs) == maximum(tws)+maximum(ts)
    @test iters ∈ mcs

    footer = map(Char, read(io,Int8,3)) |> join
    @test footer == "END"

    @test eof(io)

    close(io)
    rm("debugfile")
end

@testset "Output format" begin
    const iters = 10
    for m in 2:iters
        for L in [2,8]
            test_output(L, iters, m)
        end
    end
end

@testset "tw and t specification" begin

    function getfields(cmd)
        run(cmd);
        io = open("debugfile")
        header = map(Char, read(io,Int8,8)) |> join
        is_big_endian = read(io,Int64)
        L = read(io,Int64)
        ntws = read(io,Int64)
        tws = read(io,Int64,ntws)
        nts = read(io,Int64)
        ts = read(io,Int64,nts)
        total_meas = read(io, Int64)
        T = read(io, Float64)
        Jup = read(io,UInt64,L*L*L)
        Jright = read(io,UInt64,L*L*L)
        Jfront = read(io,UInt64,L*L*L)
        mcs = Int64[]
        for m in 1:total_meas
            mc = read(io,Int64)
            push!(mcs,mc)
            spins = read(io,UInt64,L*L*L)
        end
        return mcs,tws,ts
    end

    function glasslogspace(from,to,N)
        v = logspace(log10(from),log10(to),N)
        v = floor(Int64,v)
        v[1] = from
        v[end] = to
        return unique(v)
    end

    glasslinspace(from,to,N) = floor(Int64,linspace(from,to,N)) |> unique

    sortunique(x) = x |> unique |> sort

    for m in [10,50,150], M in m+[10,50,150]
        # Normal scale for t and tws
        mcs,tws,ts = getfields(`./exe/annealer -m 84 -i 423 -t $m -T $M -n 15 42 debugfile`)
        @test tws == glasslogspace(1,423,84)
        @test ts  == glasslinspace(m,M,15)
        tw_plus_t = [tw+t for tw in tws, t in vcat(0,ts)]
        @test sortunique(mcs) == sortunique(tw_plus_t)

        # Specify log scale for t
        mcs,tws,ts = getfields(`./exe/annealer -L -m 84 -i 423 -t $m -T $M -n 15 42 debugfile`)
        @test tws == glasslogspace(1,423,84)
        @test ts  == glasslogspace(m,M,15)
        tw_plus_t = [tw+t for tw in tws, t in vcat(0,ts)]
        @test sortunique(mcs) == sortunique(tw_plus_t)
        # Only a t
        mcs,tws,ts = getfields(`./exe/annealer -m 84 -i 423 -t $m -n 1 42 debugfile`)
        @test tws == glasslogspace(1,423,84)
        @test ts  == [m]
        tw_plus_t = [tw+t for tw in tws, t in vcat(0,ts)]
        @test sortunique(mcs) == sortunique(tw_plus_t)
        # No t
        mcs,tws,ts = getfields(`./exe/annealer -m 84 -i 423 -n 0 42 debugfile`)
        @test tws == glasslogspace(1,423,84)
        @test ts  == [0]
        @test mcs == tws
    end
end

@testset "Output format" begin
    const iters = 10
    for m in 2:iters
        for L in [2,8]
            test_output(L, iters, m)
        end
    end
end

@testset "Invalid input" begin
    # Normaly, return 0
    cmd = `./exe/annealer -i 10 -m 2 42 debugfile`
    io, proc = open(pipeline(cmd))
    wait(proc)
    @test proc.exitcode == 0
    # You can't do more measurements than iters
    cmd = `./exe/annealer -i 20 -m 21 42 debugfile`
    io, proc = open(pipeline(cmd,stderr=DevNull))
    wait(proc)
    @test proc.exitcode != 0
    # But you can measure every iter
    cmd = `./exe/annealer -i 20 -m 20 42 debugfile`
    io, proc = open(pipeline(cmd,stderr=DevNull))
    wait(proc)
    @test proc.exitcode == 0
    # You need to provide at least a temperature
    cmd = `./exe/annealer -i 20 -m 2 debugfile`
    io, proc = open(pipeline(cmd,stderr=DevNull))
    wait(proc)
    @test proc.exitcode != 0
    # You can't provide more than a temperature
    cmd = `./exe/annealer -i 20 -m 2 4 4 debugfile`
    io, proc = open(pipeline(cmd,stderr=DevNull))
    wait(proc)
    @test proc.exitcode != 0
    # It must be a positive T
    cmd = `./exe/annealer -i 20 -m 2 -2 debugfile`
    io, proc = open(pipeline(cmd,stderr=DevNull))
    wait(proc)
    @test proc.exitcode != 0
    # Don't even try to cheat argp
    cmd = `./exe/annealer -i 20 -m 2 -- -2 debugfile`
    io, proc = open(pipeline(cmd,stderr=DevNull))
    wait(proc)
    @test proc.exitcode != 0
    # measurements ≥ 1
    cmd = `./exe/annealer -i 20 -m 1 -- -2 debugfile`
    io, proc = open(pipeline(cmd,stderr=DevNull))
    wait(proc)
    @test proc.exitcode != 0
    # iters ≥ 1
    cmd = `./exe/annealer -i 1 -m 0 -- -2 debugfile`
    io, proc = open(pipeline(cmd,stderr=DevNull))
    wait(proc)
    @test proc.exitcode != 0

    rm("debugfile")
end

@testset "ferromagnet" begin
    run(`./exe/annealer -f -i 20 -m 10 42 debugfile`);
    io = open("debugfile")

    header = map(Char, read(io,Int8,8)) |> join
    is_big_endian = read(io,Int64)
    L = read(io,Int64)
    ntws = read(io,Int64)
    tws = read(io,Int64,ntws)
    nts = read(io,Int64)
    ts = read(io,Int64,nts)
    total_meas = read(io, Int64)
    T = read(io, Float64)

    Jup = read(io,UInt64,L*L*L)
    Jright = read(io,UInt64,L*L*L)
    Jfront = read(io,UInt64,L*L*L)

    # 0 001 ⋯ 001 001 001
    const ONES_BY_TRIPLETS = 0x1249249249249249

    for i in 1:L*L*L # Check J∈0,1
        @test Jup[i]    == ONES_BY_TRIPLETS
        @test Jright[i] == ONES_BY_TRIPLETS
        @test Jfront[i] == ONES_BY_TRIPLETS
    end

    close(io)
    rm("debugfile")
end

@testset "J seeding" begin

    function getJ(seed)
        run(`./exe/annealer -j $seed -i 20 -m 10 42 debugfile`);
        io = open("debugfile")

        header = map(Char, read(io,Int8,8)) |> join
        is_big_endian = read(io,Int64)
        L = read(io,Int64)
        ntws = read(io,Int64)
        tws = read(io,Int64,ntws)
        nts = read(io,Int64)
        ts = read(io,Int64,nts)
        total_meas = read(io, Int64)
        T = read(io, Float64)

        JupA = read(io,UInt64,L*L*L)
        JrightA = read(io,UInt64,L*L*L)
        JfrontA = read(io,UInt64,L*L*L)

        close(io)
        rm("debugfile")

        return JupA, JrightA, JfrontA
    end

    seed = rand(UInt32)

    @test getJ(seed) == getJ(seed) == getJ(seed) == getJ(seed)
    @test getJ(seed) != getJ(seed*2) != getJ(seed*4) != getJ(seed*5)

end

max_stars_repo_name: RedPointyJackson/tfg

max_issues_repo_name: RedPointyJackson/tfg

max_forks_repo_name: RedPointyJackson/tfg


text:
using JSON3
using ProgressBars

function distill(path::String, datafilter::RedditDataFilter)
    file = getfile(path)
    data_vector = Vector{Dict{Symbol,Any}}
    multithreaded = Threads.nthreads() > 1
    values = multithreaded ? [data_vector() for t = 1:Threads.nthreads()] : data_vector()
    decompressor_stream = get_decompressor_stream(file)
    open(decompressor_stream, file) do stream
        bar = ProgressBar(eachline(stream))
        bar.description = "File: $(basename(file)); Threads: $(Threads.nthreads())"
        for line in bar
            if multithreaded
                Threads.@spawn begin
                    if contains(line, datafilter)
                        linedata = JSON3.read(line)
                        if contains(linedata, datafilter)
                            data = convert(Dict, linedata)
                            data = filter(datafilter, data)
                            fetched_data = fetch(data)
                            push!(values[Threads.threadid()], fetched_data)
                        end
                    end
                end
            else
                if contains(line, datafilter)
                    linedata = JSON3.read(line)
                    if contains(linedata, datafilter)
                        data = convert(Dict, linedata)
                        data = filter(datafilter, data)
                        push!(values, data)
                    end
                end
            end
        end
    end
    return multithreaded ? collect(Iterators.flatten(values)) : values
end


text:
using DelayEmbeddings, DynamicalSystemsBase
using Test
using Random

@testset "Pecora" begin
# %% Generate data
lor = Systems.lorenz([0.0;1.0;0.0];ρ=60)
data = trajectory(lor, 200; Δt=0.02, Ttr = 10)
metric = Chebyshev()

@testset "Pecora univariate" begin
# %% Timeseries case

    s = data[:, 2] # input timeseries = first entry of lorenz
    optimal_τ = estimate_delay(s, "mi_min")
    Tmax = 100
    K = 14
    samplesize = 1

    τs = (0,)
    Random.seed!(123)
    es_ref, Γs = pecora(s, τs; delays = 0:Tmax, w = optimal_τ, samplesize = samplesize, K = K, metric = metric)
    (max1,_) = DelayEmbeddings.findlocalextrema(vec(es_ref))
    @test optimal_τ - 2 ≤ max1[1]-1 ≤ optimal_τ + 2
    maxi = maximum(es_ref)
    @test maxi ≤ 1.3

    Random.seed!(123)
    es_ref2, Γs = pecora(s, τs; delays = 0:Tmax, w = optimal_τ, K = K, metric = metric)
    (max2,_) = DelayEmbeddings.findlocalextrema(vec(es_ref2))
    @test optimal_τ - 2 ≤ max2[1]-1 ≤ optimal_τ + 2
    maxi2 = maximum(es_ref2)
    @test maxi - .1 ≤ maxi2 ≤ maxi + .1

    τs = (0, max1[1]-1,)
    Random.seed!(123)
    es, Γs = pecora(s, τs; delays = 0:Tmax, w = optimal_τ, samplesize = samplesize, K = K, metric = metric)
    (max2,min1) = DelayEmbeddings.findlocalextrema(vec(es))
    @test 1 ≤ max2[1]-1 ≤ 4
    @test optimal_τ - 2 ≤ min1[2]-1 ≤ optimal_τ + 2

    τs = (0, max1[1]-1, max2[1]-1)
    Random.seed!(123)
    es, Γs = pecora(s, τs; delays = 0:Tmax, w = optimal_τ, samplesize = samplesize, K = K, metric = metric)
    (_,min2) = DelayEmbeddings.findlocalextrema(vec(es))
    @test max2[1]-2 ≤ min2[1]-1 ≤ max2[1]

end

@testset "Pecora multivariate" begin
## %% Trajectory case
    s = Dataset(data)
    optimal_τ = estimate_delay(s[:,2], "mi_min")
    Tmax = 100
    K = 14
    samplesize = 1

    js = (2,)
    τs = (0,)
    Random.seed!(123)
    es_ref, _ = pecora(s[:,2], τs; delays = 0:Tmax, w = optimal_τ, samplesize = samplesize, K = K, metric = metric)
    Random.seed!(123)
    es_ref2, _ = pecora(s, τs, js; delays = 0:Tmax, w = optimal_τ, samplesize = .2, K = K, metric = metric)
    Random.seed!(123)
    es, _= pecora(s, τs, js; delays = 0:Tmax, samplesize = samplesize, w = optimal_τ, K = K, metric = metric)

    @test round.(es[:,2], digits = 4) == round.(vec(es_ref),digits = 4)
    (x_maxi,_) = DelayEmbeddings.findlocalextrema(es[:,1])
    @test 8 ≤ x_maxi[2]-1 ≤ 10
    (z_maxi,_) = DelayEmbeddings.findlocalextrema(es[:,3])
    @test 13 ≤ z_maxi[2]-1 ≤ 15
    (x_maxi_ref,_) = DelayEmbeddings.findlocalextrema(es[:,2])
    (x_maxi_ref2,_) = DelayEmbeddings.findlocalextrema(es_ref2[:,2])

    @test x_maxi_ref[1:4] == x_maxi_ref2[1:4]
    @test es[x_maxi_ref[1:4],2] .- .1 ≤ es_ref2[x_maxi_ref2[1:4],2] ≤ es[x_maxi_ref[1:4],2] .+ .1

end
end


text:
@testset "util" begin

@testset "sort_ind!" begin
    v = MArray{Tuple{2,2,2}}(rand(8))
    ind = MArray{Tuple{2,2,2}}(collect(1:8))
    @inferred(MaxwellBase.sort_ind!(ind,v))
    @test issorted(v[ind])
    @test @inferred(MaxwellBase.countdiff(ind,v)) == (8,8)
end  # @testset "sort8!"

end  # @testset "util"


text:
import PhysicsEngine

@testset "Vector2" begin
  v1 = PhysicsEngine.Vector2((5.0, 7.0))
  v2 = PhysicsEngine.Vector2((-3.5, 2.5))

  @test v1 + v2 == PhysicsEngine.Vector2((1.5, 9.5))
  @test v1 - v2 == PhysicsEngine.Vector2((8.5, 4.5))
  @test PhysicsEngine.dot(v1, v2) == -17.5 + 17.5
  @test PhysicsEngine.magnitude(v1) == sqrt(74)
end

@testset "Vector3" begin
  v1 = PhysicsEngine.Vector3((1.0, 5.0, -2.0))
  v2 = PhysicsEngine.Vector3((5.5, -4.0, -0.5))

  @test v1 + v2 == PhysicsEngine.Vector3((6.5, 1.0, -2.5))
  @test v1 - v2 == PhysicsEngine.Vector3((-4.5, 9.0, -1.5))
  @test PhysicsEngine.dot(v1, v2) == 5.5 + -20.0 + 1.0
  @test v1 * v2 == PhysicsEngine.Vector3((-10.5, -10.5, -31.5))
  @test PhysicsEngine.magnitude(v2) == sqrt(46.5)
end

@testset "Vector4" begin
  v1 = PhysicsEngine.Vector4((13, 14, 13, 17))
  v2 = PhysicsEngine.Vector4((3, 4, 5, 6))

  @test v1 + v2 == PhysicsEngine.Vector4((16, 18, 18, 23))
  @test v1 - v2 == PhysicsEngine.Vector4((10, 10, 8, 11))
  @test PhysicsEngine.dot(v1, v2) == 39 + 56 + 65 + 102
  @test PhysicsEngine.magnitude(v2) == sqrt(86)
end

@testset "Matrix" begin
  m1 = PhysicsEngine.Matrix4x4((2, 9, 6, 7, 6, 1, 4, 1, 1, 1, 4, 8, 4, 3, 8, 6))
  m2 = PhysicsEngine.Matrix4x4((6, 7, 9, 3, 4, 6, 3, 9, 3, 9, 8, 4, 8, 5, 9, 4))
  v1 = PhysicsEngine.Vector4((14.0, 24.0, 13.0, 19.0))
  v2 = PhysicsEngine.Vector4((3, 7, 1, 1))

  function isapproxmatrix(x::PhysicsEngine.Matrix4x4, y::PhysicsEngine.Matrix4x4)
    for pair in zip(x, y)
      xValue = pair[1]
      yValue = pair[2]
      if !(xValue ≈ yValue)
        return false
      end
    end
    return true
  end

  @test m1 * m2 == PhysicsEngine.Matrix4x4((122, 157, 156, 139, 60, 89, 98, 47, 86, 89, 116, 60, 108, 148, 163, 95))
  @test m1 * v1 == PhysicsEngine.Vector4((455, 179, 242, 346))
  @test PhysicsEngine.identityMatrix() == PhysicsEngine.Matrix4x4((1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1))
  @test PhysicsEngine.translationMatrix(PhysicsEngine.vec3(2.0, 3.0, -1.0)) == PhysicsEngine.Matrix4x4((1.0, 0.0, 0.0, 2.0, 0.0, 1.0, 0.0, 3.0, 0.0, 0.0, 1.0, -1.0, 0.0, 0.0, 0.0, 1.0))
  @test PhysicsEngine.translationMatrix(PhysicsEngine.vec3(2.5, -2.5, 4.0)) * v2 == PhysicsEngine.Vector4((5.5, 4.5, 5.0, 1.0))
  @test PhysicsEngine.scalingMatrix(PhysicsEngine.vec3(1.0, 2.0, 0.5)) == PhysicsEngine.Matrix4x4((1.0, 0.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 1.0))
  @test PhysicsEngine.scalingMatrix(PhysicsEngine.vec3(2.0, -2.0, 5.5)) * v2 == PhysicsEngine.Vector4((6.0, -14.0, 5.5, 1.0))
  @test isapproxmatrix(PhysicsEngine.rotationXMatrix(45.0),  PhysicsEngine.Matrix4x4((1.0, 0.0, 0.0, 0.0, 0.0, cos(-45.0), -sin(-45.0), 0.0, 0.0, sin(-45.0), cos(-45.0), 0.0, 0.0, 0.0, 0.0, 1.0)))
  @test isapproxmatrix(PhysicsEngine.rotationYMatrix(30.0), PhysicsEngine.Matrix4x4((cos(-30.0), 0.0, sin(-30.0), 0.0, 0.0, 1.0, 0.0, 0.0, -sin(-30.0), 0.0, cos(-30.0), 0.0, 0.0, 0.0, 0.0, 1.0)))
  @test isapproxmatrix(PhysicsEngine.rotationZMatrix(22.0), PhysicsEngine.Matrix4x4((cos(-22.0), -sin(-22.0), 0.0, 0.0, sin(-22.0), cos(-22.0), 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0)))

  #println(PhysicsEngine.invertMatrix(m1))
  #@test isapproxmatrix(PhysicsEngine.invertMatrix(m1))
  # @test rotating vectors
end

max_stars_repo_name: jonathanbuchanan/PhysicsEngine

max_issues_repo_name: jonathanbuchanan/PhysicsEngine

max_forks_repo_name: jonathanbuchanan/PhysicsEngine


text:
using JSON
using OrderedCollections: OrderedDict
using DataStructures
using Printf
using StaticArrays
import YAML

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

logData = false

function ReadDataSet(filename::String)
    if logData
        printstyled("     Datafile: ", color=:green)
        printstyled(filename,"\n", color=:white)
    end

    if uppercase(filename[end-3:end]) == "JSON"
        f = open(filename, "r")
        dicttxt = read(f,String)
        close(f)
        DataDict=JSON.parse(dicttxt, dicttype=DataStructures.OrderedDict)
    elseif uppercase(filename[end-3:end]) == "YAML"
        DataDict = YAML.load_file(filename)
    else
        error("AUS und AUS")
    end

    setDict = OrderedDict{String, Any}()

    # OrderedDict für Reihenfolge wie in der Datei
    isZONE = DataDict["description"]["type"]=="ZONEdata"
    for (key, value) in DataDict
        if key=="description"; continue;    end
        valueDict = OrderedDict{Symbol, Any}()
        setDict[key] = valueDict
        for (key1, value1) in value
            valueDict[Symbol(key1)] = value1
        end
        if logData
            printstyled("  data set ID: ", color=:cyan)
            printstyled(key,"  ", color=:white)
            printstyled("(Prj: ",valueDict[:Project],")\n", color=:light_black)
        end
#         println()
        if isZONE
             valueDict[:Omegarated] = valueDict[:frated]*2.0*pi
             valueDict[:Speedrad] = valueDict[:Omegarated]/valueDict[:NoPolepairs]
             valueDict[:Irated] = valueDict[:Srated]/valueDict[:Vrated]/sqrt(3.0)
             valueDict[:Zbase] = valueDict[:Vrated]*valueDict[:Vrated]/valueDict[:Srated]
             valueDict[:Tbase] = valueDict[:Srated]/valueDict[:Speedrad]
             valueDict[:Speedrpm] = valueDict[:frated]*60.0/valueDict[:NoPolepairs]
        end
        valueDict[:_class] = :Map
    end
    println()
    return setDict
end

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

function Ergebnis(tttt, liste::Vector{String})
    printstyled(lpad("[signal name]",21), color=:light_green)
    printstyled("      [1]", color=:green)
    printstyled("    [end]", color=:cyan)
    printstyled("     [min : max] \n", color=:light_yellow)
    for w in liste
        if !haskey(tttt.result_info, w) && !haskey(tttt.evaluatedParameters, Symbol(w))
            printstyled(lpad(w,21),"\n", color=:red)
            continue
        end
        daten = get_result(tttt,w)
        L1 = length(daten)
        L2 = length(daten[1])
        RES = Array{Float64,2}(zeros(L1,L2))    # für MIN/MAX notwendig.
        for i = 1:L1
            for j = 1:L2
                RES[i,j] = daten[i][j]
            end
        end

        if L2 == 1
            printstyled(lpad(w,21), color=:light_green)
            printstyled(@sprintf("  %7.4g",RES[1]), color=:green)
            printstyled(@sprintf("  %7.4g",RES[end]), color=:cyan)
            printstyled(@sprintf("  %7.4g :%7.4g\n",minimum(RES),maximum(RES)), color=:light_yellow)
        else
            for iz =1:L2
                printstyled(lpad(w*"[$iz]",21), color=:light_green)
                printstyled(@sprintf("  %7.4g",RES[1,iz]), color=:green)
                printstyled(@sprintf("  %7.4g",RES[end,iz]), color=:cyan)
                printstyled(@sprintf("  %7.4g :%7.4g\n",minimum(RES,dims=1)[iz],maximum(RES,dims=1)[iz]), color=:light_yellow)
            end

        end
    end
end

function ergebnis(tttt, w::String)
    if length(w)==0
    else
        daten = get_result(tttt,w)
        datenlg = length(daten[1])
        if datenlg == 1
            printstyled(lpad(w,21), color=:light_green)
            printstyled(@sprintf("  %7.4g",daten[1]), color=:green)
            printstyled(@sprintf("  %7.4g",daten[end]), color=:cyan)
            printstyled(@sprintf("  %7.4g :%7.4g\n",minimum(daten),maximum(daten)), color=:light_yellow)
        else
            for iz =1:datenlg
                printstyled(lpad(w*"[$iz]",21), color=:light_green)
                printstyled(@sprintf("  %7.4g",daten[1][iz]), color=:green)
                printstyled(@sprintf("  %7.4g",daten[end][iz]), color=:cyan)
                printstyled("      coming soon\n", color=:light_black)
#                 printstyled(@sprintf("  %7.4g :%7.4g\n",minimum(daten[:][iz]),maximum(daten[:][iz])), color=:light_yellow)
            end

        end
    end
end

function Speichern(mod, liste::AbstractArray= [])

    datei = "myfile.txt"
    if length(liste)==0
        printstyled("Länge der Liste = 0!!\n", color=:red)
        return
    end
    df=get_result(mod)
    open(datei, "w") do io
        show(io, MIME("text/csv"), df[!,liste])
        printstyled(@sprintf("save CSV in file: %s with %d colums, %s rows\n",datei, length(liste),nrow(df)), color=:green)
    end
end

function TemplateInit(instModel; tempName="Start.Template.jl")::Nothing

    function StartWert(wert::Float64)::String
        return string(wert)
    end

    function StartWert(wert::SVector{2, Float64})::String   #FIXME Länge 2 ist angenommen
        w = "["
        lg = length(wert.data)
        for i= 1:lg
            w = w * string(wert.data[i])
            if i<lg
                w = w * ", "
            end
        end
        w = w * "]"
        return w
    end

    Etxt = "EIN = (\n"
    Atxt = "AUS = (\n"

    for Z in instModel.equationInfo.x_info

        TXT1 = "    (\"$(Z.x_name)\", " * StartWert(Z.startOrInit) * " ),\n"
        if Z.scalar
            TXT2 = "    (\"$(Z.der_x_name)\", 0.0),\n"
        elseif Z.length==2
            TXT2 = "    (\"$(Z.der_x_name)\", [0.0, 0.0]),\n"
        else
            printstyled("ERROR - $(Z_xname)  length >2\n", color=:red)
        end
        Etxt = Etxt * TXT1
        Atxt = Atxt * TXT2
    end

    Etxt = Etxt * ")\n\n\n"
    Atxt = Atxt * ")\n"

    Ausgabe = ""

    open(tempName, "w") do io
        print(io, Etxt)
        print(io, Atxt)
    end

    nothing

end

max_stars_repo_name: johhell/EPSphasor

max_issues_repo_name: johhell/EPSphasor

max_forks_repo_name: johhell/EPSphasor