text:
using BrainFlow
using BrainFlowViz

function get_some_board_data(board_shim, nsamples)
    data = BrainFlow.get_current_board_data(nsamples, board_shim)
    data = transpose(data)
    return view(data, :, 2:9)
end

nchannels = 8
nsamples = 256

### Start streaming
BrainFlow.enable_dev_logger(BrainFlow.BOARD_CONTROLLER)
params = BrainFlowInputParams()
board_shim = BrainFlow.BoardShim(BrainFlow.SYNTHETIC_BOARD, params)
BrainFlow.prepare_session(board_shim)
BrainFlow.start_stream(board_shim)

# brief sleep
sleep(0.5)

data_func = ()->get_some_board_data(board_shim, nsamples)

# start a live plotting task
t = @task BrainFlowViz.plot_data(
    data_func,
    nsamples,
    nchannels,
    theme = :dark,
    color = :lime,
    delay = 0.02,
    )

sleep(0.5)
schedule(t)

sleep(4)

### to stop the task:
schedule(t, InterruptException(), error=true)

# this should go into a 'finally' of try/catch/finally
BrainFlow.stop_stream(board_shim)
BrainFlow.release_session(board_shim)


text:
using DelimitedFiles
using DataFrames
using Plots

gr(size = (200, 200))

function read_files()
  directory = "dados"
  dados = Array{Float64, 2}(undef,500,20)

  for i in 1:20
    fn = "dados_ "*string(i)
    dados[:,i] = readdlm(directory*"/"*fn, header=true, Float64)[1][:,2]
  end

  return dados
end

dados = read_files()
df = DataFrame(dados, :auto)
dados_label = reshape(string.("dados ", 1:20), (1, 20))

histogram(dados[:,1:10], bins=20, layout = (5,2), title = dados_label[:,1:10], legend = false)
savefig("hist_1-10.png")

ext: jl


text:
# series

# Makie

using CairoMakie, Random, LaTeXStrings
CairoMakie.activate!()
let
    Random.seed!(123)
    m = 200
    function spiral(; a = 1, n = 100, h = 0, k = 0)
        φ = LinRange(rand() + 1, 6π, n)
        h .+ rand(-1:2:1) * a * cos.(φ) ./ φ, k .+ rand(-1:2:1) * a * sin.(φ) ./ φ
    end
    x = raw"x=a\,\cos(\varphi)/\varphi,\quad "
    y = raw"y=a\,\sin(\varphi)/\varphi,\quad "
    curves = [spiral(; a = rand(), h = rand(-1:1)) for i in 1:m]

    with_theme(theme_minimal()) do
        fig = Figure(resolution = (600, 400))
        ax = Axis(fig[1, 1])
        CairoMakie.series!(curves; color = to_colormap(:Hokusai1, m), linewidth = 1.5)
        text!(latexstring(x * y * "\\varphi>0"), position = (-1, -1))
        hidedecorations!(ax; grid = false)
        display(fig)
    end
end;

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
const DATAPATH = joinpath(@__DIR__, "data")

@testset "Trajectory" begin
    @testset "Errors handling" begin
        remove_chemfiles_warning() do
            @test_throws ChemfilesError Trajectory(joinpath(DATAPATH, "not-here.xyz"))
            @test_throws ChemfilesError Trajectory(joinpath(DATAPATH, "empty.unknown"))

            trajectory = Trajectory(joinpath(DATAPATH, "water.xyz"))
            @test_throws ChemfilesError take!(trajectory)
            close(trajectory)
        end
    end

    @testset "Read frames" begin
        trajectory = Trajectory(joinpath(DATAPATH, "water.xyz"))

        @test path(trajectory) == joinpath(DATAPATH, "water.xyz")

        @test length(trajectory) == 100

        frame = read(trajectory)

        @test size(frame) == 297
        pos = positions(frame)
        @test pos[:, 1] ≈ Float64[0.417219, 8.303366, 11.737172]
        @test pos[:, 125] ≈ Float64[5.099554, -0.045104, 14.153846]

        read!(trajectory,frame)
        @test size(frame) == 297
        pos = positions(frame)
        @test pos[:, 1] ≈ Float64[0.357720, 8.306397, 11.744901]
        @test pos[:, 125] ≈ Float64[5.007032, 0.066848, 14.184831]

        topology = Topology(frame)
        @test size(frame) == 297
        @test name(Atom(frame, 0)) == "O"
        @test name(Atom(frame, 1)) == "H"

        set_cell!(trajectory, UnitCell([30.0, 30.0, 30.0]))
        frame = read_step(trajectory, 41)

        @test lengths(UnitCell(frame)) == [30.0, 30.0, 30.0]

        pos = positions(frame)
        @test pos[:, 1] ≈ Float64[0.761277, 8.106125, 10.622949]
        @test pos[:, 125] ≈ Float64[5.13242, 0.079862, 14.194161]

        read_step!(trajectory,41,frame)

        @test lengths(UnitCell(frame)) == [30.0, 30.0, 30.0]

        pos = positions(frame)
        @test pos[:, 1] ≈ Float64[0.761277, 8.106125, 10.622949]
        @test pos[:, 125] ≈ Float64[5.13242, 0.079862, 14.194161]

        topology = Topology(frame)
        @test size(topology) == 297
        @test bonds_count(topology) == 0

        guess_bonds!(frame)
        topology = Topology(frame)
        @test bonds_count(topology) == 180
        @test angles_count(topology) == 84

        topology = Topology()
        a = Atom("Cs")
        for i=1:297
            add_atom!(topology, a)
        end

        set_topology!(trajectory, topology)
        frame = read_step(trajectory, 10)
        @test name(Atom(frame, 10)) == "Cs"

        set_topology!(trajectory, joinpath(DATAPATH, "topology.xyz"))
        frame = read(trajectory)
        @test name(Atom(frame, 100)) == "Rd"
    end

    @testset "Iteration" begin
        trajectory = Trajectory(joinpath(DATAPATH, "water.xyz"))

        count = 0
        for frame in trajectory
            count += 1
        end

        @test count == length(trajectory)
    end

    @testset "Write frames" begin
        EXPECTED_CONTENT = """4
                              Properties=species:S:1:pos:R:3
                              X 1 2 3
                              X 1 2 3
                              X 1 2 3
                              X 1 2 3
                              """

        frame = Frame()
        for i=1:4
            add_atom!(frame, Atom("X"), Float64[1, 2, 3])
        end

        trajectory = Trajectory("test-tmp.xyz", 'w');
        write(trajectory, frame)
        close(trajectory)
        @test isopen(trajectory) == false

        s = open("test-tmp.xyz") do file
            read(file, String)
        end

        @test s == EXPECTED_CONTENT

        rm("test-tmp.xyz")
    end

    @testset "Function syntax" begin
        trajectory = nothing
        Trajectory(joinpath(DATAPATH, "water.xyz")) do tr
            trajectory = tr
            @test isopen(trajectory)
        end
        @test !isopen(trajectory)
    end

    @testset "In-memory reader" begin
        trajectory = MemoryTrajectory("XYZ", """3

H 0 0 0
C 1 2 0
Zn 3 2 -5""")
        frame = read(trajectory)
        expected = transpose([
            0 0 0
            1 2 0
            3 2 -5
        ])
        @test positions(frame) == expected
        @test name(Atom(frame, 0)) == "H"
        @test name(Atom(frame, 1)) == "C"
        @test name(Atom(frame, 2)) == "Zn"
    end

    @testset "In-memory writer" begin
        frame = Frame()
        add_atom!(frame, Atom("Cu"), Float64[2, 1, 4])
        add_atom!(frame, Atom("Hf"), Float64[0, 3.2, 0])
        add_atom!(frame, Atom("Li"), Float64[1, 0, -4])

        trajectory = MemoryTrajectory("PDB")
        write(trajectory, frame)
        expected = """MODEL    1
CRYST1    0.000    0.000    0.000  90.00  90.00  90.00 P 1           1
HETATM    1 Cu           1       2.000   1.000   4.000  1.00  0.00          Cu
HETATM    2 Hf           2       0.000   3.200   0.000  1.00  0.00          Hf
HETATM    3 Li           3       1.000   0.000  -4.000  1.00  0.00          Li
ENDMDL
"""

        @test String(take!(trajectory)) == expected
    end
end


text:
## SparseTensorCOO: coordinate format

"""
    SparseTensorCOO{Tv,Ti<:Integer,N} <: AbstractSparseTensor{Tv,Ti,N}

`N`-dimensional sparse tensor stored in the **COO**rdinate format.
Elements are stored as a vector of indices (using type `Ti`)
and a vector of values (of type `Tv`).

Fields:
+ `dims::Dims{N}`              : tuple of dimensions
+ `inds::Vector{NTuple{N,Ti}}` : vector of indices
+ `vals::Vector{Tv}`           : vector of values
"""
struct SparseTensorCOO{Tv,Ti<:Integer,N} <: AbstractSparseTensor{Tv,Ti,N}
    dims::Dims{N}                   # Dimensions
    inds::Vector{NTuple{N,Ti}}      # Stored indices
    vals::Vector{Tv}                # Stored values

    function SparseTensorCOO{Tv,Ti,N}(dims::Dims{N}, inds::Vector{NTuple{N,Ti}},
                            vals::Vector{Tv}) where {Tv,Ti<:Integer,N}
        check_Ti(dims, Ti)
        check_coo_buffers(inds, vals)
        check_coo_inds(dims, inds)
        return new(dims, inds, vals)
    end
end
function SparseTensorCOO(dims::Dims{N}, inds::Vector{NTuple{N,Ti}},
                        vals::Vector{Tv}) where {Tv,Ti<:Integer,N}
    if issorted(inds; by = reverse)
        _inds = inds
        _vals = vals
    else
        perm = sortperm(inds; by = reverse)
        _inds = inds[perm]
        _vals = vals[perm]
    end
    SparseTensorCOO{Tv,Ti,N}(dims, _inds, _vals)
end

"""
    SparseTensorCOO{Tv,Ti<:Integer}(undef, dims)
    SparseTensorCOO{Tv,Ti<:Integer,N}(undef, dims)

Construct an uninitialized `N`-dimensional `SparseTensorCOO`
with indices using type `Ti` and elements of type `Tv`.
Here uninitialized means it has no stored entries.

Here `undef` is the `UndefInitializer`. If `N` is supplied,
then it must match the length of `dims`.

# Examples
```julia-repl
julia> A = SparseTensorCOO{Float64, Int8, 3}(undef, (2, 3, 4)) # N given explicitly
2×3×4 SparseTensorCOO{Float64, Int8, 3} with 0 stored entries

julia> B = SparseTensorCOO{Float64, Int8}(undef, (4,)) # N determined by the input
4-element SparseTensorCOO{Float64, Int8, 1} with 0 stored entries

julia> similar(B, 2, 4, 1) # use typeof(B), and the given size
2×4×1 SparseTensorCOO{Float64, Int8, 3} with 0 stored entries
```
"""
SparseTensorCOO{Tv,Ti,N}(::UndefInitializer, dims::Dims{N}) where {Tv,Ti<:Integer,N} =
    SparseTensorCOO(dims, Vector{NTuple{N,Ti}}(), Vector{Tv}())
SparseTensorCOO{Tv,Ti}(::UndefInitializer, dims::Dims{N}) where {Tv,Ti<:Integer,N} =
    SparseTensorCOO{Tv,Ti,N}(undef, dims)

"""
    SparseTensorCOO(Ti, A::AbstractArray)

Convert an AbstractArray `A` into a `SparseTensorCOO`
with indices using type `Ti`.

# Examples
```julia-repl
julia> A = SparseTensorCOO(Int8, Float16[1.1 0.0 0.0; 2.1 0.0 2.3])
2×3 SparseTensorCOO{Float16, Int8, 2} with 3 stored entries:
  [1, 1]  =  1.1
  [2, 1]  =  2.1
  [2, 3]  =  2.3

julia> B = SparseTensorCOO(Int16, Float16[1.1 0.0 0.0; 2.1 0.0 2.3])
2×3 SparseTensorCOO{Float16, Int16, 2} with 3 stored entries:
  [1, 1]  =  1.1
  [2, 1]  =  2.1
  [2, 3]  =  2.3
```
"""
function SparseTensorCOO(Ti::Type{<:Integer}, A::AbstractArray)
    Tv, N = eltype(A), ndims(A)
    dims = size(A)
    nzidx = findall(!iszero, A)
    inds = convert(Vector{NTuple{N,Ti}}, CartesianIndices(A)[nzidx])
    vals = convert(Vector{Tv}, A[nzidx])
    SparseTensorCOO{Tv,Ti,N}(dims, inds, vals)
end

## Minimal AbstractArray interface

size(A::SparseTensorCOO) = A.dims

function getindex(A::SparseTensorCOO{Tv,<:Integer,N}, I::Vararg{Int,N}) where {Tv,N}
    @boundscheck checkbounds(A, I...)
    ptr = searchsortedfirst(A.inds, I; by = reverse)
    return (ptr == length(A.inds) + 1 || A.inds[ptr] != I) ? zero(Tv) : A.vals[ptr]
end

function setindex!(A::SparseTensorCOO{Tv,Ti,N}, v, I::Vararg{Int,N}) where {Tv,Ti<:Integer,N}
    @boundscheck checkbounds(A, I...)
    ind, val = convert(NTuple{N,Ti}, I), convert(Tv, v)
    ptr = searchsortedfirst(A.inds, ind; by = reverse)
    if ptr == length(A.inds) + 1 || A.inds[ptr] != ind
        if !iszero(val)
            insert!(A.inds, ptr, ind)
            insert!(A.vals, ptr, val)
        end
    else
        A.vals[ptr] = val
    end
    return A
end

IndexStyle(::Type{<:SparseTensorCOO}) = IndexCartesian()

## Overloads for specializing outputs

similar(::SparseTensorCOO{<:Any,Ti}, ::Type{Tv}, dims::Dims{N}) where {Tv,Ti<:Integer,N} =
    SparseTensorCOO{Tv,Ti,N}(undef, dims)

## Overloads for improving efficiency

# technically specializes the output since the state is different
function iterate(A::SparseTensorCOO{Tv}, state=((eachindex(A),),1)) where {Tv}
    idxstate, nextptr = state
    y = iterate(idxstate...)
    y === nothing && return nothing
    if nextptr > length(A.inds) || A.inds[nextptr] != Tuple(y[1])
        val = zero(Tv)
    else
        val = A.vals[nextptr]
        nextptr += 1
    end
    val, ((idxstate[1], Base.tail(y)...), nextptr)
end

## AbstractSparseTensor interface

function dropstored!(f::Function, A::SparseTensorCOO)
    ptrs = findall(f, A.vals)
    deleteat!(A.inds, ptrs)
    deleteat!(A.vals, ptrs)
    return A
end

numstored(A::SparseTensorCOO) = length(A.vals)
storedindices(A::SparseTensorCOO) = A.inds
storedvalues(A::SparseTensorCOO) = A.vals
storedpairs(A::SparseTensorCOO) = Iterators.map(Pair, A.inds, A.vals)

## AbstractSparseTensor optional interface (internal)

findall_stored(f::Function, A::SparseTensorCOO) =
    [convert(keytype(A), CartesianIndex(ind)) for (ind, val) in storedpairs(A) if f(val)]

## Utilities

"""
    check_coo_buffers(inds, vals)

Check that the `inds` and `vals` buffers are valid:
+ their lengths match (`length(inds) == length(vals)`)
If not, throw an `ArgumentError`.
"""
function check_coo_buffers(inds::Vector{NTuple{N,Ti}}, vals::Vector{Tv}) where {Tv,Ti<:Integer,N}
    length(inds) == length(vals) ||
        throw(ArgumentError("the buffer lengths (length(inds) = $(length(inds)), length(vals) = $(length(vals))) do not match"))
    return nothing
end

"""
    check_coo_inds(dims, inds)

Check that the indices in `inds` are valid:
+ each index is in bounds (`1 ≤ inds[ptr][k] ≤ dims[k]`)
+ the indices are sorted (`issorted(inds; by=CartesianIndex)`)
+ the indices are all unique (`allunique(inds`)
If not, throw an `ArgumentError`.
"""
function check_coo_inds(dims::Dims{N}, inds::Vector{NTuple{N,Ti}}) where {Ti<:Integer,N}
    # Check all the conditions in a single pass over inds for efficiency
    itr = iterate(inds)
    itr === nothing && return nothing
    prevind, state = itr
    checkbounds_dims(dims, prevind...)
    itr = iterate(inds, state)
    while itr !== nothing
        thisind, state = itr
        if reverse(prevind) < reverse(thisind)
            checkbounds_dims(dims, thisind...)
        elseif reverse(prevind) > reverse(thisind)
            throw(ArgumentError("inds are not sorted"))
        else
            throw(ArgumentError("inds are not all unique"))
        end
        prevind = thisind
        itr = iterate(inds, state)
    end
    return nothing
end


text:
module ConfusionTest

using Base.Test
using ForwardDiff

using Base.Test

# Perturbation Confusion (Issue #83) #
#------------------------------------#

D = ForwardDiff.derivative

@test D(x -> x * D(y -> x + y, 1), 1) == 1
@test ForwardDiff.gradient(v -> sum(v) * D(y -> y * norm(v), 1), [1]) == ForwardDiff.gradient(v -> sum(v) * norm(v), [1])

const A = rand(10,8)
y = rand(10)
x = rand(8)

@test A == ForwardDiff.jacobian(x) do x
    ForwardDiff.gradient(y) do y
        dot(y, A*x)
    end
end

# Issue #238                         #
#------------------------------------#

m,g = 1, 9.8
t = 1
q = [1,2]
q̇ = [3,4]
L(t,q,q̇) = m/2 * dot(q̇,q̇) - m*g*q[2]

∂L∂q̇(L, t, q, q̇) = ForwardDiff.gradient(a->L(t,q,a), q̇)
Dqq̇(L, t, q, q̇) = ForwardDiff.jacobian(a->∂L∂q̇(L,t,a,q̇), q)
@test Dqq̇(L, t, q, q̇)  == zeros(2,2)

q = [1,2]
p = [5,6]
function Legendre_transformation(F, w)
    z = zeros(w)
    M = ForwardDiff.hessian(F, z)
    b = ForwardDiff.gradient(F, z)
    v = cholfact(M)\(w-b)
    dot(w,v) - F(v)
end
function Lagrangian2Hamiltonian(Lagrangian, t, q, p)
    L = q̇ -> Lagrangian(t, q, q̇)
    Legendre_transformation(L, p)
end

Lagrangian2Hamiltonian(L, t, q, p)
@test ForwardDiff.gradient(a->Lagrangian2Hamiltonian(L, t, a, p), q) == [0.0,g]

end # module


text:
function mark_frame_perimetr!(r::Robot)
    num_vert,num_hor = ugol(r)
    for side in (Nord,Ost,Sud,West)
        while isborder(r,side)==false
            putmarkers!(r,side)
        end
    end
    moves!(r, Ost, num_hor)
    moves1!(r, Nord, num_vert)
end
function moves!(r::Robot,side::HorizonSide)
    num_steps=0
    while isborder(r,side)==false
        move!(r,side)
        num_steps+=1
    end
    return num_steps
end
function moves!(r::Robot,side::HorizonSide,num_steps::Int)
    for _ in 1:num_steps
        move!(r,side)
    end
end
function moves1!(r::Robot,side::HorizonSide,num_steps::Int)
    q=0
    while q<num_steps
        if isborder(r,Nord)==true
            while isborder(r,side)==true
                move!(r,Ost)
            end
            move!(r,Nord)
            q+=1
            if isborder(r,West)==true
                while isborder(r,West)==true
                    move!(r,Nord)
                    q+=1
                end
            end
            move!(r,West)
        end
        if q<num_steps
        move!(r,side)
        q+=1
        end
    end
end
function putmarkers!(r::Robot, side::HorizonSide)
    while isborder(r,side)==false
        move!(r,side)
        putmarker!(r)
    end
    if isborder(r,side)==true
        putmarker!(r)
    end
end
function ugol(r::Robot)
    num_vert=0
    num_hor=0
    while (isborder(r,Sud)==false || isborder(r,West)==false)
        num_vert+= moves!(r, Sud)
        num_hor+= moves!(r, West)
        repeat
    end
    return num_vert,num_hor
end


text:
using .TestUtils

@testgroup "reporting" begin

TESTCASES = Jute.collect_testobjs() do
    @testcase "multiple tests" begin
        @test 1 == 1
        @test 2 == 2
        @test 3 == 3
    end

    @testcase "returning value" begin
        @test 1 == 1
        @test_result 10
    end

    @testcase "multiple tests and one failure" begin
        @test 1 == 1
        @test 2 == 1
        @test 3 == 3
    end

    @testcase "uncaught exception" begin
        @test 1 == 1
        error("Uncaught exception")
        @test 1 == 1
    end

    @testcase "caught exception" begin
        @test 1 == 1
        @test_throws ErrorException error("Caught exception")
        @test 1 == 1
    end

    @testcase "skip test" begin
        @test 1 == 1
        @test_skip 1 == 2
        @test 1 == 1
    end

    @testcase "expected failure" begin
        @test 1 == 1
        @test_broken 1 == 2
        @test 1 == 1
    end

    @testcase "unexpected pass" begin
        @test 1 == 1
        @test_broken 1 == 1
        @test 1 == 1
    end

    @testcase "with fixtures" for x in [1], y in [2]
    end

    @testcase "test_fail" begin
        @test_fail "Failure explanation"
    end
end

TESTCASES_WITH_GROUPS = Jute.collect_testobjs() do
    @testcase "root testcase 1" begin end
    @testgroup "group 1" begin
        @testcase "group 1 testcase 1" begin end
        @testcase "group 1 testcase 2" begin end
        @testgroup "subgroup 1" begin
            @testgroup "subsubgroup 1" begin
                @testcase "subsubgroup 1 testcase 1" begin end
                @testcase "subsubgroup 1 testcase 2" begin end
            end
        end
        @testcase "group 1 testcase 3" begin end
    end
    @testgroup "group 2" begin
        @testcase "group 2 testcase 1" begin end
        @testcase "group 2 testcase 2" begin end
    end
    @testcase "root testcase 2" begin end
end

failure_template = "Test Failed at <<<path>>>"
error_template = "Error During Test at <<<path>>>"
exception_template = "Got exception outside of a @test"

@testcase "verbosity0" begin
    exitcode, output = nested_run_with_output(TESTCASES, Dict(:verbosity => 0))
    @test exitcode == 1

    template = """
        ================================================================================
        multiple tests and one failure
        $failure_template
          Expression: 2 == 1
           Evaluated: 2 == 1
        ================================================================================
        uncaught exception
        $error_template
          $exception_template
          Uncaught exception
          Stacktrace:
        <<<MULTILINE>>>
        ================================================================================
        unexpected pass
        $error_template
         Unexpected Pass
         Expression: 1 == 1
         Got correct result, please change to @test if no longer broken.

        ================================================================================
        test_fail
        Test Failed
          Failure explanation
    """

    test_match_text(template, output)
end

@testcase "verbosity1" begin
    exitcode, output = nested_run_with_output(TESTCASES, Dict(:verbosity => 1))
    @test exitcode == 1

    template = """
        Collecting testcases...
        Using 10 out of 10 testcase definitions...
        ================================================================================
        Platform: <<<platform>>>, Julia <<<julia_version>>>, Jute <<<jute_version>>>
        --------------------------------------------------------------------------------
        ....*.F..E....B..B..E..F
        --------------------------------------------------------------------------------
        returning value (<<<time>>>) [10]
        --------------------------------------------------------------------------------
        20 tests passed, 2 failed, 2 errored in <<<full_time>>> (total test time <<<test_time>>>)
        ================================================================================
        multiple tests and one failure
        $failure_template
          Expression: 2 == 1
           Evaluated: 2 == 1
        ================================================================================
        uncaught exception
        $error_template
          $exception_template
          Uncaught exception
          Stacktrace:
        <<<MULTILINE>>>
        ================================================================================
        unexpected pass
        $error_template
         Unexpected Pass
         Expression: 1 == 1
         Got correct result, please change to @test if no longer broken.

        ================================================================================
        test_fail
        Test Failed
          Failure explanation
    """

    test_match_text(template, output)
end

@testcase "verbosity2" begin
    exitcode, output = nested_run_with_output(TESTCASES, Dict(:verbosity => 2))
    @test exitcode == 1

    template = """
        Collecting testcases...
        Using 10 out of 10 testcase definitions...
        ================================================================================
        Platform: <<<platform>>>, Julia <<<julia_version>>>, Jute <<<jute_version>>>
        --------------------------------------------------------------------------------
        multiple tests (<<<time>>>) [PASS] [PASS] [PASS]
        returning value (<<<time>>>) [PASS] [10]
        multiple tests and one failure (<<<time>>>) [PASS] [FAIL] [PASS]
        uncaught exception (<<<time>>>) [PASS] [ERROR]
        caught exception (<<<time>>>) [PASS] [PASS] [PASS]
        skip test (<<<time>>>) [PASS] [BROKEN] [PASS]
        expected failure (<<<time>>>) [PASS] [BROKEN] [PASS]
        unexpected pass (<<<time>>>) [PASS] [ERROR] [PASS]
        with fixtures [1,2] (<<<time>>>) [PASS]
        test_fail (<<<time>>>) [FAIL]
        --------------------------------------------------------------------------------
        20 tests passed, 2 failed, 2 errored in <<<full_time>>> (total test time <<<test_time>>>)
        ================================================================================
        multiple tests and one failure
        $failure_template
          Expression: 2 == 1
           Evaluated: 2 == 1
        ================================================================================
        uncaught exception
        $error_template
          $exception_template
          Uncaught exception
          Stacktrace:
        <<<MULTILINE>>>
        ================================================================================
        unexpected pass
        $error_template
         Unexpected Pass
         Expression: 1 == 1
         Got correct result, please change to @test if no longer broken.

        ================================================================================
        test_fail
        Test Failed
          Failure explanation
    """

    test_match_text(template, output)
end

@testcase "verbosity1 with groups" begin
    exitcode, output = nested_run_with_output(TESTCASES_WITH_GROUPS, Dict(:verbosity => 1))
    @test exitcode == 0

    template = """
        Collecting testcases...
        Using 9 out of 9 testcase definitions...
        ================================================================================
        Platform: <<<platform>>>, Julia <<<julia_version>>>, Jute <<<jute_version>>>
        --------------------------------------------------------------------------------
        .
        group 1: ..
          subgroup 1:
            subsubgroup 1: ..
        group 1 (cont.): .
        group 2: ..
        .
        --------------------------------------------------------------------------------
        9 tests passed, 0 failed, 0 errored in <<<full_time>>> (total test time <<<test_time>>>)
    """

    test_match_text(template, output)
end

@testcase "verbosity2 with groups" begin
    exitcode, output = nested_run_with_output(TESTCASES_WITH_GROUPS, Dict(:verbosity => 2))
    @test exitcode == 0

    template = """
        Collecting testcases...
        Using 9 out of 9 testcase definitions...
        ================================================================================
        Platform: <<<platform>>>, Julia <<<julia_version>>>, Jute <<<jute_version>>>
        --------------------------------------------------------------------------------
        root testcase 1 (<<<time>>>) [PASS]
        group 1/
          group 1 testcase 1 (<<<time>>>) [PASS]
          group 1 testcase 2 (<<<time>>>) [PASS]
          subgroup 1/
            subsubgroup 1/
              subsubgroup 1 testcase 1 (<<<time>>>) [PASS]
              subsubgroup 1 testcase 2 (<<<time>>>) [PASS]
          group 1 testcase 3 (<<<time>>>) [PASS]
        group 2/
          group 2 testcase 1 (<<<time>>>) [PASS]
          group 2 testcase 2 (<<<time>>>) [PASS]
        root testcase 2 (<<<time>>>) [PASS]
        --------------------------------------------------------------------------------
        9 tests passed, 0 failed, 0 errored in <<<full_time>>> (total test time <<<test_time>>>)
    """

    test_match_text(template, output)
end

@testcase "captured output" begin
    testcases = Jute.collect_testobjs() do
        @testcase "passing testcase" begin
            println(stdout, "stdout from passing testcase")
            println(stderr, "stderr from passing testcase")
        end

        @testcase "failing testcase" begin
            println(stdout, "stdout from failing testcase")
            @test 1 == 2
            println(stderr, "stderr from failing testcase")
        end
    end

    exitcode, output = nested_run_with_output(
        testcases, Dict(:verbosity => 0, :capture_output => true))

    template = """
        ================================================================================
        failing testcase
        Captured output:
        stdout from failing testcase
        stderr from failing testcase

        $failure_template
          Expression: 1 == 2
           Evaluated: 1 == 2
    """

    test_match_text(template, output)
end

end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    UnionProcess(p₁, p₂)

Union (or superposition) of spatial point processes `p₁` and `p₂`.
"""
struct UnionProcess{P₁<:PointProcess,P₂<:PointProcess} <: PointProcess
  p₁::P₁
  p₂::P₂
end

ishomogeneous(p::UnionProcess) = ishomogeneous(p.p₁) && ishomogeneous(p.p₂)

struct UnionSampling end

default_sampling_algorithm(::UnionProcess) = UnionSampling()

function rand_single(p::UnionProcess, g::Geometry, algo::UnionSampling)
  pp₁ = rand(p.p₁, g)
  pp₂ = rand(p.p₂, g)

  X = [coordinates(pp₁[i]) for i in 1:nelements(pp₁)]
  Y = [coordinates(pp₂[i]) for i in 1:nelements(pp₂)]

  PointSet([X; Y])
end


text:
# Activate environment

using FileIO,JSON

# Grab the two arguments that are passed in
# The first is the directory containing the output of the mapper
# The second is the name of the file where the final word counts should be saved
inputdir = ARGS[1]
outputfile = ARGS[2]

# Read in the first word count file
fnames = readdir(inputdir)
allcounts = JSON.parsefile.(inputdir*"/".*fnames)

globalcounts = merge(+,allcounts...)

# Sort and print the top 5 words with their counts
globalranks = sort(collect(globalcounts), by=x->x[2], rev=true)
println(globalranks[1:5,:])

# Save to the output file
output = join([k*","*string(v) for (k,v) in globalranks[1:5,:]],'\n')*"\n"
write(outputfile*".csv",output)

max_stars_repo_name: llsc-supercloud/teaching-examples

max_issues_repo_name: llsc-supercloud/teaching-examples

max_forks_repo_name: llsc-supercloud/teaching-examples


text:
using Random: randperm

# Sum of square error between two arrays
function SSE(x::AbstractArray{T}, y::AbstractArray{T})::T where {T<:Real}
    diff = (x - y)
    return sum(diff .* diff)
end

# Sum of square error between two arrays, with weights
function SSE(x::AbstractArray{T}, y::AbstractArray{T}, w::AbstractArray{T})::T where {T<:Real}
    diff = (x - y)
    return sum(diff .* diff .* w)
end

# Mean of square error between two arrays
function MSE(x::AbstractArray{T}, y::AbstractArray{T})::T where {T<:Real}
    return SSE(x, y)/size(x)[1]
end

# Mean of square error between two arrays
function MSE(x::AbstractArray{T}, y::AbstractArray{T}, w::AbstractArray{T})::T where {T<:Real}
    return SSE(x, y, w)/sum(w)
end

# Loss function. Only MSE implemented right now. TODO
# Also need to put actual loss function in scoreFuncBatch!
function EvalLoss(tree::Node, dataset::Dataset{T}, options::Options)::T where {T<:Real}
    (prediction, completion) = evalTreeArray(tree, dataset.X, options)
    if !completion
        return convert(T, 1000000000)
    end

    if dataset.weighted
        return MSE(prediction, dataset.y, dataset.weights)
    else
        return MSE(prediction, dataset.y)
    end
end

# Score an equation
function scoreFunc(dataset::Dataset{T},
                   baseline::T, tree::Node,
                   options::Options)::T where {T<:Real}
    mse = EvalLoss(tree, dataset, options)
    return mse / baseline + countNodes(tree)*options.parsimony
end

# Score an equation with a small batch
function scoreFuncBatch(dataset::Dataset{T}, baseline::T,
                        tree::Node, options::Options)::T where {T<:Real}
    batchSize = options.batchSize
    batch_idx = randperm(dataset.n)[1:options.batchSize]
    batch_X = dataset.X[batch_idx, :]
    batch_y = dataset.y[batch_idx]
    (prediction, completion) = evalTreeArray(tree, batch_X, options)
    if !completion
        return convert(T, 1000000000)
    end

    if dataset.weighted
        mse = MSE(prediction, batch_y)
    else
        batch_w = dataset.weights[batch_idx]
        mse = MSE(prediction, batch_y, batch_w)
    end
    return mse / baseline + countNodes(tree) * options.parsimony
end


text:
using PowerSystems

##Nodes
nodes_DAIB= [Bus(1 , #number
                 "Bus 1", #Name
                 "REF" , #BusType (REF, PV, PQ)
                 0, #Angle in radians
                 1.04, #Voltage in pu
                 (min=0.94, max=1.06), #Voltage limits in pu
                 0.69),  #Base voltage in kV
                 Bus(2 , "Bus 2"  , "PV" ,  0 , 1.0 , (min=0.94, max=1.06), 0.69)]

branch_DAIB = [Line("Line1", #name
                    true, #available
                    0.0, #active power flow initial condition (from-to)
                    0.0, #reactive power flow initial condition (from-to)
                    Arc(from=nodes_DAIB[1],to=nodes_DAIB[2]), #Connection between buses
                    0.0, #resistance in pu
                    0.075, #reactance in pu
                    (from=0.0, to=0.0), #susceptance in pu
                    5.0, #rate in MW
                    1.04)]  #angle limits (-min and max)

##Inf bus
inf_gen_DAIB = StaticSource(1, #number
                    :InfBus, #name
                    nodes_DAIB[1],#bus
                    1.00, #VR
                    0.0, #VI
                    0.000005) #Xth

## Inverter

converter = AvgCnvFixedDC(690.0, #Rated Voltage
                          2.75) #Rated MVA
dc_source = FixedDCSource(600.0) #Not in the original data, guessed.

filter = LCLFilter(0.08, #Series inductance lf in pu
                   0.003, #Series resitance rf in pu
                   0.074, #Shunt capacitance cf in pu
                   0.2, #Series ractance rg to grid connection (#Step up transformer or similar)
                   0.01) #Series resistance lg to grid connection (#Step up transformer or similar)
pll = PLL(500.0, #ω_lp: Cut-off frequency for LowPass filter of PLL filter.
          0.084, #k_p: PLL proportional gain
          4.69) #k_i: PLL integral gain

virtual_H = VirtualInertia(2.0, #Ta:: VSM inertia constant
                           400.0, #kd:: VSM damping coefficient
                           20.0, #kω:: Frequency droop gain in pu
                           2*pi*50.0) #ωb:: Rated angular frequency

Q_control = ReactivePowerDroop(0.2, #kq:: Reactive power droop gain in pu
                              1000.0) #ωf:: Reactive power cut-off low pass filter frequency

outer_control = VirtualInertiaQdroop(virtual_H, Q_control)

vsc = CombinedVIwithVZ(0.59, #kpv:: Voltage controller proportional gain
                       736.0, #kiv:: Voltage controller integral gain
                       0.0, #kffv:: Binary variable enabling the voltage feed-forward in output of current controllers
                       0.0, #rv:: Virtual resistance in pu
                       0.2, #lv: Virtual inductance in pu
                       1.27, #kpc:: Current controller proportional gain
                       14.3, #kiv:: Current controller integral gain
                       0.0, #kffi:: Binary variable enabling the current feed-forward in output of current controllers
                       50.0, #ωad:: Active damping low pass filter cut-off frequency
                       0.2) #kad:: Active damping gain

Darco_Inverter = DynInverter(1, #number
                             :DARCO, #name
                             nodes_DAIB[2], #bus location
                             1.0, #ω_ref
                             1.02, #V_ref
                             0.5, #P_ref
                             0.0, #Q_ref
                             2.75, #MVABase
                             converter, #Converter
                             outer_control, #OuterControl
                             vsc, #Voltage Source Controller
                             dc_source, #DC Source
                             pll, #Frequency Estimator
                             filter) #Output Filter

loads_DAIB = [PowerLoad("Bus1", true, nodes_DAIB[1], PowerSystems.ConstantPower, 0.0, 0.0, 0.0, 0.0)]

DAIB = DynamicSystem(nodes_DAIB, branch_DAIB, [Darco_Inverter], [inf_gen_DAIB], 100.0, 50.0);

#DAIB = DynamicSystem(nodes_DAIB, branch_DAIB, [Darco_Inverter], loads_DAIB, 100.0, 50.0)


text:
#=
cbf:
- Julia version: 1.5.1
- Author: Berk
- Date: 2020-09-02
=#

# Test CBF imports
#     gams01, netmod_dol2, netmod_kar1, netmod_kar2, du-opt, du-opt5, nvs03, ex1223, ex1223a, ex1223b, gbd
filename = string("data/cblib/shortfall_20_15.cbf.gz");
# Potential examples from MINLPLib
# gams01, netmod_dol2, netmod_kar1, netmod_kar2, du-opt, du-opt5, nvs03, ex1223, ex1223a, ex1223b, gbd
model = JuMP.read_from_file(filename);
set_optimizer(model, Gurobi.Optimizer)
optimize!(model);
moi_obj = JuMP.getobjectivevalue(model);
@test moi_obj ≈ -1.0792654303

# Testing CBF import to ModelData
md = OCTHaGOn.CBF_to_ModelData(filename);
md.name = "shortfall_20_15"
find_bounds!(md, all_bounds = true);
update_bounds!(md, lbs = Dict(md.vks .=> 0.));
update_bounds!(md, ubs = Dict(md.vks .=> 1.));

# Testing constraint import.
n_samples = 200;
bbl = md.bbls[1];
X = lh_sample(bbl, n_samples=n_samples);
Y = md.bbls[1](X);

# Testing sample_and_eval for combined LH and boundary sampling.
uniform_sample_and_eval!(md);
@test size(bbl.X) == (500,length(bbl.vks))
@test all(feasibility(md) .>= 0)

learn_constraint!(md);

# Solving the model
status = optimize!(md)
OCT_vars = JuMP.getvalue.(md.vars);


text:
# Day 19: Beacon Scanner
# Not super proud of this one, but it gets the job done...
const BeaconOffset = Tuple{Int,Int,Int}
const ScannerReading = Vector{BeaconOffset}

function load(::Day{19}, path)
    scanner_readings = Vector{ScannerReading}()
    for line in readlines(path)
        if isempty(line) && continue
        elseif line[2] == '-'
            push!(scanner_readings, ScannerReading()) # add scanner
        else
            push!(last(scanner_readings), Tuple(parse.(Int, split(line, ","))))
        end
    end
    return scanner_readings
end

# Part 1 – count number of unique beacons
function distances!(r::ScannerReading)
    n = length(r)
    sort!(r) # tuples get sorted in lexicographic order
    return [r[j] .- r[i] for i in 1:(n - 1) for j in (i + 1):n]
end

# This brilliant idea was shamelessly stolen from Gunnar Farnebäck in the Julia Zulip:
# We use an second distance-metric that is invariant w.r.t. rotation and sign changes.
function invariant_distances!(r::ScannerReading)
    n = length(r)
    sort!(r)
    return [sum(abs.(r[j] .- r[i])) for i in 1:n for j in (i + 1):n]
end

function rotate_scan!(new_reading::ScannerReading, old_reading::ScannerReading)
    old_dist = distances!(old_reading)
    for fwdir in [1, 2, 3]
        for fwsign in [1, -1]
            for updir in setdiff([1, 2, 3], fwdir)
                for upsign in [1, -1]
                    # The last direction is now fully determined:
                    rightdir = first(setdiff([1, 2, 3], [fwdir, updir]))
                    if [fwdir, updir] in [[1, 2], [2, 3], [3, 1]] # even permutation
                        rightsign = fwsign * upsign
                    else
                        rightsign = -fwsign * upsign
                    end

                    # Apply the rotation
                    directions = [fwdir, updir, rightdir]
                    signs = [fwsign, upsign, rightsign]
                    rotated = [Tuple(beacon[directions] .* signs) for beacon in new_reading]

                    # and stop the routine when a match is found
                    matches = length(intersect(old_dist, distances!(rotated)))
                    if matches >= 12
                        new_reading .= rotated
                        return nothing
                    end
                end
            end
        end
    end
    return nothing
end

function align_scan!(new_reading::ScannerReading, old_reading::ScannerReading)
    diffs = [a .- b for a in old_reading for b in new_reading]
    for d in diffs
        if count(==(d), diffs) >= 12
            map!(r -> r .+ d, new_reading, new_reading)
            return d
        end
    end
    return (0, 0, 0)
end

function solve1(::Day{19}, data)
    # precompute invariant distances for each scanner:
    n = length(data)
    dists = invariant_distances!.(data)
    checked = zeros(Bool, n)

    # Use positions of beacon 1 as reference
    beacon_positions = Set(data[1])
    checked[1] = true

    while !all(checked)
        # Select scanners to match
        matches = Dict{Tuple{Int,Int},Int}()
        for candidate in findall((!).(checked))
            for c in findall(checked)
                matches[(candidate, c)] = length(intersect(dists[candidate], dists[c]))
            end
        end
        new_scanner, old_scanner = argmax(matches)

        # Shift reading in place
        rotate_scan!(data[new_scanner], data[old_scanner])
        align_scan!(data[new_scanner], data[old_scanner])

        # Add beacon positions
        union!(beacon_positions, data[new_scanner])
        checked[new_scanner] = true
    end
    return length(beacon_positions)
end

# Part 2 – largest Manhattan distance between any two scanners
function solve2(::Day{19}, data)
    # precompute invariant distances for each scanner:
    n = length(data)
    dists = invariant_distances!.(data)
    checked = zeros(Bool, n)

    # Use positions of beacon 1 as reference
    checked[1] = true
    positions = [(0, 0, 0)]

    while !all(checked)
        # Select scanners to match
        matches = Dict{Tuple{Int,Int},Int}()
        for candidate in findall((!).(checked))
            for c in findall(checked)
                matches[(candidate, c)] = length(intersect(dists[candidate], dists[c]))
            end
        end
        new_scanner, old_scanner = argmax(matches)

        # Shift reading in place
        rotate_scan!(data[new_scanner], data[old_scanner])
        push!(positions, align_scan!(data[new_scanner], data[old_scanner]))
        checked[new_scanner] = true
    end
    return maximum(sum(abs, positions[i] .- positions[j]) for i in 1:n for j in (i + 1):n)
end

testresult1(::Day{19}) = 79
testresult2(::Day{19}) = 3621


text:
function paint(r, counter, k)
    if (counter%2)==1
        while !isborder(r, West)
            move!(r, West)
            putmarker!(r)
        end
    else
        putmarker!(r)
        for _ = 1:((k-1)-counter)
            move!(r, Ost)
            putmarker!(r)
        end
    end
end

function goToHome(r, bot, right)
    while !isborder(r, West)
        move!(r, West)
    end
    for _ = 1:bot
        move!(r, Nord)
    end
    for _ = 1:right
        move!(r, Ost)
    end
end

function drawPir(r)
    right = 0

    while !isborder(r, West)
        right += 1
        move!(r, West)
    end

    k = 1
    counter = 1
    bot = 0
    if !isborder(r, Sud)
        move!(r, Sud)
        bot+=1
        putmarker!(r)
        while !isborder(r, Ost)
            k+=1
            move!(r, Ost)
            putmarker!(r)
        end
        move!(r, Sud)
        bot+=1
    end

    while !isborder(r, Sud)
        if counter==k
            break
        end

        paint(r, counter, k)

        counter+=1
        move!(r, Sud)
        bot+=1
    end

    paint(r, counter, k)

    goToHome(r, bot, right)

end

max_stars_repo_name: Trilonka/MireaAl

max_issues_repo_name: Trilonka/MireaAl

max_forks_repo_name: Trilonka/MireaAl


text:
"""
    eigROT_Distance(P, Q; edge_length = 1, α = 1.0)

computes the ROT distance matrix of P's column vectors on a graph.

# Input Argument
- `P::Matrix{Float64}`: a matrix whose columns are vector measures with the same total mass.
- `Q::SparseMatrixCSC{Int64,Int64}`: the unweighted incidence matrix of the graph.
- `edge_length::Any`: the length vector (default: 1 represents unweighted graphs)
- `α::Float64`: default is 1.0. ROT parameter.

# Output Argument
- `dis::Matrix{Float64}`: distance matrix, dis[i,j] = d_{ROT}(pᵢ, pⱼ; α).

"""
function eigROT_Distance(P::Matrix{Float64}, Q::SparseMatrixCSC{Int64,Int64};
                            edge_length::Any = 1, α::Float64 = 1.0)
    n = size(P, 2)
    total_mass = sum(P, dims = 1)[:]
    if norm(total_mass - total_mass[1] * ones(n), Inf) > 10^4 * eps()
        @error("P's column measures do not share the same total mass.")
        return
    end
    dis = zeros(n,n)
    le2 = [edge_length; edge_length]
    Q2 = [Q -Q]
    m2 = size(Q2, 2)
    for i = 1:(n - 1), j = (i + 1):n
        f = P[:, i] - P[:, j]
        md = Model(optimizer_with_attributes(Clp.Optimizer, "LogLevel" => 0))
        @variable(md, w[1:m2] >= 0.0)
        edge_length == 1 ? @objective(md, Min, sum(w)) : @objective(md, Min, sum(w .* le2))
        @constraint(md, Q2 * w .== f)
        JuMP.optimize!(md)
        wt = abs.(JuMP.value.(w))
        dis[i,j] = edge_length == 1 ? norm(wt.^α, 1) : norm((wt.^α) .* le2, 1)
    end
    return dis + dis'
end

"""
    ROT_Distance(A, B, Q; edge_length = 1, α = 1.0)

computes the ROT distance matrix from A's column vectors to B's column vectors.
If A, B are vector inputs, then it also returns the cost value and the optimal
transport plan.

# Input Argument
- `A::Any`: a vector or matrix whose columns are initial probability measures.
- `B::Any`: a vector or matrix whose columns are terminal probability measures.
- `Q::SparseMatrixCSC{Int64,Int64}`: the unweighted incidence matrix of the graph.
- `edge_length::Any`: the length vector (default: `1` represents unweighted graphs)
- `α::Float64`: default is `1.0`. ROT parameter.
- `retPlan::Bool`: an indicator if return the optimal plan (default: `false`)

# Output Argument
- `dis::Matrix{Float64}`: distance matrix, dis[i,j] = d_{ROT}(aᵢ, bⱼ; α).

"""
function ROT_Distance(A::Any, B::Any, Q::SparseMatrixCSC{Int64,Int64};
                        edge_length::Any = 1, α::Float64 = 1.0, retPlan::Bool = false)
    m = ndims(A) > 1 ? size(A, 2) : 1
    n = ndims(B) > 1 ? size(B, 2) : 1
    dis = zeros(m, n)
    le2 = [edge_length; edge_length]
    Q2 = [Q -Q]
    m2 = size(Q2, 2)
    for i = 1:m, j = 1:n
        f = (ndims(A) > 1 && ndims(B) > 1) ? B[:,j] - A[:,i] : B - A
        md = Model(optimizer_with_attributes(Clp.Optimizer, "LogLevel" => 0))
        @variable(md, w[1:m2] >= 0.0)
        edge_length == 1 ? @objective(md, Min, sum(w)) : @objective(md, Min, sum(w .* le2))
        @constraint(md, Q2 * w .== f)
        JuMP.optimize!(md)
        wt = abs.(JuMP.value.(w))
        dis[i, j] = edge_length == 1 ? norm(wt.^α, 1) : norm((wt.^α) .* le2, 1)
        if ndims(A) == 1 && ndims(B) == 1
            if retPlan
                return dis[1, 1], wt
            else
                return dis[1, 1]
            end
        end
    end
    return dis
end

using OptimalTransport
"""
    eigEMD_Distance(P, C)

EIGEMD\\_DISTANCE computes the EMD distance matrix of P's column vectors on a graph.

# Input Argument
- `P::Matrix{Float64}`: a matrix whose columns are probability measures.
- `C::Matrix{Float64}`: the cost matrix, which is the ground metric between pairwise nodes.

# Output Argument
- `dis::Matrix{Float64}`: distance matrix, dis[i,j] = d\\_EMD(pᵢ, pⱼ).

"""
function eigEMD_Distance(P, C)
    n = size(P,2)
    dis = zeros(n,n)
    for i = 1:n-1, j = i+1:n
        dis[i,j] = emd2(P[:,i], P[:,j], C)
    end
    return dis + dis'
end


text:
using ChaosTools
using DynamicalSystemsBase
using Test
using LinearAlgebra
using OrdinaryDiffEq

println("\nTesting basins of attraction...")

@testset "Basins tests" begin

@testset "Discrete map" begin
    ds = Systems.henon(zeros(2); a = 1.4, b = 0.3)
    xg = yg = range(-2.,2.,length=100)
    basin, attractors = basins_of_attraction((xg,yg), ds; show_progress = false)
    # pcolormesh(xg,yg, basin')
    @test 4260 ≤ count(basin .== 1) ≤ 4280
    @test 5700 ≤ count(basin .== -1) ≤ 5800
end

@testset "Test basin stroboscopic map" begin
    ds = Systems.duffing([0.1, 0.25]; ω = 1.0, f = 0.2, d = 0.15, β = -1)
    xg = yg = range(-2.2,2.2,length=100)
    T = 2π/1.0
    basin,attractors = basins_of_attraction((xg,yg), ds;
    T, diffeq = (;alg=Tsit5()), show_progress = false)
    # pcolormesh(xg, yg, basin')
    @test length(unique(basin)) == 2
    @test 4900 ≤ count(basin .== 1) ≤ 5100
    @test  4900 ≤  count(basin .== 2) ≤ 5100
end

@testset "Test basin poincare map" begin
    ds = Systems.thomas_cyclical(b = 0.1665)
    xg = yg = range(-6.0, 6.0; length = 100)
    pmap = poincaremap(ds, (3, 0.0), 1e6;
        rootkw = (xrtol = 1e-8, atol = 1e-8), diffeq=(reltol=1e-9,)
    )
    basin,attractors = basins_of_attraction((xg,yg), pmap; show_progress = false)
    # pcolormesh(xg,yg, basin')
    @test length(attractors) == 3
    @test  4610 ≤ count(basin .== 1) ≤ 4641
    @test  2660 ≤ count(basin .== 2) ≤ 2691
    @test  2640 ≤ count(basin .== 3) ≤ 2691
end

@testset "basins_of_attraction" begin
    ds = Systems.magnetic_pendulum(γ=1, d=0.2, α=0.2, ω=0.8, N=3)
    xg = range(-2,2,length=100)
    yg = range(-2,2,length=100)
    complete_state(y) = SVector(0.0, 0.0)
    basin, attractors = basins_of_attraction((xg,yg), ds;
    idxs=1:2, Δt = 1., complete_state, show_progress = false)
    # pcolormesh(xg,yg, basin')
    @test count(basin .== 1) == 3332
    @test count(basin .== 2) == 3332

    # Now test the zoom capability
    xg = yg = range(-2,-1.9,length=50)
    basins, att = basins_of_attraction((xg,yg), ds;
    idxs = 1:2, Δt = 1., complete_state, show_progress = false,
    attractors = attractors, mx_chk_lost = 1000, ε = 1e-3)
    @test count(basins .== 2) == 407
    @test count(basins .== 3) == 737
end

@testset "3D basins" begin
    ds = Systems.lorenz84()
    xg=yg=range(-1.,2.,length=100)
    zg=range(-1.5,1.5,length=30)
    bsn,att = basins_of_attraction((xg, yg, zg), ds; Δt = 1.0, show_progress = false)
    @test length(size(bsn)) == 3
    for i in 1:size(bsn)[3]
        # While there are 4 attractors, because system is chaotic we might
        # miss one of the 4 due to coarse state space partition
        @test sort!(unique(bsn[:, :, i])) ∈ ([1,2,3], [1,2,3,4])
    end
end

@testset "Basins for in-place system" begin
    ds = Systems.henon_iip(zeros(2); a = 1.4, b = 0.3)
    xg = yg = range(-2.,2.;length=100)
    basin, attractors = basins_of_attraction((xg,yg), ds; show_progress = false)
    # pcolormesh(xg,yg, basin')
    @test 4260 ≤ count(basin .== 1) ≤ 4280
    @test 5600 ≤ count(basin .== -1) ≤ 5800

    ds = Systems.lorenz_iip()
    xg = yg = range(-30.0, 30.0; length=10)
    zg = range(0, 50; length=10)
    basin, attractors = basins_of_attraction((xg,yg,zg), ds; show_progress=false, mx_chk_fnd_att = 1000)
    @test length(attractors) == 1
    @test unique(basin) == [1]
end

@testset "matching attractors" begin
    d, α, ω = 0.3, 0.2, 0.5
    ds = Systems.magnetic_pendulum(; d, α, ω)
    xg = yg = range(-3, 3, length = 100)
    b₋, a₋ = basins_of_attraction((xg,yg), ds; Δt=1., idxs=1:2, show_progress = false)
    @testset "method $method" for method ∈ (:overlap, :distance)
        @testset "γ3 $γ3" for γ3 ∈ [0.2, 0.1] # still 3 at 0.2, but only 2 at 0.1
            ds = Systems.magnetic_pendulum(; d, α, ω,  γs = [1, 1, γ3])
            b₊, a₊ = basins_of_attraction((xg,yg), ds; Δt=1., idxs=1:2, show_progress = false)
            match_attractors!(b₋, a₋, b₊, a₊, method)
            for k in keys(a₊)
                dist = minimum(norm(x .- y) for x ∈ a₊[k] for y ∈ a₋[k])
                @test dist < 0.2
            end
        end
    end
end

end


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

using ACE, JuLIP, BenchmarkTools

#---

fname = "/Users/ortner/Dropbox/Public/SHIPPots/Si_split_1.2_reg_test_v05.json"
D = load_dict(fname)
V = decode_dict(D)

#---

at = bulk(:Si, cubic=true) * 3;
@show length(at)

println("Energy:")
@btime energy($V, $at)
println("Forces:")
@btime forces($V, $at)

#---
print("time/atom (energy): ");
println(" $((@timed energy(V, at)).time / length(at) * 1000) ms")

print("time/atom (forces): ");
println(" $((@timed forces(V, at)).time / length(at) * 1000) ms")


text:
;(import Pkg)
;((. Pkg add) (PackageSpec (kw path "https://github.com/BenLauwens/ThinkJulia.jl")))

(using ThinkJulia)

(function (square t)
  (for ((i (: 1 4)))
    (forward t 100)
    (turn t -90)))

(= t (Turtle))
(@svg (begin
  (square t)))


text:
@testset "testpattern" begin
    tp = testpattern()
    @test typeof(tp) <: Matrix
    @test eltype(tp) <: RGBA
    @test size(tp) === (300,400)
    tp2 = augment(NoOp())
    @test tp == tp2
end

# @testset "rand_mutex" begin
#     mutex = Augmentor.rand_mutex[]
#     typeof(mutex) <: Threads.Mutex
#     # check that its not a null pointer
#     @test reinterpret(Int, mutex.handle) > 0
# end

@testset "safe_rand" begin
    num = @inferred Augmentor.safe_rand()
    @test 0 <= num <= 1
    @test typeof(num) <: Float64
    num = @inferred Augmentor.safe_rand(2)
    @test all(0 .<= num .<= 1)
    @test typeof(num) <: Vector{Float64}
end

@testset "maybe_copy" begin
    A = [1 2 3; 4 5 6; 7 8 9]
    Ao = OffsetArray(A, (-2,-1))
    @test @inferred(Augmentor.maybe_copy(A)) === A
    @test @inferred(Augmentor.maybe_copy(Ao)) === Ao
    let Ast = @SMatrix [1 2 3; 4 5 6; 7 8 9]
        @test @inferred(Augmentor.maybe_copy(Ast)) === Ast
    end
    let v = view(A, 2:3, 1:2)
        @test @inferred(Augmentor.maybe_copy(v)) == A[2:3, 1:2]
        @test typeof(Augmentor.maybe_copy(v)) <: Array
    end
    let v = view(OffsetArray(A, (-2,-1)), 0:1, 0:1)
        @test @inferred(Augmentor.maybe_copy(v)) == A[2:3, 1:2]
        @test typeof(Augmentor.maybe_copy(v)) <: Array
    end
    let v = view(A, IdentityRange(2:3), IdentityRange(1:2))
        @test @inferred(Augmentor.maybe_copy(v)) == OffsetArray(A[2:3,1:2], (1, 0))
        @test typeof(Augmentor.maybe_copy(v)) <: OffsetArray
    end
    let v = channelview(rect)
        @test @inferred(Augmentor.maybe_copy(v)) == channelview(rect)
        @test typeof(Augmentor.maybe_copy(v)) <: Array
    end
    let p = permuteddimsview(A, (2,1))
        @test @inferred(Augmentor.maybe_copy(p)) == A'
        @test typeof(Augmentor.maybe_copy(p)) <: Array
    end

    let p = permuteddimsview(Ao, (2,1))
        pc = @inferred(Augmentor.maybe_copy(p))
        #@test pc == Ao'
        @test pc == OffsetArrays.no_offset_view(Ao')
        @test typeof(pc) <: Array
    end
    let p = view(permuteddimsview(A, (2,1)), IdentityRange(2:3), IdentityRange(1:2))
        @test @inferred(Augmentor.maybe_copy(p)) == OffsetArray(A'[2:3, 1:2],1,0)
        @test typeof(Augmentor.maybe_copy(p)) <: OffsetArray
    end
    let Aa = Augmentor.prepareaffine(A)
        @test @inferred(Augmentor.maybe_copy(Aa)) == OffsetArray(A, (0,0))
        @test typeof(Augmentor.maybe_copy(Aa)) <: OffsetArray
    end
    let Ar = reshape(view(A,:,:),1,3,3)
        @test @inferred(Augmentor.maybe_copy(Ar)) == reshape(A,1,3,3)
        @test typeof(Augmentor.maybe_copy(Ar)) <: Array{Int,3}
    end
    let Ar = reshape(view(Ao,:,:),1,3,3)
        @test @inferred(Augmentor.maybe_copy(Ar)) == reshape(A,1,3,3)
        @test typeof(Augmentor.maybe_copy(Ar)) <: Array{Int,3}
    end
end

@testset "plain_array" begin
    A = [1 2 3; 4 5 6; 7 8 9]
    @test @inferred(Augmentor.plain_array(A)) === A
    @test @inferred(Augmentor.plain_array(OffsetArray(A, (-2,-1)))) === A
    let As = sparse(A)
       @test @inferred(Augmentor.plain_array(As)) == A
       @test typeof(Augmentor.plain_array(As)) <: Array
       Ar = reshape(As, 3, 3, 1)
       @test typeof(Ar) <: Base.ReshapedArray
       @test @inferred(Augmentor.plain_array(Ar)) == reshape(A,3,3,1)
       @test typeof(Augmentor.plain_array(Ar)) <: Array
    end
    let Ast = @SMatrix [1 2 3; 4 5 6; 7 8 9]
        @test @inferred(Augmentor.plain_array(Ast)) === Ast
    end
    let v = view(A, 2:3, 1:2)
        @test @inferred(Augmentor.plain_array(v)) == A[2:3, 1:2]
        @test typeof(Augmentor.plain_array(v)) <: Array
    end
    let v = view(OffsetArray(A, (-2,-1)), 0:1, 0:1)
        @test @inferred(Augmentor.plain_array(v)) == A[2:3, 1:2]
        @test typeof(Augmentor.plain_array(v)) <: Array
    end
    let v = view(A, IdentityRange(2:3), IdentityRange(1:2))
        @test @inferred(Augmentor.plain_array(v)) == A[2:3, 1:2]
        @test typeof(Augmentor.plain_array(v)) <: Array
    end
    let v = channelview(rect)
        @test @inferred(Augmentor.plain_array(v)) == channelview(rect)
        @test typeof(Augmentor.plain_array(v)) <: Array
    end
    let p = permuteddimsview(A, (2,1))
        @test @inferred(Augmentor.plain_array(p)) == A'
        @test typeof(Augmentor.plain_array(p)) <: Array
    end
    let p = view(permuteddimsview(A, (2,1)), IdentityRange(2:3), IdentityRange(1:2))
        @test @inferred(Augmentor.plain_array(p)) == A'[2:3, 1:2]
        @test typeof(Augmentor.plain_array(p)) <: Array
    end
    let Aa = Augmentor.prepareaffine(A)
        @test @inferred(Augmentor.plain_array(Aa)) == A
        @test typeof(Augmentor.plain_array(Aa)) <: Array
    end
    let Ar = reshape(view(A,:,:),1,3,3)
        @test @inferred(Augmentor.plain_array(Ar)) == reshape(A,1,3,3)
        @test typeof(Augmentor.plain_array(Ar)) <: Array{Int,3}
    end
end

@testset "plain_axes" begin
    A = [1 2 3; 4 5 6; 7 8 9]
    @test @inferred(Augmentor.plain_axes(A)) === A
    @test @inferred(Augmentor.plain_axes(OffsetArray(A, (-2,-1)))) === A
    let v = view(A, 2:3, 1:2)
        @test @inferred(Augmentor.plain_axes(v)) === v
    end
    let v = view(OffsetArray(A, (-2,-1)), 0:1, 0:1)
        @test @inferred(Augmentor.plain_axes(v)) === v
    end
    let v = view(A, IdentityRange(2:3), IdentityRange(1:2))
        @test @inferred(Augmentor.plain_axes(v)) === view(A, 2:3, 1:2)
    end
    let v = channelview(rect)
        @test @inferred(Augmentor.plain_axes(v)) === v
    end
    let p = permuteddimsview(A, (2,1))
        @test @inferred(Augmentor.plain_axes(p)) === p
    end
    let p = view(permuteddimsview(A, (2,1)), IdentityRange(2:3), IdentityRange(1:2))
        @test @inferred(Augmentor.plain_axes(p)) === view(parent(p), 2:3, 1:2)
    end
    let Aa = Augmentor.prepareaffine(A)
        @test @inferred(Augmentor.plain_axes(Aa)) === view(Aa, axes(Aa)...)
    end
    let Ar = reshape(view(A,:,:),1,3,3)
        @test @inferred(Augmentor.plain_axes(Ar)) === Ar
    end
end

@testset "match_idx" begin
    A = [1 2 3; 4 5 6; 7 8 9]
    @test @inferred(Augmentor.match_idx(A, axes(A))) === A
    let img = @inferred Augmentor.match_idx(A, (2:4, 2:4))
        @test axes(img) === Base.IdentityUnitRange.((2:4, 2:4))
        @test typeof(img) <: OffsetArray
    end
    let B = view(A,1:3,1:3)
        @test @inferred(Augmentor.match_idx(B, axes(B))) === B
    end
    let B = view(A,1:3,1:3)
        img = @inferred(Augmentor.match_idx(B, B.indices))
        @test axes(img) === Base.IdentityUnitRange.((1:3, 1:3))
        @test typeof(img) <: OffsetArray
    end
    let img = @inferred Augmentor.match_idx(view(A,1:3,1:3), (2:4,2:4))
        @test axes(img) === Base.IdentityUnitRange.((2:4, 2:4))
        @test typeof(img) <: OffsetArray
    end
    let C = Augmentor.prepareaffine(A)
        @test @inferred(Augmentor.match_idx(C, (2:4, 2:4))) === C
    end
end

@testset "direct_view" begin
    A = [1 2 3; 4 5 6; 7 8 9]
    Av = view(A, IdentityRange(2:3), IdentityRange(1:2))
    As = view(A, 3:-1:2, 1:1:2)
    @test_throws MethodError Augmentor.indirect_axes((), ())
    @test_throws MethodError Augmentor.direct_axes((), ())
    @test_throws MethodError Augmentor.direct_view(A, ())
    @test @inferred(Augmentor.direct_view(A, (2:3,1:2))) === Av
    @test @inferred(Augmentor.indirect_view(A, (2:3,1:2))) === Av
    @test @inferred(Augmentor.direct_view(A, (IdentityRange(2:3),IdentityRange(1:2)))) === Av
    @test @inferred(Augmentor.indirect_view(A, (IdentityRange(2:3),IdentityRange(1:2)))) === Av
    @test @inferred(Augmentor.direct_view(A, (3:-1:2,1:1:2))) === As
    @test @inferred(Augmentor.indirect_view(A, (3:-1:2,1:1:2))) === As
    @test @inferred(Augmentor.direct_view(Av, (3:3,1:2))) === view(Av,3:3,1:2)
    @test @inferred(Augmentor.indirect_view(Av, (2:2,1:2))) === view(Av,3:3,1:2)
    @test @inferred(Augmentor.direct_view(Av, (IdentityRange(3:3),IdentityRange(1:2)))) === view(Av,3:3,1:2)
    @test @inferred(Augmentor.indirect_view(Av, (IdentityRange(2:2),IdentityRange(1:2)))) === view(Av,3:3,1:2)
    @test_throws MethodError Augmentor.direct_view(Av, (3:-1:2,2:-1:1))
    @test @inferred(Augmentor.indirect_view(Av, (2:-1:1,2:-1:1))) === view(A,3:-1:2,2:-1:1)
end

@testset "vectorize" begin
    @test @inferred(Augmentor.vectorize(2)) === 2:2
    @test @inferred(Augmentor.vectorize(2.3)) === 2.3:2.3
    @test @inferred(Augmentor.vectorize(3:4)) === 3:4
    @test @inferred(Augmentor.vectorize(3:1:4)) === 3:1:4
    @test @inferred(Augmentor.vectorize(3.0:4)) === 3.0:4
    @test @inferred(Augmentor.vectorize(3.0:1:4)) === 3.0:1:4
    @test @inferred(Augmentor.vectorize(Base.OneTo(4))) === Base.OneTo(4)
end

@testset "round_if_float" begin
    @test @inferred(Augmentor.round_if_float(3,2)) === 3
    @test @inferred(Augmentor.round_if_float(3.1111,2)) === 3.11
    @test @inferred(Augmentor.round_if_float((3,3.1111),2)) === (3,3.11)
end

@testset "unionrange" begin
    @test_throws MethodError Augmentor.unionrange(1:1:2, 4:5)
    @test @inferred(Augmentor.unionrange(1:5, 2:6)) === 1:6
    @test @inferred(Augmentor.unionrange(2:6, 1:5)) === 1:6
    @test @inferred(Augmentor.unionrange(1:2, 4:5)) === 1:5
    @test @inferred(Augmentor.unionrange(Base.OneTo(2), 4:5)) === 1:5
    @test @inferred(Augmentor.unionrange(1:6, 2:3)) === 1:6
    @test @inferred(Augmentor.unionrange(2:3, 1:6)) === 1:6
end

@testset "_2dborder!" begin
    A = rand(2, 5, 6)
    @test @inferred(Augmentor._2dborder!(A, 0.)) === A
    s = 0.
    ndim, h, w = size(A)
    for i = 1:h, j = 1:w
        if i == 1 || j == 1 || i == h || j == w
            for d = 1:ndim
                s += abs(A[d, i, j])
            end
        end
    end
    @test s == 0.
end


text:
export LREnsemble
export trainerr, trainerrseq, testerr, testerrseq
export normcardseq

mutable struct LREnsemble{F, B<:AbBasis{F}}
    prob::Problem
    ωs::Vector{B}           # Predictor Parameters
    βs::Matrix{F}
    zs::Vector{Vector{F}}   # Predictors
    fit::Vector{F}
    λs::Vector{F}           # Covariances
    δβ::Vector{F}           # Change of dir of β
    η::Vector{F}            # Equiangular direction
    As::Vector{F}           # Equiangles
    γs::Vector{F}           # Arclengths
    Gram::Matrix{F}         # Gram Matrix
    active::BitVector
    γzerohit::Vector{F}
end

function LREnsemble(prob::Problem{MSELoss, B},
                    ω::B,
                    z::Vector{F}=calc(ω, prob.X),
                    C::F=dot(prob.y, z)) where {F, B<:AbBasis{F}}
    normz = norm(z)
    LREnsemble{F, B}(
        prob,
        [ω],
        hcat(0.),
        [z],
        zeros(prob.y),
        [C],
        [1./normz],
        z/normz,
        [normz],
        [0.],
        hcat(normz^2),
        trues(1),
        [Inf])
end

clone(e::LREnsemble) =
    LREnsemble(
        e.prob,
        deepcopy(e.ωs),
        copy(e.βs),
        deepcopy(e.zs),
        copy(e.fit),
        copy(e.λs),
        copy(e.δβ),
        copy(e.η),
        copy(e.As),
        copy(e.γs),
        copy(e.Gram),
        copy(e.active),
        copy(e.γzerohit))

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

function mindist(ωs::Vector{B}, u::Vector{F}, s::F, active::BitVector=trues(length(ωs))) where {F, B<:AbBasis{F}}
    d = Inf
    for (ω, actv) in zip(ωs, active)
        if ω.sgn==s && actv
            d = min(d, norm(u-ω.u))
        end
    end
    d
end
mindist(ωs::Vector{B}, ω::B) where {F, B<:AbBasis{F}} = mindist(ωs, ω.u, ω.sgn)
mindist(ens::LREnsemble{F, B}, u::Vector{F}, s::F) where {F, B<:AbBasis{F}} = mindist(ens.ωs, u, s)#, ens.active)
mindist(ens::LREnsemble{F, B}, ω::B) where {F, B<:AbBasis{F}} = mindist(ens.ωs, ω.u, ω.sgn)#, ens.active)
mindist(us::Vector{Vector{F}}, u::Vector{F}) where F = (length(us)==0)?Inf:minimum(norm(u-ui) for ui in us)
dists(ens::LREnsemble, na=-1.) = [ω1.sgn==ω2.sgn?norm(ω1.u-ω2.u):na for ω1 in ens.ωs, ω2 in ens.ωs]

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

function normcardseq(ens::LREnsemble)
    mapslices(β->norm(β, 1), ens.βs, 1)[:], mapslices(countnz, ens.βs, 1)[:]
end

trainerr(ens::LREnsemble, los=ens.prob.loss) = los(ens.prob.y, ens.fit)
testerr(ens::LREnsemble, X::Matrix, y::Vector, los=ens.prob.loss) = los(y, sum(ens.βs[:, end] .* calc.(ens.ωs, (X,))))

function trainerrseq(ens::LREnsemble{F}, los=ens.prob.loss) where F
    fits = hcat(ens.zs...)*ens.βs
    [los(ens.prob.y, fits[:, i]) for i in 1:size(fits)[2]]
end

function testerrseq(ens::LREnsemble{F}, X::Matrix{F}, y::Vector{F}, los=ens.prob.loss) where F
    zs = hcat(calc.(ens.ωs, (X,))...)
    fits = zs*ens.βs
    [los(y, fits[:, i]) for i in 1:size(fits)[2]]
end

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

function comparecatchup(ens::LREnsemble{F}, zold::Vector{F}, znews::Vector{Vector{F}}) where F<:AbstractFloat
    C, A = ens.λs[end], ens.As[end]
    r = ens.prob.y-ens.fit
    Cold = r'*zold
    Cnews = [r'*znew for znew in znews]
    Aold = ens.η'*zold
    Anews = [ens.η'*znew for znew in znews]
    γold = (C-Cold)/(A-Aold)
    γnews = [(C-Cnews[i])/(A-Anews[i]) for i in 1:length(znews)]
    printfmtln("\t𝒜  :C={:6.4f} A={:6.4f} ", C, A)
    printfmtln("\tOLD:C={:6.4f} A={:6.4f} γ={:6.4f} C(γ)={:6.4f}", Cold, Aold, γold, C-γold*A)
    for i in 1:length(znews)
    printfmtln("\tNEW:C={:6.4f} A={:6.4f} γ={:6.4f} C(γ)={:6.4f}", Cnews[i], Anews[i], γnews[i], C-γnews[i]*A)
    end
end

comparecatchup(ens::LREnsemble{F, B}, zold::Vector{F}, ωnews::Vector{B}) where {F, B<:AbBasis{F}} =
    comparecatchup(ens, zold, [calc(ω, ens.prob.X) for ω in ωnews])

function printens(ens::LREnsemble)
    println("LARS Ensemble")
    println("\tLength ", length(ens.ωs))
    println("\tLoss ", ens.prob.loss(ens.prob.y, ens.fit))
    println("\tλs\t", ens.λs)
    println("\tC(u)\t", [(ens.prob.y-ens.fit)'*z for z in ens.zs])
    println("\tAs\t", ens.As)
    println("\tA(u)\t", [ens.η'*z for z in ens.zs])
    println("\tγs\t", ens.γs)
    pretty("\tβs ", ens.βs)
    println("\tδβ\t", ens.δβ)
    println("\tγ0\t", ens.γzerohit)
    println("\t𝒜\t", ens.active*1)
    println("\tSigns\t", [ω.sgn>0?"+":"-" for ω in ens.ωs])
    pretty("\tGram ", ens.Gram)
    pretty("\tDists ", dists(ens, 1e-3), 3)
    printfmtln("\tClosest\t{:7.5f}", minimum([((ω1.sgn==ω2.sgn)&&(ω1.u!=ω2.u))?norm(ω1.u-ω2.u):Inf for ω1 in ens.ωs, ω2 in ens.ωs]))
end

function pretty(name, M::Matrix{F}, dig=4) where F<:AbstractFloat
    s = reprmime("text/plain", round.(Int, M*10.0^dig))
    s = s[findfirst("\n", s)[1]:end]
    s = replace(s, "\n"=>"\n\t")
    println(name*s)
end


text:
using Test #, Random

@testset "Testing solution to Exercise 2" begin

println("All tests for exercise 2 are in the notebook.")

end; # Exercise 2


text:
function _dlsym(sym::Symbol, libs::Clibrary...)
	for (ind, lib) in enumerate(libs)
		isLast = ind == length(libs)
		handle = Libdl.dlsym(lib.handle, sym, throw_error = isLast)
		isnothing(handle) || return handle
	end
	error("Libdl.dlsym returned a C_NULL handle and did not throw an error")
end

macro cbindings(exprs...) return _cbindings(__module__, exprs...) end

function _cbindings(mod::Module, exprs...)
	length(exprs) >= 1 || error("Expected @cbindings to have (optional) library expressions followed by a block of @cextern statements")

	syms = []
	libs = map(enumerate(exprs[1:end-1])) do (ind, expr)
		sym = gensym("lib")
		push!(syms, Expr(:$, sym))
		return quote
			$(sym) = $(esc(expr))
		end
	end

	function addlibs(e)
		if Base.is_expr(e, :macrocall) && length(e.args) >= 1 && e.args[1] in _externExprs
			append!(e.args, syms)
		elseif e isa Expr
			for expr in e.args
				addlibs(expr)
			end
		end
	end
	block = Base.is_expr(exprs[end], :block) ? exprs[end] : Expr(:block, exprs[end])
	addlibs(block)

	return quote
		let
			$(libs...)
			@eval $(block)
			nothing
		end
	end
end

macro cextern(exprs...) return _cextern(__module__, exprs...) end

function _cextern(mod::Module, sig::Expr, libs::Union{Symbol, Expr, Clibrary}...)
	return _cextern(mod, nothing, sig, libs...)
end

function _cextern(mod::Module, conv::Union{Symbol, Nothing}, sig::Expr, libs::Union{Symbol, Expr, Clibrary}...)
	isempty(libs) && error("Usage of @cextern requires at least 1 library to be provided")
	Base.is_expr(sig, :(::), 2) || error("Expected @cextern to have a function signature, `func(arg1::Arg1Type, ..., argN::ArgNType)::ReturnType`, or a global variable, `variable::VarType`, expression, but found `$(sig)`")
	(externName, externType) = sig.args
	return (externName isa Symbol ? _cextern_variable : _cextern_function)(mod, conv, externName, externType, libs...)
end

function _cextern_function(mod::Module, conv::Union{Symbol, Nothing}, func::Expr, retType::Union{Symbol, Expr}, libs::Union{Symbol, Expr, Clibrary}...)
	isnothing(conv) || (startswith(String(conv), "__") && endswith(String(conv), "__") && length(String(conv)) > 4) || error("Expected @cextern to have a calling convention specified as `__CONVENTION__`, such as `__cdecl__` or `__stdcall__`, but got `$(conv)`")
	conv = isnothing(conv) ? conv : Symbol(String(conv)[3:end-2])

	(Base.is_expr(func, :call) && length(func.args) >= 1) || error("Expected @cextern to have a function `func(...)` in the function signature, but got `$(func)`")
	func.args[1] isa Symbol || error("Expected @cextern to have a Symbol for a function name, but got `$(func)`")
	name = func.args[1]

	argNames = []
	argTypes = map(func.args[2:end]) do arg
		Base.is_expr(arg, :(::), 2) || Base.is_expr(arg, :(...), 1) || error("Expected @cextern to have function arguments specified as `argName::ArgType` or `argName...`, but got `$(arg)`")
		if Base.is_expr(arg, :(::), 2)
			(argName, argType) = arg.args
		else
			(argName, argType) = (arg, :Vararg)
		end
		push!(argNames, argName)
		return argType
	end

	libs = map(l -> l isa Clibrary ? l : esc(l), libs)
	sig = (
		retType,
		:(Tuple{$(argTypes...)}),
		(isnothing(conv) ? () : (:(CBinding.Cconvention{$(QuoteNode(conv))}),))...,
	)

	# NOTE: about quote block:
	#   1. inner @eval the function type in module scope in case it defines new types
	#   2. outer @eval bind the function pointer and "bake" it into the function
	return quote
		@eval $(name)($(argNames...)) = $(Expr(:$, :(
			(@eval CBinding.Cfunction{$(sig...)})($(QuoteNode(name)), $(libs...))
		)))($(argNames...))
	end
end

function _cextern_variable(mod::Module, unused::Nothing, varName::Symbol, varType::Union{Symbol, Expr}, libs::Union{Symbol, Expr, Clibrary}...)
	libs = map(l -> l isa Clibrary ? l : Expr(:$, esc(l)), libs)

	# NOTE: @eval the variable type in module scope in case it defines new types and so it can be const
	return quote
		@eval const $(varName) = CBinding.Cglobal{$(varType)}($(QuoteNode(varName)), $(libs...))
	end
end


text:
using LiteQTL
using Test
using DelimitedFiles
using LinearAlgebra

URL=joinpath(@__DIR__, "..", "data", "HC_M2_0606_R.zip")
output_dir=joinpath(@__DIR__, "..", "data","SPLEEN_CLEAN_DATA")
scan="FALSE"
r_cleaning_script = joinpath(@__DIR__, "..", "r", "cleaning.R")

@test isfile(URL)
@test isfile(r_cleaning_script)

run(`Rscript --vanilla $r_cleaning_script $URL $output_dir $scan`)
#
export_matrix="false"
output_file="julia_result.csv"
# rqtl_file is needed to find gmap.csv.
rqtl_file=URL
julia_scan_script = joinpath(@__DIR__, "..", "bin", "MyApp", "src", "MyApp.jl")

run(`julia $julia_scan_script $output_dir $output_file $rqtl_file $export_matrix`)

# readdlm(file, ','; skipstart=1)[:,2:end-1]
@info "Getting Rqtl scan results"
rqtl_result_file = joinpath(output_dir, "rqtl_result.csv")
rqtl_result = readdlm(rqtl_result_file, ','; skipstart=1)
rqtl_result_float = convert(Array{Float64,2}, rqtl_result[:,2:end])
tmp = zeros(size(rqtl_result_float)[2], size(rqtl_result_float)[1])
rqtl_max_array = find_max_idx_value(transpose!(tmp, rqtl_result_float))

@info "Getting Julia scan results"
julia_result_file = joinpath(output_dir, "julia_result.csv")
julia_result = readdlm(julia_result_file, ','; skipstart=1)[:,4:5]

# compare indexes
@info "Comparing both index and lod score"
rqtl_idx = trunc.(Int, rqtl_max_array[:,1])
julia_idx = convert(Array{Int64,1}, julia_result[:,1])

@info "Do the indices and lod scores match at least 95 percent?"
@test sum(isapprox.(rqtl_idx, julia_idx, atol=3)) >= size(rqtl_idx)[1] * 0.95
@test sum(isapprox.(rqtl_max_array[:,2], abs.(julia_result[:,2]), atol=0.00001)) >= size(rqtl_idx)[1] * 0.95


text:
@testset "Test center_position" begin
    @test (2,) == center_position(ones(3))
    @test (3,) == center_position(ones(4))
    @test (3,) == center_position(ones(5))
    @test (2,3,4) == center_position(ones(3,4,6))
end

@testset "Test ft_center_diff" begin
    @test NDTools.ft_center_diff((3,4,5,6)) == (1,2,2,3)
end

@testset "Test center_value" begin
    x = randn((1,2,3,4))
    @test center_value(x) == x[center_position(x)...]
end

@testset "Test select_sizes" begin
    @test (1, 3, 2) == select_sizes(randn((4,3,2)), (2,3))
    @test (1, 3, 1) == select_sizes(randn((4,3,2)), 2)
    @test (3, 2) == select_sizes_squeeze(randn((4,3,2)), (2,3))
    @test (1, ) == select_sizes_squeeze(randn((1,)), (1,))
    @test (1, ) == select_sizes_squeeze(randn((1,)), 1)
end

@testset "Test single_dim_size" begin
    @test single_dim_size(4, 3) == (1, 1, 1, 3)
    @test single_dim_size(4, 5) == (1, 1, 1, 5)
    @test single_dim_size(2, 5) == (1, 5)
    @test single_dim_size(Val(2), 5) == (1, 5)
    @test single_dim_size(3,5, 4) == (1, 1, 5, 1)

end

@testset "Test reorient" begin
    @test size(reorient([1,2,3,4], 3)) == (1,1,4)
    @test size(reorient([1,2,3,4], Val(3))) == (1,1,4)
    @test reorient([1,2,3,4], 2)[:] == [1,2,3,4]
end

@testset "Test collect_dim" begin
    @test (1,1,1,5) == size(collect_dim((3,4,5,6,7),4))
    @test collect_dim((3,4,5,6,7),4)[1,1,1,3] == 5
    @test collect_dim(1:5,4)[:] == [1,2,3,4,5]
    @test collect_dim(1:5, Val(4))[:] == [1,2,3,4,5]
end

@testset "Test expand_size" begin
    @test (1,2,3,7,8) == Tuple(expand_size((1,2,3), (4,5,6,7,8)))
    @test (1,2) == Tuple(expand_size((1,2,3,4,5,6), (4,5)))
end

@testset "Test expand_add" begin
    @test expand_add((1,2,3),(4,5,6,7,8,9)) == (5, 7, 9, 7, 8, 9)
    @test expand_add((1,1), (0,)) == (1,)
    @test expand_add((1,2), (0,0)) == (1,2)
    @test expand_add((1,2), (0,0, 0)) == (1,2, 0)
end

@testset "Test optional_pos_z" begin
    @test NDTools.optional_pos_z((1,2,3,4), (1,1,1,1)) == 2
    @test NDTools.optional_pos_z((5,5), (1,1)) == 1
    @test NDTools.optional_pos_z((5,), (1,)) == 1
end

@testset "Test curry" begin
    g = NDTools.curry(+,10.0)
    @test g(3) == 13
end


text:
module PandaModels
import JuMP
import InfrastructureModels
import PowerModels
import PowerModels: ids, ref, var, con, sol, nw_id_default
import Memento
import JSON
import Cbc
import Ipopt
import Juniper
const _IM = InfrastructureModels
const _PM = PowerModels
const _PdM = PandaModels
const LOGGER = Memento.getlogger(PowerModels)

export run_powermodels_pf,
    run_powermodels_opf,
    run_powermodels_ots,
    run_powermodels_tnep,
    run_pandamodels_vstab,
    run_pandamodels_multi_vstab,
    run_pandamodels_qflex,
    run_pandamodels_multi_qflex,
    run_pandamodels_vstab_test,
    run_pandamodels_qflex_test

include("input/pp_to_pm.jl")
include("input/tools.jl")
include("models/vstab.jl")
include("models/qflex.jl")
include("models/call_pandamodels.jl")
include("models/call_powermodels.jl")
include("models/run_pm_vstab_dev.jl")
include("models/run_pm_qflex_dev.jl")

end


text:
module HuggingFaceApi

using SHA
using Pkg
using Pkg.MiniProgressBars
using Downloads
import Dates

using OhMyArtifacts

using JSON3
using StructTypes

const my_artifacts = Ref{String}()

const ENDPOINT = Ref("https://huggingface.co")
const token_path = Ref(expanduser("~/.huggingface/token"))

const PYTORCH_WEIGHTS_NAME = Ref("pytorch_model.bin")
const CONFIG_NAME = Ref("config.json")

function hgf_dir()
    return dirname(hgf_artifacts())
end

function hgf_artifacts()
    global my_artifacts
    return my_artifacts[]
end

function env_is_true(key)
    e = get(ENV, key, false)
    e isa Bool && return e
    e isa String && return lowercase(e) in ("1", "on", "yes", "true") ? true : false
    return false
end

function __init__()
    global ENDPOINT, my_artifacts
    # my artifacts
    my_artifacts[] = @my_artifacts_toml!()

    # set 🤗 endpoint
    env_is_true("HUGGINGFACE_CO_STAGING") && (ENDPOINT[] = "https://moon-staging.huggingface.co")
    return
end

"add 🤗 autocomplete in REPL"
function huggingface_emoji()
    REPL = Base.require(Base.PkgId(Base.UUID((0x3fa0cd96_eef1_5676, 0x8a61_b3b8758bbffb)), "REPL"))
    if !haskey(REPL.REPLCompletions.emoji_symbols, "\\:huggingface:")
        REPL.REPLCompletions.emoji_symbols["\\:huggingface:"] = "🤗"
    end
    return :🤗
end

function get_endpoint()
    global ENDPOINT
    return ENDPOINT[]
end

function get_token_path()
    global token_path
    return abspath(token_path[])
end

export hf_hub_download, HuggingFaceURL, cached_download

include("auth.jl")
include("download.jl")
include("api.jl")
#include("datasets.jl")

end


text:
function _is_available_source(x, bus::PSY.Bus)
    return PSY.get_available(x) && x.bus == bus && !isa(x, PSY.ElectricLoad)
end

"""
Calculates the From - To complex power flow (Flow injected at the bus) of branch of type
TapTransformer
"""
function flow_val(b::PSY.TapTransformer)
    !PSY.get_available(b) && return 0.0
    Y_t = PSY.get_series_admittance(b)
    c = 1 / PSY.get_tap(b)
    arc = PSY.get_arc(b)
    V_from = arc.from.magnitude * (cos(arc.from.angle) + sin(arc.from.angle) * 1im)
    V_to = arc.to.magnitude * (cos(arc.to.angle) + sin(arc.to.angle) * 1im)
    I = (V_from * Y_t * c^2) - (V_to * Y_t * c)
    flow = V_from * conj(I)
    return flow
end

"""
Calculates the From - To complex power flow (Flow injected at the bus) of branch of type
Line
"""
function flow_val(b::PSY.ACBranch)
    !PSY.get_available(b) && return 0.0
    Y_t = PSY.get_series_admittance(b)
    arc = PSY.get_arc(b)
    V_from = arc.from.magnitude * (cos(arc.from.angle) + sin(arc.from.angle) * 1im)
    V_to = arc.to.magnitude * (cos(arc.to.angle) + sin(arc.to.angle) * 1im)
    I = V_from * (Y_t + (1im * PSY.get_b(b).from)) - V_to * Y_t
    flow = V_from * conj(I)
    return flow
end

"""
Calculates the From - To complex power flow (Flow injected at the bus) of branch of type
Transformer2W
"""
function flow_val(b::PSY.Transformer2W)
    !PSY.get_available(b) && return 0.0
    Y_t = PSY.get_series_admittance(b)
    arc = PSY.get_arc(b)
    V_from = arc.from.magnitude * (cos(arc.from.angle) + sin(arc.from.angle) * 1im)
    V_to = arc.to.magnitude * (cos(arc.to.angle) + sin(arc.to.angle) * 1im)
    I = V_from * (Y_t + (1im * PSY.get_primary_shunt(b))) - V_to * Y_t
    flow = V_from * conj(I)
    return flow
end

function flow_val(b::PSY.PhaseShiftingTransformer)
    error("Systems with PhaseShiftingTransformer not supported yet")
    return
end

"""
Calculates the From - To complex power flow using external data of voltages of branch of type
TapTransformer
"""
function flow_func(b::PSY.TapTransformer, V_from::Complex{Float64}, V_to::Complex{Float64})
    !PSY.get_available(b) && return (0.0, 0.0)
    Y_t = PSY.get_series_admittance(b)
    c = 1 / PSY.get_tap(b)
    I = (V_from * Y_t * c^2) - (V_to * Y_t * c)
    flow = V_from * conj(I)
    return real(flow), imag(flow)
end

"""
Calculates the From - To complex power flow using external data of voltages of branch of type
Line
"""
function flow_func(b::PSY.ACBranch, V_from::Complex{Float64}, V_to::Complex{Float64})
    !PSY.get_available(b) && return (0.0, 0.0)
    Y_t = PSY.get_series_admittance(b)
    I = V_from * (Y_t + (1im * PSY.get_b(b).from)) - V_to * Y_t
    flow = V_from * conj(I)
    return real(flow), imag(flow)
end

"""
Calculates the From - To complex power flow using external data of voltages of branch of type
Transformer2W
"""
function flow_func(b::PSY.Transformer2W, V_from::Complex{Float64}, V_to::Complex{Float64})
    !PSY.get_available(b) && return (0.0, 0.0)
    Y_t = PSY.get_series_admittance(b)
    I = V_from * (Y_t + (1im * PSY.get_primary_shunt(b))) - V_to * Y_t
    flow = V_from * conj(I)
    return real(flow), imag(flow)
end

function flow_func(
    b::PSY.PhaseShiftingTransformer,
    V_from::Complex{Float64},
    V_to::Complex{Float64},
)
    error("Systems with PhaseShiftingTransformer not supported yet")
    return
end

"""
Updates the flow on the branches
"""
function _update_branch_flow!(sys::PSY.System)
    for b in PSY.get_components(PSY.ACBranch, sys)
        S_flow = PSY.get_available(b) ? flow_val(b) : 0.0 + 0.0im
        PSY.set_active_power_flow!(b, real(S_flow))
        PSY.set_reactive_power_flow!(b, imag(S_flow))
    end
end

"""
Obtain total load on bus b
"""
function _get_load_data(sys::PSY.System, b::PSY.Bus)
    active_power = 0.0
    reactive_power = 0.0
    for l in PSY.get_components(PSY.ElectricLoad, sys, x -> !isa(x, PSY.FixedAdmittance))
        !PSY.get_available(l) && continue
        if (l.bus == b)
            active_power += PSY.get_active_power(l)
            reactive_power += PSY.get_reactive_power(l)
        end
    end
    return active_power, reactive_power
end

function _get_fixed_admittance_power(
    sys::PSY.System,
    b::PSY.Bus,
    result::AbstractVector,
    ix::Int,
)
    active_power = 0.0
    reactive_power = 0.0
    for l in PSY.get_components(PSY.FixedAdmittance, sys)
        !PSY.get_available(l) && continue
        if (l.bus == b)
            Vm_squared =
                b.bustype == PSY.BusTypes.PQ ? result[2 * ix - 1]^2 : PSY.get_magnitude(b)^2
            active_power += Vm_squared * real(PSY.get_Y(l))
            reactive_power += Vm_squared * imag(PSY.get_Y(l))
        end
    end
    return active_power, reactive_power
end

function _power_redistribution_ref(
    sys::PSY.System,
    P_gen::Float64,
    Q_gen::Float64,
    bus::PSY.Bus,
)
    devices_ =
        PSY.get_components(PSY.StaticInjection, sys, x -> _is_available_source(x, bus))
    if length(devices_) == 1
        device = first(devices_)
        PSY.set_active_power!(device, P_gen)
        PSY.set_reactive_power!(device, Q_gen)
        return
    elseif length(devices_) > 1
        devices = sort(collect(devices_); by=x -> PSY.get_max_active_power(x))
    else
        error("No devices in bus $(PSY.get_name(bus))")
    end

    sum_basepower = sum(PSY.get_max_active_power.(devices))
    p_residual = P_gen
    units_at_limit = Vector{Int}()
    for (ix, d) in enumerate(devices)
        p_limits = PSY.get_active_power_limits(d)
        part_factor = p_limits.max / sum_basepower
        p_frac = P_gen * part_factor
        p_set_point = clamp(p_frac, p_limits.min, p_limits.max)
        if (p_frac >= p_limits.max + BOUNDS_TOLERANCE) ||
           (p_frac <= p_limits.min - BOUNDS_TOLERANCE)
            push!(units_at_limit, ix)
            @warn "Unit $(PSY.get_name(d)) set at the limit $(p_set_point). P_max = $(p_limits.max) P_min = $(p_limits.min)"
        end
        PSY.set_active_power!(d, p_set_point)
        p_residual -= p_set_point
    end

    if !isapprox(p_residual, 0.0, atol=ISAPPROX_ZERO_TOLERANCE)
        removed_power = sum(PSY.get_max_active_power.(devices[units_at_limit]))
        reallocated_p = 0.0
        it = 0
        while !isapprox(p_residual, 0.0, atol=ISAPPROX_ZERO_TOLERANCE)
            if length(devices) == length(units_at_limit) + 1
                @warn "all devices at the active Power Limit"
                break
            end
            for (ix, d) in enumerate(devices)
                ix ∈ units_at_limit && continue
                p_limits = PSY.get_active_power_limits(d)
                part_factor = p_limits.max / (sum_basepower - removed_power)
                p_frac = p_residual * part_factor
                current_p = PSY.get_active_power(d)
                p_set_point = p_frac + current_p
                if (p_set_point >= p_limits.max + BOUNDS_TOLERANCE) ||
                   (p_set_point <= p_limits.min - BOUNDS_TOLERANCE)
                    push!(units_at_limit, ix)
                    @warn "Unit $(PSY.get_name(d)) set at the limit $(p_set_point). P_max = $(p_limits.max) P_min = $(p_limits.min)"
                end
                p_set_point = clamp(p_set_point, p_limits.min, p_limits.max)
                PSY.set_active_power!(d, p_set_point)
                reallocated_p += p_frac
            end
            p_residual -= reallocated_p
            if isapprox(p_residual, 0, atol=ISAPPROX_ZERO_TOLERANCE)
                break
            end
            it += 1
            if it > 10
                break
            end
        end
        if !isapprox(p_residual, 0.0, atol=ISAPPROX_ZERO_TOLERANCE)
            remaining_unit_index = setdiff(1:length(devices), units_at_limit)
            @assert length(remaining_unit_index) == 1 remaining_unit_index
            device = devices[remaining_unit_index[1]]
            @debug "Remaining residual $q_residual, $(PSY.get_name(bus))"
            p_set_point = PSY.get_active_power(device) + p_residual
            PSY.set_active_power!(device, p_set_point)
            p_limits = PSY.get_reactive_power_limits(device)
            if (p_set_point >= p_limits.max + BOUNDS_TOLERANCE) ||
               (p_set_point <= p_limits.min - BOUNDS_TOLERANCE)
                @error "Unit $(PSY.get_name(device)) P=$(p_set_point) above limits. P_max = $(p_limits.max) P_min = $(p_limits.min)"
            end
        end
    end
    _reactive_power_redistribution_pv(sys, Q_gen, bus)
    return
end

function _reactive_power_redistribution_pv(sys::PSY.System, Q_gen::Float64, bus::PSY.Bus)
    @debug "Reactive Power Distribution $(PSY.get_name(bus))"
    devices_ =
        PSY.get_components(PSY.StaticInjection, sys, x -> _is_available_source(x, bus))

    if length(devices_) == 1
        @debug "Only one generator in the bus"
        PSY.set_reactive_power!(first(devices_), Q_gen)
        return
    elseif length(devices_) > 1
        devices = sort(collect(devices_); by=x -> PSY.get_max_reactive_power(x))
    else
        error("No devices in bus $(PSY.get_name(bus))")
    end

    total_active_power = sum(PSY.get_active_power.(devices))

    if isapprox(total_active_power, 0.0, atol=ISAPPROX_ZERO_TOLERANCE)
        @debug "Total Active Power Output at the bus is $(total_active_power). Using Unit's Base Power"
        sum_basepower = sum(PSY.get_base_power.(devices))
        for d in devices
            part_factor = PSY.get_base_power(d) / sum_basepower
            PSY.set_reactive_power!(d, Q_gen * part_factor)
        end
        return
    end

    q_residual = Q_gen
    units_at_limit = Vector{Int}()

    for (ix, d) in enumerate(devices)
        q_limits = PSY.get_reactive_power_limits(d)
        if isapprox(q_limits.max, 0.0, atol=BOUNDS_TOLERANCE) &&
           isapprox(q_limits.min, 0.0, atol=BOUNDS_TOLERANCE)
            push!(units_at_limit, ix)
            @info "Unit $(PSY.get_name(d)) has no Q control capability. Q_max = $(q_limits.max) Q_min = $(q_limits.min)"
            continue
        end

        fraction = PSY.get_active_power(d) / total_active_power

        if fraction == 0.0
            PSY.set_reactive_power!(d, 0.0)
            continue
        else
            @assert fraction > 0
        end

        q_frac = Q_gen * fraction
        q_set_point = clamp(q_frac, q_limits.min, q_limits.max)

        if (q_frac >= q_limits.max + BOUNDS_TOLERANCE) ||
           (q_frac <= q_limits.min - BOUNDS_TOLERANCE)
            push!(units_at_limit, ix)
            @warn "Unit $(PSY.get_name(d)) set at the limit $(q_set_point). Q_max = $(q_limits.max) Q_min = $(q_limits.min)"
        end

        PSY.set_reactive_power!(d, q_set_point)
        q_residual -= q_set_point

        if isapprox(q_residual, 0.0, atol=ISAPPROX_ZERO_TOLERANCE)
            break
        end
    end

    if !isapprox(q_residual, 0.0, atol=ISAPPROX_ZERO_TOLERANCE)
        it = 0
        while !isapprox(q_residual, 0.0, atol=ISAPPROX_ZERO_TOLERANCE)
            if length(devices) == length(units_at_limit) + 1
                @debug "Only one device not at the limit in Bus"
                break
            end
            removed_power = sum(PSY.get_active_power.(devices[units_at_limit]))
            reallocated_q = 0.0
            for (ix, d) in enumerate(devices)
                ix ∈ units_at_limit && continue
                q_limits = PSY.get_reactive_power_limits(d)

                if removed_power < total_active_power
                    fraction =
                        PSY.get_active_power(d) / (total_active_power - removed_power)
                elseif isapprox(removed_power, total_active_power)
                    fraction = 1
                else
                    error("Remove power can't be larger than the total active power")
                end

                if fraction == 0.0
                    continue
                else
                    PSY.InfrastructureSystems.@assert_op fraction > 0
                end

                current_q = PSY.get_reactive_power(d)
                q_frac = q_residual * fraction
                q_set_point = clamp(q_frac + current_q, q_limits.min, q_limits.max)
                reallocated_q += q_frac
                if (q_frac >= q_limits.max + BOUNDS_TOLERANCE) ||
                   (q_frac <= q_limits.min - BOUNDS_TOLERANCE)
                    push!(units_at_limit, ix)
                    @warn "Unit $(PSY.get_name(d)) set at the limit $(q_set_point). Q_max = $(q_limits.max) Q_min = $(q_limits.min)"
                end

                PSY.set_reactive_power!(d, q_set_point)
            end
            q_residual -= reallocated_q
            if isapprox(q_residual, 0, atol=ISAPPROX_ZERO_TOLERANCE)
                break
            end
            it += 1
            if it > 1
                break
            end
        end
    end

    if !isapprox(q_residual, 0.0, atol=ISAPPROX_ZERO_TOLERANCE)
        remaining_unit_index = setdiff(1:length(devices), units_at_limit)
        @assert length(remaining_unit_index) == 1 remaining_unit_index
        device = devices[remaining_unit_index[1]]
        @debug "Remaining residual $q_residual, $(PSY.get_name(bus))"
        q_set_point = PSY.get_reactive_power(device) + q_residual
        PSY.set_reactive_power!(device, q_set_point)
        q_limits = PSY.get_reactive_power_limits(device)
        if (q_set_point >= q_limits.max + BOUNDS_TOLERANCE) ||
           (q_set_point <= q_limits.min - BOUNDS_TOLERANCE)
            @error "Unit $(PSY.get_name(device)) Q=$(q_set_point) above limits. Q_max = $(q_limits.max) Q_min = $(q_limits.min)"
        end
    end

    return
end

"""
Updates system voltages and powers with power flow results
"""
function write_powerflow_solution!(sys::PSY.System, result::Vector{Float64})
    buses = enumerate(
        sort!(collect(PSY.get_components(PSY.Bus, sys)), by=x -> PSY.get_number(x)),
    )

    for (ix, bus) in buses
        if bus.bustype == PSY.BusTypes.REF
            P_gen = result[2 * ix - 1]
            Q_gen = result[2 * ix]
            _power_redistribution_ref(sys, P_gen, Q_gen, bus)
        elseif bus.bustype == PSY.BusTypes.PV
            Q_gen = result[2 * ix - 1]
            bus.angle = result[2 * ix]
            _reactive_power_redistribution_pv(sys, Q_gen, bus)
        elseif bus.bustype == PSY.BusTypes.PQ
            Vm = result[2 * ix - 1]
            θ = result[2 * ix]
            PSY.set_magnitude!(bus, Vm)
            PSY.set_angle!(bus, θ)
        end
    end

    _update_branch_flow!(sys)
    return
end

"""
Return power flow results in dictionary of dataframes.
"""
function write_results(sys::PSY.System, result::Vector{Float64})
    @info("Voltages are exported in pu. Powers are exported in MW/MVAr.")
    buses = sort!(collect(PSY.get_components(PSY.Bus, sys)), by=x -> PSY.get_number(x))
    N_BUS = length(buses)
    bus_map = Dict(buses .=> 1:N_BUS)
    sys_basepower = PSY.get_base_power(sys)
    sources = PSY.get_components(PSY.StaticInjection, sys, d -> !isa(d, PSY.ElectricLoad))
    Vm_vect = fill(0.0, N_BUS)
    θ_vect = fill(0.0, N_BUS)
    P_gen_vect = fill(0.0, N_BUS)
    Q_gen_vect = fill(0.0, N_BUS)
    P_load_vect = fill(0.0, N_BUS)
    Q_load_vect = fill(0.0, N_BUS)

    for (ix, bus) in enumerate(buses)
        P_load_vect[ix], Q_load_vect[ix] = _get_load_data(sys, bus) .* sys_basepower
        P_admittance, Q_admittance = _get_fixed_admittance_power(sys, bus, result, ix)
        P_load_vect[ix] += P_admittance
        Q_load_vect[ix] += Q_admittance
        if bus.bustype == PSY.BusTypes.REF
            Vm_vect[ix] = PSY.get_magnitude(bus)
            θ_vect[ix] = PSY.get_angle(bus)
            P_gen_vect[ix] = result[2 * ix - 1] * sys_basepower
            Q_gen_vect[ix] = result[2 * ix] * sys_basepower
        elseif bus.bustype == PSY.BusTypes.PV
            Vm_vect[ix] = PSY.get_magnitude(bus)
            θ_vect[ix] = result[2 * ix]
            for gen in sources
                !PSY.get_available(gen) && continue
                if gen.bus == bus
                    P_gen_vect[ix] += PSY.get_active_power(gen) * sys_basepower
                end
            end
            Q_gen_vect[ix] = result[2 * ix - 1] * sys_basepower
        elseif bus.bustype == PSY.BusTypes.PQ
            Vm_vect[ix] = result[2 * ix - 1]
            θ_vect[ix] = result[2 * ix]
            for gen in sources
                !PSY.get_available(gen) && continue
                if gen.bus == bus
                    P_gen_vect[ix] += PSY.get_active_power(gen) * sys_basepower
                    Q_gen_vect[ix] += PSY.get_reactive_power(gen) * sys_basepower
                end
            end
        end
    end

    branches = PSY.get_components(PSY.ACBranch, sys)
    N_BRANCH = length(branches)
    P_from_to_vect = fill(0.0, N_BRANCH)
    Q_from_to_vect = fill(0.0, N_BRANCH)
    P_to_from_vect = fill(0.0, N_BRANCH)
    Q_to_from_vect = fill(0.0, N_BRANCH)
    for (ix, b) in enumerate(branches)
        !PSY.get_available(b) && continue
        bus_f_ix = bus_map[PSY.get_arc(b).from]
        bus_t_ix = bus_map[PSY.get_arc(b).to]
        V_from = Vm_vect[bus_f_ix] * (cos(θ_vect[bus_f_ix]) + sin(θ_vect[bus_f_ix]) * 1im)
        V_to = Vm_vect[bus_t_ix] * (cos(θ_vect[bus_t_ix]) + sin(θ_vect[bus_t_ix]) * 1im)
        P_from_to_vect[ix], Q_from_to_vect[ix] = flow_func(b, V_from, V_to) .* sys_basepower
        P_to_from_vect[ix], Q_to_from_vect[ix] = flow_func(b, V_to, V_from) .* sys_basepower
    end

    bus_df = DataFrames.DataFrame(
        bus_number=PSY.get_number.(buses),
        Vm=Vm_vect,
        θ=θ_vect,
        P_gen=P_gen_vect,
        P_load=P_load_vect,
        P_net=P_gen_vect - P_load_vect,
        Q_gen=Q_gen_vect,
        Q_load=Q_load_vect,
        Q_net=Q_gen_vect - Q_load_vect,
    )

    branch_df = DataFrames.DataFrame(
        line_name=PSY.get_name.(branches),
        bus_from=PSY.get_number.(PSY.get_from.(PSY.get_arc.(branches))),
        bus_to=PSY.get_number.(PSY.get_to.(PSY.get_arc.(branches))),
        P_from_to=P_from_to_vect,
        Q_from_to=Q_from_to_vect,
        P_to_from=P_to_from_vect,
        Q_to_from=Q_to_from_vect,
        P_losses=P_from_to_vect + P_to_from_vect,
        Q_losses=Q_from_to_vect + Q_to_from_vect,
    )
    DataFrames.sort!(branch_df, [:bus_from, :bus_to])

    return Dict("bus_results" => bus_df, "flow_results" => branch_df)
end


text:
module IO

export  InpStiffness,
        save_mesh

include(joinpath("INP", "INP.jl"))
using .INP

include("VTK.jl")
using .VTK

end