text:
# ____________________________________________________________________
# SpecLineGauss
#
mutable struct SpecLineGauss <: AbstractComponent
    norm::Parameter
    center::Parameter
    fwhm::Parameter
    voff::Parameter
    spec_res_kms::Float64
    norm_integrated::Bool

    function SpecLineGauss(center::Number)
        out = new(Parameter(1),
                  Parameter(center),
                  Parameter(3000),
                  Parameter(0),
                  0., true)

        @assert center > 0
        out.norm.low = 0
        out.center.low = 0
        out.fwhm.low = 0
        out.voff.low = 0
        out.center.fixed = true
        return out
    end
end

function prepare!(comp::SpecLineGauss, domain::Domain{1})
    return fill(NaN, length(domain))
end

function evaluate!(buffer, comp::SpecLineGauss, x::Domain{1},
                   norm, center, fwhm, voff)
    x0 = center - (voff / 3.e5) * center
    sigma_line = fwhm / 2.355      / 3.e5 * center
    sigma_spec = comp.spec_res_kms / 3.e5 * center
    sigma = sqrt(sigma_line^2 + sigma_spec^2)

    map!(x -> begin
         (abs.(x - x0) .> 4sigma)  &&  (return 0.)
         ret = norm * exp(-((x - x0) / sigma)^2 / 2)
         comp.norm_integrated  &&  (ret /= (sqrt(2pi) * sigma))
         return ret
         end, buffer, x.coords[:,1])
end

#=
    x = Domain(500:1:1500.)
    comp = QSFit.SpecLineGauss(1000.)
    comp.fwhm.val = 3e4
    ceval = GFit.CompEval(comp, x)
    GFit.evaluate_cached(ceval)
    @gp x[:] ceval.buffer ./ maximum(ceval.buffer) "w l"
=#


text:
function _add_irt!(ata_model::AbstractModel)
    try
        if ata_model.settings.irt.model == "1PL"
            ata_model.settings.irt.parameters = DataFrames.DataFrame(ata_model.settings.bank[!, [:b]])#nqp values in interval\r\n",
        elseif ata_model.settings.irt.model == "2PL"
            ata_model.settings.irt.parameters = DataFrames.DataFrame(ata_model.settings.bank[!, [:a, :b]]) #nqp values in interval\r\n",
        elseif ata_model.settings.irt.model == "3PL"
            ata_model.settings.irt.parameters = DataFrames.DataFrame(ata_model.settings.bank[!, [:a, :b, :c]]) #nqp values in interval\r\n",
        else
            error!(ata_model, "Only 1PL, 2PL and 3PL IRT models are allowed.")
        end
        CSV.write("opt/irt_parameters.csv", ata_model.settings.irt.parameters)
    catch e

    end
    success!(ata_model, "IRT item parameters loaded.")
    return nothing
end

ext: jl


text:
using DFControl, Test, LinearAlgebra

testdir = joinpath(dirname(dirname(pathof(DFControl))), "test")
testjobpath = joinpath(testdir, "testassets", "test_job")
@testset "initial creation" begin
    if ispath(testjobpath)
        rm(testjobpath, recursive=true)
    end

    name = "Ni"
    dir = testjobpath
    bin_dir = joinpath(homedir(), "Software/qe/bin")
    pw_exec = Exec("pw", "pw.x", bin_dir, :nk => 4)
    @test !isempty(pw_exec.flags)
    pseudoset = :test

    header = ["#SBATCH -N 1"]

    str = Structure(joinpath(testdir, "testassets/Ni.cif"))

    calculations = [Calculation{QE}("vcrelax", :calculation => "vc-relax", :verbosity => "high", :ion_dynamics => "bfgs", :cell_dynamics => "bfgs";
                                      exec = pw_exec,
                                      data = [InputData(:k_points, :automatic,
                                                        [6, 6, 6, 1, 1, 1])]),
                    Calculation{QE}(; name = "scf", exec = pw_exec,
                                      flags = Dict(:calculation => "scf", :verbosity => "high"),
                                      data = [InputData(:k_points, :automatic,
                                                        [4, 4, 4, 1, 1, 1])])]
    job = Job(name, str, calculations, :ecutwfc => 40.0, :occupations => "smearing", :degauss=>0.01, :conv_thr => 1e-6, :nbnd => 18;
                #kwargs
                header = header, dir = dir, server="localhost_test")

    set_pseudos!(job, :test)
    job.environment = "test_default"

    set_kpoints!(job["scf"], (6, 6, 6, 1, 1, 1))

    job.structure[element(:Ni)][1].magnetization = Vec3(0,0,0.1)
    job.structure[element(:Ni)][1].dftu.U = 4.0

    push!(job, Calculations.gencalc_bands(job["scf"], Structures.high_symmetry_kpath(job.structure, 20)))
    push!(job, Calculations.gencalc_nscf(job["scf"], (5,5,5)))

    push!(job, Calculations.gencalc_projwfc(job["nscf"], 2.0, 20.0, 0.1))
    @show isempty(job["projwfc"].exec.flags)

    @test all(values(job[:ecutwfc]) .== 40.0)
    for c in job.calculations
        pop!(c, :ecutwfc, nothing)
    end

    save(job)
    @test all(values(job[:ecutwfc]) .== 41.0)
    @test all(values(job[:ecutrho]) .== 236.0)
    @test job.version == 1
    @test length(job) == 5
    @test data(job["scf"], :k_points).data == [6,6,6,1,1,1]
    @test job["nscf"].exec == pw_exec
    @test job["projwfc"].exec.exec == "projwfc.x"
    @test job["projwfc"].exec.dir == pw_exec.dir
    @test show(job) == nothing
    job[:ecutwfc] = 40.0
    for c in job.calculations
        pop!(c, :ecutrho, nothing)
    end
    save(job)

    job2 = Job(abspath(job), "localhost_test")
    for (c1, c2) in zip(job2.calculations, job.calculations)
        @test c2 == c1
    end
    @test job2.structure == job.structure
end

refjobpath =joinpath(testdir, "testassets", "reference_job")

@testset "reference comparison" begin
    job = Job(testjobpath, "localhost_test")

    orig_job = deepcopy(job)
    job.structure = Structures.create_supercell(job.structure, 1, 0, 0, make_afm = true)

    job2 = Job(refjobpath, "localhost_test")
    @test job2.structure == job.structure

    for f in DFControl.Utils.searchdir(job2, ".out")
        cp(f, joinpath(job, splitdir(f)[2]), force=true)
    end
    for f in DFControl.Utils.searchdir(job2, "dos")
        cp(f, joinpath(job, splitdir(f)[2]), force=true)
    end

    for a in job.structure.atoms
        a.projections = [Projection("s"), Projection("p"), Projection("d")]
    end
    wanexec = Exec("wan","wannier90.x", joinpath(homedir(), "Software/wannier90"))
    append!(job, Calculations.gencalc_wan(job, 0.000011, wanexec = wanexec))
    for (c1, c2) in zip(job2.calculations, job.calculations)
        @test c2 == c1
    end
    save(job)
    @test !ispath(joinpath(job, "scf.out"))
    job = Job(testjobpath, "localhost_test")

    for (c1, c2) in zip(job2.calculations, job.calculations)
        @test c2 == c1
    end
    save(orig_job)
    for f in DFControl.Utils.searchdir(job2, ".out")
        cp(f, joinpath(job, splitdir(f)[2]), force=true)
    end
    for f in DFControl.Utils.searchdir(job2, "dos")
        cp(f, joinpath(job, splitdir(f)[2]), force=true)
    end
end


text:
examples = [
    "linearreg",
    "logireg",
    "mlogireg",
    "lasso"
]

for e in examples
    f = "$e.jl"
    println("Running $f ...")
    include(f)
    println()
end


text:
importall ReinforcementLearning
import Base.==
import Base.isequal
import Base.hash

pkgs = ["ArgParse"]
for pkg in pkgs; Pkg.installed(pkg) == nothing && Pkg.add(pkg); end

using ArgParse

function parse_commandline()
	s = ArgParseSettings()

	@add_arg_table s begin
		"--epochs"
			help = "number of epochs"
			default = 20
			arg_type=Int
		"--eps"
			help = "epsilon parameter for the random behaviour"
			default = 0.05
			arg_type = Float64
		"--gamma"
			help = "the discount factor"
			default = 0.9
			arg_type = Float64
		"--alpha"
			help = "learning step"
			default = 0.8
			arg_type = Float64
		"--lambda"
			help = "td lambda"
			default = 0.5
			arg_type = Float64

		"--dims"
			help = "dimensions of the maze"
			nargs='+'
			default = [3,3]
			arg_type = Int
		"--agent"
			help = "Agent type"
			default = "Q"
	end
	return parse_args(s)
end

#recursive backtracking algorithm
function generate_maze(h = 4, w = 4)
	maze = zeros(h, w, 4)
	unvisited = ones(h, w)

	function neighbours(r,c)
		ns = Array{Tuple{Int, Int, Int}, 1}()
		for i=1:4
			if i == 1 && (r - 1) >= 1 && unvisited[r-1, c] == 1
				push!(ns, (r-1, c, 1))
			elseif i == 2 && (c + 1) <= w && unvisited[r, c+1] == 1
				push!(ns, (r, c+1, 2))
			elseif i == 3 && (r+1) <= h && unvisited[r+1, c] == 1
				push!(ns, (r+1, c, 3))
			elseif i == 4 && (c-1) >= 1 && unvisited[r, c-1] == 1
				push!(ns, (r, c-1, 4))
			end
		end
		return shuffle(ns)
	end

	start = rand(1:h*w)
	r = div(start, h) + 1
	r = r > h ? h : r
	c = start % w
	c = c == 0 ? w : c

	stack = Array{Tuple{Int, Int}, 1}()
	curr = (r, c)
	unvisited[r,c] = 0
	while countnz(unvisited) != 0
		ns = neighbours(curr[1], curr[2])
		if length(ns) > 0
			r,c = curr
			rn, cn, d = ns[1]
			push!(stack, (r,c))
			maze[r, c, d] = 1
			dn = d - 2 <= 0 ? d + 2 : d - 2
			maze[rn, cn, dn] = 1
			curr = (rn, cn)
			unvisited[rn, cn] = 0
		elseif length(stack) != 0
			curr = pop!(stack)
		end

	end
	return maze
end

function print_maze(maze)
	h,w,_ = size(maze)
	rows = 2*h + 1
	cols = 2*w + 1

	for i=1:rows
		println("")
		for j=1:cols
			if i == 1 || i == rows || j == 1 || j == cols
				print("#")
			elseif i % 2 == 1 && j % 2 == 1
				print("#")
			elseif i % 2 == 1 && j % 2 == 0
				r = div(i - 1, 2)
				c = div(j, 2)
				if maze[r, c, 3] == 1
					print(" ")
				else
					print("#")
				end
			elseif i % 2 == 0 && j % 2 == 1
				r = div(i, 2)
				c = div(j - 1, 2)
				if maze[r,c,2] == 1
					print(" ")
				else
					print("#")
				end
			else
				print(" ")
			end

		end
	end
	print("\n")
end

#=======USAGE==============
h,w=(3,3)
maze = generate_maze(h, w)

print_maze(maze)
=========================#

mutable struct MazeState <: AbsState
	loc
end

@enum ActionEnum UP RIGHT DOWN LEFT

mutable struct MazeAction <: AbsAction; act::ActionEnum; end

==(lhs::MazeAction, rhs::MazeAction) = lhs.act == rhs.act
isequal(lhs::MazeAction, rhs::MazeAction) = lhs.act == rhs.act
hash(a::MazeAction) = hash(a.act)

==(lhs::MazeState, rhs::MazeState) = (lhs.loc[1] == rhs.loc[1] && lhs.loc[2] == rhs.loc[2])
isequal(lhs::MazeState, rhs::MazeState) = ==(lhs, rhs)
hash(s::MazeState) = hash(s.loc)

hash(t::Tuple{MazeState, MazeAction}) = hash([t[1].loc, t[2].act])

mutable struct MazeEnv <: AbsEnvironment
	dims
	maze
	start
	goal

	MazeEnv(dims) = new(dims, generate_maze(dims[1], dims[2]), (rand(1:dims[1]), 1), (rand(1:dims[1]), dims[2]))
	MazeEnv() = Maze((5, 5))
end

getInitialState(env::MazeEnv) = MazeState(env.start)
getActions(state::MazeState, env::MazeEnv) = [MazeAction(UP), MazeAction(RIGHT), MazeAction(DOWN), MazeAction(LEFT)]

function getAllStates(env::MazeEnv)
	states = MazeState[MazeState((0,0))]
	for i=1:env.dims[1]
		for j=1:env.dims[2]
			push!(states, MazeState((i, j)))
		end
	end
	return states
end

function isTerminal(state::MazeState, env::MazeEnv)
	(state.loc[1] == 0 && state.loc[2] == 0) ||
		(state.loc[1] == env.goal[1] && state.loc[2] == env.goal[2])
end

function getSuccessors(state::MazeState, a::MazeAction, env::MazeEnv)
	successors = Array{Tuple{MazeState, Float64, Float64}, 1}()
	push!(successors, transfer(env, state, a))
	return successors
end

function transfer(env::MazeEnv, state::MazeState, action::MazeAction)
	loc = state.loc
	reward = -1.0
	next_state = nothing

	if (loc[1] == 0 && loc[2] == 0)
		return (state, -10.0 * prod(size(env.maze)[1:2]), 1.0)
	elseif (loc[1] == env.goal[1] && loc[2] == env.goal[2])
		return (state, prod(size(env.maze)[1:2])*10.0, 1.0)
	end

	if action.act == UP
		orientation = 1
	elseif action.act == RIGHT
		orientation = 2
	elseif action.act == DOWN
		orientation = 3
	else
		orientation = 4
	end

	if env.maze[loc[1], loc[2], orientation] == 0
		reward = -10.0 * prod(size(env.maze)[1:2])
		next_state = MazeState((0,0))
	else
		if orientation == 1
			loc = (loc[1] - 1, loc[2])
		elseif orientation == 2
			loc = (loc[1], loc[2] + 1)
		elseif orientation == 3
			loc = (loc[1] + 1, loc[2])
		else
			loc = (loc[1], loc[2] - 1)
		end

		next_state = MazeState(loc)
		reward = (loc[1] == env.goal[1] && loc[2] == env.goal[2]) ? prod(size(env.maze)[1:2])*10.0 : -1.0
	end
	return (next_state, reward, 1.0)
end

function agent_experiment(agent, env, optimumReward, threshold, steps, rInitial)
	totalRewards = 0.0
	rewards = Any[]
	nstates = Any[]

	@time while totalRewards != optimumReward
		while totalRewards < threshold
			totalRewards, numberOfStates = playEpisode(env, agent; randomInitial=rInitial, threshold = steps)
			push!(rewards, totalRewards)
			push!(nstates, numberOfStates)
		end
		totalRewards, numberOfStates = playEpisode(env, agent; learn=false, threshold = steps)
	end

	println("Training took $(length(rewards)) episodes")

	#=
	for i=1:length(rewards)
		println("Epoch: $i, totalReward: $(rewards[i]), nsteps: $(nsteps[i])")
	end
	=#

	totalRewards, numberOfStates = playEpisode(env, agent; learn=false, threshold = steps)
	println("After training: $totalRewards, #states: $numberOfStates")

end

function main()
	args = parse_commandline()
	hm, wm = args["dims"]
	env = MazeEnv((hm, wm))

	print_maze(env.maze)
	println("Start: $(env.start)")
	println("Goal: $(env.goal)")

	ss = getAllStates(env)

	println("Value Iteration")
	@time policy, V = synchronous_value_iteration(env; Ɣ=args["gamma"], verbose=false)

	policyAgent = PolicyAgent(policy)
	optimumReward, numberOfStates = playEpisode(env, policyAgent; threshold = prod(args["dims"]) * 2, verbose=false)

	println("Optimum Reward: $optimumReward")

	#=
	for s in ss
		println("State: $(s), Value: $(V[s]), Action: $(policy.mapping[s][1][1])")
	end
	=#

	println("Policy Iteration")
	@time policy, V = policy_iteration(env::AbsEnvironment; Ɣ=args["gamma"], verbose=false)

	policyAgent = PolicyAgent(policy)
	optimumReward, numberOfStates = playEpisode(env, policyAgent; threshold = prod(args["dims"]) * 2)

	println("Optimum Reward: $optimumReward")

	#=
	for s in ss
		println("State: $(s), Value: $(V[s]), Action: $(policy.mapping[s][1][1])")
	end
	=#

	println("Gauss-Seidel Value Iteration")
	@time policy, V = gauss_seidel_value_iteration(env; Ɣ=args["gamma"], verbose=false)

	#=
	for s in ss
		println("State: $(s), Value: $(V[s]), Action: $(policy.mapping[s][1][1])")
	end
	=#

	policyAgent = PolicyAgent(policy)
	optimumReward, numberOfStates = playEpisode(env, policyAgent; threshold = prod(args["dims"]) * 2)

	println("Optimum Reward: $optimumReward")

	threshold = optimumReward - (2 * prod(args["dims"]) * args["eps"])
	steps = prod(args["dims"]) * 2

	println("\nQ Learning")
	agentQ = QLearner(env;ε=args["eps"], α=args["alpha"], Ɣ=args["gamma"])
	agent_experiment(agentQ, env, optimumReward, threshold, steps, false)

	println("\nSarsa")
	agentSarsa = SarsaLearner(env;ε=args["eps"], α=args["alpha"], Ɣ=args["gamma"])
	agent_experiment(agentSarsa, env, optimumReward, threshold, steps, false)

	println("\nSarsa Lambda")
	agentSarsaLambda = SarsaLambdaLearner(env; λ=args["lambda"], ε=args["eps"], α=args["alpha"], Ɣ=args["gamma"])
	agent_experiment(agentSarsaLambda, env, optimumReward, threshold, steps, false)

	println("\nMonte Carlo")
	agentMC = MonteCarloAgent(env)
	agent_experiment(agentMC, env, optimumReward, threshold, steps, true)
end

if VERSION >= v"0.5.0-dev+7720"
    PROGRAM_FILE=="maze.jl" && main()
else
    !isinteractive() && !isdefined(Core.Main,:load_only) && main()
end

max_stars_repo_name: ozanarkancan/ReinforcementLearning

max_issues_repo_name: ozanarkancan/ReinforcementLearning

max_forks_repo_name: ozanarkancan/ReinforcementLearning


text:
export SARTS
export SART
export SA
export RLDataSet
export fetch_gc_bucket
export fetch_gc_file

abstract type RLDataSet end

"""
`(:state, :action, :reward, :terminal, :next_state)`
type of the returned batches.
"""
const SARTS = (:state, :action, :reward, :terminal, :next_state)

"""
`(:state, :action, :reward, :terminal)`
type of the returned batches.
"""
const SART = (:state, :action, :reward, :terminal)

"""
`(:state, :action)`
type of the returned batches.
"""
const SA = (:state, :action)

# fetch_methods
"""
fetch a gc bucket from `src` to `dest`.
"""
function fetch_gc_bucket(src, dest)
    if Sys.iswindows()
        try run(`cmd /C gsutil -v`) catch x throw("gsutil not found, install gsutil to proceed further") end
        run(`cmd /C gsutil -m cp -r $src $dest`)
    else
        try run(`gsutil -v`) catch x throw("gsutil not found, install gsutil to proceed further") end
        run(`gsutil -m cp -r $src $dest`)
    end
    return dest
end

"""
fetch a gc file from `src` to `dest`.
"""
function fetch_gc_file(src, dest)
    if Sys.iswindows()
        try run(`cmd /C gsutil -v`) catch x throw("gsutil not found, install gsutil to proceed further") end
        run(`cmd /C gsutil -m cp $src $dest`)
    else
        try run(`gsutil -v`) catch x throw("gsutil not found, install gsutil to proceed further") end
        run(`gsutil -m cp $src $dest`)
    end
    return dest
end


text:
using Catalyst, DiffEqBase, ModelingToolkit, Test

@parameters t k1 k2
@variables S I R
rxs = [Reaction(k1, [S,I], [I], [1,1], [2]),
       Reaction(k2, [I], [R]) ]
rs = ReactionSystem(rxs, t, [S,I,R], [k1,k2])

specset = Set([S.op =>1, I.op => 2, R.op => 3])
@test issetequal(specset, speciesmap(rs))

pset = Set([k1.op => 1, k2.op => 2])
@test issetequal(pset, paramsmap(rs))

rxs2 = [Reaction(k2, [I], [R], [1], [1]),
        Reaction(k1, [S,I], [I], [1,1], [2])]
rs2 = ReactionSystem(rxs2, t, [R,I,S], [k2,k1])
@test rs == rs2

rs3 = make_empty_network()
@parameters k3 k4
@variables D
addspecies!(rs3, S)
addspecies!(rs3, D.op)
addparam!(rs3, k3)
addparam!(rs3, k4.op)
@test issetequal(species(rs3), [S.op, D.op])
@test issetequal(params(rs3), [k3.op, k4.op])
addreaction!(rs3, Reaction(k3, [S], [D]))
addreaction!(rs3, Reaction(k4, [S,I], [D]))
merge!(rs, rs3)
addspecies!(rs2, S)
addspecies!(rs2, D.op)
addparam!(rs2, k3)
addparam!(rs2, k4.op)
addreaction!(rs2, Reaction(k3, [S], [D]))
addreaction!(rs2, Reaction(k4, [S,I], [D]))
@test rs2 == rs

rxs = [Reaction(k1, [S,I], [I], [1,1], [2]),
       Reaction(k2, [I], [R]) ]
rs = ReactionSystem(rxs, t, [S,I,R], [k1,k2])
rs3 = make_empty_network()
addspecies!(rs3, S)
addspecies!(rs3, D.op)
addparam!(rs3, k3)
addparam!(rs3, k4.op)
addreaction!(rs3, Reaction(k3, [S], [D]))
addreaction!(rs3, Reaction(k4, [S,I], [D]))
rs4 = merge(rs, rs3)
@test rs2 == rs4

rxs = [Reaction(k1*S, [S,I], [I], [2,3], [2]),
       Reaction(k2*R, [I], [R]) ]
rs = ReactionSystem(rxs, t, [S,I,R], [k1,k2])
deps = [s.op for s in dependents(rxs[2], rs)]
@test isequal(deps, [R.op,I.op])
@test isequal(dependents(rxs[1], rs), dependants(rxs[1], rs))
addspecies!(rs, Variable(:S))
@test numspecies(rs) == 3
addspecies!(rs, Variable(:S), disablechecks=true)
@test numspecies(rs) == 4
addparam!(rs, Variable(:k1))
@test numparams(rs) == 2
addparam!(rs, Variable(:k1), disablechecks=true)
@test numparams(rs) == 3


text:
"""
  # ComplexPortraits
"""
module ComplexPortraits

using Colors, RecipesBase

"""macro for importing the *huge* set of symbols."""
macro import_huge()
  :(
    using ComplexPortraits: generate_stepfct, generate_sawfct,
                            generate_sawfctx, generate_rectfct,
                            generate_spike,
                            angle_zero_one, brightenRGB,
                            hsv_colors, nist_colors,
                            cs_c,
                            cs_m, cs_e,
                            cs_j,
                            cs_p, cs_q, cs_d,
                            cs_gridReIm_lowsat,
                            cs_gridReIm, cs_gridReIm_logabs, cs_gridReIm_abs,
                            cs_stripes, cs_angle_abs_stripes,
                            cs_useImage,
                            portrait,
                            phaseplot, phaseplot!, PhasePlot
  )
end

"""macro for importing the normal set of symbols."""
macro import_normal()
  :(
    using ComplexPortraits: hsv_colors, nist_colors,
                            cs_c,
                            cs_m, cs_e,
                            cs_j,
                            cs_p, cs_q, cs_d,
                            cs_gridReIm, cs_gridReIm_logabs, cs_gridReIm_abs,
                            cs_stripes, cs_angle_abs_stripes,
                            cs_useImage,
                            portrait,
                            phaseplot, phaseplot!, PhasePlot
  )
end

include("./utils.jl")
include("./schemes.jl")

"""
Generate a 'portrait' of a complex function.

Depending on the point_color argument different color schemes are possible.

z_upperleft and z_lowerright define the rectangle for the portrait
in the complex plane for the function f.

Additional arguments are the number of pixels and the point_color function.

Returns an Matrix with the (pixel-)colors as entries.
"""
function portrait(z_upperleft, z_lowerright, f;
                  no_pixels=(600,600), point_color=cs_j())
  if length(no_pixels) == 1
    no_pixels = (no_pixels[1], no_pixels[1])
  end
  return [
    let z = x + y*1im
      point_color(z, f(z))
    end
    for y in LinRange(imag(z_upperleft), imag(z_lowerright), no_pixels[1]),
        x in LinRange(real(z_upperleft), real(z_lowerright), no_pixels[2])]
end

include("./plotrecipes.jl")

end

# vim:syn=julia:cc=79:fdm=indent:sw=2:ts=2:


text:
####################################################

#  RRT.jl
#  nouyang 2017

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

using Plots
using DataStructures
using GeometryTypes
using Distributions

using plotfxn
using algT
using algfxn

# pts are pure Points, nodes have ID information
function rrtPlan(room, parameters, startstate, goalstate, obstaclesList)
    roomWidth, roomHeight, walls, obstacles = room.width, room.height, room.walls, room.obstacles
    numPts, connectRadius = parameters.numSamples, parameters.connectRadius

    nodeslist = Vector{algT.Node}()
    edgeslist = Vector{algT.Edge}()

    startNode = algT.Node(0, algT.Pt2D(startstate), 0) #id, state, parentid
    push!(nodeslist, startNode)

    currID = 1 #current node ID

    isPathFound = false

    nn = algT.Node(99999, algT.Pt2D(goalstate))
    newNode = nothing
    endNode = nothing

    for i in 1:numPts
            print("\n ------------------------ \n")
        randPt = algfxn.sampleFree(roomWidth, roomHeight, obstacles)

        nn = algT.Node(algfxn.nearestN(randPt, nodeslist))
        newPt = algfxn.steer(nn.state, randPt, connectRadius) #no node ID yet
        newMove = LineSegment(Point(nn.state), Point(newPt))

        if algfxn.isFreeMotion(newMove, obstacles, walls)
            newNode = algT.Node(currID, newPt, nn.id ) # include parent node id
            newEdge = algT.Edge(nn, newNode)
            currID += 1
            push!(nodeslist, newNode)
            push!(edgeslist, newEdge)
            if algfxn.inGoalRegion(newPt, goalstate)
                print("\n --goal found! $newPt, goal is $goalstate-- \n")
                endNode = newNode
                isPathFound = true
                break
            end
        end
    end

    ## Collect Paths
    print("\n\n - pathfound $isPathFound - \n\n")

    solPath = Vector{algT.Node}()

    # TODO handle this more gracefully
    if endNode != nothing
        goalNode = algT.Node(currID+1, goalstate, endNode.id)
    else
        goalNode = algT.Node(currID+1, goalstate)
    end
    push!(solPath, goalNode)

    if isPathFound
        pID = goalNode.parentID
        node = []
        while true
            node = algfxn.findNode(pID, nodeslist)
            pID = node.parentID
            push!(solPath, node)
            if node.id == 0
                break
            end
        end
    end

    #@show edgeslist
    #@show nodeslist
    @show solPath

    finalPathCost = algfxn.costPath(solPath)
    return (nodeslist, edgeslist, isPathFound, solPath, finalPathCost)
end


text:
@testset "SWIPL basics" begin
    prolog = Swipl()

    c = c_const("c")
    p = c_pred("p", 1)
    lit = p(c)

    start(prolog)

    assertz(prolog, lit)

    xvar = c_var("X")

    qlit = p(xvar)

    result = query(prolog, qlit)
    @assert length(result) == 1

    q = c_pred("q", 1)
    a = c_const("a")
    b = c_const("b")

    assertz(prolog, q(a))
    assertz(prolog, q(b))
    assertz(prolog, q(c))

    result2 = query(prolog, q(xvar))

    @assert length(result2) == 3

    rpred = c_pred("r", 2)

    assertz(prolog, rpred(a, b))
    assertz(prolog, rpred(a, c))

    yvar = c_var("Y")

    result3 = query(prolog, rpred(xvar, yvar))

    @assert length(result3) == 2

    f = c_functor("f", 1)

    tpred = c_pred("t", 2)

    assertz(prolog, tpred(a, f(b)))

    results4 = query(prolog, tpred(xvar, yvar))

    @assert length(results4) == 1
    @assert isa(results4[begin][yvar], Structure)

    n = c_pred("n", 1)

    assertz(prolog, n(1))
    assertz(prolog, n(15))
    results5 = query(prolog, n(xvar))

    @assert length(results5) == 2
    @assert isa(results5[begin][xvar], Int64) & isa(results5[begin+1][xvar], Int64)

    results6 = query(prolog, n(xvar); max_solutions=1)
    @assert length(results6) == 1

    number_list = List([1,2,3,4,5])
    member = c_pred("member", 2)

    results7 = query(prolog, member(xvar, number_list))
    @assert length(results7) == 5
    @assert isa(results7[begin][xvar], Int64)

    stop(prolog)

end

@testset "SWIPL with rules" begin
   prolog = Swipl()

   start(prolog)

   edge = c_pred("edge", 2)
   path = c_pred("path", 2)

   xvar = c_var("X")
   yvar = c_var("Y")
   zvar = c_var("Z")

   assertz(prolog, edge(1,2))
   assertz(prolog, edge(2,3))
   assertz(prolog, edge(2,4))
   assertz(prolog, edge(4,5))

   clause = path(xvar, yvar) <= edge(xvar, yvar)

   assertz(prolog, clause)

   result = query(prolog, path(xvar, yvar))

   @assert length(result) == 4

   assertz(prolog, path(xvar, yvar) <= edge(xvar, zvar) & path(zvar, yvar))

   result2 = query(prolog, path(xvar, yvar))

   @assert length(result2) == 8

   stop(prolog)

end


text:
#=
block_diag.jl
Like SparseArrays:blockdiag but slightly different name, just to be safe,
because SparseArrays are also an AbstractArray type.

Motivated by compressed sensing dynamic MRI
where each frame has different k-space sampling
2020-06-08 Jeff Fessler
=#

export block_diag
import LinearMaps
using SparseArrays: blockdiag, sparse

same_idim(As::LinearMapAX...) = all(map(A -> A._idim == As[1]._idim, As))
same_odim(As::LinearMapAX...) = all(map(A -> A._odim == As[1]._odim, As))

"""
    B = block_diag(As::LinearMapAX... ; tryop::Bool)
Make block diagonal `LinearMapAX` object from blocks.

Return a `LinearMapAM` unless `tryop` and all blocks have same `idim` and `odim`.

Default for `tryop` is `true` if all blocks are type `LinearMapAO`.
"""
block_diag(As::LinearMapAO... ; tryop::Bool = true) = block_diag(tryop, As...)

block_diag(As::LinearMapAX... ; tryop::Bool = false) = block_diag(tryop, As...)

function block_diag(tryop::Bool, As::LinearMapAX...)
    B = LinearMaps.blockdiag(map(A -> A._lmap, As)...)
    nblock = length(As)
    prop = (nblock = nblock,)
    if (tryop && nblock > 1 && same_idim(As...) && same_odim(As...)) # operator version
        return LinearMapAA(B ; prop=prop,
            idim = (As[1]._idim..., nblock),
            odim = (As[1]._odim..., nblock),
        )
    end
    return LinearMapAA(B, prop)
end


text:
function relativeerror(D::MaskedEDM{T}, D̂::EDM{T}) where T <: Real
    return norm(D.D - D̂.D)/norm(D.D) * 100
end


text:
@testset "Test optimization problem" begin
    # Define logistic map
    function logistic_f(x_next, x, k, θ)
        x_next .=  θ*x - θ*x.^2
        return
    end
    function logistic_g(y, x, k, θ)
        y .= x
        return
    end
    # Define initial conditions
    x0 = [0.5]
    x0_list = Vector{Vector{Float64}}(10)
    x0_list[1] = copy(x0)
    grid = collect(1:10:101)
    k0_list = grid[1:end-1]
    # Initialize buffer
    y = Vector{Vector{Float64}}(100)
    for i = 1:length(y)
        y[i] = zeros(1)
    end
    # Simulate
    x_final = copy(x0)
    for i = 1:length(grid)-1
        x_final = ms.simulate_space_state!(y[grid[i]:grid[i+1]-1],
                                           logistic_f, logistic_g,
                                           x0, (grid[i], grid[i+1]-1), (3.78,))
        if i != length(grid)-1
            x0 .= x_final
            x0_list[i+1] = copy(x0)
        end
    end
    # Sanity check: Check if multiple simulations corresponds to a single one
    y2 = Vector{Vector{Float64}}(100)
    for i = 1:length(y)
        y2[i] = zeros(1)
    end
    x_final2 = ms.simulate_space_state!(y2, logistic_f, logistic_g,
                                       [0.5], (1, 100), (3.78,))
    @test vcat(y...) == vcat(y2...)

    # Simulate using multiple
    θ = [3.78, 3.78, 1.00]
    # Define polynomial model
    # f(x) = θ1*x - θ2*x^2
    @everywhere function f(y, dx, dθ, x, k, θ)
        y .= θ[1]*x - θ[2]*x.^2
        dx .= θ[1] - 2*θ[2]*x
        dθ .=  [ x  -x.^2  0]
        return
    end
    # g(x) = x
    @everywhere function g(y, dx, dθ, x, k, θ)
        y .= θ[3]*x
        dx .= θ[3]
        dθ .=  [ 0 0 x ]
        return
    end
    opt = ms.OptimizationProblem(f, g, x0_list,
                                  y, k0_list, θ)

    @testset "Test cost funciton" begin
      θ_ext = [3.78; 3.78; 1.00; vcat(x0_list...)]
      @test ms.cost_function(opt, θ_ext) ≈ 0.0
      θ_ext = [3.78; 3.78; 2.00; vcat(x0_list...)]
      @test ms.cost_function(opt, θ_ext) ≈ sum(vcat(y...).^2)
    end

    @testset "Test gradient" begin
        function wrapper_cost(θ_ext)
            return ms.cost_function(opt, θ_ext)
        end

        function wrapper_gradient(θ_ext)
            return ms.gradient(opt, θ_ext)
        end

        function finite_difference_gradient(θ_ext)
            return Calculus.gradient(wrapper_cost, θ_ext, :central)
        end

        for θ in ([3.21, 3.2, 1.00],
                [3.21, 1.2, 1.54],
                [3.56, 3.07, 1.54])
            for x0_ext in (collect(linspace(0.2, 0.3, 10)),
                           collect(linspace(0.1, 0.8, 10)))
                θ_ext = [θ; x0_ext]
                @test isapprox(wrapper_gradient(θ_ext),
                               finite_difference_gradient(θ_ext),
                               rtol=1e-5)
            end
        end
    end

    @testset "Test Hessian" begin
        θ = [3.21, 3.2, 1.00]
        x0_ext = collect(linspace(0.2, 0.3, 10))
        θ_ext = [θ; x0_ext]
        hess = ms.hessian(opt, θ_ext)
        p = ones(13)
        @test length(dot(hess, p)) == 13
    end

    @testset "Test constraints" begin
        θ_ext = [3.78; 3.78; 1.00; vcat(x0_list...)]
        @test ms.constraint(opt, θ_ext) ≈ zeros(9) atol=1e-10
        θ_ext = [3.78; 3.78; 2.00; vcat(x0_list...)]
        @test ms.constraint(opt, θ_ext) ≈ zeros(9) atol=1e-10
    end

    @testset "Test constraint Jacobian θ" begin
        function wrapper_constr(θ_ext)
            return copy(ms.constraint(opt, θ_ext))
        end

        function wrapper_jacobian(θ_ext)
            J = ms.jacobian(opt, θ_ext)
            return full(sparse(J.i, J.j, J.v))
        end

        function finite_difference_jacobian(θ_ext)
            return Calculus.jacobian(wrapper_constr, θ_ext, :central)
        end

        for θ in ([3.21, 3.2, 1.00],
                [3.21, 1.2, 1.54],
                [3.56, 3.07, 1.54])
            for x0_ext in (collect(linspace(0.2, 0.3, 10)),
                           collect(linspace(0.1, 0.8, 10)))
                θ_ext = [θ; x0_ext]
                finite_difference_jacobian(θ_ext)
                @test isapprox(wrapper_jacobian(θ_ext),
                               finite_difference_jacobian(θ_ext),
                               rtol=1e-5)
            end
        end
    end
end

@testset "Test extended vector conversions" begin
    srand(1)
    Nθ = 10
    Nx = 3
    M = 5
    θ = randn(Nθ)
    x0_list = Vector{Vector{Float64}}(M)
    for i = 1:M
        x0_list[i] = randn(Nx)
    end
    θ_extended = [θ; vcat(x0_list...)]
    @testset "Test build extended vector" begin
        θ_extended_buffer = zeros(θ_extended)
        ms.build_extended_vector!(θ_extended_buffer, θ, x0_list, Nθ, M, Nx)
        @test θ_extended_buffer ≈ θ_extended
    end

    @testset "Test read extended vector" begin
        θ_buffer = zeros(θ)
        x0_list_buffer = Vector{Vector{Float64}}(M)
        for i = 1:M
            x0_list_buffer[i] = zeros(Nx)
        end
        ms.read_extended_vector!(θ_buffer, x0_list_buffer, θ_extended, Nθ, M, Nx)
        @test θ_buffer ≈ θ
        @test x0_list_buffer ≈ x0_list
    end
end

@testset "Test jacobian indices" begin
    Nθ = 3; Nx = 2; M = 3
    i, j = ms.jacobian_indices(Nθ, M, Nx)
    @test i == [1, 2, 1, 2, 1, 2,
                3, 4, 3, 4, 3, 4,
                1, 2, 1, 2,
                3, 4, 3, 4,
                1, 2,
                3, 4]
    @test j == [1, 1, 2, 2, 3, 3,
                1, 1, 2, 2, 3, 3,
                4, 4, 5, 5,
                6, 6, 7, 7,
                6, 7,
                8, 9]
end

@testset "Test Jacobian from matrices" begin
    Nθ = 3; Nx = 2; M = 3
    i, j = ms.jacobian_indices(Nθ, M, Nx)
    nnz = length(i)
    v = Vector{Float64}(nnz)
    jacobian = ms.SparseMatrixCOO{Float64,Int64}(Nx*(M-1), Nθ + M*Nx, i, j, v)
    jacobian_θ = Vector{Matrix{Float64}}(M-1)
    for i = 1:M-1
        jacobian_θ[i] = ones(Nx, Nθ)
    end
    jacobian_x0 = Vector{Matrix{Float64}}(M-1)
    for i = 1:M-1
        jacobian_x0[i] = ones(Nx, Nx)
    end
    ms.jacobian_from_matrices!(jacobian, jacobian_θ, jacobian_x0,
                            Nθ, M, Nx; reset_constants=true)
    J = full(sparse(jacobian.i, jacobian.j, jacobian.v))
    @test J == Float64[1 1 1 1 1 -1  0  0  0;
                       1 1 1 1 1  0 -1  0  0;
                       1 1 1 0 0  1  1 -1  0;
                       1 1 1 0 0  1  1  0 -1]
end

@testset "Test single shooting" begin
    # From example 2.11, "Linear System Theory and design" Chen
    # With not very realistic values
    # Define linear function
    function linear_f(x_next, x, k)
        A = [-0.005    0   -0.1;
             0         0    0.1;
             1.0      -1.0    0]
        x_next .=  A*x
        return
    end
    # g(x) = x
    function linear_g(y, x, k)
      y .= dot([1, -1, 0], x)
      return
    end
    # Define initial conditions
    x0 = [5, 0.6, 0.5]
    # Initialize buffer
    N = 1000
    y = Vector{Vector{Float64}}(N)
    for i = 1:length(y)
       y[i] = zeros(1)
    end
    # Simulate
    x_final = ms.simulate_space_state!(y, linear_f, linear_g,
                                       x0, (1, N))

    # Define initial values
    k0_list = [1]
    x0_list = [[5, 0.6, 0.5]]
    θ0 = [-0.0047, 0.097, 0.097, 0.98]

    function f(y, dx, dθ, x, k, θ)
        A = [θ[1]  0 -θ[2];
             0     0  θ[3];
             θ[4] -θ[4]  0]
         y .= A*x
         dx .= A
         dθ .= [x[1] -x[3]    0      0;
                0      0      x[3]   0;
                0      0       0    x[1]-x[2]]
         return
    end
    function g(y, dx, dθ, x, k, θ)
      y .= dot([1, -1, 0], x)
      dx .= [1 -1 0]
      dθ .=  0
      return
    end

    opt = ms.OptimizationProblem(f, g, x0_list, y, k0_list, θ0)

    res = ms.solve(opt, options=Dict("gtol" => 1e-12,
                                     "xtol" => 1e-10))

    θext_opt = [-0.005, 0.1, 0.1, 1.0, 5, 0.6, 0.5]
    final_error = norm((res["x"] - θext_opt)./θext_opt,  Inf)
    initial_error = norm(([θ0; x0_list[1]] - θext_opt)./θext_opt, Inf)

    @test final_error / initial_error < 0.01
end

@testset "Test multiple shooting" begin
    M = 10
    function linear_f(x_next, x, k)
        A = [-0.005    0   -0.1;
             0         0    0.1;
             1.0      -1.0    0]
        x_next .=  A*x
        return
    end
    function linear_g(y, x, k)
        y .= dot([1, -1, 0], x)
        return
    end

    # Define initial conditions
    x0 = [5, 0.6, 0.5]
    x0_list = Vector{Vector{Float64}}(10)
    x0_list[1] = copy(x0)
    grid = collect(1:100:1001)
    k0_list = grid[1:end-1]
    # Initialize buffer
    y = Vector{Vector{Float64}}(1000)
    for i = 1:length(y)
        y[i] = zeros(1)
    end
    # Simulate
    x_final = copy(x0)
    for i = 1:length(grid)-1
        x_final = ms.simulate_space_state!(y[grid[i]:grid[i+1]-1],
                                          linear_f, linear_g,
                                          x0, (grid[i], grid[i+1]-1))
        if i != length(grid)-1
           x0 .= x_final
           x0_list[i+1] = copy(x0)
        end
    end

    θ0 = [-0.0049, 0.097, 0.097, 1.0]
    function f(y, dx, dθ, x, k, θ)
        A = [θ[1]  0 -θ[2];
             0     0  θ[3];
             θ[4] -θ[4]  0]
         y .= A*x
         dx .= A
         dθ .= [x[1] -x[3]    0      0;
                0      0      x[3]   0;
                0      0       0    x[1]-x[2]]
         return
    end
    function g(y, dx, dθ, x, k, θ)
      y .= dot([1, -1, 0], x)
      dx .= [1 -1 0]
      dθ .=  0
      return
    end

    opt = ms.OptimizationProblem(f, g, x0_list, y, k0_list, θ0)

    res = ms.solve(opt, options=Dict("gtol" => 1e-8,
                                     "xtol" => 1e-8,
                                     "maxiter" => 200))
    θ_opt = [-0.005, 0.1, 0.1, 1.0]
    final_error = norm((res["x"][1:4] - θ_opt)./θ_opt,  Inf)
    initial_error = norm((θ0 - θ_opt)./θ_opt, Inf)

    @test final_error / initial_error < 0.1
end


text:
function determineK(L₁₁,L₂₂,Lₘ)
    K = Lₘ/√(L₁₁*L₂₂)
end

function determineCancelledCenter(L₁₁,L₂₂)
    Lₘ = L₁₁
    Lₐ = Lₘ
    Lᵦ = L₂₂-Lₘ
    Lc = L₁₁-Lₘ
    K = Lₘ/√(L₁₁*L₂₂)
    return K,Lₐ,Lᵦ,Lc
end


text:
# The , as a decimal separator isn't yet supported by DataFrames
# We replace all , with . as a workaround

import LegacyStrings

dataset_path = joinpath(dirname(@__FILE__), "acute-inflammations-2.data.orig")

f = open(dataset_path)
data = LegacyStrings.utf16(read(f)[3:end])
close(f)

data = replace(data, ",", ".")

f = open("$dataset_path.custom", "w")
write(f, data)
close(f)

max_stars_repo_name: CoventryResearch/uci-data

max_issues_repo_name: CoventryResearch/uci-data

max_forks_repo_name: CoventryResearch/uci-data


text:
module ModeEstimation

using ..Turing
using Bijectors
using Random
using SciMLBase: OptimizationFunction, OptimizationProblem, AbstractADType

using DynamicPPL
using DynamicPPL: Model, AbstractContext, VarInfo, VarName,
    _getindex, getsym, getfield, settrans!,  setorder!,
    get_and_set_val!, istrans

export  constrained_space,
        MAP,
        MLE,
        OptimLogDensity,
        OptimizationContext,
        get_parameter_bounds,
        optim_objective,
        optim_function,
        optim_problem

struct constrained_space{x} end

struct MLE end
struct MAP end

"""
    OptimizationContext{C<:AbstractContext} <: AbstractContext

The `OptimizationContext` transforms variables to their constrained space, but
does not use the density with respect to the transformation. This context is
intended to allow an optimizer to sample in R^n freely.
"""
struct OptimizationContext{C<:AbstractContext} <: AbstractContext
    context::C
end

DynamicPPL.NodeTrait(::OptimizationContext) = DynamicPPL.IsParent()
DynamicPPL.childcontext(context::OptimizationContext) = context.context
DynamicPPL.setchildcontext(::OptimizationContext, child) = OptimizationContext(child)

# assume
function DynamicPPL.tilde_assume(rng::Random.AbstractRNG, ctx::OptimizationContext, spl, dist, vn, inds, vi)
    return DynamicPPL.tilde_assume(ctx, spl, dist, vn, inds, vi)
end

function DynamicPPL.tilde_assume(ctx::OptimizationContext{<:LikelihoodContext}, spl, dist, vn, inds, vi)
    r = vi[vn]
    return r, 0
end

function DynamicPPL.tilde_assume(ctx::OptimizationContext, spl, dist, vn, inds, vi)
    r = vi[vn]
    return r, Distributions.logpdf(dist, r)
end

# dot assume
function DynamicPPL.dot_tilde_assume(rng::Random.AbstractRNG, ctx::OptimizationContext, sampler, right, left, vns, inds, vi)
    return DynamicPPL.dot_tilde_assume(ctx, sampler, right, left, vns, inds, vi)
end

function DynamicPPL.dot_tilde_assume(ctx::OptimizationContext{<:LikelihoodContext}, sampler::SampleFromPrior, right, left, vns, _, vi)
    # Values should be set and we're using `SampleFromPrior`, hence the `rng` argument shouldn't
    # affect anything.
    r = DynamicPPL.get_and_set_val!(Random.GLOBAL_RNG, vi, vns, right, sampler)
    return r, 0
end

function DynamicPPL.dot_tilde_assume(ctx::OptimizationContext, sampler::SampleFromPrior, right, left, vns, _, vi)
    # Values should be set and we're using `SampleFromPrior`, hence the `rng` argument shouldn't
    # affect anything.
    r = DynamicPPL.get_and_set_val!(Random.GLOBAL_RNG, vi, vns, right, sampler)
    return r, loglikelihood(right, r)
end

"""
    OptimLogDensity{M<:Model,C<:Context,V<:VarInfo}

A struct that stores the log density function of a `DynamicPPL` model.
"""
struct OptimLogDensity{M<:Model,C<:AbstractContext,V<:VarInfo}
    "A `DynamicPPL.Model` constructed either with the `@model` macro or manually."
    model::M
    "A `DynamicPPL.AbstractContext` used to evaluate the model. `LikelihoodContext` or `DefaultContext` are typical for MAP/MLE."
    context::C
    "A `DynamicPPL.VarInfo` struct that will be used to update model parameters."
    vi::V
end

"""
    OptimLogDensity(model::Model, context::AbstractContext)

Create a callable `OptimLogDensity` struct that evaluates a model using the given `context`.
"""
function OptimLogDensity(model::Model, context::AbstractContext)
    init = VarInfo(model)
    return OptimLogDensity(model, context, init)
end

"""
    (f::OptimLogDensity)(z)

Evaluate the log joint (with `DefaultContext`) or log likelihood (with `LikelihoodContext`)
at the array `z`.
"""
function (f::OptimLogDensity)(z)
    spl = DynamicPPL.SampleFromPrior()

    varinfo = DynamicPPL.VarInfo(f.vi, spl, z)
    f.model(varinfo, spl, f.context)
    return -DynamicPPL.getlogp(varinfo)
end

function (f::OptimLogDensity)(F, G, H, z)
    # Throw an error if a second order method was used.
    if H !== nothing
        error("Second order optimization is not yet supported.")
    end

    spl = DynamicPPL.SampleFromPrior()

    if G !== nothing
        # Calculate log joint and the gradient
        l, g = Turing.gradient_logp(
            z,
            DynamicPPL.VarInfo(f.vi, spl, z),
            f.model,
            spl,
            f.context
        )

        # Use the negative gradient because we are minimizing.
        G[:] = -g

        # If F is something, return that since we already have the
        # log joint.
        if F !== nothing
            F = -l
            return F
        end
    end

    # No gradient necessary, just return the log joint.
    if F !== nothing
        F = f(z)
        return F
    end

    return nothing
end

#################################################
# Generic optimisation objective initialisation #
#################################################

function transform!(f::OptimLogDensity)
  spl = DynamicPPL.SampleFromPrior()

  ## Check link status of vi in OptimLogDensity
  linked = DynamicPPL.islinked(f.vi, spl)

  ## transform into constrained or unconstrained space depending on current state of vi
  if !linked
    DynamicPPL.link!(f.vi, spl)
  else
    DynamicPPL.invlink!(f.vi, spl)
  end

  return nothing
end

function transform!(p::AbstractArray, vi::DynamicPPL.VarInfo, ::constrained_space{true})
    spl = DynamicPPL.SampleFromPrior()

    linked = DynamicPPL.islinked(vi, spl)

    # !linked && DynamicPPL.link!(vi, spl)
    !linked && return identity(p)
    vi[spl] = p
    DynamicPPL.invlink!(vi,spl)
    p .= vi[spl]

    linked && DynamicPPL.link!(vi,spl)

    return nothing
end

function transform!(p::AbstractArray, vi::DynamicPPL.VarInfo, ::constrained_space{false})
  spl = DynamicPPL.SampleFromPrior()

  linked = DynamicPPL.islinked(vi, spl)
  linked && DynamicPPL.invlink!(vi, spl)
  vi[spl] = p
  DynamicPPL.link!(vi, spl)
  p .= vi[spl]
  !linked && DynamicPPL.invlink!(vi, spl)

  return nothing
end

function transform(p::AbstractArray, vi::DynamicPPL.VarInfo, con::constrained_space)
    tp = copy(p)
    transform!(tp, vi, con)
  return tp
end

abstract type AbstractTransform end

struct ParameterTransform{T<:DynamicPPL.VarInfo, S<:constrained_space} <: AbstractTransform
    vi::T
    space::S
end

struct Init{T<:DynamicPPL.VarInfo, S<:constrained_space} <: AbstractTransform
    vi::T
    space::S
end

function (t::AbstractTransform)(p::AbstractArray)
    return transform(p, t.vi, t.space)
end

function (t::Init)()
    return t.vi[DynamicPPL.SampleFromPrior()]
end

 function get_parameter_bounds(model::DynamicPPL.Model)
    vi = DynamicPPL.VarInfo(model)
    spl = DynamicPPL.SampleFromPrior()

    ## Check link status of vi
    linked = DynamicPPL.islinked(vi, spl)

    ## transform into unconstrained
    !linked && DynamicPPL.link!(vi, spl)

    lb = transform(fill(-Inf,length(vi[DynamicPPL.SampleFromPrior()])), vi, constrained_space{true}())
    ub = transform(fill(Inf,length(vi[DynamicPPL.SampleFromPrior()])), vi, constrained_space{true}())

    return lb, ub
end

function _optim_objective(model::DynamicPPL.Model, ::MAP, ::constrained_space{false})
  ctx = OptimizationContext(DynamicPPL.DefaultContext())
  obj = OptimLogDensity(model, ctx)

  transform!(obj)
  init = Init(obj.vi, constrained_space{false}())
  t = ParameterTransform(obj.vi, constrained_space{true}())

  return (obj=obj, init = init, transform=t)
end

function _optim_objective(model::DynamicPPL.Model, ::MAP, ::constrained_space{true})
    ctx = OptimizationContext(DynamicPPL.DefaultContext())
    obj = OptimLogDensity(model, ctx)

    init = Init(obj.vi, constrained_space{true}())
    t = ParameterTransform(obj.vi, constrained_space{true}())

    return (obj=obj, init = init, transform=t)
  end

function _optim_objective(model::DynamicPPL.Model, ::MLE,  ::constrained_space{false})
    ctx = OptimizationContext(DynamicPPL.LikelihoodContext())
    obj = OptimLogDensity(model, ctx)

    transform!(obj)
    init = Init(obj.vi, constrained_space{false}())
    t = ParameterTransform(obj.vi, constrained_space{true}())

    return (obj=obj, init = init, transform=t)
end

function _optim_objective(model::DynamicPPL.Model, ::MLE, ::constrained_space{true})
    ctx = OptimizationContext(DynamicPPL.LikelihoodContext())
    obj = OptimLogDensity(model, ctx)

    init = Init(obj.vi, constrained_space{true}())
    t = ParameterTransform(obj.vi, constrained_space{true}())

    return (obj=obj, init = init, transform=t)
end

function optim_objective(model::DynamicPPL.Model, estimator::Union{MLE, MAP}; constrained::Bool=true)
    return _optim_objective(model, estimator, constrained_space{constrained}())
end

function optim_function(model::DynamicPPL.Model, estimator::Union{MLE, MAP}; constrained::Bool=true, autoad::Union{Nothing, AbstractADType}=nothing)
    obj, init, t = optim_objective(model, estimator; constrained=constrained)

    l(x,p) = obj(x)
    f = isa(autoad, AbstractADType) ? OptimizationFunction(l, autoad) : OptimizationFunction(l; grad = (G,x,p) -> obj(nothing, G, nothing, x), hess = (H,x,p) -> obj(nothing, nothing, H, x))

    return (func=f, init=init, transform = t)
end

function optim_problem(model::DynamicPPL.Model, estimator::Union{MAP, MLE}; constrained::Bool=true, init_theta=nothing, autoad::Union{Nothing, AbstractADType}=nothing, kwargs...)
    f = optim_function(model, estimator; constrained=constrained, autoad=autoad)

    init_theta = init_theta === nothing ? f.init() : f.init(init_theta)

    prob = OptimizationProblem(f.func, init_theta, nothing; kwargs...)

    return (prob=prob, init=f.init, transform = f.transform)
end

end


text:
abstract type API end

mutable struct OpenFigiAPI <: API
    protocol::AbstractString
    base::AbstractString
    version::AbstractString
    path::AbstractString
    headers::Vector{Pair{AbstractString, AbstractString}}
end

mutable struct OpenFigiAsset
    figi::Union{String, Nothing}
    marketSector::Union{String, Nothing}
    securityType::Union{String, Nothing}
    ticker::Union{String, Nothing}
    name::Union{String, Nothing}
    exchCode::Union{String, Nothing}
    securityDescription::Union{String, Nothing}
    securityType2::Union{String, Nothing}
    compositeFIGI::Union{String, Nothing}
    shareClassFIGI::Union{String, Nothing}
end

"""
    OpenFigiAPI(protocol::AbstractString="https",
                base::AbstractString="api.openfigi.com",
                version::AbstractString="v3",
                path::AbstractString="mapping";
                headers::Vector{Pair{String, String}}=["Content-Type"=>"application/json"],
                apikey::String="")::OpenFigiAPI

Create OpenFigiAPI.

See also: [`fetchsecuritydata`](@ref)

# Example
```jldoctest; setup = :(using FinancialSymbology)
julia> OpenFigiAPI()
OpenFigiAPI: https://api.openfigi.com/v3/mapping
```
"""
function OpenFigiAPI(protocol::AbstractString="https",
                     base::AbstractString="api.openfigi.com",
                     version::AbstractString="v3",
                     path::AbstractString="mapping";
                     headers::Vector{Pair{String, String}}=["Content-Type"=>"application/json"],
                     apikey::String="")::OpenFigiAPI

    apikey::String = haskey(ENV, APIKEYNAME) && length(apikey) == 0 ? ENV[APIKEYNAME] : apikey
    push!(headers, APIKEYNAME=>apikey)

    return OpenFigiAPI(protocol, base, version, path, headers)

end

"""
    OpenFigiAsset(; figi::Union{String, Nothing}=nothing,
                  marketSector::Union{String, Nothing}=nothing,
                  securityType::Union{String, Nothing}=nothing,
                  ticker::Union{String, Nothing}=nothing,
                  name::Union{String, Nothing}=nothing,
                  exchCode::Union{String, Nothing}=nothing,
                  securityDescription::Union{String, Nothing}=nothing,
                  securityType2::Union{String, Nothing}=nothing,
                  compositeFIGI::Union{String, Nothing}=nothing,
                  shareClassFIGI::Union{String, Nothing}=nothing)::OpenFigiAsset

Each element of the `StructArray` returned by [`fetchsecuritydata`](@ref) for an
individual [`Identifier`](@ref identifier_header) is an `OpenFigiAsset`.

# Example
```jldoctest; setup = :(using FinancialSymbology)
julia> aapl = fetchsecuritydata(Ticker("AAPL US Equity"));

julia> aapl[1]
FIGI: BBG000B9XRY4 Common Stock

julia> aapl[1].shareClassFIGI
"BBG001S5N8V8"
```
"""
function OpenFigiAsset(; figi::Union{String, Nothing}=nothing,
                       marketSector::Union{String, Nothing}=nothing,
                       securityType::Union{String, Nothing}=nothing,
                       ticker::Union{String, Nothing}=nothing,
                       name::Union{String, Nothing}=nothing,
                       exchCode::Union{String, Nothing}=nothing,
                       securityDescription::Union{String, Nothing}=nothing,
                       securityType2::Union{String, Nothing}=nothing,
                       compositeFIGI::Union{String, Nothing}=nothing,
                       shareClassFIGI::Union{String, Nothing}=nothing)::OpenFigiAsset

    return OpenFigiAsset(figi, marketSector, securityType, ticker, name, exchCode,
                         securityDescription, securityType2, compositeFIGI, shareClassFIGI)

end

keystosymbols(d::Dict{String, T} where T)::Dict{Symbol, T} where T = Dict([Symbol(k)=>v for (k, v) in d])

function OpenFigiAsset(d::Dict{String, Any})::OpenFigiAsset
    return OpenFigiAsset(; keystosymbols(d)...)
end

max_stars_repo_name: robfs/figi

max_issues_repo_name: robfs/figi

max_forks_repo_name: robfs/figi


text:
using ServiceSolicitation
# using Aqua
using Test

@testset "ServiceSolicitation.jl" begin
    # Aqua.test_all(ServiceSolicitation)

    function rangemap!(f::F, allargs, start, stop) where {F}
        dest = first(allargs)
        args = Base.tail(allargs)
        @inbounds @simd for i ∈ start:stop
            dest[i] = f(Base.unsafe_getindex.(args, i)...)
        end
        nothing
    end

    function tmap!(f::F, args::Vararg{AbstractArray,K}) where {K,F}
        dest = first(args)
        N = length(dest)
        mapfun! = (allargs, start, stop) -> rangemap!(f, allargs, start, stop)
        batch(mapfun!, (N, num_threads()), args...)
        dest
    end

    x = rand(1024); y = rand(length(x)); z = similar(x);
    foo(x,y) = exp(-0.5abs2(x-y))
    @test tmap!(foo, z, x, y) ≈ foo.(x, y)

end


text:
isdefined(Base, :__precompile__) && __precompile__()

module WoodburyMatrices

using Compat

import Base: *, \, A_ldiv_B!, convert, copy, det, full, show, similar, size

export Woodbury, SymWoodbury, liftFactor

#### Woodbury matrices ####
"""
`W = Woodbury(A, U, C, V)` creates a matrix `W` identical to `A + U*C*V` whose inverse will be calculated using
the Woodbury matrix identity.
"""
type Woodbury{T,AType,UType,VType,CType,CpType} <: AbstractMatrix{T}
    A::AType
    U::UType
    C::CType
    Cp::CpType
    V::VType
    tmpN1::Vector{T}
    tmpN2::Vector{T}
    tmpk1::Vector{T}
    tmpk2::Vector{T}
end

function Woodbury{T}(A, U::AbstractMatrix{T}, C, V::AbstractMatrix{T})
    N = size(A, 1)
    k = size(U, 2)
    if size(A, 2) != N || size(U, 1) != N || size(V, 1) != k || size(V, 2) != N
        throw(DimensionMismatch("Sizes of U ($(size(U))) and/or V ($(size(V))) are inconsistent with A ($(size(A)))"))
    end
    if k > 1
        if size(C, 1) != k || size(C, 2) != k
            throw(DimensionMismatch("C should be $(k)x$(k)"))
        end
    end
    Cp = inv(inv(C) + V*(A\U))
    # temporary space for allocation-free solver
    tmpN1 = Array(T, N)
    tmpN2 = Array(T, N)
    tmpk1 = Array(T, k)
    tmpk2 = Array(T, k)
    # don't copy A, it could be huge
    Woodbury{T,typeof(A),typeof(U),typeof(V),typeof(C),typeof(Cp)}(A, copy(U), copy(C), Cp, copy(V), tmpN1, tmpN2, tmpk1, tmpk2)
end

Woodbury{T}(A, U::Vector{T}, C, V::Matrix{T}) = Woodbury(A, reshape(U, length(U), 1), C, V)

size(W::Woodbury) = size(W.A)

function show(io::IO, W::Woodbury)
    println(io, summary(W), ":")
    print(io, "A:\n", W.A)
    print(io, "\nU:\n")
    Base.print_matrix(io, W.U)
    if isa(W.C, Matrix)
        print(io, "\nC:\n")
        Base.print_matrix(io, W.C)
    else
        print(io, "\nC: ", W.C)
    end
    print(io, "\nV:\n")
    Base.print_matrix(io, W.V)
end

full{T}(W::Woodbury{T}) = convert(Matrix{T}, W)
convert{T}(::Type{Matrix{T}}, W::Woodbury{T}) = full(W.A) + W.U*W.C*W.V

copy(W::Woodbury) = Woodbury(W.A, W.U, W.C, W.V)

## Woodbury matrix routines ##

*(W::Woodbury, B::StridedVecOrMat)=W.A*B + W.U*(W.C*(W.V*B))

function \(W::Woodbury, R::StridedVecOrMat)
    AinvR = W.A\R
    AinvR - W.A\(W.U*(W.Cp*(W.V*AinvR)))
end

det(W::Woodbury)=det(W.A)*det(W.C)/det(W.Cp)

function A_ldiv_B!(W::Woodbury, B::AbstractVector)
    length(B) == size(W, 1) || throw(DimensionMismatch("Vector length $(length(B)) must match matrix size $(size(W,1))"))
    copy!(W.tmpN1, B)
    Alu = lufact(W.A) # Note. This makes an allocation (unless A::LU). Alternative is to destroy W.A.
    A_ldiv_B!(Alu, W.tmpN1)
    A_mul_B!(W.tmpk1, W.V, W.tmpN1)
    A_mul_B!(W.tmpk2, W.Cp, W.tmpk1)
    A_mul_B!(W.tmpN2, W.U, W.tmpk2)
    A_ldiv_B!(Alu, W.tmpN2)
    for i = 1:length(W.tmpN2)
        @inbounds B[i] = W.tmpN1[i] - W.tmpN2[i]
    end
    B
end

"""
Given a type `T`, an integer `n` and `m` tuples of the form (i, j, v), build
sparse matrices `rows`, `vals`, `cols` such that the product
`out = rows * vals * cols` is equivalent to:

```julia
out = zeros(T, n, n)

for (i, j, v) in args
    out[i, j] = v
end
```

The first two components (`i` and `j`) of each tuple should be integers
whereas the third component should be of type `T`

Example:

```
julia> r, v, c = WoodburyMatrices.sparse_factors(Float64, 3,
                                                 (1, 1, 2.0),
                                                 (2, 2, 3.0),
                                                 (3, 3, 4.0));

julia> r*c*v - Diagonal([2.0, 3.0, 4.0])
3x3 sparse matrix with 0 Float64 entries:
```

"""
function sparse_factors{T}(::Type{T}, n::Int, args::@compat(Tuple{Int, Int, Any})...)
    m = length(args)
    rows = spzeros(T, n, m)
    cols = spzeros(T, m, n)
    vals = zeros(T, m, m)

    ix = 1
    for (i, (row, col, val)) in enumerate(args)
        rows[row, ix] = 1
        cols[ix, col] = 1
        vals[ix, ix] = val
        ix += 1
    end

    rows, vals, cols
end

include("SymWoodburyMatrices.jl")

end


text:
using Catlab, Catlab.CategoricalAlgebra, Catlab.Graphics, Catlab.Programs
using CombinatorialSpaces.ExteriorCalculus

@present DiffusionSpace2D(FreeExtCalc2D) begin
  X::Space
  k::Hom(Form1(X), Form1(X)) # diffusivity of space, usually constant (scalar multiplication)
  L::Hom(Form0(X)⊗Form0(X), Form0(X))
end

FicksLaw2D = @free_diagram DiffusionSpace2D begin
  C::Form0{X}
  dC::Form1{X}
  J::DualForm1{X}

  dC == d₀{X}(C)
  J == ⋆₁{X}(k(dC))
end

DiffusionConservation2D = @free_diagram DiffusionSpace2D begin
  (C, Ċ)::Form0{X}
  ϕ::DualForm1{X}
  dϕ::DualForm2{X}

  dϕ == dual_d₁{X}(ϕ)
  Ċ == ∂ₜ{Form0{X}}(C)
  Ċ == ⋆₀⁻¹{X}(dϕ)
end;

using Decapodes

compose_diffusion = @relation (C, Ċ, ϕ) begin
  ficks_law(C, ϕ)
  mass_conservation(C, Ċ, ϕ)
end

Decapodes.OpenDiagrams.draw(compose_diffusion)

using Decapodes.OpenDiagrams

composed_diffusion = oapply(compose_diffusion, [
    OpenDiagram(FicksLaw2D, [:C, :J]),
    OpenDiagram(DiffusionConservation2D, [:C, :Ċ, :ϕ]),
]);

res = diag2dwd(composed_diffusion.functor)


text:
using SimpleCG

plane_x = Plane(Vector3(0, 1, 0), 0)
plane_y = Plane(Vector3(0, 0, 1), -50)
plane_z = Plane(Vector3(1, 0, 0), -20)

light1 = PointLight(White * 1000, Vector3(30, 40, 20))
light2 = SpotLight(Red * 3000, Vector3(0, 30, 10), Vector3(0, -1, -1), π/9, π/6, 1)
light3 = SpotLight(Green * 3000, Vector3(6, 30, 20), Vector3(0, -1, -1), π/9, π/6, 1)
light4 = SpotLight(Blue * 3000, Vector3(-6, 30, 20), Vector3(0, -1, -1), π/9, π/6, 1)

canvas = Canvas(800, 800);
camera = PerspectiveCamera(Vector3(0, 40, 15), Vector3(0, -1.25, -1), Vector3(0, 1, 0), π/3)
scene = UnionGeometry([plane_x, plane_y, plane_z])
SimpleCG.inner_render!(canvas, scene, camera, SimpleCG.light_render!,
    Lights([light1, light2, light3, light4]))
imshow(canvas)


text:
using JLD2

# Each file (big-endian, 32 bits) contains 12 records of 128 x 60 values

Nφ = 128
Nλ = 60

function convert_lat_lon_bin(filepath, Nmonths)

    lat_lon_data = Array{Float32}(undef, Nλ*Nφ*Nmonths)

    read!(filepath, lat_lon_data)

    lat_lon_data = bswap.(lat_lon_data) |> Array{Float64}
    if Nmonths == 1
        lat_lon_var  = reshape(lat_lon_data, Nφ, Nλ)
    else
        lat_lon_var  = reshape(lat_lon_data, Nφ, Nλ, Nmonths)
    end

    return lat_lon_var
end

files = [("bathymetry_lat_lon_128x60",       "bathymetry",  1),
         ("sea_surface_temperature_25_128x60x12", "sst25", 12),
         ("tau_x_128x60x12",                      "tau_x", 12),
         ("tau_y_128x60x12",                      "tau_y", 12)]

for (filename, varname, Nmonths) in files

    var_data = convert_lat_lon_bin("$(filename)_FP32.bin", Nmonths)

    jldopen("$filename.jld2", "w") do file
        if Nmonths == 1
            file[varname] = var_data
        else
            for month in 1:Nmonths
                file["$varname/$month"] = var_data[:,:,month]
            end
        end
    end
end


text:
using LinearAlgebra
using DoubleFloats
using SemiLagrangian:
    AbstractInterpolation,
    Lagrange,
    B_SplineLU,
    B_SplineFFT,
    interpolate!,
    isbspline,
    getbspline,
    getprecal,
    get_allprecal,
    get_order,
    EdgeType,
    InsideEdge,
    CircEdge,
    interpolatemod!,
    getinverse,
    OpTuple

function test_interp(
    interp::AbstractInterpolation{Rational{BigInt},edge},
    dec,
    sz,
) where {edge}

    @time @testset "test interpolation  $interp dec=$dec" begin
        fct = if edge == CircEdge
            order = 3
            getbspline(big(order), 0)
        else
            order = 0
            x -> x^3 - x^2 - x // 6 + 1 // 4
        end
        mesh = collect((order + 1) * big.(0:(sz-1))) // sz
        deb = fct.(mesh)
        #       println("deb=$deb")
        ref = if edge == CircEdge
            decf = mod(dec, 1 // 1)
            decint = div(dec, 1 // 1)
            circshift(fct.(mesh .+ (order + 1) * decf // sz), (-decint,))
        else
            fct.(mesh .+ (order + 1) * dec // sz)
        end
        res = zeros(Rational{BigInt}, sz)
        interpolate!(res, deb, dec, interp)
        if edge == CircEdge
            res2 = zeros(Rational{BigInt}, sz)
            decint = convert(Int, floor(dec))
            decfloat = dec - decint
            precal = getprecal(interp, decfloat)
            # println("outside : precal=$precal decint=$decint decfloat=$decfloat")
            # println("order=$(get_order(interp)) order from precal=$(size(precal,1)-1)")
            interpolate!(res2, deb, decint, precal, interp)
            println("norm=$(float(norm(res2-ref,Inf)))")
            if isbspline(interp)
                @test res2 == ref
            else
                @test size(filter(x -> x == 0, res2 - ref), 1) == 120
            end
        end
        println("norm=$(float(norm(res-ref,Inf)))")
        diff = res - ref
        # for i=1:sz
        #     println("i=$i diff=$(convert(Float64,diff[i])), $(diff[i])")
        # end
        fl = edge != CircEdge || isbspline(interp)
        println("fl=$fl isbspline=$(isbspline(interp))")
        if fl
            @test res == ref
        else
            @test size(filter(x -> x == 0, res - ref), 1) == 120
        end

    end
end
function test_interp(interp, sz)
    for i = 0:5
        println("test_interp i=$i")
        test_interp(interp, big"1235" // 10240 + i, sz)
    end
end

function test_inv0(
    t_interp::Vector{I},
    dec::NTuple{N,T},
    sz::NTuple{N,Int},
) where {T,N,I<:AbstractInterpolation{T,CircEdge}}
    res = ntuple(x -> zeros(T, sz), N)
    res2 = ntuple(x -> zeros(T, sz), N)
    for ind in CartesianIndices(sz)
        for i = 1:N
            res[i][ind] = T(ind.I[i] - 1)
        end
    end
    decbegin = ntuple(x -> dec[x] < 0 ? -Int(floor(dec[x])) : 0, N)
    decend = ntuple(x -> dec[x] > 0 ? Int(ceil(dec[x])) : 0, N)
    for i = 1:N
        interpolatemod!(res2[i], res[i], ind -> dec, t_interp, T(sz[i]), decbegin, decend)
    end

    resinv = ntuple(x -> zeros(T, sz), N)
    resinv2 = ntuple(x -> zeros(T, sz), N)

    calinverse!(resinv, res2)
    calinverse!(resinv2, resinv)

    note = resinv .- res

    note[1] .= mod.(note[1] .+ (sz[1] / 2 + dec[1]), sz[1]) .- sz[1] / 2
    note[2] .= mod.(note[2] .+ (sz[2] / 2 + dec[2]), sz[2]) .- sz[2] / 2

    @show norm.(note)

    for i = 1:N
        @test norm(note[i]) < 1e-12
        @test norm(res2[i] - resinv2[i]) < 1e-12
    end

end

function test_interp2d(
    t_interp::Vector{I},
    coeff::T,
    sz::Tuple{Int,Int},
) where {T,N,I<:AbstractInterpolation{T,CircEdge}}
    prec = 1000 * eps(T)
    @show prec
    ref = zeros(T, sz)
    ref2 = zeros(T, sz)
    cmplxref = zeros(Complex{T}, sz)
    opref = zeros(OpTuple{2,T}, sz)
    bufdec = zeros(Complex{T}, sz)
    opbufdec = zeros(OpTuple{2,T}, sz)
    res1 = zeros(T, sz)
    res2 = zeros(T, sz)
    res3 = zeros(T, sz)
    res4 = zeros(T, sz)
    opres1 = zeros(OpTuple{2,T}, sz)
    cmplxres2 = zeros(Complex{T}, sz)
    cmplxres3 = zeros(Complex{T}, sz)
    cmplxres4 = zeros(Complex{T}, sz)

    for ind in CartesianIndices(sz)
        x = 2T(pi) * sum(ind.I ./ sz)
        bufdec[ind] = coeff * (cos(x + 1) + im * cos(x + 2))
        opbufdec[ind] = coeff * OpTuple((cos(x + 1), cos(x + 2)))
        ref[ind] = sin(x + 3) + cos(x - 1)
        ref2[ind] = sin(x + 1) + 3cos(-x + 2) / 5
        cmplxref[ind] = coeff * (sin(x + 4) + im * cos(x + 5))
        opref[ind] = coeff * OpTuple((sin(x + 4), cos(x + 5)))
    end

    @time interpolate!(res1, ref, opbufdec, t_interp)
    @time interpolate!(res2, ref, bufdec, t_interp)
    @time interpolate!(res3, ref, ind -> reim(bufdec[ind]), t_interp)
    @time interpolate!(res4, ref2, bufdec, t_interp)

    norm2 = norm(res2 - res3)
    norm1 = norm(res1 - res3)

    @test norm1 < prec
    @test norm2 < prec
    @show norm1, norm2

    @time interpolate!(opres1, opref, opbufdec, t_interp)
    @time interpolate!(cmplxres2, cmplxref, bufdec, t_interp)
    @time interpolate!(cmplxres3, cmplxref, ind -> reim(bufdec[ind]), t_interp)
    @time interpolate!(cmplxres4, ref + im * ref2, bufdec, t_interp)

    optocmplx(v::OpTuple{2,T}) where {T} = v.v[1] + im * v.v[2]

    norm1 = norm(optocmplx.(opres1) - cmplxres3)
    norm2 = norm(cmplxres2 - cmplxres3)

    norm4 = norm(res2 + im * res4 - cmplxres4)

    @test norm1 < prec
    @test norm2 < prec
    @test norm4 < prec
    @show norm1, norm2, norm4

end

function test_inv(
    t_interp::Vector{I},
    coeff::T,
    sz::NTuple{N,Int},
) where {T,N,I<:AbstractInterpolation{T,CircEdge}}
    res = ntuple(x -> zeros(T, sz), N)
    res2 = ntuple(x -> zeros(T, sz), N)
    res3 = ntuple(x -> zeros(T, sz), N)
    dec = ntuple(x -> zeros(T, sz), N)
    resnew = ntuple(x -> zeros(T, sz), N)
    for ind in CartesianIndices(sz)
        for i = 1:N
            dec[i][ind] = coeff * cos(i + 2T(pi) * sum(ind.I ./ sz))
            res[i][ind] = T(ind.I[i] - 1)
        end
    end
    extr = extrema.(dec)
    decbegin = ntuple(x -> -Int(floor(extr[x][1])), N)
    decend = ntuple(x -> Int(ceil(extr[x][2])), N)
    for i = 1:N
        interpolatemod!(
            res2[i],
            res[i],
            ind -> ntuple(x -> dec[x][ind], N),
            t_interp,
            T(sz[i]),
            decbegin,
            decend,
        )
    end
    for i = 1:N
        interpolatemod!(
            res3[i],
            res2[i],
            ind -> ntuple(x -> dec[x][ind], N),
            t_interp,
            T(sz[i]),
            decbegin,
            decend,
        )
    end

    resinv = ntuple(x -> zeros(T, sz), N)
    resinv2 = ntuple(x -> zeros(T, sz), N)

    calinverse!(resinv, res2)

    calinverse!(resinv2, resinv)

    @show Float32.(res[1]) + im * Float32.(res[2])
    @show Float32.(res2[1]) + im * Float32.(res2[2])
    @show Float32.(resinv2[1]) + im * Float32.(resinv2[2])
    @show Float32.(dec[1]) + im * Float32.(dec[2])
    @show Float32.(resinv[1]) + im * Float32.(resinv[2])
    @show Float32.(resinv[1][CartesianIndex(10, 10):CartesianIndex(12, 12)]) +
          im * Float32.(resinv[2][CartesianIndex(10, 10):CartesianIndex(12, 12)])

    resfmr =
        ntuple(x -> mod.(res2[x] - res[x] - dec[x] .+ sz[x] / 2, sz[x]) .- sz[x] / 2, 2)
    resfmr2 =
        ntuple(x -> mod.(res3[x] - res2[x] - dec[x] .+ sz[x] / 2, sz[x]) .- sz[x] / 2, 2)
    @show norm(resfmr)
    @show norm(resfmr2)

    decnew = ntuple(x -> resinv[x] - res[x], 2)
    extr = extrema.(decnew)
    decbegin = ntuple(x -> -Int(floor(extr[x][1])), N)
    decend = ntuple(x -> Int(ceil(extr[x][2])), N)
    for i = 1:N
        interpolatemod!(
            resnew[i],
            res2[i],
            ind -> ntuple(x -> decnew[x][ind], N),
            t_interp,
            T(sz[i]),
            decbegin,
            decend,
        )
    end
    @show resnew

    for i = 1:N
        @test norm(res2[i] - resinv2[i]) < 1e-12
    end

end
function test_invfmr(
    t_interp::Vector{I},
    coeff::T,
    sz::NTuple{N,Int},
) where {T,N,I<:AbstractInterpolation{T,CircEdge}}
    res = ntuple(x -> zeros(T, sz), N)
    res2 = ntuple(x -> zeros(T, sz), N)
    res3 = ntuple(x -> zeros(T, sz), N)
    dec = ntuple(x -> zeros(T, sz), N)
    dec2 = ntuple(x -> zeros(T, sz), N)
    resnew = ntuple(x -> zeros(T, sz), N)
    for ind in CartesianIndices(sz)
        for i = 1:N
            dec[i][ind] = coeff * cos(i + 2T(pi) * sum(ind.I ./ sz))
            res[i][ind] = T(ind.I[i] - 1)
        end
    end
    extr = extrema.(dec)
    decbegin = ntuple(x -> -Int(floor(extr[x][1])), N)
    decend = ntuple(x -> Int(ceil(extr[x][2])), N)
    for i = 1:N
        interpolatemod!(
            res2[i],
            res[i],
            ind -> ntuple(x -> dec[x][ind], N),
            t_interp,
            T(sz[i]),
            decbegin,
            decend,
        )
    end
    @show res
    @show res2
    extr = extrema.(ntuple(x -> -dec[x], N))
    decbegin = ntuple(x -> -Int(floor(extr[x][1])), N)
    decend = ntuple(x -> Int(ceil(extr[x][2])), N)
    for i = 1:N
        interpolatemod!(
            dec2[i],
            -dec[i],
            ind -> ntuple(x -> dec[x][ind], N),
            t_interp,
            T(sz[i]),
            decbegin,
            decend,
        )
    end
    for i = 1:N
        dec2[i] .= mod.(dec2[i] .+ sz[i] / 2, sz[i]) .- sz[i] / 2
    end
    @show dec
    @show dec2
    extr = extrema.(dec2)
    decbegin = ntuple(x -> -Int(floor(extr[x][1])), N)
    decend = ntuple(x -> Int(ceil(extr[x][2])), N)

    for i = 1:N
        interpolatemod!(
            res3[i],
            res2[i],
            ind -> ntuple(x -> dec2[x][ind], N),
            t_interp,
            T(sz[i]),
            decbegin,
            decend,
        )
    end

    @show res2
    @show res3

    res4 = ntuple(x -> res[x] - res3[x], N)
    for i = 1:N
        res4[i] .= mod.(res4[i] .+ sz[i] / 2, sz[i]) .- sz[i] / 2
    end

    @show res4

    @show norm(res4)
    res5 = ntuple(x -> zeros(T, sz), N)
    ind = 1
    while norm(res4) > eps(Double64) * length(res4[1])
        if norm(res4) > 1
            res4 = ntuple(x -> res4[x] / 2, N)
        end
        extr = extrema.(res4)
        decbegin = ntuple(x -> -Int(floor(extr[x][1])), N)
        decend = ntuple(x -> Int(ceil(extr[x][2])), N)
        for i = 1:N
            interpolatemod!(
                res5[i],
                res3[i],
                ind -> ntuple(x -> res4[x][ind], N),
                t_interp,
                T(sz[i]),
                decbegin,
                decend,
            )
        end
        res4 = ntuple(x -> res[x] - res5[x], N)
        for i = 1:N
            res4[i] .= mod.(res4[i] .+ sz[i] / 2, sz[i]) .- sz[i] / 2
            res3[i] .= res5[i]
        end

        @show res4

        @show ind, norm(res4)

        ind += 1

    end

end
function test_getinv(
    t_interp::Vector{I},
    coeff::T,
    sz::NTuple{N,Int},
) where {T,N,I<:AbstractInterpolation{T,CircEdge}}
    res = ntuple(x -> zeros(T, sz), N)
    res2 = ntuple(x -> zeros(T, sz), N)
    res3 = ntuple(x -> zeros(T, sz), N)
    dec = ntuple(x -> zeros(T, sz), N)
    dec2 = ntuple(x -> zeros(T, sz), N)
    resnew = ntuple(x -> zeros(T, sz), N)
    for ind in CartesianIndices(sz)
        for i = 1:N
            dec[i][ind] = coeff * cos(i + 2T(pi) * sum(ind.I ./ sz))
            res[i][ind] = T(ind.I[i] - 1)
        end
    end

    sz_2 = div.(sz, 2)

    for i = 1:N
        interpolatemod!(res2[i], res[i], dec, t_interp, T(sz[i]))
    end

    dec2 = getinverse(dec, t_interp)

    for i = 1:N
        interpolatemod!(res3[i], res2[i], dec2, t_interp, T(sz[i]))

        res2[i] .= mod.(res[i] .- res3[i] .+ sz_2[i], sz[i]) .- sz_2[i]
    end

    note = norm(res2)
    @show note, sqrt(eps(T))
    @test note < sqrt(eps(T))

end

@testset "test inverse" begin
    T = Double64
    test_getinv([Lagrange(11, T), Lagrange(11, T)], T(0.00911), (128, 100))
    # test_inv0([Lagrange(11,T), Lagrange(11,T)], (zero(T),zero(T)), (20,30))
    # test_inv0([Lagrange(11,T), Lagrange(11,T)], (T(pi)/10,T(pi)/9), (20,30))
    # test_inv([Lagrange(11, T), Lagrange(11, T)], T(0.011), (20, 30))
end

function test_interpfloat(
    interp::AbstractInterpolation{T,edge},
    sz,
    tol,
    nb = 100,
) where {T,edge}

    tabdec =
        T.([
            big"0.345141526199181716726626262655544",
            -big"0.3859416191876155241320011187619",
            -big"1.28561390114441619187615524132001118762519",
            -big"0.885901390114441619187615524132001118762519",
            -big"5.678513256790098898776656565545454544544545", # only for edge==CircEdge
            big"4.9876651456677809099887665655556565565656565", # only for edge==CircEdge
            big"0.186666659416191876155241320011187619",
            big"0.590999232323232323232365566787878898898",
            big"1.231098015934444444444444788888888878878",
        ])
    mesh = T.(collect(big.(0:(sz-1))) / sz)

    fct1(x) = T(cos(2big(pi) * x + big"0.25"))
    fct2(x) = T(exp(-(cos(2big(pi) * x + big"0.25") - 1)^2))
    tabfct = [fct1, fct2]
    @time @testset "test interpolation  $interp" begin
        nmax = 0

        for (i_fct, fct) in enumerate(tabfct), dec in tabdec
            deb = fct.(mesh)
            #       println("deb=$deb")
            fp = deb
            fi = zeros(T, sz)

            decint = convert(Int, floor(dec))
            value = dec - decint
            if get_order(interp) % 2 == 0 && value > 0.5
                value -= 1
                decint += 1
            end
            if edge != CircEdge && abs(decint) > 2
                continue
            end
            precal =
                edge == CircEdge ? getprecal(interp, value) :
                get_allprecal(interp, decint, value)
            for i = 1:nb
                fi .= fp
                ref = fct.(mesh .+ i * dec / sz)
                interpolate!(fp, fi, decint, precal, interp)
                nmax = max(nmax, float(norm(fp - ref, Inf)))
                if isbspline(interp) && edge != CircEdge
                    @show typeof(interp), i_fct, dec, sz, nb, nmax
                else
                    @test isapprox(fp, ref, atol = tol)
                end
            end
        end
        @show typeof(interp), sz, nb, nmax
    end
end
T = Double64
test_interp2d([Lagrange(11, T), Lagrange(11, T)], T(1.25), (128, 100))

test_interp(Lagrange(3, Rational{BigInt}; edge = InsideEdge), big"3" // 1024, 128)

test_interp(Lagrange(3, Rational{BigInt}; edge = CircEdge), 128)

test_interp(B_SplineLU(3, 128, Rational{BigInt}), 128)

test_interpfloat(Lagrange(3, BigFloat, edge = CircEdge), 128, 1e-3, 100)
test_interpfloat(Lagrange(3, Float64, edge = CircEdge), 128, 1e-3, 100)

test_interpfloat(Lagrange(7, BigFloat, edge = InsideEdge), 128, 1e-3, 3)
test_interpfloat(Lagrange(7, Float64, edge = InsideEdge), 128, 1e-3, 3)

test_interpfloat(Lagrange(21, BigFloat, edge = CircEdge), 256, 1e-20)
test_interpfloat(Lagrange(9, Float64, edge = CircEdge), 256, 1e-10)

test_interpfloat(Lagrange(4, BigFloat, edge = CircEdge), 256, 1e-5)
test_interpfloat(Lagrange(4, Float64, edge = CircEdge), 256, 1e-5)

test_interpfloat(Lagrange(22, BigFloat, edge = CircEdge), 256, 1e-20)
test_interpfloat(Lagrange(12, Float64, edge = CircEdge), 256, 1e-10)

test_interpfloat(B_SplineLU(3, 256, BigFloat), 256, 1e-5)
test_interpfloat(B_SplineLU(3, 256, Float64), 256, 1e-5)

test_interpfloat(B_SplineLU(21, 256, BigFloat), 256, 1e-30)
test_interpfloat(B_SplineLU(11, 256, Float64), 256, 1e-12)

test_interpfloat(B_SplineFFT(3, 256, BigFloat), 256, 1e-5)
test_interpfloat(B_SplineFFT(3, 256, Float64), 256, 1e-5)

test_interpfloat(B_SplineFFT(21, 256, BigFloat), 256, 1e-30)
test_interpfloat(B_SplineFFT(11, 256, Float64), 256, 1e-12)


text:
"""Passed to NMR.runjob() to specify a MapReduce job"""
struct Nmr
  jobid::Int
  inputs::Array{String}
  fun::Function
  outputfilename::String
end

"""Exception designed to be thrown by mappers when a runtime error occurs"""
struct MapperException <: Exception
  msg::String
end

"""Exception designed to be thrown by type initiations when a runtime error occurs"""
struct UDFException <: Exception
  msg::String
end