text:
struct Tile{T <: AbstractString}
    id::Int
    data::Vector{T}
    orientation::Symbol
    rotation::Int
end

function Base.show(io::IO, t::Tile)
    print(io, "$(t.id),$(t.orientation),$(t.rotation)")
end

function read_data(filename)
    blocks = split(String(read(filename)), "\n\n")
    function parse_block(block)
        lines = split(block, "\n")
        id = parse(Int, match(r"Tile (\d+):", lines[1]).captures[1])
        data = lines[2:end]
        return id, data, :original, 0
    end
    tiles = [Tile(parse_block(block)...) for block in blocks]
    return Dict(t.id => t for t in tiles)
end

sides(t::Tile) = let len = length(t.data)
    top = t.data[1]
    bottom = t.data[end]
    left = join([t.data[i][1] for i ∈ 1:len])
    right = join([t.data[i][end] for i ∈ 1:len])
    top, bottom, left, right
end

# n = 0, 1, 2, 3 (rotating left by 90 degrees)
# Pass n as -1 to use the tile's own rotation setting
function hashes(t::Tile, n::Int = -1, orientation::Symbol = :original)
    the_tile = if orientation == :original
        t
    elseif orientation == :horizontal
        hflip(t)
    else
        vflip(t)
    end
    top, bottom, left, right = sides(the_tile)
    rotation = n >= 0 ? n : the_tile.rotation
    if rotation == 0
        [hash(top), hash(bottom), hash(left), hash(right)]
    elseif rotation == 1
        [hash(right), hash(left), hash(reverse(top)), hash(reverse(bottom))]
    elseif rotation == 2
        [hash(reverse(bottom)), hash(reverse(top)), hash(reverse(right)), hash(reverse(left))]
    elseif rotation == 3
        [hash(reverse(left)), hash(reverse(right)), hash(bottom), hash(top)]
    else
        error("Incorrect usage")
    end
end

# Flip vertically
vflip(t::Tile) = Tile(
    t.id,
    t.data[end:-1:1],
    :vertical,
    t.rotation
)

# Flip horizontally
hflip(t::Tile) = Tile(
    t.id,
    [reverse(r) for r in t.data],
    :horizontal,
    t.rotation
)

# Match a tile against all other tiles.
# No need to transform current tile as long as we rotate/flip others.
function match_tile(t::Tile, tiles::Dict)
    # @show "Matching tile $(t.id)"
    side_hashes = hashes(t, 0)
    matched = Set{Int}()
    for (j, sh) in enumerate(side_hashes)
        for candidate in values(tiles)
            t == candidate && continue  # skip myself
            for x in [candidate, hflip(candidate), vflip(candidate)]
                for i in 0:3
                    csh = hashes(x, i)
                    idx = findfirst(==(sh), csh)
                    if idx !== nothing
                        push!(matched, x.id)
                    end
                end
            end
        end
    end
    return matched
end

function match_all(tiles::Dict)
    return Dict(tile.id => match_tile(tile, tiles) for tile in values(tiles))
end

# Given a tile `x`, try to fit tile `y` by rotating/flipping `y`.
# Only need to fit `y` by aligning its top or left to `x`.
# Return the orientation and rotatiion values once it's fitted.
function fit(x::Tile, y::Tile, side::Symbol)
    if side == :top
        x_bottom_hash = hashes(x, x.rotation, x.orientation)[2]
        # x.id == 2213 && @info "x (2213)" x hashes(x)
        for (orientation, this_y) in zip([:original, :horizontal, :vertical], [y, hflip(y), vflip(y)])
            for rotation in 0:3
                csh = hashes(this_y, rotation)[1]  # top hash
                # x.id == 2213 && @info "this_y ($(y.id))" hashes(this_y, i)
                if csh == x_bottom_hash
                    return orientation, rotation
                end
            end
        end
    elseif side == :left
        x_right_hash = hashes(x, x.rotation, x.orientation)[4]
        # x.id == 1021 && @info "x (1021)" x hashes(x)
        for (orientation, this_y) in zip([:original, :horizontal, :vertical], [y, hflip(y), vflip(y)])
            for rotation in 0:3
                csh = hashes(this_y, rotation)[3]  # left hash
                if csh == x_right_hash
                    return orientation, rotation
                end
            end
        end
    end
    return nothing
end

# Align `y` against `x` on a specific `side`.
# Return a new Tile object transformed from `y`.
# The tile `y` in `tiles` is mutated to the new orientation.
function align!(tiles::Dict, x::Tile, y::Tile, side::Symbol)
    result = fit(x, y, side)
    if result !== nothing
        orientation, rotation = result
        oriented_tile = Tile(y.id, y.data, orientation, rotation)
        tiles[y.id] = oriented_tile
        return oriented_tile
    end
    return nothing
end

# Assemble the puzzle!
# Since we don't know which corner piece is the top-left one,
# and our algorithm requires to start from top-left, we can
# just try all corner pieces.
function assemble(tiles::Dict, mappings::Dict)
    corners = [k for (k,v) in mappings if length(v) == 2]
    for c in corners
        image, picture = try_assemble(tiles, mappings, c)
        if count(isnothing, picture) == 0
            return image, picture
        end
    end
    error("Too bad, cannot assemble the puzzle")
end

# Bootstrap the assembly process by taking a corner piece and
# figure out how it should be oriented, and then properly
# match its adjacent pieces.
function bootstrap(tiles::Dict, mappings::Dict, corner::Int)
    tile = tiles[corner]
    n1, n2 = [tiles[id] for id in mappings[corner]]
    for t in (tile, hflip(tile), vflip(tile))
        for rotation in 0:3
            r1 = fit(tile, n1, :top)
            if r1 !== nothing  # good fit for (1,1) and (2,1)
                good_corner = Tile(tile.id, tile.data, t.orientation, rotation)
                good_bottom = Tile(n1.id, n1.data, r1[1], r1[2])
                r2 = fit(tile, n2, :left)
                good_right = Tile(n2.id, n2.data, r2[1], r2[2])
                return good_corner, good_bottom, good_right
            end
        end
    end
    return nothing
end

#=
Try to assemble the puzzle using the specified top-left `corner` piece.

The algorithm works as follows:
1. Align the top-left corner piece (1,1) and its adjacent pieces (2,1) & (1,2)
2. Assemble pieces in a diagonal manner (from lower left to upper right)
   - the left-edge piece can be aligned with its top neighbor
   - the middle pieces can be aligned with both its top & left neighbors
   - the top-edge piece can be aligned with its left neighbor

The loop index is designed to go out of bound due to the diagonal sweeping
process. The if-conditions are are used to skip out-of-bound cells.

Returns a tuple with these elements:
1. `image` - a matrix of tile id's
2. `picture` - a matrix of tiles
=#
function try_assemble(tiles::Dict, mappings::Dict, corner::Int)
    # create an image array
    sz = round(Int, sqrt(length(tiles)))
    image = fill(0, (sz, sz))

    # aligned tiles
    picture = Matrix{Tile}(undef, (sz, sz))

    # Bootstrap and align the 3 top-left pieces
    top_left_pieces = bootstrap(tiles, mappings, corner)
    picture[1,1], picture[2,1], picture[1,2] = top_left_pieces
    image[1,1], image[2,1], image[1,2] = getproperty.(top_left_pieces, :id)

    # Mutate the tiles
    tiles[image[1,1]] = top_left_pieces[1]
    tiles[image[2,1]] = top_left_pieces[2]
    tiles[image[1,2]] = top_left_pieces[3]

    # Use a Set to remember which tile has already been placed
    visited = union(corner, mappings[corner])

    for row in 3:(2sz-1) # intentionally go out of bound for the diagnoal swipe
        pos = (row, 1)
        for i in 1:row-2   # discount 2 edge tiles
            pos = (pos[1]-1, pos[2]+1) # move up-right
            if pos[1] > sz || pos[2] > sz
                continue  # skip out of bound imaginary space
            end
            left = (pos[1], pos[2]-1)  # left neighbor (col - 1)
            top  = (pos[1]-1, pos[2])  # top neighbor (row - 1)
            # @show pos left top
            left_tile = image[left[1], left[2]]
            top_tile = image[top[1], top[2]]
            me = intersect(mappings[left_tile], mappings[top_tile])
            me = setdiff(me, visited)
            @assert length(me) == 1
            me_id = pop!(me)
            image[pos[1], pos[2]] = me_id
            push!(visited, me_id)
            picture[pos[1], pos[2]] = align!(tiles, tiles[left_tile], tiles[me_id], :left)
            # println("Found tile at row=$row i=$i pos=$pos => $me_id")
        end
        # place tile at (row, 1). It has to be the only unvisited top neighbor's neighbors
        if row <= sz
            top_id = image[row - 1, 1]
            # @show "Looking for top edge title" row top_id visited
            me = setdiff(mappings[top_id], visited)
            @assert length(me) == 1
            me_id = pop!(me)
            image[row, 1] = me_id
            push!(visited, me_id)
            picture[row, 1] = align!(tiles, tiles[top_id], tiles[me_id], :top)
            # println("Found left edge tile ($(row), 1)=> $me_id")
        end
        # place tile at (1, row)
        if row <= sz
            # @show "Looking for left edge title"
            left_id = image[1, row - 1]
            me = setdiff(mappings[left_id], visited)
            @assert length(me) == 1
            me_id = pop!(me)
            image[1, row] = me_id
            push!(visited, me_id)
            picture[1, row] = align!(tiles, tiles[left_id], tiles[me_id], :left)
            # println("Found top edge tile (1, $(row)) => $me_id")
        end
    end
    return image, picture
end

function matrix(t::Tile)
    sz = length(t.data)
    return [r[j] == '#' ? 1 : 0 for (i,r) in enumerate(t.data), j in 1:sz]
end

function flipped(A::AbstractMatrix{T}, orientation::Symbol) where T
    A = copy(A)
    if orientation == :horizontal
        return A[:, 1:end] = A[:, end:-1:1]
    elseif orientation == :vertical
        return A[1:end, :] = A[end:-1:1, :]
    else
        return A
    end
end

# Generally, rotating a matrix *clockwise* involves two steps:
# 1. transpose the matrix
# 2. flip the columns
#
# Convention: n = 0, 1, 2, 3 (counter-clockwise)
function rotated(A::AbstractMatrix{T}, n::Int) where T
    # adjust number of rotations since we're doing counter-clockwise
    m = (4 - n) % 4
    A = copy(A)
    for i in 1:m
        A = transpose(A)
        A[:, 1:end] = A[:, end:-1:1]
    end
    return A
end

# Return a bitmap matrix for the tile that is flipped/rotated accordingly.
function oriented(t::Tile)
    A = matrix(t)
    A = flipped(A, t.orientation)
    A = rotated(A, t.rotation)
    return A
end

# Get rid of the border of a bitmap matrix.
remove_border(A::AbstractMatrix) = A[2:end-1,2:end-1]

# Making a bitmap for a picture (2D array of tiles)
function make_bitmap(picture)
    M = remove_border.(oriented.(picture))
    return vcat([hcat(M[r, :]...) for r in 1:size(picture,1)]...)
end

# Parse a monster stencil as a bitmap array
function monster_stencil()
    str =   """
                  #
#    ##    ##    ###
 #  #  #  #  #  #   """
    parse_line(s) = [c == ' ' ? 0 : 1 for c in s]
    vcat(map(transpose, parse_line.(split(str, "\n")))...)
end

# Check if the specific location at `A[row,col]` (top-left)
# contains a monster image. Apply the stencil using bitwise-and
# operation and see if the result contains the same number of 1's
# as in the stencil.
function is_monster(stencil, A, row, col)
    r, c = size(stencil)
    area = A[row:row+r-1, col:col+c-1]
    masked_area = stencil .& area
    return sum(masked_area) == sum(stencil)
end

# Find all monsters in a bitmap. Return the locations of those
# monsters (top-left position)
function find_monsters(stencil, bitmap)
    br, bc = size(bitmap)
    sr, sc = size(stencil)
    locations = []
    for i in 1:br-sr+1
        for j in 1:bc-sc+1
            if is_monster(stencil, bitmap, i, j)
                # @info "Found monster: bitmap=$(hash(bitmap)) i=$i j=$j"
                push!(locations, (i,j))
            end
        end
    end
    return locations
end

# Mark the monsters on the bitmap with a value of 2.
function mark_monsters!(stencil, bitmap, locations)
    for loc in locations
        for row in 1:size(stencil, 1)
            for col in 1:size(stencil, 2)
                if stencil[row, col] == 1
                    bitmap[loc[1] + row - 1, loc[2] + col - 1] = 2
                end
            end
        end
    end
end

# Flip and rotate the bitmap and see if we can locate any monster.
# If so, just return the new bitmap and the locations of monsters.
function orient_and_find_monsters(stencil, bitmap)
    for bm in (bitmap, flipped(bitmap, :horizontal), flipped(bitmap, :vertical))
        for rotation in 0:3
            bm2 = rotated(bm, rotation)
            locations = find_monsters(stencil, bm2)
            if length(locations) > 0
                # @info "Found monster in bitmap=$(hash(bm2))"
                return bm2, locations
            end
        end
    end
    error("can't find any monsters :-(")
end

function part1()
    tiles = read_data("day20.txt")
    result = match_all(tiles)
    return prod([k for (k,v) in result if length(v) == 2])
end

function part2()
    tiles = read_data("day20.txt")
    mappings = match_all(tiles)
    # @info "mappings" mappings
    puzzle_ids, puzzle_picture = assemble(tiles, mappings)

    bitmap = make_bitmap(puzzle_picture)
    stencil = monster_stencil()

    bitmap, locations = orient_and_find_monsters(stencil, bitmap)

    mark_monsters!(stencil, bitmap, locations)
    return count(==(1), bitmap)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct AdminConversationsEkmApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_adminConversationsEkmListOriginalConnectedChannelInfo(_api::AdminConversationsEkmApi; channel_ids=nothing, team_ids=nothing, limit=nothing, cursor=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/admin.conversations.ekm.listOriginalConnectedChannelInfo", ["slackAuth"])
    Swagger.set_param(_ctx.query, "channel_ids", channel_ids)  # type String
    Swagger.set_param(_ctx.query, "team_ids", team_ids)  # type String
    Swagger.set_param(_ctx.query, "limit", limit)  # type Int32
    Swagger.set_param(_ctx.query, "cursor", cursor)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

List all disconnected channels—i.e., channels that were once connected to other workspaces and then disconnected—and the corresponding original channel IDs for key revocation with EKM.
Param: channel_ids::String
Param: team_ids::String
Param: limit::Int32
Param: cursor::String
Return: Dict{String, Any}
"""
function adminConversationsEkmListOriginalConnectedChannelInfo(_api::AdminConversationsEkmApi; channel_ids=nothing, team_ids=nothing, limit=nothing, cursor=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminConversationsEkmListOriginalConnectedChannelInfo(_api; channel_ids=channel_ids, team_ids=team_ids, limit=limit, cursor=cursor, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminConversationsEkmListOriginalConnectedChannelInfo(_api::AdminConversationsEkmApi, response_stream::Channel; channel_ids=nothing, team_ids=nothing, limit=nothing, cursor=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminConversationsEkmListOriginalConnectedChannelInfo(_api; channel_ids=channel_ids, team_ids=team_ids, limit=limit, cursor=cursor, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export adminConversationsEkmListOriginalConnectedChannelInfo

ext: jl


text:
# #
# Reproduce the figures in the Appendix of
# [Duchêne](https://www.ams.org/open-math-notes/omn-view-listing?listingId=111309)
# comparing solutions to the water waves system with predictions of
# the Serre-Green-Naghdi and Whitham-Green-Naghdi systems
# #
export Integrate,Figures
using WaterWaves1D,Plots,LinearAlgebra;
#using JLD #(uncomment if using @save)

#----
"""
	Integrate(scenario;kwargs)

Integrate in time the water waves system as well as
the Serre-Green-Naghdi, Whitham-Green-Naghdi and Isobe-Kakinuma models
with an initial datum depending on the argument (`scenario`).
If `scenario==1`, the initial datum is a heap of water with zero velocity.
If `scenario==2`, the initial datum is generates at first order a right-going wave.

Other arguments are optional:
- `μ` is the shallowness parameter,
- `ϵ` in the nonlinearity parameter,
- `p` is the power in the initial datum `η = exp(-|x|^p)`
- `N` the number of collocation points,
- `L` the half-length of the mesh,
- `T` the final time of integration,
- `dt` the timestep,
- `dealias`: dealiasing with Orlicz rule `1-dealias/(dealias+2)` (default is `0`, i.e. no dealiasing),
- `iterate`: use GMRES if `true`, and LU decomposition otherwise (default is `true`),
- `precond` gives some choice in the preconditioner for GMRES (default is `true` for the default preconditioner),
- `method` and `maxiter` are the corresponding options for generating the initial data in the water waves model; see doc of `WaterWaves`
- `name`: a string used to save the figures.

Return `(Diff,problems,plt)` where `Diff` is a measure of the errors, `problems` contains all the information and `plt` is a plot of the final time solution.
"""
function Integrate(scenario;μ=.1,ϵ=.1,p=2,N=2^10,L=15,T=10,dt=0.001,dealias=0,iterate=true,precond=true,method=1,maxiter=100,name=nothing)
	#---- preparation

	if name != nothing ns=floor(Int,max(1,T/dt/100)) else ns=1 end

	if typeof(dealias)==Int64 dealias=dealias*ones(Int64,4) end

	# the set of parameters as indicated by the user
	param = ( μ  = μ, ϵ  = ϵ,
				N  = N, L  = L,
	            T  = T, dt = dt,
				ns=ns )

	# generate the initial data
	mesh=Mesh(param)
	k=mesh.k

	η = exp.(-abs.(mesh.x ).^p)
	if scenario == 1
		v = 0*η
	elseif scenario == 2
		v = 2/ϵ*(sqrt.(1 .+ ϵ*η) .- 1)
	else
		@error("the argument must be 1 or 2")
	end
	init     = Init(mesh,η,v)

	# construct the models to be solved
	modelWW	 = WaterWaves(param;dealias=dealias[1],method=method,maxiter=maxiter)
	modelSGN = SerreGreenNaghdi(param;ktol=0, gtol=1e-14, iterate=iterate, precond = precond, dealias = dealias[2])
	modelWGN = WhithamGreenNaghdi(param;SGN=false, ktol=0, gtol=1e-14, iterate=iterate, precond = precond, dealias = dealias[3])
	modelIK2 = IsobeKakinuma(param; ktol=0, gtol=1e-14, iterate=iterate, precond = precond, dealias = dealias[4])

	# construct the associated problems
	problems = []
	push!(problems, Problem(modelWW, init, param) )
	push!(problems, Problem(modelSGN, init, param) )
	push!(problems, Problem(modelWGN, init, param) )
	push!(problems, Problem(modelIK2, init, param) )

	#---- computation
	for problem in problems
		solve!( problem )
	end

	#---- analysis
	# Errors at final time
	(ηww,vww,xww,t)=solution(problems[1])
	(ηSGN,vSGN)=solution(problems[2],x=xww)
	(ηWGN,vWGN)=solution(problems[3],x=xww)
	(ηIK2,vIK2)=solution(problems[4],x=xww)

	Diff = [norm(ηww-ηSGN,Inf) norm(ηww-ηWGN,Inf) norm(ηww-ηIK2,Inf)]
	@info string("Error of SGN, WGN, IK2: ",round.(Diff, sigdigits=3))

	#---- plots
	# solutions at final time
	plt = plot_solution(problems;fourier=true)
	display(plt)

	if name != nothing
		#save(problem,name);
		#@save(name,xww,ηww,ηSGN,ηWGN,ηIK2,scenario,μ,ϵ,p,N,L,T,dt,dealias,iterate,precond,name)

		if scenario == 1
			xlims=(0,last(xww))
			ind=range(1,step=length(xww)÷32,length=32);
		else
			xlims=(first(xww),last(xww))
			ind=range(1,step=length(xww)÷16,length=16);
		end

		# Plot and save the solutions at final time
		plot_solution(problems;fourier=false)
		plot!(xlims=xlims,legend=:topleft)
		savefig(string(name,".pdf"));

		# Plot and save the differences at final time
		plot(xww, [ηww-ηSGN ηww-ηWGN ηww-ηIK2],
		 	color=[:2 :3 :4],primary=false)
		scatter!(xww[ind], [ηww[ind]-ηSGN[ind] ηww[ind]-ηWGN[ind] ηww[ind]-ηIK2[ind]],
		 	color=[:2 :3 :4],marker=([:d :h :ut], 8, Plots.stroke(1, :gray)),#linewidth=0,
		 	xlims=xlims,
		   title=string("difference (surface deformation) at t=",t),
		   label=["Serre-Green-Naghdi" "Whitham-Green-Naghdi" "Isobe-Kakinuma"])
		savefig(string(name,"-differences.pdf"));

		# Generate and save an animation of the evolution in time
		create_animation(problems;fourier=false,xlims=xlims,ylims=(-0.5,1),name=name)

	end

	return Diff,problems,plt
end

#----
"""
	Figures()

Reproduce the figures in Section I.5 of the monograph
[Many Models for Water Waves](https://www.ams.org/open-math-notes/omn-view-listing?listingId=111309)
"""
function Figures()
	MU=[1 0.1 0.01]
	EPS=[0.5 0.25 0.125]

	DiffSGN=zeros(3,3);DiffWGN=zeros(3,3);DiffIK2=zeros(3,3);

	for i in 1:3
		for j in 1:3
			mu=MU[i];eps=EPS[j]
			if mu == 1. mu = 1 end
			name=string("mu",mu,"eps",eps)
			name=replace(name, '.' => "")
			@info string("Save plots and animation with the name ",name)
			Diff,=Integrate(1;p=3,μ=mu,ϵ=eps,N=2^10,T=10,L=15,dt=0.001,dealias=0,maxiter=100,name=name);
			DiffSGN[i,j]=Diff[1];DiffWGN[i,j]=Diff[2];DiffIK2[i,j]=Diff[3];
			#@save("MM4WW",DiffSGN,DiffWGN,DiffIK2)
		end
	end
end
nothing


text:
# This file implements methods from the paper
# "Boltzmann Encoded Adversarial Machines" of Fisher et al. (2018)

function gibbssample_gamma!(particles::Particles, gbrbm::GaussianBernoulliRBM2, nsteps::Int = 10;
      autocorcoeff::Float64 = 0.9, betasd::Float64 = 0.9)

   nsamples = size(particles[1], 1)
   vs = map(i -> view(particles[1], i, :), 1:nsamples)
   hs = map(i -> view(particles[2], i, :), 1:nsamples)

   for k = 1:nsamples
      v = vs[k]
      h = hs[k]
      beta = initgammaprocess(autocorcoeff, betasd)

      for i = 1:nsteps
         beta = samplegammaprocess(beta, autocorcoeff, betasd)
         gbrbm2 = GaussianBernoulliRBM2(gbrbm.weights, gbrbm.visbias, gbrbm.hidbias,
               copy(gbrbm.sd) / sqrt(beta))
         samplevisible!(v, gbrbm2, h)
         samplehidden!(h, gbrbm, v, 1/beta)
      end
   end
   particles
end

"""
Try Temperature Driven Sampling for NegativeBinomialBernoulliRBM
"""

#function gibbssample_gamma!(particles::Particles, rbm::NegativeBinomialBernoulliRBM, nsteps::Int = 10;
#   autocorcoeff::Float64 = 0.9, betasd::Float64 = 0.9)

#nsamples = size(particles[1], 1)
#vs = map(i -> view(particles[1], i, :), 1:nsamples)
#hs = map(i -> view(particles[2], i, :), 1:nsamples)

#for k = 1:nsamples
#   v = vs[k]
#   h = hs[k]
#   beta = initgammaprocess(autocorcoeff, betasd)
#
#   for i = 1:nsteps
#      beta = samplegammaprocess(beta, autocorcoeff, betasd)
#      samplevisible!(v, rbm, h)
#      samplehidden!(h, rbm, v, 1/beta)
#   end
#end
#particles
#end

function nearestneighbourcritic(xdata::M, xmodel::M, k::Int,
      distanceweighted::Bool = true) where {M <: AbstractArray{Float64,2}}

   ndatasamples = size(xdata, 1)
   if ndatasamples != size(xmodel, 1)
      error("Number of samples in `xdata` and `xmodel` must be equal.")
   end

   x = vcat(xdata, xmodel)

   nsamples = 2 * ndatasamples

   # Distancematrix is filled with the distances
   # between all combinations of samples in x.
   distancematrix = Matrix{Float64}(nsamples, nsamples)
   for i = 1:nsamples
      for j = (i+1):nsamples
         distancematrix[i, j] = distancematrix[j, i] = norm(x[j, :] - x[i, :])
      end
      distancematrix[i, i] = 0.0
   end

   ret = Vector{Float64}(ndatasamples)
   eps = 0.1
   if distanceweighted
      for i = 1:ndatasamples
         knearestindices = selectperm(distancematrix[i, :], (1:k)+1)
         nearestdatasamplesindices = knearestindices[knearestindices .<= ndatasamples]
         datasamplesdistances = distancematrix[i, nearestdatasamplesindices]
         ret[i] = 2.0 * ( sum(1.0 ./ (datasamplesdistances .+ eps)) /
               sum(1.0 ./ (distancematrix[i, knearestindices] .+ eps))) - 1
      end
   else
      for i = 1:ndatasamples
         knearestindices = selectperm(distancematrix[i, :], (1:k)+1)
         nnearestdatasamples = sum(knearestindices .<= ndatasamples)
         ret[i] = 2 * nnearestdatasamples / k - 1
      end
   end

   ret
end

"""
Update beta with an autoregressive Gamma process.
beta_0 ~ Gamma(nu,c/(1-phi)) = Gamma(1/var, var)
h_t ~ Possion( phi/c * h_{t-1})
beta_t ~ Gamma(nu + z_t, c)
Achieves a stationary distribution with mean 1 and variance var:
Gamma(nu, var) = Gamma(1/var, var)
"""
function initgammaprocess(autocorr_coeff::Float64,betasd::Float64)

   var = betasd.^2
   nu = 1 ./ var
   c = (1-autocorr_coeff).*var

   # beta_0 ~ Gamma(nu,c/(1-phi)) = Gamma(1/var, var)
   beta = Gamma(nu,c ./ (1-autocorr_coeff))
   beta_sample = rand(beta)
end

function samplegammaprocess(beta_sample::Float64, autocorr_coeff::Float64,
      beta_std::Float64)

   var = beta_std^2

   nu = 1/var
   c = (1-autocorr_coeff)*var

   # h_t ~ Possion( phi/c * h_{t-1})
   lambda = beta_sample*autocorr_coeff/c
   z = Poisson(lambda)
   z_sample = rand(z)

   # beta_t ~ Gamma(nu + z_t, c)
   a = nu + z_sample
   b = c
   beta = Gamma(a,b)
   beta_sample = rand(beta)

   beta_sample
end


text:
using JLD2, Dates, ProgressMeter
"""
    run!(
        ρ,u,v;
        Lx,Ly,Δt,
        c1, c2, λ,
        final_time,
        bcond_x="periodic",bcond_y="periodic",
        Fx=nothing,Fy=nothing,
        method="HLLE",
        simu_name="simu",
        should_save=false,save_step=1,
        should_plot=false,plot_step=1,
        save_video=false,fps=40,
        range=(0.,5.),
        kwargs...
    )

Run a simulation of the SOH system starting from the data `(ρ,u,v)`. The domain
grid has the following structure.

```
|------- |----------------|---------------|----------------------|--------|
| ghost  |    ghost       |       ...     |    ghost             | ghost  |
|--------║================|===============|======================║--------|
| ghost  ║  (2,ncelly+1)  |      ...      | (ncellx+1,ncelly+1)  ║ ghost  |
|--------║----------------|---------------|----------------------║--------|
| ghost  ║  (2,ncelly)    |      ...      | (ncellx+1,ncelly)    ║ ghost  |
|--------║----------------|---------------|----------------------║--------|
|        ║                |               |                      ║        |
|--------║----------------|---------------|----------------------║--------|
| ghost  ║    (2,3)       |      ...      |   (ncellx+1,3)       ║ ghost  |
|--------║----------------|---------------|----------------------║--------|
| ghost  ║    (2,2)       |      ...      |   (ncellx+1,2)       ║ ghost  |
---------║================|===============|======================║--------|
| ghost  |    ghost       |      ...      |    ghost             | ghost  |
|--------|----------------|---------------|----------------------|--------|
````

The ghost cells (first and last column and rows) are used for the boundary
conditions.

# Arguments
- `ρ` -- density (matrix of size `ncellx*ncelly`)
- `u` -- x-coordinate of the velocity (matrix of size `ncellx*ncelly`)
- `v` -- y-coordinate of the velocity (matrix of size `ncellx*ncelly`)
- `Lx` -- length of the domain along the x-axis
- `Ly` -- length of the domain along the y-axis
- `Δt` -- time-step
- `c1` -- coefficient c1
- `c2` -- coefficient c1
- `λ` -- coefficient λ
- `final_time` -- final time of the simulation
- `bcond_x` -- (optional, default : `"periodic"`) boundary condition along the x-axis
- `bcond_y` -- (optional, default : `"periodic"`) boundary condition along the y-axis
- `Fx` -- (optional, default : `nothing`) if specified, x-coordinate of the exterior force
- `Fy` -- (optional, default : `nothing`) if specified, y-coordinate of the exterior force
- `method` -- (optional, default : `"HLLE`) can be either `"Roe"` or `"HLLE"`
- `simu` -- (optional, default : `"simu"`) name of the simulation directory
- `should_save` -- (optional, default : `false`) flag for saving the data (ρ,u,v)
- `save_step` -- (optional, default : `1`) number of iterations between two saved data
- `should_plot` -- (optional, default : `false`) flag for plotting and saving the heatmap
- `plot_step` -- (optional, default : `1`) number of iterations between two saved plots
- `save_video` -- (optional, default : `false`) flag for saving a video
- `fps` -- (optional, default : `40`) frame per seconds
- `range` -- (optional, default : `(0.0,5.0)`) range of the density on the heatmap
- `kwargs` -- Keyword arguments to be passed in the function `plot_rhoUV`
"""
function run!(
    ρ,u,v;
    Lx,Ly,Δt,
    c1, c2, λ,
    final_time,
    bcond_x="periodic",bcond_y="periodic",
    Fx=nothing,Fy=nothing,
    method="HLLE",
    simu_name="simu",
    should_save=false,save_step=1,
    should_plot=false,plot_step=1,
    save_video=false,fps=40,
    range=(0.,5.),
    kwargs...
)
    #=======================================================================================#
    #============ Collect and display the parameters =======================================#
    #=======================================================================================#

    ncellx = size(ρ)[1] - 2
    Δx = Lx/ncellx
    ncelly = size(ρ)[2] - 2
    Δy = Ly/ncelly
    parameters = Dict{Symbol,Any}(
        :Lx => Lx, :Ly => Ly, :ncellx => ncellx, :ncelly => ncelly,
        :Δx => Δx, :Δy => Δy, :Δt => Δt,
        :c1 => c1, :c2 => c2, :λ => λ, :CFL => c1 * Δt/Δx * sqrt(λ/(c1-c2)),
        :final_time => final_time,
        :bcond_x => bcond_x, :bcond_y => bcond_y,
        :Fx => Fx, :Fy => Fy,
        :method => method,
        :date => string(Dates.now())
    )

    println("\n************* Model parameters **************\n")
    println("c1 = $(parameters[:c1])")
    println("c2 = $(parameters[:c2])")
    println("λ = $(parameters[:λ])")
    println("\n*********************************************\n")

    println("************* Domain parameters *************\n")
    println("Lx = $(parameters[:Lx])")
    println("Ly = $(parameters[:Ly])")
    println("Boundary condition in x = $(parameters[:bcond_x])")
    println("Boundary condition in y = $(parameters[:bcond_y])")
    println("Exterior force : $(!isnothing(parameters[:Fx]) && !isnothing(parameters[:Fy]))")
    println("\n*********************************************\n")

    println("************* Numerical method **************\n")
    println("Δx = $(parameters[:Δx])")
    println("Δy = $(parameters[:Δy])")
    println("Δt = $(parameters[:Δt])")
    println("CFL = $(parameters[:CFL])")
    println("Scheme : $(parameters[:method])")
    println("Final time : $(parameters[:final_time])")
    println("\n*********************************************\n")

    println("************* Saving parameters *************\n")
    if should_save
        if save_step == 1
            println("Save all data")
        else
            println("Save data every $save_step iterations")
        end
        estimated_size = round(8*ncellx*ncelly*3*floor(Int, final_time / Δt)/save_step*1e-9, digits=2)
        println("Estimated size : $estimated_size GB")
    else
        println("Only the initial and final data will be saved")
    end
    if should_plot
        if plot_step == 1
            println("Plot all data")
        else
            println("A plot will be save every $plot_step iterations")
        end
    else
        println("No plot will be saved")
    end
    if save_video
        println("A video with $fps fps will be saved")
    else
        println("No video will be saved")
    end
    println("\n*********************************************\n")

    #=======================================================================================#
    #============ Create the simulation directory and save the initial data ================#
    #=======================================================================================#

    dir_name = make_new_dir(simu_name)
    data_dir = mkdir(joinpath(dir_name,"data"))
    data = Dict(:ρ => ρ, :u => u, :v => v)
    save(joinpath(data_dir,"data_0.jld2"),"iter_0",data)

    if should_plot || save_video
        print("Initializing plot... ")
        plot_dir = make_new_dir(joinpath(dir_name,"plots"))
        if save_video
            fig, ax, hm, arrows, stream = plot_rhoUV(ρ,u,v,Δx,Δy,range,should_plot,save_video,plot_dir,"0.png";fps=fps,kwargs...)
        else
            fig, ax, hm, arrows = plot_rhoUV(ρ,u,v,Δx,Δy,range,should_plot,save_video,plot_dir,"0.png";kwargs...)
        end
        println("Done.\n")
    end

    #=======================================================================================#
    #============ The big loop =============================================================#
    #=======================================================================================#

    flux_ρ = zeros(ncellx + 1, ncelly + 1)
    flux_u = zeros(ncellx + 1, ncelly + 1)
    flux_v = zeros(ncellx + 1, ncelly + 1)

    println("Run the simulation...\n")
    ntime = floor(Int, final_time / Δt)
    p = Progress(ntime)
    for itime in 1:ntime
        scheme_iter!(ρ,u,v,Δx,Δy,Δt,c1,c2,λ,bcond_x,bcond_y,Fx,Fy,flux_ρ,flux_u,flux_v,method)
        if should_save
            if itime%save_step == 0
                save_data!(data,ρ,u,v,data_dir,"data_$itime.jld2",key="iter_$itime")
            end
        end
        if should_plot || save_video
            if itime%plot_step == 0
                round_time = nice_float2string(itime*Δt,2)
                if save_video
                    update_plot!(fig,ax,hm,arrows,ρ,u,v,"time=$round_time",plot_dir,"$itime.png",should_plot,stream)
                else
                    update_plot!(fig,ax,hm,arrows,ρ,u,v,"time=$round_time",plot_dir,"$itime.png")
                end
            end
        end
        next!(p)
    end

    #=======================================================================================#
    #============ Collect the simulation time and save the final data ======================#
    #=======================================================================================#

    parameters[:duration] = p.tlast-p.tinit
    save(joinpath(dir_name,"parameters.jld2"), "parameters", parameters)
    if !should_save
        save_data!(data,ρ,u,v,data_dir,"data_$ntime.jld2",key="iter_$ntime")
    end

    if should_plot || save_video
        if save_video
            save(joinpath(dir_name,"video.mp4"), stream)
            return fig, ax, hm, arrows, stream
        else
            return fig, ax, hm, arrows
        end
    end

end


text:
# Tests for Revenue Additive DEA Models
@testset "RevenueAdditiveDEAModel" begin

    # Test using Book data
    X = [1; 1; 1; 1; 1; 1; 1; 1];
    Y = [7 7; 4 8; 8 4; 3 5; 3 3; 8 2; 6 4; 1.5 5];
    P = [1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1];

    # With :Ones
    dearevenueaddones = dearevenueadd(X, Y, P, :Ones)

    @test typeof(dearevenueaddones) == RevenueAdditiveDEAModel

    @test nobs(dearevenueaddones) == 8
    @test ninputs(dearevenueaddones) == 1
    @test noutputs(dearevenueaddones) == 2
    @test ismonetary(dearevenueaddones) == false

    @test efficiency(dearevenueaddones) == efficiency(dearevenueaddones, :Economic)
    @test efficiency(dearevenueaddones, :Economic)   ≈ [0.0; 2.0; 2.0; 6.0; 8.0; 4.0; 4.0; 7.5]
    @test efficiency(dearevenueaddones, :Technical)  ≈ [0.0; 0.0; 0.0; 6.0; 8.0; 2.0; 4.0; 7.5]
    @test efficiency(dearevenueaddones, :Allocative) ≈ [0.0; 2.0; 2.0; 0.0; 0.0; 2.0; 0.0; 0.0]
    @test normfactor(dearevenueaddones) ≈ [1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0]

    # With :Ones and CRS (Same result in this example)
    dearevenueaddonescrs = dearevenueadd(X, Y, P, :Ones, rts = :CRS)

    @test efficiency(dearevenueaddonescrs, :Economic)   ≈ [0.0; 2.0; 2.0; 6.0; 8.0; 4.0; 4.0; 7.5]
    @test efficiency(dearevenueaddonescrs, :Technical)  ≈ [0.0; 0.0; 0.0; 6.0; 8.0; 2.0; 4.0; 7.5]
    @test efficiency(dearevenueaddonescrs, :Allocative) ≈ [0.0; 2.0; 2.0; 0.0; 0.0; 2.0; 0.0; 0.0]

    # Default is :Ones
    @test efficiency(dearevenueadd(X, Y, P)) == efficiency(dearevenueaddones)

    # With :MIP
    dearevenueaddmip = dearevenueadd(X, Y, P, :MIP)

    @test efficiency(dearevenueaddmip, :Economic)   ≈ [0.0; 0.5; 0.5; 2.0; 8/3; 2.0; 1.0; 5.0] atol = 1e-5
    @test efficiency(dearevenueaddmip, :Technical)  ≈ [0.0; 0.0; 0.0; 26/15; 8/3; 1.0; 11/12; 37/9] atol = 1e-5
    @test efficiency(dearevenueaddmip, :Allocative) ≈ [0.0; 0.5; 0.5; 4/15; 0.0; 1.0; 1/12; 8/9] atol = 1e-5
    @test normfactor(dearevenueaddmip) ≈ [7.0; 4.0; 4.0; 3.0; 3.0; 2.0; 4.0; 1.5]

    # With :MIP and monetary = true
    dearevenueaddmipmon = dearevenueadd(X, Y, P, :MIP, monetary = true)

    @test ismonetary(dearevenueaddmipmon) == true

    @test efficiency(dearevenueaddmipmon) ≈ efficiency(dearevenueaddmip) .* normfactor(dearevenueaddmip)
    @test efficiency(dearevenueaddmipmon, :Technical) ≈ efficiency(dearevenueaddmip, :Technical) .* normfactor(dearevenueaddmip)
    @test efficiency(dearevenueaddmipmon, :Allocative) ≈ efficiency(dearevenueaddmip, :Allocative) .* normfactor(dearevenueaddmip)

    # With custom weights
    dearevenueaddcustom = dearevenueadd(X, Y, P, rhoY = 1 ./ Y)

    @test efficiency(dearevenueaddcustom) ≈ efficiency(dearevenueaddmip)

    # Print
    show(IOBuffer(), dearevenueaddones)
    show(IOBuffer(), dearevenueadd(X, Y, P, :Ones, dispos = :Weak))

    # Test errors
    @test_throws DimensionMismatch dearevenueadd([1; 2 ; 3], [4 ; 5], [1; 1; 1] , :Ones) #  Different number of observations
    @test_throws DimensionMismatch dearevenueadd([1; 2; 3], [4; 5; 6], [1; 2; 3; 4], :Ones) # Different number of observation in prices
    @test_throws DimensionMismatch dearevenueadd([1 1; 2 2; 3 3 ], [4; 5; 6], [1 1 1; 2 2 2; 3 3 3], :Ones) # Different number of output prices and Outputs
    @test_throws DimensionMismatch dearevenueadd([1; 2; 3] , [4; 5; 6], [1; 1; 1], rhoY = [1 1; 2 2; 3 3]) # Different size of output weights
    @test_throws ArgumentError dearevenueadd([1; 2; 3], [4; 5; 6], [1; 2; 3], :Ones, rts = :Error) # Invalid returns to scale
    @test_throws ArgumentError dearevenueadd([1; 2; 3], [1; 2; 3], [1; 1; 1], :Error) # Invalid inuts direction
    @test_throws ArgumentError dearevenueadd([1; 2; 3], [1; 2; 3], [1; 1; 1], :Ones, dispos = :Error) # Invalid inputs disposal
    @test_throws ArgumentError dearevenueadd([1; 2; 3], [1; 2; 3], [1; 1; 1], :Ones, rhoY= [1; 1; 1]) # Weights not allowed if model != :Custom

end


text:
"""
    Quaternion{T <: Real} = Elliptic{Binion{T}}

A quaternion is an elliptic Cayley-Dickson construct with binions.
"""
const Quaternion{T <: Real} = Elliptic{Binion{T}}

function Quaternion(z::Binion{T}) where T <: Real
    Quaternion{T}(z, zero(Binion{T}))
end

function Quaternion(a::T, b::T, c::T, d::T) where T <: Real
    Quaternion{T}(Binion{T}(a, b), Binion{T}(c, d))
end

function Quaternion(a::T, b::T, c::T) where T <: Real
    Quaternion{T}(Binion{T}(a, b), Binion{T}(c, zero(T)))
end

function Quaternion(a::T, b::T) where T <: Real
    Quaternion{T}(Binion{T}(a, b), zero(Binion{T}))
end

function Quaternion(a::T) where T <: Real
    Quaternion{T}(Binion{T}(a, zero(T)), zero(Binion{T}))
end

function Quaternion(a::Real, b::Real)
    Quaternion(promote(a, b)...)
end

function Quaternion(a::Real, b::Real, c::Real)
    Quaternion(promote(a, b, c)...)
end

function Quaternion(a::Real, b::Real, c::Real, d::Real)
    Quaternion(promote(a, b, c, d)...)
end

"""
    unreal(z::Quaternion)

The unreal part of a quaternion is a 3-vector.
"""
function unreal(z::Quaternion)
    vcat(unreal(z.l), array(z.r))
end

function show(io::IO, z::Quaternion)
    print(io, "[1: ")
    print(io, z.l.l)
    print(io, ", i: ")
    print(io, z.l.r)
    print(io, ", j: ")
    print(io, z.r.l)
    print(io, ", ij: ")
    print(io, z.r.r)
    print(io, "]")
end

function random(::Type{Quaternion{T}}) where T <: Real
    Quaternion{T}(random(Binion{T}), random(Binion{T}))
end


text:
using Documenter: DocMeta, doctest
using DevEnvConfig
using Test

using ..DevEnvConfig: Success, Warning, Error

const gitconf = DevEnvConfig.GitTools.GitConfig("loc_user", "user@foo.com", "hub_user", "main")

@testset "DevEnvConfig" begin
    @test newpkg("Pk1";dir=mktempdir(),private=false,                  testing_gitconfig=gitconf) != Error
    @test newpkg("Pk2";dir=mktempdir(),private=true,  docrepo="Docs",  testing_gitconfig=gitconf) != Error
end

# Run all doctests in the package. Also run the doctests in docs/src.
# Set manual=false to run only doctest from doc strings in the sources.
# If this is too slow, you can comment it out. The doctests will be run
# anyway by "docs/make.jl" each time the documention is built. If you
# comment this out, a doctest failure will *not* be detected by the CI
# workflow "Tests.yml" since it only executes this file.

# DocMeta.setdocmeta!(DevEnvConfig, :DocTestSetup, :(using DevEnvConfig); recursive=true)
# doctest(DevEnvConfig; manual = true)


text:
###################################################
## network/solvers/roadtypespeed.jl
## A variation on the SOCP solver with a constant speed for each road-type
###################################################

"""
    constantSpeedSolver :
    optimize travel times to minimize SOCP error from data with given paths

    `uniqueSpeed` can be
    - false: one speed per roadType
    - true : one speed for everyone
"""
function constantSpeedSolver(s::IterativeState; uniqueSpeed::Bool = false, args...)
    g = s.data.network.graph
    paths = s.paths
    tripData = s.trips
    roads = s.data.network.roads

    #Create the model (will be changed to avoid hard-coded parameters)
    # !BarConvTol needs to be changed
    m = Model(solver = MosekSolver(MSK_DPAR_OPTIMIZER_MAX_TIME=10000.; args...))

    # DECISION VARIABLES
    # inversed Road-type speeds times
    @variable(m, s[r=1:6] >= 0)

    # Absolute difference between tripData times and computed times
    @variable(m, maxRatioError[d=eachindex(tripData)] >= 0)
    @variable(m, T[d=eachindex(tripData)] >= 0)

    # OBJECTIVE
    @objective(m, Min, sum(tripData[d].weight * maxRatioError[d] for d in eachindex(tripData)))

    # CONSTRAINTS
    # big T constraints
    @constraint(m, pathTime[d=eachindex(tripData)],
        T[d] == sum(paths[d][1][edge] * s[roads[src(edge), dst(edge)].roadType] * roads[src(edge), dst(edge)].distance for edge in keys(paths[d][1])))

    # second order cone constraints (define epsilon), equal to time of first path
    @constraint(m, errorLower[d=eachindex(tripData)],
        norm([2 * sqrt(tripData[d].time), T[d] - maxRatioError[d]])
        <= T[d] + maxRatioError[d]
        )
    @constraint(m, errorUpper[d=eachindex(tripData)],
        T[d] <= maxRatioError[d] * tripData[d].time
        )

    if uniqueSpeed
        @constraint(m, equalSpeeds[i=2:6], s[i] == s[1])
    end

    # SOLVE SOCP
    status = solve(m)

    invSpeeds = getvalue(s)

    return Float64[1/iSp for iSp in invSpeeds]
end


text:
struct YeeGrid{T,nD}
	x::StructuredGrid{nD}
	y::StructuredGrid{nD}
	z::StructuredGrid{nD}
end

function YeeGrid{T}(refgrid::StructuredGrid{nD}) where {nD,T}
	YeeGrid{T,nD}(
		nodes(T,XComponent,refgrid),
		nodes(T,YComponent,refgrid),
		nodes(T,ZComponent,refgrid))
end
function YeeGrid{T}(nodes...) where T
	refgrid = StructuredGrid(nodes...)
	YeeGrid{T}(refgrid)
end

Base.size(yeegrid::YeeGrid) = (size(yeegrid.x),size(yeegrid.y),size(yeegrid.z))


text:
import Base.Docs: doc as getdocs

"""
    get_docstring(obj)

Extract and style an object's docstring.
"""
function get_docstring(obj)
    # get doc and docstring
    doc = getdocs(obj)

    docstring = if isnothing(doc)
        "no docstring"
    else
        highlight(highlight(doc), :docstring)
    end
    return doc, unescape_brackets(docstring)
end

"""
    style_super_types(info)::String

Style a vector of super types
"""
function style_super_types(info)::String
    if !isnothing(info.supertypes)
        stypes = "{dim}(supertypes):{/dim} {bold blue}$(info.name){/bold blue}"
        for sup in info.supertypes[2:end]
            abstract = isabstracttype(sup) ? "underline" : ""
            name = split(string(sup), ".")[end]
            stypes *= "{blue} <{/blue} {bold $abstract}$name{/bold $abstract}"
        end
    else
        stypes = "{dim}(supertypes): no super types{/dim}"
    end
    return stypes
end

"""
    style_sub_types(info)::String

Style a vector of sub types.
"""
function style_sub_types(info)::String
    if !isnothing(info.subtypes)
        subtypes = "  {dim}(subtypes):{/dim} "
        for sub in info.subtypes
            name = split(string(sub), ".")[end]
            subtypes *= "{bold}$name{/bold} ~ "
        end
        subtypes = subtypes[1:(end - 3)]
    else
        subtypes = "  {dim}(subtypes): no subtypes{/dim}"
    end
    return subtypes
end


text:
using Documenter, DynamicGrids

makedocs(
    modules = [DynamicGrids],
    sitename = "DynamicGrids.jl",
)

deploydocs(
    repo = "github.com/cesaraustralia/DynamicGrids.jl.git",
)


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

# What is included in the file `FMI2_int.jl` (internal functions)?
# - optional, more comfortable calls to the C-functions from the FMI-spec (example: `fmiGetReal!(c, v, a)` is bulky, `a = fmiGetReal(c, v)` is more user friendly)

"""
TODO: FMI specification reference.

Returns a string representing the header file used to compile the FMU.

Returns "default" by default.

For more information call ?fmi2GetTypesPlatform
"""
# this is basically already defined in FMI2_c.jl

"""
TODO: FMI specification reference.

Returns the version of the FMI Standard used in this FMU.

For more information call ?fmi2GetVersion
"""
# this is basically already defined in FMI2_c.jl

"""
TODO: FMI specification reference.

Set the DebugLogger for the FMU.
"""
function fmi2SetDebugLogging(c::FMU2Component)
    fmi2SetDebugLogging(c.fmu, c.compAddr, fmi2False, Unsigned(0), C_NULL)
end

"""
TODO: FMI specification reference.

Setup the simulation but without defining all of the parameters.

For more information call ?fmi2SetupExperiment (#ToDo endless recursion)
"""
function fmi2SetupExperiment(c::FMU2Component, startTime::Union{Real, Nothing} = nothing, stopTime::Union{Real, Nothing} = nothing; tolerance::Union{Real, Nothing} = nothing)

    if startTime == nothing
        startTime = fmi2GetDefaultStartTime(c.fmu.modelDescription)
        if startTime == nothing
            startTime = 0.0
        end
    end

    # default stopTime is set automatically if doing nothing
    # if stopTime == nothing
    #     stopTime = fmi2GetDefaultStopTime(c.fmu.modelDescription)
    # end

    # default tolerance is set automatically if doing nothing
    # if tolerance == nothing
    #     tolerance = fmi2GetDefaultTolerance(c.fmu.modelDescription)
    # end

    c.t = startTime

    toleranceDefined = (tolerance != nothing)
    if !toleranceDefined
        tolerance = 0.0 # dummy value, will be ignored
    end

    stopTimeDefined = (stopTime != nothing)
    if !stopTimeDefined
        stopTime = 0.0 # dummy value, will be ignored
    end

    fmi2SetupExperiment(c, fmi2Boolean(toleranceDefined), fmi2Real(tolerance), fmi2Real(startTime), fmi2Boolean(stopTimeDefined), fmi2Real(stopTime))
end

"""
TODO: FMI specification reference.

Get the values of an array of fmi2Real variables.

For more information call ?fmi2GetReal!
"""
function fmi2GetReal(c::FMU2Component, vr::fmi2ValueReferenceFormat)

    vr = prepareValueReference(c, vr)

    nvr = Csize_t(length(vr))
    values = zeros(fmi2Real, nvr)
    fmi2GetReal!(c.fmu.cGetReal, c.compAddr, vr, nvr, values)

    if length(values) == 1
        return values[1]
    else
        return values
    end
end

"""
TODO: FMI specification reference.

Get the values of an array of fmi2Real variables.

For more information call ?fmi2GetReal!
"""
function fmi2GetReal!(c::FMU2Component, vr::fmi2ValueReferenceFormat, values::AbstractArray{fmi2Real})

    vr = prepareValueReference(c, vr)
    # values = prepareValue(values)
    @assert length(vr) == length(values) "fmi2GetReal!(...): `vr` and `values` need to be the same length."

    nvr = Csize_t(length(values))
    # values[:] = fmi2Real.(values)
    fmi2GetReal!(c.fmu.cGetReal, c.compAddr, vr, nvr, values)
    nothing
end
function fmi2GetReal!(c::FMU2Component, vr::fmi2ValueReferenceFormat, values::Real)
    @assert false "fmi2GetReal! is only possible for arrays of values, please use an array instead of a scalar."
end

"""
TODO: FMI specification reference.

Set the values of an array of fmi2Real variables.

For more information call ?fmi2SetReal
"""
function fmi2SetReal(c::FMU2Component, vr::fmi2ValueReferenceFormat, values::Union{AbstractArray{<:Real}, <:Real})

    vr = prepareValueReference(c, vr)
    values = prepareValue(values)
    @assert length(vr) == length(values) "fmi2SetReal(...): `vr` ($(length(vr))) and `values` ($(length(values))) need to be the same length."

    nvr = Csize_t(length(vr))
    fmi2SetReal(c.fmu.cSetReal, c.compAddr, vr, nvr, Array{fmi2Real}(values))
end

"""
TODO: FMI specification reference.

Get the values of an array of fmi2Integer variables.

For more information call ?fmi2GetInteger!
"""
function fmi2GetInteger(c::FMU2Component, vr::fmi2ValueReferenceFormat)

    vr = prepareValueReference(c, vr)

    nvr = Csize_t(length(vr))
    values = zeros(fmi2Integer, nvr)
    fmi2GetInteger!(c.fmu.cGetInteger, c.compAddr, vr, nvr, values)

    if length(values) == 1
        return values[1]
    else
        return values
    end
end

"""
TODO: FMI specification reference.

Get the values of an array of fmi2Integer variables.

For more information call ?fmi2GetInteger!
"""
function fmi2GetInteger!(c::FMU2Component, vr::fmi2ValueReferenceFormat, values::AbstractArray{fmi2Integer})

    vr = prepareValueReference(c, vr)
    # values = prepareValue(values)
    @assert length(vr) == length(values) "fmi2GetInteger!(...): `vr` and `values` need to be the same length."

    nvr = Csize_t(length(values))
    fmi2GetInteger!(c.fmu.cGetInteger, c.compAddr, vr, nvr, values)
    nothing
end
function fmi2GetInteger!(c::FMU2Component, vr::fmi2ValueReferenceFormat, values::Integer)
    @assert false "fmi2GetInteger! is only possible for arrays of values, please use an array instead of a scalar."
end

"""
TODO: FMI specification reference.

Set the values of an array of fmi2Integer variables.

For more information call ?fmi2SetInteger
"""
function fmi2SetInteger(c::FMU2Component, vr::fmi2ValueReferenceFormat, values::Union{AbstractArray{<:Integer}, <:Integer})

    vr = prepareValueReference(c, vr)
    values = prepareValue(values)
    @assert length(vr) == length(values) "fmi2SetInteger(...): `vr` and `values` need to be the same length."

    nvr = Csize_t(length(vr))
    fmi2SetInteger(c.fmu.cSetInteger, c.compAddr, vr, nvr, Array{fmi2Integer}(values))
end

"""
TODO: FMI specification reference.

Get the values of an array of fmi2Boolean variables.

For more information call ?fmi2GetBoolean!
"""
function fmi2GetBoolean(c::FMU2Component, vr::fmi2ValueReferenceFormat)

    vr = prepareValueReference(c, vr)

    nvr = Csize_t(length(vr))
    values = Array{fmi2Boolean}(undef, nvr)
    fmi2GetBoolean!(c.fmu.cGetBoolean, c.compAddr, vr, nvr, values)

    if length(values) == 1
        return values[1]
    else
        return values
    end
end

"""
TODO: FMI specification reference.

Get the values of an array of fmi2Boolean variables.

For more information call ?fmi2GetBoolean!
"""
function fmi2GetBoolean!(c::FMU2Component, vr::fmi2ValueReferenceFormat, values::AbstractArray{fmi2Boolean})

    vr = prepareValueReference(c, vr)
    # values = prepareValue(values)
    @assert length(vr) == length(values) "fmi2GetBoolean!(...): `vr` and `values` need to be the same length."

    nvr = Csize_t(length(values))
    #values = fmi2Boolean.(values)
    fmi2GetBoolean!(c.fmu.cGetBoolean, c.compAddr, vr, nvr, values)

    nothing
end
function fmi2GetBoolean!(c::FMU2Component, vr::fmi2ValueReferenceFormat, values::Bool)
    @assert false "fmi2GetBoolean! is only possible for arrays of values, please use an array instead of a scalar."
end

"""
TODO: FMI specification reference.

Set the values of an array of fmi2Boolean variables.

For more information call ?fmi2SetBoolean
"""
function fmi2SetBoolean(c::FMU2Component, vr::fmi2ValueReferenceFormat, values::Union{AbstractArray{Bool}, Bool})

    vr = prepareValueReference(c, vr)
    values = prepareValue(values)
    @assert length(vr) == length(values) "fmi2SetBoolean(...): `vr` and `values` need to be the same length."

    nvr = Csize_t(length(vr))
    fmi2SetBoolean(c.fmu.cSetBoolean, c.compAddr, vr, nvr, Array{fmi2Boolean}(values))
end

"""
TODO: FMI specification reference.

Get the values of an array of fmi2String variables.

For more information call ?fmi2GetString!
"""
function fmi2GetString(c::FMU2Component, vr::fmi2ValueReferenceFormat)

    vr = prepareValueReference(c, vr)

    nvr = Csize_t(length(vr))
    vars = Vector{fmi2String}(undef, nvr)
    values = string.(zeros(nvr))
    fmi2GetString!(c.fmu.cGetString, c.compAddr, vr, nvr, vars)
    values[:] = unsafe_string.(vars)

    if length(values) == 1
        return values[1]
    else
        return values
    end
end

"""
TODO: FMI specification reference.

Get the values of an array of fmi2String variables.

For more information call ?fmi2GetString!
"""
function fmi2GetString!(c::FMU2Component, vr::fmi2ValueReferenceFormat, values::AbstractArray{fmi2String})

    vr = prepareValueReference(c, vr)
    @assert length(vr) == length(values) "fmi2GetString!(...): `vr` and `values` need to be the same length."

    nvr = Csize_t(length(vr))
    fmi2GetString!(c.fmu.cGetString, c.compAddr, vr, nvr, values)

    nothing
end
function fmi2GetString!(c::FMU2Component, vr::fmi2ValueReferenceFormat, values::String)
    @assert false "fmi2GetString! is only possible for arrays of values, please use an array instead of a scalar."
end

"""
TODO: FMI specification reference.

Set the values of an array of fmi2String variables.

For more information call ?fmi2SetString
"""
function fmi2SetString(c::FMU2Component, vr::fmi2ValueReferenceFormat, values::Union{AbstractArray{String}, String})

    vr = prepareValueReference(c, vr)
    values = prepareValue(values)
    @assert length(vr) == length(values) "fmi2SetReal(...): `vr` and `values` need to be the same length."

    nvr = Csize_t(length(vr))
    ptrs = pointer.(values)
    fmi2SetString(c.fmu.cSetString, c.compAddr, vr, nvr, ptrs)
end

"""
TODO: FMI specification reference.

Get the pointer to the current FMU state.

For more information call ?fmi2GetFMUstate
"""
function fmi2GetFMUstate(c::FMU2Component)
    state = fmi2FMUstate()
    stateRef = Ref(state)
    fmi2GetFMUstate!(c.fmu.cGetFMUstate, c.compAddr, stateRef)
    state = stateRef[]
    state
end

"""
TODO: FMI specification reference.

Free the allocated memory for the FMU state.

For more information call ?fmi2FreeFMUstate
"""
function fmi2FreeFMUstate!(c::FMU2Component, state::fmi2FMUstate)
    stateRef = Ref(state)
    fmi2FreeFMUstate!(c.fmu.cFreeFMUstate, c.compAddr, stateRef)
    state = stateRef[]
    return nothing
end

"""
TODO: FMI specification reference.

Returns the size of a byte vector the FMU can be stored in.

For more information call ?fmi2SerzializedFMUstateSize
"""
function fmi2SerializedFMUstateSize(c::FMU2Component, state::fmi2FMUstate)
    size = 0
    sizeRef = Ref(Csize_t(size))
    fmi2SerializedFMUstateSize!(c.fmu.cSerializedFMUstateSize, c.compAddr, state, sizeRef)
    size = sizeRef[]
end

"""
TODO: FMI specification reference.

Serialize the data in the FMU state pointer.

For more information call ?fmi2SerzializeFMUstate
"""
function fmi2SerializeFMUstate(c::FMU2Component, state::fmi2FMUstate)
    size = fmi2SerializedFMUstateSize(c, state)
    serializedState = Array{fmi2Byte}(undef, size)
    status = fmi2SerializeFMUstate!(c.fmu.cSerializeFMUstate, c.compAddr, state, serializedState, size)
    @assert status == Int(fmi2StatusOK) ["Failed with status `$status`."]
    serializedState
end

"""
TODO: FMI specification reference.

Deserialize the data in the serializedState fmi2Byte field.

For more information call ?fmi2DeSerzializeFMUstate
"""
function fmi2DeSerializeFMUstate(c::FMU2Component, serializedState::AbstractArray{fmi2Byte})
    size = length(serializedState)
    state = fmi2FMUstate()
    stateRef = Ref(state)

    status = fmi2DeSerializeFMUstate!(c.fmu.cDeSerializeFMUstate, c.compAddr, serializedState, Csize_t(size), stateRef)
    @assert status == Int(fmi2StatusOK) "Failed with status `$status`."

    state = stateRef[]
end

"""
TODO: FMI specification reference.

Computes directional derivatives.

For more information call ?fmi2GetDirectionalDerivatives
"""
function fmi2GetDirectionalDerivative(c::FMU2Component,
                                      vUnknown_ref::AbstractArray{fmi2ValueReference},
                                      vKnown_ref::AbstractArray{fmi2ValueReference},
                                      dvKnown::Union{AbstractArray{fmi2Real}, Nothing} = nothing)

    nUnknown = Csize_t(length(vUnknown_ref))

    dvUnknown = zeros(fmi2Real, nUnknown)
    status = fmi2GetDirectionalDerivative!(c, vUnknown_ref, vKnown_ref, dvUnknown, dvKnown)
    @assert status == fmi2StatusOK ["Failed with status `$status`."]

    return dvUnknown
end

"""
TODO: FMI specification reference.

Computes directional derivatives.

For more information call ?fmi2GetDirectionalDerivatives
"""
function fmi2GetDirectionalDerivative!(c::FMU2Component,
                                      vUnknown_ref::AbstractArray{fmi2ValueReference},
                                      vKnown_ref::AbstractArray{fmi2ValueReference},
                                      dvUnknown::AbstractArray, # ToDo: Data-type
                                      dvKnown::Union{AbstractArray{fmi2Real}, Nothing} = nothing)

    nKnown = Csize_t(length(vKnown_ref))
    nUnknown = Csize_t(length(vUnknown_ref))

    if dvKnown == nothing
        dvKnown = ones(fmi2Real, nKnown)
    end

    status = fmi2GetDirectionalDerivative!(c.fmu.cGetDirectionalDerivative, c.compAddr, vUnknown_ref, nUnknown, vKnown_ref, nKnown, dvKnown, dvUnknown)

    return status
end

"""
TODO: FMI specification reference.

Computes directional derivatives.

For more information call ?fmi2GetDirectionalDerivatives
"""
function fmi2GetDirectionalDerivative(c::FMU2Component,
                                      vUnknown_ref::fmi2ValueReference,
                                      vKnown_ref::fmi2ValueReference,
                                      dvKnown::fmi2Real = 1.0)

    fmi2GetDirectionalDerivative(c, [vUnknown_ref], [vKnown_ref], [dvKnown])[1]
end

# CoSimulation specific functions
"""
TODO: FMI specification reference.

Sets the n-th time derivative of real input variables.
vr defines the value references of the variables
the array order specifies the corresponding order of derivation of the variables
"""
function fmi2SetRealInputDerivatives(c::FMU2Component, vr::fmi2ValueReferenceFormat, order::AbstractArray{fmi2Integer}, values::AbstractArray{fmi2Real})
    vr = prepareValueReference(c, vr)
    order = prepareValue(order)
    values = prepareValue(values)
    nvr = Csize_t(length(vr))
    fmi2SetRealInputDerivatives(c.fmu.cSetRealInputDerivatives, c.compAddr, vr, nvr, order, values)
end

"""
TODO: FMI specification reference.

vr defines the value references of the variables
the array order specifies the corresponding order of derivation of the variables

For more information call ?fmi2GetRealOutputDerivatives
"""
function fmi2GetRealOutputDerivatives(c::FMU2Component, vr::fmi2ValueReferenceFormat, order::AbstractArray{fmi2Integer})
    vr = prepareValueReference(c, vr)
    order = prepareValue(order)
    nvr = Csize_t(length(vr))
    values = zeros(fmi2Real, nvr)
    fmi2GetRealOutputDerivatives!(c.fmu.cGetRealOutputDerivatives, c.compAddr, vr, nvr, order, values)

    if length(values) == 1
        return values[1]
    else
        return values
    end
end

"""
TODO: FMI specification reference.

The computation of a time step is started.

For more information call ?fmi2DoStep
"""
function fmi2DoStep(c::FMU2Component, communicationStepSize::Union{Real, Nothing} = nothing; currentCommunicationPoint::Union{Real, Nothing} = nothing, noSetFMUStatePriorToCurrentPoint::Bool = true)

    # skip `fmi2DoStep` if this is set (allows evaluation of a CS_NeuralFMUs at t_0)
    if c.skipNextDoStep
        c.skipNextDoStep = false
        return fmi2StatusOK
    end

    if currentCommunicationPoint == nothing
        currentCommunicationPoint = c.t
    end

    if communicationStepSize == nothing
        communicationStepSize = fmi2GetDefaultStepSize(c.fmu.modelDescription)
        if communicationStepSize == nothing
            communicationStepSize = 1e-2
        end
    end

    c.t = currentCommunicationPoint
    status = fmi2DoStep(c.fmu.cDoStep, c.compAddr, fmi2Real(currentCommunicationPoint), fmi2Real(communicationStepSize), fmi2Boolean(noSetFMUStatePriorToCurrentPoint))
    c.t += communicationStepSize

    return status
end

"""
TODO: FMI specification reference.
"""
function fmi2SetTime(c::FMU2Component, t::Real)
    status = fmi2SetTime(c.fmu.cSetTime, c.compAddr, fmi2Real(t))
    c.t = t
    return status
end

# Model Exchange specific functions

"""
TODO: FMI specification reference.

Set a new (continuous) state vector and reinitialize chaching of variables that depend on states.

For more information call ?fmi2SetContinuousStates
"""
function fmi2SetContinuousStates(c::FMU2Component, x::Union{AbstractArray{Float32}, AbstractArray{Float64}})
    nx = Csize_t(length(x))
    status = fmi2SetContinuousStates(c.fmu.cSetContinuousStates, c.compAddr, Array{fmi2Real}(x), nx)
    if status == fmi2StatusOK
        c.x = x
    end
    return status
end

"""
TODO: FMI specification reference.

Increment the super dense time in event mode.

For more information call ?fmi2NewDiscretestates
"""
function fmi2NewDiscreteStates(c::FMU2Component)
    eventInfo = fmi2EventInfo()
    fmi2NewDiscreteStates!(c, eventInfo)
    eventInfo
end

"""
TODO: FMI specification reference.

This function must be called by the environment after every completed step
If enterEventMode == fmi2True, the event mode must be entered
If terminateSimulation == fmi2True, the simulation shall be terminated

For more information call ?fmi2CompletedIntegratorStep
"""
function fmi2CompletedIntegratorStep(c::FMU2Component,
                                     noSetFMUStatePriorToCurrentPoint::fmi2Boolean)
    enterEventMode = zeros(fmi2Boolean, 1)
    terminateSimulation = zeros(fmi2Boolean, 1)

    status = fmi2CompletedIntegratorStep!(c.fmu.cCompletedIntegratorStep,
                                          c.compAddr,
                                          noSetFMUStatePriorToCurrentPoint,
                                          pointer(enterEventMode),
                                          pointer(terminateSimulation))

    return (status, enterEventMode[1], terminateSimulation[1])
end

"""
TODO: FMI specification reference.

Compute state derivatives at the current time instant and for the current states.

For more information call ?fmi2GetDerivatives
"""
function fmi2GetDerivatives(c::FMU2Component)
    nx = Csize_t(length(c.fmu.modelDescription.stateValueReferences))
    derivatives = zeros(fmi2Real, nx)
    fmi2GetDerivatives!(c, derivatives)
    return derivatives
end

"""
TODO: FMI specification reference.

Compute state derivatives at the current time instant and for the current states.

For more information call ?fmi2GetDerivatives
"""
function fmi2GetDerivatives!(c::FMU2Component, derivatives::AbstractArray{fmi2Real})
    status = fmi2GetDerivatives!(c, derivatives, Csize_t(length(derivatives)))
    if status == fmi2StatusOK
        c.ẋ = derivatives
    end
    return status
end

"""
TODO: FMI specification reference.

Returns the event indicators of the FMU.

For more information call ?fmi2GetEventIndicators
"""
function fmi2GetEventIndicators(c::FMU2Component)
    ni = Csize_t(c.fmu.modelDescription.numberOfEventIndicators)
    eventIndicators = zeros(fmi2Real, ni)
    fmi2GetEventIndicators!(c, eventIndicators)
    return eventIndicators
end

"""
TODO: FMI specification reference.

Returns the event indicators of the FMU.

For more information call ?fmi2GetEventIndicators
"""
function fmi2GetEventIndicators!(c::FMU2Component, eventIndicators::AbstractArray{fmi2Real})
    ni = Csize_t(length(eventIndicators))
    fmi2GetEventIndicators!(c.fmu.cGetEventIndicators, c.compAddr, eventIndicators, ni)
end

"""
TODO: FMI specification reference.

Return the new (continuous) state vector x.

For more information call ?fmi2GetContinuousStates
"""
function fmi2GetContinuousStates(c::FMU2Component)
    nx = Csize_t(length(c.fmu.modelDescription.stateValueReferences))
    x = zeros(fmi2Real, nx)
    fmi2GetContinuousStates!(c.fmu.cGetContinuousStates, c.compAddr, x, nx)
    x
end

"""
TODO: FMI specification reference.

Return the new (continuous) state vector x.

For more information call ?fmi2GetNominalsOfContinuousStates
"""
function fmi2GetNominalsOfContinuousStates(c::FMU2Component)
    nx = Csize_t(length(c.fmu.modelDescription.stateValueReferences))
    x = zeros(fmi2Real, nx)
    fmi2GetNominalsOfContinuousStates!(c.fmu.cGetNominalsOfContinuousStates, c.compAddr, x, nx)
    x
end

"""
ToDo
"""
function fmi2GetStatus(c::FMU2Component, s::fmi2StatusKind)
    rtype = nothing
    if s == fmi2Terminated
        rtype = fmi2Boolean
    else
        @assert false "fmi2GetStatus(_, $(s)): StatusKind $(s) not implemented yet, please open an issue."
    end
    value = zeros(rtype, 1)

    status = fmi2Error
    if rtype == fmi2Boolean
        status = fmi2GetStatus!(c.fmu.cGetStatus, c.compAddr, s, value)
    end

    status, value[1]
end

"""
Returns the start/default value for a given value reference.

TODO: Add this command in the documentation.
"""
function fmi2GetStartValue(c::FMU2Component, vrs::fmi2ValueReferenceFormat)

    vrs = prepareValueReference(c, vrs)

    starts = []

    for vr in vrs
        mvs = fmi2ModelVariablesForValueReference(c.fmu.modelDescription, vr)

        if length(mvs) == 0
            @warn "fmi2GetStartValue(...): Found no model variable with value reference $(vr)."
        end

        if mvs[1]._Real != nothing
            push!(starts, mvs[1]._Real.start)
        elseif mvs[1]._Integer != nothing
            push!(starts, mvs[1]._Integer.start)
        elseif mvs[1]._Boolean != nothing
            push!(starts, mvs[1]._Boolean.start)
        elseif mvs[1]._String != nothing
            push!(starts, mvs[1]._String.start)
        elseif mvs[1]._Enumeration != nothing
            push!(starts, mvs[1]._Enumeration.start)
        else
            @assert false "fmi2GetStartValue(...): Value reference $(vr) has no data type."
        end
    end

    if length(vrs) == 1
        return starts[1]
    else
        return starts
    end
end

# further functions

"""
This function samples the directional derivative by manipulating corresponding values (central differences).
"""
function fmi2SampleDirectionalDerivative(c::FMU2Component,
                                       vUnknown_ref::Array{fmi2ValueReference},
                                       vKnown_ref::Array{fmi2ValueReference},
                                       steps::Array{fmi2Real} = ones(fmi2Real, length(vKnown_ref)).*1e-5)

    dvUnknown = zeros(fmi2Real, length(vUnknown_ref), length(vKnown_ref))

    fmi2SampleDirectionalDerivative!(c, vUnknown_ref, vKnown_ref, dvUnknown, steps)

    dvUnknown
end

# TODO should probably in FMI2_ext.jl
"""
This function samples the directional derivative by manipulating corresponding values (central differences) and saves in-place.
"""
function fmi2SampleDirectionalDerivative!(c::FMU2Component,
                                          vUnknown_ref::Array{fmi2ValueReference},
                                          vKnown_ref::Array{fmi2ValueReference},
                                          dvUnknown::AbstractArray,
                                          steps::Array{fmi2Real} = ones(fmi2Real, length(vKnown_ref)).*1e-5)

    for i in 1:length(vKnown_ref)
        vKnown = vKnown_ref[i]
        origValue = fmi2GetReal(c, vKnown)

        fmi2SetReal(c, vKnown, origValue - steps[i]*0.5)
        negValues = fmi2GetReal(c, vUnknown_ref)

        fmi2SetReal(c, vKnown, origValue + steps[i]*0.5)
        posValues = fmi2GetReal(c, vUnknown_ref)

        fmi2SetReal(c, vKnown, origValue)

        if length(vUnknown_ref) == 1
            dvUnknown[1,i] = (posValues-negValues) ./ steps[i]
        else
            dvUnknown[:,i] = (posValues-negValues) ./ steps[i]
        end
    end

    nothing
end


text:
abstract type PSFMode end
struct ModeWidefield <: PSFMode end
struct ModeConfocal <: PSFMode end
struct Mode4Pi <: PSFMode end
struct STED <: PSFMode end

abstract type PSFMethod end
struct MethodSincR <: PSFMethod end
struct MethodPropagate <: PSFMethod end
struct MethodPropagateIterative<: PSFMethod end
struct MethodRichardsWolf <: PSFMethod end
struct MethodShell <: PSFMethod end
export Aberrations
export Zernike_Piston, Zernike_Tilt, Zernike_Tip, Zernike_ObliqueAstigmatism, Zernike_Defocus, Zernike_VerticalAstigmatism, Zernike_VerticalTrefoil
export Zernike_VerticalComa, Zernike_HorizontalComa, Zernike_ObliqueTrefoil, Zernike_ObliqueQuadrafoil, Zernike_ObliqueSecondaryAstigmatism, Zernike_Spherical
export Zernike_VerticalSecondaryAstigmatism, Zernike_VerticalQuadrafoil

const Zernike_Piston = 0
const Zernike_Tilt = 1
const Zernike_Tip = 2
const Zernike_ObliqueAstigmatism = 3
const Zernike_Defocus = 4
const Zernike_VerticalAstigmatism = 5
const Zernike_VerticalTrefoil = 6
const Zernike_VerticalComa = 7
const Zernike_HorizontalComa= 8
const Zernike_ObliqueTrefoil = 9
const Zernike_ObliqueQuadrafoil = 10
const Zernike_ObliqueSecondaryAstigmatism = 11
const Zernike_Spherical = 12
const Zernike_VerticalSecondaryAstigmatism = 13
const Zernike_VerticalQuadrafoil = 14

"""
    Aberrations
    defining Zernike phase aberrations via a list of indices and coefficients and an indexing style.
"""
struct Aberrations
    indices # zernike indices as a vector
    coefficients # corresponding zernike coefficients
    index_style
end

"""
    Aberrations(indices=[],coefficients=[];index_style = :OSA)

defining Zernike phase aberrations via a list of indices and coefficients and an indexing style.
By default, no Aberrations are defined.

Arguments:
+ indices:  Vector of indices
+ coefficients: Vector of corresponding coefficients
+ index_style: type of indexing used. By defaul :OSA is used (See: https://en.wikipedia.org/wiki/Zernike_polynomials#OSA/ANSI_standard_indices)

Here is a list of constants defining the main indices (OSA style):
`:Zernike_Piston` = 0
`:Zernike_Tilt` = 1
`:Zernike_Tip` = 2
`:Zernike_ObliqueAstigmatism` = 3
`:Zernike_Defocus` = 4
`:Zernike_VerticalAstigmatism` = 5
`:Zernike_VerticalTrefoil` = 6
`:Zernike_VerticalComa` = 7
`:Zernike_HorizontalComa` = 8
`:Zernike_ObliqueTrefoil` = 0
`:Zernike_ObliqueQuadrafoil` = 10
`:Zernike_ObliqueSecondaryAstigmatism` = 11
`:Zernike_Spherical` = 12
`:Zernike_VerticalSecondaryAstigmatism` = 13
`:Zernike_VerticalQuadrafoil` = 14

See also:
+ package `ZernikePolynomials.jl`, and https://en.wikipedia.org/wiki/Zernike_polynomials
"""
function Aberrations(indices=[],coefficients=[];index_style = :OSA)
    Aberrations(indices,coefficients,index_style)
end

"""
    PSFParams

This structure stores all the general parameters needed to calculate the point spread function.
Only pixel-pitch and image size information is handles seperately.

Structure Members:
+ λ:            Vacuum wavelength
+ NA:           numerical aperture
+ n:            refractive index of the embedding AND immersion medium
+ dtype:        real-valued data type to generate PSF for
+ mode:         microscopy mode to calculate PSF for. See the constructor PSFParams() for more details.
+ polarization: a function calculating the polarization from a given Tuple of relative-k pupil vectors
+ aplanatic:    aplanatic factor. Provided as a function of angle θ
+ method:       the method of calculation
+ FFTPlan:      information on how to calculate the FFTW plan

See also:
+ PSFParams()

Example:
```jdoctest
julia> using FFTW, PSFs

julia> ppm = PSFParams(0.580, 1.4, 1.518;pol=pol_scalar,method=PSFs.MethodSincR, aberrations= aberr, FFTPlan=FFTW.MEASURE)
PSFParams(0.58, 1.4, 1.518, Float32, ModeWidefield, PSFs.pol_scalar, PSFs.var"#42#43"(), PSFs.MethodSincR, 0x00000000, nothing, nothing)

```
"""
struct PSFParams
    λ  # Vacuum wavelength
    NA # numerical aperture
    n # refractive index of the embedding AND immersion medium
    dtype # real-valued data type to generate PSF for
    mode # microscopy mode (widefield, confocal, etc.) to calculate PSF for ::PSFMode
    polarization # a function calculating the polarization from a given Tuple of relative-k pupil vectors. Can also be pol_scalar()
    aplanatic # aplanatic factor. Provided as a function of angle θ
    method # the method of calculation (e.g. PSFs.MethodPropagate).
    FFTPlan # if not nothing this will be the plan of how to perform FFTs
    aberrations
    pixelshape # here functions can be supplied, which account for the pixelshape influence (itegration over pixel size).
end

"""
    PSFParams(λ=0.5, NA=1.2, n=1.33; pol=pol_scalar, dtype=Float32, mode=ModeWidefield,
    aplanatic = aplanatic_detection, method=MethodPropagateIterative, FFTPlan=nothing,
    aberrations=Aberrations(), pixelshape=nothing)

creates the PSFParams structure via this constructor.

Arguments:
+ λ:            Vacuum wavelength (same units as sampling, default is 0.5 µm)
+ NA:           numerical aperture
+ n:            refractive index of the embedding AND immersion medium
+ pol: a function calculating the polarization from a given Tuple of relative-k pupil vectors
+ dtype:        real-valued data type to generate PSF for
+ mode:         microscopy mode to calculate PSF for ::PSFMode. Currently only `ModeWidefield` is implemented
+ method:         microscopy mode to calculate PSF for
                valid options are currently:
                + MethodPropagate: Angulare spectrum propagation. This version does NOT account for wrap around problems yielding problems at larger out-of-focus distances
                + MethodPropagateIterativ (default): Angulare spectrum propagation accounting from wrap-around problems in each propagation step by applying a perfectly matched layer (PML).
                + MethodShell: Angulare spectrum propagation with a slightly different calculation order. This version does NOT account for wrap around problems yielding problems at larger out-of-focus distances
                + MethodSincR: Based on first calculating a SincR function in real space and applying consecutive filtering steps. It accounts for wrap around problems but requires a quite high sampling along the Z direction.
+ aplanatic:    aplanatic factor. Provided as a function of angle θ.
+ FFTPlan:      information on how to calculate the FFTW plan. Default: nothing (using FFTW.ESTIMATE)

Example:
```jdoctest
julia> using FFTW, PSFs

julia> ppm = PSFParams(0.58, 1.4, 1.518;pol=pol_scalar,method=PSFs.MethodSincR, aberrations= aberr, FFTPlan=FFTW.MEASURE)
PSFParams(0.58, 1.4, 1.518, Float32, ModeWidefield, PSFs.pol_scalar, PSFs.var"#42#43"(), PSFs.MethodSincR, 0x00000000, nothing, nothing)

```
"""
function PSFParams(my_λ=0.5, my_NA=1.2, my_n=1.33; pol=pol_scalar, dtype=Float32, mode=ModeWidefield,
                    aplanatic = aplanatic_detection, method=MethodPropagateIterative, FFTPlan=nothing,
                    aberrations=Aberrations(), pixelshape=nothing)
    PSFParams(my_λ, my_NA, my_n, dtype, mode, pol, aplanatic, method, FFTPlan, aberrations, pixelshape)
end


text:
using Documenter
using Effects

makedocs(modules=[Effects],
         authors="Beacon Biosignals, Inc.",
         repo="https://github.com/beacon-biosignals/Effects.jl/blob/{commit}{path}#{line}",
         sitename="Effects.jl",
         format=Documenter.HTML(prettyurls=get(ENV, "CI", "false") == "true",
                                canonical="https://beacon-biosignals.github.io/Effects.jl",
                                assets=String[]),
         pages=["Home" => "index.md",
                "API" => "api.md"])

deploydocs(; repo="github.com/beacon-biosignals/Effects.jl",
           devbranch="main",
           push_preview=true)