text:
using DataFrames

# Create toy data ---------------------------------
# const COUNTRIES = ["CZE", "DEU", "HUN", "IDN"]
# const PRODUCTS = ["1602", "1604", "1605", "1701", "1702", "1704", "1805", "1806", "1901", "1902"]

POP = [10474410;
       81776930;
       10000023;
       241834215]

POP = convert(Array{Float64, 1}, POP) # Some pop-models might have "half people"

EXP =  [5.66868e7 1.40287e7 1.76446e6 1.24263e8 4.21838e6 1.54845e8 3.0603e6 1.56225e8 3.14846e7 2.60488e7; # CZE
	1.231e9 6.55645e8 7.13669e7 8.70726e8 3.60293e8 8.40892e8 2.75906e8 3.22708e9 1.16899e9 1.82894e8; # DEU
	1.2379e8 1.00424e6 6.27152e6 1.3462e8 7.32589e7 7.9118e7 1.37126e6 7.98984e7 3.50975e7 3.30034e7; # HUN
	5.3835e5 2.63427e8 3.20992e8 1.91463e6 1.57506e7 9.4886e7 1.10084e8 4.17582e7 5.15586e7 1.32186e8] # IDN

# Test RCA ----------------------------------------
@testset "RCA tests" begin
	# Check inputs:
	@test_throws MethodError get_rca(DataFrame(EXP))
	@test_throws MethodError get_rca(EXP[1, :])
	@test_throws MethodError get_rca(["1" "2";
					  "3" "4"])

	# Check output types:
	@test typeof(get_rca(EXP, binary=true)) == Array{Int64, 2}
	@test typeof(get_rca(EXP, binary=false)) == Array{Float64, 2}

	# Check output result:
	@test get_rca(EXP, binary=false)[4, 3] == (EXP[4, 3] / sum(EXP[4,:])) / (sum(EXP[:,3]) / sum(EXP))
	@test get_rca(EXP, binary=false)[3, 2] == (EXP[3, 2] / sum(EXP[3,:])) / (sum(EXP[:,2]) / sum(EXP))

	@test get_rca(EXP, binary=true)[4, 3] == 1
	@test get_rca(EXP, binary=true)[3, 2] == 0
end

# Test RpcA ---------------------------------------
@testset "RPCA tests" begin
	# Check inputs:
	@test_throws MethodError get_rpca(EXP, convert(Array{Int64, 1}, POP))
	@test_throws MethodError get_rpca(DataFrame(EXP), POP)
	@test_throws MethodError get_rpca(EXP[1, :], POP)
	@test_throws MethodError get_rpca(["1" "2"
					   "3" "4"], POP)

	# Check output types:
	@test typeof(get_rpca(EXP, POP, binary=true)) == Array{Int64, 2}
	@test typeof(get_rpca(EXP, POP, binary=false)) == Array{Float64, 2}

	# Check output result:
	@test get_rpca(EXP, POP, binary=false)[4, 3] == (EXP[4, 3] / POP[4]) / (sum(EXP[:, 3]) / sum(POP))
	@test get_rpca(EXP, POP, binary=false)[3, 2] == (EXP[3, 2] / POP[3]) / (sum(EXP[:, 2]) / sum(POP))

	@test get_rpca(EXP, POP, binary=true)[4, 3] == 1
	@test get_rpca(EXP, POP, binary=true)[3, 2] == 0
end

max_stars_repo_name: sbpost/compost

max_issues_repo_name: sbpost/compost

max_forks_repo_name: sbpost/compost


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

import Zygote
using Zygote: @adjoint

using Flux, DiffEqFlux
using OrdinaryDiffEq
using DiffEqCallbacks

include("FMI2_neural.jl")

"""
Structure definition for a NeuralFMU, that runs in mode `Model Exchange` (ME).
"""
mutable struct ME_NeuralFMU
    neuralODE::NeuralODE
    solution::ODESolution
    fmu::FMU

    tspan
    saveat
    saved_values
    recordValues

    ME_NeuralFMU() = new()
end

"""
Structure definition for a NeuralFMU, that runs in mode `Co-Simulation` (CS).
"""
mutable struct CS_NeuralFMU
    model
    #simulationResult::fmi2SimulationResult
    fmu::FMU

    tspan
    saveat
    valueStack

    CS_NeuralFMU() = new()
end

NeuralFMU = ME_NeuralFMU
NeuralFMUs = Union{ME_NeuralFMU, CS_NeuralFMU}

# time caching (to set correct time in ME-NeuralFMUs)
function NeuralFMUGetCachedTime(fmu)
    Zygote.ignore() do
        return fmu.next_t
    end

    return 0.0
end
function NeuralFMUCacheTime(fmu, t)
    fmu.next_t = t
    nothing
end
function NeuralFMUCacheTime_Gradient(c̄, fmu, t)
    tuple(0.0, 0.0)
end
@adjoint NeuralFMUCacheTime(fmu, t) = NeuralFMUCacheTime(fmu, t), c̄ -> NeuralFMUCacheTime_Gradient(c̄, fmu, t)

# state caching (to set correct state in ME-NeuralFMUs)
function NeuralFMUCacheState(fmu, x)
    fmu.next_x = x
end
function NeuralFMUCacheState_Gradient(c̄, fmu, x)
    tuple(0.0, 0.0)
end
@adjoint NeuralFMUCacheState(fmu, x) = NeuralFMUCacheState(fmu, x), c̄ -> NeuralFMUCacheState_Gradient(c̄, fmu, x)

# state derivative caching (to set correct state derivative in ME-NeuralFMUs)
function NeuralFMUCacheStateDerivative(fmu, dx)
    fmu.next_dx = dx
end
function NeuralFMUCacheStateDerivative_Gradient(c̄, fmu, dx)
    tuple(0.0, 0.0)
end
@adjoint NeuralFMUCacheStateDerivative(fmu, dx) = NeuralFMUCacheStateDerivative(fmu, dx), c̄ -> NeuralFMUCacheStateDerivative_Gradient(c̄, fmu, dx)

# helper to add an additional time state later used to setup time inside the FMU
function NeuralFMUInputLayer(fmu, inputs)
    t = inputs[1]
    x = inputs[2:end]
    NeuralFMUCacheTime(fmu, t)
    x
end

# helper to add an additional time state derivative (for ME-NeuralFMUs)
function NeuralFMUOutputLayerME(inputs)
    dt = 1.0
    dx = inputs
    vcat([dt], dx) # [dt, dx...]
end

# helper to add an additional time state derivative (for CS-NeuralFMUs)
function NeuralFMUOutputLayerCS(fmu, inputs)
    out = inputs
    t = fmu.t # NeuralFMUGetCachedTime(fmu)
    vcat([t], out)
end

# function that saves periodically FMU values (in case you are interested in any)
function saveFunc(nfmu, u, t, integrator)
    values = fmi2GetReal(nfmu.fmu, nfmu.recordValues)
    return (values...,)
end

"""
Constructs a ME-NeuralFMU where the FMU is at an unknown location inside of the NN.

Arguents:
    - `saveat` time points to save the NeuralFMU output, if empty, solver step size is used (may be non-equidistant)
    - `addTopLayer` adds an input layer to add a time state (default) to set the FMU-time
    - `addBottomLayer` adds an input layer to remove the time state (default)

Keyword arguments:
    - `fixstep` forces fixed step integration
    - `recordFMUValues` additionally records internal FMU variables (currently not supported because of issues)
"""
function ME_NeuralFMU(fmu, model, tspan, alg=nothing; saveat=[], addTopLayer = true, addBottomLayer = true, recordFMUValues = [], fixstep=0.0)
    nfmu = ME_NeuralFMU()
    nfmu.fmu = fmu

    ext_model = nothing

    if addTopLayer && addBottomLayer
        ext_model = Chain(inputs -> NeuralFMUInputLayer(nfmu.fmu, inputs),
                          model.layers...,
                          inputs -> NeuralFMUOutputLayerME(inputs))
    elseif addTopLayer
        ext_model = Chain(inputs -> NeuralFMUInputLayer(nfmu.fmu, inputs),
                          model.layers...)
    elseif addBottomLayer
        ext_model = Chain(model.layers...,
                          inputs -> NeuralFMUOutputLayerME(inputs))
    else
        ext_model = Chain(model.layers...)
    end

    nfmu.recordValues = recordFMUValues
    if length(nfmu.recordValues) > 0
        @assert false "ME_NeuralFMU(...): keyword `recordFMUValues` is currently not supported (under development). Please remove keyword."

        tmp = zeros(Float64, length(nfmu.recordValues))
        type = typeof((tmp...,))

        nfmu.saved_values = SavedValues(Float64, type)
        cb = SavingCallback((u,t,integrator)->saveFunc(nfmu, u, t, integrator), nfmu.saved_values; saveat=saveat)

        if fixstep > 0.0
            nfmu.neuralODE = NeuralODE(ext_model, tspan, alg, callback=cb, dt=fixstep, adaptive=false)
        else
            nfmu.neuralODE = NeuralODE(ext_model, tspan, alg, callback=cb)
        end
    else
        if fixstep > 0.0
            nfmu.neuralODE = NeuralODE(ext_model, tspan, alg; saveat=saveat, dt=fixstep, adaptive=false)
        else
            nfmu.neuralODE = NeuralODE(ext_model, tspan, alg; saveat=saveat)
        end
    end

    nfmu.tspan = tspan
    nfmu.saveat = saveat

    nfmu
end

"""
Constructs a CS-NeuralFMU where the FMU is at an unknown location inside of the NN.

Arguents:
    - `saveat` time points to save the NeuralFMU output, if empty, solver step size is used (may be non-equidistant)
    - `addTopLayer` adds an input layer to add a time state (default) to set the FMU-time
    - `addBottomLayer` adds an input layer to remove the time state (default)
"""
function CS_NeuralFMU(fmu, model, tspan; saveat=[], addTopLayer = true, addBottomLayer = true, recordValues = [])
    nfmu = CS_NeuralFMU()
    nfmu.fmu = fmu

    if addTopLayer && addBottomLayer
        nfmu.model = Chain(inputs -> NeuralFMUInputLayer(nfmu.fmu, inputs),
                          model.layers...,
                          inputs -> NeuralFMUOutputLayerCS(nfmu.fmu, inputs))
    elseif addTopLayer
        nfmu.model = Chain(inputs -> NeuralFMUInputLayer(nfmu.fmu, inputs),
                          model.layers...)
    elseif addBottomLayer
        nfmu.model = Chain(model.layers...,
                          inputs -> NeuralFMUOutputLayerCS(nfmu.fmu, inputs))
    else
        nfmu.model = model.layers
    end

    nfmu.tspan = tspan
    nfmu.saveat = saveat

    nfmu
end

function (nfmu::ME_NeuralFMU)(x_start, t_start::Real = nfmu.tspan[1]; reset::Bool=true)

    if reset
        fmiReset(nfmu.fmu)
        fmiSetupExperiment(nfmu.fmu, t_start)
        fmiEnterInitializationMode(nfmu.fmu)
        fmiExitInitializationMode(nfmu.fmu)
    end

    x0 = [t_start, x_start...]
    nfmu.solution = nfmu.neuralODE(x0)
end

function (nfmu::CS_NeuralFMU)(t_step, t_start = nfmu.tspan[1], t_stop = nfmu.tspan[end]; inputs=[], reset::Bool=true)

    if reset
        while fmiReset(nfmu.fmu) != 0
        end
        while fmiSetupExperiment(nfmu.fmu, t_start) != 0
        end
        while fmiEnterInitializationMode(nfmu.fmu) != 0
        end
        while fmiExitInitializationMode(nfmu.fmu) != 0
        end
    end

    ts = t_start:t_step:t_stop
    modelInput = collect.(eachrow(hcat(ts, inputs)))
    nfmu.valueStack = nfmu.model.(modelInput)

    nfmu.valueStack
end

# adapting the Flux functions
function Flux.params(neuralfmu::ME_NeuralFMU)
    Flux.params(neuralfmu.neuralODE)
end

function Flux.params(neuralfmu::CS_NeuralFMU)
    Flux.params(neuralfmu.model)
end

# FMI version independent dosteps
# function fmiDoStepME(fmu::FMU2, x, t = -1.0; setValueReferences = [], setValues = [], getValueReferences = [])
#     fmi2DoStepME(fmu, x, t;
#                 setValueReferences = setValueReferences,
#                 setValues = setValues,
#                 getValueReferences = getValueReferences)
# end
function fmiDoStepME(fmu::FMU2, x, t = -1.0, setValueReferences = [], setValues = [], getValueReferences = [])
    fmi2DoStepME(fmu, x, t,
                setValueReferences,
                setValues,
                getValueReferences)
end

function fmiDoStepCS(fmu::FMU2, dt; setValueReferences = [], setValues = [], getValueReferences = [])
    fmi2DoStepCS(fmu, dt;
                setValueReferences = setValueReferences,
                setValues = setValues,
                getValueReferences = getValueReferences)
end

function fmiInputDoStepCSOutput(fmu::FMU2, dt, u)
    fmi2InputDoStepCSOutput(fmu, dt, u)
end

# define neutral gradients (=feed-trough) for ccall-functions
function neutralGradient(c̄)
    c̄
end

@adjoint fmiSetupExperiment(fmu, startTime, stopTime) = fmiSetupExperiment(fmu, startTime, stopTime), c̄ -> neutralGradient(c̄)
@adjoint fmiEnterInitializationMode(fmu) = fmiEnterInitializationMode(fmu), c̄ -> neutralGradient(c̄)
@adjoint fmiExitInitializationMode(fmu) = fmiExitInitializationMode(fmu), c̄ -> neutralGradient(c̄)
@adjoint fmiReset(fmu) = fmiReset(fmu), c̄ -> neutralGradient(c̄)

ext: jl


text:
using Pkg
current_env = Base.load_path()[1]
Pkg.activate(@__DIR__)

using RadialBasisFunctionModels
using Test
using InteractiveUtils: subtypes

#%% initialize empty machine
# should default to 64 bit precision
@testset "Initialization" begin
	mach = RBFMachine(;poly_deg = 1)
	@test mach.features isa Vector{Vector{Float64}}
	@test mach.labels isa Vector{Vector{Float64}}

	@test mach.φ isa Gaussian

	@test_throws AssertionError RadialBasisFunctionModels.fit!(mach)

	add_data!( mach, rand(1), rand(1) )
	@test_throws AssertionError RadialBasisFunctionModels.fit!(mach)

	add_data!( mach, rand(1), rand(1) )
	@test isnothing(RadialBasisFunctionModels.fit!(mach))
	# do we interpolate?
	@test mach( mach.features[1] ) ≈ mach.labels[1]

	@test mach.valid

	add_data!(mach, rand(1), rand(1))

	@test !mach.valid
end

#%%
@testset "Precision" begin
	for T in subtypes(AbstractFloat)
		features = [ rand(T, 2) for i = 1 : 4 ]
		labels = [ rand(T, 3) for i = 1 : 4 ]

		mach = RBFMachine(; features, labels )
		F = typeof(mach)
		@test F.parameters[1] == typeof(features) == Vector{Vector{T}}
		@test F.parameters[2] == typeof(labels) == Vector{Vector{T}}

		RadialBasisFunctionModels.fit!(mach)

		@test mach.valid
		T == Float16 || @test mach( mach.features[1] ) ≈ mach.labels[1]
		@test mach( features[1] ) isa Vector{T}

		for S in subtypes(AbstractFloat)
			add_data!(mach, rand(S,2), rand(S,3) )
		end
		@test !mach.valid
		@test mach.features isa Vector{Vector{T}}
		@test mach.labels isa Vector{Vector{T}}

	end
end

#%% Test if precision is kept if kernel_args have different precision
@testset "PrecisionII" begin
	for T in subtypes(AbstractFloat)
		for S in subtypes(AbstractFloat)
			features = [ rand(T, 2) for i = 1 : 4 ]
			labels = [ rand(T, 3) for i = 1 : 4 ]

			mach = RBFMachine(;
				kernel_name = :gaussian,
				kernel_args = ones( S, 1 ),
				features, labels
			)
			F = typeof(mach)
			@test F.parameters[1] == typeof(features) == Vector{Vector{T}}
			@test F.parameters[2] == typeof(labels) == Vector{Vector{T}}

			RadialBasisFunctionModels.fit!(mach)

			@test mach.valid
			# @test mach( mach.features[1] ) ≈ mach.labels[1] fails for Float16
			@test mach( features[1] ) isa Vector{T}
		end
	end
end

#%%
@testset "Kernels" begin
	features = [ rand(2) for i = 1 : 4 ]
	labels = [ rand(3) for i = 1 : 4 ]
	for kn in keys(RadialBasisFunctionModels.SymbolToRadialConstructor)
		if kn != :thin_plate_spline # default thin plate spline has cpd order = 3 => needs quadratic polynomial
			mach = RBFMachine(;
				kernel_name = kn,
				poly_deg = 1,
				features, labels
			)
			@test isnothing(RadialBasisFunctionModels.fit!(mach))
			@test mach.valid
			@test all( mach(features[i]) ≈ labels[i] for i = 1 : 4 )
		else
			@test_throws AssertionError RBFMachine(;
				kernel_name = kn,
				poly_deg = 1,
				features, labels
			)

			mach = RBFMachine(;
				kernel_name = kn,
				kernel_args = [1,],
				poly_deg = 1,
				features, labels
			)
			@test isnothing(RadialBasisFunctionModels.fit!(mach))
			@test mach.valid
			@test all( mach(features[i]) ≈ labels[i] for i = 1 : 4 )
		end
	end

	# Multiquadric has cpd order 1 => needs constant polynomial
	@test_throws AssertionError RBFMachine( kernel_name = :multiquadric, poly_deg = -1 )

	# Cubic has order 2
	@test_throws AssertionError RBFMachine( kernel_name = :cubic, poly_deg = -1 )
	@test_throws AssertionError RBFMachine( kernel_name = :cubic, poly_deg = 0 )

	for kernel_name in [:gaussian, :inv_multiquadric], poly_deg in [-1,0,1]
		@test RBFMachine(;kernel_name, poly_deg) isa RBFMachineWithKernel
	end
end

#%%
Pkg.activate(current_env)


text:
# MIT License

# Copyright (c) Microsoft Corporation.

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE

abstract type AbstractRay{T<:Real,N} end
export AbstractRay

"""
    Ray{T,N} <: AbstractRay{T,N}

Purely geometric ray, defined as `origin + alpha * direction`.

```julia
Ray(origin::SVector{N,T}, direction::SVector{N,T})
```

Has the following accessor methods:
```julia
direction(ray::Ray{T,N}) -> SVector{N,T}
origin(ray::Ray{T,N}) -> SVector{N,T}
```
"""
struct Ray{T,N} <: AbstractRay{T,N}
    origin::SVector{N,T}
    direction::SVector{N,T}

    function Ray(origin::SVector{N,T}, direction::SVector{N,T}) where {T<:Real,N}
        return new{T,N}(origin, normalize(direction))
    end

    function Ray(origin::AbstractArray{T,1}, direction::AbstractArray{T,1}) where {T<:Real}
        @assert length(origin) == length(direction)
        N = length(origin)
        return new{T,N}(SVector{N,T}(origin), normalize(SVector{N,T}(direction)))
    end
end
export Ray

direction(ray::Ray{T,N}) where {T<:Real,N} = ray.direction
origin(ray::Ray{T,N}) where {T<:Real,N} = ray.origin
export origin, direction

function Base.print(io::IO, a::Ray{T,N}) where {T,N}
    println(io, "$(rpad("Origin:", 20)) $(origin(a))")
    println(io, "$(rpad("Direction:", 20)) $(direction(a))")
end

"""
    point(ray::AbstractRay{T,N}, alpha::T) -> SVector{T, N}

Returns a point on the ray at origin + alpha * direction. Alpha must be >= 0.
"""
function point(ray::AbstractRay{T,N}, alpha::T) where {N,T<:Real}
    @assert alpha >= zero(T) "Alpha must be nonnegative. alpha value: $alpha"
    return origin(ray) + alpha * direction(ray)
end

"""
    closestpointonray(r::Ray{T,N}, point::SVector{N,T}) -> SVector{T,N

Returns the point on the ray closest to point.
"""
function closestpointonray(r::AbstractRay{T,N}, point::SVector{N,T}) where {T,N}
    #find t that gives the smallest distance between the ray and the point
    #ray = p0 + tr̂, distance = ||line - point||
    o = origin(r)
    r̂ = direction(r)
    t = dot(point .- o, r̂)

    if t >= 0.0
        return o .+ t .* r̂
    else
        return o
    end
end

"""
    distance(r::Ray{T,N}, point::SVector{N,T}) -> Union{T,Nothing}

Returns distance to the position on the ray closest to point. If t < 0 returns nothing.
"""
function distance(r::AbstractRay{T,N}, point::SVector{N,T}) where {T,N}
    #find t that gives the smallest distance between the ray and the point
    #ray = p0 + tr̂, distance = ||line - point||
    # this has redundant code with closestpointonray but it is faster to do this way than to call closestpointonray and then comput norm(result .- point)

    o = origin(r)
    r̂ = direction(r)
    t = dot(point .- o, r̂)

    if t >= 0.0
        temp = o .+ t .* r̂
        return norm(temp .- point)
    else
        return nothing
    end
end

"""
    α(ray::AbstractRay{T,N}, point::SVector{N,T}) -> T

Computes the alpha corresponding to the closest position on the ray to point
"""
α(ray::AbstractRay{T,N}, point::SVector{N,T}) where {T<:Real,N} = dot(direction(ray), (point .- origin(ray)))


text:
export DiffEqObjective, build_loss_objective

struct DiffEqObjective{F,F2} <: Function
  cost_function::F
  cost_function2::F2
end

function diffeq_sen_full(f, u0, tspan, p, t, alg; kwargs...)
  prob = ODEForwardSensitivityProblem(f,u0,tspan,p)
  sol = solve(prob,alg;kwargs...)(t)
  nvar = length(u0)
  sol[1:nvar,:], [sol[i*nvar+1:i*nvar+nvar,:] for i in 1:length(p)]
end

function diffeq_sen_l2!(res, df, u0, tspan, p, t, data, alg;kwargs...)
  prob = ODEProblem(df,u0,tspan,p)
  sol = solve(prob, alg, saveat=t; kwargs...)
  function dg(out,u,p,t,i)
    @. out = 2 * (data[:,i] - u)
  end
  res .= adjoint_sensitivities(sol,alg,dg,t,kwargs...)[2][1,:]
end

(f::DiffEqObjective)(x) = f.cost_function(x)
(f::DiffEqObjective)(x,y) = f.cost_function2(x,y)

function build_loss_objective(prob::DiffEqBase.DEProblem,alg,loss,regularization=nothing,args...;
                              prior=nothing,mpg_autodiff = false,
                              verbose_opt = false,verbose_steps = 100,
                              prob_generator = STANDARD_PROB_GENERATOR,
                              autodiff_prototype = mpg_autodiff ? zero(prob.p) : nothing,
                              autodiff_chunk = mpg_autodiff ? ForwardDiff.Chunk(autodiff_prototype) : nothing,flsa_gradient = false,adjsa_gradient= false,
                              kwargs...)
  if verbose_opt
    count = 0 # keep track of # function evaluations
  end
  cost_function = function (p)
    tmp_prob = prob_generator(prob,p)
    if typeof(loss) <: Union{L2Loss,LogLikeLoss}
      sol = solve(tmp_prob,alg,args...;saveat=loss.t,save_everystep=false,dense=false,kwargs...)
    else
      sol = solve(tmp_prob,alg,args...;kwargs...)
    end

    loss_val = loss(sol)

    if prior != nothing
      loss_val += prior_loss(prior,p)
    end
    if regularization != nothing
      loss_val += regularization(p)
    end

    if verbose_opt
      count::Int += 1
      if mod(count,verbose_steps) == 0
        println("Iteration: $count")
        println("Current Cost: $loss_val")
        println("Parameters: $p")
      end
    end
    loss_val
  end

  if mpg_autodiff
    gcfg = ForwardDiff.GradientConfig(cost_function, autodiff_prototype, autodiff_chunk)
    g! = (x, out) -> ForwardDiff.gradient!(out, cost_function, x, gcfg)
  elseif flsa_gradient
    if typeof(loss) <: L2Loss
      function g!(x,out)
        sol_,sens = diffeq_sen_full(prob.f,prob.u0,prob.tspan,x,loss.t,alg)
        l2lossgradient!(out,sol_,loss.data,sens,length(prob.p))
      end
    else
      throw("LSA gradient only for L2Loss")
    end
  elseif adjsa_gradient
    g! = (x,out) -> diffeq_sen_l2!(out,prob.f,prob.u0,prob.tspan,x,loss.t,loss.data,alg)
  else
    g! = (x, out) -> Calculus.finite_difference!(cost_function,x,out,:central)
  end

  cost_function2 = function (p,grad)
    if length(grad)>0
      g!(p,grad)
    end
    cost_function(p)
  end
  DiffEqObjective(cost_function,cost_function2)
end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

run barrier tests
=#

using Test
using Printf
import Hypatia.Cones
include(joinpath(@__DIR__, "cone.jl"))

function cone_types(T::Type{<:Real})
    cones_T = [
        Cones.Nonnegative{T},
        Cones.PosSemidefTri{T, T},
        Cones.PosSemidefTri{T, Complex{T}},
        Cones.DoublyNonnegativeTri{T},
        Cones.PosSemidefTriSparse{Cones.PSDSparseDense, T, T},
        Cones.PosSemidefTriSparse{Cones.PSDSparseDense, T, Complex{T}},
        Cones.LinMatrixIneq{T},
        Cones.EpiNormInf{T, T},
        Cones.EpiNormInf{T, Complex{T}},
        Cones.EpiNormEucl{T},
        Cones.EpiPerSquare{T},
        Cones.EpiNormSpectral{T, T},
        Cones.EpiNormSpectral{T, Complex{T}},
        Cones.MatrixEpiPerSquare{T, T},
        Cones.MatrixEpiPerSquare{T, Complex{T}},
        Cones.GeneralizedPower{T},
        Cones.HypoPowerMean{T},
        Cones.HypoGeoMean{T},
        Cones.HypoRootdetTri{T, T},
        Cones.HypoRootdetTri{T, Complex{T}},
        Cones.HypoPerLog{T},
        Cones.HypoPerLogdetTri{T, T},
        Cones.HypoPerLogdetTri{T, Complex{T}},
        Cones.EpiPerSepSpectral{Cones.VectorCSqr{T}, T},
        Cones.EpiPerSepSpectral{Cones.MatrixCSqr{T, T}, T},
        Cones.EpiPerSepSpectral{Cones.MatrixCSqr{T, Complex{T}}, T},
        Cones.EpiRelEntropy{T},
        Cones.EpiTrRelEntropyTri{T},
        Cones.WSOSInterpNonnegative{T, T},
        Cones.WSOSInterpNonnegative{T, Complex{T}},
        Cones.WSOSInterpPosSemidefTri{T},
        Cones.WSOSInterpEpiNormEucl{T},
        Cones.WSOSInterpEpiNormOne{T},
        ]

    if T <: LinearAlgebra.BlasReal
        append!(cones_T, [
            Cones.PosSemidefTriSparse{Cones.PSDSparseCholmod, T, T},
            Cones.PosSemidefTriSparse{Cones.PSDSparseCholmod, T, Complex{T}},
            ])
    end

    return cones_T
end

sep_spectral_funs = [
    Cones.InvSSF(),
    Cones.NegLogSSF(),
    Cones.NegEntropySSF(),
    Cones.Power12SSF(1.5),
    ]

@testset "cone tests" begin

println("starting oracle tests")
@testset "oracle tests" begin
real_types = [
    Float64,
    # Float32,
    # BigFloat,
    ]
@testset "$cone" for T in real_types, cone in cone_types(T)
    println("$cone")
    test_time = @elapsed test_oracles(cone)
    @printf("%8.2e seconds\n", test_time)
end
end

# println("\nstarting barrier tests")
# @testset "barrier tests" begin
# real_types = [
#     Float64,
#     # Float32,
#     # BigFloat,
#     ]
# @testset "$cone" for T in real_types, cone in cone_types(T)
#     println("$cone")
#     test_time = @elapsed test_barrier(cone)
#     @printf("%8.2e seconds\n", test_time)
# end
# end
#
# println("\nstarting time/allocation measurements")
# @testset "allocation tests" begin
# real_types = [
#     Float64,
#     # Float32,
#     # BigFloat,
#     ]
# @testset "$cone" for T in real_types, cone in cone_types(T)
#     println("\n$cone")
#     test_time = @elapsed show_time_alloc(cone)
#     @printf("%8.2e seconds\n", test_time)
# end
# println()
# end

end
;


text:
using OpenVX

include("imageio.jl")

function make_filter_graph(ctx, input::VX.vx_image, output::VX.vx_image, width, height)

    graph = VX.CreateGraph(ctx)

    yuv_image = VX.CreateVirtualImage(graph, width, height, VX.DF_IMAGE_IYUV)
    luma_image = VX.CreateVirtualImage(graph, width, height, VX.DF_IMAGE_U8)

    hyst = VX.CreateThreshold(ctx, VX.THRESHOLD_TYPE_RANGE, VX.TYPE_UINT8)

    lower = Ref{Int32}(80)
    upper = Ref{Int32}(100)

    VX.SetThresholdAttribute(hyst, VX.THRESHOLD_ATTRIBUTE_THRESHOLD_LOWER,
                             lower, sizeof(Int32))
    VX.SetThresholdAttribute(hyst, VX.THRESHOLD_ATTRIBUTE_THRESHOLD_UPPER,
                             upper, sizeof(Int32))

    gradient_size = 3

    nodes = [VX.ColorConvertNode(graph, input, yuv_image),
             VX.ChannelExtractNode(graph, yuv_image, VX.CHANNEL_Y, luma_image),
             VX.CannyEdgeDetectorNode(graph, luma_image, hyst, gradient_size,
                                      VX.NORM_L1, output)]

    VX.ReleaseImage(Ref(yuv_image))
    VX.ReleaseImage(Ref(luma_image))

    for i in 1:length(nodes)
        VX.ReleaseNode(Ref(nodes[i]))
    end

    return graph
end

ctx = VX.CreateContext()

filename = "assets/airplane.bmp"
width = 512
height = 512

image = image_from_file(filename, width, height)
output = VX.CreateImage(ctx, width, height, VX.DF_IMAGE_U8)

graph = make_filter_graph(ctx, image, output, width, height)
VX.VerifyGraph(graph)
VX.ProcessGraph(graph)
file_from_gray_image(output, "airplane_canny.bmp", width, height)


text:
#!/usr/bin/env julia

include("../src/AcceleratedDubins.jl")
using .AcceleratedDubins

using PyPlot
using PyCall

# Save figures into a directory "images"
save_figures = false
if save_figures
    mkpath("images")
end

for run in 1:1
    println("COMROB: Dubins path with multiple radii demo")
    println("Library examples")

    ##############################
    ### Vehicle settings
    minr = 1.
    maxr = 5.

    vr = 5.
    vs = 15.

    speed_params = [vr, vs, 2.6, 4]

    ##############################
    ### Position settings

    ### tested settings with visible differences
    # alpha = 0.8*pi
    # beta = 0.1*pi
    # start = [6, -7, alpha]
    # stop = [19, -5, beta]

    ### random settings
    alpha = 2*pi*rand()
    beta = 2*pi*rand()
    start = [10*rand(), -10+5*rand(), alpha]
    stop = [10+20*rand(), -10+5*rand(), beta]

    while true
        fig, ax = plt.subplots(1, 2, num="Fastest vs. shortest")
        # fig.set_size_inches(6.4*2, 4.8)
        # plt.clf()
        ax[1].axis("equal")
        ax[1].set_xlabel("x")
        ax[1].set_ylabel("y")

        ax[2].axis("equal")
        ax[2].set_ylabel("Speed [m/s]")
        ax[2].set_xlabel("Time [s]")
        # ax[2].set_ylim(vr, vs)
        # simple path computation
        # compute fastest path
        # tries all radii combinations
        _, path_fastest = fastest_path(start, stop, Array(minr:2.0:maxr), speed_params)
        # tries optimized path
        _, path_optimized = optimized_path(start, stop, minr, maxr, speed_params)
        # compute shortest path
        _, path_shortest = retrieve_path(start, stop, minr, minr, speed_params)
        # sample the path for plotting
        confx, confy = sample_path(path_fastest)
        ax[1].plot(confx, confy, color="green")
        confx, confy = sample_path(path_optimized)
        ax[1].plot(confx, confy, color="blue")
        confx, confy = sample_path(path_shortest)
        ax[1].plot(confx, confy, color="orange")

        times, speeds = speed_profile(path_fastest, speed_params)
        ax[2].plot(times, speeds, color="green")
        times, speeds = speed_profile(path_optimized, speed_params)
        ax[2].plot(times, speeds, color="blue")
        times, speeds = speed_profile(path_shortest, speed_params)
        ax[2].plot(times, speeds, color="orange")

        time_fastest = path_time(path_fastest, speed_params)
        time_fastest_no_opt = (sum(path_fastest.lengths) / speed_by_radius(minr, speed_params[1], speed_params[2]))

        time_optimized = path_time(path_optimized, speed_params)
        time_optimized_no_opt = (sum(path_optimized.lengths) / speed_by_radius(minr, speed_params[1], speed_params[2]))

        time_shortest = path_time(path_shortest, speed_params)
        time_shortest_no_opt = (sum(path_shortest.lengths) / speed_by_radius(minr, speed_params[1], speed_params[2]))

        println("[GREEN] Fastest path:    vp = ", round(time_fastest, digits=3), "    no_opt = ", round(time_fastest_no_opt, digits=3))
        println("[BLUE] Optimized path:   vp = ", round(time_optimized, digits=3), "    no_opt = ", round(time_optimized_no_opt, digits=3))
        println("[ORANGE] Shortest path:  vp = ", round(time_shortest, digits=3), "    no_opt = ", round(time_shortest_no_opt, digits=3))

        plt.pause(0.1)
        println("Press enter for next random demo, q and enter to quit")
        str = readline()
        close()
        if occursin("q", str) || occursin("Q", str)
            break
        end

        alpha = 2*pi*rand()
        beta = 2*pi*rand()
        start = [10*rand(), -10+5*rand(), alpha]
        stop = [10+20*rand(), -10+5*rand(), beta]
    end

end


text:
function mandel( c )
    z = c
    itrMax = 80
    for n in 1:itrMax
        if abs(z) > 2
            return n-1
        end
        z = z^2 + c
    end
    return itrMax
end


text:
include("demo_pulsed_deck.jl")
include("demo_steady_deck.jl")

include("arc_deck.jl")

include("act_1_deck.jl")
include("act_2_deck.jl")

include("proteus_deck.jl")
include("charybdis_deck.jl")


text:
export extract_amp, extract_vis, extract_lcamp, extract_cphase

function getvisfield(obs)
    obsamps = obs.data::PyObject
    u = deepcopy((get(obsamps, Vector{Float64}, "u")))
    v = deepcopy((get(obsamps, Vector{Float64}, "v")))
    err = deepcopy((get(obsamps, Vector{Float64}, "sigma")))
    vis = deepcopy((get(obsamps, Vector{Complex{Float64}}, "vis")))
    t1 = Symbol.(deepcopy((get(obsamps, Vector{String}, "t1"))))
    t2 = Symbol.(deepcopy((get(obsamps, Vector{String}, "t2"))))
    baseline = tuple.(t1, t2)
    time = deepcopy((get(obsamps, Vector{Float64}, "time")))
    freq = zeros(length(time))
    bw = zeros(length(time))

    return  StructArray{Comrade.EHTVisibilityDatum{Float64}}(
        visr = real.(vis),
        visi = imag.(vis),
        u = u,
        v = v,
        error = err,
        time = time,
        frequency = freq,
        bandwidth = bw,
        baseline = baseline
    )
end

function getampfield(obs)
    obsamps = obs.amp::PyObject
    uamp = deepcopy(get(obsamps, Vector{Float64}, "u"))
    vamp = deepcopy(get(obsamps, Vector{Float64}, "v"))
    erramp = deepcopy(get(obsamps, Vector{Float64}, "sigma"))
    amps = deepcopy(get(obsamps, Vector{Float64}, "amp"))
    t1 = Symbol.(deepcopy(get(obsamps, Vector{String}, "t1")))
    t2 = Symbol.(deepcopy(get(obsamps, Vector{String}, "t2")))
    baseline = tuple.(t1, t2)
    time = deepcopy(get(obsamps, Vector{Float64}, "time"))
    freq = zeros(length(time))
    bw = zeros(length(time))

    return  StructArray{Comrade.EHTVisibilityAmplitudeDatum{Float64}}(
        amp = amps,
        u = uamp,
        v = vamp,
        error = erramp,
        time = time,
        frequency = freq,
        bandwidth = bw,
        baseline = baseline
    )
end

function getcpfield(obs)
    obscp = obs.cphase::PyObject
    u1 = deepcopy((get(obscp, Vector{Float64}, "u1")))
    v1 = deepcopy((get(obscp, Vector{Float64}, "v1")))
    u2 = deepcopy((get(obscp, Vector{Float64}, "u2")))
    v2 = deepcopy((get(obscp, Vector{Float64}, "v2")))
    u3 = deepcopy((get(obscp, Vector{Float64}, "u3")))
    v3 = deepcopy((get(obscp, Vector{Float64}, "v3")))
    cp = deg2rad.(deepcopy((get(obscp, Vector{Float64}, "cphase"))))
    errcp = deg2rad.(deepcopy((get(obscp, Vector{Float64}, "sigmacp"))))

    t1 = Symbol.(deepcopy((get(obscp, Vector{String}, "t1"))))
    t2 = Symbol.(deepcopy((get(obscp, Vector{String}, "t2"))))
    t3 = Symbol.(deepcopy((get(obscp, Vector{String}, "t3"))))
    baseline = tuple.(t1, t2, t3)
    time = deepcopy(get(obscp, Vector{Float64}, "time"))
    freq = zeros(length(time))
    bw = zeros(length(time))

    return StructArray{Comrade.EHTClosurePhaseDatum{Float64}}(
        phase = cp,
        u1 = u1,
        v1 = v1,
        u2 = u2,
        v2 = v2,
        u3 = u3,
        v3 = v3,
        error = errcp,
        time = time,
        frequency = freq,
        bandwidth = bw,
        triangle = baseline
    )

end

function getlcampfield(obs)
    obslcamp = obs.logcamp::PyObject
    u1 = deepcopy((get(obslcamp, Vector{Float64}, "u1")))
    v1 = deepcopy((get(obslcamp, Vector{Float64}, "v1")))
    u2 = deepcopy((get(obslcamp, Vector{Float64}, "u2")))
    v2 = deepcopy((get(obslcamp, Vector{Float64}, "v2")))
    u3 = deepcopy((get(obslcamp, Vector{Float64}, "u3")))
    v3 = deepcopy((get(obslcamp, Vector{Float64}, "v3")))
    u4 = deepcopy((get(obslcamp, Vector{Float64}, "u4")))
    v4 = deepcopy((get(obslcamp, Vector{Float64}, "v4")))
    camp = deepcopy((get(obslcamp, Vector{Float64}, "camp")))
    errcamp = deepcopy((get(obslcamp, Vector{Float64}, "sigmaca")))

    t1 = Symbol.(deepcopy((get(obslcamp, Vector{String}, "t1"))))
    t2 = Symbol.(deepcopy((get(obslcamp, Vector{String}, "t2"))))
    t3 = Symbol.(deepcopy((get(obslcamp, Vector{String}, "t3"))))
    t4 = Symbol.(deepcopy((get(obslcamp, Vector{String}, "t4"))))
    baseline = tuple.(t1, t2, t3, t4)
    time = deepcopy(get(obslcamp, Vector{Float64}, "time"))
    freq = zeros(length(time))
    bw = zeros(length(time))

    return StructArray{Comrade.EHTLogClosureAmplitudeDatum{Float64}}(
        amp = camp,
        u1 = u1,
        v1 = v1,
        u2 = u2,
        v2 = v2,
        u3 = u3,
        v3 = v3,
        u4 = u4,
        v4 = v4,
        error = errcamp,
        time = time,
        frequency = freq,
        bandwidth = bw,
        quadrangle = baseline
    )

end

function getradec(obs)::Tuple{Float64, Float64}
    return (float(obs.ra), float(obs.dec))
end

function getmjd(obs)::Int
    return Int(obs.mjd)
end

function getsource(obs)::Symbol
    return Symbol(obs.source)
end

"""
    extract_amp(obs)
Extracts the visibility amplitudes from an ehtim observation object

Returns an EHTObservation with visibility amplitude data
"""

function extract_amp(obs; kwargs...)
    obs.add_amp(;kwargs...)
    data = getampfield(obs)
    ra, dec = getradec(obs)
    mjd = getmjd(obs)
    source = getsource(obs)
    bw = obs.bw
    rf = obs.rf

    return Comrade.EHTObservation(data = data, mjd = mjd,
                   ra = ra, dec= dec,
                   bandwidth=bw, frequency=rf,
                   source = source,
    )
end

"""
    extract_vis(obs)
Extracts the complex visibilities from an ehtim observation object

Returns an EHTObservation with complex visibility data
"""
function extract_vis(obs)
    data = getvisfield(obs)
    ra, dec = getradec(obs)
    mjd = getmjd(obs)
    source = getsource(obs)
    bw = obs.bw
    rf = obs.rf
    return Comrade.EHTObservation(
                   data = data, mjd = mjd,
                   ra = ra, dec= dec,
                   bandwidth=bw, frequency=rf,
                   source = source,
    )
end

"""
    extract_cphase(obs)
Extracts the closure phases from an ehtim observation object

Returns an EHTObservation with closure phases datums
"""
function extract_cphase(obs; kwargs...)
    obs.add_cphase(;kwargs...)
    data = getcpfield(obs)
    ra, dec = getradec(obs)
    mjd = getmjd(obs)
    source = getsource(obs)
    bw = obs.bw
    rf = obs.rf

    return Comrade.EHTObservation(data = data, mjd = mjd,
                   ra = ra, dec= dec,
                   bandwidth=bw, frequency=rf,
                   source = source,
    )

end

"""
    extract_lcamp(obs)
Extracts the log-closure amp. from an ehtim observation object

Returns an EHTObservation with closure amp. datums
"""
function extract_lcamp(obs; kwargs...)
    obs.add_logcamp(;kwargs...)
    data = getlcampfield(obs)
    ra, dec = getradec(obs)
    mjd = getmjd(obs)
    source = getsource(obs)
    bw = obs.bw
    rf = obs.rf

    return Comrade.EHTObservation(data = data, mjd = mjd,
                   ra = ra, dec= dec,
                   bandwidth=bw, frequency=rf,
                   source = source,
    )

end


text:
tol = 1e-7
################################################
# Uncertain values represented by KDE estimates
################################################

# Create a random sample
d = Normal()
some_sample = rand(d, 2000)

# Create an uncertain value from the random sample by kernel density
# estimation.
uv = UncertainValue(some_sample)

@test resample(uv) isa Float64
@test resample(uv, 10) isa Vector{Float64}

@test resample(uv, NoConstraint()) isa Float64
@test resample(uv, NoConstraint(), 10) isa Vector{Float64}

r1a = resample(uv, TruncateLowerQuantile(0.2))
r1b = resample(uv, TruncateLowerQuantile(0.2), 10)
@test r1a isa Float64
@test r1b isa Vector{Float64}
#@test quantile(uv, 0.2) <= r1a + tol
##@test all(quantile(uv, 0.2) .<= r1b .+ tol)

r1a = resample(uv, TruncateUpperQuantile(0.8))
r1b = resample(uv, TruncateUpperQuantile(0.8), 10)
@test r1a isa Float64
@test r1b isa Vector{Float64}
#@test r1a <= quantile(uv, 0.8) + tol
#@test all(r1b .<= quantile(uv, 0.8) + tol)

r1a = resample(uv, TruncateQuantiles(0.2, 0.8))
r1b = resample(uv, TruncateQuantiles(0.2, 0.8), 10)
@test r1a isa Float64
@test r1b isa Vector{Float64}
#@test quantile(uv, 0.2) - tol <= r1a <= quantile(uv, .8) + tol
#@test all(quantile(uv, 0.2) - tol .<= r1b .<= quantile(uv, .8) + tol)

r1a = resample(uv, TruncateMinimum(-0.5))
r1b = resample(uv, TruncateMinimum(-0.5), 10)
@test r1a isa Float64
@test r1b isa Vector{Float64}
#@test -0.5 <= r1a + tol
#@test all(-0.5 .<= r1b .+ tol)

r1a = resample(uv, TruncateMaximum(0.5))
r1b = resample(uv, TruncateMaximum(0.5), 10)
@test r1a isa Float64
@test r1b isa Vector{Float64}
#@test r1a <= 0.5 + tol
#@test all(r1b .<= 0.5 + tol)

r1a = resample(uv, TruncateRange(-0.5, 0.5))
r1b = resample(uv, TruncateRange(-0.5, 0.5), 10)
@test r1a isa Float64
@test r1b isa Vector{Float64}
#@test -0.5 - tol <= r1a <= 0.5 + tol
#@test all(-0.5 - tol .<= r1b .<= 0.5 + tol)

# TruncateStd will not work, but with the default fallback we should still be able
# to resample.
r1a = resample(uv, TruncateStd(1))
r1b = resample(uv, TruncateStd(1), 10)
@test r1a isa Float64
@test r1b isa Vector{Float64}


text:
__precompile__()
module StructIO

using Base: @pure
using Base.Meta
export @io, unpack, pack, fix_endian, packed_sizeof

"""
    needs_bswap(endianness::Symbol)

Returns `true` if the given endianness does not match the current host system.
"""
@pure function needs_bswap(endianness::Symbol)
    if ENDIAN_BOM == 0x01020304
        return endianness == :LittleEndian
    else
        return endianness == :BigEndian
    end
end

"""
    bswap!(ptr::Ptr{UInt8}, sz)

Byte-swap a chunk of data in-place
"""
function bswap!(ptr::Ptr{UInt8}, sz)
    # Count from outside edge to middle
    for i = 0:div(sz-1,2)
        # Swap two mirrored bytes
        ptr_hi = ptr + sz - i - 1
        ptr_lo = ptr + i
        val_hi = unsafe_load(ptr_hi)
        val_lo = unsafe_load(ptr_lo)
        unsafe_store!(ptr_hi, val_lo)
        unsafe_store!(ptr_lo, val_hi)
    end
end

"""
    fix_endian(x, endianness::Symbol)

Returns a byte-swapped version of `x` if the given endianness must be swapped
for the current host system.
"""
@pure function fix_endian(x, endianness::Symbol)
    if needs_bswap(endianness)
        return bswap(x)
    end
    return x
end

# Alignment traits
abstract type PackingStrategy end
struct Packed <: PackingStrategy; end
struct Default <: PackingStrategy; end

"""
    packing_strategy(x)

Return the packing strategy for the given type, defaults to `Default`, is
overridden by auto-generated methods for specific types from `@io` invocations.
"""
function packing_strategy(x)
    return Default
end

# Sizeof computation
@pure function packed_sizeof(T::DataType, ::Type{Default})
    return Core.sizeof(T)
end

@pure function packed_sizeof(T::DataType, ::Type{Packed})
    @assert fieldcount(T) != 0 && isbitstype(T)
    return sum(packed_sizeof, T.types)
end

@pure function packed_sizeof(T::DataType)
    if fieldcount(T) == 0
        return packed_sizeof(T, Default)
    else
        return packed_sizeof(T, packing_strategy(T))
    end
end

"""
    fieldsize(T::DataType, field_idx)

Return the size (in bytes) of a field within `T` in memory
"""
@pure function fieldsize(T::DataType, field_idx)
    @assert fieldcount(T) != 0 && isbitstype(T)
    @assert field_idx <= fieldcount(T)

    # We figure out the (padded) size of the given field by looking at the
    # offset of the next field (if it exists) or just the overall size of the
    # parent
    offset = fieldoffset(T, field_idx)
    if field_idx == fieldcount(T)
        # If there are no further fields, use the total size of the object
        return Core.sizeof(T) - offset
    else
        # If there are fields, diff this one with the next
        return fieldoffset(T, field_idx + 1) - offset
    end
end

"""
    @io <type definition>
        ...
    end

Generates `packing_strategy()` and `packed_sizeof()` methods for the type being
defined within the given type definition.  This enables usage of the `unpack`
method.
"""
macro io(typ, annotations...)
    alignment = :align_default
    if length(annotations) == 1
        ann = annotations[1]
        if isa(ann, Symbol) || haskey(alignments, ann)
            alignment = ann
        end
    end

    # Get typename, collapsing type expressions until we get the actual type
    T = typ.args[2]
    if isexpr(T, :(<:))
        T = T.args[1]
    end
    if isexpr(T, :curly)
        T = T.args[1]
    end

    ret = Expr(:toplevel, :(Base.@__doc__ $(typ)))
    strat = (alignment == :align_default ? StructIO.Default : StructIO.Packed)
    push!(ret.args, :(StructIO.packing_strategy(::Type{T}) where {T <: $T} = $strat))
    return esc(ret)
end

"""
    unsafe_unpack(io, T, target, endianness, ::Type{Default})

Unpack an object of type `T` from `io` into `target`, byte-swapping if
`endianness` dictates we should, assuming a `Default` packing strategy.  All
packed structs recurse until bitstypes objects are eventually reached, at which
point `Default` packing is the only behavior.
"""
function unsafe_unpack(io, ::Type{T}, target, endianness, ::Type{Default}) where {T}
    sz = Core.sizeof(T)

    if !needs_bswap(endianness)
        # If we don't need to bswap, just read directly into `target`
        unsafe_read(io, target, sz)
    elseif fieldcount(T) == 0
        # If this is a primitive data type, unpack it directly and bswap()
        unsafe_read(io, target, sz)

        # Special case small sizes, LLVM should turn this into a jump table
        if sz == 1
        elseif sz == 2
            ptr = Base.unsafe_convert(Ptr{T}, target)
            unsafe_store!(ptr, bswap(unsafe_load(ptr)))
        elseif sz == 4
            ptr = Base.unsafe_convert(Ptr{T}, target)
            unsafe_store!(ptr, bswap(unsafe_load(ptr)))
        elseif sz == 8
            ptr = Base.unsafe_convert(Ptr{T}, target)
            unsafe_store!(ptr, bswap(unsafe_load(ptr)))
        else
            # Otherwise, for large primitive objects, fall back to our
            # `bswap!()` method which will swap in-place
            void_ptr = Base.unsafe_convert(Ptr{Cvoid}, target)
            bswap!(Base.unsafe_convert(Ptr{UInt8}, void_ptr), sz)
        end
    else
        # If we need to bswap, but it's not a primitive type, recurse!
        target_ptr = Base.unsafe_convert(Ptr{Cvoid}, target)
        for i = 1:fieldcount(T)
            # Unpack this field into `target` at the appropriate offset
            fT = fieldtype(T, i)
            target_i = target_ptr + fieldoffset(T, i)

            # Unpack from this point in the IOStream into this field
            unsafe_unpack(io, fT, target_i, endianness, Default)

            # If bytes_read != Core.sizeof(fT), move it on forward
            skip(io, fieldsize(T, i) - Core.sizeof(fT))
        end
    end
end

"""
    unsafe_pack(io, source, endianness, ::Type{Packed/Default})

Pack `source` into `io`, byte-swapping if `endianness` dictates we should.  The
last argument is a packing strategy, used to determine the layout of the data
in memory.  All `Packed` objects recurse until bitstypes objects are eventually
reached, at which point `Default` packing is identical to `Packed` behavior.
"""
function unsafe_pack(io, source::Ref{T}, endianness, ::Type{Default}) where {T}
    sz = packed_sizeof(T)
    if !needs_bswap(endianness)
        # If we don't need to bswap, just write directly from `source`
        unsafe_write(io, source, sz)
    elseif fieldcount(T) == 0
        # Hopefully, LLVM turns this into a jump list for us
        @GC.preserve source if sz == 1
            write(io, source[])
        elseif sz == 2
            ptr = Base.unsafe_convert(Ptr{T}, source)
            write(io, bswap(unsafe_load(ptr)))
        elseif sz == 4
            ptr = Base.unsafe_convert(Ptr{T}, source)
            write(io, bswap(unsafe_load(ptr)))
        elseif sz == 8
            ptr = Base.unsafe_convert(Ptr{T}, source)
            write(io, bswap(unsafe_load(ptr)))
        else
            # If we must bswap something of unknown size, copy first so as
            # to not clobber `source`, then bswap, then write
            source_copy = Ref{T}(copy(source[]))
            @GC.preserve source_copy begin
                void_ptr = Base.unsafe_convert(Ptr{Cvoid}, source_copy)
                ptr = Base.unsafe_convert(Ptr{UInt8}, void_ptr)
                bswap!(ptr, sz)
                unsafe_write(io, ptr, sz)
            end
        end
    else
        # If we need to bswap, but it's not a primitive type, recurse!
        for i = 1:fieldcount(T)
            # Pack field `i` into `io`
            f = getfield(source, fieldname(source, i))
            unsafe_pack(io, f, endianness, Default)
        end
    end
end

# `Packed` packing strategy override for `unsafe_unpack`
function unsafe_unpack(io, T, target, endianness, ::Type{Packed})
    # If this type cannot be subdivided, packing strategy means nothing, so
    # hand it off to the `Default` packing strategy method
    if fieldcount(T) == 0
        return unsafe_unpack(io, T, target, endianness, Default)
    end

    # Otherwise, iterate over the fields, unpacking each into `target`
    target_ptr = Base.unsafe_convert(Ptr{Cvoid}, target)
    for i = 1:fieldcount(T)
        # Unpack this field into `target` at the appropriate offset
        fT = fieldtype(T, i)
        target_i = target_ptr + fieldoffset(T, i)
        unsafe_unpack(io, fT, target_i, endianness, Packed)
    end
end

# `Packed` packing strategy override for `unsafe_pack`
function unsafe_pack(io, source::Ref{T}, endianness, ::Type{Packed}) where {T}
    # If this type cannot be subdivided, packing strategy means nothing, so
    # hand it off to the `Default` packing strategy method
    if fieldcount(T) == 0
        return unsafe_pack(io, source, endianness, Default)
    end

    # Otherwise, iterate over the fields, packing each into `io`
    for i = 1:fieldcount(T)
        # Unpack this field into `target` at the appropriate offset
        fT = fieldtype(T, i)
        f = Ref{fT}(getfield(source[], fieldname(T, i)))
        unsafe_pack(io, f, endianness, Packed)
    end
end

"""
    unpack(io::IO, T::Type, endianness::Symbol = :NativeEndian)

Given an input `io`, unpack type `T`, byte-swapping according to the given
`endianness` of `io`. If `endianness` is `:NativeEndian` (the default), no
byteswapping will occur.  If `endianness` is `:LittleEndian` or `:BigEndian`,
byteswapping will occur of the endianness if the host system does not match
the endianness of `io`.
"""
function unpack(io::IO, T::Type, endianness::Symbol = :NativeEndian)
    # Create a `Ref{}` pointing to type T, we'll unpack into that
    r = Ref{T}()
    packstrat = fieldcount(T) == 0 ? Default : packing_strategy(T)
    unsafe_unpack(io, T, r, endianness, packstrat)

    # De-reference `r` and return its unpacked contents
    return r[]
end

"""
    pack(io::IO, source, endianness::Symbol = :NativeEndian)

Given an input `source`, pack it into `io`, byte-swapping according to the
given `endianness` of `io`. If `endianness` is `:NativeEndian` (the default),
no byteswapping will occur.  If `endianness` is `:LittleEndian` or
`:BigEndian`, byteswapping will occur if the endianness of the host system
does not match the endianness of `io`.
"""
function pack(io::IO, source::T, endianness::Symbol = :NativeEndian) where {T}
    r = Ref{T}(source)
    packstrat = fieldcount(T) == 0 ? Default : packing_strategy(T)
    unsafe_pack(io, r, endianness, packstrat)
    return nothing
end

end # module


text:
function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_absolute_accuracy(obj::BaseAbstractUnivariateSolver)
    return jcall(obj, "getAbsoluteAccuracy", jdouble, ())
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_evaluations(obj::BaseAbstractUnivariateSolver)
    return jcall(obj, "getEvaluations", jint, ())
end

function get_function_value_accuracy(obj::BaseAbstractUnivariateSolver)
    return jcall(obj, "getFunctionValueAccuracy", jdouble, ())
end

function get_max(obj::BaseAbstractUnivariateSolver)
    return jcall(obj, "getMax", jdouble, ())
end

function get_max_evaluations(obj::BaseAbstractUnivariateSolver)
    return jcall(obj, "getMaxEvaluations", jint, ())
end

function get_min(obj::BaseAbstractUnivariateSolver)
    return jcall(obj, "getMin", jdouble, ())
end

function get_relative_accuracy(obj::BaseAbstractUnivariateSolver)
    return jcall(obj, "getRelativeAccuracy", jdouble, ())
end

function get_start_value(obj::BaseAbstractUnivariateSolver)
    return jcall(obj, "getStartValue", jdouble, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function solve(obj::BaseAbstractUnivariateSolver, arg0::jint, arg1::UnivariateFunction, arg2::jdouble)
    return jcall(obj, "solve", jdouble, (jint, UnivariateFunction, jdouble), arg0, arg1, arg2)
end

function solve(obj::BaseAbstractUnivariateSolver, arg0::jint, arg1::UnivariateFunction, arg2::jdouble, arg3::jdouble)
    return jcall(obj, "solve", jdouble, (jint, UnivariateFunction, jdouble, jdouble), arg0, arg1, arg2, arg3)
end

function solve(obj::BaseAbstractUnivariateSolver, arg0::jint, arg1::UnivariateFunction, arg2::jdouble, arg3::jdouble, arg4::jdouble)
    return jcall(obj, "solve", jdouble, (jint, UnivariateFunction, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4)
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


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 `FMI3_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)

"""
Source: FMISpec3.0, Version D5ef1c1: 2.3.2. State: Instantiated

FMU enters Initialization mode.

For more information call ?fmi3EnterInitializationMode
"""
function fmi3EnterInitializationMode(c::FMU3Instance, startTime::Real = 0.0, stopTime::Real = startTime; tolerance::Real = 0.0)
    if c.state != fmi3InstanceStateInstantiated
        @warn "fmi3EnterInitializationMode(...): Needs to be called in state `fmi3IntanceStateInstantiated`."
    end
    c.t = startTime

    toleranceDefined = (tolerance > 0.0)
    stopTimeDefined = (stopTime > startTime)

    status = fmi3EnterInitializationMode(c.fmu.cEnterInitializationMode, c.compAddr, fmi3Boolean(toleranceDefined), fmi3Float64(tolerance), fmi3Float64(startTime), fmi3Boolean(stopTimeDefined), fmi3Float64(stopTime))
    checkStatus(c, status)
    if status == fmi3StatusOK
        c.state = fmi3InstanceStateInitializationMode
    end
    return status
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.3.1. Super State: FMU State Setable

Set the DebugLogger for the FMU.
"""
function fmi3SetDebugLogging(c::FMU3Instance)
    fmi3SetDebugLogging(c, fmi3False, Unsigned(0), C_NULL)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Get the values of an array of fmi3Float32 variables.

For more information call ?fmi3GetFloat32
"""
function fmi3GetFloat32(c::FMU3Instance, vr::fmi3ValueReferenceFormat)

    vr = prepareValueReference(c, vr)

    nvr = Csize_t(length(vr))
    values = zeros(fmi3Float32, nvr)
    fmi3GetFloat32!(c, vr, nvr, values, nvr)

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

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Get the values of an array of fmi3Float32 variables.

For more information call ?fmi3GetFloat32!
"""
function fmi3GetFloat32!(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::AbstractArray{fmi3Float32})

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

    nvr = Csize_t(length(values))
    fmi3GetFloat32!(c, vr, nvr, values, nvr)
    nothing
end
function fmi3GetFloat32!(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::fmi3Float32)
    @assert false "fmi3GetFloat32! is only possible for arrays of values, please use an array instead of a scalar."
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Set the values of an array of fmi3Float32 variables.

For more information call ?fmi3SetFloat32
"""
function fmi3SetFloat32(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::Union{AbstractArray{fmi3Float32}, fmi3Float32})

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

    nvr = Csize_t(length(vr))
    fmi3SetFloat32(c, vr, nvr, values, nvr)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Get the values of an array of fmi3Float64 variables.

For more information call ?fmi3GetFloat64
"""
function fmi3GetFloat64(c::FMU3Instance, vr::fmi3ValueReferenceFormat)

    vr = prepareValueReference(c, vr)

    nvr = Csize_t(length(vr))
    values = zeros(fmi3Float64, nvr)
    fmi3GetFloat64!(c, vr, nvr, values, nvr)

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

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Get the values of an array of fmi3Float64 variables.

For more information call ?fmi3GetFloat64!
"""
function fmi3GetFloat64!(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::AbstractArray{fmi3Float64})

    vr = prepareValueReference(c, vr)

    @assert length(vr) == length(values) "fmi3GetFloat64!(...): `vr` and `values` need to be the same length."

    nvr = Csize_t(length(values))
    fmi3GetFloat64!(c, vr, nvr, values, nvr)
    nothing
end
function fmi3GetFloat64!(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::fmi3Float64)
    @assert false "fmi3GetFloat64! is only possible for arrays of values, please use an array instead of a scalar."
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Set the values of an array of fmi3Float64 variables.

For more information call ?fmi3SetFloat64
"""
function fmi3SetFloat64(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::Union{AbstractArray{fmi3Float64}, fmi3Float64})

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

    nvr = Csize_t(length(vr))
    fmi3SetFloat64(c, vr, nvr, values, nvr)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Get the values of an array of fmi3Int8 variables.

For more information call ?fmi3GetInt8
"""
function fmi3GetInt8(c::FMU3Instance, vr::fmi3ValueReferenceFormat)

    vr = prepareValueReference(c, vr)

    nvr = Csize_t(length(vr))
    values = zeros(fmi3Int8, nvr)
    fmi3GetInt8!(c, vr, nvr, values, nvr)

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

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Get the values of an array of fmi3Int8 variables.

For more information call ?fmi3GetInt8!
"""
function fmi3GetInt8!(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::AbstractArray{fmi3Int8})

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

    nvr = Csize_t(length(values))
    fmi3GetInt8!(c, vr, nvr, values, nvr)
    nothing
end
function fmi3GetInt8!(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::fmi3Int8)
    @assert false "fmi3GetInt8! is only possible for arrays of values, please use an array instead of a scalar."
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Set the values of an array of fmi3Int8 variables.

For more information call ?fmi3SetInt8
"""
function fmi3SetInt8(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::Union{AbstractArray{fmi3Int8}, fmi3Int8})

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

    nvr = Csize_t(length(vr))
    fmi3SetInt8(c, vr, nvr, values, nvr)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Get the values of an array of fmi3UInt8 variables.

For more information call ?fmi3GetUInt8
"""
function fmi3GetUInt8(c::FMU3Instance, vr::fmi3ValueReferenceFormat)

    vr = prepareValueReference(c, vr)

    nvr = Csize_t(length(vr))
    values = zeros(fmi3UInt8, nvr)
    fmi3GetUInt8!(c, vr, nvr, values, nvr)

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

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Get the values of an array of fmi3UInt8 variables.

For more information call ?fmi3GetUInt8!
"""
function fmi3GetUInt8!(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::AbstractArray{fmi3UInt8})

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

    nvr = Csize_t(length(values))
    fmi3GetUInt8!(c, vr, nvr, values, nvr)
    nothing
end
function fmi3GetUInt8!(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::fmi3UInt8)
    @assert false "fmi3GetUInt8! is only possible for arrays of values, please use an array instead of a scalar."
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Set the values of an array of fmi3UInt8 variables.

For more information call ?fmi3SetUInt8
"""
function fmi3SetUInt8(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::Union{AbstractArray{fmi3UInt8}, fmi3UInt8})

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

    nvr = Csize_t(length(vr))
    fmi3SetUInt8(c, vr, nvr, values, nvr)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Get the values of an array of fmi3Int16 variables.

For more information call ?fmi3GetInt16
"""
function fmi3GetInt16(c::FMU3Instance, vr::fmi3ValueReferenceFormat)

    vr = prepareValueReference(c, vr)

    nvr = Csize_t(length(vr))
    values = zeros(fmi3Int16, nvr)
    fmi3GetInt16!(c, vr, nvr, values, nvr)

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

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Get the values of an array of fmi3Int16 variables.

For more information call ?fmi3GetInt16!
"""
function fmi3GetInt16!(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::AbstractArray{fmi3Int16})

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

    nvr = Csize_t(length(values))
    fmi3GetInt16!(c, vr, nvr, values, nvr)
    nothing
end
function fmi3GetInt16!(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::fmi3Int16)
    @assert false "fmi3GetInt16! is only possible for arrays of values, please use an array instead of a scalar."
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Set the values of an array of fmi3Int16 variables.

For more information call ?fmi3SetInt16
"""
function fmi3SetInt16(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::Union{AbstractArray{fmi3Int16}, fmi3Int16})

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

    nvr = Csize_t(length(vr))
    fmi3SetInt16(c, vr, nvr, values, nvr)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Get the values of an array of fmi3UInt16 variables.

For more information call ?fmi3GetUInt16
"""
function fmi3GetUInt16(c::FMU3Instance, vr::fmi3ValueReferenceFormat)

    vr = prepareValueReference(c, vr)

    nvr = Csize_t(length(vr))
    values = zeros(fmi3UInt16, nvr)
    fmi3GetUInt16!(c, vr, nvr, values, nvr)

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

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Get the values of an array of fmi3UInt16 variables.

For more information call ?fmi3GetUInt16!
"""
function fmi3GetUInt16!(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::AbstractArray{fmi3UInt16})

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

    nvr = Csize_t(length(values))
    fmi3GetUInt16!(c, vr, nvr, values, nvr)
    nothing
end
function fmi3GetUInt16!(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::fmi3UInt16)
    @assert false "fmi3GetUInt16! is only possible for arrays of values, please use an array instead of a scalar."
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Set the values of an array of fmi3UInt16 variables.

For more information call ?fmi3SetUInt16
"""
function fmi3SetUInt16(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::Union{AbstractArray{fmi3UInt16}, fmi3UInt16})

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

    nvr = Csize_t(length(vr))
    fmi3SetUInt16(c, vr, nvr, values, nvr)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Get the values of an array of fmi3Int32 variables.

For more information call ?fmi3GetInt32
"""
function fmi3GetInt32(c::FMU3Instance, vr::fmi3ValueReferenceFormat)

    vr = prepareValueReference(c, vr)

    nvr = Csize_t(length(vr))
    values = zeros(fmi3Int32, nvr)
    fmi3GetInt32!(c, vr, nvr, values, nvr)

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

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Get the values of an array of fmi3Int32 variables.

For more information call ?fmi3GetInt32!
"""
function fmi3GetInt32!(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::AbstractArray{fmi3Int32})

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

    nvr = Csize_t(length(values))
    fmi3GetInt32!(c, vr, nvr, values, nvr)
    nothing
end
function fmi3GetInt32!(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::fmi3Int32)
    @assert false "fmi3GetInt32! is only possible for arrays of values, please use an array instead of a scalar."
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Set the values of an array of fmi3Int32 variables.

For more information call ?fmi3SetInt32
"""
function fmi3SetInt32(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::Union{AbstractArray{fmi3Int32}, fmi3Int32})

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

    nvr = Csize_t(length(vr))
    fmi3SetInt32(c, vr, nvr, values, nvr)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Get the values of an array of fmi3UInt32 variables.

For more information call ?fmi3GetUInt32
"""
function fmi3GetUInt32(c::FMU3Instance, vr::fmi3ValueReferenceFormat)

    vr = prepareValueReference(c, vr)

    nvr = Csize_t(length(vr))
    values = zeros(fmi3UInt32, nvr)
    fmi3GetUInt32!(c, vr, nvr, values, nvr)

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

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Get the values of an array of fmi3UInt32 variables.

For more information call ?fmi3GetUInt32!
"""
function fmi3GetUInt32!(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::AbstractArray{fmi3UInt32})

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

    nvr = Csize_t(length(values))
    fmi3GetUInt32!(c, vr, nvr, values, nvr)
    nothing
end
function fmi3GetUInt32!(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::fmi3UInt32)
    @assert false "fmi3GetUInt32! is only possible for arrays of values, please use an array instead of a scalar."
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Set the values of an array of fmi3UInt32 variables.

For more information call ?fmi3SetUInt32
"""
function fmi3SetUInt32(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::Union{AbstractArray{fmi3UInt32}, fmi3UInt32})

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

    nvr = Csize_t(length(vr))
    fmi3SetUInt32(c, vr, nvr, values, nvr)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Get the values of an array of fmi3Int64 variables.

For more information call ?fmi3GetInt64
"""
function fmi3GetInt64(c::FMU3Instance, vr::fmi3ValueReferenceFormat)

    vr = prepareValueReference(c, vr)

    nvr = Csize_t(length(vr))
    values = zeros(fmi3Int64, nvr)
    fmi3GetInt64!(c, vr, nvr, values, nvr)

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

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Get the values of an array of fmi3Int64 variables.

For more information call ?fmi3GetInt64!
"""
function fmi3GetInt64!(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::AbstractArray{fmi3Int64})

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

    nvr = Csize_t(length(values))
    fmi3GetInt64!(c, vr, nvr, values, nvr)
    nothing
end
function fmi3GetInt64!(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::fmi3Int64)
    @assert false "fmi3GetInt64! is only possible for arrays of values, please use an array instead of a scalar."
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Set the values of an array of fmi3Int64 variables.

For more information call ?fmi3SetInt64
"""
function fmi3SetInt64(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::Union{AbstractArray{fmi3Int64}, fmi3Int64})

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

    nvr = Csize_t(length(vr))
    fmi3SetInt64(c, vr, nvr, values, nvr)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Get the values of an array of fmi3UInt64 variables.

For more information call ?fmi3GetUInt64
"""
function fmi3GetUInt64(c::FMU3Instance, vr::fmi3ValueReferenceFormat)

    vr = prepareValueReference(c, vr)

    nvr = Csize_t(length(vr))
    values = zeros(fmi3UInt64, nvr)
    fmi3GetUInt64!(c, vr, nvr, values, nvr)

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

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Get the values of an array of fmi3UInt64 variables.

For more information call ?fmi3GetUInt64!
"""
function fmi3GetUInt64!(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::AbstractArray{fmi3UInt64})

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

    nvr = Csize_t(length(values))
    fmi3GetUInt64!(c, vr, nvr, values, nvr)
    nothing
end
function fmi3GetUInt64!(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::fmi3UInt64)
    @assert false "fmi3GetUInt64! is only possible for arrays of values, please use an array instead of a scalar."
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Set the values of an array of fmi3UInt64 variables.

For more information call ?fmi3SetUInt64
"""
function fmi3SetUInt64(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::Union{AbstractArray{fmi3UInt64}, fmi3UInt64})

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

    nvr = Csize_t(length(vr))
    fmi3SetUInt64(c, vr, nvr, values, nvr)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Get the values of an array of fmi3Boolean variables.

For more information call ?fmi3GetBoolean
"""
function fmi3GetBoolean(c::FMU3Instance, vr::fmi3ValueReferenceFormat)

    vr = prepareValueReference(c, vr)

    nvr = Csize_t(length(vr))
    values = Array{fmi3Boolean}(undef, nvr)
    fmi3GetBoolean!(c, vr, nvr, values, nvr)

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

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Get the values of an array of fmi3Boolean variables.

For more information call ?fmi3GetBoolean!
"""
function fmi3GetBoolean!(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::AbstractArray{fmi3Boolean})

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

    nvr = Csize_t(length(values))
    fmi3GetBoolean!(c, vr, nvr, values, nvr)

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

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Set the values of an array of fmi3Boolean variables.

For more information call ?fmi3SetBoolean
"""
function fmi3SetBoolean(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::Union{AbstractArray{Bool}, Bool})

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

    nvr = Csize_t(length(vr))
    fmi3SetBoolean(c, vr, nvr, Array{fmi3Boolean}(values), nvr)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Get the values of an array of fmi3String variables.

For more information call ?fmi3GetString
"""
function fmi3GetString(c::FMU3Instance, vr::fmi3ValueReferenceFormat)

    vr = prepareValueReference(c, vr)

    nvr = Csize_t(length(vr))
    vars = Vector{Ptr{Cchar}}(undef, nvr)
    values = string.(zeros(nvr))
    fmi3GetString!(c, vr, nvr, vars, nvr)
    values[:] = unsafe_string.(vars)

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

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Get the values of an array of fmi3String variables.

For more information call ?fmi3GetString!
"""
function fmi3GetString!(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::AbstractArray{fmi3String})

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

    nvr = Csize_t(length(vr))
    # values = Vector{Ptr{Cchar}}.(values)
    vars = Vector{Ptr{Cchar}}(undef, nvr)
    fmi3GetString!(c, vr, nvr, vars, nvr)
    values[:] = unsafe_string.(vars)
    nothing
end
function fmi3GetString!(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::String)
    @assert false "fmi3GetString! is only possible for arrays of values, please use an array instead of a scalar."
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Set the values of an array of fmi3String variables.

For more information call ?fmi3SetString
"""
function fmi3SetString(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::Union{AbstractArray{String}, String})

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

    nvr = Csize_t(length(vr))
    ptrs = pointer.(values)
    fmi3SetString(c, vr, nvr, ptrs, nvr)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Get the values of an array of fmi3Binary variables.

For more information call ?fmi3GetBinary
"""
function fmi3GetBinary(c::FMU3Instance, vr::fmi3ValueReferenceFormat)

    vr = prepareValueReference(c, vr)

    nvr = Csize_t(length(vr))
    values = Array{fmi3Binary}(undef, nvr)
    valueSizes = Array{Csize_t}(undef, nvr)
    fmi3GetBinary!(c, vr, nvr, valueSizes, values, nvr)
    if length(values) == 1
        return values[1]
    else
        return values
    end
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Get the values of an array of fmi3Binary variables.

For more information call ?fmi3GetBinary!
"""
function fmi3GetBinary!(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::AbstractArray{fmi3Binary})

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

    nvr = Csize_t(length(vr))
    valueSizes = Array{Csize_t}(undef, nvr)
    fmi3GetBinary!(c, vr, nvr, valueSizes, values, nvr)
end
function fmi3GetBinary!(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::fmi3Binary)
    @assert false "fmi3GetBinary! is only possible for arrays of values, please use an array instead of a scalar."
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Set the values of an array of fmi3Binary variables.

For more information call ?fmi3SetBinary
"""
function fmi3SetBinary(c::FMU3Instance, vr::fmi3ValueReferenceFormat, valueSizes::Union{AbstractArray{Csize_t}, Csize_t}, values::Union{AbstractArray{fmi3Binary}, fmi3Binary})

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

    nvr = Csize_t(length(vr))
    fmi3SetBinary(c, vr, nvr, valueSizes, values, nvr)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Get the values of an array of fmi3Clock variables.

For more information call ?fmi3GetClock
"""
function fmi3GetClock(c::FMU3Instance, vr::fmi3ValueReferenceFormat)

    vr = prepareValueReference(c, vr)

    nvr = Csize_t(length(vr))
    values = Array{fmi3Clock}(undef, nvr)
    fmi3GetClock!(c, vr, nvr, values)

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

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Get the values of an array of fmi3Clock variables.

For more information call ?fmi3GetClock!
"""
function fmi3GetClock!(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::AbstractArray{fmi3Clock})

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

    nvr = Csize_t(length(vr))
    fmi3GetClock!(c, vr, nvr, values)
    nothing
end
function fmi3GetClock!(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::fmi3Clock)
    @assert false "fmi3GetClock! is only possible for arrays of values, please use an array instead of a scalar."
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Set the values of an array of fmi3Clock variables.

For more information call ?fmi3SetClock
"""
function fmi3SetClock(c::FMU3Instance, vr::fmi3ValueReferenceFormat, values::Union{AbstractArray{fmi3Clock}, fmi3Clock})

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

    nvr = Csize_t(length(vr))
    fmi3SetClock(c, vr, nvr, values)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.4. Getting and Setting the Complete FMU State

Get the pointer to the current FMU state.

For more information call ?fmi3GetFMUstate
"""
function fmi3GetFMUState(c::FMU3Instance)
    state = fmi3FMUState()
    stateRef = Ref(state)
    fmi3GetFMUState!(c, stateRef)
    state = stateRef[]
    state
end

"""
function fmi3FreeFMUState(c::FMU3Instance, FMUstate::Ref{fmi3FMUState})

Free the allocated memory for the FMU state.

For more information call ?fmi3FreeFMUstate
"""
function fmi3FreeFMUState!(c::FMU3Instance, state::fmi3FMUState)
    stateRef = Ref(state)
    fmi3FreeFMUState!(c, stateRef)
    state = stateRef[]
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.4. Getting and Setting the Complete FMU State

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

For more information call ?fmi3SerzializedFMUstateSize
"""
function fmi3SerializedFMUStateSize(c::FMU3Instance, state::fmi3FMUState)
    size = 0
    sizeRef = Ref(Csize_t(size))
    fmi3SerializedFMUStateSize!(c, state, sizeRef)
    size = sizeRef[]
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.4. Getting and Setting the Complete FMU State

Serialize the data in the FMU state pointer.

For more information call ?fmi3SerzializeFMUstate
"""
function fmi3SerializeFMUState(c::FMU3Instance, state::fmi3FMUState)
    size = fmi3SerializedFMUStateSize(c, state)
    serializedState = Array{fmi3Byte}(undef, size)
    fmi3SerializeFMUState!(c, state, serializedState, size)
    serializedState
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.4. Getting and Setting the Complete FMU State

Deserialize the data in the serializedState fmi3Byte field.

For more information call ?fmi3DeSerzializeFMUstate
"""
function fmi3DeSerializeFMUState(c::FMU3Instance, serializedState::AbstractArray{fmi3Byte})
    size = length(serializedState)
    state = fmi3FMUState()
    stateRef = Ref(state)
    fmi3DeSerializeFMUState!(c, serializedState, Csize_t(size), stateRef)
    state = stateRef[]
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.11. Getting Partial Derivatives

Computes directional derivatives.

For more information call ?fmi3GetDirectionalDerivative
"""
function fmi3GetDirectionalDerivative(c::FMU3Instance,
                                      unknowns::AbstractArray{fmi3ValueReference},
                                      knowns::AbstractArray{fmi3ValueReference},
                                      seed::AbstractArray{fmi3Float64} = Array{fmi3Float64}([]))
    sensitivity = zeros(fmi3Float64, length(unknowns))

    fmi3GetDirectionalDerivative!(c, unknowns, knowns, sensitivity, seed)

    sensitivity
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.11. Getting Partial Derivatives

Computes directional derivatives.

For more information call ?fmi3GetDirectionalDerivative
"""
function fmi3GetDirectionalDerivative!(c::FMU3Instance,
                                      unknowns::AbstractArray{fmi3ValueReference},
                                      knowns::AbstractArray{fmi3ValueReference},
                                      sensitivity::AbstractArray,
                                      seed::AbstractArray{fmi3Float64}= Array{fmi3Float64}([]))

    nKnowns = Csize_t(length(knowns))
    nUnknowns = Csize_t(length(unknowns))

    if length(seed) == 0
        seed = ones(fmi3Float64, nKnowns)
    end

    nSeed = Csize_t(length(seed))
    nSensitivity = Csize_t(length(sensitivity))

    fmi3GetDirectionalDerivative!(c, unknowns, nUnknowns, knowns, nKnowns, seed, nSeed, sensitivity, nSensitivity)

    nothing
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.11. Getting Partial Derivatives

Computes directional derivatives.

For more information call ?fmi3GetDirectionalDerivative
"""
function fmi3GetDirectionalDerivative(c::FMU3Instance,
                                      unknown::fmi3ValueReference,
                                      known::fmi3ValueReference,
                                      seed::fmi3Float64 = 1.0)

    fmi3GetDirectionalDerivative(c, [unknown], [known], [seed])[1]
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.11. Getting Partial Derivatives

Computes adjoint derivatives.

For more information call ?fmi3GetAdjointDerivative
"""
function fmi3GetAdjointDerivative(c::FMU3Instance,
                                      unknowns::fmi3ValueReference,
                                      knowns::fmi3ValueReference,
                                      seed::fmi3Float64 = 1.0)

    fmi3GetAdjointDerivative(c, [unknowns], [knowns], [seed])[1]
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.11. Getting Partial Derivatives

Computes adjoint derivatives.

For more information call ?fmi3GetAdjointDerivative
"""
function fmi3GetAdjointDerivative(c::FMU3Instance,
                                      unknowns::AbstractArray{fmi3ValueReference},
                                      knowns::AbstractArray{fmi3ValueReference},
                                      seed::AbstractArray{fmi3Float64} = Array{fmi3Float64}([]))
    sensitivity = zeros(fmi3Float64, length(unknowns))

    fmi3GetAdjointDerivative!(c, unknowns, knowns, sensitivity, seed)

    sensitivity
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.11. Getting Partial Derivatives

Computes adjoint derivatives.

For more information call ?fmi3GetAdjointDerivative
"""
function fmi3GetAdjointDerivative!(c::FMU3Instance,
                                      unknowns::AbstractArray{fmi3ValueReference},
                                      knowns::AbstractArray{fmi3ValueReference},
                                      sensitivity::AbstractArray,
                                      seed::AbstractArray{fmi3Float64}= Array{fmi3Float64}([]))

    nKnowns = Csize_t(length(knowns))
    nUnknowns = Csize_t(length(unknowns))

    if length(seed) == 0
        seed = ones(fmi3Float64, nKnowns)
    end

    nSeed = Csize_t(length(seed))
    nSensitivity = Csize_t(length(sensitivity))

    fmi3GetAdjointDerivative!(c, unknowns, nUnknowns, knowns, nKnowns, seed, nSeed, sensitivity, nSensitivity)

    nothing
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.12. Getting Derivatives of Continuous Outputs

Retrieves the n-th derivative of output values.

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

For more information call ?fmi3GetOutputDerivatives
"""
function fmi3GetOutputDerivatives(c::FMU3Instance, vr::fmi3ValueReferenceFormat, order::AbstractArray{Integer})
    vr = prepareValueReference(c, vr)
    nvr = Csize_t(length(vr))
    values = Array{fmi3Float64}(undef, nvr)
    fmi3GetOutputDerivatives!(c, vr, nvr, fmi3Int32.(order), values, nvr)
    values
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.12. Getting Derivatives of Continuous Outputs

Retrieves the n-th derivative of output values.

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

For more information call ?fmi3GetOutputDerivatives
"""
function fmi3GetOutputDerivatives(c::FMU3Instance, vr::fmi3ValueReference, order::Integer)
    nvr = Csize_t(1)
    values = Array{fmi3Float64}(undef, nvr)
    fmi3GetOutputDerivatives!(c, [vr], nvr, [order], values, nvr)
    values
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.3.2. State: Instantiated

This function returns the number of continuous states.
This function can only be called in Model Exchange.
For more information call ?fmi3GetNumberOfContinuousStates
"""
function fmi3GetNumberOfContinuousStates(c::FMU3Instance)
    size = 0
    sizeRef = Ref(Csize_t(size))
    fmi3GetNumberOfContinuousStates!(c, sizeRef)
    size = sizeRef[]
    Int32(size)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.3.2. State: Instantiated

This function returns the number of event indicators.
This function can only be called in Model Exchange.
For more information call ?fmi3GetNumberOfEventIndicators
"""
function fmi3GetNumberOfEventIndicators(c::FMU3Instance)
    size = 0
    sizeRef = Ref(Csize_t(size))
    fmi3GetNumberOfEventIndicators!(c, sizeRef)
    size = sizeRef[]
    Int32(size)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.10. Dependencies of Variables

The number of dependencies of a given variable, which may change if structural parameters are changed, can be retrieved by calling the following function:
For more information call ?fmi3GetNumberOfVariableDependencies
"""
function fmi3GetNumberOfVariableDependencies(c::FMU3Instance, vr::Union{fmi3ValueReference, String})
    if typeof(vr) == String
        vr = fmi3String2ValueReference(c.fmu.modelDescription, vr)
    end
    size = 0
    sizeRef = Ref(Csize_t(size))
    fmi3GetNumberOfVariableDependencies!(c, vr, sizeRef)
    size = sizeRef[]
    Int32(size)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.10. Dependencies of Variables

The actual dependencies (of type dependenciesKind) can be retrieved by calling the function fmi3GetVariableDependencies:
For more information call ?fmi3GetVariableDependencies
"""
function fmi3GetVariableDependencies(c::FMU3Instance, vr::Union{fmi3ValueReference, String})
    if typeof(vr) == String
        vr = fmi3String2ValueReference(c.fmu.modelDescription, vr)
    end
    nDependencies = fmi3GetNumberOfVariableDependencies(c, vr)
    elementIndiceOfDependents = Array{Csize_t}(undef, nDependencies)
    independents = Array{fmi3ValueReference}(undef, nDependencies)
    elementIndiceOfIndependents = Array{Csize_t}(undef, nDependencies)
    dependencyKinds = Array{fmi3DependencyKind}(undef, nDependencies)

    fmi3GetVariableDependencies!(c, vr, elementIndiceOfDependents, independents, elementIndiceOfIndependents, dependencyKinds, nDependencies)

    return elementIndiceOfDependents, independents, elementIndiceOfIndependents, dependencyKinds
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.3.3. State: Initialization Mode

Return the new (continuous) state vector x.

For more information call ?fmi3GetContinuousStates
"""
function fmi3GetContinuousStates(c::FMU3Instance)
    nx = Csize_t(c.fmu.modelDescription.numberOfContinuousStates)
    x = zeros(fmi3Float64, nx)
    fmi3GetContinuousStates!(c, x, nx)
    x
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.3.3. State: Initialization Mode

Return the new (continuous) state vector x.

For more information call ?fmi3GetNominalsOfContinuousStates
"""
function fmi3GetNominalsOfContinuousStates(c::FMU3Instance)
    nx = Csize_t(c.fmu.modelDescription.numberOfContinuousStates)
    x = zeros(fmi3Float64, nx)
    fmi3GetNominalsOfContinuousStates!(c, x, nx)
    x
end

"""
Source: FMISpec3.0, Version D5ef1c1: 3.2.1. State: Continuous-Time Mode

Set independent variable time and reinitialize chaching of variables that depend on time.

For more information call ?fmi3SetTime
"""
function fmi3SetTime(c::FMU3Instance, time::Real)
    fmi3SetTime(c, fmi3Float64(time))
end

"""
Source: FMISpec3.0, Version D5ef1c1: 3.2.1. State: Continuous-Time Mode

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

For more information call ?fmi3SetContinuousStates
"""
function fmi3SetContinuousStates(c::FMU3Instance, x::Union{AbstractArray{Float32}, AbstractArray{Float64}})
    nx = Csize_t(length(x))
    fmi3SetContinuousStates(c, Array{fmi3Float64}(x), nx)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 3.2.1. State: Continuous-Time Mode

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

For more information call ?fmi3GetContinuousDerivatives
"""
function  fmi3GetContinuousStateDerivatives(c::FMU3Instance)
    nx = Csize_t(c.fmu.modelDescription.numberOfContinuousStates)
    derivatives = zeros(fmi3Float64, nx)
    fmi3GetContinuousStateDerivatives!(c, derivatives, nx)
    derivatives
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.3.5. State: Event Mode

This function is called to signal a converged solution at the current super-dense time instant. fmi3UpdateDiscreteStates must be called at least once per super-dense time instant.

For more information call ?fmi3UpdateDiscreteStates
"""
function fmi3UpdateDiscreteStates(c::FMU3Instance)
    discreteStatesNeedUpdate = fmi3True
    terminateSimulation = fmi3True
    nominalsOfContinuousStatesChanged = fmi3True
    valuesOfContinuousStatesChanged = fmi3True
    nextEventTimeDefined = fmi3True
    nextEventTime = fmi3Float64(1.0)
    refdS = Ref(discreteStatesNeedUpdate)
    reftS = Ref(terminateSimulation)
    refnOCS = Ref(nominalsOfContinuousStatesChanged)
    refvOCS = Ref(valuesOfContinuousStatesChanged)
    refnETD = Ref(nextEventTimeDefined)
    refnET = Ref(nextEventTime)

    fmi3UpdateDiscreteStates!(c, refdS, reftS, refnOCS, refvOCS, refnETD, refnET)

    discreteStatesNeedUpdate = refdS[]
    terminateSimulation = reftS[]
    nominalsOfContinuousStatesChanged =refnOCS[]
    valuesOfContinuousStatesChanged = refvOCS[]
    nextEventTimeDefined = refnETD[]
    nextEventTime = refnET[]

    discreteStatesNeedUpdate, terminateSimulation, nominalsOfContinuousStatesChanged, valuesOfContinuousStatesChanged, nextEventTimeDefined, nextEventTime
end

"""
Source: FMISpec3.0, Version D5ef1c1: 3.2.1. State: Continuous-Time Mode

Returns the event indicators of the FMU.

For more information call ?fmi3GetEventIndicators
"""
function fmi3GetEventIndicators(c::FMU3Instance)
    ni = Csize_t(c.fmu.modelDescription.numberOfEventIndicators)
    eventIndicators = zeros(fmi3Float64, ni)
    fmi3GetEventIndicators!(c, eventIndicators, ni)
    eventIndicators
end

"""
Source: FMISpec3.0, Version D5ef1c1: 3.2.1. State: Continuous-Time Mode

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

For more information call ?fmi3CompletedIntegratorStep
"""
function fmi3CompletedIntegratorStep(c::FMU3Instance,
    noSetFMUStatePriorToCurrentPoint::fmi3Boolean)
    enterEventMode = fmi3Boolean(true)
    terminateSimulation = fmi3Boolean(true)
    refEventMode = Ref(enterEventMode)
    refterminateSimulation = Ref(terminateSimulation)
    status = fmi3CompletedIntegratorStep!(c,
                                        noSetFMUStatePriorToCurrentPoint,
                                        refEventMode,
                                        refterminateSimulation)
    enterEventMode = refEventMode[]
    terminateSimulation = refterminateSimulation[]

    return (status, enterEventMode, terminateSimulation)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 3.2.1. State: Continuous-Time Mode

The model enters Event Mode.

For more information call ?fmi3EnterEventMode
"""
function fmi3EnterEventMode(c::FMU3Instance, stepEvent::Bool, stateEvent::Bool, rootsFound::AbstractArray{fmi3Int32}, nEventIndicators::Csize_t, timeEvent::Bool)
    fmi3EnterEventMode(c, fmi3Boolean(stepEvent), fmi3Boolean(stateEvent), rootsFound, nEventIndicators, fmi3Boolean(timeEvent))
end

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

TODO: Add this command in the documentation.
"""
function fmi3GetStartValue(c::FMU3Instance, vrs::fmi3ValueReferenceFormat)

    vrs = prepareValueReference(c, vrs)

    starts = []

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

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

        push!(starts, mvs[1].datatype.start)
    end

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

"""
This function samples the directional derivative by manipulating corresponding values (central differences).
"""
function fmi3SampleDirectionalDerivative(c::FMU3Instance,
                                       vUnknown_ref::AbstractArray{fmi3ValueReference},
                                       vKnown_ref::AbstractArray{fmi3ValueReference},
                                       steps::AbstractArray{fmi3Float64} = ones(fmi3Float64, length(vKnown_ref)).*1e-5)

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

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

    dvUnknown
end

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

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

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

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

        fmi3SetFloat64(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