text:
"""
    sign_cadlag(x::N) where {N<:Real}

This function works like the sign function but is ``1`` for input ``0``.

### Input

- `x` -- real scalar

### Output

``1`` if ``x ≥ 0``, ``-1`` otherwise.

### Notes

This is the sign function right-continuous at zero (see
[càdlàg function](https://en.wikipedia.org/wiki/C%C3%A0dl%C3%A0g)).
It can be used with vector-valued arguments via the dot operator.

### Examples

```jldoctest
julia> LazySets.sign_cadlag.([-0.6, 1.3, 0.0])
3-element Array{Float64,1}:
 -1.0
  1.0
  1.0
```
"""
function sign_cadlag(x::N) where {N<:Real}
    return x < zero(x) ? -one(x) : one(x)
end

"""
    substitute(substitution::Dict{Int, T}, x::AbstractVector{T}) where {T}

Apply a substitution to a given vector.

### Input

- `substitution` -- substitution (a mapping from an index to a new value)
- `x`            -- vector

### Output

A fresh vector corresponding to `x` after `substitution` was applied.
"""
function substitute(substitution::Dict{Int, T}, x::AbstractVector{T}) where {T}
    return substitute!(substitution, copy(x))
end

"""
    substitute!(substitution::Dict{Int, T}, x::AbstractVector{T}) where {T}

Apply a substitution to a given vector.

### Input

- `substitution` -- substitution (a mapping from an index to a new value)
- `x`            -- vector (modified in this function)

### Output

The same (but see the Notes below) vector `x` but after `substitution` was
applied.

### Notes

The vector `x` is modified in-place if it has type `Vector` or `SparseVector`.
Otherwise, we first create a new `Vector` from it.
"""
function substitute!(substitution::Dict{Int, T}, x::AbstractVector{T}) where {T}
    return substitute!(Vector(x), substitution)
end

function substitute!(substitution::Dict{Int, T},
                     x::Union{Vector{T}, SparseVector{T}}) where {T}
    for (index, value) in substitution
        x[index] = value
    end
    return x
end

"""
    reseed(rng::AbstractRNG, seed::Union{Int, Nothing})

Reset the RNG seed if the seed argument is a number.

### Input

- `rng`  -- random number generator
- `seed` -- seed for reseeding

### Output

The input RNG if the seed is `nothing`, and a reseeded RNG otherwise.
"""
function reseed(rng::AbstractRNG, seed::Union{Int, Nothing})
    if seed != nothing
        return Random.seed!(rng, seed)
    end
    return rng
end

"""
    StrictlyIncreasingIndices

Iterator over the vectors of `m` strictly increasing indices from 1 to `n`.

### Fields

- `n` -- size of the index domain
- `m` -- number of indices to choose (resp. length of the vectors)

### Notes

The vectors are modified in-place.

The iterator ranges over ``\\binom{n}{m}`` (`n` choose `m`) possible vectors.

This implementation results in a lexicographic order with the last index growing
first.

### Examples

```jldoctest
julia> for v in LazySets.StrictlyIncreasingIndices(4, 2)
           println(v)
       end
[1, 2]
[1, 3]
[1, 4]
[2, 3]
[2, 4]
[3, 4]
```
"""
struct StrictlyIncreasingIndices
    n::Int
    m::Int

    function StrictlyIncreasingIndices(n::Int, m::Int)
        @assert n >= m > 0 "require n >= m > 0"
        new(n, m)
    end
end

Base.eltype(::Type{StrictlyIncreasingIndices}) = Vector{Int}
Base.length(sii::StrictlyIncreasingIndices) = binomial(sii.n, sii.m)

# initialization
function Base.iterate(sii::StrictlyIncreasingIndices)
    v = [1:sii.m;]
    return (v, v)
end

# normal iteration
function Base.iterate(sii::StrictlyIncreasingIndices, state::AbstractVector{Int})
    v = state
    i = sii.m
    diff = sii.n
    if i == diff
        return nothing
    end
    while v[i] == diff
        i -= 1
        diff -= 1
    end
    # update vector
    v[i] += 1
    for j in i+1:sii.m
        v[j] = v[j-1] + 1
    end
    # detect termination: first index has maximum value
    if i == 1 && v[1] == (sii.n - sii.m + 1)
        return (v, nothing)
    end
    return (v, v)
end

# termination
function Base.iterate(sii::StrictlyIncreasingIndices, state::Nothing)
    return nothing
end

"""
    subtypes(interface, concrete::Bool)

Return the concrete subtypes of a given interface.

### Input

- `interface` -- an abstract type, usually a set interface
- `concrete`  -- if `true`, seek further the inner abstract subtypes of the given
                 interface, otherwise return only the direct subtypes of `interface`

### Output

A list with the subtypes of the abstract type `interface`, sorted alphabetically.

### Examples

Consider the `AbstractPolytope` interface. If we include the abstract subtypes
of this interface,

```jldoctest subtypes
julia> using LazySets: subtypes

julia> subtypes(AbstractPolytope, false)
4-element Array{Any,1}:
 AbstractCentrallySymmetricPolytope
 AbstractPolygon
 HPolytope
 VPolytope
```

We can use this function to obtain the concrete subtypes of
`AbstractCentrallySymmetricPolytope` and `AbstractPolygon` (further until all
concrete types are obtained), using the `concrete` flag:

```jldoctest subtypes
julia> subtypes(AbstractPolytope, true)
14-element Array{Type,1}:
 Ball1
 BallInf
 HPolygon
 HPolygonOpt
 HPolytope
 Hyperrectangle
 Interval
 LineSegment
 Singleton
 SymmetricIntervalHull
 VPolygon
 VPolytope
 ZeroSet
 Zonotope
```
"""
function subtypes(interface, concrete::Bool)

    subtypes_to_test = subtypes(interface)

    # do not seek the concrete subtypes further
    if !concrete
        return sort(subtypes_to_test, by=string)
    end

    result = Vector{Type}()
    i = 0
    while i < length(subtypes_to_test)
        i += 1
        subtype = subtypes_to_test[i]
        new_subtypes = subtypes(subtype)
        if isempty(new_subtypes)
            # base type found
            push!(result, subtype)
        else
            # yet another interface layer
            append!(subtypes_to_test, new_subtypes)
        end
    end
    return sort(result, by=string)
end

"""
    implementing_sets(op::Function;
                      signature::Tuple{Vector{Type}, Int}=(Type[], 1),
                      type_args=Float64, binary::Bool=false)

Compute a dictionary containing information about availability of (unary or
binary) concrete set operations.

### Input

- `op`        -- set operation (respectively its `Function` object)
- `signature` -- (optional, default: `Type[]`) the type signature of the
                 function without the `LazySet` type(s) (see also the `index`
                 option and the `Examples` section below)
- `index`     -- (optional, default: `1`) index of the set type in the signature
                 in the unary case (see the `binary` option)
- `type_args` -- (optional, default: `Float64`) type arguments added to the
                 `LazySet`(s) when searching for available methods; valid
                 inputs are a type or `nothing`, and in the unary case (see the
                 `binary` option) it can also be a list of types
- `binary`    -- (optional, default: `false`) flag indicating whether `op` is a
                 binary function (`true`) or a unary function (`false`)

### Output

A dictionary with three keys each mapping to a list:
- `"available"` -- This list contains all set types such that there exists an
                   implementation of `op`.
- `"missing"`   -- This list contains all set types such that there does not
                   exist an implementation of `op`. Note that this is the
                   complement of the `"available"` list.
- `"specific"`  -- This list contains all set types such that there exists a
                   type-specific implementation. Note that those set types also
                   occur in the `"available"` list.

In the unary case, the lists contain set types.
In the binary case, the lists contain pairs of set types.

### Examples

```jldoctests
julia> using LazySets: implementing_sets

julia> dict = implementing_sets(tovrep);

julia> dict["available"]  # tovrep is only available for polyhedral set types
6-element Array{Type,1}:
 HPolygon
 HPolygonOpt
 HPolyhedron
 HPolytope
 VPolygon
 VPolytope

julia> dict = implementing_sets(σ; signature=Type[AbstractVector{Float64}], index=2);

julia> dict["missing"]  # every set type implements function σ
0-element Array{Type,1}

julia> N = Rational{Int};  # restriction of the number type

julia> dict = implementing_sets(σ; signature=Type[AbstractVector{N}], index=2, type_args=N);

julia> dict["missing"]  # some set types are not available with number type N
4-element Array{Type,1}:
 Ball2
 Ballp
 Bloating
 Ellipsoid

julia> dict = LazySets.implementing_sets(convex_hull; binary=true);  # binary case

julia> (HPolytope, HPolytope) ∈ dict["available"]  # dict contains pairs now
true
```
"""
function implementing_sets(op::Function;
                           signature::Vector{Type}=Type[],
                           index::Int=1,
                           type_args=Float64,
                           binary::Bool=false)
    function get_unary_dict()
        return Dict{String, Vector{Type}}("available" => Vector{Type}(),
                                          "missing" => Vector{Type}(),
                                          "specific" => Vector{Type}())
    end

    if !binary
        # unary case
        dict = get_unary_dict()
        _implementing_sets_unary!(dict, op, signature, index, type_args)
        return dict
    end

    # binary case by considering all set types as second argument
    binary_dict = Dict{String, Vector{Tuple{Type, Type}}}(
        "available" => Vector{Tuple{Type, Type}}(),
        "missing" => Vector{Tuple{Type, Type}}(),
        "specific" => Vector{Tuple{Type, Type}}())
    signature_copy = copy(signature)
    insert!(signature_copy, 1, LazySet)  # insert a dummy type
    for set_type in subtypes(LazySet, true)
        augmented_set_type = set_type
        try
            if type_args isa Type
                augmented_set_type = set_type{type_args}
            else
                @assert type_args == nothing "for binary functions, " *
                    "`type_args` must not be a list"
            end
        catch e
            # type combination not available
            push!(binary_dict["missing"], (set_type, LazySet))
            continue
        end
        signature_copy[1] = augmented_set_type  # replace dummy type
        unary_dict = get_unary_dict()
        _implementing_sets_unary!(unary_dict, op, signature_copy, 1, type_args)

        # create pairs
        for key in ["available", "missing", "specific"]
            i = 1
            for other_set_type in unary_dict[key]
                push!(binary_dict[key], (set_type, other_set_type))
                i += 1
            end
        end
    end
    return binary_dict
end

function _implementing_sets_unary!(dict, op, signature, index, type_args)
    function create_type_tuple(given_type)
        if isempty(signature)
            return Tuple{given_type}
        end
        # create dynamic tuple
        new_signature = copy(signature)
        insert!(new_signature, index, given_type)
        return Tuple{new_signature...}
    end

    for set_type in subtypes(LazySet, true)
        augmented_set_type = set_type
        try
            if type_args isa Type
                augmented_set_type = set_type{type_args}
            elseif type_args != nothing
                augmented_set_type = set_type{type_args...}
            end
        catch e
            # type combination not available
            push!(dict["missing"], set_type)
            continue
        end
        tuple_set_type = create_type_tuple(augmented_set_type)
        if !hasmethod(op, tuple_set_type)
            # implementation not available
            push!(dict["missing"], set_type)
            continue
        end
        # implementation available
        push!(dict["available"], set_type)

        # check if there is a type-specific method
        super_type = supertype(augmented_set_type)
        tuple_super_type = create_type_tuple(super_type)
        if !hasmethod(op, tuple_super_type)
            continue
        end
        method_set_type = which(op, tuple_set_type)
        method_super_type = which(op, tuple_super_type)
        if method_set_type != method_super_type
            push!(dict["specific"], set_type)
        end
    end
end


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

# What is included in the file `FMI2_c.jl`?
# - default implementations for the `fmi2CallbackFunctions`
# - julia-implementaions of the functions inside the FMI-specification
# Any c-function `f(c::fmi2Component, args...)` in the spec is implemented as `f(c::FMU2Component, args...)`.
# Any c-function `f(args...)` without a leading `fmi2Component`-arguemnt is already implented as `f(c_ptr, args...)` in FMICore, where `c_ptr` is a pointer to the c-function (inside the DLL).

# already defined in FMICore.jl:
# - fmi2Instantiate

import FMICore: fmi2Instantiate, fmi2FreeInstance!, fmi2GetTypesPlatform, fmi2GetVersion
import FMICore: fmi2SetDebugLogging, fmi2SetupExperiment, fmi2EnterInitializationMode, fmi2ExitInitializationMode, fmi2Terminate, fmi2Reset
import FMICore: fmi2GetReal!, fmi2SetReal, fmi2GetInteger!, fmi2SetInteger, fmi2GetBoolean!, fmi2SetBoolean, fmi2GetString!, fmi2SetString
import FMICore: fmi2GetFMUstate!, fmi2SetFMUstate, fmi2FreeFMUstate!, fmi2SerializedFMUstateSize!, fmi2SerializeFMUstate!, fmi2DeSerializeFMUstate!
import FMICore: fmi2GetDirectionalDerivative!, fmi2SetRealInputDerivatives, fmi2GetRealOutputDerivatives!
import FMICore: fmi2DoStep, fmi2CancelStep, fmi2GetStatus!, fmi2GetRealStatus!, fmi2GetIntegerStatus!, fmi2GetBooleanStatus!, fmi2GetStringStatus!
import FMICore: fmi2SetTime, fmi2SetContinuousStates, fmi2EnterEventMode, fmi2NewDiscreteStates!, fmi2EnterContinuousTimeMode, fmi2CompletedIntegratorStep!
import FMICore: fmi2GetDerivatives!, fmi2GetEventIndicators!, fmi2GetContinuousStates!, fmi2GetNominalsOfContinuousStates!

"""
Source: FMISpec2.0.2[p.21]: 2.1.5 Creation, Destruction and Logging of FMU Instances

Function that is called in the FMU, usually if an fmi2XXX function, does not behave as desired. If “logger” is called with “status = fmi2OK”, then the message is a pure information message. “instanceName” is the instance name of the model that calls this function. “category” is the category of the message. The meaning of “category” is defined by the modeling environment that generated the FMU. Depending on this modeling environment, none, some or all allowed values of “category” for this FMU are defined in the modelDescription.xml file via element “<fmiModelDescription><LogCategories>”, see section 2.2.4. Only messages are provided by function logger that have a category according to a call to fmi2SetDebugLogging (see below). Argument “message” is provided in the same way and with the same format control as in function “printf” from the C standard library. [Typically, this function prints the message and stores it optionally in a log file.]
"""
function fmi2CallbackLogger(_componentEnvironment::Ptr{FMU2ComponentEnvironment},
            _instanceName::Ptr{Cchar},
            _status::Cuint,
            _category::Ptr{Cchar},
            _message::Ptr{Cchar})

    message = unsafe_string(_message)
    category = unsafe_string(_category)
    status = fmi2StatusToString(_status)
    instanceName = unsafe_string(_instanceName)
    componentEnvironment = unsafe_load(_componentEnvironment)

    if status == fmi2StatusOK && componentEnvironment.logStatusOK
        @info "[$status][$category][$instanceName]: $message"
    elseif (status == fmi2StatusWarning && componentEnvironment.logStatusWarning) ||
           (status == fmi2StatusPending && componentEnvironment.logStatusPending)
        @warn "[$status][$category][$instanceName]: $message"
    elseif (status == fmi2StatusDiscard && componentEnvironment.logStatusDiscard) ||
           (status == fmi2StatusError   && componentEnvironment.logStatusError) ||
           (status == fmi2StatusFatal   && componentEnvironment.logStatusFatal)
        @error "[$status][$category][$instanceName]: $message"
    end

    return nothing
end

# (cfmi2CallbackLogger, fmi2CallbackLogger) = Cfunction{                      fmi2ComponentEnvironment,               Ptr{Cchar},         Cuint,           Ptr{Cchar},          Tuple{Ptr{Cchar}, Vararg}   }() do componentEnvironment::fmi2ComponentEnvironment, instanceName::Ptr{Cchar}, status::Cuint, category::Ptr{Cchar}, message::Tuple{Ptr{Cchar}, Vararg}
#     printf(message)
#     nothing
# end

"""
Source: FMISpec2.0.2[p.21-22]: 2.1.5 Creation, Destruction and Logging of FMU Instances

Function that is called in the FMU if memory needs to be allocated. If attribute “canNotUseMemoryManagementFunctions = true” in <fmiModelDescription><ModelExchange / CoSimulation>, then function allocateMemory is not used in the FMU and a void pointer can be provided. If this attribute has a value of “false” (which is the default), the FMU must not use malloc, calloc or other memory allocation functions. One reason is that these functions might not be available for embedded systems on the target machine. Another reason is that the environment may have optimized or specialized memory allocation functions. allocateMemory returns a pointer to space for a vector of nobj objects, each of size “size” or NULL, if the request cannot be satisfied. The space is initialized to zero bytes [(a simple implementation is to use calloc from the C standard library)].
"""
function fmi2CallbackAllocateMemory(nobj::Csize_t, size::Csize_t)
    ptr = Libc.calloc(nobj, size)
    @debug "cbAllocateMemory($(nobj), $(size)): Allocated $(nobj) x $(size) bytes at $(ptr)."
	ptr
end

"""
Source: FMISpec2.0.2[p.22]: 2.1.5 Creation, Destruction and Logging of FMU Instances

Function that must be called in the FMU if memory is freed that has been allocated with allocateMemory. If a null pointer is provided as input argument obj, the function shall perform no action [(a simple implementation is to use free from the C standard library; in ANSI C89 and C99, the null pointer handling is identical as defined here)]. If attribute “canNotUseMemoryManagementFunctions = true” in <fmiModelDescription><ModelExchange / CoSimulation>, then function freeMemory is not used in the FMU and a null pointer can be provided.
"""
function fmi2CallbackFreeMemory(obj::Ptr{Cvoid})
    @debug "cbFreeMemory($(obj)): Freeing object at $(obj)."
	Libc.free(obj)
    nothing
end

"""
Source: FMISpec2.0.2[p.22]: 2.1.5 Creation, Destruction and Logging of FMU Instances

Optional call back function to signal if the computation of a communication step of a co-simulation slave is finished. A null pointer can be provided. In this case the master must use fmiGetStatus(..) to query the status of fmi2DoStep. If a pointer to a function is provided, it must be called by the FMU after a completed communication step.
"""
function fmi2CallbackStepFinished(componentEnvironment::Ptr{Cvoid}, status::Cuint)
    @debug "cbStepFinished(_, $(status)): Step finished."
    nothing
end

# Common function for ModelExchange & CoSimulation

"""
Source: FMISpec2.0.2[p.22]: 2.1.5 Creation, Destruction and Logging of FMU Instances

Disposes the given instance, unloads the loaded model, and frees all the allocated memory and other resources that have been allocated by the functions of the FMU interface.
If a null pointer is provided for “c”, the function call is ignored (does not have an effect).

Removes the component from the FMUs component list.
"""
function fmi2FreeInstance!(c::FMU2Component; popComponent::Bool = true)

    if popComponent
        ind = findall(x -> x.compAddr==c.compAddr, c.fmu.components)
        @assert length(ind) == 1 "fmi2FreeInstance!(...): Freeing $(length(ind)) instances with one call, this is not allowed."
        deleteat!(c.fmu.components, ind)
    end

    fmi2FreeInstance!(c.fmu.cFreeInstance, c.compAddr)

    nothing
end

"""
Source: FMISpec2.0.2[p.22]: 2.1.4 Inquire Platform and Version Number of Header Files

Returns the string to uniquely identify the “fmi2TypesPlatform.h” header file used for compilation of the functions of the FMU.
The standard header file, as documented in this specification, has fmi2TypesPlatform set to “default” (so this function usually returns “default”).
"""
function fmi2GetTypesPlatform(fmu::FMU2)

    typesPlatform = fmi2GetTypesPlatform(fmu.cGetTypesPlatform)

    unsafe_string(typesPlatform)
end
# special case
function fmi2GetTypesPlatform(c::FMU2Component)
    fmi2GetTypesPlatform(c.fmu)
end

"""
Source: FMISpec2.0.2[p.22]: 2.1.4 Inquire Platform and Version Number of Header Files

Returns the version of the “fmi2Functions.h” header file which was used to compile the functions of the FMU. The function returns “fmiVersion” which is defined in this header file. The standard header file as documented in this specification has version “2.0”
"""
function fmi2GetVersion(fmu::FMU2)

    fmi2Version = fmi2GetVersion(fmu.cGetVersion)

    unsafe_string(fmi2Version)
end
# special case
function fmi2GetVersion(c::FMU2Component)
    fmi2GetVersion(c.fmu)
end

# helper
function checkStatus(c::FMU2Component, status::fmi2Status)
    @assert (status != fmi2StatusWarning) || !c.fmu.executionConfig.assertOnWarning "Assert on `fmi2StatusWarning`. See stack for errors."

    if status == fmi2StatusError
        c.state = fmi2ComponentStateError
        @assert !c.fmu.executionConfig.assertOnError "Assert on `fmi2StatusError`. See stack for errors."

    elseif status == fmi2StatusFatal
        c.state = fmi2ComponentStateFatal
        @assert false "Assert on `fmi2StatusFatal`. See stack for errors."
    end
end

"""
Source: FMISpec2.0.2[p.22]: 2.1.5 Creation, Destruction and Logging of FMU Instances

The function controls debug logging that is output via the logger function callback. If loggingOn = fmi2True, debug logging is enabled, otherwise it is switched off.
"""
function fmi2SetDebugLogging(c::FMU2Component, logginOn::fmi2Boolean, nCategories::Unsigned, categories::Ptr{Nothing})

    status = fmi2SetDebugLogging(c.fmu.cSetDebugLogging, c.compAddr, logginOn, nCategories, categories)
    checkStatus(c, status)
    return status
end

"""
Source: FMISpec2.0.2[p.23]: 2.1.6 Initialization, Termination, and Resetting an FMU

Informs the FMU to setup the experiment. This function must be called after fmi2Instantiate and before fmi2EnterInitializationMode is called.The function controls debug logging that is output via the logger function callback. If loggingOn = fmi2True, debug logging is enabled, otherwise it is switched off.
"""
function fmi2SetupExperiment(c::FMU2Component,
    toleranceDefined::fmi2Boolean,
    tolerance::fmi2Real,
    startTime::fmi2Real,
    stopTimeDefined::fmi2Boolean,
    stopTime::fmi2Real)

    if c.state != fmi2ComponentStateInstantiated
        @warn "fmi2SetupExperiment(...): Needs to be called in state `fmi2ComponentStateInstantiated`."
    end

    if startTime != 0.0
        if c.fmu.executionConfig.autoTimeShift
            #@info "fmi2SetupExperiment(...): You picked a start time which is not zero. Many FMUs don't support that feature, so all time intervals are shifted in the background to achieve the desired result. If this feature is unwanted, please use `myFMU.executionConfig.autoTimeShift=false`."
            c.t_offset = -startTime
            stopTime -= startTime
            startTime = 0.0 # equivalent to: startTime -= startTime
        end
    end

    status = fmi2SetupExperiment(c.fmu.cSetupExperiment,
                c.compAddr, toleranceDefined, tolerance, startTime, stopTimeDefined, stopTime)
    checkStatus(c, status)

    # remain in status on success, nothing to do here

    return status
end

"""
Source: FMISpec2.0.2[p.23]: 2.1.6 Initialization, Termination, and Resetting an FMU

Informs the FMU to enter Initialization Mode. Before calling this function, all variables with attribute <ScalarVariable initial = "exact" or "approx"> can be set with the “fmi2SetXXX” functions (the ScalarVariable attributes are defined in the Model Description File, see section 2.2.7). Setting other variables is not allowed. Furthermore, fmi2SetupExperiment must be called at least once before calling fmi2EnterInitializationMode, in order that startTime is defined.
"""
function fmi2EnterInitializationMode(c::FMU2Component)

    if c.state != fmi2ComponentStateInstantiated
        @warn "fmi2EnterInitializationMode(...): Needs to be called in state `fmi2ComponentStateInstantiated`."
    end
    status = fmi2EnterInitializationMode(c.fmu.cEnterInitializationMode, c.compAddr)
    checkStatus(c, status)
    if status == fmi2StatusOK
        c.state = fmi2ComponentStateInitializationMode
    end
    return status
end

"""
Source: FMISpec2.0.2[p.23]: 2.1.6 Initialization, Termination, and Resetting an FMU

Informs the FMU to exit Initialization Mode.
"""
function fmi2ExitInitializationMode(c::FMU2Component)

    if c.state != fmi2ComponentStateInitializationMode
        @warn "fmi2ExitInitializationMode(...): Needs to be called in state `fmi2ComponentStateInitializationMode`."
    end

    status = fmi2ExitInitializationMode(c.fmu.cExitInitializationMode, c.compAddr)
    checkStatus(c, status)
    if status == fmi2StatusOK
        c.state = fmi2ComponentStateEventMode
    end
    return status
end

"""
Source: FMISpec2.0.2[p.24]: 2.1.6 Initialization, Termination, and Resetting an FMU

Informs the FMU that the simulation run is terminated.
"""
function fmi2Terminate(c::FMU2Component; soft::Bool=false)
    if c.state != fmi2ComponentStateContinuousTimeMode && c.state != fmi2ComponentStateEventMode
        if soft
            return fmi2StatusOK
        else
            @warn "fmi2Terminate(_): Needs to be called in state `fmi2ComponentStateContinuousTimeMode` or `fmi2ComponentStateEventMode`."
        end
    end

    status = fmi2Terminate(c.fmu.cTerminate, c.compAddr)
    checkStatus(c, status)
    if status == fmi2StatusOK
        c.state = fmi2ComponentStateTerminated
    end
    return status
end

"""
Source: FMISpec2.0.2[p.24]: 2.1.6 Initialization, Termination, and Resetting an FMU

Is called by the environment to reset the FMU after a simulation run. The FMU goes into the same state as if fmi2Instantiate would have been called.
"""
function fmi2Reset(c::FMU2Component; soft::Bool=false)
    if c.state != fmi2ComponentStateTerminated && c.state != fmi2ComponentStateError
        if soft
            return fmi2StatusOK
        else
            @warn "fmi2Reset(_): Needs to be called in state `fmi2ComponentStateTerminated` or `fmi2ComponentStateError`."
        end
    end

    if c.fmu.cReset == C_NULL
        fmi2FreeInstance!(c.fmu.cFreeInstance, c.compAddr)
        compAddr = fmi2Instantiate(c.fmu.cInstantiate, pointer(c.fmu.instanceName), c.fmu.type, pointer(c.fmu.modelDescription.guid), pointer(c.fmu.fmuResourceLocation), Ptr{fmi2CallbackFunctions}(pointer_from_objref(c.callbackFunctions)), fmi2Boolean(false), fmi2Boolean(false))

        if compAddr == Ptr{Cvoid}(C_NULL)
            @error "fmi2Reset(...): Reinstantiation failed!"
            return fmi2StatusError
        end

        c.compAddr = compAddr
        return fmi2StatusOK
    else
        status = fmi2Reset(c.fmu.cReset, c.compAddr)
        checkStatus(c, status)
        if status == fmi2StatusOK
            c.state = fmi2ComponentStateInstantiated
        end
        return status
    end
end

"""
Source: FMISpec2.0.2[p.24]: 2.1.7 Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference
"""
function fmi2GetReal!(c::FMU2Component, vr::AbstractArray{fmi2ValueReference}, nvr::Csize_t, value::AbstractArray{fmi2Real})

    status = fmi2GetReal!(c.fmu.cGetReal,
          c.compAddr, vr, nvr, value)
    checkStatus(c, status)
    return status
end

"""
Source: FMISpec2.0.2[p.24]: 2.1.7 Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference
"""
function fmi2SetReal(c::FMU2Component, vr::AbstractArray{fmi2ValueReference}, nvr::Csize_t, value::AbstractArray{fmi2Real})

    status = fmi2SetReal(c.fmu.cSetReal,
                c.compAddr, vr, nvr, value)
    checkStatus(c, status)
    return status
end

"""
Source: FMISpec2.0.2[p.24]: 2.1.7 Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference
"""
function fmi2GetInteger!(c::FMU2Component, vr::AbstractArray{fmi2ValueReference}, nvr::Csize_t, value::AbstractArray{fmi2Integer})

    status = fmi2GetInteger!(c.fmu.cGetInteger,
                c.compAddr, vr, nvr, value)
    checkStatus(c, status)
    return status
end

"""
Source: FMISpec2.0.2[p.24]: 2.1.7 Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference
"""
function fmi2SetInteger(c::FMU2Component, vr::AbstractArray{fmi2ValueReference}, nvr::Csize_t, value::AbstractArray{fmi2Integer})

    status = fmi2SetInteger(c.fmu.cSetInteger,
                c.compAddr, vr, nvr, value)
    checkStatus(c, status)
    return status
end

"""
Source: FMISpec2.0.2[p.24]: 2.1.7 Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference
"""
function fmi2GetBoolean!(c::FMU2Component, vr::AbstractArray{fmi2ValueReference}, nvr::Csize_t, value::AbstractArray{fmi2Boolean})

    status = fmi2GetBoolean!(c.fmu.cGetBoolean,
                c.compAddr, vr, nvr, value)
    checkStatus(c, status)
    return status
end

"""
Source: FMISpec2.0.2[p.24]: 2.1.7 Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference
"""
function fmi2SetBoolean(c::FMU2Component, vr::AbstractArray{fmi2ValueReference}, nvr::Csize_t, value::AbstractArray{fmi2Boolean})

    status = fmi2SetBoolean(c.fmu.cSetBoolean,
                c.compAddr, vr, nvr, value)
    checkStatus(c, status)
    return status
end

"""
Source: FMISpec2.0.2[p.24]: 2.1.7 Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference
"""
function fmi2GetString!(c::FMU2Component, vr::AbstractArray{fmi2ValueReference}, nvr::Csize_t, value::Union{AbstractArray{Ptr{Cchar}}, AbstractArray{Ptr{UInt8}}})

    status = fmi2GetString!(c.fmu.cGetString,
                c.compAddr, vr, nvr, value)
    checkStatus(c, status)
    return status
end

"""
Source: FMISpec2.0.2[p.24]: 2.1.7 Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference
"""
function fmi2SetString(c::FMU2Component, vr::AbstractArray{fmi2ValueReference}, nvr::Csize_t, value::Union{AbstractArray{Ptr{Cchar}}, AbstractArray{Ptr{UInt8}}})

    status = fmi2SetString(c.fmu.cSetString,
                c.compAddr, vr, nvr, value)
    checkStatus(c, status)
    return status
end

"""
Source: FMISpec2.0.2[p.26]: 2.1.8 Getting and Setting the Complete FMU State

fmi2GetFMUstate makes a copy of the internal FMU state and returns a pointer to this copy
"""
function fmi2GetFMUstate!(c::FMU2Component, FMUstate::Ref{fmi2FMUstate})

    status = fmi2GetFMUstate!(c.fmu.cGetFMUstate,
                c.compAddr, FMUstate)
    checkStatus(c, status)
    return status
end

"""
Source: FMISpec2.0.2[p.26]: 2.1.8 Getting and Setting the Complete FMU State

fmi2SetFMUstate copies the content of the previously copied FMUstate back and uses it as actual new FMU state.
"""
function fmi2SetFMUstate(c::FMU2Component, FMUstate::fmi2FMUstate)

    status = fmi2SetFMUstate(c.fmu.cSetFMUstate,
                c.compAddr, FMUstate)
    checkStatus(c, status)
    return status
end

"""
Source: FMISpec2.0.2[p.26]: 2.1.8 Getting and Setting the Complete FMU State

fmi2FreeFMUstate frees all memory and other resources allocated with the fmi2GetFMUstate call for this FMUstate.
"""
function fmi2FreeFMUstate!(c::FMU2Component, FMUstate::Ref{fmi2FMUstate})

    status = fmi2FreeFMUstate!(c.fmu.cFreeFMUstate,
                c.compAddr, FMUstate)
    checkStatus(c, status)
    return status
end

"""
Source: FMISpec2.0.2[p.26]: 2.1.8 Getting and Setting the Complete FMU State

fmi2SerializedFMUstateSize returns the size of the byte vector, in order that FMUstate can be stored in it.
"""
function fmi2SerializedFMUstateSize!(c::FMU2Component, FMUstate::fmi2FMUstate, size::Ref{Csize_t})

    status = fmi2SerializedFMUstateSize!(c.fmu.cSerializedFMUstateSize,
                c.compAddr, FMUstate, size)
    checkStatus(c, status)
    return status
end

"""
Source: FMISpec2.0.2[p.26]: 2.1.8 Getting and Setting the Complete FMU State

fmi2SerializeFMUstate serializes the data which is referenced by pointer FMUstate and copies this data in to the byte vector serializedState of length size
"""
function fmi2SerializeFMUstate!(c::FMU2Component, FMUstate::fmi2FMUstate, serialzedState::AbstractArray{fmi2Byte}, size::Csize_t)

    status = fmi2SerializeFMUstate!(c.fmu.cSerializeFMUstate,
                c.compAddr, FMUstate, serialzedState, size)
    checkStatus(c, status)
    return status
end

"""
Source: FMISpec2.0.2[p.26]: 2.1.8 Getting and Setting the Complete FMU State

fmi2DeSerializeFMUstate deserializes the byte vector serializedState of length size, constructs a copy of the FMU state and returns FMUstate, the pointer to this copy.
"""
function fmi2DeSerializeFMUstate!(c::FMU2Component, serializedState::AbstractArray{fmi2Byte}, size::Csize_t, FMUstate::Ref{fmi2FMUstate})

    status = fmi2DeSerializeFMUstate!(c.fmu.cDeSerializeFMUstate,
                c.compAddr, serializedState, size, FMUstate)
    checkStatus(c, status)
    return status
end

"""
Source: FMISpec2.0.2[p.26]: 2.1.9 Getting Partial Derivatives

This function computes the directional derivatives of an FMU.
"""
function fmi2GetDirectionalDerivative!(c::FMU2Component,
                                       vUnknown_ref::AbstractArray{fmi2ValueReference},
                                       nUnknown::Csize_t,
                                       vKnown_ref::AbstractArray{fmi2ValueReference},
                                       nKnown::Csize_t,
                                       dvKnown::AbstractArray{fmi2Real},
                                       dvUnknown::AbstractArray) # ToDo: Datatype for AbstractArray
    @assert fmi2ProvidesDirectionalDerivative(c.fmu) ["fmi2GetDirectionalDerivative!(...): This FMU does not support build-in directional derivatives!"]

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

# Functions specificly for isCoSimulation
"""
Source: FMISpec2.0.2[p.104]: 4.2.1 Transfer of Input / Output Values and Parameters

Sets the n-th time derivative of real input variables.
vr defines the value references of the variables
the array order specifies the corresponding order of derivation of the variables
"""
function fmi2SetRealInputDerivatives(c::FMU2Component, vr::AbstractArray{fmi2ValueReference}, nvr::Csize_t, order::AbstractArray{fmi2Integer}, value::AbstractArray{fmi2Real})

    status = fmi2SetRealInputDerivatives(c.fmu.cSetRealInputDerivatives,
                c.compAddr, vr, nvr, order, value)
    checkStatus(c, status)
    return status
end

"""
Source: FMISpec2.0.2[p.104]: 4.2.1 Transfer of Input / Output Values and Parameters

Retrieves the n-th derivative of output values.
vr defines the value references of the variables
the array order specifies the corresponding order of derivation of the variables
"""
function fmi2GetRealOutputDerivatives!(c::FMU2Component,  vr::AbstractArray{fmi2ValueReference}, nvr::Csize_t, order::AbstractArray{fmi2Integer}, value::AbstractArray{fmi2Real})

    status = fmi2GetRealOutputDerivatives!(c.fmu.cGetRealOutputDerivatives,
                c.compAddr, vr, nvr, order, value)
    checkStatus(c, status)
    return status
end

"""
Source: FMISpec2.0.2[p.104]: 4.2.2 Computation

The computation of a time step is started.
"""
function fmi2DoStep(c::FMU2Component, currentCommunicationPoint::fmi2Real, communicationStepSize::fmi2Real, noSetFMUStatePriorToCurrentPoint::fmi2Boolean)
    @assert c.fmu.cDoStep != C_NULL ["fmi2DoStep(...): This FMU does not support fmi2DoStep, probably it's a ME-FMU with no CS-support?"]

    status = fmi2DoStep(c.fmu.cDoStep,
          c.compAddr, currentCommunicationPoint, communicationStepSize, noSetFMUStatePriorToCurrentPoint)
    checkStatus(c, status)
    return status
end

"""
Source: FMISpec2.0.2[p.105]: 4.2.2 Computation

Can be called if fmi2DoStep returned fmi2Pending in order to stop the current asynchronous execution.
"""
function fmi2CancelStep(c::FMU2Component)

    status = fmi2CancelStep(c.fmu.cCancelStep, c.compAddr)
    checkStatus(c, status)
    return status
end

"""
Source: FMISpec2.0.2[p.106]: 4.2.3 Retrieving Status Information from the Slave

Informs the master about the actual status of the simulation run. Which status information is to be returned is specified by the argument fmi2StatusKind.
"""
function fmi2GetStatus!(c::FMU2Component, s::fmi2StatusKind, value)
    rtype = nothing
    if s == fmi2Terminated
        rtype = fmi2Boolean
    else
        @assert false "fmi2GetStatus!(_, $(s), $(value)): StatusKind $(s) not implemented yet, please open an issue."
    end
    @assert typeof(value) == rtype "fmi2GetStatus!(_, $(s), $(value)): Type of value ($(typeof(value))) doesn't fit type of return type $(rtype). Change type of value to $(rtype) or change status kind."

    status = fmi2Error
    if rtype == fmi2Boolean
        status = fmi2GetStatus!(c.fmu.cGetRealStatus,
                    c.compAddr, s, Ref(value))
        checkStatus(c, status)
    end
    return status
end

"""
Source: FMISpec2.0.2[p.106]: 4.2.3 Retrieving Status Information from the Slave

Informs the master about the actual status of the simulation run. Which status information is to be returned is specified by the argument fmi2StatusKind.
"""
function fmi2GetRealStatus!(c::FMU2Component, s::fmi2StatusKind, value::fmi2Real)

    status = fmi2GetRealStatus!(c.fmu.cGetRealStatus,
                c.compAddr, s, Ref(value))
    checkStatus(c, status)
    return status
end

"""
Source: FMISpec2.0.2[p.106]: 4.2.3 Retrieving Status Information from the Slave

Informs the master about the actual status of the simulation run. Which status information is to be returned is specified by the argument fmi2StatusKind.
"""
function fmi2GetIntegerStatus!(c::FMU2Component, s::fmi2StatusKind, value::fmi2Integer)

    status = fmi2GetIntegerStatus!(c.fmu.cGetIntegerStatus,
                c.compAddr, s, Ref(value))
    checkStatus(c, status)
    return status
end

"""
Source: FMISpec2.0.2[p.106]: 4.2.3 Retrieving Status Information from the Slave

Informs the master about the actual status of the simulation run. Which status information is to be returned is specified by the argument fmi2StatusKind.
"""
function fmi2GetBooleanStatus!(c::FMU2Component, s::fmi2StatusKind, value::fmi2Boolean)

    status = fmi2GetBooleanStatus!(c.fmu.cGetBooleanStatus,
                c.compAddr, s, Ref(value))
    checkStatus(c, status)
    return status
end

"""
Source: FMISpec2.0.2[p.106]: 4.2.3 Retrieving Status Information from the Slave

Informs the master about the actual status of the simulation run. Which status information is to be returned is specified by the argument fmi2StatusKind.
"""
function fmi2GetStringStatus!(c::FMU2Component, s::fmi2StatusKind, value::fmi2String)

    status = fmi2GetStringStatus!(c.fmu.cGetStringStatus,
                c.compAddr, s, Ref(value))
    checkStatus(c, status)
    return status
end

# Model Exchange specific Functions

"""
Source: FMISpec2.0.2[p.83]: 3.2.1 Providing Independent Variables and Re-initialization of Caching

Set a new time instant and re-initialize caching of variables that depend on time, provided the newly provided time value is different to the previously set time value (variables that depend solely on constants or parameters need not to be newly computed in the sequel, but the previously computed values can be reused).
"""
function fmi2SetTime(c::FMU2Component, time::fmi2Real)

    status = fmi2SetTime(c.fmu.cSetTime,
          c.compAddr, time + c.t_offset)
    checkStatus(c, status)
    if status == fmi2StatusOK
        c.t = time
    end
    return status
end

"""
Source: FMISpec2.0.2[p.83]: 3.2.1 Providing Independent Variables and Re-initialization of Caching

Set a new (continuous) state vector and re-initialize caching of variables that depend on the states. Argument nx is the length of vector x and is provided for checking purposes
"""
function fmi2SetContinuousStates(c::FMU2Component,
                                 x::AbstractArray{fmi2Real},
                                 nx::Csize_t)

    status = fmi2SetContinuousStates(c.fmu.cSetContinuousStates, c.compAddr, x, nx)
    checkStatus(c, status)
    return status
end

"""
Source: FMISpec2.0.2[p.84]: 3.2.2 Evaluation of Model Equations

The model enters Event Mode from the Continuous-Time Mode and discrete-time equations may become active (and relations are not “frozen”).
"""
function fmi2EnterEventMode(c::FMU2Component; soft::Bool=false)

    if c.state != fmi2ComponentStateContinuousTimeMode
        if soft
            return fmi2StatusOK
        else
            @warn "fmi2EnterEventMode(...): Called at the wrong time."
        end
    end

    status = fmi2EnterEventMode(c.fmu.cEnterEventMode,
          c.compAddr)
    checkStatus(c, status)
    if status == fmi2StatusOK
        c.state = fmi2ComponentStateEventMode
    end
    return status
end

"""
Source: FMISpec2.0.2[p.84]: 3.2.2 Evaluation of Model Equations

The FMU is in Event Mode and the super dense time is incremented by this call.
"""
function fmi2NewDiscreteStates!(c::FMU2Component, eventInfo::fmi2EventInfo)

    if c.state != fmi2ComponentStateEventMode
        @warn "fmi2NewDiscreteStates(...): Needs to be called in state `fmi2ComponentStateEventMode`."
    end

    status = fmi2NewDiscreteStates!(c.fmu.cNewDiscreteStates,
                    c.compAddr, Ptr{fmi2EventInfo}(pointer_from_objref(eventInfo)) )

    if eventInfo.nextEventTimeDefined == fmi2True
        eventInfo.nextEventTime -= c.t_offset
    end

    checkStatus(c, status)
    # remain in the same mode and status (or ToDo: Meta-states)
    return status
end

"""
Source: FMISpec2.0.2[p.85]: 3.2.2 Evaluation of Model Equations

The model enters Continuous-Time Mode and all discrete-time equations become inactive and all relations are “frozen”.
This function has to be called when changing from Event Mode (after the global event iteration in Event Mode over all involved FMUs and other models has converged) into Continuous-Time Mode.
"""
function fmi2EnterContinuousTimeMode(c::FMU2Component; soft::Bool=false)

    if c.state != fmi2ComponentStateEventMode
        if soft
            return fmi2StatusOK
        else
            @warn "fmi2EnterContinuousTimeMode(...): Needs to be called in state `fmi2ComponentStateEventMode`."
        end
    end

    status = fmi2EnterContinuousTimeMode(c.fmu.cEnterContinuousTimeMode,
          c.compAddr)
    checkStatus(c, status)
    if status == fmi2StatusOK
        c.state = fmi2ComponentStateContinuousTimeMode
    end
    return status
end

"""
Source: FMISpec2.0.2[p.85]: 3.2.2 Evaluation of Model Equations

This function must be called by the environment after every completed step of the integrator provided the capability flag completedIntegratorStepNotNeeded = false.
If enterEventMode == fmi2True, the event mode must be entered
If terminateSimulation == fmi2True, the simulation shall be terminated
"""
function fmi2CompletedIntegratorStep!(c::FMU2Component,
                                      noSetFMUStatePriorToCurrentPoint::fmi2Boolean,
                                      enterEventMode::fmi2Boolean,
                                      terminateSimulation::fmi2Boolean)

    status = fmi2CompletedIntegratorStep!(c.fmu.cCompletedIntegratorStep,
          c.compAddr, noSetFMUStatePriorToCurrentPoint, Ref(enterEventMode), Ref(terminateSimulation))
    checkStatus(c, status)
    return status
end

"""
Source: FMISpec2.0.2[p.86]: 3.2.2 Evaluation of Model Equations

Compute state derivatives at the current time instant and for the current states.
"""
function fmi2GetDerivatives!(c::FMU2Component,
                            derivatives::AbstractArray{fmi2Real},
                            nx::Csize_t)

    status = fmi2GetDerivatives!(c.fmu.cGetDerivatives,
          c.compAddr, derivatives, nx)
    checkStatus(c, status)
    return status
end

"""
Source: FMISpec2.0.2[p.86]: 3.2.2 Evaluation of Model Equations

Compute event indicators at the current time instant and for the current states.
"""
function fmi2GetEventIndicators!(c::FMU2Component, eventIndicators::AbstractArray{fmi2Real}, ni::Csize_t)

    status = fmi2GetEventIndicators!(c.fmu.cGetEventIndicators,
                    c.compAddr, eventIndicators, ni)
    checkStatus(c, status)
    return status
end

"""
Source: FMISpec2.0.2[p.86]: 3.2.2 Evaluation of Model Equations

Return the new (continuous) state vector x.
"""
function fmi2GetContinuousStates!(c::FMU2Component,
                                 x::AbstractArray{fmi2Real},
                                 nx::Csize_t)

    status = fmi2GetContinuousStates!(c.fmu.cGetContinuousStates,
          c.compAddr, x, nx)
    checkStatus(c, status)
    return status
end

"""
Source: FMISpec2.0.2[p.86]: 3.2.2 Evaluation of Model Equations

Return the nominal values of the continuous states.
"""
function fmi2GetNominalsOfContinuousStates!(c::FMU2Component, x_nominal::AbstractArray{fmi2Real}, nx::Csize_t)

    status = fmi2GetNominalsOfContinuousStates!(c.fmu.cGetNominalsOfContinuousStates,
                    c.compAddr, x_nominal, nx)
    checkStatus(c, status)
    return status
end

ext: jl


text:
function copy_hrsl()
   hrsl_path = "hrsl_ls_adjusted.tif"
   admin3_path = expanduser(joinpath("~", "data", "projects", "uganda2020", "outputs",
      "uga_canonical", "uga_canonical_names.shp"))
   copy_path = "ztesting.tif"
   admin_cnt = 3
   if ispath(copy_path)
      rm(copy_path)
   end
   id_raster_ds = ArchGDAL.read(hrsl_path) do hrsl_ds
      create_fraction(hrsl_ds, copy_path, admin_cnt)
   end
   ArchGDAL.destroy(id_raster_ds)
end

copy_hrsl()
@test ispath("ztesting.tif")


text:
export getNodalGradientMatrix

function getNodalGradientMatrix(mesh::TetraMesh)
# G = getGradientMatrix(mesh)

if isempty(mesh.Grad)

edges = mesh.edges
# number of edges
ne = size(edges, 1)
nn = maximum(edges[:])

# assemble gradient as sparse matrix
i =  [1:ne
      1:ne]

j =  [edges[:, 1]
      edges[:, 2]]

g =   [ ones(ne)
       -ones(ne)]

L = getEdgeLength(mesh)

for k = 1:2*ne
   g[k] *= 1/L[i[k]]    #gradient scaled using edge lengths
end

mesh.Grad = sparse(i, j, g, ne, nn)

end

return mesh.Grad

end

max_stars_repo_name: JuliaInv/JTetra

max_issues_repo_name: JuliaInv/JTetra

max_forks_repo_name: JuliaInv/JTetra


text:
############################################################################################
"""
$(TYPEDEF)

Utility struct to help assign boundaries to parameter - keeps parameter fixed. Useful for assigning buffer values for functions of parameter.

# Fields
$(TYPEDFIELDS)
"""
struct Fixed <: AbstractConstraint
    function Fixed()
        return new()
    end
end

############################################################################################
function construct_transform(info::Fixed, val)
    return Bijectors.Identity{0}(), Bijectors.Identity{0}()
end

############################################################################################
function _check(
    _rng::Random.AbstractRNG,
    constraint::Fixed,
    val::Union{R,Array{R},AbstractArray},
) where {R<:Real}
    return true
end

############################################################################################
# Fixed constraint
function _construct_flatten(
    output::Type{T},
    flattentype::F,
    unflattentype::U,
    x
) where {T<:Real,F<:FlattenTypes,U<:UnflattenTypes}
    v = T[]
    _flatten_Fixed(x) = v
    _unflatten_Fixed(v) = x
    return _flatten_Fixed, _unflatten_Fixed
end
function construct_flatten(
    output::Type{T},
    flattentype::F,
    unflattentype::U,
    constraint::Fixed,
    x,
) where {T<:AbstractFloat,F<:FlattenTypes,U<:UnflattenTypes}
    return _construct_flatten(T, flattentype, unflattentype, x)
end

############################################################################################
#=
!NOTE: For this constraint, we use a Bijector as Transformer, so we do not need to add a new functors
1. MyTransformer <: AbstractTransformer, MyInverseTransformer <: AbstractTransformer
2. define a function construct_transform(MyConstraint, val) -> MyTransformer, MyInverseTransformer
3. overload unconstrain and log_abs_det_jac on MyTransformer.
4. overload constrain on MyInverseTransformer.
=#

############################################################################################
#Export
export
    Fixed,
    construct_flatten,
    construct_transform,
    constrain,
    unconstrain,
    log_abs_det_jac


text:
export CoreTracker, CoreTrackerResult, CoreTrackerStatus, CoreTrackerOptions,
        coretracker, coretracker_startsolutions, affine_tracking,
        track, track!, setup!, iterator,
        current_x, current_t, current_Δt, iters, status, accuracy, max_corrector_iters,
        max_step_size , refinement_accuracy, max_refinement_iters,
        set_accuracy!, set_max_corrector_iters!, set_refinement_accuracy!,
        set_max_refinement_iters!, set_max_step_size!, digits_lost, options

const coretracker_supported_keywords = [:corrector, :predictor, :patch,
    :initial_step_size, :min_step_size , :max_step_size,
    :accuracy, :refinement_accuracy, :max_corrector_iters, :max_refinement_iters,
    :max_steps, :simple_step_size_alg, :auto_scaling, :terminate_ill_conditioned,
    :log_transform, :precision, :steps_jacobian_info_update,
    :max_lost_digits]

####################
# CoreTrackerState #
####################

module CoreTrackerStatus

    @doc """
        CoreTrackerStatus.states

    The possible states the coretracker can achieve are

    * `CoreTrackerStatus.success`
    * `CoreTrackerStatus.tracking`
    * `CoreTrackerStatus.terminated_maximal_iterations`
    * `CoreTrackerStatus.terminated_invalid_startvalue`
    * `CoreTrackerStatus.terminated_step_size_too_small`
    * `CoreTrackerStatus.terminated_singularity`
    * `CoreTrackerStatus.terminated_ill_conditioned`
    """
    @enum states begin
        success
        tracking
        terminated_maximal_iterations
        terminated_invalid_startvalue
        terminated_step_size_too_small
        terminated_singularity
        terminated_ill_conditioned
    end
end

####################
# CoreTrackerResult #
####################

"""
     CoreTrackerResult{V<:AbstractVector}

Containing the result of a tracked path. The fields are
* `returncode::CoreTrackerStatus.states` If the tracking was successfull then it is `CoreTrackerStatus.success`.
* `x::V` The result.
* `t::ComplexF64` The `t` when the path tracker stopped.
* `accuracy::Float64`: The estimated accuracy of `x`.
"""
struct CoreTrackerResult{V <: AbstractVector}
     returncode::CoreTrackerStatus.states
     x::V
     t::ComplexF64
     accuracy::Float64
     accepted_steps::Int
     rejected_steps::Int
end

function CoreTrackerResult(tracker)
    state = tracker.state
     CoreTrackerResult(state.status,
          copy(state.x), state.segment[state.s],
          state.accuracy,
          state.accepted_steps,
          state.rejected_steps)
end

Base.show(io::IO, result::CoreTrackerResult) = print_fieldnames(io, result)
Base.show(io::IO, ::MIME"application/prs.juno.inline", result::CoreTrackerResult) = result

##################
## PATH TRACKER ##
##################

"""
    AutoScalingOptions(;scale_min=0.01,
                     scale_abs_min=min(scale_min^2, 200 * sqrt(eps()),
                     scale_max=1.0 / eps() / sqrt(2)

Parameters for the auto scaling of the variables.
"""
struct AutoScalingOptions
    scale_min::Float64
    scale_abs_min::Float64
    scale_max::Float64
end
function AutoScalingOptions(;scale_min=0.01,
                             scale_abs_min=min(scale_min^2, 200 * sqrt(eps())),
                             scale_max=1.0 / eps() / sqrt(2))
    AutoScalingOptions(scale_min, scale_abs_min, scale_max)
end
Base.show(io::IO, opts::AutoScalingOptions) = print_fieldnames(io, opts)
Base.show(io::IO, ::MIME"application/prs.juno.inline", opts::AutoScalingOptions) = opts

mutable struct StepSizeModel
    ω::Float64
    expected_Δx₀::Float64
end
StepSizeModel() = StepSizeModel(NaN, NaN)

function reset!(model::StepSizeModel)
    model.ω = model.expected_Δx₀ = NaN
    model
end
Base.show(io::IO, opts::StepSizeModel) = print_fieldnames(io, opts)
Base.show(io::IO, ::MIME"application/prs.juno.inline", opts::StepSizeModel) = opts

######################
# CoreTrackerOptions #
######################

"""
    CoreTrackerOptions

The set of options set for a [`CoreTracker`](@ref). See the description of [`CoreTracker`](@ref)
for all possible options.
"""
mutable struct CoreTrackerOptions
    accuracy::Float64
    max_corrector_iters::Int
    refinement_accuracy::Float64
    max_refinement_iters::Int
    max_steps::Int
    initial_step_size::Float64
    min_step_size ::Float64
    max_step_size ::Float64
    simple_step_size_alg::Bool
    update_patch::Bool
    max_lost_digits::Float64
    auto_scaling::Bool
    auto_scaling_options::AutoScalingOptions
    terminate_ill_conditioned::Bool
    precision::PrecisionOption
    steps_jacobian_info_update::Int
end

function CoreTrackerOptions(;
    accuracy=1e-7,
    auto_scaling=true,
    auto_scaling_options=AutoScalingOptions(),
    initial_step_size=0.1,
    max_corrector_iters::Int=2,
    max_refinement_iters=5,
    max_step_size=Inf,
    min_step_size=1e-14,
    parameter_homotopy=false,
    max_steps=parameter_homotopy ? 10_000 : 1_000,
    precision::PrecisionOption=PRECISION_FIXED_64,
    max_lost_digits=default_max_lost_digits(precision, accuracy),
    refinement_accuracy=1e-8,
    simple_step_size_alg=false,
    steps_jacobian_info_update::Int=3,
    terminate_ill_conditioned::Bool=true,
    update_patch=true)

    CoreTrackerOptions(accuracy, max_corrector_iters, refinement_accuracy,
            max_refinement_iters, max_steps, initial_step_size, min_step_size,
            max_step_size, simple_step_size_alg, update_patch, max_lost_digits,
            auto_scaling, auto_scaling_options, terminate_ill_conditioned,
            precision, steps_jacobian_info_update)
end

Base.show(io::IO, opts::CoreTrackerOptions) = print_fieldnames(io, opts)
Base.show(io::IO, ::MIME"application/prs.juno.inline", opts::CoreTrackerOptions) = opts

const eps_double_64 = Float64(eps(Double64))
function default_max_lost_digits(prec::PrecisionOption, accuracy::Float64)
    if prec == PRECISION_FIXED_64
        # maximal_digits_available - digits necessary - buffer
        -log10(eps()) + log10(accuracy) + 1
    else
        # if higher precision is available we will more like be constrained
        # by the fact that the jacobian cannot be too ill-conditioned
        min(-log10(eps()) - 3, -log10(eps_double_64) + log10(accuracy) - 3)
    end
end

####################
# CoreTrackerState #
####################

mutable struct CoreTrackerState{T, AV<:AbstractVector{T}, MaybePatchState <: Union{AbstractAffinePatchState, Nothing}, IP}
    x::AV # current x
    x̂::AV # last prediction
    x̄::AV # canidate for new x
    ẋ::Vector{T} # derivative at current x
    inner_product::IP
    step_size::StepSizeModel
    segment::ComplexSegment
    s::Float64 # current step length (0 ≤ s ≤ length(segment))
    Δs::Float64 # current step size
    Δs_prev::Float64 # previous step size
    accuracy::Float64
    jacobian::Jacobian{T}
    steps_jacobian_info_update::Int
    status::CoreTrackerStatus.states
    patch::MaybePatchState
    accepted_steps::Int
    rejected_steps::Int
    last_step_failed::Bool
    consecutive_successfull_steps::Int
end

function CoreTrackerState(H, x₁::AbstractVector, t₁, t₀, options::CoreTrackerOptions,
                           patch::Union{Nothing, AbstractAffinePatchState}=nothing,
                           inner_product::AbstractInnerProduct=EuclideanIP())
    if isa(x₁, SVector)
        x = Vector(x₁)
    else
        x = copy(x₁)
    end
    x̂, x̄ = copy(x), copy(x)
    ẋ = Vector(x)
    step_size_model = StepSizeModel()
    segment = ComplexSegment(promote(complex(t₁), complex(t₀))...)
    s = 0.0
    Δs = convert(Float64, min(options.initial_step_size, length(segment), options.max_step_size ))
    Δs_prev = 0.0
    accuracy = 0.0
    Jac = Jacobian(jacobian(H, x, t₁))
    steps_jacobian_info_update = 0
    accepted_steps = rejected_steps = 0
    status = CoreTrackerStatus.tracking
    last_step_failed = false
    consecutive_successfull_steps = 0
    CoreTrackerState(x, x̂, x̄, ẋ, inner_product, step_size_model, segment, s, Δs, Δs_prev, accuracy,
        Jac, steps_jacobian_info_update, status, patch,
        accepted_steps, rejected_steps, last_step_failed, consecutive_successfull_steps)
end

Base.show(io::IO, state::CoreTrackerState) = print_fieldnames(io, state)
Base.show(io::IO, ::MIME"application/prs.juno.inline", state::CoreTrackerState) = state

function reset!(state::CoreTrackerState, x₁::AbstractVector, t₁, t₀, options::CoreTrackerOptions, setup_patch::Bool, loop::Bool)
    state.segment = ComplexSegment(promote(t₁, t₀)...)
    state.s = 0.0
    state.Δs = min(options.initial_step_size, length(state.segment), options.max_step_size )
    state.Δs_prev = 0.0
    state.status = CoreTrackerStatus.tracking
    embed!(state.x, x₁)
    setup_patch && state.patch !== nothing && setup!(state.patch, state.x)
    state.last_step_failed = false
    state.consecutive_successfull_steps = 0
    if !loop
        state.accuracy = 0.0
        state.accepted_steps = state.rejected_steps = 0
        reset!(state.step_size)

        if options.auto_scaling
            init_auto_scaling!(state.inner_product, state.x, options.auto_scaling_options)
        end
        reset!(state.jacobian)
        state.steps_jacobian_info_update = 0
    end

    state
end

embed!(x::ProjectiveVectors.PVector, y) = ProjectiveVectors.embed!(x, y)
embed!(x::AbstractVector, y) = x .= y

function init_auto_scaling!(ip::WeightedIP, x::AbstractVector, opts::AutoScalingOptions)
    point_norm = euclidean_norm(x)
    w = ip.weight
    for i in 1:length(w)
        wᵢ = abs(x[i])
        if wᵢ < opts.scale_min * point_norm
            wᵢ = opts.scale_min * point_norm
        elseif wᵢ > opts.scale_max * point_norm
            wᵢ = opts.scale_max * point_norm
        end
        w[i] = max(wᵢ, opts.scale_abs_min)
    end
    nothing
end
init_auto_scaling!(ip::EuclideanIP, x::AbstractVector, opts::AutoScalingOptions) = nothing

####################
# CoreTrackerCache #
####################
mutable struct CoreTrackerCache{H<:HomotopyWithCache, P<:AbstractPredictorCache,
             C<:AbstractCorrectorCache, T}
    homotopy::H
    predictor::P
    corrector::C
    out::Vector{T}
    r::Vector{T}
end
function CoreTrackerCache(H::HomotopyWithCache, predictor, corrector, state::CoreTrackerState)
    t = state.segment[state.s]
    pcache = cache(predictor, H, state.x, state.ẋ, t)
    ccache = cache(corrector, H, state.x, t)
    out = H(state.x, t)
    r = copy(out)
    CoreTrackerCache(H, pcache, ccache, out, r)
end

###############
# CoreTracker #
###############
"""
    CoreTracker(problem::AbstractProblem, x₁; kwargs...)

Construct a `CoreTracker` from the given `problem` to track elements of type `x₁`.
The path is tracked using a predictor-corrector scheme. The recommended methods to construct
a `CoreTracker` are [`coretracker`](@ref) and [`coretracker_startsolutions`](@ref).
Note that a `CoreTracker` is also a (mutable) iterator.

## Keyword arguments
* `accuracy=1e-7`: The accuracy required during the path tracking.
* `corrector::AbstractCorrector`: The corrector used during in the predictor-corrector scheme. The default is [`NewtonCorrector`](@ref).
* `initial_step_size=0.1`: The size of the first step.
* `max_corrector_iters=2`: The maximal number of correction steps in a single step.
* `max_lost_digits::Real`: The tracking is terminated if we estimate that we loose more than `max_lost_digits` in the linear algebra steps. This threshold depends on the `precision` argument.
* `max_refinement_iters=5`: The maximal number of correction steps used to refine the final value.
* `max_steps=1_000`: The maximal number of iterations the path tracker has available. Note that this changes to `10_000` for parameter homotopies.
* `max_step_size=Inf`: The maximal step size.
* `min_step_size=1e-14`: The minimal step size.
* `precision::PrecisionOption=PRECISION_FIXED_64`: The precision used for evaluating the residual in Newton's method.
* `predictor::AbstractPredictor`: The predictor used during in the predictor-corrector scheme. The default is [`Heun`](@ref)()`.
* `refinement_accuracy=1e-8`: The precision required for the final value.
* `simple_step_size_alg=false`: Use a more simple step size algorithm.
* `steps_jacobian_info_update::Int=1`: Every n-th step a linear system will be solved using a QR factorization to obtain an estimate for the condition number of the Jacobian.
* `terminate_ill_conditioned::Bool=true`: Indicates whether the path tracking should be terminated for ill-conditioned paths. A path is considerd ill-conditioned if the condition number of the Jacobian is larger than ≈1e14 or if it is larger than 1e`max_lost_digits`.
"""
struct CoreTracker{H<:AbstractHomotopy,
    Predictor<:AbstractPredictor,
    Corrector<:AbstractCorrector,
    MaybePatch<:Union{Nothing, AbstractAffinePatch},
    S<:CoreTrackerState,
    C<:CoreTrackerCache}
    # these are fixed
    homotopy::H
    predictor::Predictor
    corrector::Corrector
    affine_patch::MaybePatch
    # these are mutable
    state::S
    options::CoreTrackerOptions
    cache::C
end

function CoreTracker(prob::AbstractProblem, x₁; kwargs...)
    CoreTracker(prob.homotopy, embed(prob, x₁), complex(1.0), complex(0.0), prob; kwargs...)
end

# Tracking in Projective Space
function CoreTracker(homotopy::AbstractHomotopy, x₁::ProjectiveVectors.PVector, t₁, t₀, prob::AbstractProblem;
    patch=has_dedicated_homvars(prob.vargroups) ? EmbeddingPatch() : OrthogonalPatch(),
    corrector::AbstractCorrector=NewtonCorrector(),
    predictor::AbstractPredictor=default_predictor(x₁),
    log_transform=false,
    simple_step_size_alg=!isa(patch, EmbeddingPatch), kwargs...)

    options = CoreTrackerOptions(;
                    parameter_homotopy=isa(homotopy, ParameterHomotopy),
                    kwargs...)

    if homotopy isa PatchedHomotopy
        error("You cannot pass a `PatchedHomotopy` to CoreTracker. Instead pass the homotopy and patch separate.")
    end

    H = log_transform ? LogHomotopy(homotopy) : homotopy
    patch_state = state(patch, x₁)
    # We close over the patch state, the homotopy and its cache
    # to be able to pass things around more easily
    HC = HomotopyWithCache(PatchedHomotopy(H, patch_state), x₁, t₁)
    if is_global_patch(patch) && options.auto_scaling
        inner_product = WeightedIP(x₁)
    else
        inner_product = EuclideanIP()
    end
    # We have to make sure that the element type of x is invariant under evaluation
    tracker_state = CoreTrackerState(HC, indempotent_x(HC, x₁, t₁), t₁, t₀, options, patch_state, inner_product)
    cache = CoreTrackerCache(HC, predictor, corrector, tracker_state)

    CoreTracker(H, predictor, corrector, patch, tracker_state, options, cache)
end

# Tracking in affine space
function CoreTracker(homotopy::AbstractHomotopy, x₁::AbstractVector, t₁, t₀, prob::AbstractProblem;
    patch=nothing,
    corrector::AbstractCorrector=NewtonCorrector(),
    predictor::AbstractPredictor=default_predictor(x₁),
    log_transform=false, kwargs...)

    if patch !== nothing
        throw(ArgumentError("You can only pass `patch=$(patch)` if `affine_tracking=false`."))
    end

    options = CoreTrackerOptions(;parameter_homotopy=isa(homotopy, ParameterHomotopy), kwargs...)

    H = log_transform ? LogHomotopy(homotopy) : homotopy
    # We close over the patch state, the homotopy and its cache
    # to be able to pass things around more easily
    HC = HomotopyWithCache(H, Vector(x₁), t₁)

    inner_product = options.auto_scaling ? WeightedIP(x₁) : EuclideanIP()
    # We have to make sure that the element type of x is invariant under evaluation
    tracker_state = CoreTrackerState(HC, indempotent_x(HC, x₁, t₁), t₁, t₀, options, nothing, inner_product)
    cache = CoreTrackerCache(HC, predictor, corrector, tracker_state)

    CoreTracker(H, predictor, corrector, nothing, tracker_state, options, cache)
end

default_predictor(x::AbstractVector) = Heun()
# Do not really understand this but Heun doesn't work that great for multi-homogeneous tracking
default_predictor(x::ProjectiveVectors.PVector{T,1}) where {T} = Heun()
default_predictor(x::ProjectiveVectors.PVector{T,N}) where {T, N} = Euler()

"""
    indempotent_x(H, x₁, t₁)

This returns a vector similar to `x₁` but with an element type which is invariant under evaluation.
"""
function indempotent_x(H, x₁, t₁)
    u = Vector{Any}(undef, size(H)[1])
    evaluate!(u, H, x₁, t₁)

    if isa(x₁, SVector)
        indem_x = Vector{promote_type(typeof(u[1]), ComplexF64)}(undef, length(x₁))
    else
        indem_x = similar(x₁, promote_type(typeof(u[1]), ComplexF64))
    end
    indem_x .= x₁
end

Base.show(io::IO, C::CoreTracker) = print(io, "CoreTracker tracking a path of type $(typeof(C.state.x))")
Base.show(io::IO, ::MIME"application/prs.juno.inline", x::CoreTracker) = x

##############
## TRACKING ##
##############

"""
    track(tracker, x₁, t₁=1.0, t₀=0.0; options...)::CoreTrackerResult

Track a value `x₁` from `t₁` to `t₀` using the given `CoreTracker` `tracker`.
This returns a `CoreTrackerResult`. This modifies `tracker`.
See [`track!`](@ref) for the possible options.
"""
function track(tracker::CoreTracker, x₁::AbstractVector, t₁=1.0, t₀=0.0; kwargs...)
     track!(tracker, x₁, t₁, t₀; kwargs...)
     CoreTrackerResult(tracker)
end

"""
     track!(tracker, x₁, t₁=1.0, t₀=0.0; setup_patch=true, checkstartvalue=true, loop::Bool=false)::CoreTrackerResult

Track a value `x₁` from `t₁` to `t₀` using the given `CoreTracker` `tracker`.
Returns one of the enum values of `CoreTrackerStatus.states` indicating the status.
If the tracking was successfull it is `CoreTrackerStatus.success`.
If `setup_patch` is `true` then [`setup!`](@ref) is called at the beginning
of the tracking.

    track!(x₀, tracker, x₁, t₁=1.0, t₀=0.0; options...)::CoreTrackerStatus.states

Additionally also stores the result in `x₀` if the tracking was successfull.
"""
function track!(x₀, tracker::CoreTracker, x₁, t₁=1.0, t₀=0.0; setup_patch::Bool=tracker.options.update_patch,
            loop::Bool=false, checkstartvalue::Bool=!loop)
     _track!(tracker, x₁, t₁, t₀, setup_patch, checkstartvalue, loop)
     retcode = status(tracker)
     if retcode == CoreTrackerStatus.success
         x₀ .= current_x(tracker)
     end
     retcode
end
@inline function track!(tracker::CoreTracker, x₁, t₁=1.0, t₀=0.0;
        setup_patch::Bool=tracker.options.update_patch,
        checkstartvalue::Bool=true,
        loop::Bool=false)
    _track!(tracker, x₁, t₁, t₀, setup_patch, checkstartvalue, loop)
end

function _track!(tracker::CoreTracker, x₁, t₁, t₀,
                    setup_patch::Bool, checkstartvalue::Bool, loop::Bool)
    if t₁ == t₀
        tracker.state.status = CoreTrackerStatus.success
        return tracker.state.status
    end

    setup!(tracker, x₁, t₁, t₀, setup_patch, checkstartvalue, loop)

    while tracker.state.status == CoreTrackerStatus.tracking
        step!(tracker)
        check_terminated!(tracker)
    end
    if tracker.state.status == CoreTrackerStatus.success
        refine!(tracker)
    end

    tracker.state.status
end

"""
    setup!(coretracker, x₁, t₁=1.0, t₀=0.0, setup_patch=coretracker.options.update_patch, checkstartvalue=true, compute_ẋ=true)

Setup `coretracker` to track `x₁` from `t₁` to `t₀`. Use this if you want to use the
coretracker as an iterator.
"""
function setup!(tracker::CoreTracker, x₁::AbstractVector, t₁=1.0, t₀=0.0, setup_patch=tracker.options.update_patch, checkstartvalue=true, loop::Bool=false)
    @unpack state, cache = tracker

    try
        reset!(state, x₁, t₁, t₀, tracker.options, setup_patch, loop)
        reset!(cache.predictor, state.x, t₁)
        checkstartvalue && checkstartvalue!(tracker)
        if !loop
            compute_ẋ!(state, cache, tracker.options)
            setup!(cache.predictor, cache.homotopy, state.x, state.ẋ, current_t(state), state.jacobian)
        end
    catch err
        if !(err isa LinearAlgebra.SingularException)
            rethrow(err)
        end
        tracker.state.status = CoreTrackerStatus.terminated_singularity
    end
    tracker
end

function checkstartvalue(tracker::CoreTracker, x, t)
    embed!(tracker.state.x̄, x)
    init_auto_scaling!(tracker.state.inner_product, tracker.state.x̄, tracker.options.auto_scaling_options)
    result = correct!(tracker.state.x̄, tracker, x, t; update_jacobian_infos=true)
    isconverged(result)
end

function checkstartvalue!(tracker::CoreTracker)
    result = correct!(tracker.state.x̄, tracker; update_jacobian_infos=true)
    if isconverged(result)
        tracker.state.x .= tracker.state.x̄
    else
        tracker.state.status = CoreTrackerStatus.terminated_invalid_startvalue
    end
    nothing
end

function compute_ẋ!(state, cache, options::CoreTrackerOptions)
    @inbounds jacobian_and_dt!(state.jacobian.J, cache.out, cache.homotopy, state.x, current_t(state))
    # apply row scaling to J and compute factorization
    updated_jacobian!(state.jacobian)

    @inbounds for i in eachindex(cache.out)
        cache.out[i] = -cache.out[i]
    end
    solve!(state.ẋ, state.jacobian, cache.out)
    nothing
end

@inline function correct!(x̄, tracker::CoreTracker,
        x=tracker.state.x,
        t=tracker.state.segment[tracker.state.s],
        norm=tracker.state.inner_product;
        accuracy::Float64=tracker.options.accuracy,
        max_iters::Int=tracker.options.max_corrector_iters,
        precision::PrecisionOption=tracker.options.precision,
        update_jacobian_infos::Bool=false,
        use_qr::Bool=false)

    correct!(x̄, tracker.corrector, tracker.cache.corrector,
             tracker.cache.homotopy, x, t, norm, tracker.state.jacobian,
             accuracy, max_iters, tracker.state.step_size;
             precision=precision,
             update_jacobian_infos=update_jacobian_infos, use_qr=use_qr)
end

function step!(tracker::CoreTracker)
    @unpack state, cache, options = tracker
    @unpack x, x̂, x̄, ẋ = state
    H = cache.homotopy

    try
        t, Δt = current_t(state), current_Δt(state)
        predict!(x̂, tracker.predictor, cache.predictor, H, x, t, Δt, ẋ, tracker.state.jacobian)
        # check if we need to update the jacobian_infos
        update_jacobian_infos =
                state.last_step_failed ||
                state.steps_jacobian_info_update ≥ options.steps_jacobian_info_update
        # reset counter
        update_jacobian_infos && (state.steps_jacobian_info_update = 0)
        result = correct!(x̄, tracker, x̂, t + Δt; update_jacobian_infos=update_jacobian_infos)
        # println("ω: ", result.ω, " t: ", real(t))
        if isconverged(result)
            # Step is accepted, assign values
            state.accepted_steps += 1
            x .= x̄
            state.s += state.Δs
            state.accuracy = result.accuracy
            # Step size change
            state.Δs_prev = state.Δs
            update_stepsize!(tracker, result)
            if state.patch !== nothing && options.update_patch
                changepatch!(state.patch, x)
            end
            options.auto_scaling && auto_scaling!(state, options)
            # update derivative
            compute_ẋ!(state, cache, options)
            # tell the predictors about the new derivative if they need to update something
            update!(cache.predictor, H, x, ẋ, t + Δt, state.jacobian)
            state.steps_jacobian_info_update += 1

            # Check termination criterion: we became too ill-conditioned
            if options.terminate_ill_conditioned && is_ill_conditioned(state, options) && state.last_step_failed
                state.status = CoreTrackerStatus.terminated_ill_conditioned
                return nothing
            end

            state.last_step_failed = false
        else
            # We have to reset the patch
            state.rejected_steps += 1
            update_rank!(state.jacobian)
            # Step failed, so we have to try with a new (smaller) step size
            update_stepsize!(tracker, result)
            Δt = current_Δt(state)
            state.last_step_failed = true
        end
    catch err
        if !(err isa LinearAlgebra.SingularException)
            rethrow(err)
        end
        tracker.state.status = CoreTrackerStatus.terminated_singularity
    end
    nothing
end

function is_ill_conditioned(state::CoreTrackerState, options::CoreTrackerOptions)
    d = max(log₁₀(cond(state)), digits_lost(state))
    d > options.max_lost_digits || state.jacobian.corank > 0
end

function is_overdetermined_tracking(tracker)
    m, n = size(tracker.cache.homotopy)
    m > n
end

function check_min_step_size!(tracker::CoreTracker{<:LogHomotopy})
    @unpack state, options = tracker
    if abs((exp(-state.Δs) - 1) * exp(-state.s)) < options.min_step_size
        state.status = CoreTrackerStatus.terminated_step_size_too_small
    end
end
function check_min_step_size!(tracker::CoreTracker)
    @unpack state, options = tracker
    if state.Δs < options.min_step_size
        state.status = CoreTrackerStatus.terminated_step_size_too_small
    end
end

g(Θ) = sqrt(1+4Θ) - 1
# Choose 0.25 instead of 1.0 due to Newton-Kantorovich theorem
δ(opts::CoreTrackerOptions, ω, μ) = @fastmath min(√(0.5ω) * τ(opts, μ), 0.25)
τ(opts::CoreTrackerOptions, μ) = @fastmath opts.accuracy^(1/(2opts.max_corrector_iters - μ))

function update_stepsize!(tracker::CoreTracker, result::CorrectorResult)
    @unpack state, options = tracker
    model = state.step_size
    ord = order(tracker.predictor)

    if options.simple_step_size_alg
        simple_step_size_alg!(tracker, result)
        return nothing
    end

    # The step size control is described in https://arxiv.org/abs/1902.02968

    # we have to handle the special case that there is only 1 iteration
    # in this case we cannot estimate ω and therefore just assume ω = 2
    # Also note ||x̂-x|| = ||Δx₀||
    if result.iters == 1
        ω = isnan(model.ω) ? 2.0 : 0.5model.ω
        d_x̂_x̄ = result.norm_Δx₀
    else
        ω = result.ω₀
        d_x̂_x̄ = distance(state.x̂, state.x, state.inner_product)
    end
    Δx₀ = result.norm_Δx₀
    if isconverged(result)
        δ_N_ω = δ(options, ω, 0.25)
        # This is to handle the edge case that g(δ_N_ω) >> (ω * d_x̂_x̄) > 0 but
        # at the same time δ_N_ω < eps(). Since then g(δ_N_ω) = 0
        δ_N_ω = max(δ_N_ω, 1e-15) #
        λ = g(δ_N_ω) / (ω * d_x̂_x̄)
        Δs = nthroot(λ, ord) * state.Δs
        if state.last_step_failed
            Δs = min(Δs, state.Δs)
        end
    else
        if isfinite(model.ω)
            ω = max(ω, model.ω)
        end
        δ_N_ω = δ(options, ω, 0.25)
        ω_η = 0.5ω * Δx₀
        if δ_N_ω < ω_η
            λ = g(δ_N_ω) / g(ω_η)
            Δs = min(nthroot(λ, ord), 0.9) * state.Δs
        else
            Δs = 0.5 * state.Δs
        end
        # If we fail consecutively reduce step size more aggressively.
        if state.last_step_failed
            Δs = min(Δs, 0.25 * state.Δs)
        end
    end
    model.expected_Δx₀ = 2δ_N_ω / ω
    model.ω = ω
    Δs = min(Δs, options.max_step_size)
    if (length(state.segment) - state.s) < Δs
        state.Δs = length(state.segment) - state.s
    else
        state.Δs = Δs
        check_min_step_size!(tracker)
    end

    nothing
end

function simple_step_size_alg!(tracker::CoreTracker, result::CorrectorResult)
    @unpack state, options = tracker
    if isconverged(result)
        state.consecutive_successfull_steps += 1
        if state.consecutive_successfull_steps == 5
            Δs′ = 2 * state.Δs
            state.consecutive_successfull_steps = 0
        else
            Δs′ = state.Δs
        end
    else
        state.consecutive_successfull_steps = 0
        Δs′ = 0.5 * state.Δs
    end

    state.Δs = min(Δs′, length(state.segment) - state.s)

    if !isconverged(result)
        check_min_step_size!(tracker)
    end
end

function auto_scaling!(state::CoreTrackerState, options::CoreTrackerOptions)
    auto_scaling!(state.inner_product, state.x, options.auto_scaling_options)
end
function auto_scaling!(ip::WeightedIP, x::AbstractVector, opts::AutoScalingOptions)
    norm_x = ip(x)
    for i in 1:length(x)
        wᵢ = (abs(x[i]) + ip.weight[i]) / 2
        if wᵢ < opts.scale_min * norm_x
            wᵢ = opts.scale_min * norm_x
        elseif wᵢ > opts.scale_max * norm_x
            wᵢ = opts.scale_max * norm_x
        end
        ip.weight[i] = max(wᵢ, opts.scale_abs_min)
    end
    nothing
end
auto_scaling!(ip::EuclideanIP, x::AbstractVector, opts::AutoScalingOptions) = nothing

function check_terminated!(tracker::CoreTracker)
    if abs(tracker.state.s - length(tracker.state.segment)) < 2eps(length(tracker.state.segment))
        tracker.state.status = CoreTrackerStatus.success
    elseif iters(tracker) ≥ tracker.options.max_steps
        tracker.state.status = CoreTrackerStatus.terminated_maximal_iterations
    end
    nothing
end

function refine!(tracker::CoreTracker, accuracy=tracker.options.refinement_accuracy)
    if tracker.state.accuracy < accuracy
        return
    end
    result = correct!(tracker.state.x̄, tracker;
        accuracy=accuracy,
        max_iters=tracker.options.max_refinement_iters)
    if isconverged(result)
        tracker.state.x .= tracker.state.x̄
        tracker.state.accuracy = result.accuracy
    end
    nothing
end

function Base.iterate(tracker::CoreTracker, state=nothing)
    state === nothing && return tracker, 1

    if tracker.state.status == CoreTrackerStatus.tracking
        step!(tracker)
        check_terminated!(tracker)

        if tracker.state.status == CoreTrackerStatus.success
            refine!(tracker)
        end

        tracker, state + 1
    else
        nothing
    end
end

"""
    affine_tracking(tracker)

Returns `true` if the path tracker tracks in affine space. If `false` then then the path tracker
tracks in some affine chart of the projective space.
"""
affine_tracking(tracker::CoreTracker) = !isa(tracker.state.x, ProjectiveVectors.PVector)

#################
## Query State ##
#################
"""
     current_t(tracker::CoreTracker)

Current `t`.
"""
current_t(tracker::CoreTracker) = current_t(tracker.state)
current_t(state::CoreTrackerState) = state.segment[state.s]

"""
     current_Δt(tracker::CoreTracker)

Current step_size `Δt`.
"""
current_Δt(tracker::CoreTracker) = current_Δt(tracker.state)
current_Δt(state::CoreTrackerState) = state.segment[state.Δs] - state.segment.start

"""
     iters(tracker::CoreTracker)

Current number of iterations.
"""
iters(tracker::CoreTracker) = iters(tracker.state)
iters(state::CoreTrackerState) = state.accepted_steps + state.rejected_steps

"""
     status(tracker::CoreTracker)

Current status.
"""
status(tracker::CoreTracker) = status(tracker.state)
status(state::CoreTrackerState) = state.status

"""
    current_x(tracker::CoreTracker)

Return the current value of `x`.
"""
current_x(tracker::CoreTracker) = current_x(tracker.state)
current_x(state::CoreTrackerState) = state.x

"""
    LinearAlgebra.cond(tracker::CoreTracker)

Returns the currently computed approximation of the condition number of the
Jacobian.
"""
cond(tracker::CoreTracker) = LinearAlgebra.cond(tracker.state)
cond(state::CoreTrackerState) = state.jacobian.cond

"""
    digits_lost(tracker::CoreTracker)

Returns the currently computed approximation of the number of digits lost
during the linear system solving in Newton's method.
"""
digits_lost(tracker::CoreTracker) = digits_lost(tracker.state)
digits_lost(state::CoreTrackerState) = unpack(state.jacobian.digits_lost, 0.0)

"""
    inner(tracker::CoreTracker)

Returns the inner product used to compute distance during the path tracking.
"""
inner(tracker::CoreTracker) = inner(tracker.state)
inner(state::CoreTrackerState) = state.inner_product

"""
    options(tracker::CoreTracker)

Returns the options used in the tracker.
"""
options(tracker::CoreTracker) = tracker.options

##################
# Modify options #
##################
"""
     accuracy(tracker::CoreTracker)

Current accuracy.
"""
accuracy(tracker::CoreTracker) = tracker.options.accuracy

"""
     set_accuracy!(tracker::CoreTracker, accuracy, update_max_lost_digits=true)

Set the current accuracy to `accuracy`. If `update_max_lost_digits` is `true` then
the setting `max_lost_digits` will be updated to the default setting.
"""
function set_accuracy!(tracker::CoreTracker, accuracy, update_max_lost_digits::Bool=true)
    @unpack options = tracker
    options.accuracy = accuracy
    if update_max_lost_digits
        options.max_lost_digits = default_max_lost_digits(options.precision, accuracy)
    end
    tracker
end

"""
     refinement_accuracy(tracker::CoreTracker)

Current refinement accuracy.
"""
refinement_accuracy(tracker::CoreTracker) = tracker.options.refinement_accuracy

"""
     set_max_refinement_iters!(tracker::CoreTracker, accuracy)

Set the current refinement accuracy to `accuracy`.
"""
set_refinement_accuracy!(T::CoreTracker, accuracy) = T.options.refinement_accuracy = accuracy

"""
     max_refinement_iters(tracker::CoreTracker)

Current refinement max_steps.
"""
max_refinement_iters(T::CoreTracker) = T.options.max_refinement_iters

"""
     set_max_refinement_iters!(tracker::CoreTracker, n)

Set the current refinement max_steps to `n`.
"""
set_max_refinement_iters!(T::CoreTracker, n) = T.options.max_refinement_iters = n

"""
     max_corrector_iters(tracker::CoreTracker)

Current correction max_steps.
"""
max_corrector_iters(T::CoreTracker) = T.options.max_corrector_iters

"""
     set_max_corrector_iters!(tracker::CoreTracker, n)

Set the correction max_steps to `n`.
"""
set_max_corrector_iters!(T::CoreTracker, n) = T.options.max_corrector_iters = n

"""
     max_step_size (tracker::CoreTracker)

Current maximal step size.
"""
max_step_size(T::CoreTracker) = T.options.max_step_size

"""
     set_max_corrector_iters!(tracker::CoreTracker, Δs)

Set the maximal step size to `Δs`.
"""
set_max_step_size!(T::CoreTracker, Δs) = T.options.max_step_size  = Δs

################
# PathIterator #
################
struct PathIterator{Tracker<:CoreTracker}
    tracker::Tracker
    t_real::Bool
end
Base.IteratorSize(::Type{<:PathIterator}) = Base.SizeUnknown()
Base.IteratorEltype(::Type{<:PathIterator}) = Base.HasEltype()

"""
    iterator(tracker::CoreTracker, x₁, t₁=1.0, t₀=0.0; affine=true)

Prepare a tracker to make it usable as a (stateful) iterator. Use this if you want to inspect a specific
path. In each iteration the tuple `(x,t)` is returned.
If `affine == true` then `x` is the affine solution (internally we compute in projective space).

## Example

Assume you have `CoreTracker` `tracker` and you wan to track `x₁` from 1.0 to 0.25:
```julia
for (x,t) in iterator(tracker, x₁, 1.0, 0.25)
    println("x at t=\$t:")
    println(x)
end
```

Note that this is a stateful iterator. You can still introspect the state of the tracker.
For example to check whether the tracker was successfull
(and did not terminate early due to some problem) you can do
```julia
println("Success: ", status(tracker) == CoreTrackerStatus.success)
```
"""
function iterator(tracker::CoreTracker, x₁, t₁=1.0, t₀=0.0; kwargs...)
    setup!(tracker, x₁, t₁, t₀; kwargs...)
    PathIterator(tracker, typeof(t₁ - t₀) <: Real)
end

function current_x_t(iter::PathIterator)
    x = current_x(iter.tracker)
    t = current_t(iter.tracker)
    (copy(x), iter.t_real ? real(t) : t)
end

function Base.iterate(iter::PathIterator, state=nothing)
    state === nothing && return current_x_t(iter), 1
    iter.tracker.state.status != CoreTrackerStatus.tracking && return nothing

    step_done = false
    while !step_done && (iter.tracker.state.status == CoreTrackerStatus.tracking)

        step!(iter.tracker)
        check_terminated!(iter.tracker)

        if iter.tracker.state.status == CoreTrackerStatus.success
            refine!(iter.tracker)
        end

        step_done = !iter.tracker.state.last_step_failed
    end
    current_x_t(iter), state + 1
end

#############################
# Convencience constructors #
#############################
"""
    coretracker_startsolutions(args...; kwargs...)

Construct a [`CoreTracker`](@ref) and `startsolutions` in the same way `solve`
does it. This also takes the same input arguments as `solve`. This is convenient if you want
to investigate single paths.
"""
function coretracker_startsolutions(args...; kwargs...)
    supported, rest = splitkwargs(kwargs,problem_startsolutions_supported_keywords)
    prob, startsolutions = problem_startsolutions(args...; supported...)
    tracker = CoreTracker(prob, start_solution_sample(startsolutions); rest...)
    (tracker=tracker, startsolutions=startsolutions)
end

"""
    coretracker(args...; kwargs...)

Construct a [`CoreTracker`](@ref) in the same way `solve`
does it. This also takes the same input arguments as `solve` with the exception that you do not need to specify startsolutions.
This is convenient if you want to investigate single paths.

## Examples

### Obtain single solution
We want to construct a path tracker to track a parameterized system `f` with parameters `p`
from the parameters `a` to `b`.
```julia
tracker = coretracker(f, parameters=p, p₁=a, p₀=b)
```
You then can obtain a single solution at `b` by using
```julia
x_b = track(tracker, x_a).x
```

### Trace a path
To trace a path you can use the [`iterator`](@ref) method.

```julia
tracker = coretracker(f, parameters=p, p₁=a, p₀=b, max_step_size =0.01)
for (x, t) in iterator(tracker, x₁)
    @show (x,t)
end
```

If we want to guarantee smooth traces we can limit the maximal step size.
```julia
tracker = coretracker(f, parameters=p, p₁=a, p₀=b, max_step_size =0.01)
for (x, t) in iterator(tracker, x₁)
    @show (x,t)
end
```
"""
function coretracker(args...; kwargs...)
    tracker, _ = coretracker_startsolutions(args...; kwargs...)
    tracker
end


text:
using NCDatasets
using Formatting

using Statistics

Re  = 6371000.0 #  m
c_p  = 1004.0    #  J / kg K
g   = 9.81      #  m / s^2
Lw  = 2.26e6    #  J / kg

t_beg = 240-10*12+1
t_end = 240
t_rng = t_beg:t_end
years = (t_end - t_beg + 1) / 12.0
nc_file_dir = "extract_nc_files"
domain_nc_file_dir = "domain_nc_files"

casenames = [
    "lowres_STD_SOM",
    "lowres_SSM_SOM",
    "lowres_SSM_NK",
    "lowres_SSM_SOM_noQflux",
]

linestyle = ["k-", "r-", "g-", "b-"]

mreplace = (x,) -> replace(x, missing=>NaN)

Dataset("$domain_nc_file_dir/domain.lnd.fv4x5_gx3v7.091218.nc", "r") do ds
    global mask
    mask = ds["mask"][:] |> mreplace

    mask[mask.!=0] .= 1
end

Dataset("$nc_file_dir/$(casenames[1]).h1.nc", "r") do ds
    global lat, lon, lev, ilev, rec

    lat = ds["lat"][:] |> mreplace
    lon = ds["lon"][:] |> mreplace

    ilev = (ds["ilev"][:] |> mreplace) * 100.0

    rec = ds["time"][:] |> mreplace

    lev = (ilev[2:end] + ilev[1:end-1]) / 2.0
end

# mass
Δp = (ilev[2:end] - ilev[1:end-1])

function calMeanTransport(v)

    # Get the zonal mean
    v = mean(v, dims=(1,))[1, :, :, :]

    # Weight with mass
    for j=1:length(lat), k=1:length(lev), l=1:length(rec)
        v[j, k, l] *= Δp[k] * cos(lat[j] |>  deg2rad)
    end

    v .*= 2 * π * Re / g

    # Get the vertical integration and time mean
    v = sum(v, dims=(2,))[:, 1, :]
    v = mean(v, dims=(2,))[:, 1]

    return v
end

function integrate(x, dydx)
    y = copy(dydx) * 0.0

    for i = 2:length(x)
        y[i] = y[i-1] + (dydx[i-1] + dydx[i]) * (x[i] - x[i-1]) / 2.0
    end

    return y
end

print("Import Pyplot... ")
using PyPlot
plt[:rcParams]["font.size"] = 20

#using PyCall
#@pyimport mpl_toolkits.basemap as basemap

println("Done. ")

#include("subcode_compare_models/global_surface_temperature.jl")
#include("subcode_compare_models/map_sea_surface_pressure.jl")
include("subcode_compare_models/map_SST.jl")

#=
### Streamfunction

include(joinpath(@__DIR__, "invert_streamfunction.jl"))

fig_compare, ax_compare = plt[:subplots](1, 1, figsize=(12,8))
ax_compare[:plot]([-90, 90], [0, 0], linewidth=2, color="#cccccc")

for i = 1:length(casenames)

    casename = casenames[i]

    fig, ax = plt[:subplots](1, 1, figsize=(12,8))

    Dataset("$nc_file_dir/$casename.h0.nc", "r") do ds

        global V = mean(replace(ds["V"][:], missing=>NaN), dims=(1, 4))[1, :, :, 1]
        ψ  = invertStreamfunction(lat * π / 180.0, ilev, V; g0=g, Re=Re) / 1e10

        println(size(ψ), ";", size(lat), ";", size(ilev))
        cs = ax[:contour](lat, ilev/100, ψ'[:,:], range(-20, step=2, stop=20), colors="k")

        ax[:clabel](cs, fmt="%d")
        ax[:set_xlabel]("Lat [deg]")
        ax[:set_ylabel]("Pressure [hPa]")

        ax[:set_title](casename)
        ax[:set_xlim](-90, 90)
        ax[:set_ylim](100, 1000)
        ax[:invert_yaxis]()

        fig[:savefig]("$(casename)_psi.png", dpi=200)

        ax_compare[:plot](lat, mean(ψ[:, 19:20], dims=(2,))[:,1], linestyle[i], label="$casename")

    end

end

ax_compare[:legend]()
ax_compare[:set_xlabel]("Lat [deg]")
ax_compare[:set_ylabel]("Streamfunction [ \$ \\times 10^{10} \$ kg / s]")

ax_compare[:set_title]("Streamfunction along 510 hPa (meridional mass transport below 510 hPa)")
ax_compare[:set_xlim](-90, 90)

fig_compare[:savefig]("compare_models_psi.png", dpi=200)

=#

#=

### Meridional Static Energy Transport

fig, ax = plt[:subplots](1, 1, figsize=(12,8))

ax[:plot]([-90, 90], [0, 0], linewidth=2, color="#cccccc")

for i = 1:length(casenames)

    casename = casenames[i]

    Dataset(casename * ".nc", "r") do ds
        global MHT, MGT, MLT
        MHT = c_p * ds["VT"][:]  |> calMeanTransport
        MGT = ds["VZ"][:]       |> calMeanTransport
        MLT = Lw * ds["VQ"][:]  |> calMeanTransport

    end

    MMSET = (MHT + MGT + MLT) / 1e15

    ax[:plot](lat, MMSET, linestyle[i], linewidth=2,  label=casename)

end

ax[:legend]()
ax[:set_xlabel]("Latitude [deg]")
ax[:set_ylabel]("[PW]")

ax[:set_title]("10 years")
ax[:set_xlim](-90, 90)
ax[:set_ylim](-15, 20)

fig[:savefig]("compare_models_MSEFLX.png", dpi=200)

### Precipitation
fig, ax = plt[:subplots](1, 1, figsize=(12,8))

ax[:plot]([-90, 90], [0, 0], linewidth=2, color="#cccccc")

for i = 1:length(casenames)

    casename = casenames[i]

    Dataset(casename * ".nc", "r") do ds
        global tot_prec = replace(ds["PRECC"][:] + ds["PRECL"][:], missing=>NaN)

        tot_prec = mean(tot_prec, dims=(1,3))[1, :, 1] # average over time and longitude
    end

    ax[:plot](lat, tot_prec * 365.0 * 86400.0 * 1000.0, linestyle[i], linewidth=2,  label=casename)
end

ax[:legend]()
ax[:set_xlabel]("Latitude [deg]")
ax[:set_ylabel]("Precipitation [mm / yr]")

ax[:set_title]("10 years")
ax[:set_xlim](-90, 90)
#ax[:set_ylim](, 20)

fig[:savefig]("compare_models_PRECIP.png", dpi=200)

### Energy convergence

fig_met, ax_met = plt[:subplots](1, 1, figsize=(12,8))
ax_met[:plot]([-90, 90], [0, 0], linewidth=2, color="#cccccc")

for i = 1:length(casenames)

    fig, ax = plt[:subplots](1, 1, figsize=(12,8))
    ax[:plot]([-90, 90], [0, 0], linewidth=2, color="#cccccc")

    casename = casenames[i]

    Dataset(casename * ".h0.nc", "r") do ds
        t_rng = (:, :,  240-5*12+1:240)
        FFLX_TOA = mean( - mean(ds["FSNT"][t_rng...], dims=(3,)) + mean(ds["FLNT"][t_rng...], dims=(3,)), dims=(1,))[1, :, 1]
        FFLX_SFC = mean( - mean(ds["FSNS"][t_rng...], dims=(3,)) + mean(ds["FLNS"][t_rng...], dims=(3,)), dims=(1,))[1, :, 1]
        HFLX_SFC = mean( mean(ds["SHFLX"][t_rng...], dims=(3,)) + mean(ds["LHFLX"][t_rng...], dims=(3,)), dims=(1,))[1, :, 1]

        FFLX_TOA .*= Re * cos.(deg2rad.(lat)) * 2π
        FFLX_SFC .*= Re * cos.(deg2rad.(lat)) * 2π
        HFLX_SFC .*= Re * cos.(deg2rad.(lat)) * 2π

        EFLX_CONV = - ( FFLX_TOA - FFLX_SFC - HFLX_SFC )

        ax[:plot](lat, FFLX_TOA, "r--", linewidth=2,  label="FTOA")
        ax[:plot](lat, FFLX_SFC, "b--", linewidth=2,  label="FSFC")
        ax[:plot](lat, HFLX_SFC, "g--", linewidth=2,  label="HFLX")
        ax[:plot](lat, EFLX_CONV, "k-", linewidth=2,  label="CONV")

        MET = integrate(Re * deg2rad.(lat), EFLX_CONV)
        ax_met[:plot](lat, MET / 1e15, linestyle[i], linewidth=2, label=casename)

    end

    ax[:legend]()
    ax[:set_xlabel]("Latitude [deg]")
    ax[:set_ylabel]("[PW / m]")

    ax[:set_title](format("10 years in case {}", casename))
    ax[:set_xlim](-90, 90)

    fig[:savefig](format("compare_models_EBM_{}.png", casename), dpi=200)

end

ax_met[:legend]()
ax_met[:set_xlabel]("Latitude [deg]")
ax_met[:set_ylabel]("[PW]")

ax_met[:set_title](format("16-20 years annual implied northward heat transport"))
ax_met[:set_xlim](-90, 90)
ax_met[:set_ylim](-6, 6)

fig_met[:savefig](format("compare_models_MET.png"), dpi=200)

#ax[:set_ylim](-15, 20)

### Temperature trend

fig, ax = plt[:subplots](1, 1, figsize=(12,8))

dθ = deg2rad(lat[2] - lat[1])
for i = 1:length(casenames)

    casename = casenames[i]

    Dataset(casename * ".h1.nc", "r") do ds

        TREFHT = sum(mean(ds["TREFHT"][:], dims=(1,)), dims=(1,) )[1, :, :]
        Tavg = zeros(Float64, length(rec))
        for l = 1:length(rec)
            TREFHT[:, l] .*= cos.(deg2rad.(lat))
            Tavg[l] = sum(TREFHT[:, l]) * dθ / 2.0
        end

        Tavg .-= 273.15

        ax[:plot](collect(Float64, 0:length(rec)-1) / 365.0, Tavg, linestyle[i], linewidth=2,  label=casename)
    end

end

ax[:legend]()
ax[:set_xlabel]("Time [month]")
ax[:set_ylabel]("Temperature")

ax[:set_title]("Global Mean surface temperature of 20 years")
ax[:set_xlim](0, 20)
ax[:set_ylim](0, 20)

fig[:savefig]("compare_models_ECNVG.png", dpi=200)

=#

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
type Mul2 <: Op; ybuf; x1; x2; dx1; dx2; Mul2()=new(); end

params(::Mul2)=Any[]
ninputs(::Mul2)=2
ysize(::Mul2,x1,x2)=size(x1)
overwrites(::Mul2)=false
back_reads_x(::Mul2)=true
back_reads_y(::Mul2)=false

# x1,x2 is a pair of similarly sized input matrices

function forw(l::Mul2, x1, x2; y=nothing, o...)
    (l.x1,l.x2) = (x1,x2)
    y = initforw(l,x1,x2,y)
    mul2!(y,x1,x2)
end

function initforw(l::Mul2, x1, x2, y)
    issimilar(x1,x2) || error("Input mismatch")
    y == nothing && (y = similar!(l, :ybuf, x1))
    issimilar0(y,x1) || error("Input output mismatch")
    return y
end

function back(l::Mul2, dy; dx=nothing, x=(l.x1,l.x2), returndx=true, o...)
    returndx || return
    (dx1,dx2) = initback(l, dy, x, dx)
    mul2!(dx1,dy,x[2])
    mul2!(dx2,dy,x[1])
    return (dx1,dx2)
end

function initback(l::Mul2, dy, x, dx)
    length(x)==2 || error("Need two inputs")
    (x1,x2) = x
    issimilar0(x1,x2) || error("Input mismatch")
    issimilar0(x1,dy) || error("Input gradient mismatch")
    issimilar0(x2,dy) || error("Input gradient mismatch")
    dx == nothing && (dx = (similar!(l,:dx1,x1), similar!(l,:dx2,x2)))
    issimilar0(x1,dx[1]) || error("Gradient mismatch")
    issimilar0(x2,dx[2]) || error("Gradient mismatch")
    dx[1] === dx[2] && error("Need two different dx")
    return dx
end

# if any of the matrices is nothing, that represents zero, we return nothing

forw(l::Mul2, ::Void, ::Void; o...)=nothing
forw(l::Mul2, x1, ::Void; o...)=nothing
forw(l::Mul2, ::Void, x2; o...)=nothing

ysize(l::Mul2, ::Void, ::Void)=nothing
ysize(l::Mul2, x1, ::Void)=nothing
ysize(l::Mul2, ::Void, x2)=nothing

mul2!(c,a::Void,b::Void)=fill!(c,0)
mul2!(c,a,b::Void)=fill!(c,0)
mul2!(c,a::Void,b)=fill!(c,0)

issimilar0(a,b)=issimilar(a,b)
issimilar0(a::Void,b::Void)=true
issimilar0(a,b::Void)=true
issimilar0(a::Void,b)=true

similar!(l,n,::Void)=nothing    # don't get rid of existing storage but return nothing