text:
# The finite volume operators that are included automatically if FV is used.

function sym_central_op(f)
    # Input will be an array of Basic. Output should be in a similar array.
    if typeof(f) <: Array
        result = copy(f);
        for i=1:length(result)
            result[i] = sym_central_op(f[i]);
        end

    elseif typeof(f) <: Number
        # If the input was just a constant, the result will just be that constant.
        result = Basic(f);

    elseif typeof(f) == Basic
        # Apply a CELLn tag to signal to the code generator which side of the face the value is from.
        side1 = apply_flag_to_all_symbols("CELL1", f);
        side2 = apply_flag_to_all_symbols("CELL2", f);
        result = Basic(0.5) .* (side1 .+ side2);
    end
    return result;
end

function sym_upwind_op(v, f, alpha=Basic(0))
    # Input will be an array of Basic. Output should be in a similar array.
    if typeof(f) <: Array
        result = copy(f);
        for i=1:length(result)
            result[i] = sym_upwind_op(v, f[i], alpha);
        end

    elseif typeof(f) == Basic
        # Apply a CELLn tag to signal to the code generator which side of the face the value is from.
        side1 = apply_flag_to_all_symbols("CELL1", f);
        side2 = apply_flag_to_all_symbols("CELL2", f);
        # The formula for an adjustable first order upwind scheme.
        # F(u) = 0.5*(side1+side2)*(v.normal) + 0.5*(side1-side2)*abs(v.normal)*(1-alpha)
        result = Basic(0.5) .* ( (sym_dot_op(v, sym_normal_op())) .* (side1 + side2) .+ abs.(sym_dot_op(v, sym_normal_op())[1]) .* (Basic(1.0)-alpha) .* (side1 - side2) );
        result = result[1]; # Since dot and normal returned arrays, result was placed in an array, but this is done above redundantly.

    elseif typeof(f) <: Number
        # If the input was just a constant, the result will just be that constant.
        result = Basic(f);
    end
    return result;
end

function sym_burgerGodunov_op(u, f)
    if typeof(f) <: Array
        result = copy(f);
        for i=1:length(result)
            result[i] = sym_burgerGodunov_op(u[i], f[i]);
        end
    elseif typeof(f) == Basic
        uside1 = apply_flag_to_all_symbols("CELL1", u);
        uside2 = apply_flag_to_all_symbols("CELL2", u);
        fside1 = apply_flag_to_all_symbols("CELL1", f);
        fside2 = apply_flag_to_all_symbols("CELL2", f);
        result = conditional(isgreaterthan(sym_normal_op()[1], Basic(0)), # if normal[1] > 0
                    conditional(isgreaterthan(uside1, uside2), # if u1 < u2
                        symbolmax(fside1, fside2), # max(f1,f2)
                        # else
                        conditional(isgreaterthan(uside1*uside2, Basic(0)), # if u1*u2 > 0 (same sign)
                            symbolmin(fside1, fside2), # min(f1,f2)
                            # else
                            Basic(0))), # 0
                    # elseif normal < 0 (just reverse 1 and 2)
                    conditional(isgreaterthan(uside2, uside1), # if u2 < u1
                        symbolmax(fside2, fside1), # max(f1,f2)
                        # else
                        conditional(isgreaterthan(uside2*uside1, Basic(0)), # if u1*u2 > 0 (same sign)
                            symbolmin(fside2, fside1), # min(f1,f2)
                            # else
                            Basic(0)))) # 0
    elseif typeof(f) <: Number
        result = Basic(f);
    end
    return result;
end

_names = [:upwind, :burgerGodunov, :central];
_handles = [sym_upwind_op, sym_burgerGodunov_op, sym_central_op];


text:
### A Pluto.jl notebook ###
# v0.14.5

using Markdown
using InteractiveUtils

# ╔═╡ d4260fcc-784e-4341-ae84-c7e8ab288c4b
using LaTeXStrings, DifferentialEquations, LinearAlgebra, SparseArrays

# ╔═╡ 82272c06-3861-4aed-a08d-bed2e84570fc
begin
	using Plots
	gr()
end

# ╔═╡ b91fd324-adf3-11eb-36e5-59bf4de5e0f9
md"# 2D Multiple Attractor System"

# ╔═╡ 41f97746-5818-410f-b3cb-e59aa9a5bc99
md"_Setup vector field equations. Can't split into jump process._"

# ╔═╡ 5a0a71b5-0b18-42be-b567-c01feb0d9af6
begin
	α = 0.1
	β = 3
	n = 4
	dfx(x, y) = 9x + 9y - (1 + 2x^3 + 2y^3)
	dfy(x, y) = 1 + 2x^3 + 11y - (11x + 2y^3)
end

# ╔═╡ f56451d1-6f67-494b-82b0-2369558d6819
begin
	n_min = -3
	n_max = 3
end

# ╔═╡ a342ae72-b363-4142-929c-78e682871d81
md"_Compute the vector field_"

# ╔═╡ 7cb95f1f-8a19-484c-a33f-a13348dce332
begin
	dq = (n_max - n_min) / 10.
	xs = n_min-1:dq:n_max+1
	ys = n_min-1:dq:n_max+1

	dxs = n_min-1:.05:n_max+1
	dys = n_min-1:.05:n_max+1

	df(x, y) = normalize([dfx(x, y), dfy(x,y)]) ./ 1.75

	xxs = [x for x in xs for y in ys]
	yys = [y for x in xs for y in ys]

	quiver(xxs, yys, quiver=df)
	contour!(dxs, dys, dfx, levels=[0], color=:orange, label="X Nullcline", colorbar = false)
	p1 = contour!(dxs, dys, dfy, levels=[0], color=:green, label = "Y Nullcline", colorbar = false)
	p1
end

# ╔═╡ 7d1e052e-c9b6-4fa3-b4aa-980e60e53961
md"_Setup the Stochastic Differential Equation model with diagonal noise_"

# ╔═╡ 601e7f8b-efa8-4710-b088-78b825b0454a
function f(du, u, p, t)
	du[1] = dfx(u...)
	du[2] = dfy(u...)
	du[1], du[2]
end

# ╔═╡ e3b32bba-b8c8-48ba-8918-c73f1f1e9db9
begin
	g(du, u, p, t) = begin
	  du[1] = 1.25
	  du[2] = 1.25
	end
	sde_prob= SDEProblem(f, g, [0., 0.], (0.0, 1000.0))
	sde_sol = solve(sde_prob, SRIW1())
end

# ╔═╡ 67d5512b-fe45-47d7-9b18-c847a8010577
pstoch_1 = plot(sde_sol,title="Stochastic Simulation", labels=[L"x" L"y"])

# ╔═╡ a1d79d06-afdd-4097-94ac-716779173e6c
pstoch_2 = plot(sde_sol[1000:end], vars = [(1,2)], label="", xlabel=L"x", ylabel = L"y", title="Time Evolution State Plot")

# ╔═╡ eef51429-8458-421a-8366-35577019e204
# pstoch = plot(pstoch_1, pstoch_2, layout=(1, 2), size=(900, 350))

# ╔═╡ 10c9d45f-d8ab-4acd-88ba-8a80a747e30f
# savefig(pstoch, "ch3_four_attractor_stoch.pdf")

# ╔═╡ 9e084100-e5d1-488c-8eaf-769711af1d95
md"_Helper function to tip towards different quadrants for each run of the simulation_"

# ╔═╡ dc5e6e4a-93d6-47c3-80ef-8d22d3483e99
function get_u0(i)
	offset = 1.0
	new_u0 = [0.0, 0.0]
	if i % 4 == 0
		new_u0 .+= [offset, offset]
	elseif i % 4 == 1
		new_u0 .+= [-1* offset, offset]
	elseif i % 4 == 2
		new_u0 .+= [-1* offset, -1* offset]
	else
		new_u0 .+= [offset,-1* offset]
	end
	new_u0
end

# ╔═╡ 3481f158-573b-4ca4-9f2f-fa0e4b68edb4
md"_Add smaller noise term for longer time simulation._

Note that this is equivalent to modifying `D` in the Fokker-Plank Equations and yeilds different probability distributions."

# ╔═╡ 25b8d3e0-4402-4e84-91a0-c5d1c4684ac4
begin
	g2(du, u, p, t) = begin
	  du[1] = 0.5
	  du[2] = 0.5
	end
end

# ╔═╡ 1870a27b-683b-463a-b5cc-9ac073110091
md"_Helper to compute the 2D histograms for unnomalized probability_"

# ╔═╡ c1b075dc-c4dc-4bd9-bbd9-bbb0119fc54a
begin
	nbins = 257
	bin_min = -10.
	bin_max = 10.
	xbins = range(bin_min, bin_max, length=nbins)
	ybins = range(bin_min, bin_max, length=nbins)
	xbins

	function bucket_idx(val, vec)
		idx = 1
		while val > vec[idx]
			idx += 1
		end
		idx - 1
	end

	function althist(x, y, dt, xedges, yedges)
	   counts = spzeros(length(xedges)-1, length(yedges)-1)

		for i=1:length(x) - 1
			r = bucket_idx(x[i], yedges)
			c = bucket_idx(y[i], xedges)
			counts[r,c] += dt[i]
		end

 		counts
	end

	drop_num = 2000
end

# ╔═╡ 351e030b-0d6c-42ea-8c5b-9cad2af5c12c
unnorm = spzeros(length(xbins)-1, length(ybins)-1)

# ╔═╡ 8b509897-8b02-4af7-a0bd-990f2cb3df87
md"_Rerun the following cell to add more data points to the simulation._"

# ╔═╡ 75bcc971-1a89-4dd4-a5ba-a1397f5db274
begin
	for i in 1:4*100
		u00 = get_u0(i)
		sde_prob_long = SDEProblem(f, g2, u00, (0.0, 1e3))
		sde_sol_long= solve(sde_prob_long, SRIW1())
		unnorm .+= althist(sde_sol_long[1, drop_num:end], sde_sol_long[2, drop_num:end], diff(sde_sol_long.t[drop_num:end]), xbins, ybins)
	end
	unnorm
end

# ╔═╡ 86019cdf-c7de-4495-ba61-de0852f2b496
md"_Normalized and visualize the probability and resulting Quasi-potential._"

# ╔═╡ 13347146-e598-4d25-9d92-4b047a9de492
prob = unnorm ./ sum(unnorm)

# ╔═╡ 6c523287-b8db-4ef7-a08a-f61bf2639278
begin
	plotly()
	p = surface(- log.(prob))
	gr()
	p
end

# ╔═╡ bc84c949-6b26-4e02-84bf-a2093508a0dc
begin
	plotly()
	pp = surface(xbins, ybins,prob,c = :grayC)
	gr()
	pp
end

# ╔═╡ 6af8e6ef-2444-4c88-ad81-a645b8971527
begin
	xplot = range(bin_min, bin_max, length=nbins-1)
	yplot = range(bin_min, bin_max, length=nbins-1)
end

# ╔═╡ b3a3b34b-ffff-4b35-952e-d844b5830ce7
 begin
	pp1 = surface(xplot, yplot,-log.(prob), c = :thermal)
	pp2 =  surface(xplot, yplot, prob, c = :grayC)
	ps = plot(pp1, pp2, layout =(2, 1), size=(500, 800))
end

# ╔═╡ 4e4f16b1-7c8a-4435-91f9-6a0d4ff54b22
# savefig(ps, "ch3_4a_ad_hoc.pdf")

# ╔═╡ c37e269f-ac80-4bd1-81d0-e5a476e85272
md"_Helper function to generate -∇ ² to solve the Helmhotz equation_"

# ╔═╡ 164f3c77-5e95-4feb-8783-54c381f4b184
begin
	import DiffEqOperators
	const DiffEqOp = DiffEqOperators

	eye(N, M) = sparse(I, N, M)
	eye(M) = eye(M, M)
	diff1(M) = [ [1.0 zeros(1, M - 1)]; diagm(1 => ones(M - 1)) - eye(M) ]
	sdiff1(M) = sparse(diff1(M))
	flatshape(X) = reshape(X, length(X))

	# make the discrete -Laplacian in 2d, with Dirichlet boundaries
	function laplacian2d(Nx, Ny, Lx, Ly)
		dx = Lx / (Nx + 1)
		dy = Ly / (Ny + 1)
		Dx = sdiff1(Nx) / dx
		Dy = sdiff1(Ny) / dy
		Ax = Dx' * Dx
		Ay = Dy' * Dy
		return Dx, Dy, kron(eye(Ny), Ax) + kron(Ay, eye(Nx))
	end
end

# ╔═╡ c0960c8a-76d0-466a-889b-2179d02518a7
# fg =  plot(p1, p2, p3, layout=(3,1), size=(750, 2250))

# ╔═╡ eacf4f60-956d-4408-b67d-58465a0411dc
# savefig(fg, "ch3_four_attractor.pdf")

# ╔═╡ 23bd055e-599b-407a-9918-d9880d257db9
md"### Appendix: Helper decomposition functions"

# ╔═╡ 3fd9492a-d75e-4603-bafc-4e48fd422385
## We are trying to solve the following equation
## ∇²Uᴴ (x⃗) = -∇ ⋅ F(x⃗)
## We are starting in 2 dimensions
function helmholtz(prob::DiscreteProblem, minval::Int64, maxval::Int64, nsample::Int64)
    dx = (maxval - minval) / (nsample + 1)
    dy = (maxval - minval) / (nsample + 1)
    xs = range(minval, maxval, length=nsample)
    ys = range(minval, maxval, length=nsample)
    ∂x(x, y) = prob.f.f([0.0, 0.0], [x, y], prob.p, nothing)[1]
    ∂y(x, y) = prob.f.f([0.0, 0.0], [x, y], prob.p, nothing)[2]
    # Let's first construct F = [f1, f2]
    f1 = [∂x(x, y) for x in xs, y in ys]
    f2 = [∂y(x, y) for x in xs, y in ys]
    _, _, Δ = laplacian2d(nsample, nsample, maxval - minval, maxval - minval)

    Dx = DiffEqOp.CenteredDifference{1}(1, 2, dx, nsample)
    Dy = DiffEqOp.CenteredDifference{2}(1, 2, dy, nsample)

    bcx = DiffEqOp.Dirichlet0BC{1}(eltype(f1), size(f1))
    bcy = DiffEqOp.Dirichlet0BC{2}(eltype(f2), size(f2))

    rhs =  Dx * bcx * f1 + Dy * bcy *  f2

    u = Δ \ flatshape(rhs)
    xs, ys, reshape(u, size(rhs)...)
end

# ╔═╡ 751ac213-1312-4f0e-bdd8-fa2a10c4554f
begin
	nsample = 512
	dprob = DiscreteProblem(f, [0.0,0.0], (0, 100.), nothing)
	x, y, u_h = helmholtz(dprob, -3, 3, nsample)
end

# ╔═╡ 8101d876-4c69-43e2-8492-59a799707151
begin
	gr()
	fg = surface(x, y, u_h, colorbar=:none,  c=:viridis, aspect_ratio=:equal, xlabel=L"x", ylabel=L"y", zlabel=L"U^H")
end

# ╔═╡ b6ce865e-e1d0-47cd-bbfc-e144459d92a2
p3 = heatmap(x, y, u_h, c=:viridis)

# ╔═╡ 75902190-09c8-41aa-b3ab-ceeaee4ac821

function make_inplace(f::Function)
    function inplace_f(out, x)
        out .= f(x)
    end
    return inplace_f
end

# ╔═╡ b57411c3-a92d-4ce6-82e0-ad0cc4b56f3f
function normal(prob::DiscreteProblem, minval::Int64, maxval::Int64, nsample::Int64, u0::Matrix{Float64})
    dx = (maxval - minval) / (nsample + 1)
    dy = (maxval - minval) / (nsample + 1)

    xs = range(minval, maxval, length=nsample)
    ys = range(minval, maxval, length=nsample)
    ∂x(x, y) = prob.f.f([0.0, 0.0], [x, y], prob.p, nothing)[1]
    ∂y(x, y) = prob.f.f([0.0, 0.0], [x, y], prob.p, nothing)[2]

    f1 = [∂x(x, y) for x in xs, y in ys]
    f2 = [∂y(x, y) for x in xs, y in ys]

    Dx = DiffEqOp.UpwindDifference{1}(1, 2, dx, nsample, 1)
    Dy = DiffEqOp.UpwindDifference{2}(1, 2, dy, nsample, 1)

    bcx = DiffEqOp.Dirichlet0BC{1}(eltype(u0), size(u0))
    bcy = DiffEqOp.Dirichlet0BC{2}(eltype(u0), size(u0))

    function find_my_zero(u)
        dudx = Dx * bcx * u
        dudy = Dy * bcy * u
        return dot(dudx, f1 .+ dudx) + dot(dudy, f2 .+ dudy)
    end

    fmz! = make_inplace(find_my_zero)

    sol = nlsolve(fmz!, u0, show_trace=true, iterations=3000)
    sol.zero
end

# ╔═╡ Cell order:
# ╟─b91fd324-adf3-11eb-36e5-59bf4de5e0f9
# ╠═d4260fcc-784e-4341-ae84-c7e8ab288c4b
# ╠═82272c06-3861-4aed-a08d-bed2e84570fc
# ╟─41f97746-5818-410f-b3cb-e59aa9a5bc99
# ╠═5a0a71b5-0b18-42be-b567-c01feb0d9af6
# ╠═f56451d1-6f67-494b-82b0-2369558d6819
# ╟─a342ae72-b363-4142-929c-78e682871d81
# ╠═7cb95f1f-8a19-484c-a33f-a13348dce332
# ╟─7d1e052e-c9b6-4fa3-b4aa-980e60e53961
# ╠═601e7f8b-efa8-4710-b088-78b825b0454a
# ╠═e3b32bba-b8c8-48ba-8918-c73f1f1e9db9
# ╠═67d5512b-fe45-47d7-9b18-c847a8010577
# ╠═a1d79d06-afdd-4097-94ac-716779173e6c
# ╠═eef51429-8458-421a-8366-35577019e204
# ╠═10c9d45f-d8ab-4acd-88ba-8a80a747e30f
# ╟─9e084100-e5d1-488c-8eaf-769711af1d95
# ╠═dc5e6e4a-93d6-47c3-80ef-8d22d3483e99
# ╟─3481f158-573b-4ca4-9f2f-fa0e4b68edb4
# ╠═25b8d3e0-4402-4e84-91a0-c5d1c4684ac4
# ╟─1870a27b-683b-463a-b5cc-9ac073110091
# ╠═c1b075dc-c4dc-4bd9-bbd9-bbb0119fc54a
# ╠═351e030b-0d6c-42ea-8c5b-9cad2af5c12c
# ╟─8b509897-8b02-4af7-a0bd-990f2cb3df87
# ╠═75bcc971-1a89-4dd4-a5ba-a1397f5db274
# ╟─86019cdf-c7de-4495-ba61-de0852f2b496
# ╠═13347146-e598-4d25-9d92-4b047a9de492
# ╠═6c523287-b8db-4ef7-a08a-f61bf2639278
# ╠═bc84c949-6b26-4e02-84bf-a2093508a0dc
# ╠═6af8e6ef-2444-4c88-ad81-a645b8971527
# ╠═b3a3b34b-ffff-4b35-952e-d844b5830ce7
# ╠═4e4f16b1-7c8a-4435-91f9-6a0d4ff54b22
# ╟─c37e269f-ac80-4bd1-81d0-e5a476e85272
# ╠═164f3c77-5e95-4feb-8783-54c381f4b184
# ╠═751ac213-1312-4f0e-bdd8-fa2a10c4554f
# ╠═8101d876-4c69-43e2-8492-59a799707151
# ╠═b6ce865e-e1d0-47cd-bbfc-e144459d92a2
# ╠═c0960c8a-76d0-466a-889b-2179d02518a7
# ╠═eacf4f60-956d-4408-b67d-58465a0411dc
# ╟─23bd055e-599b-407a-9918-d9880d257db9
# ╟─3fd9492a-d75e-4603-bafc-4e48fd422385
# ╟─75902190-09c8-41aa-b3ab-ceeaee4ac821
# ╟─b57411c3-a92d-4ce6-82e0-ad0cc4b56f3f

ext: jl


text:
using BlockSparseMatrices
using Base.Test

include("test_SparseMatrixBCS.jl")


text:
export AbstractSolverType
export DiscreteSplittingType, HEVISplitting
export getdtmodel

"""
    AbstractSolverType

This is an abstract type representing a generic solver. By
a "solver," we mean an ODE solver together with any potential
implicit solver (linear solvers).
"""
abstract type AbstractSolverType end

"""
    DiscreteSplittingType

This is an abstract type representing a temporal splitting
in the discrete equations. For example, HEVI
(horizontally explicit, vertically implicit) type methods.
"""
abstract type DiscreteSplittingType end

"""
    SlowFastSplitting

The tendency is treated using a standard slow-fast splitting,
where the fast processes are purely related to acoustic/gravity
waves (ideal for a typical multirate or MIS method).
This splitting does not take into account the geometric
direction (vertical vs horizontal).
"""
struct SlowFastSplitting <: DiscreteSplittingType end

"""
    HEVISplitting

HEVI (horizontally explicit, vertically implicit) type method,
where vertical acoustic waves are treated implicitly. All other
dynamics are treated explicitly.

Note: Can potentially imagine several different types of
HEVI splittings (for example, include vertical momentum and/or
diffusion)
"""
struct HEVISplitting <: DiscreteSplittingType end

"""
    getdtmodel(ode_solver::AbstractSolverType, bl)

A function which returns a model representing the dynamics
with the most restrictive time-stepping requirements.
"""
getdtmodel(ode_solver::AbstractSolverType, bl) =
    throw(MethodError(solversetup, (ode_solver, bl)),)

"""
    solversetup(
        ode_solver::AbstractSolverType,
        dg,
        Q,
        dt,
        t0,
        diffusion_direction,
    )

A function which returns a configured ODE solver.
"""
solversetup(
    ode_solver::AbstractSolverType,
    dg,
    Q,
    dt,
    t0,
    diffusion_direction,
) = throw(MethodError(
    solversetup,
    (ode_solver, dg, Q, dt, t0, diffusion_direction),
))

include("ExplicitSolverType.jl")
include("ImplicitSolverType.jl")
include("IMEXSolverType.jl")
include("HEVISolverType.jl")
include("MultirateSolverType.jl")
include("MISSolverType.jl")
include("SplitExplicitSolverType.jl")

DefaultSolverType = IMEXSolverType
export DefaultSolverType


text:
# # Plot Recipes

@nospecialize

"""
    _process_plotrecipes!(plt, kw_list)

Grab the first in line to be processed and pass it through `apply_recipe` to generate a
list of `RecipeData` objects.
If we applied a "plot recipe" without error, then add the returned datalist's KWs,
otherwise we just add the original KW.
"""
function _process_plotrecipes!(plt, kw_list)
    still_to_process = kw_list
    kw_list = KW[]
    while !isempty(still_to_process)
        next_kw = popfirst!(still_to_process)
        _process_plotrecipe(plt, next_kw, kw_list, still_to_process)
    end
    return kw_list
end

function _process_plotrecipe(plt, kw, kw_list, still_to_process)
    if !isa(get(kw, :seriestype, nothing), Symbol)
        # seriestype was never set, or it's not a Symbol, so it can't be a plot recipe
        push!(kw_list, kw)
        return
    end
    try
        st = kw[:seriestype]
        st = kw[:seriestype] = type_alias(plt, st)
        datalist = RecipesBase.apply_recipe(kw, Val{st}, plt)
        warn_on_recipe_aliases!(plt, datalist, :plot, st)
        for data in datalist
            preprocess_attributes!(plt, data.plotattributes)
            if data.plotattributes[:seriestype] == st
                error("Plot recipe $st returned the same seriestype: $(data.plotattributes)")
            end
            push!(still_to_process, data.plotattributes)
        end
    catch err
        if err isa MethodError && occursin("plot recipe", err.args)
            push!(kw_list, kw)
        else
            rethrow()
        end
    end
    return
end

@specialize


text:
function get_additional_state(obj::AdditionalStateProvider, arg0::SpacecraftState)
    return jcall(obj, "getAdditionalState", Vector{jdouble}, (SpacecraftState,), arg0)
end

function get_name(obj::AdditionalStateProvider)
    return jcall(obj, "getName", JString, ())
end


text:
# targets for phase in / phase out of policies

pinft = 0.4.*pinf_o # infection prob at meeting

socialmaxyyt = 2 # max no. people met outside firm young_young
socialmaxoyt = 1 # max no. people met outside firm old_young
socialmaxoot = 0.5 # max no. people met outside firm old_old
socialmaxyot = 0.5 # max no. people met outside firm old_old

phomeofficet = 1

pshopt = pshopt = [[1,1,1],[1,1,1]] - alpha_lockdown*[[0.15,0.5,0],[0.15,0.5,0]]
pshop = pshopt = [[1,1,1],[1,1,1]] - alpha_lockdown*[[0.15,0.5,0],[0.15,0.5,0]]

workmeett = [4,5,8,2] # max number of people met at work: man / ser / food /public
shopmeett = [5,20,10] # max number of people met at each shopping instant  man / ser / food

# direct change in variable through policy
tauadj = 0

genhomeoffice = true

ptravel = 0

bailoutprogram = true # initially no bailouts

shorttimeprogram = true # initially no short term program


text:
# code to calculate the first derivative of a function using finite differences
# function to be evaluated is a sin function
#Initial Parameters
using Statistics
xmax=10.0
nx=200
dx=xmax/(nx-1)
x=LinRange(0,xmax,nx)
#Initialization of the sin function
l=20*dx #wavelength
k=2*pi/l
f=sin.(k.*x)
using PyPlot
figure(1)
plot(x,f)
xlabel("X")
ylabel("Amplitude")
title("Sin function")
#Initialize the analytical and numerical derivatives
ader=zeros(nx)
nder=zeros(nx)
#Numerical derivative
for i in 2:nx-1
    nder[i]=(f[i+1]-f[i-1])/(2*dx)
end
#Analytical derivative
ader=k*cos.(k.*x)
ader[1]=0
ader[nx]=0
#Calculate RMS error
t1=(ader.-nder)
t2=t1.^2
t3=mean(t2)
rms=sqrt.(t3)
difference=ader-nder
figure(2)
plot(x,ader)
plot(x,nder)
plot(x,difference)
xlabel("X")
ylabel("Amplitude")
title("Derivatives rms error: $rms")
gcf()


text:
using Pkg; Pkg.activate(".")

import Counter: Domain

using .Domain

item_name = "ar"

item_nbr = 1000

item = Item(item_name, item_nbr)

@info item.name * " has been created on " * string(item.created) * " with the value " * string(item.nbr)


text:
module GameOfLife

import REPL
using REPL.TerminalMenus
export Gol, game, game2, glidergun

mutable struct Gol
    state::Matrix{Int8}
    isperiodic::Bool
    generations::Int
end

function Gol(n::Int,
             m::Int = n;
             isperiodic::Bool = true,
             generations::Int = 100)
    state = rand(Int8[0,1], n, m)

    return Gol(state, isperiodic, generations)
end

"count the number of live neighbors"
function determine_neighbor_state(d::Gol, i::Integer, j::Integer)
    n, m = size(d.state)
    neighbor = [(p,q) for p in -1:1 for q in -1:1 if !(p == 0 && q == 0)]
    total = Int8(0)
    if d.isperiodic # periodic
        for idx in 1:8
            total += d.state[mod1(i + neighbor[idx][1], n), mod1(j + neighbor[idx][2], m)]
        end
    else # non-periodic
        for idx in 1:8
            if !(i + neighbor[idx][1] <= 0 || j + neighbor[idx][2] <= 0 || i + neighbor[idx][1] > n || j + neighbor[idx][2] > m)
                total += d.state[i + neighbor[idx][1], j + neighbor[idx][2]]
            end
        end
    end

    return total
end

"update cells"
function next_generation(d::Gol)
    n, m = size(d.state)
    next_state = zeros(Int8, n, m)
    for col in 1:m
        for row in 1:n
            if d.state[row, col] == one(d.state[row, col]) # live
                if determine_neighbor_state(d, row, col) in Int8[2, 3]
                    next_state[row, col] = Int8(1)
                end
            else # dead
                if determine_neighbor_state(d, row, col) == Int8(3)
                    next_state[row, col] = Int8(1)
                end
            end
        end
    end
    return next_state
end

"show state"
function Base.show(d::Gol)
    n, m = size(d.state)
    for i in 1:n
        for j in 1:m
            if d.state[i,j] == zero(d.state[i,j])
                print(" ")
            else
                printstyled("X", color=:light_blue, bold=true)
            end
        end
        println()
    end
end

"""
    game(d::Gol; tstep=0.1)
simulation start
"""
function game(d::Gol; tstep=0.1)
    sleep(0.1)
    show(d);
    n = length(string(d.generations))
    for iter in 1:d.generations
       d.state = next_generation(d)
       sleep(tstep)
       run(`clear`)
       show(d)
       print("\nGeneration: ", lpad(iter, n))
       sum(d.state) == 0 && break
    end
end

"""
    game(;t=0.1, gen=100)
"""
function game(;t=0.1, gen=100)
    row, col = displaysize(stdout)
    d = Gol(row, col, isperiodic=true, generations=gen)
    game(d, tstep=t)
end

"""
    glidergun(n::Integer = 20, m::Integer = 36; generations=250)
glider gun sample
"""
function glidergun(n::Integer = 20, m::Integer = 36; generations=250)
    d = Gol(n, m, isperiodic=false, generations=generations)
    d.state = zeros(Int8, n, m)
    gliderpos = Tuple{Int8,Int8}[(5,1),(6,1),(5,2),(6,2),(5,11),(6,11),(7,11),(4,12),(8,12),(3,13),(9,13),(3,14),(9,14),(6,15),(4,16),(8,16),(5,17),(6,17),(7,17),(6,18),(3,21),(4,21),(5,21),(3,22),(4,22),(5,22),(2,23),(6,23),(1,25),(2,25),(6,25),(7,25),(3,35),(4,35),(3,36),(4,36)]
    for pos in gliderpos
        d.state[pos...] = Int8(1)
    end
    return d
end

boardsize = ["(20, 20)", "(20, 40)", "(30, 30)", "(30, 40)", "(40, 40)", "(40, 80)", "(100, 200)"]
boardsize_dict = Dict(zip(1:length(boardsize), eval.(Meta.parse.(boardsize))))
menu = RadioMenu(boardsize, pagesize=4)
function game2()
    choice = request("Choose your favorite field size (row, column):", menu)
    col, row = boardsize_dict[choice]
    println()
    periodic_condition = request("Choose boundary condition:",
                                RadioMenu(["Periodic", "Non periodic"],
                                pagesize=2))
    condition = ifelse(isone(periodic_condition), true, false)

    print("\nThe number of generations: ")
    num_generation = Meta.parse(readline())
    @assert isinteger(num_generation)

    game(Gol(col, row, isperiodic=condition, generations=num_generation))
    return nothing
end

end # module


text:
# Plot stuff
import Plots: plot!, plot3d!
using RegionTrees

function plot!(octree::Octree)
    # plt = plot(legend=nothing)
    plt = scatter(octree.X[1,:], octree.X[2,:], markersize=2, marker=:x, legend=nothing)
    for leaf in allleaves(octree.root)
       v = hcat(collect(vertices(leaf.boundary))...)
       plot!(plt, v[1,[1,2,4,3,1]], v[2,[1,2,4,3,1]], color=:black)
   end
   plt
end

function plot3d!(octree::Octree)
    plt = plot(legend=nothing)
    for leaf in allleaves(octree.root)
       v = hcat(collect(vertices(leaf.boundary))...)
       facets = [[1,2,4,3,1], [1,5,7,3], [2,6,8,4], [7,8], [5,6]]
       for facet in facets
           plot!(plt, v[1, facet], v[2, facet], v[3, facet])
       end
   end
   plt
end


text:
function mac(v1,v2)

	# This funcion compute the Modal Assurance Criterion (MAC) between
	# two vectors (real or complex vectors)

	# first we convert v1 and v2 in row matrices (1,n)

	n1 = length(v1)
	v1 = reshape(v1,1,n1)
	n2 = length(v2)
	v2 = reshape(v2,1,n2)

	# julia transpose is hermitic: transpose + conjugate
	m = (abs.(v1*v2'))^2/((v1*v1')*(v2*v2')) # m is an array with one complex element
	m = Real(m[1])

	return(m)

end


text:
# Autogenerated wrapper script for simdjson_jll for powerpc64le-linux-gnu-cxx03
export libsimdjson

JLLWrappers.@generate_wrapper_header("simdjson")
JLLWrappers.@declare_library_product(libsimdjson, "libsimdjson.so.8")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libsimdjson,
        "lib/libsimdjson.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# # Custom Cone Constraint
#
# This example demonstrate how the user can extend `COSMO` with his own custom convex cones and use them in constraints. To make things easy, we will implement the simple cone of **Nonpositives**, i.e. $\mathbb{R}_{-}^n = \{x \in   \mathbb{R}^n \mid x_i \leq 0\}$.
# ## Cone definition and projection function
# We start by creating a concrete subtype of our cone interface-type: `COSMO.AbstractConvexCone{T}`, where `T` is a parameter for the floating-point precision. The only field that is strictly required for our cone definition is `dim`, which is used by the solver to store information about the cone dimension. However, you can also use the object to store other information, or allocate the workspace necessary for the projection step.
using COSMO, LinearAlgebra, SparseArrays, Test

## define new cone type
struct Nonpositives{T} <: COSMO.AbstractConvexSet{T}
    dim::Int64
end
#-
# Next, we define a projection method for `Nonpositives{T}` that takes a vector `x` and projects it onto our custom cone. For the cone of nonpositive vectors the projection simply sets all positive elements of `x` to zero:

function COSMO.project!(x::AbstractVector{T}, C::Nonpositives{T}) where {T <: AbstractFloat}
    x .= min.(x, zero(T))
end

#-
# This is all that is required to get a basic constraint working. Let's test our new cone/constraint by solving the following LP:
# $$
# \begin{array}{ll} \text{maximize} & x_1 + x_2 + x_3\\
# \text{subject to} &  x_1 \leq 3 \\
#                   &  x_2 \leq 2 \\
#                   &  x_1 + x_3 =  5
# \end{array}
# $$
# We define our decision vector $x = (x_1, x_2, x_3)$ and setup the problem:
n = 3
P = spzeros(n, n)
q = -ones(n)

## A1 * x + b1 <= 0
A1 = diagm(0 => ones(2))
b1 = [-3.; -2.]
c1 = COSMO.Constraint(A1, b1, Nonpositives, n, 1:2); #here we use our new cone

## x_1 + x_3 == 5
A2 = [1. 0 1.]
b2 = [-5.]
c2 = COSMO.Constraint(A2, b2, COSMO.ZeroSet);

## assemple and solve
model = COSMO.Model();
assemble!(model, P, q, [c1; c2]);
res = COSMO.optimize!(model);

#-
res.x

#-
obj_val = -dot(q, res.x)
# The problem was solved using constraints involving our new cone and yields the expected result. Next, we want to add the ability to detect infeasible problems.
#
# ## Support infeasibility detection
# If no further information about the new cone is provided, the infeasibility detection is disabled. However, by defining the two additional methods `in_dual` and `in_pol_recc` we can support infeasibility detection. More information on how infeasible problems are detected in `COSMO` can be found here \[1\]. We will have to add the ability for the solver to check whether a vector is in the [dual cone](https://en.wikipedia.org/wiki/Dual_cone_and_polar_cone) $\mathcal{K}^*$ of our new cone $\mathcal{K}$ and whether a vector is in the recession cone of the [polar cone](https://en.wikipedia.org/wiki/Dual_cone_and_polar_cone) ${\mathcal{K}^\circ}^\infty$ of our cone. The first function is used to check for primal infeasibility, while the second function  is used to check for dual infeasibility. Luckily, for the Nonpositives-cone this is straightforward:
#
# $$
# {\mathbb{R}_{-}^n}^* = \mathbb{R}_{-}^n, \quad {{\mathbb{R}_{-}^n}^\circ}^\infty = \mathbb{R}_{+}^n.
# $$

function COSMO.in_dual(x::AbstractVector{T}, C::Nonpositives{T}, tol::T) where {T <: AbstractFloat}
	return !any( x-> (x > -tol), x)
end

function COSMO.in_pol_recc(x::AbstractVector{T}, C::Nonpositives{T}, tol::T) where {T <: AbstractFloat}
	return !any( x-> (x < tol), x)
end
# Notice that for numerical reasons, we give the check-functions a bit of slack using the tolerance parameter `tol`.
# To test the infeasibility detection, we will attempt to use our new cone to solve the clearly dual infeasible problem:
# $$
# \begin{array}{ll} \text{minimize} & x\\
# \text{subject to} &  x\leq 3
# \end{array}
# $$
n = 1
P = spzeros(n, n)
q = [1.]

##  x <= 3 <=> x - 3 in Nonpositives
Ai = [1.]
bi = [-3.]
ci = COSMO.Constraint(Ai, bi, Nonpositives);

model = COSMO.Model();
assemble!(model, P, q, [ci]);
res = COSMO.optimize!(model);
#-
@test res.status == :Dual_infeasible

# which is what we would expect.
# ## Using new cone with JuMP
# The nice composibility of Julia-code allows us to use our new cone definition even when the problem is modelled with `JuMP`. For this to work, we define a concrete subtype of `MOI.AbstractSet`, that JuMP can use to build a constraint.
using MathOptInterface, JuMP
import Base.copy
const MOI = MathOptInterface

struct NonPos <: MOI.AbstractVectorSet
	dimension::Int
end
Base.copy(set::NonPos) = set #this is needed for MOI
# Next, we tell `COSMO` that whenever `JuMP` wants to solve a problem with a `NonPos`-cone it should translate it into a `Nonpositives` constraint and we also tell `MOI` that we now support constraints of this new mysterious cone `NonPos`.
function COSMO.processSet!(b::Vector{T}, rows::UnitRange{Int}, cs, s::NonPos) where {T <: AbstractFloat}
    push!(cs, Nonpositives{T}(length(rows)))
    nothing
end
#tell MOI, that COSMO supports constraints with this new cone
MOI.supports_constraint(optimizer::COSMO.Optimizer, ::Type{<:Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{Float64}}}, ::Type{NonPos}) = true

# We should now be able to model the LP from above in `JuMP` and solve it internally using our `Nonpositives` cone and projection function.
model = JuMP.Model(COSMO.Optimizer);
@variable(model, x[1:3]);
@objective(model, Max, x[1] + x[2] + x[3]);
@constraint(model, A1 * x[1:2] .+ b1 in NonPos(2));
@constraint(model, x[1] + x[3] == 5);
JuMP.optimize!(model)

#-
x_opt = JuMP.value.(x)
# You can see in the solver output that indeed a set `Nonpositives` of `dim: 2` was used.
#
# The discussed cone of Nonpositives is of course trivial, but the ability to define new cones and constraints for `COSMO` can be very powerful to model complex problems.

# ## References
# [1] Garstka et al. - COSMO: A conic operator splitting method for convex conic problems (2020)


text:
abstract type Perturbation end

struct ThreePhaseFault
    time::Float64
    Ybus::SparseMatrixCSC{Complex{Float64}, Int64}

end


text:
# ===============================================================
# Discretize using the correction hull of the matrix exponential
# ===============================================================

"""
    CorrectionHull{EM} <: AbstractApproximationModel

Discretization using the correction hull of the matrix exponential.

### Fields

- `exp`   -- exponentiation method
- `order` -- order of the Taylor series expansion of the matrix exponential

### Algorithm

For the homogeneous case, this method implements the transformation:

```math
Ω_0 = CH(X_0, e^{Aδ}  X_0) ⊕ FX_0
```
where ``F`` is the correction (interval) matrix.

For the inhomogeneous case, ``x' = Ax + u``,  ``x ∈ X, u ∈ U``, implements
``Ω_0 = CH(X_0, exp(Aδ)  X0) ⊕ FX0`` where ``F`` is the correction (interval) matrix.

In both cases, if ``A`` is an interval matrix, the exponential is overapproximated
using methods from `IntervalMatrices.jl`.
"""
struct CorrectionHull{EM} <: AbstractApproximationModel
   order::Int
   exp::EM
end

# convenience constructor using symbols
function CorrectionHull(; order::Int=10, exp=IntervalExpAlg(order))
    return CorrectionHull(order, _alias(exp))
end

function Base.show(io::IO, alg::CorrectionHull)
    print(io, "`CorrectionHull` approximation model with: \n")
    print(io, "    - exponentiation method: $(alg.exp) \n")
    print(io, "    - order: $(alg.order)\n")
end

Base.show(io::IO, m::MIME"text/plain", alg::CorrectionHull) = print(io, alg)

# -----------------------------------------------------------------
# Correction hull: homogeneous case x' = Ax, x in X
# -----------------------------------------------------------------

function discretize(ivp::IVP{<:CLCS, <:LazySet}, δ, alg::CorrectionHull)
    A = state_matrix(ivp)
    X0 = initial_state(ivp)
    X = stateset(ivp)

    # compute exp(A*δ) * X0
    Φ = _exp(A, δ, alg.exp)

    # compute Ω0 = CH(X0, ΦX0) ⊕ FX0
    Ω0 = _discretize_chull(A, Φ, X0, δ, alg)

    Sdis = ConstrainedLinearDiscreteSystem(Φ, X)
    return InitialValueProblem(Sdis, Ω0)
end

function _discretize_chull(A, Φ::IntervalMatrix, X0, δ, alg, P=nothing)
    X0z = _convert_or_overapproximate(X0, Zonotope)
    Y = _overapproximate(Φ * X0z, Zonotope)
    if !isnothing(P)
        Y = minkowski_sum(Y, P)
    end

    H = overapproximate(CH(X0z, Y), Zonotope)
    F = correction_hull(A, δ, alg.order)
    R = _overapproximate(F * X0z, Zonotope)
    Ω0 = minkowski_sum(H, R)
    return Ω0
end

# F(δ) without the E(δ) correction term
function _correction_hull_without_E(A, δ, p)
    timeint(δ, i) = interval((i^(-i / (i-1)) - i^(-1 / (i-1))) * δ^i, 0)
    F = sum(map(x -> timeint(δ, i) * x, A^i / factorial(i)) for i in 2:p)
    return IntervalMatrix(F)
end

function _discretize_chull(A, Φ::AbstractMatrix, X0, δ, alg)
    X0z = _convert_or_overapproximate(X0, Zonotope)
    Y = linear_map(Φ, X0z)

    H = overapproximate(CH(X0z, Y), Zonotope)
    F = _correction_hull_without_E(A, δ, alg.order)
    R = _overapproximate(F * X0z, Zonotope)

    Ω0 = minkowski_sum(H, R)
    return Ω0
end

# -----------------------------------------------------------------
# Correction hull: inhomogeneous case x' = Ax + u, x in X, u ∈ U
# -----------------------------------------------------------------
function discretize(ivp::IVP{<:CLCCS, <:LazySet}, δ, alg::CorrectionHull)
    A = state_matrix(ivp)
    X0 = initial_state(ivp)
    X = stateset(ivp)
    U = next_set(inputset(ivp), 1) # inputset(ivp)

    # here U is an interval matrix map of a lazyset, TODO refactor / dispatch
    if isa(U, LinearMap)
        Uz = _convert_or_overapproximate(Zonotope, LazySets.set(U))
        B = matrix(U) # TODO remove this case, since preprocessing normalize(ivp) would make B the identity
        Uz = isinterval(B) ? _overapproximate(B * Uz, Zonotope) : linear_map(B, Uz)
    else # LazySet
        Uz = _convert_or_overapproximate(Zonotope, U)
    end

    Φ = _exp(A, δ, alg.exp)

    A_interval = _interval_matrix(A)

    origin_not_contained_in_U = zeros(dim(U)) ∉ Uz

    if origin_not_contained_in_U
        # shift U to origin
        u = center(Uz)
        Uz = Zonotope(zeros(dim(U)), genmat(Uz))
    end

    # note: this evaluates Φ₁(A, δ) using Taylor expansion of the interval matrix
    Cδ = _Cδ(A_interval, δ, alg.order)

    if origin_not_contained_in_U
        # compute C(δ) * u
        Pu = _convert_intervals_to_box_with_vectors(Cδ * u)
    else
        Pu = nothing
    end

    Ω0 = _discretize_chull(A_interval, Φ, X0, δ, alg, Pu)

    if origin_not_contained_in_U
        # apply another correction hull for the shifted U
        F = input_correction(A_interval, δ, alg.order)
        Fu = _convert_intervals_to_box_with_vectors(F * u)
        Ω0 = minkowski_sum(Ω0, Fu)
    end
    # Ω0 = _apply_setops(Ω0, alg.setops) # TODO requires to add `setops` field to the struct

    # compute C(δ) * U
    Ud = _overapproximate(Cδ * Uz, Zonotope)

    # add inputs to Ω0
    Ω0 = minkowski_sum(Ω0, Ud)

    B = Matrix(IdentityMultiple(one(eltype(A)), size(A, 1)))
    Sdis = ConstrainedLinearControlDiscreteSystem(Φ, B, X, Ud)
    return InitialValueProblem(Sdis, Ω0)
end

# convert a vector of intervals to a hyperrectangle with normal vector type
function _convert_intervals_to_box_with_vectors(intervals)
    # convert to Hyperrectangle
    H = convert(Hyperrectangle, IntervalBox(intervals))
    # remove static vectors (they scale poorly)
    H = Hyperrectangle(Vector(center(H)), Vector(radius_hyperrectangle(H)))
    return H
end


text:
# Author: Mathias Louboutin, mlouboutin3@gatech.edu
# Date: December 2020

module SlimOptim

using Printf, LinearAlgebra, LineSearches

import LineSearches: BackTracking, HagerZhang, Static, MoreThuente, StrongWolfe
export BackTracking, HagerZhang, Static, MoreThuente, StrongWolfe

export pqn, pqn_options
export spg, spg_options
export bregman, bregman_options
#############################################################################
# Optimization algorithms
include("utils.jl") # common functions
include("linesearches.jl")   # LineSearches.jl
include("SPGSlim.jl")   # minConf_SPG
include("PQNSlim.jl")   # minConf_PQN
include("bregman.jl")   # Linearized bregman iterations

function checkls(ls)
    is_ls = (typeof(ls) <: Union{BackTracking, HagerZhang, Static, MoreThuente, StrongWolfe})
    !is_ls && throw(ArgumentError("Unrecognized linsearch input, only LineSearches.jl's[BackTracking, HagerZhang, Static, MoreThuente, StrongWolfe] are supported"))
end

end # module


text:
# common Base functionality

allequal(x) = true
allequal(x, y, z...) = x == y && allequal(y, z...)
function Base.map!(f, y::AbstractGPUArray, xs::AbstractGPUArray...)
    @assert allequal(size.((y, xs...))...)
    return y .= f.(xs...)
end
function Base.map(f, y::AbstractGPUArray, xs::AbstractGPUArray...)
    @assert allequal(size.((y, xs...))...)
    return f.(y, xs...)
end

# Break ambiguities with base
Base.map!(f, y::AbstractGPUArray) =
    invoke(map!, Tuple{Any,AbstractGPUArray,Vararg{AbstractGPUArray}}, f, y)
Base.map!(f, y::AbstractGPUArray, x::AbstractGPUArray) =
    invoke(map!, Tuple{Any,AbstractGPUArray, Vararg{AbstractGPUArray}}, f, y, x)
Base.map!(f, y::AbstractGPUArray, x1::AbstractGPUArray, x2::AbstractGPUArray) =
    invoke(map!, Tuple{Any,AbstractGPUArray, Vararg{AbstractGPUArray}}, f, y, x1, x2)

# Base functions that are sadly not fit for the the GPU yet (they only work for Int64)
Base.@pure @inline function gpu_ind2sub(A::AbstractArray, ind::T) where T
    _ind2sub(size(A), ind - T(1))
end
Base.@pure @inline function gpu_ind2sub(dims::NTuple{N}, ind::T) where {N, T}
    _ind2sub(NTuple{N, T}(dims), ind - T(1))
end
Base.@pure @inline _ind2sub(::Tuple{}, ind::T) where {T} = (ind + T(1),)
Base.@pure @inline function _ind2sub(indslast::NTuple{1}, ind::T) where T
    ((ind + T(1)),)
end
Base.@pure @inline function _ind2sub(inds, ind::T) where T
    r1 = inds[1]
    indnext = div(ind, r1)
    f = T(1); l = r1
    (ind-l*indnext+f, _ind2sub(Base.tail(inds), indnext)...)
end

Base.@pure function gpu_sub2ind(dims::NTuple{N}, I::NTuple{N2, T}) where {N, N2, T}
    Base.@_inline_meta
    _sub2ind(NTuple{N, T}(dims), T(1), T(1), I...)
end
_sub2ind(x, L, ind) = ind
function _sub2ind(::Tuple{}, L, ind, i::T, I::T...) where T
    Base.@_inline_meta
    ind + (i - T(1)) * L
end
function _sub2ind(inds, L, ind, i::IT, I::IT...) where IT
    Base.@_inline_meta
    r1 = inds[1]
    _sub2ind(Base.tail(inds), L * r1, ind + (i - IT(1)) * L, I...)
end

# This is pretty ugly, but I feel bad to add those to device arrays, since
# we're never bound checking... So getindex(a::AbstractGPUVector, 10, 10) would silently go unnoticed
# we need this here for easier implementation of repeat
@inline Base.@propagate_inbounds getidx_2d1d(x::AbstractVector, i, j) = x[i]
@inline Base.@propagate_inbounds getidx_2d1d(x::AbstractMatrix, i, j) = x[i, j]

function Base.repeat(a::AbstractGPUVecOrMat, m::Int, n::Int = 1)
    o, p = size(a, 1), size(a, 2)
    b = similar(a, o*m, p*n)
    gpu_call(a, (b, a, o, p, m, n), n) do state, b, a, o, p, m, n
        j = linear_index(state)
        j > n && return
        d = (j - 1) * p + 1
        @inbounds for i in 1:m
            c = (i - 1) * o + 1
            for r in 1:p
                for k in 1:o
                    b[k - 1 + c, r - 1 + d] = getidx_2d1d(a, k, r)
                end
            end
        end
        return
    end
    return b
end

function Base.repeat(a::AbstractGPUVector, m::Int)
    o = length(a)
    b = similar(a, o*m)
    gpu_call(a, (b, a, o, m), m) do state, b, a, o, m
        i = linear_index(state)
        i > m && return
        c = (i - 1)*o + 1
        @inbounds for i in 1:o
            b[c + i - 1] = a[i]
        end
        return
    end
    return b
end


text:
using PyPlot
using JLD
include("model_tyc3559_norv.jl")
include("model_tyc3559_v2.jl")
include("transit_orb.jl")
include("occultquad.jl")
include("kepler.jl")
include("compute_likelihood_save.jl")
include("lorentz_likelihood_hermitian_band_init.jl")
include("lorentz_likelihood_hermitian_band_save.jl")
include("bandec.jl")
include("banbks.jl")

#function fit_tyc3559_norv()

data = readdlm("tyc3559_pdc.dat")
time_phot  = vec(data[:,1])-54900.0
fflat = vec(data[:,2])
ntphot = 65266
#ntphot = 100
#time_phot = time_phot[1:ntphot]
#fflat = fflat[1:ntphot]

#param=[pbest,87.0,.023.0,.857.0*pbest,.5,0.0,1.0+6e-5,5.0,0.0,0.0,0e-5,0e-5,-1.3e-4,0.0,2450.0,-1334.0]
param_phot = vec(readdlm("param_best_norv.dat"))
pname1 = ["Period [day]","Inc [deg]","R_2/R_1","t0 [JD-2459000]","u1","u2","f0","a/R_2","ecc","omega","a1c","a1s","a2c","a2s","f_3/(f_1+f_2)","f_2/(f_1+f_2)"]
param_phot1=zeros(16)
for i=1:16
  param_phot1[i]=param_phot[i]
end

ntime=length(time_phot)
t_first = minimum(time_phot)
t_last = maximum(time_phot)
p0 = param_phot[1]
t00 = param_phot[4]
nfirst = floor((t_first-t00)/p0)
nlast = ceil((t_last-t00)/p0)
tn = nlast*p0+t00
aonr = param_phot[8]
u1 = param_phot[5]
u2 = param_phot[6]
r2onr1 = param_phot[3]
inc = param_phot[2]*pi/180.0
b = cos(inc)*aonr
tdur = 1.0/aonr*p0/pi*sqrt(1-b^2)
x = (1-sqrt(1-b^2))
depth = (r2onr1)^2*(1.0-u1*x-u2*x^2)/(1.0-u1/3-u2/6)
q1 = (u1+u2)^2
q2 = u1/(2*(u1+u2))
f0 = param_phot[7]

# Convert to new set of variables:
pname = ["t0 [JD - 24590000]","tN [JD-2459000]","b","depth","T_dur","q1","q2","f0","ecc","omega","a1c","a1s","a2c","a2s","f_3/(f_1+f_2)","f_2/(f_1+f_2)"]
param_phot[1] = t00
param_phot[2] = tn
param_phot[3] = b
param_phot[4] = depth
param_phot[5] = tdur
param_phot[6] = q1
param_phot[7] = q2
param_phot[8] = f0

#time_phot=reform(phasebin[0,*])
dt = median(time_phot[2:ntime]-time_phot[1:ntime-1])
flux_phot=fflat[1:ntime]

scattermin = 0.00040507201
weight_phot=ones(ntime)./scattermin^2
weight=weight_phot
flux=flux_phot
model1 = model_tyc3559_norv(time_phot,param_phot1)
model2 = model_tyc3559_v2(time_phot,param_phot)
@time model = model_tyc3559_v2(time_phot,param_phot)
chisquare=sum((flux-model).^2.*weight)
println("Chi-square: ",chisquare)
println("length(model1): ",length(model1))
println("length(model2): ",length(model2))
println("max(abs(model1-model2)): ",maximum(abs(model1-model2)))
plot(mod(time_phot,p0),model1)
plot(mod(time_phot,p0),model2)
#read(STDIN,Char)

# Now compute correlated likelihood:
w = 0.03027 * 1.6467e-7
kernel_param = [w,1.0428542*1.6467e-7, -0.38361831*1.6467e-7, 0.30345984/2*1.6467e-7, 0.1229159/dt, 0.48922908/dt, 0.09086397/dt, 2pi/12.203317/dt]
# Set kernel to a delta function:
#kernel_param = [w,1.0428542*1.6467e-7, -0.38361831*1.6467e-7, 0.30345984/2*1.6467e-7, 1e5 , 1e5 , 1e5 , 0.]
param_tot = [param_phot;kernel_param]
w = param_tot[17]
alpha = [param_tot[18],param_tot[19],param_tot[20]/2,param_tot[20]/2]
p=length(alpha)
beta = [complex(param_tot[21],0.0),complex(param_tot[22],0.0),complex(param_tot[23],param_tot[24]),complex(param_tot[23],-param_tot[24])]
aex,al_small,indx,logdeta = lorentz_likelihood_hermitian_band_init(alpha,beta,w,time_phot)
#Profile.clear()
#for i=1:100
#  @profile log_like = compute_likelihood_save(time_phot,param_tot,flux,aex,al_small,indx,logdeta,p)
#end
#Profile.print()
#read(STDIN,Char)
log_like = compute_likelihood_save(time_phot,param_tot,flux,aex,al_small,indx,logdeta,p)

nparam = 24
ivary = [1,2,3,4,5,6,7,8,11,12,13,15,16]
nvary = length(ivary)

# Run a Markov chain:
#errors = [0.0001,0.2,0.001,0.01,0.001,0.1,1e-6,0.1,0.,0.,1e-7,1e-7,1e-7,1e-7,1e-3,1e-7,0.,0.,0.,0.,0.,0.,0.,0.]
errors = [0.001,0.001,0.01,0.00001,0.01,0.1,0.1,1e-6,0.,0.,1e-7,1e-7,1e-7,0.,1e-2,1e-7,0.,0.,0.,0.,0.,0.,0.,0.]
#pname = ["t0 [JD - 24590000]","tN [JD-2459000]","b","depth","T_dur","q1","q2","f0","ecc","omega","a1c","a1s","a2c","a2s","f_3/(f_1+f_2)","f_2/(f_1+f_2)"]

nwalkers = nvary * 2
nsteps = 10000
#nsteps = 100
# Set up arrays to hold the results:
par_mcmc = zeros(nwalkers,nsteps,nparam)
log_like_mcmc = zeros(nwalkers,nsteps)
#logprior_mcmc = zeros(nwalkers,nsteps)
# Initialize walkers:
par_trial = param_tot
log_like_best = log_like
for j=1:nwalkers
# Select from within uncertainties:
  log_like_trial = -1e100
# Only initiate models with reasonable chi-square values:
  while log_like_trial < log_like_best - 1000.0
    par_trial = param_tot + errors.*randn(nparam)
    aonr2 = (1./par_trial[5]*p0/pi)^2*(1-par_trial[3]^2)
    if (par_trial[6] < 0) || (par_trial[6] > 1) || (par_trial[7] < 0) || (par_trial[7] > 1) || (par_trial[4] < 0) || (par_trial[5] < 0) || (abs(par_trial[3]) > 1) || (par_trial[3]^2 > aonr2) then
      log_like_trial = -1e100
   else
      log_like_trial = compute_likelihood_save(time_phot,par_trial,flux,aex,al_small,indx,logdeta,p)
   end
#    println("log_like_trial: ",log_like_trial," params: ",par_trial)
  end
  log_like_mcmc[j,1]=log_like_trial
  par_mcmc[j,1,:]=par_trial
#  logprior_mcmc[j,1] = log_prior+log_det
  println("Success: ",par_trial,log_like_trial)
end

# Initialize scale length & acceptance counter:
ascale = 2.0
accept = 0
tic()
# Next, loop over steps in markov chain:
for i=2:nsteps
  for j=1:nwalkers
    ipartner = j
# Choose another walker to 'breed' a trial step with:
    while ipartner == j
      ipartner = ceil(Int,rand()*nwalkers)
    end
# Now, choose a new set of parameters for the trial step:
    z=(rand()*(sqrt(ascale)-1.0/sqrt(ascale))+1.0/sqrt(ascale))^2
    par_trial = vec(par_mcmc[j,i-1,:])
    for itrial=1:nvary
      par_trial[ivary[itrial]]=z*par_mcmc[j,i-1,ivary[itrial]]+(1.0-z)*par_mcmc[ipartner,i-1,ivary[itrial]]
    end
# Compute model & chi-square:
#    log_prior,log_det,log_like_trial = compute_likelihood(time_phot,par_trial,flux)
#    log_like_trial = compute_likelihood(time_phot,par_trial,flux)
    aonr2 = (1./par_trial[5]*p0/pi)^2*(1-par_trial[3]^2)
    if (par_trial[6] < 0) || (par_trial[6] > 1) || (par_trial[7] < 0) || (par_trial[7] > 1) || (par_trial[4] < 0) || (par_trial[5] < 0) || (abs(par_trial[3]) > 1) || (par_trial[3]^2 > aonr2) then
      log_like_trial = -1e100
    else
      log_like_trial = compute_likelihood_save(time_phot,par_trial,flux,aex,al_small,indx,logdeta,p)
    end
# Next, determine whether to accept this trial step:
    alp = z^(nvary-1)*exp(log_like_trial - log_like_mcmc[j,i-1])
#    if rand() < 0.0001
#      println("Step: ",i," Walker: ",j," Chi-square: ",log_like_trial," Prob: ",alp," Frac: ",accept/(mod(i-1,1000)*nwalkers+j))
#    end
    if alp >= rand()
# If step is accepted, add it to the chains!
      par_mcmc[j,i,:] = par_trial
      log_like_mcmc[j,i,:] = log_like_trial
#      logprior_mcmc[j,i,:] = log_prior + log_det
      accept = accept + 1
    else
# If step is rejected, then copy last step:
      par_mcmc[j,i,:] = par_mcmc[j,i-1,:]
      log_like_mcmc[j,i,:] = log_like_mcmc[j,i-1]
#      logprior_mcmc[j,i,:] = logprior_mcmc[j,i-1,:]
    end
  end
  if mod(i,100) == 0
    eltime = toq()
    println("Number of steps: ",i," acceptance rate: ",accept/(100*nwalkers)," elapsed time [hr]: ",eltime/3600.)
    accept = 0
    tic()
  end
end

for i=1:nvary
  for j=1:nwalkers
    plot(vec(par_mcmc[j,1:nsteps,ivary[i]]))
  end
  xlabel("MCMC step")
  ylabel(pname[ivary[i]])
  println("Hit return to continue")
  read(STDIN,Char)
  clf()
end

@save "tyc3559_mcmc_norv_fixedkernel_v02.jld"


text:
import Pkg; Pkg.add("JuMP")
import Pkg; Pkg.add("Ipopt")

max_stars_repo_name: joaopedromattos/Robot-dance

max_issues_repo_name: joaopedromattos/Robot-dance

max_forks_repo_name: joaopedromattos/Robot-dance


text:
module TestPromotionConversion

using Test
using DiffCheckNumbers

constant = DiffCheck{Float64, ()}(1.0, NamedTuple{(), NTuple{0, Float64}}(( )))
x = DiffCheck{Float64, (:x, )}(1.0, NamedTuple{(:x, ), NTuple{1, Float64}}((1.0, )))
xy = DiffCheck{Float64, (:y, :x)}(1.0, NamedTuple{(:y, :x), NTuple{2, Float64}}((1.0, 1.0)))

z = DiffCheck(Complex{Float64}(1.0), :z)
zw = DiffCheck(1.0+im*0, [:z, :w])

#########
#Conversions
#########

# Testing Base.convert(::Type{DiffCheck{T, S, N}}, z::DiffCheck{T, S, N})
@test convert(DiffCheck{Float64, (), 0}, constant) === constant
@test convert(DiffCheck{Float64, (:x, ), 1}, x) === x
@test convert(DiffCheck{Float64, (:y, :x), 2}, xy) === xy
@test convert(DiffCheck{Complex{Float64}, (:z, ), 1}, z) === z
@test convert(DiffCheck{Complex{Float64}, (:w, :z), 2}, zw) === zw

# Testing Base.convert(::Type{DiffCheck{T, S, N}}, z::DiffCheck{R, S, N})
@test convert(DiffCheck{Complex{Float64}, (), 0}, constant) === DiffCheck{Complex{Float64}, ()}(value(constant), StaticTuple{Complex{Float64}, (), 0}(()))
@test convert(DiffCheck{Complex{Float64}, (:x, ), 1}, x) === DiffCheck{Complex{Float64}, (:x, )}(value(constant), StaticTuple{Complex{Float64}, (:x, ), 1}((1, )))
@test convert(DiffCheck{Complex{Float64}, (:y, :x), 2}, xy) === DiffCheck{Complex{Float64}, (:y, :x)}(value(constant), StaticTuple{Complex{Float64}, (:y, :x), 2}((1, 1)))

# convert a complex with no Imaginary part to Real
@test convert(DiffCheck{Float64, (:z, ), 1}, z) === DiffCheck{Float64, (:z, )}(1.0, StaticTuple{Complex{Float64}, (:z, ), 1}((1, )))
# throws the right exception when it is not possible
@test try
		convert(DiffCheck{Float64, (:z, ), 1}, DiffCheck(1.0+im, :z))
	  catch e
		isa(e, InexactError)
	  end

# Test convert(::Type{DiffCheck{T, (), 0}}, x::T)
@test convert(DiffCheck{Float64, (), 0}, 1.0) === constant

# Test convert(::Type{DiffCheck{T, S, N}}, x::T)
@test convert(DiffCheck{Float64, (:y, :x), 2}, 1.0) === DiffCheck{Float64, (:y, :x)}(1.0, NamedTuple{(:y, :x), NTuple{2, Float64}}((0.0, 0.0)))

# Test convert(::Type{T}, z::DiffCheck), converts only constants
@test convert(Float64, constant) === 1.0
@test try
		convert(Float64, x)
	  catch e
		isa(e, InexactError)
	  end

# Test convert(::Type{DiffCheck{T₁, S₁, N₁}}, x::DiffCheck{T₂, S₂, N₂}), only works upwards
@test convert(DiffCheck{Float64, (:y, :x), 2}, x) === DiffCheck{Float64, (:y, :x)}(1.0, NamedTuple{(:y, :x), NTuple{2, Float64}}((0.0, 1.0)))

@test try
		convert(DiffCheck{Float64, (:x, ), 1}, xy)
	  catch e
		isa(e, InexactError)
	  end

########
#Promotions
########

@test promote_type(DiffCheck{Float64, (:y, :x), 2}, DiffCheck{Complex{Float64}, (:y, :x), 2}) === DiffCheck{Complex{Float64}, (:y, :x), 2}
@test promote_type(DiffCheck{Float64, (:y, :x), 2}, Complex{Float64}) === DiffCheck{Complex{Float64}, (:y, :x), 2}
@test promote_type(DiffCheck{Float64, (:y, :x), 2}, Float64) === DiffCheck{Float64, (:y, :x), 2}
@test promote_type(DiffCheck{Float64, (), 0}, Float64) === DiffCheck{Float64, (), 0}
@test widen(DiffCheck{Float64, (:y, :x), 2}) === DiffCheck{BigFloat, (:y, :x), 2}

end


text:
using ModelingToolkit, OrdinaryDiffEq, Test, NonlinearSolve
using ModelingToolkit: topsort_equations

@variables t x(t) y(t) z(t) k(t)
eqs = [
       x ~ y + z
       z ~ 2
       y ~ 2z + k
      ]

sorted_eq = topsort_equations(eqs, [x, y, z, k])

ref_eq = [
          z ~ 2
          y ~ 2z + k
          x ~ y + z
         ]
@test ref_eq == sorted_eq

@test_throws ArgumentError topsort_equations([
                                             x ~ y + z
                                             z ~ 2
                                             y ~ 2z + x
                                            ], [x, y, z, k])

@parameters t σ ρ β
@variables x(t) y(t) z(t) a(t) u(t) F(t)
D = Differential(t)

test_equal(a, b) = @test isequal(simplify(a, polynorm=true), simplify(b, polynorm=true))

eqs = [
       D(x) ~ σ*(y-x)
       D(y) ~ x*(ρ-z)-y + β
       0 ~ z - x + y
       0 ~ a + z
       u ~ z + a
      ]

lorenz1 = ODESystem(eqs,t,name=:lorenz1)

lorenz1_aliased = structural_simplify(lorenz1)
io = IOBuffer(); show(io, MIME("text/plain"), lorenz1_aliased); str = String(take!(io))
@test all(s->occursin(s, str), ["lorenz1", "States (2)", "Parameters (3)"])
reduced_eqs = [
               D(x) ~ σ*(y - x)
               D(y) ~ β + x*(ρ - (x - y)) - y
              ]
test_equal.(equations(lorenz1_aliased), reduced_eqs)
@test isempty(setdiff(states(lorenz1_aliased), [x, y, z]))
test_equal.(observed(lorenz1_aliased), [
                                        u ~ 0
                                        z ~ x - y
                                        a ~ -z
                                       ])

# Multi-System Reduction

@variables s(t)
eqs1 = [
        D(x) ~ σ*(y-x) + F,
        D(y) ~ x*(ρ-z)-u,
        D(z) ~ x*y - β*z,
        u ~ x + y - z,
       ]

lorenz = name -> ODESystem(eqs1,t,name=name)
lorenz1 = lorenz(:lorenz1)
ss = ModelingToolkit.get_structure(initialize_system_structure(lorenz1))
@test isequal(ss.fullvars, [x, y, z, D(x), D(y), D(z), F, u])
lorenz2 = lorenz(:lorenz2)

connected = ODESystem([s ~ a + lorenz1.x
                       lorenz2.y ~ s
                       lorenz1.F ~ lorenz2.u
                       lorenz2.F ~ lorenz1.u],t,systems=[lorenz1,lorenz2])
@test length(Base.propertynames(connected)) == 10
@test isequal((@nonamespace connected.lorenz1.x), x)

# Reduced Flattened System

reduced_system = structural_simplify(connected)

@test setdiff(states(reduced_system), [
        s
        a
        lorenz1.x
        lorenz1.y
        lorenz1.z
        lorenz1.u
        lorenz2.x
        lorenz2.y
        lorenz2.z
        lorenz2.u
       ]) |> isempty

@test setdiff(parameters(reduced_system), [
        lorenz1.σ
        lorenz1.ρ
        lorenz1.β
        lorenz2.σ
        lorenz2.ρ
        lorenz2.β
       ]) |> isempty

reduced_eqs = [
               D(lorenz1.x) ~ lorenz1.σ*((lorenz1.y) - (lorenz1.x)) - ((lorenz2.z) - (lorenz2.x) - (lorenz2.y))
               D(lorenz1.y) ~ lorenz1.z + lorenz1.x*(lorenz1.ρ - (lorenz1.z)) - (lorenz1.x) - (lorenz1.y)
               D(lorenz1.z) ~ lorenz1.x*lorenz1.y - (lorenz1.β*(lorenz1.z))
               D(lorenz2.x) ~ lorenz2.σ*((lorenz2.y) - (lorenz2.x)) - ((lorenz1.z) - (lorenz1.x) - (lorenz1.y))
               D(lorenz2.y) ~ lorenz2.z + lorenz2.x*(lorenz2.ρ - (lorenz2.z)) - (lorenz2.x) - (lorenz2.y)
               D(lorenz2.z) ~ lorenz2.x*lorenz2.y - (lorenz2.β*(lorenz2.z))
              ]

test_equal.(equations(reduced_system), reduced_eqs)

observed_eqs = [
                s ~ lorenz2.y
                lorenz2.u ~ -((lorenz2.z) - (lorenz2.x) - (lorenz2.y))
                lorenz1.u ~ -((lorenz1.z) - (lorenz1.x) - (lorenz1.y))
                a ~ s - (lorenz1.x)
                lorenz1.F ~ lorenz2.u
                lorenz2.F ~ lorenz1.u
               ]
test_equal.(observed(reduced_system), observed_eqs)

pp = [
      lorenz1.σ => 10
      lorenz1.ρ => 28
      lorenz1.β => 8/3
      lorenz2.σ => 10
      lorenz2.ρ => 28
      lorenz2.β => 8/3
     ]
u0 = [
      lorenz1.x => 1.0
      lorenz1.y => 0.0
      lorenz1.z => 0.0
      lorenz2.x => 1.0
      lorenz2.y => 0.0
      lorenz2.z => 0.0
     ]
prob1 = ODEProblem(reduced_system, u0, (0.0, 100.0), pp)
solve(prob1, Rodas5())

# issue #724 and #716
let
    @parameters t
    D = Differential(t)
    @variables x(t) u(t) y(t)
    @parameters a b c d
    ol = ODESystem([D(x) ~ a * x + b * u; y ~ c * x + d * u], t, name=:ol)
    @variables u_c(t) y_c(t)
    @parameters k_P
    pc = ODESystem(Equation[u_c ~ k_P * y_c], t, name=:pc)
    connections = [
        ol.u ~ pc.u_c
        pc.y_c ~ ol.y
    ]
    connected = ODESystem(connections, t, systems=[ol, pc])
    @test equations(connected) isa Vector{Equation}
    reduced_sys = structural_simplify(connected)
    ref_eqs = [
               D(ol.x) ~ ol.a*ol.x + ol.b*pc.u_c
               0 ~ -pc.u_c - (pc.k_P*((-ol.c*(ol.x)) - (ol.d*(pc.u_c))))
              ]
    @test ref_eqs == equations(reduced_sys)
end

# issue #889
let
    @parameters t
    D = Differential(t)
    @variables x(t)
    @named sys = ODESystem([0 ~ D(x) + x], t, [x], [])
    sys = structural_simplify(sys)
    @test_throws ModelingToolkit.InvalidSystemException ODEProblem(sys, [1.0], (0, 10.0))
end

# NonlinearSystem
@parameters t
@variables u1(t) u2(t) u3(t)
@parameters p
eqs = [
       u1 ~ u2
       u3 ~ u1 + u2 + p
       u3 ~ hypot(u1, u2) * p
      ]
sys = NonlinearSystem(eqs, [u1, u2, u3], [p])
reducedsys = structural_simplify(sys)
@test observed(reducedsys) == [u2 ~ 0.5(u3 - p); u1 ~ u2]

u0 = [
      u1 => 1
      u2 => 1
      u3 => 0.3
     ]
pp = [2]
nlprob = NonlinearProblem(reducedsys, u0, pp)
reducedsol = solve(nlprob, NewtonRaphson())
residual = fill(100.0, length(states(reducedsys)))
nlprob.f(residual, reducedsol.u, pp)
@test all(x->abs(x) < 1e-5, residual)

N = 5
@variables xs[1:N]
A = reshape(1:N^2, N, N)
eqs = xs .~ A * xs
sys′ = NonlinearSystem(eqs, xs, [])
sys = structural_simplify(sys′)


text:
using CryptoTools
using Test
using Random

@testset "OneTimePad" begin
    plaintext = "I am plaintext"
    key = rand(RandomDevice(), UInt8, length(plaintext))
    cipher = OneTimePad(key)
    @test decrypt(cipher, encrypt(cipher, plaintext)) == plaintext

    plaintext = [0x00, 0x01, 0x02, 0x03, 0x04, 0x05]
    key = rand(RandomDevice(), UInt8, length(plaintext))
    cipher = OneTimePad(key)
    @test decrypt(cipher, encrypt(cipher, plaintext)) == plaintext

end

@testset "AES" begin
    # manual test
    plaintext = hex2bytes("00112233445566778899aabbccddeeff")
    key = hex2bytes("000102030405060708090a0b0c0d0e0f")
    cipher = AES(key)
    @test encrypt(cipher, plaintext, false) == hex2bytes("69c4e0d86a7b0430d8cdb78070b4c55a")

    # manual test
    plaintext = hex2bytes("00000000000000000000000000000000")
    key = hex2bytes("2c6202f9a582668aa96d511862d8a279")
    cipher = AES(key)
    @test encrypt(cipher, plaintext, false) == hex2bytes("12b620bb5eddcde9a07523e59292a6d7")

    # random test
    plaintext = rand(RandomDevice(), UInt8, 16)
    key = rand(RandomDevice(), UInt8, 16)
    cipher = AES(key)
    @test decrypt(cipher, encrypt(cipher, plaintext, false), false) == plaintext

    # check encryption capability of a plaintext that is shorter than 16 bytes
    @test decrypt(cipher, encrypt(cipher, [0x00])) == [0x00]

    plaintext = repeat(hex2bytes("00000000000000000000000000000000"), 2)
    @test decrypt(cipher, encrypt(cipher, plaintext)) == plaintext
end


text:
using Polyhedra

include("simplex.jl")
include("permutahedron.jl")
include("board.jl")

myeq(x::Real, y::Real) = myeq(promote(x, y)...)
myeq{T<:Real}(x::T, y::T) = x == y
myeq{T<:AbstractFloat}(x::T, y::T) = y < x+1024*eps(T) && x < y+1024*eps(T)
myeq{S<:Real,T<:Real}(x::Vector{S}, y::Vector{T}) = myeq(promote(x, y)...)
myeq{T<:Real}(x::Vector{T}, y::Vector{T}) = x == y
myeq{T<:AbstractFloat}(x::Vector{T}, y::Vector{T}) = myeq(norm(x - y), zero(T))
myeqzero{T<:Real}(x::T) = myeq(x, zero(T))

tomatrix(M::Matrix) = M
function tomatrix(v::Vector)
  M = Matrix{eltype(v)}(length(v), 1)
  M[:,1] = v
  M
end

function inlinspace(x, L)
  for i in 1:size(L, 1)
    y = vec(L[i,:])
    # remove component
    x = x * dot(y, y) - y * dot(y, x)
  end
  myeqzero(norm(x))
end

function inequality_fulltest(p::Polyhedron, A, b, linset)
  A = tomatrix(A)
  detecthlinearities!(p)
  removeredundantinequalities!(p)
  ine = SimpleHRepresentation(getinequalities(p))
  @test size(ine.A) == size(A)
  @test length(ine.linset) == length(linset)

  aff = SimpleHRepresentation(getinequalities(affinehull(p)))
  affAb = [aff.b aff.A]
  inaff(x) = inlinspace(x, affAb)

  for i in 1:size(A, 1)
    found = false
    for j in 1:size(ine.A, 1)
      # vec for julia 0.4
      if !((i in linset) $ (j in ine.linset)) && inaff([b[i]-ine.b[j];vec(A[i,:]-ine.A[j,:])])
        found = true
        break
      end
    end
    @test found
  end
end
function generator_fulltest(p::Polyhedron, V, R=Matrix{eltype(V)}(0, size(V, 2)), Vlinset = IntSet(), Rlinset = IntSet())
  V = tomatrix(V)
  R = tomatrix(R)
  detectvlinearities!(p)
  removeredundantgenerators!(p)
  ext = SimpleVRepresentation(getgenerators(p))
  @test size(ext.V) == size(V)
  @test size(ext.R) == size(R)
  @test length(ext.Vlinset) == length(Vlinset)
  @test length(ext.Rlinset) == length(Rlinset)
  for i in 1:size(V, 1)
    found = false
    for j in 1:size(ext.V, 1)
      if myeq(vec(V[i, :]), vec(ext.V[j, :]))
        found = true
        break
      end
    end
    @test found
  end
  linspace = ext.R[collect(ext.Rlinset),:]
  inlin(x) = inlinspace(vec(x), linspace)
  for i in 1:size(R, 1)
    found = false
    for j in 1:size(ext.R, 1)
      if !((i in Rlinset) $ (j in ext.Rlinset)) && inlin(R[i,:]-ext.R[j,:])
      #if parallel(vec(R[i, :]), vec(ext.R[j, :]), (i in Rlinset) || (j in ext.Rlinset))
        found = true
        break
      end
    end
    @test found
  end
end
#generator_fulltest(p::Polyhedron, V) = generator_fulltest(p, V, Matrix{eltype(V)}(0, size(V, 2)))

function alltests{Lib<:PolyhedraLibrary}(lib::Lib)
  simplextest(lib)
  permutahedrontest(lib)
  boardtest(lib)
end


text:
#import CuArrays: _cuview, ViewIndex, NonContiguous
#using GPUArrays: genperm
# fallback to SubArray when the view is not contiguous

#=
function LinearAlgebra.permutedims!(dest::GPUArray, src::GPUArray, perm) where N
    perm isa Tuple || (perm = Tuple(perm))
    gpu_call(dest, (dest, src, perm)) do state, dest, src, perm
        I = @cartesianidx src state
        @inbounds dest[genperm(I, perm)...] = src[I...]
        return
    end
    return dest
end
=#

import CUDAnative: pow, abs, angle
for (RT, CT) in [(:Float64, :ComplexF64), (:Float32, :ComplexF32)]
    @eval cp2c(d::$RT, a::$RT) = CUDAnative.ComplexF64(d*CUDAnative.cos(a), d*CUDAnative.sin(a))
    for NT in [RT, :Int32]
        @eval CUDAnative.pow(z::$CT, n::$NT) = CUDAnative.ComplexF64((CUDAnative.pow(CUDAnative.abs(z), n)*CUDAnative.cos(n*CUDAnative.angle(z))), (CUDAnative.pow(CUDAnative.abs(z), n)*CUDAnative.sin(n*CUDAnative.angle(z))))
    end
end

@inline function bit_count(x::UInt32)
    x = ((x >> 1) & 0b01010101010101010101010101010101) + (x & 0b01010101010101010101010101010101)
    x = ((x >> 2) & 0b00110011001100110011001100110011) + (x & 0b00110011001100110011001100110011)
    x = ((x >> 4) & 0b00001111000011110000111100001111) + (x & 0b00001111000011110000111100001111)
    x = ((x >> 8) & 0b00000000111111110000000011111111) + (x & 0b00000000111111110000000011111111)
    x = ((x >> 16)& 0b00000000000000001111111111111111) + (x & 0b00000000000000001111111111111111)
    return x
end

@inline function bit_count(x::Int32)
    x = ((x >> 1) & Int32(0b01010101010101010101010101010101)) + (x & Int32(0b01010101010101010101010101010101))
    x = ((x >> 2) & Int32(0b00110011001100110011001100110011)) + (x & Int32(0b00110011001100110011001100110011))
    x = ((x >> 4) & Int32(0b00001111000011110000111100001111)) + (x & Int32(0b00001111000011110000111100001111))
    x = ((x >> 8) & Int32(0b00000000111111110000000011111111)) + (x & Int32(0b00000000111111110000000011111111))
    x = ((x >> 16)& Int32(0b00000000000000001111111111111111)) + (x & Int32(0b00000000000000001111111111111111))
    return x
end

bit_count(UInt32(0b11111))

# TODO
# support norm(view(reshape(A, m, n), :, 1))
using LinearAlgebra
import LinearAlgebra: norm
const CuSubArr{T, N} = Union{CuArray{T, N}, SubArray{T, N, <:CuArray}}
norm2(A::CuSubArr; dims=1) = mapreduce(abs2, +, A, dims=dims) .|> CUDAnative.sqrt

export piecewise, cudiv
@inline function cudiv(x::Int)
    max_threads = 256
    threads_x = min(max_threads, x)
    threads_x, ceil(Int, x/threads_x)
end

# NOTE: the maximum block size is 65535
@inline function cudiv(x::Int, y::Int)
    max_threads = 256
    threads_x = min(max_threads, x)
    threads_y = min(max_threads ÷ threads_x, y)
    threads = (threads_x, threads_y)
    blocks = ceil.(Int, (x, y) ./ threads)
    threads, blocks
end

fix_cudiv(A::AbstractVector, firstdim::Int) = cudiv(firstdim)
fix_cudiv(A::AbstractMatrix, firstdim::Int) = cudiv(firstdim, size(A,2))

piecewise(state::AbstractVector, inds) = state
piecewise(state::AbstractMatrix, inds) = @inbounds view(state,:,inds[2])

import Base: kron, getindex
function kron(A::Union{CuArray{T1}, Adjoint{<:Any, <:CuArray{T1}}}, B::Union{CuArray{T2}, Adjoint{<:Any, <:CuArray{T2}}}) where {T1, T2}
    res = CuArrays.zeros(promote_type(T1,T2), (size(A).*size(B))...)
    @inline function kernel(res, A, B)
        state = (blockIdx().x-1) * blockDim().x + threadIdx().x
        inds = GPUArrays.gpu_ind2sub(res, state)
        inds_A = (inds.-1) .÷ size(B) .+ 1
        inds_B = (inds.-1) .% size(B) .+ 1
        state <= length(res) && (@inbounds res[state] = A[inds_A...]*B[inds_B...])
        return
    end

    X, Y = cudiv(length(res))
    @cuda threads=X blocks=Y kernel(res, A, B)
    res
end

function getindex(A::CuVector{T}, B::CuArray{<:Integer}) where T
    res = CuArrays.zeros(T, size(B)...)
    @inline function kernel(res, A, B)
        state = (blockIdx().x-1) * blockDim().x + threadIdx().x
        state <= length(res) && (@inbounds res[state] = A[B[state]])
        return
    end

    X, Y = cudiv(length(B))
    @cuda threads=X blocks=Y kernel(res, A, B)
    res
end


text:
using MultipleViewGeometry, Base.Test
using MultipleViewGeometry.ModuleCostFunction
using MultipleViewGeometry.ModuleTypes
using MultipleViewGeometry.ModuleConstraints
using MultipleViewGeometry.ModuleConstruct
using MultipleViewGeometry.ModuleDraw
using BenchmarkTools, Compat
using StaticArrays
using  MAT, Plots

# Load MATLAB matrices that represent a pair of images and that contain
# a set of manually matched corresponding points.
file = matopen("./data/teapot.mat")
X = read(file,"pts3D")
X = X[:,1:10:end]
close(file)

# Fix random seed.
srand(1234)
plotlyjs()

#𝒳 = [Point3DH(x,y,z,1.0) for x=-1:0.5:10 for y=-1:0.5:10 for z=2:-0.1:1]

𝒳 = [Point3DH(X[1,i],X[2,i],X[3,i],1.0) for i = 1:size(X,2)]

#X = reinterpret(Float64,map(SVector{4,Float64},𝒳),(4,length(𝒳)))
#Z = reinterpret(SVector{4,Float64},(3000,4))
# reinterpret(SVector{4,Float64}, X, (size(X,2),))#
# reinterpret(SVector{4,Float64}, Z, (size(Z,2),))#
# reinterpret(SVector{4,Float64}, X, (3000,1))

# Intrinsic and extrinsic parameters of camera one.
𝐊₁ = eye(3)
𝐑₁ = eye(3)
𝐭₁ = [-250.0, 0.0, 2500.0]
𝐜₁ = [0.0, 0.0, 0.0]
𝐄₁ = [𝐑₁ 𝐭₁]

# Intrinsic and extrinsic parameters of camera two.
𝐊₂ = eye(3)
𝐑₂ =  eye(3)
𝐭₂ = [250.0,   0.0, 2500.0]
𝐜₂ = [0.0, 0.0, 0.0]
𝐄₂ = [𝐑₂ 𝐭₂]

# Camera projection matrices.
𝐏₁ = construct(ProjectionMatrix(),𝐊₁,𝐑₁,𝐭₁)
𝐏₂ = construct(ProjectionMatrix(),𝐊₂,𝐑₂,𝐭₂)

# Set of corresponding points.
ℳ = project(Pinhole(),𝐏₁,𝒳)
ℳʹ = project(Pinhole(),𝐏₂,𝒳)

# Estimate of the fundamental matrix and the true fundamental matrix.
𝐅 = estimate(FundamentalMatrix(), DirectLinearTransform(), (ℳ, ℳʹ))

# Convert the array of Point3DH into a 4 x N matrix to simplify
# the plotting of the data points.
#X = reinterpret(Float64,map(SVector{4,Float64},𝒳),(4,length(𝒳)))

# Visualise the data points
p1 = Plots.plot(X[1,:],X[2,:],X[3,:],seriestype = :scatter, ms=1,grid = false, box = :none, legend = false)
draw!(WorldCoordinateSystem3D(), 450, p1)
draw!(Camera3D(), 𝐊₁, 𝐑₁, 𝐭₁, 250, p1)
draw!(Camera3D(), 𝐊₂, 𝐑₂, 𝐭₂, 250, p1)

# Plot the projections of the point cloud in the image pair.
p2 =  Plots.plot();
for n = 1:length(ℳ)
    m = ℳ[n]
    Plots.plot!([m[1]],[m[2]], grid = false, box = :none, legend = false,
                    seriestype = :scatter, ms = 2, markercolor=:Black,
                    aspect_ratio = :equal)
end

p3 =  Plots.plot();
for n = 1:length(ℳʹ)
    mʹ = ℳʹ[n]
    Plots.plot!([mʹ[1]],[mʹ[2]], grid = false, box = :none, legend = false,
                    seriestype = :scatter, ms = 2, markercolor=:Black,
                    aspect_ratio = :equal)
end

# Visualise the 3D point cloud, as well as the projected images.
l = @layout [ a; [b c] ]
p4 = Plots.plot(p1,p2, p3, layout = l)

𝒴 = triangulate(DirectLinearTransform(),𝐏₁,𝐏₂,(ℳ,ℳʹ))

Y = reinterpret(Float64,map(SVector{4,Float64},𝒴),(4,length(𝒴)))
# Visualise the data points
p5 = Plots.plot(Y[1,:],Y[2,:],Y[3,:],seriestype = :scatter, ms=1,grid = false, box = :none, legend = false)
draw!(WorldCoordinateSystem3D(), 450, p1)
draw!(Camera3D(), 𝐊₁, 𝐑₁, 𝐭₁, 250, p1)
draw!(Camera3D(), 𝐊₂, 𝐑₂, 𝐭₂, 250, p1)

# 𝐅 = construct(FundamentalMatrix(),𝐊₁,𝐑₁,𝐭₁,𝐊₂,𝐑₂,𝐭₂)
#
# # To triangulate the corresponding points using the Fundamental matrix, we first
# # have to factorise the Fundamental matrix into a pair of Camera matrices. Due
# # to projective ambiguity, the camera matrices are not unique, and so the
# # triangulated 3D points will most probably not match the original 3D points.
# # However, when working with noiseless data, the projections of the triangulated
# # points should satisfy the epipolar constraint. We can use this fact to
# # validate that the triangulation is correctly implemented.
# 𝒴 = triangulate(DirectLinearTransform(),𝐅,(ℳ,ℳʹ))
# Y = reinterpret(Float64,map(SVector{4,Float64},𝒴),(4,length(𝒴)))
# p6 = Plots.plot(Y[1,:],Y[2,:],Y[3,:],seriestype = :scatter, ms=1,grid = false, box = :none, legend = false)
#
# 𝐐₁, 𝐐₂ = construct(ProjectionMatrix(),𝐅)


text:
# define a special case of exponential Kernel, because Claudia said so

####
# Population and connection constructors

# the taus are in the trace, now. No unittype or kernel object
# this is pretty much a wrapper around neural traces
struct PopulationStateExpKernel{N}  <: AbstractPopulationState
  label::Symbol
  n::Int64
  traces::NTuple{N,Trace}
end
nneurons(ps::PopulationStateExpKernel) = ps.n

function PopulationStateExpKernel(n::Int64,traces...;
    label::Union{String,Nothing}=nothing)
  label = something(label,rand_label())
  PopulationStateExpKernel(label,n,traces)
end
function reset!(ps::PopulationStateExpKernel)
  for tra in ps.traces
    reset!(tra)
  end
  return nothing
end

struct PopulationExpKernel{N,PS<:PopulationStateExpKernel,
    TC<:NTuple{N,Connection},
    TP<:NTuple{N,PopulationStateExpKernel},
    NL<:AbstractNonlinearity} <: AbstractPopulation
  state::PS
  connections::TC
  pre_states::TP
  input::Vector{Float64} # input only used for external currents
  nonlinearity::NL # nonlinearity if present
  spike_proposals::Vector{Float64} # memory allocation
end
function PopulationExpKernel(state::PopulationStateExpKernel,input::Vector{Float64},
    (conn_pre::Tuple{C,PS} where {C<:Connection,PS<:PopulationStateExpKernel})...;
      nonlinearity::AbstractNonlinearity=NLRelu())
  connections = Tuple(getindex.(conn_pre,1))
  typeassert.(connections,ConnectionExpKernel)
  pre_states = Tuple(getindex.(conn_pre,2))
  spike_proposals = fill(Inf,nneurons(state))
  return PopulationExpKernel(state,connections,pre_states,input,nonlinearity,spike_proposals)
end
# one population only!
function PopulationExpKernel(state::PopulationStateExpKernel,conn::Connection,input::Vector{Float64};
    nonlinearity::AbstractNonlinearity=NLRelu())
  return PopulationExpKernel(state,input,(conn,state);nonlinearity=nonlinearity)
end

struct ConnectionExpKernel{N,PL<:NTuple{N,PlasticityRule}} <: Connection
  weights::Matrix{Float64}
  pre_trace::Trace
  trace_proposal::Vector{Float64}
  plasticities::PL
end
function ConnectionExpKernel(weights::Matrix{Float64},pre_trace::Trace)
  plast = NTuple{0,NoPlasticity}()
  _,npre = size(weights)
  @assert npre == nneurons(pre_trace)
  tra_prop = Vector{Float64}(undef,npre)
  return ConnectionExpKernel(weights,pre_trace,tra_prop,plast)
end
function ConnectionExpKernel(weights::Matrix{Float64},pre_trace::Trace,
     (plasticity_rules::PL where PL<:PlasticityRule)...)
  _,npre = size(weights)
  @assert npre == nneurons(pre_trace)
  tra_prop = Vector{Float64}(undef,npre)
  return ConnectionExpKernel(weights,pre_trace,tra_prop,plasticity_rules)
end

@inline function trace_proposal!(t_now::Real,conn::ConnectionExpKernel)
  return trace_proposal!(conn.trace_proposal,t_now,conn.pre_trace)
end

abstract type Recorder end

struct RecurrentNetworkExpKernel{N,TP<:NTuple{N,AbstractPopulation},TR<:NTuple{N,Recorder}}
  populations::TP
  recorders::TR
end
# one population constructor
function RecurrentNetworkExpKernel(pop::AbstractPopulation,recorders...)
  if isempty(recorders)
    #recorders=NTuple{0,RecNothing}()
    recorders=(RecNothing(),)
    #@show typeof(recorders)
  end
  return RecurrentNetworkExpKernel((pop,),recorders)
end

# reset
function reset!(rn::RecurrentNetworkExpKernel)
  # reset internal variables
  for pop in rn.populations
    reset!(pop.state)
    fill!(pop.spike_proposals,Inf)
    # reset plasticity operators if present
    for conn in pop.connections
      for plast in conn.plasticities
        reset!(plast)
      end
    end
  end
  # reset recorders
  for rec in rn.recorders
    reset!(rec)
  end
  return nothing
end
@inline function npopulations(ntw::RecurrentNetworkExpKernel)
  return length(ntw.populations)
end

function burn_spike!(t_spike::Real,ps::PopulationStateExpKernel,idx_update::Integer)
  # take care of traces
  for tra in ps.traces
    propagate_for_dynamics!(t_spike,tra) # update full trace to t_spike
    update_for_dynamics!(tra,idx_update) # add pulse at index
  end
  return nothing
end

function compute_rate(t_now::Real,external_input::Real,
    pop::PopulationExpKernel, idxneu::Integer)
  ret = external_input
  for conn in pop.connections
    tra_tnow = trace_proposal!(t_now,conn)
    wij_all = view(conn.weights,idxneu,:)
    ret += dot(wij_all,tra_tnow)
  end
  return apply_nonlinearity(ret,pop.nonlinearity)
end

# Thinning algorith, e.g.  Laub,Taimre,Pollet 2015
function compute_next_spike(t_now::Real,pop::PopulationExpKernel,ineu::Integer;Tmax::Real=100.0)
  t_start = t_now
  t = t_now
  ext_inp = pop.input[ineu]
  freq(t) = compute_rate(t,ext_inp,pop,ineu)
  while (t-t_start)<Tmax
    (M::Float64) = freq(t)
    Δt =  -log(rand())/M # rand(Exponential())/M
    t = t+Δt
    u = rand()*M # random between 0 and M
    (u_up::Float64) = freq(t)
    if u <= u_up
      return t
    end
  end
  return Tmax + t_start
end

###
# recorders methods
struct RecNothing <: Recorder end

function reset!(::RecNothing)
  return nothing
end

# general signature: record_stuff!(rec,tfire,popfire,neufire,label_fire,ntw)
# @inline function record_stuff!(::RecNothing,::Real,
#     ::Integer,::Integer,::Symbol,
#     ::RecurrentNetworkExpKernel)
#   return nothing
# end
@inline function record_stuff!(::RecNothing,::Real,
    whatevs...)
  return nothing
end

struct RecFullTrain{N} <: Recorder
  nrec::Int64
  timesneurons::NTuple{N,Tuple{Vector{Float64},Vector{Int64}}}
  k_rec::Vector{Int64}
end
function RecFullTrain(nrec::Integer,(npops::Integer)=1)
  timesneurons=ntuple(_-> (fill(NaN,nrec),fill(-1,nrec)),npops)
  k_rec = fill(0,npops)
  RecFullTrain(nrec,timesneurons,k_rec)
end
function reset!(rec::RecFullTrain)
  fill!(rec.k_rec,0)
  for (spktim,spkneu) in rec.timesneurons
    fill!(spktim,NaN)
    fill!(spkneu,-1)
  end
  return nothing
end

# general signature: record_stuff!(rec,tfire,popfire,neufire,label_fire,ntw)

function record_stuff!(rec::RecFullTrain,tfire::Real,
    popfire::Integer,neufire::Integer,::Symbol,
    ::RecurrentNetworkExpKernel)
  k = rec.k_rec[popfire]+1
  spiketimes,spikeneurons = rec.timesneurons[popfire]
  if !checkbounds(Bool,spiketimes,k)
    @error "recorder full!"
    return nothing
  end
  @inbounds spiketimes[k] = tfire
  @inbounds spikeneurons[k] = neufire
  rec.k_rec[popfire] = k
  return nothing
end

function dynamics_step!(t_now::Real,ntw::RecurrentNetworkExpKernel)
  npops = npopulations(ntw)
  proposals_best = Vector{Float64}(undef,npops)
  neuron_best = Vector{Int64}(undef,npops)
  # for each postsynaptic network, compute spike proposals
  for (kn,pop) in enumerate(ntw.populations)
    # update proposed next spike for each postsynaptic neuron
    nneu = nneurons(pop)
    for ineu in 1:nneu
      pop.spike_proposals[ineu] = compute_next_spike(t_now,pop,ineu)
    end
    # best candidate and neuron that fired it for one input network
    (proposals_best[kn], neuron_best[kn]) = findmin(pop.spike_proposals)
  end
  # select next spike (best across all input_networks)
  (tfire,popfire) = findmin(proposals_best)
  neufire = neuron_best[popfire]
  psfire = ntw.populations[popfire].state
  label_fire = psfire.label
  # update stuff for that specific neuron/population state :
  burn_spike!(tfire,psfire,neufire)
  # apply plasticity rules ! Each rule in each connection in each population
  for pop in ntw.populations
    post = pop.state
    for (conn,pre) in zip(pop.connections,pop.pre_states)
      for plast in conn.plasticities
        plasticity_update!(tfire,label_fire,neufire,post,conn,pre,plast)
      end
    end
  end
  # now trigger recorders.
  # Recorder objects will take care of which population to target
  for rec in ntw.recorders
    record_stuff!(rec,tfire,popfire,neufire,label_fire,ntw)
  end
  # update t_now
  return tfire
end

function dynamics_step_singlepopulation!(t_now::Real,ntw::RecurrentNetworkExpKernel)
  # update proposed next spike for each postsynaptic neuron
  popfire = 1
  pop = ntw.populations[1]
  nneu = nneurons(pop)
  for ineu in 1:nneu
    pop.spike_proposals[ineu] = compute_next_spike(t_now,pop,ineu)
  end
  (tfire, neufire) = findmin(pop.spike_proposals)
  psfire = pop.state
  label_fire = psfire.label
  # update stuff for that specific neuron/population state :
  burn_spike!(tfire,psfire,neufire)
  # apply plasticity rules ! Each rule in each connection in each population
  post = psfire
  for (conn,pre) in zip(pop.connections,pop.pre_states)
    for plast in conn.plasticities
      plasticity_update!(tfire,label_fire,neufire,post,conn,pre,plast)
    end
  end
  # now trigger recorders.
  # Recorder objects will take care of which population to target
  for rec in ntw.recorders
    record_stuff!(rec,tfire,popfire,neufire,label_fire,ntw)
  end
  # update t_now
  return tfire
end


text:
module Altro

import TrajectoryOptimization
import RobotDynamics
using StaticArrays
using Parameters
using DocStringExtensions
using BenchmarkTools
using Interpolations
using UnsafeArrays
using SolverLogging
using Crayons

using SparseArrays
using LinearAlgebra
using Logging
using Statistics
using TimerOutputs
using FiniteDiff

const TO = TrajectoryOptimization
const RD = RobotDynamics

using TrajectoryOptimization:
    integration, num_constraints, get_trajectory

import TrajectoryOptimization: rollout!, get_constraints, get_model, get_objective
import RobotDynamics: discrete_jacobian!, discrete_dynamics, dynamics

using TrajectoryOptimization:
    Problem,
    ConstraintList,
    AbstractObjective, Objective, #QuadraticObjective,
    AbstractTrajectory,
    DynamicsExpansion, # TODO: Move to ALTRO
    # ALConstraintSet,
    DynamicsConstraint,
    Traj,
    states, controls,
    Equality, Inequality, SecondOrderCone

using RobotDynamics:
    AbstractModel,
    QuadratureRule, Implicit, Explicit,
    AbstractKnotPoint, KnotPoint, StaticKnotPoint,
    state, control

# types
export
    ALTROSolver,
    # iLQRSolver,
    # AugmentedLagrangianSolver,
    SolverStats,
    SolverOptions

export
    solve!,
    benchmark_solve!,
    iterations,
    set_options!,
    status

# modules
export
    Problems

include("utils.jl")
include("infeasible_model.jl")
include("solvers.jl")
include("solver_opts.jl")

include("ilqr/ilqr.jl")
include("ilqr/ilqr_solve.jl")
include("ilqr/backwardpass.jl")
include("ilqr/rollout.jl")
# include("augmented_lagrangian/conic_penalties.jl")
include("augmented_lagrangian/alconval.jl")
include("augmented_lagrangian/ALconset.jl")
include("augmented_lagrangian/alcosts.jl")
include("augmented_lagrangian/al_solver.jl")
include("augmented_lagrangian/al_objective.jl")
include("augmented_lagrangian/al_methods.jl")
include("direct/primals.jl")
include("direct/pn.jl")
include("direct/pn_methods.jl")
include("altro/altro_solver.jl")

include("direct/copy_blocks.jl")
include("direct/direct_constraints.jl")

include("problems.jl")
# include("deprecated.jl")

end # module


text:
using Requests
using DataFrames
################################################################################
# this function returns a DataFrame of articles (pmid, title, abstract)
function fetchBreastCancerArticles(cancerType="breast cancer", researchType="diagnosis",
                                    minDate=1990,maxDate=2030, retmax=1000)

    search_terms = "\"$cancerType\"[TIAB] and \"$researchType\"[TIAB]"

    # Retrieve number of results - from ESearch results
    # define base search query for eutils
    base_search_query = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi"

    # run esearch to get pmids
    search_result = readstring(post(base_search_query;
    data = Dict("db" => "pubmed", "term" => "$search_terms", "retmax" => retmax,
                 "mindate"=>minDate, "maxdate"=>maxDate)))
    #print(search_result)

    # retrieve the pmids into an array
    pmids = []
    for result_line in split(search_result, "\n")
      # get list of pmids
      pmid = match(r"<Id>(\d+)<\/Id>", result_line)
      if pmid != nothing
        push!(pmids,pmid[1])
      end
    end

    # concatenate pmids into a single comma separated string
    id_string = join(collect(pmids), ",")

    # define base fetch query for eutils
    base_fetch_query = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi"

    # retrieve metadata for pmid Set
    fetch_result = readstring(post(base_fetch_query; data = Dict("db" => "pubmed",
    "id" => id_string, "rettype" => "medline", "retmode" => "text")))

    # save metadata in a file
    output_file = open("output/$researchType.txt", "w")
    write(output_file,fetch_result)
    close(output_file)

    # field set
    pmid_array = [] # title set
    title_array = [] # title set
    abstract_array=[] # article abstract
    date_created_array=[] # Create Date (CRDT).
    # read through each line
    for fetch_article in split(fetch_result, "PMID- ")
        # clear all field
        pmid=""
        full_title=""
        full_abstract=""
        date_created=""

        # extract PMID
        pmid_str = match(r"([0-9]+)", fetch_article)
        if pmid_str != nothing
          pmid=pmid_str[1]
        end

        for fetch_line in split(fetch_article, "\n")
          # get Title
          title = match(r"TI  - ([\w\W \r\n]+)", fetch_line)
          if title != nothing
              full_title=title[1]
          end

          # get abstract
          abstract_str = match(r"AB  - ([\w\W \r\n]+)", fetch_line)
          if abstract_str != nothing
              full_abstract=abstract_str[1]
          end

          # date created
          date_created_str = match(r"CRDT- ([\w\W \r\n]+)", fetch_line)
          if date_created_str != nothing
              date_created=date_created_str[1]
          end
        end

        # push everything
        if pmid != ""
            push!(pmid_array, pmid)
            push!(title_array, full_title)
            push!(abstract_array, full_abstract)
            push!(date_created_array, date_created)
        end
    end

    # return dictionary of pmid and title
    df = DataFrame(pmid=pmid_array,
                    title=title_array,
                    date_created=date_created_array,
                    main_abstract=abstract_array
                    )

    return(df, fetch_result)

end

# usage
#df, output_text = fetchBreastCancerArticles("breast cancer", "prevention",2018,2018)

# ref
#https://www.nlm.nih.gov/bsd/mms/medlineelements.html#ab

max_stars_repo_name: jmulloka/Breast-Cancer-Text-Mining

max_issues_repo_name: jmulloka/Breast-Cancer-Text-Mining

max_forks_repo_name: jmulloka/Breast-Cancer-Text-Mining


text:
import Pkg
Pkg.activate(".")
push!(LOAD_PATH, "/home/guilhermebodin/Documents/Github/LPIP.jl/src")
using LPIP, Test, GLPK

using MathOptInterface, JuMP

const MOI = MathOptInterface

A = [
    2. 1
    1 2
]
b = [4.;4]
c = [-4.; -3]
d = 3.

lpip_pb = LPIPLinearProblem{Float64}(A, b, c, d, 2)
params = LPIP.Params(;rho = 0.01, verbose = false)
lp = LPIP.interior_points(lpip_pb, params)
rho = 0.01
verbose = false
using BenchmarkTools
@benchmark begin
    lpip_pb = LPIPLinearProblem{Float64}($A, $b, $c, $d, 2)
    params = LPIP.Params(;rho = $rho, verbose = $verbose)
    lp = LPIP.interior_points(lpip_pb, params)
end

model = Model(with_optimizer(LPIP.Optimizer, rho = 0.01, verbose = true))
# model = Model(with_optimizer(GLPK.Optimizer))
@variable(model, x[1:2])
@constraint(model, A * x - b in MOI.Zeros(2))
@objective(model, Min, c'x + d)
@time optimize!(model)
JuMP.value.(x)
JuMP.objective_value(model)
JuMP.termination_status(model)


text:
using HorizonSideRobots
include("BackPath.jl")
include("Coord.jl")

import HorizonSideRobots: move!

function print_chess(r::Robot, n::Int)
    back_path = BackPath(r)
    coord = Coord(0, 0)
    side = Ost
    nord_border = false
    side_border = false
    while !nord_border
        while !side_border
            if (coord.x ÷ n + coord.y ÷ n) % 2 == 0
                putmarker!(r)
            end
            if !isborder(r, side)
                move!(coord, side)
                move!(r, side)
            else
                side_border = true
            end
        end
        if !isborder(r, Nord)
            move!(coord, Nord)
            move!(r, Nord)
            side = inverse(side)
            side_border = false
        else
            nord_border = true
        end
    end
    to_sud_west_corner = BackPath(r)
    back!(r, back_path)
end

max_stars_repo_name: Mananhina/mirea-progs

max_issues_repo_name: Mananhina/mirea-progs

max_forks_repo_name: Mananhina/mirea-progs


text:
module StackViews

export StackView

using OffsetArrays

const SlicesType = Union{Tuple, AbstractArray{<:AbstractArray}}

"""
    StackView(slices...; [dims])
    StackView(slices, [dims])
    StackView{T}(args...; kwargs...)

Stack/concatenate a list of arrays `slices` along dimension `dims` without copying the data.

If not specified, `dims` is defined as `ndims(first(slices))+1`, i.e., a new dimension in the tail.
This works better for normal Julia arrays when their memory layout is column-major.

# Example

`StackView` works very similar to `cat` and its `vcat`/`hcat`/`hvcat` variants.

```jldoctest; setup=:(using StackViews)
julia> A = reshape(collect(1:6), 2, 3);

julia> B = reshape(collect(7:12), 2, 3);

julia> StackView(A, B, dims=3) # mostly equivalent to `cat(A, B, dims=3)`
2×3×2 StackView{$Int, 3, 3, $(Tuple{Matrix{Int}, Matrix{Int}})}:
[:, :, 1] =
 1  3  5
 2  4  6

[:, :, 2] =
 7   9  11
 8  10  12
```

The difference comes to when `dim <= ndims(first(slices))`, which is `dim <= 2` for this example.
`StackView` always creats new dimension while `cat`s do not.

```jldoctest; setup=:(using StackViews; A = reshape(collect(1:6), 2, 3); B = reshape(collect(7:12), 2, 3);)
julia> StackView(A, B, dims=1) # `cat(A, B, dims=1)` outputs 4×3 Matrix
2×2×3 StackView{$Int, 3, 1, $(Tuple{Matrix{Int}, Matrix{Int}})}:
[:, :, 1] =
 1  2
 7  8

[:, :, 2] =
 3   4
 9  10

[:, :, 3] =
  5   6
 11  12

julia> StackView(A, B, dims=2) # `cat(A, B, dims=2)` outputs 2×6 Matrix
2×2×3 StackView{$Int, 3, 2, $(Tuple{Matrix{Int}, Matrix{Int}})}:
[:, :, 1] =
 1  7
 2  8

[:, :, 2] =
 3   9
 4  10

[:, :, 3] =
 5  11
 6  12
```

!!! tip
    For type-stability, you can pass `Val` as dim, e.g., `StackView(slices, Val(3))`.
"""
struct StackView{T, N, D, A} <: AbstractArray{T, N}
    slices::A
    function StackView{T, N, D}(slices) where {T, N, D}
        all(A->axes(A) == axes(first(slices)), slices) || throw(ArgumentError("all slices should be of the same axes."))
        new{T, N, D, typeof(slices)}(slices)
    end
    StackView{T, N, 0}(::A) where {T, N, A} = throw(ArgumentError("`dims=Val(0)` is not supported, do you mean `dims=Val(1)`?"))
end
StackView(slices::AbstractArray...; dims=Val(_default_dims(slices))) = StackView(slices, dims)
StackView{T}(slices::AbstractArray...; dims=Val(_default_dims(slices))) where T = StackView{T}(slices, dims)
StackView(slices, dims::Int) = StackView(slices, Val(dims)) # type-unstable

function StackView(slices::SlicesType, dims::Val = Val(_default_dims(slices)))
   return StackView{_default_eltype(slices)}(slices, dims)
end
function StackView{T}(slices::SlicesType, dims::Val = Val(_default_dims(slices))) where T
    N = _max(dims, Val(_default_dims(slices)))
    # unify all the axes to 1-based ranges
    slices = map(OffsetArrays.no_offset_view, slices)
    return StackView{T, _value(N), _value(dims)}(slices)
end

@inline _default_dims(slices) = ndims(first(slices)) + 1
@inline function _default_eltype(slices)
    T = mapreduce(eltype, promote_type, slices)
    _isconcretetype(T) || throw(ArgumentError("Input arrays should be homogenous."))
    return T
end
function _isconcretetype(T)
    # We relax the restriction and allow `Union`
    # This is particularily useful for arrays with `missing` and `nothing`
    isconcretetype(T) && return true
    isa(T, Union) && return isconcretetype(T.a) && _isconcretetype(T.b)
    return false
end

function Base.size(A::StackView{T,N,D}) where {T,N,D}
    frame_size = size(first(A.slices))
    prev, post = Base.IteratorsMD.split(frame_size, Val(D-1))
    return (_append_tuple(prev, Val(D-1))..., length(A.slices), post...)
end

function Base.axes(A::StackView{T,N,D}) where {T,N,D}
    frame_axes = axes(first(A.slices))
    prev, post = Base.IteratorsMD.split(frame_axes, Val(D-1))

    # use homogenous range to make _append_tuple happy
    fill_range = convert(typeof(first(frame_axes)), Base.OneTo(1))
    return (_append_tuple(prev, Val(D-1), fill_range)...,
            Base.OneTo(length(A.slices)),
            post...)
end

@inline function Base.getindex(A::StackView{T,N,D}, inds::Vararg{Int,N}) where {T,N,D}
    @boundscheck checkbounds(A, inds...)
    prev, post = Base.IteratorsMD.split(inds, Val(D-1))
    idx, post = first(post), Base.tail(post)
    return @inbounds A.slices[idx][prev..., post...]
end

@inline function Base.setindex!(A::StackView{T,N,D}, x, inds::Vararg{Int, N}) where {T,N,D}
    @boundscheck checkbounds(A, inds...)
    prev, post = Base.IteratorsMD.split(inds, Val(D-1))
    idx, post = first(post), Base.tail(post)
    @inbounds A.slices[idx][prev..., post...] = x
end

# utils

# For type stability
@inline _max(::Val{x}, ::Val{y}) where {x, y} = Val(max(x, y))
@inline _value(::Val{N}) where N = N
@inline _append_tuple(t::NTuple{N1}, ::Val{N1}, x=1) where N1 = t
@inline _append_tuple(t::NTuple{N1}, ::Val{N2}, x=1) where {N1, N2} = _append_tuple((t..., x), Val(N2), x)
end


text:
function getIndexes(i,j,d)
    if i == j
        return [i]
    elseif i > j && i > 2
        return 0
    else
        niz1 = []
        niz2 = collect(i:j-1)
        for k in 1:d
            if !(k in niz2)
                push!(niz1,k)
            end
        end
        return length(niz1) == 0 ? 0 : niz1
    end
end

function getSelectedColumns(selectedColumns,X)
    localX = []
    if selectedColumns != 0
        for i in selectedColumns
            push!(localX,X[:,i])
        end
        return localX
    else
        return 0
    end
end

X = [1 2 3 4 5;
     1 2 3 4 5;
     1 2 3 4 5;
     1 2 3 4 5]

d = 10

for i in 1:d-1
    for j in 1:d
        temp = getIndexes(i,j,d)
        println(temp == 0 ? "" : temp)
    end
end

max_stars_repo_path: Primer Predavanje/testing.jl

max_stars_repo_name: pikacxe/Primenjeni-Algoritmi

max_issues_repo_path: Primer Predavanje/testing.jl

max_issues_repo_name: pikacxe/Primenjeni-Algoritmi

max_forks_repo_path: Primer Predavanje/testing.jl

max_forks_repo_name: pikacxe/Primenjeni-Algoritmi


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

"""
    abstract type AbstractPosteriorDensity <: AbstractDensity end

Abstract type for posterior probability densities.
"""
abstract type AbstractPosteriorDensity <: AbstractDensity end
export AbstractPosteriorDensity

"""
    getlikelihood(posterior::AbstractPosteriorDensity)::AbstractDensity

*BAT-internal, not part of stable public API.*

The likelihood density of `posterior`. The likelihood may or may not be
normalized.
"""
function getlikelihood end

"""
    getprior(posterior::AbstractPosteriorDensity)::AbstractDensity

*BAT-internal, not part of stable public API.*

The prior density of `posterior`. The prior may or may not be normalized.
"""
function getprior end

function eval_logval(density::AbstractPosteriorDensity, v::Any, T::Type{<:Real})
    v_shaped = fixup_variate(varshape(density), v)
    R = density_logval_type(v_shaped, T)

    prior_logval = eval_logval(getprior(density), v_shaped, R)

    # Don't evaluate likelihood if prior probability is zero. Prevents
    # failures when algorithms try to explore parameter space outside of
    # definition of likelihood (as long as prior is chosen correctly).
    if !is_log_zero(prior_logval, R)
        likelihood_logval = eval_logval(getlikelihood(density), v_shaped, R)
        convert(R, likelihood_logval + prior_logval)
    else
        log_zero_density(R)
    end
end

eval_logval_unchecked(density::AbstractPosteriorDensity, v::Any) = eval_logval(density, v, default_dlt())

function var_bounds(density::AbstractPosteriorDensity)
    li_bounds = var_bounds(getlikelihood(density))
    pr_bounds = var_bounds(getprior(density))
    if ismissing(li_bounds)
        pr_bounds
    else
        li_bounds ∩ pr_bounds
    end
end

"""
    struct PosteriorDensity{
        Li<:AbstractDensity,
        Pr<:DistLikeDensity,
        ...
    } <: AbstractPosteriorDensity

A representation of a PosteriorDensity, based a likelihood and prior.
Likelihood and prior be accessed via

```julia
getlikelihood(posterior::PosteriorDensity)::Li
getprior(posterior::PosteriorDensity)::Pr
```

Constructors:

* ```PosteriorDensity(likelihood, prior)```

`likelihood` and `prior` must be convertible to an [`AbstractDensity`](@ref).

Fields:

$(TYPEDFIELDS)

!!! note

    Fields `parbounds` and `parbounds` do not form part of the stable public
    API and are subject to change without deprecation.
"""
struct PosteriorDensity{
    L<:AbstractDensity,
    P<:AbstractDensity,
    B<:AbstractVarBounds,
    S<:AbstractValueShape
} <: AbstractPosteriorDensity
    likelihood::L
    prior::P
    parbounds::B
    parshapes::S
end

export PosteriorDensity

function PosteriorDensity(likelihood::Any, prior::Any)
    li = convert(AbstractDensity, likelihood)
    pr = convert(AbstractDensity, prior)

    parbounds = _posterior_parbounds(var_bounds(li), var_bounds(pr))

    li_shape = varshape(li)
    pr_shape = varshape(pr)
    parshapes = _posterior_parshapes(li_shape, pr_shape)

    li_with_shape = _density_with_shape(li, parshapes, li_shape)

    PosteriorDensity(li_with_shape, pr, parbounds, parshapes)
end

getlikelihood(posterior::PosteriorDensity) = posterior.likelihood

getprior(posterior::PosteriorDensity) = posterior.prior

var_bounds(posterior::PosteriorDensity) = posterior.parbounds

ValueShapes.varshape(posterior::PosteriorDensity) = posterior.parshapes

ValueShapes.unshaped(density::PosteriorDensity) = PosteriorDensity(unshaped(density.likelihood), unshaped(density.prior))

(shape::AbstractValueShape)(density::PosteriorDensity) = PosteriorDensity(shape(density.likelihood), shape(density.prior))

function _posterior_parshapes(li_ps::AbstractValueShape, pr_ps::AbstractValueShape)
    if li_ps == pr_ps
        li_ps
    else
        throw(ArgumentError("Variable shapes of likelihood and prior are incompatible"))
    end
end

function _posterior_parshapes(li_ps::AbstractValueShape, pr_ps::ArrayShape{T,1}) where T
    n = totalndof(li_ps)
    if n == totalndof(pr_ps)
        @assert size(pr_ps) == (n,)
        li_ps
    else
        throw(ArgumentError("Likelihood and prior have different number of free parameters"))
    end
end

_posterior_parshapes(li_ps::ArrayShape{T,1}, pr_ps::AbstractValueShape) where T =
    _posterior_parshapes(pr_ps, li_ps)

function _posterior_parshapes(li_ps::ArrayShape{T,1}, pr_ps::ArrayShape{U,1}) where {T,U}
    if size(li_ps) == size(pr_ps)
        li_ps
    else
        throw(ArgumentError("Likelihood and prior have different number of free parameters"))
    end
end

_posterior_parshapes(li_ps::Missing, pr_ps::AbstractValueShape) = pr_ps

_posterior_parshapes(li_ps::Union{AbstractValueShape,Missing}, pr_ps::Missing) =
    throw(ArgumentError("Parameter shapes of prior must not be missing"))

_posterior_parbounds(li_bounds::AbstractVarBounds, pr_bounds::AbstractVarBounds) =
     li_bounds ∩ pr_bounds

_posterior_parbounds(li_bounds::Missing, pr_bounds::AbstractVarBounds) = pr_bounds

function _density_with_shape(density::AbstractDensity, requested_shape::AbstractValueShape, orig_shape::AbstractValueShape)
    if requested_shape == orig_shape
        density
    else
        throw(ArgumentError("Original and requested variable shape are incompatible"))
    end
end

function _density_with_shape(density::AbstractDensity, requested_shape::AbstractValueShape, orig_shape::Missing)
    DensityWithShape(density, requested_shape)
end

function example_posterior()
    prior = NamedTupleDist(
        a = Exponential(),
        b = [4.2, 3.3],
        c = Normal(1, 3),
        d = [Weibull(), Weibull()],
        e = Beta(),
        f = MvNormal([0.3, -2.9], [1.7 0.5; 0.5 2.3])
    )
    n = totalndof(varshape(prior))
    likelihood = varshape(prior)(MvNormal(float(I(n))))
    PosteriorDensity(likelihood, prior)
end