text:
export getBlockFGMRESmem,BlockFGMRES_relaxation,BlockFGMRESmem,isempty
export BlockFGMRES

function BlockFGMRES(AT::SparseMatrixCSC,r0::ArrayTypes,x0::ArrayTypes,inner::Int64,prec::Function,TOL::Float64,
                 verbose::Bool,flexible::Bool,numCores::Int64, mem::FGMRESmem =  getEmptyFGMRESmem())
Az = zeros(eltype(r0),size(r0));
function Afun(z::ArrayTypes)
	SpMatMul(AT,z,Az,numCores);
	return Az;
end
return BlockFGMRES(Afun,r0,x0,inner,prec,TOL,verbose,flexible,numCores,mem);
end

function BlockFGMRES(Afun::Function,R0::ArrayTypes,X0::ArrayTypes,inner::Int64,prec::Function,TOL::Float64,
                 verbose::Bool,flexible::Bool,numCores::Int64, mem::FGMRESmem = getEmptyFGMRESmem())
(n,m) = size(R0);
if isempty(mem)
	warn("Allocating memory in Block GMRES")
	mem = getFGMRESmem(n,flexible,eltype(R0),inner,false);
else
    resetMem(mem);
	if size(mem.V,2)!= inner*m
		error("BlockFGMRES_relaxation: size of Krylov subspace is different than inner times #rhs");
	end
end

Zbig = mem.Z;
Vbig = mem.V;

rnorm0 = vecnorm(R0);
Q = qrfact!(R0); Betta = Q[:R]; R0 = full(Q[:Q],thin=true); # this code is equivalent to (W,Betta) = qr(R0);
W = 0;
Z = 0;

TYPE = eltype(R0);
H = zeros(TYPE,(inner+1)*m,inner*m);
xi = zeros(TYPE,(inner+1)*m,m);
T = zeros(TYPE,inner*m,m);
xi[1:m,:] = Betta;

constOne = one(TYPE);
constZero = zero(TYPE);

rnorms = zeros(inner);

for j = 1:inner

	colSet_j = (((j-1)*m)+1):j*m
	if j==1
		Vbig[:,colSet_j] = R0;  # no memory problem with this line....
		Z = prec(R0); # here unfortunately we have an allocation....
	else
		Vbig[:,colSet_j] = W;  # no memory problem with this line....
		Z = prec(W);
	end
	if flexible
		Zbig[:,colSet_j] = Z;
	end

	W = Afun(Z); # w = A'*z;

	BLAS.gemm!('C','N', constOne, Vbig, W,constZero,T); # t = V'*w;
	H[1:(inner*m),colSet_j] = T;
	BLAS.gemm!('N','N', -constOne, Vbig, T,constOne,W); # w = w - V*t

	Q = qrfact!(W); Betta = Q[:R]; W = full(Q[:Q],thin=true); # this code is equivalent to (W,Betta) = qr(W);

	H[((j*m)+1):(j+1)*m,colSet_j] = Betta;

	y = H[1:(j+1)*m,1:j*m]\xi[1:(j+1)*m,:];
    rnorms[j] = vecnorm(H[1:(j+1)*m,1:j*m]*y - xi[1:(j+1)*m,:])
	if verbose
		if j==1
			resprev = rnorm0;
		else
			resprev = rnorms[j-1];
		end
        # println(string("Inner iter ",j,": BFGMRES: residual norm: ",rnorms[j],", relres: ",rnorms[j]/rnorm0,", factor: ",rnorms[j]/resprev));
    end
	if rnorms[j] < TOL
        rnorms = rnorms[1:j];
        break;
    end
end
y = pinv(H)*xi;
# y = H\xi;

if flexible
	BLAS.gemm!('N','N', constOne, Zbig, y,constZero,W); # w = Z*y  #This is the correction that corresponds to the residual.
else
	# W = Vbig*y;
	BLAS.gemm!('N','N', constOne, Vbig, y, constZero, W);
	Z = prec(W);
	W = Z;
end
addVectors(constOne,W,X0);
return X0,rnorms,mem; # w is the correction vector;
end

type BlockFGMRESmem
	Zbig    			::Array{ArrayTypes}
	AZbig				::Array{ArrayTypes}
end

function resetMem(mem::BlockFGMRESmem)
for k = 1:length(mem.AZbig)
	mem.AZbig[k][:] = 0.0;
end
for k = 1:length(mem.Zbig)
	mem.Zbig[k][:] = 0.0;
end
end

function getEmptyBlockFGMRESmem(n::Int64,m::Int64,flexible::Bool,T::Type,k::Int64)
Zbig  = Array(Array{T,2},0);
AZbig = Array(Array{T,2},0);
return BlockFGMRESmem(Zbig,AZbig);
end

import Base.isempty
function isempty(mem::BlockFGMRESmem)
return length(mem.AZbig)==0;
end

function getBlockFGMRESmem(n::Int64,m::Int64,flexible::Bool,T::Type,inner::Int64)
if m <= 1
	error("Block GMRES is used for 1 rhs. Use GMRES");
end
if flexible
	Zbig = Array(Array{T,2},inner);
	for kk=1:inner
		Zbig[kk] = zeros(T,n,m);
	end
	AZbig = Array(Array{T,2},inner);
	for kk=1:inner
		AZbig[kk] = zeros(T,n,m);
	end
else
	Zbig = Array(Array{T,2},0);
	# Zbig[1] = zeros(T,n,m);
	AZbig = Array(Array{T,2},inner);
	for kk=1:inner
		AZbig[kk] = zeros(T,n,m);
	end
end
return BlockFGMRESmem(Zbig,AZbig);
end

function BlockFGMRES_relaxation(AT::SparseMatrixCSC,r0::ArrayTypes,x0::ArrayTypes,inner::Int64,prec::Function,TOL::Float64,
                 verbose::Bool,flexible::Bool,numCores::Int64, mem::BlockFGMRESmem =  getEmptyBlockFGMRESmem())
Az = zeros(eltype(r0),size(r0));
function Afun(z::ArrayTypes)
	SpMatMul(AT,z,Az,numCores);
	return Az;
end
return BlockFGMRES_relaxation(Afun,r0,x0,inner,prec,TOL,verbose,flexible,numCores,mem);
end

function BlockFGMRES_relaxation(Afun::Function,R0::ArrayTypes,X0::ArrayTypes,inner::Int64,prec::Function,TOL::Float64,
                 verbose::Bool,flexible::Bool,numCores::Int64, mem::BlockFGMRESmem = getEmptyBlockFGMRESmem())
warning("This method is not numerically stable");
n = size(R0,1);
m = size(R0,2);
if isempty(mem)
	warn("Allocating memory in Block GMRES")
	mem = getBlockFGMRESmem(n,m,flexible,eltype(R0),inner);
else
    resetMem(mem);
	if length(mem.AZbig)!= inner
		error("BlockFGMRES_relaxation: size of Krylov subspace is different than inner");
	end
end
Zbig = mem.Zbig; # if flexible
AZbig = mem.AZbig;

rnorm0 = vecnorm(R0);

TYPE = eltype(R0);

H = zeros(TYPE,inner*m,inner*m);
Xi = zeros(TYPE,inner*m,m);
T = zeros(TYPE,inner*m,m);
Told = copy(T);
S = zeros(TYPE,m,m);

constOne = one(TYPE);
constZero = zero(TYPE);

rnorms = zeros(inner);
Z = 0;
W = 0;

for j = 1:inner
	# println("---------------");
	colSet_j = (((j-1)*m)+1):j*m
	if j==1
		Z = prec(R0);
	else
		Z = prec(AZbig[j-1]);
	end
	if flexible
		Zbig[j][:] = Z;
	end
	W = Afun(Z); # w = A'*z;
	AZbig[j][:] = W;

	BLAS.gemm!('C','N', constOne, AZbig[j], R0,constZero,S); # t = V'*w;
	Xi[colSet_j,:] = S;
	# Xi[colSet_j,:] = dot(AZbig[j],R0)
	for jj = 1:j
		BLAS.gemm!('C','N', constOne, AZbig[jj], AZbig[j],constZero,S);
		# S = dot( AZbig[jj], AZbig[j]);
		if j==jj
			S = 0.5*S+0.5*S';
		end
		colSet_jj = (((jj-1)*m)+1):jj*m
		T[colSet_jj,:] = S;
	end

	H[:,colSet_j] = T;
	H[colSet_j,:] = T';
	if j>1
		ee = eig(H)[1];
		ee = ee[ee.!=0.0];
		if ee!=[]
			if minimum(ee)/maximum(ee) < 1e-14
				if verbose
					println("Breaking GMRES primitive due to lack of accuracy in LS solver");
				end
				# T = Told;
				rnorms = rnorms[1:j-1];
				break;
			end
		end
	end

	T = pinv(H)*Xi;

	# println(real(trace(T'*H*T-2.0*T'*Xi))+rnorm0^2)
	rnorms[j] = sqrt(abs(real(trace(T'*H*T-2.0*T'*Xi))+rnorm0^2));
	if verbose
		ee = eig(H)[1];
		ee = ee[ee.!=0.0];
		# println("Condition Number H: ",minimum(ee)/maximum(ee));
		if j==1
			resprev = rnorm0;
		else
			resprev = rnorms[j-1];
		end
        println(string("Inner iter ",j,": FGMRES: residual norm: ",rnorms[j],", relres: ",rnorms[j]/rnorm0,", factor: ",rnorms[j]/resprev));
    end
    if rnorms[j] < TOL
        rnorms = rnorms[1:j];
        break;
    end
end
if flexible
	for jj = 1:length(rnorms)
		colSet_jj = (((jj-1)*m)+1):jj*m
		BLAS.gemm!('N','N', constOne, Zbig[jj], T[colSet_jj,:],constOne,X0);
		# addVectors(T[jj],Zbig[jj],X0);
	end
else
	jj = 1;
	colSet_jj = (((jj-1)*m)+1):jj*m
	BLAS.gemm!('N','N', constOne, R0, T[colSet_jj,:],constZero,Z);
	for jj = 2:length(rnorms)
		colSet_jj = (((jj-1)*m)+1):jj*m
		BLAS.gemm!('N','N', constOne, AZbig[jj-1], T[colSet_jj,:],constOne,Z);
	end
	AZbig[1] = prec(Z);
	addVectors(constOne,AZbig[1],X0);
end

return X0,rnorms,mem; # w is the correction vector;
end


text:
# ins,ans = generate(4)
# open("test.dat","w") do io
#     save(io,ins)
# end

@show generate(10)

ext: jl

max_stars_repo_name: Roger-luo/openQC

max_issues_repo_name: Roger-luo/openQC

max_forks_repo_name: Roger-luo/openQC


text:
using ArgParse
using Dates
using Distributed
@everywhere include("./mining.jl")
@everywhere using Distributed
@everywhere using Dates
@everywhere using .Mining

ccall(:jl_exit_on_sigint, Cvoid, (Cint,), 0)

function parse_cmd()
    s = ArgParseSettings()

    @add_arg_table! s begin
            "data_path"
                help = "data folder relative to data directory"
                required = true
            "--sort_heuristic"
                default = "plain" # alternative is "uncovered_nodes" but this doesn't work as well
                arg_type = String
            "--surprise_threshold"
                default = 10.
                arg_type = Float64
            "--size_threshold"
                default = 10.
                arg_type = Float64
            "--stop_threshold"
                default = Inf
                arg_type = Float64
            "--stop_after_n_structures"
                default = Inf
                arg_type = Float64
            "--overwrite", "-o"
                default = false
                arg_type = Bool
            "--single_process", "-s"
                default = false
                arg_type = Bool
            "--debug", "-d"
                default = false
                arg_type = Bool
    end
    return parse_args(s)
end

@everywhere function build_model_for_dataset(dataset::String, surprise_threshold::Float64, size_threshold::Float64, stop_threshold::Float64, stop_after_n_structures::Float64, sort_heuristic::String; overwrite::Bool=false, debug::Bool=false)
    println("$(now()) - Starting to compute model for $(dataset) with size threshold $(size_threshold), stop threshold $(stop_threshold), and at most $(stop_after_n_structures) structures...")
    flush(stdout)
    Mining.compute_anytime_model(dataset, surprise_threshold, size_threshold, stop_threshold, stop_after_n_structures, sort_heuristic; overwrite=overwrite, debug=debug)
end

function main()
    parsed_args = parse_cmd()
    dsbase = parsed_args["data_path"]
    datasets = [dsbase * "/" * x for x in readdir(dsbase) if (endswith(x, ".txt") || endswith(x, ".csv")) && !occursin("nodemap",x)]
    params = [parsed_args["surprise_threshold"], parsed_args["size_threshold"], parsed_args["stop_threshold"], parsed_args["stop_after_n_structures"], parsed_args["sort_heuristic"]]

    if parsed_args["single_process"]
        for dataset in datasets
            build_model_for_dataset(dataset, params...; overwrite=parsed_args["overwrite"], debug=parsed_args["debug"])
        end
    else
        pmap(x->build_model_for_dataset(x, params...; overwrite=parsed_args["overwrite"], debug=parsed_args["debug"]), datasets)
    end
end

main()

max_stars_repo_name: dataspider/momo

max_issues_repo_name: dataspider/momo

max_forks_repo_name: dataspider/momo


text:
using InteractiveUtils
versioninfo()
#-
using Pkg
Pkg.activate(".")

#-
## Pkg.update()

#-
Pkg.instantiate()

#-
Pkg.precompile()

#-
Pkg.status()

#-
ENV["GKS_ENCODING"]="utf-8" # for scientific formatter
using Plots;
## pyplot()
gr()

# # Standard Plots

include("line_plot_2d_1.jl")
include("line_plot_2d_2.jl")
include("datetime_plot.jl")

include("simple_fn.jl")
include("parametric2d.jl")
include("parametric3d.jl")

include("semilogy.jl")
include("semilogx.jl")
include("loglog.jl")

include("vertical_bar_plot.jl")
include("stacked_bar_plot.jl")
include("horizontal_bar_plot.jl")

include("histogram_plot.jl")
include("categorical_histogram_plot.jl")
include("histogram_2d.jl")

include("boxplot.jl")
include("boxplot2.jl")

include("heatmap.jl")

include("pie.jl")

include("contour_plot.jl")
include("fn_contour_plot.jl")
include("polar.jl")
include("rose.jl")
include("scatter2d.jl")
include("scatter3d.jl")
include("stem.jl")
include("stairstep.jl")
include("error_bar.jl")
include("two_y_axis.jl")
include("subplot1.jl")
include("subplot2.jl")
include("subplot3.jl")
include("quiver.jl")
include("surface.jl")
include("surface2.jl")
include("surface3.jl")

# # Customizing plots
include("standard_line_colors.jl")
include("line_style.jl")
include("latex1.jl")
include("latex2.jl")
include("axis.jl")
include("no_axis.jl")
include("ticks.jl")
include("formatter.jl")
include("axis_range.jl")
include("aspect_ratio.jl")
include("axis_flip.jl")
include("font.jl")


text:
function L2norm(x)
    local l2::eltype(x) = 0
    for xᵢ in x
        l2 += xᵢ^2
    end
    l2
end

function l2norm(x)
    sqrt(L2norm(x))
end


text:
abstract Lazy{T}

type LazyF{T} <: Lazy{T}
    f::Function
    evaluated::Bool
    deps::Vector{LazyF}
    value::T
    LazyF(f::Function) = new(f, false, LazyF[])
end
LazyF(f::Function, T::Type=Any) = LazyF{T}(f)

type LazyV{T} <: Lazy{T}
    deps::Vector{LazyF}
    value::T
    LazyV(value) = new(LazyF[], value)
end
LazyV{T}(value::T) = LazyV{T}(value)

value(l::LazyV) = l.value
value(l::LazyF) = begin
    if !l.evaluated
        l.value = l.f()
        l.evaluated = true
    end
    l.value
end

for op in [+, -, *, /]
    op(x::Lazy, y::Lazy) = begin
        res = LazyF(() -> op(value(x), value(y)))
        push!(x.deps, res)
        push!(y.deps, res)
        res
    end
    op(x, y::Lazy) = begin
        res = LazyF(() -> op(value(x), value(y)))
        push!(y.deps, res)
        res
    end
    op(x::Lazy, y) = begin
        res = LazyF(() -> op(value(x), value(y)))
        push!(x.deps, res)
        res
    end
end

unevaluate(lt::Lazy) = begin
    for l in lt.deps
        l.evaluated=false
        unevaluate(l)
    end
    nothing
end
update!(l::Lazy, value) = begin
    unevaluate(l)
    l.value = value
end

If(cond::Bool, val_true::Function, val_false::Function) =
    cond ? val_true() : val_false

#Testing
x = LazyV(3)
y = LazyV(4)
xpy = x+y
@assert 7 == value(xpy)
update!(x, 5)
@assert xpy.value == 7
@assert 9 == value(xpy)

c = LazyV(true)
v = If(c, val_true, val_false)

#Output types for functions.
immutable Symb{s}
end
Symb(s::Symbol) = Symb{s}()
f_type(s::Symb, args::Any...) = Any


text:
using DataFrames, FITSIO, Unitful, UnitfulAstro, CSV

export PHOENIX,
       update!,
       wave,
       params,
       ALL_PHOENIX_MODELS,
       download_PHOENIX_model

const PHOENIX_REGEX = r"lte(\d{5})-(\d\.\d{2})([\+|-]\d\.\d)(?:\.Alpha=([\+|-]\d\.\d{2}))?\.PHOENIX-ACES-AGSS-COND-2011-HiRes.fits"
const PHOENIX_WAVEFILE = "WAVE_PHOENIX-ACES-AGSS-COND-2011.fits"
const PHOENIX_PARAMETERS = ("T", "logg", "Z", "α")
const PHOENIX_WAVEUNITS = u"angstrom"
const PHOENIX_FLUXUNITS = u"erg/s/cm^3"

# Base definitions

mutable struct PHOENIX <: AbstractSpectralLibrary
    path::String
    entries::DataFrame
    wave::AbstractArray
end

"""
    PHOENIX(path)

Stellar Atmosphere Spectra from Husser et al. (2013)

This interfaces the PHOENIX ACES models as computed by Husser et al. (2013). It is parametrized by the following:

| Parameter |     Range     | Description                                     |
|----------:|:-------------:|:------------------------------------------------|
|       `T` | [2300, 12000] | Teff - The effective temperature in Kelvin      |
|    `logg` |   [0.0, 6.0]  | log(g) - The surface gravity in log-solar units |
|       `Z` |  [-4.0, 1.0]  | [Fe/H] - The iron fraction in log-solar units   |
|       `α` |  [-0.2, 1.2]  | [α/M] - The helium concentration fraction       |

The files must be organized the same way they are organized on the Göttinghen servers, as this constructs a table mapping the parameters to each discovered filename.

# References
[[1]](https://ui.adsabs.harvard.edu/abs/2013A%26A...553A...6H) Husser et al. (2013)

# Examples
```julia-repl
julia> grid = PHOENIX("PHOENIX")
PHOENIX ACES Stellar Atmosphere Spectral Library (Husser et al. 2013)
path           = PHOENIX
parameters     = T, logg, Z, α
units          = wave (Å), flux (erg cm^-3 s^-1)
num. available = 659

# Indexing
julia> grid[1]
(T = 5700, logg = 3.5, Z = 0.5, α = 0.0) => Float32[2.3576384e-16, 2.3964725e-16, 2.4359302e-16, 2.4760205e-16, 2.516754e-16, 2.5581404e-16, 2.60019e-16, 2.642913e-16, 2.6863203e-16, 2.7304218e-16  …  1.06697366e12, 1.146122e12, 1.1551366e12, 1.1556633e12, 1.1557706e12, 1.1557903e12, 1.1554768e12, 1.1539062e12, 1.1525517e12, 1.1532526e12]

# We can iterate over the library
julia> [println(model.T) for (model, flux) in grid[[1, 100, 500]]];
5700
7200
7400

# Accessing certain parameters
julia> grid.logg[1:10]
10-element Array{Float64,1}:
 3.5
 4.0
 4.5
 5.0
 5.5
 6.0
 3.5
 4.0
 4.5
 5.0

```
"""
function PHOENIX(path::String)
    entries = load_PHOENIX_table(path)
    wavefile = joinpath(path, PHOENIX_WAVEFILE)
    # Check if wavefile exists gracefully
    !isfile(wavefile) && error("Could not find wavelength file $(PHOENIX_WAVEFILE) in $(path)")
    wavelengths = read(FITS(wavefile)[1])
    return PHOENIX(path, entries, wavelengths)
end

"""
    (::PHOENIX)(T, logg, Z, α=0; use_units=false)

Return the flux density of the model at the given parameters. Raises an error if the model cannot be found. If `use_units` is true, will return the flux density as a `Unitful.Quantity` with the appropriate units.

# Examples
```julia-repl
julia> grid = PHOENIX("PHOENIX")
PHOENIX ACES Stellar Atmosphere Spectral Library (Husser et al. 2013)
path           = PHOENIX
parameters     = T, logg, Z, α
units          = wave (Å), flux (erg cm^-3 s^-1)
num. available = 659

julia> grid(6000, 4.5, 0)[1:10]
10-element Array{Float32,1}:
 9.125983e-13
 9.261124e-13
 9.398209e-13
 9.537266e-13
 9.67832e-13
 9.821403e-13
 9.966539e-13
 1.011376e-12
 1.0263094e-12
 1.0414569e-12

julia> grid(6000, 4.5, 0, use_units=true)[1:10]
10-element Array{Unitful.Quantity{Float32,𝐌*𝐋^-1*𝐓^-3,Unitful.FreeUnits{(erg, cm^-3, s^-1),𝐌*𝐋^-1*𝐓^-3,nothing}},1}:
  9.125983f-13 erg cm^-3 s^-1
  9.261124f-13 erg cm^-3 s^-1
  9.398209f-13 erg cm^-3 s^-1
  9.537266f-13 erg cm^-3 s^-1
   9.67832f-13 erg cm^-3 s^-1
  9.821403f-13 erg cm^-3 s^-1
  9.966539f-13 erg cm^-3 s^-1
  1.011376f-12 erg cm^-3 s^-1
 1.0263094f-12 erg cm^-3 s^-1
 1.0414569f-12 erg cm^-3 s^-1

julia> grid(6000, 4.5, 2)[1:10]
ERROR: Could not find spectrum with parameters (6000, 4.5, 2, 0)
[...]

julia> grid.(6000:100:6200, 4.5, 0)
3-element Array{Array{Float32,1},1}:
 [9.125983e-13, 9.261124e-13, 9.398209e-13, 9.537266e-13, 9.67832e-13, 9.821403e-13, 9.966539e-13, 1.011376e-12, 1.0263094e-12, 1.0414569e-12  …  1.200507e12, 1.2106235e12, 1.2115639e12, 1.211697e12, 1.211749e12, 1.2117668e12, 1.2116947e12, 1.2112217e12, 1.2107722e12, 1.211132e12]
 [3.4244803e-12, 3.4742658e-12, 3.5247538e-12, 3.5759542e-12, 3.6278769e-12, 3.6805316e-12, 3.733928e-12, 3.788077e-12, 3.842988e-12, 3.898672e-12  …  1.2244243e12, 1.2318009e12, 1.2325055e12, 1.2325888e12, 1.2326135e12, 1.2326139e12, 1.2325272e12, 1.2320685e12, 1.2316515e12, 1.231986e12]
 [1.1416081e-11, 1.1579233e-11, 1.1744649e-11, 1.19123565e-11, 1.2082389e-11, 1.2254776e-11, 1.2429551e-11, 1.2606744e-11, 1.2786388e-11, 1.2968517e-11  …  1.2485481e12, 1.2531918e12, 1.2535282e12, 1.253531e12, 1.2535174e12, 1.2534951e12, 1.2534e12, 1.2529466e12, 1.2525302e12, 1.2528731e12]

```
"""
function (p::PHOENIX)(T, logg, Z, α = 0; use_units = false)
    query(p) = p.T == T && p.logg == logg && p.Z == Z && p.α == α
    filenames = filter(query, p).filename
    length(filenames) == 0 && error("Could not find spectrum with parameters ($T, $logg, $Z, $α)")
    flux = read(FITS(first(filenames))[1])
    return use_units ? flux * PHOENIX_FLUXUNITS : flux
end

# Custom functinos

function load_PHOENIX_table(path::String)
    filenames = DataFrame(T = Int[], logg = Float64[], Z = Float64[], α = Float64[], filename = [])
    # Walk through the root path and search for PHOENIX ACES files
    for (root, dirs, files) in walkdir(path)
        for file in files
            m = match(PHOENIX_REGEX, file)
            # If it looks like a PHOENIX file, parse the parameters from the regex string
            if !isnothing(m)
                T = parse(Int, m.captures[1])
                logg = parse(Float64, m.captures[2])
                Z = parse(Float64, m.captures[3])
                a = isnothing(m.captures[4]) ? 0.0 : parse(Float64, m.captures[4])
                push!(filenames, (T, logg, Z, a, joinpath(root, file)))
            end
        end
    end
    return filenames
end

"""
    update!(::PHOENIX)
    update!(::PHOENIX, path::String)

Goes through the given path (or the current path) for the library and recreates the table of models.
"""
function update!(p::PHOENIX, path::String)
    p.entries = load_PHOENIX_table(path)
    p.path = path
    return p
end

update!(p::PHOENIX) = update!(p, p.path)

"""
    wave(::PHOENIX; use_units=false)

Returns the wavelength array for the library. If `use_units` is true, will return a `Unitful.Quantity` array with the appropriate units.

# Examples
```julia-repl
julia> grid = PHOENIX("PHOENIX")
PHOENIX ACES Stellar Atmosphere Spectral Library (Husser et al. 2013)
path           = PHOENIX
parameters     = T, logg, Z, α
units          = wave (Å), flux (erg cm^-3 s^-1)
num. available = 659

julia> wave(grid)[1:10]
10-element Array{Float64,1}:
 500.0
 500.1
 500.2
 500.3
 500.4
 500.5
 500.6
 500.7
 500.8
 500.9

julia> wave(grid, use_units=true)[1:10]
10-element Array{Unitful.Quantity{Float64,𝐋,Unitful.FreeUnits{(Å,),𝐋,nothing}},1}:
 500.0 Å
 500.1 Å
 500.2 Å
 500.3 Å
 500.4 Å
 500.5 Å
 500.6 Å
 500.7 Å
 500.8 Å
 500.9 Å

```
"""
wave(p::PHOENIX; use_units = false) = use_units ? p.wave * PHOENIX_WAVEUNITS : p.wave

"""
    params(::PHOENIX)

Returns a dataframe of the parameters of each model.

# Examples
```julia-repl
julia> grid = PHOENIX("PHOENIX")
PHOENIX ACES Stellar Atmosphere Spectral Library (Husser et al. 2013)
path           = PHOENIX
parameters     = T, logg, Z, α
units          = wave (Å), flux (erg cm^-3 s^-1)
num. available = 659

julia> ps = params(grid)[1:10, :]
10×4 DataFrames.DataFrame
│ Row │ T     │ logg    │ Z       │ α       │
│     │ Int64 │ Float64 │ Float64 │ Float64 │
├─────┼───────┼─────────┼─────────┼─────────┤
│ 1   │ 5700  │ 3.5     │ 0.5     │ 0.0     │
│ 2   │ 5700  │ 4.0     │ 0.5     │ 0.0     │
│ 3   │ 5700  │ 4.5     │ 0.5     │ 0.0     │
│ 4   │ 5700  │ 5.0     │ 0.5     │ 0.0     │
│ 5   │ 5700  │ 5.5     │ 0.5     │ 0.0     │
│ 6   │ 5700  │ 6.0     │ 0.5     │ 0.0     │
│ 7   │ 5800  │ 3.5     │ 0.5     │ 0.0     │
│ 8   │ 5800  │ 4.0     │ 0.5     │ 0.0     │
│ 9   │ 5800  │ 4.5     │ 0.5     │ 0.0     │
│ 10  │ 5800  │ 5.0     │ 0.5     │ 0.0     │

julia> convert(Matrix, ps)
10×4 Array{Float64,2}:
 5700.0  3.5  0.5  0.0
 5700.0  4.0  0.5  0.0
 5700.0  4.5  0.5  0.0
 5700.0  5.0  0.5  0.0
 5700.0  5.5  0.5  0.0
 5700.0  6.0  0.5  0.0
 5800.0  3.5  0.5  0.0
 5800.0  4.0  0.5  0.0
 5800.0  4.5  0.5  0.0
 5800.0  5.0  0.5  0.0

```
"""
params(p::PHOENIX) = select(p.entries, Not(:filename))

"""
    ALL_PHOENIX_MODELS

`DataFrames.DataFrame` of every single available PHOENIX ACES model.
"""
const ALL_PHOENIX_MODELS = CSV.read(joinpath(@__DIR__, "..", "data", "PHOENIX_MODELS.csv"))

"""
    download_PHOENIX_model(::DataFrame   ; path="PHOENIX")
    download_PHOENIX_model(::DataFrameRow; path="PHOENIX")

Given a row or subframe from the [`ALL_PHOENIX_MODELS`](@ref) dataframe, download into the given path making all directories on the way. Files are downloaded from the Goettingen HTTP server and should not be abused. This will also download the wavelength file if it is not already present.

!!! warning
    Do not abuse the Goettingen HTTP server by spamming this function!

# Examples
```julia-repl
julia> size(ALL_PHOENIX_MODELS, 1)
27605

julia> subset = filter(r -> 6000 ≤ r.T ≤ 6300 && 4.0 ≤ r.logg ≤ 4.5 && r.Z == 0 && r.α == 0, a);

julia> size(subset, 1)
8

julia> download_PHOENIX_model(subset)
[ Info: Downloading wavelength file to PHOENIX/WAVE_PHOENIX-ACES-AGSS-COND-2011.fits
[ Info: Downloading model to PHOENIX/Z-0.0/lte06200-4.00-0.0.PHOENIX-ACES-AGSS-COND-2011-HiRes.fits
[ Info: Downloading model to PHOENIX/Z-0.0/lte06300-4.00-0.0.PHOENIX-ACES-AGSS-COND-2011-HiRes.fits
[ Info: Downloading model to PHOENIX/Z-0.0/lte06100-4.00-0.0.PHOENIX-ACES-AGSS-COND-2011-HiRes.fits
[ Info: Downloading model to PHOENIX/Z-0.0/lte06000-4.00-0.0.PHOENIX-ACES-AGSS-COND-2011-HiRes.fits
[ Info: Downloading model to PHOENIX/Z-0.0/lte06200-4.50-0.0.PHOENIX-ACES-AGSS-COND-2011-HiRes.fits
[ Info: Downloading model to PHOENIX/Z-0.0/lte06300-4.50-0.0.PHOENIX-ACES-AGSS-COND-2011-HiRes.fits
[ Info: Downloading model to PHOENIX/Z-0.0/lte06100-4.50-0.0.PHOENIX-ACES-AGSS-COND-2011-HiRes.fits
[ Info: Downloading model to PHOENIX/Z-0.0/lte06000-4.50-0.0.PHOENIX-ACES-AGSS-COND-2011-HiRes.fits

```
"""
function download_PHOENIX_model(model::DataFrameRow; path = "PHOENIX")
    # Set up folders
    folder = joinpath(path, model.folder)
    !isdir(folder) && mkpath(folder)

    # Download wavelength file
    if !isfile(joinpath(path, PHOENIX_WAVEFILE))
        front = "http://phoenix.astro.physik.uni-goettingen.de/data/HiResFITS"
        url = join([front, PHOENIX_WAVEFILE], "/")
        filename = joinpath(path, PHOENIX_WAVEFILE)
        @info "Downloading wavelength file to $filename"
        download(url, filename)
    end

    # Download model
    front = "http://phoenix.astro.physik.uni-goettingen.de/data/HiResFITS/PHOENIX-ACES-AGSS-COND-2011"
    url = join([front, model.folder, model.filename], "/")
    filename = joinpath(folder, model.filename)
    if !isfile(filename)
        @info "Downloading model to $filename"
        download(url, filename)
    end
end

download_PHOENIX_model(models::DataFrame; path = "PHOENIX") = download_PHOENIX_model.(eachrow(models), path = path)

# Extended functions

function Base.show(io::IO, s::PHOENIX)
    println(io, "PHOENIX ACES Stellar Atmosphere Spectral Library (Husser et al. 2013)")
println(io, "path           = $(s.path)")
    pstring = join(PHOENIX_PARAMETERS, ", ")
    println(io, "parameters     = $pstring")
    println(io, "units          = wave ($PHOENIX_WAVEUNITS), flux ($PHOENIX_FLUXUNITS)")
    println(io, "num. available = $(length(s))")
end

"""
    length(::PHOENIX)

Returns the number of models in the library.
"""
Base.length(p::PHOENIX) = size(p.entries, 1)

function Base.iterate(p::PHOENIX, row::Int)
    row > length(p) && return nothing
    return p[row], row + 1
end

Base.iterate(p::PHOENIX) = iterate(p, 1)

function Base.getproperty(p::PHOENIX, name::Symbol)
    if name in propertynames(p)
        return getfield(p, name)
    elseif name in propertynames(p.entries)
        return getproperty(p.entries, name)
    end
end

function Base.getindex(p::PHOENIX, idx::Int)
    entry = p.entries[idx, :]
    params = copy(entry[Not(:filename)])
    filename = entry[:filename]
    flux = read(FITS(filename)[1])
    return params => flux
end

function Base.getindex(p::PHOENIX, idxs)
    return getindex.(Ref(p), idxs)
end

"""
    filter!(function, ::PHOENIX)

In-place version of [`filter`](@ref)
"""
function Base.filter!(fn, p::PHOENIX)
    p.entries = filter(fn, p.entries)
    return p
end

"""
    filter(function, ::PHOENIX)

Returns a library with parameters filtered by the function.

# Examples
```julia-repl
julia> grid = PHOENIX("PHOENIX")
PHOENIX ACES Stellar Atmosphere Spectral Library (Husser et al. 2013)
path           = PHOENIX
parameters     = T, logg, Z, α
units          = wave (Å), flux (erg cm^-3 s^-1)
num. available = 659

julia> query(model) = 6000 < model.T < 7000 &&
               4.0 < model.logg < 5.0 &&
               -1.0 < model.Z < 1.0 &&
               model.α == 0
query (generic function with 1 method)

julia> new_grid = filter(query, grid)
PHOENIX ACES Stellar Atmosphere Spectral Library (Husser et al. 2013)
path           = PHOENIX
parameters     = T, logg, Z, α
units          = wave (Å), flux (erg cm^-3 s^-1)
num. available = 27

```
"""
function Base.filter(fn, p::PHOENIX)
    return filter!(fn, deepcopy(p))
end

"""
    sort!(p::PHOENIX, cols;
          alg::Union{Algorithm, Nothing}=nothing, lt=isless, by=identity,
          rev::Bool=false, order::Ordering=Forward)

In-place version of [`sort`](@ref)
"""
function Base.sort!(p::PHOENIX, cols;
    alg::Union{DataFrames.Algorithm,Nothing} = nothing, lt = isless, by = identity,
    rev::Bool = false, order::DataFrames.Ordering = DataFrames.Forward)
    p.entries = sort!(p.entries, cols, alg = alg, lt = lt, by = by, rev = rev, order = order)
    return p
end

"""
    sort(p::PHOENIX, cols;
         alg::Union{Algorithm, Nothing}=nothing, lt=isless, by=identity,
         rev::Bool=false, order::Ordering=Forward)

Return a copy of the model library sorted by parameter(s) cols. cols can be either a Symbol or Integer column index, or a tuple or vector of such indices.

If alg is nothing (the default), the most appropriate algorithm is chosen automatically among TimSort, MergeSort and RadixSort depending on the type of the sorting columns and on the number of rows in df. If rev is true, reverse sorting is performed. To enable reverse sorting only for some columns, pass order(c, rev=true) in cols, with c the corresponding column index (see example below). See sort! for a description of other keyword arguments.

# Examples
```julia-repl
julia> grid = PHOENIX("PHOENIX")
PHOENIX ACES Stellar Atmosphere Spectral Library (Husser et al. 2013)
path           = PHOENIX
parameters     = T, logg, Z, α
units          = wave (Å), flux (erg cm^-3 s^-1)
num. available = 659

julia> params(grid)[1:10, :]
10×4 DataFrames.DataFrame
│ Row │ T     │ logg    │ Z       │ α       │
│     │ Int64 │ Float64 │ Float64 │ Float64 │
├─────┼───────┼─────────┼─────────┼─────────┤
│ 1   │ 5700  │ 3.5     │ 0.5     │ 0.0     │
│ 2   │ 5700  │ 4.0     │ 0.5     │ 0.0     │
│ 3   │ 5700  │ 4.5     │ 0.5     │ 0.0     │
│ 4   │ 5700  │ 5.0     │ 0.5     │ 0.0     │
│ 5   │ 5700  │ 5.5     │ 0.5     │ 0.0     │
│ 6   │ 5700  │ 6.0     │ 0.5     │ 0.0     │
│ 7   │ 5800  │ 3.5     │ 0.5     │ 0.0     │
│ 8   │ 5800  │ 4.0     │ 0.5     │ 0.0     │
│ 9   │ 5800  │ 4.5     │ 0.5     │ 0.0     │
│ 10  │ 5800  │ 5.0     │ 0.5     │ 0.0     │

julia> logg_sorted = sort(p, :logg, rev=true)
PHOENIX ACES Stellar Atmosphere Spectral Library (Husser et al. 2013)
path           = /Users/miles/dev/starfish/examples/PHOENIX
parameters     = T, logg, Z, α
units          = wave (Å), flux (erg cm^-3 s^-1)
num. available = 659

julia> params(logg_sorted)[1:10, :]
10×4 DataFrames.DataFrame
│ Row │ T     │ logg    │ Z       │ α       │
│     │ Int64 │ Float64 │ Float64 │ Float64 │
├─────┼───────┼─────────┼─────────┼─────────┤
│ 1   │ 5700  │ 6.0     │ 0.5     │ 0.0     │
│ 2   │ 5800  │ 6.0     │ 0.5     │ 0.0     │
│ 3   │ 5900  │ 6.0     │ 0.5     │ 0.0     │
│ 4   │ 6000  │ 6.0     │ 0.5     │ 0.0     │
│ 5   │ 6100  │ 6.0     │ 0.5     │ 0.0     │
│ 6   │ 6200  │ 6.0     │ 0.5     │ 0.0     │
│ 7   │ 6300  │ 6.0     │ 0.5     │ 0.0     │
│ 8   │ 6400  │ 6.0     │ 0.5     │ 0.0     │
│ 9   │ 6500  │ 6.0     │ 0.5     │ 0.0     │
│ 10  │ 6600  │ 6.0     │ 0.5     │ 0.0     │

```
"""
function Base.sort(p::PHOENIX, cols;
    alg::Union{DataFrames.Algorithm,Nothing} = nothing, lt = isless, by = identity,
    rev::Bool = false, order::DataFrames.Ordering = DataFrames.Forward)
    return sort!(deepcopy(p), cols, alg = alg, lt = lt, by = by, rev = rev, order = order)
end


text:
# SPDX-License-Identifier: X11
# 2020-12-22
# Day 6, Part 1

input = split.(
  split(readchomp(stdin), "\n\n"; keepempty=false),
  "\n";
  keepempty=false)

(println ∘ sum ∘ map)(length ∘ unique ∘ (x -> reduce(*, x)), input)

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
include("guiding_center_4d_settings.jl")

set_config(:ls_solver, :julia)
set_config(:quasi_newton_refactorize, 1)
set_config(:nls_atol, Dec128(1E-32))
set_config(:nls_rtol, Dec128(1E-32))
set_config(:nls_stol, Dec128(1E-32))
set_config(:nls_nmax, 20)


text:
import PyCall: hasproperty # Base.hasproperty only in Julia 1.2

struct LazyHelp
    o
    keys::Tuple{Vararg{String}}
    LazyHelp(o) = new(o, ())
    LazyHelp(o, k::AbstractString) = new(o, (k,))
end
function show(io::IO, st, h::LazyHelp)
    o = h.o
    for k in h.keys
        o = o[k]
    end
    print(io, _generate_docstring(o))
end
Base.show(io::IO, h::LazyHelp) = show(io, "text/plain", h)
function Base.Docs.catdoc(hs::LazyHelp...)
    Base.Docs.Text() do io
        for h in hs
            show(io, MIME"text/plain"(), h)
        end
    end
end

function _generate_docstring(f)
    py"""
    import inspect

    def get_signature(f):
        try:
            return str(inspect.signature(f))
        except ValueError:
            return ""
    """

    fname = hasproperty(f, "__name__") ? f.__name__ : ""
    sig = hasproperty(f, "__call__") ? py"get_signature"(f) : ""
    fdoc = hasproperty(f, "__doc__") ? f.__doc__ : ""

    return """
        $(fname)$(sig)
    $(fdoc)
    """
end


text:
using Convex, FirstOrderSolvers, SCS, Test, ProximalOperators

ys = [-0.0064709 -0.22443;
     -0.22443    -1.02411]
y = Variable((2, 2))

#Solve with SCS
p = minimize(norm(vec(y-ys)), isposdef(y))
solve!(p, SCSSolver(eps=1e-8, verbose=0))

ysol  = copy(y.value)

#Solve with projection
Y = Symmetric(ys)
X = similar(Y)
prox!(X, IndPSD(), Y)

#Testing SCS and projection equal
@test X.data ≈ ysol atol=1e-8

#Solve with DR
p = minimize(norm(vec(y-ys)), isposdef(y))
solve!(p, DR(eps=1e-8, verbose=0))

@test y.value ≈ ysol atol=1e-8


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct Melt <: Edit
        on::Array{Symbol,1}
        var::Symbol
        val::Symbol
    end

Normalize the dataframe by 'melting' columns into rows, lengthening the dataframe by duplicating values in the column `on` into new rows and defining 2 new columns: 1. `var` with header names from the original dataframe. 2. `val` with column values from the original dataframe. This operation can only be performed once per dataframe.

# Arguments
- `on::Array{Symbol,1}`: name of column(s) NOT included in melt
- `var::Symbol`: name of column containing header NAMES from the original dataframe
- `val::Symbol`: name of column containing VALUES from the original dataframe
"""
mutable struct Melt <: Edit
    "name of column(s) NOT included in melt"
    on::Array{Symbol,1}
    "name of column containing header NAMES from the original dataframe"
    var::Symbol
    "name of column containing VALUES from the original dataframe"
    val::Symbol
end

function Melt(; on, var, val, )
    Melt(on, var, val, )
end

"""Get [`Melt`](@ref) `on`."""
get_on(value::Melt) = value.on
"""Get [`Melt`](@ref) `var`."""
get_var(value::Melt) = value.var
"""Get [`Melt`](@ref) `val`."""
get_val(value::Melt) = value.val

"""Set [`Melt`](@ref) `on`."""
set_on!(value::Melt, val) = value.on = val
"""Set [`Melt`](@ref) `var`."""
set_var!(value::Melt, val) = value.var = val
"""Set [`Melt`](@ref) `val`."""
set_val!(value::Melt, val) = value.val = val

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
@testset "629.k-inverse-pairs-array" begin
    @test k_inverse_pairs(1000, 1000) == 663677020
    @test k_inverse_pairs(3, 0) == 1
    @test k_inverse_pairs(3, 1) == 2
end


text:
export EpisodicCompactSARTSATrajectory

using MacroTools: @forward

"""
    EpisodicCompactSARTSATrajectory(; state_type = Int, action_type = Int, reward_type = Float32, terminal_type = Bool)

Exactly the same with [`VectorialCompactSARTSATrajectory`](@ref). It only exists for multiple dispatch purpose.

!!! warning
    The `EpisodicCompactSARTSATrajectory` will not be automatically emptified when reaching the end of an episode.
"""
struct EpisodicCompactSARTSATrajectory{types,trace_types} <:
       AbstractTrajectory{SARTSA,types}
    trajectories::VectorialCompactSARTSATrajectory{types,trace_types}
end

EpisodicCompactSARTSATrajectory(; kwargs...) =
    EpisodicCompactSARTSATrajectory(VectorialCompactSARTSATrajectory(; kwargs...))

@forward EpisodicCompactSARTSATrajectory.trajectories Base.length,
Base.isempty,
Base.empty!,
Base.push!,
Base.pop!

# avoid method ambiguous
RLBase.get_trace(t::EpisodicCompactSARTSATrajectory, s::Symbol) =
    get_trace(t.trajectories, s)
Base.getindex(t::EpisodicCompactSARTSATrajectory, i::Int) = getindex(t.trajectories, i)
Base.pop!(t::EpisodicCompactSARTSATrajectory, s::Symbol...) = pop!(t.trajectories, s...)


text:
module Karger
using Base.Threads
using SparseArrays
include("./utilities.jl")
struct Graph
    n :: Int
    edges :: Matrix{Int}
    weights :: Vector{Float64}
end

struct SeedData
    uf :: UnionFind
    is_fixed :: Vector{Bool}
    seed_list :: Vector{Int}
end

function init_seeds(seeds :: Vector{Int})
    uf = UnionFind(length(seeds))
    # True if a seed is at that position
    is_fixed = (seeds .> 0)
    # Dictionary label => union-find index
    seed_roots = Dict{Int, Int}()
    # precontract all seeds
    for i in 1:length(seeds)
        if is_fixed[i]
            if haskey(seed_roots, seeds[i])
                uf[i] = seed_roots[seeds[i]]
            else
                seed_roots[seeds[i]] = i
            end
        end
    end
    return SeedData(uf, is_fixed, collect(keys(seed_roots)))
end

function karger(g :: Graph, s :: Int, t :: Int)
    seeds = zeros(Int, g.n)
    seeds[s] = 1
    seeds[t] = 2
    karger(g, seeds)
end

function karger(g :: Graph, seeds :: Vector{Int})
    karger(g, seeds, init_seeds(seeds))
end

function karger_initialized(g :: Graph, seeds :: Vector{Int}, sd :: SeedData)
    m = size(g.edges, 2)
    karger!(g, seeds, sd, Vector{Float64}(undef, m))
end

function karger_initialized!(g :: Graph, seeds :: Vector{Int}, sd :: SeedData, scores :: Vector{Float64})
    # randexp samples from the exponential distribution with scale 1.
    # We want to sample from p(score > t) = exp(-wt). This means we have to
    # _divide by_ w! See e.g. http://www.math.wm.edu/~leemis/chart/UDR/PDFs/ExponentialS.pdf
    randexp!(scores)
    scores .= scores ./ g.weights
    # this would be argsort in numpy
    perm = sortperm(scores)

    uf, is_fixed = copy(sd.uf), copy(sd.is_fixed)

    num_clusters = g.n
    for i in 1:length(g.weights)
        u, v = g.edges[:, perm[i]]
        u_root, v_root = find(uf, u), find(uf, v)

        if is_fixed[u_root]
            if !is_fixed[v_root]
                # u has a fixed label, so merge v into u
                merged = union!(uf, u, v)
            else
                # u and v's labels are both fixed and they can't be contracted
                continue
            end
        else
            # merge u into v
            merged = union!(uf, v, u)
        end

        if merged
            num_clusters -= 1
            if num_clusters == length(sd.seed_list)
                break
            end
        end
    end
    return uf
end

function potential(g :: Graph, N :: Int, seeds :: Vector{Int})
    sd = init_seeds(seeds)
    probs = Dict(seed => zeros(g.n, nthreads()) for seed in sd.seed_list)
    scores = [Vector{Float64}(undef, size(g.edges, 2)) for i in 1:nthreads()]
    @threads for i in 1:N
        uf = karger_initialized!(g, seeds, sd, scores[threadid()])
        for j in 1:g.n
            probs[seeds[find(uf, j)]][j, threadid()] += 1/N
        end
    end
    return Dict(seed => sum(probs[seed], dims=2) for seed in sd.seed_list)
end

function watershed(g :: Graph, seeds :: Vector{Int})
    m = size(g.edges, 2)
    sd = init_seeds(seeds)

    # We want to merge edges in order of descending weight
    perm = sortperm(g.weights, rev=true)

    uf, is_fixed = sd.uf, sd.is_fixed

    num_clusters = g.n
    for i in 1:length(g.weights)
        u, v = g.edges[:, perm[i]]
        u_root, v_root = find(uf, u), find(uf, v)

        if is_fixed[u_root]
            if !is_fixed[v_root]
                # u has a fixed label, so merge v into u
                merged = union!(uf, u, v)
            else
                # u and v's labels are both fixed and they can't be contracted
                continue
            end
        else
            # merge u into v
            merged = union!(uf, v, u)
        end

        if merged
            num_clusters -= 1
            if num_clusters == length(sd.seed_list)
                break
            end
        end
    end
    segmentation = Vector{Float64}(undef, g.n)
    for i in 1:g.n
        segmentation[i] = seeds[find(uf, i)]
    end
    return segmentation
end

function power_watershed(g :: Graph, seeds :: Vector{Int})
    m = size(g.edges, 2)
    sd = init_seeds(seeds)

    # We want to merge edges in order of descending weight
    perm = sortperm(g.weights, rev=true)
    weights = g.weights[perm]
    edges = g.edges[:, perm]

    uf, is_fixed = sd.uf, sd.is_fixed

    num_unfixed_clusters = g.n

    potential = Vector{Float64}(undef, g.n)
    for i in 1:g.n
        if is_fixed[i]
            num_unfixed_clusters -= 1
            # seeds are 1 and 2, we want to make those 0 or 1
            potential[i] = convert(Float64, seeds[i]) - 1
        end
    end

    # we pre-allocate arrays for the local seeds, is_fixed values and edges
    # in each plateau. These are too large but it's still a small speed-up
    # compared to re-allocating an array for each plateau
    seeds_ = Vector{Float64}(undef, g.n)
    is_fixed_ = falses(g.n)
    edges_ = Matrix{Int}(undef, 2, length(edges))
    # we also pre-allocate buffers needed internally by find_plateaus
    label_buffer = Vector{Int}(undef, g.n)
    uf_buffer = Vector{Int}(undef, g.n)

    i = 0
    # we use a while loop so that we can skip a chunk of iterations later
    # by modifying i
    while i < m
        if num_unfixed_clusters == 0
            break
        end

        i += 1

        if i == m || weights[i + 1] < weights[i]
            # no plateau (only one edge), so just do the watershed procedure
            # to save time

            u, v = edges[:, i]
            u_root, v_root = find(uf, u), find(uf, v)

            if is_fixed[u_root] && is_fixed[v_root]
                # u and v's potentials are both fixed and they can't be contracted
                continue
            end
            if is_fixed[u_root]
                # u has a fixed potential, so merge v into u
                merged = union!(uf, u, v)
            else
                # merge u into v
                merged = union!(uf, v, u)
            end
            if merged
                num_unfixed_clusters -= 1
            end

            continue
        end

        indices = i:m
        plateau_nodes, plateau_edges, num_edges = find_plateaus!(weights[indices], edges[:, indices], uf, label_buffer, uf_buffer)
        # solve each plateau individually.
        for (key, nodes_) in plateau_nodes
            # translate the node ids to 1:length(nodes_):
            old_to_new_id = Dict{Int, Int}()
            fill!(is_fixed_, false)
            found_potentials = Set{Float64}()
            n_ = length(nodes_)
            m_ = length(plateau_edges[key])
            for (j, node) in enumerate(nodes_)
                old_to_new_id[node] = j
                if is_fixed[find(uf, node)]
                    seeds_[j] = potential[find(uf, node)]
                    is_fixed_[j] = true
                    push!(found_potentials, seeds_[j])
                end
            end

            if length(found_potentials) == 0
                # no fixed nodes found in this plateau, just merge the nodes
                for node in nodes_
                    if union!(uf, node, nodes_[1])
                        num_unfixed_clusters -= 1
                    end
                end
                continue
            elseif length(found_potentials) == 1
                # Only one potential occurs in this plateau, assign that to every node.
                # This isn't necessary but it's much faster than constructing the Laplacian
                # and solving the Random Walker for the plateau
                pot = first(found_potentials)
                for node in nodes_
                    if !is_fixed[find(uf, node)]
                        num_unfixed_clusters -= 1
                        is_fixed[find(uf, node)] = true
                        potential[find(uf, node)] = pot
                    end
                end
                continue
            end

            # Create a list of the edges themselves from the list of edge indices.
            # Additionally, we use old_to_new_id to transform the node identifiers
            # (we want edge entries to be wrt the local node ids for the plateau).
            for (j, edge_idx) in enumerate(plateau_edges[key])
                edge_idx += i - 1
                edges_[1, j] = old_to_new_id[find(uf, edges[1, edge_idx])]
                edges_[2, j] = old_to_new_id[find(uf, edges[2, edge_idx])]
            end

            # Now we construct the Laplacian for the plateau.
            # We need to be careful about the fact that we are working
            # with a multi-graph: there might be self-loops (which we ignore)
            # and multiple edges between the same nodes (in which case we add up the weights)

            # if the plateau is relatively small, we just construct and solve
            # a dense system
            if n_ <= 500
                lap = zeros(n_, n_)
                for j in 1:m_
                    # ignore self-loops
                    if edges_[1, j] == edges_[2, j]
                        continue
                    end
                    # Fill the Laplacian with the negative adjacency matrix.
                    # Note that there may be multiple edges between two nodes,
                    # so we need -= 1 instead of = -1. We can use 1 instead of
                    # the actual weight because all edges inside the plateau
                    # have the same weight.
                    lap[edges_[1, j], edges_[2, j]] -= 1
                    lap[edges_[2, j], edges_[1, j]] -= 1
                end
                # add the diagonal, which contains the degrees of each node
                degrees = -sum(lap, dims=1)
                for j in 1:length(nodes_)
                    lap[j, j] = degrees[j]
                end

            # if the plateau is large, it is likely very sparse and we
            # use a sparse Laplacian matrix instead
            else
                row_indices = Vector{Int}()
                col_indices = Vector{Int}()
                values = Vector{Float64}()
                for j in 1:m_
                    # ignore self-loops
                    if edges_[1, j] == edges_[2, j]
                        continue
                    end
                    # we don't need to worry about multi-edges,
                    # they will be combined anyway when forming the sparse
                    # array
                    push!(row_indices, edges_[1, j])
                    push!(col_indices, edges_[2, j])
                    push!(values, -1)

                    # also add the edge in the other direction
                    push!(row_indices, edges_[2, j])
                    push!(col_indices, edges_[1, j])
                    push!(values, -1)

                    # Finally, we add a 1 on the diagonal for both nodes
                    # (again, they are added up later automatically)
                    push!(row_indices, edges_[1, j])
                    push!(col_indices, edges_[1, j])
                    push!(values, 1)

                    push!(row_indices, edges_[2, j])
                    push!(col_indices, edges_[2, j])
                    push!(values, 1)
                end
                lap = sparse(row_indices, col_indices, values, length(nodes_), length(nodes_))
            end

            # then, we construct the linear system
            # TODO: this indexing is slow for sparse matrices,
            # is there a way to do this already during construction?
            # is_fixed_ is too long, we only need the first n_ entries
            is_fixed_view = is_fixed_[1:n_]
            lap_unseeded = lap[.!is_fixed_view, .!is_fixed_view]
            BT = lap[.!is_fixed_view, is_fixed_view]
            # seeds_ is just as long as is_fixed, so use everything here
            b = - BT * seeds_[is_fixed_]
            # TODO: should we tell Julia that the Laplacian is symmetric
            # and positive definite, to improve performance?
            x = lap_unseeded \ b

            # We now have the solution for the unseeded nodes of the
            # current plateau in x. We just need to fill them into
            # the right places.
            # k is the index for x, i.e. it enumerates the unseeded nodes
            # in the plateau. j enumerates all nodes in the plateau.
            # node gives the global index of the node j in the plateau.
            k = 1
            for (j, node) in enumerate(nodes_)
                if is_fixed_[j]
                    continue
                end
                potential[find(uf, node)] = x[k]
                is_fixed[find(uf, node)] = true
                num_unfixed_clusters -= 1
                k += 1
            end
        end
        # now we skip the remaining edges in the plateau
        i += num_edges - 1
    end

    for j in 1:g.n
        if !is_fixed[find(uf, j)]
            error("Something went wrong, not all nodes are fixed yet")
        end
        potential[j] = potential[find(uf, j)]
    end
    return potential
end

function power_watershed_multi(g :: Graph, seeds :: Vector{Int})
    m = size(g.edges, 2)
    sd = init_seeds(seeds)

    # We want to merge edges in order of descending weight
    perm = sortperm(g.weights, rev=true)
    weights = g.weights[perm]
    edges = g.edges[:, perm]

    uf, is_fixed, seed_list = sd.uf, sd.is_fixed, sd.seed_list

    num_unfixed_clusters = g.n
    num_seeds = length(seed_list)

    potentials = zeros(Float64, g.n, num_seeds)
    for i in 1:g.n
        if is_fixed[i]
            num_unfixed_clusters -= 1
            potentials[i, seeds[i]] = 1
        end
    end

    # we pre-allocate arrays for the local seeds, is_fixed values and edges
    # in each plateau. These are too large but it's still a small speed-up
    # compared to re-allocating an array for each plateau
    seeds_ = Matrix{Float64}(undef, g.n, num_seeds)
    is_fixed_ = falses(g.n)
    edges_ = Matrix{Int}(undef, 2, length(edges))
    # we also pre-allocate buffers needed internally by find_plateaus
    label_buffer = Vector{Int}(undef, g.n)
    uf_buffer = Vector{Int}(undef, g.n)

    i = 0
    # we use a while loop so that we can skip a chunk of iterations later
    # by modifying i
    while i < m
        if num_unfixed_clusters == 0
            break
        end

        i += 1

        if i == m || weights[i + 1] < weights[i]
            # no plateau (only one edge), so just do the watershed procedure
            # to save time

            u, v = edges[:, i]
            u_root, v_root = find(uf, u), find(uf, v)

            if is_fixed[u_root] && is_fixed[v_root]
                # u and v's potentials are both fixed and they can't be contracted
                continue
            end
            if is_fixed[u_root]
                # u has a fixed potential, so merge v into u
                merged = union!(uf, u, v)
            else
                # merge u into v
                merged = union!(uf, v, u)
            end
            if merged
                num_unfixed_clusters -= 1
            end

            continue
        end

        indices = i:m
        plateau_nodes, plateau_edges, num_edges = find_plateaus!(weights[indices], edges[:, indices], uf, label_buffer, uf_buffer)
        # solve each plateau individually.
        for (key, nodes_) in plateau_nodes
            # translate the node ids to 1:length(nodes_):
            old_to_new_id = Dict{Int, Int}()
            fill!(is_fixed_, false)
            fill!(seeds_, false)
            found_potentials = false
            n_ = length(nodes_)
            m_ = length(plateau_edges[key])
            for (j, node) in enumerate(nodes_)
                old_to_new_id[node] = j
                if is_fixed[find(uf, node)]
                    seeds_[j, :] .= potentials[find(uf, node), :]
                    is_fixed_[j] = true
                    found_potentials = true
                end
            end

            if !found_potentials
                # no fixed nodes found in this plateau, just merge the nodes
                for node in nodes_
                    if union!(uf, node, nodes_[1])
                        num_unfixed_clusters -= 1
                    end
                end
                continue
            end

            # Create a list of the edges themselves from the list of edge indices.
            # Additionally, we use old_to_new_id to transform the node identifiers
            # (we want edge entries to be wrt the local node ids for the plateau).
            for (j, edge_idx) in enumerate(plateau_edges[key])
                edge_idx += i - 1
                edges_[1, j] = old_to_new_id[find(uf, edges[1, edge_idx])]
                edges_[2, j] = old_to_new_id[find(uf, edges[2, edge_idx])]
            end

            # Now we construct the Laplacian for the plateau.
            # We need to be careful about the fact that we are working
            # with a multi-graph: there might be self-loops (which we ignore)
            # and multiple edges between the same nodes (in which case we add up the weights)

            # if the plateau is relatively small, we just construct and solve
            # a dense system
            if n_ <= 500
                lap = zeros(n_, n_)
                for j in 1:m_
                    # ignore self-loops
                    if edges_[1, j] == edges_[2, j]
                        continue
                    end
                    # Fill the Laplacian with the negative adjacency matrix.
                    # Note that there may be multiple edges between two nodes,
                    # so we need -= 1 instead of = -1. We can use 1 instead of
                    # the actual weight because all edges inside the plateau
                    # have the same weight.
                    lap[edges_[1, j], edges_[2, j]] -= 1
                    lap[edges_[2, j], edges_[1, j]] -= 1
                end
                # add the diagonal, which contains the degrees of each node
                degrees = -sum(lap, dims=1)
                for j in 1:length(nodes_)
                    lap[j, j] = degrees[j]
                end

            # if the plateau is large, it is likely very sparse and we
            # use a sparse Laplacian matrix instead
            else
                row_indices = Vector{Int}()
                col_indices = Vector{Int}()
                values = Vector{Float64}()
                for j in 1:m_
                    # ignore self-loops
                    if edges_[1, j] == edges_[2, j]
                        continue
                    end
                    # we don't need to worry about multi-edges,
                    # they will be combined anyway when forming the sparse
                    # array
                    push!(row_indices, edges_[1, j])
                    push!(col_indices, edges_[2, j])
                    push!(values, -1)

                    # also add the edge in the other direction
                    push!(row_indices, edges_[2, j])
                    push!(col_indices, edges_[1, j])
                    push!(values, -1)

                    # Finally, we add a 1 on the diagonal for both nodes
                    # (again, they are added up later automatically)
                    push!(row_indices, edges_[1, j])
                    push!(col_indices, edges_[1, j])
                    push!(values, 1)

                    push!(row_indices, edges_[2, j])
                    push!(col_indices, edges_[2, j])
                    push!(values, 1)
                end
                lap = sparse(row_indices, col_indices, values, length(nodes_), length(nodes_))
            end

            # then, we construct the linear system
            # TODO: this indexing is slow for sparse matrices,
            # is there a way to do this already during construction?
            # is_fixed_ is too long, we only need the first n_ entries
            is_fixed_view = is_fixed_[1:n_]
            lap_unseeded = lap[.!is_fixed_view, .!is_fixed_view]
            BT = lap[.!is_fixed_view, is_fixed_view]
            # seeds_ is just as long as is_fixed, so use everything here
            b = - BT * seeds_[is_fixed_, :]
            # TODO: should we tell Julia that the Laplacian is symmetric
            # and positive definite, to improve performance?
            x = lap_unseeded \ b

            # We now have the solution for the unseeded nodes of the
            # current plateau in x. We just need to fill them into
            # the right places.
            # k is the index for x, i.e. it enumerates the unseeded nodes
            # in the plateau. j enumerates all nodes in the plateau.
            # node gives the global index of the node j in the plateau.
            k = 1
            for (j, node) in enumerate(nodes_)
                if is_fixed_[j]
                    continue
                end
                potentials[find(uf, node), :] .= x[k, :]
                is_fixed[find(uf, node)] = true
                num_unfixed_clusters -= 1
                k += 1
            end
        end
        # now we skip the remaining edges in the plateau
        i += num_edges - 1
    end

    for j in 1:g.n
        if !is_fixed[find(uf, j)]
            error("Something went wrong, not all nodes are fixed yet")
        end
        potentials[j, :] .= potentials[find(uf, j), :]
    end

    return Dict(seed => potentials[:, seed] for seed in seed_list)
end

export Graph, karger, sample_cuts, potential, watershed, power_watershed, power_watershed_multi
end


text:
putAdj(A::Assoc,AA::AbstractSparseMatrix) =  Assoc(copy(A.row),copy(A.col),copy(A.val),AA)
putAdj(A::Assoc,AA::Array{Int,2}) =  putAdj(A::Assoc,sparse(AA))

StringOrNum = Union{AbstractString,Number}
UnionArray = Array{Union{AbstractString,Number}}

# put for new array
putRow(A::Assoc,ARow::UnionArray) =  Assoc(ARow,copy(A.col),copy(A.val),copy(A.A))
putCol(A::Assoc,ACol::UnionArray) =  Assoc(copy(A.row),ACol,copy(A.val),copy(A.A))
putVal(A::Assoc,AVal::UnionArray) =  Assoc(copy(A.row),copy(A.col),AVal,copy(A.A))

# put for new scalar
function putRow(A::Assoc,ARow::StringOrNum)
    r,c,v = find(A)
    Assoc(ARow,c,v)
end

function putCol(A::Assoc,ACol::StringOrNum)
    r,c,v = find(A)
    Assoc(r,ACol,v)
end

function putVal(A::Assoc,AVal::StringOrNum)
    r,c,v = find(A)
    Assoc(r,c,AVal)
end

#=
nocol : remove columns mapping
=#

function nocol(A::Assoc)
    return condense(Assoc(A.row,promote(ones(size(A.col)),A.col)[1],A.val,A.A))
end

#=
norow : remove rows mapping
=#

function norow(A::Assoc)
    return condense(Assoc(promote(ones(size(A.row)),A.row)[1],A.col,A.val,A.A))
end


text:
module MachineLearning

end


text:
module Inputs

include("PWscf.jl")
# include("CP.jl")

end


text:
module SemiLagrangian

# import Base.Threads: @sync, @spawn, nthreads, threadid

using Polynomials: StandardBasisPolynomial
using Polynomials
using FFTW
using Base.Threads

# include("mpoly/MultiPoly.jl")
include("mpiinterface.jl")

include("util.jl")
include("cplxlagrange.jl")
include("fftbig.jl")
include("mesh.jl")
include("interpolation.jl")
include("lagrange.jl")
include("hermite.jl")
include("spline.jl")
include("bspline.jl")
include("bsplinelu.jl")
include("bsplinefft.jl")
include("advection.jl")
include("util_poisson.jl")
include("poisson.jl")
include("rotation.jl")
include("translation.jl")
include("quasigeostrophic.jl")

export UniformMesh, start, stop, AbstractInterpolation, get_order
export Advection, AdvectionData
export nosplit, standardsplit, strangsplit, triplejumpsplit, order6split, hamsplit_3_11, ymsplit, table2split
export Lagrange, Hermite, B_SplineLU, B_SplineFFT, interpolate!
export compute_charge!,
    compute_elfield!, compute_elfield, compute_ee, compute_ke, advection!
export dotprod, getpoissonvar, getrotationvar, gettranslationvar
export TimeOptimization, NoTimeOpt, SimpleThreadsOpt, SplitThreadsOpt, MPIOpt, TimeAlgorithm, NoTimeAlg, ABTimeAlg_ip, ABTimeAlg_init, ABTimeAlg_new
export TypePoisson, StdPoisson, StdPoisson2d, StdOrder2_1, StdOrder2_2, StdAB, StdAB2, StdRK4, StdABinit, StdABp
export get_type, sizeall, getdata, OpTuple, getgeovar, initdata!, getenergyall

end


text:
@testset "gc" begin

try
    throw(nothing)
catch
    # clear exception_in_transit, which can hold a reference to objects
    # (this works around JuliaLang/julia#20784)
end

# test outstanding contexts
@test length(CUDAdrv.context_instances) == 1
destroy!(ctx)
for i in 1:50
    GC.gc()
end
if length(CUDAdrv.context_instances) > 0
    for (handle, object) in CUDAdrv.context_instances
        warn("CUDA context $handle has outstanding object at $object")
    end
end
@test length(CUDAdrv.context_instances) == 0

end


text:
function _unitrange(start,stop)
  s = 0.0
  for i = start:stop
      s += 2.1
  end
  return s
end

function _unitrange2(start,stop)
  s = 0.0
  for i = start:stop
      s = s + 2.1 + float(i)
  end
  return s
end

function _steprange(start,step,stop)
  s = 0.0
  for i = start:step:stop
      s += 2.1
  end
  return s
end

function _while(n)
  s = 0.0
  i=1
  while i<n
      s += 2.1
      i+=1
  end
  return s
end

function exports()
  _unitrange(1,2)
  _unitrange2(1,2)
  _steprange(1,2,3)
  _while(3)
end

exports()


text:
struct ELM{T₁ <: Number, T₂ <: ActivationFunction}
    input_weights::Matrix{T₁}
    output_weights::Matrix{T₁}
    activation_functions::Vector{T₂}

    function ELM{T₁}(samples::T₂,
                     targets::T₃,
                     input_weights::T₄,
                     activation_functions::Vector{T₅};
                     batch_size::Int = 1000) where {T₁ <: Number,
                                                    T₂ <: AbstractMatrix,
                                                    T₃ <: AbstractVector,
                                                    T₄ <: AbstractMatrix,
                                                    T₅ <: ActivationFunction}
        X = samples
        t = targets
        W = input_weights
        L = first(size(W))
        fs = activation_functions
        B = train_output_weights(T₁, X, t, W, fs, batch_size)
        new{T₁,eltype(fs)}(W, B, fs)
    end

    function ELM{T₁}(samples::T₂,
                     targets::T₃,
                     n_neurons::Int,
                     activation_function::T₄;
                     batch_size::Int = 1000) where {T₁ <: Number,
                                                    T₂ <: AbstractMatrix,
                                                    T₃ <: AbstractVector,
                                                    T₄ <: NaiveActivationFunction}
        X = samples
        t = targets
        L = n_neurons
        f = activation_function
        D = first(size(X))
        W = gaussian_projection_matrix(T₁, L, D)
        ELM{T₁}(X, t, W, f, batch_size = batch_size)
    end
end

function train_output_weights(::Type{T₁},
                              samples::T₂,
                              targets::T₃,
                              input_weights::T₄,
                              activation_functions::Vector{T₅},
                              batch_size::Int) where {T₁ <: Number,
                                                      T₂ <: AbstractMatrix,
                                                      T₃ <: AbstractVector,
                                                      T₄ <: AbstractMatrix,
                                                      T₅ <: ActivationFunction}
    X = samples
    t = targets
    ψ = calculate_sample_weights(t)
    W = input_weights
    fs = activation_functions
    L = length(fs)
    D, N = size(X)

    HH = zeros(T₁, L + 1, L + 1)
    TH = zeros(T₁, 1, L + 1)

    for batch in partition_range(1:N, batch_size)
        X₀ = @view X[:,batch]
        t₀ = @view t[batch]
        ψ₀ = @view ψ[batch]
        train_on_batch!(HH, TH, X₀, t₀, ψ₀, W, fs)
    end

    output_weights = TH * LinearAlgebra.pinv(HH)
end

function train_on_batch!(HH::T₁,
                         TH::T₂,
                         samples::T₃,
                         targets::T₄,
                         sample_weights::T₅,
                         input_weights::T₆,
                         activation_functions::Vector{T₇}) where {T₁ <: AbstractMatrix,
                                                                  T₂ <: AbstractMatrix,
                                                                  T₃ <: AbstractMatrix,
                                                                  T₄ <: AbstractVector,
                                                                  T₅ <: AbstractVector,
                                                                  T₆ <: AbstractMatrix,
                                                                  T₇ <: ActivationFunction}
    X = samples
    T = reshape(targets, 1, :)
    Ψ = LinearAlgebra.Diagonal(sample_weights)
    W = input_weights
    fs = activation_functions
    N = last(size(X))

    H = project(X, W, fs)
    H = [H; ones(eltype(H), 1, N)]

    H = H * Ψ
    T = T * Ψ

    HH .+= (H * H')
    TH .+= (T * H')

    HH, TH
end

function predict(elm::T₁,
                 samples::T₂) where {T₁ <: ELM,
                                     T₂ <: AbstractMatrix}
    X = samples
    W = elm.input_weights
    B = elm.output_weights
    fs = elm.activation_functions
    N = last(size(X))

    H = project(X, W, fs)
    H = [H; ones(eltype(H), 1, N)]

    y = vec(B * H)
end

function predict(elm::T₁,
                 sample::T₂) where {T₁ <: ELM,
                                    T₂ <: AbstractVector}
    y = predict(elm, reshape(sample, :, 1))
    first(y)
end


text:
#ДАНО: Робот - рядом с горизонтальной перегородкой (под ней),
#бесконечно продолжающейся в обе стороны, в которой имеется
#проход шириной в одну клетку.
#РЕЗУЛЬТАТ: Робот - в клетке под проходом
function find_passage(r)
    n=0; side=Ost
    while isborder(r,Nord)==true # прохода сверху нет
        n+=1
        move!(r,side,n)
        side=inverse(side)
    end
end

move!(r,side,num_steps)=for _ in 1:num_steps move!(r,side) end

inverse(side::HorizonSide)=HorizonSide(mod(Int(side)+2,4))


text:
# This file defines all of the mlpack types.

struct AdaBoostModel
  ptr::Ptr{Nothing}
end

struct ApproxKFNModel
  ptr::Ptr{Nothing}
end

struct BayesianLinearRegression
  ptr::Ptr{Nothing}
end

struct CFModel
  ptr::Ptr{Nothing}
end

struct DSModel
  ptr::Ptr{Nothing}
end

struct DecisionTreeModel
  ptr::Ptr{Nothing}
end

struct DTree
  ptr::Ptr{Nothing}
end

struct FastMKSModel
  ptr::Ptr{Nothing}
end

struct GMM
  ptr::Ptr{Nothing}
end

struct HMMModel
  ptr::Ptr{Nothing}
end

struct HoeffdingTreeModel
  ptr::Ptr{Nothing}
end

struct KDEModel
  ptr::Ptr{Nothing}
end

struct LARS
  ptr::Ptr{Nothing}
end

struct LinearRegression
  ptr::Ptr{Nothing}
end

struct LinearSVMModel
  ptr::Ptr{Nothing}
end

struct LocalCoordinateCoding
  ptr::Ptr{Nothing}
end

struct LogisticRegression
  ptr::Ptr{Nothing}
end

struct LSHSearch
  ptr::Ptr{Nothing}
end

struct NBCModel
  ptr::Ptr{Nothing}
end

struct KNNModel
  ptr::Ptr{Nothing}
end

struct KFNModel
  ptr::Ptr{Nothing}
end

struct PerceptronModel
  ptr::Ptr{Nothing}
end

struct ScalingModel
  ptr::Ptr{Nothing}
end

struct RandomForestModel
  ptr::Ptr{Nothing}
end

struct RANNModel
  ptr::Ptr{Nothing}
end

struct SoftmaxRegression
  ptr::Ptr{Nothing}
end

struct SparseCoding
  ptr::Ptr{Nothing}
end


text:
using Test

include("../src/runtests.jl")


text:
@testset "inplace" begin
    # R^N → R
    f1 = OnceDifferentiable(exponential, rand(2), 0.0)
    fia1 = OnceDifferentiable(exponential, rand(2); inplace = false)
    fi1 = OnceDifferentiable(exponential, exponential_gradient, rand(2); inplace = false)
    fi2 = OnceDifferentiable(exponential, exponential_gradient, rand(2); inplace = false)
    fi3 = OnceDifferentiable(exponential, exponential_gradient, exponential_value_gradient,
                       rand(2); inplace = false)
    xr = rand(2)
    @test value!(f1, xr) ≈ value!(fia1, xr) ≈ value!(fi1, xr) ≈ value!(fi2, xr) ≈ value!(fi3, xr)
    @test gradient!(f1, xr) ≈ gradient!(fia1, xr)
    @test gradient!(fia1, xr) ≈ gradient!(fi1, xr) ≈ gradient!(fi2, xr) ≈ gradient!(fi3, xr)
    vg1 = value_gradient!(fia1, xr)
    vg2 = value_gradient!(fi1, xr)
    vg3 = value_gradient!(fi2, xr)
    vg4 = value_gradient!(fi3, xr)
    @test vg1[1] ≈ vg2[1] ≈ vg3[1] ≈ vg4[1]
    @test vg1[2] ≈ vg2[2] ≈ vg3[2] ≈ vg4[2]

    ft1 = TwiceDifferentiable(exponential, rand(2), 0.0)
    ftia1 = TwiceDifferentiable(exponential, rand(2); inplace = false)
    fti1 = TwiceDifferentiable(exponential, exponential_gradient, rand(2); inplace = false)
    fti2 = TwiceDifferentiable(exponential, exponential_gradient, rand(2); inplace = false)
    fti3 = TwiceDifferentiable(exponential, exponential_gradient, exponential_hessian,
                               rand(2); inplace = false)
    fti4 = TwiceDifferentiable(exponential, exponential_gradient, exponential_value_gradient,
                               exponential_hessian, rand(2); inplace = false)

    @test value!(ft1, xr) ≈ value!(ftia1, xr) ≈ value!(fti1, xr) ≈ value!(fti2, xr) ≈ value!(fti3, xr) ≈ value!(fti4, xr)
    @test gradient!(ft1, xr) ≈ gradient!(ftia1, xr)
    @test gradient!(ftia1, xr) ≈ gradient!(fti1, xr) ≈ gradient!(fti2, xr) ≈ gradient!(fti3, xr) ≈ gradient!(fti4, xr)
    vg1 = value_gradient!(ftia1, xr)
    vg2 = value_gradient!(fti1, xr)
    vg3 = value_gradient!(fti2, xr)
    vg4 = value_gradient!(fti3, xr)
    vg5 = value_gradient!(fti4, xr)
    @test vg1[1] ≈ vg2[1] ≈ vg3[1] ≈ vg4[1] ≈ vg5[1]
    @test vg1[2] ≈ vg2[2] ≈ vg3[2] ≈ vg4[2] ≈ vg5[2]
    @test hessian!(ft1, xr) ≈ hessian!(ftia1, xr)
    @test hessian!(fti1, xr) ≈ hessian!(fti2, xr)
    @test hessian!(fti3, xr) ≈ hessian!(fti4, xr)

    # R^N → R^N
    f1 = OnceDifferentiable(exponential_gradient!, rand(2), rand(2))
    fia1 = OnceDifferentiable(exponential_gradient, rand(2), rand(2); inplace = false)
    fia2 = OnceDifferentiable(exponential_gradient, rand(2), rand(2); inplace = false, autodiff = :forward)
    fi1 = OnceDifferentiable(exponential_gradient, exponential_hessian, rand(2), rand(2); inplace = false)
    fi2 = OnceDifferentiable(exponential_gradient, exponential_hessian, rand(2), rand(2); inplace = false)
    fi3 = OnceDifferentiable(exponential_gradient, exponential_hessian, exponential_gradient_hessian,
    rand(2), rand(2); inplace = false)
    xr = fill(2.0, 2)
    @test value!(f1, xr) ≈ value!(fia1, xr) ≈ value!(fia2, xr) ≈ value!(fi1, xr) ≈ value!(fi2, xr) ≈ value!(fi3, xr)
    @test jacobian!(f1, xr) ≈ jacobian!(fia1, xr) ≈ jacobian!(fia2, xr) ≈ jacobian!(fi1, xr) ≈ jacobian!(fi2, xr) ≈ jacobian!(fi3, xr)
    vj1 = value_jacobian!(f1, xr)
    vj2 = value_jacobian!(fia1, xr)
    vj3 = value_jacobian!(fia2, xr)
    vj4 = value_jacobian!(fi1, xr)
    vj5 = value_jacobian!(fi2, xr)
    vj6 = value_jacobian!(fi3, xr)
    @test vj1[1] ≈ vj2[1] ≈ vj3[1] ≈ vj4[1] ≈ vj5[1] ≈ vj6[1]
    @test vj1[2] ≈ vj2[2] ≈ vj3[2] ≈ vj4[2] ≈ vj5[2] ≈ vj6[2]
    # Finite doesn't cut it here, so had to split it up
    xr = [1.0, 3.4]
    @test value!(f1, xr) ≈ value!(fia1, xr) ≈ value!(fi1, xr) ≈ value!(fi2, xr) ≈ value!(fi3, xr)
    # need to split this up as the example is hard to approximate with finite differences
    @test jacobian!(f1, xr) ≈ jacobian!(fia1, xr) ≈ jacobian!(fia2, xr)
    @test jacobian!(fi1, xr) ≈ jacobian!(fi2, xr) ≈ jacobian!(fi3, xr)
end

@testset "autodiff" begin
    # base line api
    f1 = OnceDifferentiable(exponential, exponential_gradient!, exponential_value_gradient!,
                            rand(2), 0.0, rand(2))
    # default autodiff R^N → R
    fa1 = OnceDifferentiable(exponential, rand(2))
    # specific autodiff R^N → R
    fa2 = OnceDifferentiable(exponential, rand(2); autodiff = :finite)
    fa3 = OnceDifferentiable(exponential, rand(2); autodiff = :forward)
    # random input
    xr = rand(2)
    # test that values are the same approximately
    @test value!(f1, xr) ≈ value!(fa1, xr) ≈ value!(fa2, xr) ≈ value!(fa3, xr)
    @test gradient!(f1, xr) ≈ gradient!(fa1, xr) ≈ gradient!(fa2, xr) ≈ gradient!(fa3, xr)
end

@testset "inplace/autodiff" begin
    # R^N → R^N
    fa1 = OnceDifferentiable(exponential_gradient!, rand(2), rand(2); inplace = true)
    fia1 = OnceDifferentiable(exponential_gradient, rand(2), rand(2); inplace = false)
    fia2 = OnceDifferentiable(exponential_gradient, rand(2), rand(2); inplace = false, autodiff = :finite)
    fia3 = OnceDifferentiable(exponential_gradient, rand(2), rand(2); inplace = false, autodiff = :forward)
    xr = rand(2)
    @test value!(fa1, xr) ≈ value!(fia1, xr) ≈ value!(fia2, xr) ≈ value!(fia3, xr)
    @test jacobian!(fa1, xr) ≈ jacobian!(fia1, xr) ≈ jacobian!(fia2, xr) ≈ jacobian!(fia3, xr)
end


text:
using CSV, DataFrames, Statistics

results_dir = ARGS[1]
results = CSV.read("$results_dir/results.csv")

by(results, [:sampler, :Nd], :a0_ess=>mean, :a1_ess=>mean, :a0_sig_ess=>mean, :tree_depth=>mean, :epsilon=>mean) |> println


text:
module computational_funcs
# some computational functions for calculating stuff

using DifferentialEquations:ODEProblem, solve, Tsit5
using LinearAlgebra:Matrix, I, mul!
using StaticArrays

export propagator, full_propagator, fast_propagator, full_fast_propagator, expectation_value, track_expectation_value

# in place von Neumann equation
function comp_U!(du, u, p, t)
    mul!(du, -1.0im * p.func(p, t), u)
end

# von Neumann equation for use with static arrays
function comp_U(u, p, t)
    -1.0im * p.func(p, t) * u
end

# computes the closed system propagator U
function propagator(p, td, u0::T)::T where T
    # p is a named tuple containing at least a field called func which computes the Hamiltonian
    tspan = (0.0, td)
    prob = ODEProblem(comp_U, u0, tspan, p)
    sol = solve(prob, abstol = 1e-6, reltol = 1e-6, save_everystep=false, alg = Tsit5())
    return sol.u[end]
end

# computes the propagator but returns the ODE solution object
function full_propagator(p, td, u0)
    tspan = (0.0, td)
    prob = ODEProblem(comp_U, u0, tspan, p)
    sol = solve(prob, abstol = 1e-6, reltol = 1e-6, save_everystep=true, alg = Tsit5())
    return sol
end

# this solves the ODE much faster *but* it doesn't specify the tolerances so you might see some deviation
function fast_propagator(p, td, u0::T)::T where T
    tspan = (0.0, td)
    prob = ODEProblem(comp_U, u0, tspan, p)
    sol = solve(prob, save_everystep=false, alg = Tsit5())
    return sol.u[end]
end

function full_fast_propagator(p, td, u0)
    tspan = (0.0, td)
    prob = ODEProblem(comp_U, u0, tspan, p)
    sol = solve(prob, save_everystep=true, alg = Tsit5())
    return sol
end

# computes the expectation value of an operator and a state psi
function expectation_value(psi, operator)::AbstractFloat
    real(psi' * operator * psi)
end

# track the expectation value of an operator as a function of time
function track_expectation_value(sol, psi, operator, res = 1000)
    z = @MVector zeros(res)
    for i = 1:res
        t = (i-1)/res * sol.t[end]
        @inbounds z[i] = expectation_value(sol(t) * psi, operator)
    end
    return z
end

end


text:
@testset "print" begin

println(Seq())
println(e"a" + E"b" | p"c" > tuple)

end


text:
struct FormatLogger <: AbstractLogger
    f::Function
    stream::IO
    always_flush::Bool
end

"""
    FormatLogger(f::Function, io::IO=stderr; always_flush=true)

Logger sink that formats the message and finally writes to `io`.
The formatting function should be of the form `f(io::IOContext, log_args::NamedTuple)`
where `log_args` has the following fields:
`(level, message, _module, group, id, file, line, kwargs)`.
See [`LoggingExtras.handle_message_args`](@ref) for more information on what field is.

# Examples
```julia-repl
julia> using Logging, LoggingExtras

julia> logger = FormatLogger() do io, args
           println(io, args._module, " | ", "[", args.level, "] ", args.message)
       end;

julia> with_logger(logger) do
           @info "This is an informational message."
           @warn "This is a warning, should take a look."
       end
Main | [Info] This is an informational message.
Main | [Warn] This is a warning, should take a look.
```
"""
function FormatLogger(f::Function, io::IO=stderr; always_flush=true)
    return FormatLogger(f, io, always_flush)
end

function handle_message(logger::FormatLogger, args...; kwargs...)
    log_args = handle_message_args(args...; kwargs...)
    # We help the user by passing an IOBuffer to the formatting function
    # to make sure that everything writes to the logger io in one go.
    iob = IOBuffer()
    ioc = IOContext(iob, logger.stream)
    logger.f(ioc, log_args)
    write(logger.stream, take!(iob))
    logger.always_flush && flush(logger.stream)
    return nothing
end
shouldlog(logger::FormatLogger, arg...) = true
min_enabled_level(logger::FormatLogger) = BelowMinLevel
catch_exceptions(logger::FormatLogger) = true # Or false? SimpleLogger doesn't, ConsoleLogger does.


text:
using Test
using TimeseriesPrediction
using StaticArrays
using LinearAlgebra
import Statistics:var
println("\nTesting local models...")

import Random
Random.seed!(42)

ds = Systems.roessler()
data = trajectory(ds, 5000; dt=0.1)
N_train = 45000
s_train = data[1:N_train, 1]
s_test  = data[N_train:end,1]

@testset "ALM localmodel_tsp" begin
    @testset "γ=$γ and τ=$τ" for γ ∈ [2,3], τ ∈ [14,15]
        p = 50
        method = AverageLocalModel()
        ntype = FixedMassNeighborhood(2)
        stepsize = 1
        s_pred = localmodel_tsp(s_train, γ, τ, p;
        method=method, ntype=ntype, stepsize=stepsize)
        @test length(s_pred) == p+1
        @test norm(s_test[1:p+1] - s_pred)/p < 5e-2

        #Repeat with reconstruction given
        R = reconstruct(s_train, γ, τ)
        s_pred = localmodel_tsp(R, p; method=method, ntype=ntype, stepsize=stepsize)[:,γ+1]
        @test length(s_pred) == p+1
        @test norm(s_test[1:p+1] - s_pred)/p < 5e-2
    end
end

@testset "LinearLocalModel Ridge Reg" begin
    @testset "γ=$γ and τ=$τ" for γ ∈ [3,4], τ ∈ [14,15]
        p = 50
        method = LinearLocalModel()
        ntype = FixedMassNeighborhood(3)
        stepsize = 1
        s_pred = localmodel_tsp(s_train, γ, τ, p;
        method=method, ntype=ntype, stepsize=stepsize)
        @test length(s_pred) == p+1
        @test norm(s_test[1:p+1] - s_pred)/p < 6e-2

        #Repeat with reconstruction given
        R = reconstruct(s_train, γ, τ)
        s_pred = localmodel_tsp(R, p; method=method, ntype=ntype, stepsize=stepsize)[:,γ+1]
        @test length(s_pred) == p+1
        @test norm(s_test[1:p+1] - s_pred)/p < 6e-2
    end
end
@testset "LinearLocalModel McNames Reg" begin
    @testset "γ=$γ and τ=$τ" for γ ∈ [3,4], τ ∈ [14,15]
        p = 50
        method = LinearLocalModel(TimeseriesPrediction.ω_safe, 0.1,1.)
        ntype = FixedMassNeighborhood(3)
        stepsize = 1
        s_pred = localmodel_tsp(s_train, γ, τ, p;
        method=method, ntype=ntype, stepsize=stepsize)
        @test length(s_pred) == p+1
        @test norm(s_test[1:p+1] - s_pred)/p < 5e-2

        #Repeat with reconstruction given
        R = reconstruct(s_train, γ, τ)
        s_pred = localmodel_tsp(R, p; method=method, ntype=ntype, stepsize=stepsize)[:,γ+1]
        @test length(s_pred) == p+1
        @test norm(s_test[1:p+1] - s_pred)/p < 5e-2
    end
end

@testset "FixedSizeNeighborhood" begin
    @testset "γ=$γ and τ=$τ" for γ ∈ [4], τ ∈ [14,15]
        p = 25
        method = AverageLocalModel()
        ntype = FixedSizeNeighborhood(0.5)
        stepsize = 1
        s_pred = localmodel_tsp(s_train, γ, τ, p;
        method=method, ntype=ntype, stepsize=stepsize)
        @test length(s_pred) == p+1
        @test norm(s_test[1:p+1] - s_pred)/p < 0.2
    end
end

@testset "Multivariate Input predict" begin
    sm_train = data[1:N_train,SVector(1,2)]
    @testset "γ=$γ and τ=$τ" for γ ∈ [3,4], τ=15
        R = reconstruct(sm_train,γ,τ)
        num_points = 50
        method = AverageLocalModel()
        ntype = FixedMassNeighborhood(2)
        stepsize = 1
        sind = SVector{2, Int}(((γ+1)*2 - i for i in 2-1:-1:0)...)
        pred = localmodel_tsp(R,num_points; method=method,ntype=ntype,stepsize=stepsize)[:,sind]
        @test size(pred) == (num_points+1, 2)
        @test norm(s_test[1:num_points+1] - pred[:, 1])/num_points < 5e-2
        @test norm(data[N_train:N_train+num_points, 2] - pred[:, 2])/num_points < 5e-2

        pred2 = localmodel_tsp(sm_train, γ,τ,num_points; method=method,ntype=ntype,stepsize=stepsize)
        @test norm(pred.data-pred2.data) < 1e-10
    end
    @testset "γ=3 and multi τ" begin
    sm_train = data[1:N_train,SVector(1,2)]

        γ = 3;
        τ=[15 15; 30 29; 45 45]#[14 15; 29 30; 45 47]
        R = reconstruct(sm_train,γ,τ)
        num_points = 25
        method = AverageLocalModel()
        ntype = FixedMassNeighborhood(2)
        stepsize = 1
        svind = SVector{2, Int}(7,8)
        pred = localmodel_tsp(R,num_points; method=method,ntype=ntype,stepsize=stepsize)[:,svind]
        @test size(pred) == (num_points+1, 2)
        @test norm(s_test[1:num_points+1] - pred[:, 1])/num_points < 1e-1
        @test norm(data[N_train:N_train+num_points, 2] - pred[:, 2])/num_points < 1e-1

        pred2 = localmodel_tsp(sm_train,γ,τ,num_points; method=method,ntype=ntype,stepsize=stepsize)
        @test norm(pred.data-pred2.data) < 1e-10
    end
end

@testset "MSE" begin
    @testset "p=$p" for p ∈ [50,100]
        γ = 3; τ = 15;
        R = reconstruct(s_train,γ,τ)
        R_test = reconstruct(s_test[end-γ*τ-p-50:end],γ,τ)
        method = AverageLocalModel()
        ntype = FixedMassNeighborhood(2)
        stepsize = 1

        @test MSEp(R,R_test,p; method=method,ntype=ntype,stepsize=stepsize)/p < 5e-2
        @test MSEp(R,R_test,1; method=method,ntype=ntype,stepsize=stepsize) < 5e-2
    end
end

@testset "Maps" begin
    ds = Systems.standardmap()
    data = trajectory(ds,50000; dt=1)
    N_train = 49900
    p = 25
    stepsize = 1
    s = data[1:N_train,1]
    s_real = data[N_train:N_train+p*stepsize,1]
    method = AverageLocalModel()
    ntype = FixedMassNeighborhood(2)
    s_pred = localmodel_tsp(s,2,1,p; method=method, ntype=ntype)
    @test norm(s_pred -s_real)/var(s)/p < 0.2
end

#####################################################################################
#                                  Cross Prediction                                   #
#####################################################################################

ds = Systems.roessler()
data = trajectory(ds, 5000; dt=0.1)
N_train = 45000
N_test = 1000
source_train = data[1:N_train, 1]
target_train = data[1:N_train, 2]
source_pred = data[N_train+1:N_train+N_test, 1]
target_test = data[N_train+1:N_train+N_test, 2]

@testset "ALM localmodel_cp" begin
    @testset "γ=$γ and τ=$τ" for γ ∈ [3,4], τ ∈ [14,15]
        method = AverageLocalModel()
        ntype = FixedMassNeighborhood(2)
        target_pred = localmodel_cp(    source_train,
                                        target_train,
                                        source_pred, γ, τ;
                                        method=method, ntype=ntype)
        @test length(target_pred) == N_test - γ*τ
        @test norm(target_test[1+γ*τ:end] .- target_pred)/N_test < 5e-3

        #Repeat with reconstruction given
        R      = reconstruct(source_train, γ, τ)
        R_pred = reconstruct(source_pred, γ, τ)
        target_pred = localmodel_cp(
                        R, target_train, R_pred;
                        method=method, ntype=ntype, y_idx_shift=γ*τ)
        @test length(target_pred) == N_test - γ*τ
        @test norm(target_test[1+γ*τ:end] .- target_pred)/N_test < 5e-3
    end
end

@testset "LinearLocalModel Ridge Reg" begin
    @testset "γ=$γ and τ=$τ" for γ ∈ [3,4], τ ∈ [14,15]
        method = LinearLocalModel()
        ntype = FixedMassNeighborhood(3)
        target_pred = localmodel_cp(    source_train,
                                        target_train,
                                        source_pred, γ, τ;
                                        method=method,
                                        ntype=ntype)
        @test length(target_pred) == N_test - γ*τ
        @test norm(target_test[1+γ*τ:end] .- target_pred)/N_test < 5e-3

        #Repeat with reconstruction given
        R      = reconstruct(source_train, γ, τ)
        R_pred = reconstruct(source_pred, γ, τ)
        target_pred = localmodel_cp(
                        R, target_train, R_pred;
                        method=method, ntype=ntype, y_idx_shift=γ*τ)
        @test length(target_pred) == N_test - γ*τ
        @test norm(target_test[1+γ*τ:end] .- target_pred)/N_test < 5e-3
    end
end
@testset "LinearLocalModel McNames Reg" begin
    @testset "γ=$γ and τ=$τ" for γ ∈ [3,4], τ ∈ [14,15]
        p = 50
        method = LinearLocalModel(TimeseriesPrediction.ω_safe, 0.1,1.)
        ntype = FixedMassNeighborhood(3)
        stepsize = 1
        target_pred = localmodel_cp(    source_train,
                                        target_train,
                                        source_pred, γ, τ;
                                        method=method,
                                        ntype=ntype)
        @test length(target_pred) == N_test - γ*τ
        @test norm(target_test[1+γ*τ:end] .- target_pred)/N_test < 5e-3

        #Repeat with reconstruction given
        R      = reconstruct(source_train, γ, τ)
        R_pred = reconstruct(source_pred, γ, τ)
        target_pred = localmodel_cp(
                        R, target_train, R_pred;
                        method=method, ntype=ntype, y_idx_shift=γ*τ)
        @test length(target_pred) == N_test - γ*τ
        @test norm(target_test[1+γ*τ:end] .- target_pred)/N_test < 5e-3
    end
end

@testset "FixedSizeNeighborhood" begin
    @testset "γ=$γ and τ=$τ" for γ ∈ [4], τ ∈ [14,15]
        method = AverageLocalModel()
        ntype = FixedSizeNeighborhood(0.5)
        target_pred = localmodel_cp(    source_train,
                                        target_train,
                                        source_pred, γ, τ;
                                        method=method,
                                        ntype=ntype)
        @test length(target_pred) == N_test - γ*τ
        @test norm(target_test[1+γ*τ:end] .- target_pred)/N_test < 5e-3
    end
end

@testset "Multivariate Input predict" begin
    sm_train = data[1:N_train,SVector(1,2)]
    sm_pred = data[1+N_train:N_train+N_test,SVector(1,2)]
    @testset "γ=$γ and τ=$τ" for γ ∈ [3,4], τ=15
        source_train = reconstruct(sm_train,γ,τ)
        source_pred = reconstruct(sm_pred,γ,τ)
        method = AverageLocalModel()
        ntype = FixedMassNeighborhood(2)
        target_pred = localmodel_cp(    source_train,
                                        target_train,
                                        source_pred;
                                        method=method,
                                        ntype=ntype,
                                        y_idx_shift=γ*τ)
        @test length(target_pred) == N_test - γ*τ
        @test norm(target_test[1+γ*τ:end] .- target_pred)/N_test < 5e-3
    end
    @testset "γ=3 and multi τ" begin
        γ = 3;
        τ=[15 15; 30 29; 45 45]
        source_train = reconstruct(sm_train,γ,τ)
        source_pred = reconstruct(sm_pred,γ,τ)
        method = AverageLocalModel()
        ntype = FixedMassNeighborhood(2)
        target_pred = localmodel_cp(    source_train,
                                        target_train,
                                        source_pred;
                                        method=method,
                                        ntype=ntype,
                                        y_idx_shift=maximum(τ))
        @test length(target_pred) == N_test - maximum(τ)
        @test norm(target_test[1+maximum(τ):end] .- target_pred)/N_test < 5e-3
    end
end


text:
using IsoPkg
using Test,Pkg

group_tmp="temp_"*basename(tempname())
IsoPkg.switch(group_tmp) #Switch the current project group to a temp empty group

IsoPkg.add("Glob") #Install Glob
IsoPkg.add("Glob@1.2.0") #Install Glob v1.2.0 and pin the version

#Notice: If the following `using` statements are in the same session, only the first one does the real loading.
@iso using Glob #Load Glob
@iso using Glob "1.2.0" #Load Glob v1.2.0
@iso "Glob1" pkg"add Glob@1.3.0" #Add Glob v1.3.0 as name Glob1
@iso "Glob1" using Glob #Load Glob v1.3.0

@test readdir(IsoPkg.env_path())==["Glob","Glob1","Glob@1.2.0"]

IsoPkg.pin("Glob1")
IsoPkg.free("Glob1")
IsoPkg.pin("Glob") #Pin Glob version (will automatically change its name to match the version)
IsoPkg.free("Glob@1.2.0") #Free Glob v1.2.0 version (will automatically remove the version in its name)

IsoPkg.status() #Show status of all packages
IsoPkg.update() #Update all packages

IsoPkg.rm("Glob1") #Remove Glob v1.3.0
IsoPkg.rm("Glob") #Remove Glob
IsoPkg.rm("Glob") #Because there is only one version left, the version number can be omitted

@test readdir(IsoPkg.env_path())==[]

Base.rm(IsoPkg.env_path())
IsoPkg.switch()


text:
include("../../RawFileofflineProc/dataStruct.jl")
include("WindIDen.jl")
using  Base.Threads
function CalDynmicTh(CumSumDataobj::spaceCumSumData, vPos)
    segID = getAbsPos(CumSumDataobj ,vPos)
    CalTimePeriod,calCFactor  = getSegDynParaValue( glbSegAlarmParasDictObj,segID )
    # curAlarmsinSEG[idx].DynTh    =  typemax(UInt ) # getSegStaticThPar
    vsc=1
    avg = calRectAvg(CumSumDataobj, vPos ,vsc,CalTimePeriod  )
    if calCFactor != typemax( UInt )
        curAlarmsinSEG[idx].DynTh    = avg * calCFactor
    else
        curAlarmsinSEG[idx].DynTh    = typemax( UInt )
    end
end

function postProcFeatureData(dataObject::spaceCumSumData ,data::Array{UInt16,1} )
    cumSumProc( dataObject,  data  )
    # println("          <?>                          ",length(data ))
    vsc = 4 #dataObject.defaultSpaceCnt
    vTcNum = 20 #
    # vTcNum = dataObject.defaultTimeCnt
    ######### calRectSum( CumSumDataobj::spaceCumSumData, vPos ,vsc ,vTcNum   )
    arng =getSpaceLoopRng( dataObject )
    # println( "---------------",arng )

     @threads for vpos = arng
        # println("vpos ,vsc ,vTcNum :$vpos ,$vsc ,$vTcNum   "  )
        calRectSum( dataObject, vpos ,vsc ,vTcNum   )
        # CalDynmicTh(CumSumDataobj, vPos)
    end

    AlarmCal(dataObject)
end

function AlarmCal( CumSumDataobj::spaceCumSumData,threshold = 3000)
    chnID = CumSumDataobj.startID > 1 ? 2 : 1
    # println("-CHnID:-$chnID, $(CumSumDataobj.startID)----")
    return  doesWinding( CumSumDataobj) # in windIDen.jl
    # return  AlarmCalAna( CumSumDataobj,threshold)
end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
λ_params = Dict(
    "A" => (a=-11.4, b=0.10),
    "B" => (a=-26.0, b=0.17),
    "C" => (a=-123.0, b=0.30),
    #"D" => Monin-Obhukov length is infinite
    "E" => (a=123.0, b=0.30),
    "F" => (a=26.0, b=0.17)
)

"""
    monin_obukhov(stability::String, roughness::Number)
returns the Monin-Obukhov length for a given Pasquill-Gifford stability class
and surface roughness (in meters)
Curve fit from
    Pasquill, F., *Atmospheric Diffusion, 2nd Ed.*, Halstead Press, New York, 1974.
"""
function monin_obukhov(stability::String, roughness::Number)
    if stability ∈ Set(["A","B","C","E","F"])
        a, b = λ_params[stability]
        λ = a*roughness^b
    elseif stability == "D"
        λ = Inf
    else
        err = "$stability is not a valid Pasquill-Gifford stability class"
        error(err)
    end

    return λ
end


text:
using SBML
using Catalyst
using EzXML
using MathML
using ModelingToolkit
using OrdinaryDiffEq
using Symbolics
using Test

DATA = joinpath(@__DIR__, "data")
SBML_FILE = joinpath(DATA, "case01.xml")
SBML_DOC = readxml(SBML_FILE)
SYSDICT = sbml_to_sysinfo(SBML_FILE)

t = Variable(:t)
A = Num(Variable(:A))
B = Num(Variable(:B))
k1 = Num(Sym{ModelingToolkit.Parameter{Real}}(:k1))
k2 = Num(Sym{ModelingToolkit.Parameter{Real}}(:k2))
ConvA_k2 = Num(Sym{ModelingToolkit.Parameter{Real}}(:ConvA_k2))
comp1 = Num(Sym{ModelingToolkit.Parameter{Real}}(:comp1))

# test build_par_map
trueparmap = [k1 => 0.8]
parmap = build_par_map(SYSDICT["listOfParameters"])
@test isequal(parmap, trueparmap)

# test build_comp_map
truecompmap = [comp1 => 1.0]
compmap = build_comp_map(SYSDICT["listOfCompartments"])
@test isequal(compmap, truecompmap)

# test build_spec_map
truespecmap = [A => (1.0, comp1),
               B => (0.0, comp1)]
specmap = build_spec_map(SYSDICT["listOfSpecies"])
@test isequal(specmap, truespecmap)

# test promote_local_parameters
doc = promotelocalparameters(readxml(SBML_FILE))
localpar = lastelement(lastelement(lastelement(lastelement(lastelement(firstelement(doc.root))))))
@test isequal(localpar["id"], "ConvA_k2")

# test getkineticlaw
reaction = firstelement(lastelement(firstelement(readxml(SBML_FILE).root)))
kineticlaw = getkineticlaw(reaction)
@test isequal(nodename(kineticlaw),"kineticLaw")

# test getmath
kineticlaw = lastelement(reaction)
math = getmath(kineticlaw)
@test isequal(nodename(math),"math")

# test _getlistofspeciesreference
reactants, r_stoich = SBML._getlistofspeciesreference(reaction,"Reactant")
products, p_stoich = SBML._getlistofspeciesreference(reaction,"Product")
@test isequal(reactants, [A])
@test isequal(r_stoich, [1])
@test isequal(products, [B])
@test isequal(p_stoich, [1])

# test build_reactions()
truesubs = Dict(Num(Variable(:comp1)) => Num(Sym{ModelingToolkit.Parameter{Real}}(:comp1)),
                Num(Variable(:k1)) => Num(Sym{ModelingToolkit.Parameter{Real}}(:k1)),
                Num(Variable(:k2)) => Num(Sym{ModelingToolkit.Parameter{Real}}(:k2)))
reactions = build_reactions([reaction],truesubs)
kineticlaw = comp1*(A*k1 - (B*k2))
truereactions = [(kineticlaw,[A],[B],[1],[1])]
@test isequal(reactions[1][1],truereactions[1][1])

# test _process_doc
model = SBML._process_doc(readxml(SBML_FILE))
@test isequal(model.parameters, [k1 => 0.8, ConvA_k2 => 0.6])
@test isequal(model.compartments, truecompmap)
@test isequal(model.species, truespecmap)
kineticlaw = comp1*(A*k1 - (B*ConvA_k2))
truereactions = [(kineticlaw,[A],[B],[1],[1])]
@test isequal(model.reactions, truereactions)

# test _mathml_substitutions
subs = SBML._mathml_substitutions(SYSDICT)
@test isequal(subs, truesubs)

# test _par_subs
truesubs = [Num(Variable(:k1)) => Num(Sym{ModelingToolkit.Parameter{Real}}(:k1))]
subs = SBML._par_subs(SYSDICT["listOfParameters"])
@test isequal(subs, truesubs)

# test _var_subs

# test ReactionSystem
model = SbmlModel([k1=>1.],[comp1=>1.],[A=>1.],[(comp1*k1*A, [A], nothing)])
rs = ReactionSystem(model)
@test isequal(rs.eqs, Reaction[Reaction(comp1*k1*A, [A], nothing; use_only_rate=true)])
@test isequal(rs.iv, t)
@test isequal(rs.states, [A])
@test isequal(rs.ps, [k1,comp1])

model = SbmlModel(SBML_DOC)
rs = ReactionSystem(model)
@test isequal(rs.eqs, Reaction[Reaction(comp1*(k1*A - (ConvA_k2*B)), [A], [B]; use_only_rate=true)])
@test isequal(rs.iv, t)
@test isequal(rs.states, [A, B])
@test isequal(rs.ps, [k1,ConvA_k2,comp1])

model = SbmlModel(SBML_FILE)
rs = ReactionSystem(model)
@test isequal(rs.eqs, Reaction[Reaction(comp1*(k1*A - (ConvA_k2*B)), [A], [B]; use_only_rate=true)])
@test isequal(rs.iv, t)
@test isequal(rs.states, [A, B])
@test isequal(rs.ps, [k1,ConvA_k2,comp1])

# test ODESystem
odesys = ODESystem(model)  # Todo: Add test cases
trueeqs = Equation[Differential(t)(A) ~ -comp1*(A*k1 - (B*ConvA_k2)), Differential(t)(B) ~ comp1*(A*k1 - (B*ConvA_k2))]
@test isequal(odesys.eqs, trueeqs)
@test isequal(odesys.iv, t)
@test isequal(odesys.states, [A, B])
# isequal(Sym{ModelingToolkit.Parameter{Real},Nothing}[k1, ConvA_k2, comp1, comp1], [k1, ConvA_k2, comp1])
@test isequal(odesys.ps, [k1, ConvA_k2, comp1])
u0 = [A => 1., B => 0.]
par = [k1 => 0.8, ConvA_k2 => 0.6, comp1 => 1.]
@test isequal(odesys.defaults, Dict(append!(u0, par)))

odesys = ODESystem(SBML_DOC)
e1 = odesys.eqs[1]
println(e1)
println(typeof(e1))
println(e1.lhs)
println(e1.rhs)
@test isequal(odesys.eqs, trueeqs)
@test isequal(odesys.iv, t)
@test isequal(odesys.states, [A, B])
# isequal(Sym{ModelingToolkit.Parameter{Real},Nothing}[k1, ConvA_k2, comp1, comp1], [k1, ConvA_k2, comp1])
@test isequal(odesys.ps, [k1, ConvA_k2, comp1])
u0 = [A => 1., B => 0.]
par = [k1 => 0.8, ConvA_k2 => 0.6, comp1 => 1.]
@test isequal(odesys.defaults, Dict(append!(u0, par)))

odesys = ODESystem(SBML_FILE)
@test isequal(odesys.eqs, trueeqs)
@test isequal(odesys.iv, t)
@test isequal(odesys.states, [A, B])
# isequal(Sym{ModelingToolkit.Parameter{Real},Nothing}[k1, ConvA_k2, comp1, comp1], [k1, ConvA_k2, comp1])
@test isequal(odesys.ps, [k1, ConvA_k2, comp1])
u0 = [A => 1., B => 0.]
par = [k1 => 0.8, ConvA_k2 => 0.6, comp1 => 1.]
@test isequal(odesys.defaults, Dict(append!(u0, par)))

# test ODEProblem
oprob = ODEProblem(model, [0., 10.])
sol = solve(oprob, Tsit5())
println("soln")
println(sol.u)
println(typeof(sol.u))
@test isapprox(sol.u[end], [0.42857284277223784, 0.5714271572277619])

@test_nowarn ODEProblem(SBML_DOC, [0., 10.])

@test_nowarn ODEProblem(SBML_FILE, [0., 10.])

# test _substitute_par!
doc = deepcopy(SBML_DOC)
kineticlaw = findfirst("//x:kineticLaw", doc.root, ["x"=>namespace(doc.root)])
SBML._substitute_par!(kineticlaw, "k2", "ConvA")
k1 = Num(Variable(:k1))
ConvA_k2 = Num(Variable(:ConvA_k2))
comp1 = Num(Variable(:comp1))
@test isequal(parse_node(firstelement(kineticlaw))[1], comp1*(k1*A - (ConvA_k2*B)))

#=DATA = "./test/data"
SBML_FILE = joinpath(DATA, "case01.xml")
SBML_DOC = readxml(SBML_FILE)
node = SBML_DOC.root
ns = namespace(node)
list_name = "listOfParameters"
list_name = "listOfReactions/x:reaction/x:kineticLaw/x:listOfParameters"
list_name = "listOfReactions/x:reaction/x:kineticLaw/x:listOfLocalParameters"
str = "/x:sbml/x:model/x:$(list_name)"
findall(str, node, ["x" => ns])

using ModelingToolkit
@parameters β γ t
@variables S(t) I(t) R(t)

rxs = [Reaction(β, [S,I], [I], [1,1], [2])
       Reaction(γ, [I], [R])]
rs  = ReactionSystem(rxs, t, [S,I,R], [β,γ])

u₀map    = [S => 999.0, I => 1.0, R => 0.0]
parammap = [β => 1/10000, γ => 0.01]
tspan    = (0.0, 250.0)

# solve as ODEs
odesys = ODESystem(rs; defaults=merge(u₀map,parammap))=#