text:
#======================================
Methods using semidefinite programming
======================================#
using SumOfSquares, DynamicPolynomials, SemialgebraicSets, SDPA

"""
    new_sos(backend; kwargs...)

Return a new (empty) sum-of-squares optimization problem for the given backend.

### Input

- `backend` -- the backend (also called JuMP `solver`)
- `kwargs`  -- additional keyword arguments

### Output

An instance of `SOSModel` for the given backend and options.
"""
@inline function new_sos(backend, kwargs...)
    𝑂 = Dict(kwargs)

    if VERSION < v"1.6"
        if :QUIET ∈ keys(𝑂)
            # for mosek solver
            SOSModel(with_optimizer(backend, QUIET=𝑂[:QUIET]))
        else
            SOSModel(with_optimizer(backend))
        end
    else
        if :QUIET ∈ keys(𝑂)
            # for mosek solver
            SOSModel(backend, QUIET=𝑂[:QUIET])
        else
            SOSModel(backend)
        end
    end
end

"""
    enclose_SumOfSquares(f::Function, dom::Interval_or_IntervalBox,
                         order::Int=5, backend=SDPA.Optimizer; kwargs...)

Compute a range enclosure using sum-of-squares optimization.

### Input

- `f`       -- function
- `dom`     -- hyperrectangular domain, either a unidimensional `Interval` or
               a multidimensional `IntervalBox`
- `order`   -- (optional, default: `5`) maximum degree of the SDP relaxation
- `backend` -- (optional, default: `SDPA.Optimizer`) the optimization backend
               (aka JuMP solver)
- `kwargs`  -- additional keyword arguments

### Output

An interval representing the range enclosure (minimum and maximum) of `f` over
its domain `dom`.

### Algorithm

The range enclosure is computed using polynomial optimization methods.
We refer to the documentation and examples of `SumOfSquares.jl` for details.

### Notes

Use the `backend` keyword argument to pass an SDP solver backend of your choice;
additional arguments to the solver can be passed in `kwargs`.

For instance, to use `SDPA` (default choice), use it as:

```julia
julia> backend = SDPA.Optimizer

julia> enclose_SumOfSquares(f, dom, order, backend=backend)
...
```

To use `MOSEK` in non-verbose mode, let

```julia
julia> using MosekTools

julia> backend = MosekTools.Mosek.Optimizer;

julia> enclose_SumOfSquares(f, dom, order, backend=backend, QUIET=true)
...
```

To get the runtime, use `MOI.get(model, MOI.SolveTime())`.
"""
function enclose_SumOfSquares(f::Function, dom::Interval_or_IntervalBox;
                              order::Int=5, backend=SDPA.Optimizer, kwargs...)
    _enclose_SumOfSquares(f, dom, order, backend; kwargs...)
end

function _enclose_SumOfSquares(f::Function, dom::Interval, order::Int,
                               backend; kwargs...)

    # polynomial variables
    @polyvar x
    p = f(x)

    # box constraints
    B = @set inf(dom) <= x && x <= sup(dom)

    # ============
    # Upper bound
    # ============
    model = new_sos(backend, kwargs...)
    @variable(model, γ) # JuMP decision variable
    @constraint(model, p <= γ, domain=B, maxdegree=order)
    @objective(model, Min, γ)
    optimize!(model)
    upper_bound = objective_value(model)

    # ============
    # Lower bound
    # ============
    model = new_sos(backend, kwargs...)
    @variable(model, γ) # JuMP decision variable
    @constraint(model, p >= γ, domain=B, maxdegree=order)
    @objective(model, Max, γ)
    optimize!(model)
    lower_bound = objective_value(model)

    return Interval(lower_bound, upper_bound)
end

function _enclose_SumOfSquares(f::Function, dom::IntervalBox{N}, order::Int,
                               backend; kwargs...) where {N}

    # polynomial variables
    @polyvar x[1:N]
    p = f(x...)

    # box constraints
    Bi =[@set inf(dom[i]) <= x[i] && x[i] <= sup(dom[i]) for i in 1:N]
    B = reduce(intersect, Bi)

    # ============
    # Upper bound
    # ============
    model = new_sos(backend, kwargs...)
    @variable(model, γ) # JuMP decision variable
    @constraint(model, p <= γ, domain=B, maxdegree=order)
    @objective(model, Min, γ)
    optimize!(model)
    upper_bound = objective_value(model)

    # ============
    # Lower bound
    # ============
    model = new_sos(backend, kwargs...)
    @variable(model, γ) # JuMP decision variable
    @constraint(model, p >= γ, domain=B, maxdegree=order)
    @objective(model, Max, γ)
    optimize!(model)
    lower_bound = objective_value(model)

    return Interval(lower_bound, upper_bound)
end


text:
export free_cell_solver, display_solution

"""
    free_cell_solver(T::Tableau; depth::Int=1, verbose::Int=0)

Solve a free_cell game `T` using the `guided_path_finder` function from
`ImplicitGraphs`.

* `depth` is the amount of look ahead
* `verbose` sets reporting frequency (0 for no reporting)
"""
function free_cell_solver(S::Tableau; depth::Int = 1, verbose::Int = 0)
    G = free_cell_graph()
    T = Victory()
    X = guided_path_finder(G, S, T, score = score, depth = depth, verbose = verbose)
    return X
end

"""
    display_solution(X::Vector{Tableau}, pause=0.5)

Print out a list of free_cell positions (presumably output from `free_cell_solver`)
with `pause` seconds of delay between each step.
"""
function display_solution(X::Vector{Tableau}, pause = 0.5)
    n = length(X)
    for j = 1:n
        println("Step $j of $n")
        println(X[j])
        println()
        sleep(pause)
    end
    nothing
end

function solver_test(T::Tableau, max_depth::Int = 2)
    XX = []
    for d = 1:max_depth
        println("Solving at depth = $d")
        @time X = free_cell_solver(T, depth = d)
        println("Solution has $(length(X)) moves\n")
        push!(XX, X)
    end
    return XX
end

ext: jl


text:
using Combinatorics

function parse_line(l)
    lhs, rhs = split(l ," | ")
    digits = split(lhs)
    output = split(rhs)
    digits, output
end

const data = parse_line.(readlines())

# Part One

sum(data) do (_, out)
    count(∈((2, 3, 4, 7)), length.(out))
end |> println

# Mapping from segments on to corresponding digit
const segments2digit = Dict(
    [1, 2, 3, 5, 6, 7]    => 0,
    [3, 6]                => 1,
    [1, 3, 4, 5, 7]       => 2,
    [1, 3, 4, 6, 7]       => 3,
    [2, 3, 4, 6]          => 4,
    [1, 2, 4, 6, 7]       => 5,
    [1, 2, 4, 5, 6, 7]    => 6,
    [1, 3, 6]             => 7,
    [1, 2, 3, 4, 5, 6, 7] => 8,
    [1, 2, 3, 4, 6, 7]    => 9
)

is_decode_key(perm::Vector, digits) = all(
            sort([findfirst(==(seg), perm) for seg in digit]) ∈ keys(segments2digit)
            for digit in digits
        )

decode(digit, perm::Vector) = get(segments2digit, sort([findfirst(==(seg), perm) for seg in digit]), nothing)

function decode(digits::Vector)
    for perm in permutations("abcdefg")
        if is_decode_key(perm, digits)
            return perm
        end
    end
    nothing
end

# Part Two

sum(data) do (digits, output)
    perm_key = decode(digits)
    sum(decode(d, perm_key)*10^i for (i, d) in zip(3:-1:0, output))
end |> println

max_stars_repo_name: pLOPeGG/JuliaOfCode

max_issues_repo_name: pLOPeGG/JuliaOfCode

max_forks_repo_name: pLOPeGG/JuliaOfCode


text:
using ModelingToolkit
using Test

@Param t
@Unknown x(t)
@Unknown y(t)
@Unknown z(t)
x1 = Unknown(:x ,dependents = [t])
y1 = Unknown(:y ,dependents = [t])
z1 = Unknown(:z ,dependents = [t])
@test x1 == x
@test y1 == y
@test z1 == z
@test convert(Expr, x) == :x
@test convert(Expr, y) == :y
@test convert(Expr, z) == :z

@Param begin
    t
    s
end
t1 = Parameter(:t)
s1 = Parameter(:s)
@test t1 == t
@test s1 == s
@test convert(Expr, t) == :t
@test convert(Expr, s) == :s
@test convert(Expr, cos(t + sin(s))) == :(cos(t + sin(s)))

@Deriv D''~t
D1 = Differential(t, 2)
@test D1 == D
@test convert(Expr, D) == D

@Const c=0 v=2
c1 = Constant(0)
v1 = Constant(2)
@test c1 == c
@test v1 == v
@test convert(Expr, c) == 0
@test convert(Expr, v) == 2


text:
include("znccUtils.jl")

#=
	2D IMPLEMENTATION
=#

#  Lewis, J.P.: Fast normalized cross-correlation. Ind. Light Magic10(2001)
function crossCorrelation!( ::ZNCC, ccr::A{T,2}, num::A{T,2},
                                    padF::A{C{T},2}, padG::A{C{T},2}, plan, iplan,
                                    sumF::A{U,2}, sumF2::A{U,2}, meanF, meanF2,
                                    stdG, sizeF, sizeG, shifts
                          ) where {T<:Real,U<:Real}

    # Computing numerator with FFT.
    crossCorrelation!( FFT(), padF, padG, ccr, plan, iplan );
    Base.circshift!( num, ccr, shifts )

    # Computing denominator through summed-area tables.
    N = prod(sizeG);
    N2F, stdI = N*meanF2, stdG*sqrt(N)
    @inbounds for col in 1:size(ccr,2)
        for row in 1:size(ccr,1)
            F, F2 = znccSumValues( sumF, sumF2, row, col, sizeF, sizeG )
            ccr[ row, col ] = num[ row, col ]/( sqrt( F2 + N2F - 2*F*meanF )*stdG );
        end
    end
end

function PIV_2D( ::ZNCC, img1::A{T,2}, img2::A{T,2},
                         IA::III, SM::III, overlap::III, mpass::I, width::I,
                         peak::S, sigNoise::S, filtFun::Function, threshold::F,
                         mpFun::Function, reps::Array{I,1};
                         corrType=Float32
			   ) where {T<:Real}

	# Calculating vector field size.
    h, w   = size( img1 );
    step   = IA .- overlap;
    VFsize = 1 .+ floor.( Int64, ((h,w,0).-IA)./step );
    U  = zeros( corrType, VFsize[1:2] );
    V  = zeros( corrType, VFsize[1:2] );
    SN = zeros( corrType, VFsize[1:2] );

	ignoreSN = sigNoise == "";

    @inbounds for mp in mpass:-1:1

        # Scaling IA, SM and step (overlap) to multi-pass iteration.
        IA_mp   = mpFun.(  IA , mp );
        SM_mp   = mpFun.(  SM , mp );
        step_mp = mpFun.( step, mp );
        last_mp = ( mp == 1 );

        ssize   = IA_mp .+ 2 .* SM_mp
        search  = zeros( T, ssize[1:2] );
        sumS    = zeros( T, ssize[1:2] .+ 1 );
        sumS2   = zeros( T, ssize[1:2] .+ 1 );

        csize   = 2 .* ( IA_mp .+ SM_mp  );
        cmatrix = zeros( corrType, csize[1:2] );
        num     = zeros( corrType, csize[1:2] );
		shifts  = div.( csize[1:2], 2 ) .+ SM_mp[1:2] .- 1;
        pads    = zeros( Complex{corrType}, csize[1:2] );
        padi    = zeros( Complex{corrType}, csize[1:2] );
        plan    =  FFTW.plan_fft!( padi );
        iplan   = FFTW.plan_ifft!( padi );

        # Constructing grid coordinates for each interrogation volume,
        grid = constructGrid( (h, w, IA_mp[3]), IA_mp, step_mp, mp )

        for g in grid
            #  g[1],  g[2] are start and stop rows of interrogation area
            #  g[3],  g[4] are start and stop columns of interrogation area
            #  g[7],  g[8] are start and stop rows inside UV arrays
            #  g[9], g[10] are start and stop columns inside UV arrays

            interr = view( img1, g[1]:g[2], g[3]:g[4] )

            if threshold > 0 && filtFun( interr ) < threshold
                continue;
            end

            # 1-. Shifting the searchArea ranges by previously calculated displacements
            offH = round(Int64, U[ g[7], g[9] ]);
            offW = round(Int64, V[ g[7], g[9] ]);

            # 1.2-. Search volume coordinates after shifting.
            scoords, so = searchGrid( g, SM, offH, offW, 0, h, w, 0 )

            setTo0!( search )
            putWithinSearch!( search, img2, so, scoords )

            setTo0!( sumS  );
            setTo0!( sumS2 );
            znccSumArrays!( search, sumS, sumS2 );

            meanI = Statistics.mean( interr );
            meanS = Statistics.mean( search );

            setTo0!( pads )
            setTo0!( padi )
            putWithinPadded!( padi, interr, meanI );
            putWithinPadded!( pads, search, meanS );

            meanS2 = meanS*meanS;
            meanI2 = meanI*meanI;
            stdI   = Statistics.std( interr, mean=meanI );

            # 2-. Cross-Correlation
            crossCorrelation!( ZNCC(), cmatrix, num, pads, padi, plan, iplan,
                               sumS, sumS2, meanS, meanS2, stdI, ssize, IA_mp, shifts );

            # 3-. Calculation of displacement
            ( r, c ) = approxTranslation( cmatrix, peak, Val(last_mp) )
			if ( !ignoreSN && last_mp )
            	SN[ g[7], g[9] ] = SNRatio( cmatrix, sigNoise, width )
			end

            # 4-. Updating U, V matrices
            t1 =  g[7]:g[8];
            t2 =  g[9]:g[10];
            U[ t1, t2 ] .= offH - r;
            V[ t1, t2 ] .= offW - c;
        end
    end # Multi-pass loop

    return U, V, SN
end

#=
	3D IMPLEMENTATION
=#

function crossCorrelation!( ::ZNCC, ccr::A{T,3}, num::A{T,3},
                                    padF::A{C{T},3}, padG::A{C{T},3}, plan, iplan,
                                    sumF::A{U,3}, sumF2::A{U,3}, meanF, meanF2,
                                    stdG, sizeF, sizeG, shifts
						  ) where {T<:Real, U<:Real}

    # Computing numerator with FFT.
    crossCorrelation!( FFT(), padF, padG, ccr, plan, iplan );
    Base.circshift!( num, ccr, shifts )

    # Computing denominator through summed-area tables.
    N = prod(sizeG);
    N2F, stdG  = N*meanF2, stdG*sqrt(N);
    @inbounds for z in 1:size(ccr,3)
        for c in 1:size(ccr,2)
            for r in 1:size(ccr,1)
                F, F2 = znccSumValues( sumF, sumF2, r, c, z, sizeF, sizeG )
                ccr[r,c,z] = num[r,c,z]/( sqrt( F2 + N2F - 2*F*meanF )*stdG );
            end
        end
    end
end

function PIV_3D( ::ZNCC, vol1::A{T,3}, vol2::A{T,3},
                         IA::III, SM::III, overlap::III, mpass::I, width::I,
                         peak::S, sigNoise::S, filtFun::Function, threshold::F,
                         mpFun::Function, reps::Array{I,1};
                         corrType=Float32
			   ) where {T<:Real}

	# Calculating vector field size.
    h, w, d = size( vol1 );
    step   = IA .- overlap;
    VFsize = 1 .+ floor.(Int64,(size(vol1).-IA)./step);
    U  = zeros( corrType, VFsize );
    V  = zeros( corrType, VFsize );
    W  = zeros( corrType, VFsize );
    SN = zeros( corrType, VFsize );

	ignoreSN = sigNoise == "";

    @inbounds for mp in mpass:-1:1

        # Scaling IA, SM and step (overlap) to multi-pass iteration.
        IA_mp   = mpFun.(  IA , mp );
        SM_mp   = mpFun.(  SM , mp );
        step_mp = mpFun.( step, mp );
        last_mp = ( mp == 1 );

        ssize   = IA_mp .+ 2 .* SM_mp
        search  = zeros( T, ssize );
        sumS    = zeros( T, ssize .+ 1 );
        sumS2   = zeros( T, ssize .+ 1 );

        # Initialize cross-correlation variables once per multi-pass iteration
        csize   = 2 .* ( IA_mp .+ SM_mp  );
        cmatrix = zeros( corrType, csize );
        num     = zeros( corrType, csize );
        shifts  = div.(csize,2) .+ SM_mp .- 1;
        pads    = zeros( Complex{corrType}, csize );
        padi    = zeros( Complex{corrType}, csize );
        plan    =  FFTW.plan_fft!( pads );
        iplan   = FFTW.plan_ifft!( pads );

        # Constructing grid coordinates for each interrogation volume.
        grid = constructGrid( size(vol1), IA_mp, step_mp, mp )

        for g in grid
            #  g[1],  g[2] are start and stop rows of interrogation volume
            #  g[3],  g[4] are start and stop columns of interrogation volume
            #  g[5],  g[6] are start and stop depths of interrogation volume
            #  g[7],  g[8] are start and stop rows inside UVW arrays
            #  g[9], g[10] are start and stop columns inside UVW arrays
            # g[11], g[12] are start and stop depths insie UVW arrays

            interr = view( vol1, g[1]:g[2], g[3]:g[4], g[5]:g[6] )

            if threshold > 0 && filtFun( interr ) < threshold
                continue;
            end

            # 1-. Shifting the searchArea ranges by previously calculated displacements
            offH = round(Int64, U[ g[7], g[9], g[11] ]);
            offW = round(Int64, V[ g[7], g[9], g[11] ]);
            offD = round(Int64, W[ g[7], g[9], g[11] ]);

            # 1.2-. Search volume coordinates after shifting.
            scoords, so = searchGrid( g, SM, offH, offW, offD, h, w, d )

            setTo0!( search )
            putWithinSearch!( search, vol2, so, scoords )

            setTo0!( sumS  );
            setTo0!( sumS2 );
            znccSumArrays!( search, sumS, sumS2 );

            meanI = Statistics.mean( interr );
            meanS = Statistics.mean( search );

            setTo0!( pads )
            setTo0!( padi )
            putWithinPadded!( padi, interr, meanI );
            putWithinPadded!( pads, search, meanS );

            meanS2 = meanS*meanS;
            meanI2 = meanI*meanI;
            stdI   = Statistics.std( interr, mean=meanI );

            # 2-. Cross-Correlation
            crossCorrelation!( ZNCC(), cmatrix, num, pads, padi, plan, iplan,
                               sumS, sumS2, meanS, meanS2, stdI, ssize, IA_mp, shifts );

            # 3-. Calculation of displacement
            ( r, c, z ) = approxTranslation( cmatrix, peak, Val(last_mp) )
			if ( !ignoreSN && last_mp )
            	SN[ g[7], g[9], g[11] ] = SNRatio( cmatrix, sigNoise, width )
			end

            # 4-. Updating U, V, W matrices
            t1 =  g[7]:g[8];
            t2 =  g[9]:g[10];
            t3 = g[11]:g[12];
            U[ t1, t2, t3 ] .= offH - r;
            V[ t1, t2, t3 ] .= offW - c;
            W[ t1, t2, t3 ] .= offD - z;
        end
    end # Multi-pass loop

    return U, V, W, SN
end

""" ZNCC out of place CROSS CORRELATION """
function crossCorrelation( ::ZNCC, f::A{T,3}, g::A{T,3}; typ=Float32 ) where {T<:Real, U<:Real}

    sizef  = size( f );
    sizeg  = size( g );

    sumf   = zeros( typ, sizef .+ 1 );
    sumf2  = zeros( typ, sizef .+ 1 );
    znccSumArrays!( f, sumf, sumf2 );

    csize = sizef .+ sizeg;
    corr  = zeros( typ, csize );
    num   = zeros( typ, csize );
    padf  = zeros( Complex{typ}, csize );
    padg  = zeros( Complex{typ}, csize );
    plan  =  plan_fft!( padf );
    iplan = plan_ifft!( padf );
	shifts = div.( csize, 2 ) .- 1;

    meanf = Statistics.mean( f );
    meang = Statistics.mean( g );

    putWithinPadded!( padf, f, meanf );
    putWithinPadded!( padg, g, meang );

    meanf2 = meanf^2;
    meang2 = meang^2;
    stdg   = Statistics.std( g, mean=meang );

    crossCorrelation!( ZNCC(), corr, num, padf, padg, plan, iplan,
                       sumf, sumf2, meanf, meanf2, stdg, sizef, sizeg, shifts );

    GC.gc();
    return corr;
end

function crossCorrelation( ::ZNCC, f::A{T,2}, g::A{T,2}; typ=Float32 ) where {T<:Real, U<:Real}

    sizef  = size( f );
    sizeg  = size( g );

    sumf   = zeros( typ, sizef .+ 1 );
    sumf2  = zeros( typ, sizef .+ 1 );
    znccSumArrays!( f, sumf, sumf2, (0,0) );

    csize = sizef .+ sizeg;
    corr  = zeros( typ, csize );
    num   = zeros( typ, csize );
    padf  = zeros( Complex{typ}, csize );
    padg  = zeros( Complex{typ}, csize );
    plan  =  plan_fft!( padf );
    iplan = plan_ifft!( padf );
	shifts = div.( csize[1:2], 1 ) .- 1;

    meanf = Statistics.mean( f );
    meang = Statistics.mean( g );

    putWithinPadded!( padf, f, meanf );
    putWithinPadded!( padg, g, meang );

    meanf2 = meanf^2;
    meang2 = meang^2;
    stdg   = Statistics.std( g, mean=meang );

    crossCorrelation!( ZNCC(), corr, num, padf, padg, plan, iplan,
                       sumf, sumf2, meanf, meanf2, stdg, sizef, sizeg, shifts);

    GC.gc();
    return corr;
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule SYMPHONY_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("SYMPHONY")
JLLWrappers.@generate_main_file("SYMPHONY", UUID("fad677ff-a6e1-547e-a877-6c265eeefd90"))
end  # module SYMPHONY_jll


text:
@testset "kernelproduct" begin
    rng = MersenneTwister(123456)
    x = rand(rng)*2
    v1 = rand(rng, 3)
    v2 = rand(rng, 3)

    k1 = LinearKernel()
    k2 = SqExponentialKernel()
    k3 = RationalQuadraticKernel()
    X = rand(rng, 2,2)

    k = KernelProduct(k1,k2)
    ks1 = 2.0*k1
    ks2 = 0.5*k2
    @test length(k) == 2
    @test string(k) == (
        "Product of 2 kernels:\n\tLinear Kernel (c = 0.0)\n\tSquared " *
        "Exponential Kernel"
    )
    @test k(v1, v2) == (k1 * k2)(v1, v2)
    @test (k * k3)(v1,v2) ≈ (k3 * k)(v1, v2)
    @test (k1 * k2)(v1, v2) == KernelProduct(k1, k2)(v1, v2)
    @test (k * ks1)(v1, v2) ≈ (ks1 * k)(v1, v2)
    @test (k * k)(v1, v2) == KernelProduct([k1, k2, k1, k2])(v1, v2)
    @test KernelProduct([k1, k2]) == KernelProduct((k1, k2)) == k1 * k2

    @test (KernelProduct([k1, k2]) * KernelProduct([k2, k1])).kernels == [k1, k2, k2, k1]
    @test (KernelProduct([k1, k2]) * k3).kernels == [k1, k2, k3]
    @test (k3 * KernelProduct([k1, k2])).kernels == [k3, k1, k2]

    @test (KernelProduct((k1, k2)) * KernelProduct((k2, k1))).kernels == (k1, k2, k2, k1)
    @test (KernelProduct((k1, k2)) * k3).kernels == (k1, k2, k3)
    @test (k3 * KernelProduct((k1, k2))).kernels == (k3, k1, k2)

    # Standardised tests.
    TestUtils.test_interface(k, Float64)
    test_ADs(
        x->SqExponentialKernel() * LinearKernel(c= x[1]), rand(1);
        ADs = [:ForwardDiff, :ReverseDiff, :Zygote],
    )

    test_params(k1 * k2, (k1, k2))
end


text:
using Test, MyAwesomePackage

out = plusTwo(3)

@test out == 5

max_stars_repo_name: sylvaticus/MyAwesomePackage

max_issues_repo_name: sylvaticus/MyAwesomePackage

max_forks_repo_name: sylvaticus/MyAwesomePackage


text:
using SpikingNeuralNetworks
SNN = SpikingNeuralNetworks
SNN.@load_units
using SpikeNetOpt
SNO = SpikeNetOpt

##
# smoke tests code mostly borrowed from
# https://github.com/AStupidBear/SpikingNeuralNetworks.jl/blob/master/test/iz_neuron.jl
##
using Plots
RS = SNN.IZ(; N = 1, param = SNN.IZParameter(; a = 0.02, b = 0.2, c = -65, d = 8))
IB = SNN.IZ(; N = 1, param = SNN.IZParameter(; a = 0.02, b = 0.2, c = -55, d = 4))
CH = SNN.IZ(; N = 1, param = SNN.IZParameter(; a = 0.02, b = 0.2, c = -50, d = 2))
FS = SNN.IZ(; N = 1, param = SNN.IZParameter(; a = 0.1, b = 0.2, c = -65, d = 2))
TC1 = SNN.IZ(; N = 1, param = SNN.IZParameter(; a = 0.02, b = 0.25, c = -65, d = 0.05))
TC2 = SNN.IZ(; N = 1, param = SNN.IZParameter(; a = 0.02, b = 0.25, c = -65, d = 0.05))
RZ = SNN.IZ(; N = 1, param = SNN.IZParameter(; a = 0.1, b = 0.26, c = -65, d = 2))
LTS = SNN.IZ(; N = 1, param = SNN.IZParameter(; a = 0.1, b = 0.25, c = -65, d = 2))
P = [RS, IB, CH, FS, TC1, TC2, RZ, LTS]

SNN.monitor(P, [:v])
T = 2second
for t = 0:T
    for p in [RS, IB, CH, FS, LTS]
        p.I = [10]
    end
    TC1.I = [(t < 0.2T) ? 0mV : 2mV]
    TC2.I = [(t < 0.2T) ? -30mV : 0mV]
    RZ.I = [(0.5T < t < 0.6T) ? 10mV : 0mV]
    SNN.sim!(P, [], 0.1ms)
end
for p in P
    v = SNN.vecplot(p, :v)
    v |> display
end
cell_type = "IZHI"
ngt_spikes = 10
current_ = SNO.current_search(cell_type, P[1], ngt_spikes)
P[1].I = [current_ * nA]
@show(current_)
SNN.monitor(P[1], [:v])
ALLEN_DURATION = 2000 * ms
ALLEN_DELAY = 1000 * ms

SNN.sim!(
    [P[1]];
    dt = 1 * ms,
    delay = ALLEN_DELAY,
    stimulus_duration = ALLEN_DURATION,
    simulation_duration = ALLEN_DURATION + ALLEN_DELAY + 443ms,
)

v = SNN.vecplot(P[1], :v)
v |> display
cell_type = "ADEXP"
ngt_spikes = 10
current_ = SNO.current_search(cell_type, P[1], ngt_spikes)
P[1].I = [current_ * nA]
@show(current_)
SNN.monitor(P[1], [:v])
ALLEN_DURATION = 2000 * ms
ALLEN_DELAY = 1000 * ms

SNN.sim!(
    [P[1]];
    dt = 1 * ms,
    delay = ALLEN_DELAY,
    stimulus_duration = ALLEN_DURATION,
    simulation_duration = ALLEN_DURATION + ALLEN_DELAY + 443ms,
)

v = SNN.vecplot(P[1], :v)
v |> display


text:
using JSON
using PipeLayout
import PipeLayout: ploss_coeff_nice, randdemand, select_covering_subset

# set random seed for reproducible results
srand(23)

# create topologies for ground structures
#  [x] square meshes of increasing sizes
#  [x] other meshes (with crossing diagonals? cf. Jakob's thesis)
#  [ ] random (or sobolev) points with voronoi triangulation
#  [ ] grids with holes
ground_structures = Dict{String, Topology}()
for args in [(6, 7, 18), (9, 9, 12), (13, 13, 9), (19, 19, 6)]
    m, n, w = args
    for deg in [4, 8]
        key = @sprintf "m%02d_n%02d_d%02d_l%02d" n m deg w
        ground_structures[key] = gridtopology(RegGrid{m,n,deg}(w))
    end
end

# for each topology: select subsets of nodes for terminals
#  [x] randomly
#  [ ] greedy farthest neighbor
terminals = Dict{String, Vector{Node}}()
for (gskey, topo) in ground_structures
    for nterminals in [3, 5, 7]
        key = @sprintf "%s.t%02d" gskey nterminals
        terminals[key] = select_covering_subset(topo.nodes, nterminals)
    end
end

# for each set of terminals (and related topo), create demand situations
#  [x] random distribution (source vs sink)
#  [x] same, but single source
#  [ ] hard distribution: maximize transport momentum
#  [ ] for each normalized demand: several flow scalings
demands = Dict{String, Vector{Float64}}()
for (tkey, nodes) in terminals
    nnodes = length(nodes)
    nodeidx = collect(1:nnodes)

    # collect assignments to source/sink role
    sourcesinks = []

    # single source, at first position
    push!(sourcesinks, (nodeidx[1:1], nodeidx[2:end]))

    # half and half (less sources)
    if nnodes > 3
        half = div(nnodes, 2)
        push!(sourcesinks, (nodeidx[1:half], nodeidx[half+1:end]))
    end

    for (sources, sinks) in sourcesinks
        for variant in ["a", "b"]
            dist = randdemand(nnodes, sources, sinks, 1.0)
            for flow in [25, 50, 100, 200]
                key = @sprintf "%s.%02d_%02d_%s_%03d" tkey length(sources) length(sinks) variant flow
                demands[key] = flow * dist
            end
        end
    end
end

# uniform pressure bounds
#  [x] one for each boundary node
pressure_bounds = Bounds(40.0, 80.0)

# diameter data
#  [x] full data and smaller subset
_diams = [0.4,  0.6,  0.8, 1.0,  1.2]
_costs = [0.68, 0.91, 1.2, 1.55, 1.96]
diameters = Dict{String, Vector{Diameter}}()
diameters["full"] = [Diameter(d,c) for (d,c) in zip(_diams, _costs)]
diameters["some"] = diameters["full"][1:2:end]

# create instances by combination
#  [x] come up with uniq naming
#  [x] write instances to files (json)
instances = String[]
for (key, demand) in demands
    gs, trm, dmd = split(key, ".")
    topology = ground_structures[gs]
    nodes = terminals["$gs.$trm"]
    pressure = [pressure_bounds for n in nodes]

    for (diamkey, diams) in diameters
        instkey = "$key.$diamkey"
        instance = Instance(nodes, demand, pressure, diams, ploss_coeff_nice)
        push!(instances, instkey)

        open("$instkey.instance.json", "w") do f
            JSON.print(f, instance, 2)
        end
        open("$instkey.topology.json", "w") do f
            JSON.print(f, topology, 2)
        end
    end
end

sort!(instances)
open("instances.txt", "w") do f
    write(f, join(instances, "\n"))
end


text:
### A Pluto.jl notebook ###
# v0.16.4

using Markdown
using InteractiveUtils

# ╔═╡ 029fefe6-56bf-4e5c-9198-7a2bb8a00889
begin
	using Parameters
	using UUIDs
	using PlutoUI
end

# ╔═╡ 9827859f-ffd9-4c55-abd8-7969dcaa90c9
number_of_iterations = 100

# ╔═╡ 317bfe52-fc27-4c2c-a58f-a7b55ce1973f
number_of_agents = 10

# ╔═╡ 4ed1ab97-435d-4928-8a36-f93916a89c28
number_of_assets = 10

# ╔═╡ 07de0609-6f9b-406d-9fbf-cdfeaacdd335

# ╔═╡ 8f4a5b5e-1319-4df1-8fa4-b8a1e3b2d327
uuid4()

# ╔═╡ 886205dc-3f0f-4b5c-928e-4030b3cfd812
A = rand(100,2)

# ╔═╡ 6f3c3296-69f0-441e-8d01-814dd996f143
size(A,2)

# ╔═╡ b71b93a7-8106-4207-917c-4c7f907de970
rand(1:number_of_assets,number_of_assets)

# ╔═╡ 59cea89b-d09f-46be-ac26-5ce6befebf7f
d = Dict("A"=>rand(10),"B"=>2)

# ╔═╡ 9f7c16e9-6cf4-4d95-9d6c-d0fd9e0ffd91
typeof(d)

# ╔═╡ 0f860b70-7f65-4d04-a27d-088e752c496f
function ingredients(path::String; module_name::Symbol = :model)

	# this is from the Julia source code (evalfile in base/loading.jl)
	# but with the modification that it returns the module instead of the last object

	# JV change ... want to
	# name = Symbol(basename(path))
	name = module_name

	# original -
	m = Module(name)
	Core.eval(m,
        Expr(:toplevel,
             :(eval(x) = $(Expr(:core, :eval))($name, x)),
             :(include(x) = $(Expr(:top, :include))($name, x)),
             :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
             :(include($path))))
	m
end

# ╔═╡ df0780f6-2854-11ec-21a7-8b1cd061e361
begin

	# load my code into module M -
	M = ingredients("../Include.jl"; module_name = :ABM)

	# alias my types that are declared -
	VLAgentModel = M.VLAgentModel
	VLGameWorld = M.VLGameWorld
	VLOrderModel = M.VLOrderModel
	ActionType = M.ActionType
	OrderType = M.OrderType

	# return -
	nothing
end

# ╔═╡ 914359d9-bfeb-4ab3-8448-041a3f7444f7
function my_technical_agent_trade_logic(agent_model, asset_price_array, iteration_index)

	with_terminal() do
		println("$(asset_price_array)")
	end

	# initialize -
	order_model_array = Array{VLOrderModel,1}()
	agent_position_array = agent_model.agent_position_array

	# get the parameters from this agent -
	parameters_dict = agent_model.agent_parameters

	# get the c-parameter array -
	c_parameter_array = parameters_dict["c_parameter_array"]
	price_lag_array = parameters_dict["price_lag_array"]

	# compute the setpoint array -
	(number_of_iterations, number_of_assets) = size(asset_price_array)
	for asset_index = 1:number_of_assets

		# what is the lag?
		asset_lag_index = price_lag_array[asset_index]

		# what is the lagged price for this asset?
		price_setpoint_value = 1.0
		if (asset_lag_index<iteration_index)
			price_setpoint_value = asset_price_array[(iteration_index - asset_lag_index),asset_index]
		end

		# build the order model -
		c_value = c_parameter_array[asset_index]
		current_price = asset_price_array[asset_index]
		desired_position_size = c_value*sign(current_price - price_setpoint_value)
		Δω = agent_position_array[iteration_index,asset_index] - desired_position_size

		# build the order -
		# an order has this data -
		local_asset_index = asset_index
    	action = (Δω < 0) ? -1 : 1
    	order_type = 1
    	quantity = 1
    	price = current_price
    	agent_id = agent_model.agent_id
		order_model = VLOrderModel(local_asset_index,action,order_type,quantity,price,agent_id);
		push!(order_model_array,order_model)
	end

	# return -
	return order_model_array
end

# ╔═╡ 76280c54-afd3-4655-a991-54b957de8b5f
begin

	# build the game world -
	game_world = M.build(VLGameWorld; number_of_iterations=number_of_iterations, number_of_assets=number_of_assets);

	# build agent array -
	my_agent_array = Array{VLAgentModel,1}()
	for agent_index = 1:number_of_agents

		# let's use the defaults for the various agent settings -
		agent_model = M.build(VLAgentModel; number_of_iterations=number_of_iterations, number_of_assets=number_of_assets)

		# setup parameters for this agent -
		c_array = rand(1:number_of_assets,number_of_assets)
		lag_array = rand(1:10,number_of_assets)
		p = Dict("c_parameter_array"=>c_array,"price_lag_array"=>lag_array)
		agent_model.agent_parameters = p

		# let's have all my agents be technical traders -
		agent_model.agent_trade_logic = my_technical_agent_trade_logic

		# grab -
		push!(my_agent_array, agent_model)
	end

	# attach the agents to the world -
	game_world.agent_array = my_agent_array

	# return -
	nothing
end

# ╔═╡ 92dbc919-0f2a-4abd-8cac-191ca51c2e26
R = M.run(game_world)

# ╔═╡ 0ee60583-af5a-48dd-a5e0-7fd1311bafc7
instances(ActionType)

# ╔═╡ ce1ba111-1cba-4c4a-945e-104f3b8acee0
html"""
<style>
main {
    max-width: 1200px;
    margin: auto;
	width: 90%;
    font-family: "Roboto, monospace";
}

a {
    color: blue;
    text-decoration: none;
}

.H1 {
    padding: 0px 30px;
}
</style>
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[compat]
Parameters = "~0.12.3"
PlutoUI = "~0.7.16"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[HypertextLiteral]]
git-tree-sha1 = "f6532909bf3d40b308a0f360b6a0e626c0e263a8"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.1"

[[IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[Parameters]]
deps = ["OrderedCollections", "UnPack"]
git-tree-sha1 = "34c0e9ad262e5f7fc75b10a9952ca7692cfc5fbe"
uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
version = "0.12.3"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "f19e978f81eca5fd7620650d7dbea58f825802ee"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.1.0"

[[PlutoUI]]
deps = ["Base64", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "4c8a7d080daca18545c56f1cac28710c362478f3"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.16"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[UnPack]]
git-tree-sha1 = "387c1f73762231e86e0c9c5443ce3b4a0a9a0c2b"
uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
version = "1.0.2"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
"""

# ╔═╡ Cell order:
# ╠═029fefe6-56bf-4e5c-9198-7a2bb8a00889
# ╠═9827859f-ffd9-4c55-abd8-7969dcaa90c9
# ╠═317bfe52-fc27-4c2c-a58f-a7b55ce1973f
# ╠═4ed1ab97-435d-4928-8a36-f93916a89c28
# ╠═76280c54-afd3-4655-a991-54b957de8b5f
# ╠═df0780f6-2854-11ec-21a7-8b1cd061e361
# ╠═92dbc919-0f2a-4abd-8cac-191ca51c2e26
# ╠═07de0609-6f9b-406d-9fbf-cdfeaacdd335
# ╠═914359d9-bfeb-4ab3-8448-041a3f7444f7
# ╠═0ee60583-af5a-48dd-a5e0-7fd1311bafc7
# ╠═8f4a5b5e-1319-4df1-8fa4-b8a1e3b2d327
# ╠═886205dc-3f0f-4b5c-928e-4030b3cfd812
# ╠═6f3c3296-69f0-441e-8d01-814dd996f143
# ╠═b71b93a7-8106-4207-917c-4c7f907de970
# ╠═59cea89b-d09f-46be-ac26-5ce6befebf7f
# ╠═9f7c16e9-6cf4-4d95-9d6c-d0fd9e0ffd91
# ╟─0f860b70-7f65-4d04-a27d-088e752c496f
# ╟─ce1ba111-1cba-4c4a-945e-104f3b8acee0
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using Luna

a = 13e-6
gas = :Ar
pres = 5
flength = 15e-2

τfwhm = 30e-15
λ0 = 800e-9
energy=1e-6

modes = (
    Capillary.MarcatiliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=0.0, loss=false),
    Capillary.MarcatiliMode(a, gas, pres, n=1, m=2, kind=:HE, ϕ=0.0, loss=false),
)
nmodes = length(modes)

grid = Grid.RealGrid(flength, λ0, (160e-9, 3000e-9), 1e-12)

energyfun, energyfunω = Fields.energyfuncs(grid)
normfun = NonlinearRHS.norm_modal(grid.ω)

densityfun = let dens0=PhysData.density(gas, pres)
    z -> dens0
end

ionpot = PhysData.ionisation_potential(gas)
ionrate = Ionisation.ionrate_fun!_ADK(ionpot)

responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),
             Nonlinear.PlasmaCumtrapz(grid.to, grid.to, ionrate, ionpot))

inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy)

Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs, modes,
                              :y; full=true)

statsfun = Stats.default(grid, Eω, modes, linop, transform; gas=gas, windows=((150e-9, 300e-9),))
output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
linop = LinearOps.make_const_linop(grid, modes, λ0)

Luna.run(Eω, grid, linop, transform, FT, output)

Plotting.pygui(true)
Plotting.stats(output)
Plotting.prop_2D(output)


text:
struct WangCombinedStrategy <: SVDD.InitializationStrategyCombined
    solver
    gamma_search_range
    C
    scoring_function
end

fnr(x::ConfusionMatrix) = 1-tpr(x)

function SVDD.calculate_gamma(model, strategy::WangCombinedStrategy)
    m = deepcopy(model)
    data_target, data_outliers = SVDD.generate_binary_data_for_tuning(m.data)
    artificial_data = hcat(data_target, data_outliers)
    ground_truth = vcat(fill(:inlier, size(data_target, 2)),
                    fill(:outlier, size(data_outliers, 2)))

    @info("[Gamma Search] Searching for parameters gamma and C.")
    best_gamma = 1.0
    best_C = 1.0
    best_score = Inf
    for gamma in strategy.gamma_search_range

        init_strategy = FixedParameterInitialization(GaussianKernel(gamma), 1.0)
        initialize!(m, init_strategy)
        set_adjust_K!(m, true)
        for C in strategy.C
            try
                @info("[Gamma Search] Testing C = $C, gamma = $gamma")
                set_param!(m, Dict(:C => C))
                SVDD.fit!(m, strategy.solver)
            catch e
                @info("[Gamma Search] Fitting failed for C = $C, gamma $gamma.")
                println(e)
                continue
            end
            cm = ConfusionMatrix(SVDD.classify.(SVDD.predict(m, artificial_data)), ground_truth)

            score = 0.5 * fpr(cm) + 0.5 * fnr(cm)
            if best_score > score
                @info("[Gamma Search] New best found with gamma = $gamma and score = $score.")
                best_gamma = gamma
                best_C = C
                best_score = score
            end
        end
    end
    @info("Best paremeters: C=$best_C, gamma=$best_gamma")
    return (best_C, MLKernels.GaussianKernel(best_gamma))
end

function SVDD.get_parameters(model, strategy::WangCombinedStrategy)
    SVDD.calculate_gamma(model, strategy)
end

max_stars_repo_name: holtri/lama

max_issues_repo_name: holtri/lama

max_forks_repo_name: holtri/lama


text:
# calculate standard deviation of sensor data

module CalcStdDev
    using DataFrames, CSV, Statistics

    function main()
        data_path = joinpath(split(@__FILE__, "src")[1], "data/sensor_data_200.txt")
        df_200_mm = CSV.read(data_path, DataFrame,
                             header=["date", "time", "ir", "lidar"],
                             delim=' ')

        # calculate from definition
        zs = df_200_mm.lidar # observation
        mean_def = sum(zs) / length(zs) # mean
        diff_square = [(z - mean_def)^2 for z in zs]

        sampling_var = sum(diff_square) / length(zs) # sampling variance
        unbiased_var = sum(diff_square) / (length(zs) - 1) # unbiased variance

        sampling_std_dev = sqrt(sampling_var)
        unbiased_std_dev = sqrt(unbiased_var)

        # calculate by Statistics
        stats_std_dev = Statistics.std(zs, corrected=false)

        println("Sampling standard deviation = $sampling_std_dev")
        println("Unbiased standard deviation = $unbiased_std_dev")
        println("Standard deviation by Statistics = $stats_std_dev")
    end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
using MusicXML
using Test

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


text:
#Reference: https://jump.dev/JuMP.jl/0.17/refexpr.html
using CSV, DelimitedFiles, DataFrames
using JuMP, CPLEX

#Preparing optimization model
model = Model(CPLEX.Optimizer)
set_optimizer_attribute(model, "CPX_PARAM_EPINT", 1e-5)

const N_GEN = 145
const N_HRS = 24

dfGenerator = CSV.read(".//inputs//data_generators.csv", DataFrame)
dfDemand = CSV.read(".//inputs//demand_reserves.csv", DataFrame)

@variable(model, genOnOff[1:N_GEN, 1:N_HRS], Bin)
@variable(model, genOut[1:N_GEN, 1:N_HRS] >= 0)
@variable(model, genStartUp[1:N_GEN, 1:N_HRS], Bin)
@variable(model, genShutDown[1:N_GEN, 1:N_HRS], Bin)

#Setting the objective function
@objective(model, Min, sum(sum(dfGenerator.FixedCost[g]*genOnOff[g,t] + dfGenerator.VariableCost[g]*genOut[g,t] + dfGenerator.StartUpCost[g]*genStartUp[g,t] + dfGenerator.ShutdownCost[g]*genShutDown[g, t] for g in 1:N_GEN) for t in 1:N_HRS) )

#Initial constraints
#Generators should have an initial and defined value t=0
@constraint(model, constInit[g=1:N_GEN], genOnOff[g,1] == dfGenerator.Uini[g])

#Initial values of Generator shutdown and startup values should be defined
@constraint(model, constShutDn[g=1:N_GEN], genShutDown[g,1] == 0)
@constraint(model, constShutUp[g=1:N_GEN], genStartUp[g,1] == 0)

#Demand constraint
#for t in 1:N_HRS
@constraint(model, constDemand[t=1:N_HRS], sum(genOut[g, t] for g=1:N_GEN) == dfDemand.Demand[t])
#end

#Power bounds
@constraint(model, constMaxPow[t=1:N_HRS, g=1:N_GEN],  genOut[g,t] <= dfGenerator.MaxPowerOut[g]*genOnOff[g,t] )
@constraint(model, constMinPow[t=1:N_HRS, g=1:N_GEN],  genOut[g,t] >= dfGenerator.MinPowerOut[g]*genOnOff[g,t] )

#Logical Conditions
@constraint(model, constLog1[t=2:N_HRS, g=1:N_GEN], (genStartUp[g,t] - genShutDown[g,t]) == (genOnOff[g,t] - genOnOff[g,t-1]))
@constraint(model, constLog2[t=1:N_HRS, g=1:N_GEN], genStartUp[g,t] + genShutDown[g, t] <= 1)

#Ramp limits
@constraint(model, constRampUp[t=2:N_HRS, g=1:N_GEN], genOut[g,t] - genOut[g,t-1] <= dfGenerator.RampUpLimit[g]*genOnOff[g, t-1] + dfGenerator.RampStartUpLimit[g]*genStartUp[g,t])
@constraint(model, constRampDown[t=2:N_HRS, g=1:N_GEN], (genOut[g,t-1] - genOut[g,t]) <= (dfGenerator.RampDownLimit[g]*genOnOff[g,t] + dfGenerator.RampShutDownLimit[g]*genShutDown[g,t]))

#Reserves
@constraint(model, constReserves[t=1:N_HRS], sum(dfGenerator.MaxPowerOut[g]*genOnOff[g,t] for g=1:N_GEN) >= (dfDemand.Demand[t] + dfDemand.Reserves[t]))

mystatus = optimize!(model)

println("------------------------------------")
println("------- TERMINATION STATUS -------")
if termination_status(model) == MOI.OPTIMAL
    optimal_solution = value.(genOnOff)
    optimal_objective = objective_value(model)
    println("Optimal result")
elseif termination_status(model) == MOI.TIME_LIMIT && has_values(model)
    suboptimal_solution = value.(genOnOff)
    suboptimal_objective = objective_value(model)
    println("Time limit was exceeded")
else
    error("The model was not solved correctly.")
end

#Printing general results
println("------------------------------------")
println("------- OBJECTIVE VALUE -------")
println("Objective value: ", JuMP.objective_value(model))

println("------------------------------------")
println("------- PRIMAL STATUS -------")
println(primal_status(model))

println("------------------------------------")
println("------- DUAL STATUS -------")
println(JuMP.dual_status(model))

#Storing decision variables results
#Store whole model (TBD)
#write_to_file(model, ".//outputs//my_model.mps")
typeFileAccess = "w+"

open(".//outputs//resultsGenOut.csv", typeFileAccess) do io
           #writedlm(io, [x y], ',')
           writedlm(io,  JuMP.value.(genOut), ',')
    end;

open(".//outputs//resultsGenCommit.csv", typeFileAccess) do io
               #writedlm(io, [x y], ',')
               writedlm(io, JuMP.value.(genOnOff), ',')
         end;

open(".//outputs//resultsGenStartUp.csv", typeFileAccess) do io
              #writedlm(io, [x y], ',')
              writedlm(io, JuMP.value.(genStartUp), ',')
        end;

open(".//outputs//resultsGenShutDown.csv", typeFileAccess) do io
              #writedlm(io, [x y], ',')
              writedlm(io, JuMP.value.(genShutDown), ',')
        end;

println("------------------------------------")
println("------- OPTIMAL SOLUTION -------")
num_print_hrs = min(10, N_HRS)
num_print_gen = min(20, N_GEN)

for t=1:num_print_hrs
  print("Hour ", t-1, "\t")
end
println("")
for g=1:num_print_gen
  for t=1:num_print_hrs
    print(round(JuMP.value(genOut[g,t]); digits=2 ),"\t")
  #  println(round(pi; digits = 3))

  end
  println("\t Generator: $g")
end

println("------------------------------------")
println("------- GENERATORS COMMITED --------")
for t=1:num_print_hrs
  print("Hour ", t-1, "\t")
end

for g=1:num_print_gen
  for t=1:num_print_hrs
    #print(JuMP.value(genOnOff[g,t]),"\t")
    print(convert(Int32, JuMP.value(genOnOff[g,t])),"\t")
  end
  println("\t Generator: $g")
end

println("-------------------------------")
println("----------LAST DAY GEN OFF/ON---------")
print(JuMP.value.(genOnOff[:,end]))

println("-------------------------------")
#Clean memory
GC.gc()

max_stars_repo_name: power-systems-op/grace-arpae

max_issues_repo_name: power-systems-op/grace-arpae

max_forks_repo_name: power-systems-op/grace-arpae


text:
using Distributions, Plots, LaTeXStrings; pyplot()

cUnif = Uniform(0,2π)
xGrid, N = 0:0.1:2π, 10^6

stephist( rand(N)*2π, bins=xGrid,
	normed=:true, c=:blue,
	label="MC Estimate")
plot!( xGrid, pdf.(cUnif,xGrid),
	c=:red,ylims=(0,0.2),label="PDF", ylabel="Density",xticks=([0:π/2:2π;],
	["0", L"\dfrac{\pi}{2}", L"\pi", L"\dfrac{3\pi}{2}", L"2\pi"]))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
# Содержание
# Функция вычисления правой части прямой задачи     :   |directRP|
# Функция вычисления трехдиагонального
#                       якобиана правой части       :   |DRP_y|
# Обертка `DRP_y`                                   :   |∂DRP_∂y|
# Функция вычисления якобиана правой части          :   |∂directRP_∂y|
# Функция нахождения решения прямой задачи          :   |solve|

@doc raw"""
    directRP(y::Vector, m::Int,
             Xₙ::Vector, N::Int,
             ε::Real,
             ulₘ::Vector, urₘ::Vector,
             qₙ::Vector)

Вычисляет вектор правой части с помощью конечно-разностной
аппроксимации пространственных производных.

# Arguments
- `y::Vector`:          Вектор размера `N-1`, решение в текущий момент без ГТ.
- `m::Int`:             Номер шага в сетке по времени.
- `Xₙ::Vector`:         Вектор размера `N+1` сетки по ``X``, вместе с ГТ.
- `N::Int`:             Число **интервалов** в полной сетке по ``X``.
- `ε::Real`:            Малый параметр при старшей производной.
- `ulₘ::Vector`:        Вектор сеточных значений левого ГУ.
- `urₘ::Vector`:        Вектор сеточных значений правого ГУ.
- `qₙ::Vector`:         Вектор размера `N-1` сеточные значений  неоднородности без ГТ.

!!! warn
    Функция не входит в публичный API, поэтому размерность
    входных векторов отличается от аналогичных в [`solve`](@ref).
     - `qₙ`, `y` размера `N-1`!
     - `Xₙ` размера `N+1`.

"""
function directRP(y::Vector, m::Int,
                  Xₙ::Vector, N::Int,
                  ε::Real,
                  ulₘ::Vector, urₘ::Vector,
                  qₙ::Vector)

    @assert length(Xₙ) == N+1
    @assert length(qₙ) == N-1
    @assert length(y)  == N-1
    @assert m < length(ulₘ)

    RP = zero(y)                    # Создаем нулевой вектор того же типа и размера

    # Здесь нужно применить сдвиг индексов в `Xₙ` на +1!
    # Xₙ[1] соответствует первому узлу нашей сетки — ``x_0``
    # Xₙ[2] — ``x_1``
    # Xₙ[3] — ``x_2``
    #
    # Для удобства, формулы записаны в следующем виде
    # d[1] = ... Xₙ[1+1+1] ...,
    # первая единица соответствует единице в d[1],
    # вторая — берем следующий элемент,
    # третья — совершаем сдвиг индексов для Xₙ
    RP[1]     = 2ε / (Xₙ[1+1+1] - Xₙ[1-1+1]) * ( (y[2] - y[1]) /  (Xₙ[1+1+1] - Xₙ[1+1]) - (y[1] - ulₘ[m])/(Xₙ[1+1] - Xₙ[1-1+1])) +
    + y[1]   * (y[2]   - ulₘ[m])/(Xₙ[1+1+1] - Xₙ[1-1+1]) - y[1] * qₙ[1]

    # Xₙ[n  ] — соответствует ``x_{n-1}``
    # Xₙ[n+1] — соответствует ``x_{n  }``
    # Xₙ[n+2] — соответствует ``x_{n+1}``
    #
    # Для удобства, формулы записаны в следующем виде
    # d[n] = ... Xₙ[n-1+1],
    # первое слагаемое соответствует d[n],
    # второе — берем предыдущий элемент,
    # третье — совершаем сдвиг индексов для Xₙ
    for n in 2:N-2
        RP[n] = 2ε / (Xₙ[n+1+1] - Xₙ[n-1+1]) * ( (y[n+1] - y[n]) /  (Xₙ[n+1+1] - Xₙ[n+1]) - (y[n] - y[n-1])/(Xₙ[n+1] - Xₙ[n-1+1])) +
        + y[n]   * (y[n+1] - y[n-1])/(Xₙ[n+1+1] - Xₙ[n-1+1]) - y[n]   * qₙ[n]
    end

    RP[N-1]   =  2ε / (Xₙ[N+1] - Xₙ[N-2+1]) * ( (urₘ[m] - y[N-1]) /  (Xₙ[N+1] - Xₙ[N-1+1]) - (y[N-1] - y[N-2])/(Xₙ[N-1+1] - Xₙ[N-2+1])) +
    + y[N-1] * (urₘ[m] - y[N-2])/(Xₙ[N+1] - Xₙ[N-2+1]) - y[N-1] * qₙ[N-1]

    return RP
end

@doc raw""""
    DRP(y::Int, m::Int,
        Xₙ::Vector, N::Int,
        ε::Real,
        ulₘ::Vector, urₘ::Vector,
        qₙ::Vector) -> (::Vector, ::Vector, ::Vector)

Возращает три диагонали якоибана правой части прямой задачи,
чтобы после сформировать трехдигональную матрицу `Tridiagonal`.

!!! warn
    Функция не входит в публичный API, поэтому размерность
    входных векторов отличается от аналогичных в [`solve`](@ref).
     - `qₙ`, `y` размера `N-1`!
     - `Xₙ` размера `N+1`.

# Return
Возвращает три вектора `dl`, `d`, `dl` элементов якоибана.
Нижняя и верхняя диагонали длины ``N-2``, главная — ``N-1``.
```
    [ d[1]  du[1]                         ]
    [ dl[1] d[2]  du[2]                   ]
    [ 0     dl[2] d[3] du[3]              ]
    [           ...  ...  ...             ]
    [                ...  ...     du[N-2] ]
    [                     dl[N-2] d[N-1]  ]
```

# Example
```
julia> # Данные необохдимо подготовить, например, как в examples/example_direct.jl

julia> # Подготовим массивы, выбросив граничные точки

julia> # ведь тестируемая функция — для внутреннего использования

julia> qq = qₙ[2:N];

julia> y = y₀[2:N];

julia> dl, d, du = NonLinearReactionAdvectionDiffusionWithFrontData.DRP_y(y, 1, Xₙ, N, ε, ulₘ, urₘ, qq);

julia> Tridiagonal( dl, d, du )
49×49 Tridiagonal{Float64,Array{Float64,1}}:
 -991.038   305.462      ⋅         ⋅         ⋅         ⋅         ⋅     …       ⋅         ⋅         ⋅         ⋅         ⋅        ⋅
[...]
```
"""
function DRP_y(y::Vector, m::Int,
             Xₙ::Vector, N::Int,
             ε::Real,
             ulₘ::Vector, urₘ::Vector,
             qₙ::Vector)

    @assert length(Xₙ) == N+1
    @assert length(qₙ) == N-1
    @assert length(y) == N-1
    @assert m < length(ulₘ)

    dl = zeros(N-2);        # Поддиагональные элементы
    d = zeros(N-1);         #    Диагональные элементы
    du = zeros(N-2);        # Наддиагональные элементы

    # Здесь нужно применить сдвиг индексов на +1 для Xₙ!
    # Xₙ[1] соответствует ПЕРВОМУ узлу нашей сетки — ``x_0``
    # Xₙ[2] — ``x_1``
    # Xₙ[3] — ``x_2``
    #
    # Для удобства, формулы записаны в следующем виде
    # d[1] = ... Xₙ[1+1+1] ...,
    # первая единица соответствует единице в d[1],
    # вторая — берем следующий элемент,
    # третья — совершаем сдвиг индексов для Xₙ
    d[1] = 2ε / ( Xₙ[1+1+1] - Xₙ[1-1+1]) * ( -1 / (Xₙ[1+1+1] - Xₙ[1+1]) - 1 / (Xₙ[1+1] - Xₙ[1-1+1])) +
            + (y[2] - ulₘ[m]) / ( Xₙ[1+1+1] - Xₙ[1-1+1])  - qₙ[1]

    du[1] = 2ε / ( Xₙ[1+1+1] - Xₙ[1-1+1]) * ( 1/ (Xₙ[1+1+1] - Xₙ[1+1]) ) +
            + y[1] / (Xₙ[1+1+1] - Xₙ[1-1+1])

    # Xₙ[n  ] — соответствует ``x_{n-1}``
    # Xₙ[n+1] — соответствует ``x_{n  }``
    # Xₙ[n+2] — соответствует ``x_{n+1}``
    #
    # Для удобства, формулы записаны в следующем виде
    # d[n] = ... Xₙ[n-1+1],
    # первое слагаемое соответствует d[n],
    # второе — берем предыдущий элемент,
    # третье — совершаем сдвиг индексов для Xₙ
    for n in 2:N-2
        # Сдвиг индексов для поддиагональных элементов `dl`
        # см # Return документации функции
        dl[n - 1] = 2ε / ( Xₙ[n+1+1] - Xₙ[n-1+1]) * ( 1 / (Xₙ[n+1] - Xₙ[n-1+1])) -
                    y[n] / ( Xₙ[n+1+1] - Xₙ[n-1+1])

        d[n]      = 2ε / ( Xₙ[n+1+1] - Xₙ[n-1+1]) * ( -1 / (Xₙ[n+1+1] - Xₙ[n+1]) - 1 / (Xₙ[n+1] - Xₙ[n-1+1])) +
                + (y[n+1] - y[n-1]) / ( Xₙ[n+1+1] - Xₙ[n-1+1])  - qₙ[n]

        du[n]     = 2ε / ( Xₙ[n+1+1] - Xₙ[n-1+1]) * ( 1/ (Xₙ[n+1+1] - Xₙ[n+1]) ) +
                + y[n] / (Xₙ[n+1+1] - Xₙ[n-1+1])
    end

    dl[N-2] = 2ε / ( Xₙ[N+1] - Xₙ[N-2+1]) * ( 1 / (Xₙ[N-1+1] - Xₙ[N-2+1])) -
                 y[N-1] / ( Xₙ[N+1] - Xₙ[N-2+1])
    d[N-1] =  2ε / ( Xₙ[N+1] - Xₙ[N-2+1]) * ( -1 / (Xₙ[N+1] - Xₙ[N-1+1]) - 1 / (Xₙ[N-1+1] - Xₙ[N-2+1]))             +
                + (urₘ[m] - y[N-2]) / ( Xₙ[N+1] - Xₙ[N-2+1])  - qₙ[N-1]

    return dl, d, du
end

@doc raw"""
    ∂DRP_∂y(y::Vector, m::Int,
            Xₙ::Vector, N::Int,
            ε::Real,
            ulₘ::Vector, urₘ::Vector,
            qₙ::Vector) -> Tridiagonal

Обертка фнукции [`DRP_y`](@ref), которая возвращает `Tridiagonal( DRP_y(...))`
трехдиагональную матрицу из векторов, которые возвращает `DRP_y(...)`.

!!! warn
    Функция не входит в публичный API, поэтому размерность
    входных векторов отличается от аналогичных в [`solve`](@ref).
     - `qₙ`, `y` размера `N-1`!
     - `Xₙ` размера `N+1`.
"""
function ∂DRP_∂y(y::Vector, m::Int,
               Xₙ::Vector, N::Int,
               ε::Real,
               ulₘ::Vector, urₘ::Vector,
               qₙ::Vector)
    return Tridiagonal( (DRP_y(y, m, Xₙ, N, ε, ulₘ, urₘ, qₙ))... )
end

@doc raw"""
    ∂directRP_∂y(y::Vector, m::Int,
                  Xₙ::Vector, N::Int,
                  ε::Real,
                  ulₘ::Vector, urₘ::Vector,
                  qₙ::Vector)

Функция якобиана для `adjointRP`.
Размерности входных векторов такие же, как и у [`directRP`](@ref).

!!! warn
    Функция не входит в публичный API, поэтому размерность
    входных векторов отличается от аналогичных в [`solve`](@ref).
     - `qₙ`, `y` размера `N-1`!
     - `Xₙ` размера `N+1`.
"""
function ∂directRP_∂y(y::Vector, m::Int,
                  Xₙ::Vector, N::Int,
                  ε::Real,
                  ulₘ::Vector, urₘ::Vector,
                  qₙ::Vector)
    return ForwardDiff.jacobian( z -> directRP(z, m, Xₙ, N, ε, ulₘ, urₘ, qₙ), y)
end

@doc raw"""
    solve(y₀::Vector, Xₙ::Vector, N::Int,
          Tₘ::Vector, M::Int,
          ε::Real, ulₘ::Vector, urₘ::Vector,
          qₙ::Vector,
          RP::Function = directRP,
          jac::Function = ∂directRP_∂y;
          α::Complex = complex(0.5, 0.5),
          create_mesh::Function = x -> [NaN, NaN]) -> Matrix, Matrix, Vector

Функция, которая находит решение прямой задачи с помощью
одностадийной схемы Розенброка с комплексным коэффициентом.

# Arguments
- `y₀::Vector`:         Сеточные значения начального условия.
- `Xₙ::Vector`:         Пространственная сетка по ``X``.
- `N::Int`:             Число **интервалов** сетки по ``X``.
- `Tₘ::Vector`:         Пространственная сетка по ``t``.
- `M::Int`:             Число **интервалов** сетки по ``t``.
- `ε::Real`:            Малый параметр при старшей производной.
- `ulₘ::Function`:      Сеточные значения левого ГУ.
- `urₘ::Function`:      Сеточные значения правого ГУ.
- `qₙ::Vector`:         Сеточные значения "неоднородности", см. постановку задачи.
- `RP::Function`:       Функция вычисления вектора правой части — `directRP`.
- `jac::Function`:      Якобиан правой части по вектору `y` — `∂DRP_∂y`.
- `α::Complex`:         Коэффициент схемы. При `α = 0` — схема Эйлера, при `α = complex(0.5, 0.5)` — схема Розенброка с комплексным коэффициентом.
- `create_mesh::Function`: Принимает один аргумент `x_tp` положение переходного слоя и формирует подходяющую сетку.

# Return
Тройку:
    - Матрицу размера ``(N+1, M+1)``, содержащую значения искомой функции на сетках ``X\_n, T\_m``,.
    - Матрицу размера ``(N+1, M+1)``, содержащую значения стеки ``X\_n`` на каждом моменте времени.
    - Вектора размера ``M+1``, содержащий координаты положения переходного слоя в каждый момент времени.

!!! info
     - `Xₙ`, `qₙ`, `y₀`     векторы размера ``N+1``!
     - `Tₘ`, `ulₘ`, `urₘ`   векторы размера ``M+1``!
"""
function solve(y₀::Vector, Xₙ::Vector, N::Int,
               Tₘ::Vector, M::Int,
               ε::Real, ulₘ::Vector, urₘ::Vector,
               qₙ::Vector,
               RP::Function = directRP,
               jac::Function = ∂DRP_∂y;
               α::Complex = complex(0.5, 0.5),
               create_mesh::Function = x -> [NaN, NaN])

    if length(Xₙ) != N+1
        throw(ArgumentError("""length(Xₙ) == $(length(Xₙ)), N == $(N)
        Массив Xₙ должен иметь размерность N+1."""))
    end
    if length(qₙ) != N+1
        throw(ArgumentError("""length(qₙ) == $(length(qₙ)), N == $(N)
        Массив qₙ должен иметь размерность N+1."""))
    end
    if length(y₀) != N+1
        throw(ArgumentError("""length(y) == $(length(y)), N == $(N)
        Массив `y` должен иметь размерность N+1."""))
    end

    if length(Tₘ) != M+1
        throw(ArgumentError("""length(Tₘ) == $(length(Tₘ)), M == $(M)
        Массив `Tₘ` должен иметь размерность M+1"""))
    end
    if length(ulₘ) != M+1
        throw(ArgumentError("""length(ulₘ) == $(length(ulₘ)), M == $(M)
        Массив `ulₘ` должен иметь размерность M+1"""))
    end
    if length(urₘ) != M+1
        throw(ArgumentError("""length(y) == $(length(y)), M == $(M)
        Массив `ulₘ` должен иметь размерность M+1"""))
    end

    u = zeros(N+1, M+1);                # Решение на всей сетке ``Xₙ x Tₘ``
    u[1  , :]   = ulₘ;                  # Запись левого  ГУ на каждом шаге
    u[N+1, :]   = urₘ;                  # Запись правого ГУ на каждом шаге
    u[:, 1] = y₀;                       # Запись начального условия, это перепишет ГУ.
                                        #
    Xₜₚ = zeros(M+1);                   # Массив для хранения положения переходного слоя
                                        # Определяется апостриорным способом
    XX = zeros(N+1, M+1);               # Для хранения сеток: переходный слой определяется
                                        # на каждом шаге и под него формируется сетка
    XX[:, 1] = Xₙ;                      # Сетка на первом шаге задается автоматом.

    # Векторы, которые подлежат изменению внутри цикла
    # Дальше, в вычислениях должны использоваться только они!
    y = strip_borderPoints(y₀, N);      # Вектор содержащий решение на текущем шаге
    q = strip_borderPoints(qₙ, N);      # Сетка без граничных точек
    X = XX[:, 1];                       # Сетка по X на первом временном шаге

    # Для определения положения переходного слоя
    # Найдем полусумму вырожденных корней на всех шагах по времени
    # используя оригинальную сетку.
    ϕl = phidetermination(qₙ, ulₘ, Xₙ, N, Tₘ, M);                   # Левый вырожденный корень
    ϕr = phidetermination(qₙ, urₘ, Xₙ, N, Tₘ, M, reverseX = true);  # Правый вырожденный корень
    ϕ  = Φ(ϕl, ϕr, N, M);                                           # Серединный корень

    x_tp = find_f_zeros(u[:, 1] - ϕ[:, 1], X);         # Найдем положение переходного слоя
    Xₜₚ[1] = x_tp;

    shouldupdate_mesh = all(!, isnan.(create_mesh(x_tp)));  # Проверяем, передали ли нам функцию

    for m in 1:M
        X = zeros(N+1);                         # Создаем новый массив, чтобы точно все обнулилось
        X = XX[:, m];                           # Используем соответствующую текущему шагу по времени сетку
        τ = (Tₘ[m+1] - Tₘ[m]);                  # Шаг по времени

        j = jac(y, m, X, N, ε, ulₘ, urₘ, q);    # Якобиан от текущего решения
        rp = RP(y, m, X, N, ε, ulₘ, urₘ, q);    # Правая часть от текущего решения

        W = (I - α * τ * j) \ rp;
        y = y .+ τ * real(W);

        u[2:N, m+1] = y       # Сохраним вектор решения на следующем шаге

        if shouldupdate_mesh
            y1 = u[:, m+1];       # Решение на следующем шаге вместе с Граничными Точками

            # Теперь нам нужно определить положение переходного слоя в **следующий** момент времени.
            # Напоминаю, что оно определяет как пересечение решения с полусуммой вырожденных решений.
            # Их мы нашли на оригинальной сетке для всех моментов времени.
            # Найдем `ϕ[:, m+1]` на текущей сетке.
            ϕsp = Spline1D(X, ϕ[:, m+1]);
            ϕₘ =  ϕsp(X);

            try
                x_tp = find_f_zeros(y1 - ϕₘ, X);        # Найдем положение переходного слоя
            catch e
                @error "Прерываем вычисления на шаге $(m), но возвращаем уже посчитанное." e
                return u, XX, Xₜₚ;
            end
            Xₜₚ[m+1] = x_tp;                        # Сохраним его

            XX[:, m+1] = create_mesh(x_tp);         # Обновим сетку

            spl = Spline1D(X, y1)                   # Задаем аппроксимирующий объект: решение на следующем шаге на текущей сетке
            y = spl(XX[2:end-1, m+1])               # Аппроксимируем решением на следующем шаге на сетку следующего шага

            qsp = Spline1D(X[2:end-1], q)           # Задаем аппроксимирующий объект: `q` на текущей сетке
            q = qsp(XX[2:end-1, m+1])               # Аппроксимируем `q` на сетку следующего шага

        else
            XX[:, m+1] = X                          # Если сетку обновлять не нужно, примем за следующую — текущую
        end

    end

    return u, XX, Xₜₚ;
end


text:
######################### Specification for NRL-TB ############################

using DelimitedFiles: readdlm

function default_orbitals(s)
   s = string(s)
   # if s == "Al"  return "spd"
   # if s == "Si"  return "sp"
   if s == "C" || s == "Si"
      return :sp
   else
      return :spd
   end
   error("unkown species in `NRLParams`")
end

# This function reads the ASCII format data files from the NRL server.
#
# The website has been lost, but a functional version remains on the Internet Archive:
# https://web.archive.org/web/20121003160812/http://cst-www.nrl.navy.mil/bind/
#
function NRLHamiltonian(s; orbitals=default_orbitals(s), cutoff=:forceshift)
#   s = string(s)
#   orbitals = string(orbitals)
   if s == :C && orbitals == :sp
      H = C_sp
   elseif s == :Si && orbitals == :sp
      H = Si_sp
   elseif s == :Si && orbitals == :spd
      H = Si_spd
   elseif s == :Al && orbitals == :spd
      H = Al_spd
   else
        # Relatively universal!
        data_dir = joinpath(dirname(Base.find_package("SKTB")), "..", "nrl_data")
        cd(data_dir)
        fname = NRLFILENAME[s] # Lookup species -> filename

        M = readdlm(fname,skipstart=1)
        H =  NRLHamiltonian{9, Function}(9, 10,	   # norbital, nbond
                           M[3,1], M[3,2], cutoff_NRL,			# Rc, lc
                           M[7,1],		# λ
                           [M[8,1],   M[12,1],  M[12,1],  M[12,1],
                            M[16,1],  M[16,1],  M[16,1],  M[16,1], M[16,1]],        #a
                           [M[9,1],   M[13,1],  M[13,1],  M[13,1],
                            M[17,1],  M[17,1],  M[17,1],  M[17,1], M[17,1]],       	#b
                           [M[10,1],  M[14,1],  M[14,1],  M[14,1],
                            M[18,1],  M[18,1],  M[18,1],  M[18,1], M[18,1]],  	   #c
                           [M[11,1],  M[15,1],  M[15,1],  M[15,1],
                            M[19,1],  M[19,1],  M[19,1],  M[19,1], M[19,1]],       	#d

                           [M[24,1],   M[28,1],  M[32,1],  M[36,1],  M[40,1],
                            M[44,1],   M[48,1],  M[52,1],  M[56,1],  M[60,1]],		#e
                           [M[25,1],   M[29,1],  M[33,1],  M[37,1],  M[41,1],
                            M[45,1],   M[49,1],  M[53,1],  M[57,1],  M[61,1]],      #f
                           [M[26,1],   M[30,1],  M[34,1],  M[38,1],  M[42,1],
                            M[46,1],   M[50,1],  M[54,1],  M[58,1],  M[62,1]],      #g
                           [M[27,1],   M[31,1],  M[35,1],  M[39,1],  M[43,1],
                            M[47,1],   M[51,1],  M[55,1],  M[59,1],  M[63,1]],      #h

                           [M[64,1],   M[68,1],   M[72,1],   M[76,1],   M[80,1],
                            M[84,1],   M[88,1],   M[92,1],   M[96,1],   M[100,1]],  #p
                           [M[65,1],   M[69,1],   M[73,1],   M[77,1],   M[81,1],
                            M[85,1],   M[89,1],   M[93,1],   M[97,1],   M[101,1]],  #q
                           [M[66,1],   M[70,1],   M[74,1],   M[78,1],   M[82,1],
                            M[86,1],   M[90,1],   M[94,1],   M[98,1],   M[102,1]],  #r
                           [0.0,   0.0,   0.0,   0.0,   0.0,
                            0.0,   0.0,   0.0,   0.0,   0.0],     						#s
                           [M[67,1],   M[71,1],   M[75,1],   M[79,1],   M[83,1],
                            M[87,1],   M[91,1],   M[95,1],   M[99,1],   M[103,1]],  #t
                           )
   #   error("unkown species / orbitals combination in `NRLParams`")
   end

   if cutoff == :original
      H.fcut = cutoff_NRL_original
   elseif cutoff == :energyshift
      H.fcut = cutoff_NRL_Eshift
   elseif cutoff == :forceshift
      H.fcut = cutoff_NRL_Fshift
   else
      error("unknown cut-off type")
   end
   return H
end

# Lookup table between species and NRL data filename
NRLFILENAME = Dict(
   :Ag => "ag_par",
   :Al => "al_par",
   (:Al, :spd) => "al_par",
#    : => " al_par_t2g_eq_eg",
    :Au => "au_par",
    :Ba => "ba_par",
    :Ca => "ca_par",
#    :Ca => "ca_par_315",
#    :Co => "co_ferro_par",
    :Co => "co_par",
#    : => "co_para_par",
    :C => "c_par",
#    :C => "c_par.105",
    :Cr => "cr_par",
    :Cu => "cu_par",
#    :Cu => "cu_par_99",
#    :Fe => "fe_ferro_par",
#    :Fe => "fe_par",
    :Fe => "fe_para_par",
    :Ga => "ga_par",
    (:Ge, :sp)  => "ge_par.sp.125 ",
    (:Ge, :spd) => "ge_par.spd.125 ",
    :Hf => "hf_par",
    :In => "in_par",
    :Ir => "ir_par",
    :Mg => "mg_par",
    :Mn => "mn_par",
    :Mo => "mo_par",
    :Nb => "nb_par",
    :Ni => "ni_par",
    :Os => "os_par",
    :Pb => "pb_par",
    :Pd => "pd_par",
#    :Pd => "pd_par.105",
    :Pt => "pt_par",
    :Re => "re_par",
    :Rh => "rh_par",
    :Ru => "ru_par",
    :Sc => "sc_par",
    (:Si, :sp) => "si_par",
#    :Si => " si_par.125",
    (:Si, :spd) => "si_par.spd",
    :Sn => "sn_case1_par",
#    :Sn => "sn_case2_par",
    :Sr => "sr_par",
    :Ta => "ta_par",
    :Tc => "tc_par",
#    :Ti => "ti_gga_par",
#    :Ti => "ti_par",
    :Ti => "ti_par_01",
    :V => "v_par",
    :W => "w_par",
    :Y => "y_par",
    :Zr => "zr_par"
)

# Directly entered we have parameters for Si_sp; Si_spd, C_sp and Al_spd

# SILICON
# 'Si' : silicon with s&p orbitals
# reduce Rc = 12.5 to 7.5
Si_sp  =  NRLHamiltonian{4, Function}(4, 4,			    # norbital, nbond
                    12.5, 0.5, cutoff_NRL,			# Rc, lc
                    1.10356625153,		# λ
                    [-0.053233461902,  0.357859715265,  0.357859715265,  0.357859715265],  	#a
                    [-0.907642743185,  0.303647693101,  0.303647693101,  0.303647693101],   	#b
                    [-8.83084913674,   7.09222903560,   7.09222903560,   7.09222903560],     	#c
                    [56.5661321469,    -77.4785508399,  -77.4785508399,  -77.4785508399],     	#d
                    [219.560813651,    10.127687621,    -22.959028107,   10.265449263],  		#e
                    [-16.2132459618,   -4.4036811240,   1.7207707741,   4.6718241428],   		#f
                    [-15.5048968097,   0.2266767834,   1.4191307713,   -2.2161562721],   		#g
                    [1.26439940008,   0.92267194054,   1.03136916513,   1.11134828469],  		#h
                    [1.0, 		     0.0,     		 1.0, 			 1.0],			    	#p
                    [5.157587186,     8.873646665,     11.250489009,   -692.184231145],    	#q
                    [0.660009308,     -16.240770475,  -1.1701322929,    396.153248956],    	#r
                    [-0.0815441307,   5.1822969049,   -1.0591485021,   -13.8172106270],        #s
                    [1.10814448800,   1.24065238343,   1.13762861032,   1.57248559510],  		#t
                   )

# SILICON
# 'Si' : silicon with s&p&d orbitals  :BUT: ignore d-d orbital interactions
# reduce Rc = 12.5 to 7.5
Si_spd  =  NRLHamiltonian{9, Function}(9, 10,			# norbital, nbond
                     12.5, 0.5, cutoff_NRL,			# Rc, lc
                     1.1108,		# λ
                     [-0.0555,  0.4127,   0.4127,   0.4127,
                      0.9691,   0.9691,   0.9691,   0.9691,  0.9691],         	#a
                     [-1.1131,  -0.0907,  -0.0907,  -0.0907,
                      -0.9151,  -0.9151,  -0.9151,  -0.9151, -0.9151],       	#b
                     [-7.3201,  5.3155,   5.3155,   5.3155,
                      -5.9743,  -5.9743,  -5.9743,  -5.9743, -5.9743],  	    #c
                     [74.8905,  -44.0417, -44.0417, -44.0417,
                      602.0289, 602.0289, 602.0289, 602.0289, 602.0289],       	#d

                     [234.6937,   9.5555,  -22.6782,   -1.5942,  -7571.4416,
                      -1.8087,   0.8933,  0.0,   0.0,  0.0],				    #e
                     [-18.6013,   -4.1279,   1.3611,   4.7914,   2.2354,
                      -3.4695,   0.1058,  0.0,  0.0,  0.0],                     #f
                     [-15.0266,   0.2499,    1.3879,   -1.5693,   7.0122,
                      -7.7637,   -0.0224, 0.0,  0.0,  0.0],                     #g
                     [1.2502,    0.8761,    1.01655,    1.1030,   1.6234,
                      1.6294,   0.8217,   0.0,  0.0,  0.0],                     #h

                     [1.0,   0.0,  1.0,   1.0,   0.0,
                      0.0,   0.0,  1.0,   1.0,   1.0],		                    #p
                     [2.4394,   -12.0027,  13.9608,   188.0012,   11.4724,
                      -0.6071,  -2.1340,   0.0,  0.0,  0.0],                    #q
                     [0.9091,   -14.6860,  -1.1961,  -143.3625,  -0.4454,
                      0.05789,  -0.5209,   0.0,  0.0,  0.0],                    #r
                     [-0.0749,   6.1856,   -1.2606,   33.5043,   -0.5838,
                      0.0221,   -0.0948,   0.0,  0.0,  0.0],                    #s
                     [1.0590,   1.2218,    1.1118,    1.4340,   1.0598,
                      0.8130,   1.0580,    0.0,  0.0,  0.0],                    #t
                    )

# CARBON
# 'C' : carbon with s&p orbitals
# reduce Rc = 10.5 to 6.0
C_sp  =  NRLHamiltonian{4, Function}( 4, 4,			    # norbital, nbond
                    10.5, 0.5, cutoff_NRL,			# Rc, lc
                    1.59901905594,		# λ
                    [-0.102789972814,  0.542619178314,  0.542619178314,  0.542619178314],  	#a
                    [-1.62604640052,   2.73454062799,   2.73454062799,   2.73454062799],   	#b
                    [-178.884826119,  -67.139709883,   -67.139709883,  -67.139709883],     	#c
                    [4516.11342028,    438.52883145,   438.52883145,   438.52883145],      	#d
                    [74.0837449667,   -7.9172955767,   -5.7016933899,   24.9104111573],  		#e
                    [-18.3225697598,   3.6163510241,   1.0450894823,   -5.0603652530],   		#f
                    [-12.5253007169,   1.0416715714,   1.5062731505,   -3.6844386855],   		#g
                    [1.41100521808,   1.16878908431,   1.13627440135,   1.36548919302],  		#h
                    [1.0,             0.0,             1.0,             1.0],              	#p
                    [0.18525064246,   1.85250642463,   -1.29666913067,   0.74092406925],   	#q
                    [1.56010486948,   -2.50183774417,   0.28270660019,   -0.07310263856],  	#r
                    [-0.308751658739,   0.178540723033,   -0.022234235553,   0.016694077196],  #s
                    [1.13700564649,   1.12900344616,   0.76177690688,   1.02148246334],  		#t
                   )

# ALUMINIUM
# 'Al' : Aluminum with s&p&d orbitals
# reduce Rc = 16.5 to 9.5
Al_spd  =  NRLHamiltonian{9, Function}(9, 10,			        # norbital, nbond
                     16.5, 0.5, cutoff_NRL,			# Rc, lc
                     1.08515601511,		# λ
                     [-0.0368090795665,  0.394060871550,  0.394060871550,  0.394060871550,
                      1.03732517161,  1.03732517161,  1.03732517161, 1.03732517161, 1.03732517161],         	#a
                     [1.41121601477,  0.996479629379,   0.996479629379,  0.996479629379,
                      2.25910876474,  2.25910876474,  2.25910876474,  2.25910876474,  2.25910876474],       	#b
                     [13.7933378593,   7.02388797078,    7.02388797078,   7.02388797078,
                      -34.3716752929,  -34.3716752929,  -34.3716752929,  -34.3716752929,  -34.3716752929],  	#c
                     [-150.317796096,   -77.7996182049,  -77.7996182049,  -77.7996182049,
                      293.811629762,  293.811629762,  293.811629762,  293.811629762,  293.811629762],       	#d

                     [-45.1824404773,   11.1265443008,  -27.7049616756,   7.48992761254,  -29.0814831341,
                      0.843008479887,   35.6686973234,  -8939.68482560,   -55.7867097600,  41.7418125111],      #e
                     [19.0441568385,   -4.74833564810,   1.14888504976,   3.01675751048,   12.2929753319,
                      -1.52618018997,  -8.20900836372,   730.518353338,   0.853972348751,  -12.0915149851],     #f
                     [-2.81748968422,   0.273179395549,    1.33493438322,   -1.27543114992,   -1.75865704211,
                      0.378014000015,   -0.777295830901,   282.319390496,   2.30786449075,    0.905794614748],  #g
                     [1.05012810116,    0.880778921750,    0.983097680652,    1.01546352470,    1.03433851149,
                      0.964916606335,   1.08355725127,    1.35770591641,    0.997222907112,    0.898850379448], #h

                     [-20.9441522845,   -13.0267838833,   -13.7830192613,   560.641345191,   -103.077140819,
                      20.8403415336,   23.9108771944,   -295.728688028,   -17.2027621079,   -42.9299886946],    #p
                     [17.5240112799,   7.92017585690,   5.15785376743,   -215.309856616,   33.5869977281,
                      -6.65151229760,   -5.86527866863,   80.7470264719,   -6.54916621891,   23.2260645871],    #q
                     [-1.33002925172,   -0.523366384472,   -1.08061004452,   24.2658321255,   -1.86799882707,
                      0.195368101148,   0.725698443913,   -2.93711938026,   1.11096322734,   -0.538315401068],  #r
                     [0.0,   0.0,   0.0,   0.0,   0.0,
                      0.0,   0.0,   0.0,   0.0,   0.0],     													#s
                     [1.06584516722,   0.943623371617,   0.915429215594,   1.17753799190,   0.988337965212,
                      0.873041790591,   0.999293973116,   1.02005972107,   1.01466433826,   1.14341718458],     #t
                    )


text:
"""
$(SIGNATURES)

Convert a judoc html string into a html string (i.e. replace `{{ ... }}` blocks).
"""
function convert_html(hs::AS, allvars::PageVars; isoptim::Bool=false)::String
    isempty(hs) && return hs
    # Tokenize
    tokens = find_tokens(hs, HTML_TOKENS, HTML_1C_TOKENS)

    # Find hblocks ({{ ... }})
    hblocks, tokens = find_all_ocblocks(tokens, HTML_OCB)
    filter!(hb -> hb.name != :COMMENT, hblocks)

    # Find qblocks (qualify the hblocks)
    qblocks = qualify_html_hblocks(hblocks)

    fhs = process_html_qblocks(hs, allvars, qblocks)

    # See issue #204, basically not all markdown links are processed  as
    # per common mark with the JuliaMarkdown, so this is a patch that kind
    # of does
    if haskey(allvars, "reflinks") && allvars["reflinks"].first
        fhs = find_and_fix_md_links(fhs)
    end
    # if it ends with </p>\n but doesn't start with <p>, chop it off
    # this may happen if the first element parsed is an ocblock (not text)
    δ = ifelse(endswith(fhs, "</p>\n") && !startswith(fhs, "<p>"), 5, 0)

    isempty(fhs) && return ""

    if !isempty(GLOBAL_PAGE_VARS["prepath"].first) && isoptim
        fhs = fix_links(fhs)
    end

    return String(chop(fhs, tail=δ))
end

"""
$SIGNATURES

Return the HTML corresponding to a JuDoc-Markdown string.
"""
function jd2html(st::AbstractString; internal::Bool=false, dir::String="")::String
    isempty(st) && return st
    if !internal
        FOLDER_PATH[] = isempty(dir) ? mktempdir() : dir
        set_paths!()
        CUR_PATH[] = "index.md"
        def_GLOBAL_LXDEFS!()
        def_GLOBAL_PAGE_VARS!()
    end
    m, v = convert_md(st * EOS, collect(values(GLOBAL_LXDEFS)); isinternal=internal)
    h = convert_html(m, v)
    return h
end

"""
$SIGNATURES

Take a qualified html block stack and go through it, with recursive calling.
"""
function process_html_qblocks(hs::AS, allvars::PageVars, qblocks::Vector{AbstractBlock},
                              head::Int=1, tail::Int=lastindex(hs))::String
    htmls = IOBuffer()
    head  = head # (sub)string index
    i     = 1    # qualified block index
    while i ≤ length(qblocks)
        β = qblocks[i]
        # write what's before the block
        fromβ = from(β)
        (head < fromβ) && write(htmls, subs(hs, head, prevind(hs, fromβ)))

        if β isa HTML_OPEN_COND
            content, head, i = process_html_cond(hs, allvars, qblocks, i)
            write(htmls, content)
        # should not see an HEnd by itself --> error
        elseif β isa HEnd
            throw(HTMLBlockError("I found a lonely {{end}}."))
        # it's a function block, process it
        else
            write(htmls, convert_html_fblock(β, allvars))
            head = nextind(hs, to(β))
        end
        i += 1
    end
    # write whatever is left after the last block
    head ≤ tail && write(htmls, subs(hs, head, tail))
    return String(take!(htmls))
end

"""
$SIGNATURES

Recursively process a conditional block from an opening HTML_COND_OPEN to a {{end}}.
"""
function process_html_cond(hs::AS, allvars::PageVars, qblocks::Vector{AbstractBlock},
                           i::Int)::Tuple{String,Int,Int}
    if i == length(qblocks)
        throw(HTMLBlockError("Could not close the conditional block " *
                             "starting with '$(qblocks[i].ss)'."))
    end

    init_idx      = i
    else_idx      = 0
    elseif_idx    = Vector{Int}()
    accept_elseif = true # false as soon as we see an else block
    content       = ""

    # inbalance keeps track of whether we've managed to find a
    # matching {{end}}. It increases if it sees other opening {{if..}}
    # and decreases if it sees a {{end}}
    inbalance = 1
    probe = qblocks[i] # just for allocation
    while i < length(qblocks) && inbalance > 0
        i    += 1
        probe = qblocks[i]
        # if we have a {{elseif ...}} or {{else}} with the right level,
        # keep track of them
        if inbalance == 1
            if probe isa HElseIf
                accept_elseif || throw(HTMLBlockError("Saw a {{elseif ...}} after a "*
                                                      "{{else}} block."))
                push!(elseif_idx, i)
            elseif probe isa HElse
                else_idx      = i
                accept_elseif = false
            end
        end
        inbalance += hbalance(probe)
    end
    # we've exhausted the candidate qblocks and not found an appropriate {{end}}
    if inbalance > 0
        throw(HTMLBlockError("Could not close the conditional block " *
                             "starting with '$(qblocks[init_idx].ss)'."))
    end

    # we've found the closing {{end}} at index `i`
    β_close = probe
    i_close = i

    # We now have a complete conditional block with possibly
    # several else-if blocks and possibly an else block.
    # Check which one applies and discard the rest

    k   = 0 # index of the first verified condition
    lag = 0 # aux var to distinguish 1st block case

    # initial block may be {{if..}} or {{isdef..}} or...
    # if its not just a {{if...}}, need to act appropriately
    # and if the condition is verified then k=1
    βi = qblocks[init_idx]
    if βi isa Union{HIsDef,HIsNotDef,HIsPage,HIsNotPage}
        lag = 1
        if βi isa HIsDef
            k = Int(haskey(allvars, βi.vname))
        elseif βi isa HIsNotDef
            k = Int(!haskey(allvars, βi.vname))
        else
            # current path is relative to /src/ for instance
            # /src/pages/blah.md -> pages/blah
            # if starts with `pages/`, replaces by `pub/`: pages/blah => pub/blah
            rpath = splitext(unixify(CUR_PATH[]))[1]
            rpath = replace(rpath, Regex("^pages") => "pub")
            # compare with β.pages
            inpage = any(p -> splitext(p)[1] == rpath, βi.pages)

            if βi isa HIsPage
                k = Int(inpage)
            else
                k = Int(!inpage)
            end
        end
    end

    # If we've not yet found a verified condition, keep looking
    # Now all cond blocks ahead are {{if ...}} {{elseif  ...}}
    if iszero(k)
        if lag == 0
            conds = [βi.vname, (qblocks[c].vname for c in elseif_idx)...]
        else
            conds = [qblocks[c].vname for c in elseif_idx]
        end
        known = [haskey(allvars, c) for c in conds]

        if !all(known)
            idx = findfirst(.!known)
            throw(HTMLBlockError("At least one of the condition variable could not be found: " *
                                 "couldn't find '$(conds[idx])'."))
        end
        # check that all these variables are bool
        bools = [isa(allvars[c].first, Bool) for c  in conds]
        if !all(bools)
            idx = findfirst(.!bools)
            throw(HTMLBlockError("At least one of the condition variable is not a Bool: " *
                                 "'$(conds[idx])' is not a bool."))
        end
        # which one is verified?
        u = findfirst(c -> allvars[c].first, conds)
        k = isnothing(u) ? 0 : u + lag
    end

    if iszero(k)
        # if we still haven't found a verified condition, use the else if there is one
        if !iszero(else_idx)
            # use elseblock, the content is from it to the β_close
            head = nextind(hs, to(qblocks[else_idx]))
            tail = prevind(hs, from(β_close))
            # get the stack of blocks in there
            action_qblocks = qblocks[else_idx+1:i_close-1]
            # process
            content = process_html_qblocks(hs, allvars, action_qblocks, head, tail)
        end
    else
        # determine the span of blocks
        oidx = 0 # opening conditional
        cidx = 0 # closing conditional
        if k == 1
            oidx = init_idx
        else
            oidx = elseif_idx[k-1]
        end
        if k ≤ length(elseif_idx)
            cidx = elseif_idx[k]
        else
            cidx = ifelse(iszero(else_idx), i_close, else_idx)
        end
        # Get the content
        head = nextind(hs, to(qblocks[oidx]))
        tail = prevind(hs, from(qblocks[cidx]))
        # Get the stack of blocks in there
        action_qblocks = qblocks[oidx+1:cidx-1]
        # process
        content = process_html_qblocks(hs, allvars, action_qblocks, head, tail)
    end
    # move the head after the final {{end}}
    head = nextind(hs, to(β_close))

    return content, head, i_close
end


text:
const T_MAX = Float64x2

function compute_svd(A::AbstractMatrix{T_MAX}; n_sv_hint=nothing, strategy=:default)
    if !isnothing(n_sv_hint)
        @info "n_sv_hint is set but will not be used in the current implementation!"
    end
    if strategy ≠ :default
        @info "strategy is set but will not be used in the current implementation!"
    end

    return tsvd(A)
end

function compute_svd(A::AbstractMatrix; n_sv_hint=nothing, strategy=:default)
    if !isnothing(n_sv_hint)
        @info "n_sv_hint is set but will not be used in the current implementation!"
    end
    # m, n = size(A)
    # isnothing(n_sv_hint) && (n_sv_hint = min(m, n))
    # n_sv_hint = min(m, n, n_sv_hint)

    if strategy == :default
        u, s, v = svd(A)
    elseif strategy == :accurate
        u, s, v = svd(A; alg=QRIteration())
    else
        throw(DomainError(strategy, "unknown strategy"))
    end

    return u, s, v
end


text:
using ReinforcementLearning
using StableRNGs
using Flux
using Flux.Losses
using Nature
using ArgParse
# using Serialization
using TensorBoardLogger, Logging
using Infiltrator
using Dates
using Serialization

N_STARTING_PLAYERS = 36
UPDATE_FREQ = 128
WORLD_SIZE = (64, 64)
CLIP = 0.1f0
seed = 123

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--episode-len"
            help = "Maximum number of steps to run"
            arg_type = Int
            default = 500
        "--num-frames"
            help = "Maximum number of steps to run"
            arg_type = Int
            default = 1
        "--max-steps"
            help = "Maximum number of steps to run"
            arg_type = Int
            default = 1_000_000
        "--tb"
            help = "Log to Tensorboard"
            action = :store_true
        "name"
            help = "Name of Run"
            arg_type = String
            required = true
    end

    return parse_args(s)
end

args = parse_commandline()
trial_id = args["name"]

rng = StableRNG(seed)
env = NatureEnv(num_starting_players=N_STARTING_PLAYERS,
                world_size=WORLD_SIZE,
                episode_len=args["episode-len"],
                window=3,
                num_frames=args["num-frames"],
                vocab_size=5,
                food_generators=[
                    FoodGen([15,15],[60,60]),
                    FoodGen([45,45],[60,60]),
                   ])
RLBase.reset!(env)

hook = NatureHook(env, trial_id, args["max-steps"])
stop_condition = StopAfterStep(args["max-steps"], is_show_progress=!haskey(ENV, "CI"))

if args["tb"]
    println("Tensorboard Logging enabled.")
    disable_logging(Logging.Debug)
    lg=TBLogger("tensorboard_logs/$trial_id", tb_append)
    global_logger(lg)
else
    println("Tensorboard Logging not enabled.")
    disable_logging(Logging.Info)
    global_logger(ConsoleLogger())
end

if ispath("checkpoints/$trial_id") && length(readdir("checkpoints/$trial_id/")) > 0
    dir_roots = map(x->x[1:end-4], readdir("checkpoints/$trial_id/"))
    last_step = sort(parse.(Int, dir_roots))[end]
    println("Resuming from checkpoint $last_step")
    last_checkpoint = deserialize("checkpoints/$trial_id/$last_step.jls")
    global_logger().global_step = last_checkpoint[:tb_step]
    agents = last_checkpoint[:policy]
    stop_condition.cur = last_checkpoint[:hook_step]
    hook.step = last_checkpoint[:hook_step]
else
    ispath("checkpoints/$trial_id") || mkdir("checkpoints/$trial_id")
    agents = build_MultiPPOManager(env, N_STARTING_PLAYERS, UPDATE_FREQ, CLIP)
end

ex = Experiment(agents, env, stop_condition, hook, "# PPO with Nature")

run(ex)
serialize("checkpoints/$trial_id.jls", ex.policy)


text:
# We have n particles here.
# For each particle i (i ∈ N⁺)  :
# We have Mₘᵢₙ ≤ mᵢ ≤ Mₘₐₓ
# We have an three dimentional vector v⃗ᵢ for i's velocity.
# And an three dimentional point Pᵢ for i's locationn.
# We have Gravitational Constant G ≈ 6.67408*10⁻¹¹ (m³kg⁻¹s⁻²) .
# And a small time step Δt→0 (actually, it's a very small value).
# At the beginning, every particle will spread in an area randomly.
# At the beginning, every particle will have a velocity which direction and size are all given randomly but in a range.
#
# ⁽ᶜᵒᵈᵉ ᵇʸ ᴿᵉˢᵇⁱ⁾

# Modules
using LinearAlgebra

# Constants
const n = 1024*16
const G = 6.67408e-11
const Δt = 60*12
const Δt² = Δt^2
const Δtₕ = Δt/2
const c = 299792458 # m/s
const Δt⁻²c⁻² = (c*Δt)^(-2)
const Mₘₐₓ = 1e22
const Mₘᵢₙ = 1e20
const Areaₘᵢₙ = [-1e11 -1e11 -1e11]
const Areaₘₐₓ = [1e11 1e11 1e11]
const Areaₜₑₘₚ = Areaₘₐₓ - Areaₘᵢₙ
const v⃗ₘᵢₙ  = [-5e2 -5e2 -5e2]
const v⃗ₘₐₓ  = [5e2 5e2 5e2]
const v⃗ₜₑₘₚ  = v⃗ₘₐₓ  - v⃗ₘᵢₙ
const SavePer = 5
const Total = 1024*16

# Generate
MGΔtΔtₕ = (rand(n)*(Mₘₐₓ-Mₘᵢₙ) + ones(n)*Mₘᵢₙ)*G*Δt*Δtₕ
P₀ = zeros(3, n)
P₁ = zeros(3, n)
Δtv⃗ₕ  = zeros(3, n)
#Δt²a⃗ₕ = zeros(3, n)
for i=1:3
  P₀[i, :] = rand(n)*Areaₜₑₘₚ[i] + ones(n)*Areaₘᵢₙ[i]
  Δtv⃗ₕ[i,  :] = (rand(n)*v⃗ₜₑₘₚ[i]  + ones(n)*v⃗ₘᵢₙ[i])*Δtₕ
end

#Set a blackhole
P₀[:, 1] = [0.0 0.0 0.0]
Δtv⃗ₕ[:, 1] = [0.0 0.0 0.0]
const MGΔtΔtₕ[1] = 1e30*G*Δt*Δtₕ

# Functions
function update(SYNC_ARRAY, P₀, P₁, Δtv⃗ₕ)
  @inbounds Threads.@threads for i=1:n
  #for i=1:n
    Δt²a⃗ₕ  = [0.0, 0.0, 0.0]
    #γ = sqrt(1 - (norm(Δtv⃗ₕ[i])^2)*Δt⁻²c⁻²)
    for j=1:(i-1)
      Δt²a⃗ₕ  += MGΔtΔtₕ[j]/(norm(P₀[:, j] - P₀[:, i])^3)*(P₀[:, j] - P₀[:, i])
    end
    for j=(i+1):n
      Δt²a⃗ₕ  += MGΔtΔtₕ[j]/(norm(P₀[:, j] - P₀[:, i])^3)*(P₀[:, j] - P₀[:, i])
    end
    P₁[:, i] = P₀[:, i] + 2*Δtv⃗ₕ[:, i] + Δt²a⃗ₕ #*γ
    Δtv⃗ₕ[:, i] = Δtv⃗ₕ[:, i] + Δt²a⃗ₕ #*γ
    # Mark
    #SYNC_ARRAY[i] = 1
  end
end

# Run
for j=0:Total-1
  SYNC_ARRAY = zeros(n)
  if mod(j, 2)==0
    @time update(SYNC_ARRAY, P₀, P₁, Δtv⃗ₕ)
  else
    @time update(SYNC_ARRAY, P₁, P₀, Δtv⃗ₕ)
  end
  # Wait
  #while sum(SYNC_ARRAY) < n
  #  pass
  #end
  if mod(j, SavePer)==0
    println(Δtv⃗ₕ[:, 1]/Δtₕ)
    println(P₀[:, 1])
  end
end

max_stars_repo_name: RESBI/N-Body-in-Julia

max_issues_repo_name: RESBI/N-Body-in-Julia

max_forks_repo_name: RESBI/N-Body-in-Julia


text:
struct CHD{K, V} <:AbstractDict{K, V}
  slots ::Vector{UInt8}
  keys ::Vector{K}
  vals ::Vector{V}
  count ::Int

  r ::Vector{UInt64}
  indices ::Vector{UInt16}
end

"""
one-based. negative means not found
"""
function chd_keyindex(c ::CHD{K, V}, key ::K) ::Int where {K, V}
  r0 = c.r[1]
  h = hash(key) ⊻ r0
  i = h % UInt64(length(c.indices))
  @inbounds ri = c.indices[i+1]

  ri >= UInt16(length(c.r)) && return -1

  @inbounds r = c.r[ri+1]
  ti = (h ⊻ r) % UInt64(length(c.keys))

  @inbounds (c.slots[ti+1] == 0x0) && return -1
  @inbounds (c.keys[ti+1] != key) && return -1

  return ti+1
end

@inline function Base.getindex(c ::CHD{K, V}, key) where {K, V}
  ti = chd_keyindex(c, K(key))
  ti <= 0 && throw(KeyError("key $key not found"))
  @inbounds return c.vals[ti]
end

@inline function Base.length(c ::CHD) ::Int
  return c.count
end

@inline function Base.iterate(c ::CHD{K, V}, count::Int=1) ::Union{Nothing, Tuple{Pair{K, V}, Int}} where {K, V}
  # skip empty slots
  while (count <= length(c.keys)) && (c.slots[count] == 0x0)
    count += 1
  end
  count > length(c.keys) && return nothing
  return ((c.keys[count] => c.vals[count]), count+1)
end

@inline function Base.in(k, keys ::Base.KeySet{K, CHD{K,V}}) where {K, V}
  ti = chd_keyindex(keys.dict, K(k))
  return ti > 0
end

@inline function Base.get(h::CHD{K,V}, key, default) where {K, V}
    index = chd_keyindex(h, key)
    @inbounds return (index <= 0) ? default : h.vals[index]::V
end

@inline function Base.getkey(h::CHD{K,V}, key, default) where {K, V}
    index = chd_keyindex(h, key)
    @inbounds return (index <= 0) ? default : h.keys[index]::K
end

max_stars_repo_name: kyungminlee/MinimalPerfectHashing

max_issues_repo_name: kyungminlee/MinimalPerfectHashing

max_forks_repo_name: kyungminlee/MinimalPerfectHashing


text:
"""
Description of the symbolic regression problem.
"""
function symbolic()
  def = algorithm.genetic() do alg
    alg.population = population.simple() do pop
      pop.size = 100
      pop.species = species.simple() do sp
        sp.fitness = fitness.scalar()
        sp.representation = koza.tree() do t
          t.min_depth = 1
          t.max_depth = 18
          t.inputs = ["x::Float64"]
          t.terminals = ["x::Float64"]
          t.non_terminals = [
            "add(x::Float64, y::Float64)::Float64 = x + y"
            "sub(x::Float64, y::Float64)::Float64 = x - y"
            "mul(x::Float64, y::Float64)::Float64 = x * y"
          ]
        end
      end
      pop.breeder = breeder.flat() do br
        br.selection = selection.tournament(2)
        br.mutation = koza.subtree_mutation() do mut
          mut.rate = 0.1
        end
        br.crossover = koza.subtree_crossover(0.9)
      end
    end
    alg.evaluator = evaluator.regression()
    alg.termination = [criterion.generations(1000)]
  end
  algorithm.compose!(def)
end


text:
# In this example we show how to pass a DMRGObserver to
# the dmrg function which allows tracking energy convergence and
# convergence of local operators.
using ITensors

"""
  Get MPO of transverse field Ising model Hamiltonian with field strength h
"""
function tfimMPO(sites, h::Float64)
  # Input operator terms which define a Hamiltonian
  N = length(sites)
  ampo = OpSum()
  for j in 1:(N - 1)
    ampo += -1, "Z", j, "Z", j + 1
  end
  for j in 1:N
    ampo += h, "X", j
  end
  # Convert these terms to an MPO tensor network
  return MPO(ampo, sites)
end

N = 100
sites = siteinds("S=1/2", N)
psi0 = randomMPS(sites; linkdims=16)

# define parameters for DMRG sweeps
sweeps = Sweeps(15)
setmaxdim!(sweeps, 10, 20, 100, 100, 200)
setcutoff!(sweeps, 1E-8)

H = tfimMPO(sites, 1.0)
energy, psi = dmrg(H, psi0, sweeps)

As = map(psi.data) do t
    T = reshape(t.tensor.storage.data, size(t))
    if ndims(T) == 2
      return T
    else
      return permutedims(T, (3, 1, 2))
    end
end

size.(As)
maximum(map(x->maximum(size(x)), As))
length(Ranks(As).middle)
Ranks(As)

using RenyiEntropyOfMagic

ranks=Ranks(
    left_boundary=8,
    right_boundary=8,
    middle=[(8, 8) for _ in 1:98]
)

replicated_expectation(As; ranks)

nsites = length(As)
Os = [Consts.O for _ in 1:length(As)]
lA, rA = first(As), last(As)
lO, rO = first(Os), last(Os)

lK, lU = boundary_decomposition(lA, lO; rank=ranks.left_boundary)
rK, rU = boundary_decomposition(rA, rO; rank=ranks.right_boundary)

K, prevU = rK, rU
idx = length(As)-1
A, O = As[idx], Os[idx]
left, right = ranks.middle[idx-1]
C, U1, U2 = middle_decomposition(A, O; left, right)

M = transpose(U2) * prevU
K = C * kron(M, M) * K