text:
# Example for Variable Projection:
# Author: Ziyi Yin, Rafael Orozco, Philipp Witte, Mathias Louboutin, Gabrio Rizzuti, Felix J. Herrmann
# Date: April, 2020
#
using DrWatson
@quickactivate "ExtSrcImg"

include("gen_geometry.jl")
include("modeling_extended_source_spg_raf_function.jl")

using PyPlot, FFTW, Images
using JUDI, JOLI, JUDI4Flux
using LinearAlgebra, PyPlot, Random, Statistics
using IterativeSolvers
using JLD
using MAT
using Random

function weight(n, d, xs, zs, nsrc,mode)
	if mode == "medical"
		zs_ave = zeros(1,n[2],nsrc); #hacky way
		zs_ave .= mean(zs)

		xs_ave = zeros(1,n[2],nsrc); #hacky way
		xs_ave .= mean(zs)

		x = zeros(n[1],1,nsrc);
	    x .= reshape((collect(1:n[1]).-1)*d[1], :, 1)
		z = zeros(1,n[2],nsrc);
	    z .= reshape((collect(1:n[2]).-1)*d[2], 1, :)

		W = zeros(Float32,n[1],n[2],nsrc)
		#W .= sqrt.((z.-zs_ave).^2f0)
		 W .= sqrt.((x.-zs_ave).^2f0.+(z.-zs_ave).^2f0)

		# x = zeros(n[1],1,nsrc);
	 #    x .= reshape((collect(1:n[1]).-1)*d[1], :, 1)
		# z = zeros(1,n[2],nsrc);
	 #    z .= reshape((collect(1:n[2]).-1)*d[2], 1, :)
		# W = zeros(Float32,n[1],n[2],nsrc)
		# W .= sqrt.((x.-xs).^2f0.+(z.-zs).^2f0)

	elseif mode == "seismic"
		x = zeros(n[1],1,nsrc);
	    x .= reshape((collect(1:n[1]).-1)*d[1], :, 1)
		z = zeros(1,n[2],nsrc);
	    z .= reshape((collect(1:n[2]).-1)*d[2], 1, :)
		W = zeros(Float32,n[1],n[2],nsrc)
		W .= sqrt.((x.-xs).^2f0.+(z.-zs).^2f0)
		#W = sqrt.(r.^2f0.+δ^2f0)./δ
	end

    return reshape(W,n[1],n[2],1,nsrc)
end

function CG(d_obs, q0, m0, W; maxiter=5,λ=1f2)
	# use CG to solve damped least square normal equations
	r1 = d_obs-𝒢(q0, m0);
	r = 𝒢_T(r1,m0)-λ*W.^2f0.*q0;
	p = r;
	r_vec = zeros(maxiter+1);
	println("residual: $(norm(r1)) before updating q")
	for i = 1:maxiter

		if norm(r1)<1f-3
			break
		end

		temp = sum(sum(r.*r,dims=1),dims=2);
		temp1 = 𝒢(p, m0);
		temp2 = 𝒢_T(temp1,m0)+λ*W.^2f0.*p;
		alpha = temp./sum(sum(p.*temp2,dims=1),dims=2);
		q0 = q0 + alpha.* p;
		r1 = r1 - alpha.* temp1;
		r = r - alpha .* temp2;
		beta = sum(sum(r.*r,dims=1),dims=2)./temp;
		p = r + beta.*p;
		r_vec[i+1]=norm(r1)
		println("residual: $(norm(r1)) after updating q")

	end
	return q0, r1, r_vec
end

# Objective function: estimate source and compute residual + function value
function objective(d_obs, q0, m0, W; λ=1f2)

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

	# Function value
	fval = .5f0 *sum(r_data.^2) + 0.5f0*λ*sum((W.*r_source).^2)

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

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

	num_simultaneous_sources = size(xsrc_index)[1]
	#check if Simultaneous
	if num_simultaneous_sources > 1 && nsrc == 1
		#for Simultaneous we need to put all the weights in the same slice.
		for j=1:num_simultaneous_sources

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

    return weights
end

function GenSimData(q₀; model=model)   # pass q₀ directly instead of z
    # Simultaneous observed data
    d_obs = F*vec(q₀)
    d_obs = reshape(d_obs, recGeometry.nt[1], nrec, 1, nsrc)  # reshape into tensor
    return d_obs
end

function objective_function_spg(x)
	#cFact = Float32(1 ./model.d[1]^2); #scale model up to make gradient have a normal magnitude

	r = d_obs-𝒢(x, m0);
  	f = .5f0*norm(r)^2;
  	g = -vec(𝒢_T(r,m0));

  	return f, (g/norm(g, Inf));
end

############################################ Models ###############################################

M = load("../data/brain.jld")
m = M["m"]
m = m./1f6
N = load("../data/mute_sp_phantom.jld")
mute = N["mute"]
m0 = zeros(size(m))
m0 .= m[1,1]
#m0 = imfilter(deepcopy(m),Kernel.gaussian(5))
n = size(m)
dx = 50f0 #from hauptmann
dz = 50f0   # in micro meters (μm)
#dx = 687f0 #from hauptmann
#dz = 687f0   # in micro meters (μm)

d            = (dx,dz);
o            = (0.,0.);
#the only difference between the models is the slowness

model = Model(n, d, o, m,nb=120)
model0 = Model(n, d, o, m0,nb=120)

extentx = n[1]*d[1] / 1000 #in mm
extentz = n[2]*d[2] / 1000 #in mm

# receiver sampling and recording time
dtR = 0.02f0  # micro seconds (μs)
timeR = 25
#dtR = 2f0  # micro seconds (μs)
#timeR = 1500f0

timeS = timeR
dtS = dtR

#Mode = "seismic"
mode = "medical"

if mode == "seismic"
	# Set up receiver structure
	xrec,yrec,zrec,nrec,nsrc =GenRecGeometrySeismic();
	recGeometry = Geometry(xrec, yrec, zrec; dt=dtR, t=timeR, nsrc=nsrc)

	# Set up source structure
	xsrc_index,zsrc_index = GenSrcGeometrySeismic(nsrc);

	#generate tensor with sources
	q       = GenSimSourceMulti(xsrc_index, zsrc_index, nsrc, model.n);

	#get muting region
	mute_mask = GenMuteSeismic()

elseif mode == "medical"

	xrec, yrec, zrec, nrec, nsrc, xsrc_index, zsrc_index, nsrc_simult = Gen_Exp_Medical_Skull_Phantom_Source()

	#generate tensor with sources
	q           = GenSimSourceMulti(xsrc_index, zsrc_index, nsrc, model.n);

	#get muting region
	mute_mask = GenMuteMedical_brain()
else
	error("Please choose a correct mode.")
end

recGeometry = Geometry(xrec, yrec, zrec; dt=dtR, t=timeR, nsrc=nsrc)

# Setup wavelet
#f0 = 0.0035  # dominant frequency in [kHz]
f0 = 5  # dominant frequency in [MHz]
#f0 = 0.0030  # dominant frequency in [kHz]
wavelet = ricker_wavelet(timeS,dtS,f0)

######################################### Set up JUDI operators ###################################

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

# Write shots as array
opt = Options(return_array = true, subsampling_factor=10, isic=true)
#opt = Options(return_array = true, isic=true)

# Modeling operators
Pr  = judiProjection(info, recGeometry)
F0  = judiModeling(info, model0; options=opt)	# modeling operator true model
F  = judiModeling(info, model; options=opt)	# modeling operator initial model
Pw = judiLRWF(info, wavelet)

# Combined operators
F = Pr*F*adjoint(Pw)
F0 = Pr*F0*adjoint(Pw)

# Extended modeling CNN layers
𝒢 = ExtendedQForward(F)
𝒢_T = ExtendedQAdjoint(F)

######################################### Run Forward simulations ###################################
# Generate shots
d_obs   = GenSimData(q); # Data in true model

######################################### Generate 3 gradient images ###################################

#run cg to invert for estimated source
q_iter = 25

q_init = zeros(Float32, n[1], n[2]);
q_inv = spg_tv_norm(d_obs, objective_function_spg, sum(q,dims=4)[:,:,1,1], q_init, model0; maxiter=q_iter, tv_norm=0.5)

#get jacobian using estimated source and use it to get gradient
J = judiJacobian(F0, q_inv)
grad_extend=reshape(adjoint(J)* vec(d_obs-𝒢(q_inv,m0)),n[1],n[2])

figure(2);imshow(grad_extend'/norm(grad_extend',Inf),cmap="gray",extent=(0,extentx,extentz,0),vmin=-0.05,vmax=0.05);title("Extended Source Imaging - Inversion");
xlabel("Lateral position [mm]");ylabel("Depth [mm]");
cb = colorbar(fraction=0.038, pad=0.04);
fig = PyPlot.gcf(); fig.set_size_inches(12.0, 7.0)

#### 3.) Extended source gradient - adjoint scheme
q_adj = 𝒢_T(d_obs,m0)

print(vars)
######################################### Visualize results   ###################################

#show models used

v = reshape(sqrt.(1f0./(m.*1f6)),n[1], n[2]);
using PyCall
@pyimport numpy.ma as ma
v1 = pycall(ma.masked_less, Any, v', 2)
v2 = pycall(ma.masked_greater, Any, v', 2)

PyPlot.rc("font", family="serif"); PyPlot.rc("xtick", labelsize=8); PyPlot.rc("ytick", labelsize=8)
fig,ax = PyPlot.subplots()
p2 = ax.imshow(v2,extent=(0,extentx,extentz,0),interpolation="nearest",cmap="winter",vmin=1.5,vmax=1.55);
cb2 = colorbar(p2,fraction=0.037, pad=0.08);
p1 = ax.imshow(v1,extent=(0,extentx,extentz,0),interpolation="nearest",cmap="Reds",vmin=2.45,vmax=2.55)
cb1 = colorbar(p1,fraction=0.043, pad=0.08);
PyPlot.scatter(xrec/1000,zrec/1000,marker=".",color="yellow",label="receivers",s=2);
PyPlot.scatter(xsrc_index*d[1]/1000,zsrc_index*d[2]/1000,label="sources",color="white",s=1);
legend(loc=1,fontsize=7);
xlabel("Lateral position [mm]", fontsize=8);ylabel("Horizontal position [mm]", fontsize=8);
fig = PyPlot.gcf(); fig.set_size_inches(6, 4.5)
savefig("phantom_setup.png")

if mode == "seismic"	# one single extended source in the middle

	q_inv_stack=sum(q_inv,dims=4)[:,:,1,1];
	figure(8);imshow(q_inv_stack',vmin=-2f-3,vmax=2f-3,extent=(0,extentx,extentz,0));title("Extended Source from Inversion");
	xlabel("Lateral position [samples]");ylabel("Depth [samples]");
	colorbar(fraction=0.028, pad=0.04);
	fig = PyPlot.gcf(); fig.set_size_inches(12.0, 7.0)

	q_adj_stack=sum(q_adj,dims=4)[:,:,1,1];
	figure(9);imshow(q_adj_stack',vmin=-40000,vmax=40000,extent=(0,extentx,extentz,0));title("Extended Source from Adjoint");
	xlabel("Lateral position [samples]");ylabel("Depth [samples]");
	colorbar(fraction=0.028, pad=0.04);
	fig = PyPlot.gcf(); fig.set_size_inches(12.0, 7.0)

	if nsrc != 1
		mid_loc = Int(nsrc/2)
		q_inv_single=q_inv[:,:,1,mid_loc];
		figure(18);imshow(q_inv_single',vmin=-2f-3,vmax=2f-3,extent=(0,extentx,extentz,0));title("Extended Source from Inversion");
		xlabel("Lateral position [samples]");ylabel("Depth [samples]");
		colorbar(fraction=0.028, pad=0.04);
		fig = PyPlot.gcf(); fig.set_size_inches(12.0, 7.0)

		q_adj_single=q_adj[:,:,1,mid_loc];
		figure(19);imshow(q_adj_single',vmin=-40000,vmax=40000,extent=(0,extentx,extentz,0));title("Extended Source from Adjoint");
		xlabel("Lateral position [samples]");ylabel("Depth [samples]");
		colorbar(fraction=0.028, pad=0.04);
		fig = PyPlot.gcf(); fig.set_size_inches(12.0, 7.0)
	end

elseif mode == "medical"
	PyPlot.rc("font", family="serif"); PyPlot.rc("xtick", labelsize=8); PyPlot.rc("ytick", labelsize=8)
	figure(); imshow(reshape(q_adj,n[1], n[2])'.*mute',extent=(0,extentx,extentz,0),vmin=-0.8*norm(q_adj,Inf),vmax=0.8*norm(q_adj,Inf));
	xlabel("Lateral position [mm]", fontsize=8);ylabel("Horizontal position [mm]", fontsize=8);
	fig = PyPlot.gcf(); fig.set_size_inches(6, 4.5)
	savefig("phantom_adj_rand.png")

	PyPlot.rc("font", family="serif"); PyPlot.rc("xtick", labelsize=8); PyPlot.rc("ytick", labelsize=8)
	figure(); imshow(reshape(q_inv,n[1], n[2])'.*mute',extent=(0,extentx,extentz,0),vmin=-0.8*norm(q_inv,Inf),vmax=0.8*norm(q_inv,Inf));
	xlabel("Lateral position [mm]", fontsize=8);ylabel("Horizontal position [mm]", fontsize=8);
	fig = PyPlot.gcf(); fig.set_size_inches(6, 4.5)
	savefig("phantom_inv_water.png")

	PyPlot.rc("font", family="serif"); PyPlot.rc("xtick", labelsize=8); PyPlot.rc("ytick", labelsize=8)
	figure(); imshow(reshape(q,n[1], n[2])'.*mute',extent=(0,extentx,extentz,0));
	cb = colorbar(fraction=0.037, pad=0.04);
	title("True Source", fontsize=8);
	xlabel("Lateral position [mm]", fontsize=8);ylabel("Horizontal position [mm]", fontsize=8);
	fig = PyPlot.gcf(); fig.set_size_inches(6.66, 5)
	savefig("phantom_true.png", bbox_inches="tight")
end

######################################### Save results   ########################################
#@save "tv_brain_water.jld" q q_inv q_adj m m0 d n xrec yrec zrec xsrc_index zsrc_index grad_extend mute d_obs


text:
function _register_job_def(config::AWSConfig, input::AbstractArray, expected::AbstractArray)
    @test input == expected
    return REGISTER_JOB_DEF_RESP
end

function _submit_job(config::AWSConfig, input::AbstractArray, expected::AbstractArray)
    @test input == expected
    return SUBMIT_JOB_RESP
end

@testset "run_batch" begin

    @testset "Defaults" begin
        withenv("AWS_BATCH_JOB_ID" => nothing) do
            @test_throws AWSBatch.BatchEnvironmentError run_batch()
        end
    end

    @testset "From Job Definition" begin
        expected_job = [
                "jobName" => "example",
                "jobQueue" => "HighPriority",
                "jobDefinition" => "arn:aws:batch:us-east-1:012345678910:job-definition/sleep60:1",
                "parameters" => Dict{String,String}(),
                "containerOverrides" => Dict(
                    "command" => ["sleep", "60"],
                    "memory" => 128,
                    "vcpus" => 1,
                ),
            ]

        patches = [
            @patch describe_job_definitions(args...) = DESCRIBE_JOBS_DEF_RESP
            @patch submit_job(config, input) = _submit_job(config, input, expected_job)
        ]

        apply(patches) do
            job = run_batch(; name="example", definition="sleep60", queue="HighPriority")
            @test job.id == "24fa2d7a-64c4-49d2-8b47-f8da4fbde8e9"
        end
    end

    @testset "From Current Job" begin
        withenv(BATCH_ENVS...) do
            expected_job = [
                "jobName" => "example",
                "jobQueue" => "HighPriority",
                "jobDefinition" => "arn:aws:batch:us-east-1:012345678910:job-definition/sleep60:1",
                "parameters" => Dict{String,String}(),
                "containerOverrides" => Dict(
                    "command" => ["sleep", "60"],
                    "memory" => 128,
                    "vcpus" => 1,
                ),
            ]

            expected_job_def = [
                "type" => "container",
                "parameters" => Dict{String,String}(),
                "containerProperties" => [
                    "image" => "busybox",
                    "vcpus" => 1,
                    "memory" => 128,
                    "command" => ["sleep", "60"],
                    "jobRoleArn" => "arn:aws:iam::012345678910:role/sleep60",
                ],
                "jobDefinitionName" => "sleep60",
            ]

            patches = [
                @patch instance_region() = "us-east-1"
                @patch describe_jobs(args...) = DESCRIBE_JOBS_RESP
                @patch describe_job_definitions(args...) = Dict("jobDefinitions" => Dict())
                @patch register_job_definition(config, input) =
                    _register_job_def(config, input, expected_job_def)
                @patch submit_job(config, input) = _submit_job(config, input, expected_job)
            ]

            apply(patches) do
                job = run_batch()
                @test job.id == "24fa2d7a-64c4-49d2-8b47-f8da4fbde8e9"
            end
        end
    end

    @testset "Using a Job Definition" begin
        withenv(BATCH_ENVS...) do
            expected_job = [
                "jobName" => "example",
                "jobQueue" => "HighPriority",
                "jobDefinition" => "arn:aws:batch:us-east-1:012345678910:job-definition/sleep60:1",
                "parameters" => Dict{String,String}(),
                "containerOverrides" => Dict(
                    "command" => ["sleep", "60"],
                    "memory" => 128,
                    "vcpus" => 1,
                ),
            ]

            patches = [
                @patch instance_region() = "us-east-1"
                @patch describe_jobs(args...) = DESCRIBE_JOBS_RESP
                @patch describe_job_definitions(args...) = Dict("jobDefinitions" => Dict())
                @patch submit_job(config, input) = _submit_job(config, input, expected_job)
            ]

            apply(patches) do
                definition = JobDefinition("arn:aws:batch:us-east-1:012345678910:job-definition/sleep60:1")
                job = run_batch(definition=definition)
                @test job.id == "24fa2d7a-64c4-49d2-8b47-f8da4fbde8e9"
            end
        end
    end
end

ext: jl


text:
function get_occ_ocl_free(cloud, camera, depth_image, resolution)
    pixels = GL.point_cloud_to_pixel_coordinates(cloud, camera)

    THRESHOLD = 0.01

    occluded = fill(true, size(pixels)[2])
    occupied = fill(false, size(pixels)[2])

    for (i,(x,y)) in enumerate(eachcol(pixels))
        if 1 <= x <= camera.width && 1 <= y <= camera.height
            if abs(cloud[3,i] - depth_image[y,x]) <= resolution
                occupied[i] = true
            end
            if cloud[3, i] < depth_image[y,x] + 0.0
                occluded[i] = false
            end
        end
    end
    occluded = occluded .& (.!(occupied));
    free = .!(occluded) .& .!(occupied)
    occupied, occluded, free
end

export get_occ_ocl_free


text:
"""
	struct PseudoSpectralSolver <: AbstractPropagator

An iteration scheme used for solving modified diffusion equations using pseudo-spectral Fourier methods.
The `method` parameter specifies whether to use a second-order (:RK2) or fourth-order (:RQM4) time-stepping scheme.
"""
struct PseudoSpectralSolver<: AbstractPropagator
	dims   :: NTuple{3,Int}
	method :: Symbol

	# Operator grids for solving MDE
	LW1    :: FieldGrid{Float64}
	LW2    :: FieldGrid{Float64}
	LD1    :: FieldGrid{Float64}
	LD2    :: FieldGrid{Float64}
end

function PseudoSpectralSolver(dims::NTuple{3,<:Integer}; method::Symbol = :RK2)
	if method != :RK2 && method != :RQM4
		throw(ArgumentError("Invalid propagator method `$(method)` for PseudoSpectralSolver."))
	end

	# Operator grids
	LW1 = zeros(Float64, dims)
	LW2 = zeros(Float64, dims)
	LD1 = zeros(Float64, ksize(dims))
	LD2 = zeros(Float64, ksize(dims))

	return PseudoSpectralSolver(dims, method, LW1, LW2, LD1, LD2)
end

PseudoSpectralSolver(Nx::Integer, Ny::Integer, Nz::Integer; kwargs...) = PseudoSpectralSolver((Nx, Ny, Nz); kwargs...)

#==============================================================================#

Base.show(io::IO, solver::PseudoSpectralSolver) = @printf(io, "PseudoSpectralSolver(dims = %s)", solver.dims)

"""
	update!(solver, omega, ksq, vol, b, ds)

Update differential operators for a given set of `omega` and `ksq` grids.
"""
function update!(solver::PseudoSpectralSolver, omega::FieldGrid, ksq::FieldGrid, vol::Real, b::Real, ds::Real)
	@assert size(omega) == solver.dims
	@assert size(ksq) == size(solver.LD1)

	lw_coeff = ds / 2.0
	ld_coeff = ds * b^2 / 6.0

	@. solver.LW1 = exp(-lw_coeff * vol * omega)
	@. solver.LW2 = exp(-lw_coeff * vol * omega / 2.0)

	@. solver.LD1 = exp(-ld_coeff * ksq)
	@. solver.LD2 = exp(-ld_coeff * ksq / 2.0)

	return nothing
end

"""
	propagate!(solver, fft, qin, qout)

Take one step forward along the chain propagator from input `qin` to output `qout`.
Utilizes pre-allocated fast-Fourier transforms of the correct dimension from the provided `FFTHolder`.
"""
function propagate!(solver::PseudoSpectralSolver, plan::FFTHolder, qin::FieldGrid, qout::FieldGrid)
	@assert solver.dims == plan.dims

	# Offload params
	FT, iFT = plan.FT, plan.iFT
	LW1, LW2, LD1, LD2 = solver.LW1, solver.LW2, solver.LD1, solver.LD2

	# Pre-allocated temp grids from FFT
	q1, q2, qr = rgrid(plan, 1), rgrid(plan, 2), rgrid(plan, 3)
	qk = kgrid(plan, 1)

	if solver.method == :RK2
		@. qr = LW1 * qin
		mul!(qk, FT, qr)
		@. qk *= LD1
		mul!(qr, iFT, qk)
		@. qout = LW1 * qr

	elseif solver.method == :RQM4
		# Full step
		@. qr = LW1 * qin
		mul!(qk, FT, qr)
		@. qk *= LD1
		mul!(qr, iFT, qk)
		@. q1 = LW1 * qr

		# Two half-steps
		@. qr = LW2 * qin
		mul!(qk, FT, qr)
		@. qk *= LD2
		mul!(qr, iFT, qk)
		@. qr *= LW1
		mul!(qk, FT, qr)
		@. qk *= LD2
		mul!(qr, iFT, qk)
		@. q2 = LW2 * qr

		# Richardson extrapolation
		@. qout = (4.0*q2 - q1)/3.0
	end

	return nothing
end


text:
## Graph iterators ##

get_node_id(::Integer, ptr::NodePtr, ::Bool, ::Type{R}) where {R<:DirParam} =
        (@_inline_meta; ptr.node)
get_edge_id(::Integer, ptr::NodePtr, ::Bool, ::Type{R}) where {R<:DirParam} =
        (@_inline_meta; ptr.id)
get_edge(n::Integer, ptr::NodePtr, isdir::Bool, ::Type{Forward}) =
        (@_inline_meta; Edge(n, ptr.node, ptr.id, isdir))
get_edge(n::Integer, ptr::NodePtr, isdir::Bool, ::Type{Reverse}) =
        (@_inline_meta; Edge(ptr.node, n, ptr.id, isdir))

# TODO: refactor Edge system for efficiency
# TODO: 1. create edges with native types (N,E)
# TODO: 2. add getters source() and target() that convert ids to Int
# TODO: 3. throughout Junet, convert e.source to source(e)

"""
    PtrView(g::Graph, i::Integer, ::Type{R}, fun::F)

A view used to browse through `NodePtr` objects associated with node `i`.
`R` shows the requested direction: `Forward` for outgoing ones, `Reverse`—for incoming,
and `Both` for both. Each `NodePtr` is then fed into the function `fun`, for which
there are 3 choices: `get_node_id`, `get_edge_id`, and `get_edge`.

# Examples

To get outgoing edges of node `1`, run:
```julia
    PtrView(g, 1, Forward, get_edge)
```

Similarly, for in-neighbors:
```julia
    PtrView(g, 1, Reverse, get_node_id)
```

# Type parameters

This view class is very flexible and can be used to produce 27 different kinds of views.
Under the hood, it uses multiple dispatch to construct and run the correct methods.
It relies on a number of parameters:

* `O` — shows which vector of `NodePtrs` associated with a node to use
* `D` — same as `D` in `Graph`, signifies directedness
* `R` — same as `R` passed in the constructor
* `F` — type of function used to transform the elements
* `N` and `E` are the graph's node and edge ID types

"""
struct PtrView{O<:DirParam,D<:DirParam,R<:DirParam,F,N,E}
    src  :: N
    node :: Node{N,E}
    fun  :: F
end

# const NodeIDView = PtrView{O,D,R,<:typeof(get_node_id)} where {O,D,R}
# const EdgeIDView = PtrView{O,D,R,<:typeof(get_edge_id)} where {O,D,R}
# const EdgeView = PtrView{O,D,R,<:typeof(get_edge_id)} where {O,D,R}

PtrView(g::Graph{N,E,D}, i::Integer, ::Type{R}, fun::F) where {N,E,D,R<:DirParam,F<:Function} =
            PtrView{dir_xor(D,R),D,R,F,N,E}(N(i), g.nodes[i], fun)

eltype(::Type{PtrView{O,D,R,F,N,E}}) where {O,D,R,F<:typeof(get_node_id),N,E} = N
eltype(::Type{PtrView{O,D,R,F,N,E}}) where {O,D,R,F<:typeof(get_edge_id),N,E} = E
eltype(::Type{PtrView{O,D,R,F,N,E}}) where {O,D,R,F<:typeof(get_edge),N,E} = Edge{N,E}
IteratorEltype(::Type{PtrView}) = HasEltype()

length(x::PtrView{O}) where {O<:Forward} = length(x.node.forward)
length(x::PtrView{O}) where {O<:Reverse} = length(x.node.reverse)
length(x::PtrView{O}) where {O<:Both} = length(x.node.forward) + length(x.node.reverse)
IteratorSize(::Type{PtrView}) = HasLength()

function iterate(x::PtrView, i=1)
    @_inline_meta
    if i < length(x)
        return x[i], i + 1
    end
end

getindex(x::PtrView{O,D,R}, i::Integer) where {O<:Forward,D,R} =
    (@_inline_meta; @_propagate_inbounds_meta; x.fun(x.src, x.node.forward[i], true, R))
getindex(x::PtrView{O,D,R}, i::Integer) where {O<:Reverse,D,R} =
    (@_inline_meta; @_propagate_inbounds_meta; x.fun(x.src, x.node.reverse[i], true, R))
function getindex(x::PtrView{O,D,R}, i::Integer) where {O<:Both,D<:Both,R}  # R <: Forward | Both
    @_inline_meta; @_propagate_inbounds_meta
    let l = length(x.node.reverse)
        if i <= l
            return x.fun(x.src, x.node.reverse[i], false, Forward)
        else
            return x.fun(x.src, x.node.forward[i - l], false, Forward)
        end
    end
end
function getindex(x::PtrView{O,D,R}, i::Integer) where {O<:Both,D<:Both,R<:Reverse}  # D <: Both
    @_inline_meta; @_propagate_inbounds_meta
    let l = length(x.node.reverse)
        if i <= l
            return x.fun(x.src, x.node.reverse[i], false, Reverse)
        else
            return x.fun(x.src, x.node.forward[i - l], false, Reverse)
        end
    end
end
function getindex(x::PtrView{O,D}, i::Integer) where {O<:Both,D<:Forward}  # R <: Both
    @_inline_meta; @_propagate_inbounds_meta
    @inbounds let l = length(x.node.reverse)
        if i <= l
            return x.fun(x.src, x.node.reverse[i], true, Reverse)
        else
            return x.fun(x.src, x.node.forward[i - l], true, Forward)
        end
    end
end
function getindex(x::PtrView{O,D}, i::Integer) where {O<:Both,D<:Reverse}  # R <: Both
    @_inline_meta; @_propagate_inbounds_meta
    let l = length(x.node.forward)
        if i <= l
            return x.fun(x.src, x.node.forward[i], true, Reverse)
        else
            return x.fun(x.src, x.node.reverse[i - l], true, Forward)
        end
    end
end

# TODO: implement custom bounds checking / strip extra bounds checks
# TODO: make it a child of AbstractArray, overload the necessary methods to make it work properly
# TODO: play around with lambdas or parent-scope counter to shake off the extra comparisons

"""
    EdgeIter(g::Graph)

Iterator over all edges in the graph.

Under the hood, it uses a 3-tuple of type {N, Int, Bool} to hold the state
during the iteration.
"""
struct EdgeIter{N,E,D}
    graph :: Graph{N,E,D}
end

eltype(::Type{EdgeIter{N,E,D}}) where {N,E,D} = Edge{N,E}
IteratorEltype(::Type{EdgeIter}) = HasEltype()

length(x::EdgeIter) = edgecount(x.graph)
IteratorSize(::Type{EdgeIter}) = HasLength()

function iterate(x::EdgeIter{N,E,D}, state=(1, 1)) where {N,E,D}
    node_ind, ptr_ind = state
    @inbounds while true
        ptrs = fwd_ptrs(x.graph.nodes[node_ind], D)
        if ptr_ind <= length(ptrs)
            edge = get_edge(N(node_ind), ptrs[ptr_ind], D<:Directed, Forward)
            return edge, (node_ind, ptr_ind + 1)
        else
            node_ind += 1
            if node_ind > nodecount(x.graph)
                return nothing
            end
            ptr_ind = 1
        end
    end
end


text:
# Test that example workflow from documentation works

module TestExamples

using OpenStreetMap
import Winston

const MAP_FILENAME = "tech_square.osm"

# Read data from an OSM XML file
nodesLLA, highways, buildings, features = getOSMData( MAP_FILENAME )
println("Number of nodes: $(length(nodesLLA))")
println("Number of highways: $(length(highways))")
println("Number of buildings: $(length(buildings))")
println("Number of features: $(length(features))")

# Get bounds from OSM file
# (This is not always available, depending on how OSM data was exported.)
boundsLLA = getBounds(parseMapXML(MAP_FILENAME))

# Convert to ENU coordinates
lla_reference = center(boundsLLA) # Manual reference point for coordinate transform (optional)
nodes = ENU( nodesLLA, lla_reference )
bounds = ENU( boundsLLA, lla_reference )

# Crop map to bounds
cropMap!(nodes, bounds, highways=highways, buildings=buildings, features=features, delete_nodes=false)

# Extract highway classes (note that OpenStreetMap calls all paths “highways”
roads = roadways(highways)
peds = walkways(highways)
cycles = cycleways(highways)
bldg_classes = classify(buildings)
feat_classes = classify(features)

# Find all highway intersections
intersections = findIntersections(highways)

# Segment only specific levels of roadways
# (e.g., freeways through residential streets, levels 1-6)
segments = segmentHighways(nodes, highways, intersections, roads, Set(1:6))

# Create the routing network
network = createGraph(segments, intersections)

# Compute the shortest and fastest routes from point A to B
loc_start = ENU(-5000, 5500, 0)
loc_end = ENU(5500, -4000, 0)

node0 = nearestNode(nodes, loc_start, network)
node1 = nearestNode(nodes, loc_end, network)
shortest_route, shortest_distance = shortestRoute(network, node0, node1)

fastest_route, fastest_time = fastestRoute(network, node0, node1)
fastest_distance = distance(nodes, fastest_route)

println("Shortest route: $(shortest_distance) m  (Nodes: $(length(shortest_route)))")
println("Fastest route: $(fastest_distance) m  Time: $(fastest_time/60) min  (Nodes: $(length(fastest_route)))")

# Display the shortest and fastest routes
fignum_shortest = plotMap(nodes, highways=highways, bounds=bounds, roadways=roads, route=shortest_route)

fignum_fastest = plotMap(nodes, highways=highways, bounds=bounds, roadways=roads, route=fastest_route)

# Extract Nodes near to (within range) our route's starting location:
loc0 = nodes[node0]
filteredENU = filter((k,v)->haskey(network.v,k), nodes)
local_indices = nodesWithinRange(filteredENU, loc0, 100.0)

# Identify Driving Catchment Areas (within limit):
start_index = nearestNode(filteredENU, loc0)
node_indices, distances = nodesWithinDrivingDistance(network, local_indices, 300.0)

node_indices, distances = nodesWithinDrivingTime(network, local_indices, 50.0)

fignum = plotMap(nodes,
                 highways=highways,
                 buildings=buildings,
                 features=features,
                 bounds=bounds,
                 width=500,
                 feature_classes=feat_classes,
                 building_classes=bldg_classes,
                 roadways=roads)

Winston.savefig("osm_map.png")

end # module TestExamples


text:
"""
    FPlane{FT} <: AbstractRotation

A parameter object for constant rotation around a vertical axis.
"""
struct FPlane{FT} <: AbstractRotation
    f :: FT
end

"""
    FPlane([FT=Float64;] f=nothing, rotation_rate=Ω_Earth, latitude=nothing)

Returns a parameter object for constant rotation at the angular frequency
`f/2`, and therefore with background vorticity `f`, around a vertical axis.
If `f` is not specified, it is calculated from `rotation_rate` and
`latitude` according to the relation `f = 2*rotation_rate*sind(latitude).

By default, `rotation_rate` is assumed to be Earth's.

Also called `FPlane`, after the "f-plane" approximation for the local effect of
a planet's rotation in a planar coordinate system tangent to the planet's surface.
"""
function FPlane(FT::DataType=Float64; f=nothing, rotation_rate=Ω_Earth, latitude=nothing)

    use_f = !isnothing(f)
    use_planet_parameters = !isnothing(latitude)

    if !xor(use_f, use_planet_parameters)
        throw(ArgumentError("Either both keywords rotation_rate and latitude must be " *
                            "specified, *or* only f must be specified."))
    end

    if use_f
        return FPlane{FT}(f)
    elseif use_planet_parameters
        return FPlane{FT}(2rotation_rate*sind(latitude))
    end
end

@inline x_f_cross_U(i, j, k, grid, coriolis::FPlane, U) = - coriolis.f * ℑxyᶠᶜᵃ(i, j, k, grid, U[2])
@inline y_f_cross_U(i, j, k, grid, coriolis::FPlane, U) =   coriolis.f * ℑxyᶜᶠᵃ(i, j, k, grid, U[1])

@inline z_f_cross_U(i, j, k, grid::AbstractGrid{FT}, coriolis::FPlane, U) where FT = zero(FT)

Base.show(io::IO, f_plane::FPlane{FT}) where FT =
    print(io, "FPlane{$FT}: f = ", @sprintf("%.2e", f_plane.f))


text:
__precompile__()

module PiecewiseLinearOpt

using JuMP
import MathOptInterface
const MOI = MathOptInterface
using LinearAlgebra
using Random

export PWLFunction, UnivariatePWLFunction, BivariatePWLFunction, piecewiselinear

include("types.jl")
include("jump.jl")

end # module


text:
""" Functions that solve 1OKP problems """

# All functions are public

" struct to give to the Combo C library "
struct Combo_item
    p::Clonglong
    w::Clonglong
    x::Cint
    i::Cint
end

"""
Returns the optimal solution for the scalarized exact problem computed by the Combo solver.

`λ` is the scalarization vector.

`assignment` is the current assignment to take into account.
"""
function SOLVE1OKP_combo(prob::BiOKP, λ::Vector{Float64}, assignment::Assignment)

    # Handling extreme senarios
	if prob.nbVar == assignment.lastAssigned # no free variables
        returnedSol = Sol(Float64.(assignment.assign), assignment.profit, assignment.weight, true)

        CONFIG.debug && DEBUG_feasibleBinarySolution(prob,returnedSol)

		return returnedSol
	else
        # check if we are in the case where all the objects can be added or if no object can be added at all

		objectsCanBeAdded = false
        iter = assignment.lastAssigned+1
		while !objectsCanBeAdded && iter <= prob.nbVar
			objectsCanBeAdded = objectsCanBeAdded || prob.weights[iter] <= (prob.maxWeight - assignment.weight)
            iter += 1
		end

		if !objectsCanBeAdded # we can't add anymore objects

            sol = Sol(append!(Float64.(assignment.assign[1:assignment.lastAssigned]), zeros(Float64, prob.nbVar-assignment.lastAssigned)), assignment.profit, assignment.weight, true)

            CONFIG.debug && DEBUG_feasibleBinarySolution(prob,sol)

			return sol

		elseif sum(prob.weights[(assignment.lastAssigned+1):end]) <= (prob.maxWeight - assignment.weight) # we can add all the objects

            sol = Sol(append!(Float64.(assignment.assign[1:assignment.lastAssigned]), ones(Float64, prob.nbVar-assignment.lastAssigned)), assignment.profit + [sum(prob.profits[iter, (assignment.lastAssigned+1):end]) for iter=1:prob.nbObj], assignment.weight + sum(prob.weights[(assignment.lastAssigned+1):end]), true)

            CONFIG.debug && DEBUG_feasibleBinarySolution(prob,sol)

			return sol

		end
	end

    weightedProfits = Vector{Float64}(undef,prob.nbVar)
    for var in 1:prob.nbVar
        weightedProfits[var] = (λ[1]*prob.profits[1,var] + λ[2]*prob.profits[2,var])
    end

	items = [Combo_item(weightedProfits[iter], prob.weights[iter], 0, iter) for iter = (assignment.lastAssigned+1):prob.nbVar]

	z = ccall((:solve, comboPath),Clonglong,(Ref{Combo_item}, Cint, Clonglong, Clonglong, Clonglong),items, (prob.nbVar - assignment.lastAssigned), (prob.maxWeight - assignment.weight), 0, 0)

    if z == 0
        z = dot(weightedProfits[(assignment.lastAssigned+1):end], broadcast(it->it.x, items))
    end

    x = Float64.(append!(assignment.assign[1:assignment.lastAssigned], falses(prob.nbVar-assignment.lastAssigned)))
    for it in items
    	x[it.i] = Float64(it.x)
    end

    returnedSol = Sol(x, assignment.profit + [sum(x[(assignment.lastAssigned+1):end] .* prob.profits[iter, (assignment.lastAssigned+1):end]) for iter = 1:prob.nbObj],assignment.weight + sum(prob.weights[(assignment.lastAssigned+1):end] .* x[(assignment.lastAssigned+1):end]), true)

    CONFIG.debug && DEBUG_feasibleBinarySolution(prob,returnedSol)

    return returnedSol
end

"""
Returns the optimal solution for the scalarized and relaxed problem (continuous relaxation).

`λ` is the scalarization vector.

`assignment` is the current assignment to take into account.
"""
function SOLVE1OKP_linear(prob::BiOKP, λ::Vector{Float64}, assignment::Assignment)

    utilities = Vector{Float64}(undef,prob.nbVar-assignment.lastAssigned)

    for i in 1:prob.nbVar-assignment.lastAssigned
        utilities[i] = (λ[1] * prob.profits[1,i+assignment.lastAssigned] + λ[2] * prob.profits[2,i+assignment.lastAssigned]) ./ prob.weights[i+assignment.lastAssigned]
    end

    permList = sortperm(utilities, rev = true)

	nbUnassignedVars = prob.nbVar-assignment.lastAssigned

	weight = assignment.weight
	profit = assignment.profit[1:end]

    x = Float64.(assignment.assign[1:end])

    for iter = 1:nbUnassignedVars
        permList[iter] += assignment.lastAssigned
        x[assignment.lastAssigned + iter] = 0.
    end

    iter = 1
    isBinary = true
    while iter <= nbUnassignedVars && weight + prob.weights[permList[iter]] <= prob.maxWeight
        weight += prob.weights[permList[iter]]
        x[permList[iter]] = 1.
        profit[1] += prob.profits[1,permList[iter]]
        profit[2] += prob.profits[2,permList[iter]]

        iter += 1
    end

    if iter <= nbUnassignedVars && (prob.maxWeight - weight) > 0.
        x[permList[iter]] = (prob.maxWeight - weight) / prob.weights[permList[iter]]
        profit += x[permList[iter]] * prob.profits[1:end, permList[iter]]
        weight += x[permList[iter]] * prob.weights[permList[iter]]
        isBinary = false
    end

    sol = Sol(x, profit, weight, isBinary)

    CONFIG.debug && DEBUG_feasibleSolution(prob,sol)

	return sol
end

"""
Returns the optimal solution for the scalarized exact problem computed by the solver GLPK (using JuMP).

`λ` is the scalarization vector.

`assignment` is the current assignment to take into account.
"""
function SOLVE1OKP_GLPK(prob::BiOKP, λ::Vector{Float64}, assignment::Assignment)

	model = Model(GLPK.Optimizer)
	x = @variable(model, x[1:(prob.nbVar-assignment.lastAssigned)], Bin)
	@constraint(model, Weights, sum(x .* prob.weights[(assignment.lastAssigned+1):end]) + assignment.weight <= prob.maxWeight)
	@objective(model, Max, λ[1] * sum(x .* prob.profits[1, (assignment.lastAssigned+1):end]) + λ[2] * sum(x .* prob.profits[2, (assignment.lastAssigned+1):end]))

	optimize!(model)

	X = append!(assignment.assign[1:assignment.lastAssigned],(Float64).(value.(x)))

	termStatus = termination_status(model)

	if termStatus == MOI.OPTIMAL
        sol = Sol(X, [sum(X .* prob.profits[1,1:end]), sum(X .* prob.profits[2,1:end])], sum(X .* prob.weights), true)

        CONFIG.debug && DEBUG_feasibleBinarySolution(prob,sol)

        return sol

	elseif termStatus == MOI.INFEASIBLE
        error("GLPK is stating that the problem is infeasible")
	else
		error("MOI or GLPK aren't working as they should be")
	end
end


text:
using HOODESolver: filtred_f, PrepareU0
using Test

include("henon_heiles.jl")
include("dataprep_u0.jl")

# these 3 functions come from old version with sin and cos computation
# to test the new version
function reffltfct( u, s, c )
    a = c*u[1] + s*u[3]
    b = 2 * u[2] * a
    return [ s*b, u[4], -c*b, -u[2] - a^2 + u[2]^2, 1 ]
end
function reffltfctgen2( u::Array{Array{T,1},1}, s, c, n_tau ) where T <: Number
    f = reffltfct.( u, s, c)
    # convert an array of vector to a matrix
    f = collect(transpose(reshape(collect(Iterators.flatten(f)), 5, n_tau)))
    return f
end
function reffltfctgen( uMat::Array{T,2}, s, c, n_tau ) where T <: Number
    return reffltfctgen2(reshape(mapslices(x->[x],uMat,dims = (2,)), n_tau), s, c, n_tau)
end

function test_tau_A()
    epsilon = big"0.5"
    n_tau = 64
    A = [0 0 1 0; 0 0 0 0; -1 0 0 0; 0 0 0 0]
    parphi = HOODESolver.PreparePhi(epsilon, n_tau, A, henon_heiles)
    tol= 1e-77
    @testset "preparephi tau_A" begin
        for i = 1:n_tau
            s,c = setprecision(precision(BigFloat)+32) do
                v = (i-1)*2big(pi)/n_tau
                sin(v), cos(v)
            end

            res_ref = one(zeros(BigFloat,5,5))
            res_ref[1,1] = res_ref[3,3] = c
            res_ref[1,3] = s
            res_ref[3,1]= -res_ref[1,3]
            @test isapprox(res_ref, convert(Matrix, parphi.tau_Ap[i]), atol=tol)
            res_ref[3,1], res_ref[1,3] = res_ref[1,3], res_ref[3,1]
            @test isapprox(res_ref, convert(Matrix, parphi.tau_Ap_inv[i]), atol=tol)
        end
    end
end
function test_fct()
    epsilon = big"0.5"
    n_tau = 64
    A = [0 0 1 0; 0 0 0 0; -1 0 0 0; 0 0 0 0]
    parphi = HOODESolver.PreparePhi(epsilon, n_tau, A, henon_heiles)
    tol= 1e-76
    @testset "preparephi fct" begin
        for i = 1:n_tau
            s,c = setprecision(precision(BigFloat)+32) do
                v = (i-1)*2big(pi)/n_tau
                sin(v), cos(v)
            end
            u0 = rand(BigFloat,5)
            res_ref = reffltfct(u0, s, c)
            res = filtred_f(u0,parphi.tau_Ap_inv[i], parphi.tau_Ap[i], parphi)
            @test isapprox(res_ref, res, atol=tol)
        end
    end
end

function testpreparephi0()

    epsilon = big"0.5"
    n_tau = 64
    tol=1e-76
    @time @testset "preparephi and function for all tau" begin
        parphi = HOODESolver.PreparePhi(
    epsilon,
    n_tau,
    [0 0 1 0; 0 0 0 0; -1 0 0 0; 0 0 0 0],
    henon_heiles
)
        u_mat = rand(BigFloat, 5, n_tau)
        u_mat_tr = collect(transpose(u_mat))
        tau = LinRange(zero(BigFloat), one(BigFloat), n_tau + 1)[1:end - 1]
        s = sinpi.(2tau)
        c = cospi.(2tau)
        old_m = collect(transpose(reffltfctgen(u_mat_tr, s, c, n_tau)))
        new_m = HOODESolver.filtredfct(parphi,u_mat)
        println("norm=$(norm(old_m-new_m))")
        @test isapprox(old_m, new_m, atol=tol)
    end
end

function testpreparephi()
    @time @testset "phi and prepare u0" begin
        tab_ref, tab_eps = setprecision(1024) do
            get_prepare_u0_for_test()
        end
        for i_prec=1:3
            prec = 256*2^(i_prec-1)
            tol=max((1e-77)^i_prec, 1e-200)
            setprecision(prec) do
                for i_eps=1:size(tab_eps,1)
                   epsilon = BigFloat(tab_eps[i_eps])
                   n_tau = 32
                   parphi = HOODESolver.PreparePhi(
    epsilon,
    n_tau,
    [0 0 1 0; 0 0 0 0; -1 0 0 0; 0 0 0 0],
    henon_heiles
)
                    u0 = BigFloat.([0.125, 0.140625, 0.15625, 0.171875])
                    for ord=2:10
            #     paru0 = PrepareU0(parphi, ord, u0, 1024)
                        paru0 = PrepareU0(parphi, ord, u0)
 #                   println("prec=$prec ord=$ord norm=$(norm(tab_ref[:, :, ord]- paru0.ut0))")
                        @test isapprox(tab_ref[:, :, ord, i_eps], paru0.ut0[1:4,:], atol=tol, rtol=tol)
                    end
                end
            end
        end
    end
end

test_tau_A()
test_fct()
testpreparephi0()
testpreparephi()


text:
using MPIFiles, Test

@testset "Testing DatasetStore submodule" begin

# Bruker store

storeA = BrukerDatasetStore(joinpath(datadir, "BrukerStore"))
@test readonly(storeA) == true

studiesA = getStudies(storeA)
@test [s.name for s in studiesA] == ["SF", "Wuerfelphantom", "Size Distributions", "EasyAxisContrast"]

# MDF store

storeB = MDFDatasetStore(joinpath(tmpdir, "MDFStoreA"))
@test readonly(storeB) == false
empty!(storeB)

# First add studies and remove them
studyA = Study(storeB,"StudyA")
studyB = Study(storeB,"StudyB")

addStudy(storeB, studyA)
addStudy(storeB, studyB)

studies = getStudies(storeB)
@test studies == [studyA, studyB]

remove(studyB)
studies = getStudies(storeB)
@test studies == [studyA]

empty!(storeB)

# Now lets copy over the Bruker studies
exportData(storeA, storeB, keepExpNum=true)
storeC = MDFDatasetStore(joinpath(tmpdir, "MDFStoreB"))
empty!(storeC)
exportData(storeB, storeC, SNRThresh=4.0)
createArtifact(storeB, "https://")

@test validate(storeA)
@test validate(storeB)
@test validate(storeC)

studiesB = getStudies(storeB)
studiesC = getStudies(storeC)
@test getExperiments(studiesB[2])[1].num == 18
@test getExperiments(studiesC[2])[1].num == 1

# export into new/existing Study
newStudy = Study(storeB, "NewStudy")
exportData(studiesB[2], storeB, newStudy)
@info [s.name for s in getStudies(storeB)]
@test getExperiments(newStudy)[1].num == 1

# Next we test changing parameters and study names.
s = studiesC[2]
@test validate(s)
changeStudy(s, "NewStudyName")
studiesD = getStudies(storeC)
sNew = studiesD[end]
@test validate(sNew)
@test sNew.name == "NewStudyName"
enforceStudy(sNew)
@test validate(sNew)

# Experiment handling

#exps = getExperiments(storeB, studyA)
#@test exps == Experiment[]
#@info getNewExperimentNum(storeB,studyA) == 1

# Calibration handling
@test getNewCalibNum(storeB) == 3

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
    AbstractChannel{T}

Representation of a channel passing objects of type `T`.
"""
abstract type AbstractChannel{T} end

push!(c::AbstractChannel, v) = (put!(c, v); c)
popfirst!(c::AbstractChannel) = take!(c)

"""
    Channel{T=Any}(size::Int=0)

Constructs a `Channel` with an internal buffer that can hold a maximum of `size` objects
of type `T`.
[`put!`](@ref) calls on a full channel block until an object is removed with [`take!`](@ref).

`Channel(0)` constructs an unbuffered channel. `put!` blocks until a matching `take!` is called.
And vice-versa.

Other constructors:

* `Channel()`: default constructor, equivalent to `Channel{Any}(0)`
* `Channel(Inf)`: equivalent to `Channel{Any}(typemax(Int))`
* `Channel(sz)`: equivalent to `Channel{Any}(sz)`

!!! compat "Julia 1.3"
    The default constructor `Channel()` and default `size=0` were added in Julia 1.3.
"""
mutable struct Channel{T} <: AbstractChannel{T}
    cond_take::Threads.Condition                 # waiting for data to become available
    cond_wait::Threads.Condition                 # waiting for data to become maybe available
    cond_put::Threads.Condition                  # waiting for a writeable slot
    @atomic state::Symbol
    excp::Union{Exception, Nothing}      # exception to be thrown when state !== :open

    data::Vector{T}
    @atomic n_avail_items::Int           # Available items for taking, can be read without lock
    sz_max::Int                          # maximum size of channel

    function Channel{T}(sz::Integer = 0) where T
        if sz < 0
            throw(ArgumentError("Channel size must be either 0, a positive integer or Inf"))
        end
        lock = ReentrantLock()
        cond_put, cond_take = Threads.Condition(lock), Threads.Condition(lock)
        cond_wait = (sz == 0 ? Threads.Condition(lock) : cond_take) # wait is distinct from take iff unbuffered
        return new(cond_take, cond_wait, cond_put, :open, nothing, Vector{T}(), 0, sz)
    end
end

function Channel{T}(sz::Float64) where T
    sz = (sz == Inf ? typemax(Int) : convert(Int, sz))
    return Channel{T}(sz)
end
Channel(sz=0) = Channel{Any}(sz)

# special constructors
"""
    Channel{T=Any}(func::Function, size=0; taskref=nothing, spawn=false)

Create a new task from `func`, bind it to a new channel of type
`T` and size `size`, and schedule the task, all in a single call.

`func` must accept the bound channel as its only argument.

If you need a reference to the created task, pass a `Ref{Task}` object via
the keyword argument `taskref`.

If `spawn = true`, the Task created for `func` may be scheduled on another thread
in parallel, equivalent to creating a task via [`Threads.@spawn`](@ref).

Return a `Channel`.

# Examples
```jldoctest
julia> chnl = Channel() do ch
           foreach(i -> put!(ch, i), 1:4)
       end;

julia> typeof(chnl)
Channel{Any}

julia> for i in chnl
           @show i
       end;
i = 1
i = 2
i = 3
i = 4
```

Referencing the created task:

```jldoctest
julia> taskref = Ref{Task}();

julia> chnl = Channel(taskref=taskref) do ch
           println(take!(ch))
       end;

julia> istaskdone(taskref[])
false

julia> put!(chnl, "Hello");
Hello

julia> istaskdone(taskref[])
true
```

!!! compat "Julia 1.3"
    The `spawn=` parameter was added in Julia 1.3. This constructor was added in Julia 1.3.
    In earlier versions of Julia, Channel used keyword arguments to set `size` and `T`, but
    those constructors are deprecated.

```jldoctest
julia> chnl = Channel{Char}(1, spawn=true) do ch
           for c in "hello world"
               put!(ch, c)
           end
       end
Channel{Char}(1) (2 items available)

julia> String(collect(chnl))
"hello world"
```
"""
function Channel{T}(func::Function, size=0; taskref=nothing, spawn=false) where T
    chnl = Channel{T}(size)
    task = Task(() -> func(chnl))
    task.sticky = !spawn
    bind(chnl, task)
    if spawn
        schedule(task) # start it on (potentially) another thread
    else
        yield(task) # immediately start it, yielding the current thread
    end
    isa(taskref, Ref{Task}) && (taskref[] = task)
    return chnl
end
Channel(func::Function, args...; kwargs...) = Channel{Any}(func, args...; kwargs...)

# This constructor is deprecated as of Julia v1.3, and should not be used.
# (Note that this constructor also matches `Channel(::Function)` w/out any kwargs, which is
# of course not deprecated.)
# We use `nothing` default values to check which arguments were set in order to throw the
# deprecation warning if users try to use `spawn=` with `ctype=` or `csize=`.
function Channel(func::Function; ctype=nothing, csize=nothing, taskref=nothing, spawn=nothing)
    # The spawn= keyword argument was added in Julia v1.3, and cannot be used with the
    # deprecated keyword arguments `ctype=` or `csize=`.
    if (ctype !== nothing || csize !== nothing) && spawn !== nothing
        throw(ArgumentError("Cannot set `spawn=` in the deprecated constructor `Channel(f; ctype=Any, csize=0)`. Please use `Channel{T=Any}(f, size=0; taskref=nothing, spawn=false)` instead!"))
    end
    # Set the actual default values for the arguments.
    ctype === nothing && (ctype = Any)
    csize === nothing && (csize = 0)
    spawn === nothing && (spawn = false)
    return Channel{ctype}(func, csize; taskref=taskref, spawn=spawn)
end

closed_exception() = InvalidStateException("Channel is closed.", :closed)

isbuffered(c::Channel) = c.sz_max==0 ? false : true

function check_channel_state(c::Channel)
    if !isopen(c)
        # if the monotonic load succeed, now do an acquire fence
        (@atomic :acquire c.state) === :open && concurrency_violation()
        excp = c.excp
        excp !== nothing && throw(excp)
        throw(closed_exception())
    end
end
"""
    close(c::Channel[, excp::Exception])

Close a channel. An exception (optionally given by `excp`), is thrown by:

* [`put!`](@ref) on a closed channel.
* [`take!`](@ref) and [`fetch`](@ref) on an empty, closed channel.
"""
function close(c::Channel, excp::Exception=closed_exception())
    lock(c)
    try
        c.excp = excp
        @atomic :release c.state = :closed
        notify_error(c.cond_take, excp)
        notify_error(c.cond_wait, excp)
        notify_error(c.cond_put, excp)
    finally
        unlock(c)
    end
    nothing
end
isopen(c::Channel) = ((@atomic :monotonic c.state) === :open)

"""
    bind(chnl::Channel, task::Task)

Associate the lifetime of `chnl` with a task.
`Channel` `chnl` is automatically closed when the task terminates.
Any uncaught exception in the task is propagated to all waiters on `chnl`.

The `chnl` object can be explicitly closed independent of task termination.
Terminating tasks have no effect on already closed `Channel` objects.

When a channel is bound to multiple tasks, the first task to terminate will
close the channel. When multiple channels are bound to the same task,
termination of the task will close all of the bound channels.

# Examples
```jldoctest
julia> c = Channel(0);

julia> task = @async foreach(i->put!(c, i), 1:4);

julia> bind(c,task);

julia> for i in c
           @show i
       end;
i = 1
i = 2
i = 3
i = 4

julia> isopen(c)
false
```

```jldoctest
julia> c = Channel(0);

julia> task = @async (put!(c, 1); error("foo"));

julia> bind(c, task);

julia> take!(c)
1

julia> put!(c, 1);
ERROR: TaskFailedException
Stacktrace:
[...]
    nested task error: foo
[...]
```
"""
function bind(c::Channel, task::Task)
    T = Task(() -> close_chnl_on_taskdone(task, c))
    _wait2(task, T)
    return c
end

"""
    channeled_tasks(n::Int, funcs...; ctypes=fill(Any,n), csizes=fill(0,n))

A convenience method to create `n` channels and bind them to tasks started
from the provided functions in a single call. Each `func` must accept `n` arguments
which are the created channels. Channel types and sizes may be specified via
keyword arguments `ctypes` and `csizes` respectively. If unspecified, all channels are
of type `Channel{Any}(0)`.

Returns a tuple, `(Array{Channel}, Array{Task})`, of the created channels and tasks.
"""
function channeled_tasks(n::Int, funcs...; ctypes=fill(Any,n), csizes=fill(0,n))
    @assert length(csizes) == n
    @assert length(ctypes) == n

    chnls = map(i -> Channel{ctypes[i]}(csizes[i]), 1:n)
    tasks = Task[ Task(() -> f(chnls...)) for f in funcs ]

    # bind all tasks to all channels and schedule them
    foreach(t -> foreach(c -> bind(c, t), chnls), tasks)
    foreach(schedule, tasks)
    yield() # Allow scheduled tasks to run

    return (chnls, tasks)
end

function close_chnl_on_taskdone(t::Task, c::Channel)
    isopen(c) || return
    lock(c)
    try
        isopen(c) || return
        if istaskfailed(t)
            close(c, TaskFailedException(t))
            return
        end
        close(c)
    finally
        unlock(c)
    end
    nothing
end

struct InvalidStateException <: Exception
    msg::String
    state::Symbol
end
showerror(io::IO, ex::InvalidStateException) = print(io, "InvalidStateException: ", ex.msg)

"""
    put!(c::Channel, v)

Append an item `v` to the channel `c`. Blocks if the channel is full.

For unbuffered channels, blocks until a [`take!`](@ref) is performed by a different
task.

!!! compat "Julia 1.1"
    `v` now gets converted to the channel's type with [`convert`](@ref) as `put!` is called.
"""
function put!(c::Channel{T}, v) where T
    check_channel_state(c)
    v = convert(T, v)
    return isbuffered(c) ? put_buffered(c, v) : put_unbuffered(c, v)
end

# Atomically update channel n_avail, *assuming* we hold the channel lock.
function _increment_n_avail(c, inc)
    # We hold the channel lock so it's safe to non-atomically read and
    # increment c.n_avail_items
    newlen = c.n_avail_items + inc
    # Atomically store c.n_avail_items to prevent data races with other threads
    # reading this outside the lock.
    @atomic :monotonic c.n_avail_items = newlen
end

function put_buffered(c::Channel, v)
    lock(c)
    did_buffer = false
    try
        # Increment channel n_avail eagerly (before push!) to count data in the
        # buffer as well as offers from tasks which are blocked in wait().
        _increment_n_avail(c, 1)
        while length(c.data) == c.sz_max
            check_channel_state(c)
            wait(c.cond_put)
        end
        check_channel_state(c)
        push!(c.data, v)
        did_buffer = true
        # notify all, since some of the waiters may be on a "fetch" call.
        notify(c.cond_take, nothing, true, false)
    finally
        # Decrement the available items if this task had an exception before pushing the
        # item to the buffer (e.g., during `wait(c.cond_put)`):
        did_buffer || _increment_n_avail(c, -1)
        unlock(c)
    end
    return v
end

function put_unbuffered(c::Channel, v)
    lock(c)
    taker = try
        _increment_n_avail(c, 1)
        while isempty(c.cond_take.waitq)
            check_channel_state(c)
            notify(c.cond_wait)
            wait(c.cond_put)
        end
        check_channel_state(c)
        # unfair scheduled version of: notify(c.cond_take, v, false, false); yield()
        popfirst!(c.cond_take.waitq)
    finally
        _increment_n_avail(c, -1)
        unlock(c)
    end
    schedule(taker, v)
    yield()  # immediately give taker a chance to run, but don't block the current task
    return v
end

"""
    fetch(c::Channel)

Wait for and get the first available item from the channel. Does not
remove the item. `fetch` is unsupported on an unbuffered (0-size) channel.
"""
fetch(c::Channel) = isbuffered(c) ? fetch_buffered(c) : fetch_unbuffered(c)
function fetch_buffered(c::Channel)
    lock(c)
    try
        while isempty(c.data)
            check_channel_state(c)
            wait(c.cond_take)
        end
        return c.data[1]
    finally
        unlock(c)
    end
end
fetch_unbuffered(c::Channel) = throw(ErrorException("`fetch` is not supported on an unbuffered Channel."))

"""
    take!(c::Channel)

Remove and return a value from a [`Channel`](@ref). Blocks until data is available.

For unbuffered channels, blocks until a [`put!`](@ref) is performed by a different
task.
"""
take!(c::Channel) = isbuffered(c) ? take_buffered(c) : take_unbuffered(c)
function take_buffered(c::Channel)
    lock(c)
    try
        while isempty(c.data)
            check_channel_state(c)
            wait(c.cond_take)
        end
        v = popfirst!(c.data)
        _increment_n_avail(c, -1)
        notify(c.cond_put, nothing, false, false) # notify only one, since only one slot has become available for a put!.
        return v
    finally
        unlock(c)
    end
end

# 0-size channel
function take_unbuffered(c::Channel{T}) where T
    lock(c)
    try
        check_channel_state(c)
        notify(c.cond_put, nothing, false, false)
        return wait(c.cond_take)::T
    finally
        unlock(c)
    end
end

"""
    isready(c::Channel)

Determine whether a [`Channel`](@ref) has a value stored to it. Returns
immediately, does not block.

For unbuffered channels returns `true` if there are tasks waiting
on a [`put!`](@ref).
"""
isready(c::Channel) = n_avail(c) > 0
isempty(c::Channel) = n_avail(c) == 0
function n_avail(c::Channel)
    # Lock-free equivalent to `length(c.data) + length(c.cond_put.waitq)`
    @atomic :monotonic c.n_avail_items
end

lock(c::Channel) = lock(c.cond_take)
lock(f, c::Channel) = lock(f, c.cond_take)
unlock(c::Channel) = unlock(c.cond_take)
trylock(c::Channel) = trylock(c.cond_take)

function wait(c::Channel)
    isready(c) && return
    lock(c)
    try
        while !isready(c)
            check_channel_state(c)
            wait(c.cond_wait)
        end
    finally
        unlock(c)
    end
    nothing
end

eltype(::Type{Channel{T}}) where {T} = T

show(io::IO, c::Channel) = print(io, typeof(c), "(", c.sz_max, ")")

function show(io::IO, ::MIME"text/plain", c::Channel)
    show(io, c)
    if !(get(io, :compact, false)::Bool)
        if !isopen(c)
            print(io, " (closed)")
        else
            n = n_avail(c)
            if n == 0
                print(io, " (empty)")
            else
                s = n == 1 ? "" : "s"
                print(io, " (", n, " item$s available)")
            end
        end
    end
end

function iterate(c::Channel, state=nothing)
    try
        return (take!(c), nothing)
    catch e
        if isa(e, InvalidStateException) && e.state === :closed
            return nothing
        else
            rethrow()
        end
    end
end

IteratorSize(::Type{<:Channel}) = SizeUnknown()

max_stars_repo_name: BSnelling/julia

max_issues_repo_name: BSnelling/julia

max_forks_repo_name: BSnelling/julia


text:
module NoisyQuantumComputerSimulator

export Circuit, exec, reset_state!
export Gate, ID, X, Y, Z, H, PHASE, S, T, CZ, RX, RY, RZ, CNOT, CCNOT, SWAP, CSWAP, CCSWAP, controlled, noisify, damp_amplitude
export damping_residual_kraus, damping_kraus, damping_kraus_map

include("FullRegisterGate.jl")
include("Gates.jl")

using LinearAlgebra
using .FullRegisterGate
using .Gates

mutable struct Circuit
	qubit_count::Int
	density_matrix::Matrix{Complex{Float64}}
	commands::Array{Function}

	# Circuit creates quantum register with the given size and sets it's value to |0⟩.
	function Circuit(size::Integer, gates...)
		state_size = 2^size
		density_matrix = zeros(Complex{Float64}, state_size, state_size)
		density_matrix[1,1] = 1
		circuit = new(Int(size), density_matrix, Array{Function}[])
		for gate ∈ gates
			circuit += gate
		end
		return circuit
	end
end

# reset_state! sets quantum register state to |0⟩.
function reset_state!(c::Circuit)
	c.density_matrix .= complex(0.0, 0.0)
	c.density_matrix[1,1] = 1
end

function add_gate(c::Circuit, g::Gate)::Circuit
	if !isempty(g.control_bit_indexes) && !isempty(g.kraus_operators)
		error("applying Kraus operators within a gate with control bit(s) is not supported; instead, apply Kraus operators to control and target qubits individually")
	end
	for control ∈ g.control_bit_indexes
		if control < 0 || control >= c.qubit_count
			error("control bit at index $control is outside of bounds of quantum register [0, $(c.qubit_count - 1)]")
		end
	end
	if g.swap_bit_index_a ≠ -1 && g.swap_bit_index_b ≠ -1
		if g.swap_bit_index_a < 0 || g.swap_bit_index_a >= c.qubit_count
			error("bit to SWAP at index $(g.swap_bit_index_a) is outside of bounds of quantum register [0, $(c.qubit_count - 1)]")
		end
		if g.swap_bit_index_b < 0 || g.swap_bit_index_b >= c.qubit_count
			error("bit to SWAP at index $(g.swap_bit_index_b) is outside of bounds of quantum register [0, $(c.qubit_count - 1)]")
		end
		# build SWAP as a sequence of 3 CNOTs.
		not_matrix = Matrix{Complex{Float64}}([0 1; 1 0])
		controls_a = vcat(g.swap_bit_index_a, g.control_bit_indexes)
		controls_b = vcat(g.swap_bit_index_b, g.control_bit_indexes)
		add_unitary_gate!(c, not_matrix, ()->(), g.swap_bit_index_a, controls_b)
		add_unitary_gate!(c, not_matrix, ()->(), g.swap_bit_index_b, controls_a)
		add_unitary_gate!(c, not_matrix, ()->(), g.swap_bit_index_a, controls_b)
	else
		add_unitary_gate!(c, g.matrix, g.matrix_func, g.qubit_index, g.control_bit_indexes)
	end
	if size(g.kraus_operators)[1] > 0
		add_noise!(c, g.kraus_operators, g.qubit_index)
	end
	return c
end

# Defining + allows adding gates like this: Circuit += gate
Base.:+(c::Circuit, g::Gate) = add_gate(c, g)

function exec(c::Circuit, params::Dict{String,Float64} = Dict{String,Float64}())::Array{Complex{Float64}}
	for command ∈ c.commands
		command(params)
	end
	return c.density_matrix
end

# add_unitary_gate! adds the given unitary gate with optional control bits to the given Circuit.
function add_unitary_gate!(c::Circuit,
	matrix::AbstractMatrix{Complex{Float64}},
	matrix_func::Function,
	gate_lowest_index::Integer,
	control_bit_indexes::AbstractArray{Int64,1} = Array{Int64,1}([]))

	push!(c.commands, (params::Dict{String,Float64})->begin
		small_matrix = !isempty(matrix) ? matrix : matrix_func(params)
		matrix_size = size(small_matrix)[1]
		if matrix_size == 0
			error("matrix is empty")
		end
		gate_size = Int(log2(matrix_size))
		if gate_lowest_index < 0 || gate_lowest_index + gate_size > c.qubit_count
			error("quantum register of size $(c.qubit_count) cannot fit given gate of size $(gate_size) at index $(gate_lowest_index)")
		end
		big_matrix = expandGateToFullRegister(c.qubit_count, small_matrix, gate_lowest_index, control_bit_indexes)
		c.density_matrix = big_matrix * c.density_matrix * big_matrix'
		return
	end)
	return
end

# add_noise! adds the specified via Kraus operators noise gate to the given Circuit.
function add_noise!(c::Circuit,
	kraus_operators::AbstractVector{Matrix{Complex{Float64}}},
	gate_lowest_index::Integer)

	push!(c.commands, (params::Dict{String,Float64})->begin
		kraus_probabilities = Vector{Float64}([])
		big_kraus_operators = Vector{Matrix{Complex{Float64}}}([])
		for kraus ∈ kraus_operators
			big_kraus = expandGateToFullRegister(c.qubit_count, kraus, gate_lowest_index)
			trace = tr(big_kraus * c.density_matrix * big_kraus')
			if imag(trace) ≉ 0
				error("expected trace be a real value since it's a probability, got $trace")
			end
			kraus_probability = real(tr(big_kraus * c.density_matrix * big_kraus'))
			push!(big_kraus_operators, big_kraus)
			push!(kraus_probabilities, kraus_probability)
		end
		random_kraus_index = random_index(kraus_probabilities)
		big_kraus = big_kraus_operators[random_kraus_index]
		big_kraus_probability = kraus_probabilities[random_kraus_index]
		c.density_matrix = big_kraus * c.density_matrix * big_kraus' / big_kraus_probability
		return
	end)
	return
end

# random_index returns random index according to the given probabilities.
# Min returned value is 1; max return value is the length of the given array of probabilities.
function random_index(probabilities::AbstractVector{Float64})
	random = rand()
	return value_index_by_probabilities(random, probabilities)
end

# value_index_by_probabilities returns index of the given value according to the given probabilities.
# Min returned value is 1; max return value is the length of the given array of probabilities.
function value_index_by_probabilities(value::Float64, probabilities::AbstractVector{Float64})
	total_probability = zero(Float64)
	for (index, probability) ∈ enumerate(probabilities)
		total_probability += probability
		if value < total_probability
			return index
		end
	end
	return length(probabilities) # Unlikely, but control might reach here in case of some rounding issues.
end

end # module


text:
function nextEdges(edges::Array{Int64,1})
	theNexts = Array{Int64,1}()
	n = length(edges)
	for k in 1:n
		push!(theNexts, k==n ? edges[1] : edges[k+1] )
	end
	theNexts
end

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

	end
	forward
end

function cycleBasis(V,EV)
	I = Int64[];
	J = Int64[];
	Val = Int64[];
	for e in 1:size(EV,2)
		i,j = EV[:,e]
		push!(I,i); push!(I,j);
		push!(J,e); push!(J,e);
		push!(Val,1); push!(Val,1);
	end
	S = sparse(J,I,Val)
	rows = rowvals(S)
	VE = Array{Array{Int64,1},1}()
	for k=1:size(S,2)
		# extract column
		col = [rows[h] for h in nzrange(S,k)]
		push!(VE, col)
	end
	forward = edgeSlopeOrdering(VE,V,EV)
	for (v,(edges,nexts)) in enumerate(forward)
		for (h,e) in enumerate(edges)
			S[edges[h],v] = nexts[h]
		end
	end
	return S
end

S = cycleBasis(W,EW)

# return the first index with positive value in Val. If none return zero.
function getFirstPositive(Val)
	for k=1:length(Val)
		if Val[k] > 0  return k end
	end
	return 0
end

# Extract 1-cycles from the sparse matrix EV evaluated on next edges
function edgeCycles(S,V,EV)
	cycles = []
	I,J,Val = findnz(S)
	k = getFirstPositive(Val)
	while k != 0
		v = J[k]
		firstEdge = I[k]
		cycle = [firstEdge]
		next_edge = abs(S[abs(firstEdge),v])
		if S[abs(firstEdge),v] > 0
			S[abs(firstEdge),v] = -S[abs(firstEdge),v]  # sign as visited
		end
		while next_edge != firstEdge
			push!(cycle, next_edge)
			if S[abs(next_edge),v] > 0
				S[abs(next_edge),v] = -S[abs(next_edge),v]
			end
			nextv = pop!(setdiff(Set(EV[abs(next_edge)]), Set(v)))
			next_edge = S[abs(next_edge),nextv]
			v = nextv
		end
		push!(cycles,cycle)
		I,J,Val = findnz(S)
		k = getFirstPositive(Val)
	end
	cycles
end

edgeCycles(S,W,EW)


text:
"""
    sayhi()

A function that says hi.

# Arguments
None.

# Notes
None.

# Examples
```julia
julia> sayhi()
"Hi!"
```
"""
function sayhi()
    return "Hi!"
end


text:
using VLKeggSDK

# test -
test_compound = "C21H29N7O14P2"

# cut -
d = extract_atom_dictionary(test_compound)


text:
using Documenter, QuantExQASM

makedocs(
    modules = [QuantExQASM],
    clean = false,
    sitename = "QuantExQASM.jl",
    pages = Any[
        "Home" => "index.md",
        #"Tutorial" => "tutorial.md",
        "Manual" => Any[
            "Gates" => "gates_circuits/GateOps.md",
            "Circuits" => "gates_circuits/circuits.md",
            "Algorithms" => Any[ "NCU" => "algo/ncu.md", "Grover" => "algo/grover.md" ],
        ],
    ]
)
deploydocs(
    repo = "github.com/ICHEC/QuantExQASM.jl.git",
)


text:
"""
    struct Paraboloid <: Shape

Paraboloid shape model.

"""
struct Paraboloid <: Shape
    id::Int64
    h0::Float64
    h1::Float64
    radius::Float64
    inverted::Bool
    obj_to_world::Transformation
    world_to_obj::Transformation
end

Paraboloid() = Paraboloid(1, 0.0, 1.0, 1.0, false, Transformation(), Transformation())
Paraboloid(T::Transformation) = Paraboloid(1, 0.0, 1.0, 1.0, false, T, inv(T))

"""
    can_intersect(s::Paraboloid) = true

Is a paraboloid intersectable?

TODO: should this function be replaced by an abstract type like `IntersectableShape`?

"""
can_intersect(s::Paraboloid) = true

"""
    obj_bounds(P::Paraboloid)

The bounding box or a paraboloid in object space.

"""
function obj_bounds(P::Paraboloid)
    v0 = Point3(-P.radius, -P.radius, P.h0)
    v1 = Point3(P.radius, P.radius, P.h1)
    return BoundingBox(v0, v1)
end

"""
    world_bounds(P::Paraboloid)

The bounding box or a paraboloid transformed to world space.

"""
function world_bounds(P::Paraboloid)
    P.obj_to_world(obj_bounds(P))
end

"""
    (T::Transformation)(P::Paraboloid)

Apply a `Transformation` to a `Paraboloid` object, return a new object.

"""
function (T::Transformation)(P::Paraboloid)
    T2 = T * P.obj_to_world
    Paraboloid(P.id, P.h0, P.h1, P.radius, P.inverted, T2, inv(T2))
end

function paraboloid_coefs(ray::Ray, Par::Paraboloid)
    A = Par.h1 * (ray.direction.x^2 + ray.direction.y^2)
    B = 2*Par.h1 * (ray.origin.x * ray.direction.x + ray.origin.y * ray.direction.y) - Par.radius^2 * ray.direction.z
    C = Par.h1 * (ray.origin.x^2 + ray.origin.y^2) - Par.radius^2 * ray.origin.z
    return A, B, C
end

function shape_intersect(r::Ray, Par::Paraboloid)
    ray = Par.world_to_obj(r)

    hit, tnear, tfar = quadratic(paraboloid_coefs(ray, Par)...)

    # Not hit or both hits out of bounds
    if !hit || tnear > r.tmax || tfar < r.tmin || (tnear < r.tmin && tfar > r.tmax)
        return nothing, NaN, NaN
    end

    # Compute hit point, check if it's on our paraboloid.
    t = tnear > r.tmin ? tnear : tfar
    P = ray(t)
    if P.z > Par.h1 || P.z < Par.h0
        return nothing, NaN, NaN
    end

    r = sqrt(P.x^2 + P.y^2)

    n = normalize(Normal3(2r, 0.0, -1.0))
    s = normalize(Vector3(2r, 0.0, 1.0))

    n = flip_normal(n, ray.direction)

    DG = DifferentialGeometry(Par.obj_to_world(P), Par.obj_to_world(n), Par.obj_to_world(s))

    return DG, t, 5e-4 * t

end

"""
    intersectP(r::Ray, Par::Paraboloid)

"Quick" true/false intersect test between given ray and paraboloid.

"""
function intersectP(r::Ray, Par::Paraboloid)
    ray = Par.world_to_obj(r)

    hit, tnear, tfar = quadratic(paraboloid_coefs(ray, Par)...)

    # Not hit or both hits out of bounds
    if !hit || tnear > r.tmax || tfar < r.tmin || (tnear < r.tmin && tfar > r.tmax)
        return false
    end

    # Compute hit point, check if it's on our paraboloid.
    t = tnear > r.tmin ? tnear : tfar
    P = ray(t)
    if P.z > Par.h1 || P.z < Par.h0
        return false
    else
        return true
    end
end

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
## IO for element

# this is to avoid copy and pasting in print below, don't use this
# anywhere else. gnarly hack, the variable names here have to
# correspond to those in print
macro writeandcheck(line)
    esc(quote
        if pretty
            write(io, repeat("  ",depth)*$line*"\n")
        else
            write(io, $line)
        end
        written += 1
        if written == maxlines return end
    end)
end

function Base.print(io::IO, elem::HTMLElement{T},
                    maxlines=Inf, depth=0, written=0; pretty=false) where {T}
    opentag = isempty(elem.attributes) ? "<$T" : "<$T "
    for (name,value) in elem.attributes
        opentag *= "$name=\"$value\""
    end
    opentag *= ">"
    closetag = "</$T>"
    # TODO make inline elements printed all on one line
    if isempty(children(elem))
        @writeandcheck(opentag * closetag)
    else
        @writeandcheck(opentag)
        for child in elem.children
            print(io,child,maxlines,depth+1,written,pretty=pretty)
        end
        @writeandcheck(closetag)
    end
end

prettyprint(io::IO, elem::HTMLElement) = print(io, elem, Inf, pretty=true)
prettyprint(elem::HTMLElement) = print(stdout, elem, pretty=true)

# TODO maybe query tty_cols for a default?
function Base.show(io::IO, elem::HTMLElement)
    write(io,summary(elem)*":\n")
    if get(io, :compact, false)
        write(io, summary(elem))
    elseif get(io, :limit, false)
        print(io, elem, 20, pretty=true)
    else
        print(io, elem, Inf, pretty=true)
    end
end

### IO for Text

function Base.show(io::IO, t::HTMLText)
    write(io,"HTML Text: $(t.text)")
end

function Base.print(io::IO, node::HTMLText,
                    maxlines=Inf, depth=0, written=0; pretty=false)
    if pretty
        for line in split(strip(text(node)), "\n")
            @writeandcheck(line)
        end
    else
        @writeandcheck(text(node))
    end
end

### io for Document

function Base.show(io::IO, doc::HTMLDocument)
    write(io, "HTML Document:\n")
    write(io, "<!DOCTYPE $(doc.doctype)>\n")
    Base.print(io, doc.root, pretty=true)
end

function Base.print(io::IO, doc::HTMLDocument; pretty=false)
    write(io, "<!DOCTYPE $(doc.doctype)>")
    Base.print(io, doc.root, pretty=pretty)
end

prettyprint(io::IO, doc::HTMLDocument) = Base.print(io, doc, pretty=true)
prettyprint(doc::HTMLDocument) = prettyprint(stdout, doc)


text:
using MDCDL
using LinearAlgebra

using TestImages: testimage
using Plots: plot, png

using Base.Filesystem
using Dates: format, now

using ImageFiltering.Algorithm: FIR, FFT
using ComputationalResources

function main()
    ########## Configurations #########
    # choose NSOLT type: (Rnsolt | Cnsolt)
    Nsolt = Rnsolt
    # TP := eltype(Nsolt) (<:AbstractFloat)
    TP = Float64
    # decimation factor: (<:NTuple{D,Int} where D is #dims)
    df = (2,2)
    # polyphase order: (<:NTuple{D,Int} where D)
    ord = (4,4)
    # number of channels: (<:Union{Integer,Tuple{Int,Int}} for Rnsolt)
    #                     (<:Integer for Cnsolt)
    nch = 8

    # size of minibatches (<:NTuple{D,Int})
    szx = (16,16)
    # number of minibatches (<:Integer)
    nSubData = 4

    # path of log files (do nothing if isa(logdir, Nothing))
    # logdir = begin
    #     tm = format(now(), "yyyymmdd_HH_MM_SS_sss")
    #     joinpath(@__DIR__, "results", tm)
    # end
    logdir = nothing
    # save minibatches?
    do_save_trainingset = false

    addresource(CPUThreads)
    @show haveresource(CPUThreads)
    # options for sparse coding
    sc_options = ( iterations = 1000, sparsity = 0.5, filter_domain=:convolution, resource=CPUThreads(FIR()),)
    # options for dictionary update
    du_options = ( iterations = 100, stepsize = 1e-3,)

    # general options of dictionary learning
    options = ( epochs  = 100,
                verbose = :none, # :none, :standard, :specified, :loquacious
                sc_options = sc_options,
                du_options = du_options,
                logdir = logdir,)
    ####################################
    logdir != nothing && !isdir(logdir) && mkpath(logdir)
    do_save_trainingset = do_save_trainingset && logdir != nothing

    # original image
    orgImg = TP.(testimage("cameraman"))

    # generate minibatches of orgImage
    trainingSet = map(1:nSubData) do nsd
        pos = rand.(UnitRange.(0,size(orgImg) .- szx))
        orgImg[UnitRange.(1 .+ pos, szx .+ pos)...]
    end

    if do_save_trainingset
        datadir = joinpath(logdir, "data").
        !isdir(datadir) && mkpath(datadir)
        map(idx->save(joinpath(datadir, "$idx.png"), trainingSet[idx]), 1:nSubData)
    end

    # create NSOLT instance
    nsolt = Nsolt(TP, df, ord, nch)
    # set random orthonormal matrices to the initial matrices.
    MDCDL.rand!(nsolt, isPropMat = false, isPropAng = false, isSymmetry = false)

    # dictionary learning
    MDCDL.train!(nsolt, trainingSet; options...)

    #plot(nsolt)
    return nsolt
end

main()

max_stars_repo_name: nixir/MDCDL

max_issues_repo_name: nixir/MDCDL

max_forks_repo_name: nixir/MDCDL


text:
#MAIN FILE
#---------------------------------------------------------------------

#Cores Files and Variables that are necessary for the software to run.
CoreFiles = ["Directories.txt","TOOLS.jl","CAL.jl","CMP.jl"]
#Packages needed.
packages_list = []

#List of dir that allow the soft to know where are the usefull things
#Format : "Thing" => "Where it is"
Dir = Dict{String,String}()

#---------------------------------------------------------------------
"""
	\tCheckFiles()\n
	\t\tReturn Boolean
Internal Function that check Core Files.\n
"""
function CheckFiles()
	#Check if the Cores Files are inside the current directory
	OK_array = map(x->(isfile(x)||isdir(x)),CoreFiles)
	#Total True == Number of Files ?
	if(sum(OK_array) == length(CoreFiles))
		println("Yes. ",sum(OK_array),"/",length(CoreFiles)," core files in the directory")
	else
		println("Warning, missing files !")
		for i in 1:length(OK_array)
			if(OK_array[i] == false)
				println("\t-> "*CoreFiles[i])
			end
		end
		return false
	end
	return true
end
#---------------------------------------------------------------------
function LoadDir()
	#Take the paths inside Directories.txt and put them into Dir
	open(CoreFiles[1],"r") do d_file
		while ((!eof(d_file))&&isreadable(d_file))
			(k,v) = split(readline(d_file))
			if(isdir(v))
				Dir[k] = v
				println(k*" @ "*v)
			else
				println("\t->"*v*" is not a valid Directory")
				return false
			end
		end
	end
	return true
end
#---------------------------------------------------------------------

# INIT 							~~~~~~~~~~~~~~~~~~
println("\n---  P R O G R A M    S T A R T S  ---")
println("\nEveryones here ? ")
while(secure(CheckFiles()))end

include("TOOLS.jl")
include("CAL.jl")

println("\nLoading Packages...")
while(secure(CheckPkg(packages_list)))end

println("\nLoading DirectoryFile...")
while(secure(LoadDir()))end

println("\nLoading CalibrationFile...")
while(secure(LoadCal("Calibration.txt")))end

# COMPONENTS DATABASE BUILDING 	~~~~~~~~~~~~~~~~~~

# HARDWARE CONNECTIONS 			~~~~~~~~~~~~~~~~~~

# SESSION BUILD 				~~~~~~~~~~~~~~~~~~

# CALIBRATE 					~~~~~~~~~~~~~~~~~~

# EXTRACTION CHOICES 			~~~~~~~~~~~~~~~~~~

# IMAGE GATHERING 				~~~~~~~~~~~~~~~~~~

# PRE-PROCESSING				~~~~~~~~~~~~~~~~~~

# PROCESSING (each comp, H&V)	~~~~~~~~~~~~~~~~~~

# POSITIONS REVIEW				~~~~~~~~~~~~~~~~~~

# SEND HARDWARE TASK			~~~~~~~~~~~~~~~~~~

# GATHER HARDWARE LOG			~~~~~~~~~~~~~~~~~~

# CLOSING 						~~~~~~~~~~~~~~~~~~
println("\n---  P R O G R A M    E N D I N G  ---")


text:
using VoronoiCells
using Test

include("Points.jl")
include("Rectangle.jl")
include("Clipping.jl")
include("Sort.jl")
include("Cells.jl")
include("Area.jl")
include("Plot.jl")


text:
using Test
using PETScX

@testset "mat tests" begin
    for petsclib in PETScX.petsclibs
        PETScX.Initialize(petsclib)
        PetscScalar = PETScX.scalartype(petsclib)
        PetscInt = PETScX.inttype(petsclib)

        # Create a matrix
        m = 10
        n = m
        nz = 3
        A = PETScX.MatAIJ(petsclib, m, n, nz)

        # Check some dimension information
        row_rng = PETScX.ownershiprange(A)
        @test row_rng == 0:(m - 1)
        @test size(A) == (m, n)
        @test length(A) == m * n

        vals = PetscScalar.([1, -1, -1, 1])
        for r in row_rng
            rows = cols = PetscInt.([r - 1, r])
            PETScX.setvalues!(A, rows, cols, vals, PETScX.ADD_VALUES)
        end

        # Fix up the last entry
        rows = cols = PetscInt.([m - 1])
        PETScX.setvalues!(A, rows, cols, vals, PETScX.ADD_VALUES)

        # Do assembly
        PETScX.assemblybegin!(A)
        PETScX.assemblyend!(A)

        # Check the the values we inserted are correct
        vals = zeros(PetscScalar, 3)
        for r in row_rng
            rows = PetscInt.([r])
            if r == 0
                cols = PetscInt.([r, r + 1])
                PETScX.getvalues!(vals, A, rows, cols)
                @test PetscScalar.([2, -1]) == vals[1:2]
            elseif r == m - 1
                cols = PetscInt.([r - 1, r])
                PETScX.getvalues!(vals, A, rows, cols)
                @test PetscScalar.([-1, 2]) == vals[1:2]
            else
                cols = PetscInt.([r - 1, r, r + 1])
                PETScX.getvalues!(vals, A, rows, cols)
                @test PetscScalar.([-1, 2, -1]) == vals
            end
        end

        PETScX.destroy(A)

        PETScX.Finalize(petsclib)
    end
end


text:
module PkgPage

import Franklin
import CommonMark
import NodeJS
using Pkg: TOML
using FranklinUtils

export newpage
export serve

export lx_lead, lx_figure, lx_alert
export env_column, env_table, env_section

export hfun_col

const F = Franklin

include("instantiate.jl")

include("latex/basic.jl")
include("latex/environments.jl")

end # module


text:
export read_quake

@doc """
    Ev = read_quake(fmt, file [, keywords])

Read data in file format `fmt` from `file` into SeisEvent object `Ev`.
* Formats: suds, qml, uw
* Keywords: full, v

Note: because earthquake data are usually discrete, self-contained files, no "in-place" version of `read_quake` exists, and  `read_quake` doesn't support wildcards in the file string.

See also: `read_data`, `get_data`, `read_meta`, `UW.readuwevt`
""" read_quake
function read_quake(fmt::String, fname::String;
  full    ::Bool    = false,              # full header
  v       ::Integer = KW.v                # verbosity level
  )

  if fmt == "suds"
    Ev = SUDS.readsudsevt(fname, full=full, v=v)
  elseif fmt == "uw"
    Ev = UW.readuwevt(fname, full=full, v=v)
  elseif fmt in ("qml", "quakeml")
    hdr, source = read_qml(fname)
    Ev = SeisEvent(hdr = hdr[1], source = source[1])
  end
  return Ev
end


text:
using SimpleML
using Base.Test

using HDF5
using JLD
using Calculus

f=jldopen(joinpath(Pkg.dir(),"SimpleML/test/mnist.jld"), "r")

function testLogisticOnevsAll()

    X=read(f, "X")
    y=read(f, "y")
    m = size(X, 1); #number of training set data points
    n = size(X, 2); #number of features
    X = [ones(m, 1) X]

    all_theta = learnLogisticAll(X, y, .1, 10)

    p=predictLogisticAll(all_theta, X)
    acc=mean(int(p.==y))
    println("Accuracy of training set for Logistic: $acc")
    @assert acc > .95

end

function testNN2()
    input_sz  = 400
    hidden_sz = 25
    num_labels = 10

    #Read input data and saved coefficients
    X=read(f, "X")
    y=read(f, "y")
    m = size(X, 1); #number of training set data points
    n = size(X, 2); #number of features
    fwts=jldopen(joinpath(Pkg.dir(),"SimpleML/test/wts.jld"), "r")
    theta1 = read(fwts, "theta1")
    theta2 = read(fwts, "theta2")

    #Test prediction function
    p=SimpleML.predict2nn(theta1, theta2, X)
    acc=mean(int(p.==y))
    @test_approx_eq acc 0.9758

    #Test Cost Function
    params = [theta1[:] ; theta2[:] ]
    theta1_grad = zeros(size(theta1))
    theta2_grad = zeros(size(theta2))
    grad = [theta1_grad[:] ; theta2_grad[:]]

    lambda=0
    @test_approx_eq_eps SimpleML.nn2Cost(params, X, y, input_sz,  hidden_sz, num_labels,  lambda, grad) 0.2382217 1e-6

    lambda=1
    @test_approx_eq_eps SimpleML.nn2Cost(params, X, y, input_sz,  hidden_sz, num_labels, lambda, grad) 0.4602118 1e-6

    #Test learning
    theta1, theta2 = learn2nn(X, y, input_sz, hidden_sz, num_labels, lambda)
    p=predict2nn(theta1, theta2, X)
    acc=mean(int(p.==y))
    println("Accuracy of training set for NN: $acc")
    @assert acc > .97

end

#test gradient calculations for dummy data
function testNNGradient(lambda)
    lambda=0
    input_sz = 3
    hidden_sz = 5
    num_labels = 3
    m = 5
    theta1 = debugInitialWeights(hidden_sz, input_sz)
    theta2 = debugInitialWeights(num_labels, hidden_sz)

    X  = debugInitialWeights(m, input_sz - 1)
    y  = 1 + mod([1:m], num_labels)'
    params = [theta1[:] ; theta2[:]]
    grad = zeros(size(params))

    J=SimpleML.nn2Cost(params, X, y, input_sz, hidden_sz, num_labels, lambda, grad)
    costFunc = (p) -> SimpleML.nn2Cost(p, X, y, input_sz, hidden_sz, num_labels, lambda, nothing)
    numgrad = derivative(costFunc, params, :central)
    diff = norm(numgrad-grad)/norm(numgrad+grad)
    @assert diff < 1e-9
end

#generate a well distributed, but deterministic values for initial weights, for testing
function debugInitialWeights(l_out, l_in)
    W = zeros(l_out, 1 + l_in)
    W = reshape(sin(1:length(W)), size(W)) / 10
    return W
end

#Run the tests
testLogisticOnevsAll()
@test_approx_eq_eps SimpleML.sigmoidGradient([1 -0.5 0 0.5 1])  [0.196612 0.235004 0.250000 0.235004 0.196612] 1e-6
testNNGradient(0)
testNNGradient(3)
testNN2()


text:
module Bridges

using LinearAlgebra

using MathOptInterface
using ConstraintProgrammingExtensions
const CP = ConstraintProgrammingExtensions
const MOI = MathOptInterface
const MOIU = MOI.Utilities
const MOIB = MOI.Bridges
const MOIBC = MOIB.Constraint

include("helpers.jl")
include("Constraint/Constraint.jl")

end


text:
using Base.Test
using Autoreload

m1_before = """
module M
type MyType
  x::Int
end

f(var::MyType) = 1
end
"""

m1_after = """
module M
type MyType
  x::Int
end

f(var::MyType) = 2
end
"""

open(joinpath(dirname(@__FILE__), "M.jl"), "w") do f
    write(f, m1_before)
end

using Autoreload
arequire("M")
my_var = M.MyType(5)
M.f(my_var)

open(joinpath(dirname(@__FILE__), "M.jl"), "w") do f
    write(f, m1_after)
end

areload()
@test M.f(my_var) == 2

rm(joinpath(dirname(@__FILE__), "M.jl"))


text:
module Multicombinations
using Base

import Base: start, next, done, eltype, length

export multichoose, multicombinations, integersums, is2mc

# multichoose

function multichoose(n, k)
    binomial(n + k - 1, k)
end

# multicombinations

immutable MState{T}
    a::T
    k::Int
end

eltype(c::MState) = typeof(c.a)

eltype{T}(c::MState{UnitRange{T}}) = Array{T, 1}

eltype{T}(c::MState{Range{T}}) = Array{T, 1}

length(c::MState) = multichoose(length(c.a), c.k)

multicombinations(a, k::Integer) = MState(a, k)

function nextmulticombination(c::MState, s)
    s = copy(s)
    for i = c.k:-1:1
        if s[i] < endof(c.a)
            s[i] = nextind(c.a, s[i])
            for j = i:c.k-1
                s[j+1] = s[j]
            end
            return s
        end
    end
    return [0]
end

done(c::MState, s) = length(s) > 0 && s[1] < 1

function start(c::MState)
    if c.k > 0 && length(c.a) > 0
        ones(Int, c.k)
    elseif c.k == 0
        Array(Int, 0)
    else
        [0]
    end
end

function next(c::MState, s)
    r = c.a[s]
    s = nextmulticombination(c, s)
    (r, s)
end

# integer solutions to x_1 + x_2 + ... + x_n = k

immutable IState
    n::Int
    k::Int
end

eltype(c::IState) = Array{Int, 1}

length(c::IState) = multichoose(c.n, c.k)

integersums(n::Integer, k::Integer) = IState(n, k)

done(c::IState, s) = length(s) > 0 && s[1] < 0

function start(c::IState)
    if c.k >= 0 && c.n > 0
        [c.k, zeros(Int, c.n - 1)]
    elseif c.k == 0 && c.n == 0
        Array(Int, 0)
    else
        [-1]
    end
end

function nextintegersum(c::IState, s)
    if length(s) == 0 || s[end] == c.k
        return [-1]
    end
    s = copy(s)
    sum = 0
    i = 1
    while i < length(s)
        sum += s[i]
        if sum == c.k
            s[i] -= 1
            s[i+1] += 1
            return s
        end
        i += 1
    end
    i -= 1
    while s[i] == 0
        i -= 1
    end
    s[i] -= 1
    last = s[end]
    s[end] = 0
    s[i+1] = last + 1
    return s
end

next(c::IState, s) = (s, nextintegersum(c, s))

# conversion

function is2mc(s)
    ret = Array(Int, sum(s))
    k = 1
    for i = 1:length(s)
        for j = 1:s[i]
            ret[k] = i
            k += 1
        end
    end
    ret
end

end


text:
module net

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
function FixedNumberInterpolationGrid(arg0::jint)
    return FixedNumberInterpolationGrid((jint,), arg0)
end

function get_grid_points(obj::FixedNumberInterpolationGrid, arg0::jdouble, arg1::jdouble)
    return jcall(obj, "getGridPoints", Vector{jdouble}, (jdouble, jdouble), arg0, arg1)
end


text:
module TestTaxon
    using Test
    using NCBITaxonomy
    using StringDistances

    # Match with a known ID
    @test taxon(1234567) == NCBITaxon("Exiligada punctata", 1234567)

    # Match with a case insensitive search
    @test_throws NameHasNoDirectMatch taxon("exiligada punctata")
    @test taxon("exiligada punctata"; casesensitive=false) == NCBITaxon("Exiligada punctata", 1234567)

    # Strict matching with a scientific name
    bos = taxon("Bos taurus")
    @test typeof(bos) == NCBITaxon
    @test bos.name == "Bos taurus"
    @test bos.id == 9913

    # Strict matching with a vernacular name
    cow = taxon("cow")
    @test typeof(cow) == NCBITaxon
    @test cow.name == "Bos taurus"
    @test cow.id == 9913

    # Test show method
    @test sprint(io -> show(io, cow)) == "Bos taurus (ncbi:9913)"

    # Fuzzy matching with a scientific name
    box = taxon("Box taurus"; strict=false)
    @test typeof(box) == NCBITaxon
    @test box.name == "Bos taurus"
    @test box.id == 9913

    # Fuzzy matching with a scientific name and a custom distance
    box = taxon("Box taurus"; strict=false, dist=StringDistances.DamerauLevenshtein)
    @test typeof(box) == NCBITaxon
    @test box.name == "Bos taurus"
    @test box.id == 9913

    #Vernacular name
    chub = vernacular(ncbi"Leuciscus cephalus")
    @test "European chub" in chub
    @test "chub" in chub

    # Vernacular missing
    @test isnothing(vernacular(ncbi"Lamellodiscus elegans"))

    # Synonyms
    @test "Bos bovis" in synonyms(ncbi"Bos taurus")
    @test isnothing(synonyms(ncbi"Lamellodiscus elegans"))

    # Limit by rank
    @test_throws AssertionError taxon("Lamellodiscus"; rank=:specs)
    @test_throws NameHasNoDirectMatch taxon("Lamellodiscus"; rank=:species)
    @test typeof(taxon("Lamellodiscus"; rank=:genus)) <: NCBITaxon

    # Check with multiple names
    @test taxon(10088) in alternativetaxa("Mus")
    @test taxon(862507) in alternativetaxa("Mus")

    # Check the pairs
    @test Pair(taxon(10088), 1.0) in similarnames("mouse")
    @test Pair(taxon(10090), 1.0) in similarnames("mouse")

end


text:
using SQLite
using DataFrames
using LibPQ
using BenchmarkTools

function main()
    sqlite_db = SQLite.DB("fakedata.db")
    pg_conn = LibPQ.Connection("dbname=postgres user=postgres password=postgres host=localhost")
    user_df = DBInterface.execute(sqlite_db, "SELECT * FROM users") |> DataFrame

    execute(pg_conn, "DROP TABLE IF EXISTS users")
    execute(
        pg_conn,
        """
        CREATE TABLE users (
            user_id INTEGER PRIMARY KEY,
            user_name TEXT,
            user_age INTEGER,
            user_address TEXT
        )
        """
    )

    execute(pg_conn, "BEGIN;")

    LibPQ.load!(
        user_df,
        pg_conn,
        """
        INSERT INTO users (user_id, user_name, user_age, user_address)
        VALUES            (    \$1,       \$2,      \$3,          \$4)
        """
    )

    execute(pg_conn, "COMMIT;")

    close(pg_conn)
end

@benchmark main()

max_stars_repo_name: renzmann/postgres-loading

max_issues_repo_name: renzmann/postgres-loading

max_forks_repo_name: renzmann/postgres-loading


text:
using CSV
using DataFrames
using LinearAlgebra
using Plots

function ansatz(r1,r2,α)
    norm_r1 = norm(r1)
    norm_r2 = norm(r2)
    r12 = norm(r1-r2)
    wf = exp(-2*norm_r1)*exp(-2*norm_r2)*exp(r12/(2*(1+α*r12)))
    return wf
end

function prob_density(r1,r2,α)
    return ansatz(r1,r2,α)^2
end

function E_local(r1,r2,α)
    norm_r1 = norm(r1)
    norm_r2 = norm(r2)
    r12 = norm(r1-r2)
    dot_product = dot(r1/norm_r1-r2/norm_r2,r1-r2)
    energy = -4+dot_product/(r12*(1+α*r12)^2)-1/(r12*(1+α*r12)^3)-1/(4*(1+α*r12)^4)+1/r12
    return energy
end

function metropolis(N, α)

    L = 1
    r1 = rand(3)*2*L .- L
    r2 = rand(3)*2*L .- L #random number from -L to L
    E = 0
    E2 = 0
    Eln_average = 0
    ln_average = 0
    rejection_ratio = 0
    step = 0
    max_steps = -1

    for i in 1:N
        chose = rand()
        step = step + 1
        if chose < 0.5
            r1_trial = r1 + 0.5*(rand(3)*2*L .- L)
            r2_trial = r2
        else
            r2_trial = r2 + 0.5*(rand(3)*2*L .- L)
            r1_trial = r1
        end
        if prob_density(r1_trial,r2_trial,α) >= prob_density(r1,r2,α)
            r1 = r1_trial
            r2 = r2_trial
        else
            dummy = rand()
            if dummy < prob_density(r1_trial,r2_trial,α)/prob_density(r1,r2,α)
                r1 = r1_trial
                r2 = r2_trial
            else
                rejection_ratio += 1/N
            end
        end

        if step > max_steps
            E += E_local(r1,r2,α)/(N-max_steps)
            E2 += E_local(r1,r2,α)^2/(N-max_steps)
            r12 = norm(r1-r2)
            Eln_average += (E_local(r1,r2,α)*-r12^2/(2*(1+α*r12)^2))/(N-max_steps)
            ln_average += -r12^2/(2*(1+α*r12)^2)/(N-max_steps)
        end
    end
    return E, E2, Eln_average, ln_average, rejection_ratio
end

function iterate_alpha()
	α = 0

	alpha_iterations = 15
	N_metropolis = 5000
	random_walkers = 200
	γ = 0.5

	energy_plot = []
	alpha_plot = []
	variance_plot = []

	for i in 1:alpha_iterations
	    E = 0
	    E2 = 0
	    dE_dalpha = 0
	    Eln = 0
	    ln = 0
	    rejection_ratio = 0

	    for j in 1:random_walkers
	        E_met, E2_met, Eln_met, ln_met, rejections_met = metropolis(N_metropolis,α)
	        E += E_met/random_walkers
	        E2 += E2_met/random_walkers
	        Eln += Eln_met/random_walkers
	        ln += ln_met/random_walkers
	        rejection_ratio += rejections_met/random_walkers
	    end

		α = α + 0.05

	    energy_plot = append!(energy_plot, E)
	    alpha_plot = append!(alpha_plot, α)
	    variance_plot = append!(variance_plot, E2-E^2)
	end
	return energy_plot,alpha_plot,variance_plot
end

energy,alpha,variance = iterate_alpha();
CSV.write(datadir("sims/helium_atom.csv"),(Energy = energy, Variance = variance, Alpha = alpha))

energy_plot = plot()
plot!(energy_plot,energy,
	label = false,
	title = "Energy expenctation Plot",
	linecolor = :blue,
	ylabel = "<E>",
	xlabel = "timestep"
)
safesave(plotsdir("energy_helium_atom.png"),energy_plot)

variance_plot = plot()
plot!(variance_plot,variance,
	title = "Variance in energy",
	label = false,
	linecolor = :blue,
	ylabel = "Var(E)",
	xlabel = "timestep"
)
safesave(plotsdir("variance_helium_atom.png"),variance_plot)

alpha_plot = plot()
plot!(alpha_plot,alpha,
	title = "Variation Parameter",
	label = false,
	linecolor = :blue,
	ylabel = "α",
	xlabel = "timestep"
)
safesave(plotsdir("alpha_helium_atom.png"),alpha_plot)

alpha_vs_energy_plot = plot()
scatter!(alpha_vs_energy_plot,alpha,
	energy,
	title = "Energy vs α",
	xlabel = "α",
	ylabel = "<E>",
	label = false,
	markercolor = :blue
 )
safesave(plotsdir("alpha_vs_energy_helium_atom.png"),alpha_vs_energy_plot)

max_stars_repo_name: Ved-Mahajan/Variational-Monto-Carlo

max_issues_repo_name: Ved-Mahajan/Variational-Monto-Carlo

max_forks_repo_name: Ved-Mahajan/Variational-Monto-Carlo


text:
# Single-line boxes show the number of bytes, double-lined have
# a variable number of bytes. E.g. here: 1, 1, 4, N, 4.
# +---+---+===============+---+---+---+---+
# |CMF|FLG|COMPRESSED DATA|     ADLER32   |
# +---+---+===============+---+---+---+---+

"""
    zlib_decompress!(
        ::Decompressor, output, input, [n_out::Integer]
    )::Union{LibDeflateError, Int}

Zlib decompress from `input` to `output`.
If the precise number of decompressed bytes is known, pass it in as `n_out`
for added performance. If it is wrong, the function will return an error.

Return the number of bytes written, or a `LibDeflateError`.

See also: [`unsafe_zlib_decompress!`](@ref)
"""
function zlib_decompress! end

function zlib_decompress!(
    decompressor::Decompressor, output, input
)::Union{LibDeflateError,Int}
    GC.@preserve output input begin
        write = WriteableMemory(output)
        read = ReadableMemory(input)
        unsafe_zlib_decompress!(
            Base.SizeUnknown(),
            decompressor,
            pointer(write),
            sizeof(write),
            pointer(read),
            sizeof(read),
        )
    end
end

function zlib_decompress!(
    decompressor::Decompressor, output, input, n_out::Integer
)::Union{LibDeflateError,Int}
    GC.@preserve output input begin
        write = WriteableMemory(output)
        read = ReadableMemory(input)
        n_out > sizeof(write) && return LibDeflateErrors.deflate_insufficient_space
        unsafe_zlib_decompress!(
            Base.HasLength(),
            decompressor,
            pointer(write),
            n_out,
            pointer(read),
            sizeof(read),
        )
    end
end

"""
    unsafe_zlib_decompress!(
        size::Union{Base.SizeUnknown, Base.HasLength},
        ::Decompressor,
        out_ptr::Ptr, n_out::Integer,
        in_ptr::Ptr, len::Integer
    )::Union{LibDeflateError, Int}

Zlib decompress data beginning at `in_ptr` and `len` bytes onwards, into `out_ptr`.
Return the number of bytes written, or a `LibDeflateError`.
`size`` can be `SizeUnknown` or `HasLength`. If the former, `n_out` tells how much space
is available at the output. If the latter, `n_out` is the exact number of bytes
that the payload decompresses to. The latter is faster.

See also: [`zlib_decompress!`](@ref)
"""
function unsafe_zlib_decompress!(
    size::Union{Base.SizeUnknown,Base.HasLength},
    decompressor::Decompressor,
    out_ptr::Ptr,
    n_out::Integer,
    in_ptr::Ptr,
    len::Integer,
)::Union{LibDeflateError,Int}
    # Must be at least 6 bytes in length + compressed
    len < 6 && return LibDeflateErrors.zlib_input_too_short

    ptr = in_ptr
    header = ltoh(unsafe_load(Ptr{UInt16}(ptr)))
    ptr += 2

    # Parse CMF: First 4 bits must be 0x8 = DEFLATE algorithm
    header & 0x000f != 0x0008 && return LibDeflateErrors.zlib_not_deflate

    # Next 4 bits must be 7, as the window size in libdeflate is hardcoded to 32 KiB.
    header & 0x00f0 != 0x0070 && return LibDeflateErrors.zlib_wrong_window_size

    # libdeflate does not support a custom decompression dict, I think
    header & 0x2000 != 0x0000 && return LibDeflateErrors.zlib_needs_compression_dict

    # This is ntoh, because the header checksum is interpreted as a big-endian integer.
    iszero(mod(ntoh(header), UInt16(31))) || return LibDeflateErrors.zlib_bad_header_check

    # Decompress payload
    nbytes = unsafe_decompress!(size, decompressor, out_ptr, n_out, ptr, len - 6)
    nbytes isa LibDeflateError && return nbytes

    # Then check adler32, also stored as big-endian
    exp_adler32 = ntoh(unsafe_load(Ptr{UInt32}(in_ptr) + len - 4))
    obs_adler32 = unsafe_adler32(out_ptr, nbytes)
    obs_adler32 == exp_adler32 || return LibDeflateErrors.zlib_bad_adler32

    return nbytes
end

"""
    zlib_compress!(
        ::Compressor, output, input
    )::Union{LibDeflateError, Int}

Zlib compress from `input` to `output`.
Return the number of bytes written, or a `LibDeflateError`.

See also: [`unsafe_zlib_compress!`](@ref)
"""
function zlib_compress!(compressor::Compressor, output, input)::Union{LibDeflateError,Int}
    GC.@preserve output input begin
        write = WriteableMemory(output)
        read = ReadableMemory(input)
        unsafe_zlib_compress!(
            compressor, pointer(write), sizeof(write), pointer(read), sizeof(read)
        )
    end
end

"""
    unsafe_zlib_compress!(
        ::Compressor,
        out_ptr::Ptr, max_outlen::Integer,
        in_ptr::Ptr, len::Integer
    )::Union{LibDeflateError, Int}

Zlib compress data beginning at `in_ptr` and `len` bytes onwards, into `out_ptr`.
Return the number of bytes written, or a `LibDeflateError`.

See also: [`zlib_compress!`](@ref)
"""
function unsafe_zlib_compress!(
    compressor::Compressor, out_ptr::Ptr, max_outlen::Integer, in_ptr::Ptr, len::Integer
)::Union{LibDeflateError,Int}
    max_outlen < 6 && return LibDeflateErrors.zlib_insufficient_space
    # Be sure to check that all these 4 results are zero mod 31 when big-endian ordered
    header = if compressor.level == 1
        0x0178
    elseif compressor.level == DEFAULT_COMPRESSION_LEVEL
        0x5e78
    elseif compressor.level == 12
        0xda78
    else
        0x9c78
    end
    unsafe_store!(Ptr{UInt16}(out_ptr), htol(header))
    out_ptr += 2
    n_bytes = unsafe_compress!(compressor, out_ptr, max_outlen - 6, in_ptr, len)
    n_bytes isa LibDeflateError && return n_bytes
    out_ptr += n_bytes
    # Alder32 is also big-endian stored
    unsafe_store!(Ptr{UInt32}(out_ptr), hton(unsafe_adler32(in_ptr, len)))
    return n_bytes + 6
end


text:
# Autogenerated wrapper script for yq_jll for i686-linux-gnu
export yq

JLLWrappers.@generate_wrapper_header("yq")
JLLWrappers.@declare_executable_product(yq)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        yq,
        "bin/yq",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()