text:
using LAR
include("src/inters.jl")

# Compute the codomain of functions (next,prev) for the edge indices
function nextPrev(edges::Array{Int64,1})
	theNexts,thePrevs = Array{Int64,1}(), Array{Int64,1}()
	n = length(edges)
	for k in 1:n
		push!(theNexts, k==n ? edges[1] : edges[k+1] )
		push!(thePrevs, k==1 ? edges[n] : edges[k-1] )
	end
	theNexts,thePrevs
end

# Circular ordering of edges around vertices
function edgeSlopeOrdering(VE,V,EV)
	VE_sorted = Array{Array{Int64,1},1}()
	forwardBackward = Tuple{Array{Int64,1},Array{Int64,1},Array{Int64,1}}[]
	for (v,ve) in enumerate(VE)
		ve_angle = []
		if ve != []
			for edge in ve
				w = pop!(setdiff(Set(EV[edge]), Set(v)))
				(x,y) = V[:,w] - V[:,v]
				angle = atan2(y,x)
				push!(ve_angle, 180*angle/pi)
			end
		end
		pairs = sort(collect(zip(ve_angle,ve)))
		sortedEdges = [pair[2] for pair in pairs]
		push!(VE_sorted, sortedEdges)
		theNexts,thePrevs = nextPrev(sortedEdges)
		push!(forwardBackward, (sortedEdges,theNexts,thePrevs))

	end
	forwardBackward
end

function theSign(EV,edge,v)
	w = pop!(setdiff(Set(EV[edge]), Set(v)))
	v < w ? 1 : -1
end

# Return the sparse matrix codifying (next,prev) edges for each pair [v,e]
function cycleBasis(V,EV)
	I = Int64[];
	J = Int64[];
	Vnext, Vprev = Int64[], Int64[]
	for (e,(i,j)) in enumerate(EV)
		push!(I,i); push!(I,j);
		push!(J,e); push!(J,e);
		push!(Vnext,0); push!(Vnext,0);
		push!(Vprev,0); push!(Vprev,0);
	end
	Snext = sparse(J,I,Vnext)
	Sprev = sparse(J,I,Vprev)
	rows = rowvals(Snext)
	VE = Array{Array{Int64,1},1}()
	for k=1:size(Snext,2)
		# extract column
		col = [rows[h] for h in nzrange(Snext,k)]
		push!(VE, col)
	end
	forwardBackward = edgeSlopeOrdering(VE,V,EV)
	for (v,(edges,nexts,prevs)) in enumerate(forwardBackward)
		for (h,e) in enumerate(edges)
			esign = theSign(EV,edges[h],v)
			Snext[edges[h],v] = nexts[h]*esign
			Sprev[edges[h],v] = prevs[h]*esign
		end
	end
	return Snext,Sprev
end

# Return the sparse matrix of the signed ∂–1 operator
function boundary_1(EV)
	we = cellComplex(EV)';
	rows = rowvals(we);
	for e in 1:size(we,2)
		i,j = collect(nzrange(we,e))
		h,k = rows[i], rows[j]
		if h < k
			we[h,e] = -1
		else
			we[k,e] = -1
		end
	end;
	we
end

datafile = readcsv("test/svg/test1.lines");
V = reshape(datafile',(size(datafile',1)÷2,size(datafile',2)*2));
len = length(datafile);
EV = collect(reshape(1:(len÷2), 2,(len÷4)));
lineArray = (V,EV);
W,EW = lines2lar(lineArray);
lar = Lar();
lar.EV = [EW[:,k] for k=1:size(EW,2)];
W,newcells,oldindex = larvalidate(W,lar,10^2);
EW = newcells.EV;
viewLarIndices(W,EW,0.75)

using Iterators

function colProd(op,chain1)
	chain0 = op * chain1
	chain0_cells = []
	I,J,V = findnz(chain0)
	@itr for (i,v) in zip(I,V)
    	push!( chain0_cells, i*v )
    end
	chain0_cells
end

function rowProd(chain0,op)
	chain1 = chain0' * op
	chain1_cells = []
	I,J,V = findnz(chain1)
	@itr for (j,v) in zip(J,V)
    	push!( chain1_cells, j*v )
    end
	chain1_cells
end

function extract_2cell(seed, Snext,Sprev, useCounts, B_1, EW)
	chain1 = map(Int64, spzeros(size(B_1,2),1));
	if useCounts[seed] == 0
		chain1[seed,1] = +1
	elseif useCounts[seed] == 1
		chain1[seed,1] = -1
	end
	chain1_cells = [chain1[seed,1] * seed]
	chain0_cells = colProd(B_1, chain1)
	useCounts[seed] += 1

	while chain0_cells != []
		chain1_cells = rowProd(chain0, B_1)

		next1_cells = [ S[abs(cell1),abs(cell0)] for cell0 in chain0_cells,
			cell1 in chain1_cells if Set(abs(cell0)) < Set(EW[cell1])]

		direction_test = all([ useCounts[cell]<2 for cell in next1_cells ])
		if direction_test==true
			for cell in Set(next1_cells)
				useCounts[cell] += 1
				chain1[cell,1] = 1
				push!(chain1_cells, cell)
			end
		else
			return []
		end
		chain0 = B_1 * chain1
		cells0,val = findnz(chain0[:,1])
		chain0_cells = [cells0[k] for (k,v) in enumerate(val) if abs(val[k])==1]
	end
	chain1_cells
end

Snext,Sprev = cycleBasis(W,EW)
B_1 = boundary_1(EW)
S = Snext
useCounts = [0 for k in 1:length(EW)]
extract_2cell(1, S, useCounts, B_1, EW)


text:
module PiMC

using LinearAlgebra
using Base.Threads

include("estimators.jl")
export serial_pi_estimation, pi_estimation, atomic_pi_estimation

end # module

ext: jl


text:
#=
#define FMPR_RND_DOWN  0    RoundingMode{:Down}()
#define FMPR_RND_UP    1
#define FMPR_RND_FLOOR 2
#define FMPR_RND_CEIL  3
#define FMPR_RND_NEAR  4
=#

import Base.Rounding: rounding_raw, setrounding_raw, rounding, setrounding;

const ARB_ROUNDING_MODE = Ref{Cint}(0)

to_arb(::RoundingMode{:Nearest}) = Cint(0)
to_arb(::RoundingMode{:ToZero}) = Cint(1)
to_arb(::RoundingMode{:Up}) = Cint(2)
to_arb(::RoundingMode{:Down}) = Cint(3)
to_arb(::RoundingMode{:FromZero}) = Cint(4)

function from_arb(c::Integer)
    if c == 0
        return RoundNearest
    elseif c == 1
        return RoundToZero
    elseif c == 2
        return RoundUp
    elseif c == 3
        return RoundDown
    elseif c == 4
        return RoundFromZero
    else
        throw(ArgumentError("invalid MPFR rounding mode code: $c"))
    end
    RoundingMode(c)
end

rounding_raw(::Type{ArbFloat}) = ARB_ROUNDING_MODE[]
setrounding_raw(::Type{ArbFloat},i::Integer) = ARB_ROUNDING_MODE[] = (Cint)(i)

rounding(::Type{ArbFloat}) = from_arb(rounding_raw(ArbFloat))
setrounding(::Type{ArbFloat},r::RoundingMode) = setrounding_raw(ArbFloat,to_arb(r))

rounding_raw(::Type{ArbFloat{P}}) where {P} = ARB_ROUNDING_MODE[]
setrounding_raw(::Type{ArbFloat{P}},i::Integer) where {P} = ARB_ROUNDING_MODE[] = (Cint)(i)

rounding(::Type{ArbFloat{P}}) where {P} = from_arb(rounding_raw(ArbFloat))
setrounding(::Type{ArbFloat{P}},r::RoundingMode) where {P} = setrounding_raw(ArbFloat,to_arb(r))


text:
export AbstractModel
"""
Abstract type whose subtypes are the models from which initial-value problems can be built,
through `Problem( model :: AbstractModel, initial :: InitialData, param :: NamedTuple )`
"""
abstract type AbstractModel end
show(io::IO, m::AbstractModel) =
    try print(io,m.info)
    catch
        print(io,"Model: $(m.label)")
    end

function Base.:(==)(m1::AbstractModel, m2::AbstractModel)
    x=rand(2)
    init=Random(x)
    U1=m1.mapto(init);U2=m2.mapto(init);
    data1=m1.mapfro(U1);data2=m2.mapfro(U2);
    m1.f!(U1);m2.f!(U2);

    U1 == U2 && data1 == data2 && m1.label == m2.label
end


text:
export CuStream

type CuStream
  ptr::CUstream

  function CuStream(ptr)
    s = new(ptr)
    finalizer(s, cuStreamDestroy)
    s
  end
end

function CuStream()
  p = CUstream[0]
  cuStreamCreate(p, 0)
  CuStream(p[1])
end

Base.unsafe_convert(::Type{CUstream}, s::CuStream) = s.ptr

function test_stream(f::Function, num::Int)
  while length(custreams) < num
    push!(custreams, CuStream())
  end
  results = []
  @sync begin
    for i = 1:num
      @async begin
        while true
          s = custreams[i]
          results[i] = f()
        end
      end
    end
  end
end


text:
function vcatsort(v1, v2)
  n2 = length(v2)
  n2 == 0 && return v1
  n1 = length(v1)
  n1 == 0 && return v2

  n = n1 + n2
  res = similar(v1, n)
  c1, c2 = 1, 1
  @inbounds for i = 1:n
    if c2 == n2 + 1
      res[i] = v1[c1]
      c1 += 1
    elseif c1 == n1 + 1
      res[i] = v2[c2]
      c2 += 1
    else
      if v1[c1] < v2[c2]
        res[i] = v1[c1]
        c1 += 1
      else
        res[i] = v2[c2]
        c2 += 1
      end
    end
  end

  return res
end

function sparse_dropzeros(rows, cols, vals::Vector, nrows, ncols)
  M = sparse(rows, cols, vals, ncols, nrows)
  dropzeros!(M)
  return M
end

function get_QM_data(QM::AbstractQuadraticModel{T, S}, uplo::Symbol) where {T <: Real, S}
  # constructs A and Q transposed so we can create K upper triangular.
  # Q is lower triangular in QuadraticModels.jl
  if uplo == :U # A is Aᵀ of QuadraticModel QM
    A = sparse_dropzeros(QM.data.Acols, QM.data.Arows, QM.data.Avals, QM.meta.ncon, QM.meta.nvar)
    Q = sparse_dropzeros(QM.data.Hcols, QM.data.Hrows, QM.data.Hvals, QM.meta.nvar, QM.meta.nvar)
  else
    A = sparse_dropzeros(QM.data.Arows, QM.data.Acols, QM.data.Avals, QM.meta.nvar, QM.meta.ncon)
    Q = sparse_dropzeros(QM.data.Hrows, QM.data.Hcols, QM.data.Hvals, QM.meta.nvar, QM.meta.nvar)
  end

  id = QM_IntData(
    vcatsort(QM.meta.ilow, QM.meta.irng),
    vcatsort(QM.meta.iupp, QM.meta.irng),
    QM.meta.irng,
    QM.meta.ifree,
    QM.meta.ifix,
    QM.meta.ncon,
    QM.meta.nvar,
    0,
    0,
  )
  id.nlow, id.nupp = length(id.ilow), length(id.iupp) # number of finite constraints
  @assert QM.meta.lcon == QM.meta.ucon # equality constraint (Ax=b)
  @assert length(QM.meta.lvar) == length(QM.meta.uvar) == QM.meta.nvar
  fd = QM_FloatData(Q, A, QM.meta.lcon, QM.data.c, QM.data.c0, QM.meta.lvar, QM.meta.uvar, uplo)
  return fd, id
end

function allocate_workspace(
  QM::AbstractQuadraticModel,
  iconf::InputConfig,
  itol::InputTol,
  start_time,
  T0::DataType,
)
  sc = StopCrit(false, false, false, itol.max_iter, itol.max_time, start_time, 0.0)

  # save inital IntData to compute multipliers at the end of the algorithm
  idi = IntDataInit(
    QM.meta.nvar,
    QM.meta.ncon,
    QM.meta.ilow,
    QM.meta.iupp,
    QM.meta.irng,
    QM.meta.ifix,
    QM.meta.jlow,
    QM.meta.jupp,
    QM.meta.jrng,
    QM.meta.jfix,
  )

  if !QM.meta.minimize
    QM.data.Hvals .= .-QM.data.Hvals
    QM.data.c .= .-QM.data.c
    QM.data.c0 = -QM.data.c0
  end

  QM = SlackModel(QM)
  if QM.meta.ncon == length(QM.meta.jfix) && !iconf.presolve && iconf.scaling
    QM = deepcopy(QM) # if not modified by SlackModel and presolve
  end

  uplo = iconf.sp.uplo
  fd_T0, id = get_QM_data(QM, uplo) # apply presolve at the same time

  T = T0 # T0 is the data type, in mode :multi T will gradually increase to T0
  ϵ = Tolerances(
    T(itol.ϵ_pdd),
    T(itol.ϵ_rb),
    T(itol.ϵ_rc),
    one(T),
    one(T),
    T(itol.ϵ_μ),
    T(itol.ϵ_Δx),
    iconf.normalize_rtol,
  )

  S0 = typeof(fd_T0.c)
  if iconf.scaling
    if uplo == :U
      m, n = id.nvar, id.ncon
    else
      m, n = id.ncon, id.nvar
    end
    sd = ScaleData{T0, S0}(
      fill!(S0(undef, id.ncon), one(T0)),
      fill!(S0(undef, id.nvar), one(T0)),
      fill!(S0(undef, id.nvar), one(T0)),
      S0(undef, m),
      S0(undef, n),
    )
  else
    empty_v = S0(undef, 0)
    sd = ScaleData{T0, S0}(empty_v, empty_v, empty_v, empty_v, empty_v)
  end

  # allocate data for iterations
  if iconf.mode == :multi
    T = Float32
  end
  S = S0.name.wrapper{T, 1}

  res = init_residuals(S(undef, id.ncon), S(undef, id.nvar), zero(T), zero(T), iconf, id)

  itd = IterData(
    S(undef, id.nvar + id.ncon), # Δxy
    S(undef, id.nlow), # Δs_l
    S(undef, id.nupp), # Δs_u
    S(undef, id.nlow), # x_m_lvar
    S(undef, id.nupp), # uvar_m_x
    S(undef, id.nvar), # init Qx
    S(undef, id.nvar), # init ATy
    S(undef, id.ncon), # Ax
    zero(T), #xTQx
    zero(T), #cTx
    zero(T), #pri_obj
    zero(T), #dual_obj
    zero(T), #μ
    zero(T),#pdd
    zeros(T, 6), #l_pdd
    one(T), #mean_pdd
    nnz(fd_T0.Q) > 0,
    QM.meta.minimize,
  )

  dda_type = Symbol(:DescentDirectionAllocs, iconf.solve_method)
  dda = eval(dda_type)(id, S)

  cnts = Counters(zero(Int), zero(Int), 0, 0, iconf.kc, iconf.max_ref, zero(Int), iconf.w)

  pt = Point(S(undef, id.nvar), S(undef, id.ncon), S(undef, id.nlow), S(undef, id.nupp))

  spd = StartingPointData{T, S}(S(undef, id.nvar), S(undef, id.nlow), S(undef, id.nupp))

  #####
  return sc, idi, fd_T0, id, ϵ, res, itd, dda, pt, sd, spd, cnts, T
end

function allocate_extra_workspace_32(itol::InputTol, iconf::InputConfig, fd_T0::QM_FloatData)
  T = Float32
  ϵ32 = Tolerances(
    T(itol.ϵ_pdd32),
    T(itol.ϵ_rb32),
    T(itol.ϵ_rc32),
    one(T),
    one(T),
    T(itol.ϵ_μ),
    T(itol.ϵ_Δx),
    iconf.normalize_rtol,
  )
  fd32 = convert_FloatData(T, fd_T0)
  return fd32, ϵ32, T
end

function allocate_extra_workspace_64(itol::InputTol, iconf::InputConfig, fd_T0::QM_FloatData)
  T = Float64
  fd64 = convert_FloatData(T, fd_T0)
  ϵ64 = Tolerances(
    T(itol.ϵ_pdd64),
    T(itol.ϵ_rb64),
    T(itol.ϵ_rc64),
    one(T),
    one(T),
    T(itol.ϵ_μ),
    T(itol.ϵ_Δx),
    iconf.normalize_rtol,
  )
  T = Float32
  return fd64, ϵ64, T
end

function initialize!(
  fd::QM_FloatData{T},
  id::QM_IntData,
  res::AbstractResiduals{T},
  itd::IterData{T},
  dda::DescentDirectionAllocs{T},
  pt::Point{T},
  spd::StartingPointData{T},
  ϵ::Tolerances{T},
  sc::StopCrit{Tc},
  iconf::InputConfig{Tconf},
  cnts::Counters,
  T0::DataType,
) where {T <: Real, Tc <: Real, Tconf <: Real}
  pad = PreallocatedData(iconf.sp, fd, id, itd, pt, iconf)

  # init system
  # solve [-Q-D    A' ] [x] = [0]  to initialize (x, y, s_l, s_u)
  #       [  A     0  ] [y] = [b]
  itd.Δxy[1:(id.nvar)] .= 0
  itd.Δxy[(id.nvar + 1):end] = fd.b
  if typeof(pad) <: PreallocatedDataNewton
    itd.Δs_l .= zero(T)
    itd.Δs_u .= zero(T)
    pt.s_l .= one(T)
    pt.s_u .= one(T)
    itd.x_m_lvar .= one(T)
    itd.uvar_m_x .= one(T)
  end
  out = solver!(itd.Δxy, pad, dda, pt, itd, fd, id, res, cnts, T0, :init)
  pt.x .= itd.Δxy[1:(id.nvar)]
  pt.y .= itd.Δxy[(id.nvar + 1):(id.nvar + id.ncon)]

  starting_points!(pt, fd, id, itd, spd)

  # stopping criterion
  #     rcNorm, rbNorm = norm(rc), norm(rb)
  #     optimal = pdd < ϵ_pdd && rbNorm < ϵ_rb && rcNorm < ϵ_rc
  res.rb .= itd.Ax .- fd.b
  res.rc .= itd.ATy .- itd.Qx .- fd.c
  res.rc[id.ilow] .+= pt.s_l
  res.rc[id.iupp] .-= pt.s_u
  res.rcNorm, res.rbNorm = norm(res.rc, Inf), norm(res.rb, Inf)
  typeof(res) <: ResidualsHistory && push_history_residuals!(res, itd, pad, id)
  set_tol_residuals!(ϵ, res.rbNorm, res.rcNorm)

  sc.optimal = itd.pdd < ϵ.pdd && res.rbNorm < ϵ.tol_rb && res.rcNorm < ϵ.tol_rc
  sc.small_μ = itd.μ < ϵ.μ

  return pad
end

function set_tol_residuals!(ϵ::Tolerances{T}, rbNorm::T, rcNorm::T) where {T <: Real}
  if ϵ.normalize_rtol == true
    ϵ.tol_rb, ϵ.tol_rc = ϵ.rb * (one(T) + rbNorm), ϵ.rc * (one(T) + rcNorm)
  else
    ϵ.tol_rb, ϵ.tol_rc = ϵ.rb, ϵ.rc
  end
end

############ tools for sparse matrices ##############

function get_diag_sparseCSC(M_colptr, n; tri = :U)
  # get diagonal index of M.nzval
  # we assume all columns of M are non empty, and M triangular (:L or :U)
  @assert tri == :U || tri == :L
  if tri == :U
    diagind = M_colptr[2:end] .- one(Int)
  else
    diagind = M_colptr[1:(end - 1)]
  end
  return diagind
end

function get_diag_Q(Q_colptr, Q_rowval, Q_nzval::Vector{T}, n) where {T <: Real}
  diagval = spzeros(T, n)
  @inbounds @simd for j = 1:n
    for k = Q_colptr[j]:(Q_colptr[j + 1] - 1)
      if j == Q_rowval[k]
        diagval[j] = Q_nzval[k]
      end
    end
  end
  return diagval
end

function get_diag_Q_dense(Q::SparseMatrixCSC{T, Int}) where {T <: Real}
  n = size(Q, 1)
  diagval = zeros(T, n)
  fill_diag_Q_dense!(Q.colptr, Q.rowval, Q.nzval, diagval, n)
  return diagval
end

function fill_diag_Q_dense!(
  Q_colptr,
  Q_rowval,
  Q_nzval::Vector{T},
  diagval::Vector{T},
  n,
) where {T <: Real}
  for j = 1:n
    k = Q_colptr[j + 1] - 1
    if k > 0
      i = Q_rowval[k]
      if j == i
        diagval[j] = Q_nzval[k]
      end
    end
  end
end


text:
module RandomForestBehaviors

using Reexport

include("MvDecisionTrees.jl")
using .MvDecisionTrees

include("GindeleRandomForestBehaviors.jl")
@reexport using .GindeleRandomForestBehaviors

include("DynamicForestBehaviors.jl")
@reexport using .DynamicForestBehaviors

include("GaussianMixtureRegressors.jl")
@reexport using .GaussianMixtureRegressors

end # module


text:
using Statistics
using DataFrames
using CSV
using JLD

# %%

data_dir = "/home/user/Downloads/raw_data/csv";
data_files = readdir(data_dir);
output_file = "/home/user/thesis/data/scaling/Network_stochastic_edge_acceptance.csv";

# %%

dfs_raw = [];
for data_file in data_files
	df_raw = CSV.read(joinpath(data_dir, data_file)) |> DataFrame
	push!(dfs_raw, df_raw)
end

df = dfs_raw[1];
for df_raw in dfs_raw[2:end]
	append!(df, df_raw)
end

# %%

df[!, :Δr]  = df[:, :r_1] .- df[:, :r_0];
df[!, :Δm]  = df[:, :m_1] .- df[:, :m_0];
df[!, :t_0] = Int.(df[:, :r_0] .* df[:, :N]);
df[!, :t_1] = Int.(df[:, :r_1] .* df[:, :N]);

# t_data = Dict();
# for row in 1:size(df)[1]
# 	t_data[(df[row, :N], df[row, :seed])] = (df[row, :t_0], df[row, :t_1])
# end
# t_data
# save("/home/user/thesis/code/scaling/t_data.jld", "t_data", t_data)

df[!, :N] = Int.(df[:, :N]);
df[!, :t] = Int.(df[:, :t]);
df[!, :seed] = Int.(df[:, :seed]);
sort!(df, [:N, :seed]);

# %%

df |> CSV.write(output_file)

max_stars_repo_name: cameronperot/explosive-percolation

max_issues_repo_name: cameronperot/explosive-percolation

max_forks_repo_name: cameronperot/explosive-percolation


text:
# ---
# title: 787. Cheapest Flights Within K Stops
# id: problem787
# author: Indigo
# date: 2021-01-28
# difficulty: Medium
# categories: Dynamic Programming, Heap, Breadth-first Search
# link: <https://leetcode.com/problems/cheapest-flights-within-k-stops/description/>
# hidden: true
# ---
#
# There are `n` cities connected by `m` flights. Each flight starts from city
# `u` and arrives at `v` with a price `w`.
#
# Now given all the cities and flights, together with starting city `src` and
# the destination `dst`, your task is to find the cheapest price from `src` to
# `dst` with up to `k` stops. If there is no such route, output `-1`.
#
#
#
#     **Example 1:**
#     Input:
#     n = 3, edges = [[0,1,100],[1,2,100],[0,2,500]]
#     src = 0, dst = 2, k = 1
#     Output: 200
#     Explanation:
#     The graph looks like this:
#     ![](https://s3-lc-upload.s3.amazonaws.com/uploads/2018/02/16/995.png)
#
#     The cheapest price from city 0 to city 2 with at most 1 stop costs 200, as marked red in the picture.
#
#
#     **Example 2:**
#     Input:
#     n = 3, edges = [[0,1,100],[1,2,100],[0,2,500]]
#     src = 0, dst = 2, k = 0
#     Output: 500
#     Explanation:
#     The graph looks like this:
#     ![](https://s3-lc-upload.s3.amazonaws.com/uploads/2018/02/16/995.png)
#
#     The cheapest price from city 0 to city 2 with at most 0 stop costs 500, as marked blue in the picture.
#
#
#
#
# **Constraints:**
#
#   * The number of nodes `n` will be in range `[1, 100]`, with nodes labeled from `0` to `n`` - 1`.
#   * The size of `flights` will be in range `[0, n * (n - 1) / 2]`.
#   * The format of each flight will be `(src, ``dst``, price)`.
#   * The price of each flight will be in the range `[1, 10000]`.
#   * `k` is in the range of `[0, n - 1]`.
#   * There will not be any duplicated flights or self cycles.
#
#
## @lc code=start
using LeetCode

function find_cheapest_price(n::Int, flights::Vector{Vector{Int}}, src::Int, dst::Int, K::Int)
    graph = [Dict{Int, Int}() for i in 1:n]
    for flight in flights
        graph[flight[1]][flight[2]] = flight[3]
    end
    cost = fill(typemax(Int) >> 1, n)
    cost[src] = 0
    changed = [src]
    for i in 1:K+1
        tmp = Int[]
        for new_s in changed
            for (s_to_t, st_cost) in graph[new_s]
                if cost[s_to_t] > st_cost + cost[new_s]
                    cost[s_to_t] = st_cost + cost[new_s]
                    push!(tmp, s_to_t)
                end
            end
        end
        changed = tmp
    end
    return (cost[dst] == typemax(Int) >> 1) ? -1 : cost[dst]
end
## @lc code=end


text:
using JLD
using Plots
using LaTeXStrings
using Colors
#reset past font scaling
Plots.scalefontsizes()
#scale fonts
Plots.scalefontsizes(1.75)

orange = colorant"rgb(85%,55%,13%)"

#Extracting stored varibles
ld = load("./out/jld/TestCGConvergence.jld")
for key in keys(ld)
  symb = Symbol(key)
  data = ld[key]
  @eval $(symb) = $(data)
end

outPlot = plot(xlabel=L"\mathrm{Iterations}",
            ylabel=L"\log_{10}(\max_{\rho, \sigma, \nu} \|\mu - \mu_{\mathrm{app}}\| / \|\mu\|)")

plot!(outPlot, log10.(maxError), legend=false, linestyle=:dash, markershape=:dtriangle, markersize=15, linewidth=5, color=orange)

savefig(outPlot, "./out/plots/TestCGConvergence.pdf")


text:
@testset "utility" begin
    p1 = rand(SVector{3,Float64})
    p2 = rand(SVector{3,Float64})
    p3 = rand(SVector{3,Float64})
    p4 = rand(SVector{3,Float64})

    m3 = asMat(p1, p2, p3)
    m4 = asMat(p1, p2, p3, p4)
    m3_pad = asMatOnePad(p1, p2, p3)
    m4_pad = asMatOnePad(p1, p2, p3, p4)
    @test all(m3 .== hcat(p1, p2, p3))
    @test all(m4 .== hcat(p1, p2, p3, p4))
    @test all(m3_pad .== hcat(onePad(p1), onePad(p2), onePad(p3)))
    @test all(m4_pad .== hcat(onePad(p1), onePad(p2), onePad(p3), onePad(p4)))

    m3 = asMat( (p1, p2, p3) )
    m4 = asMat( (p1, p2, p3, p4) )
    m3_pad = asMatOnePad( (p1, p2, p3) )
    m4_pad = asMatOnePad( (p1, p2, p3, p4) )
    @test all(m3 .== hcat(p1, p2, p3))
    @test all(m4 .== hcat(p1, p2, p3, p4))
    @test all(m3_pad .== hcat(onePad(p1), onePad(p2), onePad(p3)))
    @test all(m4_pad .== hcat(onePad(p1), onePad(p2), onePad(p3), onePad(p4)))

    m3 = asMat( SVector(p1, p2, p3) )
    m4 = asMat( SVector(p1, p2, p3, p4) )
    m3_pad = asMatOnePad( SVector(p1, p2, p3) )
    m4_pad = asMatOnePad( SVector(p1, p2, p3, p4) )
    @test all(m3 .== hcat(p1, p2, p3))
    @test all(m4 .== hcat(p1, p2, p3, p4))
    @test all(m3_pad .== hcat(onePad(p1), onePad(p2), onePad(p3)))
    @test all(m4_pad .== hcat(onePad(p1), onePad(p2), onePad(p3), onePad(p4)))
end


text:
using Base.Test, Distributions, DataFrames
using UnobservableQueue
const UQ = UnobservableQueue

@testset "Unobservable Queue" begin
    @testset "Queue Simulation" begin
        # M/M/2 queue
        c = 2 # number of servers
        μ = 1 / c # service rate (mean service time is 1/μ) (single server)
        λ = 0.9 # arrival rate (mean interarrival time is 1/λ)
        adist = Exponential(1/λ) # arrival distribution
        sdist = Exponential(1/μ) # service distribution (single server)
        ncust = 10_000 # total number of customers generated
        timelimit = 1_000 # time limit for simulation
        seed = 8710 # rng seed
        s1 = UQ.Sim(ncust, timelimit, seed)
        q1 = UQ.Queue(c, adist, sdist)
        df1 = UQ.runsim(s1, q1)
        @test mean(df1[:dtime] - df1[:stime]) ≈ (1 / μ) atol = 0.1
        @test df1[:atime] |> diff |> mean ≈ (1 / λ) atol = 0.1
        @test all(diff(df1[:atime]) .> 0)
        @test all(diff(sort(df1, :dorder)[:dtime]) .> 0)

        # G/G/15 queue
        c = 15 # number of servers
        μa = (1 / 2) * 0.99
        μs = 1 / 2
        σ = 0.9
        adist = LogNormal(exp(1/μa), σ) # arrival distribution
        sdist = LogNormal(exp(1/μs), σ) # service distribution (single server)
        ncust = 10_000 # total number of customers generated
        timelimit = 10_000 # time limit for simulation
        seed = 3001 # rng seed
        s2 = UQ.Sim(ncust, timelimit, seed)
        q2 = UQ.Queue(c, adist, sdist)
        df2 = UQ.runsim(s2, q2)
        @test mean(df2[:dtime] - df2[:stime]) ≈ mean(rand(sdist, 10_000)) atol = 100
        @test df2[:atime] |> diff |> mean ≈ mean(rand(adist, 10_000)) atol = 100
        @test all(diff(df2[:atime]) .> 0)
        @test all(diff(sort(df2, :dorder)[:dtime]) .> 0)
    end

    @testset "Server Estimation" begin
        # Introduce observation error into true data
        c = 2 # number of servers
        μ = 1 / c # service rate (mean service time is 1/μ) (single server)
        λ = 0.9 # arrival rate (mean interarrival time is 1/λ)
        adist = Exponential(1/λ) # arrival distribution
        sdist = Exponential(1/μ) # service distribution (single server)
        ncust = 10_000 # total number of customers generated
        timelimit = 1_000 # time limit for simulation
        seed = 8710 # rng seed
        s1 = UQ.Sim(ncust, timelimit, seed)
        q1 = UQ.Queue(c, adist, sdist)
        df1 = UQ.runsim(s1, q1)
        noiseout = UQ.disorder(df1, 2)
        @test count(noiseout[:DepartOrderErr] - noiseout[:DepartOrder] .== 0.0) / size(noiseout,1) ≈ 0.9 atol = 0.1
        @test all(sort(noiseout, :DepartMeas)[:DepartOrder] .== collect(1:size(noiseout,1)))
        @test all(sort(noiseout, :DepartMeas)[:DepartOrder] .== sort(noiseout, :Depart)[:DepartOrder])

        # Deltamax algorithm
        o1 = [1,2,3,5,4]
        o2 = [1,3,6,5,4]
        @test UQ.c_deltamax(o1) == 2 && UQ.c_deltamax(o2) == 3

        # Order-based algorithm
        c = 15 # number of servers
        μa = (1 / c) * 0.99
        μs = 1 / c
        σ = 0.9
        adist = LogNormal(exp(1/μ), σ) # arrival distribution
        sdist = LogNormal(exp(1/μ), σ) # service distribution (single server)
        ncust = 10_000 # total number of customers generated
        timelimit = 10_000 # time limit for simulation
        seed = 3001 # rng seed
        s2 = UQ.Sim(ncust, timelimit, seed)
        q2 = UQ.Queue(c, adist, sdist)
        df2 = UQ.runsim(s2, q2)
        outorder1 = sort(df1, :dorder)[:aorder]
        outorder2 = sort(df2, :dorder)[:aorder]
        outorder3 = randperm(1_000)
        @test UQ.c_order(outorder3) == UQ.c_order_slow(outorder3)
        @test UQ.c_order([1,3,5]) == 3
        @test UQ.c_order(outorder1) == 2

        # Variance and Uninformed algorithms
        cmax = 19
        A1, A2 = df1[:atime][1:1000], df2[:atime][1:1000]
        D1, D2 = df1[:dtime][1:1000], df2[:dtime][1:1000]
        (cvar1, cunf1, VS1) = UQ.c_var_unf(A1, D1, cmax)
        (cvar2, cunf2, VS2) = UQ.c_var_unf(A2, D2, cmax)
        @test cvar1 == 2 && cunf1 == 2

        # Build distributions
        (adist, sdist) = UQ.builddist(2, "Uniform", "LogNormal", 0.99)
        (adist2, sdist2) = UQ.builddist(15, "Exponential", "Beta", 0.9)
        @test (1/mean(adist)) / (2 * (1/mean(sdist))) ≈ 0.99 atol = 0.0001
        @test (1/mean(adist2)) / (15 * (1/mean(sdist2))) ≈ 0.9 atol = 0.0001
    end

    @testset "Metrics" begin
        # Convergence estimation
        settings = ["Order" "Exponential" "Exponential" 400 9 0.9;
                    "Order" "Uniform" "LogNormal" 400 15 0.9]
        n_runs = 10 # limit runs for debugging and timing
        n_methods = 3 # number of methods to compare
        max_servers = 19 # max number of servers
        obs_max = 1_000 # max observations available
        time_limit = 10_000 # max simulation time
        ϵ = 0.05 # convergence quality
        window = 20 # observation window for convergence estimate
        window_detail = 50 # observation window for error estimate
        step = 20 # how many observations to skip while calculating convergence
        seed = 8710
        param_inf = Paraminf(settings, n_runs, n_methods, max_servers, obs_max, time_limit, ϵ, window, window_detail, step, seed)
        (conv, conv_meas, ests, ests_meas) = UQ.convergence(param_inf,1)
        @test (conv[1] > 0) && (conv[1] < 400) && ((conv[2] == 0) || conv[2] > conv[1])
        @test (conv[1] == conv_meas[1]) && (conv_meas[2] == 0)

        # Error estimation
        settings = ["Order" "Exponential" "Exponential" 400 9 0.9;
                    "Order" "Uniform" "LogNormal" 400 15 0.9;
                    "Order" "LogNormal" "LogNormal" 60 9 0.9]
        n_runs = 10 # limit runs for debugging and timing
        n_methods = 3 # number of methods to compare
        max_servers = 19 # max number of servers
        obs_max = 1_000 # max observations available
        time_limit = 10_000 # max simulation time
        ϵ = 0.05 # convergence quality
        window = 20 # observation window for convergence estimate
        window_detail = 50 # observation window for error estimate
        step = 20 # how many observations to skip while calculating convergence
        seed = 8710
        param_inf = Paraminf(settings, n_runs, n_methods, max_servers, obs_max, time_limit, ϵ, window, window_detail, step, seed)
        (err, err_meas, detail, detail_meas) = UQ.esterror(param_inf, 3)
        @test err[1] == 0 && err[2] >= 1
        @test err_meas[1] == 0 && err_meas[2] >= err[2]

        # Comparison
        settings = ["Order" "Exponential" "Exponential" 400 9 0.9;
                    "Order" "Uniform" "LogNormal" 400 15 0.9;
                    "Order" "LogNormal" "LogNormal" 60 9 0.9]
        n_runs = size(settings,1) # experimental runs
        n_methods = 3 # number of methods to compare
        max_servers = 19 # max number of servers
        obs_max = 1_000 # max observations available
        time_limit = 10_000 # max simulation time
        ϵ = 0.05 # convergence quality
        window = 20 # observation window for convergence estimate
        window_detail = 50 # observation window for error estimate
        step = 20 # how many observations to skip while calculating convergence
        seed = 8710
        param_inf = Paraminf(settings, n_runs, n_methods, max_servers, obs_max, time_limit, ϵ, window, window_detail, step, seed)
        (rerr, rconv, rraw) = infer(param_inf)
        @test mean(rerr[1][:,1]) < mean(rerr[1][:,2])
        @test mean(rerr[2][:,1]) < mean(rerr[2][:,2])
        @test (rconv[1][1,1] < rconv[1][1,2]) || (rconv[1][1,2] == 0)
        @test (rconv[2][1,1] < rconv[2][1,2]) || (rconv[2][1,2] == 0)
    end
end

@testset "LCFS" begin
    # LCFS Order-based algorithm
    c = 7 # number of servers
    μ = 1 / c # service rate (mean service time is 1/μ) (single server)
    λ = 0.9 # arrival rate (mean interarrival time is 1/λ)
    adist = Exponential(1/λ) # arrival distribution
    sdist = Exponential(1/μ) # service distribution (single server)
    ncust = 10_000 # total number of customers generated
    timelimit = 10_000 # time limit for simulation
    seed = 3001 # rng seed
    s2 = UQ.Sim(ncust, timelimit, seed)
    q2 = UQ.Queue(c, adist, sdist)
    df = UQ.runsimLCFS(s2, q2)
    lcfsorder, dfl = UQ.lcfs(df, :dtime)
    chat = UQ.c_order_LCFS(lcfsorder)
    @test chat == c

    # LCFS Variance, Uninformed algorithm
    c = 7 # number of servers
    λ = 0.99 # arrival rate (mean interarrival time is 1/λ)
    μ = 1 / c # service rate (mean service time is 1/μ) (single server)
    adist = Exponential(1/λ) # arrival distribution
    sdist = Exponential(1/μ) # service distribution (single server)
    ncust = 500 # total number of customers generated
    timelimit = 10_000 # time limit for simulation
    seed = 1302 # rng seed
    cmax = 19
    s2 = UQ.Sim(ncust, timelimit, seed)
    q2 = UQ.Queue(c, adist, sdist)
    df1 = UQ.runsimLCFS(s2, q2)
    A1, D1 = df1[:atime], df1[:dtime]
    (cvar1, cunf1, VS1) = UQ.c_var_unf_LCFS(df1, A1, D1, cmax)
    @test cvar1 == c
    @test cunf1 == c

    # LCFS first depature function
    df = DataFrame()
    df[:dtime] = [1.1,1.2,1.3,1.4,1.5,1.6]
    df[:nd] = [1,2,3,4,5,6]
    @test UQ.firstdepart(df, 5, 1.2) == 1.5
    @test UQ.firstdepart(df, 5, 1.6) == 1.6

    # LCFS Introduce observation error into true data
    c = 2 # number of servers
    λ = 0.9 # arrival rate (mean interarrival time is 1/λ)
    μ = 1 / c # service rate (mean service time is 1/μ) (single server)
    adist = Exponential(1/λ) # arrival distribution
    sdist = Exponential(1/μ) # service distribution (single server)
    ncust = 10_000 # total number of customers generated
    timelimit = 1_000 # time limit for simulation
    seed = 8710 # rng seed
    s1 = UQ.Sim(ncust, timelimit, seed)
    q1 = UQ.Queue(c, adist, sdist)
    df1 = UQ.runsimLCFS(s1, q1)
    noiseout = UQ.disorderLCFS(df1)
    @test all(sort(noiseout, :dmeas)[:dorder] .== collect(1:size(noiseout,1)))
    @test all(sort(noiseout, :dmeas)[:dorder] .== sort(noiseout, :dtime)[:dorder])

    # LCFS Convergence estimation
    settings = ["Order" "Exponential" "Exponential" 400 9 0.9;
                "Order" "Uniform" "LogNormal" 400 15 0.9]
    n_runs = 10 # limit runs for debugging and timing
    n_methods = 3 # number of methods to compare
    max_servers = 19 # max number of servers
    obs_max = 1_000 # max observations available
    time_limit = 10_000 # max simulation time
    ϵ = 0.05 # convergence quality
    window = 20 # observation window for convergence estimate
    window_detail = 50 # observation window for error estimate
    step = 20 # how many observations to skip while calculating convergence
    seed = 8710
    param_inf = Paraminf(settings, n_runs, n_methods, max_servers, obs_max, time_limit, ϵ, window, window_detail, step, seed)
    @time (conv, conv_meas, ests, ests_meas) = UQ.convergenceLCFS(param_inf,1)
    @test (conv[1] > 0) && (conv[1] < 400) && ((conv[2] == 0) || conv[2] >= conv[1])
    @test (conv[1] == conv_meas[1]) && (conv_meas[2] >= conv_meas[1])

    # LCFS Error estimation
    settings = ["Order" "Exponential" "Exponential" 400 9 0.9;
                "Order" "Uniform" "LogNormal" 400 15 0.9;
                "Order" "LogNormal" "LogNormal" 60 9 0.9]
    n_runs = 10 # limit runs for debugging and timing
    n_methods = 3 # number of methods to compare
    max_servers = 19 # max number of servers
    obs_max = 1_000 # max observations available
    time_limit = 10_000 # max simulation time
    ϵ = 0.05 # convergence quality
    window = 20 # observation window for convergence estimate
    window_detail = 50 # observation window for error estimate
    step = 20 # how many observations to skip while calculating convergence
    seed = 8710
    param_inf = Paraminf(settings, n_runs, n_methods, max_servers, obs_max, time_limit, ϵ, window, window_detail, step, seed)
    @time (err, err_meas, detail, detail_meas) = UQ.esterrorLCFS(param_inf, 3)
    @test err[1] < err[2]
    @test err_meas[1] < err_meas[2]

    # LCFS Comparison
    settings = ["Order" "Exponential" "Exponential" 400 9 0.9;
                "Order" "Uniform" "LogNormal" 400 15 0.9;
                "Order" "LogNormal" "LogNormal" 60 9 0.9]
    n_runs = size(settings,1) # experimental runs
    n_methods = 3 # number of methods to compare
    max_servers = 19 # max number of servers
    obs_max = 1_000 # max observations available
    time_limit = 10_000 # max simulation time
    ϵ = 0.05 # convergence quality
    window = 20 # observation window for convergence estimate
    window_detail = 50 # observation window for error estimate
    step = 20 # how many observations to skip while calculating convergence
    seed = 8710
    param_inf = Paraminf(settings, n_runs, n_methods, max_servers, obs_max, time_limit, ϵ, window, window_detail, step, seed)
    @time (rerr, rconv, rraw) = inferLCFS(param_inf)
    @test mean(rerr[1][:,1]) < mean(rerr[1][:,2])
    @test mean(rerr[2][:,1]) < mean(rerr[2][:,2])
    @test (rconv[1][1,1] <= rconv[1][1,2]) || (rconv[1][1,2] == 0)
    @test (rconv[2][1,1] <= rconv[2][1,2]) || (rconv[2][1,2] == 0)
end


text:
# ------------------------------------------------------------------------------------------
# ## Create a caesar cipher
#
# ### Caesar ciphers
#
# A caesar cipher is an encryption scheme that shifts all letters in the alphabet by some
# specified offset to other letters in the alphabet.
#
# For example, a shift of 1 would turn the letter "A" into the letter "B" and the letter "M"
# to the letter "N".
#
# ### Goal
#
# We want to add a method to the `+` operator such that we can add together a string and an
# integer shift to encrypt a message. For example,
#
# ```julia
# 4 + "hello" == "lipps"
# ```
#
# We'll walk you through the steps to do this!
#
# ### Test it out
#
# Once you think you have it working, try to decrypt the following string by adding a shift
# of -7.
# ```julia
# "Kv'uv{'tlkksl'pu'{ol'hmmhpyz'vm'kyhnvuz'mvy'\u80v|'hyl'jy|ujo\u80'huk'{hz{l'nvvk'~p{o'rl{
# jo|w5"
# ```
#
# ### Let's get started!
#
# First, we want a way to convert between characters and integers. Actually, under the hood,
# all of our characters are being represented as numbers via their *ASCII representation*.
#
# You can start to get a feel for how this works by running the following lines of code.
#
# ```julia
# convert(Int, 'a')
# convert(Int, 'b')
# convert(Char, 97)
# convert(Char, 98)
# ```
# ------------------------------------------------------------------------------------------

convert(Int, 'a')

convert(Int, 'b')

convert(Char, 97)

convert(Char, 98)

# ------------------------------------------------------------------------------------------
# What happens when you try to add an integer to a character? (Note that the difference
# between `Char`s and `String`s is important here!)
# ------------------------------------------------------------------------------------------

'a' + 1

'm' + 1

# Basically our caesar cipher already works for characters.
# We still need to add a method to `+` that makes it work for strings!

# ------------------------------------------------------------------------------------------
# When we treat a string elementwise, what is the type of (`typeof`) each element?
# ------------------------------------------------------------------------------------------

[typeof(x) for x in "Mew"]

map(x -> typeof(x), collect("Mew"))

# ------------------------------------------------------------------------------------------
# Try to write a function called `caesar(shift, stringin)` that encodes its input string,
# `stringin`, by shifting all letters in the alphabet by `shift`.
#
# One way to do this is to use the `map` or `broadcast` function!
# ------------------------------------------------------------------------------------------

caesar(shift, stringin) = map(x -> x + shift, stringin)

# ------------------------------------------------------------------------------------------
# If you think you have this working, try out
# ```julia
# caesar(-4, "lipps")
# ```
# ------------------------------------------------------------------------------------------

caesar(-4, "lipps")

# ------------------------------------------------------------------------------------------
# Last, we want to extend the `+` operator to include a way to apply this cipher.
#
# The `+` operator lives in a place called "Base". Everything that lives in Base is
# accessible to us as users by default, but we need a special incantation to modify the
# things that live in Base. If we want to modify `+`, our incantation is
#
# ```julia
# import Base: +
# ```
# ------------------------------------------------------------------------------------------

import Base: +

# ------------------------------------------------------------------------------------------
# After you have imported `+` from Base, we are ready to modify it.
#
# If you're not sure how to add a method for `+`, let's go through an example first.
#
# Let's say we want to add a method for `+` that allows us to concatenate strings so that we
# can say
# ```julia
# "hello " + "world"
# ```
# and receive the output string ```hello world```.
#
# To do this, we can write a method for `+` like this
# ------------------------------------------------------------------------------------------

+(x::String, y::String) = string(x, y)

# ------------------------------------------------------------------------------------------
# Note that I've added type information about `x` and `y` with a double colon and the word
# `String` to say that `x` and `y` are both strings.
#
# To test that this works, try out
# ```julia
# "hello " + "world"
# ```
# ------------------------------------------------------------------------------------------

"hello " + "world"

# ------------------------------------------------------------------------------------------
# And now that you've extended `+` once, let's add another method for `+` that calls the
# `caesar` function we've written.
# ------------------------------------------------------------------------------------------

+(shift::Int, stringin::String) = caesar(shift, stringin)

# ------------------------------------------------------------------------------------------
# Now try it out via
#
# ```julia
# -7 + "Kv'uv{'tlkksl'pu'{ol'hmmhpyz'vm'kyhnvuz'mvy'\u80v|'hyl'jy|ujo\u80'huk'{hz{l'nvvk'~p{
# o'rl{jo|w5"
# ```
# ------------------------------------------------------------------------------------------

-7 + "Kv'uv{'tlkksl'pu'{ol'hmmhpyz'vm'kyhnvuz'mvy'\u80v|'hyl'jy|ujo\u80'huk'{hz{l'nvvk'~p{o'rl{jo|w5"

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
# Run a  test from the GALLERY, tests or doc/scripts
#
# gmtest(test) run an example from the gallery were TEST is is either 'ex01' ... 'ex45'
#
# gmtest(test, test_dir) runs a test from the tests suit. TEST_DIR is the directory name where TEST live
#		example: gmtest("poldecimate","gmtspatial")
#
# gmtest(test, test_dir, family) runs a test from a certainly family of tests.
#	Currently only family = "scripts" is implemented (means run test from the 'scripts' dir. TEST_DIR is ignored
#		example: gmtest("GMT_insert", "", "scripts")

using Compat
import Compat.String

global g_root_dir, out_path, GM
include("gallery.jl")

# Edit these for your own needs
g_root_dir = "C:/progs_cygw/GMTdev/gmt5/trunk/"
out_path   = "V:/"		# Set this if you want to save the PS files in a prticular place
GM         = "C:/programs/GraphicsMagick/gm"

function run_tests(what::String)
	if (what == "all_examples")
		for k = 1:46
			ex = @sprintf("ex%02d", k)
			gmtest(ex)
		end
	elseif (what == "all_tests")
		do_tests()
	elseif (what[1] == 'e' && what[2] == 'x')
		gmtest(what)
	end
end

# ----------------------------------------------------------------------------------------
gmtest(test, test_dir::String, family::String) = gmtest(test, test_dir, family, 3)
function gmtest(test, test_dir="", family="", nargin::Int=1)
# Run an example or test and print its tested status, i.e. if it PASS or FAIL

	global g_root_dir, out_path, GM

	if ((nargin == 1) && (length(test) == 4) && startswith(test, "ex"))	# Run example from gallery
		ps, orig_path = gallery(test, g_root_dir, out_path)
	else
		ps, orig_path = call_test(test, test_dir, g_root_dir, out_path, family)
		if (isa(ps, Bool))
			if (ps)		println("Test " * test * " PASS")
			else 		println("Test " * test * " FAIL")
			end
			return
		end
	end

	if (isempty(ps))
		println("    Test -> " * test * " <- does not exist or some other error occurred")
		return
	end
	pato, fname = fileparts(ps)
	png_name = pato * fname * ".png"
	ps_orig  = orig_path * fname * ".ps"

	# Compare the ps file with its original.
	cm = `$GM compare -density 200 -maximum-error 0.001 -highlight-color magenta -highlight-style
		 assign -metric rmse -file $png_name $ps_orig $ps`

	run(pipeline(ignorestatus(cm), stdout=DevNull, stderr="errs.txt"))
	t = readstring("errs.txt")
	rm("errs.txt")
	if (isempty(t))
		println("    Test " * test * " PASS")
		rm(png_name)
	else
		ind = searchindex(t, "image")
		if (ind == 0)
			println("    Test " * test * " FAIL with:\n " * t)		# Quite likely a Ghostscript error
		else
			println("    Test " * test * " FAIL with: " * t[ind:end-1])
		end
	end
end

# ---------------------------------------------------------------------------------------
function call_test(test, test_dir, g_root_dir, out_path, family)
# Here PS hold the full name of the created PS file and ORIG_PATH must
# contain the path to where the original postscript file, against which the PS will be compared

	ps = [];	orig_path = [];					# Defaults for the case of error
	if (isempty(family))						# A test from the tests suit
		pato = g_root_dir * "test/" * test_dir * "/"
	elseif (family == "scripts")				# A test for the 'scripts' directory
		pato = g_root_dir * "doc/scripts/ml/"
	else
		@sprintf("Family %s is unknown or not yet implemented", family)
		return
	end

	try
		#ps, orig_path = evalfile(pato * test * ".jl")(out_path)	# Load & run the test file. Errors because UndefVarError: fileparts not defined
		include(pato * test * ".jl")					# Load the test file
		ps, orig_path = eval(Symbol(test))(out_path)	# and now run it
		if (family == "scripts")
			orig_path = orig_path[1:end-3]				# Because original PS lieve in a subdir below
		end
	catch
		println("Error executing test script")
	end
	return ps, orig_path
end

# ---------------------------------------------------------------------------------------
function do_tests()
	tests = Any["poldecimate" "gmtspatial";
	            "spheres" "potential";
	            "measure" "gmtspatial"]

	for k = 1:size(tests,1)
		gmtest(tests[k,1], tests[k,2])
	end
end

# ---------------------------------------------------------------------------------------
function fileparts(file)
	path = dirname(file)
	f = basename(file)
	fname, ext = splitext(f)
	return path, fname, ext
end

# ---------------------------------------------------------------------------------------
function addpath(path)
	push!(LOAD_PATH, path)
end

# ---------------------------------------------------------------------------------------
function mfilename(arg)
# If ARG == "fullpath" return the full path of the script being run, otherwise only its name.
# But the problem with this function is that it WONT WORK because it always return the path to gmtest.jl
	t = @__FILE__
	path, fname = fileparts(t)
	if (arg == "fullpath")
		return path * fname
	else
		return fname
	end
end

mfilename() = mfilename("treta")	# Output only the script name

# ---------------------------------------------------------------------------------------
feval(fn_str, args...) = eval(parse(fn_str))(args...)	# This is not a good one.

# ---------------------------------------------------------------------------------------
function strfind(str::AbstractString, sub_str)
# To mimic Matalab's function except that it returns [0] instead of [] in case of failure
# and it always returns a vector.
	ind = [search(str, sub_str)]
	ind_n = ind[1]
	while (ind_n != 0)
		ind_n = searchindex(str, sub_str, ind_n+1)
		if (ind_n != 0)		# Found another
			push!(ind, ind_n)
		end
	end
	return ind
end


text:
struct InformationRewardBeliefMDP{P<:POMDP, R, S, B, A} <: MDP{InformationBelief{B}, A}
    pomdp::P
    resample::R
    im::InformationMeasure
    discount_information_gain::Bool
    max_frac_replaced::Float64
    rng::AbstractRNG
    _pm::Vector{S}
    _rm::Vector{Float64}
    _wm::Vector{Float64}
end

function InformationRewardBeliefMDP(pomdp::P, resampler, im, discount_information_gain, max_frac_replaced, rng=Random.GLOBAL_RNG) where {P<:POMDP}
    S = statetype(pomdp)
    b0 = resample(resampler, initialstate_distribution(pomdp), rng)
    InformationRewardBeliefMDP{P, typeof(resampler), S, typeof(b0), actiontype(pomdp)}(
        pomdp,
        resampler,
        im,
        discount_information_gain,
        max_frac_replaced,
        rng,
        S[],
        Float64[],
        Float64[]
        )
end

function POMDPs.generate_sr(bmdp::InformationRewardBeliefMDP{P,R,S,B,A}, ib::InformationBelief{B}, a, rng::AbstractRNG) where {P,R,S,B,A}
    b = ib.b

    if n_particles(b) > bmdp.resample.n
        b = resample(bmdp.resample, b, rng)
    end
    s_o = rand(rng, Base.filter(s->!isterminal(bmdp.pomdp, s), particles(b)))
    sp_o, o = generate_so(bmdp.pomdp, s_o, a, rng)
    while true
        sp_o, o = generate_so(bmdp.pomdp, s_o, a, rng)
        obs_weight(bmdp.pomdp, s_o, a, sp_o, o) > eps(1.0/bmdp.resample.n) && break
        # if obs_weight(bmdp.pomdp, s_o, a, sp_o, o) > eps(1.0/bmdp.resample.n)
        #     break
        # end
    end
    obs_w = obs_weight(bmdp.pomdp, s_o, a, sp_o, o)

    ps = particles(b)
    pm = bmdp._pm
    rm = bmdp._rm
    wm = bmdp._wm
    resize!(pm, 0)
    resize!(rm, 0)
    resize!(wm, 0)

    all_terminal = true
    for i in 1:n_particles(b)
        s = ps[i]
        if !isterminal(bmdp.pomdp, s)
            all_terminal = false
            sp, r = generate_sr(bmdp.pomdp, s, a, rng)
            if sp isa Real && !isfinite(sp)
                @warn "Generated sp is not finite"
                @show sp s a transition(bmdp.pomdp,s,a)
            end
            push!(pm, sp)
            push!(rm, r)
            push!(wm, obs_weight(bmdp.pomdp, s, a, sp, o))
        end
    end
    ws = sum(wm)
    if all_terminal
        @warn("All states in particle collection were terminal.")
        @debug b length(particles(b)) s_o a sp_o o pm wm ws obs_w
        # return (resample(bmdp.resample, reset_distribution(bmdp.pomdp, b, a, o), rng), ib.i), 0.0
        return InformationBelief{B}(resample(bmdp.resample, initialstate_distribution(bmdp.pomdp), rng), ib.i), [0.0, 0.0]
    elseif ws < eps(1.0/length(wm))
        # @warn("Observation weights are low, ws=",ws)
        @debug b length(particles(b)) s_o a sp_o o pm wm ws obs_w
        # return (resample(bmdp.resample, reset_distribution(bmdp.pomdp, b, a, o), rng), ib.i), 0.0
        wm = ones(length(wm))
        ws = sum(wm)
    end

    # Information Measure calculation
    bpred = ParticleCollection(pm)
    wbp = WeightedParticleBelief(pm, wm/ws)
    pc = resample(bmdp.resample, wbp, rng)

    ip = @inferred information(bmdp.im, bp=pc, wbp=wbp, b=b, bpred=bpred, a=a, o=o)
    @assert isfinite(ip)

    ps = particles(pc)
    mpw = max_possible_weight(bmdp.pomdp, a, o)
    frac_replaced = bmdp.max_frac_replaced*max(0.0, 1.0 - maximum(wm)/mpw)
    n_replaced = floor(Int, frac_replaced*length(ps))
    is = randperm(rng, length(ps))[1:n_replaced]
    for i in is
        ps[i] = new_particle(bmdp.pomdp, b, a, o, rng)
    end

    d = bmdp.discount_information_gain ? discount(bmdp) : 1
    InformationBelief{B}(pc, ip), Array{Float64}([mean(rm,weights(wm)), d * ip - ib.i])
end

# function initialstate(bmdp::InformationRewardBeliefMDP, rng::AbstractRNG)
#     return resample(bmdp.resample, initialstate_distribution(bmdp.pomdp), rng)
# end

POMDPs.actions(bmdp::InformationRewardBeliefMDP{P,R,S,B,A}, ib::InformationBelief{B}) where {P,R,S,B,A} = actions(bmdp.pomdp, ib.b)
POMDPs.actions(bmdp::InformationRewardBeliefMDP) = actions(bmdp.pomdp)

POMDPs.isterminal_act(bmdp::InformationRewardBeliefMDP{P,R,S,B,A}, a::A) where {P,R,S,B,A} = isterminal_act(bmdp.pomdp, a)
POMDPs.isterminal(bmdp::InformationRewardBeliefMDP{P,R,S,B,A}, ib::InformationBelief{B}) where {P,R,S,B,A} = all(isterminal(bmdp.pomdp, s) for s in particles(ib.b))
# POMDPs.isterminal(bmdp::InformationRewardBeliefMDP, ib) = all(isterminal(bmdp.pomdp, s) for s in particles(ib.b))

POMDPs.discount(bmdp::InformationRewardBeliefMDP) = discount(bmdp.pomdp)

function POMDPSimulators.simulate(sim::RolloutSimulator, mdp::InformationRewardBeliefMDP{P,R,S,B,A}, policy::Policy, initialstate::InformationBelief{B}) where {P,R,S,B,A}
    if isterminal(mdp, initialstate)
        return [0.0, 0.0]
    end

    if sim.eps == nothing
        eps = 0.0
    else
        eps = sim.eps
    end

    if sim.max_steps == nothing
        max_steps = typemax(Int)
    else
        max_steps = sim.max_steps
    end

    s = initialstate

    disc = 1.0
    r_total = [0.0, 0.0]
    step = 1

    while disc > eps && step <= max_steps
        a = action(policy, s.b)

        sp, r = @inferred generate_sr(mdp, s, a, sim.rng)

        r_total += disc*r

        s = sp

        disc *= discount(mdp)
        step += 1

        if isterminal(mdp, s, a)
            break
        end
    end

    return r_total
end


text:
using hypergalaopt,Optim
using Test

@testset "hyper-gala-opt.jl" begin
        arg_item = 10
        arg_optimize = Optim.SimulatedAnnealing()
        #arg_optimize = [Optim.SimulatedAnnealing(), Optim.BFGS()]
        arg_function = f(a;c=10) = sum(@. 100 + (a-3)^2 + (c-100)^2)
        arg_label=1
        hypergalaop(arg_item,arg_optimize,arg_function,arg_label)
end


text:
"""
AbstractSystem

An abstract type representing a system
"""
abstract type AbstractSystem end

"""
state_dim(s::AbstractSystem)

Return the state dimension of a system
"""
state_dim(s::AbstractSystem) = s.state_dim

"""
state_set(s::AbstractSystem)

Return the state dimension of a system
"""
state_set(s::AbstractSystem) = s.X

"""
has_input(s::AbstractSystem)

Return the input condition of a system
"""
has_input(s::AbstractSystem) = false

"""
has_uncertainty(s::AbstractSystem)

Return the uncertainty condition of a system
"""
has_uncertainty(s::AbstractSystem) = false

"""
AbstractUncertainSystem

An abstract type representing an uncertain system
"""
abstract type AbstractUncertainSystem <:AbstractSystem end

"""
uncertainty_dim(s::AbstractUncertainSystem)

Return the uncertainty dimension of a system
"""
uncertainty_dim(s::AbstractUncertainSystem) = s.uncertainty_dim

"""
has_uncertainty(s::AbstractUncertainSystem)

Return the uncertainty condition of a system
"""
has_uncertainty(s::AbstractUncertainSystem) = true

"""
uncertainty_set(s::AbstractUncertainSystem)

Return the uncertainty set of a system
"""
uncertainty_set(s::AbstractUncertainSystem) = s.W

"""
AbstractControlSystem

An abstract type representing a control system
"""
abstract type AbstractControlSystem <:AbstractSystem end

"""
input_dim(s::AbstractControlSystem)

Return the input dimension of a control system
"""
input_dim(s::AbstractControlSystem) = s.input_dim

"""
has_input(s::AbstractControlSystem)

Return the input condition of a system
"""
has_input(s::AbstractControlSystem) = true

"""
input_set(s::AbstractControlSystem)

Return the input constraint of a control system
"""
input_set(s::AbstractControlSystem) = s.U

"""
AbstractUncertainControlSystem

An abstract type representing a control system with uncertainties
"""
abstract type AbstractUncertainControlSystem <:AbstractUncertainSystem end

"""
input_dim(s::AbstractUncertainControlSystem)

Return the input dimension of a control system
"""
input_dim(s::AbstractUncertainControlSystem) = s.input_dim

"""
has_input(s::AbstractUncertainControlSystem)

Return the input condition of a system
"""
has_input(s::AbstractUncertainControlSystem) = true

"""
input_set(s::AbstractUncertainControlSystem)

Return the input constraint of a control system
"""
input_set(s::AbstractUncertainControlSystem) = s.U

"""
AbstractSolver

An abstract solver object
"""
abstract type AbstractSolver end

"""
AbstractOptions

An abstract options object
"""
abstract type AbstractOptions end

"""
AbstractSolution

An abstract solution object
"""
abstract type AbstractSolution end


text:
#=
Credits: https://discourse.julialang.org/t/swap-array-contents/7774/12?page=2
Professor Steven @stevengj (physics) MIT lecturer, & an active julia member

=#
#---------------------------------
function swap!(x::Float64, y::Float64)
    x, y = y, x
end
#----------------------------------
function swap!(x::Int64, y::Int64)
    x, y = y, x
end
#-------------------------------------
"""
A wrapper
around Array{T,N}
but in application: it could be anything, simply put

```input:
WrapArray{T,N}: inherits defintions from  `AbstractArray{T,N}`, has a type template `T`, & size `N`
````

```output:
 a::Array{T,N} : which has a type of template `T`, & Size N
```
"""
mutable struct WrapArray{T,N} <: AbstractArray{T,N}
    a::Array{T,N}
end
T = Float64;
N = 1;
#------------
"""
```input:
    x: original x value
    y: original y value
```
```output:
    x: swapped x value
    y: swapped y value
```
"""

function oldschoolSwap(x, y)
    tmp = x
    x = y
    y = tmp
    return x, y
end

#---------------------------------------------------
function swap!(x::Array{T,N}, y::Array{T,N}) where {T,N}
    x.a, y.a = y.a, x.a # a is for any arbitrary item
    return x, y
end
#--------------------------------------------------
function swap!(x::VecOrMat, y::VecOrMat) #where {T,N}
    x.a, y.a = y.a, x.a # a is for any arbitrary item
    return x, y
end
#----------------------------------------------
#=
function swap!(x::WrapArray{T,N}, y::WrapArray{T,N}) where {T,N}
    x.a, y.a = y.a, x.a # a is for any arbitrary item
    return x, y
end

#footnote:
Since this is a mutable type, the changes to the contents of x and y would be reflected in all references to those objects elsewhere. As you point out, above, though, this adds an extra indirection to array accesses.

=#
""" A type-Matching function chooses the right function based on given type

```input:
x:  x generic , could be of type Float64 or Array{T,N}
y:
```

```output:
swap!(x,y) function
```

Returns `swap!(x::Float64, y::Float64)` or  swap!(x::Array{T,N}, y::Array{T,N})
"""
function swapping(x, y)
    if typeof(x) == typeof(y) == Float64 # && typeof(y) == Float64
        return swap!(x::Float64, y::Float64) #
    elseif typeof(x) == Array{T,N} == typeof(y) # && typeof(y) = Array{T,N}
        return swap!(x::Array{T,N}, y::Array{T,N})
    else
        println("ERROR: faulty input argiments, or unidentified ERROR")
    end
end

#= Redundant: reason: function definiton above includes the definition below

function swap!(x::WrapArray{T,1}, y::WrapArray{T,1}) where {T,N}
    x.a, y.a = y.a, x.a # a is for any arbitrary item
    return x, y
end

=#

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
module AbstractNumbers

using SpecialFunctions

# Yeah yeah, how ironic
abstract type AbstractNumber{T} <: Number end

# convert(Number, my_type) & is the only interface one needs to overload
@inline basetype(T::Type{<: AbstractNumber}) = error("AbstractNumbers.basetype not implemented for $T")
@inline basetype(x::T) where T <: AbstractNumber = basetype(T)
@inline number(x::AbstractNumber) = error("AbstractNumbers.number not implemented for $(typeof(x))")

@inline function Base.convert(::Type{AN}, x::AbstractNumber{T2}) where {AN <: AbstractNumber{T}, T2} where T
    AN(T(number(x)))
end

@inline function Base.convert(::Type{AN}, x::Number) where AN <: AbstractNumber{T} where T
    AN(x)
end
@inline function Base.convert(::Type{T}, x::AbstractNumber) where T <: Number
    if isabstracttype(T)
        return x
    else
        return convert(T, number(x))
    end
end

"""
    like(num::AbstractNumber, x::T)

Creates a number from `x` like the first argument. It discards the eltype of `num`
and uses the type of `x` instead.
usage:
    ```
    like(MyAbstractNumber(1f0), 22) === MyAbstractNumber{Int}(22)
    ```
"""
@inline like(num::AN, x::T) where {AN <: AbstractNumber, T} = like(AN, x)
@inline like(::Type{AN}, x::T) where {AN <: AbstractNumber, T} = basetype(AN){T}(x)
# Not sure if this is the most performant and correct way, but this enables
# that == and friends return a bool. Returning an AbstractNumber{Bool} seems a bit odd
@inline like(::Type{AN}, x::Tuple) where AN <: AbstractNumber = like.(AN, x)

@inline Base.typemax(::Type{Num}) where Num <: AbstractNumber{T} where T = like(Num, typemax(T))
@inline Base.typemin(::Type{Num}) where Num <: AbstractNumber{T} where T = like(Num, typemin(T))

@inline Base.typemax(::T) where T <: AbstractNumber = typemax(T)
@inline Base.typemin(::T) where T <: AbstractNumber = typemax(T)

@inline Base.eltype(x::T) where T <: AbstractNumber = T
@inline Base.eltype(::Type{T}) where T <: AbstractNumber = T

include("overloads.jl")

rem(x::AbstractNumber, y::AbstractNumber, r::RoundingMode) = like(x, rem(number(x), number(y), r))

# Overload ambiguities
Base.:^(a::Irrational{:ℯ}, b::AbstractNumber) = like(a, Base.:^(a, number(b)))
Base.:^(a::AbstractNumber, b::Rational) = like(a, Base.:^(number(a), b))
Base.:^(a::AbstractNumber, b::Integer) = like(a, Base.:^(number(a), b))
Base.log(a::Irrational{:ℯ}, b::AbstractNumber) = like(a, Base.:^(a, number(b)))

export AbstractNumber

end # module


text:
module Footings

export
    AbstractElement,
    Parameter,
    Temporary,
    Result,
    Step,
    Model,
    add_parameter!,
    add_temporary!,
    add_result!,
    add_step!,
    parse_expr,
    @model,
    build_documentation,
    build_runner,
    build_auditor

include("model.jl")
# inlcude("audit_runner.jl")
# include("model_runner.jl")

end


text:
# functions specific to creating "batched" fields, which are fields that
# simultaneously store and operate on several maps at a time. this is mainly
# useful on GPUs where operating on batches of fields is sometimes no slower
# than a single field.

"""
    batch(fs::FlatField...)
    batch(fs::Vector{<:FlatField})
    batch(fs::TUple{<:FlatField})

Turn a length-N array of `FlatField`'s into a single batch-length-N `FlatField`.
For the inverse operation, see [`unbatch`](@ref).
"""
batch(f::F) where {N,θ,∂m,F<:FlatS0{<:Flat{N,θ,∂m}}} = f
batch(fs::F...) where {N,θ,∂m,F<:FlatS0{<:Flat{N,θ,∂m}}} =
    basetype(F){Flat{N,θ,∂m,length(fs)}}(cat(map(firstfield,fs)..., dims=3))
batch(fs::F...) where {F<:Union{FlatS2,FlatS02}} =
    FieldTuple{basis(F)}(map(batch, map(firstfield,fs)...))
batch(fs::Union{Vector{<:FlatField},Tuple{<:FlatField}}) = batch(fs...)

"""
    batch(f::FlatField, D::Int)

Construct a batch-length-`D` `FlatField` from an unbatched `FlatField` which
will broadcast as if it were `D` copies of `f` (without actually making `D`
copies of the data in `f`)
"""
batch(f::F, D::Int) where {N,θ,∂m,D′,F<:FlatS0{Flat{N,θ,∂m,D′}}} =
    (D′==D || D′==1) ? basetype(F){Flat{N,θ,∂m,D}}(firstfield(f)) : error("Can't change batch-length from $(D′) to $D.")
batch(f::F, D::Int) where {F<:Union{FlatS2,FlatS02}} = FieldTuple{basis(F)}(map(f->batch(f,D), f.fs))
batch(x::T, D::Int) where {T} = D==1 ? x : error("batching $T not implemented.")
batch(x, ::Nothing) = x

"""
    batch_promote!(to, f)

Promote `f` to the same batch size as `to` by replication. If both are already
the same batch size, no copy is made and `f` is returned. If promotion needs to
happen, the answer is stored in-place in `to` and returned.
"""
batch_promote!(to, f) = batchsize(to)==batchsize(f) ? f : (to .= f)

"""
    unbatch(f::FlatField)

If `f` is a batch-length-`D` field, return length-`D` vector of each batch
component, otherwise just return `f`. For the inverse operation, see
[`batch`](@ref).
"""
unbatch(f::FlatField{<:Flat{<:Any,<:Any,<:Any,1}}) = f
unbatch(f::FlatField{<:Flat{<:Any,<:Any,<:Any,D}}) where {D} = [batchindex(f,i) for i=1:D]

"""
    batchindex(f::FlatField, I)

Get the `I`th indexed batch (`I` can be a slice).
"""
batchindex(f::F, I) where {N,θ,∂mode,P<:Flat{N,θ,∂mode},F<:FlatS0{P}} =
    basetype(F){Flat{N,θ,∂mode,length(I)}}(f[:,:,I])
batchindex(f::FlatField, I) =
    FieldTuple{basis(f)}(map(f->batchindex(f, I), f.fs))

"""
    batchsize(f::FlatField)

The number of batches of in this object.
"""
batchsize(::FlatField{<:Flat{<:Any,<:Any,<:Any,D}}) where {D} = D
batchsize(x) = 1

struct BatchedVals{T,D,V<:AbstractVector{T}} <: Real
    vals :: V
    BatchedVals(v::V) where {T,V<:AbstractVector{T}} = new{T,length(v),V}(v)
end

"""
    BatchedReal(::Vector{<:Real}) <: Real

Holds a vector of real numbers and broadcasts algebraic operations over them,
as well as broadcasting with batched `FlatField`s, but is itself a `Real`.
"""
const BatchedReal{D,V,T<:Real} = BatchedVals{T,D,V}

batch(r::Real) = r
batch(r::Real, D::Int) = batch(collect(repeated(r,D)))
batch(v::AbstractVector) = BatchedVals(v)
batch(rs::Real...) = BatchedVals(collect(rs))
batchindex(br::BatchedVals, I) = getindex(br.vals,I)
batchsize(::BatchedVals{<:Any,D}) where {D} = D
struct BatchedRealStyle{D} <: AbstractArrayStyle{0} end
BroadcastStyle(::Type{<:BatchedReal{D}}) where {D} = BatchedRealStyle{D}()
BroadcastStyle(::FlatS0Style{F,M}, ::BatchedRealStyle{D′}) where {D′,N,θ,∂m,D,M,F<:FlatS0{Flat{N,θ,∂m,D}}} =
    (D==1 || D′==1 || D==D′) ? FlatS0Style{basetype(F){Flat{N,θ,∂m,max(D,D′)}},M}() : Broadcast.Unknown
BroadcastStyle(::FieldTupleStyle{B,Names,FS}, S2::BatchedRealStyle) where {B,Names,FS} =
    FieldTupleStyle{B,Names,Tuple{map_tupleargs(S1->typeof(Broadcast.result_style(S1(),S2)), FS)...}}()
for op in [:+, :-, :*, :/, :<, :<=, :&, :|, :(==)]
    @eval begin
        ($op)(a::BatchedReal, b::BatchedReal) = batch(broadcast(($op), a.vals, b.vals))
        ($op)(a::BatchedReal, b::Real)        = batch(broadcast(($op), a.vals, b))
        ($op)(a::Real,        b::BatchedReal) = batch(broadcast(($op), a,      b.vals))
    end
end
for op in [:-, :sqrt, :one, :zero, :isfinite, :eps]
    @eval ($op)(br::BatchedReal) = batch(broadcast(($op),br.vals))
end
for op in [:any, :all]
    @eval ($op)(br::BatchedReal) = ($op)(br.vals)
end
eltype(::BatchedVals{T}) where {T} = T
broadcastable(::Type{<:FlatS0{<:Flat,T}}, br::BatchedReal) where {T} = reshape(T.(br.vals),1,1,length(br.vals))
unbatch(br::BatchedVals) = br.vals
unbatch(r::Real) = r
Base.show(io::IO, br::BatchedReal) = print(io, "Batched", br.vals)
(::Type{T})(br::BatchedReal) where {T<:Real} = batch(T.(br.vals))
convert(::Type{<:BatchedVals{T,N}}, v::Bool) where {T,N} = batch(T(v),N)
Base.hash(bv::BatchedVals, h::UInt) = hash(bv.vals,hash(typeof(bv),h))

batch(L::Diagonal{<:Any,<:FlatField}, D::Int) = Diagonal(batch(diag(L), D))

"""
    batchmap(f, args...)

map function `f` over `args`, unbatching them first if they are batched, and
then batching the result.
"""
function batchmap(f, args...)
    _unbatch(x) = batchsize(x)==1 ? Ref(x) : unbatch(x)
    batch(broadcast(f, map(_unbatch, args)...)...)
end


text:
import NbodyGradient: init_nbody

@testset "init_nbody" begin

    elements = "elements.txt"
    H = [3,1,1]
    t0 = 7257.93115525 - 7300.0
    init = ElementsIC(t0, H, elements)
    init.elements[:,3] .-= 7300.0
    x, v, jac_init = init_nbody(init)
    elements_big = big.(init.elements)
    t0big = big(t0)
    init_big = ElementsIC(t0big, H, elements_big)
    xbig, vbig, jac_init_big = init_nbody(init_big)

    # elements = readdlm("elements.txt",',')

    # elements[2:3,1] /= 100.0
    # n_body = 4
    n_body = 3

    # jac_init     = zeros(Float64,7*n_body,7*n_body)
    # jac_init_big = zeros(BigFloat,7*n_body,7*n_body)
    jac_init_num = zeros(BigFloat, 7 * n_body, 7 * n_body)
    # x,v = init_nbody(elements,t0,n_body,jac_init)
    # elements_big = big.(elements)
    # t0big = big(t0)
    # xbig,vbig = init_nbody(elements_big,t0big,n_body,jac_init_big)
    elements0 = copy(init.elements)
    # dq = big.([1e-10,1e-5,1e-6,1e-6,1e-6,1e-5,1e-5])
    # dq = big.([1e-10,1e-8,1e-8,1e-8,1e-8,1e-8,1e-8])
    dq = big.([1e-12,1e-10,1e-10,1e-10,1e-10,1e-10,1e-10])

    # dq = big.([1e-15,1e-15,1e-15,1e-15,1e-15,1e-15,1e-15])
    # t0big = big(t0)
    # Now, compute derivatives numerically:
    for j = 1:n_body
        for k = 1:7
            elementsbig = big.(elements0)
            dq0 = dq[k]; if j == 1 && k == 1 ; dq0 = big(1e-15); end
            elementsbig[j,k] += dq0
            initp = ElementsIC(t0big, H, elementsbig)
            xp, vp, _ = init_nbody(initp)
            elementsbig[j,k] -= 2dq0
            initm = ElementsIC(t0big, H, elementsbig)
            xm, vm, _ = init_nbody(initm)
            for l in 1:n_body, p in 1:3
                i1 = (l - 1) * 7 + p
                if k == 1
                    j1 = j * 7
                else
                    j1 = (j - 1) * 7 + k - 1
                end
                jac_init_num[i1,  j1] = (xp[p,l] - xm[p,l]) / dq0 * .5
                jac1 = jac_init[i1,j1]; jac2 = jac_init_num[i1,j1]
                # if abs(jac1-jac2) > 1e-4*abs(jac1+jac2) && abs(jac1+jac2) > 1e-14
                #  println(l," ",p," ",j," ",k," ",jac_init_num[i1,j1]," ",jac_init[i1,j1]," ",jac_init_num[i1,j1]/jac_init[i1,j1])
                # end
                jac_init_num[i1 + 3,j1] = (vp[p,l] - vm[p,l]) / dq0 * .5
                jac1 = jac_init[i1 + 3,j1]; jac2 = jac_init_num[i1 + 3,j1]
                # if abs(jac1-jac2) > 1e-4*abs(jac1+jac2) && abs(jac1+jac2) > 1e-14
                # println(l," ",p+3," ",j," ",k," ",jac_init_num[i1+3,j1]," ",jac_init[i1+3,j1]," ",jac_init_num[i1+3,j1]/jac_init[i1+3,j1])
                # end
            end
        end
        jac_init_num[j * 7,j * 7] = 1.0
    end

    # println("Maximum jac_init-jac_init_num: ",maximum(abs.(jac_init-jac_init_num)))
    # println("Maximum jac_init-jac_init_num: ",maximum(abs.(asinh.(jac_init)-asinh.(jac_init_num))))
    # println("Maximum jac_init-jac_init_big: ",maximum(abs.(asinh.(jac_init)-asinh.(jac_init_big))))
    # println("Maximum jac_init/jac_init_big-1: ",maximum(abs.(jac_init[jac_init_big .!= 0.0]./jac_init_big[jac_init_big .!= 0.0].-1)))
    # println("Maximum jac_init_big-jac_init_num: ",maximum(abs.(asinh.(jac_init_num)-asinh.(jac_init_big))))
    # println(convert(Array{Float64,2},jac_init_big-jac_init_num))
    # @test isapprox(jac_init_num,jac_init)
    @test isapprox(jac_init_num, jac_init;norm=maxabs)
    @test isapprox(jac_init, convert(Array{Float64,2}, jac_init_big);norm=maxabs)
    # end
    for i in 1:21, j in 1:21
        if jac_init_big[i,j] != 0.0
            dij = abs(jac_init[i,j] / jac_init_big[i,j] - 1)
            if dij > 1e-12
                jac_max = dij
                println(i, " ", j, " ", convert(Float64, jac_max), " ", jac_init[i,j], " ", convert(Float64, jac_init_big[i,j]), " ", convert(Float64, jac_init[i,j] - jac_init_big[i,j]))
            end
        end
    end
end


text:
# Convert full BATSRUS AMR output to VTK.

using FortranFiles, Batsrus, WriteVTK

##

# Number of cells per block in each direction.
# These values are set by the Config.pl script.
# Set 1 for ignored directions!
const nI, nJ, nK = 4, 4, 4

# Maximum number of ghost cells set by Config.pl script.
# Valid values are 0,1,2,3,4,5
const nG = 2

# Refinement ratios in the 3 dimensions. Either 1 or 2.
# The values are set by the Config.pl script.
const iRatio, jRatio, kRatio = min(2, nI), min(2, nJ), min(2, nK)

# Number of dimensions in which grid adaptation is done
const nDimAmr = iRatio + jRatio + kRatio - 3

# Number of children per node
const nChild = 2^nDimAmr

# Possible values for the status variable
const Unset_     = -100 # index for unset values (that are otherwise larger)
const Unused_      = -1 # unused block (not a leaf)
const Refine_      = -2 # parent block to be refined
const DontCoarsen  = -3 # block not to be coarsened
const Coarsen_     = -4 # child block to be coarsened
const Used_        =  1 # currently used block (leaf)
const RefineNew_   =  2 # child block to be refined
const Refined_     =  3 # refined child block
const CoarsenNew_  =  4 # parent block to be coarsened
const Coarsened_   =  5 # coarsened parent block

# Deepest AMR level relative to root nodes (limited by 32 bit integers)
const MaxLevel = 30

# Named indexes of iTree_IA
const Status_   =  1
const Level_    =  2 # grid level
const Proc_     =  3 # processor index
const Block_    =  4 # block index
const MinLevel_ =  5 # minimum level allowed
const MaxLevel_ =  6 # maximum level allowed
const Coord0_   =  6 # equal to Coord1_-1
const Coord1_   =  7 # coordinate of node in 1st dimension
const Coord2_   =  8 # coordinate of node in 2nd dimension
const Coord3_   =  9 # coordinate of node in 3rd dimension
const CoordLast_=  9 # Coord0_ + MaxDim (?)
const Parent_   = 10 # Parent_ must be equal to Child0_
const Child0_   = 10 #
const Child1_   = Child0_ + 1
const ChildLast_= Child0_ + nChild

#=
the status of the node (used, unused, to be refined, to be coarsened, etc.);
the current, the maximum allowed and minimum allowed AMR levels for this node;
the three integer coordinates with respect to the whole grid;
the index of the parent node (if any);
the indexes of the children nodes (if any);
the processor index where the block is stored for active nodes;
the local block index for active nodes.
=#

filetag = "3d__ful_2_t00000141_n00000142"

## Loading AMR tree
f = FortranFile(filetag*".tree")

nDim, nInfo, nNode = read(f, Int32, Int32, Int32)
iRatio_D = read(f, (Int32,nDim)) # Array of refinement ratios
nRoot_D = read(f, (Int32,nDim)) # The number of root nodes in all dimension
iTree_IA = read(f, (Int32,(nInfo,nNode)))

close(f)

# Get all the used blocks
blockused_ = [iNode for iNode in 1:nNode if iTree_IA[Status_,iNode] == 1]

## If you know how blocks are organized, you can get VTK!
data = readdata(filetag*".out")

# vtkMultiBlockDataSet
vtmfile = vtk_multiblock("vtm/my_vtm_file")

for (i,iblock) in enumerate(blockused_)
   x = data.x[64*(i-1)+1:64*(i-1)+4,:,:,1][:]
   y = data.x[64*(i-1)+1:4:64*(i-1)+13,:,:,2][:]
   z = data.x[64*(i-1)+1:16:64*(i-1)+49,:,:,3][:]
   p = data.w[64*(i-1)+1:64*(i-1)+64,:,:,8][:]

   vtkfile = vtk_grid(vtmfile, x, y, z)
   vtkfile["Pressure",VTKPointData()] = p
end

outfiles = vtk_save(vtmfile)

# vtkHierarchicalBoxDataSet


text:
module GalewskyShallowWaterThetaMethod

using Gridap
using GridapDistributed
using GridapPETSc
using GridapP4est
using PartitionedArrays
const PArrays=PartitionedArrays
using GridapGeosciences
using SparseMatricesCSR
using FileIO

gamg = """
       -ksp_type gmres -ksp_rtol 1.0e-06 -ksp_atol 0.0 -ksp_gmres_restart 100
       -ksp_monitor -pc_type gamg -pc_gamg_type agg
       -mg_levels_esteig_ksp_type gmres -mg_coarse_sub_pc_type lu
       -mg_coarse_sub_pc_factor_mat_ordering_type nd -pc_gamg_process_eq_limit 50
       -pc_gamg_square_graph 9 pc_gamg_agg_nsmooths 1
       """

mumps = """
        -ksp_type preonly -ksp_error_if_not_converged true
        -pc_type lu -pc_factor_mat_solver_type mumps
        -mat_mumps_icntl_1 4
        -mat_mumps_icntl_14 100
        -mat_mumps_icntl_28 1
        -mat_mumps_icntl_29 2
        -mat_mumps_cntl_3 1.0e-6
        """

options = """
          -snes_type newtonls
          -snes_linesearch_type basic
          -snes_linesearch_damping 1.0
          -snes_rtol 1.0e-8
          -snes_atol 0.0
          -snes_monitor
          -snes_converged_reason
          -mm_ksp_type cg
          -mm_ksp_monitor
          -mm_ksp_rtol 1.0e-14
          -mm_pc_type jacobi
          """

  options_pcasm = """
          -snes_type newtonls
          -snes_linesearch_type basic
          -snes_linesearch_damping 1.0
          -snes_rtol 1.0e-6
          -snes_atol 1.0e-9
          -snes_monitor
          -snes_converged_reason
          -ksp_type gmres
          -ksp_rtol 1.0e-4
          -ksp_atol 1.0e-14
          -ksp_gmres_restart 5000
          -ksp_monitor
          -pc_type asm
          -sub_ksp_type preonly
          -sub_pc_type lu
          -mm_ksp_type cg
          -mm_ksp_monitor
          -mm_ksp_rtol 1.0e-4
          -mm_pc_type jacobi
          """

    options_fas =
      """
      -snes_view
      -snes_monitor
      -snes_converged_reason
      -snes_type ngmres
      -npc_snes_type fas
      -npc_snes_fas_levels 3
      -npc_snes_max_it 1
      -npc_fas_levels_snes_type gs
      -npc_fas_levels_snes_gs_sweeps 3
      -npc_snes_fas_smoothup 6
      -npc_snes_fas_smoothdown 6
     """

  function main_galewsky(parts,ir,
                        np,numrefs,dt,τ,
                        write_solution,write_solution_freq,title,order,degree,
                        verbose,mumps_relaxation,nstep)

    function mysnessetup(snes)
      ksp      = Ref{GridapPETSc.PETSC.KSP}()
      pc       = Ref{GridapPETSc.PETSC.PC}()
      mumpsmat = Ref{GridapPETSc.PETSC.Mat}()
      @check_error_code GridapPETSc.PETSC.SNESSetFromOptions(snes[])
      @check_error_code GridapPETSc.PETSC.SNESGetKSP(snes[],ksp)
      @check_error_code GridapPETSc.PETSC.KSPView(ksp[],C_NULL)
      @check_error_code GridapPETSc.PETSC.KSPSetType(ksp[],GridapPETSc.PETSC.KSPPREONLY)
      @check_error_code GridapPETSc.PETSC.KSPGetPC(ksp[],pc)
      @check_error_code GridapPETSc.PETSC.PCSetType(pc[],GridapPETSc.PETSC.PCLU)
      @check_error_code GridapPETSc.PETSC.PCFactorSetMatSolverType(pc[],GridapPETSc.PETSC.MATSOLVERMUMPS)
      @check_error_code GridapPETSc.PETSC.PCFactorSetUpMatSolverType(pc[])
      @check_error_code GridapPETSc.PETSC.PCFactorGetMatrix(pc[],mumpsmat)
      @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[],  4, 4)
      # percentage increase in the estimated working space
      @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[],  14, mumps_relaxation)
      @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[], 28, 2)
      @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[], 29, 2)
      @check_error_code GridapPETSc.PETSC.MatMumpsSetCntl(mumpsmat[], 3, 1.0e-6)
    end

    function set_ksp_mm(ksp)
      @check_error_code GridapPETSc.PETSC.KSPSetOptionsPrefix(ksp[],"mm_")
      @check_error_code GridapPETSc.PETSC.KSPSetFromOptions(ksp[])
    end

    t = PArrays.PTimer(parts,verbose=true)
    PArrays.tic!(t,barrier=true)
    ngcells, ngdofs = galewsky(parts,numrefs,dt,τ,
                              write_solution,write_solution_freq,title,order,degree,
                              verbose,mumps_relaxation,mysnessetup,set_ksp_mm,nstep)
    PArrays.toc!(t,"Simulation")
    display(t)
    map_main(t.data) do data
      out = Dict{String,Any}()
      merge!(out,data)
      out["ngdofs"] = ngdofs
      out["ngcells"] = ngcells
      out["np"] = np
      out["order"] = order
      out["dt"] = dt
      out["τ"] = τ
      out["degree"] = degree
      out["ir"] = ir
      save("$title.bson",out)
    end
  end

  ########

  function main(;
    np::Integer,
    nr::Integer,
    numrefs::Integer,
    dt::Float64,
    τ::Float64=dt/4,
    write_solution=false,
    write_solution_freq=4,
    title::AbstractString,
    k::Integer=1,
    degree::Integer=4,
    verbose::Bool=true,
    mumps_relaxation=1000,
    nstep::Integer=Int(24*60^2*20/dt)) # 20 days

    numrefs>=1 || throw(ArgumentError("numrefs should be larger or equal than 1"))

    prun(mpi,np) do parts
      for ir=1:nr
        GridapPETSc.with(args=split(options_pcasm)) do
          str_r   = lpad(ir,ceil(Int,log10(nr)),'0')
          title_r = "$(title)_ir$(str_r)"
          main_galewsky(parts,ir,np,numrefs,dt,τ,
              write_solution,write_solution_freq,title_r,
              k,degree,
              verbose,mumps_relaxation,nstep)
        end
      end
    end
  end

  include("../../../sequential/GalewskyInitialConditions.jl")

  function galewsky(parts,numrefs,dt,τ,
                    write_solution,write_solution_freq,title,order,degree,
                    verbose,mumps_relaxation,mysnessetup,set_ksp_mm,nstep)
      T      = dt*nstep
      θ      = 0.5
      model  = CubedSphereDiscreteModel(parts, numrefs; radius=rₑ)

      nls    = PETScNonlinearSolver()
      mmls   = PETScLinearSolver(set_ksp_mm)

      function ts()
        shallow_water_theta_method_full_newton_time_stepper(nls, model, order, degree,
        h₀, u₀, f, topography, g, θ, T, nstep, τ;
        mass_matrix_solver=mmls,
        am_i_root=PArrays.get_part_id(parts)==1,
        write_solution=write_solution,
        write_solution_freq=write_solution_freq,
        write_diagnostics=true,
        write_diagnostics_freq=1,
        dump_diagnostics_on_screen=true,
        output_dir=title)
      end

      if (PArrays.get_part_id(parts)==1)
        _,_,ndofs = ts()
      else
        _,_,ndofs = ts()
      end

      num_cells(model), ndofs
  end
end

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
Folds.collect(itr; kwargs...) = Folds.collect(itr, parallel_executor(itr; kwargs...))

function Folds.collect(itr, ex::Executor)
    acc =
        transduce(Map(SingletonVector), wheninit(collector, append!!), collector(), itr, ex)
    return reshape_as(finish!(unreduced(acc)), itr)
end

Folds.collect(itr, ex::PreferParallel) = Folds.collect(itr, executor_for(itr, ex))
Folds.collect(itr, ex::SequentialEx) =
    reshape_as(collect(extract_transducer(itr)...; ex.kwargs...), itr)
Folds.collect(itr, ex::ThreadedEx) =
    reshape_as(tcollect(extract_transducer(itr)...; ex.kwargs...), itr)
Folds.collect(itr, ex::DistributedEx) =
    reshape_as(dcollect(extract_transducer(itr)...; ex.kwargs...), itr)

reshape_as(ys, xs) = reshape_as(ys, xs, IteratorSize(xs))
reshape_as(ys, _, ::IteratorSize) = ys
reshape_as(ys, xs, ::HasShape) = reshape(ys, size(xs))
reshape_as(::Empty{T}, xs, isize::HasShape) where {T<:AbstractVector} =
    reshape_as(T(undef, length(xs)), xs, isize)

function as_copy_args(T, itr)
    xf, coll = extract_transducer(itr)
    return (xf, T, coll)
end

Folds.copy(T, itr; kwargs...) = Folds.copy(T, itr, parallel_executor(itr; kwargs...))

Folds.copy(T, itr, ex::SequentialEx) = copy(as_copy_args(T, itr)...; ex.kwargs...)

# TODO: implement `copy(itr, ex)`; move materializer to here

Folds.copy(T, itr, ex::Executor) = copy_default(T, itr, ex)

copy_default(T, itr, ex) =
    Folds.reduce(append!!, itr |> Map(SingletonVector), ex; init = Empty(T))

struct SingletonShimOf{T} end
(::SingletonShimOf{T})(x) where {T} = singletonshim(T, x)

function copy_default(T::Type{<:AbstractSet}, itr, ex)
    xf, array = extract_transducer(itr)
    if xf isa Union{Map{identity},IdentityTransducer}
        if array isa PartitionableArray
            basesize =
                something(get_basesize(ex), max(1, length(array) ÷ Threads.nthreads()))
            return Folds.reduce(
                union!!,
                Iterators.partition(array, basesize),
                set_basesize(ex, 1);
                init = Empty(T),
            )
        end
    end
    return Folds.reduce(union!!, itr |> Map(SingletonVector), ex; init = Empty(T))
end

function copy_default(T::Type{<:AbstractDict}, itr, ex)
    xf, array = extract_transducer(itr)
    if xf isa Union{Map{identity},IdentityTransducer}
        if array isa PartitionableArray
            basesize =
                something(get_basesize(ex), max(1, length(array) ÷ Threads.nthreads()))
            return Folds.reduce(
                merge!!,
                Iterators.partition(array, basesize),
                set_basesize(ex, 1);
                init = emptyshim(T),
            )
        end
    end
    return Folds.reduce(merge!!, itr |> Map(SingletonShimOf{T}()), ex; init = emptyshim(T))
end

Folds.map(f, itr; kwargs...) = Folds.collect(itr |> Map(f), parallel_executor(itr; kwargs...))
Folds.map(f, itr, ex::Executor) = Folds.collect(itr |> Map(f), ex)

map_check_no_kwargs(::NamedTuple{(),Tuple{}}) = nothing
@noinline function map_check_no_kwargs(kwargs)
    error(
        "`map(f, itrs..., executor)` does not accept any keyword arguments" *
        "\ngot:" *
        string(kwargs),
    )
end

function Folds.map(f, itr, itrs...; kwargs...)
    args0, ex0 = de_snoc(itr, itrs...)
    if ex0 isa Executor
        xs = zip(args0...)
        ex = ex0
        map_check_no_kwargs((; kwargs...))
    else
        xs = zip(itr, itrs...)
        ex = parallel_executor(bottom_foldable(xs); kwargs...)
    end
    return Folds.collect(xs |> MapSplat(f), ex)
end

# TODO: use `Folds.copy(AbstractSet, ...)`?
Folds.set(itr; kwargs...) = Folds.copy(Set, itr; kwargs...)
Folds.set(itr, ex; kwargs...) = Folds.copy(Set, itr, ex; kwargs...)

# TODO: use `Folds.copy(AbstractDict, ...)`?
Folds.dict(itr; kwargs...) = Folds.copy(Dict, itr; kwargs...)
Folds.dict(itr, ex; kwargs...) = Folds.copy(Dict, itr, ex; kwargs...)


text:
# ---
# title: 1305. All Elements in Two Binary Search Trees
# id: problem1305
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Sort, Tree
# link: <https://leetcode.com/problems/all-elements-in-two-binary-search-trees/description/>
# hidden: true
# ---
#
# Given two binary search trees `root1` and `root2`.
#
# Return a list containing _all the integers_ from _both trees_ sorted in
# **ascending** order.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2019/12/18/q2-e1.png)
#
#
#
#     Input: root1 = [2,1,4], root2 = [1,0,3]
#     Output: [0,1,1,2,3,4]
#
#
# **Example 2:**
#
#
#
#     Input: root1 = [0,-10,10], root2 = [5,1,7,0,2]
#     Output: [-10,0,0,1,2,5,7,10]
#
#
# **Example 3:**
#
#
#
#     Input: root1 = [], root2 = [5,1,7,0,2]
#     Output: [0,1,2,5,7]
#
#
# **Example 4:**
#
#
#
#     Input: root1 = [0,-10,10], root2 = []
#     Output: [-10,0,10]
#
#
# **Example 5:**
#
# ![](https://assets.leetcode.com/uploads/2019/12/18/q2-e5-.png)
#
#
#
#     Input: root1 = [1,null,8], root2 = [8,1]
#     Output: [1,1,8,8]
#
#
#
#
# **Constraints:**
#
#   * Each tree has at most `5000` nodes.
#   * Each node's value is between `[-10^5, 10^5]`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
mutable struct Population
  all::Dict{String, Int}
  size::Int
end

Population() = Population(Dict(), 0)
Base.keys(p::Population) = keys(p.all)

function binarize(p::Population, subpop)
  subpop_char_list = [string(Int(s)) for s in subpop]
  return join(subpop_char_list, "")
end

function label(p::Population, subpop)
  b = binarize(p, subpop)
  if !(b in keys(p))
    p.size += 1
    p.all[b] = p.size
  end
  return p.all[b]
end

# Example usage:
# p = Population()
# label(p, [1,0,1])  # 1
# label(p, [1,0,1])  # 1
# label(p, [1,1,1])  # 2
# label(p, [1,1,1,1]) # 3


text:
#!/usr/bin/env julia

Springpot =  RheoModelClass(
        # Model name
        name="springpot",
        # Model parameters,
        p = [:cᵦ, :β],
        # Relaxation modulus
        G = quote
              cᵦ*t^(-β)/gamma(1 - β)
            end,
        # Creep modulus
        J = quote
              (t^β)/(cᵦ*gamma(1 + β))
            end,
        # Storage modulus
        Gp = quote
                cᵦ*(ω^β)*cos(π*β/2)
             end,
        # Loss modulus
        Gpp = quote
                cᵦ*(ω^β)*sin(π*β/2)
              end,
        # Constraints
        constraint = quote
                (0<β<1)
                end,
        # Network
        info= "
                ____ ╱╲ ____
                     ╲╱  cᵦ, β
                "
        )

Spring =  RheoModelClass(
        # Model name
        name="spring",
        # Model parameters,
        p = [:k],
        # Relaxation modulus
        G = quote
              k
            end,
        # Creep modulus
        J = quote
              1/k
            end,
        # Storage modulus
        Gp = quote
                k
             end,
        # Loss modulus
        Gpp = quote
                0.0
              end,
        # Network
        info= "
                ___╱╲  ╱╲  ╱╲  ________
                     ╲╱  ╲╱  ╲╱  k
                "
        )

Dashpot =  RheoModelClass(
        # Model name
        name="dashpot",
        # Model parameters,
        p = [:η],
        # Creep modulus
        J = quote
              t/η
            end,
        # Storage modulus
        Gp = quote
                0.0
             end,
        # Loss modulus
        Gpp = quote
                η*ω
              end,
        # Network
        info= "
                 ___
             _____| |_____
                 _|_|
                     η
                "
        )


text:
@testset "Bouncing ball" begin
    H = readsxmodel("../examples/bball/bball_flattened.xml")
    @test nmodes(H) == 1
    @test ntransitions(H) == 1

    flow, invariant = H.modes[1][1], H.modes[1][2]
    @test flow == [:(x' = v), :(v' = -1)]
    @test invariant == [:(x >= 0)]

    assignment, guard = H.resetmaps[1][1], H.resetmaps[1][2]
    @test assignment == [:(x' = x), :(v' = -0.75v)]
    @test guard == [:(x = 0), :(v < 0)]
end

@testset "Van der Pol" begin
    H = readsxmodel("../examples/van_der_pol_fourth_quadrant/van_der_pol_fourth_quadrant.xml")
    @test nmodes(H) == 1
    @test ntransitions(H) == 0

    flow, invariant = H.modes[1][1], H.modes[1][2]
    @test flow == [:(x' = y), :(y' = -x + (1 - x ^ 2) * y)]
    @test invariant == [:(x <= 0), :(y >= 0)]
end

@testset "Lotka-Volterra" begin
    H = readsxmodel("../examples/lotka_volterra_fourth_quadrant/lotka_volterra_fourth_quadrant.xml")
    @test nmodes(H) == 1
    @test ntransitions(H) == 0

    flow, invariant = H.modes[1][1], H.modes[1][2]
    @test flow == [:(x' = x * (1 - y)), :(y' = -((1 - x) * y))]
    @test invariant == [:(x >= 0), :(y >= 0)]
end

@testset "Circle" begin
    H = readsxmodel("../examples/circle/circle_flattened.xml")
    @test nmodes(H) == 2
    @test ntransitions(H) == 2
    for (i, ti) in enumerate(transitions(H))
        if i == 1
            @test source(H, ti) == 1 && target(H, ti) == 2
        elseif i == 2
            @test source(H, ti) == 2 && target(H, ti) == 1
        end
    end

    m = H.modes[1]
    flow, invariant = m[1], m[2]
    @test flow == [:(x' = -y), :(y' = x)] && invariant == [:(y >= 0)]
    m = H.modes[2]
    flow, invariant = m[1], m[2]
    @test flow == [:(x' = -y), :(y' = x)] && invariant == [:(y <= 0)]

    # assignments in each mode 1, 2
    @test (H.resetmaps[1][1], H.resetmaps[2][1]) == ([:(x' = x), :(y' = y)], [:(x' = x), :(y' = y)])

    # guards in each mode 1, 2
    @test (H.resetmaps[1][2], H.resetmaps[2][2]) == ([:(y = 0)], [:(y >= 0)])
end


text:
__precompile__()

module ComputationalResources

export
    # Resources
    AbstractResource,
    AbstractCPU,
    ArrayFireLibs,
    CPU1,
    CPUThreads,
    CUDALibs,
    OpenCLLibs,
    # Settings
    TileSize,
    # functions
    addresource,
    haveresource,
    rmresource

### Resources

"""
    AbstractResource

The abstract supertype of all computational resources.
"""
abstract AbstractResource{T}

"""
    AbstractCPU

An abstract type indicating that computation should be performed using
the CPU.
"""
abstract AbstractCPU{T} <: AbstractResource{T}

"""
    CPU1()
    CPU1(settings)

Indicate that a computation should be performed using the CPU in single-threaded mode.
Optionally pass in an object specifying algorithmic settings.

# Examples:
```julia
filter(CPU1(), image, kernel)
filter(CPU1(TileSize(64,8)), image, kernel)
```
"""
immutable CPU1{T} <: AbstractCPU{T}
    settings::T
end
CPU1() = CPU1(nothing)
CPU1(r::AbstractResource) = CPU1(r.settings)

"""
    CPUThreads()
    CPUThreads(settings)

Indicate that a computation should be performed using the CPU in multi-threaded mode.
Optionally pass in an object specifying algorithmic settings.

# Examples:
```julia
filter(CPUThreads(), image, kernel)
filter(CPUThreads(TileSize(64,8)), image, kernel)
```
"""
immutable CPUThreads{T} <: AbstractCPU{T}
    settings::T
end
CPUThreads() = CPUThreads(nothing)
CPUThreads(r::AbstractResource) = CPUThreads(r.settings)

"""
    ArrayFireLibs()
    ArrayFireLibs(settings)

Indicate that computation should be performing using the ArrayFire libraries.
Optionally pass in an object specifying algorithmic settings.

# Examples:
```julia
filter(ArrayFireLibs(), image, kernel)
filter(ArrayFireLibs(backend), image, kernel)
```
"""
immutable ArrayFireLibs{T} <: AbstractResource{T}
    settings::T
end
ArrayFireLibs() = ArrayFireLibs(nothing)
ArrayFireLibs(r::AbstractResource) = ArrayFireLibs(r.settings)

"""
    CUDALibs()
    CUDALibs(settings)

Indicate that computation should be performing using the CUDA libraries.
Optionally pass in an object specifying algorithmic settings.

# Examples:
```julia
filter(CUDALibs(), image, kernel)
filter(CUDALibs(backend), image, kernel)
```
"""
immutable CUDALibs{T} <: AbstractResource{T}
    settings::T
end
CUDALibs() = CUDALibs(nothing)
CUDALibs(r::AbstractResource) = CUDALibs(r.settings)

"""
    OpenCLLibs()
    OpenCLLibs(settings)

Indicate that computation should be performing using the OpenCL libraries.
Optionally pass in an object specifying algorithmic settings.

# Examples:
```julia
filter(OpenCLLibs(), image, kernel)
filter(OpenCLLibs(backend), image, kernel)
```
"""
immutable OpenCLLibs{T} <: AbstractResource{T}
    settings::T
end
OpenCLLibs() = OpenCLLibs(nothing)
OpenCLLibs(r::AbstractResource) = OpenCLLibs(r.settings)

### Settings

"""
    TileSize(dims)

Request that an array computation be performed using tiles (blocks) of size `dims`.
"""
immutable TileSize{N}
    dims::NTuple{N,Int}
end

# Hold the available resources, allows packages to control loading of code
const resources = Set{Type}([CPU1])

"""
    addresource(T)

Add `T` to the list of available resources. For example,
`addresource(OpenCLLibs)` would indicate that you have a GPU and the
OpenCL libraries installed.
"""
addresource{T<:AbstractResource}(::Type{T}) = push!(resources, T)

"""
    rmresource(T)

Remove `T` from the list of available resources. For example,
`rmresource(OpenCLLibs)` would indicate that any future package loads
should avoid loading their specializations for OpenCL.
"""
rmresource{T<:AbstractResource}(::Type{T}) = delete!(resources, T)

"""
    haveresource(T)

Returns `true` if `T` is an available resource. For example,
`haveresource(OpenCLLibs)` tests whether the `OpenCLLibs` have been
added as an available resource. This function is typically used inside
a module's `__init__` function.

# Example:

```julia
# The __init__ function for MyPackage:
function __init__()
    ...  # other initialization code, possibly setting the LOAD_PATH
    if haveresource(OpenCLLibs)
        @eval using MyPackageCL  # a separate module containing OpenCL implementations
    end
    # Put additional resource checks here
    ...  # other initialization code, possibly cleaning up the LOAD_PATH
end
```
"""
haveresource{T<:AbstractResource}(::Type{T}) = T ∈ resources

"""
ComputationalResources makes it possible to dispatch to different methods that employ different computational resources. The exported resources are:

- `CPU1` (single-threaded computation)
- `CPUThreads` (multi-threaded computation)
- `ArrayFireLibs` (using the [ArrayFire package](https://github.com/JuliaComputing/ArrayFire.jl)
- `CUDALibs` (GPU computation using NVIDIA's CUDA libraries)
- `OpenCLLibs` (GPU computation using the OpenCL libraries)

There are also functions that interact with package initialization
machinery to control the availability of supported methods:

- `addresource(T)`: request newly-loaded packages to support resource `T`
- `rmresource(T)`: stop asking newly-loaded packages to support resource `T`
- `haveresource(T)`: test whether resource `T` has been requested
"""
ComputationalResources

end # module


text:
using Remez
using Test

import Remez: ratfn_leastsquares, ratfn_eval, ratfn_minimax, goldensection, winnow_extrema

# Test Gaussian elimination.
@testset "Gaussian test 1" begin
    m = BigFloat[1 1 2; 3 5 8; 13 34 21]
    v = BigFloat[1, -1, 2]
    ret = m \ v

    @test ret[1] ≈ big(109)/big(26)
    @test ret[2] ≈ big(-105)/big(130)
    @test ret[3] ≈ big(-31)/big(26)
end

# Test leastsquares rational functions.
@testset "Leastsquares test 1" begin
    n = 10000
    a = Array{BigFloat}(undef, n+1)
    for i = 0:1:n
        a[1+i] = i/BigFloat(n)
    end
    (nc, dc) = ratfn_leastsquares(exp, a, 2, 2)

    for x in a
        @test isapprox(exp(x), ratfn_eval(nc, dc, x); rtol=0,atol=1e-4)
    end
end

# Test golden section search.
@testset "Golden section test 1" begin
    x, y = goldensection(sin, big"0.0", big"0.1", big"4.0")
    @test isapprox(x, asin(BigFloat(1)))
    @test isapprox(y, 1)
end

# Test extrema-winnowing algorithm.
@testset "Winnow test 1" begin
    extrema = [(x, sin(20*x)*sin(197*x))  for x in big"0.0":big"0.001":big"1.0"]
    winnowed = winnow_extrema(extrema, 12)
    prevx, prevy = -1, 0
    for (x,y) in winnowed
        @test x > prevx
        @test y != 0
        @test prevy * y <= 0 # tolerates initial prevx having no sign
        @test abs(y) > 0.9
        prevx, prevy = x, y
    end
end

# Test actual minimax approximation.
@testset "Minimax test 1 (polynomial)" begin
    (nc, dc, e, x) = ratfn_minimax(exp, (0, 1), 4, 0)
    @test 0 < e < 1e-3
    for x = big"0.0":big"0.001":big"1.0"
        @test abs(ratfn_eval(nc, dc, x) - exp(x)) <= e * 1.0000001
    end
end
@testset "Minimax test 2 (rational)" begin
    (nc, dc, e, x) = ratfn_minimax(exp, (0, 1), 2, 2)
    @test 0 < e < 1e-3
    for x = big"0.0":big"0.001":big"1.0"
        @test abs(ratfn_eval(nc, dc, x) - exp(x)) <= e * 1.0000001
    end
end
@testset "Minimax test 3 (polynomial, weighted)" begin
    (nc, dc, e, x) = ratfn_minimax(exp, (0, 1), 4, 0,
                                   (x,y)->1/y)
    @test 0 < e < 1e-3
    for x = big"0.0":big"0.001":big"1.0"
        @test abs(ratfn_eval(nc, dc, x) - exp(x))/exp(x) <= e * 1.0000001
    end
end
@testset "Minimax test 4 (rational, weighted)" begin
    (nc, dc, e, x) = ratfn_minimax(exp, (0, 1), 2, 2,
                                   (x,y)->1/y)
    @test 0 < e < 1e-3
    for x = big"0.0":big"0.001":big"1.0"
        @test abs(ratfn_eval(nc, dc, x) - exp(x))/exp(x) <= e * 1.0000001
    end
end

@testset "Minimax test 5 (rational, weighted, odd degree)" begin

    (nc, dc, e, x) = ratfn_minimax(exp, (0, 1), 2, 1,
                                   (x,y)->1/y)
    @test 0 < e < 1e-3
    for x = big"0.0":big"0.001":big"1.0"
        @test abs(ratfn_eval(nc, dc, x) - exp(x))/exp(x) <= e * 1.0000001
    end
end


text:
const Array34{T} = Union{AbstractArray{T,3}, AbstractArray{T,4}}

function interp!(ca::AbstractArray{TO,3}, img1::AbstractArray{T,3}, img2::AbstractArray{T,3}, coefs::Vector) where {TO,T}
    for i = 1:size(ca,3)
         interp!(view(ca,:,:,i), view(img1,:,:,i), view(img2,:,:,i), coefs[i])
    end
    return ca
end

function interp!(ca::AbstractArray{TO,2}, img1::AbstractArray{T,2}, img2::AbstractArray{T,2}, coef) where {TO,T}
    for i in eachindex(ca)
        ca[i] = img1[i]*coef + img2[i]*(1.0-coef)
    end
    return ca
end

function check01(coef)
    if coef<0.0 || coef >1.0
        error("All coefficients must lie in the interval (0.0,1.0)")
    end
end


text:
using jInv.LinearSolvers
using Test
using ConvDiffMIPDECO
using jInv.ForwardShare
using jInv.Mesh
using MUMPSjInv

domain = [0 1. 0 1]

N = [16;32;64;128]
v = randn(2)
v = v/norm(v)

u  = X -> cos.(pi*X[:,1]).*cos.(pi*X[:,2])
ux = X -> -pi*sin.(pi*X[:,1]).*cos.(pi*X[:,2])
uy = X -> -pi*cos.(pi*X[:,1]).*sin.(pi*X[:,2])
sig = .01;

f  = X ->  2*sig*pi^2*u(X)+v[1]*ux(X)+v[2]*uy(X)

err = zeros(length(N))

times = zeros(length(N),3)

for j=1:2
	for k=1:length(N)
		Ainvs = (getJuliaSolver(), getMUMPSsolver())

		nk      = [N[k];N[k]]
		Mk  = getRegularMesh(domain,nk)
		xck = getCellCenteredGrid(Mk)
		xnk = getNodalGrid(Mk)
		pFork = getConvDiffFEMParam(Mk,sig=sig,v=v,bc=(:dir,:neu,:neu,:neu),gd=u,Ainv=Ainvs[j])

		dobsk,pFork = getData(f(xck), pFork)
		# pFor.Fields -= mean(pFor.Fields)
		utruek = u(xnk)
		# utrue -= mean(utrue)
		times[k,j] = @elapsed begin
			err[k] = norm(pFork.Fields - utruek,Inf)/norm(utruek,Inf)
		end
		println(err[k])
end
@test all(abs.(diff(log2.(err))).>1.2)
end

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO