text:
struct ODEInterpolatingAdjointSensitivityFunction{C<:AdjointDiffCache,Alg<:InterpolatingAdjoint,
                                                  uType,SType,CPS,pType,fType<:DiffEqBase.AbstractDiffEqFunction} <: SensitivityFunction
  diffcache::C
  sensealg::Alg
  discrete::Bool
  y::uType
  sol::SType
  checkpoint_sol::CPS
  prob::pType
  f::fType
  noiseterm::Bool
end

mutable struct CheckpointSolution{S,I,T}
  cpsol::S # solution in a checkpoint interval
  intervals::I # checkpoint intervals
  cursor::Int # sol.prob.tspan = intervals[cursor]
  tols::T
end

function ODEInterpolatingAdjointSensitivityFunction(g,sensealg,discrete,sol,dg,f,checkpoints,tols;noiseterm=false)
  tspan = reverse(sol.prob.tspan)
  checkpointing = ischeckpointing(sensealg, sol)
  (checkpointing && checkpoints === nothing) && error("checkpoints must be passed when checkpointing is enabled.")

  checkpoint_sol = if checkpointing
    intervals = map(tuple, @view(checkpoints[1:end-1]), @view(checkpoints[2:end]))
    interval_end = intervals[end][end]
    tspan[1] > interval_end && push!(intervals, (interval_end, tspan[1]))
    cursor = lastindex(intervals)
    interval = intervals[cursor]

    if typeof(sol.prob) <: SDEProblem
      # replicated noise
      _sol = deepcopy(sol)
      sol.W.save_everystep = false
      _sol.W.save_everystep = false
      idx1 = searchsortedfirst(_sol.t, interval[1]-1000eps(interval[1]))
      #idx2 = searchsortedfirst(_sol.t, interval[2])
      #forwardnoise = DiffEqNoiseProcess.NoiseGrid(_sol.t[idx1:idx2], _sol.W.W[idx1:idx2])
      forwardnoise = DiffEqNoiseProcess.NoiseWrapper(_sol.W, indx=idx1)
      dt = abs(_sol.W.dt)
      if dt < 1000eps(_sol.t[end])
        dt = interval[2] - interval[1]
      end
      cpsol = solve(remake(sol.prob, tspan=interval, u0=sol(interval[1]), noise=forwardnoise), sol.alg, save_noise=false; dt=dt, tstops=_sol.t[idx1:end] ,tols...)
    else
      cpsol = solve(remake(sol.prob, tspan=interval, u0=sol(interval[1])), sol.alg; tols...)
    end
    CheckpointSolution(cpsol, intervals, cursor, tols)
  else
    nothing
  end

  diffcache, y = adjointdiffcache(g,sensealg,discrete,sol,dg,f;quad=false,noiseterm=noiseterm)

  return ODEInterpolatingAdjointSensitivityFunction(diffcache,sensealg,
                                                    discrete,y,sol,
                                                    checkpoint_sol,sol.prob,f,noiseterm)
end

function findcursor(intervals, t)
  # equivalent with `findfirst(x->x[1] <= t <= x[2], intervals)`
  lt(x, t) = <(x[2], t)
  return searchsortedfirst(intervals, t, lt=lt)
end

# u = λ'
# add tstop on all the checkpoints
function (S::ODEInterpolatingAdjointSensitivityFunction)(du,u,p,t)
  @unpack sol, y, checkpoint_sol, discrete, prob, f = S
  idx = length(y)

  if checkpoint_sol === nothing
    if typeof(t) <: ForwardDiff.Dual && eltype(S.y) <: AbstractFloat
      y = sol(t)
    else
      sol(y,t)
    end
  else
    intervals = checkpoint_sol.intervals
    interval = intervals[checkpoint_sol.cursor]
    if !(interval[1] <= t <= interval[2])
      cursor′ = findcursor(intervals, t)
      interval = intervals[cursor′]
      cpsol_t = checkpoint_sol.cpsol.t
      if typeof(t) <: ForwardDiff.Dual && eltype(S.y) <: AbstractFloat
        y = sol(interval[1])
      else
        sol(y, interval[1])
      end
      if typeof(sol.prob) <: SDEProblem
        #idx1 = searchsortedfirst(sol.t, interval[1])
        _sol = deepcopy(sol)
        _sol.W.save_everystep = false
        idx1 = searchsortedfirst(_sol.t, interval[1]-100eps(interval[1]))
        idx2 = searchsortedfirst(_sol.t, interval[2])
        #forwardnoise = DiffEqNoiseProcess.NoiseGrid(_sol.t[idx1:idx2], _sol.W.W[idx1:idx2])
        forwardnoise = DiffEqNoiseProcess.NoiseWrapper(_sol.W, indx=idx1)
        prob′ = remake(prob, tspan=intervals[cursor′], u0=y, noise=forwardnoise)
        dt = abs(cpsol_t[end]-cpsol_t[end-1])
        if dt < 10000eps(cpsol_t[end])
          dt = interval[2] - interval[1]
        end
        cpsol′ = solve(prob′, sol.alg, noise=forwardnoise, save_noise=false; dt=dt, tstops=_sol.t[idx1:idx2], checkpoint_sol.tols...)
      else
        prob′ = remake(prob, tspan=intervals[cursor′], u0=y)
        cpsol′ = solve(prob′, sol.alg; dt=abs(cpsol_t[end] - cpsol_t[end-1]), checkpoint_sol.tols...)
      end
      checkpoint_sol.cpsol = cpsol′
      checkpoint_sol.cursor = cursor′
    end
    checkpoint_sol.cpsol(y, t)
  end

  λ     = @view u[1:idx]
  grad  = @view u[idx+1:end]

  if length(u) == length(du)
    dλ    = @view du[1:idx]
    dgrad = @view du[idx+1:end]

  elseif length(u) != length(du) &&  StochasticDiffEq.is_diagonal_noise(prob) && !isnoisemixing(S.sensealg)
    idx1 = [length(u)*(i-1)+i for i in 1:idx] # for diagonal indices of [1:idx,1:idx]

    dλ    = @view du[idx1]
    dgrad = @view du[idx+1:end,1:idx]

  else
    # non-diagonal noise and noise mixing case
    dλ    = @view du[1:idx,1:idx]
    dgrad = @view du[idx+1:end,1:idx]
  end

  if S.noiseterm
    if length(u) == length(du)
      vecjacobian!(dλ, y, λ, p, t, S, dgrad=dgrad)
    elseif length(u) != length(du) &&  StochasticDiffEq.is_diagonal_noise(prob) && !isnoisemixing(S.sensealg)
      vecjacobian!(dλ, y, λ, p, t, S)
      jacNoise!(λ, y, p, t, S, dgrad=dgrad)
    else
      jacNoise!(λ, y, p, t, S, dgrad=dgrad, dλ=dλ)
    end
  else
    vecjacobian!(dλ, y, λ, p, t, S, dgrad=dgrad)
  end

  dλ .*= -one(eltype(λ))

  discrete || accumulate_cost!(dλ, y, p, t, S, dgrad)
  return nothing
end

# g is either g(t,u,p) or discrete g(t,u,i)
@noinline function ODEAdjointProblem(sol,sensealg::InterpolatingAdjoint,
                                     g,t=nothing,dg=nothing;
                                     checkpoints=sol.t,
                                     callback=CallbackSet(),
                                     reltol=nothing, abstol=nothing,
                                     kwargs...)
  @unpack f, p, u0, tspan = sol.prob
  tspan = reverse(tspan)
  discrete = t != nothing

  numstates = length(u0)
  numparams = p === nothing || p === DiffEqBase.NullParameters() ? 0 : length(p)

  len = numstates+numparams

  λ = p === nothing || p === DiffEqBase.NullParameters() ? similar(u0) : one(eltype(u0)) .* similar(p, len)
  λ .= false

  sense = ODEInterpolatingAdjointSensitivityFunction(g,sensealg,discrete,sol,dg,f,
                                                     checkpoints,
                                                     (reltol=reltol,abstol=abstol))

  init_cb = t !== nothing && tspan[1] == t[end]
  cb = generate_callbacks(sense, g, λ, t, callback, init_cb)
  z0 = vec(zero(λ))
  original_mm = sol.prob.f.mass_matrix
  if original_mm === I || original_mm === (I,I)
    mm = I
  else
    adjmm = copy(sol.prob.f.mass_matrix')
    zzz = similar(adjmm, numstates, numparams)
    fill!(zzz, zero(eltype(zzz)))
    # using concrate I is slightly more efficient
    II = Diagonal(I, numparams)
    mm = [adjmm       zzz
          copy(zzz')   II]
  end

  jac_prototype = sol.prob.f.jac_prototype
  if !sense.discrete || jac_prototype === nothing
    adjoint_jac_prototype = nothing
  else
    _adjoint_jac_prototype = copy(jac_prototype')
    zzz = similar(_adjoint_jac_prototype, numstates, numparams)
    fill!(zzz, zero(eltype(zzz)))
    II = Diagonal(I, numparams)
    adjoint_jac_prototype = [_adjoint_jac_prototype zzz
                             copy(zzz')             II]
  end

  odefun = ODEFunction(sense, mass_matrix=mm, jac_prototype=adjoint_jac_prototype)
  return ODEProblem(odefun,z0,tspan,p,callback=cb)
end

@noinline function SDEAdjointProblem(sol,sensealg::InterpolatingAdjoint,
                                     g,t=nothing,dg=nothing;
                                     checkpoints=sol.t,
                                     callback=CallbackSet(),
                                     reltol=nothing, abstol=nothing,
                                     diffusion_jac=nothing, diffusion_paramjac=nothing,
                                     kwargs...)
  @unpack f, p, u0, tspan = sol.prob
  tspan = reverse(tspan)
  discrete = t != nothing

  if length(unique(round.(checkpoints, digits=13))) != length(checkpoints)
    @warn "The given checkpoints are not unique. To avoid issues in the interpolation the checkpoints were redefined. You may want to check sol.t if default checkpoints were used."
    checkpoints = unique(round.(checkpoints, digits=13))
  end

  p === DiffEqBase.NullParameters() && error("Your model does not have parameters, and thus it is impossible to calculate the derivative of the solution with respect to the parameters. Your model must have parameters to use parameter sensitivity calculations!")
  numstates = length(u0)
  numparams = length(p)

  len = numstates+numparams

  λ = one(eltype(u0)) .* similar(p, len)
  λ .= false

  sense_drift = ODEInterpolatingAdjointSensitivityFunction(g,sensealg,discrete,sol,dg,sol.prob.f,
                                                     checkpoints,(reltol=reltol,abstol=abstol))

  diffusion_function = ODEFunction(sol.prob.g, jac=diffusion_jac, paramjac=diffusion_paramjac)
  sense_diffusion = ODEInterpolatingAdjointSensitivityFunction(g,sensealg,discrete,sol,dg,diffusion_function,
                                                     checkpoints,(reltol=reltol,abstol=abstol);noiseterm=true)

  init_cb = t !== nothing && tspan[1] == t[end]
  cb = generate_callbacks(sense_drift, g, λ, t, callback, init_cb)
  z0 = vec(zero(λ))
  original_mm = sol.prob.f.mass_matrix
  if original_mm === I || original_mm === (I,I)
    mm = I
  else
    adjmm = copy(sol.prob.f.mass_matrix')
    zzz = similar(adjmm, numstates, numparams)
    fill!(zzz, zero(eltype(zzz)))
    # using concrate I is slightly more efficient
    II = Diagonal(I, numparams)
    mm = [adjmm       zzz
          copy(zzz')   II]
  end

  jac_prototype = sol.prob.f.jac_prototype
  if !sense_drift.discrete || jac_prototype === nothing
    adjoint_jac_prototype = nothing
  else
    _adjoint_jac_prototype = copy(jac_prototype')
    zzz = similar(_adjoint_jac_prototype, numstates, numparams)
    fill!(zzz, zero(eltype(zzz)))
    II = Diagonal(I, numparams)
    adjoint_jac_prototype = [_adjoint_jac_prototype zzz
                             copy(zzz')             II]
  end

  sdefun = SDEFunction(sense_drift,sense_diffusion,mass_matrix=mm,jac_prototype=adjoint_jac_prototype)

  # replicated noise
  _sol = deepcopy(sol)
  backwardnoise = DiffEqNoiseProcess.NoiseWrapper(_sol.W, reverse=true)

  if StochasticDiffEq.is_diagonal_noise(sol.prob) && typeof(sol.W[end])<:Number
    # scalar noise case
    noise_matrix = nothing
  else
    noise_matrix = similar(z0,length(z0),numstates)
    noise_matrix .= false
  end

  return SDEProblem(sdefun,sense_diffusion,z0,tspan,p,
    callback=cb,
    noise=backwardnoise,
    noise_rate_prototype = noise_matrix
    )
end


text:
using Documenter, Compose

makedocs(
    modules = [Compose],
    clean = true,
    format = :html,
    sitename = "Compose.jl",
    pages = Any[
        "Home" => "index.md"
        "Library" => "library.md"
    ]
)

deploydocs(
    repo   = "github.com/GiovineItalia/Compose.jl.git",
    julia  = "0.5",
    osname = "linux",
    deps = nothing,
    make = nothing,
    target = "build"
)

ext: jl


text:
# Autogenerated wrapper script for SYMPHONY_jll for x86_64-apple-darwin-libgfortran5
export libosisym, libsym, symphony

using CoinUtils_jll
using Osi_jll
using Clp_jll
using Cgl_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("SYMPHONY")
JLLWrappers.@declare_library_product(libosisym, "@rpath/libOsiSym.3.dylib")
JLLWrappers.@declare_library_product(libsym, "@rpath/libSym.3.dylib")
JLLWrappers.@declare_executable_product(symphony)
function __init__()
    JLLWrappers.@generate_init_header(CoinUtils_jll, Osi_jll, Clp_jll, Cgl_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libosisym,
        "lib/libOsiSym.3.6.17.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsym,
        "lib/libSym.3.6.17.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        symphony,
        "bin/symphony",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using CollegeStratCollege, ModelObjectsLH, ModelParams, Test

cs = CollegeStratCollege;

function learning_test(switches)
    h0 = 1.2;
    hHs = 1.5;
    @testset "Learning" begin
        println(switches);
        objId = ObjectId(:learning);
        ls = cs.init_learning_set(objId, switches);
        pt = ModelParams.param_table(ls, true);
        println(pt);
        for isTwoYear ∈ (false, true)
            l = cs.init_learning(ls, isTwoYear);
            fracH = cs.frac_h(l);
            @test 0.0 <= fracH <= 1.0
            for dh ∈ ([0.5 0.6; 0.7 0.8], 0.6);
                hM, hHsM = cs.h_stocks(l, dh, h0, hHs);
                @test size(hM) == size(hHsM) == size(dh)
                @test all(hM .>= h0)
                @test all(hHsM .>= hHs)
            end
        end
	end
end

@testset "Learning" begin
    for switches ∈ (cs.LearnCollegeOnlySwitches(),
        LearnTwoFourSwitches(0.3, 0.5, true, true, true))

        learning_test(switches);
    end
end

# ----------------


text:
module Logging
    using Logging
    using Colors

    function start_logger()
        printstyled("Logger started!", bold=true, color=colorant"blue")
    end
end

max_stars_repo_name: dustyirwin/PlaymatSimulator

max_issues_repo_name: dustyirwin/PlaymatSimulator

max_forks_repo_name: dustyirwin/PlaymatSimulator


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

# What is included in the file `FMI2_ext.jl` (external/additional functions)?
# - new functions, that are useful, but not part of the FMI-spec (example: `fmi2Load`, `fmi2SampleDirectionalDerivative`)

using Libdl
using ZipFile
import Downloads

"""
Create a copy of the .fmu file as a .zip folder and unzips it.
Returns the paths to the zipped and unzipped folders.

Via optional argument ```unpackPath```, a path to unpack the FMU can be specified (default: system temporary directory).
"""
function fmi2Unzip(pathToFMU::String; unpackPath=nothing)

    fileNameExt = basename(pathToFMU)
    (fileName, fileExt) = splitext(fileNameExt)

    if unpackPath == nothing
        # cleanup=true leads to issues with automatic testing on linux server.
        unpackPath = mktempdir(; prefix="fmijl_", cleanup=false)
    end

    zipPath = joinpath(unpackPath, fileName * ".zip")
    unzippedPath = joinpath(unpackPath, fileName)

    # only copy ZIP if not already there
    if !isfile(zipPath)
        cp(pathToFMU, zipPath; force=true)
    end

    @assert isfile(zipPath) ["fmi2Unzip(...): ZIP-Archive couldn't be copied to `$zipPath`."]

    zipAbsPath = isabspath(zipPath) ?  zipPath : joinpath(pwd(), zipPath)
    unzippedAbsPath = isabspath(unzippedPath) ? unzippedPath : joinpath(pwd(), unzippedPath)

    @assert isfile(zipAbsPath) ["fmi2Unzip(...): Can't deploy ZIP-Archive at `$(zipAbsPath)`."]

    numFiles = 0

    # only unzip if not already done
    if !isdir(unzippedAbsPath)
        mkpath(unzippedAbsPath)

        zarchive = ZipFile.Reader(zipAbsPath)
        for f in zarchive.files
            fileAbsPath = normpath(joinpath(unzippedAbsPath, f.name))

            if endswith(f.name,"/") || endswith(f.name,"\\")
                mkpath(fileAbsPath) # mkdir(fileAbsPath)

                @assert isdir(fileAbsPath) ["fmi2Unzip(...): Can't create directory `$(f.name)` at `$(fileAbsPath)`."]
            else
                # create directory if not forced by zip file folder
                mkpath(dirname(fileAbsPath))

                numBytes = write(fileAbsPath, read(f))

                if numBytes == 0
                    @info "fmi2Unzip(...): Written file `$(f.name)`, but file is empty."
                end

                @assert isfile(fileAbsPath) ["fmi2Unzip(...): Can't unzip file `$(f.name)` at `$(fileAbsPath)`."]
                numFiles += 1
            end
        end
        close(zarchive)
    end

    @assert isdir(unzippedAbsPath) ["fmi2Unzip(...): ZIP-Archive couldn't be unzipped at `$(unzippedPath)`."]
    @info "fmi2Unzip(...): Successfully unzipped $numFiles files at `$unzippedAbsPath`."

    (unzippedAbsPath, zipAbsPath)
end

# Checks with dlsym for available function in library.
# Prints an info text and returns C_NULL if not (soft-check).
# TODO used in FMI3_ext.jl too other spot to put it?
function dlsym_opt(libHandle, symbol)
    addr = dlsym(libHandle, symbol; throw_error=false)
    if addr == nothing
        @info "This FMU does not support function '$symbol'."
        addr = Ptr{Cvoid}(C_NULL)
    end
    addr
end

"""
Sets the properties of the fmu by reading the modelDescription.xml.
Retrieves all the pointers of binary functions.

Returns the instance of the FMU struct.

Via optional argument ```unpackPath```, a path to unpack the FMU can be specified (default: system temporary directory).
"""
function fmi2Load(pathToFMU::String; unpackPath=nothing, type=nothing)
    # Create uninitialized FMU
    fmu = FMU2()

    if startswith(pathToFMU, "http")
        @info "Downloading FMU from `$(pathToFMU)`."
        pathToFMU = Downloads.download(pathToFMU)
    end

    pathToFMU = normpath(pathToFMU)

    # set paths for fmu handling
    (fmu.path, fmu.zipPath) = fmi2Unzip(pathToFMU; unpackPath=unpackPath)

    # set paths for modelExchangeScripting and binary
    tmpName = splitpath(fmu.path)
    pathToModelDescription = joinpath(fmu.path, "modelDescription.xml")

    # parse modelDescription.xml
    fmu.modelDescription = fmi2LoadModelDescription(pathToModelDescription)
    fmu.modelName = fmu.modelDescription.modelName
    fmu.instanceName = fmu.modelDescription.modelName

    if (fmi2IsCoSimulation(fmu.modelDescription) && fmi2IsModelExchange(fmu.modelDescription) && type==:CS)
        fmu.type = fmi2TypeCoSimulation::fmi2Type
    elseif (fmi2IsCoSimulation(fmu.modelDescription) && fmi2IsModelExchange(fmu.modelDescription) && type==:ME)
        fmu.type = fmi2TypeModelExchange::fmi2Type
    elseif fmi2IsCoSimulation(fmu.modelDescription) && (type===nothing || type==:CS)
        fmu.type = fmi2TypeCoSimulation::fmi2Type
    elseif fmi2IsModelExchange(fmu.modelDescription) && (type===nothing || type==:ME)
        fmu.type = fmi2TypeModelExchange::fmi2Type
    else
        error(unknownFMUType)
    end

    fmuName = fmi2GetModelIdentifier(fmu.modelDescription; type=fmu.type) # tmpName[length(tmpName)]

    directoryBinary = ""
    pathToBinary = ""

    directories = []

    fmuExt = ""
    osStr = ""

    juliaArch = Sys.WORD_SIZE
    @assert (juliaArch == 64 || juliaArch == 32) "fmi2Load(...): Unknown Julia Architecture with $(juliaArch)-bit, must be 64- or 32-bit."

    if Sys.iswindows()
        if juliaArch == 64
            directories = [joinpath("binaries", "win64"), joinpath("binaries","x86_64-windows")]
        else
            directories = [joinpath("binaries", "win32"), joinpath("binaries","i686-windows")]
        end
        osStr = "Windows"
        fmuExt = "dll"
    elseif Sys.islinux()
        if juliaArch == 64
            directories = [joinpath("binaries", "linux64"), joinpath("binaries", "x86_64-linux")]
        else
            directories = []
        end
        osStr = "Linux"
        fmuExt = "so"
    elseif Sys.isapple()
        if juliaArch == 64
            directories = [joinpath("binaries", "darwin64"), joinpath("binaries", "x86_64-darwin")]
        else
            directories = []
        end
        osStr = "Mac"
        fmuExt = "dylib"
    else
        @assert false "fmi2Load(...): Unsupported target platform. Supporting Windows, Linux and Mac. Please open an issue if you want to use another OS."
    end

    @assert (length(directories) > 0) "fmi2Load(...): Unsupported architecture. Supporting Julia for Windows (64- and 32-bit), Linux (64-bit) and Mac (64-bit). Please open an issue if you want to use another architecture."
    for directory in directories
        directoryBinary = joinpath(fmu.path, directory)
        if isdir(directoryBinary)
            pathToBinary = joinpath(directoryBinary, "$(fmuName).$(fmuExt)")
            break
        end
    end
    @assert isfile(pathToBinary) "fmi2Load(...): Target platform is $(osStr), but can't find valid FMU binary at `$(pathToBinary)` for path `$(fmu.path)`."

    # make URI ressource location
    tmpResourceLocation = string("file:///", fmu.path)
    tmpResourceLocation = joinpath(tmpResourceLocation, "resources")
    fmu.fmuResourceLocation = replace(tmpResourceLocation, "\\" => "/") # URIs.escapeuri(tmpResourceLocation)

    @info "fmi2Load(...): FMU resources location is `$(fmu.fmuResourceLocation)`"

    if fmi2IsCoSimulation(fmu.modelDescription) && fmi2IsModelExchange(fmu.modelDescription)
        @info "fmi2Load(...): FMU supports both CS and ME, using CS as default if nothing specified."
    end

    fmu.binaryPath = pathToBinary
    loadBinary(fmu)

    # dependency matrix
    # fmu.dependencies

    fmu
end

function loadBinary(fmu::FMU2)
    lastDirectory = pwd()
    cd(dirname(fmu.binaryPath))

    # set FMU binary handler
    fmu.libHandle = dlopen(fmu.binaryPath) # , RTLD_NOW|RTLD_GLOBAL

    cd(lastDirectory)

    # retrieve functions
    fmu.cInstantiate                  = dlsym(fmu.libHandle, :fmi2Instantiate)
    fmu.cGetTypesPlatform             = dlsym(fmu.libHandle, :fmi2GetTypesPlatform)
    fmu.cGetVersion                   = dlsym(fmu.libHandle, :fmi2GetVersion)
    fmu.cFreeInstance                 = dlsym(fmu.libHandle, :fmi2FreeInstance)
    fmu.cSetDebugLogging              = dlsym(fmu.libHandle, :fmi2SetDebugLogging)
    fmu.cSetupExperiment              = dlsym(fmu.libHandle, :fmi2SetupExperiment)
    fmu.cEnterInitializationMode      = dlsym(fmu.libHandle, :fmi2EnterInitializationMode)
    fmu.cExitInitializationMode       = dlsym(fmu.libHandle, :fmi2ExitInitializationMode)
    fmu.cTerminate                    = dlsym(fmu.libHandle, :fmi2Terminate)
    fmu.cReset                        = dlsym(fmu.libHandle, :fmi2Reset)
    fmu.cGetReal                      = dlsym(fmu.libHandle, :fmi2GetReal)
    fmu.cSetReal                      = dlsym(fmu.libHandle, :fmi2SetReal)
    fmu.cGetInteger                   = dlsym(fmu.libHandle, :fmi2GetInteger)
    fmu.cSetInteger                   = dlsym(fmu.libHandle, :fmi2SetInteger)
    fmu.cGetBoolean                   = dlsym(fmu.libHandle, :fmi2GetBoolean)
    fmu.cSetBoolean                   = dlsym(fmu.libHandle, :fmi2SetBoolean)

    fmu.cGetString                    = dlsym_opt(fmu.libHandle, :fmi2GetString)
    fmu.cSetString                    = dlsym_opt(fmu.libHandle, :fmi2SetString)

    if fmi2CanGetSetState(fmu.modelDescription)
        fmu.cGetFMUstate                  = dlsym_opt(fmu.libHandle, :fmi2GetFMUstate)
        fmu.cSetFMUstate                  = dlsym_opt(fmu.libHandle, :fmi2SetFMUstate)
        fmu.cFreeFMUstate                 = dlsym_opt(fmu.libHandle, :fmi2FreeFMUstate)
    end

    if fmi2CanSerializeFMUstate(fmu.modelDescription)
        fmu.cSerializedFMUstateSize       = dlsym_opt(fmu.libHandle, :fmi2SerializedFMUstateSize)
        fmu.cSerializeFMUstate            = dlsym_opt(fmu.libHandle, :fmi2SerializeFMUstate)
        fmu.cDeSerializeFMUstate          = dlsym_opt(fmu.libHandle, :fmi2DeSerializeFMUstate)
    end

    if fmi2ProvidesDirectionalDerivative(fmu.modelDescription)
        fmu.cGetDirectionalDerivative     = dlsym_opt(fmu.libHandle, :fmi2GetDirectionalDerivative)
    end

    # CS specific function calls
    if fmi2IsCoSimulation(fmu.modelDescription)
        fmu.cSetRealInputDerivatives      = dlsym(fmu.libHandle, :fmi2SetRealInputDerivatives)
        fmu.cGetRealOutputDerivatives     = dlsym(fmu.libHandle, :fmi2GetRealOutputDerivatives)
        fmu.cDoStep                       = dlsym(fmu.libHandle, :fmi2DoStep)
        fmu.cCancelStep                   = dlsym(fmu.libHandle, :fmi2CancelStep)
        fmu.cGetStatus                    = dlsym(fmu.libHandle, :fmi2GetStatus)
        fmu.cGetRealStatus                = dlsym(fmu.libHandle, :fmi2GetRealStatus)
        fmu.cGetIntegerStatus             = dlsym(fmu.libHandle, :fmi2GetIntegerStatus)
        fmu.cGetBooleanStatus             = dlsym(fmu.libHandle, :fmi2GetBooleanStatus)
        fmu.cGetStringStatus              = dlsym(fmu.libHandle, :fmi2GetStringStatus)
    end

    # ME specific function calls
    if fmi2IsModelExchange(fmu.modelDescription)
        fmu.cEnterContinuousTimeMode      = dlsym(fmu.libHandle, :fmi2EnterContinuousTimeMode)
        fmu.cGetContinuousStates          = dlsym(fmu.libHandle, :fmi2GetContinuousStates)
        fmu.cGetDerivatives               = dlsym(fmu.libHandle, :fmi2GetDerivatives)
        fmu.cSetTime                      = dlsym(fmu.libHandle, :fmi2SetTime)
        fmu.cSetContinuousStates          = dlsym(fmu.libHandle, :fmi2SetContinuousStates)
        fmu.cCompletedIntegratorStep      = dlsym(fmu.libHandle, :fmi2CompletedIntegratorStep)
        fmu.cEnterEventMode               = dlsym(fmu.libHandle, :fmi2EnterEventMode)
        fmu.cNewDiscreteStates            = dlsym(fmu.libHandle, :fmi2NewDiscreteStates)
        fmu.cGetEventIndicators           = dlsym(fmu.libHandle, :fmi2GetEventIndicators)
        fmu.cGetNominalsOfContinuousStates= dlsym(fmu.libHandle, :fmi2GetNominalsOfContinuousStates)
    end
end

"""
TODO: FMI specification reference.

Create a new instance of the given fmu, adds a logger if logginOn == true.

Returns the instance of a new FMU component.

For more information call ?fmi2Instantiate

# Keywords
- `visible` if the FMU should be started with graphic interface, if supported (default=`false`)
- `loggingOn` if the FMU should log and display function calls (default=`false`)
- `externalCallbacks` if an external DLL should be used for the fmi2CallbackFunctions, this may improve readability of logging messages (default=`false`)
- `logStatusOK` whether to log status of kind `fmi2OK` (default=`true`)
- `logStatusWarning whether to log status of kind `fmi2Warning` (default=`true`)
- `logStatusDiscard whether to log status of kind `fmi2Discard` (default=`true`)
- `logStatusError whether to log status of kind `fmi2Error` (default=`true`)
- `logStatusFatal whether to log status of kind `fmi2Fatal` (default=`true`)
- `logStatusPending whether to log status of kind `fmi2Pending` (default=`true`)
"""
function fmi2Instantiate!(fmu::FMU2; pushComponents::Bool = true, visible::Bool = false, loggingOn::Bool = fmu.executionConfig.loggingOn, externalCallbacks::Bool = fmu.executionConfig.externalCallbacks,
                          logStatusOK::Bool=true, logStatusWarning::Bool=true, logStatusDiscard::Bool=true, logStatusError::Bool=true, logStatusFatal::Bool=true, logStatusPending::Bool=true,
                          instanceName::String=fmu.modelName)

    compEnv = FMU2ComponentEnvironment()
    compEnv.logStatusOK = logStatusOK
    compEnv.logStatusWarning = logStatusWarning
    compEnv.logStatusDiscard = logStatusDiscard
    compEnv.logStatusError = logStatusError
    compEnv.logStatusFatal = logStatusFatal
    compEnv.logStatusPending = logStatusPending

    ptrLogger = @cfunction(fmi2CallbackLogger, Cvoid, (Ptr{FMU2ComponentEnvironment}, Ptr{Cchar}, Cuint, Ptr{Cchar}, Ptr{Cchar}))
    if externalCallbacks
        if fmu.callbackLibHandle == C_NULL
            @assert Sys.iswindows() && Sys.WORD_SIZE == 64 "`externalCallbacks=true` is only supported for Windows 64-bit."

            cbLibPath = joinpath(dirname(@__FILE__), "callbackFunctions", "binaries", "win64", "callbackFunctions.dll")

            # check permission to execute the DLL
            perm = filemode(cbLibPath)
            permRWX = 16895
            if perm != permRWX
                chmod(cbLibPath, permRWX; recursive=true)
            end

            fmu.callbackLibHandle = dlopen(cbLibPath)
        end
        ptrLogger = dlsym(fmu.callbackLibHandle, :logger)
    end
    ptrAllocateMemory = @cfunction(fmi2CallbackAllocateMemory, Ptr{Cvoid}, (Csize_t, Csize_t))
    ptrFreeMemory = @cfunction(fmi2CallbackFreeMemory, Cvoid, (Ptr{Cvoid},))
    ptrStepFinished = C_NULL # ToDo
    ptrComponentEnvironment = Ptr{FMU2ComponentEnvironment}(pointer_from_objref(compEnv))
    callbackFunctions = fmi2CallbackFunctions(ptrLogger, ptrAllocateMemory, ptrFreeMemory, ptrStepFinished, ptrComponentEnvironment)

    guidStr = "$(fmu.modelDescription.guid)"

    compAddr = fmi2Instantiate(fmu.cInstantiate, pointer(instanceName), fmu.type, pointer(guidStr), pointer(fmu.fmuResourceLocation), Ptr{fmi2CallbackFunctions}(pointer_from_objref(callbackFunctions)), fmi2Boolean(visible), fmi2Boolean(loggingOn))

    if compAddr == Ptr{Cvoid}(C_NULL)
        @error "fmi2Instantiate!(...): Instantiation failed!"
        return nothing
    end

    component = nothing

    # check if address is already inside of the components (this may be in FMIExport.jl)
    for c in fmu.components
        if c.compAddr == compAddr
            component = c
            break
        end
    end

    if component != nothing
        @info "fmi2Instantiate!(...): This component was already registered. This may be because you created the FMU by yourself with FMIExport.jl."
    else
        component = FMU2Component(compAddr, fmu)
        component.jacobianUpdate! = fmi2GetJacobian!
        component.componentEnvironment = compEnv
        component.callbackFunctions = callbackFunctions

        if pushComponents
            push!(fmu.components, component)
        end
    end

    component
end

"""
Reloads the FMU-binary. This is useful, if the FMU does not support a clean reset implementation.
"""
function fmi2Reload(fmu::FMU2)
    dlclose(fmu.libHandle)
    loadBinary(fmu)
end

"""
Unload a FMU.

Free the allocated memory, close the binaries and remove temporary zip and unziped FMU model description.
"""
function fmi2Unload(fmu::FMU2, cleanUp::Bool = true)

    while length(fmu.components) > 0
        fmi2FreeInstance!(fmu.components[end])
    end

    dlclose(fmu.libHandle)

    # the components are removed from the component list via call to fmi2FreeInstance!
    @assert length(fmu.components) == 0 "fmi2Unload(...): Failure during deleting components, $(length(fmu.components)) remaining in stack."

    if cleanUp
        try
            rm(fmu.path; recursive = true, force = true)
            rm(fmu.zipPath; recursive = true, force = true)
        catch e
            @warn "Cannot delete unpacked data on disc. Maybe some files are opened in another application."
        end
    end
end

"""
This function samples the directional derivative by manipulating corresponding values (central differences).
"""
function fmi2SampleDirectionalDerivative(c::FMU2Component,
                                       vUnknown_ref::AbstractArray{fmi2ValueReference},
                                       vKnown_ref::AbstractArray{fmi2ValueReference},
                                       steps::Union{AbstractArray{fmi2Real}, Nothing} = nothing)

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

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

    dvUnknown
end

"""
This function samples the directional derivative by manipulating corresponding values (central differences) and saves in-place.
"""
function fmi2SampleDirectionalDerivative!(c::FMU2Component,
                                          vUnknown_ref::AbstractArray{fmi2ValueReference},
                                          vKnown_ref::AbstractArray{fmi2ValueReference},
                                          dvUnknown::AbstractArray, # ToDo: datatype
                                          steps::Union{AbstractArray{fmi2Real}, Nothing} = nothing)

    step = 0.0

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

        if steps === nothing
            # smaller than 1e-6 leads to issues
            step = max(2.0 * eps(Float32(origValue)), 1e-6)
        else
            step = steps[i]
        end

        fmi2SetReal(c, vKnown, origValue - step)
        negValues = fmi2GetReal(c, vUnknown_ref)

        fmi2SetReal(c, vKnown, origValue + step)
        posValues = fmi2GetReal(c, vUnknown_ref)

        fmi2SetReal(c, vKnown, origValue)

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

    nothing
end

"""
Builds the jacobian over the FMU `fmu` for FMU value references `rdx` and `rx`, so that the function returns the jacobian ∂rdx / ∂rx.

If FMI built-in directional derivatives are supported, they are used.
As fallback, directional derivatives will be sampled with central differences.
For optimization, if the FMU's model description has the optional entry 'dependencies', only dependent variables are sampled/retrieved. This drastically boosts performance for systems with large variable count (like CFD).

If sampling is used, sampling step size can be set (for each direction individually) using optional argument `steps`.
"""
function fmi2GetJacobian(comp::FMU2Component,
                         rdx::AbstractArray{fmi2ValueReference},
                         rx::AbstractArray{fmi2ValueReference};
                         steps::Union{AbstractArray{fmi2Real}, Nothing} = nothing)
    mat = zeros(fmi2Real, length(rdx), length(rx))
    fmi2GetJacobian!(mat, comp, rdx, rx; steps=steps)
    return mat
end

"""
Fills the jacobian over the FMU `fmu` for FMU value references `rdx` and `rx`, so that the function returns the jacobian ∂rdx / ∂rx.

If FMI built-in directional derivatives are supported, they are used.
As fallback, directional derivatives will be sampled with central differences.
For optimization, if the FMU's model description has the optional entry 'dependencies', only dependent variables are sampled/retrieved. This drastically boosts performance for systems with large variable count (like CFD).

If sampling is used, sampling step size can be set (for each direction individually) using optional argument `steps`.
"""
function fmi2GetJacobian!(jac::AbstractMatrix{fmi2Real},
                          comp::FMU2Component,
                          rdx::AbstractArray{fmi2ValueReference},
                          rx::AbstractArray{fmi2ValueReference};
                          steps::Union{AbstractArray{fmi2Real}, Nothing} = nothing)

    @assert size(jac) == (length(rdx), length(rx)) ["fmi2GetJacobian!: Dimension missmatch between `jac` $(size(jac)), `rdx` ($length(rdx)) and `rx` ($length(rx))."]

    if length(rdx) == 0 || length(rx) == 0
        jac = zeros(length(rdx), length(rx))
        return nothing
    end

    ddsupported = fmi2ProvidesDirectionalDerivative(comp.fmu)

    # ToDo: Pick entries based on dependency matrix!
    #depMtx = fmi2GetDependencies(fmu)
    rdx_inds = collect(comp.fmu.modelDescription.valueReferenceIndicies[vr] for vr in rdx)
    rx_inds  = collect(comp.fmu.modelDescription.valueReferenceIndicies[vr] for vr in rx)

    for i in 1:length(rx)

        sensitive_rdx_inds = 1:length(rdx)
        sensitive_rdx = rdx

        # sensitive_rdx_inds = Int64[]
        # sensitive_rdx = fmi2ValueReference[]

        # for j in 1:length(rdx)
        #     if depMtx[rdx_inds[j], rx_inds[i]] != fmi2DependencyIndependent
        #         push!(sensitive_rdx_inds, j)
        #         push!(sensitive_rdx, rdx[j])
        #     end
        # end

        if length(sensitive_rdx) > 0
            if ddsupported
                # doesn't work because indexed-views can`t be passed by reference (to ccalls)
                #try
                    fmi2GetDirectionalDerivative!(comp, sensitive_rdx, [rx[i]], view(jac, sensitive_rdx_inds, i))
                #catch e
                #    jac[sensitive_rdx_inds, i] = fmi2GetDirectionalDerivative(comp, sensitive_rdx, [rx[i]])
                #end
            else
                # doesn't work because indexed-views can`t be passed by reference (to ccalls)
                #try
                    fmi2SampleDirectionalDerivative!(comp, sensitive_rdx, [rx[i]], view(jac, sensitive_rdx_inds, i))
                #catch e
                #    jac[sensitive_rdx_inds, i] = fmi2SampleDirectionalDerivative(comp, sensitive_rdx, [rx[i]], steps)
                #end
            end
        end
    end

    return nothing
end

"""
Builds the jacobian over the FMU `fmu` for FMU value references `rdx` and `rx`, so that the function returns the jacobian ∂rdx / ∂rx.

If FMI built-in directional derivatives are supported, they are used.
As fallback, directional derivatives will be sampled with central differences.
No performance optimization, for an optimized version use `fmi2GetJacobian`.

If sampling is used, sampling step size can be set (for each direction individually) using optional argument `steps`.
"""
function fmi2GetFullJacobian(comp::FMU2Component,
                             rdx::AbstractArray{fmi2ValueReference},
                             rx::AbstractArray{fmi2ValueReference};
                             steps::Union{AbstractArray{fmi2Real}, Nothing} = nothing)
    mat = zeros(fmi2Real, length(rdx), length(rx))
    fmi2GetFullJacobian!(mat, comp, rdx, rx; steps=steps)
    return mat
end

"""
Fills the jacobian over the FMU `fmu` for FMU value references `rdx` and `rx`, so that the function returns the jacobian ∂rdx / ∂rx.

If FMI built-in directional derivatives are supported, they are used.
As fallback, directional derivatives will be sampled with central differences.
No performance optimization, for an optimized version use `fmi2GetJacobian!`.

If sampling is used, sampling step size can be set (for each direction individually) using optional argument `steps`.
"""
function fmi2GetFullJacobian!(jac::AbstractMatrix{fmi2Real},
                              comp::FMU2Component,
                              rdx::AbstractArray{fmi2ValueReference},
                              rx::AbstractArray{fmi2ValueReference};
                              steps::Union{AbstractArray{fmi2Real}, Nothing} = nothing)
    @assert size(jac) == (length(rdx),length(rx)) "fmi2GetFullJacobian!: Dimension missmatch between `jac` $(size(jac)), `rdx` ($length(rdx)) and `rx` ($length(rx))."

    @warn "`fmi2GetFullJacobian!` is for benchmarking only, please use `fmi2GetJacobian`."

    if length(rdx) == 0 || length(rx) == 0
        jac = zeros(length(rdx), length(rx))
        return nothing
    end

    if fmi2ProvidesDirectionalDerivative(comp.fmu)
        for i in 1:length(rx)
            jac[:,i] = fmi2GetDirectionalDerivative(comp, rdx, [rx[i]])
        end
    else
        jac = fmi2SampleDirectionalDerivative(comp, rdx, rx)
    end

    return nothing
end

function fmi2Get!(comp::FMU2Component, vrs::fmi2ValueReferenceFormat, dstArray::AbstractArray)
    vrs = prepareValueReference(comp, vrs)

    @assert length(vrs) == length(dstArray) "fmi2Get!(...): Number of value references doesn't match number of `dstArray` elements."

    retcodes = zeros(fmi2Status, length(vrs)) # fmi2StatusOK

    for i in 1:length(vrs)
        vr = vrs[i]
        mv = fmi2ModelVariablesForValueReference(comp.fmu.modelDescription, vr)
        mv = mv[1]

        if mv._Real != nothing
            #@assert isa(dstArray[i], Real) "fmi2Get!(...): Unknown data type for value reference `$(vr)` at index $(i), should be `Real`, is `$(typeof(dstArray[i]))`."
            dstArray[i] = fmi2GetReal(comp, vr)
        elseif mv._Integer != nothing
            #@assert isa(dstArray[i], Union{Real, Integer}) "fmi2Get!(...): Unknown data type for value reference `$(vr)` at index $(i), should be `Integer`, is `$(typeof(dstArray[i]))`."
            dstArray[i] = fmi2GetInteger(comp, vr)
        elseif mv._Boolean != nothing
            #@assert isa(dstArray[i], Union{Real, Bool}) "fmi2Get!(...): Unknown data type for value reference `$(vr)` at index $(i), should be `Bool`, is `$(typeof(dstArray[i]))`."
            dstArray[i] = fmi2GetBoolean(comp, vr)
        elseif mv._String != nothing
            #@assert isa(dstArray[i], String) "fmi2Get!(...): Unknown data type for value reference `$(vr)` at index $(i), should be `String`, is `$(typeof(dstArray[i]))`."
            dstArray[i] = fmi2GetString(comp, vr)
        elseif mv._Enumeration != nothing
            @warn "fmi2Get!(...): Currently not implemented for fmi2Enum."
        else
            @assert isa(dstArray[i], Real) "fmi2Get!(...): Unknown data type for value reference `$(vr)` at index $(i), is `$(mv.datatype.datatype)`."
        end
    end

    return retcodes
end

function fmi2Get(comp::FMU2Component, vrs::fmi2ValueReferenceFormat)
    vrs = prepareValueReference(comp, vrs)
    dstArray = Array{Any,1}(undef, length(vrs))
    fmi2Get!(comp, vrs, dstArray)
    return dstArray
end

function fmi2Set(comp::FMU2Component, vrs::fmi2ValueReferenceFormat, srcArray::AbstractArray)
    vrs = prepareValueReference(comp, vrs)

    @assert length(vrs) == length(srcArray) "fmi2Set(...): Number of value references doesn't match number of `srcArray` elements."

    retcodes = zeros(fmi2Status, length(vrs)) # fmi2StatusOK

    for i in 1:length(vrs)
        vr = vrs[i]
        mv = fmi2ModelVariablesForValueReference(comp.fmu.modelDescription, vr)
        mv = mv[1]

        if mv._Real != nothing
            @assert isa(srcArray[i], Real) "fmi2Set(...): Unknown data type for value reference `$(vr)` at index $(i), should be `Real`, is `$(typeof(srcArray[i]))`."
            retcodes[i] = fmi2SetReal(comp, vr, srcArray[i])
        elseif mv._Integer != nothing
            @assert isa(srcArray[i], Union{Real, Integer}) "fmi2Set(...): Unknown data type for value reference `$(vr)` at index $(i), should be `Integer`, is `$(typeof(srcArray[i]))`."
            retcodes[i] = fmi2SetInteger(comp, vr, Integer(srcArray[i]))
        elseif mv._Boolean != nothing
            @assert isa(srcArray[i], Union{Real, Bool}) "fmi2Set(...): Unknown data type for value reference `$(vr)` at index $(i), should be `Bool`, is `$(typeof(srcArray[i]))`."
            retcodes[i] = fmi2SetBoolean(comp, vr, Bool(srcArray[i]))
        elseif mv._String != nothing
            @assert isa(srcArray[i], String) "fmi2Set(...): Unknown data type for value reference `$(vr)` at index $(i), should be `String`, is `$(typeof(srcArray[i]))`."
            retcodes[i] = fmi2SetString(comp, vr, srcArray[i])
        elseif mv._Enumeration != nothing
            @warn "fmi2Set(...): Currently not implemented for fmi2Enum."
        else
            @assert false "fmi2Set(...): Unknown data type for value reference `$(vr)` at index $(i), is `$(mv.datatype.datatype)`."
        end
    end

    return retcodes
end

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

TODO: Add this command in the documentation.
"""
function fmi2GetStartValue(md::fmi2ModelDescription, vrs::fmi2ValueReferenceFormat = md.valueReferences)

    vrs = prepareValueReference(md, vrs)

    starts = []

    for vr in vrs
        mvs = fmi2ModelVariablesForValueReference(md, vr)

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

        push!(starts, fmi2GetStartValue(mvs[1]) )
    end

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

function fmi2GetStartValue(fmu::FMU2, vrs::fmi2ValueReferenceFormat = fmu.modelDescription.valueReferences)
    fmi2GetStartValue(fmu.modelDescription, vrs)
end

function fmi2GetStartValue(c::FMU2Component, vrs::fmi2ValueReferenceFormat = c.fmu.modelDescription.valueReferences)
    fmi2GetStartValue(c.fmu, vrs)
end

function fmi2GetStartValue(mv::fmi2ScalarVariable)
    if mv._Real != nothing
        return mv._Real.start
    elseif mv._Integer != nothing
        return mv._Integer.start
    elseif mv._Boolean != nothing
        return mv._Boolean.start
    elseif mv._String != nothing
        return mv._String.start
    elseif mv._Enumeration != nothing
        return mv._Enumeration.start
    else
        @assert false "fmi2GetStartValue(...): Variable $(mv) has no data type."
    end
end

"""
Returns the `unit` entry of the corresponding model variable.

ToDo: update docstring format.
"""
function fmi2GetUnit(mv::fmi2ScalarVariable)
    if mv._Real != nothing
        return mv._Real.unit
    else
        return nothing
    end
end

"""
Returns the `inital` entry of the corresponding model variable.

ToDo: update docstring format.
"""
function fmi2GetInitial(mv::fmi2ScalarVariable)
    return mv.initial
end


text:
# ---
# title: 237. Delete Node in a Linked List
# id: problem237
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Linked List
# link: <https://leetcode.com/problems/delete-node-in-a-linked-list/description/>
# hidden: true
# ---
#
# Write a function to **delete a node** in a singly-linked list. You will
# **not** be given access to the `head` of the list, instead you will be given
# access to **the node to be deleted** directly.
#
# It is **guaranteed** that the node to be deleted is **not a tail node** in the
# list.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/01/node1.jpg)
#
#
#
#     Input: head = [4,5,1,9], node = 5
#     Output: [4,1,9]
#     Explanation: You are given the second node with value 5, the linked list should become 4 -> 1 -> 9 after calling your function.
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/01/node2.jpg)
#
#
#
#     Input: head = [4,5,1,9], node = 1
#     Output: [4,5,9]
#     Explanation: You are given the third node with value 1, the linked list should become 4 -> 5 -> 9 after calling your function.
#
#
# **Example 3:**
#
#
#
#     Input: head = [1,2,3,4], node = 3
#     Output: [1,2,4]
#
#
# **Example 4:**
#
#
#
#     Input: head = [0,1], node = 0
#     Output: [1]
#
#
# **Example 5:**
#
#
#
#     Input: head = [-3,5,-99], node = -3
#     Output: [5,-99]
#
#
#
#
# **Constraints:**
#
#   * The number of the nodes in the given list is in the range `[2, 1000]`.
#   * `-1000 <= Node.val <= 1000`
#   * The value of each node in the list is **unique**.
#   * The `node` to be deleted is **in the list** and is **not a tail** node
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
@rule AR(:y, Marginalisation) (m_x::NormalDistributionsFamily, q_θ::NormalDistributionsFamily, q_γ::GammaShapeRate, meta::ARMeta) = begin
    mθ, Vθ = mean_cov(q_θ)
    mx, Wx = mean_invcov(m_x)

    mγ = mean(q_γ)

    mA = as_companion_matrix(mθ)
    mV = ar_transition(getvform(meta), getorder(meta), mγ)

    D = Wx + mγ * Vθ
    C = mA * inv(D)

    my = C * Wx * mx
    Vy = add_transition!(C * mA', mV)

    return convert(promote_variate_type(getvform(meta), NormalMeanVariance), my, Vy)
end


text:
function random_dist(n::Int64)
    point = [0.0, 0.0]
    f = open("random.dat", "w")

    for i = 1:n
        point = [rand(),rand()]
        println(f, point[1], '\t', point[2],"\n\n")
    end

    close(f)
end

function barnsley_tree(n::Int64)
    point = [0.0, 0.0]
    f = open("barnsley_tree.dat", "w")
    for i = 1:n
        rnd = rand()
        if (rnd <= 0.02)
            point = [0.03 0; 0 0.1]*point
        elseif(rnd > 0.02 && rnd <= 0.62)
            point = [0.85 0.0; 0.0 0.85]*point + [0, 1.50]
        elseif(rnd > 0.62 && rnd <= 0.72)
            point = [0.8 0; 0 0.8]*point + [0, 1.50]
        elseif(rnd > 0.72 && rnd <= 0.79)
            point = [0.2 -0.08; 0.15 0.22]*point + [0, 0.85]
        elseif(rnd > 0.79 && rnd <= 0.86)
            point = [-0.2 0.08; 0.15 0.22]*point + [0, 0.85]
        elseif(rnd > 0.86 && rnd <= 0.93)
            point = [0.25 -0.1; 0.12 0.25]*point + [0, 0.3]
        else
            point = [-0.2 0.1; 0.12 0.2]*point + [0, 0.4]
        end

        println(f, point[1], '\t', point[2])

        if (i%10000==0)
            println(f,"\n\n")
        end
    end
    close(f)
end

function barnsley_fern(n::Int64)
    point = [0.0, 0.0]
    f = open("barnsley.dat", "w")
    for i = 1:n
        rnd = rand()
        if (rnd <= 0.01)
            point = [0 0; 0 0.16]*point
        elseif(rnd > 0.01 && rnd <= 0.86)
            point = [0.85 0.04; -0.04 0.85]*point + [0, 1.60]
        elseif(rnd > 0.86 && rnd <= 0.93)
            point = [0.2 -0.26; 0.23 0.22]*point + [0, 1.60]
        else
            point = [-0.15 0.28; 0.26 0.24]*point + [0, 0.44]
        end

        println(f, point[1], '\t', point[2])

        if (i%10000==0)
            println(f,"\n\n")
        end
    end
    close(f)
end

function sierpensky(n::Int64)
    element_x = rand()
    element_y = rand()
    point = [element_x element_y]

    f = open("sierpensky.dat", "w")

    for i = 1:n
        if (i > 20)
            println(f, point[1], '\t', point[2])
        end
        rnd = rand()

        if (rnd <= 0.33)
            point = 0.5*point
        elseif(rnd > 0.33 && rnd <= 0.66)
            point[1] = (point[1]+0.5)/2
            point[2] = (point[2]+1)/2
        else
            point[1] = (point[1]+1)/2
            point[2] = (point[2])/2
        end
    end

    close(f)
end

#random_dist(1000)
#barnsley_fern(1000000)
barnsley_tree(1000000)
#sierpensky(10000)


text:
@testset "Compressed array string" begin
    @test NaiveNASlib.compressed_string([1,2,3,4]) == "[1, 2, 3, 4]"
    @test NaiveNASlib.compressed_string(1:23) == "[1,…, 23]"
    @test NaiveNASlib.compressed_string([1,2,3,5,6,8,9,10:35...]) == "[1, 2, 3, 5, 6, 8,…, 35]"
    @test NaiveNASlib.compressed_string(-ones(Int, 24)) == "[-1×24]"
    @test NaiveNASlib.compressed_string([-1, 5:9...,-ones(Int, 6)...,23,25,99,100,102,23:32...,34]) ==  "[-1, 5,…, 9, -1×6, 23, 25, 99, 100, 102, 23,…, 32, 34]"
    @test NaiveNASlib.compressed_string([1,2,4:13..., 10:20...,-1]) == "[1, 2, 4,…, 13, 10,…, 20, -1]"
end


text:
using RobotDynamics
using StaticArrays
using BenchmarkTools
const RD = RobotDynamics

nonlinear_model = Cartpole()

N = 5
dt = 0.01
n = state_dim(nonlinear_model)
m = control_dim(nonlinear_model)

Z = RD.Traj(n, m, dt, N)
times = RD.get_times(Z)

for knot_point in Z
    knot_point.z = @SVector rand(n+m)
end

linearized_model = LinearizedModel(nonlinear_model, Z, dt=dt, integration=RK3, is_affine=true)

## Test that the linearization matches exactly at the knot points
for k=1:N-1
    local z = Z[k]
    @test discrete_dynamics(RK3, nonlinear_model, z) ≈
        discrete_dynamics(PassThrough, linearized_model, z) atol=1e-4

    local F1 = RD.DynamicsJacobian(nonlinear_model)
    local F2 = RD.DynamicsJacobian(linearized_model)
    discrete_jacobian!(RK3, F1, nonlinear_model, z)
    discrete_jacobian!(PassThrough, F2, linearized_model, z)
    @test F1 ≈ F2
end

## Test with updated trajectory
for z in Z
    z.z = @SVector rand(n+m)
end

update_trajectory!(linearized_model, Z)

for i=1:N-1
    local z = Z[i]
    @test discrete_dynamics(RK3, nonlinear_model, z) ≈
        discrete_dynamics(PassThrough, linearized_model, z) atol=1e-4

    local F1 = RD.DynamicsJacobian(nonlinear_model)
    local F2 = RD.DynamicsJacobian(linearized_model)
    discrete_jacobian!(RK3, F1, nonlinear_model, z)
    discrete_jacobian!(PassThrough, F2, linearized_model, z)
    @test F1 ≈ F2
end

@test (@ballocated update_trajectory!($linearized_model, $Z) samples=2 evals=2) == 0
@test RD.integration(linearized_model) == RK3

## LinearModel
n,m = 10,5
A,B = gencontrollable(n,m)
d = rand(n)
model = LinearModel(A,B)
model_affine = LinearModel(A,B,d)

x,u = rand(model)
dynamics(model_affine, x, u)

# continuous -> continuous
linmodel = LinearizedModel(model)
@test RD.is_discrete(linmodel) == false
@test RD.get_A(linmodel) ≈ A
@test RD.get_B(linmodel) ≈ B

# continuous -> discrete
dt = 0.01
linmodel = LinearizedModel(model, integration=RK2, dt=dt)
@test linmodel.Z[1].dt == dt
@test RD.is_discrete(linmodel)
@test RD.is_timevarying(linmodel) == false
@test RD.is_affine(linmodel) == false

@test RD.get_A(linmodel) ≈ I + A'*(I + 0.5*A*dt)*dt
@test RD.get_B(linmodel) ≈ (0.5*A*dt * B + B)*dt

# discrete -> continuous
model = LinearModel(A, B, dt=dt)
@test RD.is_discrete(model)
@test_throws AssertionError LinearizedModel(model)

# discrete -> discrete
linmodel = LinearizedModel(model, dt=dt, integration=PassThrough)
@test RD.get_A(linmodel) == A
@test RD.get_B(linmodel) == B

## Exponential
model = LinearModel(A, B)
@test RD.is_discrete(model) == false

linmodel = LinearizedModel(model, integration=Exponential, dt=dt)
RD.get_A(linmodel)
Ec = [A B; zeros(m,n+m)]
Ed = exp(Ec*dt)
@test Ed[1:n,1:n] ≈ RD.get_A(linmodel)
@test Ed[1:n,n .+ (1:m)] ≈ RD.get_B(linmodel)

## Try a large system
# dt = 0.01
# n,m = 150,50
# A,B = gencontrollable(n,m)
# e = n+m
# E = SizedMatrix{e,e}(zeros(e,e))
# @time Ed = RD.matrix_exponential!(E, A, B, dt)
# model = LinearModel(A,B)
# @test model.use_static == false
# linmodel = LinearizedModel(model, dt=dt, integration=Exponential)

# @btime RD.linearize!($linmodel)

###################################
## Exponential Discretization Tests
###################################

# using ControlSystems
# using StaticArrays

# dt = 0.01
# C_111 = ss([-5], [2], [3], [0])
# D_111 = c2d(C_111, dt, :zoh)
# exp([-5 2; 0 0]*dt)

# C_111_rd = LinearModel(-5*ones(1,1), 2*ones(1,1))
# D_111_rd = LinearizedModel(C_111_rd, dt=dt, integration=Exponential)

# @test D_111[1].A ≈ RD.get_A(D_111_rd)
# @test D_111[1].B ≈ RD.get_B(D_111_rd)

# ####
# C_212 = ss([-5 -3; 2 -9], [1; 2], [1 0; 0 1], [0; 0])
# D_212 = c2d(C_212, dt, :zoh)
# C_212_rd = LinearModel([-5 -3; 2 -9], reshape([1; 2], 2, 1))
# D_212_rd = LinearizedModel(C_212_rd, dt=dt, integration=Exponential)

# @test D_212[1].A ≈ RD.get_A(D_212_rd)
# @test D_212[1].B ≈ RD.get_B(D_212_rd)

# ###
# C_221 = ss([-5 -3; 2 -9], [1 0; 0 2], [1 0], [0 0])
# D_221 = c2d(C_221, dt, :zoh)
# C_221_rd = LinearModel([-5 -3; 2 -9], [1 0; 0 2])
# D_221_rd = LinearizedModel(C_221_rd, dt=dt, integration=Exponential)

# @test D_221[1].A ≈ RD.get_A(D_221_rd)
# @test D_221[1].B ≈ RD.get_B(D_221_rd)

# Constant Discrete Jacobian Tests:
n,m = 5, 5
A, B = gencontinuous(n,m)
random_model = LinearModel(A, B)
random_model_discrete = LinearModel(n, m; dt=dt)
x,u = rand(random_model)
z = KnotPoint(x, u, dt)

struct RandomLinear{n,m,T} <: AbstractModel
    A::SMatrix{n,n,T}
    B::SMatrix{n,m,T}
end
RD.state_dim(::RandomLinear{n}) = n
RD.control_dim(::RandomLinear{<:Any, m}) = m
RD.dynamics(model::RandomLinear, x, u) = model.A*x + model.B*u

random_model_test = RandomLinear{n,m,Float64}(A,B)
F = RD.DynamicsJacobian(n,m)

for RK in (RK2, RK3, RK4)
    local linmodel = LinearizedModel(random_model_test, integration=RK, dt=dt)
    discrete_jacobian!(RK, F, random_model_test, z)
    @test RD.get_A(F) ≈ RD.get_A(linmodel)
    @test RD.get_B(F) ≈ RD.get_B(linmodel)
end


text:
export list_voice_regions

"""
    list_voice_regions(c::Client) -> Vector{VoiceRegion}

Get a list of the [`VoiceRegion`](@ref)s that can be used when creating [`Guild`](@ref)s.
"""
function list_voice_regions(c::Client)
    return Response{Vector{VoiceRegion}}(c, :GET, "/voice/regions")
end


text:
using Pkg # hideall
Pkg.activate("_literate/EX-telco/Project.toml")
Pkg.update()
macro OUTPUT()
    return isdefined(Main, :Franklin) ? Franklin.OUT_PATH[] : "/tmp/"
end;

# An application of the [MLJ
# toolbox](https://alan-turing-institute.github.io/MLJ.jl/dev/) to the
# Telco Customer Churn dataset, aimed at practicing data scientists
# new to MLJ (Machine Learning in Julia). This tutorial does not
# cover exploratory data analysis.

# MLJ is a *multi-paradigm* machine learning toolbox (i.e., not just
# deep-learning).

# For other MLJ learning resources see the [Learning
# MLJ](https://alan-turing-institute.github.io/MLJ.jl/dev/learning_mlj/)
# section of the
# [manual](https://alan-turing-institute.github.io/MLJ.jl/dev/).

# **Topics covered**: Grabbing and preparing a dataset, basic
# fit/predict workflow, constructing a pipeline to include data
# pre-processing, estimating performance metrics, ROC curves, confusion
# matrices, feature importance, basic feature selection, controlling iterative
# models, hyper-parameter optimization (tuning).

# **Prerequisites for this tutorial.** Previous experience building,
# evaluating, and optimizing machine learning models using
# scikit-learn, caret, MLR, weka, or similar tool. No previous
# experience with MLJ. Only fairly basic familiarity with Julia is
# required. Uses
# [DataFrames.jl](https://dataframes.juliadata.org/stable/) but in a
# minimal way ([this
# cheatsheet](https://ahsmart.com/pub/data-wrangling-with-data-frames-jl-cheat-sheet/index.html)
# may help).

# **Time.** Between two and three hours, first time through.

# ## Summary of methods and types introduced

# |code   | purpose|
# |:-------|:-------------------------------------------------------|
# | `OpenML.load(id)` | grab a dataset from [OpenML.org](https://www.openml.org)|
# | `scitype(X)`      | inspect the scientific type (scitype) of object `X`|
# | `schema(X)`       | inspect the column scitypes (scientific types) of a table `X`|
# | `coerce(X, ...)`   | fix column encodings to get appropriate scitypes|
# | `partition(data, frac1, frac2, ...; rng=...)` | vertically split `data`, which can be a table, vector or matrix|
# | `unpack(table, f1, f2, ...)` | horizontally split `table` based on conditions `f1`, `f2`, ..., applied to column names|
# | `@load ModelType pkg=...`           | load code defining a model type|
# | `input_scitype(model)` | inspect the scitype that a model requires for features (inputs)|
# | `target_scitype(model)`| inspect the scitype that a model requires for the target (labels)|
# | `ContinuousEncoder`   | built-in model type for re-encoding all features as `Continuous`|
# | `model1 ∣> model2 ∣> ...` | combine multiple models into a pipeline|
# | `measures("under curve")` | list all measures (metrics) with string "under curve" in documentation|
# | `accuracy(yhat, y)` | compute accuracy of predictions `yhat` against ground truth observations `y`|
# | `auc(yhat, y), brier_loss(yhat, y)` | evaluate two probabilistic measures (`yhat` a vector of probability distributions)|
# | `machine(model, X, y)` | bind `model` to training data `X` (features) and `y` (target)|
# | `fit!(mach, rows=...)` | train machine using specified rows (observation indices)|
# | `predict(mach, rows=...)`, | make in-sample model predictions given specified rows|
# | `predict(mach, Xnew)` | make predictions given new features `Xnew`|
# | `fitted_params(mach)` | inspect learned parameters|
# | `report(mach)`        | inspect other outcomes of training|
# | `confmat(yhat, y)`    | confusion matrix for predictions `yhat` and ground truth `y`|
# | `roc(yhat, y)` | compute points on the receiver-operator Characteristic|
# | `StratifiedCV(nfolds=6)` | 6-fold stratified cross-validation resampling strategy|
# | `Holdout(fraction_train=0.7)` | holdout resampling strategy|
# | `evaluate(model, X, y; resampling=..., options...)` | estimate performance metrics for `model` using the data `X`, `y`|
# | `FeatureSelector()` | transformer for selecting features|
# | `Step(3)` | iteration control for stepping 3 iterations|
# | `NumberSinceBest(6)`, `TimeLimit(60/5), InvalidValue()` | stopping criterion iteration controls|
# | `IteratedModel(model=..., controls=..., options...)` | wrap an iterative `model` in controls|
# | `range(model,  :some_hyperparam, lower=..., upper=...)` | define a numeric range|
# | `RandomSearch()` | random search tuning strategy|
# | `TunedModel(model=..., tuning=..., options...)` | wrap the supervised `model` in specified `tuning` strategy|

# ## Warm up: Building a model for the iris dataset

# Before turning to the Telco Customer Churn dataset, we very quickly
# build a predictive model for Fisher's well-known iris data set, as way of
# introducing the main actors in any MLJ workflow. Details that you
# don't fully grasp should become clearer in the Telco study.

# This section is a condensed adaption of the [Getting Started
# example](https://alan-turing-institute.github.io/MLJ.jl/dev/getting_started/#Fit-and-predict)
# in the MLJ documentation.

# First, using the built-in iris dataset, we load and inspect the features
# `X_iris` (a table) and target variable `y_iris` (a vector):

using MLJ

#-

X_iris, y_iris = @load_iris;
schema(X_iris)

#-

y_iris[1:4]

#

levels(y_iris)

# We load a decision tree model, from the package DecisionTree.jl:

DecisionTree = @load DecisionTreeClassifier pkg=DecisionTree # model type
model = DecisionTree(min_samples_split=5)                    # model instance

# In MLJ, a *model* is just a container for hyper-parameters of
# some learning algorithm. It does not store learned parameters.

# Next, we bind the model together with the available data in what's
# called a *machine*:

mach = machine(model, X_iris, y_iris)

# A machine is essentially just a model (ie, hyper-parameters) plus data, but
# it additionally stores *learned parameters* (the tree) once it is
# trained on some view of the data:

train_rows = vcat(1:60, 91:150); # some row indices (observations are rows not columns)
fit!(mach, rows=train_rows)
fitted_params(mach)

# A machine stores some other information enabling [warm
# restart](https://alan-turing-institute.github.io/MLJ.jl/dev/machines/#Warm-restarts)
# for some models, but we won't go into that here. You are allowed to
# access and mutate the `model` parameter:

mach.model.min_samples_split  = 10
fit!(mach, rows=train_rows) # re-train with new hyper-parameter

# Now we can make predictions on some other view of the data, as in

predict(mach, rows=71:73)

# or on completely new data, as in

Xnew = (sepal_length = [5.1, 6.3],
        sepal_width = [3.0, 2.5],
        petal_length = [1.4, 4.9],
        petal_width = [0.3, 1.5])
yhat = predict(mach, Xnew)

# These are probabilistic predictions which can be manipulated using a
# widely adopted interface defined in the Distributions.jl
# package. For example, we can get raw probabilities like this:

pdf.(yhat, "virginica")

# A single prediction is displayed like this:

yhat[2]

# We now turn to the Telco dataset.

# ## Getting the Telco data

import DataFrames

#-

data = OpenML.load(42178) # data set from OpenML.org
df0 = DataFrames.DataFrame(data)
first(df0, 4)

# The object of this tutorial is to build and evaluate supervised
# learning models to predict the `Churn` variable, a binary variable
# measuring customer retention, based on other variables that are
# relevant.

# In the table, observations correspond to rows, and features to
# columns, which is the convention for representing all
# two-dimensional data in MLJ.

# ## Type coercion

# *Introduces:* `scitype`, `schema`, `coerce`

# A ["scientific
# type"](https://juliaai.github.io/ScientificTypes.jl/dev/) or
# *scitype* indicates how MLJ will *interpret* data. For example,
# `typeof(3.14) == Float64`, while `scitype(3.14) == Continuous` and
# also `scitype(3.14f0) == Continuous`. In MLJ, model data
# requirements are articulated using scitypes.

# Here are common "scalar" scitypes:

# \figalt{scalar scitypesb}{./scitypes.svg}

# There are also container scitypes. For example, the scitype of any
# `N`-dimensional array is `AbstractArray{S, N}`, where `S` is the scitype of the
# elements:

scitype(["cat", "mouse", "dog"])

# The `schema` operator summarizes the column scitypes of a table:

schema(df0) |> DataFrames.DataFrame  # converted to DataFrame for better display

# All of the fields being interpreted as `Textual` are really
# something else, either `Multiclass` or, in the case of
# `TotalCharges`, `Continuous`. In fact, `TotalCharges` is
# mostly floats wrapped as strings. However, it needs special
# treatment because some elements consist of a single space, " ",
# which we'll treat as "0.0".

fix_blanks(v) = map(v) do x
    if x == " "
        return "0.0"
    else
        return x
    end
end

df0.TotalCharges = fix_blanks(df0.TotalCharges);

# Coercing the `TotalCharges` type to ensure a `Continuous` scitype:

coerce!(df0, :TotalCharges => Continuous);

# Coercing all remaining `Textual` data to `Multiclass`:

coerce!(df0, Textual => Multiclass);

# Finally, we'll coerce our target variable `Churn` to be
# `OrderedFactor`, rather than `Multiclass`, to enable a reliable
# interpretation of metrics like "true positive rate".  By convention,
# the first class is the negative one:

coerce!(df0, :Churn => OrderedFactor)
levels(df0.Churn) # to check order

# Re-inspecting the scitypes:

schema(df0) |> DataFrames.DataFrame

# ## Preparing a holdout set for final testing

# *Introduces:* `partition`

# To reduce training times for the purposes of this tutorial, we're
# going to dump 90% of observations (after shuffling) and split off
# 30% of the remainder for use as a lock-and-throw-away-the-key
# holdout set:

df, df_test, df_dumped = partition(df0, 0.07, 0.03, # in ratios 7:3:90
                                   stratify=df0.Churn,
                                   rng=123);

# The reader interested in including all data can instead do
# `df, df_test = partition(df0, 0.7, rng=123)`.

# ## Splitting data into target and features

# *Introduces:* `unpack`

# In the following call, the column with name `Churn` is copied over
# to a vector `y`, and every remaining column, except `customerID`
# (which contains no useful information) goes into a table `X`. Here
# `Churn` is the target variable for which we seek predictions, given
# new versions of the features `X`.

y, X = unpack(df, ==(:Churn), !=(:customerID));
schema(X).names

#-

intersect([:Churn, :customerID], schema(X).names)

# We'll do the same for the holdout data:

ytest, Xtest = unpack(df_test, ==(:Churn), !=(:customerID));

# ## Loading a model and checking type requirements

# *Introduces:* `@load`, `input_scitype`, `target_scitype`

# For tools helping us to identify suitable models, see the [Model
# Search](https://alan-turing-institute.github.io/MLJ.jl/dev/model_search/#model_search)
# section of the manual. We will build a gradient tree-boosting model,
# a popular first choice for structured data like we have here. Model
# code is contained in a third-party package called
# [EvoTrees.jl](https://github.com/Evovest/EvoTrees.jl) which is
# loaded as follows:

Booster = @load EvoTreeClassifier pkg=EvoTrees

# Recall that a *model* is just a container for some algorithm's
# hyper-parameters. Let's create a `Booster` with default values for
# the hyper-parameters:

booster = Booster()

# This model is appropriate for the kind of target variable we have because of
# the following passing test:

scitype(y) <: target_scitype(booster)

# However, our features `X` cannot be directly used with `booster`:

scitype(X) <: input_scitype(booster)

# As it turns out, this is because `booster`, like the majority of MLJ
# supervised models, expects the features to be `Continuous`. (With
# some experience, this can be gleaned from `input_scitype(booster)`.)
# So we need categorical feature encoding, discussed next.

# ## Building a model pipeline to incorporate feature encoding

# *Introduces:* `ContinuousEncoder`, pipeline operator `|>`

# The built-in `ContinuousEncoder` model transforms an arbitrary table
# to a table whose features are all `Continuous` (dropping any fields
# it does not know how to encode). In particular, all `Multiclass`
# features are one-hot encoded.

# A *pipeline* is a stand-alone model that internally combines one or
# more models in a linear (non-branching) pipeline. Here's a pipeline
# that adds the `ContinuousEncoder` as a pre-processor to the
# gradient tree-boosting model above:

pipe = ContinuousEncoder() |> booster

# Note that the component models appear as hyper-parameters of
# `pipe`. Pipelines are an implementation of a more general [model
# composition](https://alan-turing-institute.github.io/MLJ.jl/dev/composing_models/#Composing-Models)
# interface provided by MLJ that advanced users may want to learn about.

# From the above display, we see that component model hyper-parameters
# are now *nested*, but they are still accessible (important in hyper-parameter
# optimization):

pipe.evo_tree_classifier.max_depth

# ## Evaluating the pipeline model's performance

# *Introduces:* `measures` (function), **measures:** `brier_loss`,
#  `auc`, `accuracy`; `machine`, `fit!`, `predict`, `fitted_params`,
#  `report`, `roc`, **resampling strategy** `StratifiedCV`,
#  `evaluate`, `FeatureSelector`

# Without touching our test set `Xtest`, `ytest`, we will estimate the
# performance of our pipeline model, with default hyper-parameters, in
# two different ways:

# **Evaluating by hand.** First, we'll do this "by hand" using the `fit!` and `predict`
# workflow illustrated for the iris data set above, using a
# holdout resampling strategy. At the same time we'll see how to
# generate a **confusion matrix**, **ROC curve**, and inspect
# **feature importances**.

# **Automated performance evaluation.** Next we'll apply the more
# typical and convenient `evaluate` workflow, but using `StratifiedCV`
# (stratified cross-validation) which is more informative.

# In any case, we need to choose some measures (metrics) to quantify
# the performance of our model. For a complete list of measures, one
# does `measures()`. Or we also can do:

measures("Brier")

# We will be primarily using `brier_loss`, but also `auc` (area under
# the ROC curve) and `accuracy`.

# ### Evaluating by hand (with a holdout set)

# Our pipeline model can be trained just like the decision tree model
# we built for the iris data set. Binding all non-test data to the
# pipeline model:

mach_pipe = machine(pipe, X, y)

# We already encountered the `partition` method above. Here we apply
# it to row indices, instead of data containers, as `fit!` and
# `predict` only need a *view* of the data to work.

train, validation = partition(1:length(y), 0.7)
fit!(mach_pipe, rows=train)

# We note in passing that we can access two kinds of information from a trained machine:

# - The **learned parameters** (eg, coefficients of a linear model): We use `fitted_params(mach_pipe)`
# - Other **by-products of training** (eg, feature importances): We use `report(mach_pipe)`

fp = fitted_params(mach_pipe);
keys(fp)

# For example, we can check that the encoder did not actually drop any features:

Set(fp.continuous_encoder.features_to_keep) == Set(schema(X).names)

# And, from the report, extract feature importances:

rpt = report(mach_pipe)
keys(rpt.evo_tree_classifier)

#-

fi = rpt.evo_tree_classifier.feature_importances
feature_importance_table =
    (feature=Symbol.(first.(fi)), importance=last.(fi)) |> DataFrames.DataFrame

# For models not reporting feature importances, we recommend the
# [Shapley.jl](https://expandingman.gitlab.io/Shapley.jl/) package.

# Returning to predictions and evaluations of our measures:

ŷ = predict(mach_pipe, rows=validation);
print(
    "Measurements:\n",
    "  brier loss: ", brier_loss(ŷ, y[validation]) |> mean, "\n",
    "  auc:        ", auc(ŷ, y[validation]),                "\n",
    "  accuracy:   ", accuracy(mode.(ŷ), y[validation])
)

# Note that we need `mode` in the last case because `accuracy` expects
# point predictions, not probabilistic ones. (One can alternatively
# use `predict_mode` to generate the predictions.)

# While we're here, lets also generate a **confusion matrix** and
# [receiver-operator
# characteristic](https://en.wikipedia.org/wiki/Receiver_operating_characteristic)
# (ROC):

confmat(mode.(ŷ), y[validation])

# Note: Importing the plotting package and calling the plotting
# functions for the first time can take a minute or so.

using Plots

#-

roc_curve = roc(ŷ, y[validation])
plt = scatter(roc_curve, legend=false)
plot!(plt, xlab="false positive rate", ylab="true positive rate")
plot!([0, 1], [0, 1], linewidth=2, linestyle=:dash, color=:black)

savefig(joinpath(@OUTPUT, "EX-telco-roc.svg")) # hide

# \fig{EX-telco-roc.svg}

# ### Automated performance evaluation (more typical workflow)

# We can also get performance estimates with a single call to the
# `evaluate` function, which also allows for more complicated
# resampling - in this case stratified cross-validation. To make this
# more comprehensive, we set `repeats=3` below to make our
# cross-validation "Monte Carlo" (3 random size-6 partitions of the
# observation space, for a total of 18 folds) and set
# `acceleration=CPUThreads()` to parallelize the computation.

# We choose a `StratifiedCV` resampling strategy; the complete list of options is
# [here](https://alan-turing-institute.github.io/MLJ.jl/dev/evaluating_model_performance/#Built-in-resampling-strategies).

e_pipe = evaluate(pipe, X, y,
                  resampling=StratifiedCV(nfolds=6, rng=123),
                  measures=[brier_loss, auc, accuracy],
                  repeats=3,
                  acceleration=CPUThreads())

# (There is also a version of `evaluate` for machines. Query the
# `evaluate` and `evaluate!` doc-strings to learn more about these
# functions and what the `PerformanceEvaluation` object `e_pipe` records.)

# While [less than ideal](https://arxiv.org/abs/2104.00673), let's
# adopt the common practice of using the standard error of a
# cross-validation score as an estimate of the uncertainty of a
# performance measure's expected value. Here's a utility function to
# calculate 95% confidence intervals for our performance estimates
# based on this practice, and it's application to the current
# evaluation:

using Measurements

#-

function confidence_intervals(e)
    factor = 2.0 # to get level of 95%
    measure = e.measure
    nfolds = length(e.per_fold[1])
    measurement = [e.measurement[j] ± factor*std(e.per_fold[j])/sqrt(nfolds - 1)
                   for j in eachindex(measure)]
    table = (measure=measure, measurement=measurement)
    return DataFrames.DataFrame(table)
end

confidence_intervals_basic_model = confidence_intervals(e_pipe)

# ## Filtering out unimportant features

# *Introduces:* `FeatureSelector`

# Before continuing, we'll modify our pipeline to drop those features
# with low feature importance, to speed up later optimization:

unimportant_features = filter(:importance => <(0.005), feature_importance_table).feature

pipe2 = ContinuousEncoder() |>
    FeatureSelector(features=unimportant_features, ignore=true) |> booster

# ## Wrapping our iterative model in control strategies

# *Introduces:* **control strategies:** `Step`, `NumberSinceBest`, `TimeLimit`,
# `InvalidValue`, **model wrapper** `IteratedModel`, **resampling strategy:** `Holdout`

# We want to optimize the hyper-parameters of our model. Since our
# model is iterative, these parameters include the (nested) iteration
# parameter `pipe.evo_tree_classifier.nrounds`. Sometimes this
# parameter is optimized first, fixed, and then maybe optimized again
# after the other parameters. Here we take a more principled approach,
# **wrapping our model in a control strategy** that makes it
# "self-iterating". The strategy applies a stopping criterion to
# *out-of-sample* estimates of the model performance, constructed
# using an internally constructed holdout set. In this way, we avoid
# some data hygiene issues, and, when we subsequently optimize other
# parameters, we will always being using an optimal number of
# iterations.

# Note that this approach can be applied to any iterative MLJ model,
# eg, the neural network models provided by
# [MLJFlux.jl](https://github.com/FluxML/MLJFlux.jl).

# First, we select appropriate controls from [this
# list](https://alan-turing-institute.github.io/MLJ.jl/dev/controlling_iterative_models/#Controls-provided):

controls = [
    Step(1),              # to increment iteration parameter (`pipe.nrounds`)
    NumberSinceBest(4),   # main stopping criterion
    TimeLimit(2/3600),    # never train more than 2 sec
    InvalidValue()        # stop if NaN or ±Inf encountered
]

# Now we wrap our pipeline model using the `IteratedModel` wrapper,
# being sure to specify the `measure` on which internal estimates of
# the out-of-sample performance will be based:

iterated_pipe = IteratedModel(model=pipe2,
                              controls=controls,
                              measure=brier_loss,
                              resampling=Holdout(fraction_train=0.7))

# We've set `resampling=Holdout(fraction_train=0.7)` to arrange that
# data attached to our model should be internally split into a train
# set (70%) and a holdout set (30%) for determining the out-of-sample
# estimate of the Brier loss.

# For demonstration purposes, let's bind `iterated_model` to all data
# not in our don't-touch holdout set, and train on all of that data:

mach_iterated_pipe = machine(iterated_pipe, X, y)
fit!(mach_iterated_pipe);

# To recap, internally this training is split into two separate steps:

# - A controlled iteration step, training on the holdout set, with the total number of iterations determined by the specified stopping criteria (based on the out-of-sample performance estimates)
# - A final step that trains the atomic model on *all* available
#   data using the number of iterations determined in the first step. Calling `predict` on `mach_iterated_pipe` means using the learned parameters of the second step.

# ## Hyper-parameter optimization (model tuning)

# *Introduces:* `range`, **model wrapper** `TunedModel`, `RandomSearch`

# We now turn to hyper-parameter optimization. A tool not discussed
# here is the `learning_curve` function, which can be useful when
# wanting to visualize the effect of changes to a *single*
# hyper-parameter (which could be an iteration parameter). See, for
# example, [this section of the
# manual](https://alan-turing-institute.github.io/MLJ.jl/dev/learning_curves/)
# or [this
# tutorial](https://github.com/ablaom/MLJTutorial.jl/blob/dev/notebooks/04_tuning/notebook.ipynb).

# Fine tuning the hyper-parameters of a gradient booster can be
# somewhat involved. Here we settle for simultaneously optimizing two
# key parameters: `max_depth` and `η` (learning_rate).

# Like iteration control, **model optimization in MLJ is implemented as
# a model wrapper**, called `TunedModel`. After wrapping a model in a
# tuning strategy and binding the wrapped model to data in a machine
# called `mach`, calling `fit!(mach)` instigates a search for optimal
# model hyperparameters, within a specified range, and then uses all
# supplied data to train the best model. To predict using that model,
# one then calls `predict(mach, Xnew)`. In this way the wrapped model
# may be viewed as a "self-tuning" version of the unwrapped
# model. That is, wrapping the model simply transforms certain
# hyper-parameters into learned parameters (just as `IteratedModel`
# does for an iteration parameter).

# To start with, we define ranges for the parameters of
# interest. Since these parameters are nested, let's force a
# display of our model to a larger depth:

show(iterated_pipe, 2)

#-

p1 = :(model.evo_tree_classifier.η)
p2 = :(model.evo_tree_classifier.max_depth)

r1 = range(iterated_pipe, p1, lower=-2, upper=-0.5, scale=x->10^x)
r2 = range(iterated_pipe, p2, lower=2, upper=6)

# Nominal ranges are defined by specifying `values` instead of `lower`
# and `upper`.

# Next, we choose an optimization strategy from [this
# list](https://alan-turing-institute.github.io/MLJ.jl/dev/tuning_models/#Tuning-Models):

tuning = RandomSearch(rng=123)

# Then we wrap the model, specifying a `resampling` strategy and a
# `measure`, as we did for `IteratedModel`.  In fact, we can include a
# battery of `measures`; by default, optimization is with respect to
# performance estimates based on the first measure, but estimates for
# all measures can be accessed from the model's `report`.

# The keyword `n` specifies the total number of models (sets of
# hyper-parameters) to evaluate.

tuned_iterated_pipe = TunedModel(model=iterated_pipe,
                                 range=[r1, r2],
                                 tuning=tuning,
                                 measures=[brier_loss, auc, accuracy],
                                 resampling=StratifiedCV(nfolds=6, rng=123),
                                 acceleration=CPUThreads(),
                                 n=40)

# To save time, we skip the `repeats` here.

# Binding our final model to data and training:

mach_tuned_iterated_pipe = machine(tuned_iterated_pipe, X, y)
fit!(mach_tuned_iterated_pipe)

# As explained above, the training we have just performed was split
# internally into two separate steps:

# - A step to determine the parameter values that optimize the aggregated cross-validation scores
# - A final step that trains the optimal model on *all* available data. Future predictions `predict(mach_tuned_iterated_pipe, ...)` are based on this final training step.

# From `report(mach_tuned_iterated_pipe)` we can extract details about
# the optimization procedure. For example:

rpt2 = report(mach_tuned_iterated_pipe);
best_booster = rpt2.best_model.model.evo_tree_classifier

#-

print(
    "Optimal hyper-parameters: \n",
    "  max_depth: ", best_booster.max_depth, "\n",
    "  η:         ", best_booster.η
)

# Using the `confidence_intervals` function we defined earlier:

e_best = rpt2.best_history_entry
confidence_intervals(e_best)

# Digging a little deeper, we can learn what stopping criterion was
# applied in the case of the optimal model, and how many iterations
# were required:

rpt2.best_report.controls |> collect

# Finally, we can visualize the optimization results:

plot(mach_tuned_iterated_pipe, size=(600,450))

savefig(joinpath(@OUTPUT, "EX-telco-tuning.svg")) # hide

# \fig{EX-telco-tuning.svg}

# ## Saving our model

# *Introduces:* `MLJ.save`

# Here's how to serialize our final, trained self-iterating,
# self-tuning pipeline machine using Julia's native serializer (see
# [the
# manual](https://alan-turing-institute.github.io/MLJ.jl/dev/machines/#Saving-machines)
# for more options):

MLJ.save("tuned_iterated_pipe.jls", mach_tuned_iterated_pipe)

# We'll deserialize this in "Testing the final model" below.

# ## Final performance estimate

# Finally, to get an even more accurate estimate of performance, we
# can evaluate our model using stratified cross-validation and all the
# data attached to our machine. Because this evaluation implies
# [nested
# resampling](https://mlr.mlr-org.com/articles/tutorial/nested_resampling.html),
# this computation takes quite a bit longer than the previous one
# (which is being repeated six times, using 5/6th of the data each
# time):

e_tuned_iterated_pipe = evaluate(tuned_iterated_pipe, X, y,
                                 resampling=StratifiedCV(nfolds=6, rng=123),
                                 measures=[brier_loss, auc, accuracy])

#-

confidence_intervals(e_tuned_iterated_pipe)

# For comparison, here are the confidence intervals for the basic
# pipeline model (no feature selection and default hyperparameters):

confidence_intervals_basic_model

# As each pair of intervals overlap, it's doubtful the small changes
# here can be assigned statistical significance. Default `booster`
# hyper-parameters do a pretty good job.

# ## Testing the final model

# We now determine the performance of our model on our
# lock-and-throw-away-the-key holdout set. To demonstrate
# deserialization, we'll pretend we're in a new Julia session (but
# have called import/using on the same packages). Then the
# following should suffice to recover our model trained under
# "Hyper-parameter optimization" above:

mach_restored = machine("tuned_iterated_pipe.jls")

# We compute predictions on the holdout set:

ŷ_tuned = predict(mach_restored, Xtest);
ŷ_tuned[1]

# And can compute the final performance measures:

print(
    "Tuned model measurements on test:\n",
    "  brier loss: ", brier_loss(ŷ_tuned, ytest) |> mean, "\n",
    "  auc:        ", auc(ŷ_tuned, ytest),                "\n",
    "  accuracy:   ", accuracy(mode.(ŷ_tuned), ytest)
)

# For comparison, here's the performance for the basic pipeline model

mach_basic = machine(pipe, X, y)
fit!(mach_basic, verbosity=0)

ŷ_basic = predict(mach_basic, Xtest);

print(
    "Basic model measurements on test set:\n",
    "  brier loss: ", brier_loss(ŷ_basic, ytest) |> mean, "\n",
    "  auc:        ", auc(ŷ_basic, ytest),                "\n",
    "  accuracy:   ", accuracy(mode.(ŷ_basic), ytest)
)

rm("tuned_iterated_pipe.jls") # hide

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
using OffsetArrays

include("ReadDataFile.jl")
include("SetProb.jl")

using .ReadDataFile
using .SetProb

include("qinit.jl")
include("out1.jl")
include("src1.jl")
include("bc1.jl")
include("rp1_acoustics.jl")
include("b4step1.jl")
include("claw1.jl")

#
#
#
#     ===============================================================
function claw1ez()    # No arguments
#     ===============================================================
#
#     An easy-to-use clawpack driver routine for simple applications
#     Documentation is available at
#                 http://www.amath.washington.edu/~claw/doc.html
#
#     Authors: Randall J. LeVeque, Grady I. Lemoine
#
#      implicit double precision (a-h,o-z)
#      external bc1,rp1,src1,b4step1

#      double precision, dimension(:,:), allocatable :: q, aux
#      double precision, dimension(:), allocatable :: work, tout
#      integer, dimension(:), allocatable :: mthlim, iout_q, iout_aux

#      dimension method(7),dtv(5),cflv(4),nv(2),mthbc(2)
    method = Vector{Int}(undef, 7)
    dtv    = Vector{Float64}(undef, 5)
    cflv   = Vector{Float64}(undef, 4)
    nv     = Vector{Int}(undef, 2)
    mthbc  = Vector{Int}(undef, 2)

###       integer :: allocate_status, outstyle
###       logical :: outaux_init_only, use_fwaves, output_t0
###       logical :: outaux_always, dt_variable

    outaux_init_only = false
    outaux_always = false

###       character*12 fname
### #
###       open(10,file='fort.info',status='unknown',form='formatted')

    unit10 = open("fort.info", "w")
###
### #
### #     # Read the input in standard form from claw.data:
### #     # For a description of input parameters see the documentation at
### #                 http://www.amath.washington.edu/~claw
###
### #     ## New in 4.4:   Input file name changed from claw1ez.data
### #     ## Open file and skip over leading lines with # comments:
###       fname = 'claw.data'
###       call opendatafile(55,fname)
###
### #
### #     # Read the input in standard form from claw.data:
### #
###
    r = Reader("claw.data")

    ndim     = read_int(r)
    xlower   = read_flt(r)
    xupper   = read_flt(r)
    mx       = read_int(r)
    meqn     = read_int(r)
    mwaves   = read_int(r)
    maux     = read_int(r)
    t0       = read_flt(r)

    println("read from claw.data : ndim: $ndim xlower: $xlower xupper: $xupper mx: $mx meqn: $meqn mwaves: $mwaves maux: $maux t0: $t0")

    outstyle = read_int(r)

    tout = Vector{Float64}
    if     outstyle == 1
        nout      = read_int(r)
        tfinal    = read_flt(r)
        output_t0 = read_bool(r)       # Not currently used
        nstepout  = 1

        println("nout: $nout tfinal: $tfinal output_t0: $output_t0")

    elseif outstyle == 2
        nout      = read_int(r)
        tout      = read_flt_array(r, nout)
        println("tout: $tout")
        nstepout  = 1

    elseif outstyle == 3
        nstepout  = read_int(r)
        nstop     = read_int(r)
        output_t0 = read_bool(r)       # Not currently used
        nout      = nstop

    else
        println("nout: $nout tfinal: $tfinal output_t0: $output_t0")
        println("*** Unrecognized output style $outstyle")
        error("*** Exiting claw1ez")
    end

    output_format = read_int(r)    # Not used yet

    iout_q = read_flt_array(r, meqn)
    println("iout_q: $iout_q")

#        iout_aux

    iout_aux = Vector{Int}
    if maux > 0
        iout_aux = read_int_array(r, maux)
        outaux_init_only = read_bool(r)

        # TODO
    else
        outaux_always = false
        outaux_init_only = false    # Just to initialize
    end

    dtv[1]  = read_flt(r)    # Initial dt
    dtv[2]  = read_flt(r)    # Max dt
    cflv[1] = read_flt(r)    # Max CFL number
    cflv[2] = read_flt(r)    # Desired CFL number
    nv[1]   = read_int(r)    # Maximum number of steps

    println("dtv: $dtv cflv: $cflv nv: $nv")

    dt_variable::Bool = read_bool(r)
    method[1] = dt_variable ? 1 : 0

    method[2] = read_int(r) # Order

    # method(3) (transverse order) not used in 1D
    # No dimensional splitting in 1D

    method[4] = read_int(r) # Verbosity
    method[5] = read_int(r) # Source term splitting style
    method[6] = read_int(r) # Index into aux of capacity function
    method[7] = maux       # Number of aux variables

    println("method: $method")

    use_fwaves::Bool = read_bool(r)

    mthlim::Array{Int, 1} = read_int_array(r, mwaves)
    println("mthlim: $mthlim")

    mbc      = read_int(r)
    mthbc[1] = read_int(r)
    mthbc[2] = read_int(r)

    println("mbc: $mbc mthbc: $mthbc")

    # No restart in 1D, and there's nothing after the restart info, so
    # just close the file now.
    close_reader(r)

    # Check consistency for periodic BCs
    if (mthbc[1] == 2 && mthbc[2] != 2) ||
       (mthbc[2] == 2 && mthbc[1] != 2)
          error("*** ERROR ***  periodic boundary conditions require mthbc(1) and mthbc(2) BOTH be set to 2")
    end

    println("running...")
    println(" ")

    # grid spacing
    dx = (xupper - xlower) / float(mx)

    # time increments between outputing solution:
    if outstyle == 1
        dtout = (tfinal - t0)/float(nout)
    end

    # call user's routine setprob to set any specific parameters
    # or other initialization required.
    #
    prob_data = CParam()
    read_prob_data(prob_data)

    # Allocate aux

    aux = (maux > 0) ? OffsetArray{Float64}(undef, 1:maux, 1-mbc:mx+mbc) : OffsetArray{Float64}(undef, 1:1, 1:1)

    q = OffsetArray{Float64}(undef, 1:meqn, 1-mbc:mx+mbc)

    # set initial conditions:

    qinit(prob_data, meqn, mbc, mx, xlower, dx, q, maux, aux)

    # output initial data
    out1(meqn,mbc,mx,xlower,dx,q,t0,0,aux,maux,
         (outaux_init_only || outaux_always) )
#
#     ----------
#     Main loop:
#     ----------
#
    tend = t0
    for n=1:nout
        tstart = tend
        if outstyle == 1
            tend = tstart + dtout
        elseif outstyle == 2
            tend = tout[n]
        elseif outstyle == 3
            tend = tstart - 1.0  # single-step mode
        end

        info::Int =
        claw1(prob_data,meqn,mwaves,maux,mbc,mx,
              q,aux,xlower,dx,tstart,tend,dtv,cflv,nv,method,mthlim,
              mthbc,
              use_fwaves,
              bc1,rp1,src1,b4step1)

        # check to see if an error occured:
        if info != 0
            println("*** ERROR in claw1 ***  info = $info")
            if info == 1
                println("***   either mx > maxmx or mbc < 2")
            elseif info == 2
                println("***   dt does not divide (tend - tstart)")
                println("***   and dt is fixed since method(1)=0")
            elseif info == 3
                println("***   method(1)=1 and cflv(2) > cflv(1)")
            elseif info == 4
                println("***   mwork is too small")
            elseif info == 11
                println("***   Too many times steps, n > nv(1)")
            elseif info == 12
                println("***   The Courant number is greater than cflv(1)")
                println("***   and dt is fixed since method(1)=0")
            end

            return

        end

        dtv[1] = dtv[5]  # use final dt as starting value on next call

        # output solution at this time
        # ------------------------------

        # if outstyle=1 or 2, then nstepout=1 and we output every time
        # we reach this point, since claw1 was called for the entire time
        # increment between outputs.

        # if outstyle=3 then we only output if we have taken nstepout
        # time steps since the last output.

        # iframe is the frame number used to form file names in out1
        iframe::Int = convert(Int, n/nstepout)
        if iframe*nstepout == n
            out1(meqn,mbc,mx,xlower,dx,q,tend,iframe,
                 aux,maux,outaux_always)

            @printf("CLAW1EZ: Frame %4d output files done at time t = %12.4e\n", iframe,tend)
            @printf(unit10,"tend = %15.4e \ninfo =%5d\nsmallest dt =%15.4e\nlargest dt =%15.4e\n",
                    tend,info,dtv[3],dtv[4])
            @printf(unit10,"last dt =%15.4e\nlargest cfl =%15.4e\nlast cfl =%15.4e\nsteps taken =%4d\n",
                    dtv[5],cflv[3],cflv[4],nv[2])

        end
    end
    Base.close(unit10)

end # claw1ez