text:
function _recon2d_slices_tvrdart!(u, A, b0, niter, w_tv, sigmas, tau)
    At = sparse(A')
    H, W = size(u)

    b = vec(b0)

    return u
end

function obj_param(A, p, u, w_tv, epsilon)

end

function der_param()

end

function _estimate_param(A, p, u, param0, w_tv, epsilon)
    obj = obj_param(A, p, u, w_tv, epsilon)
    objder = der_param()
    # residual = minimize()..

    param = res.x
    return param
end

"""
    recon2d_slices_tvrdart!(u::Array{T, 2}, A, b, niter::Int, w_tv::T, c=10.0)

Reconstruct 3D slice by slice by TVR-DART

# Args
u : Initial guess of image
A : Forward opeartor
b : Projection data (2 dimension or 1 dimension)
niter: number of iterations (<50)
w_tv: weight for TV term
"""
function recon2d_slices_tvrdart!(u::Array{T, 3}, A, b, niter::Int, w_tv::T, nmaterials, epsilon=1e-4, K=4) where {T <: AbstractFloat}

    tau = c / sum(ops_norm)
    println("@ step sizes sigmas: ", sigmas, ", tau: $tau")

    recon2d_slices_sirt!(u, A, b, 100)
    u_max = maximum(u)
    b .= b ./ u_max
    u .= u ./ u_max

    # initial params
    K = 4*ones(nmaterials-1)
    gv = LinRange(0, 1, nmaterials)

    params = [gv[2:], K]

    # estimate optimal params by small subset
    sum(u, dims=[1,3])

    # solve gray values and thresholds
    yhalf = Int(floor(size(u,2)/2))
    yrange = [ yhalf, yhalf+2 ]
    u_slice = u[:, yhalf:yhalf+2, :]
    data_slice = b[:, :, yhalf:yhalf+2]
    p_slice = reshape(data_slice)
    grays = range(0, stop=1, length=nmaterials)[2:end-1]
    param0 = [grays[2:end]..., w_tv]
    param = _estimate_param(A, p_slice, u_slice, param0, w_tv, epsilon)

    return _recon2d_slices_tvrdart!(u, A, b, niter, w_tv, param)
end

def gv2param(gv,K):
    param = np.hstack((gv[1:],K))
    return param;

max_stars_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_issues_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_forks_repo_name: JuliaTomo/DynamicTomographyWithFlow


text:
using Reexport

@reexport module TriangulationPartitions
    # Necessary modules
    include("Simplices/Simplices.jl")
    include("DelaunayTriangulations/DelaunayTriangulations.jl")

	dimension(DT::DelaunayTriangulation) = DelaunayTriangulations.dimension(DT)

    # Abstract types and functions that work on all subtypes
    include("AbstractTriangulationPartition.jl")

    # Composite types for different types of data
    include("composite_types/DatasetTriangulationPartition.jl")
    include("composite_types/EmbeddingTriangulationPartition.jl")
    include("composite_types/PointTriangulationPartition.jl")
	include("composite_types/MatrixTriangulationPartition.jl")

    # Functions to create triangulations from data
    include("triangulate.jl")

    # Plot triangulations
    include("plot_recipes/plot_recipes.jl")

end

"""
    TriangulationPartitions

A module handling partitioned formed by triangulations.
"""
TriangulationPartitions

ext: jl


text:
if VERSION >= v"1.5.0"
    using Pkg

    Pkg.develop(PackageSpec(url="https://github.com/carlobaldassi/GaussDCA.jl"))

    msa = map(Residue, rand(1:21,100,20))
    dca = gaussdca(msa, min_separation=2)

    @test  isnan(dca[1,1])
    @test  isnan(dca[1,2])
    @test !isnan(dca[1,3])
end


text:
################################################################################
#
# More interesting search: Enumeration and other things on the Minkowski side
#
################################################################################

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

function enum_ctx_from_ideal(A::NfOrdIdl,
                v::fmpz_mat; prec::Int = 100, limit::Int = 0, Tx::DataType = Int, TU::DataType = Float64, TC::DataType = Float64)

  l, t = lll(A, v, prec = prec)
  OK = order(A)
  K = nf(OK)
  temp = FakeFmpqMat(basis_matrix(A, copy = false))*basis_matrix(OK, copy = false)
  b = temp.num
  b_den = temp.den

  if limit == 0
    limit = nrows(l)
  end

  E = enum_ctx_from_gram(l, Tx = Tx, TC = TC, TU = TU,
                  limit = limit)::enum_ctx{Tx, TC, TU}
  E.t = t*b
  E.t_den = b_den
  ## we want to find x sth. norm(x) <= sqrt(|disc|)*norm(A)
  ## |N(x)^2|^(1/n) <= T_2(x)/n
  ## so if T_2(x) <= n * D^(1/n)
  ## then |N(x)| <= D^(1/2)
  #d = abs(discriminant(order(A))) * norm(A)^2
  #due to limit, the lattice is smaller and the disc as well.
  d = fmpz(ceil(abs(prod(TC[TC(E.C[i,i]) for i=1:E.limit]))))
  ## but we don't want to overshoot too much the length of the last
  ## basis element.
  den = basis_matrix(OK, copy = false).den ## we ignore the den above, but this
                                ## changes the discriminant!!!
  b = min(den^2 * (root(d, E.limit)+1)*E.limit * E.d, E.G[E.limit, E.limit]*E.limit)
  @v_do :ClassGroup 3 println("T_2 from disc ", (root(d, E.limit)+1)*E.limit * E.d)
  @v_do :ClassGroup 3 println("    from Gram ", E.G[E.limit, E.limit]*E.limit)
  @v_do :ClassGroup 3 println(" using ", b)
  enum_ctx_start(E, b)
  return E
end

_start = 0.0
function class_group_small_real_elements_relation_start(clg::ClassGrpCtx,
                A::NfOrdIdl; prec::Int = 200, val::Int = 0,
                limit::Int = 0)
  global _start
  @v_do :ClassGroup_time 2 rt = time_ns()
  while true
    try
      I = IdealRelationsCtx{Int, Float64, Float64}(clg, A, prec = prec, val = val, limit = limit)
      @v_do :ClassGroup_time 2 _start += time_ns()-rt
      return I
    catch e
      if isa(e, LowPrecisionCholesky)
        print_with_color(:red, "prec too low in cholesky,")
        prec = Int(ceil(1.2*prec))
        println(" increasing to ", prec)
        if prec > 1000
          _start = A
          error("1:too much prec")
        end
      elseif isa(e, LowPrecisionLLL)
        print_with_color(:red, "prec too low in LLL,")
        prec = Int(ceil(1.2*prec))
        println(" increasing to ", prec)
        if prec > 1000
          error("2:too much prec")
        end
      else
        rethrow(e)
      end
    end
  end
end

global _elt = UInt(0)

function class_group_small_real_elements_relation_next(I::IdealRelationsCtx)
  global _elt, _next
  K = nf(order(I.A))
  while true
    if enum_ctx_next(I.E)
#       println("elt is", I.E.x)
      # should we (again) add content_is_one()?
      if !isone(content(I.E.x))
        continue
      end
      @v_do :ClassGroup_time 2 rt = time_ns()
#      I.M = I.E.x * I.E.t
      ccall((:fmpz_mat_mul, libflint), Nothing, (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{fmpz_mat}), I.M, I.E.x, I.E.t)
      q = elem_from_mat_row(K, I.M, 1, I.E.t_den)
#      println("found ", q, " norm ", norm(q)//norm(I.A))
      @v_do :ClassGroup_time 2 _elt += time_ns()- rt
      return q
    end
    @v_do :ClassGroup 2 print_with_color(:red, "restart after ")
    @v_do :ClassGroup 2 print(I.E.cnt)
    @v_do :ClassGroup 3 println(" for ", I.A, I.E.c)
    @v_do :ClassGroup 2 println(" length now ", I.E.c*2)
#    throw(NoElements())
    I.restart += 1
    if I.restart > 100
      _elt = I
      error("too much restarting");
    end
    enum_ctx_start(I.E, I.E.c*2)
  end
end

function show(io::IO, I::IdealRelationsCtx)
  println(io, "IdealRelationCtx for ", I.A)
  println(io, "  current length bound ", I.c, " stats: ", I.cnt, " and ", I.bad)
end


text:
"""
based on  https://github.com/Visceral-Project/EvaluateSegmentation/blob/master/source/VolumeSimilarityCoefficient.h

"""
module VolumeMetric
using ..BasicStructs, Parameters

"""
calculate volume metric  based on precalulated constants
TnTpFpFns - list of basic metrics - in the order [tn,tp,fp,fn ]
res - ResultMetrics struct holding result of all metrics calculated for given run
return modified ResultMetrics with added volume metric
"""
function calculateVolumeMetric(TnTpFpFns::Vector{Float64} , res::ResultMetrics) ::ResultMetrics
    tn = TnTpFpFns[1]
    tp = TnTpFpFns[2]
    fp = TnTpFpFns[3]
    fn = TnTpFpFns[4]
return setproperties(res, (vol= getVolumMetric(tp,fp, fn)  ))

end #calculateVolumeMetric

"""
calculates volume metric
"""
function getVolumMetric(tp,fp, fn )
    return 1- abs(fn-fp)/(2*tp + fn + fp)
end#getVolumMetric

end#VolumeMetric


text:
using sbl

using HTTP
using JSON3
using DataFrames

function download_json(cid, record_type="3d")
    pc_uri = "https://pubchem.ncbi.nlm.nih.gov/rest/pug/compound/cid/$cid/JSON?record_type=$record_type"

    HTTP.request("GET", pc_uri).body
end

function read_json(filename)
    pb_json_txt = ""

    open(filename, "r") do pb_json_f
        pb_json_txt = read(pb_json_f, String)
    end

    pb_json_txt
end

pb_json_txt = read_json("sbl/test/data/aspirin_pug.json");

r = JSON3.read(pb_json_txt, PCResult)

# pb_json = JSON.parse(String(pb_json_response.body))

# for (i, compound) in enumerate(pb_json["PC_Compounds"])

# end

# pb_json = jsontable(String(pb_json_response.body))

# pb_df = DataFrame(pb_json)

# print(pb_df)


text:
global const marmousi_folder=joinpath(dirname(pathof(GeoPhyInv)), "media", "marmousi2")
global const overthrust_folder=joinpath(dirname(pathof(GeoPhyInv)), "media", "overthrust")

"""
```julia
mod=Medium(attrib)
```
As of now, only seismic models are predefined in this package. Choose `attrib::Symbol`

* `=:acou_homo2D` : a test homogeneous acoustic model
* `=:acou_homo2` : a test homogeneous acoustic model, but with coarser spatial sampling (faster testing)
* `=:marmousi2` : marmousi2 model with lower resolution; useful for surface seismic experiments
* `=:marmousi2_small` : a smaller section of marmousi2
* `=:pizza` : `:acou_homo2` with some perturbations
"""
function Medium(attrib::Symbol, δ::Real = 0.0; verbose = false)
    bfrac = 0.1
    δ = Float64(δ)
    vp0 = [1500.0, 3500.0] # some bounds for vp
    rho0 = [1500.0, 3500.0] # some density bounds
    vs0 = [1000, 2000] # bounds for vs
    if ((attrib == :acou_homo2D))
        mgrid = repeat([range(-1000.0, stop = 1000.0, length = 201)], 2)
        model = Medium(mgrid, [:vp, :rho])
        update!(model, [:vp, :rho], [vp0, rho0])
        fill!(model)
    elseif ((attrib == :acou_homo3D))
        mgrid = repeat([range(-1000.0, stop = 1000.0, length = 201)], 3)
        model = Medium(mgrid, [:vp, :rho])
        update!(model, [:vp, :rho], [vp0, rho0])
        fill!(model)
    elseif ((attrib == :elastic_homo2D))
        mgrid = fill(range(-1000, stop = 1000, length = 255), 2)
        model = Medium(mgrid, [:vp, :vs, :rho])
        update!(model, [:vp, :vs, :rho], [vp0, vs0, rho0])
        fill!(model)
    elseif ((attrib == :elastic_homo3D))
        mgrid = fill(range(-1000, stop = 1000, length = 255), 3)
        model = Medium(mgrid, [:vp, :vs, :rho])
        update!(model, [:vp, :vs, :rho], [vp0, vs0, rho0])
        fill!(model)
    elseif ((attrib == :pizza))
        vp0 = [1700.0, 2300.0] # bounds for vp
        rho0 = [1700.0, 2300.0] # density bounds
        mgrid = repeat([range(-1000.0, stop = 1000.0, length = 51)], 2)
        nz, nx = length.(mgrid)
        model = Medium(mgrid, [:vp, :rho])
        update!(model, [:vp, :rho], [vp0, rho0])
        fill!(model)

        # add some noise to starting model
        update!(model, [:vp, :rho], randn_perc = 0.5)

        # add perturbations
        for ellip_loc in [[500.0, 0.0], [0.0, 500.0], [-500.0, 0.0], [0.0, -500.0]]
            update!(
                model,
                [:vp, :rho],
                ellip_rad = 50.0,
                ellip_loc = ellip_loc,
                ellip_pert = 100.0,
            )
        end

    elseif ((attrib == :acou_homo2))
        vp0 = [1700.0, 2300.0] # bounds for vp
        rho0 = [1700.0, 2300.0] # density bounds
        mgrid = repeat([range(-1000.0, stop = 1000.0, length = 51)], 2)
        nz, nx = length.(mgrid)
        model = Medium(mgrid, [:vp, :rho])
        update!(model, [:vp, :rho], [vp0, rho0])
        fill!(model)

    elseif (attrib == :marmousi2)
        c = h5open(joinpath(marmousi_folder, "marmousi2.h5"), "r") do file
            vp = read(file, "vp")
            vs = read(file, "vs")
            rho = read(file, "rho")
            xgrid = read(file, "xgrid")
            zgrid = read(file, "zgrid")
            mgrid = [
                range(zgrid[1], stop = zgrid[end], length = size(vp, 1)),
                range(xgrid[1], stop = xgrid[end], length = size(vp, 2)),
            ]
            model = Medium(mgrid, [:vp, :rho, :vs])
            copyto!(model[:vp], vp)
            copyto!(model[:rho], rho)
            copyto!(model[:vs], vs)
            update!(model, bfrac)
        end
    elseif (attrib == :marmousi2_small)
        fmodel = Medium(:marmousi2)
        fmgrid = fmodel.mgrid
        mgrid = [
            range(500, stop = 3500, step = step(fmgrid[1])),
            range(4000, step = step(fmgrid[2]), stop = 13000),
        ]
        model = update(fmodel, mgrid)

    elseif (attrib == :overthrust)
        vp = []
        for i in range(1, 8, step = 1)
            file_ =
                h5open(joinpath(overthrust_folder, "overthrust_" * string(i) * ".h5"), "r")
            vp_sect = read(file_, "vp")
            if i == 1
                vp = deepcopy(vp_sect)
            else
                vp = cat(vp, vp_sect, dims = 3)
            end
            close(file_)
        end
        file_ = h5open(joinpath(overthrust_folder, "overthrust_mgrid.h5"), "r")
        xgrid = read(file_, "xgrid")
        ygrid = read(file_, "ygrid")
        zgrid = read(file_, "zgrid")
        close(file_)
        mgrid = [
            range(zgrid[1], stop = zgrid[end], length = size(vp, 1)),
            range(ygrid[1], stop = ygrid[end], length = size(vp, 2)),
            range(xgrid[1], stop = xgrid[end], length = size(vp, 3)),
        ]
        model = Medium(mgrid, [:vp, :rho])
        rho0 = [1500.0, 2500.0]
        update!(model, [:rho], [rho0])
        fill!(model)
        copyto!(model[:vp], vp)
        update!(model, bfrac)

        #=
        	elseif(attrib == :seismic_marmousi2_surf)
        		model=Medium_trun(Seismic(:seismic_marmousi2_high_res),
        				     xmin=6000., xmax=12000.,)
        		update!(model, bfrac) # adjust bounds just inside the bounds
        	elseif(attrib == :seismic_marmousi2_downhole)
        		model=Medium_trun(Seismic(:seismic_marmousi2_high_res),
        				     xmin=9025., xmax=9125., zmin=1400., zmax=1600.,)
        		update!(model, bfrac) # adjust bounds just inside the bounds
        	elseif(attrib == :seismic_marmousi2_rvsp)
        		model=Medium_trun(Seismic(:seismic_marmousi2_high_res),
        				     xmin=8000., xmax=10000., zmax=1700.,zmin=500.)
        		update!(model, bfrac) # adjust bounds just inside the bounds
        		=#

    else
        error("invalid attrib")
    end
    if (δ == 0.0)
        verbose && Models.print(model, string(attrib))
        return model
    elseif (δ > 0.0)
        mgrid_out = broadcast(x -> range(x[1], stop = x[end], step = δ), model.mgrid)
        model_out = update(model, mgrid_out)
        verbose && Models.print(model_out, string(attrib))
        return model_out
    else
        error("invalid δ")
    end

end


text:
@testset "solve_adbc.jl" begin
    a = 0.0
    b = 1.0
    domain = ClosedIntervalDomain(a, b)
    g_a(t) = 0.0
    g_b(t) = 0.0
    f(ξ, t) = (1.0 - ξ)ξ
    u₀(ξ) = 0.0
    T = 1.0
    N = 100
    Δx = (b - a) / N
    h(x) = 3.1Δx
    σ = Δx / 2
    filter = GaussianFilter(h, σ)
    equation = DiffusionEquation(domain, filter, f, g_a, g_b)
    ū_adbc = solve_adbc(equation, u₀, (0.0, T), N, T / 1000)
end


text:
using TimerOutputs, Printf
using Oceananigans

include("benchmark_utils.jl")

#####
##### Benchmark setup and parameters
#####

const timer = TimerOutput()

Ni = 2   # Number of iterations before benchmarking starts.
Nt = 10  # Number of iterations to use for benchmarking time stepping.

# Run benchmark across these parameters.
            Ns = [(128, 128, 128)]
   float_types = [Float64]  # Float types to benchmark.
         archs = [CPU()]    # Architectures to benchmark on.
@hascuda archs = [GPU()]    # Benchmark GPU on systems with CUDA-enabled GPUs.

#####
##### Forcing function definitions
#####

@inline function Fu_params(i, j, k, grid, time, U, C, params)
    if k == 1
        return @inbounds -2*params.K/grid.Δz^2 * (U.u[i, j, 1] - 0)
    elseif k == grid.Nz
        return @inbounds -2*params.K/grid.Δz^2 * (U.u[i, j, grid.Nz] - 0)
    else
        return 0
    end
end

const K = 0.1
@inline function Fu_consts(i, j, k, grid, time, U, C, params)
    if k == 1
        return @inbounds -2*K/grid.Δz^2 * (U.u[i, j, 1] - 0)
    elseif k == grid.Nz
        return @inbounds -2*K/grid.Δz^2 * (U.u[i, j, grid.Nz] - 0)
    else
        return 0
    end
end

@inline FT_params(i, j, k, grid, time, U, C, params) = @inbounds ifelse(k == 1, -params.λ * (C.T[i, j, 1] - 0), 0)

const λ = 1e-4
@inline FT_consts(i, j, k, grid, time, U, C, params) = @inbounds ifelse(k == 1, -λ * (C.T[i, j, 1] - 0), 0)

#####
##### Run benchmarks
#####

for arch in archs, FT in float_types, N in Ns
    Nx, Ny, Nz = N
    Lx, Ly, Lz = 1, 1, 1

    forced_model_params = Model(architecture = arch, float_type = FT,
		                grid = RegularCartesianGrid(size=(Nx, Ny, Nz), length=(Lx, Ly, Lz)),
                                forcing=ModelForcing(Fu=Fu_params, FT=FT_params), parameters=(K=0.1, λ=1e-4))

    time_step!(forced_model_params, Ni, 1)  # First 1-2 iterations usually slower.

    bn =  benchmark_name(N, "with forcing (params)", arch, FT)
    @printf("Running benchmark: %s...\n", bn)
    for i in 1:Nt
        @timeit timer bn time_step!(forced_model_params, 1, 1)
    end

    forced_model_consts = Model(architecture = arch, float_type = FT,
		                grid = RegularCartesianGrid(size=(Nx, Ny, Nz), length=(Lx, Ly, Lz)),
                                forcing=ModelForcing(Fu=Fu_consts, FT=FT_consts))

    time_step!(forced_model_consts, Ni, 1)  # First 1-2 iterations usually slower.

    bn =  benchmark_name(N, "with forcing (consts)", arch, FT)
    @printf("Running benchmark: %s...\n", bn)
    for i in 1:Nt
        @timeit timer bn time_step!(forced_model_consts, 1, 1)
    end

    unforced_model = Model(architecture = arch, float_type = FT,
			   grid = RegularCartesianGrid(size=(Nx, Ny, Nz), length=(Lx, Ly, Lz)))

    time_step!(unforced_model, Ni, 1)  # First 1-2 iterations usually slower.

    bn =  benchmark_name(N, "  no forcing         ", arch, FT)
    @printf("Running benchmark: %s...\n", bn)
    for i in 1:Nt
        @timeit timer bn time_step!(unforced_model, 1, 1)
    end
end

#####
##### Print benchmark results
#####

println()
print_benchmark_info()
print_timer(timer, title="Forcing function benchmarks")
println()


text:
mutable struct MultiStepDynamicBuffer{O,A} <: AbstractBuffer
    transitions::Vector{Transition{O,A}}

    MultiStepDynamicBuffer{O,A}() where {O,A} = new(Transition{O,A}[])
end

function add!(buffer::MultiStepDynamicBuffer{O,A}, obs::O, action::A, reward::Float64) where {O,A}
    push!(buffer.transitions, Transition{O,A}(obs, action, reward))
end

reset!(buffer::MultiStepDynamicBuffer) = empty!(buffer.transitions)

function Base.iterate(buffer::MultiStepDynamicBuffer, state=1)
    state > length(buffer.transitions) && return nothing

    return buffer.transitions[state], state + 1
end

function Base.iterate(rbuffer::Base.Iterators.Reverse{MultiStepDynamicBuffer{O,A}}, state=length(rbuffer.itr.transitions)) where {O,A}
    state < 1 && return nothing

    return rbuffer.itr.transitions[state], state - 1
end

Base.length(buffer::MultiStepDynamicBuffer) = length(buffer.transitions)

Base.eltype(::MultiStepDynamicBuffer{O,A}) where {O,A} = Transition{O,A}

Base.popfirst!(buffer::MultiStepDynamicBuffer) = popfirst!(buffer.transitions)


text:
@testset "exports" begin
    # Ensure that every exported name is actually defined
    for name in names(RadauBig)
        @test isdefined(RadauBig, name)
        !isdefined(RadauBig, name) && println("MISSING: ", name)
    end
end


text:
include("entity.jl")

function generatePopulation(n, genesLength)
    data = []
    for i in 1:n
        entity = generateEntity(genesLength)
        push!(data,entity)
    end
    return data
end

function calculatePopulationFitness!(data,values, weights, maxCapacity)
    for i in 1:length(data)
        calculateFitness!(data[i],values,weights,maxCapacity)
    end
    sort!(data, by = d -> d.fitness, rev = false)
end

function printPopulation(data)
    for i in 1:length(data)
        printEntity(data[i])
    end
end

function selectPopulation(data,n)
    return copy(data[1:n])
end

function crossoverPopulation(data,crossoverPoint1,crossoverPoint2)
    newData = []
    for i in 1:length(data)
        for j in 1:length(data)
            entity1 = deepcopy(data[i])
            entity2 = deepcopy(data[j])
            crossover!(entity1,entity2,crossoverPoint1,crossoverPoint2)
            push!(newData,entity1)
            push!(newData,entity2)
        end
    end
    return newData
end

function mutatePopulation!(data, mutatePercentage)
    for i in 1:length(data)
        mutate!(data[i],mutatePercentage)
    end
end

max_stars_repo_name: pikacxe/Primenjeni-Algoritmi

max_issues_repo_name: pikacxe/Primenjeni-Algoritmi

max_forks_repo_name: pikacxe/Primenjeni-Algoritmi


text:
function createButton(bar_name::String, on_change, userdata, type::Int32 = 0, initial_button_state::Bool = false)
    func =  (x)->on_change(x, userdata)
    CxxWrap.gcprotect(userdata)
    CxxWrap.gcprotect(func)
    CxxWrap.gcprotect(on_change)
    return jl_cpp_cv2.createButton(bar_name,func, type, initial_button_state)
end

function setMouseCallback(winname::String, onMouse, userdata)
    func =  (event, x, y, flags)->onMouse(event, x, y, flags, userdata)
    CxxWrap.gcprotect(userdata)
    CxxWrap.gcprotect(func)
    CxxWrap.gcprotect(onMouse)
    return jl_cpp_cv2.setMouseCallback(winname,func)
end

function createTrackbar(trackbarname::String, winname::String, value::Ref{Int32}, count::Int32, onChange, userdata)
    func =  (x)->onChange(x, userdata)
    CxxWrap.gcprotect(userdata)
    CxxWrap.gcprotect(func)
    CxxWrap.gcprotect(onChange)
    return jl_cpp_cv2.createTrackbar(trackbarname, winname, value, count, func)
end

function CascadeClassifier(filename::String)
	return cpp_to_julia(jlopencv_cv_cv_CascadeClassifier_cv_CascadeClassifier_CascadeClassifier(julia_to_cpp(filename)))
end

function detect(cobj::cv_Ptr{T}, image::InputArray, mask::InputArray) where {T <: Feature2D}
	return cpp_to_julia(jlopencv_cv_cv_Feature2D_cv_Feature2D_detect(julia_to_cpp(cobj),julia_to_cpp(image),julia_to_cpp(mask)))
end
detect(cobj::cv_Ptr{T}, image::InputArray; mask::InputArray = (CxxMat())) where {T <: Feature2D} = detect(cobj, image, mask)

function detectMultiScale(cobj::CascadeClassifier, image::InputArray, scaleFactor::Float64, minNeighbors::Int32, flags::Int32, minSize::Size{Int32}, maxSize::Size{Int32})
	return cpp_to_julia(jlopencv_cv_cv_CascadeClassifier_cv_CascadeClassifier_detectMultiScale(julia_to_cpp(cobj),julia_to_cpp(image),julia_to_cpp(scaleFactor),julia_to_cpp(minNeighbors),julia_to_cpp(flags),julia_to_cpp(minSize),julia_to_cpp(maxSize)))
end
detectMultiScale(cobj::CascadeClassifier, image::InputArray; scaleFactor::Float64 = Float64(1.1), minNeighbors::Int32 = Int32(3), flags::Int32 = Int32(0), minSize::Size{Int32} = (Size{Int32}(0,0)), maxSize::Size{Int32} = (Size{Int32}(0,0))) = detectMultiScale(cobj, image, scaleFactor, minNeighbors, flags, minSize, maxSize)

function empty(cobj::CascadeClassifier)
	return cpp_to_julia(jlopencv_cv_cv_CascadeClassifier_cv_CascadeClassifier_empty(julia_to_cpp(cobj)))
end

function SimpleBlobDetector_create(parameters::SimpleBlobDetector_Params)
	return cpp_to_julia(jlopencv_cv_cv_SimpleBlobDetector_create(julia_to_cpp(parameters)))
end
SimpleBlobDetector_create(; parameters::SimpleBlobDetector_Params = (SimpleBlobDetector_Params())) = SimpleBlobDetector_create(parameters)


text:
np = 4
nt = 500

if nworkers() < np
    addprocs(np - nworkers())
    @everywhere using MessageUtils
end
pids = workers()

ts = tspace()

#Launch processes on each worker that simulates CPU intensive / blocking stuff for a random amount of time
@everywhere begin
    function simworker(ts)
        if iseven(myid())
            key = string(myid())
        else
            key = r".*"
        end

        while true
            r = take!(ts, key)
#            lsecs = rand()
            lsecs = 0.1
            t1 = time()
            while (time() - t1) < 0.1
                # empty loop simulating CPU intensive code
            end

            put!(r[2], "Finished $(r[1]) in $lsecs seconds @ worker $(myid())")
        end
    end
end

for p in pids
    @spawnat p simworker(ts)
end

# Start n tasks concurrently - each adding its request to ts and waiting for a result
tic()
@sync begin
    for i in 1:nt
        @async begin
            rr = RemoteRef()
            put!(ts, ("$((i%(np))+2)", rr))
            result = take!(rr)
#            println(result)
#             println("jobs left to be taken : $(length(ts))")
        end
    end
end
toc()


text:
""" Computing in the category of finite sets and Petri cospans
"""
module AlgebraicPetri

export TheoryPetriNet, PetriNet, OpenPetriNetOb, AbstractPetriNet, ns, nt, ni, no,
  os, ot, is, it,
  add_species!, add_transition!, add_transitions!,
  add_input!, add_inputs!, add_output!, add_outputs!, inputs, outputs,
  TransitionMatrices, vectorfield,
  TheoryLabelledPetriNet, LabelledPetriNet, AbstractLabelledPetriNet, sname, tname, snames, tnames,
  TheoryReactionNet, ReactionNet, AbstractReactionNet, concentration, concentrations, rate, rates,
  TheoryLabelledReactionNet, LabelledReactionNet, AbstractLabelledReactionNet,
  Open, OpenPetriNet, OpenLabelledPetriNet, OpenReactionNet, OpenLabelledReactionNet,
  OpenPetriNetOb, OpenLabelledPetriNetOb, OpenReactionNetOb, OpenLabelledReactionNetOb

using Catlab
using Catlab.CategoricalAlgebra
using Catlab.CategoricalAlgebra.FinSets
using Catlab.Present
using Catlab.Theories
using LabelledArrays
using LinearAlgebra: mul!

vectorify(n::Vector) = n
vectorify(n::Tuple) = length(n) == 1 ? [n] : n
vectorify(n) = [n]

state_dict(n) = Dict(s=>i for (i, s) in enumerate(n))

# Petri Nets
############

""" ACSet definition for a Petri net.

See Catlab.jl documentation for description of the @present syntax.
"""
@present TheoryPetriNet(FreeSchema) begin
  T::Ob
  S::Ob
  I::Ob
  O::Ob

  it::Hom(I,T)
  is::Hom(I,S)
  ot::Hom(O,T)
  os::Hom(O,S)
end

@abstract_acset_type AbstractPetriNet
@acset_type PetriNet(TheoryPetriNet,index=[:it,:is,:ot,:os]) <: AbstractPetriNet
const OpenPetriNetOb, OpenPetriNet = OpenCSetTypes(PetriNet,:S)

""" Open(p::AbstractPetriNet)

Converts a PetriNet to an OpenPetriNet where each state is exposed as a leg of
the cospan. The OpenPetriNet can be composed over an undirected wiring diagram
(see this
[blog post](https://www.algebraicjulia.org/blog/post/2020/11/structured-cospans-2/)
for a description of this compositional tooling)
"""
Open(p::AbstractPetriNet) = OpenPetriNet(p, map(x->FinFunction([x], ns(p)), 1:ns(p))...)

""" Open(p::AbstractPetriNet, legs...)

Generates on OpenPetriNet with legs bundled as described by `legs`
"""
Open(p::AbstractPetriNet, legs...) = OpenPetriNet(p, map(l->FinFunction(l, ns(p)), legs)...)

""" Open(n, p::AbstractPetriNet, m)

Generates on OpenPetriNet with two legs, `n` and `m`
"""
Open(n, p::AbstractPetriNet, m) = Open(p, n, m)

""" PetriNet(n::Int, ts::Vararg{Union{Pair,Tuple}})

Constructs a PetriNet object with `n` states and transitions described by `ts`.
Transitions are given as `(input_states)=>(output_states)`.

A PetriNet modelling the SIR model with 3 states and 2 transitions can be
constructed as follows:
```@example
PetriNet(3, (1,2)=>(2,2), 2=>3)
```
"""
PetriNet(n::Int, ts::Vararg{Union{Pair,Tuple}}) = begin
  p = PetriNet()
  add_species!(p, n)
  add_transitions!(p, length(ts))
  for (i,(ins,outs)) in enumerate(ts)
    ins = vectorify(ins)
    outs = vectorify(outs)
    add_inputs!(p, length(ins), repeat([i], length(ins)), collect(ins))
    add_outputs!(p, length(outs), repeat([i], length(outs)), collect(outs))
  end
  p
end

""" Number of states in a Petri net
"""
ns(p::AbstractPetriNet) = nparts(p,:S)

""" Number of transitions in a Petri net
"""
nt(p::AbstractPetriNet) = nparts(p,:T)

""" Number of input relationships in a Petri net
"""
ni(p::AbstractPetriNet) = nparts(p,:I)

""" Number of output relationships in a Petri net
"""
no(p::AbstractPetriNet) = nparts(p,:O)

is(p::AbstractPetriNet, args...) = subpart(p, args..., :is)
os(p::AbstractPetriNet, args...) = subpart(p, args..., :os)
it(p::AbstractPetriNet, args...) = subpart(p, args..., :it)
ot(p::AbstractPetriNet, args...) = subpart(p, args..., :ot)

""" Add a species to the Petri net. Label and concentration can be provided
depending on the kind of Petri net.

Returns the ID of the species
"""
add_species!(p::AbstractPetriNet;kw...) = add_part!(p,:S;kw...)

""" Add `n` species to the Petri net. Label and concentration can be provided
depending on the kind of Petri net.

Returns the ID of the species
"""
add_species!(p::AbstractPetriNet,n;kw...) = add_parts!(p,:S,n;kw...)

""" Add a transition to the Petri net. Label and rate can be provided
depending on the kind of Petri net.

Returns the ID of the transition
"""
add_transition!(p::AbstractPetriNet;kw...) = add_part!(p,:T;kw...)

""" Add `n` transitions to the Petri net. Label and rate can be provided
depending on the kind of Petri net.

Returns the ID of the transition
"""
add_transitions!(p::AbstractPetriNet,n;kw...) = add_parts!(p,:T,n;kw...)

""" add_input!(p::AbstractPetriNet,t,s;kw...)

Add an input relationship to the Petri net between the transition `t` and species `s`.

Returns the ID of the input relationship
"""
add_input!(p::AbstractPetriNet,t,s;kw...) = add_part!(p,:I;it=t,is=s,kw...)

""" add_inputs!(p::AbstractPetriNet,n,t,s;kw...)

Add input relationships to the Petri net between the transitions `t` and species `s`.

Returns the ID of the input relationship
"""
add_inputs!(p::AbstractPetriNet,n,t,s;kw...) = add_parts!(p,:I,n;it=t,is=s,kw...)

""" add_output!(p::AbstractPetriNet,t,s;kw...)

Add an output relationship to the Petri net between the transition `t` and species `s`.

Returns the ID of the input relationship
"""
add_output!(p::AbstractPetriNet,t,s;kw...) = add_part!(p,:O;ot=t,os=s,kw...)

""" add_outputs!(p::AbstractPetriNet,n,t,s;kw...)

Add output relationships to the Petri net between the transitions `t` and species `s`.

Returns the ID of the input relationship
"""
add_outputs!(p::AbstractPetriNet,n,t,s;kw...) = add_parts!(p,:O,n;ot=t,os=s,kw...)

""" Name of species

Note that this returns an index if labels are not present in the PetriNet
"""
sname(p::AbstractPetriNet,s) = (1:ns(p))[s]

""" Name of transition

Note that this returns an index if labels are not present in the PetriNet
"""
tname(p::AbstractPetriNet,t) = (1:nt(p))[t]

""" Names of species in  a Petri net

Note that this returns indices if labels are not present in the PetriNet
"""
snames(p::AbstractPetriNet) = map(s->sname(p, s), 1:ns(p))

""" Names of transitions in  a Petri net

Note that this returns indices if labels are not present in the PetriNet
"""
tnames(p::AbstractPetriNet) = map(t->tname(p, t), 1:nt(p))

# Note: although indexing makes this pretty fast, it is often faster to bulk-convert
# the PetriNet net into a transition matrix, if you are working with all of the transitions
""" Input relationships for a transition
"""
inputs(p::AbstractPetriNet,t) = subpart(p,incident(p,t,:it),:is)

""" Output relationships for a transition
"""
outputs(p::AbstractPetriNet,t) = subpart(p,incident(p,t,:ot),:os)

""" TransitionMatrices

This data structure stores the transition matrix of an AbstractPetriNet object.
This is primarily used for constructing the vectorfield representation of the
Petri net.
"""
struct TransitionMatrices
  input::Matrix{Int}
  output::Matrix{Int}
  TransitionMatrices(p::AbstractPetriNet) = begin
    input, output = zeros(Int,(nt(p),ns(p))), zeros(Int,(nt(p),ns(p)))
    for i in 1:ni(p)
      input[subpart(p,i,:it),subpart(p,i,:is)] += 1
    end
    for o in 1:no(p)
      output[subpart(p,o,:ot),subpart(p,o,:os)] += 1
    end
    new(input,output)
  end
end

""" PetriNet(tm::TransitionMatrices)

Constructs a PetriNet from its TransitionMatrices representation.
"""
PetriNet(tm::TransitionMatrices) = begin
  (m,n) = size(tm.input)
  p = PetriNet()
  add_species!(p,n)
  add_transitions!(p,m)
  for i in 1:m
    for j in 1:n
      add_inputs!(p,tm.input[i,j],i,j)
      add_outputs!(p,tm.output[i,j],i,j)
    end
  end
  p
end

valueat(x::Number, u, t) = x
valueat(f::Function, u, t) = try f(u,t) catch e f(t) end

""" vectorfield(pn::AbstractPetriNet)

Generates a Julia function which calculates the vectorfield of the Petri net
being simulated under the law of mass action.

The resulting function has a signature of the form `f!(du, u, p, t)` and can be
passed to the DifferentialEquations.jl solver package.
"""
vectorfield(pn::AbstractPetriNet) = begin
  tm = TransitionMatrices(pn)
  dt = tm.output - tm.input
  f(du,u,p,t) = begin
    rates = zeros(eltype(du),nt(pn))
    u_m = [u[sname(pn, i)] for i in 1:ns(pn)]
    p_m = [p[tname(pn, i)] for i in 1:nt(pn)]
    for i in 1:nt(pn)
      rates[i] = valueat(p_m[i],u,t) * prod(u_m[j] ^ tm.input[i,j] for j in 1:ns(pn))
    end
    for j in 1:ns(pn)
      du[sname(pn, j)] = sum(rates[i] * dt[i,j] for i in 1:nt(pn))
    end
    return du
  end
  return f
end

""" ACSet definition for a Petri net with labels on transitions and states.

See Catlab.jl documentation for description of the @present syntax.
"""
@present TheoryLabelledPetriNet <: TheoryPetriNet begin
  Name::AttrType

  tname::Attr(T, Name)
  sname::Attr(S, Name)
end

@abstract_acset_type AbstractLabelledPetriNet <: AbstractPetriNet
@acset_type LabelledPetriNetUntyped(TheoryLabelledPetriNet, index=[:it,:is,:ot,:os]) <: AbstractLabelledPetriNet
const LabelledPetriNet = LabelledPetriNetUntyped{Symbol}
const OpenLabelledPetriNetObUntyped, OpenLabelledPetriNetUntyped = OpenACSetTypes(LabelledPetriNetUntyped,:S)
const OpenLabelledPetriNetOb, OpenLabelledPetriNet = OpenLabelledPetriNetObUntyped{Symbol}, OpenLabelledPetriNetUntyped{Symbol}

Open(p::AbstractLabelledPetriNet) = OpenLabelledPetriNet(p, map(x->FinFunction([x], ns(p)), 1:ns(p))...)
Open(p::AbstractLabelledPetriNet, legs...) = begin
  s_idx = Dict(sname(p, s)=>s for s in 1:ns(p))
  OpenLabelledPetriNet(p, map(l->FinFunction(map(i->s_idx[i], l), ns(p)),legs)...)
end
Open(n, p::AbstractLabelledPetriNet, m) = Open(p, n, m)

""" LabelledPetriNet(n, ts::Vararg{Union{Pair,Tuple}})

Constructs a LabelledPetriNet object with state names as elements of `n` and
labelled transitions described by `ts`.
Transitions are given as `transition_name=>((input_states)=>(output_states))`.

A LabelledPetriNet modelling the SIR model with 3 states and 2 transitions can be
constructed as follows:
```@example
LabelledPetriNet([:S, :I, :R], :inf=>((:S,:I)=>(:I,:I)), :rec=>(:I=>:R))
```
"""
LabelledPetriNet(n, ts::Vararg{Union{Pair,Tuple}}) = begin
  p = LabelledPetriNet()
  n = vectorify(n)
  state_idx = state_dict(n)
  add_species!(p, length(n), sname=n)
  for (name,(ins,outs)) in ts
    i = add_transition!(p, tname=name)
    ins = vectorify(ins)
    outs = vectorify(outs)
    add_inputs!(p, length(ins), repeat([i], length(ins)), map(x->state_idx[x], collect(ins)))
    add_outputs!(p, length(outs), repeat([i], length(outs)), map(x->state_idx[x], collect(outs)))
  end
  p
end

# Reaction Nets
###############

""" ACSet definition for a Petri net with rates on transitions and
concentrations on states.

See Catlab.jl documentation for description of the @present syntax.
"""
@present TheoryReactionNet <: TheoryPetriNet begin
  Rate::AttrType
  Concentration::AttrType

  rate::Attr(T, Rate)
  concentration::Attr(S, Concentration)
end

@abstract_acset_type AbstractReactionNet <: AbstractPetriNet
@acset_type ReactionNet(TheoryReactionNet, index=[:it,:is,:ot,:os]) <: AbstractReactionNet
const OpenReactionNetOb, OpenReactionNet = OpenACSetTypes(ReactionNet,:S)

Open(p::ReactionNet{R,C}) where {R,C} = OpenReactionNet{R,C}(p, map(x->FinFunction([x], ns(p)), 1:ns(p))...)
Open(p::ReactionNet{R,C}, legs...) where {R,C} = OpenReactionNet{R,C}(p, map(l->FinFunction(l, ns(p)), legs)...)
Open(n, p::AbstractReactionNet, m) = Open(p, n, m)

""" ReactionNet{R,C}(n, ts::Vararg{Union{Pair,Tuple}}) where {R,C}

Constructs a ReactionNet object with state concentrations as elements of `n`
and transitions described by `ts`. `R` is the data type used to store rates and
`C` is the data type used to store concentrations.

Transitions are given as `transition_rate=>((input_states)=>(output_states))`.

A ReactionNet modelling the SIR model with 3 states and 2 transitions, an
initial population of 10 susceptible, 1 infected, 0 recovered and an infection
rate of 0.5 and recovery rate of 0.1 can be
constructed as follows:
```@example
ReactionNet{Float64, Float64}([10,1,0], 0.5=>((1,2)=>(2,2)), 0.1=>(2=>3))
```
"""
ReactionNet{R,C}(n, ts::Vararg{Union{Pair,Tuple}}) where {R,C} = begin
  p = ReactionNet{R,C}()
  add_species!(p, length(n), concentration=n)
  for (i, (rate,(ins,outs))) in enumerate(ts)
    i = add_transition!(p, rate=rate)
    ins = vectorify(ins)
    outs = vectorify(outs)
    add_inputs!(p, length(ins), repeat([i], length(ins)), collect(ins))
    add_outputs!(p, length(outs), repeat([i], length(outs)), collect(outs))
  end
  p
end

concentration(p::AbstractReactionNet,s) = subpart(p,s,:concentration)
rate(p::AbstractReactionNet,t) = subpart(p,t,:rate)

concentrations(p::AbstractReactionNet) = map(s->concentration(p, s), 1:ns(p))
rates(p::AbstractReactionNet) = map(t->rate(p, t), 1:nt(p))

""" ACSet definition for a ReactionNet with labels on transitions and states.

See Catlab.jl documentation for description of the @present syntax.
"""
@present TheoryLabelledReactionNet <: TheoryReactionNet begin
  Name::AttrType

  tname::Attr(T, Name)
  sname::Attr(S, Name)
end

@abstract_acset_type AbstractLabelledReactionNet <: AbstractPetriNet
@acset_type LabelledReactionNetUntyped(TheoryLabelledReactionNet, index=[:it,:is,:ot,:os]) <: AbstractLabelledReactionNet
const LabelledReactionNet{R,C} = LabelledReactionNetUntyped{R,C,Symbol}
const OpenLabelledReactionNetObUntyped, OpenLabelledReactionNetUntyped = OpenACSetTypes(LabelledReactionNetUntyped,:S)
const OpenLabelledReactionNetOb{R,C} = OpenLabelledReactionNetObUntyped{R,C,Symbol}
const OpenLabelledReactionNet{R,C} = OpenLabelledReactionNetUntyped{R,C,Symbol}

Open(p::LabelledReactionNet{R,C}) where {R,C} = OpenLabelledReactionNet{R,C}(p, map(x->FinFunction([x], ns(p)), 1:ns(p))...)
Open(p::LabelledReactionNet{R,C}, legs...) where {R,C} = begin
  s_idx = Dict(sname(p, s)=>s for s in 1:ns(p))
  OpenLabelledReactionNet{R,C}(p, map(l->FinFunction(map(i->s_idx[i], l), ns(p)), legs)...)
end
Open(n, p::AbstractLabelledReactionNet, m) = Open(p, n, m)

# Ex. LabelledReactionNet{Number, Int}((:S=>990, :I=>10, :R=>0), (:inf, .3/1000)=>((:S, :I)=>(:I,:I)), (:rec, .2)=>(:I=>:R))
""" LabelledReactionNet{R,C}(n, ts::Vararg{Union{Pair,Tuple}}) where {R,C}

Constructs a LabelledReactionNet object with labelled state concentrations as elements of `n`
and labelled transitions described by `ts`. `R` is the data type used to store rates and
`C` is the data type used to store concentrations.

Transitions are given as `(t_name=>t_rate)=>((input_states)=>(output_states))`.

A LabelledReactionNet modelling the SIR model with 3 states and 2 transitions, an
initial population of 10 susceptible, 1 infected, 0 recovered and an infection
rate of 0.5 and recovery rate of 0.1 can be
constructed as follows:
```@example
ReactionNet{Float64, Float64}([:S=>10,:I=>1,:R=>0], (:inf=>0.5)=>((1,2)=>(2,2)), (:rec=>0.1)=>(2=>3))
```
"""
LabelledReactionNet{R,C}(n, ts::Vararg{Union{Pair,Tuple}}) where {R,C} = begin
  p = LabelledReactionNet{R,C}()
  n = vectorify(n)
  states = map(first, collect(n))
  concentrations = map(last, collect(n))
  state_idx = state_dict(states)
  add_species!(p, length(states), concentration=concentrations, sname=states)
  for (i, ((name,rate),(ins,outs))) in enumerate(ts)
    i = add_transition!(p,rate=rate, tname=name)
    ins = vectorify(ins)
    outs = vectorify(outs)
    add_inputs!(p, length(ins), repeat([i], length(ins)), map(x->state_idx[x], collect(ins)))
    add_outputs!(p, length(outs), repeat([i], length(outs)), map(x->state_idx[x], collect(outs)))
  end
  p
end

sname(p::Union{AbstractLabelledPetriNet, AbstractLabelledReactionNet},s) = subpart(p,s,:sname)
tname(p::Union{AbstractLabelledPetriNet, AbstractLabelledReactionNet},t) = subpart(p,t,:tname)

# Interoperability between different types

PetriNet(pn::AbstractPetriNet) = begin
  pn′ = PetriNet()
  copy_parts!(pn′, pn)
  pn′
end

LabelledPetriNet(pn::Union{AbstractLabelledPetriNet, AbstractLabelledReactionNet}) = begin
  pn′ = LabelledPetriNet()
  copy_parts!(pn′, pn)
  pn′
end

LabelledPetriNet(pn::AbstractPetriNet, snames, tnames) = begin
  pn′ = LabelledPetriNet()
  copy_parts!(pn′, pn)
  map(k->set_subpart!(pn′, k, :sname, snames[k]), keys(snames))
  map(k->set_subpart!(pn′, k, :tname, tnames[k]), keys(tnames))
  pn′
end

ReactionNet{R,C}(pn::Union{AbstractReactionNet, AbstractLabelledReactionNet}) where {R, C} = begin
  pn′ = ReactionNet{R,C}()
  copy_parts!(pn′, pn)
  pn′
end

ReactionNet{R,C}(pn::AbstractPetriNet, concentrations, rates) where {R, C} = begin
  pn′ = ReactionNet{R,C}()
  copy_parts!(pn′, pn)
  map(k->set_subpart!(pn′, k, :concentration, concentrations[k]), keys(concentrations))
  map(k->set_subpart!(pn′, k, :rate, rates[k]), keys(rates))
  pn′
end

LabelledReactionNet{R,C}(pn::AbstractLabelledReactionNet) where {R, C} = begin
  pn′ = LabelledReactionNet{R,C}()
  copy_parts!(pn′, pn)
  pn′
end

LabelledReactionNet{R,C}(pn::AbstractPetriNet, s_labels, t_labels, concentrations, rates) where {R, C} = begin
  pn′ = LabelledReactionNet{R,C}()
  copy_parts!(pn′, pn)
  map(k->set_subpart!(pn′, k, :sname, s_labels[k]), keys(s_labels))
  map(k->set_subpart!(pn′, k, :tname, t_labels[k]), keys(t_labels))
  map(k->set_subpart!(pn′, k, :concentration, concentrations[k]), keys(concentrations))
  map(k->set_subpart!(pn′, k, :rate, rates[k]), keys(rates))
  pn′
end

LabelledReactionNet{R,C}(pn::Union{AbstractPetriNet}, states, transitions) where {R, C} = begin
  pn′ = LabelledReactionNet{R,C}()
  copy_parts!(pn′, pn)
  for (i, (k, v)) in enumerate(states)
    set_subpart!(pn′, i, :sname, k)
    set_subpart!(pn′, i, :concentration, v)
  end
  for (i, (k, v)) in enumerate(transitions)
    set_subpart!(pn′, i, :tname, k)
    set_subpart!(pn′, i, :rate, v)
  end
  pn′
end

""" Concentration of a ReactionNet
"""
concentration(p::AbstractLabelledReactionNet,s) = subpart(p,s,:concentration)

""" Rate of a RectionNet
"""
rate(p::AbstractLabelledReactionNet,t) = subpart(p,t,:rate)

""" All concentrations of a ReactionNet
"""
concentrations(p::AbstractLabelledReactionNet) = begin
  snames = [sname(p, s) for s in 1:ns(p)]
  LVector(;[(snames[s]=>concentration(p, s)) for s in 1:ns(p)]...)
end

""" All rates of a ReactionNet
"""
rates(p::AbstractLabelledReactionNet) = begin
  tnames = [tname(p, s) for s in 1:nt(p)]
  LVector(;[(tnames[t]=>rate(p, t)) for t in 1:nt(p)]...)
end

include("interoperability.jl")
include("visualization.jl")
include("Epidemiology.jl")
include("BilayerNetworks.jl")
include("ModelComparison.jl")
end

max_stars_repo_name: AlgebraicJulia/AlgebraicPetri

max_issues_repo_name: AlgebraicJulia/AlgebraicPetri

max_forks_repo_name: AlgebraicJulia/AlgebraicPetri


text:
push!(LOAD_PATH, "../src")
push!(LOAD_PATH, "../Distances")

using Distances, Aircraft_module, Constants_module

function check_NMAC(All_AC::Array{Aircraft,1})
	All_Coords = Matrix{Float64}(2, 0)

	R = Matrix{Float64}(length(All_AC), length(All_AC))
	for ac in All_AC
		All_Coords = hcat(All_Coords, [ac.dyn.x; ac.dyn.y])
	end

	pairwise!(R, Euclidean(), All_Coords) #inline

	# iterate through the upper right side of R (not including diag)
	# so that not double count
	TOTAL_num_colli = 0
	for col = 2 : length(All_AC)
		for row = 1 : col - 1
			if R[row, col] <= collision_range
				TOTAL_num_colli = TOTAL_num_colli + 1
			end
		end
	end

	return TOTAL_num_colli
end

function check_desti(ac::Aircraft)
	dist = sqrt((ac.dyn.x - ac.x_dest) ^ 2 + (ac.dyn.y - ac.y_dest) ^ 2)
	if dist <= dest_range
		return true
	else
		return false
	end
end

function check_in_airspace(ac::Aircraft)
	if ac.dyn.x >= 0 && ac.dyn.x <= airspace_dim  &&
		ac.dyn.y >= 0 && ac.dyn.y <= airspace_dim
		return true
	else
		return false
	end
end

function check_colli(new_ac::Aircraft, All_AC::Array{Aircraft,1})
	for ac in All_AC
		dist = sqrt((ac.dyn.x - new_ac.dyn.x) ^ 2 + (ac.dyn.y - new_ac.dyn.y) ^ 2)
		if dist <= collision_range
			break
			return true
		end
	end
end

max_stars_repo_name: parachutel/DeepCorrectionCAS


text:
"A stateful browser counter."
mutable struct Countbrowser;value::Int;end
(c::Countbrowser)() =COUNTBROWSER.value += 1
"For next value: COUNTBROWSER(). For current value: COUNTBROWSER.value"
COUNTBROWSER = Countbrowser(0)

"Get application path for developer applications"
function fwhich(s)
    fi = ""
    if Sys.iswindows()
        try
            fi = split(read(`where.exe $s`, String), "\r\n")[1]
            if !isfile(fi)
                fi = ""
            end
        catch
            fi =""
        end
    else
        try
            fi = readchomp(`which $s`)
        catch
            fi =""
        end
    end
    fi
end
function browser_path_unix_apple(shortname)
    trypath = ""
    if shortname == "chrome"
        if Sys.isapple()
            return "Google Chrome"
        else
            return "google-chrome"
        end
    end
    if shortname == "firefox"
        return "firefox"
    end
    if shortname == "safari"
        if Sys.isapple()
            return "safari"
        else
            return ""
        end
    end
    if shortname == "phantomjs"
        return fwhich(shortname)
    end
    return ""
end
function browser_path_windows(shortname)
    # windows accepts English paths anywhere.
    # Forward slash is acceptable too.
    trypath = ""
    homdr = ENV["HOMEDRIVE"]
    path32 = homdr * "/Program Files (x86)/"
    path64 = homdr * "/Program Files/"
    if shortname == "chrome"
        trypath = path64 * "Google/Chrome/Application/chrome.exe"
        isfile(trypath) && return trypath
        trypath = path32 * "Google/Chrome/Application/chrome.exe"
        isfile(trypath) && return trypath
    end
    if shortname == "firefox"
        trypath = path64 * "Mozilla Firefox/firefox.exe"
        isfile(trypath) && return trypath
        trypath = path32 * "Mozilla Firefox/firefox.exe"
        isfile(trypath) && return trypath
    end
    if shortname == "safari"
        trypath = path64 * "Safari/Safari.exe"
        isfile(trypath) && return trypath
        trypath = path32 * "Safari/Safari.exe"
        isfile(trypath) && return trypath
    end
    if shortname == "iexplore"
        trypath = path64 * "Internet Explorer/iexplore.exe"
        isfile(trypath) && return trypath
        trypath = path32 * "Internet Explorer/iexplore.exe"
        isfile(trypath) && return trypath
    end
    if shortname == "phantomjs"
        return fwhich(shortname)
    end
    return ""
end
"Constructs launch command"
function launch_command(shortbrowsername; url = DEFAULT_REDIRECT_URI, privatemode = false)
    if Sys.iswindows()
        pt = browser_path_windows(shortbrowsername)
    else
        pt = browser_path_unix_apple(shortbrowsername)
    end
    pt == "" && return ``
    if shortbrowsername == "iexplore"
        prsw = "-private"
    else
        prsw = "--incognito"
    end
    if shortbrowsername == "phantomjs"
        if isdefined(:SRCPATH)
            script = joinpath(SRCPATH, "phantom.js")
        else
            script = "phantom.js"
        end
        return Cmd(`$pt $script $url`)
    else
        if Sys.iswindows()
            if privatemode
                return Cmd( [ pt, url, prsw])
            else
                return Cmd( [ pt, url])
            end
        else
            if Sys.isapple()
                if privatemode
                    return Cmd(`open --fresh -n $url -a $pt --args $prsw`)
                else
                    return Cmd(`open --fresh -n $url -a $pt`)
                end
            elseif Sys.islinux() || Sys.isbsd()
                return Cmd(`xdg-open $(url) $pt`)
            end
        end
    end
end

function open_on_page(shortbrowsername; url = DEFAULT_REDIRECT_URI)
    id = "open_on_page "
    dmc = launch_command(shortbrowsername; url)
    if dmc == ``
        @debug id, "Could not find ", shortbrowsername
        return false
    else
        try
            if shortbrowsername == "phantomjs"
                # Run enables text output of phantom messages in the REPL. In Windows
                # standalone REPL, run will freeze the main thread if not run async.
                @async run(dmc)
            else
                run(dmc, wait = false)
            end
        catch
            @debug id, "Failed to spawn ", shortbrowsername
            return false
        end
    end
    return true
end

"Try to open one browser from BROWSERS.
In some cases we expect an immediate indication
of failure, for example when the corresponding browser
is not found on the system. In other cases, we will
just wait in vain. In those cases,
call this function again after a reasonable timeout.
The function remembers which browsers were tried before.
"
function open_a_browser(;url = DEFAULT_REDIRECT_URI)
    id = "open_a_browser "
    if COUNTBROWSER.value > length(BROWSERS)
        return false
    end
    success = false
    b = ""
    while COUNTBROWSER.value < length(BROWSERS) && !success
        b = BROWSERS[COUNTBROWSER()]
        @debug id, "Trying to launch browser no. ", COUNTBROWSER.value, ", ",  b
        success = open_on_page(b, url = url)
    end
    if success
        @debug id, b, " on ", url, " seems to work."
    end
    success, b
end
function open_browser(b;url = DEFAULTURL)
    if b ∉ BROWSERS
        throw(ArgumentError("$b not in $BROWSERS"))
    end
    open_on_page(b, url = url)
end


text:
@testset "MPC Methods" begin

    N = 10
    dt = 0.1
    p = 2
    model = UnicycleGame(p=p)
    opts = Options(inner_print=0.0, outer_print=0.0)
    roadway_opts = HighwayRoadwayOptions()
    roadway = build_roadway(roadway_opts)
    player = Vector{Player}(undef,p)
    player[1] = Player(model, roadway.lane[1])
    player[2] = Player(model, roadway.lane[2])
    sce = Scenario(model, roadway, player)
    vis = Visualizer()
    prob = GameProblem(N,dt,sce,opts)
    m_stats = MPCStatistics()
    m_opts = MPCOptions(M=5)
    x0 = get_state(sce)
    simulate_MPC!(vis, prob, x0, sce, m_stats, m_opts)

end


text:
const BlockTriPertToeplitz{T} = BlockMatrix{T,Tridiagonal{Matrix{T},Vcat{Matrix{T},1,Tuple{Vector{Matrix{T}},Fill{Matrix{T},1,Tuple{OneToInf{Int}}}}}},
                                        NTuple{2,BlockedUnitRange{Vcat{Int,1,Tuple{Vector{Int},InfStepRange{Int,Int}}}}}}

const BlockTridiagonalToeplitzLayout = BlockLayout{TridiagonalToeplitzLayout,DenseColumnMajor}

function BlockTridiagonal(adjA::Adjoint{T,BlockTriPertToeplitz{T}}) where T
    A = parent(adjA)
    BlockTridiagonal(Matrix.(adjoint.(A.blocks.du)),
                     Matrix.(adjoint.(A.blocks.d)),
                     Matrix.(adjoint.(A.blocks.dl)))
end

for op in (:-, :+)
    @eval begin
        function $op(A::BlockTriPertToeplitz{T}, λ::UniformScaling) where T
            TV = promote_type(T,eltype(λ))
            BlockTridiagonal(Vcat(convert.(AbstractVector{Matrix{TV}}, A.blocks.dl.args)...),
                             Vcat(convert.(AbstractVector{Matrix{TV}}, broadcast($op, A.blocks.d, Ref(λ)).args)...),
                             Vcat(convert.(AbstractVector{Matrix{TV}}, A.blocks.du.args)...))
        end
        function $op(λ::UniformScaling, A::BlockTriPertToeplitz{V}) where V
            TV = promote_type(eltype(λ),V)
            BlockTridiagonal(Vcat(convert.(AbstractVector{Matrix{TV}}, broadcast($op, A.blocks.dl.args))...),
                             Vcat(convert.(AbstractVector{Matrix{TV}}, broadcast($op, Ref(λ), A.blocks.d).args)...),
                             Vcat(convert.(AbstractVector{Matrix{TV}}, broadcast($op, A.blocks.du.args))...))
        end
        $op(adjA::Adjoint{T,BlockTriPertToeplitz{T}}, λ::UniformScaling) where T = $op(BlockTridiagonal(adjA), λ)
        $op(λ::UniformScaling, adjA::Adjoint{T,BlockTriPertToeplitz{T}}) where T = $op(λ, BlockTridiagonal(adjA))
    end
end

*(a::AbstractVector, b::AbstractFill{<:Any,2,Tuple{OneTo{Int},OneToInf{Int}}}) = ApplyArray(*,a,b)

sizes_from_blocks(A::Diagonal, ::NTuple{2,OneToInf{Int}}) = size.(A.diag, 1), size.(A.diag,2)
sizes_from_blocks(A::Tridiagonal, ::NTuple{2,OneToInf{Int}}) = size.(A.d, 1), size.(A.d,2)
sizes_from_blocks(A::LazyBandedMatrices.Tridiagonal, ::NTuple{2,OneToInf{Int}}) = size.(A.d, 1), size.(A.d,2)
sizes_from_blocks(A::Bidiagonal, ::NTuple{2,OneToInf{Int}}) = size.(A.dv, 1), size.(A.dv,2)
sizes_from_blocks(A::LazyBandedMatrices.Bidiagonal, ::NTuple{2,OneToInf{Int}}) = size.(A.dv, 1), size.(A.dv,2)

axes_print_matrix_row(_, io, X, A, i, ::AbstractVector{<:PosInfinity}, sep) = nothing
axes_print_matrix_row(::NTuple{2,BlockedUnitRange}, io, X, A, i, ::AbstractVector{<:PosInfinity}, sep) = nothing

function BlockSkylineSizes(A::BlockTriPertToeplitz, (l,u)::NTuple{2,Int})
    N = max(length(A.blocks.du.args[1])+1,length(A.blocks.d.args[1]),length(A.blocks.dl.args[1]))
    block_sizes = Vector{Int}(undef, N) # assume square
    block_starts = BandedMatrix{Int}(undef, (N+l,N),  (l,u))
    block_strides = Vector{Int}(undef, N)
    for J=1:N
        block_starts[max(1,J-u),J] = J == 1 ? 1 :
                            block_starts[max(1,J-1-u),J-1]+block_sizes[J-1]*block_strides[J-1]

        for K=max(1,J-u)+1:J+l
            block_starts[K,J] = block_starts[K-1,J]+size(A[Block(K-1,J)],1)
        end
        block_strides[J] = block_starts[J+l,J] + size(A[Block(J+l,J)],1) - block_starts[max(1,J-u),J]
        block_sizes[J] = size(A[Block(J,J)],2)
    end

    block_stride∞ = 0
    for K=max(1,N+1-u):N+1+l
        block_stride∞ += size(A[Block(K,N+1)],1)
    end
    block_size∞ = size(A[Block(N+1,N+1)],2)

    bs∞ = fill(block_starts[max(1,N-u),N]+block_strides[N]*size(A[Block(N,N)],2):block_stride∞*block_size∞:∞, l+u+1)
    for k=2:l+u+1
        bs∞[k] = bs∞[k-1] .+ size(A[Block(N+1-u+k-1,N+1)],1)
    end

    BlockSkylineSizes(axes(A),
                        _BandedMatrix(Hcat(block_starts.data, Vcat(adjoint.(bs∞)...)), ℵ₀, l, u),
                        Vcat(block_strides, Fill(block_stride∞,∞)),
                        Fill(l,∞),Fill(u,∞))
end

function BlockBandedMatrix(A::BlockTriPertToeplitz{T}, (l,u)::NTuple{2,Int}) where T
    data = T[]
    append!(data,vec(A[Block.(1:1+l),Block(1)]))
    N = max(length(A.blocks.du.args[1])+1,length(A.blocks.d.args[1]),length(A.blocks.dl.args[1]))
    for J=2:N
        append!(data, vec(A[Block.(max(1,J-u):J+l),Block(J)]))
    end
    tl = mortar(Fill(vec(A[Block.(max(1,N+1-u):N+1+l),Block(N+1)]),∞))

    B = _BlockSkylineMatrix(Vcat(data,tl), BlockSkylineSizes(A, (l,u)))
end

BlockBandedMatrix(A::BlockTriPertToeplitz) = BlockBandedMatrix(A, blockbandwidths(A))


text:
# ---------------------------------------------------------------------------- #
#
#   mesherBAMG.jl
#
#   BAMG-specific read and write
#
#   λυτέος
#   Fall 2017
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

"""
    MesherBAMG

MesherBAMG type:
Type for BAMG meshes.
"""
struct MesherBAMG <: Mesher

end

"""
    writeMsh( mesh::Mesh2D, mesher::MesherBAMG, flname::String )

Writes 2D mesh in `mesh` to `flname` in format used by BAMG.
"""
function writeMsh( mesh::Mesh2D, mesher::MesherBAMG, flname::String )

    np = size(mesh.p,1)

    fid = open( flname, "w" )

    @printf( fid, "MeshVersionFormatted 0\n\n" )

    @printf(fid, "Dimension 2\n\n")

    # Vertices
    @printf(fid, "Vertices\n%i\n", np)
    for jj = 1:size(mesh.p,1)
        @printf(fid, "%16.15e\t%16.15e\t%i\n", mesh.p[jj,1], mesh.p[jj,2], 1)
    end
    @printf(fid,"\n")

    # Edges
    nbound = size(mesh.fb,1)
    @printf(fid, "Edges\n%i\n",nbound)
    for ii = 1:nbound
        indf = mesh.fb[ii,1]
        @printf(fid, "%i\t%i\t%i\n", mesh.f[indf,1], mesh.f[indf,2], mesh.fb[ii,2])
    end
    @printf(fid,"\n")

    # Triangles
    @printf(fid, "Triangles\n%i\n",size(mesh.t,1))
    for jj = 1:size(mesh.t,1)
        @printf(fid, "%i\t%i\t%i\t%i\n", mesh.t[jj,1], mesh.t[jj,2], mesh.t[jj,3], 1)
    end
    @printf(fid,"\n")

    close(fid)

    # TODO: Corners

end

"""
    readBAMG( flname::String )

Reads (2D) BAMG mesh from `flname`.
"""
function readBAMG( flname::String )

    dim = 2

    # Open file
    fid = open( flname, "r" )

    # Check dimension
    for il in 1:3
        x = readline(fid)
    end
    x = readline(fid)
    if parse(Int64, x) != dim
        error(" BAMG read-in: dimension should be 2")
    end

    # Read vertices
    while isempty(x) || x[1] != 'V'
        x = readline(fid)
    end
    x  = readline(fid)
    np = parse( Int64, x )
    p  = fill( 0.0, np, 2 )
    for il in 1:np
        x  = readline(fid)
        xx = split(x)
        p[il,1] = parse( Float64, xx[1] )
        p[il,2] = parse( Float64, xx[2] )
    end

    # Read boundary
    while isempty(x) || x[1] != 'E'
        x = readline(fid)
    end
    x   = readline(fid)
    nb  = parse( Int64, x )
    bel = fill( 0, nb, 3 )
    for il in 1:nb
        x  = readline(fid)
        xx = split(x)
        bel[il,1] = parse( Int64, xx[1] )
        bel[il,2] = parse( Int64, xx[2] )
        bel[il,3] = parse( Int64, xx[3] )
    end

    # Read triangles
    while isempty(x) || x[1] != 'T'
        x = readline(fid)
    end
    x  = readline(fid)
    nt = parse( Int64, x )
    t  = fill( 0, nt, 3 )
    for il in 1:nt
        x  = readline(fid)
        xx = split(x)
        t[il,1] = parse( Int64, xx[1] )
        t[il,2] = parse( Int64, xx[2] )
        t[il,3] = parse( Int64, xx[3] )
    end

    close(fid)

    # Add boundary tags
    nbounds = maximum( bel[:,3] )
    tags = Vector{String}( nbounds )
    for bb in 1:nbounds
        tags[bb] = string("bound", bb)
    end

    return( p, t, bel, tags )

end

"""
    runMesher( mesher::MesherBAMG, flMesh::String, flMetric::String, flMeshNew::String )

Runs BAMG mesher.
"""
function runMesher( mesher::MesherBAMG, flMesh::String, flMetric::String, flMeshNew::String )

    run( `bamg -b $flMesh -M $flMetric -o $flMeshNew` )

end


text:
using DelimitedFiles

function solve(input,n_generation)
    raw = readdlm(input,'\n')
    initial_state = parse(raw[1])
    spread = Dict{String,String}(map(parse, raw[2:end]))
    state = string("...",initial_state,"...")
    state_sum = Int64(0)
    zero_index=4
    sum_diffs = zeros(Int16,42)
    current_gen = 1
    for generation in 1:n_generation
        next_state = Array{String,1}()
        for pot in 3:length(state)-2
            note = state[pot-2:pot+2]
            if note in keys(spread)
                push!(next_state,spread[note])
            else
                push!(next_state,".")
            end
        end
        next_state = string("...",join(next_state),"...")
        zero_index+=1
        next_state_sum = sum((findall(x->x=='#',next_state) .-zero_index))
        sum_diff = next_state_sum-state_sum
        popfirst!(sum_diffs)
        push!(sum_diffs,sum_diff)
        state = next_state
        state_sum = next_state_sum
        current_gen += 1
        if length(unique(sum_diffs))==1
            break
        end
    end
    if length(unique(sum_diffs))==1
        current_gen -= 1
        state_sum += sum_diffs[1]*(n_generation-current_gen)
    end
    @show state_sum
end

function parse(s)
    m = match(r"(initial state: (?<initialstate>.+)|(?<left>.+) => (?<right>.+))", s)
    initial_state = m[:initialstate]
    if initial_state!= nothing
        return initial_state
    else
        return (m[:left],m[:right])
    end
end

solve("input_test.txt",20)
solve("input.txt",20)
solve("input.txt",50000000000)


text:
const Double = Float64
# For no of courses
const ncInt = UInt8
# For time variables
const TimeInt = UInt8
# For indexing types
const TypeInt = UInt16
# For indexing colleges
const CollInt = UInt8
# For indexing school levels
const SchoolInt = UInt8
# For indexing grid points
const GridInt = UInt16

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


text:
using BenchmarkTools

const BENCHMARKS_DIR = "../benchmark/"
cp(joinpath(BENCHMARKS_DIR, "benchmark_utils.jl"), joinpath(@__DIR__, "benchmark_utils.jl"), force=true)

function run_benchmark(replace_strings, benchmark_name, module_suffix="")
    benchmark_filepath = joinpath(BENCHMARKS_DIR, "benchmark_" * benchmark_name * ".jl")
    txt = read(benchmark_filepath, String)

    for strs in replace_strings
        txt = replace(txt, strs[1] => strs[2])
    end

    test_script_filepath = benchmark_name * "_benchmark_test.jl"

    open(test_script_filepath, "w") do f
        write(f, "module Test_$benchmark_name" * "_$module_suffix\n")
        write(f, txt)
        write(f, "\nend # module")
    end

    try
        include(test_script_filepath)
    catch err
        @error sprint(showerror, err)
        rm(test_script_filepath)
        return false
    end

    rm(test_script_filepath)
    return true
end

@testset "Performance benchmarks" begin
    @info "Testing performance benchmarks..."

    @testset "Performance benchmark scripts" begin
        @info "  Testing performance benchmark scripts..."

        @testset "Static ocean benchmark" begin
            @info "    Running static ocean benchmarks..."

            replace_strings = [
                ("Ns = [(32, 32, 32), (64, 64, 64), (128, 128, 128), (256, 256, 256)]",
                 "Ns = [(16, 16, 16)]")
            ]

            @test run_benchmark(replace_strings, "static_ocean")
        end

        @testset "Channel benchmark" begin
            @info "    Running channel benchmarks..."

            replace_strings = [
                ("Ns = [(32, 32, 32), (64, 64, 64), (128, 128, 128), (256, 256, 256)]",
                 "Ns = [(16, 16, 16)]")
            ]

            @test run_benchmark(replace_strings, "channel")
        end

        @testset "Turbulence closures benchmark" begin
            @info "    Running turbulence closures benchmark..."

            replace_strings = [
                ("Ns = [(32, 32, 32), (128, 128, 128)]",
                 "Ns = [(16, 16, 16)]")
            ]

            @test run_benchmark(replace_strings, "turbulence_closures")
        end

        @testset "Tracers benchmark" begin
            @info "    Running tracers benchmark..."

            replace_strings = [
                ("(32, 32, 32)", "(16, 16, 16)"),
                ("(256, 256, 256)", "(16, 16, 16)"),
                ("test_cases = [(0, 0), (0, 1), (0, 2), (1, 0), (2, 0), (2, 3), (2, 5), (2, 10)]",
                 "test_cases = [(0, 0), (2, 0), (2, 3)]")
            ]

            @test run_benchmark(replace_strings, "tracers")
        end
    end

    @testset "Selected performance benchmarks" begin
        @info "  Running selected performance benchmarks..."

        for arch in archs
            sizes = [(16, 16, 16), (32, 32, 32)]
            for sz in sizes
                grid = RegularCartesianGrid(size=sz, extent=(1, 1, 1))
                model = IncompressibleModel(architecture=arch, grid=grid)

                @info "Benchmarking [$arch, Float64, $sz]..."
                b = @benchmark time_step!($model, 1)
                display(b)
                println()

                @test model isa IncompressibleModel
            end
        end
    end
end

rm(joinpath(@__DIR__, "benchmark_utils.jl"))


text:
module TestUtilities

using Test
using MLJBase

@test MLJBase.finaltypes(Union{Missing,Int}) == [Union{Missing,Int64}]
@test MLJBase.finaltypes(Float64) == [Float64]

abstract type Foo end
struct Bar <: Foo end
struct Baz <: Foo end

@test MLJBase.finaltypes(Foo) == [Bar, Baz]

end # module
true


text:
using Test, LinearAlgebra, SparseArrays, Random
using QDLDL
rng = Random.MersenneTwister(1312)

@testset "update values" begin

    # create random KKT system
    nz = 100
    nc = 70
    H = sprand(nz,nz,0.05);
    H = H'*H + I
    b = randn(nz + nc)
    A1 = sprand(nc,nz,0.8)
    K1 = [H A1';A1 -1e-3*I]

    # get triu of K1 and create QDLDL struct
    triuK1 = triu(K1);
    F = qdldl(triuK1)

    # compare with backslash
    @test norm(K1\b - F\b,Inf) < 1e-12

    # create a new KKT system with the same sparsity pattern
    A2 = copy(A1)
    A2.nzval .= randn(length(A2.nzval))
    K2 = [H A2';A2 -1e-7*I]
    triuK2 = triu(K2)

    # update factorization of F in place (non allocating)
    update_values!(F,1:length(triuK2.nzval),triuK2.nzval)
    refactor!(F)

    # compare with backslash using the refactored F
    @test norm(K2\b - F\b,Inf)<1e-12

end


text:
mutable struct Impact{T} <: Contact{T}
    ainv3::Adjoint{T,SVector{3,T}}
    p::SVector{3,T}
    offset::SVector{3,T}

    function Impact(body::Body{T}, normal::AbstractVector; p::AbstractVector = zeros(3),  offset::AbstractVector = zeros(3)) where T
        # Derived from plane equation a*v1 + b*v2 + distance*v3 = p - offset
        V1, V2, V3 = orthogonalcols(normal) # gives two plane vectors and the original normal axis
        A = [V1 V2 V3]
        Ainv = inv(A)
        ainv3 = Ainv[3,SA[1; 2; 3]]'

        new{T}(ainv3, p, offset), body.id
    end
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, constraint::Impact{T}) where {T}
    summary(io, constraint)
    println(io,"")
    println(io, " ainv3:  "*string(constraint.ainv3))
    println(io, " offset: "*string(constraint.offset))
end

## Additional force for friction default
@inline additionalforce(bound::Impact{T}, x::AbstractVector, q::UnitQuaternion) where T = szeros(T, 6)


text:
using IrisFlower
using Test

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


text:
# Gradient tests for adjoint extended modeling
# Author: Ziyi Yin, ziyi.yin@gatech.edu
# Date: Aug 2020
#
using JUDI4Flux
using JUDI, SegyIO, JOLI, Flux
using Test, ImageFiltering, LinearAlgebra, Random

Random.seed!(11)

# Set up model structure
n = (120, 100)   # (x,y,z) or (x,z)
d = (10., 10.)
o = (0., 0.)

# Velocity [km/s]
v = ones(Float32,n) .+ 0.4f0
v[:, Int(round(end/2)):end] .= 4f0

# Slowness squared [s^2/km^2]
m = (1f0 ./ v).^2
m0 = imfilter(m, Float32.(Kernel.gaussian(10)))

# Setup info and model structure
nsrc = 4	# number of sources
model = Model(n, d, o, m)
model0 = Model(n, d, o, m0)
m = reshape(m,n[1],n[2],1,1)
m0 = reshape(m0,n[1],n[2],1,1)

# Set up receiver geometry
nxrec = 120
xrec = range(50f0, stop=1150f0, length=nxrec)
yrec = 0f0
zrec = range(50f0, stop=50f0, length=nxrec)

# receiver sampling and recording time
time = 1000f0   # receiver recording time [ms]
dt = 1f0    # receiver sampling interval [ms]
nt = Int(time/dt)+1
# Set up receiver structure
recGeometry = Geometry(xrec, yrec, zrec; dt=dt, t=time, nsrc=nsrc)

# setup wavelet
f0 = 0.01f0     # MHz
wavelet = ricker_wavelet(time, dt, f0)

# Set up info structure for linear operators
ntComp = get_computational_nt(recGeometry, model)
info = Info(prod(n), nsrc, ntComp)

# return linearized data as Julia array
opt = Options(return_array=true, sum_padding=true)

# Linear operators
Pr = judiProjection(info, recGeometry)
A_inv = judiModeling(info, model; options=opt)
A_inv0 = judiModeling(info, model0; options=opt)
Pw = judiLRWF(info, wavelet)
F = Pr*A_inv*adjoint(Pw)
F0 = Pr*A_inv0*adjoint(Pw)

function GenSimSourceMulti(xsrc_index, zsrc_index, nsrc, n)
	weights = zeros(Float32, n[1], n[2], 1, nsrc)

	for j=1:nsrc
        weights[xsrc_index[j], zsrc_index[j], 1, j] = 1f0
    end

    return weights
end

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

G = ExtendedQForward(F)

function misfit_objective(d_obs, q0, m0, G)

	# Data residual and source residual
  	r_data   = G(q0, m0) - d_obs

	# Function value
	fval = .5f0 *sum(r_data.^2)

	print("Function value: ", fval, "\n")
	return fval
end

if nsrc == 1
    xsrc_index = Int.(round.((range(50,stop=50,length=nsrc))))   # indices where sim source weights are 1
    zsrc_index = Int.(round.((range(10,stop=10,length=nsrc))))
else
    xsrc_index = Int.(round.((range(1,stop=n[1]-1,length=nsrc))))   # indices where sim source weights are 1
    zsrc_index = Int.(round.((range(10,stop=10,length=nsrc))))
end

q = GenSimSourceMulti(xsrc_index, zsrc_index, nsrc, model.n);
q0 = deepcopy(q)
for i = 1:nsrc
    q0[:,:,1,i] = imfilter(q[:,:,1,i], Float32.(Kernel.gaussian(10)))
end

d_obs = G(q,m)

J = judiJacobian(F0,q0)
gradient_m = adjoint(J)*vec(G(q0,m0)-d_obs)

gs_inv = gradient(x -> misfit_objective(d_obs, q0, x, G),m0)

g1 = vec(gradient_m)
g2 = vec(gs_inv[1])

@test isapprox(norm(g1-g2) / norm(g1), 0f0; atol=1f-1)
@test isapprox(dot(g1,g2)/norm(g1)/norm(g2),1f0;rtol=1f-1)


text: const LIB_OPEN_SPIEL = "$(dirname(@__FILE__))/../../../build/julia/libspieljl"


text:
# Create example data (continent population in million people)
continents = ["Africa", "America", "Asia", "Europe", "Oceania"]
population_1960 = [285, 425, 1687, 606, 16]
population_2010 = [1044, 944, 4170, 735, 36]
population_matrix = [population_1960 population_2010]
# Plot with respect to 500 millions of people
barplot(continents, population_matrix, baseline=500)
legend("1960", "2010") # add legend


text:
# test to see that deleteHybridizationUpdate undoes all attributes
# prompted by Cecile finding cases when containRoot was not updated
# Claudia December 2015

PhyloNetworks.CHECKNET || error("need CHECKNET==true in PhyloNetworks to test snaq in test_correctLik.jl")

@testset "test: delete hybridization" begin
global seed, currT0, besttree, net, success,hybrid,flag,nocycle,flag2,flag3

#seed = 2738
seed = 56326

currT0 = readTopologyLevel1("(((((((1,2),3),4),5),(6,7)),(8,9)),10);");
## printEdges(currT0)
## printNodes(currT0)
## writeTopologyLevel1(currT0)
checkNet(currT0)
# warning: the random number generator has a local scope:
# with subsets of tests, the same seed would be re-used over and over.
Random.seed!(seed);
besttree = deepcopy(currT0);

# ===== first hybridization ==========================
# @testset "first hybridization" begin
success,hybrid,flag,nocycle,flag2,flag3 = addHybridizationUpdate!(besttree);
@test success
#printEdges(besttree)
#printNodes(besttree)
@test_logs writeTopologyLevel1(besttree)
# "(1:1.0,2:1.0,(3:1.0,(4:1.0,(5:1.0,((6:1.0,(7:0.4918119005933492,#H11:0.0::0.2227763931131359):0.5081880994066508):1.0,((8:1.0,(9:0.7497438380865815)#H11:0.2502561619134185::0.7772236068868641):1.0,10:1.0):1.0):1.0):1.0):1.0):1.0);"
net = deepcopy(besttree);
# test contain root
@test !net.edge[15].containRoot
@test !net.edge[19].containRoot
@test !net.edge[20].containRoot

# test inCycle
@test net.edge[12].inCycle == 11
@test net.edge[13].inCycle == 11
@test net.edge[16].inCycle == 11
@test net.edge[18].inCycle == 11
@test net.edge[19].inCycle == 11
@test net.edge[20].inCycle == 11

@test net.node[12].inCycle == 11
@test net.node[13].inCycle == 11
@test net.node[16].inCycle == 11
@test net.node[17].inCycle == 11
@test net.node[19].inCycle == 11
@test net.node[20].inCycle == 11

# test partition
@test length(net.partition) == 6
@test [n.number for n in net.partition[1].edges] == [15]
@test [n.number for n in net.partition[2].edges] == [11]
@test [n.number for n in net.partition[3].edges] == [10]
@test [n.number for n in net.partition[4].edges] == [9,7,5,3,1,2,4,6,8]
@test [n.number for n in net.partition[5].edges] == [17]
@test [n.number for n in net.partition[6].edges] == [14]
#end

# ===== second hybridization ==========================
#@testset "second hybridization" begin
success = false
success,hybrid,flag,nocycle,flag2,flag3 = addHybridizationUpdate!(besttree);
@test success
#printEdges(besttree)
#printNodes(besttree)
@test_logs writeTopologyLevel1(besttree,true)
# "(3,(4,(5,(((6,(7,#H11:::0.2227763931131359):0.5081880994066508):1.0,((8,(9)#H11:::0.7772236068868641):1.0,10):1.0):0.9561820134835957,#H13:0.0::0.28555825607592755):0.04381798651640434):1.0):1.0,((1,2):0.011459257289162528)#H13:0.9885407427108375::0.7144417439240724);"
net = deepcopy(besttree);

# test contain root
@test !net.edge[15].containRoot
@test !net.edge[19].containRoot
@test !net.edge[20].containRoot
@test !net.edge[3].containRoot
@test !net.edge[1].containRoot
@test !net.edge[2].containRoot
@test !net.edge[23].containRoot
@test !net.edge[22].containRoot

# test inCycle
@test net.edge[12].inCycle == 11
@test net.edge[13].inCycle == 11
@test net.edge[16].inCycle == 11
@test net.edge[18].inCycle == 11
@test net.edge[19].inCycle == 11
@test net.edge[20].inCycle == 11

@test net.node[12].inCycle == 11
@test net.node[13].inCycle == 11
@test net.node[16].inCycle == 11
@test net.node[17].inCycle == 11
@test net.node[19].inCycle == 11
@test net.node[20].inCycle == 11

@test net.edge[9].inCycle == 13
@test net.edge[7].inCycle == 13
@test net.edge[5].inCycle == 13
@test net.edge[22].inCycle == 13
@test net.edge[23].inCycle == 13

@test net.node[5].inCycle == 13
@test net.node[7].inCycle == 13
@test net.node[9].inCycle == 13
@test net.node[21].inCycle == 13
@test net.node[22].inCycle == 13

# test partition
@test length(net.partition) == 10
@test [n.number for n in net.partition[1].edges] == [15]
@test [n.number for n in net.partition[2].edges] == [11]
@test [n.number for n in net.partition[3].edges] == [10]
@test [n.number for n in net.partition[4].edges] == [17]
@test [n.number for n in net.partition[5].edges] == [14]
@test [n.number for n in net.partition[6].edges] == [3,1,2]
@test [n.number for n in net.partition[7].edges] == [21]
@test [n.number for n in net.partition[8].edges] == [8]
@test [n.number for n in net.partition[9].edges] == [6]
@test [n.number for n in net.partition[10].edges] == [4]
#end

## # ===== identify containRoot for net.node[21]
## net0=deepcopy(net);
## hybrid = net.node[21];
## hybrid.isBadDiamondI
## nocycle, edgesInCycle, nodesInCycle = identifyInCycle(net,hybrid);
## [n.number for n in edgesInCycle] == [23,9,7,5,22] || error("wrong identified cycle")
## [n.number for n in nodesInCycle] == [13,14,-4,-5,-6] || error("wrong identified cycle")
## edgesRoot = identifyContainRoot(net,hybrid);
## [n.number for n in edgesRoot] == [3,1,2] || error("wrong identified contain root")
## [e.containRoot for e in edgesRoot] == [false,false,false] || error("edges root should be false contain root")
## edges = hybridEdges(hybrid);
## [e.number for e in edges] == [22,23,3] || error("wrong identified edges")
## undoGammaz!(hybrid,net);
## othermaj = getOtherNode(edges[1],hybrid);
## othermaj.number == -6 || error("wrong othermaj")
## edgesmaj = hybridEdges(othermaj);
## [e.number for e in edgesmaj] == [22,5,4] || error("wrong identified edges")
## edgesmaj[3].containRoot || error("wrong edgesmaj[3] contain root")
## undoContainRoot!(edgesRoot);
## [e.containRoot for e in edgesRoot] == [true,true,true] || error("edges root should be false contain root")
## deleteHybrid!(hybrid,net,true, false)
## [e.containRoot for e in edgesRoot] == [true,true,true] || error("edges root should be false contain root")
## printEdges(net)

# ================= delete second hybridization =============================
@testset "delete 2nd hybridization" begin
deleteHybridizationUpdate!(net,net.node[21], false,false);
@test length(net.partition) == 6
# 15,11,10,[9,7,5,3,1,2,4,6,8],17,14
@test [n.number for n in net.partition[1].edges] == [15]
@test [n.number for n in net.partition[2].edges] == [11]
@test [n.number for n in net.partition[3].edges] == [10]
@test [n.number for n in net.partition[4].edges] == [17]
@test [n.number for n in net.partition[5].edges] == [14]
@test [n.number for n in net.partition[6].edges] == [3,1,2,8,6,4,9,7,5]
#printNodes(net)
#printEdges(net)
#printPartitions(net)
@test_logs writeTopologyLevel1(net)
# "(1:1.0,2:1.0,(3:1.0,(4:1.0,(5:1.0,((6:1.0,(7:0.4918119005933492,#H11:0.0::0.2227763931131359):0.5081880994066508):1.0,((8:1.0,(9:0.7497438380865815)#H11:0.2502561619134185::0.7772236068868641):1.0,10:1.0):1.0):1.0):1.0):1.0):1.0);"

# test contain root
@test !net.edge[15].containRoot
@test !net.edge[19].containRoot
@test !net.edge[20].containRoot
@test net.edge[1].containRoot
@test net.edge[2].containRoot
@test net.edge[3].containRoot

# test inCycle
@test net.edge[12].inCycle == 11
@test net.edge[13].inCycle == 11
@test net.edge[16].inCycle == 11
@test net.edge[18].inCycle == 11
@test net.edge[19].inCycle == 11
@test net.edge[20].inCycle == 11

@test net.node[12].inCycle == 11
@test net.node[13].inCycle == 11
@test net.node[16].inCycle == 11
@test net.node[17].inCycle == 11
@test net.node[19].inCycle == 11
@test net.node[20].inCycle == 11
end

# =============== delete first hybridization ===================
@testset "delete 1st hybridization" begin
deleteHybridizationUpdate!(net,net.node[19]);
checkNet(net)
@test length(net.partition) == 0
#printEdges(net)
end
end


text:
add_format(format"PLY_ASCII", "ply\nformat ascii 1.0", ".ply")
add_format(format"PLY_BINARY", "ply\nformat binary_little_endian 1.0", ".ply")

function type(t)
    if t == "float" return Float32
    else throw(ArgumentError("unexpected property type $t")) end
end

parseheader(io::IO, ::format"PLY_BINARY") = parseheader(io, format"PLY_ASCII"())
function parseheader(io::IO, ::format"PLY_ASCII")
    types, properties  = Type[], Symbol[]
    vertices, faces = 0, 0

    vertexmode = false
    line = readline(io)

    while !startswith(line, "end_header")
        if startswith(line, "e") # element *
            tmp = split(line)
            vertexmode = tmp[2] == "vertex"
            if vertexmode vertices = parse(Int, tmp[3])
            else faces = parse(Int, tmp[3]) end
        elseif startswith(line, "pr") && vertexmode # vertex property
            tmp = split(line)
            if !(tmp[3] in ["x", "y", "z"])
                push!(types, type(tmp[2]))
                push!(properties, Symbol(tmp[3]))
            end
        end
        line = readline(io)
    end

    if length(properties) == 0
        (vertices, faces, [], Tuple{}(), Point{3, Float32})
    else
        (vertices, faces, types, Tuple(properties),
        MetaPoint{3, Float32, Tuple(properties), Tuple{types...}})
    end
end

function load(file::File{format"PLY_ASCII"}; mesh=Mesh)
  io = file.io
  # read the header
  vertices, n_faces, types, properties, vertex = parseheader(io, format"PLY_ASCII"())

  points = Array{vertex}(undef, vertices)
  faces = Array{Triangle}(undef, n_faces)

  # read the data
  if length(types) == 0
    for i in 1:vertices
        points[i] = vertex(Tuple(parse.(Float32, split(readline(io)))))
    end
  else
    for i in 1:vertices
        tmp = split(readline(io))
        points[i] = vertex(
            Tuple(parse.(Float32, tmp[1:3])),
            NamedTuple{properties}(parse.(types, tmp[4:end])) # todo simplify meshes with metadata
        )
    end
  end

  for i = 1:n_faces
      line = split(readline(io))
      N = parse(Int, popfirst!(line))
      @assert N == 3
      faces[i] = Polygon{3}(parse.(UInt32, line))
  end
  mesh(points, faces)
end

function load(file::File{format"PLY_BINARY"}; mesh=Mesh)
    io = file.io
    # read the header
    vertices, n_faces, types, properties, vertex = parseheader(io, format"PLY_BINARY"())

    # read vertex data
    raw_vertices = read(io, sizeof(vertex) * vertices)

    # read face data
    faces = Array{Triangle}(undef, n_faces)
    for i = 1:n_faces
        N = read(io, UInt8)
        @assert N == 3
        faces[i] = read(io, Triangle)
    end

    mesh(unsafe_wrap(Vector{vertex}, unsafe_convert(Ptr{vertex}, raw_vertices), vertices), faces)
end


text:
using Whirl: runtests

runtests()

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
function make_global_state_map(inputs::SimulationInputs)
    dic = inputs.global_state_map
    if !isempty(dic)
        return dic
    end
    device_wrappers = get_dynamic_injectors(inputs)
    branches_wrappers = get_dynamic_branches(inputs)
    buses_diffs = get_voltage_buses_ix(inputs)
    n_buses = get_bus_count(inputs)
    for d in device_wrappers
        dic[PSY.get_name(d)] = OrderedDict{Symbol, Int}(get_global_index(d)...)
    end
    if !isempty(branches_wrappers)
        for br in branches_wrappers
            dic[PSY.get_name(br)] = OrderedDict{Symbol, Int}(get_global_index(br)...)
        end
    end
    if !isempty(buses_diffs)
        for ix in buses_diffs
            dic["V_$(ix)"] = OrderedDict{Symbol, Int}(:R => ix, :I => ix + n_buses)
        end
    end
    return dic
end

function get_state_from_ix(global_index::MAPPING_DICT, idx::Int)
    for (name, device_ix) in global_index
        if idx ∈ values(device_ix)
            state = [state for (state, number) in device_ix if number == idx]
            IS.@assert_op length(state) == 1
            return name, state[1]
        end
    end
    error("State with index $(idx) not found in the global index")
end


text:
# Load modules
using Plots, LaTeXStrings
pyplot()
using ThinFilmsTools

let
    λ = [632.8]
    n = [RIdb.air(λ) RIdb.silicon(λ)]
    p = LinRange(0,1,100)
    neff_looyenga = [RI.looyenga([p[i] 1.0-p[i]], n)[1] for i in eachindex(p)]
    neff_bruggeman = [RI.bruggeman(p[i], n)[1] for i in eachindex(p)]
    neff_maxwell = [RI.maxwell([p[i] 1.0-p[i]], n)[1] for i in eachindex(p)]
    neff_monecke = [RI.monecke(p[i], n)[1] for i in eachindex(p)]
    plot(p, real.(neff_looyenga), label="Looyenga", line=(:solid))
    plot!(p, real.(neff_bruggeman), label="Bruggeman", line=(:dash))
    plot!(p, real.(neff_maxwell), label="Maxwell", line=(:dashdot))
    plot!(p, real.(neff_monecke), label="Monecke", line=(3.0, :solid))
    xaxis!("Porosity, %")
    yaxis!("Effective dielectric permittivity")
    gui()
end

## lorentzlorenz EMA
# D. Schwarz et al. / Thin Solid Films 519 (2011) 2994-2997
let
    plt = plot();
    f = [0.0,0.25,0.5,0.75,1.0]
    ϵg = complex.(1.0:0.005:1.8)
    ϵₘ = 2.1*ones(length(ϵg))
    n = sqrt.([ϵg ϵₘ])
    for i in eachindex(f)
        neff = RI.lorentz_lorenz(f[i], n)
        plot!(plt, real.(ϵg), real.(neff.^2), label="f = $(f[i])")
    end
    xaxis!("Guest dielectric constant, ϵg")
    yaxis!("Lorentz-Lorenz effective permittivity")
    gui()
end


text:
function calcsort!(A,values::UnitRange{Int})
    store = zeros(Int, size(values))
    for i in 1:length(A)
        for j in 1:length(values)
            if A[i]==values[j]
                store[j]+=1
            end
        end
    end
    k=1
    for i in 1:length(values)
        while store[i]>0
            A[k]=values[i]
            k+=1
            store[i]-=1
        end
    end
    return A
end

function calcsort!(A,values::Vector{Int})
    store = zeros(Int, size(values))
    for i in 1:length(A)
        for j in 1:length(values)
            if A[i]==values[j]
                store[j]+=1
            end
        end
    end
    k=1
    for i in 1:length(values)
        while store[i]>0
            A[k]=values[i]
            k+=1
            store[i]-=1
        end
    end
    return A
end


text:
import Compat: Libc, @compat

function Base.connect(::Type{Postgres},
                      host::AbstractString="",
                      user::AbstractString="",
                      passwd::AbstractString="",
                      db::AbstractString="",
                      port::AbstractString="")
    conn = PQsetdbLogin(host, port, C_NULL, C_NULL, db, user, passwd)
    status = PQstatus(conn)

    if status != CONNECTION_OK
        errmsg = bytestring(PQerrorMessage(conn))
        PQfinish(conn)
        error(errmsg)
    end

    conn = PostgresDatabaseHandle(conn, status)
    finalizer(conn, DBI.disconnect)
    return conn
end

function Base.connect(::Type{Postgres},
                      host::AbstractString,
                      user::AbstractString,
                      passwd::AbstractString,
                      db::AbstractString,
                      port::Integer)
    Base.connect(Postgres, host, user, passwd, db, string(port))
end

# Note that for some reason, `do conn` notation
# doesn't work using this version of the function
function Base.connect(::Type{Postgres};
                      dsn::AbstractString="")
    conn = PQconnectdb(dsn)
    status = PQstatus(conn)
    if status != CONNECTION_OK
        errmsg = bytestring(PQerrorMessage(conn))
        PQfinish(conn)
        error(errmsg)
    end
    conn = PostgresDatabaseHandle(conn, status)
    finalizer(conn, DBI.disconnect)
    return conn
end

function DBI.disconnect(db::PostgresDatabaseHandle)
    if db.closed
        return
    else
        PQfinish(db.ptr)
        db.closed = true
        return
    end
end

function DBI.errcode(db::PostgresDatabaseHandle)
    return db.status = PQstatus(db.ptr)
end

function DBI.errstring(db::PostgresDatabaseHandle)
    return bytestring(PQerrorMessage(db.ptr))
end

function DBI.errcode(res::PostgresResultHandle)
    return PQresultStatus(res.ptr)
end

function DBI.errstring(res::PostgresResultHandle)
    return bytestring(PQresultErrorMessage(res.ptr))
end

DBI.errcode(stmt::PostgresStatementHandle) = DBI.errcode(stmt.result)
DBI.errstring(stmt::PostgresStatementHandle) = DBI.errstring(stmt.result)

function checkerrclear(result::Ptr{PGresult})
    status = PQresultStatus(result)

    try
        if status == PGRES_FATAL_ERROR
            statustext = bytestring(PQresStatus(status))
            errmsg = bytestring(PQresultErrorMessage(result))
            error("$statustext: $errmsg")
        end
    finally
        PQclear(result)
    end
end

escapeliteral(db::PostgresDatabaseHandle, value) = value
escapeliteral(db::PostgresDatabaseHandle, value::AbstractString) = escapeliteral(db, bytestring(value))

function escapeliteral(db::PostgresDatabaseHandle, value::Union{ASCIIString, UTF8String})
    strptr = PQescapeLiteral(db.ptr, value, sizeof(value))
    str = bytestring(strptr)
    PQfreemem(strptr)
    return str
end

function escapeidentifier(db::PostgresDatabaseHandle, value::Union{ASCIIString, UTF8String})
    strptr = PQescapeIdentifier(db.ptr, value, sizeof(value))
    str = bytestring(strptr)
    PQfreemem(strptr)
    return str
end

Base.run(db::PostgresDatabaseHandle, sql::AbstractString) = checkerrclear(PQexec(db.ptr, sql))

function checkcopyreturnval(db::PostgresDatabaseHandle, returnval::Int32)
    if returnval == -1
        errcode = bytestring(DBI.errcode(db))
        errmsg = bytestring(DBI.errmsg(db))
        error("Error $errcode: $errmsg")
    end
end

function copy_from(db::PostgresDatabaseHandle, table::AbstractString,
                   filename::AbstractString, format::AbstractString)
    f = open(filename)
    try
        Base.run(db, string("COPY ", escapeidentifier(db, table), " FROM STDIN ", format))
        for row in eachline(f)
            # send row to postgres
            checkcopyreturnval(db, PQputCopyData(db.ptr, row, length(row)))
        end
        checkcopyreturnval(db, PQputCopyEnd(db.ptr, C_NULL))
    finally
        close(f)
    end
    return checkerrclear(PQgetResult(db.ptr))
end

hashsql(sql::AbstractString) = bytestring(string("__", hash(sql), "__"))

function getparamtypes(result::Ptr{PGresult})
    nparams = PQnparams(result)
    return @compat [pgtype(OID{Int(PQparamtype(result, i-1))}) for i = 1:nparams]
end

LIBC = @windows ? "msvcrt.dll" : :libc
strlen(ptr::Ptr{UInt8}) = ccall((:strlen, LIBC), Csize_t, (Ptr{UInt8},), ptr)

function getparams!(ptrs::Vector{Ptr{UInt8}}, params, types, sizes, lengths::Vector{Int32}, nulls)
    fill!(nulls, false)
    for i = 1:length(ptrs)
        if params[i] === nothing || params[i] === NA || params[i] === Union{}
            nulls[i] = true
        else
            ptrs[i] = pgdata(types[i], ptrs[i], params[i])
            if sizes[i] < 0
                warn("Calling strlen--this should be factored out.")
                lengths[i] = strlen(ptrs[i]) + 1
            end
        end
    end
    return
end

function cleanupparams(ptrs::Vector{Ptr{UInt8}})
    for ptr in ptrs
        Libc.free(ptr)
    end
end

DBI.prepare(db::PostgresDatabaseHandle, sql::AbstractString) = PostgresStatementHandle(db, sql)

DBI.finish(stmt::PostgresStatementHandle) = nothing

function DBI.execute(stmt::PostgresStatementHandle)
    result = PQexec(stmt.db.ptr, stmt.stmt)
    return stmt.result = PostgresResultHandle(result)
end

function DBI.execute(stmt::PostgresStatementHandle, params::Vector)
    nparams = length(params)

    if nparams > 0 && isempty(stmt.paramtypes)
        paramtypes = [pgtype(typeof(p)) for p in params]
    else
        paramtypes = stmt.paramtypes
    end

    if nparams != length(paramtypes)
        error("Number of parameters in statement ($(length(stmt.paramtypes))) does not match number of " *
            "parameter values ($nparams).")
    end

    sizes = zeros(Int64, nparams)
    lengths = zeros(Cint, nparams)
    param_ptrs = fill(convert(Ptr{UInt8}, 0), nparams)
    nulls = falses(nparams)
    for i = 1:nparams
        if paramtypes[i] === nothing
            paramtypes[i] = pgtype(typeof(params[i]))
        end

        sizes[i] = sizeof(paramtypes[i])

        if sizes[i] > 0
            lengths[i] = sizes[i]
        end
    end
    formats = fill(PGF_TEXT, nparams)

    getparams!(param_ptrs, params, paramtypes, sizes, lengths, nulls)

    oids = Oid[convert(Oid, oid(p)) for p in paramtypes]

    result = PQexecParams(stmt.db.ptr, stmt.stmt, nparams,
        oids,
        [convert(Ptr{UInt8}, nulls[i] ? C_NULL : param_ptrs[i]) for i = 1:nparams],
        pointer(lengths), pointer(formats), PGF_TEXT)

    cleanupparams(param_ptrs)

    return stmt.result = PostgresResultHandle(result)
end

function executemany{T<:AbstractVector}(stmt::PostgresStatementHandle,
        params::Union{DataFrame,AbstractVector{T}})
    nparams = isa(params, DataFrame) ? ncol(params) : length(params[1])

    if nparams > 0 && isempty(stmt.paramtypes)
        if isa(params, DataFrame)
            paramtypes = collect(PostgresType, eltypes(params))
        else
            paramtypes = [pgtype(typeof(p)) for p in params[1]]
        end
    else
        paramtypes = stmt.paramtypes
    end

    if nparams != length(paramtypes)
        error("Number of parameters in statement ($(length(stmt.paramtypes))) does not match number of " *
            "parameter values ($nparams).")
    end

    sizes = zeros(Int64, nparams)
    lengths = zeros(Cint, nparams)
    param_ptrs = fill(convert(Ptr{UInt8}, 0), nparams)
    nulls = falses(nparams)
    for i = 1:nparams
        if paramtypes[i] === nothing
            paramtypes[i] = pgtype(typeof(params[1][i]))
        end

        if sizes[i] > 0
            lengths[i] = sizes[i]
        end
    end
    formats = fill(PGF_TEXT, nparams)

    result = C_NULL
    rowiter = isa(params, DataFrame) ? eachrow(params) : params
    for paramvec in rowiter
        getparams!(param_ptrs, paramvec, paramtypes, sizes, lengths, nulls)

        oids = Oid[convert(Oid, oid(p)) for p in paramtypes]

        result = PQexecParams(stmt.db.ptr, stmt.stmt, nparams,
            oids,
            [convert(Ptr{UInt8}, nulls[i] ? C_NULL : param_ptrs[i]) for i = 1:nparams],
            pointer(lengths), pointer(formats), PGF_TEXT)
    end

    cleanupparams(param_ptrs)

    return stmt.result = PostgresResultHandle(result)
end

function DBI.fetchrow(stmt::PostgresStatementHandle)
    error("DBI API not fully implemented")
end

# Assumes the row exists and has the structure described in PostgresResultHandle
function unsafe_fetchrow(result::PostgresResultHandle, rownum::Integer)
    return Any[PQgetisnull(result.ptr, rownum, i-1) == 1 ? Union{} :
               jldata(datatype, PQgetvalue(result.ptr, rownum, i-1))
               for (i, datatype) in enumerate(result.types)]
end

function unsafe_fetchcol_dataarray(result::PostgresResultHandle, colnum::Integer)
    return @data([PQgetisnull(result.ptr, i, colnum) == 1 ? NA :
            jldata(result.types[colnum+1], PQgetvalue(result.ptr, i, colnum))
            for i = 0:(PQntuples(result.ptr)-1)])
end

function DBI.fetchall(result::PostgresResultHandle)
    return Vector{Any}[row for row in result]
end

function DBI.fetchdf(result::PostgresResultHandle)
    df = DataFrame()
    for i = 0:(length(result.types)-1)
        df[symbol(bytestring(PQfname(result.ptr, i)))] = unsafe_fetchcol_dataarray(result, i)
    end

    return df
end

function Base.length(result::PostgresResultHandle)
    return PQntuples(result.ptr)
end

function Base.start(result::PostgresResultHandle)
    return 0
end

function Base.next(result::PostgresResultHandle, state)
    return (unsafe_fetchrow(result, state), state + 1)
end

function Base.done(result::PostgresResultHandle, state)
    return state >= result.nrows
end

# delegate statement iteration to result
Base.start(stmt::PostgresStatementHandle) = Base.start(stmt.result)
Base.next(stmt::PostgresStatementHandle, state) = Base.next(stmt.result, state)
Base.done(stmt::PostgresStatementHandle, state) = Base.done(stmt.result, state)

max_stars_repo_name: beffaCo/updatePostgresql

max_forks_repo_name: beffaCo/updatePostgresql


text:
using UnitTestDesign

entries = [
    (args = 5, vals = 4),
    (args = 5, vals = 8),
    (args = 10, vals = 4),
    (args = 10, vals = 8),
    (args = 40, vals = 4),
    (args = 40, vals = 8),
]
for e in entries
    vals = fill(1:e[:vals], e[:args])
    c = e[:vals]^e[:args]
    p = size(all_pairs(vals...))[1]
    t = size(all_triples(vals...))[1]
    ve = size(values_excursion(vals...))[1]
    pe = size(pairs_excursion(vals...))[1]
    # This is the power of 10 for the arguments.
    # val^arg = 10^(log10(val^arg)) = 10^(arg * log10(val))
    pow10 = e[:args] * log10(e[:vals])
    println("$(e[:vals]) $(e[:args]) $c $(pow10) $p $t $ve $pe")
end

params = fill(1:4, 40)
ts = all_pairs(params...;
    wayness = Dict(3 => [collect(3:6)])
    )


text:
export pivotedqr

function pivotedqr(ψ::AbstractMatrix{T}, u::AbstractVector{T}; invert::Bool=true, verbose::Bool = true, ϵrel::Float64 = 1e-1, maxterms::Int64=typemax(Int64)) where {T}

    m, n = size(ψ)

    # Compute a pivoted QR decomposition of ψ, use the pivot for the sensor placement

    F = qr(ψ, Val(true))

    residue = copy(u)
    idxset = Int64[]
    dict = collect(1:n)
    ϵu = norm(u)
    ϵhist = [ϵu]
    ϵrel *= ϵu

    @inbounds for k=1:n
        # Update set of selected basis
        push!(idxset, F.p[k])

        # Update candidate dictionary
        filter!(x-> x != F.p[k], dict)

        # Compute residual
        c = view(ψ,:,idxset)\u
        residue = u - view(ψ,:,idxset)*c

        # Calculate stopping critera
        ϵ = norm(residue)
        if verbose == true
            push!(ϵhist, copy(ϵ))
        end
        if ϵ < ϵrel || k == maxterms
            break
        end

    end

    # Solve the system with the set of indices
    if verbose == true
        if invert == true
            c = zeros(length(idxset))
            c .= view(ψ,:, idxset)\u
            return idxset, c, ϵhist
        else
            return idxset, ϵhist
        end
    else
        if invert == true
            c = zeros(length(idxset))
            c .= view(ψ,:, idxset)\u
            return idxset, c
        else
            return idxset
        end
    end
end

max_forks_repo_name: mleprovost/QROMP