text:
using LowRankApprox: pqrfact, psvdfact, LRAOptions

export pQR_solver, pSVD_solver, rSVD_solver
LinearAlgebra.ishermitian(::DictionaryOperator) = false
function pQR_solver(A::DictionaryOperator{T}; threshold=default_threshold(A), lraoptions::LRAOptions=LRAOptions(T), verbose=false, options...) where {T}
    lraoptions.rtol = threshold
    lraoptions.verb = verbose
    verbose && println("Calculating partial QR factorization (LowRankApprox.jl) with `LRAOptions` $lraoptions")
    pqr = pqrfact(A,lraoptions)
    verbose && println("\t Solver truncated at R = ", size(pqr.Q,2), " dof out of ",size(A,2))
    BasisFunctions.GenericSolverOperator(A, pqr)
end
function pSVD_solver(A::DictionaryOperator{T}; threshold=default_threshold(A), lraoptions::LRAOptions=LRAOptions(T), verbose=false, options...) where {T}
    lraoptions.rtol = threshold
    lraoptions.verb = verbose
    verbose && println("Calculating partial SVD factorization (LowRankApprox.jl) `LRAOptions` $lraoptions")
    psvd = psvdfact(A,lraoptions)
    verbose && println("\t Solver truncated at R = ", size(psvd.U,2), " dof out of ",size(A,2))
    BasisFunctions.GenericSolverOperator(A, psvd)
end
rSVD_solver(A::DictionaryOperator; options...) =
    RandomizedSvdSolver(A; options...)

"""
A Randomized Truncated SVD solver storing a decomposition of an operator A.
If A is MxN, and of rank R, the cost of constructing this solver is MR^2.

For tensor product operators it returns a decomposition of the linearized system
"""
struct RandomizedSvdSolver{T} <: DictionarySolverOperator{T}
    # Keep the original operator
    A          ::  DictionaryOperator
    # Random matrix
    W           ::  BasisFunctions.DenseMatrixOperator{T}
    # Decomposition
    Ut          ::  Array{T,2}
    Sinv        ::  Array{T,1}
    V           ::  Array{T,2}
    # For storing intermediate results when applying
    y           ::  Array{T,1}
    sy          ::  Array{T,1}
    scratch_src ::  Array{T,1}
    scratch_dest::  Array{T,1}

    smallcoefficients   :: Bool
    smallcoefficients_atol            :: Float64
    smallcoefficients_rtol            :: Float64

    function RandomizedSvdSolver{T}(A::DictionaryOperator;
                threshold = default_threshold(A), rankestimate = 5, verbose = false, smallcoefficients = false,
                smallcoefficients_atol = NaN, smallcoefficients_rtol = NaN, options...) where {T}

        W, U, S, Vt = randomizedsvd(A, threshold = threshold, rankestimate = rankestimate, verbose = verbose)
        # Small coefficients #
        # One of the tolerances can not be NaN
        # If we know the norm of the function, we can wet the atol as rtol*\|f\| and rtol to NaN to simulate a relative tolerance using an absolute tolerace.
        smallcoefficients && (@assert !isnan(smallcoefficients_atol) || !isnan(smallcoefficients_rtol))
        if verbose && smallcoefficients
             !isnan(smallcoefficients_atol) && println("Coefficients kept smaller than: $(smallcoefficients_atol)")
             !isnan(smallcoefficients_rtol) && println("Coefficients kept smaller than: $(smallcoefficients_rtol)||b||")
        end

        Wsrc = DiscreteVectorDictionary{T}(size(W,2))
        Wdest = src(A)
        WM = ArrayOperator(W, Wsrc, Wdest)
        scratch_src = zeros(T, length(dest(A)))
        scratch_dest = zeros(T, length(src(A)))
        y = zeros(T, size(Vt,2))
        sy = zeros(T, length(S))
        new(A, WM, U', S.^(-1), Vt', y, sy, scratch_src, scratch_dest, smallcoefficients, smallcoefficients_atol, smallcoefficients_rtol)
    end
end

operator(s::RandomizedSvdSolver) = s.A

RandomizedSvdSolver(A::DictionaryOperator; options...) =
    RandomizedSvdSolver{eltype(A)}(A::DictionaryOperator; options...)

apply!(s::RandomizedSvdSolver, coef_dest, coef_src) = _apply!(s, coef_dest, coef_src,
    s.W, s.Ut, s.Sinv, s.V, s.y, s.sy, s.scratch_src, s.scratch_dest,
    s.smallcoefficients, s.smallcoefficients_atol, s.smallcoefficients_rtol)

# We don't need to (de)linearize coefficients when they are already vectors
function _apply!(s::RandomizedSvdSolver, coef_dest::Vector, coef_src::Vector,
    W, Ut, Sinv, V, y, sy, scratch_src, scratch_dest, smallcoefficients, smallcoefficients_atol, smallcoefficients_rtol)

    # Applying the truncated SVD to P*Rhs
    mul!(sy, Ut, coef_src)
    for i in 1:length(sy)
        sy[i] = sy[i]*Sinv[i]
    end
    if smallcoefficients
        if !isnan(smallcoefficients_rtol)
            M = maximum(abs.(coef_src))
            sy[abs.(sy) .> smallcoefficients_rtol*M] .= 0
        end
        if !isnan(smallcoefficients_atol)
            sy[abs.(sy) .> smallcoefficients_atol] .= 0
        end
    end
    mul!(y, V, sy)
    apply!(W, coef_dest, y)
end

function _apply!(s::RandomizedSvdSolver, coef_dest, coef_src,
    W, Ut, Sinv, V, y, sy, scratch_src, scratch_dest, smallcoefficients, smallcoefficients_atol, smallcoefficients_rtol)

    linearize_coefficients!(src(s), scratch_src, coef_src)
    # Call the implementation above for vectors
    _apply!(s, scratch_dest, scratch_src, W, Ut, Sinv, V, y, sy,
        scratch_src, scratch_dest, smallcoefficients, smallcoefficients_atol, smallcoefficients_rtol)
    delinearize_coefficients!(dest(s), coef_dest, scratch_dest)
end

"""
A Truncated SVD solver storing a decomposition of an operator A.
The cost is equal to the computation of the SVD of A.

For tensor product operators it returns a decomposition of the linearized system
"""
struct ExactSvdSolver{T} <: DictionarySolverOperator{T}
    # Keep the original operator
    A          ::  DictionaryOperator
    # Decomposition
    Ut          ::  Array{T,2}
    Sinv        ::  Array{T,1}
    V           ::  Array{T,2}
    # For storing intermediate results when applying
    y           ::  Array{T,1}
    sy          ::  Array{T,1}
    scratch_src ::  Array{T,1}

    function ExactSvdSolver{T}(A::DictionaryOperator; threshold = default_threshold(A), verbose = false,options...) where {T}
        C = matrix(A)
        m = maximum(abs.(C))
        USV = LAPACK.gesdd!('S',C)
        S = USV[2]
        maxind = findlast(S.>threshold*m)
        verbose && println("Solver truncated at singular value equal to ", S[maxind], " threshold was ", threshold, " and norm of A ", m)
        Sinv = 1 ./ S[1:maxind]
        y = zeros(T, size(USV[3],1))
        sy = zeros(T, maxind)

        scratch_src = zeros(T, length(dest(A)))
        new(A, USV[1][:,1:maxind]',Sinv,USV[3][1:maxind,:]',y,sy, scratch_src)
    end
end

operator(s::ExactSvdSolver) = s.A

ExactSvdSolver(A::DictionaryOperator; options...) =
    ExactSvdSolver{eltype(A)}(A::DictionaryOperator; options...)

# We don't need to (de)linearize coefficients when they are already vectors
function apply!(s::ExactSvdSolver, coef_dest::Vector, coef_src::Vector)
    # Applying the truncated SVD to P*Rhs
    mul!(s.sy, s.Ut, coef_src)
    for i =1:length(s.sy)
        s.sy[i]=s.sy[i]*s.Sinv[i]
    end
    mul!(coef_dest, s.V, s.sy)
end

function apply!(s::ExactSvdSolver, coef_dest, coef_src)
    linearize_coefficients!(src(s), s.scratch_src, coef_src)
    mul!(s.sy, s.Ut, s.scratch_src)
    for i =1:length(s.sy)
        s.sy[i]=s.sy[i]*s.Sinv[i]
    end
    mul!(s.y, s.V, s.sy)
    delinearize_coefficients!(dest(s), coef_dest, s.y)
end

# Function with equal functionality, but allocating memory
function truncatedsvd_solve(b::Vector, A::DictionaryOperator;
            threshold = default_threshold(A), R = 5, growth_factor = 2, verbose = false,
            smallcoefficients=false,smallcoefficients_atol=10, options...)
    finished=false
    T = eltype(A)
    R = min(R, size(A,2))
    random_matrix = rand(T, size(A,2), R)
    C = apply_multiple(A, random_matrix)
    c = cond(C)
    cold = threshold
    while (c < 1/threshold) && (R<size(A,2)) && (c>cold*10)
        verbose && println("c : $c\t cold : $cold\t threshold : $threshold")
        verbose && println("Solver truncated at R = ", R, " dof out of ",size(A,2))
        R0 = R
        R = min(round(Int,growth_factor*R),size(A,2))
        extra_random_matrix = rand(T, size(A,2), R-R0)
        Cextra = apply_multiple(A, extra_random_matrix)
        random_matrix = [random_matrix extra_random_matrix]
        # Extra condition: condition number has to increase by a significant amount each step, otherwise, possibly well conditioned.
        # cold = c
        C = [C Cextra]
        c = cond(C)
    end
    verbose && println("c : $c\t cold : $cold\t threshold : $threshold")
    verbose && println("Solver truncated at R = ", R, " dof out of ",size(A,2))

    y, rr = LAPACK.gelsy!(C,b[:],threshold)
    x1 = random_matrix*y;
    reshape(x1, size(src(A))...)
end

# Make pqrfact available for BigFloat
function LowRankApprox.geqp3_adap_main!(
    A::AbstractMatrix{T}, jpvt::Vector{Int}, tau::Vector{T},
    opts::LRAOptions) where T
  LowRankApprox.chkstride1(A)
  lda = stride(A, 2)
  n   = length(jpvt)
  m   = size(A,1)
  k   = length(tau)

  maxnrm = norm(view(A,:,LinearAlgebra.indmaxcolumn(A)))

  # set pivot threshold
  ptol = max(opts.atol, opts.rtol*maxnrm)

  # unblocked factorization
  for j = 1:k

    # Find column with maximum norm in trailing submatrix
    jm = LinearAlgebra.indmaxcolumn(view(A, j:m, j:n)) + j - 1

    if jm != j
        # Flip elements in pivoting vector
        tmpp = jpvt[jm]
        jpvt[jm] = jpvt[j]
        jpvt[j] = tmpp

        # Update matrix with
        for i = 1:m
            tmp = A[i,jm]
            A[i,jm] = A[i,j]
            A[i,j] = tmp
        end
    end

    # Compute reflector of columns j
    x = view(A, j:m, j)
    τj = LinearAlgebra.reflector!(x)
    tau[j] = τj

    # Update trailing submatrix with reflector
    LinearAlgebra.reflectorApply!(x, τj, view(A, j:m, j+1:n))

    # check for rank termination
    if abs(A[j,j]) <= ptol
        return j
    end
    j += 1
  end
  k
end

LowRankApprox.crandn(::Type{T}, dims::Integer...) where {T<:Complex} =
  reshape(reinterpret(T, vec(LowRankApprox.crandn(real(T), 2*dims[1], dims[2:end]...))),
          (dims...,))
LowRankApprox.crandn(::Type{T}, dims::Integer...) where {T<:BigFloat} =
    error("No `randn` sampling for BigFloat")
LowRankApprox.crandn(::Type{T}, dims::Integer...) where {T<:Complex{BigFloat}} =
    error("No `randn` sampling for BigFloat")


text:
module ProjectX
using Pkg: TOML

function __init__()
    project = Base.active_project()
    if !isfile(project)
        @warn "$(project) not a file"
        return
    end
    toml = TOML.parsefile(project)
    haskey(toml, "environment") && update_environment!(dirname(project), toml["environment"])
    haskey(toml, "juliaenv") && update_juliaenv!(dirname(project), toml["juliaenv"])
    return
end

function update_environment!(root, d)
    for (k, v) in d
        if v isa String
            @debug "Setting env. \$$(k)=`$(v)`"
            ENV[k] = v
        elseif v isa Dict{String,T} where T
            if !haskey(v, "method")
                @error "env. \$$k missing `method`. Skipping" v
                continue
            end
            if v["method"] == "abspath"
                if !haskey(v, "path")
                    @error "env. \$$k (method `abspath`) missing `path`. Skipping" v
                    continue
                end
                ENV[k] = normpath(joinpath(root, v["path"]))
                @debug "Setting \$$(k)=`$(ENV[k])`"
            else
                @error "Unrecognized method for env. \$$k: $(v["method"]) Skipping" v
                continue
            end
        else
            @error "Invalid value type for env. \$$k: $(typeof(v)). Skipping." v
            continue
        end
    end
end

function update_juliaenv!(root, d)
    if haskey(d, "load_path")
        loadpath = map(Base.parse_load_path(d["load_path"])) do path
            startswith(path, "@") && return path
            joinpath(root, path)
        end
        empty!(LOAD_PATH)
        append!(LOAD_PATH, loadpath)
        @debug "LOAD_PATH updated based in Project.toml" LOAD_PATH
    end
    if haskey(d, "startup")
        isa(d["startup"], String) || error("juliaenv.startup is not a string")
        @debug "Running a project startup script:\n$(d["startup"])"
        for (ex, str) in parseblock(d["startup"])
            @debug "Evaluating in Main:\n$(str)"
            Main.eval(ex)
        end
    end
end

# Parsing multi-expression code snippets, borrowed with modificatoins from
# Documenter.jl (https://github.com/JuliaDocs/Documenter.jl)
# Copyright (c) 2016: Michael Hatherly
function parseblock(code::AbstractString; skip = 0, keywords = true, raise=true)
    # Drop `skip` leading lines from the code block. Needed for deprecated `{docs}` syntax.
    code = string(code, '\n')
    code = last(split(code, '\n', limit = skip + 1))
    endofstr = lastindex(code)
    results = []
    cursor = 1
    while cursor < endofstr
        # Check for keywords first since they will throw parse errors if we `parse` them.
        line = match(r"^(.*)\r?\n"m, SubString(code, cursor)).match
        keyword = Symbol(strip(line))
        (ex, ncursor) =
            # TODO: On 0.7 Symbol("") is in Docs.keywords, remove that check when dropping 0.6
            if keywords && (haskey(Docs.keywords, keyword) || keyword == Symbol(""))
                (QuoteNode(keyword), cursor + lastindex(line))
            else
                try
                    Meta.parse(code, cursor; raise=raise)
                catch err
                    @warn "ProjectX: failed to parse an expression" exception = err
                    break
                end
            end
        str = SubString(code, cursor, prevind(code, ncursor))
        if !isempty(strip(str))
            push!(results, (ex, str))
        end
        cursor = ncursor
    end
    results
end

end

ext: jl


text:
using Rotations
using GaussianProcesses
using ConstrainedDynamics
using GPR

function predictdynamics(mechanism::Mechanism, gps::Vector{<:GPE}, startobservation::CState, steps::Integer, getvω::Function; regularizer::Real = 0.)
    projectionerror = 0
    oldstates = startobservation
    setstates!(mechanism, oldstates)
    for _ in 1:steps
        obs = reshape(oldstates.state, :, 1)
        μ = [predict_y(gp, obs)[1][1] for gp in gps]
        vcurr, ωcurr = getvω(μ)
        vconst, ωconst = projectv!(vcurr, ωcurr, mechanism, regularizer=regularizer)
        projectionerror += norm(vcat(reduce(vcat, vconst .- vcurr), reduce(vcat, ωconst .-ωcurr)))
        foreachactive(updatestate!, mechanism.bodies, mechanism.Δt)  # Now at xcurr, vcurr
        oldstates = CState(mechanism)
    end
    foreachactive(updatestate!, mechanism.bodies, mechanism.Δt)  # Now at xnew, undef
    return CState(mechanism), projectionerror/steps
end

function predictdynamics(mechanism::Mechanism, startobservation::CState, steps::Int)
    setstates!(mechanism, startobservation)
    ConstrainedDynamics.simulate!(mechanism, 1:steps+1, Storage{Float64}())
    return CState(mechanism)
end

function predictdynamicsmin(mechanism::Mechanism, etype::String, gps::Vector{<:GPE}, startobservation::Vector{Float64}, steps::Integer; usesin::Bool = false)
    etype == "P1" && return _predictdynamicsp1min(mechanism, gps, startobservation, steps, usesin)
    etype == "P2" && return _predictdynamicsp2min(mechanism, gps, startobservation, steps, usesin)
    etype == "CP" && return _predictdynamicscpmin(mechanism, gps, startobservation, steps, usesin)
    etype == "FB" && return _predictdynamicsfbmin(mechanism, gps, startobservation, steps, usesin)
    throw(ArgumentError("Experiment $etype not supported!"))
end

function _predictdynamicsp1min(mechanism, gps, startobservation, steps, usesin)
    l = mechanism.bodies[1].shape.rh[2]
    θold, ωold = startobservation
    θcurr = θold + mechanism.Δt*ωold
    for _ in 1:steps
        obs = usesin ? reshape([sin(θold), cos(θold), ωold], :, 1) : reshape([θold, ωold], :, 1)
        ωcurr = predict_y(gps[1], obs)[1][1]
        θold, ωold = θcurr, ωcurr
        θcurr = θcurr + ωcurr*mechanism.Δt
    end
    q = UnitQuaternion(RotX(θcurr))
    return CState([0, 0.5l*sin(θcurr), -0.5l*cos(θcurr), q2vec(q)..., zeros(6)...])
end

function _predictdynamicsp2min(mechanism, gps, startobservation, steps, usesin)
    l1, l2 = mechanism.bodies[1].shape.xyz[3], mechanism.bodies[2].shape.xyz[3]
    θ1old, ω1old, θ2old, ω2old = startobservation
    θ1curr, θ2curr = θ1old + mechanism.Δt*ω1old, θ2old + mechanism.Δt*ω2old
    for _ in 1:steps
        obs = usesin ? reshape([sin(θ1old), cos(θ1old), ω1old, sin(θ2old), cos(θ2old), ω2old], :, 1) : reshape([θ1old, ω1old, θ2old, ω2old], :, 1)
        ω1curr, ω2curr = [predict_y(gp, obs)[1][1] for gp in gps]
        θ1old, ω1old, θ2old, ω2old = θ1curr, ω1curr, θ2curr, ω2curr
        θ1curr = θ1curr + ω1curr*mechanism.Δt
        θ2curr = θ2curr + ω2curr*mechanism.Δt
    end
    q1, q2 = UnitQuaternion(RotX(θ1curr)), UnitQuaternion(RotX(θ1curr + θ2curr))
    x1 = [0, 0.5l1*sin(θ1curr), -0.5l1*cos(θ1curr)]
    x2 = [0, l1*sin(θ1curr) + 0.5l2*sin(θ1curr+θ2curr), -l1*cos(θ1curr) - 0.5l2*cos(θ1curr + θ2curr)]
    return CState([x1..., q2vec(q1)..., zeros(6)..., x2..., q2vec(q2)..., zeros(6)...])
end

function _predictdynamicscpmin(mechanism, gps, startobservation, steps, usesin)
    l = mechanism.bodies[2].shape.rh[2]
    xold, vold, θold, ωold = startobservation
    xcurr, θcurr = xold + vold*mechanism.Δt, θold + ωold*mechanism.Δt
    for _ in 1:steps
        obs = usesin ? reshape([xold, vold, sin(θold), cos(θold), ωold], :, 1) : reshape([xold, vold, θold, ωold], :, 1)
        vcurr, ωcurr = [predict_y(gp, obs)[1][1] for gp in gps]
        xold, vold, θold, ωold = xcurr, vcurr, θcurr, ωcurr
        xcurr = xcurr + vcurr*mechanism.Δt
        θcurr = θcurr + ωcurr*mechanism.Δt
    end
    q = UnitQuaternion(RotX(θcurr))
    return CState([0, xcurr, 0, 1, zeros(10)..., 0.5l*sin(θcurr)+xcurr, -0.5l*cos(θcurr), q2vec(q)..., zeros(6)...])
end

function _predictdynamicsfbmin(mechanism, gps, startobservation, steps, usesin)
    l = mechanism.bodies[1].shape.xyz[3]
    θ1old, ω1old, θ2old, ω2old = startobservation
    θ1curr, θ2curr = θ1old + ω1old*mechanism.Δt, θ2old + ω2old*mechanism.Δt
    for _ in 1:steps
        obs = usesin ? reshape([sin(θ1old), cos(θ1old), ω1old, sin(θ2old), cos(θ2old), ω2old], :, 1) : reshape([θ1old, ω1old, θ2old, ω2old], :, 1)
        ω1curr, ω2curr = [predict_y(gp, obs)[1][1] for gp in gps]
        θ1old, ω1old, θ2old, ω2old = θ1curr, ω1curr, θ2curr, ω2curr
        θ1curr = θ1curr + ω1curr*mechanism.Δt
        θ2curr = θ2curr + ω2curr*mechanism.Δt
    end
    x1 = [0, .5sin(θ1curr)l, -.5cos(θ1curr)l]
    x2 = [0, sin(θ1curr)l + .5sin(θ2curr)l, -cos(θ1curr)l - .5cos(θ2curr)l]
    x3 = [0, .5sin(θ2curr)l, -.5cos(θ2curr)l]
    x4 = [0, sin(θ2curr)l + 0.5sin(θ1curr)l, -cos(θ2curr)l - .5cos(θ1curr)l]
    q1 = UnitQuaternion(RotX(θ1curr))
    q2 = UnitQuaternion(RotX(θ2curr))
    return CState([x1..., q2vec(q1)..., zeros(6)..., x2..., q2vec(q2)..., zeros(6)..., x3..., q2vec(q2)..., zeros(6)..., x4..., q2vec(q1)..., zeros(6)...])
end


text:
# 24 bytes -> sizeof dpiDataBuffer on 64bit arch
primitive type OraDataBuffer 24 * 8 end

struct OraData
    is_null::Int32 # Specifies if the value refers to a null value (1) or not (0).
    value::OraDataBuffer
end

struct OraTimestamp
    year::Int16
    month::UInt8
    day::UInt8
    hour::UInt8
    minute::UInt8
    second::UInt8
    fsecond::UInt32 # Specifies the fractional seconds for the timestamp, in nanoseconds.
    tzHourOffset::Int8 # Specifies the hours offset from UTC. This value is only used for timestamp with time zone and timestamp with local time zone columns.
    tzMinuteOffset::Int8 # Specifies the minutes offset from UTC. This value is only used for timestamp with time zone and timestamp with local time zone columns. Note that this value will be negative in the western hemisphere. For example, when the timezone is -03:30, tzHourOffset will be -3 and tzMinuteOffset will be -30.

    function OraTimestamp(year::Integer, month::Integer, day::Integer,
                          hour::Integer, minute::Integer, second::Integer, fsecond::Integer,
                          tz_hour_offset::Integer, tz_minute_offset::Integer)

        return new(Int16(year), UInt8(month), UInt8(day),
                   UInt8(hour), UInt8(minute), UInt8(second), UInt32(fsecond),
                   Int8(tz_hour_offset), Int8(tz_minute_offset)
               )
    end
end

struct OraErrorInfo <: Exception
    code::Int32 # The OCI error code if an OCI error has taken place. If no OCI error has taken place the value is 0.
    offset::UInt16 # The parse error offset (in bytes) when executing a statement or the row offset when fetching batch error information. If neither of these cases are true, the value is 0.
    message::Ptr{UInt8} # The error message as a byte string in the encoding specified by the OraErrorInfo.encoding member.
    message_length::UInt32 # The length of the OraErrorInfo.message member, in bytes.
    encoding::Cstring # The encoding in which the error message is encoded as a null-terminated string. For OCI errors this is the CHAR encoding used when the connection was created. For ODPI-C specific errors this is UTF-8.
    fn_name::Cstring # The public ODPI-C function name which was called in which the error took place. This is a null-terminated ASCII string.
    action::Cstring # The internal action that was being performed when the error took place. This is a null-terminated ASCII string.
    sql_state::Cstring # The SQLSTATE code associated with the error. This is a 5 character null-terminated string.
    is_recoverable::Int32 # A boolean value indicating if the error is recoverable. This member always has a value of 0 unless both client and server are at release 12.1 or higher.
end

Base.showerror(io::IO, err::OraErrorInfo) = print(io, unsafe_string(err.message, err.message_length))

struct OraVersionInfo
    version::Int32 # Specifies the major version of the Oracle Client or Database.
    release::Int32 # Specifies the release version of the Oracle Client or Database.
    update::Int32 # Specifies the update version of the Oracle Client or Database.
    port_release::Int32 # Specifies the port specific release version of the Oracle Client or Database.
    port_update::Int32 # Specifies the port specific update version of the Oracle Client or Database.
    full_version::Int32 # Specifies the full version (all five components) as a number that is suitable for comparison with the result of the macro DPI_ORACLE_VERSION_TO_NUMBER.
end

mutable struct OraCommonCreateParams
    create_mode::OraCreateMode # Specifies the mode used for creating connections. It is expected to be one or more of the values from the enumeration OraCreateMode, OR’ed together. The default value is DPI_MODE_CREATE_DEFAULT.
    encoding::Cstring # Specifies the encoding to use for CHAR data, as a null-terminated ASCII string. Either an IANA or Oracle specific character set name is expected. NULL is also acceptable which implies the use of the NLS_LANG environment variable (or ASCII, if the NLS_LANG environment variable is not set). The default value is NULL.
    nencoding::Cstring # Specifies the encoding to use for NCHAR data, as a null-terminated ASCII string. Either an IANA or Oracle specific character set name is expected. NULL is also acceptable which implies the use of the NLS_NCHAR environment variable (or the same value as the OraCommonCreateParams.encoding member if the NLS_NCHAR environment variable is not set). The default value is NULL.
    edition::Ptr{UInt8} # Specifies the edition to be used when creating a standalone connection. It is expected to be NULL (meaning that no edition is set) or a byte string in the encoding specified by the OraCommonCreateParams.encoding member. The default value is NULL.
    edition_length::UInt32 # Specifies the length of the OraCommonCreateParams.edition member, in bytes. The default value is 0.
    driver_name::Ptr{UInt8} # Specifies the name of the driver that is being used. It is expected to be NULL or a byte string in the encoding specified by the OraCommonCreateParams.encoding member. The default value is NULL.
    driver_name_length::UInt32 # Specifies the length of the OraCommonCreateParams.driverName member, in bytes. The default value is 0.
end

struct OraAppContext
    namespace_name::Ptr{UInt8} # Specifies the value of the “namespace” parameter to sys_context(). It is expected to be a byte string in the encoding specified in the OraConnCreateParams structure and must not be NULL.
    namespace_name_length::UInt32 # Specifies the length of the OraAppContext.namespaceName member, in bytes.
    name::Ptr{UInt8} # Specifies the value of the “parameter” parameter to sys_context(). It is expected to be a byte string in the encoding specified in the OraConnCreateParams structure and must not be NULL.
    name_length::UInt32 # Specifies the length of the OraAppContext.name member, in bytes.
    value::Ptr{UInt8} # Specifies the value that will be returned from sys_context(). It is expected to be a byte string in the encoding specified in the OraConnCreateParams structure and must not be NULL.
    value_length::UInt32 # Specifies the length of the OraAppContext.value member, in bytes.
end

mutable struct OraConnCreateParams
    auth_mode::OraAuthMode
    connection_class::Ptr{UInt8}
    connection_class_length::UInt32
    purity::OraPurity
    new_password::Ptr{UInt8}
    new_password_length::UInt32
    app_context::Ptr{OraAppContext} # Specifies the application context that will be set when the connection is created. This value is only used when creating standalone connections. It is expected to be NULL or an array of OraAppContext structures. The context specified here can be used in logon triggers, for example. The default value is NULL.
    num_app_context::UInt32 # Specifies the number of elements found in the OraConnCreateParams.appContext member. The default value is 0.
    external_auth::Int32 # Specifies whether external authentication should be used to create the connection. If this value is 0, the user name and password values must be specified in the call to OraConn_create(); otherwise, the user name and password values must be zero length or NULL. The default value is 0.
    external_handle::Ptr{Cvoid} # Specifies an OCI service context handle created externally that will be used instead of creating a connection. The default value is NULL.
    pool_handle::Ptr{Cvoid} # Specifies the session pool from which to acquire a connection or NULL if a standalone connection should be created. The default value is NULL.
    tag::Ptr{UInt8}
    tag_length::UInt32
    match_any_tag::Int32
    out_tag::Ptr{UInt8}
    out_tag_length::UInt32
    out_tag_found::Int32
    sharding_key_columns::Ptr{Cvoid} # TODO
    num_sharding_key_columns::UInt8 # TODO
    super_sharding_key_columns::Ptr{Cvoid} # TODO
    num_super_sharding_key_columns::UInt8 # TODO
end

mutable struct OraPoolCreateParams
    min_sessions::UInt32 # Specifies the minimum number of sessions to be created by the session pool. This value is ignored if the OraPoolCreateParams.homogeneous member has a value of 0. The default value is 1.
    max_sessions::UInt32 # Specifies the maximum number of sessions that can be created by the session pool. Values of 1 and higher are acceptable. The default value is 1.
    session_increment::UInt32 # Specifies the number of sessions that will be created by the session pool when more sessions are required and the number of sessions is less than the maximum allowed. This value is ignored if the OraPoolCreateParams.homogeneous member has a value of 0. This value added to the OraPoolCreateParams.minSessions member value must not exceed the OraPoolCreateParams.maxSessions member value. The default value is 0.
    ping_interval::Int32 # Specifies the number of seconds since a connection has last been used before a ping will be performed to verify that the connection is still valid. A negative value disables this check. The default value is 60.
    ping_timeout::Int32 # Specifies the number of milliseconds to wait when performing a ping to verify the connection is still valid before the connection is considered invalid and is dropped. The default value is 5000 (5 seconds). This value is ignored in clients 12.2 and later since a much faster internal check is done by the Oracle client.
    homogeneous::Int32 # Specifies whether the pool is homogeneous or not. In a homogeneous pool all connections use the same credentials whereas in a heterogeneous pool other credentials are permitted. The default value is 1.
    external_auth::Int32 # Specifies whether external authentication should be used to create the sessions in the pool. If this value is 0, the user name and password values must be specified in the call to dpiPool_create(); otherwise, the user name and password values must be zero length or NULL. The default value is 0. External authentication cannot be used with homogeneous pools.
    get_mode::OraPoolGetMode # Specifies the mode to use when sessions are acquired from the pool. It is expected to be one of the values from the enumeration OraPoolGetMode. The default value is DPI_MODE_POOL_GET_NOWAIT. This value can be set after the pool has been created using the function dpiPool_setGetMode() and acquired using the function dpiPool_getGetMode().
    out_pool_name::Ptr{UInt8} # This member is populated upon successful creation of a pool using the function dpiPool_create(). It is a byte string in the encoding used for CHAR data. Any value specified prior to creating the session pool is ignored.
    out_pool_name_length::UInt32 # This member is populated upon successful creation of a pool using the function dpiPool_create(). It is the length of the OraPoolCreateParams.outPoolName member, in bytes. Any value specified prior to creating the session pool is ignored.
    timeout::UInt32 # Specifies the length of time (in seconds) after which idle sessions in the pool are terminated. Note that termination only occurs when the pool is accessed. The default value is 0 which means that no idle sessions are terminated. This value can be set after the pool has been created using the function dpiPool_setTimeout() and acquired using the function dpiPool_getTimeout().
    wait_timeout::UInt32 # Specifies the length of time (in milliseconds) that the caller should wait for a session to become available in the pool before returning with an error. This value is only used if the OraPoolCreateParams.getMode member is set to the value DPI_MODE_POOL_GET_TIMEDWAIT. The default value is 0. This value can be set after the pool has been created using the function dpiPool_setWaitTimeout() and acquired using the function dpiPool_getWaitTimeout().
    max_lifetime_session::UInt32 # Specifies the maximum length of time (in seconds) a pooled session may exist. Sessions in use will not be closed. They become candidates for termination only when they are released back to the pool and have existed for longer than maxLifetimeSession seconds. Session termination only occurs when the pool is accessed. The default value is 0 which means that there is no maximum length of time that a pooled session may exist. This value can be set after the pool has been created using the function dpiPool_setMaxLifetimeSession() and acquired using the function dpiPool_getMaxLifetimeSession().
    plsql_fixup_callback::Ptr{UInt8}
    plsql_fixup_callback_length::UInt32
end

struct OraDataTypeInfo
    oracle_type_num::OraOracleTypeNum # Specifies the type of the data. It will be one of the values from the enumeration OraOracleTypeNum, or 0 if the type is not supported by ODPI-C.
    default_native_type_num::OraNativeTypeNum # Specifies the default native type for the data. It will be one of the values from the enumeration OraNativeTypeNum, or 0 if the type is not supported by ODPI-C.
    oci_type_code::UInt16 # Specifies the OCI type code for the data, which can be useful if the type is not supported by ODPI-C.
    db_size_in_bytes::UInt32 # Specifies the size in bytes (from the database’s perspective) of the data. This value is only populated for strings and binary data. For all other data the value is zero.
    client_size_in_bytes::UInt32 # Specifies the size in bytes (from the client’s perspective) of the data. This value is only populated for strings and binary data. For all other data the value is zero.
    size_in_chars::UInt32 # Specifies the size in characters of the data. This value is only populated for string data. For all other data the value is zero.
    precision::Int16 # Specifies the precision of the data. This value is only populated for numeric and interval data. For all other data the value is zero.
    scale::Int8 # Specifies the scale of the data. This value is only populated for numeric data. For all other data the value is zero.
    fs_precision::Int16 # Specifies the fractional seconds precision of the data. This value is only populated for timestamp and interval day to second data. For all other data the value is zero.
    object_type_handle::Ptr{Cvoid} # Specifies a reference to the type of the object. This value is only populated for named type data. For all other data the value is NULL. This reference is owned by the object attribute, object type or statement and a call to OraObjectType_addRef() must be made if the reference is going to be used beyond the lifetime of the owning object.
end

struct OraQueryInfo
    name::Ptr{UInt8} # Specifies the name of the column which is being queried, as a byte string in the encoding used for CHAR data.
    name_length::UInt32 # Specifies the length of the OraQueryInfo.name member, in bytes.
    type_info::OraDataTypeInfo # Specifies the type of data of the column that is being queried. It is a structure of type OraDataTypeInfo.
    null_ok::Int32 # Specifies if the data that is being queried may return null values (1) or not (0).
end

struct OraStmtInfo
    is_query::Int32
    is_PLSQL::Int32
    is_DDL::Int32
    is_DML::Int32
    statement_type::OraStatementType
    is_returning::Int32
end

"High-level version for OraStmtInfo using Bool Julia type."
struct StmtInfo
    is_query::Bool
    is_PLSQL::Bool
    is_DDL::Bool
    is_DML::Bool
    statement_type::OraStatementType
    is_returning::Bool

    function StmtInfo(ora_stmt_info::OraStmtInfo)

        return new(
                Bool(ora_stmt_info.is_query),
                Bool(ora_stmt_info.is_PLSQL),
                Bool(ora_stmt_info.is_DDL),
                Bool(ora_stmt_info.is_DML),
                ora_stmt_info.statement_type,
                Bool(ora_stmt_info.is_returning)
            )
    end
end

struct OraBytes
    ptr::Ptr{UInt8}
    length::UInt32
    encoding::Cstring
end

Base.show(io::IO, ptr::Ptr{OraBytes}) = show(io, unsafe_load(ptr))

function Base.show(io::IO, ora_str::OraBytes)
    str = unsafe_string(ora_str.ptr, ora_str.length)
    enc = unsafe_string(ora_str.encoding)
    print(io, "OraBytes(", str, ", ", enc, ")")
end

#=
This structure is used for transferring encoding information from ODPI-C.
All of the information here remains valid as long as a reference is held
to the standalone connection or session pool from which the information was taken.
=#
struct OraEncodingInfo
    encoding::Cstring # The encoding used for CHAR data, as a null-terminated ASCII string.
    max_bytes_per_character::Int32 # The maximum number of bytes required for each character in the encoding used for CHAR data. This value is used when calculating the size of buffers required when lengths in characters are provided.
    nencoding::Cstring # The encoding used for NCHAR data, as a null-terminated ASCII string.
    nmax_bytes_per_character::Int32 # The maximum number of bytes required for each character in the encoding used for NCHAR data. Since this information is not directly available from Oracle it is only accurate if the encodings used for CHAR and NCHAR data are identical or one of ASCII or UTF-8; otherwise a value of 4 is assumed. This value is used when calculating the size of buffers required when lengths in characters are provided.
end

"Mirrors ODPI-C's OraEncodingInfo struct, but using Julia types."
struct EncodingInfo
    encoding::String
    max_bytes_per_character::Int
    nencoding::String
    nmax_bytes_per_character::Int

    function EncodingInfo(encoding_info_ref::Ref{OraEncodingInfo})
        encoding_info = encoding_info_ref[]
        return new(
            unsafe_string(encoding_info.encoding),
            Int(encoding_info.max_bytes_per_character),
            unsafe_string(encoding_info.nencoding),
            Int(encoding_info.nmax_bytes_per_character)
        )
    end
end

mutable struct Context
    handle::Ptr{Cvoid}

    function Context(handle::Ptr{Cvoid})
        new_context = new(handle)
        @compat finalizer(destroy!, new_context)
        return new_context
    end
end

function destroy!(ctx::Context)
    if ctx.handle != C_NULL
        result = dpiContext_destroy(ctx.handle)
        error_check(ctx, result)
        ctx.handle = C_NULL
    end
    nothing
end

mutable struct Pool
    context::Context
    handle::Ptr{Cvoid}
    name::String

    function Pool(context::Context, handle::Ptr{Cvoid}, name::String)
        new_pool = new(context, handle, name)
        @compat finalizer(destroy!, new_pool)
        return new_pool
    end
end

function destroy!(pool::Pool)
    if pool.handle != C_NULL
        result = dpiPool_release(pool.handle)
        error_check(context(pool), result)
        pool.handle = C_NULL
    end
    nothing
end

mutable struct Connection
    context::Context
    handle::Ptr{Cvoid}
    encoding_info::EncodingInfo
    pool::Union{Nothing, Pool}

    function Connection(context::Context, handle::Ptr{Cvoid}, pool::Union{Nothing, Pool})

        # this driver currently only supports UTF-8 encoding
        function check_supported_encoding(ei::EncodingInfo)
            @assert ei.encoding ∈ SUPPORTED_CONNECTION_ENCODINGS "Unsupported encoding for CHARS: $(ei.encoding). Currently, Oracle.jl supports only $(SUPPORTED_CONNECTION_ENCODINGS)."
            @assert ei.nencoding ∈ SUPPORTED_CONNECTION_ENCODINGS  "Unsupported encoding for CHARS: $(ei.nencoding). Currently, Oracle.jl supports only $(SUPPORTED_CONNECTION_ENCODINGS)."
        end

        ei = EncodingInfo(context, handle)
        check_supported_encoding(ei)

        new_connection = new(context, handle, ei, pool)
        @compat finalizer(destroy!, new_connection)
        return new_connection
    end
end

function EncodingInfo(context::Context, connection_handle::Ptr{Cvoid})
    encoding_info_ref = Ref{OraEncodingInfo}()
    result = dpiConn_getEncodingInfo(connection_handle, encoding_info_ref)
    error_check(context, result)
    return EncodingInfo(encoding_info_ref)
end

function destroy!(conn::Connection)
    if conn.handle != C_NULL
        result = dpiConn_release(conn.handle)
        error_check(context(conn), result)
        conn.handle = C_NULL
        conn.pool = nothing
    end
    nothing
end

mutable struct Stmt{statement_type}
    connection::Connection
    handle::Ptr{Cvoid}
    scrollable::Bool
    info::StmtInfo
    bind_count::UInt32
    bind_names::Vector{String}
    bind_names_index::Dict{String, UInt32} # maps bind_name to bind position
    is_open::Bool
    columns_info::Union{Nothing, Vector{OraQueryInfo}}
end

const QueryStmt = Stmt{ORA_STMT_TYPE_SELECT}

function destroy!(stmt::Stmt)
    if stmt.handle != C_NULL
        result = dpiStmt_release(stmt.handle)
        error_check(context(stmt), result)
        stmt.handle = C_NULL
    end
    nothing
end

"Holds a 1-indexed vector of OraData."
abstract type AbstractOracleValue{O,N} end

"Wraps a OraData handle managed by extern ODPI-C. 1-indexed."
struct ExternOracleValue{O,N,P} <: AbstractOracleValue{O,N}
    parent::P
    data_handle::Ptr{OraData}
    use_add_ref::Bool # triggers a call to add_ref, if data_handle is a LOB, statement, object or rowid that is owned by the statement
end

function ExternOracleValue(parent::P, oracle_type::OraOracleTypeNum, native_type::OraNativeTypeNum, handle::Ptr{OraData}; use_add_ref::Bool=false) where {P}
    return ExternOracleValue{oracle_type, native_type, P}(parent, handle, use_add_ref)
end

"Wraps a OraData handle managed by Julia. 1-indexed."
struct JuliaOracleValue{O,N,T} <: AbstractOracleValue{O,N}
    buffer::Vector{T}
end

function JuliaOracleValue(oracle_type::OraOracleTypeNum, native_type::OraNativeTypeNum, ::Type{T}, capacity::Integer=1) where {T}
    buffer = undef_vector(T, capacity)
    return JuliaOracleValue{oracle_type,native_type,T}(buffer)
end

struct FetchResult
    found::Bool
    buffer_row_index::UInt32
end

Base.show(io::IO, result::FetchResult) = print(io, "FetchResult(", result.found, ", ", Int(result.buffer_row_index), ")")

struct FetchRowsResult
    buffer_row_index::UInt32
    num_rows_fetched::UInt32
    more_rows::Int32
end

Base.show(io::IO, result::FetchRowsResult) = print(io, "FetchRowsResult(", Int(result.buffer_row_index), ", " ,Int(result.num_rows_fetched), ", ",Int(result.more_rows), ")")

struct CursorSchema
    column_query_info::Vector{OraQueryInfo}
    column_names_index::Dict{String, Int}
end

struct Cursor
    stmt::QueryStmt
    schema::CursorSchema
end

struct ResultSetRow
    schema::CursorSchema
    data::Vector{Any}
end

struct ResultSet
    schema::CursorSchema
    rows::Vector{ResultSetRow}
end

"Safe version of OraCommonCreateParams"
mutable struct CommonCreateParams
    create_mode::Union{Nothing, OraCreateMode}
    encoding::String
    nencoding::String
    edition::Union{Nothing, String}
    driver_name::Union{Nothing, String}
end

"Safe version of OraConnCreateParams"
mutable struct ConnCreateParams
    auth_mode::OraAuthMode
    pool::Union{Nothing, Pool}
end

mutable struct Variable{T}
    connection::Connection
    handle::Ptr{Cvoid}
    oracle_type::OraOracleTypeNum
    native_type::OraNativeTypeNum
    max_byte_string_size_in_bytes::UInt32 # size. the maximum size of the buffer used for transferring data to/from Oracle. This value is only used for variables transferred as byte strings.
    is_PLSQL_array::Int32 # boolean value indicating if the variable refers to a PL/SQL array or simply to buffers used for binding or fetching data.
    obj_type_handle::Ptr{Cvoid}
    buffer_handle::Ptr{OraData} # a pointer to an array of dpiData structures that are used to transfer data to/from the variable. These are allocated when the variable is created and the number of structures corresponds to the maxArraySize.
    buffer_capacity::UInt32 # maxArraySize. the maximum number of rows that can be fetched or bound at one time from the database, or the maximum number of elements that can be stored in a PL/SQL array.

    function Variable(
            connection::Connection,
            ::Type{T},
            handle::Ptr{Cvoid},
            oracle_type::OraOracleTypeNum,
            native_type::OraNativeTypeNum,
            max_byte_string_size_in_bytes::UInt32,
            is_PLSQL_array::Int32,
            obj_type_handle::Ptr{Cvoid},
            buffer_handle::Ptr{OraData},
            buffer_capacity::UInt32) where {T}

        new_ora_variable = new{T}(
            connection,
            handle,
            oracle_type,
            native_type,
            max_byte_string_size_in_bytes,
            is_PLSQL_array,
            obj_type_handle,
            buffer_handle,
            buffer_capacity)

        @compat finalizer(destroy!, new_ora_variable)

        return new_ora_variable
    end
end

function destroy!(v::Variable)
    if v.handle != C_NULL
        result = dpiVar_release(v.handle)
        error_check(context(v), result)
        v.handle = C_NULL
        v.obj_type_handle = C_NULL
        v.buffer_handle = C_NULL
    end
    nothing
end

mutable struct Lob{ORATYPE,T}
    parent::T
    handle::Ptr{Cvoid}
    is_open::Bool

    function Lob(p::T, handle::Ptr{Cvoid}, oracle_type::OraOracleTypeNum; use_add_ref::Bool=false) where {T}
        check_valid_lob_oracle_type_num(oracle_type)
        new_lob = new{oracle_type, T}(p, handle, true)
        @compat finalizer(destroy!, new_lob)

        if use_add_ref
            add_ref(new_lob)
        end

        return new_lob
    end
end

function destroy!(lob::Lob)
    if lob.handle != C_NULL
        result = dpiLob_release(lob.handle)
        error_check(context(lob), result)
        lob.handle = C_NULL
    end
    nothing
end


text:
#Ensemble Kalman Inversion: specific structures and function definitions

"""
    Inversion <: Process

An ensemble Kalman Inversion process
"""
struct Inversion <: Process end

"""
    find_ekp_stepsize(ekp::EnsembleKalmanProcess{FT, IT, Inversion}, g::Array{FT, 2}; cov_threshold::FT=0.01) where {FT}

Find largest stepsize for the EK solver that leads to a reduction of the determinant of the sample
covariance matrix no greater than cov_threshold.
"""
function find_ekp_stepsize(
    ekp::EnsembleKalmanProcess{FT, IT, Inversion},
    g::Array{FT, 2};
    cov_threshold::FT = 0.01,
) where {FT, IT}
    accept_stepsize = false
    if !isempty(ekp.Δt)
        Δt = deepcopy(ekp.Δt[end])
    else
        Δt = FT(1)
    end
    # final_params [N_par × N_ens]
    cov_init = cov(get_u_final(ekp), dims = 2)
    while accept_stepsize == false
        ekp_copy = deepcopy(ekp)
        update_ensemble!(ekp_copy, g, Δt_new = Δt)
        cov_new = cov(get_u_final(ekp_copy), dims = 2)
        if det(cov_new) > cov_threshold * det(cov_init)
            accept_stepsize = true
        else
            Δt = Δt / 2
        end
    end

    return Δt

end

"""
    update_ensemble!(ekp::EnsembleKalmanProcess{FT, IT, <:Inversion}, g::Array{FT,2} cov_threshold::FT=0.01, Δt_new=nothing) where {FT, IT}

Updates the ensemble according to which type of Process we have. Model outputs `g` need to be a `N_obs × N_ens` array (i.e data are columms).
"""
function update_ensemble!(
    ekp::EnsembleKalmanProcess{FT, IT, Inversion},
    g::Array{FT, 2};
    cov_threshold::FT = 0.01,
    Δt_new = nothing,
    deterministic_forward_map = true,
) where {FT, IT}

    # Update follows eqns. (4) and (5) of Schillings and Stuart (2017)

    #catch works when g non-square
    if !(size(g)[2] == ekp.N_ens)
        throw(DimensionMismatch("ensemble size in EnsembleKalmanProcess and g do not match, try transposing g or check ensemble size"))
    end

    # u: N_par × N_ens
    # g: N_obs × N_ens
    u = get_u_final(ekp)
    N_obs = size(g, 1)

    cov_init = cov(u, dims = 2)
    cov_ug = cov(u, g, dims = 2, corrected = false) # [N_par × N_obs]
    cov_gg = cov(g, g, dims = 2, corrected = false) # [N_obs × N_obs]

    if !isnothing(Δt_new)
        push!(ekp.Δt, Δt_new)
    elseif isnothing(Δt_new) && isempty(ekp.Δt)
        push!(ekp.Δt, FT(1))
    else
        push!(ekp.Δt, ekp.Δt[end])
    end

    # Scale noise using Δt
    scaled_obs_noise_cov = ekp.obs_noise_cov / ekp.Δt[end]
    noise = rand(MvNormal(zeros(N_obs), scaled_obs_noise_cov), ekp.N_ens)

    # Add obs_mean (N_obs) to each column of noise (N_obs × N_ens) if
    # G is deterministic
    y = deterministic_forward_map ? (ekp.obs_mean .+ noise) : (ekp.obs_mean .+ zero(noise))

    # N_obs × N_obs \ [N_obs × N_ens]
    # --> tmp is [N_obs × N_ens]
    tmp = (cov_gg + scaled_obs_noise_cov) \ (y - g)
    u += (cov_ug * tmp) # [N_par × N_ens]

    # store new parameters (and model outputs)
    push!(ekp.u, DataContainer(u, data_are_columns = true))
    push!(ekp.g, DataContainer(g, data_are_columns = true))

    # Store error
    compute_error!(ekp)

    # Check convergence
    cov_new = cov(get_u_final(ekp), dims = 2)
    cov_ratio = det(cov_new) / det(cov_init)
    if cov_ratio < cov_threshold
        @warn string(
            "New ensemble covariance determinant is less than ",
            cov_threshold,
            " times its previous value.",
            "\nConsider reducing the EK time step.",
        )
    end
end


text:
using CGPrism3D
using Test
using LinearAlgebra

J,K,alpha = 2,2,0.5
dataT = dataTet(J,K,alpha)
@show unique(dataT[2])
#dependence of tensorGlueTet3D on J,K,alpha is for dimension factors
#tet2 = tensorGlueTet3DN(dataT,dataT,[1,2,3],[1,2,3],dataT)
tb2 = tensorBlock(dataT,[],[2,3,4,5,6],[1],[0.,0.5,0.5])[1]
#U,s,V =svd(tb2)
#@show s
#@time prA = tensorGlueTet3DN(tet2,dataT,[2,7,9],[1,2,3],dataT)
#@show length(prA[1])

#@time prA2 = tensorGlueTet3DN2(tet2,dataT,[2,7,9],[1,2,3],dataT)
#@show length(prA[1]), length(prA2[1])
#@show sum(prA[2]-prA2[2])
function testA(J::Int64,K::Int64,alpha::Float64,loop)
    dataT = dataTet(J,K,alpha)
    #dependence of tensorGlueTet3D on J,K,alpha is for dimension factors
    tet2 = tensorGlueTet3DN(dataT,dataT,[1,2,3],[1,2,3],dataT)
    prA = tensorGlueTet3DN(tet2,dataT,[2,7,9],[1,2,3],dataT)

    #Using only Prism A to coarse grain
    for l in 1:loop
        println("loop=$l")
        #dataT,dataV = fullSplitTet3DN(prA,[4,5,6],[7,8,9,10,11,12],[1,2,3],dataT,J,K,alpha)
        #dataT = tensorPermuteN(dataT,[4,5,6,1,2,3])
        #prA =
        pruvA = tensorGlueTet3DN(prA,prA,[1,5,6,8,9],[2,4,6,12,11],dataT)
        pruvAA = tensorSumTetO(pruvA,[1],dataT)
        #@time pruvAAn = tensorSumTet3D(pruvA,[1],dataT)
        #choose Fsymb for higher level to do 2-2 move
        #println("Starting 2-2 move on face 1...")
        #fA0F = tensor22moveT(pruvAA,[18,15,6,2,7],dataT)
        fA0F = tensor22moveF(pruvAA,[18,15,6,2,7],K,K,alpha)
        #fA0F = tensor22moveT2(pruvAA,[18,15,6,2,7],dataT,[7,1,17],K)
        #println("2-2 move on face 1 done")
        #println("Starting 2-2 move on face 2...")
        #fA01F = tensor22moveT(fA0F,[17,18,10,9,8],dataT)
        fA01F = tensor22moveF(fA0F,[17,18,10,9,8],K,K,alpha)
        #fA01F = tensor22moveT2(fA0F,[17,18,10,9,8],dataT,[8,7,11],K)
        #println("2-2 move on face 2 done")
        #println("Starting 2-2 move on face 3...")
        #fA012F = tensor22moveT(fA01F,[14,13,3,2,4],dataT)
        fA012F = tensor22moveF(fA01F,[14,13,3,2,4],K,K,alpha)
        #fA012F = tensor22moveT2(fA01F,[14,13,3,2,4],dataT,[4,7,16],K)
        #println("2-2 move on face 3 done")
        #println("Starting 3-1 move on tet1...")
        f133A = fullSplitTet3DN(fA012F,[9,6,18],[1,2,3,4,5,10,12,13,14,15,16,17],[7,8,11],dataT,J,K,alpha)
        #println("3-1 move on tet1 done")
        #println("Starting 3-1 move on tet2...")
        CGPrsmA1 = fullSplitTet3DN(f133A[2],[8,10,2],[1,3,5,6,7,9,12,14,15],[13,4,11],dataT,J,K,alpha)[2]
        tperm = tensorPermuteN(CGPrsmA1,[5,1,11,2,6,3,10,12,7,8,4,9])
        prA = tensorProjectPrism(tperm,K)
        dataT = tensorPermuteN(fullSplitTet3DN(prA,[7,8,9,10,11,12],[4,5,6],[1,2,3],dataT,J,K,alpha)[2],[4,5,6,1,2,3])
        @show length(prA[1])
        #@show println(dataT1[1])
        #@show println(dataT[1])
        #@show sort(dataT1[1]) == sort(dataT[1])
        tb = tensorBlock(prA,[4,5,6],[7,8,9,10,11,12],[1,2,3],[0.,0.5,0.5])[1]
        @show sum(prA[2])
        U,s,V = svd(tb)
        @show s
    end
    return prA,dataT
end

#for i in 1:3
@time ff = testA(J,K,alpha,30)
#fullSplitTet3DN(ff,)
#@time ff1 = testA(J,K,alpha)
#tb = tensorBlock(ff[1],[4,5,6],[7,8,9,10,11,12],[1,2,3],[0.,0.5,0.5])[1]
#@show sum(ff[2][2])
#U,s,V = svd(tb)
#@show s
#end

#@show length(ff[1])


text:
module LTMSim

using SimpleHypergraphs
using DataFrames
using Random
using StatsBase
using DataStructures
using LightGraphs
using CSV
using Random
using LaTeXStrings

export generateHypergraph
export simulate!
export randMetaV, randMetaE
export proportionalMetaV, proportionalMetaE

export randomH, randomHkuniform, randomHduniform, randomHpreferential
export dual

export bisect, greedy_tss, greedy_tss_2section
export sub_tss
export optimize_seed_set

include("tools.jl")

include("models/random_models.jl")
include("models/dual.jl")

include("heuristics/additive.jl")
include("heuristics/subtractive.jl")
include("heuristics/optimization.jl")

end # module


text:
using MakieGallery

# load the database.  TODO this is a global and should be changed.
database = MakieGallery.load_database();

# where is the refimage repo?
repo = joinpath(get(ENV, "MAKIEGALLERY_REFIMG_PATH", joinpath(homedir(), "ReferenceImages")), "gallery")

recordings = joinpath(@__DIR__, "test_recordings")

# MakieGallery.generate_preview(recordings, joinpath(recordings, "index.html"))

# generate `thumb.jpg` for every directory in `recordings`
MakieGallery.generate_thumbnails(recordings)

# generate HTML pages for the Gallery and
MakieGallery.gallery_from_recordings(
    recordings,
    joinpath(recordings, "index.html");
    hltheme = MakieGallery.Highlights.Themes.DefaultTheme
)

# move this content to the repo
cp(recordings, repo, force = true)


text:
export to_bcoords, from_bcoords, arcintervals, totallength
export boundarymap
#######################################################################################
## Arclengths
#######################################################################################
"""
    totallength(o::Obstacle)
Return the total boundary length of `o`.
"""
@inline totallength(o::Wall) = norm(o.ep - o.sp)
@inline totallength(o::Semicircle) = π*o.r
@inline totallength(o::Circular) = 2π*o.r
@inline totallength(e::Ellipse) = 4e.arc

@inline totallength(bd::Billiard) = sum(totallength(x) for x in bd.obstacles)

"""
    arcintervals(bd::Billiard) -> s
Generate a vector `s`, with entries being the delimiters of the
arclengths of the obstacles of the billiard.
The arclength from `s[i]` to `s[i+1]` is the arclength spanned by
the `i`th obstacle.

`s` is used to transform an arc-coordinate `ξ` from local to global and
vice-versa. A local `ξ` becomes global by adding `s[i]` (where `i` is the
index of current obstacle). A global `ξ` becomes local by subtracting `s[i]`.

See also [`boundarymap`](@ref), [`to_bcoords`](@ref), [`from_bcoords`](@ref).
"""
function arcintervals(bd::Billiard{T, D}) where {T, D}
    intervals = SVector{D+1,T}(0, map(x->totallength(x), bd.obstacles)...)
    return cumsum(intervals)
end

################################################################################
## Coordinate systems
################################################################################
"""
    to_bcoords(pos, vel, o::Obstacle) -> ξ, sφ
Convert the real coordinates `pos, vel` to
boundary coordinates (also known as Birkhoff coordinates)
`ξ, sφ`, assuming that `pos` is on the obstacle.

`ξ` is the arc-coordinate, i.e. it parameterizes the arclength of the
obstacle. `sφ` is the
sine of the angle between the velocity vector and the vector normal
to the obstacle.

The arc-coordinate `ξ` is measured as:
* the distance from start point to end point in `Wall`s
* the arc length measured counterclockwise from the open face in `Semicircle`s
* the arc length measured counterclockwise from the rightmost point
  in `Circular`/`Ellipse`s

Notice that this function returns the *local* arclength. To get the global
arclength parameterizing an entire billiard, simply do
`ξ += arcintervals(bd)[i]` if the index of obstacle `o` is `i`.

See also [`from_bcoords`](@ref), which is the inverse
function.
"""
to_bcoords(p::AbstractParticle, o::Obstacle) = to_bcoords(p.pos, p.vel, o)

function to_bcoords(pos::SV, vel::SV, o::Obstacle)
    n = normalvec(o, pos)
    sinφ = cross2D(vel, n)
    ξ = _ξ(pos, o)
    return ξ, sinφ
end

# Internal function ξ is simply returning the boundary "arc-coordinate"
_ξ(pos::SV, o::Wall) = norm(pos - o.sp)

function _ξ(pos::SV{T}, o::Semicircle{T}) where {T<:AbstractFloat}
    #project pos on open face
    chrd = SV{T}(-o.facedir[2],o.facedir[1])
    d = (pos - o.c)/o.r
    x = dot(d, chrd)
    r =  acos(clamp(x, -1, 1))*o.r
    return r
end

function _ξ(pos::SV{T}, o::Circular{T}) where {T<:AbstractFloat}
    d = (pos - o.c)/o.r
    if d[2] > 0
        r = acos(clamp(d[1], -1, 1))
    else
        r = π + acos(-clamp(d[1], -1, 1))
    end
    return r*o.r
end

function _ξ(pos::SV{T}, o::Ellipse{T}) where {T<:AbstractFloat}
    θ = atan(pos[2] - o.c[2], pos[1] - o.c[1]) + π
    return ellipse_arclength(θ, o)
end

"""
    from_bcoords(ξ, sφ, o::Obstacle) -> pos, vel
Convert the boundary coordinates `ξ, sφ` on the obstacle to
real coordinates `pos, vel`.

Note that `vel` always points away from the obstacle.

This function is the inverse of [`to_bcoords`](@ref).
"""
function from_bcoords(ξ::T, sφ::T, o::Obstacle{T}) where {T}
    pos = real_pos(ξ, o)
    cφ = sqrt(1-sφ^2) # = cos(asin(sφ))
    n = normalvec(o, pos)
    vel = SV{T}(n[1]*cφ + n[2]*sφ, -n[1]*sφ + n[2]*cφ)

    return pos, vel
end

"""
    real_pos(ξ, o::Obstacle)
Converts the arclength coordinate `ξ` relative to the obstacle `o` into a real space
position vector.
"""
real_pos(ξ, o::Wall) = o.sp + ξ*normalize(o.ep - o.sp)

function real_pos(ξ, o::Semicircle{T}) where T
    sξ, cξ = sincos(ξ/o.r)
    chrd = SV{T}(-o.facedir[2], o.facedir[1])
    return o.c - o.r*(sξ*o.facedir - cξ*chrd)
end

real_pos(ξ, o::Circular{T}) where T = o.c .+ o.r * SV{T}(cossin(ξ/o.r))

function real_pos(ξ, o::Ellipse{T}) where T
    error("`from_bcoords` is not implemented for Ellipse. If you know of a way "*
    "to find a point (2-vector) on an ellipse given an arclength ξ, please let us know!")
    return nothing
end

"""
    from_bcoords(ξ, sφ, bd::Billiard, intervals = arcintervals(bd))
Same as above, but now `ξ` is considered to be the global arclength,
parameterizing the entire billiard, instead of a single obstacle.
"""
function from_bcoords(ξ, sφ, bd::Billiard{T}, intervals = arcintervals(bd)) where T

    for (i, obst) ∈ enumerate(bd)
        if ξ <= intervals[i+1]
            pos, vel  = from_bcoords(ξ - intervals[i], sφ, obst)
            return pos, vel, i
        end
    end
    throw(DomainError(ξ ,"ξ is too large for this billiard!"))
end

#######################################################################################
## Boundary Map
#######################################################################################
"""
    boundarymap(p, bd, t [,intervals]) → bmap, arclengths
Compute the boundary map of the particle `p` in the billiard `bd` by evolving
the particle for total amount `t` (either float for
time or integer for collision number).

Return a vector of 2-vectors `bmap` and also `arclengths(bd)`. The first entry of each
element of `bmap` is the arc-coordinate at collisions ``\\xi``, while the second  is the
sine of incidence angle ``\\sin(\\phi_n)``.

The measurement direction of the arclengths of the individual obstacles
is dictated by their order in `bd`. The sine of the angle is computed
*after* specular reflection has taken place.

The returned values of this function can be used in conjuction with the
function [`plot_boundarymap`](@ref) (requires `using PyPlot`) to plot the boundary map
in an intuitive way.

*Notice* - this function only works for normal specular reflection. Random reflections
or ray-splitting will give unexpected results.

See also [`to_bcoords`](@ref), [`boundarymap_portion`](@ref).
See [`parallelize`](@ref) for a parallelized version.
"""
function boundarymap(par::AbstractParticle{T}, bd::Billiard{T}, t,
                     intervals = arcintervals(bd)) where {T}

    p = copy(par)
    bmap = SV{T}[]
    if typeof(t) == Int
        sizehint!(bmap, t)
    end
    count = zero(t); t_to_write = zero(T)

    while count < t
        i, tmin = bounce!(p,bd)
        t_to_write += tmin

        if isperiodic(i, bd)
            continue # do not write output if collision with with PeriodicWall
        else
            ξ, sφ = to_bcoords(p, bd[i])
            @inbounds push!(bmap, SV{T}(ξ + intervals[i], sφ))
            # set counter
            count += increment_counter(t, t_to_write)
            t_to_write = zero(T)
        end
    end #time, or collision number, loop
    return bmap, intervals
end


text:
function χ²(y, y_fit, σ)
	return sum(((y .- y_fit) ./ σ).^2)
end

function χ²_fit(x, y, σ)
	num = sum(y ./ σ.^2) * sum(x.^2 ./ σ.^2) - sum(x ./ σ.^2) * sum(y .* x ./ σ.^2)
	den = sum(1 ./ σ.^2) * sum(x.^2 ./ σ.^2) - sum(x ./ σ.^2)^2
	a = num / den

	num = sum(y .* x ./ σ.^2) * sum(1 ./ σ.^2) - sum(x ./ σ.^2) * sum(y ./ σ.^2)
	den = sum(1 ./ σ.^2) * sum(x.^2 ./ σ.^2) - sum(x ./ σ.^2)^2
	b = num / den

	y_fit = a .+ b .* x

	return Dict(
		"y_fit"  => y_fit,
		"a"      => a,
		"b"      => b,
		"χ²/dof" => χ²(y, y_fit, σ) / (length(x)-2)
	)
end

max_stars_repo_name: cameronperot/explosive-percolation

max_issues_repo_name: cameronperot/explosive-percolation

max_forks_repo_name: cameronperot/explosive-percolation


text:
using BinDeps

@BinDeps.setup

libproj_min_ver = v"4.9.0" # First verison with supporting the geodesic API

# Validate the version of the C library.  This is somewhat duplicated
# code, but has to happen here in the setup stage before any of the other code
# is loaded.
function validate_proj_version(libname, handle)
    pj_get_release = Libdl.dlsym_e(handle, :pj_get_release)
    pj_get_release != C_NULL || return false
    verstr = bytestring(ccall(pj_get_release, Cstring, ()))
    m = match(r"(\d+).(\d+).(\d+),.+", verstr)
    m !== nothing || return false
    ver = VersionNumber(parse(Int, m[1]), parse(Int, m[2]), parse(Int, m[3]))
    ver >= libproj_min_ver
end

libproj = library_dependency("libproj", validate=validate_proj_version)

# Provide a BuildProcess for linux, since the package manager provided libproj
# can be woefully old.
src_url = "http://download.osgeo.org/proj"
libproj_src_ver = v"4.9.2"
libproj_name = "proj-$(libproj_src_ver).tar.gz"
provides(Sources, URI("$src_url/$libproj_name"), libproj, os = :Unix)
prefix   = joinpath(BinDeps.depsdir(libproj), "usr")
srcdir   = joinpath(BinDeps.depsdir(libproj), "src", "proj-$(libproj_src_ver)")
builddir = joinpath(BinDeps.depsdir(libproj), "builds", "libproj-$(libproj_src_ver)_version")
libdir   = joinpath(prefix, "lib")
provides(BuildProcess,
    (@build_steps begin
        GetSources(libproj)
        CreateDirectory(builddir)
        @build_steps begin
            ChangeDirectory(builddir)
            FileRule(joinpath(libdir, "libproj.so"),
                     @build_steps begin
                         ChangeDirectory(srcdir)
                         `./configure --prefix=$(prefix)`
                         `make`
                         `make install`
                     end
            )
        end
    end),
    libproj,
    os = :Unix,
)

@osx_only begin
    using Homebrew
    provides(Homebrew.HB, "proj", libproj, os = :Darwin)
end

# Create deps/deps.jl
# TODO(chris.foster): BinDeps.@load_dependencies doesn't work as advertised, so
# we're using @install instead.  See
# https://github.com/JuliaLang/BinDeps.jl/issues/196
@BinDeps.install Dict(:libproj => :libproj)


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

module x2

greet() = print("Hello World!")

end # module

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
"""
	gmtbegin(name::String=""; fmt)

Start a GMT session in modern mode (GMT >= 6).
'name' contains the figure name with or without extension. If an extension is used
(e.g. "map.pdf") it is used to select the image format.

As an alternative use 'fmt' as a string or symbol containing the format (ps, pdf, png, PNG, tif, jpg, eps).

By default name="GMTplot" and fmt="ps"
"""
function gmtbegin(name::String=""; fmt=nothing, verbose=nothing)
	cmd = "begin"       # Default name (GMTplot.ps) is set in gmt_main()
	if (name != "")  cmd *= " " * get_format(name, fmt)  end
	if (verbose !== nothing)  cmd *= " -V" * string(verbose)  end
	gmt("destroy")		# Always start with a clean session
	gmt(cmd)
	return nothing
end

"""
	gmtend(show=false, verbose=nothing)

Ends a GMT session in modern mode (GMT >= 6) and optionaly shows the figure
"""
function gmtend(show=nothing; verbose=nothing)
	cmd = "end"
	if (show !== nothing)  cmd *= " show"  end
	if (verbose !== nothing)  cmd *= " -V" * string(verbose)  end
	gmt(cmd)
	gmt("destroy")		# Lieve it in a clean state
	return nothing
end

"""
	function gmtfig(name::String; fmt=nothing, opts="")

Set attributes for the current modern mode session figure.
'name' name of the new (or resumed) figure. It may contain an extension.
'fmt'  figures graphics format (or formats, e.g. fmt="eps,pdf"). Not needed if 'name' has extension
'opts' Sets one or more comma-separated options (and possibly arguments) that can be passed to psconvert when preparing this figure.
"""
function gmtfig(name::String; fmt=nothing, opts="")
	global IamModern
	if (!IamModern)  error("Not in modern mode. Must run 'gmtbegin' first")  end

	cmd = "figure"       # Default name (GMTplot.ps) is set in gmt_main()
	if (name == "")  error("figure name cannot be empty")   end
	cmd *= " " * get_format(name, fmt)
	if (opts != "")  cmd *= " " * opts  end
	gmt(cmd)
	return nothing
end

function inset(fim=nothing; stop=false, kwargs...)

	global IamModern
	if (!IamModern)  error("Not in modern mode. Must run 'gmtbegin' first")  end

	d = KW(kwargs)
	cmd = parse_common_opts(d, "", [:c :F :V_params], true)
	cmd = parse_these_opts(cmd, d, [[:M :margins], [:N :no_clip]])
	cmd = parse_type_anchor(d, cmd, [[:D :inset :inset_box]])

	do_show = false
	if (fim !== nothing)
		t = lowercase(string(fim))
		if (t == "end" || t == "stop" || t == "show")  stop = true  end
		if (t == "show")  do_show = true  end
	end

	if (!stop)
		if (dbg_print_cmd(d, cmd) !== nothing)  return cmd  end		# Vd=2 cause this return
		gmt("inset begin " * cmd);
	else
		gmt("inset end");
		if (do_show || haskey(d, :show))  gmt("end" * show);  end
	end
end

function get_format(name, fmt=nothing, d=nothing)
	# Get the fig name and format. If format not specified, default to FMT (ps)
	# NAME is supposed to always exist (otherwise, errors)
	fname, ext = splitext(string(name))
	if (ext != "")
		fname *= " " * ext[2:end]
	elseif (fmt !== nothing)
		fname *= " " * string(fmt)      # No checking
	elseif (d !== nothing)
		if (haskey(d, :fmt))  fname *= " " * string(d[:fmt])
		else                  fname *= " " * FMT		# Then use default format
		end
	else
		fname *= " " * FMT
	end
	return fname
end


text:
Base.@irrational Ωe 0.00007292115 big(0.00007292115)
Ωe, const Omegae = Ωe
gsw_f(latitude::T) where {T<:Real} = T(2) * Ωe * sind(latitude)


text:
## Code for RAYS #########################################################################################################################

"""
Ray( origin, dir, tmin, tmax, depth)

It creates a **Ray**. When not specified tmin = 1e-5, tmax = +∞ and depth = 0.f0
"""
struct Ray
    origin::Point
    dir::Vec
    tmin::Float32
    tmax::Float32
    depth::Int16

    Ray( origin, dir, tmin=1e-5, tmax=Inf, depth=0) = new(origin, dir, tmin, tmax, depth)

end

Base.:*(T::Transformation, R::Ray) = Ray(T*R.origin, T*R.dir, R.tmin, R.tmax, R.depth )

Base.:isapprox(ray1::Ray, ray2::Ray) = isapprox(ray1.origin, ray2.origin) && isapprox(ray1.dir, ray2.dir)

"""
    at(ray, t)

It calculates the position of the ray at the instant *t*.
"""
at(ray::Ray, t::Number) = ray.origin + ray.dir*t

## Code for SHAPES #########################################################################################

abstract type Shape
end
Base.:≈(S1::Shape,S2::Shape) = S1.transformation ≈ S2.transformation && S1.material ≈ S2.material

"""
    Sphere(T)

It creates a **Sphere**, where T is a generic ``Transformation`` applied to the unit sphere centered in the origin, M is the ``Material`` of the sphere.
"""
struct Sphere <: Shape
    transformation::Transformation
    material::Material

    Sphere(transformation::Transformation = Transformation(),
           material::Material = Material()) =
           new(transformation, material)
#    Sphere() = new(Tansformation(), Material())
end

## Hidden methods for sphere
function _sphere_point_to_uv(point::Point)
    u = Float32(atan(point.y, point.x) / (2.0f0 * pi))
    if u >= 0.f0
        u = u
    else
        u = u + 1.f0
    end
    if abs(point.z)>1.f0
        v=1.f0
    else
        v=Float32(acos(point.z) / pi)
    end
    return Vec2D( u , v )
end

function _sphere_normal(point::Point, ray_dir::Vec)
    result = Normal(point.x, point.y, point.z)
    if toVec(point)*ray_dir < 0.f0
        return result
    else
        return Normal(-1.f0*result.x,-1.f0*result.y,-1.f0*result.z)
    end
end

"""
    Plane(T,M)

It creates a **Plane**, where T is a generic ``Transformation`` applied to the XY palane, M is the ``Material`` of the plane.
"""
struct Plane <: Shape
    transformation::Transformation
    material::Material

    Plane(transformation::Transformation=Transformation(),
          material::Material=Material()) =
          new(transformation, material)
end

## Hidden methods for plane
function _plane_point_to_uv(point::Point)
    u = Float32(point.x - floor(point.x))
    v = Float32(point.y - floor(point.y))
    return Vec2D(u,v)
end

function _plane_normal(point::Point, ray_dir::Vec)
    result = Normal(0.f0, 0.f0, 1.f0)
    Vec(0.f0, 0.f0, 1.f0) * ray_dir < 0.f0 ? nothing : result = Normal(0.f0, 0.f0, -1.f0)
    return result
end

"""
    AAB(T,M)

It creates a **axis-aligned box**, where T is a generic ``Transformation`` and  M is the ``Material`` of the cube.
By default the AAB is generated with the rear-bottom left vertex in (0,0,0), while the front-top right vertex is in (1,1,1).
"""
struct AAB <: Shape
    transformation::Transformation
    material::Material

    AAB(transformation::Transformation=Transformation(),
        material::Material=Material()) =
        new(transformation, material)
end

## Hidden methods for AAB
function _cube_point_to_uv(point::Point)
    if point.x == 0
        u = point.z
        v = 2.f0 - point.y
    elseif point.x == 1
        u = 3.f0 - point.z
        v = 2.f0 - point.y
    elseif point.z == 0
        u = 4.f0 - point.x
        v = 2.f0 - point.y
    elseif point.z == 1
        u = 1.f0 + point.x
        v = 2.f0 - point.y
    elseif point.y == 0
        u = 1.f0 + point.x
        v = 3.f0 - point.z
    elseif point.y == 1
        u = 1.f0 + point.x
        v = point.z
    else
        u = 12.f0
        v = 12.f0
    end
    return Vec2D(u/4.f0, v/3.f0)
end

function _cube_normal(point::Point, ray_dir::Vec)
    if point.x == 0 || point.x == 1
        result = Vec(1.f0, 0.f0, 0.f0)
    elseif point.y == 0 || point.y == 1
        result = Vec(0.f0, 1.f0, 0.f0)
    elseif point.z == 0 || point.z == 1
        result = Vec(0.f0, 0.f0, 1.f0)
    else
        result = Vec(9,9,9)############!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    end
    result * ray_dir < 0.f0 ? nothing : result = result * -1.f0
    return Normal(result.x, result.y, result.z)
end

"""
    LightPoint(P,M)

It creates a **point light source** in the ``Position`` P. M is the ``Material`` of the object, the only important fild is the ``emettide_radiance``.
By default the source is placed in (0.,0.,0.) and has an ``emettide_radiance`` RGB(1.,1,.1.).
"""

struct LightPoint <: Shape
    point::Point
    material::Material

    LightPoint(point::Point=Point(0.,0.,0.),
           material::Material = Material()) =
           new(point, material)
end
Base.:≈(S1::LightPoint,S2::LightPoint) = S1.point ≈ S2.point && S1.material ≈ S2.material

## Code for HITRECORD ###########################################################################################################################

"""
    HitRecord(world_point, normal, surface_point, t, ray, shape)

## Arguments:
- world point;
- normal;
- surface point (u & v coordinates);
- t (distance covered by the ray);
- ray;
- shape.
"""
struct HitRecord
    world_point::Point
    normal::Normal
    surface_point::Vec2D
    t::Float32
    ray::Ray
    shape::Shape
end

Base.:≈(H1::HitRecord,H2::HitRecord) = H1.world_point≈H2.world_point && H1.normal≈H2.normal && H1.surface_point≈H2.surface_point && H1.t ≈ H2.t && H1.ray≈H2.ray && H1.shape ≈ H2.shape
#Base.:≈(H1::HitRecord,H2::HitRecord) = H1.world_point≈H2.world_point && H1.normal≈H2.normal && H1.surface_point≈H2.surface_point && isapprox(H1.t,H2.t, atol=5*10^-2) && H1.ray≈H2.ray && H1.shape ≈ H2.shape
Base.:≈(::Nothing,H2::HitRecord) = false

## Definition of WORLD ############################################################################################################################################

"""
    World()

It creates a **World** with an array of shapes.
"""
struct World
    shapes::Array{Shape}
    World() = new([])
end

# Adding shapes method
function add_shape(world::World, shape::Shape)
    push!(world.shapes, shape)
end

## RAY INTERSECTION ###############################################################################################################################

"""
    ray_intersection(shape, ray)

Evaluates the intersection between a shape (or a World) and a ray, returning a HitRecord.
"""
struct ShapeError <: Exception
    msg::String
end

function ray_intersection(shape::Shape, ray::Ray)
    throw(ShapeError("Ray intersection whit shape $shape are not supported"))
end

## Overloading for ray_intersection with the struct World
function ray_intersection(world::World, ray::Ray)
    closest = nothing
    for i ∈ 1:length(world.shapes)
        intersection = ray_intersection(world.shapes[i], ray)

        if intersection === nothing
            continue
        elseif closest === nothing  || (intersection.t < closest.t)
            closest = intersection
        end

    end
    return closest

end

function ray_intersection(sphere::Sphere, ray::Ray)
    inverse_ray= inverse(sphere.transformation) * ray
    origin_vec = toVec(inverse_ray.origin)
    a = Float32(squared_norm(inverse_ray.dir))
    b = Float32(2.0 * origin_vec * inverse_ray.dir)
    c = Float32(squared_norm(origin_vec) - 1.f0)
    Δ = Float32(b * b - 4 * a * c)

    if Δ<0
        return nothing
    else
        t_1 = Float32(( -b - sqrt(Δ) ) / (2.0 * a))
        t_2 = Float32(( -b + sqrt(Δ) ) / (2.0 * a))
        if t_1 > inverse_ray.tmin && t_1 < inverse_ray.tmax
            first_hit_t = t_1
        elseif t_2 > inverse_ray.tmin && t_2 < inverse_ray.tmax
            first_hit_t = t_2
        else
            return nothing
        end

        hit_point = at(inverse_ray, first_hit_t)
    end

    return HitRecord(sphere.transformation * hit_point, sphere.transformation * _sphere_normal(hit_point, inverse_ray.dir), _sphere_point_to_uv(hit_point), first_hit_t, ray, sphere )
end

function ray_intersection(plane::Plane, ray::Ray)
    inverse_ray= inverse(plane.transformation) * ray
    origin_vec = toVec(inverse_ray.origin)

    if inverse_ray.dir.z ≈ 0
        return nothing
    else
        t = - inverse_ray.origin.z / inverse_ray.dir.z
        if inverse_ray.tmin < t < inverse_ray.tmax
            hit_point = at(inverse_ray, t)
        else
            return nothing
        end
    end
    return HitRecord(plane.transformation * hit_point, plane.transformation * _plane_normal(hit_point, inverse_ray.dir), _plane_point_to_uv(hit_point), t, ray, plane )
end

function ray_intersection(cube::AAB, ray::Ray)
    inverse_ray= inverse(cube.transformation) * ray
    origin_vec = toVec(inverse_ray.origin)
    t_=zeros(2)

    t_min ,t_max = sort( [- inverse_ray.origin.x / inverse_ray.dir.x , (1. - inverse_ray.origin.x) / inverse_ray.dir.x] )
    t_ymin ,t_ymax = sort( [- inverse_ray.origin.y / inverse_ray.dir.y , (1. - inverse_ray.origin.y) / inverse_ray.dir.y] )

    if (t_min > t_ymax) || (t_ymin > t_max)
        return nothing
    end
    t_min = max(t_min, t_ymin)
    t_max = min(t_max, t_ymax)

    t_zmin ,t_zmax = sort( [- inverse_ray.origin.z / inverse_ray.dir.z , (1. - inverse_ray.origin.z) / inverse_ray.dir.z] )

    if (t_min > t_zmax) || (t_zmin > t_max)
        return nothing
    end
    t_min = max(t_min, t_zmin)
    t_max = min(t_max, t_zmax)

    if (inverse_ray.tmin ≤ t_min ≤ inverse_ray.tmax)
        hit_point = at(inverse_ray, t_min)
        return HitRecord(cube.transformation * hit_point, cube.transformation * _cube_normal(hit_point, inverse_ray.dir), _cube_point_to_uv(hit_point), t_min, ray, cube)
    elseif (inverse_ray.tmin ≤ t_max ≤ inverse_ray.tmax)
        hit_point = at(inverse_ray, t_max)
        return HitRecord(cube.transformation * hit_point, cube.transformation * _cube_normal(hit_point, inverse_ray.dir), _cube_point_to_uv(hit_point), t_max, ray, cube)
    else
        return nothing
    end
end

function ray_intersection(light::LightPoint, ray::Ray)
    return nothing
end


text:
# Checks if the solution is a solution that can be checked
# * right number of elements
# * values in [1,data.nbMachines]
function checkVector(data::DataGoogle, solution::SolutionGoogle, verbose::Bool)
    ### Checking that the solution vector is the right one
    if verbose
        println("Checking the size of the vector")
    end

    if size(solution.assignment,1) != data.nbProcess
        println("Vector solution does not have the correct size")
        println("Size = ", size(solution,1), ", nb process = ", data.nbProcess)
        return false
    else
        if verbose
            println("Test passed")
        end
    end

    println("Checking that the values in the solution are machine indices")
    for i in 1:size(solution.assignment,1)
        if solution.assignment[i] > data.nbMachines || solution.assignment[i] < 0
            println("The machine id is wrong for index ", i, " : value =", solution.assignment[i])
            return false
        end
    end

    if verbose
        println("Test passed")
    end

    return true
end

function checkCapacity(data::DataGoogle, solution::SolutionGoogle, verbose::Bool)
    if verbose
        println("Checking capacity constraint (without transient)")
    end

    # create a vector with 0 consumption for each machine * resource
    resourceConsumption = zeros(data.nbMachines, data.nbResources)

    ok = true
    for r in 1:data.nbResources # for each resource

        # computing the resource consumption of the processes on each machine
        for p in 1:data.nbProcess   # for each process
            # add its resource consumption for the machine to which it is assigned
            resourceConsumption[solution.assignment[p],r] += data.processReq[p,r]
        end

        # checking hard resource consumption for each machine
        for m in 1:data.nbMachines
            if resourceConsumption[m,r] > data.hardCapacities[m,r]
                println("Capacity violation for machine ", m, " and resource ", r)
                println(resourceConsumption[m,r], " > ", data.hardCapacities[m,r])
                ok = false
            end
        end
    end

    if verbose && ok
        println("Test passed")
    end

    if verbose
        println("Checking capacity constraint (with transient)")
    end
    for r in 1:data.nbResources # for each resource
        if data.transientStatus[r]==1    # if it is transient
            # compute the resource consumption (including processes that were originally assigned to the machine )
            for p in 1:data.nbProcess
                if solution.assignment[p] != data.initialAssignment[p] # if the process has moved
                    resourceConsumption[data.initialAssignment[p],r] += data.processReq[p,r] # add its resource consumption
                end
            end

            # checking hard resource consumption for each machine (including transient processes)
            for m in 1:data.nbMachines
                if resourceConsumption[m,r] > data.hardCapacities[m,r]
                    println("When transient resources considered, capacity violation for machine ", m, " and resource ", r)
                    println(resourceConsumption[m,r], " > ", data.hardCapacities[m,r])
                    ok = false
                end
            end
        end
    end

    if verbose && ok
        println("Test passed")
    end

    return ok
end

function checkDisjunctions(data::DataGoogle, solution::SolutionGoogle, verbose::Bool)

    if verbose
        println("Checking disjunctions between processes of the same service")
    end

    ok = true
    for i in 1:data.nbProcess-1, j in i+1:data.nbProcess # for each process j ≠ i
        # if i and j are in the same service and on the same machine there is a problem
        if data.servicesProcess[i] == data.servicesProcess[j]  && solution.assignment[i] == solution.assignment[j]
            println("Disjunction constraint violation between processes ", i, " and ", j)
            println("Both belong to service",  data.servicesProcess[i], " and are assigned to machine ", solution.assignment[i])
            ok = false
        end
    end

    if verbose && ok
        println("Test passed")
    end

    return ok
end

function checkSpread(data::DataGoogle, solution::SolutionGoogle, verbose::Bool)
    if verbose
        println("Checking spread constraints ")
    end
    ok = true

    for s in 1:data.nbServices # for each service
        if data.spreadMin[s] != 0 # if it has a spread min

            locationsWithService = zeros(data.nbLocations)  # table used to know which locations have service s
            for p in 1:data.nbProcess   # for each process
                locationsWithService[data.locations[solution.assignment[p]]] = 1    # note that its current location has its service
            end
            nbLocationsWithService = sum(locationsWithService)  # number of locations with this service

            if nbLocationsWithService < data.spreadMin[s]   # should be larger than the spread
                println("Spread constraint violated for service ", s)
                println("Present in ", nbLocationsWithService, " locations < ", data.spreadMin[s])
                ok = false
            end
        end
    end

    if verbose && ok
        println("Test passed")
    end

    return ok
end

function checkDependencies(data::DataGoogle, solution::SolutionGoogle, verbose::Bool)
    println("Checking dependency constraints")
    ok = true

    # table indicating for each neighborhood if it currently has the service
    servicesInNeighborhood = zeros(data.nbServices,data.nbNeighborhoods)

    for p in 1:data.nbProcess
        servicesInNeighborhood[data.servicesProcess[p],data.neighborhoods[solution.assignment[p]]] = 1
    end

    for s in 1:data.nbServices, n in 1:data.nbNeighborhoods, t in data.dependences[s]
        # there cannot be s and not t
        if servicesInNeighborhood[s,n] == 1 && servicesInNeighborhood[t,n] == 0
            println("Dependency constraint violation")
            println("Service ", s, " is in neighbordhood ", n, " while service ", t, " is not")
            ok = false
        end
    end

    if verbose && ok
        println("Test passed")
    end

    return ok
end

function computeLoadCost(data::DataGoogle, solution::SolutionGoogle, verbose::Bool)
    val = 0

    # resource consumption for each machine and each resource
    resourceConsumption = zeros(data.nbMachines, data.nbResources)
    for r in 1:data.nbResources,  p in 1:data.nbProcess
            resourceConsumption[solution.assignment[p],r] += data.processReq[p,r]
    end

    # compute the load cost by summing the cost over machines and resources
    for m in 1:data.nbMachines, r in 1:data.nbResources
        # the cost is accounted only if is larger than 0
        val += data.weightLoadCost[r] * max(0, resourceConsumption[m,r] - data.softCapacities[m,r])
    end

    return val
end

function computeBalanceCost(data::DataGoogle, solution::SolutionGoogle, verbose::Bool)
    totalBalanceCost = 0

    # for each balance cost data
    for b in data.balanceCostDataList
        # table that will contain the remaining capacity for resource r1 on each machine
        slack_r1 = data.hardCapacities[:,b.resource1]
        # table that will contain the remaining capacity for resource r2 on each machine
        slack_r2 = data.hardCapacities[:,b.resource2]

        for p in 1:data.nbProcess   # for each process, remove its resource requirement from the remaining capacity
            slack_r1[solution.assignment[p]] -= data.processReq[p,b.resource1]
            slack_r2[solution.assignment[p]] -= data.processReq[p,b.resource2]
        end

        localBalanceCost = 0    # balance cost for the current value of b
        for m in 1:data.nbMachines  # sum over all machines of the balance costs
            localBalanceCost += max(0, b.target * slack_r1[m] - slack_r2[m])
        end

        totalBalanceCost += b.weight*localBalanceCost
    end

    return totalBalanceCost
end

function computeProcessMoveCost(data::DataGoogle, solution::SolutionGoogle, verbose::Bool)
    val = 0
    for p in 1:data.nbProcess   # count the number of processses that are not in the same machine after optimizing
        if data.initialAssignment[p] != solution.assignment[p]
            val += data.processMoveCost[p]
        end
    end

    return val
end

function computeMachineMoveCost(data::DataGoogle, solution::SolutionGoogle, verbose::Bool)
    val = 0
    for p in 1:data.nbProcess
        val += data.machineMoveCosts[data.initialAssignment[p],solution.assignment[p]]
    end

    return val
end

function computeServiceMoveCost(data::DataGoogle, solution::SolutionGoogle, verbose::Bool)
    nbMovesInService = zeros(data.nbServices)   # number of processes that moved inside each service

    for p in 1:data.nbServices  # for each process
        if data.initialAssignment[p] != solution.assignment[p]  # if it moved
            nbMovesInService[data.servicesProcess[p]] += 1  # increase the corresponding service counter
        end
    end

    return maximum(nbMovesInService)    # return the maximum among all services
end

function checkAndComputeCost(data::DataGoogle, solution::SolutionGoogle, verbose::Bool)

    ### Printing the arguments
    if verbose
        println("Checking the solution ")
        print("Solution received is ")
        println(solution.assignment)
        println("Cost : ", solution.cost)
    end

    if ! checkVector(data,solution,verbose)
        println("Cannot check the solution ")
        return
    end

    ok = true

    ok = checkCapacity(data,solution,verbose) && ok
    ok = checkDisjunctions(data, solution, verbose) && ok
    ok = checkSpread(data, solution, verbose) && ok
    ok = checkDependencies(data, solution, verbose) && ok

    if ok
        println("All tests passed: solution feasible ")
    end

    loadCost        = computeLoadCost(data, solution, verbose)
    println("Load cost = ", loadCost)
    balanceCost     = computeBalanceCost(data, solution, verbose)
    println("Balance cost = ", balanceCost)
    processMoveCost = computeProcessMoveCost(data, solution, verbose)
    println("Process move cost = ", processMoveCost)
    machineMoveCost = computeMachineMoveCost(data, solution, verbose)
    println("Machine move cost = ", machineMoveCost)
    serviceMoveCost = computeServiceMoveCost(data, solution, verbose)
    println("Service move cost = ", serviceMoveCost)

    totalCost = loadCost + balanceCost + data.processMoveWeight * processMoveCost + data.serviceMoveWeight * serviceMoveCost + data.machineMoveWeight * machineMoveCost

    println("Objective function = ", totalCost)

    println("Objective recorded in the solution ", solution.cost)

    if totalCost != solution.cost
        println("The two costs do not match")
    end

end


text:
# general helper functions for integrators.

_default_solver() = SimpleSolvers.get_config(:nls_solver)

"""
Create nonlinear solver object for a system of `N` equations with data type `DT`.
The function ``f(x)=0`` to be solved for is determined by a julia function
`function_stages!(x, b, params)`, where `x` is the current solution and `b` is
the output vector, s.th. ``b = f(x)``. `params` are a set of parameters depending
on the equation and integrator that is used.
The solver type is obtained from the config dictionary (`:nls_solver`).
"""
function create_nonlinear_solver(DT, N, params, caches, solver::Type{<:NonlinearSolver}=_default_solver(), F::Function=function_stages!)
    # create solution vector for nonlinear solver
    x = zeros(DT, N)

    # create wrapper function f!(x,b) that calls `F(x, b, params)`
    # with the appropriate `params`
    f! = (x,b) -> F(x, b, params, caches)

    # create nonlinear solver with solver type obtained from config dictionary
    s = solver(x, f!)
end

function create_nonlinear_solver(DT, N, params, caches, i::Int, solver::Type{<:NonlinearSolver}=_default_solver(), F::Function=function_stages!)
    # create solution vector for nonlinear solver
    x = zeros(DT, N)

    # create wrapper function f(x,b) that calls `F(x, b, params)`
    # with the appropriate `params`
    f = (x,b) -> F(x, b, params, caches, i)

    # create nonlinear solver with solver type obtained from config dictionary
    s = solver(x, f)
end

function create_nonlinear_solver_with_jacobian(DT, N, params, caches, solver::Type{<:NonlinearSolver}=_default_solver(), F=function_stages!, J=jacobian!)
    # create solution vector for nonlinear solver
    x = zeros(DT, N)

    # create wrapper function f!(x,b) that calls `F(x, b, params)`
    # with the appropriate `params`
    f! = (x,b) -> F(x, b, params, caches)

    # create wrapper function j!(x,df) that calls `J(x, df, params)`
    # with the appropriate `params`
    j! = (x,df) -> J(x, df, caches[DT], params)

    # create nonlinear solver with solver type obtained from config dictionary
    s = solver(x, f!; J! = j!)
end

"""
Create a vector of S solution vectors of type DT to store the solution of S
internal stages for a problem with `D` dimensions.
"""
function create_internal_stage_vector(DT, D, S)
    [zeros(DT,D) for i in 1:S]
end

"""
Create a vector of S solution matrices of type DT to store the solution of S
internal stages for a problem with `DxM` dimensions.
"""
function create_internal_stage_matrix(DT, D, M, S)
    [zeros(DT,D,M) for i in 1:S]
end

"""
Create a vector of S+1 solution vectors of type DT to store the solution of S
internal stages and the solution of the previous timestep for a problem with `D`
    dimensions.
"""
function create_internal_stage_vector_with_zero(DT, D, S)
    a = OffsetArray{Vector{DT}}(undef, 0:S)

    for i in 0:S
        a[i] = zeros(DT,D)
    end

    return a
end

"""
Create a vector of (S,M+1) solution vectors of type DT to store the solution of S
internal stages and M random processes for a problem with `D` dimensions.
"""
function create_internal_stage_vector_with_zero(DT, D, M, S)
    a = OffsetArray{Vector{DT}}(undef, 0:M, 1:S)

    for i in 0:M
        for j in 1:S
            a[i,j] = zeros(DT,D)
        end
    end

    return a
end

function update_solution!(x::Vector{T}, xₑᵣᵣ::Vector{T}, ẋ::Matrix{T}, b::Vector{T}, Δt::T) where {T}
    @assert length(x) == length(xₑᵣᵣ)
    @assert length(x) == size(ẋ, 1)
    @assert length(b) == size(ẋ, 2)

    for k in axes(ẋ, 1)
        for i in axes(ẋ, 2)
            x[k], xₑᵣᵣ[k] = compensated_summation(Δt * b[i] * ẋ[k,i], x[k], xₑᵣᵣ[k])
        end
    end
end

function update_solution!(x::Vector{T}, xₑᵣᵣ::Vector{T}, ẋ::Vector{Vector{T}}, b::Vector{T}, Δt::T) where {T}
    @assert length(b) == length(ẋ)
    @assert length(x) == length(ẋ[1])
    @assert length(x) == length(xₑᵣᵣ)

    for i in eachindex(ẋ)
        for k in eachindex(ẋ[i])
            x[k], xₑᵣᵣ[k] = compensated_summation(Δt * b[i] * ẋ[i][k], x[k], xₑᵣᵣ[k])
        end
    end
end

function update_solution!(x::SolutionVector{T}, ẋ::Matrix{T}, b::Vector{T}, Δt::T) where {T}
    @assert length(x) == size(ẋ, 1)
    @assert length(b) == size(ẋ, 2)

    local Δx::eltype(x)

    for k in axes(ẋ, 1)
        Δx = 0
        for i in axes(ẋ, 2)
            Δx += b[i] * ẋ[k,i]
        end
        x[k] += Δt * Δx
    end
end

function update_solution!(x::SolutionVector{T}, ẋ::Vector{Vector{T}}, b::Vector{T}, Δt::T) where {T}
    @assert length(b) == length(ẋ)
    @assert length(x) == length(ẋ[1])

    local Δx::T

    for k in eachindex(x)
        Δx = 0
        for i in eachindex(ẋ)
            Δx += b[i] * ẋ[i][k]
        end
        x[k] += Δt * Δx
    end
end

function update_solution!(x::Vector{T}, xₑᵣᵣ::Vector{T}, ẋ::Union{Matrix{T},Vector{Vector{T}}}, b::Vector{T}, b̂::Vector, Δt::T) where {T}
    update_solution!(x, xₑᵣᵣ, ẋ, b, Δt)
    update_solution!(x, xₑᵣᵣ, ẋ, b̂, Δt)
end

function update_solution!(x::SolutionVector{T}, ẋ::Union{Matrix{T},Vector{Vector{T}}}, b::Vector{T}, b̂::Vector, Δt::T) where {T}
    update_solution!(x, ẋ, b, Δt)
    update_solution!(x, ẋ, b̂, Δt)
end

function update_multiplier!(λ::SolutionVector{T}, Λ::Vector{Vector{T}}, b::Vector{T}) where {T}
    for Λᵢ in Λ
        @assert length(λ) == length(Λᵢ)
    end
    local t::T
    @inbounds for i in eachindex(λ)
        t = zero(T)
        for j in eachindex(b,Λ)
            t += b[j] * Λ[j][i]
        end
        λ[i] = t
    end
    nothing
end

function GeometricBase.cut_periodic_solution!(x::SolutionVector{T}, periodicity::Vector{T}) where {T}
    @assert length(x) == length(periodicity)

    for k in eachindex(x, periodicity)
        if periodicity[k] ≠ 0
            while x[k] < 0
                x[k] += periodicity[k]
            end
            while x[k] ≥ periodicity[k]
                x[k] -= periodicity[k]
            end
        end
    end
end

function GeometricBase.cut_periodic_solution!(x::SolutionVector{T}, periodicity::Vector{T}, shift::Vector{T}) where {T}
    @assert length(x) == length(periodicity)
    shift .= 0
    for k in eachindex(x, periodicity, shift)
        if periodicity[k] ≠ 0
            while x[k] + shift[k] < 0
                shift[k] += periodicity[k]
            end
            while x[k] + shift[k] ≥ periodicity[k]
                shift[k] -= periodicity[k]
            end
        end
    end
end


text:
module MacrosTests

using ReverseDiff, ForwardDiff, Test, StaticArrays
using ForwardDiff: Dual, Partials, partials

include(joinpath(dirname(@__FILE__), "utils.jl"))

tp = InstructionTape()
x, a, b = rand(3)

############
# @forward #
############

f0(x) = 1. / (1. + exp(-x))
f0(a, b) = sqrt(a^2 + b^2)

ReverseDiff.@forward f1(x::T) where {T<:Real} = 1. / (1. + exp(-x))
ReverseDiff.@forward f1(a::A, b::B) where {A,B<:Real} = sqrt(a^2 + b^2)

ReverseDiff.@forward f2(x) = 1. / (1. + exp(-x))
ReverseDiff.@forward f2(a, b) = sqrt(a^2 + b^2)

ReverseDiff.@forward function f3(x::T) where T<:Real
    return 1. / (1. + exp(-x))
end

ReverseDiff.@forward function f3(a::A, b::B) where {A,B<:Real}
    return sqrt(a^2 + b^2)
end

ReverseDiff.@forward function f4(x)
    return 1. / (1. + exp(-x))
end

ReverseDiff.@forward function f4(a, b)
    return sqrt(a^2 + b^2)
end

function test_forward(f, x, tp)
    xt = ReverseDiff.TrackedReal(x, zero(x), tp)

    y = f(x)
    @test isempty(tp)

    yt = f(xt)
    @test yt == y
    dual = f(Dual(x, one(x)))
    @test length(tp) == 1
    instruction = first(tp)
    @test typeof(instruction) <: ReverseDiff.ScalarInstruction
    @test instruction.input === xt
    @test instruction.output === yt
    @test instruction.cache[] === partials(dual, 1)
    empty!(tp)
end

function test_forward(f, a, b, tp)
    at = ReverseDiff.TrackedReal(a, zero(a), tp)
    bt = ReverseDiff.TrackedReal(b, zero(b), tp)

    c = f(a, b)
    dual = f(Dual(a, one(a), zero(a)), Dual(b, zero(b), one(b)))
    @test isempty(tp)

    tc = f(at, b)
    @test tc == c
    @test length(tp) == 1
    instruction = first(tp)
    @test typeof(instruction) <: ReverseDiff.ScalarInstruction
    @test instruction.input === (at, b)
    @test instruction.output === tc
    @test instruction.cache[] === SVector(partials(dual, 1), partials(dual, 1))
    empty!(tp)

    tc = f(a, bt)
    @test tc == c
    @test length(tp) == 1
    instruction = first(tp)
    @test typeof(instruction) <: ReverseDiff.ScalarInstruction
    @test instruction.input === (a, bt)
    @test instruction.output === tc
    @test instruction.cache[] === SVector(partials(dual, 2), partials(dual, 2))
    empty!(tp)

    tc = f(at, bt)
    @test tc == c
    @test length(tp) == 1
    instruction = first(tp)
    @test typeof(instruction) <: ReverseDiff.ScalarInstruction
    @test instruction.input === (at, bt)
    @test instruction.output === tc
    @test instruction.cache[] === SVector(partials(dual)...)
    empty!(tp)
end

for f in (ReverseDiff.@forward(f0), f1, f2, f3, f4, ReverseDiff.@forward(-))
    test_println("@forward named functions", f)
    test_forward(f, x, tp)
    test_forward(f, a, b, tp)
end

ReverseDiff.@forward f5 = (x) -> 1. / (1. + exp(-x))
test_println("@forward anonymous functions", f5)
test_forward(f5, x, tp)

ReverseDiff.@forward f6 = (a, b) -> sqrt(a^2 + b^2)
test_println("@forward anonymous functions", f6)
test_forward(f6, a, b, tp)

#########
# @skip #
#########

g0 = f0

ReverseDiff.@skip g1(x::T) where {T<:Real} = 1. / (1. + exp(-x))
ReverseDiff.@skip g1(a::A, b::B) where {A,B<:Real} = sqrt(a^2 + b^2)

ReverseDiff.@skip g2(x) = 1. / (1. + exp(-x))
ReverseDiff.@skip g2(a, b) = sqrt(a^2 + b^2)

ReverseDiff.@skip function g3(x::T) where T<:Real
    return 1. / (1. + exp(-x))
end

ReverseDiff.@skip function g3(a::A, b::B) where {A,B<:Real}
    return sqrt(a^2 + b^2)
end

ReverseDiff.@skip function g4(x)
    return 1. / (1. + exp(-x))
end

ReverseDiff.@skip function g4(a, b)
    return sqrt(a^2 + b^2)
end

function test_skip(g, x, tp)
    xt = ReverseDiff.TrackedReal(x, zero(x), tp)

    y = g(x)
    @test isempty(tp)

    yt = g(xt)
    @test yt === y
    @test isempty(tp)
end

function test_skip(g, a, b, tp)
    at = ReverseDiff.TrackedReal(a, zero(a), tp)
    bt = ReverseDiff.TrackedReal(b, zero(b), tp)

    c = g(a, b)
    @test isempty(tp)

    tc = g(at, b)
    @test tc === c
    @test isempty(tp)

    tc = g(a, bt)
    @test tc === c
    @test isempty(tp)

    tc = g(at, bt)
    @test tc === c
    @test isempty(tp)
end

for g in (ReverseDiff.@skip(g0), g1, g2, g3, g4)
    test_println("@skip named functions", g)
    test_skip(g, x, tp)
    test_skip(g, a, b, tp)
end

ReverseDiff.@skip g5 = (x) -> 1. / (1. + exp(-x))
test_println("@skip anonymous functions", g5)
test_skip(g5, x, tp)

ReverseDiff.@skip g6 = (a, b) -> sqrt(a^2 + b^2)
test_println("@skip anonymous functions", g6)
test_skip(g6, a, b, tp)

end # module


text:
function updatediscretization!(system::LegSystem)
    # determine time shift
    system.solverparams.tshift = system.t[end]-
        sum(system.solverparams.th*system.solverparams.numbeats);

    # update total number of time steps
    ntoadd = ceil((system.solverparams.th-
        system.solverparams.tshift)/system.solverparams.h);
    system.solverparams.numsteps+=ntoadd;

    # update discrete times
    ttoadd = [0:1:ntoadd;]*system.solverparams.h + system.t[end];
    append!(system.t,ttoadd[2:end]);

    # allocate additional space for solution variables
    for i = 1:length(system.branches.ID)
        system.branches.A[i] = [system.branches.A[i] zeros(system.solverparams.JL,length(ttoadd[2:end]))];
        system.branches.Q[i] = [system.branches.Q[i] zeros(system.solverparams.JL,length(ttoadd[2:end]))];
        system.branches.P[i] = [system.branches.P[i] zeros(system.solverparams.JL,length(ttoadd[2:end]))];
        if isempty(system.branches.children[i])
            system.branches.term[i].P = [system.branches.term[i].P;zeros(length(ttoadd[2:end]),1)];
            system.branches.term[i].Q = [system.branches.term[i].Q;zeros(length(ttoadd[2:end]),1)];
        end
    end

end

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics