text:
## Utilities

iv2r(iv, fs) = round(Int, minimum(iv)*fs) : round(Int, maximum(iv)*fs)

islo(trace, idx, lohi) = trace[idx] < lohi[1]
ishi(trace, idx, lohi) = trace[idx] > lohi[2]

function advance_to(isval::Function, trace, idxstart::Integer, idxstop::Integer, lohi)
    i = idxstart
    while i < idxstop && !isval(trace, i, lohi)
        i += 1
    end
    return i
end
advance_to(isval::Function, trace, startstop::AbstractUnitRange{Int}, lohi) =
    advance_to(isval, trace, first(startstop), last(startstop), lohi)

advance_to_lo(args...) = advance_to(islo, args...)
advance_to_hi(args...) = advance_to(ishi, args...)

function carrierfreq(lick, fs; twindow=4s, fwindow=2)
    # Check for the lickometer carrier signal
    sg = Periodograms.spectrogram(lick, nextpow(2, convert(AbstractFloat, twindow*fs)); fs=fs/(1Hz))
    pwr = vec(sum(sg.power; dims=2))[begin+1:end]   # drop the zero-frequency bin
    frq = sg.freq[begin+1:end]
    idx = argmax(pwr)
    pwrc = copy(pwr)
    pwrc[max(idx-fwindow, begin):min(idx+fwindow, end)] .= 0
    return frq[idx] * Hz, pwr[idx] / maximum(pwrc)
end

function checkcarrier(lick, fs, fcarrierlick=60Hz; frtol=0.05, powthresh=5)
    fpeak, prat = carrierfreq(lick, fs)
    abs(fpeak - fcarrierlick) < frtol * fcarrierlick || error("peak frequency $fpeak does not match expected $fcarrierlick")
    prat > powthresh || error("power ratio $prat does not meet SNR requirements")
    return nothing
end

function smooth(s, σ)
    σ2 = σ/sqrt(2)    # since we're using filtfilt
    f = [1/(sqrt(2π)*σ2) * exp(-i^2/(2σ2^2)) for i = floor(Int, -3σ2):ceil(Int, 3σ2)]
    return filtfilt(f, s)
end
smooth(s, σ::Quantity; fs::Quantity) = smooth(s, convert(AbstractFloat, σ*fs))

function slope(y)
    dydx = similar(y)
    for i in firstindex(y)+1:lastindex(y)-1
        dydx[i] = (y[i+1] - y[i-1])/2
    end
    dydx[begin] = y[begin+1] - y[begin]
    dydx[end] = y[end] - y[end-1]
    return dydx
end

function getindex_fraction(sortedlist, frac::Real)
    n = round(Int, (length(sortedlist)-1)*frac)
    return sortedlist[firstindex(sortedlist) + n]
end


text:
"""
  HttpErrors

Http-themed errors. Very sexy.

At present this module doesn't implement all the errors indicated in
https://tools.ietf.org/html/rfc7231. That's because I haven't needed
them yet. You're encouraged to submit a pull request.
"""
module HttpErrors

export Content,
  HttpError,
  bad_request,
  unauthorized,
  forbidden,
  not_found,
  method_not_allowed,
  conflict,
  imateapo,
  unprocessable_entity,
  internal_server_error,
  not_implemented

# HTTPError

"""
    Content

`HttpError` content-type. Alias for `Union{String, Exception, Nothing`.
"""
const Content = Union{String, Exception, Nothing}

"""
    HttpError

Http-themed error.

While nothing's stopping you from building one of these yourself, the value proposition of
this module comes from using the constructors (e.g. `bad_request`, `not_found`, &c.).

Instances come with fields `status::Int`, `name::String`, and `content::Content`.
"""
struct HttpError <: Exception
  status::Int
  name::String
  content::String
end

"""
    HttpError(status::Int, name::String, ::Nothing)

Produce an `HttpError` with no content.
"""
HttpError(status::Int, name::String, ::Nothing) = HttpError(status, name, "")

"""
    HttpError(status::Int, name::String, err::Exception)

Produce an `HttpError` with a serialized `err` as its content.
"""
HttpError(status::Int, name::String, err::Exception) =
  HttpError(status, name, sprint(showerror, err))

Base.showerror(io::IO, e::HttpError) = string(
  "$(e.status) $(e.name)",
  isempty(e.content) ? "" : " - $(e.content)"
) |> x -> write(io, x)

# bad_request

"""
    bad_request(content=nothing::Content)

400

The server cannot or will not process the request due to something that is perceived to
be a client error (e.g., malformed request syntax, invalid request content framing, or
deceptive request routing).
"""
bad_request(content=nothing::Content) =
  HttpError(400, "Bad Request", content) |> throw

# unauthorized

"""
  unauthorized(content=nothing::Content)

401

The request has not been applied because it lacks valid authentication credentials
for the target resource.
"""
unauthorized(content=nothing::Content) =
  HttpError(401, "Unauthorized", content) |> throw

# forbidden

"""
  forbidden(content=nothing::Content)

403

The server understood the request but refuses to authorize it.
"""
forbidden(content=nothing::Content) =
  HttpError(403, "Forbidden", content) |> throw

# not_found

"""
  not_found(content=nothing::Content)

404

The origin server did not find a current representation for the target resource or is
not willing to disclose that one exists.
"""
not_found(content=nothing::Content) =
  HttpError(404, "Not Found", content) |> throw

# method_not_allowed

"""
  method_not_allowed(content=nothing::Content)

405

The method received in the request-line is known by the origin server but not supported
by the target resource.
"""
method_not_allowed(content=nothing::Content) =
  HttpError(405, "Method Not Allowed", content) |> throw

# conflict

"""
  conflict(content=nothing::Content)

409

The request could not be completed due to a conflict with the current state of the target
resource. This code is used in situations where the user might be able to resolve the
conflict and resubmit the request.
"""
conflict(content=nothing::Content) =
  HttpError(409, "Conflict", content) |> throw

# im_a_teapot

"""
  im_a_teapot(content=nothing::Content)

418

Any attempt to brew coffee with a teapot should result in the error code "418 I'm a
teapot". The resulting entity body MAY be short and stout.
"""
imateapot(content=nothing::Content) =
  HttpError(418, "I'm A Teapot", content) |> throw

# unprocessable_entity

"""
  unprocessable_entity(content=nothing::Content)

422

The server understands the content type of the request entity (hence a 415
Unsupported Media Type status code is inappropriate), and the syntax of the request
entity is correct (thus a 400 Bad Request status code is inappropriate) but was
unable to process the contained instructions.
"""
unprocessable_entity(content=nothing::Content) =
  HttpError(422, "Unprocessable Entity", content) |> throw

# internal_server_error

"""
  internal_server_error(content=nothing::Content)

500

The server encountered an unexpected condition that prevented it from fulfilling the
request.
"""
internal_server_error(content=nothing::Content) =
  HttpError(500, "Internal Server Error", content) |> throw

# not_implemented

"""
  not_implemented(content=nothing::Content)

501

The server does not support the functionality required to fulfill the request.
"""
not_implemented(content=nothing::Content) =
  HttpError(501, "Not Implemented", content) |> throw

end

ext: jl


text:
# This file should contain site-specific commands to be executed on Julia startup
# Users should store their own personal commands in homedir(), in a file named .juliarc.jl


text:
#===============================================================================
 Robust NonMonotone Line search
 From DF-Sane paper https://www.ams.org/journals/mcom/2006-75-255/S0025-5718-06-01840-0/S0025-5718-06-01840-0.pdf
 This line search doesn't fit the rest of the API so I need to come up with an
 interface for nonmonotone line searches (it has to accept an fbar (or Q in Hager
 Zhang's notion) and a forcing term)
===============================================================================#
struct RNMS end

function find_steplength(rnms::RNMS, φ, fbar, ηk::T, τmin, τmax) where T
  φ0 = φ.φ0
  α₊ = T(1)
  α₋ = -α₊
  γ = T(1)/4
  for k = 1:100
    φα₊ = φ(α₊)
    if φα₊ ≤ fbar + ηk - γ*α₊^2*φ0
      return α₊, φα₊
    end
    φα₋ = φ(α₋)
    if φα₋ ≤ fbar + ηk - γ*α₋^2*φ0
      return α₋, φα₋
    end
    αt = α₊^2*φ0/(φα₊ + (2*α₊-1)*φ0)
    if αt < τmin*α₊
      α₊ = τmin*α₊
    elseif αt > τmax*α₊
      α₊ = τmax*α₊
    else
      α₊ = αt
    end

    # This needs to be safe guarded!
    αt = α₋^2*φ0/(φα₋ + (2*α₋-1)*φ0)
    if αt < τmin*α₋
      α₋ = τmin*α₋
    elseif αt > τmax*α₋
      α₋ = τmax*α₋
    else
      α₋ = αt
    end
  end
  throw(ErrorException("Line search failed."))
end

function safeguard_σ(σ::T, σmin, σmax, F) where T
  if abs(σ) < σmin || abs(σ) > σmax
    normF = norm(F)
    if normF > T(1)
      σ = T(1)
    elseif T(1)/10^5 ≤ normF ≤ T(1)
      σ = inv(normF)
    elseif normF < T(1)/10^5
      σ = T(10)^5
    end
  end
  return σ
end

struct DFSANE end
function nlsolve(F::OnceDiffed, x0, scheme::DFSANE)
  T = eltype(x0)

  σmin, σmax = 1e-10, 1e10
  τmin, τmax = T(1)/10, T(1)/2

  nexp = 2

  fvals = @SVector(T[])

  x = copy(x0)
  Fx = F(x)
  y = copy(Fx)
  normFx0 = norm(Fx)
  fx = normFx0^nexp
  fvals = push(fvals, fx)
  σ = T(1)
  σ = safeguard_σ(σ, σmin, σmax, Fx)
  for j = 1:10^5
    fbar = maximum(x->fvals[x], max(1, length(fvals)-5+1):length(fvals))
    d = -σ*Fx
    ηk = normFx0/(1+j)^2
    φ = LineObjective(OnceDiffed((_n,x)->norm(F(x))^nexp), x, d, fx, 0*fx)
    α, φα = find_steplength(RNMS(), φ, fbar, ηk, τmin, τmax)
    # LineObjective(x) should return z as well!
    s = α*d
    x .= x.+s
    y .= Fx
    Fx = F(x)
    y .-= Fx
    fx = norm(Fx)^nexp
    fvals = push(fvals, fx)
    σ = norm(s)^2/dot(s, y)
    σ = safeguard_σ(σ, σmin, σmax, Fx)
    if sqrt(φα) < 1e-8
      break
    end
  end
  x, Fx
end


text:
function ABCDQR_em1(y,u,nx::Int;max_iter::Int=100,tol::Float64=1e-6,txo::Bool=true)
	#=
	estimate A, C, Q, R, m1, P1 using the EM algorithm
	and the subspace algorithm for the starting point
	for model

	x_{t+1} = A*x_{t} + B*u_{t} + w_{t}
	y_{t}   = C*x_{t} + D*u_{t} + v_{t}

	cov(w_{t},v_{t}) = [Q 0;0 R]
	x1 -> N(m1,P1)

	javier.cara@upm.es, 2017-01

	ss=true : estationary algorithm

	=#

	# data as a matrix and by rows
	y,ny,nt = byrow(y)
	u,nu,aux = byrow(u)

	# starting values from the ssi algorithm
	i = nx+1
	Ai,Ci,Qi,Ri,Si = ACQRS_sub(y,nx,i)
	Bi = zeros(nx,nu)
	Di = zeros(ny,nu)

	# x11 = C^{-1}(y_1 - v_1)
	m1i = Ci\y[:,1:1]
	P1i = zeros(nx,nx)

	A,B,C,D,Q,R,m1,P1,loglikv,aic = ABCDQR_em(y,u,Ai,Bi,Ci,Di,Qi,Ri,m1i,P1i,max_iter,tol,txo)

	# output
	return ABCDQR(A,B,C,D,Q,R,m1,P1,loglikv[end],aic)

end


text:
using Artifacts
using LazyArtifacts
import JLD2
using ElectronicStructure: MolecularData

## TODO: encapsulate this somehow
"""
    get_molecular_data()

Return a `Dict` containing a few examples of PySCF electronic structure
calculations. They are retrieved in `JLD2` format using the `Artifacts`
framework. The `JLD2` file is read and the result is returned.

The data returned was generated using the function `generate_pyscf_molecules`
in *compare_formats.jl*.
"""
function get_molecular_data()
    rootpath = artifact"molecular_data"
    filename = joinpath(rootpath, "molecular_data_examples.jld2")
    data = JLD2.load(filename)
    return data
end

@testset "molecular data sets" begin
    data_top = get_molecular_data()
    @test collect(keys(data_top)) == ["header", "data"]
    data = data_top["data"]
    for mol_name in keys(data)
        mol1 = data[mol_name]
        @test isa(mol1, MolecularData)
        @test find_index_order(mol1.two_body_integrals) == :chemist
        iop = InteractionOperator(mol1)
        @test find_index_order(iop.two_body_tensor) == :physicist
        @test find_index_order(iop) == :physicist
        iop2 = InteractionOperator(mol1; index_order=:chemist)
        @test find_index_order(iop2.two_body_tensor) == :chemist
        iop2 = InteractionOperator(mol1; index_order=:intermediate)
        @test find_index_order(iop2.two_body_tensor) == :intermediate
        for spin_order in (:block, :alternating, :none)
            @test find_index_order(InteractionOperator(mol1; spin_order=spin_order).two_body_tensor) == :physicist
        end
    end
end


text:
using ITensors, Test

@testset "basic DMRG" begin
  @testset "Spin-one Heisenberg" begin
    N = 100
    sites = spinOneSites(N)

    ampo = AutoMPO(sites)
    for j=1:N-1
      add!(ampo,"Sz",j,"Sz",j+1)
      add!(ampo,0.5,"S+",j,"S-",j+1)
      add!(ampo,0.5,"S-",j,"S+",j+1)
    end
    H = toMPO(ampo)

    psi = randomMPS(sites)

    sweeps = Sweeps(5)
    @test length(sweeps) == 5
    maxdim!(sweeps,10,20,100,100)
    mindim!(sweeps,1,10,20,20)
    cutoff!(sweeps,1E-11)
    str = split(sprint(show, sweeps), '\n')
    @test str[1] == "Sweeps"
    @test str[2] == "1 cutoff=1.0E-11, maxdim=10, mindim=1"
    @test str[3] == "2 cutoff=1.0E-11, maxdim=20, mindim=10"
    @test str[4] == "3 cutoff=1.0E-11, maxdim=100, mindim=20"
    @test str[5] == "4 cutoff=1.0E-11, maxdim=100, mindim=20"
    @test str[6] == "5 cutoff=1.0E-11, maxdim=100, mindim=20"
    energy,psi = dmrg(H,psi,sweeps,maxiter=2,quiet=true)
    @test energy ≈ -138.94 rtol=1e-3
    # test with SVD too!
    psi = randomMPS(sites)
    energy,psi = dmrg(H,psi,sweeps,maxiter=2,quiet=true,which_factorization="svd")
    @test energy ≈ -138.94 rtol=1e-3
  end

  @testset "Transverse field Ising" begin
    N = 32
    sites = spinHalfSites(N)
    ψ0 = randomMPS(sites)

    ampo = AutoMPO(sites)
    for j = 1:N
      j < N && add!(ampo,-1.0,"Sz",j,"Sz",j+1)
      add!(ampo,-0.5,"Sx",j)
    end
    H = toMPO(ampo)

    sweeps = Sweeps(3)
    maxdim!(sweeps,10,20)
    cutoff!(sweeps,1E-12)
    energy,ψ = dmrg(H,ψ0,sweeps,quiet=true)

    # Exact energy for transverse field Ising model
    # with open boundary conditions at criticality
    energy_exact = 0.25 - 0.25/sin(π/(2*(2*N+1)))
    energy ≈ energy_exact
  end
end


text:
################################################################################
#
# Routines for synthetic data generation
#
################################################################################

export gendata

### Wrapper for JUDI data types

function gendata(
	m::Model, wav::judiVector, rcv_geom::Geometry;
	opt = Options(),
	dt_comp::Union{Nothing, Array{Any, 1}} = nothing,
	save_wavefield::Bool = false)

	return gendata_rawinput(m, wav.geometry, rcv_geom, wav.data, 1:wav.nsrc; opt = opt, dt_comp = dt_comp, save_wavefield = save_wavefield)

end

### Objective wrapper for array input

function gendata(
	m::Array{R, 2}, n::NTuple{2, Int64}, d::NTuple{2, R}, o::NTuple{2, R},
	wav::judiVector, rcv_geom::Geometry;
	opt = Options(),
	dt_comp::Union{Nothing, Array{Any, 1}} = nothing,
	save_wavefield::Bool = false)

	return gendata(Model(n, d, o, m), wav, rcv_geom; opt = opt, dt_comp = dt_comp, save_wavefield = save_wavefield)

end

### Serial/parallel implementation

## Serial implementation

function gendata_rawinput(
	model_full::Model,
	src_geom::Geometry, rcv_geom::Geometry,
	src_data::Array{R, 2},;
	opt = Options(),
	dt_comp::Union{Nothing, R} = nothing,
	save_wavefield::Bool = false)

	# Setting pre-defined absorbing layer size
	model_full.nb = NB

	# Load full geometry for out-of-core geometry containers
	typeof(rcv_geom) == GeometryOOC && (rcv_geom = Geometry(rcv_geom))
	typeof(src_geom) == GeometryOOC && (src_geom = Geometry(src_geom))
	length(model_full.n) == 3 ? dims = [3, 2, 1] : dims = [2, 1] # model dimensions for Python are (z,y,x) and (z,x)

	# Limit model to area with sources/receivers
	if opt.limit_m == true
		model = deepcopy(model_full)
		model = limit_model_to_receiver_area(src_geom, rcv_geom, model, opt.buffer_size)
	else
		model = model_full
	end

	# Set up Python model structure
	model_py = devito_model(model, "F", 1, opt, 0)

	# Remove receivers outside the modeling domain (otherwise leads to segmentation faults)
	rcv_geom = remove_out_of_bounds_receivers(rcv_geom, model)

	# Call to objective with Julia/Devito interface function
	return gendata_jldevito(model_py, model.o, src_geom, rcv_geom, src_data; opt = opt, dt_comp = dt_comp, save_wavefield = save_wavefield)

end

## Parallel implementation

function gendata_rawinput(
	model::Model,
	src_geom::Geometry, rcv_geom::Geometry,
	src_data,
	src_idx::UnitRange{Int64};
	opt = Options(),
	dt_comp::Union{Nothing, Array{Any, 1}} = nothing,
	save_wavefield::Bool = false)

	# Setup parallelization
	p = default_worker_pool()
	time_modeling_par = remote(gendata_rawinput)
	time_modeling = retry(time_modeling_par)

	# Initialize output
	nsrc = length(src_idx)
	results = Array{Any}(undef, nsrc)

	# Process shots from source channel asynchronously
	@sync begin
		for j = 1:nsrc

			# Local geometry for current position
			opt_loc = subsample(opt, j)
			src_geom_loc = subsample(src_geom, j)
			rcv_geom_loc = subsample(rcv_geom, j)

			# Selecting variables for current shot index
			src_data_loc = src_data[j]
			if dt_comp == nothing
				dt_comp_loc = nothing
			else
				dt_comp_loc = dt_comp[j]
			end

			# Local result
			results[j] = @spawn gendata_rawinput(model, src_geom_loc, rcv_geom_loc, src_data_loc, opt = opt_loc; dt_comp = dt_comp_loc, save_wavefield = save_wavefield)

		end
	end

	# Aggregating results
	argout = Array{Array, 1}(undef, nsrc)
	for j = 1:nsrc
		argout[j] = fetch(results[j])
	end
	if save_wavefield
		return argout
	else
		return judiVector(rcv_geom, argout)
	end

end

### Julia/devito implementation

function gendata_jldevito(
	model_py::PyCall.PyObject,
	origin,
	src_geom::Geometry, rcv_geom::Geometry,
	src_data::Array{R, 2};
	opt = Options(),
	dt_comp::Union{Nothing, R} = nothing,
	save_wavefield::Bool = false)

	# Loading python modules for devito implementation of objTWRIdual
	devitomod = load_codegen_TWRIdual()

	# Interpolate input data to computational grid
	if dt_comp == nothing
		dt_comp = R(model_py.critical_dt)
	end
	q_in = time_resample(src_data, src_geom, dt_comp)[1]
	nt_comp = size(q_in, 1)
	nt_rcv = Int(trunc(rcv_geom.t[1]/dt_comp+1))

	# Set up coordinates with devito dimensions
	src_coords = setup_grid(src_geom, model_py.shape, origin)
	rcv_coords = setup_grid(rcv_geom, model_py.shape, origin)

	# Computing output
	if !save_wavefield
		argout = pycall(devitomod.forward,
						Tuple{Array{R, 2}, Nothing},
						model_py,
						PyReverseDims(copy(transpose(src_coords))), PyReverseDims(copy(transpose(rcv_coords))),
						PyReverseDims(copy(transpose(q_in))),
						dt_comp, SPACE_ORDER, save_wavefield)
		return argout[1]
	else
		argout = pycall(devitomod.forward,
						Tuple{Array{R, 2}, PyObject},
						model_py,
						PyReverseDims(copy(transpose(src_coords))), PyReverseDims(copy(transpose(rcv_coords))),
						PyReverseDims(copy(transpose(q_in))),
						dt_comp, SPACE_ORDER, save_wavefield)
		return argout[2].data[:, model_py.nbpml+1:end-model_py.nbpml, model_py.nbpml+1:end-model_py.nbpml]
	end

end


text:
using Test
import DifferentiableBackwardEuler
import DPFEHM
import Zygote

doplot = false
if doplot == true
	import PyPlot
end

#Theis solution
function W(u)
	if u <= 1
		return -log(u) + -0.57721566 + 0.99999193u^1 + -0.24991055u^2 + 0.05519968u^3 + -0.00976004u^4 + 0.00107857u^5
	else
		return (u^2 + 2.334733u + 0.250621) / (u^2 + 3.330657u + 1.681534) * exp(-u) / u
	end
end
function theis(t, r, T, S, Q)
	return Q * W(r^2 * S / (4 * T * t)) / (4 * pi * T)
end
#utitlity for checking gradients
function checkgradientquickly(f, x0, gradf, n; delta::Float64=1e-8, kwargs...)
	indicestocheck = sort(collect(1:length(x0)), by=i->abs(gradf[i]), rev=true)[1:n]
	#indicestocheck = [indicestocheck; rand(1:length(x0), n)]
	f0 = f(x0)
	for i in indicestocheck
		x = copy(x0)
		x[i] += delta
		fval = f(x)
		grad_f_i = (fval - f0) / delta
		@test isapprox(gradf[i], grad_f_i; kwargs...)
	end
end

#test Theis solution against groundwater model
steadyhead = 1e3
sidelength = 50.0
thickness = 10.0
n = 101
ns = [n, n]
coords, neighbors, areasoverlengths, volumes = DPFEHM.regulargrid2d([-sidelength, -sidelength], [sidelength, sidelength], ns, thickness)
k = 1e-5
Ks = fill(k, length(neighbors))
T = thickness * k
Q = 1e-3
Ss = 0.1
specificstorage = fill(Ss, size(coords, 2))
S = Ss * thickness
Qs = zeros(size(coords, 2))
Qs[ns[2] * (div(ns[1] + 1, 2) - 1) + div(ns[2] + 1, 2)] = -Q#put a fluid source in the middle
dirichletnodes = Int[]
dirichleths = zeros(size(coords, 2))
for i = 1:size(coords, 2)
	if abs(coords[1, i]) == sidelength || abs(coords[2, i]) == sidelength
		push!(dirichletnodes, i)
		dirichleths[i] = steadyhead
	end
end
isfreenode, nodei2freenodei, freenodei2nodei = DPFEHM.getfreenodes(prod(ns), dirichletnodes)
function unpack(p)
	@assert length(p) == length(neighbors)
	Ks = p[1:length(neighbors)]
	return Ks
end
function f_gw(u, p, t)
	Ks = unpack(p)
	return -DPFEHM.groundwater_residuals(u, Ks, neighbors, areasoverlengths, dirichletnodes, dirichleths, Qs, specificstorage, volumes)
end
function f_gw_u(u, p, t)
	Ks = unpack(p)
	return -DPFEHM.groundwater_h(u, Ks, neighbors, areasoverlengths, dirichletnodes, dirichleths, Qs, specificstorage, volumes)
end
function f_gw_p(u, p, t)
	Ks = unpack(p)
	return -DPFEHM.groundwater_Ks(u, Ks, neighbors, areasoverlengths, dirichletnodes, dirichleths, Qs, specificstorage, volumes)
end
f_gw_t(u, p, t) = zeros(length(u))
t0 = 0.0
tfinal = 60 * 60 * 24 * 1e1
h0 = fill(steadyhead, size(coords, 2) - length(dirichletnodes))
p = Ks
print("groundwater forward")
@time h_gw = DifferentiableBackwardEuler.steps_diffeq(h0, f_gw, f_gw_u, f_gw_p, f_gw_t, p, t0, tfinal; abstol=1e-6, reltol=1e-6)
dpfehm_solution = map(h->DPFEHM.addboundaryconditions(h, dirichletnodes, dirichleths, isfreenode, nodei2freenodei), h_gw.u)
r0 = 0.1
goodnodes = collect(filter(i->coords[2, i] == 0 && coords[1, i] > r0 && coords[1, i] <= sidelength / 2, 1:size(coords, 2)))
rs = coords[1, goodnodes]
theis_drawdowns = theis.(h_gw.t[end], rs, T, S, Q)
gw_drawdowns = -dpfehm_solution[end][goodnodes] .+ steadyhead
if doplot
	fig, ax = PyPlot.subplots()
	ax.plot(rs, theis_drawdowns, "r.", ms=20, label="Theis")
	ax.plot(rs, gw_drawdowns, "k", linewidth=3, label="DPFEHM groundwater")
	ax.set_xlabel("x [m]")
	ax.set_ylabel("drawdown [m]")
	ax.legend()
	display(fig)
	println()
	PyPlot.close(fig)
end
@test isapprox(theis_drawdowns, gw_drawdowns; atol=1e-1)
g_gw(p) = DifferentiableBackwardEuler.steps(h0, f_gw, f_gw_u, f_gw_p, f_gw_t, p, t0, tfinal; abstol=1e-6, reltol=1e-6)[nodei2freenodei[goodnodes[round(Int, 0.25 * end)]], end]
print("groundwater gradient")
@time grad_gw_zygote = Zygote.gradient(g_gw, p)[1]
checkgradientquickly(g_gw, p, grad_gw_zygote, 3; delta=1e-8, rtol=1e-1)

#test Theis solution against richards equation model
coords_richards = vcat(coords, zeros(size(coords, 2))')
alphas = fill(0.5, length(neighbors))
Ns = fill(1.25, length(neighbors))
function f_richards(u, p, t)
	Ks = unpack(p)
	return DPFEHM.richards_residuals(u, Ks, neighbors, areasoverlengths, dirichletnodes, dirichleths, coords_richards, alphas, Ns, Qs, specificstorage, volumes)
end
function f_richards_u(u, p, t)
	Ks = unpack(p)
	return DPFEHM.richards_psi(u, Ks, neighbors, areasoverlengths, dirichletnodes, dirichleths, coords_richards, alphas, Ns, Qs, specificstorage, volumes)
end
function f_richards_p(u, p, t)
	Ks = unpack(p)
	return DPFEHM.richards_Ks(u, Ks, neighbors, areasoverlengths, dirichletnodes, dirichleths, coords_richards, alphas, Ns, Qs, specificstorage, volumes)
end
f_richards_t(u, p, t) = zeros(length(u))
print("richards forward")
h0 = fill(steadyhead, size(coords, 2))
@time h_richards = DifferentiableBackwardEuler.steps_diffeq(h0, f_richards, f_richards_u, f_richards_p, f_richards_t, p, t0, tfinal; abstol=1e-6, reltol=1e-6)
richards_drawdowns = -h_richards[end][goodnodes] .+ steadyhead
if doplot
	fig, ax = PyPlot.subplots()
	ax.plot(rs, theis_drawdowns, "r.", ms=20, label="Theis")
	ax.plot(rs, richards_drawdowns, "k", linewidth=3, label="DPFEHM richards")
	ax.set_xlabel("x [m]")
	ax.set_ylabel("drawdown [m]")
	ax.legend()
	display(fig)
	println()
	PyPlot.close(fig)
end
@test isapprox(theis_drawdowns, richards_drawdowns; atol=1e-1)
#@test isapprox(h_richards[isfreenode, :], h_gw[:, :])#make sure richards and groundwater are giving the same thing -- commenting this out for now because the recent changes in gw mean they shouldn't be exactly the same -- uncomment when those changes get propagated into richards as well
g_richards(p) = DifferentiableBackwardEuler.steps(h0, f_richards, f_richards_u, f_richards_p, f_richards_t, p, t0, tfinal; abstol=1e-6, reltol=1e-6)[goodnodes[round(Int, 0.25 * end)], end]
print("richards gradient")
@time grad_richards_zygote = Zygote.gradient(g_richards, p)[1]
checkgradientquickly(g_richards, p, grad_richards_zygote, 3; delta=1e-8, rtol=1e-1)


text:
using Documenter, AbstractPlutoDingetjes

makedocs(;

sitename = "AbstractPlutoDingetjes",
)


text:
#=

Olympic puzzle in Julia ConstraintSolver.jl

Benchmark for Prolog (BProlog)
"""
  File   : olympic.pl
  Author : Neng-Fa ZHOU
  Date   : 1993

  Purpose: solve a puzzle taken from Olympic Arithmetic Contest

   Given ten variables with the following configuration:

               X7   X8   X9   X10

                  X4   X5   X6

                     X2   X3

                        X1

  We already know that X1 is equal to 3 and want to assign each variable
  with a different integer from {1,2,...,10} such that for any three
  variables
                      Xi   Xj

                         Xk
  the following constraint is satisfied:

                    |Xi-Xj| = Xk
"""

Model created by Hakan Kjellerstrand, hakank@gmail.com
See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
using Printf
const CS = ConstraintSolver
include("constraints_utils.jl")

function olympic(print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    n = 10
    @variable(model, 1 <= x[1:n] <= n, Int)
    x1,x2,x3,x4,x5,x6,x7,x8,x9,x10 = x
    @constraint(model, x in CS.AllDifferentSet())

    @constraint(model,x1 == 3)
    my_abs(model,x2,x3,x1)
    my_abs(model,x4,x5,x2)
    my_abs(model,x5,x6,x3)
    my_abs(model,x7,x8,x4)
    my_abs(model,x8,x9,x5)
    my_abs(model,x9,x10,x6)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                xx = convert.(Integer,JuMP.value.(x; result=sol))
                # println("x:$xx")
                @printf "%d  %d  %d %d\n"  xx[7] xx[8] xx[9] xx[10]
                @printf "  %d  %d  %d\n" xx[4] xx[5] xx[6]
                @printf "    %d  %d\n" xx[2] xx[3]
                @printf "      %d\n" xx[1]
                println()

            end
        end
    else
        println("status:$status")
    end

    return status
end

@time olympic()

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
#=
Main program for two phase flow inversion
=#

include("args.jl")
include("main_fwi_sepinv.jl")

function sw_p_to_lambda_den(sw, p)
    sw = tf.reshape(sw, (1, m, n, 1))
    p = tf.reshape(p, (1, m, n, 1))
    sw = tf.image.resize_bilinear(sw, (nz, nx))
    p = tf.image.resize_bilinear(p, (nz, nx))
    sw = cast(sw, Float64)
    p = cast(p, Float64)
    sw = squeeze(sw)
    p = squeeze(p)
    # tran_lambda, tran_den = Gassman(sw)
    # tran_lambda, tran_den = RockLinear(sw) # test linear relationship
    tran_lambda, tran_den = Patchy(sw)
    tran_lambda_pad =  tf.pad(tran_lambda, [nPml (nPml+nPad); nPml nPml], constant_values=3500.0^2*2200.0/3.0) /1e6
    tran_den_pad = tf.pad(tran_den, [nPml (nPml+nPad); nPml nPml], constant_values=2200.0)
    return tran_lambda_pad, tran_den_pad
end

# NOTE Generate Data
if args["generate_data"]
    println("Generate Test Data...")
    K = 20.0 .* ones(m,n) # millidarcy
    K[8:10,:] .= 120.0
    # K[17:21,:] .= 100.0
    # for i = 1:m
    #     for j = 1:n
    #         if i <= (14 - 24)/(30 - 1)*(j-1) + 24 && i >= (12 - 18)/(30 - 1)*(j-1) + 18
    #             K[i,j] = 100.0
    #         end
    #     end
    # end
    tfCtxTrue = tfCtxGen(m,n,h,NT,Δt,Z,X,ρw,ρo,μw,μo,K,g,ϕ,qw,qo, sw0, true)
    out_sw_true, out_p_true = imseq(tfCtxTrue)
    lambdas = Array{PyObject}(undef, n_survey)
    dens = Array{PyObject}(undef, n_survey)
    for i = 1:n_survey
        sw = out_sw_true[survey_indices[i]]
        p = out_p_true[survey_indices[i]]
        lambdas[i], dens[i] = sw_p_to_lambda_den(sw, p)
    end

    misfit = Array{PyObject}(undef, n_survey)
    for i = 1:n_survey
        if !isdir("./$(args["version"])/Data$i")
            mkdir("./$(args["version"])/Data$i")
        end
        para_fname = "./$(args["version"])/para_file$i.json"
        survey_fname = "./$(args["version"])/survey_file$i.json"
        paraGen(nz_pad, nx_pad, dz, dx, nSteps, dt, f0, nPml, nPad, para_fname, survey_fname, "./$(args["version"])/Data$i/")
        # shot_inds = collect(1:3:length(z_src)) .+ mod(i-1,3) # 5src rotation
        # shot_inds = i # 1src rotation
        shot_inds = collect(1:length(z_src)) # all sources
        surveyGen(z_src[shot_inds], x_src[shot_inds], z_rec, x_rec, survey_fname)
        tf_shot_ids0 = constant(collect(0:length(shot_inds)-1), dtype=Int32)
        misfit[i] = fwi_obs_op(lambdas[i], tf_shear_pad, dens[i], tf_stf, tf_gpu_id0, tf_shot_ids0, para_fname)
    end
    config = tf.ConfigProto()
    config.intra_op_parallelism_threads = 24
    config.inter_op_parallelism_threads = 24
    sess = Session(config=config); init(sess);
    run(sess, misfit)
    error("Generate Data: Stop")
end

for iSur = 2:n_survey
    fwi_sep(iSur)
end


text:
using Plots, Test, Dates

@testset "Limits" begin
    y=[1.0*i*i for i in 1:10]
    x=[Date(2019,11,i) for i in 1:10]

    rx=[x[3],x[5]]

    p = plot(x,y, widen = false)
    vspan!(p, rx, label="", alpha=0.2)

    ref_ylims = (y[1], y[end])
    ref_xlims = (x[1].instant.periods.value, x[end].instant.periods.value)
    @test Plots.ylims(p) == ref_ylims
    @test Plots.xlims(p) == ref_xlims
end # testset


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

JLLWrappers.@generate_main_file_header("MongoC")
JLLWrappers.@generate_main_file("MongoC", UUID("90100e71-7732-535a-9be7-2e9affd1cfc1"))
end  # module MongoC_jll


text:
using Yao
using YaoBlocks.AD: apply_back
using ChainRulesCore
using Flux

# helper function with home-grown derivative
function apply_circuit(reg, params; circ::AbstractBlock)
    circ = deepcopy(circ)
    dispatch!(circ, params)
    ψ = ArrayReg(copy(reg))
    ψ |> circ
    state(ψ)
end

# derivative of apply_circuit
function ChainRulesCore.rrule(::typeof(apply_circuit),
        reg, params; circ::AbstractBlock)
    circ = deepcopy(circ)
    out = apply_circuit(reg, params; circ)
    out, function apply_circuit_pullback(outδ)
        dispatch!(circ, params)
        (_, regδ), paramsδ = apply_back((ArrayReg(copy(out)), ArrayReg(copy(outδ))), circ)
        return (NoTangent, state(regδ), paramsδ)
    end
end

# empty circuit
empty_circ = chain()
# zero state
zero_ket(nqb) = reshape([convert(ComplexF64, i == 1) for i = 1:2^nqb], 2^nqb, 1)

struct Circuit
    θ ::Vector
    circ ::AbstractBlock
    enccirc ::AbstractBlock
end

# constructors
Circuit(circ ::AbstractBlock) = Circuit(
    zeros(nparameters(circ)),
    circ,
    empty_circ(nqubits(circ))
)
Circuit(circ ::AbstractBlock, enccirc ::AbstractBlock) = Circuit(
    zeros(nparameters(circ)),
    circ,
    enccirc
)

# model call
function (m::Circuit)(x)
    init = zero_ket(nqubits(m.enccirc))
    enc = apply_circuit(init, x, circ=m.enccirc)
    ret = apply_circuit(enc, m.θ, circ=m.circ)
end

# trainable parameters
Flux.trainable(m::Circuit) = (m.θ,)

# other overloads
Yao.nparameters(m::Circuit) = nparameters(m.circ)

max_stars_repo_name: lmondada/YaoFlux

max_issues_repo_name: lmondada/YaoFlux

max_forks_repo_name: lmondada/YaoFlux


text:
# Neighbor finders

export
    NoNeighborFinder,
    find_neighbors,
    DistanceNeighborFinder,
    DistanceVecNeighborFinder,
    TreeNeighborFinder,
    CellListMapNeighborFinder

"""
    NoNeighborFinder()

Placeholder neighbor finder that returns no neighbors.
When using this neighbor finder, ensure that `nl_only` for the interactions is
set to `false`.
"""
struct NoNeighborFinder <: NeighborFinder end

"""
    find_neighbors(system, neighbor_finder, current_neighbors=nothing, step_n=0; parallel=true)

Obtain a list of close atoms in a system.
Custom neighbor finders should implement this function.
"""
function find_neighbors(s::System,
                        nf::NoNeighborFinder,
                        current_neighbors=nothing,
                        step_n::Integer=0;
                        kwargs...)
    return nothing
end

Base.show(io::IO, neighbor_finder::NoNeighborFinder) = print(io, typeof(neighbor_finder))

"""
    DistanceNeighborFinder(; nb_matrix, matrix_14, n_steps, dist_cutoff)

Find close atoms by distance.
"""
struct DistanceNeighborFinder{D} <: NeighborFinder
    nb_matrix::BitArray{2}
    matrix_14::BitArray{2}
    n_steps::Int
    dist_cutoff::D
end

function DistanceNeighborFinder(;
                                nb_matrix,
                                matrix_14=falses(size(nb_matrix)),
                                n_steps=10,
                                dist_cutoff)
    return DistanceNeighborFinder{typeof(dist_cutoff)}(nb_matrix, matrix_14, n_steps, dist_cutoff)
end

function find_neighbors(s::System,
                        nf::DistanceNeighborFinder,
                        current_neighbors=nothing,
                        step_n::Integer=0;
                        parallel::Bool=true)
    !iszero(step_n % nf.n_steps) && return current_neighbors

    if isnothing(current_neighbors)
        neighbors = NeighborList()
    else
        neighbors = current_neighbors
    end
    empty!(neighbors)

    sqdist_cutoff = nf.dist_cutoff ^ 2

    if parallel && nthreads() > 1
        nl_threads = [Tuple{Int, Int, Bool}[] for i in 1:nthreads()]

        @threads for i in 1:length(s)
            nl = nl_threads[threadid()]
            ci = s.coords[i]
            nbi = @view nf.nb_matrix[:, i]
            w14i = @view nf.matrix_14[:, i]
            for j in 1:(i - 1)
                r2 = sum(abs2, vector(ci, s.coords[j], s.box_size))
                if r2 <= sqdist_cutoff && nbi[j]
                    push!(nl, (i, j, w14i[j]))
                end
            end
        end

        for nl in nl_threads
            append!(neighbors, nl)
        end
    else
        for i in 1:length(s)
            ci = s.coords[i]
            nbi = @view nf.nb_matrix[:, i]
            w14i = @view nf.matrix_14[:, i]
            for j in 1:(i - 1)
                r2 = sum(abs2, vector(ci, s.coords[j], s.box_size))
                if r2 <= sqdist_cutoff && nbi[j]
                    push!(neighbors, (i, j, w14i[j]))
                end
            end
        end
    end

    return neighbors
end

"""
    DistanceVecNeighborFinder(; nb_matrix, matrix_14, n_steps, dist_cutoff)

Find close atoms by distance in a GPU and differentiable safe manner.
"""
struct DistanceVecNeighborFinder{D, B, I} <: NeighborFinder
    nb_matrix::B
    matrix_14::B
    n_steps::Int
    dist_cutoff::D
    is::I
    js::I
end

function DistanceVecNeighborFinder(;
                                nb_matrix,
                                matrix_14=falses(size(nb_matrix)),
                                n_steps=10,
                                dist_cutoff)
    n_atoms = size(nb_matrix, 1)
    if isa(nb_matrix, CuArray)
        is = cu(hcat([collect(1:n_atoms) for i in 1:n_atoms]...))
        js = cu(permutedims(is, (2, 1)))
        m14 = cu(matrix_14)
    else
        is = hcat([collect(1:n_atoms) for i in 1:n_atoms]...)
        js = permutedims(is, (2, 1))
        m14 = matrix_14
    end
    return DistanceVecNeighborFinder{typeof(dist_cutoff), typeof(nb_matrix), typeof(is)}(
            nb_matrix, m14, n_steps, dist_cutoff, is, js)
end

# Find the boundaries of an ordered list of integers
function findboundaries(nbs_ord, n_atoms)
    inds = zeros(Int, n_atoms)
    atom_i = 1
    for (nb_i, nb_ai) in enumerate(nbs_ord)
        while atom_i < nb_ai
            inds[atom_i] = nb_i
            atom_i += 1
        end
    end
    while atom_i < (n_atoms + 1)
        inds[atom_i] = length(nbs_ord) + 1
        atom_i += 1
    end
    return inds
end

function find_neighbors(s::System,
                        nf::DistanceVecNeighborFinder,
                        current_neighbors=nothing,
                        step_n::Integer=0;
                        kwargs...)
    !iszero(step_n % nf.n_steps) && return current_neighbors

    n_atoms = length(s)
    sqdist_cutoff = nf.dist_cutoff ^ 2
    sqdists = sqdistance.(nf.is, nf.js, (s.coords,), (s.box_size,))

    close = sqdists .< sqdist_cutoff
    close_nb = close .* nf.nb_matrix
    eligible = tril(close_nb, -1)

    fa = Array(findall(!iszero, eligible))
    nbsi, nbsj = getindex.(fa, 1), getindex.(fa, 2)
    order_i = sortperm(nbsi)
    weights_14 = @view nf.matrix_14[fa]

    nbsi_ordi, nbsj_ordi = nbsi[order_i], nbsj[order_i]
    sortperm_j = sortperm(nbsj_ordi)
    weights_14_ordi = @view weights_14[order_i]
    atom_bounds_i = findboundaries(nbsi_ordi, n_atoms)
    atom_bounds_j = findboundaries(view(nbsj_ordi, sortperm_j), n_atoms)

    return NeighborListVec(nbsi_ordi, nbsj_ordi, atom_bounds_i, atom_bounds_j,
                            sortperm_j, weights_14_ordi)
end

function allneighbors(n_atoms)
    nbs_all = findall(!iszero, tril(ones(Bool, n_atoms, n_atoms), -1))
    nbsi, nbsj = getindex.(nbs_all, 1), getindex.(nbs_all, 2)
    order_i = sortperm(nbsi)

    nbsi_ordi, nbsj_ordi = nbsi[order_i], nbsj[order_i]
    sortperm_j = sortperm(nbsj_ordi)
    atom_bounds_i = findboundaries(nbsi_ordi, n_atoms)
    atom_bounds_j = findboundaries(view(nbsj_ordi, sortperm_j), n_atoms)
    weights_14 = nothing

    return NeighborListVec(nbsi_ordi, nbsj_ordi, atom_bounds_i, atom_bounds_j,
                            sortperm_j, weights_14)
end

"""
    TreeNeighborFinder(; nb_matrix, matrix_14, n_steps, dist_cutoff)

Find close atoms by distance using a tree search.
"""
struct TreeNeighborFinder{D} <: NeighborFinder
    nb_matrix::BitArray{2}
    matrix_14::BitArray{2}
    n_steps::Int
    dist_cutoff::D
end

function TreeNeighborFinder(;
                            nb_matrix,
                            matrix_14=falses(size(nb_matrix)),
                            n_steps=10,
                            dist_cutoff)
    return TreeNeighborFinder{typeof(dist_cutoff)}(nb_matrix, matrix_14, n_steps, dist_cutoff)
end

function find_neighbors(s::System,
                        nf::TreeNeighborFinder,
                        current_neighbors=nothing,
                        step_n::Integer=0;
                        parallel::Bool=true)
    !iszero(step_n % nf.n_steps) && return current_neighbors

    if isnothing(current_neighbors)
        neighbors = NeighborList()
    else
        neighbors = current_neighbors
    end
    empty!(neighbors)

    dist_unit = unit(first(first(s.coords)))
    bv = ustrip.(dist_unit, s.box_size)
    btree = BallTree(ustripvec.(s.coords), PeriodicEuclidean(bv))
    dist_cutoff = ustrip(dist_unit, nf.dist_cutoff)

    if parallel && nthreads() > 1
        nl_threads = [Tuple{Int, Int, Bool}[] for i in 1:nthreads()]

        @threads for i in 1:length(s)
            nl = nl_threads[threadid()]
            ci = ustrip.(s.coords[i])
            nbi = @view nf.nb_matrix[:, i]
            w14i = @view nf.matrix_14[:, i]
            idxs = inrange(btree, ci, dist_cutoff, true)
            for j in idxs
                if nbi[j] && i > j
                    push!(nl, (i, j, w14i[j]))
                end
            end
        end

        for nl in nl_threads
            append!(neighbors, nl)
        end
    else
        for i in 1:length(s)
            ci = ustrip.(s.coords[i])
            nbi = @view nf.nb_matrix[:, i]
            w14i = @view nf.matrix_14[:, i]
            idxs = inrange(btree, ci, dist_cutoff, true)
            for j in idxs
                if nbi[j] && i > j
                    push!(neighbors, (i, j, w14i[j]))
                end
            end
        end
    end

    return neighbors
end

# Find neighbor lists using CellListMap.jl
"""
    CellListMapNeighborFinder(; nb_matrix, matrix_14, n_steps, dist_cutoff, x0, unit_cell)

Find close atoms by distance, and store auxiliary arrays for in-place threading. `x0` and `unit_cell`
are optional initial coordinates and system unit cell that improve the first approximation of the
cell list structure. The unit cell can be provided as a three-component vector of box sides on each
direction, in which case the unit cell is considered `OrthorhombicCell`, or as a unit cell matrix,
in which case the cell is considered a general `TriclinicCell` by the cell list algorithm.

### Example

```julia-repl
julia> coords
15954-element Vector{SVector{3, Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}}}:
 [2.5193063341012127 nm, 3.907448346081021 nm, 4.694954671434135 nm]
 [2.4173958848835233 nm, 3.916034913604175 nm, 4.699661024574953 nm]
 ⋮
 [1.818842280373283 nm, 5.592152965227421 nm, 4.992100424805031 nm]
 [1.7261366568663976 nm, 5.610326185704369 nm, 5.084523386833478 nm]

julia> box_size
3-element SVector{3, Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}} with indices SOneTo(3):
              5.676 nm
             5.6627 nm
             6.2963 nm

julia> neighbor_finder = CellListMapNeighborFinder(
           nb_matrix=s.neighbor_finder.nb_matrix, matrix_14=s.neighbor_finder.matrix_14,
           n_steps=10, dist_cutoff=1.2u"nm",
           x0 = coords, unit_cell = box_size
       )
CellListMapNeighborFinder{Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}, 3, Float64}
  Size of nb_matrix = (15954, 15954)
  n_steps = 10
  dist_cutoff = 1.2 nm

```
"""
mutable struct CellListMapNeighborFinder{N, T} <: NeighborFinder
    nb_matrix::BitArray{2}
    matrix_14::BitArray{2}
    n_steps::Int
    dist_cutoff::T
    # Auxiliary arrays for multi-threaded in-place updating of the lists
    cl::CellListMap.CellList{N, T}
    aux::CellListMap.AuxThreaded{N, T}
    neighbors_threaded::Vector{NeighborList}
end

# This function sets up the box structure for CellListMap. It uses the unit cell
# if it is given, or guesses a box size from the number of particles, assuming
# that the atomic density is similar to that of liquid water at ambient conditions.
function CellListMapNeighborFinder(;
                                   nb_matrix,
                                   matrix_14=falses(size(nb_matrix)),
                                   n_steps=10,
                                   x0=nothing,
                                   unit_cell=nothing,
                                   number_of_batches=(0, 0), # (0, 0): use default heuristic
                                   dist_cutoff::T) where T
    np = size(nb_matrix, 1)
    if isnothing(unit_cell)
        if unit(dist_cutoff) == NoUnits
            side = max(2 * dist_cutoff, (np * 0.01) ^ (1 / 3))
        else
            side = max(2 * dist_cutoff, uconvert(unit(dist_cutoff), (np * 0.01u"nm^3") ^ (1 / 3)))
        end
        sides = SVector(side, side, side)
        box = CellListMap.Box(sides, dist_cutoff)
    else
        box = CellListMap.Box(unit_cell, dist_cutoff)
    end
    if isnothing(x0)
        x = [ustrip.(box.unit_cell_max) .* rand(SVector{3, T}) for _ in 1:np]
    else
        x = x0
    end
    # Construct the cell list for the first time, to allocate
    cl = CellList(x, box; parallel=true, nbatches=number_of_batches)
    return CellListMapNeighborFinder{3, T}(
        nb_matrix, matrix_14, n_steps, dist_cutoff,
        cl, CellListMap.AuxThreaded(cl),
        [NeighborList(0, [(0, 0, false)]) for _ in 1:CellListMap.nbatches(cl)])
end

"""
    push_pair!(neighbor::NeighborList, i, j, nb_matrix, matrix_14)

Add pair to pair list. If the buffer size is large enough, update element, otherwise
push new element to `neighbor.list`.
"""
function push_pair!(neighbors::NeighborList, i, j, nb_matrix, matrix_14)
    if nb_matrix[i, j]
        push!(neighbors, (Int(i), Int(j), matrix_14[i, j]))
    end
    return neighbors
end

# This is only called in the parallel case
function reduce_pairs(neighbors::NeighborList, neighbors_threaded::Vector{NeighborList})
    neighbors.n = 0
    for i in 1:length(neighbors_threaded)
        append!(neighbors, neighbors_threaded[i])
    end
    return neighbors
end

"""
    find_neighbors(s::System,
                    nf::CellListMapNeighborFinder,
                    current_neighbors=nothing,
                    step_n::Integer=0;
                    parallel::Bool=true)

Find neighbors using `CellListMap`, without in-place updating. Should be called only
the first time the cell lists are built. Modifies the mutable `nf` structure.
"""
function find_neighbors(s::System,
                        nf::CellListMapNeighborFinder,
                        current_neighbors=nothing,
                        step_n::Integer=0;
                        parallel::Bool=true)
    !iszero(step_n % nf.n_steps) && return current_neighbors

    if isnothing(current_neighbors)
        neighbors = NeighborList()
    else
        neighbors = current_neighbors
    end
    aux = nf.aux
    cl = nf.cl
    neighbors.n = 0
    neighbors_threaded = nf.neighbors_threaded
    if parallel
        for i in 1:length(neighbors_threaded)
            neighbors_threaded[i].n = 0
        end
    else
        neighbors_threaded[1].n = 0
    end

    box = CellListMap.Box(s.box_size, nf.dist_cutoff; lcell=1)
    cl = UpdateCellList!(s.coords, box, cl, aux; parallel=parallel)

    map_pairwise!(
        (x, y, i, j, d2, pairs) -> push_pair!(pairs, i, j, nf.nb_matrix, nf.matrix_14),
        neighbors, box, cl;
        reduce=reduce_pairs,
        output_threaded=neighbors_threaded,
        parallel=parallel
    )

    nf.cl = cl
    return neighbors
end

function Base.show(io::IO, neighbor_finder::NeighborFinder)
    println(io, typeof(neighbor_finder))
    println(io,"  Size of nb_matrix = " , size(neighbor_finder.nb_matrix))
    println(io,"  n_steps = " , neighbor_finder.n_steps)
    print(io,"  dist_cutoff = ", neighbor_finder.dist_cutoff)
end


text:
abstract CodingScheme

################################################################################
#               ____ _  _ ____ ____ ___  _ _  _ ____                           #
#               |___ |\ | |    |  | |  \ | |\ | | __                           #
#               |___ | \| |___ |__| |__/ | | \| |__]                           #
################################################################################

type Gray <: CodingScheme end

function encode( ::Type{Gray}, n::Integer )
    n $ (n >> 1)
end

function decode( ::Type{Gray}, n::Integer )
    p = n
     while (n >>= 1) != 0
         p $= n
     end
     return p
end

encode{T<:CodingScheme}( ::Type{T}, N::AbstractVector ) = [ encode( T, n) for n in N ]
decode{T<:CodingScheme}( ::Type{T}, N::AbstractVector ) = [ decode( T, n) for n in N ]


text:
# Copyright 2017 The dm_control Authors.
# Copyright (c) 2019 Colin Summers, The Contributors of LyceumMuJoCo
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#    http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or  implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# ============================================================================

"""
    $(TYPEDEF)

Swing up and balance an unactuated pole by applying forces to a cart at its base. The
physical model conforms to [Neuronlike adaptive elements that can solve difficult learning
control problems](https://ieeexplore.ieee.org/document/6313077) (Barto et al., 1983).

# Spaces

* **State: (7, )**
* **Action: (1, )**
* **Observation: (5, )**
"""
struct CartpoleSwingup{Sim,OSpace} <: AbstractMuJoCoEnvironment
    sim::Sim
    observationspace::OSpace
    function CartpoleSwingup(sim::MJSim)
        observationspace = MultiShape(
            pos = MultiShape(
                cart = ScalarShape(Float64),
                pole_zz = ScalarShape(Float64),
                pole_xz = ScalarShape(Float64),
            ),
            vel = MultiShape(
                cart = ScalarShape(Float64),
                pole_ang = ScalarShape(Float64),
            ),
        )
        env = new{typeof(sim), typeof(observationspace)}(sim, observationspace)
        reset!(env)
    end
end

CartpoleSwingup() = first(tconstruct(CartpoleSwingup, 1))

function LyceumBase.tconstruct(::Type{CartpoleSwingup}, n::Integer)
    modelpath = joinpath(@__DIR__, "cartpole.xml")
    [CartpoleSwingup(s) for s in tconstruct(MJSim, n, modelpath)]
end

@inline getsim(env::CartpoleSwingup) = env.sim

@inline LyceumBase.observationspace(env::CartpoleSwingup) = env.observationspace

@inline function LyceumBase.getobservation!(obs, env::CartpoleSwingup)
    @boundscheck checkaxes(observationspace(env), obs)

    @uviews obs begin
        sobs = observationspace(env)(obs)
        sobs.pos.cart = env.sim.dn.qpos[:slider]
        sobs.pos.pole_zz = env.sim.dn.xmat[:z, :z, :pole_1]
        sobs.pos.pole_xz = env.sim.dn.xmat[:x, :z, :pole_1]
        copyto!(sobs.vel, env.sim.d.qvel)
    end
    obs
end

@inline function LyceumBase.getreward(state, action, obs, env::CartpoleSwingup)
    @boundscheck begin
        checkaxes(statespace(env), state)
        checkaxes(actionspace(env), action)
        checkaxes(observationspace(env), obs)
    end

    sobs = observationspace(env)(obs)

    upright = (sobs.pos.pole_zz + 1) / 2

    centered = tolerance(sobs.pos.cart, margin = 2)
    centered = (1 + centered) / 2

    small_control = tolerance(first(action), margin = 1, value_at_margin = 0, sigmoid = quadratic)
    small_control = (4 + small_control) / 5

    small_velocity = min(tolerance(sobs.vel.pole_ang, margin = 5))
    small_velocity = (1 + small_velocity) / 2

    mean(upright) * small_control * small_velocity * centered
end

function LyceumBase.reset!(env::CartpoleSwingup)
    reset_nofwd!(env.sim)

    qpos = env.sim.dn.qpos
    @uviews qpos begin
        qpos[:hinge_1] = pi
    end

    forward!(env.sim)

    env
end

function LyceumBase.randreset!(rng::Random.AbstractRNG, env::CartpoleSwingup)
    reset_nofwd!(env.sim)

    qpos = env.sim.dn.qpos
    @uviews qpos begin
        qpos[:slider] = 0.01 * randn(rng)
        qpos[:hinge_1] = pi + 0.01*randn(rng)
    end

    randn!(rng, env.sim.d.qvel)
    env.sim.d.qvel .*= 0.01

    forward!(env.sim)

    env
end


text:
"""
    manifest!"..."

Create and activate a throwaway environment based on the contents of a `Manifest.toml`.
The environment is deleted at exit.

Use this macro to paste a manifest file into a Julia script,
ensuring that it will always run with the exact set of desired dependencies.
"""
macro manifest!_str(s::AbstractString)
    quote
        env = _env_dir()
        write(joinpath(env, "Manifest.toml"), $s)
        deps = Dict(name => pkgs[1]["uuid"] for (name, pkgs) in TOML.parse($s))
        project = Dict("deps" => deps)
        open(io -> TOML.print(io, project), joinpath(env, "Project.toml"), "w")
        _activate(env)
    end
end

"""
    project!"..."

Create and activate a throwaway environment based on the contents of a `Project.toml`.
The environment is deleted at exit.

Use this macro to paste a project file into a Julia script,
ensuring that it will always run with the defined set of desired dependencies.
"""
macro project!_str(s::AbstractString)
    quote
        env = _env_dir()
        write(joinpath(env, "Project.toml"), $s)
        _activate(env)
    end
end

function _env_dir()
    dir = mktempdir()
    atexit(() -> rm(dir; recursive=true, force=true))
    return dir
end

function _activate(env::AbstractString)
    Pkg.activate(env)
    Pkg.instantiate()
end


text:
using PGFPlots

using CSV
using DataFrames
using Missings
using Query

function only(a)
    @assert length(a) == 1
    return first(a)
end

# filename = Pkg.dir("Multilane", "data/cor_trend_Monday_2_Apr_19_49.csv")
# filename = Pkg.dir("Multilane", "data/cor_trend_Tuesday_3_Apr_07_20.csv")
# filename = Pkg.dir("Multilane", "data/cor_trend_Tuesday_3_Apr_20_35.csv")
filename = Pkg.dir("Multilane", "data/cor_trend_Friday_27_Apr_03_39.csv")

data = CSV.read(filename)

data[:time_to_lane] = data[:steps_to_lane].*data[:dt]
data[:brakes_per_step] = data[:nb_brakes]./data[:n_steps]

data = @from i in data begin
    @where i.solver != "outcome"
    @select i
    @collect DataFrame
end

combined = by(data, [:cor, :solver]) do df
    return DataFrame(reached_lane=sum(df[:terminal].=="lane")/nrow(df),
                     n=nrow(df))
end

@show sort(combined, cols=[:solver])

plts = Plots.Plot[]
for df in groupby(combined, [:solver])
    df = sort(df, cols=[:cor])
    name = only(unique(df[:solver]))
    # 95% confidence region from Hoeffding Bound
    confidence = 0.5
    n = only(unique(df[:n]))
    confidence_radius = sqrt(log((1.0-confidence)/2)/(-2*n))
    plt = Plots.Linear(convert(Vector{Float64}, df[:cor]),
                       df[:reached_lane],
                       errorBars=ErrorBars(y=fill(confidence_radius, nrow(df))),
                       legendentry=name)
    push!(plts, plt)
end
a = Axis(plts,
         ylabel="Fraction Successful",
         xlabel="Correlation",
         legendPos="south west"
        )

file = tempname()*".pdf"
save(file, a)
run(`xdg-open $file`)


text:
using DifferentialEquations
function f(du,u,p,t)
    du[1] = -u[1]
end
u0 = [10.0]
const V = 1
prob = ODEProblem(f,u0,(0.0,10.0))

sol = solve(prob,Tsit5())
using Plots; gr()
plot(sol)

condition(u,t,integrator) = t==4 && u[1]/V<4
affect!(integrator) = integrator.u[1] += 10
cb = DiscreteCallback(condition,affect!)

sol = solve(prob,Tsit5(),tstops=[4.0],callback=cb)
using Plots; gr()
plot(sol)

println(sol(4.00000))
println(sol(4.000000000001))

function f(du,u,p,t)
    du[1] = -u[1]/6
end
u0 = [10.0]
const V = 1
prob = ODEProblem(f,u0,(0.0,10.0))

sol = solve(prob,Tsit5())
using Plots; gr()
plot(sol)

sol = solve(prob,Tsit5(),tstops=[4.0],callback=cb)
using Plots; gr()
plot(sol)

using SciMLTutorials
SciMLTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])

max_stars_repo_name: SciML/SciMLTutorialsOutput

max_issues_repo_name: SciML/SciMLTutorialsOutput

max_forks_repo_name: SciML/SciMLTutorialsOutput


text:
abstract type AbstractConsolidation end
abstract type AbstractConsolidator end
"""
    RawConsolidationParameter

An optimizer attribute used for raw parameters of the consolidator. Defers to `RawParameter`.
"""
struct RawConsolidationParameter <: ConsolidationParameter
    name::Any
end

function MOI.get(consolidator::AbstractConsolidator, param::RawConsolidationParameter)
    name = Symbol(param.name)
    if !(name in fieldnames(typeof(consolidator)))
        error("Unrecognized parameter name: $(name) for consolidator $(typeof(consolidator)).")
    end
    return getfield(consolidator, name)
end

function MOI.set(consolidator::AbstractConsolidator, param::RawConsolidationParameter, value)
    name = Symbol(param.name)
    if !(name in fieldnames(typeof(consolidator)))
        error("Unrecognized parameter name: $(name) for consolidator $(typeof(consolidator)).")
    end
    setfield!(consolidator, name, value)
    return nothing
end

# No consolidation
# ------------------------------------------------------------
"""
    NoRegularization

Empty functor object for running the L-shaped algorithm without consolidation.

"""
struct NoConsolidation <: AbstractConsolidation end

function consolidate!(::AbstractLShaped, ::NoConsolidation)
    return nothing
end

function add_cut!(::AbstractLShaped, ::NoConsolidation, ::AbstractHyperPlane)
    return nothing
end

function add_cut!(::AbstractLShaped, ::NoConsolidation, ::Integer, ::AbstractHyperPlane)
    return nothing
end
# Consolidation
# ------------------------------------------------------------
"""
    Consolidation

Functor object for using consolidation in an L-shaped algorithm. Create by supplying a [`Consolidate`](@ref) object through `consolidate ` in the `LShapedSolver` factory function and then pass to a `StochasticPrograms.jl` model.

...
# Algorithm parameters
- `tresh::T` = 0.95: Relative amount of redundant cuts in a former iteration required to consider the iteration redundant
- `at::Int = 5.0`: Number of times an iteration can be redundant before consolidation is triggered
- `rebuild::Function = at_tolerance()`: Function deciding when the master model should be rebuilt according to performed consolidations
...
"""
struct Consolidation{T <: AbstractFloat} <: AbstractConsolidation
    cuts::Vector{Vector{AnySparseOptimalityCut{T}}}
    feasibility_cuts::Vector{Vector{SparseFeasibilityCut{T}}}
    consolidated::Vector{Bool}
    redundance_count::Vector{Int}
    redundance_treshold::T
    consolidation_trigger::Int
    rebuild::Function

   function Consolidation(::Type{T}, tresh::T, at::Int, rebuild::Function) where T <: AbstractFloat
        cuts = Vector{Vector{AnySparseOptimalityCut{T}}}()
        push!(cuts, Vector{AnySparseOptimalityCut{T}}())
        feasibility_cuts = Vector{Vector{SparseFeasibilityCut{T}}}()
        push!(feasibility_cuts, Vector{SparseFeasibilityCut{T}}())
        return new{T}(cuts, feasibility_cuts, [false], [0], tresh, at, rebuild)
    end
end

function consolidate!(lshaped::AbstractLShaped, consolidation::Consolidation{T}) where T <: AbstractFloat
    for i in findall(map(!,consolidation.consolidated))
        nc = sum(num_subproblems.(consolidation.cuts[i]))
        nc == 0 && continue
        count = 0
        for cut in consolidation.cuts[i]
            # Check for redundance
            if !active(lshaped, cut)
                count += num_subproblems(cut)
            end
        end
        if count/nc >= consolidation.redundance_treshold
            consolidation.redundance_count[i] += 1
        end
        if consolidation.redundance_count[i] == consolidation.consolidation_trigger
            consolidated_cut = aggregate(consolidation.cuts[i])
            empty!(consolidation.cuts[i])
            push!(consolidation.cuts[i], consolidated_cut)
            consolidation.consolidated[i] = true
        end
    end
    for i in findall(c->length(c)>0, consolidation.feasibility_cuts)
        inactive = Vector{Int}()
        for (j,cut) in enumerate(consolidation.feasibility_cuts[i])
            # Check for redundance
            if !active(lshaped, cut)
                push!(inactive, j)
            end
        end
        deleteat!(consolidation.feasibility_cuts[i], inactive)
    end
    if minimum_requirements(lshaped,consolidation) && consolidation.rebuild(lshaped, consolidation)
        # Rebuild sparingly
        rebuild_master!(lshaped, consolidation)
        lshaped.data.consolidations += 1
    end
    # Prepare memory for next iteration
    allocate!(consolidation)
    return nothing
end

function allocate!(consolidation::Consolidation{T}) where T <: AbstractFloat
    push!(consolidation.cuts, Vector{AnySparseOptimalityCut{T}}())
    push!(consolidation.feasibility_cuts, Vector{SparseFeasibilityCut{T}}())
    push!(consolidation.consolidated, false)
    push!(consolidation.redundance_count, 0)
    return nothing
end

function add_cut!(::AbstractLShaped, ::Consolidation, ::AbstractHyperPlane)
    return nothing
end

function add_cut!(lshaped::AbstractLShaped, consolidation::Consolidation{T}, cut::HyperPlane{FeasibilityCut}) where T <: AbstractFloat
    push!(consolidation.feasibility_cuts[timestamp(lshaped)], cut)
    return nothing
end

function add_cut!(lshaped::AbstractLShaped, consolidation::Consolidation{T}, cut::AnySparseOptimalityCut) where T <: AbstractFloat
    push!(consolidation.cuts[timestamp(lshaped)], cut)
    return nothing
end

function rebuild_master!(lshaped::AbstractLShaped, consolidation::Consolidation)
    # Remove all cutsa
    remove_cut_constraints!(lshaped)
    # Readd consolidated cuts
    readd_cuts!(lshaped, consolidation)
    return nothing
end

num_optimalitycuts(consolidation::Consolidation) = sum(length.(consolidation.cuts))
num_feasibilitycuts(consolidation::Consolidation) = sum(length.(consolidation.feasibility_cuts))
function num_cutconstraints(consolidation::Consolidation)
    return num_optimalitycuts(consolidation) + num_feasibilitycuts(consolidation)
end

function num_cutconstraints(lshaped::AbstractLShaped)
    return length(lshaped.cut_constraints)
end
# Rebuild functions
# ------------------------------------------------------------
function minimum_requirements(lshaped::AbstractLShaped, consolidation::Consolidation)
    return (num_optimalitycuts(consolidation) > 0 &&
        count(consolidation.consolidated) > 0) ||
        (num_feasibilitycuts(consolidation) > 0 &&
         num_feasibilitycuts(consolidation)/num_cutconstraints(lshaped) <= 0.1)
end
"""
    at_tolerance(τ = 0.4, miniter = 0)

Rebuild master when at least nconsolidations*`miniter` iterations has passed and the ratio of number of cuts in the consolidated collection and the number of cuts in the master model has decreased below `τ`.

"""
at_tolerance() = at_tolerance(0.4, 0)
function at_tolerance(τ, miniter)
    return (lshaped, consolidation) -> begin
        return lshaped.data.iterations >= (lshaped.data.consolidations+1)*miniter && num_cutconstraints(consolidation)/num_cutconstraints(lshaped) <= τ
    end
end

# API
# ------------------------------------------------------------
"""
    DontConsolidate

Factory object for [`NoConsolidation`](@ref). Passed by default to `consolidate` in the `LShapedSolver` factory function.

"""
struct DontConsolidate <: AbstractConsolidator end

function (::DontConsolidate)(::Type{T} where T <: AbstractFloat)
    return NoConsolidation()
end

"""
    Consolidate

Factory object for [`Consolidation`](@ref). Pass to `consolidate` in the `LShapedSolver` factory function. See ?Consolidation for parameter descriptions.

"""
mutable struct Consolidate <: AbstractConsolidator
    redundance_treshold::Float64
    consolidation_trigger::Int
    rebuild::Function
end
Consolidate(; tresh::Float64 = 0.95, at::Int = 5, rebuild = at_tolerance()) = Consolidate(tresh, at, rebuild)

function (consolidator::Consolidate)(::Type{T}) where T <: AbstractFloat
    return Consolidation(T,
                         consolidator.redundance_treshold,
                         consolidator.consolidation_trigger,
                         consolidator.rebuild)
end


text:
module VTKDataTypes

using WriteVTK
using ResumableFunctions
using GeometryTypes
using Colors
using IterTools
using StaticArrays

import GeometryTypes: decompose

include("vtkcelltypes.jl")
include("types.jl")
include("utils.jl")
include("validation.jl")
include("typeconversion.jl")
include("decompose.jl")
include("triangulate.jl")
include("GLMesh.jl")
#include("surface.jl")

export  AbstractVTKData,
        AbstractStaticVTKData,
        AbstractTimeSeriesVTKData,
        AbstractVTKMultiblockData,
        AbstractVTKSimpleData,
        AbstractVTKUnstructuredData,
        AbstractVTKStructuredData,
        AbstractVTKRectilinearData,
        VTKUnstructuredData,
        VTKPolyData,
        VTKStructuredData,
        VTKRectilinearData,
        VTKUniformRectilinearData,
        VTKImageData,
        VTKMultiblockData,
        VTKTimeSeriesData,
        VTK_CELL_TYPE,
        endof,
        same_geometry,
        same_ordered_geometry,
        same_geometry_shape,
        same_data_shape,
        coherent, extents,
        cell_extents,
        dim,
        num_of_points,
        num_of_cells,
        num_of_point_vars,
        num_of_cell_vars,
        cell_type,
        cell_connectivity,
        has_var,
        var_dim,
        is_homogeneous,
        filter_cells!,
        keep_volume_cells_only!,
        keep_face_cells_only!,
        get_cell_ids,
        get_lowest_index,
        get_highest_index,
        is_valid_cell,
        add_new_cell!,
        remove_cell!,
        add_point_id_offset!,
        append,
        num_of_blocks,
        insert_new_block!,
        remove_block!,
        timespan,
        num_of_timesteps,
        insert_timed_data!,
        remove_timed_data!,
        is_valid,
        simple_block_generator,
        increase_resolution!,
        triangular,
        triangulate,
        decompose,
        GLMesh,
        bb,
        pseudo_center
end


text:
using Serialization

function get_trials()
    unique([String(split(t, "_")[1]) for t in readdir(LOGGER_FOLDER) if isfile(LOGGER_FOLDER * t)])
end

function get_snapshot(snapshot_id)
    for filename in readdir(SNAPSHOTS_FOLDER)
        _, file_snapshot_id = split(filename, "_")
        if file_snapshot_id != snapshot_id continue end

        snapshot = deserialize(SNAPSHOTS_FOLDER * filename)
        return deepcopy(snapshot)
    end
end

function get_snapshot_ids()
    snapshot_ids = []

    for filename in readdir(SNAPSHOTS_FOLDER)
        _, snapshot_id = split(filename, "_")
        push!(snapshot_ids, string(snapshot_id))
    end

    return unique(snapshot_ids)
end

function get_snapshot_ids(trial_id::String)
    snapshot_ids = []

    for filename in readdir(SNAPSHOTS_FOLDER)
        current_trial_id, snapshot_id = split(filename, "_")

        if current_trial_id == trial_id
            push!(snapshot_ids, string(snapshot_id))
        end
    end

    return snapshot_ids
end

function save_offspring_log()
    offspring_parents = Dict()

    for trial_id in get_trials()
        @info trial_id
        i = 0

        for file in readdir(LOGGER_FOLDER)
            current_trial_id, _ = split(file, "_")
            if current_trial_id != trial_id continue end

            logged_organisms = deserialize(LOGGER_FOLDER * file)

            for (_, organism) in logged_organisms
                for parent_genotype in organism.parent_genotypes
                    key = (parent_genotype, organism.genotype)

                    if haskey(offspring_parents, key)
                        offspring_parents[key] += 1
                    else
                        offspring_parents[key] = 1
                    end
                end
            end

            i += 1

            @info "$i \t $file"
        end

        serialize(CALCULATED_FOLDER * "parent_offspring_ocurrances", offspring_parents)
    end
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
data = readline()

function group(xs)
    ys = []
    push!(ys, [xs[1]])
    prev_y = xs[1]
    for x in xs[2:end]
        if x == prev_y
            push!(ys[end], x)
        else
            push!(ys, [x])
            prev_y = x
        end
    end
    ys
end

function step_audioactive(s)
    groups = group(s)
    join(["$(length(g))$(g[1])" for g in groups], "")
end

function solve_1()
    s = data
    for i in 1:50
        s = step_audioactive(s)
    end
    s
end

solve_1() |> length |> println

max_stars_repo_name: pLOPeGG/JuliaOfCode

max_issues_repo_name: pLOPeGG/JuliaOfCode

max_forks_repo_name: pLOPeGG/JuliaOfCode


text:
using JuMP
using Cbc
using Dates

dt1=now(Dates.UTC)
println(dt1)

function readData()
    arr = Array{Int, 2}(undef, 0,0 )
    if isfile("instances/instance_30_1.txt")
        myFile = open("instances/instance_30_1.txt")
        data = readlines(myFile)
        n = length(data)
        arr = fill(0,(n - 3, n -10))

        i = 0
        for line in data
            if line == ""
                continue
            end
            i += 1
            var = split(line, " ")
            for j in 1:length(var)
                arr[i, j] = parse(Int, var[j])
            end
        end
    end
    return arr
end

function square(n::Int64)
    return n*n
end

#Import data
function generate_aerodrones(array_of_aerodrone, map_limit)

    aerodrones = []
    for aerodrone_counter in 1:array_of_aerodrone[1,1]
        push!(aerodrones,
            Dict(
                "id" => aerodrone_counter,
                "x" => array_of_aerodrone[aerodrone_counter+7,1],
                "y" => array_of_aerodrone[aerodrone_counter+7,2],
                "region" => array_of_aerodrone[6,aerodrone_counter]
            )
        )
    end
    println("Cities Generated:", size(aerodrones)[1])
    return aerodrones
end

function calculate_distance_between_two_points(point1, point2)
    return sqrt((((point2[1] - point1[1]))^2) + (((point2[2] - point1[2]))^2))
end

#Calculate traversal paths (measure fitness)
function calculate_chromosome_travel_distance(chromosome)
    travel_distance = 0
    for geneId in 1:length(chromosome) - 1
        point1 = (
            aerodrones[chromosome[geneId]]["x"],
            aerodrones[chromosome[geneId]]["y"]
            )
        point2 = (
            aerodrones[chromosome[geneId + 1]]["x"],
            aerodrones[chromosome[geneId + 1]]["y"]
            )
        travel_distance += calculate_distance_between_two_points(point1, point2)
    end
    return travel_distance
end

#Determine whether the limit length is exceeded
function is_transcend_limit_distance(point1,point2)
    travel_distance = calculate_distance_between_two_points(point1, point2)
    if(travel_distance>array_of_aerodrone[7][1])
        return true
    else
        return false
    end
end

# We shuffle the chromosome here
function shuffle_chromosome(chromosome)
    for i in 1:size(chromosome)[1]
        random_point = rand(1:array_of_aerodrone[1,1]-2, 1)[1]
        chromosome[i], chromosome[random_point] = chromosome[random_point], chromosome[i]
    end
    return chromosome
end

#Generates a specified number of random chromosomes
function generate_initial_population(initial_population_size)
    chromosomes = []
    for population_counter in 1:initial_population_size
        #Cut out the beginning and the end
        initial_chromosome = [1:length(aerodrones);]
        loc_start = findfirst(el -> el == array_of_aerodrone[2,1], initial_chromosome)
        splice!(initial_chromosome, loc_start)
        loc_destination = findfirst(el -> el == array_of_aerodrone[3,1], initial_chromosome)
        splice!(initial_chromosome, loc_destination)
        #Shuffle
        upper_chromosome = rand(array_of_aerodrone[4,1]-2:array_of_aerodrone[1,1]-2, 1)[1]
        chromosome = shuffle_chromosome(initial_chromosome)
        chromosome = chromosome[1:upper_chromosome]
        #Joining together
        chromosome = vcat(array_of_aerodrone[2,1],chromosome,array_of_aerodrone[3,1])
        travel_distance = calculate_chromosome_travel_distance(chromosome)
        push!(chromosomes,
            Dict(
                "id_generation" => 0,
                "id_initial" => population_counter,
                "chromosome" => chromosome,
                "distance" => travel_distance
            )
        )
    end
    return chromosomes
end

#Crossover
function crossover(parent_one_chromosome, parent_two_chromosome, crossover_point)
    offspring_part_one = parent_one_chromosome[1:crossover_point]
    for gene in offspring_part_one
        if gene in parent_two_chromosome
            gene_loc = findfirst(el -> el == gene, parent_two_chromosome)
            splice!(parent_two_chromosome, gene_loc)
        end
    end
    return vcat(offspring_part_one, parent_two_chromosome)
end

#Random point mutation
function mutate(offspring)
    random_mutation_point1 = rand(2:length(offspring)-1)
    random_mutation_point2 = rand(2:length(offspring)-1)
    offspring[random_mutation_point1], offspring[random_mutation_point2] = offspring[random_mutation_point2], offspring[random_mutation_point1]
    return offspring
end

#Determine whether to pass all areas
function is_traverse_all_region(chromosome)
    region_number=array_of_aerodrone[5,1]
    region_array=zeros(Int8,1,region_number)
    for geneID in 1:length(chromosome)
         aerodrone_ID = chromosome[geneID]
         aerodrone_region = array_of_aerodrone[6,aerodrone_ID]
         if(aerodrone_region != 0 && aerodrone_region != region_array[aerodrone_region])
            region_array[aerodrone_region] = aerodrone_region
        end
    end
    for geneID in 1:length(region_array)
        if(region_array[geneID]==0)
            #print("Not traverse all region. ")
            return false
        end
    end
    return true
end

function is_route_transcended(chromosome)
    for geneId in 1:length(chromosome) - 1
        point1 = (
            aerodrones[chromosome[geneId]]["x"],
            aerodrones[chromosome[geneId]]["y"]
            )
        point2 = (
            aerodrones[chromosome[geneId + 1]]["x"],
            aerodrones[chromosome[geneId + 1]]["y"]
            )
        if is_transcend_limit_distance(point2,point1)
            #print("Transcend limit distance.")
            return true
        end
    end
    return false
end

#evolution
function evolve(generation_count, offsprings_count, crossover_point)
    println("Start generating:")
    for generation in 1:generation_count

        for offspring_count in 1:offsprings_count
            #Mating of two chromosomes (crossover + mutation)
            random_parent_one_id = rand(1:size(chromosomes)[1])
            random_parent_two_id = rand(1:size(chromosomes)[1])
            random_parent_one = copy(chromosomes[random_parent_one_id]["chromosome"])
            random_parent_two = copy(chromosomes[random_parent_two_id]["chromosome"])
            offspring = crossover(random_parent_one, random_parent_two, crossover_point)
            offspring = mutate(offspring)
            #Ensure that new offspring pass through all regions
            while(!is_traverse_all_region(offspring) || is_route_transcended(offspring))
                #println()
                #println("illegal route, regenerate offspring")
                random_parent_one_id = rand(1:size(chromosomes)[1])
                random_parent_one = copy(chromosomes[random_parent_one_id]["chromosome"])
                new_chromosome = generate_initial_population(1)
                random_parent_two = copy(new_chromosome[1]["chromosome"])
                #println(random_parent_two)
                offspring = crossover(random_parent_one, random_parent_two, crossover_point)
                offspring = mutate(offspring)
            end
            push!(chromosomes,
                Dict(
                    "id_generation" => generation,
                    "id_offspring" => offspring_count,
                    "chromosome" => offspring,
                    "distance" => calculate_chromosome_travel_distance(offspring)
                    )
            )
        end
        sort!(chromosomes, by=x -> x["distance"], rev=false)
        splice!(chromosomes, offsprings_count÷2+1:size(chromosomes)[1])
        println("generation: ", generation, " ", "Best route: ", chromosomes[1]["distance"])
    end
end

# creating 10 cities randomly
generation_count = 120
offsprings_count = 100000
crossover_point = 2
array_of_aerodrone = readData()
aerodrones = generate_aerodrones(array_of_aerodrone, 500)
# generating 10 initial chromosomes
chromosomes = generate_initial_population(offsprings_count)

# we are running GA for
# 5 generations
# 5 offsprings per generation
# random crossover point as 2

evolve(generation_count, offsprings_count, crossover_point)
println("--------------------------------------------------------")
println("Optimal route:", chromosomes[1]["chromosome"])
println("travel_distance:", chromosomes[1]["distance"])
if(is_route_transcended(chromosomes[1]["chromosome"]))
    println("Not avaliable since transcend the limit distance")
end
if(is_traverse_all_region(chromosomes[1]["chromosome"]))
    println("Congratulations! This route has traversed all the region")
end

print("For drawing: ")
for i in 2:length(chromosomes[1]["chromosome"])
    print(chromosomes[1]["chromosome"][i-1]," ",chromosomes[1]["chromosome"][i]," ")
end
println()

dt2 = now(Dates.UTC)
println(dt2-dt1)


text:
module NewtonRaphson

export newton_raphson

function newton_raphson(f::Function, x0::Number, fprime::Function, args::Tuple=();
                tol::AbstractFloat=1e-8, maxiter::Integer=50, eps::AbstractFloat=1e-10)
    for _ in 1:maxiter
        yprime = fprime(x0, args...)
        if abs(yprime) < eps
           warn("First derivative is zero")
           return x0
        end
        y = f(x0, args...)
        x1 = x0 - y/yprime
        if abs(x1-x0) < tol
          return x1
        end
        x0 = x1
    end
    error("Max iteration exceeded")
end

end


text:
# Before running this, please make sure to activate and instantiate the environment
# corresponding to [this `Project.toml`](https://raw.githubusercontent.com/alan-turing-institute/DataScienceTutorials.jl/master/Project.toml) and [this `Manifest.toml`](https://raw.githubusercontent.com/alan-turing-institute/DataScienceTutorials.jl/master/Manifest.toml)
# so that you get an environment which matches the one used to generate the tutorials:
#
# ```julia
# cd("DataScienceTutorials") # cd to folder with the *.toml
# using Pkg; Pkg.activate("."); Pkg.instantiate()
# ```

using MLJ
import RDatasets: dataset
using PrettyPrinting
using Random

Random.seed!(3203)
X = randn(20, 2)
y = vcat(-ones(10), ones(10))

using PyPlot

figure(figsize=(8,6))

ym1 = y .== -1
ym2 = .!ym1
plot(X[ym1, 1], X[ym1, 2], ls="none", marker="o")
plot(X[ym2, 1], X[ym2, 2], ls="none", marker="x")

X = MLJ.table(X)
y = categorical(y);

@load SVC pkg=LIBSVM

svc_mdl = SVC()
svc = machine(svc_mdl, X, y)

fit!(svc);

ypred = predict(svc, X)
misclassification_rate(ypred, y)

rc = range(svc_mdl, :cost, lower=0.1, upper=5)
tm = TunedModel(model=svc_mdl, ranges=[rc], tuning=Grid(resolution=10),
                resampling=CV(nfolds=3, rng=33), measure=misclassification_rate)
mtm = machine(tm, X, y)

fit!(mtm)

ypred = predict(mtm, X)
misclassification_rate(ypred, y)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
module LiterateDollop
export optfun

ρ = (sqrt(5) - 1)/2

fw1(a, b) = a + (1 - ρ)*(b - a)
fw2(a, b) = b - (1 - ρ)*(b - a)

function optfun(f, a, b; ϵ=1e-5)
    @assert b > a
    k = 1
    w1 = fw1(a, b)
    w2 = fw2(a, b)

    while abs(b - a) > ϵ
        if f(w1) < f(w2)
            b = w2
            w2 = w1
            w1 = fw1(a, b)
        else
            a = w1
            w1 = w2
            w2 = fw2(a, b)
        end
    end
    return a
end

end # module


text:
# The main workhorse
using Plots, DifferentialEquations
include("ecme-dox.jl")

u0 = collect(u0_tup)

param = Params(
  pDOX = DOXParams(DOX=0.0),
  pCK=CKParams(V_ATPASE_CYTO= 1E-5),
  pC1=C1Params(SCALE=1e3 * 50),
  pC3=C3Params(SCALE=60e3 * 15),
  pC4=C4Params(SCALE=60e3 * 15),
  pC5 = C5Params(ρF1=1.5),
  pANT = ANTParams(VMAX=5E-3),
  pSODi = SODParams(ET=1.43E-3 * 0.65),
  BCL=1000.0)

using JLD2
JLD2.@load "pacing.jld2" solPacing

@unpack A_CAP_V_MYO_F = param

tspan = (0.0, 10000.0)
prob = ODEProblem(rhs, u0, tspan, param)

integrator = init(prob, Rodas5(); reltol=1e-9, abstol=1e-9, dt=0.01, progress=true, tstops=tspan[1]:1000.0:tspan[end])
for (u,t) in tuples(integrator)
    (vm, dpsi, m_na, h_na, j_na, x_k, na_i, k_i, ca_i, ca_jsr, ca_nsr, ca_ss, ltr_ca, htr_ca,
     po1, po2, pc2, c0, c1, c2, c3, c4, o, cca0, cca1, cca2, cca3, y_ca, p0, p1, p2, p3, n1,
     adp_i, adp_ic, crp_i, crp_ic, ca_m, adp_m, nadh, isoc, akg, scoa, suc, fum, mal, oaa,
     sox_m, sox_i, h2o2_i, gsh_i, Q_n, Qdot_n, QH2_n, QH2_p, Qdot_p, Q_p, b1, b2, b3, b4,
     fes_ox, cytc1_ox, cytc_ox) = u
    @show t, vm, dpsi
    @show adp_i, adp_m
    @show Q_n, QH2_n
end

solPacing = solve(prob, Rodas5(); reltol=1e-9, abstol=1e-9, dt=0.01, progress=true, tstops=tspan[1]:1000.0:tspan[end])
pVolt = plot(solPacing, vars=(0, nameLUT[:vm]), label="Membrane Potential", lw=1)
plot!(pVolt, solPacing, vars=(0, nameLUT[:dpsi]), label="Mitochondrial Potential", lw=1, xlabel="time(ms)", ylabel = "Voltage (mV)")
savefig(pVolt, "voltage.png")

plot(solPacing, vars=(0, [nameLUT[:adp_i], nameLUT[:adp_m]]), label=["ADPi" "ADPm"], lw=1)
plot(solPacing, vars=(0, nameLUT[:nadh]), label="[NADH]", lw=1)
plot(solPacing, vars=(0, nameLUT[:ca_m]), label="[Ca]m", lw=1)
plot(solPacing, vars=(0, nameLUT[:isoc]), label="[ISOC]", lw=1)
plot(solPacing, vars=(0, nameLUT[:akg]), label="[AKG]", lw=1)
plot(solPacing, vars=(0, nameLUT[:scoa]), label="[SCoA]", lw=1)
plot(solPacing, vars=(0, nameLUT[:suc]), label="[SUC]", lw=1)
plot(solPacing, vars=(0, nameLUT[:fum]), label="[FUM]", lw=1)
plot(solPacing, vars=(0, nameLUT[:mal]), label="[MAL]", lw=1)
plot(solPacing, vars=(0, nameLUT[:oaa]), label="[OAA]", lw=1)
plot(solPacing, vars=(0, [nameLUT[:Q_n], nameLUT[:QH2_n]]), lw=1, label=["Q_n" "QH2_n"])
plot(solPacing, vars=(0, [nameLUT[:b1], nameLUT[:b2], nameLUT[:b3], nameLUT[:b4]]), lw=1)
plot(solPacing, vars=(0, [nameLUT[:fes_ox], nameLUT[:cytc1_ox], nameLUT[:cytc_ox]]), lw=1, label=["fes_ox" "cytc1_ox" "cytc_ox"])
plot(solPacing, vars=(0, [nameLUT[:sox_i], nameLUT[:sox_m]]), lw=1, label=["sox_i" "sox_m"])
plot(solPacing, vars=(0, nameLUT[:sox_i]), lw=1, label=["sox_i"])
plot(solPacing, vars=(0, nameLUT[:gsh_i]), lw=1)

plot()
plot(t-> ecme_dox(solPacing(t), param, t)[rateMap[:iKatp]], solPacing.t[1], solPacing.t[end], label="iKatp")
plot!(solPacing.t, 8 .- solPacing[nameLUT[:adp_i],:], label="ATPi")
plot!(solPacing, vars=(0, nameLUT[:adp_i]), lw=1, label="ADPi", xlabel="Time (ms)", ylabel=" uA/uF or mM")
savefig("katp.png")

plot!(t-> ecme_dox(solPacing(t), param, t)[rateMap[:vANT]], solPacing.t[1], solPacing.t[end], label="vANT")
plot!(t-> ecme_dox(solPacing(t), param, t)[rateMap[:vATPase]], solPacing.t[1], solPacing.t[end], label="vATPase")
plot!(t-> ecme_dox(solPacing(t), param, t)[rateMap[:vSL]], solPacing.t[1], solPacing.t[end], label="vSL")
plot!(t-> ecme_dox(solPacing(t), param, t)[rateMap[:jUp]], solPacing.t[1], solPacing.t[end], label="jUp")
plot!(t-> A_CAP_V_MYO_F * ecme_dox(solPacing(t), param, t)[rateMap[:iPCa]], solPacing.t[1], solPacing.t[end], label="iPCa")
plot!(t-> A_CAP_V_MYO_F * ecme_dox(solPacing(t), param, t)[rateMap[:iNaK]], solPacing.t[1], solPacing.t[end], label="iNaK")
plot!(t-> ecme_dox(solPacing(t), param, t)[rateMap[:vAM]], solPacing.t[1], solPacing.t[end], label="vAM")
plot!(t-> ecme_dox(solPacing(t), param, t)[rateMap[:d_adp_i]], solPacing.t[1], solPacing.t[end], label="d_adp_i")

plot(sol, vars=(0, nameLUT[:dpsi]), label="Mitochondrial Potential", lw=1)
plot(sol, vars=(0, nameLUT[:nadh]), label="[NADH]", lw=1)
plot(sol, vars=(0, nameLUT[:adp_m]), label="[ADP]m", lw=1)
plot(sol, vars=(0, nameLUT[:isoc]), label="[ISOC]", lw=1)
plot(sol, vars=(0, nameLUT[:akg]), label="[AKG]", lw=1)
plot(sol, vars=(0, nameLUT[:scoa]), label="[SCoA]", lw=1)
plot(sol, vars=(0, nameLUT[:suc]), label="[SUC]", lw=1)
plot(sol, vars=(0, nameLUT[:fum]), label="[FUM]", lw=1)
plot(sol, vars=(0, nameLUT[:mal]), label="[MAL]", lw=1)
plot(sol, vars=(0, nameLUT[:oaa]), label="[OAA]", lw=1)
plot(sol, vars=(0, [nameLUT[:Q_n], nameLUT[:QH2_n]]), lw=1, label=["Q_n" "QH2_n"])
plot(sol, vars=(0, [nameLUT[:b1], nameLUT[:b2], nameLUT[:b3], nameLUT[:b4]]), lw=1)
plot(sol, vars=(0, [nameLUT[:fes_ox], nameLUT[:cytc1_ox], nameLUT[:cytc_ox]]), lw=1, label=["fes_ox" "cytc1_ox" "cytc_ox"])
plot(sol, vars=(0, [nameLUT[:sox_i], nameLUT[:sox_m]]), lw=1, label=["sox_i" "sox_m"])
plot(sol, vars=(0, nameLUT[:sox_i]), lw=1, label=["sox_i"])
plot(sol, vars=(0, nameLUT[:gsh_i]), lw=1)

using JLD2
@save "ROS-oscillation.jld2" sol

# Get the rates
plot(t-> ecme_dox(sol(t), param, t)[end-15], sol.t[1], sol.t[end])
plot!(t-> ecme_dox(sol(t), param, t)[end-22], sol.t[1], sol.t[end])
plot(t-> ecme_dox(sol(t), param, t)[end], tspan[1], tspan[end])


text:
using ArrayInterfaceCore
using ArrayInterfaceBandedMatrices
using BandedMatrices
using Test

B=BandedMatrix(Ones(5,5), (-1,2))
B[band(1)].=[1,2,3,4]
B[band(2)].=[5,6,7]
@test ArrayInterfaceCore.has_sparsestruct(B)
rowind,colind=ArrayInterfaceCore.findstructralnz(B)
@test [B[rowind[i],colind[i]] for i in 1:length(rowind)]==[5,6,7,1,2,3,4]
B=BandedMatrix(Ones(4,6), (-1,2))
B[band(1)].=[1,2,3,4]
B[band(2)].=[5,6,7,8]
rowind,colind=ArrayInterfaceCore.findstructralnz(B)
@test [B[rowind[i],colind[i]] for i in 1:length(rowind)]==[5,6,7,8,1,2,3,4]
@test ArrayInterfaceCore.isstructured(typeof(B))
@test ArrayInterfaceCore.fast_matrix_colors(typeof(B))


text:
@testset "test common" begin

networks_opf = [build_pm_opf_model(case) for case in cases]
opf_ac_objective = [14676.9, 27564.91]
@testset "opf ac - $(i)" for (i,network) in enumerate(networks_opf)

    result = run_ac_opf(network, nlp_solver)

    @test isapprox(result["termination_status"], LOCALLY_SOLVED)
    @test isapprox(result["objective"], opf_ac_objective[i]; atol = 1e0)
end

@testset "set_start_values - $(i)" for (i,network) in enumerate(networks_opf)
    network = deepcopy(network)

    set_start_values!(network)
    result = run_ac_opf(network, nlp_solver)

    @test isapprox(result["termination_status"], LOCALLY_SOLVED)
    @test isapprox(result["objective"], opf_ac_objective[i]; atol = 1e0)
end

opf_ac_tight_objective = [14676.9, 28347.12]
@testset "opf ac tight - $(i)" for (i,network) in enumerate(networks_opf)
    network = deepcopy(network)

    tighten_constraints!(network)
    result = run_ac_opf(network, nlp_solver)

    @test isapprox(result["termination_status"], LOCALLY_SOLVED)
    @test isapprox(result["objective"], opf_ac_tight_objective[i]; atol = 1e0)
end

opf_ac_no_rate_a_objective = [14676.9, 27546.46]
@testset "opf ac no flow limits - $(i)" for (i,network) in enumerate(networks_opf)
    network = deepcopy(network)

    deactivate_rate_a!(network)
    result = run_ac_opf(network, nlp_solver)

    @test isapprox(result["termination_status"], LOCALLY_SOLVED)
    @test isapprox(result["objective"], opf_ac_no_rate_a_objective[i]; atol = 1e0)
end

mn_size = [2, 20]
@testset "build multinetwork - $(i)" for (i,network) in enumerate(networks)
    network = deepcopy(network)

    network["gen_contingencies"] = network["gen_contingencies"][1:min(10, length(network["gen_contingencies"]))]
    network["branch_contingencies"] = network["branch_contingencies"][1:min(10, length(network["branch_contingencies"]))]

    mn_case = build_scopf_multinetwork(network)

    @test isapprox(length(mn_case["nw"]), mn_size[i]; atol = 1e0)
end

@testset "rate_a deactivate / activate - $(i)" for (i,network) in enumerate(networks)
    network = deepcopy(network)

    deactivate_rate_a!(network)
    @test all(!haskey(branch, "rate_a") for (i,branch) in network["branch"])

    activate_rate_a!(network)
    @test all(haskey(branch, "rate_a") for (i,branch) in network["branch"])
end

first_cont_id = [3, 203]
@testset "contingency_order - $(i)" for (i,network) in enumerate(networks)
    order = contingency_order(network)
    @test isapprox(order[i].idx, first_cont_id[i]; atol = 1e0)
end

first_gen_cont_id = [3, 60]
@testset "contingency_order, gen only - $(i)" for (i,network) in enumerate(networks)
    network = deepcopy(network)
    network["branch_contingencies"] = []

    order = contingency_order(network)
    @test isapprox(order[i].idx, first_gen_cont_id[i]; atol = 1e0)
end

first_branch_cont_id = [9, 61]
@testset "contingency_order, branch only - $(i)" for (i,network) in enumerate(networks)
    network = deepcopy(network)
    network["gen_contingencies"] = []

    order = contingency_order(network)
    @test isapprox(order[i].idx, first_branch_cont_id[i]; atol = 1e0)
end

ref_bus_id = ["1", "17"]
@testset "correct_voltage_angles - $(i)" for (i,network) in enumerate(networks)
    network = deepcopy(network)
    update_data!(network, solutions[i])

    correct_voltage_angles!(network)

    @test isapprox(network["bus"][ref_bus_id[i]]["va"], 0.0; atol = 1e-2)
end

opf_p_delta_abs_max = [0.26746780838927353, 0.07392932526496376]
opf_q_delta_abs_max = [0.41571449226280965, 1.1036712386853447]
solution1_lines = [25,594]
@testset "write_solution1 - $(i)" for (i,network) in enumerate(networks)
    network = deepcopy(network)

    #deactivate_rate_a!(network)

    result = run_dc_opf(network, lp_solver)
    #result = run_dc_opf(network, lp_solver)
    @test isapprox(result["termination_status"], OPTIMAL)

    update_active_power_data!(network, result["solution"])

    balance = compute_power_balance_deltas!(network)
    @test isapprox(balance.p_delta_abs_max, opf_p_delta_abs_max[i])
    @test isapprox(balance.q_delta_abs_max, opf_q_delta_abs_max[i])

    correct_network_solution!(network)
    write_solution1(network, solution_file="solution1-tmp.txt")

    open("solution1-tmp.txt", "r") do file
        lines = countlines(file)
        @test lines == solution1_lines[i]
    end
    rm("solution1-tmp.txt")
end

end # close test group


text:
struct ConfusionMatrix
    pos_class::Symbol
    neg_class::Symbol

    tp::Int
    fp::Int
    tn::Int
    fn::Int
end

get_positive(x::ConfusionMatrix) = x.tp + x.fn

get_negative(x::ConfusionMatrix) = x.fp + x.tn

tpr(x::ConfusionMatrix) = x.tp / (x.tp + x.fn)
sensitivity = tpr
recall = tpr

fpr(x::ConfusionMatrix) = x.fp / (x.fp + x.tn)
specificity(x) = 1 - fpr(x)
tnr = specificity

get_n(x::ConfusionMatrix) = get_positive(x) + get_negative(x)

f1_score(x::ConfusionMatrix) = 2 * x.tp / (2 * x.tp + x.fp + x.fn)

ConfusionMatrix(tp, fp, tn, fn; pos_class = :outlier, neg_class = :inlier) = ConfusionMatrix(pos_class, neg_class, tp, fp, tn, fn)

Base.convert(x::Type{ConfusionMatrix}, d::Dict) = ConfusionMatrix(d)

function ConfusionMatrix(d::Dict)
    return ConfusionMatrix(convert(Int64, d["tp"]),
                           convert(Int64, d["fp"]),
                           convert(Int64, d["tn"]),
                           convert(Int64, d["fn"]))
end

function ConfusionMatrix(classification, ground_truth; pos_class = :outlier, neg_class = :inlier)
    @assert islabelenc(classification, LABEL_ENCODING)
    @assert islabelenc(ground_truth, LABEL_ENCODING)
    @assert length(classification) == length(ground_truth)

    tp = sum((classification .== pos_class) .& (ground_truth .== pos_class))
    fp = sum((classification .== pos_class) .& (ground_truth .== neg_class))
    tn = sum((classification .== neg_class) .& (ground_truth .== neg_class))
    fn = sum((classification .== neg_class) .& (ground_truth .== pos_class))

    return ConfusionMatrix(pos_class, neg_class, tp, fp, tn, fn)
end

function cohens_kappa(x::ConfusionMatrix)
    p_0 = (x.tp + x.tn) / (get_negative(x) + get_positive(x))
    p_e = ((x.tp + x.fp) * get_positive(x) + (x.tn + x.fn) * get_negative(x)) / (get_negative(x) + get_positive(x))^2
    return (p_0 - p_e) / (1 - p_e)
end

function matthews_corr(x::ConfusionMatrix)::Float64
    d = sqrt((x.tp + x.fp) * (x.tp + x.fn) * (x.tn + x.fp) * (x.tn + x.fn))
    return d == 0 ? 0.0 : (x.tp * x.tn - x.fp * x.fn) / d
end

function roc_auc(predictions::Vector{T}, ground_truth::Vector{Symbol}; fpr = 1.0, normalize = false) where T <: Real
    @assert length(predictions) == length(ground_truth)
    @assert islabelenc(ground_truth, LABEL_ENCODING)
    lm = labelmap(ground_truth)
    r = ROCAnalysis.roc(predictions[lm[:outlier]], predictions[lm[:inlier]])
    return AUC(r, pfa=fpr, normalize = normalize)
end


text:
#=
  Routines for finding out what cholmod would do when Cholesky factorizing a matrix

  These are calling cholmod's analyze routinem which is overkill.
  They could be made faster with more direct access to Cholmod.

  By Daniel Spielman.
=#

"""
  nnzL, flops = ask_cholmod(mat)

Estimate the number of nonzeros in the cholfact factorization of mat,
along with the number of flops needed to compute it.
Does this through a call to the analyze routine of cholmod.
Note that this is much faster than actually computing the factorization
"""
function ask_cholmod(sdd)
    ba = SuiteSparse.CHOLMOD
    #cm = ba.common()
    cm = ba.defaults(ba.common_struct[Threads.threadid()])

    anal = ba.analyze(ba.Sparse(lap(sdd)), cm);

    # s_anal = unsafe_load(get(anal.p))
    s_anal = unsafe_load(pointer(anal))

    n = Int(s_anal.n)

    nnzL = 0
    flops = 0

    for i in 1:n
        nzl = Int(unsafe_load(s_anal.ColCount,i))
        nnzL += nzl
        flops += nzl^2
    end

    return nnzL, flops
end

"""
  p = cholmod_perm(mat)

Return the permutation that cholmod would apply.
"""
function cholmod_perm(sdd)
    ba = SuiteSparse.CHOLMOD
    #cm = ba.common()
    cm = ba.defaults(ba.common_struct[Threads.threadid()])

    anal = ba.analyze(ba.Sparse(lap(sdd)), cm);

    #s_anal = unsafe_load(get(anal.p))
    s_anal = unsafe_load(pointer(anal))

    n = Int(s_anal.n)

    p = zeros(Int,n)
    for i in 1:n
        p[i] = unsafe_load(s_anal.Perm,i)+1
    end

    return p
end


text:
using CSV, DataFrames,StatsBase
using DiseaseOutbreak
using Plots

include("data_dir.jl")
# include(joinpath(@__DIR__, "..", "src", "distancing.jl"))

memphis = CSV.read(joinpath(memphis_data_dir,"Memphis_residences.csv"),type=Float64,
                    missingstring="NA")

names(memphis)
memphispos = convert(Matrix{Float64},dropmissing(memphis[:,3:4]))

n = size(memphispos,1);
idx = sample(1:n,20000)
sampleMemphispos = memphispos[idx,:]
status,d,pos = makeCitizens(sampleMemphispos)
infectCitizens!(status,10)
followCitizens!(status,d,100,0.5,14,0.000015)

uninfected(status)

newInfections(status)

cumulativeInfections(status)

gif(infectionSpread(status,pos,80),joinpath(@__DIR__, "..", "plots", "memphis.gif"),fps=15)