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 this file:
# - fmi3xxx-structs and -functions from the FMI3-specification
# - helper-structs that are not part of the FMI3-specification, but help to parse the model description nicely (e.g. `fmi3ModelDescriptionDefaultExperiment`, ...)

"""
Source: FMISpec3.0-dev, Version D5ef1c1:2.2.2. Platform Dependent Definitions

To simplify porting, no C types are used in the function interfaces, but the alias types are defined in this section.
All definitions in this section are provided in the header file fmi3PlatformTypes.h. It is required to use this definition for all binary FMUs.
"""
const fmi3Float32 = Cfloat
const fmi3Float64 = Cdouble
const fmi3Int8 = Cchar
const fmi3UInt8 = Cuchar
const fmi3Int16 = Cshort
const fmi3UInt16 = Cushort
const fmi3Int32 = Cint
const fmi3UInt32 = Cuint
const fmi3Int64 = Clonglong
const fmi3UInt64 = Culonglong
const fmi3Boolean = Cuchar
const fmi3Char = Cchar
const fmi3String = Ptr{fmi2Char}
const fmi3Byte = Cuchar
const fmi3Binary = Ptr{fmi3Byte}
const fmi3ValueReference = Cuint
const fmi3FMUState = Ptr{Cvoid}
const fmi3Instance = Ptr{Cvoid}
const fmi3InstanceEnvironment = Ptr{Cvoid}
const fmi3Enum = Array{Array{String}} # TODO: correct it
const fmi3Clock = Cint

const fmi3False = fmi3Boolean(false)
const fmi3True = fmi3Boolean(true)

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.3. Status Returned by Functions
Defines the status flag (an enumeration of type fmi3Status defined in file fmi3FunctionTypes.h) that is returned by functions to indicate the success of the function call:
The status has the following meaning:
fmi3OK - The call was successful. The output argument values are defined.

fmi3Warning - A non-critical problem was detected, but the computation can continue. The output argument values are defined. Function logMessage should be called by the FMU with further information before returning this status, respecting the current logging settings.
[In certain applications, e.g. in a prototyping environment, warnings may be acceptable. For production environments warnings should be treated like errors unless they can be safely ignored.]

fmi3Discard - The call was not successful and the FMU is in the same state as before the call. The output argument values are not defined, but the computation can continue. Function logMessage should be called by the FMU with further information before returning this status, respecting the current logging settings. Advanced simulation algorithms can try alternative approaches to drive the simulation by calling the function with different arguments or calling another function. Otherwise the simulation algorithm has to treat this return code like fmi3Error and has to terminate the simulation.
[Examples for usage of fmi3Discard are handling of min/max violation, or signal numerical problems during model evaluation forcing smaller step sizes.]

fmi3Error - The call failed. The output argument values are undefined and the simulation cannot be continued. Function logMessage should be called by the FMU with further information before returning this status, respecting the current logging settings. If a function returns fmi3Error, it is possible to restore a previously retrieved FMU state by calling fmi3SetFMUState. Otherwise fmi3FreeInstance or fmi3Reset must be called. When detecting illegal arguments or a function call not allowed in the current state according to the respective state machine, the FMU must return fmi3Error. Other instances of this FMU are not affected by the error.

fmi3Fatal - The state of all instances of the model is irreparably corrupted. [For example, due to a runtime exception such as access violation or integer division by zero during the execution of an FMI function.] Function logMessage should be called by the FMU with further information before returning this status, respecting the current logging settings, if still possible. It is not allowed to call any other function for any instance of the FMU.
"""
const fmi3Status          = Cuint
const fmi3StatusOK        = Cuint(0)
const fmi3StatusWarning   = Cuint(1)
const fmi3StatusDiscard   = Cuint(2)
const fmi3StatusError     = Cuint(3)
const fmi3StatusFatal     = Cuint(4)

"""
Source: FMISpec3.0, Version D5ef1c1: 2.4.7.4. Variable Attributes
Enumeration that defines the causality of the variable. Allowed values of this enumeration:

parameter - A data value that is constant during the simulation (except for tunable parameters, see there) and is provided by the environment and cannot be used in connections, except for parameter propagation in terminals as described in Section 2.4.9.2.6. variability must be fixed or tunable. These parameters can be changed independently, unlike calculated parameters. initial must be exact or not present (meaning exact).

calculatedParameter - A data value that is constant during the simulation and is computed during initialization or when tunable parameters change. variability must be fixed or tunable. initial must be approx, calculated or not present (meaning calculated).

input - The variable value can be provided by the importer. [For example, the importer could forward the output of another FMU into this input.]

output -  The variable value can be used by the importer. [For example, this value can be forwarded to an input of another FMU.]
The algebraic relationship to the inputs can be defined via the dependencies attribute of <fmiModelDescription><ModelStructure><Output>.

local -  Local variables are:
- continuous states and their ContinuousStateDerivatives, ClockedStates, EventIndicators or InitialUnknowns. These variables are listed in the <fmiModelDescription><ModelStructure>.
- internal, intermediate variables or local clocks which can be read for debugging purposes and are not listed in the <fmiModelDescription><ModelStructure>.
Setting of local variables:
- In Initialization Mode and before, local variables need to be set if they do have start values or are listed as <InitialUnknown>.
- In super state Initialized, fmi3Set{VariableType} must not be called on any of the local variables. Only in Model Exchange, continuous states can be set with fmi3SetContinuousStates. Local variable values must not be used as input to another model or FMU.

independent - The independent variable (usually time [but could also be, for example, angle]). All variables are a function of this independent variable. variability must be continuous. Exactly one variable of an FMU must be defined as independent.
For Model Exchange the value is the last value set by fmi3SetTime. For Co-Simulation the value of the independent variable is lastSuccessfulTime return by the last call to fmi3DoStep or the value of argument intermediateUpdateTime of fmi3CallbackIntermediateUpdate. For Scheduled Execution the value of the independent variable is not defined. [The main purpose of this variable in Scheduled Execution is to define a quantity and unit for the independent variable.] The initial value of the independent variable is the value of the argument startTime of fmi3EnterInitializationMode for both Co-Simulation and Model Exchange.
If the unit for the independent variable is not defined, it is implicitly s (seconds). If one variable is defined as independent, it must be defined with a floating point type without a start attribute. It is not allowed to call function fmi3Set{VariableType} on an independent variable. Instead, its value is initialized with fmi3EnterInitializationMode and after initialization set by fmi3SetTime for Model Exchange and by arguments currentCommunicationPoint and communicationStepSize of fmi3DoStep for Co-Simulation FMUs. [The actual value can be inquired with fmi3Get{VariableType}.]

structuralParameter - The variable value can only be changed in Configuration Mode or Reconfiguration Mode. The variability attribute must be fixed or tunable. The initial attribute must be exact or not present (meaning exact). The start attribute is mandatory. A structural parameter must not have a <Dimension> element. A structural parameter may be referenced in <Dimension> elements. If a structural parameters is referenced in <Dimension> elements, it must be of type <UInt64> and its start attribute must be larger than 0. The min attribute might still be 0.

The default of causality is local.
A continuous-time state or an event indicator must have causality = local or output, see also Section 2.4.8.

[causality = calculatedParameter and causality = local with variability = fixed or tunable are similar. The difference is that a calculatedParameter can be used in another model or FMU, whereas a local variable cannot. For example, when importing an FMU in a Modelica environment, a calculatedParameter should be imported in a public section as final parameter, whereas a local variable should be imported in a protected section of the model.]

The causality of variables of type Clock must be either input or output.

Added prefix "fmi3" to help with redefinition of constans in enums.
"""
const fmi3Causality                     = Cuint
const fmi3CausalityParameter            = Cuint(0)
const fmi3CausalityCalculatedParameter  = Cuint(1)
const fmi3CausalityInput                = Cuint(2)
const fmi3CausalityOutput               = Cuint(3)
const fmi3CausalityLocal                = Cuint(4)
const fmi3CausalityIndependent          = Cuint(5)
const fmi3CausalityStructuralParameter  = Cuint(6)

"""
Source: FMISpec3.0, Version D5ef1c1: 2.4.7.4. Variable Attributes
Enumeration that defines the time dependency of the variable, in other words, it defines the time instants when a variable can change its value. [The purpose of this attribute is to define when a result value needs to be inquired and to be stored. For example, discrete variables change their values only at event instants (ME) or at a communication point (CS and SE) and it is therefore only necessary to inquire them with fmi3Get{VariableType} and store them at event times.] Allowed values of this enumeration:
constant - The value of the variable never changes.

fixed - The value of the variable is fixed after initialization, in other words, after fmi3ExitInitializationMode was called the variable value does not change anymore.

tunable - The value of the variable is constant between events (ME) and between communication points (CS and SE) due to changing variables with causality = parameter and variability = tunable. Whenever a parameter with variability = tunable changes, an event is triggered externally (ME or CS if events are supported), or the change is performed at the next communication point (CS and SE) and the variables with variability = tunable and causality = calculatedParameter or output must be newly computed. [tunable inputs are not allowed, see Table 18.]

discrete -
Model Exchange: The value of the variable is constant between external and internal events (= time, state, step events defined implicitly in the FMU).
Co-Simulation: By convention, the variable is from a real sampled data system and its value is only changed at communication points (including event handling). During intermediateUpdate, discrete variables are not allowed to change. [If the simulation algorithm notices a change in a discrete variable during intermediateUpdate, the simulation algorithm will delay the change, raise an event with earlyReturnRequested == fmi3True and during the communication point it can change the discrete variable, followed by event handling.]

continuous - Only a variable of type == fmi3GetFloat32 or type == fmi3GetFloat64 can be continuous.
Model Exchange: No restrictions on value changes (see Section 4.1.1).

The default is continuous for variables of type <Float32> and <Float64>, and discrete for all other types.

For variables of type Clock and clocked variables the variability is always discrete or tunable.

[Note that the information about continuous states is defined with elements <ContinuousStateDerivative> in <ModelStructure>.]

Added prefix "fmi3" to help with redefinition of constans in enums.
"""
const fmi3Variability             = Cuint
const fmi3VariabilityConstant     = Cuint(0)
const fmi3VariabilityFixed        = Cuint(1)
const fmi3VariabilityTunable      = Cuint(2)
const fmi3VariabilityDiscrete     = Cuint(3)
const fmi3VariabilityContinuous   = Cuint(4)

"""
Source: FMISpec3.0, Version D5ef1c1:2.4.7.5. Type specific properties
Enumeration that defines how the variable is initialized, i.e. if a fmi3Set{VariableType} is allowed and how the FMU internally treats this value in Instantiated and Initialization Mode.
For the variable with causality = independent, the attribute initial must not be provided, because its start value is set with the startTime parameter of fmi3EnterInitializationMode.

The attribute initial for other variables can have the following values and meanings:

exact - The variable is initialized with the start value (provided under the variable type element).

approx - The start value provides an approximation that may be modified during initialization, e.g., if the FMU is part of an algebraic loop where the variable might be an iteration variable and start value is taken as initial value for an iterative solution process.

calculated - The variable is calculated from other variables during initialization. It is not allowed to provide a start value.

If initial is not present, it is defined by Table 22 based on causality and variability. If initial = exact or approx, or causality = input, a start value must be provided. If initial = calculated, or causality = independent, it is not allowed to provide a start value.

[The environment decides when to use the start value of a variable with causality = input. Examples: * Automatic tests of FMUs are performed, and the FMU is tested by providing the start value as constant input. * For a Model Exchange FMU, the FMU might be part of an algebraic loop. If the input variable is iteration variable of this algebraic loop, then initialization starts with its start value.]

If fmi3Set{VariableType} is not called on a variable with causality = input, then the FMU must use the start value as value of this input.

Added prefix "fmi3" to help with redefinition of constans in enums.
"""
const fmi3Initial           = Cuint
const fmi3InitialExact      = Cuint(0)
const fmi3InitialApprox     = Cuint(1)
const fmi3InitialCalculated = Cuint(2)

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

Argument fmuType defines the type of the FMU:
- fmi3ModelExchange: FMU with initialization and events; between events simulation of continuous systems is performed with external integrators from the environment.
- fmi3CoSimulation: Black box interface for co-simulation.
- fmi3ScheduledExecution: Concurrent computation of model partitions on a single computational resource (e.g. CPU-core)
"""
const fmi3Type                      = Cuint
const fmi3TypeModelExchange         = Cuint(0)
const fmi3TypeCoSimulation          = Cuint(1)
const fmi3TypeScheduledExecution    = Cuint(2)

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.9.4. Scheduled Execution
Enumeration that defines the IntervalQualifiers which describe how to treat the intervals and intervalCounters arguments. They have the following meaning:
fmi3IntervalNotYetKnown -  is returned for a countdown aperiodic Clock for which the next interval is not yet known. This qualifier value can only be returned directly after the Clock was active and previous calls to fmi3GetInterval never returned fmi3IntervalChanged (nor fmi3IntervalUnchanged). In Scheduled Execution this return value means that the corresponding model partition cannot be scheduled yet.

fmi3IntervalUnchanged - is returned if a previous call to fmi3GetInterval already returned a value qualified with fmi3IntervalChanged which has not changed since. In Scheduled Execution this means the corresponding model partition has already been scheduled.

fmi3IntervalChanged - is returned to indicate that the value for the interval has changed for this Clock. Any previously returned intervals (if any) are overwritten with the current value. The new Clock interval is relative to the time of the current Event Mode or Clock Update Mode in contrast to the interval of a periodic Clock, where the interval is defined as the time between consecutive Clock ticks. In Scheduled Execution this means that the corresponding model partition has to be scheduled or re-scheduled (if a previous call to fmi3GetInterval returned fmi3IntervalChanged).
"""
const fmi3IntervalQualifier                     = Cuint
const fmi3IntervalQualifierIntervalNotYetKnown  = Cuint(0)
const fmi3IntervalQualifierIntervalUnchanged    = Cuint(1)
const fmi3IntervalQualifierIntervalChanged      = Cuint(2)

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

Enumeration that defines the dependencies a single unknown variable vunknown can have in relation to a known variable vknown. They have the following meaning:
dependent - no particular structure, f(.., v_{known,i}, ..)

Only for floating point type unknowns v_{unknown}:

constant - constant factor, c ⋅ v_{known,i} where c is an expression that is evaluated before fmi3EnterInitializationMode is called.

Only for floating point type unknowns v_{unknown} in Event and Continuous-Time Mode (ME) and at communication points (CS and SE), and not for <InitialUnknown> for Initialization Mode:

fixed - fixed factor, p⋅v_{known,i} where p is an expression that is evaluated before fmi3ExitInitializationMode is called.

tunable - tunable factor, p⋅v_{known,i} where p is an expression that is evaluated before fmi3ExitInitializationMode is called and in Event Mode due to event handling (ME) or at a communication point (CS and SE)

discrete - discrete factor, d⋅v_{known,i} where d is an expression that is evaluated before fmi3ExitInitializationMode is called and in Event Mode due to an external or internal event or at a communication point (CS and SE).
"""
const fmi3DependencyKind            = Cuint
const fmi3DependencyKindIndependent = Cuint(0)
const fmi3DependencyKindConstant    = Cuint(1)
const fmi3DependencyKindFixed       = Cuint(2)
const fmi3DependencyKindTunable     = Cuint(3)
const fmi3DependencyKindDiscrete    = Cuint(4)
const fmi3DependencyKindDependent   = Cuint(5)

"""
Source: FMISpec3.0, Version D5ef1c1: 2.4.7. Definition of Model Variables
                                     2.4.4. Definition of Types

"""
mutable struct fmi3DatatypeVariable
    # mandatory
    # TODO clock
    datatype::Union{Nothing, Type{fmi3String}, Type{fmi3Float64}, Type{fmi3Float32}, Type{fmi3Int8}, Type{fmi3UInt8}, Type{fmi3Int16}, Type{fmi3UInt16}, Type{fmi3Int32}, Type{fmi3UInt32}, Type{fmi3Int64}, Type{fmi3UInt64}, Type{fmi3Boolean}, Type{fmi3Binary}, Type{fmi3Char}, Type{fmi3Byte}, Type{fmi3Enum}}

    # Optional
    canHandleMultipleSet::Union{fmi3Boolean, Nothing}
    intermediateUpdate::fmi3Boolean
    previous::Union{fmi3UInt32, Nothing}
    clocks
    declaredType::Union{String, Nothing}
    start::Union{String, fmi3String, fmi3Float32, fmi3Float64, fmi3Int8, fmi3UInt8, fmi3Int16, fmi3UInt16, fmi3Int32, fmi3UInt32, fmi3Int64, fmi3UInt64, fmi3Boolean, fmi3Binary, fmi3Char, fmi3Byte, fmi3Enum, Array{fmi3Float32}, Array{fmi3Float64}, Array{fmi3Int32}, Array{fmi3UInt32}, Array{fmi3Int64}, Array{fmi3UInt64},  Nothing}
    min::Union{fmi3Float64,fmi3Int32, fmi3UInt32, fmi3Int64, Nothing}
    max::Union{fmi3Float64,fmi3Int32, fmi3UInt32, fmi3Int64, Nothing}
    initial::Union{fmi3Initial, Nothing}
    quantity::Union{String, Nothing}
    unit::Union{String, Nothing}
    displayUnit::Union{String, Nothing}
    relativeQuantity::Union{fmi3Boolean, Nothing}
    nominal::Union{fmi3Float64, Nothing}
    unbounded::Union{fmi3Boolean, Nothing}
    derivative::Union{fmi3UInt32, Nothing}
    reinit::Union{fmi3Boolean, Nothing}
    mimeType::Union{String, Nothing}
    maxSize::Union{fmi3UInt32, Nothing}

    # # used by Clocks TODO
    # canBeDeactivated::Union{fmi3Boolean, Nothing}
    # priority::Union{fmi3UInt32, Nothing}
    # intervall
    # intervallDecimal::Union{fmi3Float32, Nothing}
    # shiftDecimal::Union{fmi3Float32, Nothing}
    # supportsFraction::Union{fmi3Boolean, Nothing}
    # resolution::Union{fmi3UInt64, Nothing}
    # intervallCounter::Union{fmi3UInt64, Nothing}
    # shitftCounter::Union{fmi3Int32, Nothing}

    # additional (not in spec)
    unknownIndex::Integer
    dependencies::Array{Integer}
    dependenciesValueReferences::Array{fmi3ValueReference}

    # Constructor
    fmi3DatatypeVariable() = new()
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.4.7. Definition of Model Variables

A fmi3ModelVariable describes the the type, name, valueRefence and optional information for every variable in the Modeldescription.
"""
mutable struct fmi3ModelVariable
    #mandatory
    name::String
    valueReference::fmi3ValueReference
    datatype::fmi3DatatypeVariable

    # Optional
    description::String

    causality::fmi3Causality
    variability::fmi3Variability
    # initial::fmi3Initial ist in fmi3 optional

    # dependencies
    dependencies #::Array{fmi3Int32}
    dependenciesKind #::Array{fmi3String}

    # Constructor for not further specified Model variable
    function fmi3ModelVariable(name::String, valueReference::fmi3ValueReference)
        new(name, valueReference, fmi3DatatypeVariable(), "", fmi3CausalityLocal::fmi3Causality, fmi3VariabilityContinuous::fmi3Variability)
    end

    # Constructor for fully specified Model Variable
    function fmi3ModelVariable(name::String, valueReference::fmi3ValueReference, type, description, causalityString, variabilityString, dependencies, dependenciesKind)
        var = fmi3VariabilityDiscrete::fmi3Variability
        if type.datatype == fmi3Float32 || type.datatype == fmi3Float64
            var = fmi3VariabilityContinuous::fmi3Variability
        end
        cau = fmi3CausalityLocal::fmi3Causality

        # if !occursin("fmi3" * variabilityString, string(instances(fmi3Variability)))
        #     display("Error: variability not known")
        # else
        #     for i in 0:(length(instances(fmi3Variability))-1)
        #         if "fmi3" * variabilityString == string(fmi3Variability(i))
        #             var = fmi3Variability(i)
        #         end
        #     end
        # end

        # if !occursin("fmi3" * causalityString, string(instances(fmi3Causality)))
        #     display("Error: causalitiy not known")
        # else
        #     for i in 0:(length(instances(fmi3Causality))-1)
        #         if "fmi3" * causalityString == string(fmi3Causality(i))
        #             cau = fmi3Causality(i)
        #         end
        #     end
        # end

        new(name, valueReference, type, description, cau, var, dependencies, dependenciesKind)
    end
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.4.1. Definition of an FMU

The central FMU data structure defining all variables of the FMU that are visible/accessible via the FMU functions.
"""
mutable struct fmi3ModelDescription
    # FMI model description
    fmiVersion::String
    modelName::String
    generationTool::String
    generationDateAndTime::String
    variableNamingConvention::String
    instantiationToken::String  # replaces GUID

    CSmodelIdentifier::String
    CSneedsExecutionTool::Bool
    CScanBeInstantiatedOnlyOncePerProcess::Bool
    CScanGetAndSetFMUstate::Bool
    CScanSerializeFMUstate::Bool
    CSprovidesDirectionalDerivatives::Bool
    CSproivdesAdjointDerivatives::Bool
    CSprovidesPerElementDependencies::Bool
    CScanHandleVariableCommunicationStepSize::Bool
    CSmaxOutputDerivativeOrder::UInt
    CSprovidesIntermediateUpdate::Bool
    CSrecommendedIntermediateInputSmoothness::Int
    CScanReturnEarlyAfterIntermediateUpdate::Bool
    CShasEventMode::Bool
    CSprovidesEvaluateDiscreteStates::Bool

    MEmodelIdentifier::String
    MEneedsExecutionTool::Bool
    MEcanBeInstantiatedOnlyOncePerProcess::Bool
    MEcanGetAndSetFMUstate::Bool
    MEcanSerializeFMUstate::Bool
    MEprovidesDirectionalDerivatives::Bool
    MEprovidesAdjointDerivatives::Bool
    MEprovidesPerElementDependencies::Bool

    SEmodelIdentifier::String
    SEneedsExecutionTool::Bool
    SEcanBeInstantiatedOnlyOncePerProcess::Bool
    SEcanGetAndSetFMUstate::Bool
    SEcanSerializeFMUstate::Bool
    SEprovidesDirectionalDerivatives::Bool
    SEprovidesAdjointDerivatives::Bool
    SEprovidesPerElementDependencies::Bool

    # helpers
    isCoSimulation::Bool
    isModelExchange::Bool
    isScheduledExecution::Bool

    description::String

    # Model variables
    modelVariables::Array{fmi3ModelVariable,1}

    # additionals
    valueReferences::Array{fmi3ValueReference}

    inputValueReferences::Array{fmi3ValueReference}
    outputValueReferences::Array{fmi3ValueReference}
    stateValueReferences::Array{fmi3ValueReference}
    derivativeValueReferences::Array{fmi3ValueReference}
    intermediateUpdateValueReferences::Array{fmi3ValueReference}

    numberOfContinuousStates::Int
    numberOfEventIndicators::Int
    enumerations::fmi3Enum

    stringValueReferences

    defaultStartTime::fmi3Float64
    defaultStopTime::fmi3Float64
    defaultTolerance::fmi3Float64

    # Constructor for uninitialized struct
    fmi3ModelDescription() = new()

    # additional fields (non-FMI-specific)
    valueReferenceIndicies::Dict{Integer,Integer}
end

""" Overload the Base.show() function for custom printing of the fmi3ModelDescription"""
Base.show(io::IO, desc::fmi3ModelDescription) = print(io,
"Model name:     $(desc.modelName)
FMI version:     $(desc.fmiVersion)
GUID:            $(desc.instantiationToken)
Description:     $(desc.description)
Co-Simulation:   $(desc.isCoSimulation)
Model-Exchange:  $(desc.isModelExchange)
Model variables: $(desc.modelVariables)"
)


text:
function readCompoundNodes( mpSim::MPsim, xf::XF, catXF::XF )

    if !XLSX.hassheet( xf, "Compound Nodes" )
        @warn "Configuration file is missing a 'Compound Nodes' sheet, cannot configure compound nodes."
        return
    end  # if !XLSX.hassheet( xf, "Compound Nodes" )

    ws = xf["Compound Nodes"]
    clearSimulationCompoundNodes!( mpSim )

    if XLSX.hassheet( catXF, "States" )
        catalogue = catXF["States"]
        nodes = generateHierarchyNodes( mpSim, ws, catalogue )
        append!( nodes, readCatalogueNodes( mpSim, ws, catalogue ) )
    else
        nodes = generateHierarchyNodes( mpSim, ws )
    end  # if XLSX.hassheet( catXF, "States" )

    addSimulationCompoundNode!( mpSim, nodes... )
    readCustomNodes!( mpSim, ws )

end  # readCompoundNodes( mpSim, xf, catXF )

function generateHierarchyNodes( mpSim::MPsim, ws::WS, catalogue::WS )

    nodes = generateHierarchyNodes( mpSim, ws )

    if isempty( nodes )
        return nodes
    end  # if isempty( nodes )

    # Find names for compound nodes in catalogue.
    names = readCatalogueNodeNames( catalogue::WS )
    catSize = length( names )
    combs = Int.( catalogue[CR( 2, 7, catSize + 1, 7 )] )[:]

    for node in nodes
        attrVals = split.( split( node.name, "/" ), ":" )
        nReqs = length( attrVals )
        attrVals = string.( collect( Iterators.flatten( attrVals ) ) )
        eligibleLines = findall( combs .== nReqs )

        isMatch = map( eligibleLines ) do lineNr
            comb = string.( catalogue[CR( lineNr + 1, 8, lineNr + 1,
                nReqs * 2 + 7 )] )[:]
            return all( comb .== attrVals )
        end  # map( eligibleLines ) do lineNr

        matchNr = findfirst( isMatch )

        if !( matchNr isa Nothing )
            setCompoundNodeName!( node, names[eligibleLines[matchNr]] )
        end  # if !( matchNr isa Nothing )
    end  # for node in nodes

    return nodes

end  # generateHierarchyNodes( mpSim, ws, catalogue )

function generateHierarchyNodes( mpSim::MPsim, ws::WS )

    nLevels = ws["F4"] isa Integer ? ws["F4"] : 0

    if nLevels <= 0
        return Vector{CompoundNode}()
    end  # if nLevels <= 0

    # Get the list of attributes to generate the hierarchy.
    attributeNames = string.( ws[CR( 7, 5, 6 + nLevels, 5 )] )
    attributeNames = attributeNames[haskey.( Ref( mpSim.attributeList ),
        attributeNames )]
    attributes = get.( Ref( mpSim.attributeList ), attributeNames, missing )
    attributes = attributes[length.( getfield.( attributes,
        :possibleValues ) ) .> 1]

    if isempty( attributes )
        return Vector{CompoundNode}()
    end  # if isempty( attributes )

    unique!( attributes )

    # Get the list of all base nodes on each level and create compound nodes
    #   where needed.
    baseNodes = get.( Ref( mpSim.baseNodeList ),
        collect( keys( mpSim.baseNodeList ) ), missing )
    baseNodesByLevel = Vector{Array{Vector{BaseNode}}}( undef,
        length( attributes ) + 1 )
    compoundNodes = Vector{CompoundNode}()

    for ii in reverse( eachindex( baseNodesByLevel ) )
        # Generate base nodes by level.
        tmpAttrs = attributes[1:( ii > length( attributes ) ? end : ii )]
        op = ii > length( attributes ) ? Base.:(>=) : Base.:(==)
        tmpNodes = filter( node -> op( length( node.requirements ), ii ),
            baseNodes )
        attrCombs = collect( Iterators.product( getfield.( tmpAttrs,
            :possibleValues )... ) )

        baseNodesByLevel[ii] = map( attrCombs ) do attrComb
            return filter( tmpNodes ) do node
                return all( haskey.( Ref( node.requirements ),
                        attributeNames ) ) &&
                    all( get.( Ref( node.requirements ), attributeNames,
                        missing ) .== attrComb )
            end  # filter( baseNodes ) do node
        end  # map( attrCombs ) do attrComb

        if ii <= length( attributes )
            # Generate compound nodes.
            makeCompound = map( eachindex( baseNodesByLevel[ii] ) ) do jj
                return ( length( baseNodesByLevel[ii][jj] ) > 1 ) ||
                    ( length( baseNodesByLevel[ii + 1][jj] ) > 0 )
            end  # map( eachindex( baseNodesByLevel ) ) do jj
            makeCompound = findall( makeCompound )

            baseNodesByLevel[ii] = vcat.( baseNodesByLevel[ii],
                baseNodesByLevel[ii + 1] )

            compoundNames = map( makeCompound ) do jj
                return join( string.( attributeNames[1:ii], ":",
                    attrCombs[jj] ), "/" )
            end  # map( makeCompound ) do jj
            tmpCompounds = CompoundNode.( compoundNames )
            setCompoundNodeComponents!.( tmpCompounds,
                map( nodeList -> getfield.( nodeList, :name ),
                baseNodesByLevel[ii][makeCompound] ) )
            append!( compoundNodes, tmpCompounds )

            # Collapse base node array over last dimension.
            baseNodesByLevel[ii] = vcat.( map( jj -> getindex(
                baseNodesByLevel[ii], fill( Colon(), ii - 1 )..., jj ),
                axes( baseNodesByLevel[ii], ii ) )... )
            # ? This line creates a vector of dim-1 arrays of the length of the
            #   final dimension, and then performs an element-wise
            #   concatenation.
        end  # if ii < length( attributes )
    end  # for ii in reverse( eachindex( baseNodesByLevel ) )

    return compoundNodes

end  # generateHierarchyNodes( mpSim, ws )

function readCatalogueNodeNames( catalogue::WS )

    # Find names for compound nodes in catalogue.
    nStates = 2

    while !( catalogue[nStates, 1] isa Missing )
        nStates += 1
    end  #  while !( catalogue[nStates, 1] isa Missing )

    nStates -= 2
    return string.( catalogue[CR( 2, 1, nStates + 1, 1 )] )[:]

end  # readCatalogueNodeNames( catalogue )

function readCatalogueNodes( mpSim::MPsim, ws::WS, catalogue::WS )

    names = readCatalogueNodeNames( catalogue::WS )
    catSize = length( names )

    # Find the compound node definitions in the catalogue.
    nNodes = ws["B4"]

    if nNodes == 0
        return Vector{CompoundNode}()
    end  # if nNodes == 0

    nodeNames = string.( ws[CR( 7, 1, 6 + nNodes, 1 )] )[:]
    catLines = map( name -> findfirst( names .== name ), nodeNames )
    isInCat = .!isa.( catLines, Nothing )
    nodeNames = nodeNames[isInCat]
    catLines = catLines[isInCat] .+ 1
    nodes = Vector{CompoundNode}( undef, length( nodeNames ) )

    for ii in eachindex( nodeNames )
        jj = catLines[ii]
        nReqs = catalogue[jj, 7]
        reqs = Dict( catalogue[jj, kk * 2 + 6] =>
            catalogue[jj, kk * 2 + 7] for kk in 1:nReqs )

        baseNodes = filter( collect( keys( mpSim.baseNodeList ) ) ) do name
            node = mpSim.baseNodeList[name]
            return all( map( collect( keys( reqs ) ) ) do attr
                return haskey( node.requirements, attr ) &&
                    ( node.requirements[attr] == reqs[attr] )
            end  # map( collect( keys( reqs ) ) ) do attr
            )
        end  # filter( collect( keys( mpSim.baseNodeList ) ) ) do name

        if !isempty( baseNodes )
            nodes[ii] = CompoundNode( nodeNames[ii] )
            setCompoundNodeComponents!( nodes[ii], baseNodes )
        end  # if !isempty( baseNodes )
    end  # for ii in eachindex( nodeNames )

    return nodes[map( ii -> isassigned( nodes, ii ), 1:length( nodes ) )]

end  # readCatalogueNodes( mpSim, ws, catalogue )

function readCustomNodes!( mpSim::MPsim, ws::WS )

    nNodes = ws["I4"]

    if nNodes == 0
        return
    end  # if nNodes == 0

    nodeNames = string.( ws[CR( 7, 8, 6 + nNodes, 8 )] )[:]
    sLine = 6

    for name in nodeNames
        sLine += 1
        nComps = ws[sLine, 10]

        if !( nComps isa Integer ) || ( nComps <= 0 )
            continue
        end  # if !( nComps isa Integer ) || ...

        # Resolve all components to the base nodes.
        comps = string.( ws[CR( sLine, 11, sLine, 10 + nComps )] )[:]
        compNodes = filter( comp -> haskey( mpSim.baseNodeList, comp ), comps )
        otherNodes = filter( comp -> haskey( mpSim.compoundNodeList, comp ),
            comps )
        append!( compNodes,
            vcat( map( cNode -> mpSim.compoundNodeList[cNode].baseNodeList,
            otherNodes )... ) )

        # Add compound node if the list of base nodes isn't empty.
        if !isempty( compNodes )
            node = CompoundNode( name )
            setCompoundNodeComponents!( node, compNodes )
            addSimulationCompoundNode!( mpSim, node )
        end  # if !isempty( compNodes )
    end  # for ii in eachindex( nodeNames )

end  # readCustomNodes!( mpSim::MPsim, ws::WS )

function readCompoundNodes( mpSim::MPsim, configDB::SQLite.DB,
    configName::String )

    nodePars = DataFrame( DBInterface.execute( configDB,
        string( "SELECT * FROM `", configName,
        "` WHERE parType IS 'Compound Node'" ) ) )

    if isempty( nodePars )
        return
    end  # if isempty( nodePars )

    nodes = map( eachindex( nodePars[:, :parName] ) ) do ii
        node = CompoundNode( nodePars[ii, :parName] )
        pars = split( nodePars[ii, :parValue], ";" )

        if length( pars[1] ) > 2
            setCompoundNodeComponents!( node, string.(
                split( pars[1][2:(end-1)], "," ) ) )
        end  # if length( pars[1] ) > 2

        if tryparse( Int, pars[2] ) isa Int
            setCompoundNodeTarget!( node, parse( Int, pars[2] ) )
        end  # if tryparse( Int, pars[2] ) isa Int

        return node
    end  # map( eachindex( nodePars[:parName] ) ) do ii

    setSimulationCompoundNodes!( mpSim, nodes )

end  # readCompoundNodes( mpSim, configDB, configName )

ext: jl


text:
using MinimumVolumeEllipsoids
using Plots

X = [
    -1 -1 1 2
    1 -1 -1 2
]

ϵ = minimum_volume_ellipsoid(X)

U = rand(ϵ, 5000)

p = scatter(U[1, :], U[2, :]; legend=:none)
scatter!(p, X[1, :], X[2, :]; legend=:none)


text:
module FrostHelperIncrementBoosterModule

using ..Ahorn, Maple

@mapdef Entity "FrostHelper/IncrementBooster" FHIncrementalBooster(x::Integer, y::Integer, boostTime::Number=0.25, respawnTime::Number=1.0, particleColor::String="93bd40", directory="objects/FrostHelper/dashIncrementBooster/", reappearSfx::String="event:/game/04_cliffside/greenbooster_reappear", enterSfx::String="event:/game/04_cliffside/greenbooster_enter", boostSfx::String="event:/game/04_cliffside/greenbooster_dash", releaseSfx::String="event:/game/04_cliffside/greenbooster_end", red::Bool=false, dashCap::Integer = -1, dashes::Integer=1, refillBeforeIncrementing::Bool=false, preserveSpeed::Bool=false)

const colors = sort(collect(keys(Ahorn.XNAColors.colors)))

const placements = Ahorn.PlacementDict(
	"Dash Increment Booster (Frost Helper)" => Ahorn.EntityPlacement(
		FHIncrementalBooster
	),
	"Dash Increment Booster (Red, Frost Helper)" => Ahorn.EntityPlacement(
		FHIncrementalBooster,
        "point",
        Dict{String, Any}(),
        function(entity)
			entity.data["reappearSfx"] = "event:/game/05_mirror_temple/redbooster_reappear"
			entity.data["boostSfx"] = "event:/game/05_mirror_temple/redbooster_dash"
			entity.data["enterSfx"] = "event:/game/05_mirror_temple/redbooster_enter"
			entity.data["releaseSfx"] = "event:/game/05_mirror_temple/redbooster_end"
			entity.data["red"] = true
			entity.data["directory"] = "objects/FrostHelper/dashIncrementBoosterRed/"
			entity.data["particleColor"] = "c268d1"
			entity.data["dashes"] = 2
        end
	)
)

Ahorn.editingOptions(entity::FHIncrementalBooster) = Dict{String, Any}(
    "particleColor" => colors
)

sprite = "objects/booster/booster00"

function Ahorn.selection(entity::FHIncrementalBooster)
    x, y = Ahorn.position(entity)
    return Ahorn.getSpriteRectangle(sprite, x, y)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::FHIncrementalBooster, room::Maple.Room)
	texture = strip(get(entity.data, "directory", "objects/FrostHelper/blueBooster/"), '/')
	realTexture = "$texture/booster00.png"

	Ahorn.drawSprite(ctx, realTexture, 0, 0)
end

end

max_stars_repo_name: JaThePlayer/FrostHelper

max_issues_repo_name: JaThePlayer/FrostHelper

max_forks_repo_name: JaThePlayer/FrostHelper


text:
#useful macro for conditionally running things in parallel
macro mpar(cond, expr)
  quote
    if $(esc(cond))
        :($(Threads.@threads($expr)))
    else
        :($($expr))
    end
  end
end

const spinlock = Threads.SpinLock()

using Dates, DataFrames

function hypersort!(df::AbstractDataFrame, sorts::Vector{Symbol}, hyperidx::Vector{Int},
  threaded::Bool=false)

  #assume the other levels are sorted
  if length(sorts) > 1
    hypersort!(df, sorts[2:end], hsidx, threaded)
  end

  target::Symbol = sorts[1]
  T::Type = eltype(df[!,target])
  gdf::Vector{SubDataFrame} = collect(groupby(df, sorts[1]))
  targetv::Vector{Float64} = (sdf::SubDataFrame->sdf[1,target]).(gdf)
  Ntarget::Int = length(targetv)

  #map the values to the index array
  cumgroups::Vector{Int} = (sdf::SubDataFrame->size(sdf,1)).(gdf)
  cumsum!(cumgroups) #this is the upperbound on the index rows
  shyperidxs::Vector{SubArray{Int}} = (i->
    view(hyperidx, (i==1 ? 1 : cumgroups[i-1]):cumgroups[i])).(1:Ntarget)
  rowmapping::Dict{T, SubArray{Int}} = Dict{T, SubArray{Int}}(zip(targetv, hsidxs))

end
function hypersort(df::AbstractDataFrame, sorts::Vector{Symbol}; threaded::Bool=false)
  hyperidx::Vector{Int} = collect(1:size(df,1))
  hypersort!(df, sorts, hyperidx, threaded)
end

#------------------------------------------
#=function sortdf3!(df::DataFrame, sorts::Vector{Symbol}; threaded::Bool = false)

  firstsort::Symbol = sorts[1]
  sort!(df, firstsort)
  (length(sorts) == 1) && (return df)

  finalsorts::Vector{Symbol} = sorts[2:end]

  dfout::DataFrame = DataFrame(view(df, falses(size(df, 1)),:))
  sdfs = groupby(df, firstsort)
  dfouts::Vector{DataFrame} = (i->deepcopy(dfout)).(1:length(sdfs))
  @mpar threaded for i ∈ 1:length(sdfs)
    sdf::SubDataFrame = sdfs[i]
    dfouts[i] = sort(sdf, finalsorts)
  end

  (d->append!(dfout,d)).(dfouts)

  return dfout
end

using DataFrames, Dates, Finometrics
function testsorts(N::Int = 10_000_000)
  local df::DataFrame
  local dfs1::DataFrame
  local dfs2::DataFrame
  local dfs3::DataFrame

  df = DataFrame(ones(N,100))
  df.rid = rand(1:25000,N)
  df.dt = rand(Date(1980,1,1):Day(1):Date(2018,12,31),N)
  df.rid2 = rand(1:1000,N)

  dfs1 = deepcopy(df)
  print("\ntime sort1: ")
  @time sort!(dfs1, [:rid, :dt, :rid2])

  Base.GC.gc()

  dfs2 = deepcopy(df)
  print("time sort2: ")
  @time begin
    sort!(dfs2, :rid)
    dfs2s = groupby(dfs2, :rid)
    Threads.@threads for i ∈ 1:length(dfs2s)
      sdfs2::SubDataFrame = dfs2s[i]
      sdfs2 .= sort(sdfs2, [:dt, :rid2])
    end
  end

  dfs3 = deepcopy(df)
  print("time sort3: ")
  @time dfs3 = sortdf3!(dfs3, [:rid, :dt, :rid2], threaded=true)

  dfs4 = deepcopy(df)
  print("time sort4: ")
  @time dfs4 .= sort(dfs4, [:rid, :dt, :rid2])

  @assert sum(dfs1.rid .== dfs2.rid)==N
  @assert sum(dfs1.dt .== dfs2.dt)==N
  @assert sum(dfs1.rid2 .== dfs2.rid2)==N

  @assert sum(dfs1.rid .== dfs3.rid)==N
  @assert sum(dfs1.dt .== dfs3.dt)==N
  @assert sum(dfs1.rid2 .== dfs3.rid2)==N

  @assert sum(dfs1.rid .== dfs4.rid)==N
  @assert sum(dfs1.dt .== dfs4.dt)==N
  @assert sum(dfs1.rid2 .== dfs4.rid2)==N

  return nothing

end

testsorts(100_000)=#

using StatsModels, DataFrames, GLM
implicit_intercept(::Type{<:Any}) = true

#=function mwe(N=20)
  df = DataFrame(x = rand(N), z = (i->Symbol(:s, i)).(rand(1:3, N)))
  LHS::Nothing, RHS = nothing, :(x+z)

  #f = term(0) ~ StatsModels.terms!(StatsModels.parse!(:(x+y)))
  f = @eval(@formula(x ~ $RHS))
  #f = @eval(@formula(0 ~ $RHS))
  f = apply_schema(f, schema(f,df), StatisticalModel)
  m = modelcols(f.rhs, df)
  #mf::ModelFrame = ModelFrame(f, df)
  #m = ModelMatrix(mf).m
  display(m)
end

mwe()=#

function mwe(N=10)
  df = DataFrame(x = rand(N),
    zm = Vector{Union{Missing, Int}}(rand(1:4,N)),
    z = Vector{Int}(rand(1:4,N)))

  println("W/out missings")
  f = @eval(@formula(x ~ z))
  f = apply_schema(f, schema(f,df), StatisticalModel)
  m = modelcols(f.rhs, df)
  display(m)

  println("With missings")
  f = @eval(@formula(x ~ zm))
  f = apply_schema(f, schema(f,df), StatisticalModel)
  m = modelcols(f.rhs, df)
  display(m)
end

mwe()

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
function glavnaia(r::Robot)
    side=Nord
    for _ in 1:4
        change=1
        while change !=0
            change=move_full_detour(r,side)
            putmarker!(r)
        end
        while ismarker(r)
            move_full_detour(r,inverse(side))
        end
        side=inverse(next(side))
    end
    putmarker!(r)
end

function move_full_detour(r::Robot, side::HorizonSide)
    num_steps=0
    while isborder(r,side)&&!isborder(r,next(side))
        move!(r,next(side))
        num_steps+=1
    end
    change=0
    if !isborder(r,side)
        move!(r,side)
        change+=1
    end
    if num_steps !=0
        while isborder(r,inverse(next(side)))
            move!(r,side)
            change+=1
        end
        for _ in 1:num_steps
            move!(r,inverse(next(side)))
        end
    end
    return change
end

next(side::HorizonSide)=HorizonSide(mod(Int(side)+1,4))

inverse(side::HorizonSide)=HorizonSide(mod(Int(side)+2,4))

max_stars_repo_name: Aminulik/mirea-progs

max_issues_repo_name: Aminulik/mirea-progs

max_forks_repo_name: Aminulik/mirea-progs


text:
# Tensor network type hierarchy
#-------------------------------
abstract AbstractTensorNetwork{S<:IndexSpace,P<:TensorSpace,T}
# Any implementation of AbstractTensor should have method definitions for the
# same set of methods which are defined for the dense implementation Tensor
# defined in tensor.jl.

# tensor characteristics
spacetype{S}(::AbstractTensorNetwork{S})=S
spacetype{S<:IndexSpace}(::Type{AbstractTensorNetwork{S}})=S
spacetype{S<:IndexSpace,P<:TensorSpace}(::Type{AbstractTensorNetwork{S,P}})=S
spacetype{S<:IndexSpace,P<:TensorSpace,T}(::Type{AbstractTensorNetwork{S,P,T}})=S
spacetype{TN<:AbstractTensorNetwork}(::Type{TN})=spacetype(super(TN))

tensortype{S,P}(::AbstractTensorNetwork{S,P})=P
tensortype{S<:IndexSpace,P<:TensorSpace}(::Type{AbstractTensorNetwork{S,P}})=P
tensortype{S<:IndexSpace,P<:TensorSpace,T}(::Type{AbstractTensorNetwork{S,P,T}})=P
tensortype{TN<:AbstractTensorNetwork}(::Type{TN})=spacetype(super(TN))

Base.eltype{S,P,T}(::AbstractTensorNetwork{S,P,T}) = T
Base.eltype{S<:IndexSpace,P<:TensorSpace,T}(::Type{AbstractTensorNetwork{S,P,T}})=T
Base.eltype{TN<:AbstractTensorNetwork}(::Type{TN}) = eltype(super(TN))


text:
## These methods are described in Appendix C of
## Lee, A., Doucet, A. and Łatuszyński, K., 2014. Perfect simulation using
## atomic regeneration with application to Sequential Monte Carlo. arXiv:1407.5770

function _signRejectionSample(μ::F1, φ::F2, c::Float64, rng::RNG) where
  {F1<:Function, F2<:Function, RNG<:AbstractRNG}
  while true
    v = φ(μ())
    if rand(rng) < abs(v)/c
      return sign(v) < 0
    end
  end
end

function _annotateμ(μ::F) where F<:Function
  calls = Ref(0)
  function ν()
    calls.x +=1
    return μ()
  end
  function νCalls()
    return calls.x
  end
  return ν, νCalls
end

function _signedEstimate(μ::F1, φ::F2, c::Float64, δ::Float64, n::Int64,
  rng::RNG) where {F1<:Function, F2<:Function, RNG<:AbstractRNG}
  ν, νCalls = _annotateμ(μ)
  q() = _signRejectionSample(ν, φ, c, rng)
  Y, flips = linear(q, 2.0, δ/c)
  if Y
    return 0.0, flips, νCalls()
  else
    v = 0.0
    for i in 1:n
      v += abs(φ(μ()))
    end
    return v/n, flips, νCalls()+n
  end
end

function _signedEstimate(μ::F1, φ::F2, a::Float64, b::Float64, δ::Float64,
  c::Float64, n::Int64, rng::RNG) where {F1<:Function, F2<:Function,
  RNG<:AbstractRNG}
  φnew(x) = φ(x) - b
  cnew = max(b-a, c-b)
  estimate, flips, calls = signedEstimate(μ, φnew, cnew, δ-b, n, rng)
  return estimate + b, flips, calls+n
end


text:
using jInv.Mesh;
using Helmholtz
using Helmholtz.PointSourceADR
using Multigrid

plotting = true;

if plotting
	using PyPlot;
	close("all");
end

include("getModels.jl");

NeumannAtFirstDim = false;
Sommerfeld = true;
gamma_val = 0.01;

model = "const";
domain = [0.0,1.0,0.0,1.0];
n_cells = [512,512];
src = div.(n_cells,2).+1;
# src[2] = 1;

if model=="const"
	NeumannAtFirstDim = false;
	Sommerfeld = true;
	f = 40.0;
	src = div.(n_cells,2).+1;
	m_coarse = getConstModel(n_cells.+1);
	pad_cells = [20;20];
else
	error("unknown model");
end
Minv = getRegularMesh(domain,n_cells);

w = 2*pi*f

if plotting
	figure()
	imshow(m_coarse'); colorbar();
	title(string("The ",model," model"))
end

n_tup = tuple((n_cells.+1)...);
q_coarse = zeros(ComplexF64,n_tup)
src2 = div.(n_cells,2).+1;
q_coarse[src2[1],src2[2]] = 1/(Minv.h[2]^2);

# src2 = div(3*n_cells,4)+1;
# q_coarse[src2[1],src2[2]] = 1/(Minv.h[2]^2);
# src = src2;

println("omega*h:");
println(w*Minv.h*sqrt(maximum(m_coarse)));

maxOmega = getMaximalFrequency(m_coarse,Minv);
ABLamp = maxOmega;

m = m_coarse;

gamma = (w*gamma_val)*ones(Float32,size(m)) + getABL(Minv.n.+1,NeumannAtFirstDim,pad_cells.+1,ABLamp);

if plotting
	figure();
	imshow(reshape(gamma,n_tup));
end

# getHelmholtzADR(central::Bool, Mesh::RegularMesh, mNodal::Array{Float64},omega::Float64, gamma::Array,T,G,LT,
									# NeumannAtFirstDim::Bool,src::Array{Int64},Sommerfeld::Bool,single::Bool = true)

X1,X2 = ndgrid(collect(0:Minv.n[1])*Minv.h[1],collect(0:Minv.n[2])*Minv.h[2]);

T  = X1;
Tx = ones(size(X1));
Ty = zeros(size(X1));
# G = [Tx,Ty];
G =  Array{Array{Float64}}(undef, Minv.dim);
G[1] = Tx;
G[2] = Ty;
LT = zeros(size(X1));
(ADR_long,ADR_short,T) = getHelmholtzADR(true,Minv, m, w, gamma,T,G,LT, NeumannAtFirstDim,[3,3],Sommerfeld);

# ADR_short = 0.5*ADR_long + 0.5*ADR_short;

if plotting
	figure();
	imshow(real(exp.(1im*w*reshape(T,n_tup)))');
end

H, = GetHelmholtzOperator(Minv, m, w, gamma_val*ones(Float32,size(m)),NeumannAtFirstDim,pad_cells.+1,ABLamp,Sommerfeld);
S =  GetHelmholtzShiftOP(m, w,0.2);

q = q_coarse;
xh = (H\q[:])
q = q[:].*exp.(1im*w*T[:]);
xs = ((ADR_short )\q[:]).*exp.(-1im*w*T[:]);
us = (H+S)\q[:];

# if plotting
# figure();
# subplot(1,4,1);
# imshow(real(I2'));
# subplot(1,4,2);
# imshow(reshape(T,n_tup)');
# subplot(1,4,3);
# imshow(reshape(real(alu),n_tup)');
# subplot(1,4,4);
# imshow(reshape(real(exp.(-1im*w*T)),n_tup)');

if plotting
	figure();
	subplot(1,3,1);
	imshow(reshape(real(xh),n_tup)',cmap = "jet"); title("Full Helmholtz solution"); colorbar();
	subplot(1,3,2);
	imshow(reshape(real(-us),n_tup)',cmap = "jet"); title("Shifted Helmholtz solution"); colorbar();
	subplot(1,3,3);
	imshow(reshape(real(xs),n_tup)',cmap = "jet"); title("First order upwind ADR solution");colorbar();
end
a = 1;

levels      = 6;
numCores 	= 2;
maxIter     = 10;
relativeTol = 1e-10;
relaxType   = "Jac-GMRES";
relaxParam  = 0.7;
relaxPre 	= 3;
relaxPost   = 3;
cycleType   ='W';
coarseSolveType = "NoMUMPS";

MG = getMGparam(ComplexF64,Int64,levels,numCores,maxIter,relativeTol,relaxType,relaxParam,relaxPre,relaxPost,cycleType,coarseSolveType,0.5,0.0);
nrhs = 1;
MGsetup(sparse(ADR_short'),Minv,MG,nrhs,true);

println("****************************** Stand-alone GMG RAP for ADR 1st: ******************************")
a = zeros(ComplexF64, size(q[:]));
solveMG(MG,q,a,true);

println("****************************** GMRES GMG RAP for ADR 1st: ******************************")
F = getAfun(sparse(ADR_short'),zeros(eltype(q),size(q)),1)
MG.maxOuterIter = 2;
solveGMRES_MG(F,MG,q,a,true,10,true);

z = 1;

# ADRparam = getADRparam(Minv,gamma,m,w,T,G,LT,src,NeumannAtFirstDim,true,M);
	# getOperator = (mesh,ADRparam)->((ADR2,ADR1,~) = getHelmholtzADR(true,ADRparam.Minv, ADRparam.m, ADRparam.w, ADRparam.gamma,
									# ADRparam.T,ADRparam.G,ADRparam.LT,ADRparam.NeumannAtFirstDim,ADRparam.src,ADRparam.Sommerfeld);
									# H = ADR2 + GetHelmholtzShiftOP(ADRparam.m, ADRparam.w,shift);
									# # M = exp(-1im*w*T);
									# # M = ADRparam.scaling;
									# # rescaleMatrix(H,M);
									# return H);

	# MG = getMGparam(levels,numCores,maxIter,relativeTol/relres,relaxType,relaxParam,
				# relaxPre,relaxPost,cycleType,coarseSolveType,0.5,0.0,"FullWeighting");
	# MGsetup(getMultilevelOperatorConstructor(ADRparam,getOperator,restrictADRParam),ADRparam.Minv,MG,eltype(q),1,true);


text:
"""
  `generateVTK(filename, points; lines, cells, point_data, path, num, time)`
Generates a vtk file with the given data. Written by Eduardo Alvarez.
  **Arguments**
  * `points::Array{Array{Float64,1},1}`     : Points to output.
  **Optional Arguments**
  * `lines::Array{Array{Int64,1},1}`  : line definitions. lines[i] contains the
                            indices of points in the i-th line.
  * `cells::Array{Array{Int64,1},1}`  : VTK polygons definiton. cells[i]
                            contains the indices of points in the i-th polygon.
  * `data::Array{Dict{String,Any},1}` : Collection of data point fields in the
                            following format:
                              data[i] = Dict(
                                "field_name" => field_name::String
                                "field_type" => "scalar" or "vector"
                                "field_data" => point_data
                              )
                            where point_data[i] is the data at the i-th point.
See `examples.jl` for an example on how to use this function.
"""
function generateVTK(filename::String, points;
                    lines::Array{Array{Int64,1},1}=Array{Int64,1}[],
                    cells::Array{Array{Int64,1},1}=Array{Int64,1}[],
                    point_data=nothing, cell_data=nothing,
                    num=nothing, time=nothing,
                    path="", comments="", _griddims::Int64=-1,
                    keep_points::Bool=false,
                    override_cell_type::Int64=-1)
  aux = num!=nothing ? ".$num" : ""
  ext = aux*".vtk"
  if path !=""
    _path = string(path, (path[end]!="/" ? "/" : ""))
  else
    _path = path
  end
  f = open(string(_path, filename, ext), "w")
  # HEADER
  header = "# vtk DataFile Version 4.0" # File version and identifier
  header = string(header, "\n", " ", comments) # Title
  header = string(header, "\n", "ASCII") # File format
  header = string(header, "\n", "DATASET UNSTRUCTURED_GRID")
  write(f, header)
  # TIME
  if time!=nothing
    line0 = "\nFIELD FieldData 1"
    line1 = "\nSIM_TIME 1 1 double"
    line2 = "\n$(time)"
    write(f, line0*line1*line2)
  end
  np = size(points)[1]
  nl = size(lines)[1]
  nc = size(cells)[1]
  _keep_points = keep_points || (nl==0 && nc==0)
  # POINTS
  write(f, string("\n", "POINTS ", np, " float"))
  for i in 1:np
    print(f, "\n", points[i][1], " ", points[i][2], " ", points[i][3])
  end
  # We do this to avoid outputting points as cells if outputting a Grid
  # or if we simply want to ignore points
  if _griddims!=-1 || !_keep_points
    auxnp = np
    np = 0
  end
  # CELLS
  auxl = size(lines)[1]
  for line in lines
    auxl += size(line)[1]
  end
  auxc = size(cells)[1]
  for cell in cells
    auxc += size(cell)[1]
  end
  write(f, "\n\nCELLS $(np+nl+nc) $(2*np+auxl+auxc)")
  for i in 1:np+nl+nc
    if i<=np
      pts = [i-1]
    elseif i<=np+nl
      pts = lines[i-np]
    else
      pts = cells[i-(nl+np)]
    end
    print(f, "\n", size(pts,1))
    for pt in pts
      print(f, " ", pt)
    end
  end
  write(f, "\n\nCELL_TYPES $(np+nl+nc)")
  for i in 1:np+nl+nc
    if i<=np
      tpe = 1
    elseif i<=np+nl
      tpe = 4
    else
      if override_cell_type==-1
        if _griddims!=-1
          if _griddims==1
            tpe = 3
          elseif _griddims==2
            tpe = 9
          elseif _griddims==3
            tpe = 12
          else
            error("Generation of VTK cells of $_griddims dimensions not implemented")
          end
        else
          tpe = 7
        end
      else
        tpe = override_cell_type
      end
    end
    print(f, "\n", tpe)
  end
  if _griddims!=-1 || !_keep_points
    np = auxnp
  end
  # POINT DATA
  if point_data!=nothing
      write(f, "\n\nPOINT_DATA $np")
  end
  _p_data = point_data!=nothing ? point_data : []
  for field in _p_data
    field_name = field["field_name"]
    field_type = field["field_type"]
    data = field["field_data"]
    if size(data)[1]!=np
      warn("Corrupted field $(field_name)! Field size != number of points.")
    end
    if field_type=="scalar"
      write(f, "\n\nSCALARS $field_name float\nLOOKUP_TABLE default")
      for entry in data
        print(f, "\n", entry)
      end
    elseif field_type=="vector"
      write(f, "\n\nVECTORS $field_name float")
      for entry in data
        print(f, "\n", entry[1], " ", entry[2], " ", entry[3])
      end
    else
      error("Unknown field type $(field_type).")
    end
  end
    # CELL DATA
    if cell_data!=nothing
        write(f, "\n\nCELL_DATA $nc")
    end
    _c_data = cell_data!=nothing ? cell_data : []
    for field in _c_data
      field_name = field["field_name"]
      field_type = field["field_type"]
      data = field["field_data"]
      if size(data)[1]!=nc
        warn("Corrupted field $(field_name)! Field size != number of cells.")
      end
      if field_type=="scalar"
        write(f, "\n\nSCALARS $field_name float\nLOOKUP_TABLE default")
        for entry in data
          print(f, "\n", entry)
        end
      elseif field_type=="vector"
        write(f, "\n\nVECTORS $field_name float")
        for entry in data
          print(f, "\n", entry[1], " ", entry[2], " ", entry[3])
        end
      else
        error("Unknown field type $(field_type).")
      end
    end
  close(f)
  return filename*ext*";"
end


text:
using JGE:unsafe_pointer, @safe_contain

function test_function_call(args...)
    return args
end

@noinline function test_case_nomacro()
    x = Ref(Cint(1))
    y = Ref(Ptr{Cvoid}())
    test_function_call(unsafe_pointer(x), unsafe_pointer(y))
end

@noinline function test_case()
    x = Ref(Cint(1))
    y = Ref(Ptr{Cvoid}())
    @safe_contain test_function_call(unsafe_pointer(x), unsafe_pointer(y))
end

@noinline function test_impl(ptrs, expected_values, fail)
    for (ptrs, expected) in zip(ptrs, expected_values)
        comp = unsafe_load(ptrs) == expected[]
        if fail
            @test !comp
        else
            @test comp
        end
    end
end

function test_feature_nomacro()
    # pointed values are not kept by Julia, so unsafe_load does not return the expected value
    test_impl(test_case_nomacro(), [Ref(Cint(1)), Ref(Ptr{Cvoid}())], true)
end

function test_feature()
    # pointed values kept in the deps field of the container object, so unsafe_load returns the expected value
    test_impl(test_case()[], [Ref(Cint(1)), Ref(Ptr{Cvoid}())], false)
end

test_feature_nomacro()
test_feature()


text:
function add_friends_neighbors_of_neighbors!(
    state::Tuple{AbstractGraph, AbstractArray},
    agent_idx::Integer,
    config::Config,
    new_friends_count::Int64,
    rng::MersenneTwister
)
    graph, agent_list = state
    this_agent = agent_list[agent_idx]
    # neighbors of neighbors
    friend_candidates = Int64[]

    for neighbor in inneighbors(graph, agent_idx)
        append!(
            friend_candidates,
            setdiff(
                all_neighbors(graph, neighbor),
                inneighbors(graph, agent_idx),
                agent_idx
            )
        )

    end

    shuffle!(rng, friend_candidates)
    # order neighbors by frequency of occurence in input_candidates descending
    if length(friend_candidates) > 0
        friends_queue = first.(
            sort(collect(countmap(friend_candidates)), by=last, rev=true)
        )

        if (length(friends_queue) - config.network.new_follows) < 0
            new_friends_count = length(friends_queue)

        for _ in 1:new_friends_count
            new_neighbor = popfirst!(friends_queue)
            add_edge!(graph, new_neighbor, agent_idx)
            end
        end
    end

    return state
end

function add_friends_random!(
    state::Tuple{AbstractGraph, AbstractArray},
    agent_idx::Integer,
    config::Config,
    new_friends_count::Int64
    )

    graph, agent_list = state
    this_agent = agent_list[agent_idx]

    friends_queue = Array{Tuple{Int64,Int64}, 1}()
    not_neighbors = setdiff(
        [1:(agent_idx - 1); (agent_idx + 1):nv(graph)],
        inneighbors(graph, agent_idx)
    )

        for candidate in not_neighbors
            if (
                abs(this_agent.opinion - agent_list[candidate].opinion)
                < config.opinion_threshs.befriend
            )
                push!(friends_queue, (candidate, outdegree(graph, candidate)))
        end
    end

    friends_queue = first.(sort(friends_queue, by=last, rev=true))

    if (length(friends_queue) - config.network.new_follows) < 0
        new_friends_count = length(friends_queue)

        for _ in 1:new_friends_count
            new_neighbor = popfirst!(friends_queue)
            add_edge!(graph, new_neighbor, agent_idx)
        end
    end

    return state
end

# suppress output of include()
;


text:
###############################################################################
#
# leaf level simulations of PAR and APAR
#
###############################################################################
"""
    leaf_fluxes(leaf::LeafBios{FT},
                in_rad::IncomingRadiation{FT},
                wls::WaveLengths{FT},
                rt_con::RTCache{FT}
    ) where {FT<:AbstractFloat}

Return leaf PAR and APAR, given
- `leaf` [`LeafBio`](@ref) type struct
- `in_rad` [`IncomingRadiation`](@ref) type struct
- `wls` [`WaveLengths`](@ref) type struct
- `rt_con` [`RTCache`](@ref) type cache

Note that `in_rad` assumes direct light with zenith angle of 0, and a zenith
    angle correction needs to be made before passing it to this function.
"""
function leaf_fluxes(
            leaf::LeafBios{FT},
            in_rad::IncomingRadiation{FT},
            wls::WaveLengths{FT},
            rt_con::RTCache{FT}
) where {FT<:AbstractFloat}
    # unpack variables
    @unpack dWL_iPAR, iPAR, WL_iPAR = wls;
    cf_con = rt_con.cf_con;
    cf_con.kChlrel .= view(leaf.kChlrel, iPAR);

    # PAR energy from direct light
    cf_con.E_iPAR .= view(in_rad.E_direct, iPAR) .* view(leaf.α_SW, iPAR);
    e2phot!(WL_iPAR, cf_con.E_iPAR, cf_con.PAR_dir);

    # PAR energy from diffuse light
    cf_con.E_iPAR .= view(in_rad.E_diffuse, iPAR) .* view(leaf.α_SW, iPAR);
    e2phot!(WL_iPAR, cf_con.E_iPAR, cf_con.PAR_diff);

    # absorbed PAR energy from direct and diffuse light
    cf_con.PAR_diffCab .= cf_con.kChlrel .* cf_con.PAR_diff;
    cf_con.PAR_dirCab  .= cf_con.kChlrel .* cf_con.PAR_dir;

    # toral PAR and APAR
    _dif    = numerical∫(cf_con.PAR_diff   , dWL_iPAR);
    _dir    = numerical∫(cf_con.PAR_dir    , dWL_iPAR);
    _difCab = numerical∫(cf_con.PAR_diffCab, dWL_iPAR);
    _dirCab = numerical∫(cf_con.PAR_dirCab , dWL_iPAR);

    # return PAR and APAR in μmol m⁻² s⁻¹
    return 1000 * (_dir + _dif), 1000 * (_dirCab + _difCab)
end


text:
import Simplex: nrows, ncols

@testset "Artificial.jl" begin
    """
    Canonical form:

    min  -10 x1 - 12 x2 - 12 x3
    s.t. - x1 - 2 x2 - 2 x3 - x4                == -20
        2 x1 +   x2 + 2 x3      + x5           ==  20
        2 x1 + 2 x2 +   x3           - x6      ==  10
        - x1 -   x2                       - x7 == -10
        x1, x2, x4, x5, x7 >= 0
        -1 <= x3 <= 10
        0 <= x6 <= 10
    """
    A = sparse(
        [1; 2; 3; 4; 1; 2; 3; 4; 1; 2; 3; 1; 2; 3; 4],
        [1; 1; 1; 1; 2; 2; 2; 2; 3; 3; 3; 4; 5; 6; 7],
        Float64[-1; 2; 2; -1; -2; 1; 2; -1; -2; 2; 1; -1; 1; -1; -1])
    b = Float64[-20; 20; 10; -10]
    c = Float64[-10; -12; -12; 0; 0; 0; 0]
    xl = Float64[0; 0; -1; 0; 0; 0; 0]
    xu = Float64[Inf; Inf; 10; Inf; Inf; 10; Inf]
    senses = fill(MatOI.EQUAL_TO, 4)
    lp = MatOI.LPSolverForm(MOI.MIN_SENSE, c, A, b, senses, xl, xu)

    @testset "reformulate" begin
        """
        Phase-one form:
        min                                           x8 + x9 + x10 + x11
        s.t. - x1 - 2 x2 - 2 x3 - x4                - x8                  ==  -20
            2 x1 +   x2 + 2 x3      + x5                + x9             ==  20
            2 x1 + 2 x2 +   x3           - x6                + x10       ==  10
            - x1 -   x2                       - x7                 - x11 ==  -10
            x1, x2, x4, x5, x7, x8, x9, x10, x11 >= 0
            -1 <= x3 <= 10
            0 <= x6 <= 10
        """
        p1, artif_col_idx = Simplex.PhaseOne.Artificial.reformulate(lp)
        @test nrows(p1) == 4
        @test ncols(p1) == 11
        @test p1.A == sparse(
            [1; 2; 3; 4; 1; 2; 3; 4; 1; 2; 3; 1; 2; 3; 4; 1; 2; 3; 4],
            [1; 1; 1; 1; 2; 2; 2; 2; 3; 3; 3; 4; 5; 6; 7; 8; 9; 10; 11],
            Float64[-1; 2; 2; -1; -2; 1; 2; -1; -2; 2; 1; -1; 1; -1; -1; -1; 1; 1; -1])
        @test p1.b == [-20; 20; 10; -10]
        @test p1.c == [0; 0; 0; 0; 0; 0; 0; 1; 1; 1; 1]
        @test p1.v_lb == [0; 0; -1; 0; 0; 0; 0; 0; 0; 0; 0]
        @test p1.v_ub == [Inf; Inf; 10; Inf; Inf; 10; Inf; Inf; Inf; Inf; Inf]

        # The original problems should not be modified.
        @test nrows(lp) == 4
        @test ncols(lp) == 7
        @test lp.A == A
        @test lp.b == b
        @test lp.c == c
        @test lp.v_lb == xl
        @test lp.v_ub == xu
    end

    lp = nothing
end


text:
module CellMapApplyTests

using Test
using CellwiseValues
using TensorValues

using ..CellValuesMocks
using ..MapsMocks

l = 10

a = VectorValue(10,10)
b = VectorValue(15,20)
p1 = VectorValue(1,1)
p2 = VectorValue(2,2)
p3 = VectorValue(3,3)
p = [p1,p2,p3]
m = MockMap(a)
r = evaluate(m,p)

cm = TestIterCellValue(m,l)
cp = TestIterCellValue(p,l)
rm = [ CachedArray(r) for i in 1:l]

test_iter_cell_map(cm,cp,rm)

cm2 = apply(-,cm,broadcast=true)
rm = [ CachedArray(-r) for i in 1:l]

test_iter_cell_map(cm2,cp,rm)

cm = TestIndexCellValue(m,l)
cp = TestIndexCellValue(p,l)
rm = [ CachedArray(r) for i in 1:l]

test_index_cell_map_with_index_arg(cm,cp,rm)

ca2 = evaluate(cm,cp)

test_index_cell_array(ca2,rm)

end # module


text:
module INIBuilder

using REPL
using REPL.TerminalMenus
using Circuitscape

include("filepicker.jl")
include("run.jl")

end


text:
abstract type AbstractInterval end

struct NoBoundsInterval <: AbstractInterval
end

struct LowerAndUpperBoundInterval <: AbstractInterval
    left::String
    right::String
    function LowerAndUpperBoundInterval(
            left::String,
            right::String,
            )::LowerAndUpperBoundInterval
        correct_left = strip(left)
        correct_right = strip(right)
        result::LowerAndUpperBoundInterval = new(
            correct_left,
            correct_right,
            )
        return result
    end
end

struct LowerBoundOnlyInterval <: AbstractInterval
    left::String
    function LowerBoundOnlyInterval(
            left::String,
            )::LowerBoundOnlyInterval
        correct_left = strip(left)
        result::LowerBoundOnlyInterval = new(
            correct_left,
            )
        return result
    end
end

struct UpperBoundOnlyInterval <: AbstractInterval
    right::String
    function UpperBoundOnlyInterval(
            right::String,
            )::UpperBoundOnlyInterval
        correct_right = strip(right)
        result::UpperBoundOnlyInterval = new(
            correct_right,
            )
        return result
    end
end

function is_interval(x::String)::Bool
    if is_no_bounds_interval(x)
        return true
    elseif is_lower_bound_only_interval(x)
        return true
    elseif is_upper_bound_only_interval(x)
        return true
    elseif is_lower_and_upper_bound_interval(x)
        return true
    else
        return false
    end
end

function get_lower_and_upper_bound_interval_regex()::Regex
    lower_and_upper_bound_interval_regex::Regex =
        r"\[(\w\w*?)\,(\w\w*?)\)"
    return lower_and_upper_bound_interval_regex
end

function get_lower_bound_only_interval_regex()::Regex
    lower_bound_only_interval_regex::Regex =
        r"\[(\w\w*?)\,\)"
    return lower_bound_only_interval_regex
end

function get_upper_bound_only_interval_regex()::Regex
    upper_bound_only_interval_regex::Regex =
        r"\[\,(\w\w*?)\)"
    return upper_bound_only_interval_regex
end

function get_no_bounds_interval_regex()::Regex
    no_bounds_interval_regex::Regex =
        r"\[\,\)"
    return no_bounds_interval_regex
end

function is_no_bounds_interval(x::String)::Bool
    result::Bool = occursin(
        get_no_bounds_interval_regex(),
        x,
        )
    return result
end

function is_lower_and_upper_bound_interval(x::String)::Bool
    result::Bool = occursin(
        get_lower_and_upper_bound_interval_regex(),
        x,
        )
    return result
end

function is_lower_bound_only_interval(x::String)::Bool
    result::Bool = occursin(
        get_lower_bound_only_interval_regex(),
        x,
        )
    return result
end

function is_upper_bound_only_interval(x::String)::Bool
    result::Bool = occursin(
        get_upper_bound_only_interval_regex(),
        x,
        )
    return result
end

function construct_interval(x::String)::AbstractInterval
    if is_no_bounds_interval(x)
        result = NoBoundsInterval()
    elseif is_lower_bound_only_interval(x)
        loweronly_regexmatch::RegexMatch = match(
            get_lower_bound_only_interval_regex(),
            x,
            )
        loweronly_left::String = strip(
            convert(String, loweronly_regexmatch[1])
            )
        result = LowerBoundOnlyInterval(loweronly_left)
    elseif is_upper_bound_only_interval(x)
        upperonly_regexmatch::RegexMatch = match(
            get_upper_bound_only_interval_regex(),
            x,
            )
        upperonly_right::String = strip(
            convert(String, upperonly_regexmatch[1])
            )
        result = UpperBoundOnlyInterval(upperonly_right)
    elseif is_lower_and_upper_bound_interval(x)
        lowerandupper_regexmatch::RegexMatch = match(
            get_lower_and_upper_bound_interval_regex(),
            x,
            )
        lowerandupper_left::String = strip(
            convert(String, lowerandupper_regexmatch[1])
            )
        lowerandupper_right::String = strip(
            convert(String, lowerandupper_regexmatch[2])
            )
        result = LowerAndUpperBoundInterval(
            lowerandupper_left,
            lowerandupper_right,
            )
    else
        error("argument is not a valid interval")
    end
    return result
end

function interval_contains_x(
        interval::NoBoundsInterval,
        x::AbstractString,
        )::Bool
    result::Bool = true
    return result
end

function interval_contains_x(
        interval::LowerAndUpperBoundInterval,
        x::AbstractString,
        )::Bool
    x_stripped::String = strip(convert(String, x))
    left::String = strip(interval.left)
    right::String = strip(interval.right)
    result::Bool = (left <= x_stripped) && (x_stripped < right)
    return result
end

function interval_contains_x(
        interval::LowerBoundOnlyInterval,
        x::AbstractString,
        )::Bool
    x_stripped::String = strip(convert(String, x))
    left::String = strip(interval.left)
    result::Bool = left <= x_stripped
    return result
end

function interval_contains_x(
        interval::UpperBoundOnlyInterval,
        x::AbstractString,
        )
    x_stripped::String = strip(convert(String, x))
    right::String = strip(interval.right)
    result::Bool = x_stripped < right
    return result
end


text:
# memory operations

function Base.unsafe_copyto!(ctx::ZeContext, dev::ZeDevice, dst::Union{Ptr{T},ZePtr{T}},
                             src::Union{Ptr{T},ZePtr{T}}, N::Integer) where T
    bytes = N*sizeof(T)
    bytes==0 && return
    execute!(global_queue(ctx, dev)) do list
        append_copy!(list, dst, src, bytes)
    end
end

function unsafe_fill!(ctx::ZeContext, dev::ZeDevice, ptr::Union{Ptr{T},ZePtr{T}},
                      pattern::Union{Ptr{T},ZePtr{T}}, N::Integer) where T
    bytes = N*sizeof(T)
    bytes==0 && return
    execute!(global_queue(ctx, dev)) do list
        append_fill!(list, ptr, pattern, sizeof(T), bytes)
    end
end


text:
include("core.jl")
include("time.jl")
include("params.jl")
include("defs.jl")
include("instances.jl")
include("model.jl")


text:
using Pkg
Pkg.activate(pwd())

# # Useful packages
# ## Statistics

using Statistics

x = rand(10);

mean(x)
var(x)
std(x)

# ## LinearAlgebra

using LinearAlgebra

A = [-4.0 -17.0; 2.0 2.0]

det(A)
inv(A)
norm(A)
eigvals(A)
eigvecs(A)

#+

D = Diagonal([1,2,3])

#+

D + I
D - I

# ## Random

using Random
using Random: seed!

seed!(1234);
rand(2)

#+

seed!(1234);
rand(2)

#+

randperm(4)

#+

v = [1,2,3,4]
shuffle(v)


text:
#################### Flat Distribution ####################

struct Flat <: ContinuousUnivariateDistribution end

minimum(d::Flat) = -Inf
maximum(d::Flat) = Inf
insupport(d::Flat, x::Real) = true

logpdf(d::Flat, x::Real) = 0.0

function Truncated(d::Flat, l::Float64, u::Float64)
  Truncated{Flat, Continuous}(d, l, u, 0.0, 1.0, 1.0, 0.0)
end

#################### BDiagNormal Distribution ####################

const BDiagNormal = MvNormal{Float64, PBDiagMat, Vector{Float64}}

function BDiagNormal(μ::Vector{Float64}, Σ::Matrix{Float64})
  n = div(length(μ), size(Σ, 1))
  BDiagNormal(μ, PBDiagMat(Σ, n))
end

function BDiagNormal(μ::Vector{Float64}, Σ::Vector{Matrix{Float64}})
  BDiagNormal(μ, PBDiagMat(Σ))
end

BDiagNormal(μ::AbstractVector, Σ::AbstractMatrix) =
  BDiagNormal(convert(Vector{Float64}, μ), convert(Matrix{Float64}, Σ))

BDiagNormal(μ::AbstractVector, Σ::AbstractVector) =
  BDiagNormal(convert(Vector{Float64}, μ), Matrix{Float64}[Σ...])

#################### Null Distribution ####################

struct NullUnivariateDistribution <: UnivariateDistribution{ValueSupport} end

#################### SymUniform Distribution ####################

struct SymUniform <: ContinuousUnivariateDistribution
  SymUniform() = Uniform(-1.0, 1.0)
  SymUniform(μ::Real, σ::Real) = Uniform(μ - σ, μ + σ)
end

#################### Type Aliases ####################

const SymDistributionType =
          Union{Type{Biweight}, Type{Cosine}, Type{Epanechnikov}, Type{Normal},
                Type{SymTriangularDist}, Type{Triweight}, Type{SymUniform}}

const KernelDensityType = SymDistributionType


text:
abstract type AbstractDecisionRule{S<:Grid,T<:Grid,nx} end

function Base.show(io::IO, dr::AbstractDecisionRule)
    println(io, typeof(dr))
end

outdim(dr::AbstractDecisionRule{<:Grid,<:Grid,nx}) where nx = nx

import Base

function MSM(v::AbstractVector{<:AbstractVector{T}}) where T
    x = cat(v...;dims=1)
    sizes = [length(e) for e in v]
    offset = 0
    coords = Tuple{Int64, Int64}[]
    for s in sizes
        push!(coords, (offset+1, offset+s ))
        offset += s
    end
    views = [view(x, c[1]:c[2]) for c in coords]
    MSM{T}(x,sizes,views)
end

function MSM{T}(data::Vector{T}, sizes) where T
    # WARNING: this holds  a reference on data (rename?)
    offset = 0
    coords = Tuple{Int64, Int64}[]

    for s in sizes
        push!(coords, (offset+1, offset+s ))
        offset += s
    end
    views = [view(data, c[1]:c[2]) for c in coords]
    MSM{T}(data, sizes, views)

end

function MSM(data::Vector{T}, sizes) where T
    # WARNING: this holds  a reference on data (rename?)
    offset = 0
    coords = Tuple{Int64, Int64}[]

    for s in sizes
        push!(coords, (offset+1, offset+s ))
        offset += s
    end
    views = [view(data, c[1]:c[2]) for c in coords]
    MSM{T}(data, sizes, views)

end

function zeros_like(m::MSM{T}) where T
    data = m.data*0
    sizes = m.sizes
    return MSM{T}(data, sizes)
end

Base.getindex(m::MSM, I::Vararg{Int, 2}) = m.views[I[1]][I[2]]
Base.length(m::MSM) = length(m.views)
vecvec(m::MSM) = [copy(e) for e in m.views]

function reset!(m)
    m.data[:] .*= 0.0
end

# ---------------------- #
# Constant decision rule #
# ---------------------- #

mutable struct ConstantDecisionRule{nx} <: AbstractDecisionRule{EmptyGrid,EmptyGrid,nx}
    constants::SVector{nx,Float64}
end

function ConstantDecisionRule(constants::Vector{Float64})
    nx = length(constants)
    ConstantDecisionRule{nx}(SVector{nx,Float64}(constants...))
end

function ConstantDecisionRule(model::Model)
    ConstantDecisionRule(model.calibration[:controls])
end

(dr::ConstantDecisionRule)(x::Point) = dr.constants
(dr::ConstantDecisionRule)(x::Vector{Point{d}}) where d = [dr.constants for n=1:length(x)]
(dr::ConstantDecisionRule)(i::Int, x::Union{Point{d},Vector{Point{d}}}) where d = dr(x)
(dr::ConstantDecisionRule)(i::Int, j::Int, x::Union{Point{d},Vector{Point{d}}}) where d = dr(x)

# (dr::ConstantDecisionRule)(x::AbstractVector, y::AbstractVector) = dr.constants
# (dr::ConstantDecisionRule)(x::AbstractVector, y::AbstractMatrix) = repeat(dr.constants', size(y, 1), 1)
# (dr::ConstantDecisionRule)(x::AbstractMatrix, y::AbstractVector) = repeat(dr.constants', size(x, 1), 1)
# (dr::ConstantDecisionRule)(x::AbstractMatrix, y::AbstractMatrix) = repeat(dr.constants', size(x, 1), 1)
# (dr::ConstantDecisionRule)(i::Int, x::Union{AbstractVector,AbstractMatrix}) = dr(x)
# (dr::ConstantDecisionRule)(i::Int, j::Int, x::Union{AbstractVector,AbstractMatrix}) = dr(x)

# ------------------------------ #
# 2-dimensional Taylor Expansion #
# ------------------------------ #

mutable struct BiTaylorExpansion{nx} <: AbstractDecisionRule{EmptyGrid,EmptyGrid,nx}
    m0::Vector{Float64}
    s0::Vector{Float64}
    x0::Vector{Float64}
    x_m::Matrix{Float64}
    x_s::Matrix{Float64}
end

(dr::BiTaylorExpansion)(m::Point, s::Point) = SVector( (dr.x0 + dr.x_m*(m-dr.m0) + dr.x_s*(s-dr.s0))... )

(dr::BiTaylorExpansion)(m::AbstractVector{<:Point}, s::AbstractVector{<:Point}) = [dr(m[i], s[i]) for i=1:length(s)]

(dr::BiTaylorExpansion)(m::AbstractVector, s::AbstractVector) = Array(dr(SVector(m...), SVector(s...)))
(dr::BiTaylorExpansion)(m::AbstractMatrix, s::AbstractVector) = vcat([(dr(m[i, :], s))' for i=1:size(m, 1) ]...)
(dr::BiTaylorExpansion)(m::AbstractVector, s::AbstractMatrix) = vcat([(dr(m, s[i, :]))' for i=1:size(s, 1) ]...)
(dr::BiTaylorExpansion)(m::AbstractMatrix, s::AbstractMatrix) = vcat([(dr(m[i, :], s[i, :]))' for i=1:size(m, 1) ]...)

# User defined functions
struct CFunDR{S,T,nx} <: AbstractDecisionRule{S,T,nx}
    fun::Function
    dprocess::AbstractDiscretizedProcess
end
CFunDR(fun::Function, dprocess::AbstractDiscretizedProcess, n_x::Int) = CFunDR{typeof(dprocess.grid), EmptyGrid, n_x}(fun, dprocess)
(cfdr::CFunDR)(i::Int, x::Point{d}) where d = cfdr.fun(node(Point, cfdr.dprocess,i),x)
(cfdr::CFunDR)(i::Int, x::Vector{Point{d}}) where d = [cfdr(i,e) for e in x]


text:
module ConvergencePlots

export  ConvergencePlot,
        addpoint!,
        updateplot!,
        closeplot!

using PyPlot, Parameters, Random

mutable struct ConvergencePlot
    fig
    axs
    names::Vector{String}
    history::Dict
    options::Dict
    npoints::Int
end
function ConvergencePlot(
    npoints::Int = 100000;
    names = ["Residual"],
    options = nothing,
    show = false,
)
    history = Dict(n => Float64[] for n in names)
    options = getoptions(names, options)
    @assert keys(history) == keys(options)
    fig, axs = emptyplot(length(names))
    if show
        PyPlot.show(block=false)
    end
    return ConvergencePlot(fig, axs, names, history, options, npoints)
end

function addpoint!(plot::ConvergencePlot, y::Real; show = true)
    @assert length(keys(plot.history)) == 1
    k = collect(keys(plot.history))[1]
    return addpoint!(plot, Dict(k => y), show = show)
end
function addpoint!(plot::ConvergencePlot, y::Dict; show = true)
    @unpack history, npoints = plot
    for k in keys(y)
        if length(history[k]) < npoints
            push!(history[k], y[k])
        else
            history[k][1:npoints-1] .= history[k][2:npoints]
            history[k][npoints] = y[k]
        end
    end
    updateplot!(plot; show = show)
    return plot
end

function updateplot!(plot::ConvergencePlot; show = true)
    @unpack fig, axs, names, history, options = plot
    keepinteractive() do
        for ax in axs
            ax.clear()
        end
        for i in 1:length(names)
            ax = axs[i]
            n = names[i]
            ax.plot(1:length(history[n]), history[n]; options[n]...)
            if mod(i, 2) == 1
                ax.set_xlabel("Iteration")
            end
            ax.set_ylabel(n)
            ax.tick_params(axis="y", labelcolor=options[n].color)
        end
        fig.canvas.draw()
        fig.legend(loc=2)
        if show
            PyPlot.show(block=false)
        end
    end
    return plot
end

function closeplot!(plot::ConvergencePlot)
    PyPlot.close(plot.fig)
    fig, axs = emptyplot(length(plot.names))
    @pack! plot = fig, axs
    return plot
end

function emptyplot(ncurves)
    if ncurves % 2 == 0
        nsubplots = (ncurves ÷ 2)
    else
        nsubplots = (ncurves ÷ 2) + 1
    end
    local fig, axs
    keepinteractive() do
        fig, axs1 = PyPlot.subplots(nsubplots)
        axs = []
        for i in 1:nsubplots
            ax1 = nsubplots == 1 ? axs1 : axs1[i]
            push!(axs, ax1)
            if 2i <= ncurves
                ax2 = ax1.twinx()
                push!(axs, ax2)
            end
        end
    end
    return fig, axs
end

randcolor(rng) = (rand(rng), rand(), rand()) .* 0.7

function getoptions(names, old_options)
    rng = Random.MersenneTwister()
    Random.seed!(rng, 12344321)

    options = Dict()
    if old_options === nothing
        for n in names
            options[n] = (label = n,)
        end
    elseif old_options isa Dict
        for n in names
            if !(haskey(old_options[n], :label))
                options[n] = merge(old_options[n], (label = n,))
            end
        end
    else # Same options for all names
        for n in names
            options[n] = merge(old_options, (label = n,))
        end
    end

    # Fix the same color problem when using 2 y-axes
    for i in 1:(length(names) ÷ 2)
        n1, n2 = names[2i - 1], names[2i]
        if !(haskey(options[n2], :color))
            color = randcolor(rng)
            options[n2] = merge(options[n2], (color = color,))
        end
        if !(haskey(options[n1], :color))
            color = randcolor(rng)
            options[n1] = merge(options[n1], (color = color,))
        end
    end
    if mod(length(names), 2) == 1
        n = names[end]
        if !(haskey(options[n], :color))
            color = randcolor(rng)
            options[n] = merge(options[n], (color = color,))
        end
    end
    return options
end

function keepinteractive(f)
    was_interactive = PyPlot.isinteractive()
    was_interactive && PyPlot.ioff()
    f()
    was_interactive && PyPlot.ion()
end

end


text:
#!/usr/bin/env julia

using PencilArrays
using PencilArrays.MPITopologies

using MPI

using BenchmarkTools
using LinearAlgebra: transpose!
using Random
using Test

include("include/MPITools.jl")
using .MPITools

const BENCHMARK_ARRAYS = "--benchmark" in ARGS
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 1.0

Indexation(::Type{IndexLinear}) = LinearIndices
Indexation(::Type{IndexCartesian}) = CartesianIndices

# For testing `similar` for PencilArray.
struct DummyArray{T,N} <: AbstractArray{T,N}
    dims :: Dims{N}
    DummyArray{T}(::UndefInitializer, dims) where {T} = new{T,length(dims)}(dims)
end
Base.size(x::DummyArray) = x.dims
Base.getindex(::DummyArray{T}, ind...) where {T} = zero(T)
Base.similar(x::DummyArray, ::Type{S}, dims::Dims) where {S} = DummyArray{S}(undef, dims)

function benchmark_fill!(::Type{T}, u, val) where {T <: IndexStyle}
    indices = Indexation(T)(u)
    @inbounds for I in indices
        u[I] = val
    end
    u
end

function test_array_wrappers(p::Pencil, ::Type{T} = Float64) where {T}
    u = PencilArray{T}(undef, p)

    @test match(r"PencilArray{.*}\(::Pencil{.*}\)", summary(u)) !== nothing

    perm = permutation(u)
    @test perm === permutation(typeof(u))
    let topo = topology(p)
        @test topo === topology(u)
        I = coords_local(topo)
        @test range_remote(u, I) === range_remote(p, I)
    end

    @test eltype(u) === eltype(u.data) === T
    @test length.(axes(u)) === size(u)
    @test sizeof_global(u) == sizeof(T) * prod(size_global(u))
    @test sizeof_global((u, u)) == 2 * sizeof_global(u)
    let umat = [u for i = 1:2, j = 1:3]
        @test sizeof_global(umat) == 6 * sizeof_global(u)
    end

    let
        A = PermutedDimsArray(parent(u), perm)
        @test strides(A) === strides(u)
    end

    randn!(u)
    @test check_iteration_order(u)

    @inferred global_view(u)
    ug = global_view(u)
    @test check_iteration_order(ug)

    if BENCHMARK_ARRAYS
        for S in (IndexLinear, IndexCartesian)
            @info("Filling arrays using $S (Array, PencilArray, GlobalPencilArray)",
                  permutation(p))
            for v in (parent(u), u, ug)
                val = 3 * oneunit(eltype(v))
                @btime benchmark_fill!($S, $v, $val)
            end
            println()
        end
    end

    @testset "similar" begin
        let v = @inferred similar(u)
            @test typeof(v) === typeof(u)
            @test size(v) == size(u)
            @test pencil(v) === pencil(u)
        end

        let v = @inferred similar(u, Int)
            @test v isa PencilArray
            @test size(v) == size(u)
            @test eltype(v) === Int
            @test pencil(v) === pencil(u)
        end

        let v = @inferred similar(u, (3, 4))
            @test v isa Matrix
            @test size(v) == (3, 4)
            @test eltype(v) === eltype(u)
        end

        let v = @inferred similar(u, Int, (3, 4))
            @test v isa Matrix
            @test size(v) == (3, 4)
            @test eltype(v) === Int
        end

        let A = DummyArray{Int}(undef, size_local(p, MemoryOrder()))
            local u = @inferred PencilArray(p, A)
            @test parent(u) === A

            v = @inferred similar(u)
            @test typeof(v) === typeof(u)
            @test size(v) == size(u)

            w = @inferred similar(u, (4, 2))
            @test w isa DummyArray{Int,2}
            @test size(w) == (4, 2)

            z = @inferred similar(u, Float32, (3, 4, 2))
            @test z isa DummyArray{Float32,3}
            @test size(z) == (3, 4, 2)
        end
    end

    let v = similar(u)
        @test typeof(v) === typeof(u)

        psize = size_local(p, LogicalOrder())
        @test psize === size(v) === size(u)
        @test psize === size_local(u, LogicalOrder()) === size_local(v, LogicalOrder())

        vp = parent(v)
        randn!(vp)
        I = size(v) .>> 1  # non-permuted indices
        J = perm * I
        @test v[I...] == vp[J...]  # the parent takes permuted indices
    end

    let psize = size_local(p, MemoryOrder())
        a = zeros(T, psize)
        u = PencilArray(p, a)
        @test u.data === a
        @test IndexStyle(typeof(u)) === IndexStyle(typeof(a)) === IndexLinear()

        b = zeros(T, psize .+ 2)
        @test_throws DimensionMismatch PencilArray(p, b)
        @test_throws DimensionMismatch PencilArray(p, zeros(T, 3, psize...))

        # This is allowed.
        w = PencilArray(p, zeros(T, psize..., 3))
        @test size_global(w) === (size_global(p)..., 3)

        @inferred PencilArray(p, zeros(T, psize..., 3))
        @inferred size_global(w)
    end

    nothing
end

function test_multiarrays(pencils::Vararg{Pencil,M};
                          element_type::Type{T} = Float64) where {M,T}
    @assert M >= 3
    @inferred ManyPencilArray{T}(undef, pencils...)

    A = ManyPencilArray{T}(undef, pencils...)

    @test ndims(A) === ndims(first(pencils))
    @test eltype(A) === T
    @test length(A) === M

    @inferred first(A)
    @inferred last(A)
    @inferred A[Val(2)]
    @inferred A[Val(M)]

    @test_throws ErrorException @inferred A[2]  # type not inferred

    @test A[Val(1)] === first(A) === A[Val(UInt8(1))] === A[1]
    @test A[Val(2)] === A[2] === A.arrays[2] === A[Val(Int32(2))]
    @test A[Val(M)] === last(A)

    @test_throws BoundsError A[Val(0)]
    @test_throws BoundsError A[Val(M + 1)]

    @testset "In-place extra dimensions" begin
        e = (3, 2)
        @inferred ManyPencilArray{T}(undef, pencils...; extra_dims=e)
        A = ManyPencilArray{T}(undef, pencils...; extra_dims=e)
        @test extra_dims(first(A)) === extra_dims(last(A)) === e
        @test ndims_extra(first(A)) == ndims_extra(last(A)) == length(e)
    end

    @testset "In-place transpose" begin
        u = A[Val(1)]
        v = A[Val(2)]
        w = A[Val(3)]

        randn!(u)
        u_orig = copy(u)

        transpose!(v, u)  # this also modifies `u`!
        @test compare_distributed_arrays(u_orig, v)

        # In the 1D decomposition case, this is a local transpose, since v and w
        # only differ in the permutation.
        transpose!(w, v)
        @test compare_distributed_arrays(u_orig, w)
    end

    nothing
end

function check_iteration_order(u::Union{PencilArray,GlobalPencilArray})
    p = parent(parent(u)) :: Array  # two `parent` are needed for GlobalPencilArray
    cart = CartesianIndices(u)
    lin = LinearIndices(u)

    # Check that the behaviour of `cart` is consistent with that of
    # CartesianIndices.
    @assert size(CartesianIndices(p)) == size(p)
    @test size(cart) == size(u)

    # Same for `lin`.
    @assert size(LinearIndices(p)) == size(p)
    @test size(lin) == size(u)

    # Check that Cartesian indices iterate in memory order.
    for (n, I) in enumerate(cart)
        l = lin[I]
        @assert l == n
        u[n] == p[n] == u[I] == u[l] || return false
    end

    # Also test iteration on LinearIndices and their conversion to Cartesian
    # indices.
    for (n, l) in enumerate(lin)
        @assert l == n
        # Convert linear to Cartesian index.
        I = cart[l]  # this is relatively slow, don't do it in real code!
        u[n] == p[n] == u[I] == u[l] || return false
    end

    true
end

function compare_distributed_arrays(u_local::PencilArray, v_local::PencilArray)
    comm = get_comm(u_local)
    root = 0
    myrank = MPI.Comm_rank(comm)

    u = gather(u_local, root)
    v = gather(v_local, root)

    same = Ref(false)
    if u !== nothing && v !== nothing
        @assert myrank == root
        same[] = u == v
    end
    MPI.Bcast!(same, root, comm)

    same[]
end

function main()
    Nxyz = (16, 21, 41)
    comm = MPI.COMM_WORLD
    Nproc = MPI.Comm_size(comm)
    myrank = MPI.Comm_rank(comm)

    silence_stdout(comm)

    rng = MersenneTwister(42 + myrank)

    # Let MPI_Dims_create choose the values of (P1, P2).
    proc_dims = MPITopologies.dims_create(comm, Val(2))

    # Note that using dims_create is the default in MPITopology
    @test MPITopology(comm, proc_dims) == MPITopology(comm, Val(2))

    @test_throws ArgumentError MPITopology(comm, proc_dims .- 1)
    @test_throws ArgumentError MPITopology(comm, proc_dims .+ 1)
    topo = MPITopology(comm, proc_dims)
    @test match(
        r"^MPI topology: 2D decomposition \(\d+×\d+ processes\)$",
        string(topo),
    ) !== nothing
    @test ndims(topo) == length(proc_dims) == 2

    pen1 = @inferred Pencil(topo, Nxyz)
    let p = @inferred Pencil(topo, Nxyz, (2, 3))  # this is the default decomposition
        @test decomposition(p) === decomposition(pen1)
    end
    pen2 = @inferred Pencil(pen1, decomp_dims=(1, 3), permute=Permutation(2, 3, 1))
    pen3 = @inferred Pencil(pen2, decomp_dims=(1, 2), permute=Permutation(3, 2, 1))

    @test match(r"Pencil{3, 2, NoPermutation}", summary(pen1)) !== nothing
    @test match(r"Pencil{3, 2, Permutation{.*}}", summary(pen2)) !== nothing

    println("Pencil 1: ", pen1, "\n")
    println("Pencil 2: ", pen2, "\n")
    println("Pencil 3: ", pen3, "\n")

    @testset "Pencil constructors" begin
        p = @inferred Pencil((5, 4, 4, 3), MPI.COMM_WORLD)
        @test decomposition(p) == (2, 3, 4)
        @test ndims(topology(p)) == 3
        @test permutation(p) == NoPermutation()

        p = @inferred Pencil((5, 4, 4, 3), (2, 3), MPI.COMM_WORLD;
                             permute = Permutation(2, 3, 4, 1))
        @test decomposition(p) == (2, 3)
        @test ndims(topology(p)) == 2
        @test permutation(p) == Permutation(2, 3, 4, 1)
    end

    @testset "ManyPencilArray" begin
        test_multiarrays(pen1, pen2, pen3)
    end

    @testset "Topology" begin
        @test topology(pen2) === topo
        @test range_remote(pen2, coords_local(topo)) == range_local(pen2)
        @test eachindex(topo) isa LinearIndices
        for (n, I) in zip(eachindex(topo), CartesianIndices(topo))
            for order in (MemoryOrder(), LogicalOrder())
                @test range_remote(pen2, Tuple(I), order) ==
                    range_remote(pen2, n, order)
            end
        end
    end

    # Note: the permutation of pen2 was chosen such that the inverse permutation
    # is different.
    @assert permutation(pen2) != inv(permutation(pen2))

    @testset "Pencil constructor checks" begin
        # Too many decomposed directions
        @test_throws ArgumentError Pencil(
            MPITopology(comm, (Nproc, 1, 1)), Nxyz, (1, 2, 3))

        # Invalid permutations
        @test_throws TypeError Pencil(
            topo, Nxyz, (1, 2), permute=(2, 3, 1))
        @test_throws ArgumentError Pencil(
            topo, Nxyz, (1, 2), permute=Permutation(0, 3, 15))

        # Decomposed dimensions may not be repeated.
        @test_throws ArgumentError Pencil(topo, Nxyz, (2, 2))

        # Decomposed dimensions must be in 1:N = 1:3.
        @test_throws ArgumentError Pencil(topo, Nxyz, (1, 4))
        @test_throws ArgumentError Pencil(topo, Nxyz, (0, 2))

        @test Pencils.complete_dims(Val(5), (2, 3), (42, 12)) ===
            (1, 42, 12, 1, 1)
    end

    @testset "PencilArray" begin
        test_array_wrappers(pen1)
        test_array_wrappers(pen2)
        test_array_wrappers(pen3)
    end

    transpose_methods = (Transpositions.PointToPoint(),
                         Transpositions.Alltoallv())

    @testset "transpose! $method" for method in transpose_methods
        T = Float64
        u1 = PencilArray{T}(undef, pen1)
        u2 = PencilArray{T}(undef, pen2)
        u3 = PencilArray{T}(undef, pen3)

        # Set initial random data.
        randn!(rng, u1)
        u1 .+= 10 * myrank
        u1_orig = copy(u1)

        # Direct u1 -> u3 transposition is not possible!
        @test_throws ArgumentError transpose!(u3, u1, method=method)

        # Transpose back and forth between different pencil configurations
        transpose!(u2, u1, method=method)
        @test compare_distributed_arrays(u1, u2)

        transpose!(u3, u2, method=method)
        @test compare_distributed_arrays(u2, u3)

        transpose!(u2, u3, method=method)
        @test compare_distributed_arrays(u2, u3)

        transpose!(u1, u2, method=method)
        @test compare_distributed_arrays(u1, u2)

        @test u1_orig == u1

        # Test transpositions without permutations.
        let pen2 = Pencil(pen1, decomp_dims=(1, 3))
            u2 = PencilArray{T}(undef, pen2)
            transpose!(u2, u1, method=method)
            @test compare_distributed_arrays(u1, u2)
        end

    end

    # Test arrays with extra dimensions.
    @testset "extra dimensions" begin
        T = Float32
        u1 = PencilArray{T}(undef, pen1, 3, 4)
        u2 = PencilArray{T}(undef, pen2, 3, 4)
        u3 = PencilArray{T}(undef, pen3, 3, 4)
        @test range_local(u2) ===
            (range_local(pen2)..., Base.OneTo.((3, 4))...)
        @test range_remote(u2, 1) ===
            (range_remote(pen2, 1)..., Base.OneTo.((3, 4))...)
        randn!(rng, u1)
        transpose!(u2, u1)
        @test compare_distributed_arrays(u1, u2)
        transpose!(u3, u2)
        @test compare_distributed_arrays(u2, u3)

        for v in (u1, u2, u3)
            @test check_iteration_order(v)
        end

        @inferred global_view(u1)
    end

    # Test slab (1D) decomposition.
    @testset "1D decomposition" begin
        T = Float32
        topo = MPITopology(comm, (Nproc, ))
        @test ndims(topo) == 1

        pen1 = Pencil(topo, Nxyz, (1, ))
        pen2 = Pencil(pen1, decomp_dims=(2, ))

        # Same decomposed dimension as pen2, but different permutation.
        pen3 = Pencil(pen2, permute=Permutation(3, 2, 1))

        u1 = PencilArray{T}(undef, pen1)
        u2 = PencilArray{T}(undef, pen2)
        u3 = PencilArray{T}(undef, pen3)

        randn!(rng, u1)
        transpose!(u2, u1)
        @test compare_distributed_arrays(u1, u2)

        transpose!(u3, u2)
        @test compare_distributed_arrays(u1, u3)
        @test check_iteration_order(u3)

        # Test transposition between two identical configurations.
        transpose!(u2, u2)
        @test compare_distributed_arrays(u1, u2)

        let v = similar(u2)
            @test u2.pencil === v.pencil
            transpose!(v, u2)
            @test compare_distributed_arrays(u1, v)
        end

        test_multiarrays(pen1, pen2, pen3)
    end

    begin
        periods = zeros(Int, length(proc_dims))
        comm_cart = MPI.Cart_create(comm, collect(proc_dims), periods, false)
        @inferred MPITopologies.create_subcomms(Val(2), comm_cart)
        @test_throws ArgumentError MPITopology{3}(comm_cart)  # wrong dimensionality
        @inferred MPITopology{2}(comm_cart)
        @inferred MPITopologies.get_cart_ranks_subcomm(pen1.topology.subcomms[1])

        @inferred PencilArrays.to_local(pen2, (1:2, 1:2, 1:2), MemoryOrder())
        @inferred PencilArrays.to_local(pen2, (1:2, 1:2, 1:2), LogicalOrder())

        @inferred PencilArrays.size_local(pen2, MemoryOrder())

        T = Int
        @inferred PencilArray{T}(undef, pen2)
        @inferred PencilArray{T}(undef, pen2, 3, 4)

        u1 = PencilArray{T}(undef, pen1)
        u2 = PencilArray{T}(undef, pen2)

        @inferred Nothing gather(u2)
        @inferred transpose!(u2, u1)
        @inferred Transpositions.get_remote_indices(1, (2, 3), 8)
    end

    nothing
end

MPI.Init()
main()