text:
module ConversionWrapper

using JavaCall

const AbsoluteDate = @jimport org.orekit.time.AbsoluteDate
const AbstractIntegrator = @jimport org.hipparchus.ode.AbstractIntegrator
const AbstractPropagatorBuilder = @jimport org.orekit.propagation.conversion.AbstractPropagatorBuilder
const AbstractPropagatorConverter = @jimport org.orekit.propagation.conversion.AbstractPropagatorConverter
const AdamsBashforthIntegratorBuilder = @jimport org.orekit.propagation.conversion.AdamsBashforthIntegratorBuilder
const AdamsMoultonIntegratorBuilder = @jimport org.orekit.propagation.conversion.AdamsMoultonIntegratorBuilder
const AdditionalEquations = @jimport org.orekit.propagation.integration.AdditionalEquations
const AttitudeProvider = @jimport org.orekit.attitudes.AttitudeProvider
const BatchLSModel = @jimport org.orekit.estimation.leastsquares.BatchLSModel
const BatchLSODModel = @jimport org.orekit.estimation.leastsquares.BatchLSODModel
const Class = @jimport java.lang.Class
const ClassicalRungeKuttaIntegratorBuilder = @jimport org.orekit.propagation.conversion.ClassicalRungeKuttaIntegratorBuilder
const DSSTBatchLSModel = @jimport org.orekit.estimation.leastsquares.DSSTBatchLSModel
const DSSTForceModel = @jimport org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel
const DSSTKalmanModel = @jimport org.orekit.estimation.sequential.DSSTKalmanModel
const DSSTPropagator = @jimport org.orekit.propagation.semianalytical.dsst.DSSTPropagator
const DSSTPropagatorBuilder = @jimport org.orekit.propagation.conversion.DSSTPropagatorBuilder
const DataContext = @jimport org.orekit.data.DataContext
const DormandPrince54IntegratorBuilder = @jimport org.orekit.propagation.conversion.DormandPrince54IntegratorBuilder
const DormandPrince853IntegratorBuilder = @jimport org.orekit.propagation.conversion.DormandPrince853IntegratorBuilder
const EcksteinHechlerPropagatorBuilder = @jimport org.orekit.propagation.conversion.EcksteinHechlerPropagatorBuilder
const EulerIntegratorBuilder = @jimport org.orekit.propagation.conversion.EulerIntegratorBuilder
const FiniteDifferencePropagatorConverter = @jimport org.orekit.propagation.conversion.FiniteDifferencePropagatorConverter
const ForceModel = @jimport org.orekit.forces.ForceModel
const Frame = @jimport org.orekit.frames.Frame
const GillIntegratorBuilder = @jimport org.orekit.propagation.conversion.GillIntegratorBuilder
const GraggBulirschStoerIntegratorBuilder = @jimport org.orekit.propagation.conversion.GraggBulirschStoerIntegratorBuilder
const HighamHall54IntegratorBuilder = @jimport org.orekit.propagation.conversion.HighamHall54IntegratorBuilder
const IntegratedPropagatorBuilder = @jimport org.orekit.propagation.conversion.IntegratedPropagatorBuilder
const JString = @jimport java.lang.String
const JacobianPropagatorConverter = @jimport org.orekit.propagation.conversion.JacobianPropagatorConverter
const KalmanModel = @jimport org.orekit.estimation.sequential.KalmanModel
const KalmanODModel = @jimport org.orekit.estimation.sequential.KalmanODModel
const KeplerianPropagatorBuilder = @jimport org.orekit.propagation.conversion.KeplerianPropagatorBuilder
const List = @jimport java.util.List
const LutherIntegratorBuilder = @jimport org.orekit.propagation.conversion.LutherIntegratorBuilder
const MidpointIntegratorBuilder = @jimport org.orekit.propagation.conversion.MidpointIntegratorBuilder
const ModelObserver = @jimport org.orekit.estimation.leastsquares.ModelObserver
const NumericalPropagator = @jimport org.orekit.propagation.numerical.NumericalPropagator
const NumericalPropagatorBuilder = @jimport org.orekit.propagation.conversion.NumericalPropagatorBuilder
const ODEIntegratorBuilder = @jimport org.orekit.propagation.conversion.ODEIntegratorBuilder
const Object = @jimport java.lang.Object
const Orbit = @jimport org.orekit.orbits.Orbit
const OrbitType = @jimport org.orekit.orbits.OrbitType
const OsculatingToMeanElementsConverter = @jimport org.orekit.propagation.conversion.OsculatingToMeanElementsConverter
const ParameterDriversList = @jimport org.orekit.utils.ParameterDriversList
const PositionAngle = @jimport org.orekit.orbits.PositionAngle
const PropagationType = @jimport org.orekit.propagation.PropagationType
const Propagator = @jimport org.orekit.propagation.Propagator
const PropagatorBuilder = @jimport org.orekit.propagation.conversion.PropagatorBuilder
const PropagatorConverter = @jimport org.orekit.propagation.conversion.PropagatorConverter
const SpacecraftState = @jimport org.orekit.propagation.SpacecraftState
const TLE = @jimport org.orekit.propagation.analytical.tle.TLE
const TLEPropagatorBuilder = @jimport org.orekit.propagation.conversion.TLEPropagatorBuilder
const ThreeEighthesIntegratorBuilder = @jimport org.orekit.propagation.conversion.ThreeEighthesIntegratorBuilder
const TideSystem = @jimport org.orekit.forces.gravity.potential.TideSystem
const UnnormalizedSphericalHarmonicsProvider = @jimport org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider

export AbstractPropagatorBuilder
export AbstractPropagatorConverter
export AdamsBashforthIntegratorBuilder
export AdamsMoultonIntegratorBuilder
export ClassicalRungeKuttaIntegratorBuilder
export DSSTPropagatorBuilder
export DormandPrince54IntegratorBuilder
export DormandPrince853IntegratorBuilder
export EcksteinHechlerPropagatorBuilder
export EulerIntegratorBuilder
export FiniteDifferencePropagatorConverter
export GillIntegratorBuilder
export GraggBulirschStoerIntegratorBuilder
export HighamHall54IntegratorBuilder
export IntegratedPropagatorBuilder
export JacobianPropagatorConverter
export KeplerianPropagatorBuilder
export LutherIntegratorBuilder
export MidpointIntegratorBuilder
export NumericalPropagatorBuilder
export ODEIntegratorBuilder
export OsculatingToMeanElementsConverter
export PropagatorBuilder
export PropagatorConverter
export TLEPropagatorBuilder
export ThreeEighthesIntegratorBuilder
export add_additional_equations
export add_force_model
export build_integrator
export build_kalman_model
export build_ls_model
export build_propagator
export convert
export copy
export equals
export get_adapted_propagator
export get_all_force_models
export get_attitude_provider
export get_class
export get_evaluations
export get_frame
export get_initial_orbit_date
export get_integrator_builder
export get_mass
export get_mu
export get_orbit_type
export get_orbital_parameters_drivers
export get_position_angle
export get_position_scale
export get_propagation_parameters_drivers
export get_rms
export get_selected_normalized_parameters
export hash_code
export notify
export notify_all
export reset_orbit
export set_attitude_provider
export set_mass
export to_string
export wait

include("abstract_propagator_builder.jl")
include("abstract_propagator_converter.jl")
include("adams_bashforth_integrator_builder.jl")
include("adams_moulton_integrator_builder.jl")
include("classical_runge_kutta_integrator_builder.jl")
include("dormand_prince54_integrator_builder.jl")
include("dormand_prince853_integrator_builder.jl")
include("dsst_propagator_builder.jl")
include("eckstein_hechler_propagator_builder.jl")
include("euler_integrator_builder.jl")
include("finite_difference_propagator_converter.jl")
include("gill_integrator_builder.jl")
include("gragg_bulirsch_stoer_integrator_builder.jl")
include("higham_hall54_integrator_builder.jl")
include("integrated_propagator_builder.jl")
include("jacobian_propagator_converter.jl")
include("keplerian_propagator_builder.jl")
include("luther_integrator_builder.jl")
include("midpoint_integrator_builder.jl")
include("numerical_propagator_builder.jl")
include("ode_integrator_builder.jl")
include("osculating_to_mean_elements_converter.jl")
include("propagator_builder.jl")
include("propagator_converter.jl")
include("three_eighthes_integrator_builder.jl")
include("tle_propagator_builder.jl")

end


text:
using RandomNumbers.Xorshifts

@testset "PMDPg.jl" begin

    # initialize
    pp = simple_pp()
    rng = Xorshift128Plus(1)
    mg = PMDPs.PMDPg(pp)

    @test PMDPs.empty_product(mg) == PMDPs.Product(falses(PMDPs.n_resources(mg)), PMDPs.selling_period_end(mg))

    # Test methods
    @test 1<=PMDPs.sample_request(mg, 1, rng) <= PMDPs.n_products(mg)+1

    s = PMDPs.State(pp.c₀, 1, 1)
    @test typeof(PMDPs.sample_customer_budget(mg, s, rng)) == Float64
    @test minimum(pp.B[1])<=PMDPs.sample_customer_budget(mg, s, rng)<=maximum(pp.B[1])

    sₑ = PMDPs.State(pp.c₀, 1, mg.empty_product_id)
    @test PMDPs.sample_customer_budget(mg, sₑ, rng)==PMDPs.EMPTY_PRODUCT_USER_BUDGET

    s₀ = PMDPs.State(pp.c₀, 1, 1)
    a = 10.
    rng = Xorshift128Plus(1)
    sp, r, info = POMDPs.gen(mg, s, a, rng)
    @test sp.t==2
    @test r>=0.
    @test info.b>=0.

end

ext: jl


text:
module Progress
using Logging

using Atom: msg

const progs = Set()

"""
    JunoProgressLogger(previous_logger<:AbstractLogger)

Logger that only handles messages with `progress` argument set and passes all
others throught to the previously defined logger.

If the `progress` argument is set then a corresponding progress bar will be shown
in the Juno frontend with it's name set to the logging message.

Possible keyword arguments to logging messages consumed by this logger are:
  - `progress`:
    - `0 <= progress < 1`: create or update progress bar
    - `progress == nothing || progress = NaN`: set progress bar to indeterminate progress
    - `progress > 1 || progress == "done"`: destroy progress bar
  - `right_text`: Shown instead of a approximation of remaining time.
  - `message`: Shown in a tooltip.
  - `_id`: The progress bars' ID. Should be set (to a symbol or string) in most cases.
"""
struct JunoProgressLogger <: AbstractLogger
  previous_logger
end

function Logging.handle_message(j::JunoProgressLogger, level, message, _module,
                                group, id, file, line; kwargs...)
  if haskey(kwargs, :progress)
    if !(id in progs)
      push!(progs, id)
      msg("progress", "add", id)
    end
    msg("progress", "leftText", id, message)

    prog = kwargs[:progress]
    if prog isa Real
      if 0 <= prog < 1
        msg("progress", "progress", id, prog)
      elseif prog >= 1
        msg("progress", "delete", id)
        delete!(progs, id)
      end
    elseif prog isa Nothing
      msg("progress", "progress", id, NaN)
    elseif prog isa AbstractString && prog == "done"
      msg("progress", "delete", id)
      delete!(progs, id)
    else
      msg("progress", "progress", id, NaN)
    end

    if haskey(kwargs, :right_text)
      msg("progress", "rightText", id, kwargs[:right_text])
    end

    if haskey(kwargs, :message)
      msg("progress", "message", id, kwargs[:message])
    end
  else
    if (Logging.min_enabled_level(j.previous_logger) <= Logging.LogLevel(level) || Base.CoreLogging.env_override_minlevel(group, _module)) &&
        Logging.shouldlog(j.previous_logger, level, _module, group, id)
      Logging.handle_message(j.previous_logger, level, message, _module,
                             group, id, file, line; kwargs...)
    end
  end
end

Logging.shouldlog(::JunoProgressLogger, level, _module, group, id) = true

Logging.catch_exceptions(::JunoProgressLogger) = true

Logging.min_enabled_level(j::JunoProgressLogger) = min(Logging.min_enabled_level(j.previous_logger), Logging.LogLevel(-1))
end


text:
module TestExamplesEmptyResultHandling
using LiterateTest.AssertAsTest: @assert
include("../examples/empty_result_handling.jl")
end  # module


text:
# TODO, NOTE
# We should be doing each f1xfm **AFTER** each trial as been cycle averaged

export f1xfm, cycle_mean
# ============================================================================ #
@inline calcf1(cs::Vector{C}, d::Matrix{T}) where {C<:Complex,T<:Real} =
    vec(cs' * d) .* (2.0 / length(cs))
# ---------------------------------------------------------------------------- #
@inline calcf1(cs::Vector{C}, d::Vector{T}) where {C<:Complex,T<:Real} =
    dot(cs, d) .* (2.0 / length(cs))
# ---------------------------------------------------------------------------- #
@inline f1basis(npt::Int) = exp.(-im*2.0*pi*range(0.0, stop=1.0, length=npt))
# ============================================================================ #
function f1xfm(d::Array{T}, f1::AbstractFloat, dur::AbstractFloat) where {T<:Real}
    bpc = floor(Int, size(d, 1) / (dur * f1))
    tmp = cycle_mean(d, bpc)
    return calcf1(f1basis(bpc), tmp)
end
# ============================================================================ #
function f1xfm(d::Vector{Matrix{T}}, f1::AbstractFloat,
    dur::AbstractFloat) where {T<:Real}

    siz = size(d[1], 1)
    check_size(x::Matrix{T}) = size(x, 1) == siz

    out = [Vector{Float64}(undef, size(d[k], 2)) for k in 1:length(d)]

    if !all(check_size, d[2:end])
        # psth matricies have different number of bins so the f1basis must be
        # caclulated for each matrix
        @inbounds for k = 1:length(d)
            out[k] = f1xfm(d[k], f1, dur)
        end
    else
        # our sinusodial basis can be reused so pre-allocate
        bpc = floor(Int, size(d[1], 1) / (dur * f1))
        cs = f1basis(bpc)

        @inbounds for k = 1:length(d)
            out[k] = calcf1(cs, cycle_mean(d[k], bpc))
        end
    end

    return out
end
# ============================================================================ #
function f1xfm(d::Vector{Matrix{T}}, f1::Vector{F},
    dur::AbstractFloat) where {T<:Real, F<:AbstractFloat}

    length(f1) != length(d) && error("Number of temporal frequencies and number of trial groups *MUST* match")

    out = Vector{Vector{Float64}}(undef, length(d))

    @inbounds for k = 1:length(d)
        out[k] = f1xfm(d[k], f1[k], dur)
    end

    return out
end
# ============================================================================ #
function cycle_pad(npt::Int, bpc::Integer)
    ncycle = floor(Int, npt / bpc)
    if npt % bpc > 0
        ncycle += 1
        npad = (bpc * ncycle) - npt
    else
        npad = 0
    end
    return npad, ncycle
end
# ---------------------------------------------------------------------------- #
function cycle_mean(d::Vector{T}, bpc::Integer) where {T<:Real}
    pad, ncycle = cycle_pad(length(d), bpc)
    den = [fill(ncycle, bpc-pad); fill(ncycle-1, pad)]
    return vec(
                sum(
                    reshape(cat(d, zeros(T, pad), dims=1), bpc, ncycle),
                    dims=2
                )
            ) ./ den
end
# ---------------------------------------------------------------------------- #
function cycle_mean(d::Matrix{T}, bpc::Integer) where {T<:Real}
    pad, ncycle = cycle_pad(size(d, 1), bpc)
    ntrial = size(d, 2)

    den = [fill(ncycle, bpc-pad); fill(ncycle-1, pad)]
    tmp = cat(d, zeros(pad, ntrial), dims=1)

    out = Matrix{T}(undef, bpc, ntrial)
    @inbounds @fastmath for k in 1:ntrial
        out[:,k] = sum(reshape(tmp[:,k], bpc, ncycle), dims=2) ./ den
    end

    return out
end
# ============================================================================ #


text:
# This file is for functions that explictly mess with the dimensions of a NameDimsArray

"""
    rename(nda::NamedDimsArray, names)

Returns a new `NameDimsArray` with the given dimension `names`.
`rename` outright replaces the names; while still wrapping the same backing array.
Unlike the constructor, it does not require that new names are compatible
with the old names (though you do still need to match the number of dimensions).
"""
rename(nda::NamedDimsArray, names) = NamedDimsArray(parent(nda), names)

function Base.dropdims(nda::NamedDimsArray; dims)
    numerical_dims = dim(nda, dims)
    data = dropdims(parent(nda); dims=numerical_dims)
    L = remaining_dimnames_after_dropping(names(nda), numerical_dims)
    return NamedDimsArray{L}(data)
end

function Base.permutedims(nda::NamedDimsArray{L}, perm) where {L}
    numerical_perm = dim(nda, perm)
    new_names = permute_dimnames(L, numerical_perm)

    return NamedDimsArray{new_names}(permutedims(parent(nda), numerical_perm))
end

# We use CoVector to workout if we are taking the tranpose of a tranpose etc
const CoVector = Union{Adjoint{<:Any, <:AbstractVector}, Transpose{<:Any, <:AbstractVector}}
for f in (
    :(Base.transpose),
    :(Base.adjoint),
    :(Base.permutedims),
    :(LinearAlgebra.pinv)
)
    # Vector
    @eval function $f(nda::NamedDimsArray{L,T,1}) where {L,T}
        new_names = (:_, first(L))
        return NamedDimsArray{new_names}($f(parent(nda)))
    end

    # Vector Double Transpose
    if f !== :permutedims
        @eval function $f(nda::NamedDimsArray{L,T,2,A}) where {L,T,A<:CoVector}
            new_names = (last(L),)  # drop the name of the first dimensions
            return NamedDimsArray{new_names}($f(parent(nda)))
        end
    end

    # Matrix
    @eval function $f(nda::NamedDimsArray{L,T,2}) where {L,T}
        new_names = (last(L), first(L))
        return NamedDimsArray{new_names}($f(parent(nda)))
    end
end

# reshape
# For now we only implement the version that drops dimension names
Base.reshape(nda::NamedDimsArray, dims::Int...) = reshape(parent(nda), dims...)


text:
#
function find_exe_paths(prog)
    # 8f1kRCu
    env_pathext = get(ENV, "PATHEXT", "")

    # 4fpQ2RB
    if env_pathext == ""
        # 9dqlPRg
        return []
    end

    # 6qhHTHF
    # Split into a list of extensions
    ext_s = split(env_pathext, ";")

    # 2pGJrMW
    # Strip
    ext_s = [strip(ext) for ext=ext_s]

    # 2gqeHHl
    # Remove empty.
    # Must be done after the stripping at 2pGJrMW.
    ext_s = filter(x -> x != "", ext_s)

    # 2zdGM8W
    # Convert to lowercase
    ext_s = [lowercase(ext) for ext=ext_s]

    # 2fT8aRB
    # Uniquify
    ext_s = unique(ext_s)

    # 4ysaQVN
    env_path = get(ENV, "PATH", "")

    # 5gGwKZL
    if env_path == ""
        # 7bVmOKe
        # Go ahead with "dir_path_s" being empty
        dir_path_s = []
    else
        # 6mPI0lg
        # Split into a list of dir paths
        dir_path_s = split(env_path, ";")
    end

    # 5rT49zI
    # Insert empty dir path to the beginning.
    #
    # Empty dir handles the case that "prog" is a path, either relative or
    # absolute. See code 7rO7NIN.
    unshift!(dir_path_s, "")

    # 2klTv20
    # Uniquify
    dir_path_s = unique(dir_path_s)

    # 9gTU1rI
    # Check if "prog" ends with one of the file extension in "ext_s".
    #
    # "ext_s" are all in lowercase, ensured at 2zdGM8W.
    prog_lc = lowercase(prog)

    prog_has_ext = (findfirst(ext -> endswith(prog_lc, ext), ext_s)) != 0
    # Index value "0" means not found.

    # 6bFwhbv
    exe_path_s = String[]

    for dir_path=dir_path_s
        # 7rO7NIN
        # Synthesize a path
        if dir_path == ""
            path = prog
        else
            path = dir_path * "\\" *  prog
        end

        # 6kZa5cq
        # If "prog" ends with executable file extension
        if prog_has_ext
            # 3whKebE
            if isfile(path)
                # 2ffmxRF
                push!(exe_path_s, path)
            end
        end

        # 2sJhhEV
        # Assume user has omitted the file extension
        for ext=ext_s
            # 6k9X6GP
            # Synthesize a path with one of the file extensions in PATHEXT
            path_2 = path * ext

            # 6kabzQg
            if isfile(path_2)
                # 7dui4cD
                push!(exe_path_s, path_2)
            end
        end
    end

    # 8swW6Av
    # Uniquify
    exe_path_s = unique(exe_path_s)

    # 7y3JlnS
    return exe_path_s
end

#
function main()
    # 9mlJlKg
    # If not exactly one command argument is given
    if length(ARGS) != 1
        usage = """Usage: aoikwinwhich PROG

#/ PROG can be either name or path
aoikwinwhich notepad.exe
aoikwinwhich C:\\Windows\\notepad.exe

#/ PROG can be either absolute or relative
aoikwinwhich C:\\Windows\\notepad.exe
aoikwinwhich Windows\\notepad.exe

#/ PROG can be either with or without extension
aoikwinwhich notepad.exe
aoikwinwhich notepad
aoikwinwhich C:\\Windows\\notepad.exe
aoikwinwhich C:\\Windows\\notepad"""

        # 7rOUXFo
        # Print program usage
        print(usage)

        # 3nqHnP7
        # Exit
        return 1
    end

    # 9m5B08H
    # Get executable name or path
    prog = ARGS[1]

    # 8ulvPXM
    # Find executable paths
    exe_path_s = find_exe_paths(prog)

    # 5fWrcaF
    # If has found none
    if length(exe_path_s) == 0
        # 3uswpx0
        # Exit
        return 2
    # If has found some
    else
        # 9xPCWuS
        # Print result
        print(join(exe_path_s, "\n") * "\n")

        # 4s1yY1b
        # Exit
        return 0
    end
end

# 4zKrqsC
# Program entry
exit(main())

max_stars_repo_name: AoiKuiyuyou/AoikWinWhich-Julia

max_issues_repo_name: AoiKuiyuyou/AoikWinWhich-Julia

max_forks_repo_name: AoiKuiyuyou/AoikWinWhich-Julia


text:
using BinDeps

@BinDeps.setup

world = library_dependency("libworld", aliases=["libworld", "world-0"])

const version = "0.3.1"

# TODO
if Sys.iswindows() && Sys.WORD_SIZE == 32
    error("Your platform isn't supported yet.")
end

github_root = "https://github.com/r9y9/World-cmake"
arch = Sys.WORD_SIZE == 64 ? "x86_64" : "i686"
major = version[1]
provides(Binaries,
         URI("$(github_root)/releases/download/v$(version)/WORLD-$(major)_mingw$(Sys.WORD_SIZE)_$(arch).zip"),
         world, unpacked_dir = "usr/lib", os = :Windows)

provides(Sources,
         URI("$(github_root)/archive/v$(version).tar.gz"),
         world,
         unpacked_dir="World-cmake-$(version)")

prefix = joinpath(BinDeps.depsdir(world), "usr")
srcdir = joinpath(BinDeps.depsdir(world), "src", "World-cmake-$(version)")
builddir = joinpath(srcdir, "build")

provides(SimpleBuild,
          (@build_steps begin
              GetSources(world)
              @build_steps begin
                  ChangeDirectory(builddir)
                  `cmake -DCMAKE_BUILD_TYPE=Release -DBUILD_SHARED_LIBS=ON -DCMAKE_INSTALL_PREFIX=$prefix ..`
                  `make`
                  `make install`
              end
           end), world, os = :Unix)

@BinDeps.install Dict(:libworld => :libworld)


text:
using AbstractTrees
using AbstractTrees: repr_tree
using Test
import Base: ==

if VERSION >= v"1.1.0-DEV.838" # requires https://github.com/JuliaLang/julia/pull/30291
    @testset "Ambiguities" begin
        @test isempty(detect_ambiguities(AbstractTrees, Base, Core))
    end
end

@testset "Builtins" begin include("builtins.jl") end
@testset "Custom tree types" begin include("trees.jl") end
@testset "Printing" begin include("printing.jl") end
@testset "Examples" begin include("examples.jl") end


text:
mutable struct pointcloud2_t <: LCMType
    utime::Int64
    seq::Int32
    frame_id::String
    height::Int32
    width::Int32

    nfields::Int32
    fields::Vector{pointfield_t}

    is_bigendian::Bool
    point_step::Int32
    row_step::Int32
    data_nbytes::Int32
    data::Vector{UInt8}
    is_dense::Bool
end

@lcmtypesetup(pointcloud2_t,
    fields => (nfields,),
    data => (data_nbytes,),
)


text:
# general contracted DZ basis [3s1p/2s1p] for H₂
#     exponents    contract-coeff
# S   6.0          0.7               0.4
#     2.0          0.6               0.3
#     0.8          0.5               0.2
# P   0.9          1.

function gto_norm(n, a)
   # normalization factor of function rⁿ exp(-ar²)
    s = 2^(2n+3) * factorial(n+1) * (2a)^(n+1.5) / (factorial(2n+2) * √π)
    sqrt(s)
end

CHARGE_OF  = 1
PTR_COORD  = 2
NUC_MOD_OF = 3
PTR_ZETA   = 4
ATM_SLOTS  = 6

ATOM_OF    = 1
ANG_OF     = 2
NPRIM_OF   = 3
NCTR_OF    = 4
KAPPA_OF   = 5
PTR_EXP    = 6
PTR_COEFF  = 7
BAS_SLOTS  = 8

ptr_env = 20
env = fill(0.0, ptr_env)

i = 0
#                     CHARGE_OF,PTR_COORD
atm =           Int32[1         ptr_env  0 0 0 0]
ptr_env += 3
atm = vcat(atm, Int32[1         ptr_env  0 0 0 0])
ptr_env += 3
#             x  y   z (Bohr)
append!(env, [0, 0, -0.8])
append!(env, [0, 0,  0.8])

# basis for atom #0
# 3s -> 2s
append!(env, [6., 2., .8])
append!(env, [.7 * gto_norm(0, 6.),
              .6 * gto_norm(0, 2.),
              .5 * gto_norm(0, .8),
              .4 * gto_norm(0, 6.),
              .3 * gto_norm(0, 2.),
              .2 * gto_norm(0, .8)])
#                     ATOM_OF, ANG_OF, NPRIM_OF, NCTR_OF, KAPPA_OF, PTR_EXP, PTR_COEFF
bas =           Int32[0        0       3         2         0        ptr_env  ptr_env+3  0]
ptr_env += 9
bas = vcat(bas, Int32[0        1       1         1        0         ptr_env  ptr_env+1  0])
append!(env, [.9])
append!(env, [1. * gto_norm(1, 0.9)])
ptr_env += 2

# basis functions for atom #1, they are the same to thoes of atom #0
bas = vcat(bas, bas)

shls = Int32[0, 1]

natm = length(atm)
nbas = length(bas)

@info "Input" atm bas #env

bas = permutedims(bas, (2,1))   # Column (Julia) -> Row (C) order
atm = permutedims(atm, (2,1))   # Column (Julia) -> Row (C) order

# Libcint API
const LIBCINT = joinpath(@__DIR__, "path/to/libcint.so/dylib/dll")

function CINTcgtos_spheric(bas_id, bas)
   @static if VERSION >= v"1.5.0" # clearer syntax
      @ccall LIBCINT.CINTcgtos_spheric(bas_id::Cint, bas::Ptr{Cint})::Cint
   else # works anyway
      ccall( (:CINTcgtos_spheric, LIBCINT), Cint, (Cint, Ptr{Cint},), bas_id, bas)
   end
end

function cint1e_ipnuc_sph(buf, shls, atm, natm, bas, nbas, env)
   @static if VERSION >= v"1.5.0" # clearer syntax
      @ccall LIBCINT.cint1e_ipnuc_sph(
                                      buf  :: Ptr{Cdouble},
                                      shls :: Ptr{Cint},
                                      atm  :: Ptr{Cint},
                                      natm :: Cint,
                                      bas  :: Ptr{Cint},
                                      nbas :: Cint,
                                      env  :: Ptr{Cdouble}
                                     )::Cvoid
   else # works anyway
      ccall( (:cint1e_ipnuc_sph, LIBCINT), Cvoid, (
                                      Ptr{Cdouble},
                                      Ptr{Cint},
                                      Ptr{Cint},
                                      Cint,
                                      Ptr{Cint},
                                      Cint,
                                      Ptr{Cdouble}, ),
                                      buf  ,
                                      shls ,
                                      atm  ,
                                      natm ,
                                      bas  ,
                                      nbas ,
                                      env  ,
                                     )
   end
end

di = CINTcgtos_spheric(0, bas) # C uses 0-indexing!
dj = CINTcgtos_spheric(1, bas) # C uses 0-indexing!
@assert di == (2bas[ANG_OF,1]+1) * bas[NCTR_OF,1]
@assert dj == (2bas[ANG_OF,2]+1) * bas[NCTR_OF,2]

buf = zeros(3, dj, di) # Column (Julia) -> Row (C) order

cint1e_ipnuc_sph(buf, shls, atm, natm, bas, nbas, env)
@info "∇⟨i|nuc|j⟩ = " permutedims(buf, (3,2,1))

max_stars_repo_name: gustavojra/libcint

max_issues_repo_name: gustavojra/libcint

max_forks_repo_name: gustavojra/libcint


text:
module Test_afdisyn

using FaultDetectionTools
using DescriptorSystems
using LinearAlgebra
using Polynomials
using Random
using Test

# Testing examples for AFDSYN
println("Test_afdisyn")
@testset "afdisyn" begin
rand()

## Examples without control and disturbance inputs
p = 3; mf = 2;
sysf = fdimodset(rss(1,p,mf),faults = 1:mf);

# solve an AFDP using the nullspace based approach
@time Q, Rf, info = afdisyn(sysf; separate = false, atol = 1.e-7, minimal = false, rdim = [2]); info
R = fdIFeval(Q, sysf; minimal = true, atol = 1.e-7)
@test iszero(Rf.sys[1][:,Rf.faults]-R.sys[1][:,R.faults],atol=1.e-7) &&
      info.HDesign2[1] == [0.0 1.0 0.0; 0.0 0.0 1.0] &&
      all(isstable.(Q.sys)) && info.gap == [Inf] == fdif2ngap(R,fditspec(R))[1]

@time Q, Rf, info = afdisyn(sysf; atol3 = 1.e-7, minimal = false, rdim = [2], separate = true); info
R = fdIFeval(Q, sysf; minimal = true, atol = 1.e-7)
@test iszero(Rf.sys[1][:,Rf.faults]-R.sys[1][:,R.faults],atol=1.e-7) &&
      info.HDesign2[1] == [0.0 1.0 0.0; 0.0 0.0 1.0] && all(isstable.(Q.sys)) &&
      all(isstable.(Q.sys)) && info.gap == [Inf] == fdif2ngap(R,fditspec(R))[1]

p = 3; mf = 2; mw = 3;
sysf = fdimodset(rss(1,p,mf+mw),faults = 1:mf, noise = mf .+ (1:mw));

# solve an AFDP using the nullspace based approach
@time Q, Rf, info = afdisyn(sysf; separate = false, atol3 = 1.e-7, minimal = false, rdim = [2]); info
R = fdIFeval(Q, sysf; minimal = true, atol = 1.e-7)
@test iszero(Rf.sys[1][:,Rf.faults]-R.sys[1][:,R.faults],atol=1.e-7) && info.HDesign[1] == [0.0 1.0 0.0; 0.0 0.0 1.0] &&
all(isstable.(Q.sys)) && info.gap ≈ fdif2ngap(R,fditspec(R))[1]

@time Q, Rf, info = afdisyn(sysf; atol3 = 1.e-7, minimal = false, rdim = [2], separate = true); info
R = fdIFeval(Q, sysf; minimal = true, atol = 1.e-7)
@test iszero(Rf.sys[1][:,Rf.faults]-R.sys[1][:,R.faults],atol=1.e-7) && info.HDesign[1] == [0.0 1.0 0.0; 0.0 0.0 1.0]

## Example without disturbance inputs
p = 3; mf = 2; mu = 2
sysf = fdimodset(rss(1,p,mf+mu),controls = 1:mu, faults = mu+1:mu+mf);

# solve an EFDP using the nullspace based approach
@time Q, Rf, info = afdisyn(sysf; FDtol = 0.00001, atol = 1.e-7, minimal = false, rdim = [2]); info
R = fdIFeval(Q, sysf; minimal = true, atol = 1.e-7);
@test iszero(Rf.sys[1][:,Rf.faults]-R.sys[1][:,R.faults],atol=1.e-7) && iszero(R.sys[1][:,R.controls],atol=1.e-7) &&
      info.HDesign2[1] == [0.0 1.0 0.0; 0.0 0.0 1.0]

# solve using observer based nullspace
@time Q, Rf, info = afdisyn(sysf; nullspace = false, atol3 = 1.e-7, minimal = false, rdim = [2]); info
R = fdIFeval(Q, sysf; minimal = true, atol = 1.e-7);
@test iszero(Rf.sys[1][:,Rf.faults]-R.sys[1][:,R.faults],atol=1.e-7) && iszero(R.sys[1][:,R.controls],atol=1.e-7) &&
      info.HDesign2[1] == [0.0 1.0 0.0; 0.0 0.0 1.0]

@time Q, Rf, info = afdisyn(sysf; nullspace = false, atol3 = 1.e-7, minimal = true, rdim = [1]); info
R = fdIFeval(Q, sysf; minimal = true, atol = 1.e-7);
@test iszero(Rf.sys[1][:,Rf.faults]-R.sys[1][:,R.faults],atol=1.e-7) && iszero(R.sys[1][:,R.controls],atol=1.e-7) &&
      info.HDesign2[1] == [0.0 1.0 0.0]

SFDI = fdigenspec(sysf, atol = 1.e-7)
@time Q, Rf, info = afdisyn(sysf, SFDI; nullspace = false, atol3 = 1.e-7, minimal = true, rdim = ones(Int,3)); info
R = fdIFeval(Q, sysf; minimal = true, atol = 1.e-7);
@test iszero(vcat(Rf.sys...)-vcat(R.sys...)[:,R.faults],atol=1.e-7) && iszero(vcat(R.sys...)[:,R.controls],atol=1.e-7) &&
      info.HDesign2 == [[0.0 1.0 0.0], [0.0 1.0], [0.0 1.0]] &&
      fditspec(Rf, FDtol = 1.e-6) == SFDI  && fdisspec(Rf, FDGainTol = 1.e-3) == SFDI

SFDI = fdigenspec(sysf, atol = 1.e-7)
@time Q, Rf, info = afdisyn(sysf, SFDI[1:1,:]; nullspace = false, atol3 = 1.e-7, minimal = true, rdim = ones(Int,1),HDesign = [[0.03964565583901385 0.8670124660000067 -1.7144387091257782]]); info
#@time Q, Rf, info = efdisyn(sysf, SFDI[1:1,:]; nullspace = false, atol3 = 1.e-7, minimal = true, rdim = ones(Int,1),HDesign = [randn(1,3)]); info
@test fdiscond(Rf,SFDI[1:1,:])[2] ≈ ones(1) && fdif2ngap(Rf,SFDI[1:1,:])[2] ≈ ones(1)

# system with noise inputs
p = 3; mf = 2; mu = 2; mw = 3;
sysf = fdimodset(rss(1,p,mf+mu+mw),controls = 1:mu, faults = mu+1:mu+mf, noise = (mu+mf+1) .+ Vector(1:mf));
SFDI = fdigenspec(sysf, atol = 1.e-7)
@time Q, Rfw, info = afdisyn(sysf, SFDI; nullspace = false, atol3 = 1.e-7, minimal = true, rdim = ones(Int,3)); info
R = fdIFeval(Q, sysf; minimal = true, atol = 1.e-7);
@test iszero(vcat(Rfw.sys...)[:,Rfw.faults]-vcat(R.sys...)[:,R.faults],atol=1.e-7) &&
      iszero(vcat(Rfw.sys...)[:,Rfw.noise]-vcat(R.sys...)[:,R.noise],atol=1.e-7) &&
      iszero(vcat(R.sys...)[:,R.controls],atol=1.e-7) &&
      info.HDesign == [[0.0 1.0 0.0], [0.0 1.0], [0.0 1.0]] &&
      fditspec(Rfw, FDtol = 1.e-6) == SFDI  && fdisspec(Rfw, FDGainTol = 1.e-3) == SFDI &&
      fdiscond(Rfw,SFDI)[1] ≈ fdiscond(R,SFDI)[1] && fdif2ngap(R,SFDI)[1] ≈ info.gap &&
      any(isinf.(fdimmperf(Rfw,2))) && all(isfinite.(fdimmperf(Rfw,Inf)))

# Example 5.11c - Solution of an AFDIP

# define s as an improper transfer function
s = rtf('s');
# define Gu(s), Gf(s) and Gw(s)
Gu = [(s+1)/(s+2); (s+2)/(s+3)];     # enter Gu(s)
Gw = [1/(s+2); 0];                   # enter Gw(s)
Gf = [(s+1)/(s+2) 0; 0 1];           # enter Gf(s)
p = 2; mu = 1; mw = 1; mf = 2;       # enter dimensions
SFDI = eye(mf) .> 0

# build the synthesis model with additive faults
sysf = fdimodset(dss([Gu Gf Gw]), c = 1:mu,f = mu .+ (1:mf), n = (mu+mf) .+ (1:mw));

Q, R, info = afdisyn(sysf, SFDI; smarg =-3, sdeg = -3)

# check synthesis conditions: Qt*[Gu Gd;I 0] = 0 and Qt*[Gf; 0] = Rf
Rt = fdIFeval(Q,sysf) # form Qt*[Gu Gd Gf;I 0 0];
@test iszero(vcat(Rt.sys...)[:,[Rt.controls;Rt.disturbances]],atol=1.e-7) &&
      iszero(vcat(Rt.sys...)[:,[Rt.faults;Rt.noise]]-vcat(R.sys...),atol=1.e-7) &&
      fdif2ngap(R,SFDI)[1] ≈ info.gap

# check weak and strong fault detectability
@test fditspec(R) == fdisspec(R) == SFDI

## Model Niemann 1998, Optim. Appl. Meth.
n = 5; mu = 2; md = 1; mf = 2; m = mu+md+mf; p = 5;
a = [
-0.0782 0.2939 0.0220 0.0208 -0.0291
0.0077 -0.0278 0.0015 -0.0015 0.0026
-1.3251 5.32608 -0.5263 0.2214 -0.4777
1.0809 -4.4452 0.3770 -0.4631 0.4032
2.1532 -8.6386 0.7811 -0.5745 0.7816]*1000;
bd = [
    0.0443
   -0.0042
    0.7910
   -0.6598
   -1.2881]*1000;
bu = [
    -0.007197 0.003005
0.003473 0.000284
1.218502 -0.032866
1.322502 0.020147
-0.082265 0.024388 ];
bf = bu;
c = eye(n); dd = zeros(n,md); du = zeros(n,mu); df = zeros(n,mf);

sys = dss(a,[bu bd bf],c,[du dd df]);
sysf = fdimodset(sys, c = 1:mu, d = mu .+ (1:md) , f = (mu+md) .+ (1:mf));

SFDI = fdigenspec(sysf) # determine achievable structure matrix
nb = size(SFDI,1)
@time Q, Rf = afdisyn(sysf, SFDI; atol = 1.e-7, rdim = 1);

# check synthesis conditions: Q*[Gu Gd;I 0] = 0 and Q*[Gf; 0] = Rf
R = fdIFeval(Q,sysf); # form Q*[Gu Gd Gf;I 0 0];
@test iszero(vcat(R.sys...)[:,[R.controls;R.disturbances]],atol=1.e-7) &&
      iszero(vcat(R.sys...)[:,R.faults]-vcat(Rf.sys...),atol=1.e-7) &&
      order.(Q.sys) == ones(Int,nb)

# check weak and strong fault detectability
@test fditspec(Rf) == fdisspec(Rf) == SFDI

## Example 1, Xiong & Saif 2000, Int. J. Robust Nonlinear Control
n = 4; mu = 2; md = 2; mf = 2; m = mu+md+mf; p = 4;
a = [ -9.9477 -0.7476 0.2632 -5.0337
      52.1659 2.7452 5.5532 -24.4221
      26.0922 2.6361 -4.1975 -19.2774
      0.0 0.0 1.0 0.0];
bu = [ 0.4422 0.1761
       3.5446 -7.5922
        -5.52 4.49
        0.0 0.0];
bd = [ 0 0; 0 1 ; 1 0; 0 0 ];
bf = [bu[:,1] zeros(n,1)];
c = [ 1 0 0 0; 0 1 0 0; 0 0 1 0; 0 1 1 1];
dd = zeros(n,md); du = zeros(n,mu); df = [zeros(n,mf-1) [0 0 1 0]'];

sys = dss(a,[bu bd bf],c,[du dd df]);
sysf = fdimodset(sys, c = 1:mu, d = mu .+ (1:md) , f = (mu+md) .+ (1:mf));

SFDI = fdigenspec(sysf) # determine achievable structure matrix
nb = size(SFDI,1)
@time Q, Rf = afdisyn(sysf, SFDI; atol = 1.e-7, rdim = 1);

# check synthesis conditions: Q*[Gu Gd;I 0] = 0 and Q*[Gf; 0] = Rf
R = fdIFeval(Q,sysf); # form Q*[Gu Gd Gf;I 0 0];
@test iszero(vcat(R.sys...)[:,[R.controls;R.disturbances]],atol=1.e-7) &&
      iszero(vcat(R.sys...)[:,R.faults]-vcat(Rf.sys...),atol=1.e-7) &&
      order.(Q.sys) == [1, 1, 2]

# check weak and strong fault detectability
@test fditspec(Rf) == fdisspec(Rf) == SFDI

## Example 2, Xiong & Saif 2000, Int. J. Robust Nonlinear Control
n = 5; mu = 3; md = 2; mf = 2; m = mu+md+mf; p = 3;
a = [
0 0 -0.0034 0 0
0 -0.041 0.0013 0 0
0 0 -1.1471 0 0
0 0 -0.0036 0 0
0 0.094 0.0057 0 -0.051];
bu = [
    -1 0 0
     0 0 0
     0 0 0.948
     0.916 -1 0
    -0.598 0 0];
bd = [
     0 1
     0.062 -0.132
     0 -7.189
     0 0
     0 0 ];
bf = bu[:,1:2];
c = [
1 0 0 0 0
0 0 1 0 0
0 0 0 1 0];
dd = zeros(p,md); du = zeros(p,mu); df = zeros(p,mu-1);

sys = dss(a,[bu bd bf],c,[du dd df]);

sysf = fdimodset(sys, c = 1:mu, d = mu .+ (1:md) , f = (mu+md) .+ (1:mf));

SFDI = fdigenspec(sysf) # determine achievable structure matrix
nb = size(SFDI,1)
@time Q, Rf = afdisyn(sysf, SFDI; atol = 1.e-7, rdim = 1);

# check synthesis conditions: Q*[Gu Gd;I 0] = 0 and Q*[Gf; 0] = Rf
R = fdIFeval(Q,sysf); # form Q*[Gu Gd Gf;I 0 0];
@test iszero(vcat(R.sys...)[:,[R.controls;R.disturbances]],atol=1.e-7) &&
      iszero(vcat(R.sys...)[:,R.faults]-vcat(Rf.sys...),atol=1.e-7) &&
      order.(Q.sys) == [1, 1, 1]

# check weak and strong fault detectability
@test fditspec(Rf) == fdisspec(Rf) == SFDI

## Yuan et al. IJC (1997)  strong synthesis
p = 3; mu = 1; md = 0; mf = 8; n = 4;
A = [
-1 1 0 0
1 -2 1 0
0 1 -2 1
0 0 1 -2  ];
Bu = [1 0 0 0]';
Bf = [
1 0 0 0 1 0 0 0
0 1 0 0 -1 1 0 0
0 0 1 0 0 -1 1 0
0 0 0 1 0 0 -1 1];
C = [
 1 0 0 0
0 0 1 0
0 0 0 1];
Du = zeros(p,mu); Df = zeros(p,mf);

sys = dss(A,[Bu Bf],C,[Du Df]);
sysf = fdimodset(sys, c = 1:mu, d = mu .+ (1:md) , f = (mu+md) .+ (1:mf));

SFDI = fdigenspec(sysf; FDfreq = 0, atol = 1.e-7) # determine achievable strong structure matrix
nb = size(SFDI,1)

@time Q, Rf = afdisyn(sysf, SFDI; FDfreq = 0, atol = 1.e-7, rdim = 1);

# check synthesis conditions: Q*[Gu Gd;I 0] = 0 and Q*[Gf; 0] = Rf
R = fdIFeval(Q,sysf); # form Q*[Gu Gd Gf;I 0 0];
@test iszero(vcat(R.sys...)[:,[R.controls;R.disturbances]],atol=1.e-7) &&
      iszero(vcat(R.sys...)[:,R.faults]-vcat(Rf.sys...),atol=1.e-7) &&
      order.(Q.sys) == [2, 2, 1, 2, 2, 1, 1, 1, 2, 2, 2, 2]

# check strong fault detectability
@test fdisspec(Rf) == SFDI

## Yuan et al. IJC (1997)  weak synthesis with stabillization
p = 3; mu = 1; md = 0; mf = 8; n = 4;
A = [
-1 1 0 0
1 -2 1 0
0 1 -2 1
0 0 1 -2  ];
Bu = [1 0 0 0]';
Bf = [
1 0 0 0 1 0 0 0
0 1 0 0 -1 1 0 0
0 0 1 0 0 -1 1 0
0 0 0 1 0 0 -1 1];
C = [
 1 0 0 0
0 0 1 0
0 0 0 1];
Du = zeros(p,mu); Df = zeros(p,mf);

sys = dss(A,[Bu Bf],C,[Du Df]);
sysf = fdimodset(sys, c = 1:mu, d = mu .+ (1:md) , f = (mu+md) .+ (1:mf));

SFDI = fdigenspec(sysf; atol = 1.e-7) # determine achievable weak structure matrix
nb = size(SFDI,1)

@time Q, Rf = afdisyn(sysf, SFDI; atol = 1.e-7, rdim = 1);

# check synthesis conditions: Q*[Gu Gd;I 0] = 0 and Q*[Gf; 0] = Rf
R = fdIFeval(Q,sysf); # form Q*[Gu Gd Gf;I 0 0];
@test iszero(vcat(R.sys...)[:,[R.controls;R.disturbances]],atol=1.e-7) &&
      iszero(vcat(R.sys...)[:,R.faults]-vcat(Rf.sys...),atol=1.e-7) &&
      order.(Q.sys) == [2, 2, 1, 2, 2, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2] &&
      order(vcat(Q.sys...)) == 32

# check weak  fault detectability
@test fditspec(Rf; atol = 1.e-6, FDtol = 0.0001) == SFDI

## Yuan et al. IJC (1997)  weak synthesis with pole assignment and
#  enforcing least global order of 6!
p = 3; mu = 1; md = 0; mf = 8; n = 4;
A = [
-1 1 0 0
1 -2 1 0
0 1 -2 1
0 0 1 -2  ];
Bu = [1 0 0 0]';
Bf = [
1 0 0 0 1 0 0 0
0 1 0 0 -1 1 0 0
0 0 1 0 0 -1 1 0
0 0 0 1 0 0 -1 1];
C = [
 1 0 0 0
0 0 1 0
0 0 0 1];
Du = zeros(p,mu); Df = zeros(p,mf);

sys = dss(A,[Bu Bf],C,[Du Df]);
sysf = fdimodset(sys, c = 1:mu, d = mu .+ (1:md) , f = (mu+md) .+ (1:mf));

SFDI = fdigenspec(sysf; atol = 1.e-7) # determine achievable weak structure matrix
nb = size(SFDI,1)

@time Q, Rf = afdisyn(sysf, SFDI; smarg = -5, sdeg = -5, poles = [-5, -5], atol = 1.e-7, rdim = 1);

# check synthesis conditions: Q*[Gu Gd;I 0] = 0 and Q*[Gf; 0] = Rf
R = fdIFeval(Q,sysf); # form Q*[Gu Gd Gf;I 0 0];
@test iszero(vcat(R.sys...)[:,[R.controls;R.disturbances]],atol=1.e-7) &&
      iszero(vcat(R.sys...)[:,R.faults]-vcat(Rf.sys...),atol=1.e-7) &&
      order.(Q.sys) == [2, 2, 1, 2, 2, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2] &&
      order(vcat(Q.sys...)) == 32 && order(gbalmr(vcat(Q.sys...),atol=1.e-7)[1]) == 6 &&
      count(ghanorm(vcat(Q.sys...))[2] .> 1.e-7) == 6 &&
      fditspec(Rf) == SFDI

@test all(iszero.(gbalmr(Q; balance = true).sys - Q.sys, atol=1.e-7) )

## Example with solvable strong synthesis
#
s = rtf('s');
# define 2x2 Gf(s) without zeros
Gf = [s/(s+1) 0; (s^2+4)/(s+3)^2 1]; # enter Gf(s)
sysf = fdimodset(dss(Gf; minimal = true, atol = 1.e-7),faults = 1:2);

# weak achievable specifications
SFDI = fdigenspec(sysf)

@time Q, Rf, info = afdisyn(sysf, SFDI; sdeg = -3, rdim = 1); info

# check synthesis conditions: Q*[Gu Gd;I 0] = 0 and Q*[Gf; 0] = Rf
R = fdIFeval(Q,sysf); # form Q*[Gu Gd Gf;I 0 0];
@test iszero(vcat(R.sys...)[:,[R.controls;R.disturbances]],atol=1.e-7) &&
      iszero(vcat(R.sys...)[:,R.faults]-vcat(Rf.sys...),atol=1.e-7) &&
      fditspec(Rf) == SFDI

# the same results  must result using the resulting design matrix H
Q1, Rf1, info1 = afdisyn(sysf, SFDI; sdeg = -3, rdim = 1, HDesign = info.HDesign, HDesign2 = info.HDesign2); info1

@test all(iszero.(Q.sys .- Q1.sys,atol=1.e-7)) && all(iszero.(Rf.sys .- Rf1.sys,atol=1.e-7))

# build minimal realization of the model with faults
sysf = fdimodset(dss(Gf; minimal = true, atol = 1.e-7),faults = 1:2);

# strong achievable specifications
SFDI = fdigenspec(sysf, FDfreq = [0, 2], FDGainTol = 1.e-3)

@time Q, Rf, info = afdisyn(sysf, SFDI; sdeg = -3, rdim = 1, FDfreq = [0, 2]); info

# check synthesis results
R = fdIFeval(Q, sysf; minimal = true, atol = 1.e-7);
@test iszero(vcat(R.sys...)[:,[R.controls;R.disturbances]],atol=1.e-7) &&
      iszero(vcat(R.sys...)[:,R.faults]-vcat(Rf.sys...),atol=1.e-7) &&
      fdisspec(Rf,[0,2]) == SFDI

# the same results  must result using the resulting design matrix H
Q1, Rf1, info1 = afdisyn(sysf, SFDI; sdeg = -3, rdim = 1, FDfreq = [0, 2], FDGainTol = 0.001, HDesign = info.HDesign,  HDesign2 = info.HDesign2); info1

@test all(iszero.(Q.sys .- Q1.sys,atol=1.e-7)) && all(iszero.(Rf.sys .- Rf1.sys,atol=1.e-7))

## Example with solvable strong synthesis (frequencies are also zeros)
#  Using simple basis based synthesis
#
s = rtf('s');
# define 2x1 Gf(s) without zeros
Gf = [s/(s+1) 0; (s^2+4)/(s+3)^2 1]; # enter Gf(s)

# build minimal realzation of the model with faults
sysf = fdimodset(dss([Gf I]; minimal = true, atol = 1.e-7),faults = 1:2, noise = 3:4);

# strong achievable specifications
SFDI = fdigenspec(sysf, FDfreq = [0, 2], FDGainTol = 1.e-3)

@time Q, Rf, info = afdisyn(sysf, SFDI; sdeg = -3, rdim = 1, FDfreq = [0, 2], simple = true); info

# check synthesis results
R = fdIFeval(Q, sysf; minimal = true, atol = 1.e-7);
@test iszero(vcat(R.sys...)[:,[R.controls;R.disturbances]],atol=1.e-7) &&
      iszero(vcat(R.sys...)-vcat(Rf.sys...),atol=1.e-7) &&
      fdisspec(Rf,[0,2]) == SFDI

Q1, Rf1, info1 = afdisyn(sysf, SFDI; sdeg = -3, rdim = 1, FDfreq = [0, 2], simple = true, HDesign = info.HDesign, HDesign2 = info.HDesign2); info1

@test all(iszero.(Q.sys .- Q1.sys,atol=1.e-7)) && all(iszero.(Rf.sys .- Rf1.sys,atol=1.e-7))

## Example with solvable strong synthesis (frequencies are also poles)
#
s = rtf('s');
# define 2x1 Gf(s) without zeros
Gf = [1/s; 1/(s^2+4)]; # enter Gf(s)

# build minimal realzation of the model with faults
sysf = fdimodset(dss([Gf eye(2)]; minimal = true, atol = 1.e-7),faults = 1, noise = 2:3);

# strong achievable specifications
SFDI = fdigenspec(sysf, FDfreq = [0, 2], FDGainTol = 1.e-3);
SFDI = [SFDI;SFDI];

@time Q, Rf, info = afdisyn(sysf, SFDI; sdeg = -3, rdim = [1,2], FDfreq = [0, 2], simple = true); info

# check synthesis results
R = fdIFeval(Q, sysf; minimal = true, atol = 1.e-7);
@test iszero(vcat(R.sys...)[:,[R.controls;R.disturbances]],atol=1.e-7) &&
      iszero(vcat(R.sys...)[:,R.faults]-vcat(Rf.sys...)[:,Rf.faults],atol=1.e-7) &&
      iszero(vcat(R.sys...)[:,R.noise]-vcat(Rf.sys...)[:,Rf.noise],atol=1.e-7) &&
      fdisspec(Rf,[0,2]) == SFDI

Q1, Rf1, info1 = afdisyn(sysf, SFDI; sdeg = -3, rdim = [1,2], FDfreq = [0, 2], simple = true, HDesign = info.HDesign); info1

@test all(iszero.(Q.sys .- Q1.sys,atol=1.e-7)) && all(iszero.(Rf.sys .- Rf1.sys,atol=1.e-7))

## Example with solvable strong synthesis (frequencies are also poles)
#  Using simple basis based synthesis
s = rtf('s');
# define 2x1 Gf(s) without zeros
Gf = [1/s; 1/(s^2+4)]; # enter Gf(s)

# build minimal realzation of the model with faults
sysf = fdimodset(dss(Gf; minimal = true, atol = 1.e-7),faults = 1);

# strong achievable specifications
SFDI = fdigenspec(sysf, FDfreq = [0, 2], FDGainTol = 1.e-3);
SFDI = [SFDI;SFDI];

@time Q, Rf, info = afdisyn(sysf, SFDI; sdeg = -3, rdim = [1,2], FDfreq = [0, 2], simple = true); info

# check synthesis results
R = fdIFeval(Q, sysf; minimal = true, atol = 1.e-7);
@test iszero(vcat(R.sys...)[:,[R.controls;R.disturbances]],atol=1.e-7) &&
      iszero(vcat(R.sys...)[:,R.faults]-vcat(Rf.sys...),atol=1.e-7) &&
      fdisspec(Rf,[0,2]) == SFDI

Q1, Rf1, info1 = afdisyn(sysf, SFDI; sdeg = -3, rdim = [1,2], FDfreq = [0, 2], simple = true, HDesign2 = info.HDesign2); info1

@test all(iszero.(Q.sys .- Q1.sys,atol=1.e-7)) && all(iszero.(Rf.sys .- Rf1.sys,atol=1.e-7))

## Example with solvable strong synthesis (frequencies are also zeros)

s = rtf('s');
# define 2x1 Gf(s) without zeros
Gf = [s; s^2+4]; # enter Gf(s)

# build minimal realzation of the model with faults
sysf = fdimodset(dss(Gf; minimal = true, atol = 1.e-7),faults = 1);

# strong achievable specifications
SFDI = fdigenspec(sysf, FDfreq = [0, 2], FDGainTol = 1.e-3);
SFDI = [SFDI;SFDI];

separate = true
@time Q, Rf, info = afdisyn(sysf, SFDI; atol = 1.e-7, separate, smarg = -3, sdeg = -3, rdim = [1,2], FDGainTol=1.e-5, FDfreq = [0, 2], simple = true); info

# check synthesis results
R = fdIFeval(Q, sysf; minimal = true, atol = 1.e-7);
@test iszero(vcat(R.sys...)[:,[R.controls;R.disturbances]],atol=1.e-7) &&
      iszero(vcat(R.sys...)[:,R.faults]-vcat(Rf.sys...),atol=1.e-7) &&
      fdisspec(Rf,[0,2]; FDGainTol=1.e-5) == SFDI

Q1, Rf1, info1 = afdisyn(sysf, SFDI; atol = 1.e-7, separate, smarg = -3, sdeg = -3, rdim = [1,2], FDGainTol=1.e-5, FDfreq = [0, 2], simple = true, HDesign2 = info.HDesign2); info1
#Q1, Rf1, info1 = efdisyn(sysf, SFDI; atol = 1.e-7, smarg = -3, sdeg = -3, rdim = [1,2], FDGainTol=1.e-5, FDfreq = [0, 2], simple = false, HDesign = info.HDesign); info1
R1 = fdIFeval(Q1, sysf; minimal = true, atol = 1.e-7);
@test iszero(vcat(R1.sys...)[:,[R1.controls;R1.disturbances]],atol=1.e-7) &&
      iszero(vcat(R1.sys...)[:,R1.faults]-vcat(Rf1.sys...),atol=1.e-7) &&
      fdisspec(Rf1,[0,2]; FDGainTol=1.e-5) == SFDI

@test all(iszero.(Q.sys .- Q1.sys,atol=1.e-7)) && all(iszero.(Rf.sys .- Rf1.sys,atol=1.e-7))

# example with a weak AFDIP
s = rtf('s');
# define 1x2 Gf(s) without finite zeros
Gf = [1/s 1/(s^2+4)]; # enter Gf(s)

# build minimal realzation of the model with faults
sysf = fdimodset(dss(Gf; minimal = true, atol = 1.e-7), faults = 1:2);
SFDI = eye(2) .> 0

@time Q, Rf, info = afdisyn(sysf, SFDI; separate = true, sdeg = -3, rdim = 1); info
# check synthesis results
R = fdIFeval(Q, sysf; minimal = true, atol = 1.e-7);

@test iszero(vcat(R.sys...)[:,[R.controls;R.disturbances]],atol=1.e-7) &&
      iszero(vcat(R.sys...)[:,R.faults]-vcat(Rf.sys...),atol=1.e-7) &&
      fditspec(Rf) == (ones(2,2) .> 0) &&
      fdif2ngap(Rf,SFDI)[1] ≈ info.gap ≈ fdif2ngap(R,SFDI)[1] &&
      fdimmperf(Rf) == zeros(2) && fdimmperf(Rf,SFDI) ≈ ones(2)

Q1, Rf1, info1 = afdisyn(sysf, SFDI; separate = true, sdeg = -3, rdim = 1, HDesign = info.HDesign); info1

@test all(iszero.(Q.sys .- Q1.sys,atol=1.e-7)) && all(iszero.(Rf.sys.-Rf1.sys,atol=1.e-7))

@time Q, Rf, info = afdisyn(sysf, SFDI; separate = false, sdeg = -3, rdim = 1); info
# check synthesis results
R = fdIFeval(Q, sysf; minimal = true, atol = 1.e-7);

@test iszero(vcat(R.sys...)[:,[R.controls;R.disturbances]],atol=1.e-7) &&
      iszero(vcat(R.sys...)[:,R.faults]-vcat(Rf.sys...),atol=1.e-7) &&
      fditspec(Rf) == (ones(2,2) .> 0) &&
      fdif2ngap(Rf,SFDI)[1] ≈ info.gap ≈ fdif2ngap(R,SFDI)[1] &&
      fdimmperf(Rf) == zeros(2) && fdimmperf(Rf,SFDI) ≈ ones(2)

Q1, Rf1, info1 = afdisyn(sysf, SFDI; separate = false, sdeg = -3, rdim = 1, HDesign = info.HDesign); info1

@test all(iszero.(Q.sys .- Q1.sys,atol=1.e-7)) && all(iszero.(Rf.sys.-Rf1.sys,atol=1.e-7))

# Example - Solution of an approximate fault detection problem (AFDIP)
# Example 5.3 of (V,2017) (modified)
s = rtf('s'); # define the Laplace variable s
mu = 1; mw = 1; p = 2; mf = mu+p; # set dimensions
Gu = [(s+1)/(s+2); (s+2)/(s+3)];
# enter Gu(s)
Gw = [(s-1)/(s+2); 0];
# enter Gw(s)
# build the model with additive faults having Gf(s) = [Gu(s) eye(p)];
sysf = fdimodset(dss([Gu Gw]),c=1,f = 1,fs = 1:2,n = 2);
# select SFDI
S = fdigenspec(sysf); SFDI = S[(sum(S,dims = 2).==2)[:],:]
nb = size(SFDI,1);
# perform synthesis with AFDISYN
Q, Rfw, info = afdisyn(sysf, SFDI; atol = 1.e-7, smarg = -3, sdeg = -3)

# check synthesis results
R = fdIFeval(Q, sysf; minimal = true, atol = 1.e-7);

@test iszero(vcat(R.sys...)[:,[R.controls;R.disturbances]],atol=1.e-7) &&
      iszero(vcat(R.sys...)[:,R.faults]-vcat(Rfw.sys...)[:,Rfw.faults],atol=1.e-7) &&
      iszero(vcat(R.sys...)[:,R.noise]-vcat(Rfw.sys...)[:,Rfw.noise],atol=1.e-7) &&
      fditspec(Rfw) == SFDI &&
      fdif2ngap(Rfw,SFDI)[1] ≈ info.gap ≈ fdif2ngap(R,SFDI)[1]  &&
      norm(fdimmperf(Rfw) - fdimmperf(Rfw,SFDI)) < 1.e-7 &&
      all(fdimmperf(R,Rfw) .< 1.e-7)

end # test fdisyn

end # module

max_stars_repo_name: andreasvarga/FaultDetectionTools

max_issues_repo_name: andreasvarga/FaultDetectionTools

max_forks_repo_name: andreasvarga/FaultDetectionTools


text:
"""
    AnovaValue <: AnovaEffect

A set of values for an Anova item for which a mean square is not required.

`name` - the name of this value

`ss` - sum of squares

`df` - degrees of freedom
"""
struct AnovaValue <: AnovaEffect
    name::String
    ss::Float64
    df::Float64
end

AnovaValue(ss, df) = AnovaValue("", ss, df)

import Base.isapprox
isapprox(x::AnovaValue, y::AnovaValue) =
    x.name == y.name &&
    x.ss ≈ y.ss &&
    x.df == y.df

export AnovaValue


text:
using Test: @testset, @test, @test_throws
using Dates: year, month, day

@eval Main TRACK_TODOS = true
using Todo

@testset "Todo" begin

@testset "TodoTask" begin
	task = TodoTask("x (A) 2019-05-09 2003-01-23 +PROJ1 this is the description +PROJ2 @context/file.jl:23 tag:a-value tag2:another-value not-tag: @ctx_2 :not-tag +PROJ3")
	@test isComplete(task) === true
	@test priority(task) == 'A'
	@test year(completed(task)) == 2019 && month(completed(task)) == 5 && day(completed(task)) == 9
	@test year(created(task)) == 2003 && month(created(task)) == 1 && day(created(task)) == 23
	@test description(task) == "+PROJ1 this is the description +PROJ2 @context/file.jl:23 tag:a-value tag2:another-value not-tag: @ctx_2 :not-tag +PROJ3"
	@test projects(task) == ["PROJ1", "PROJ2", "PROJ3"]
	@test contexts(task) == ["context/file.jl:23", "ctx_2"]
	@test tags(task) == ["tag" => "a-value", "tag2" => "another-value"]

	@test_throws ErrorException TodoTask("x x this is a bad description")
	@test_throws ErrorException TodoTask("(A) (B) this is another bad description")
	@test_throws ErrorException TodoTask("x 2000-01-02 2000-01-01 2000-01-01 this is also a bad description")
	@test_throws ErrorException TodoTask("2000-01-01 2000-01-01 this is not marked complete, but has a completion date")
end

@testset "load_todos/save_todos" begin
	tasks = [
		TodoTask("(A) 2003-01-23 my first task in +PROJ1 @context/file.jl:1 tag1:value1"),
		TodoTask("x (B) 2019-05-13 2003-01-23 my second task in +PROJ2 @/context/file.jl:2 tag2:value2"),
		TodoTask("(Z) my last task in +PROJ1 +PROJ2 @../another/context/file.jl:3 tag3:value3"),
	]

	mktempdir() do dir
		todotxt = joinpath(dir, "todo.txt")
		todomd  = joinpath(dir, "todo.md")

		save_todos(todotxt, tasks)
		saved = String(read(todotxt))
		@test occursin("(A) 2003-01-23 my first task in +PROJ1 @context/file.jl:1 tag1:value1", saved)
		@test occursin("x (B) 2019-05-13 2003-01-23 my second task in +PROJ2 @/context/file.jl:2 tag2:value2", saved)
		@test occursin("(Z) my last task in +PROJ1 +PROJ2 @../another/context/file.jl:3 tag3:value3", saved)

		loadedTasks = load_todos(todotxt)
		@test tasks == loadedTasks

		save_todos(todotxt, loadedTasks)
		@test String(read(todotxt)) == saved

		save_todos(todomd, tasks, format = :markdown)
		saved = String(read(todomd))
		@test occursin("- [ ] my first task in +PROJ1 @context/file.jl:1 tag1:value1", saved)
		@test occursin("- [x] my second task in +PROJ2 @/context/file.jl:2 tag2:value2", saved)
		@test occursin("- [ ] my last task in +PROJ1 +PROJ2 @../another/context/file.jl:3 tag3:value3", saved)
	end
end

@testset "@todo_str" begin
	before = length(todos())
	@eval Main module TodoTest
		using Todo

		todo"provide more module contents"
		todo"""
			provide a
			multi-line
			todo string
			"""

		function f()
			todo"implement this function"
			todo"""
				test
				multi-line
				in
				this function
				"""
		end

		macro m()
			todo"implement this macro"
			return :(todo"test macro generated code")
		end

		@generated function g(x)
			todo"implement this generated function with value"
			return :(todo"test generated function code")
		end
	end
	after = length(todos())

	expected = 0
	@testset "module" begin
		expected += 2  # 2 in module
		@test sum(values(todo_hits())) == expected
	end

	@testset "+ macro" begin
		@eval Main TodoTest.@m()
		expected += 1+1  # 1 in the macro, 1 in the expansion of the macro
		@test sum(values(todo_hits())) == expected
	end

	@testset "+ repeat macro" begin
		@eval Main TodoTest.@m()
		expected += 1+1  # another 1 in the macro and 1 in the macro expansion
		@test sum(values(todo_hits())) == expected
	end

	@testset "+ function" begin
		@eval Main TodoTest.f()
		expected += 2  # 2 in the function
		@test sum(values(todo_hits())) == expected
	end

	@testset "+ repeat function" begin
		@eval Main TodoTest.f()
		expected += 2  # 2 in the function again
		@test sum(values(todo_hits())) == expected
	end

	@testset "+ generated" begin
		@eval Main TodoTest.g(1)
		expected += 1+1  # 1 in generating function, 1 in generated function
		@test sum(values(todo_hits())) == expected
	end

	@testset "+ repeat generated" begin
		@eval Main TodoTest.g(1.0)
		expected += 1+1  # another 1 in generating function, 1 in generated function
		@test sum(values(todo_hits())) == expected
	end

	@testset "todos() + todo_hits()" begin
		@test before == 0
		@test after == 6
		@test length(todos()) == 8
		@test length(todo_hits()) == 8

		reset_hits()
		@test length(todos()) == 8
		@test length(todo_hits()) == 0
		@test sum(values(todo_hits())) == 0
	end
end

end


text:
using MAT
using JuMP
using CPLEX

"""
    opti()

Solve the multi-time-steps optimization problem for ANM of DSO
"""
function opti()
	data = matread("./data/data_4nodes.mat")

	tic()

	# Get data
	r = data["r"]
	x = data["x"]
	I_up = data["I_up"]
	Vmin = data["Vmin"]
	Vmax = data["Vmax"]
	pDemand = data["pDemand"]
	pSolar = data["pSolarMax"]
	Ccurt = 300
	Cunfeas = 8300

	solver = CplexSolver()
	m = Model(solver=solver)

	@variable(m, v[1:5,1:169]>=0)
	# for i = 1:169
	# 	JuMP.fix(v[1,i], 0.0529)
	# end
	@variable(m, l[2:5,1:169]>=0)
	@variable(m, P[1:5,1:169])
	@variable(m, Q[1:5,1:169])
	@variable(m, q[1:5,1:169])
	@variable(m, p[1:5,1:169])
	@variable(m, p_prod[2:5,1:169]>=0)
	for i = 1:5
		for j = 1:169
			setlowerbound(q[i,j], -0.01)
			setupperbound(q[i,j], 0.01)
		end
	end

	@constraint(m, [t=1:169], P[1,t] == 0)
	@constraint(m, [t=1:169], Q[1,t] == 0)
	@constraint(m, [i=1:5,t=1:169], P[i,t] == p[i,t] + sum(P[j,t]-r[j-1]*l[j,t] for j = i+1:5))
	@constraint(m, [i=1:5,t=1:169], Q[i,t] == q[i,t] + sum(Q[j,t]-x[j-1]*l[j,t] for j = i+1:5))
	@constraint(m, [i=2:5,t=1:169], v[i-1,t] == v[i,t] -2*(r[i-1]*P[i,t]+x[i-1]*Q[i,t])+(x[i-1]^2+r[i-1]^2)*l[i,t])
	@constraint(m, [i=2:5,t=1:169], v[i,t]*l[i,t] >= P[i,t]^2 + Q[i,t]^2)
	@constraint(m, [i=2:5,t=1:169], p_prod[i,t] <= pSolar[t,i-1])
	@constraint(m, [i=2:5,t=1:169], p[i,t] == p_prod[i,t] - pDemand[t,i-1])
	@constraint(m, [i=1:5,t=1:169], v[i,t] <= Vmax[i])
	@constraint(m, [i=1:5,t=1:169], v[i,t] >= Vmin[i])
	@constraint(m, [i=2:5,t=1:169], l[i,t] <= I_up[i-1])

	@objective(m, Min, Ccurt*sum(pSolar[t,i-1]-p_prod[i,t] for i=2:5,t=1:169))

	status = solve(m)

	p_prod_value = getvalue(p_prod)
	print(sum(pSolar[t,i-1]-p_prod_value[i,t] for i=2:5,t=1:169))

	toc()
end

opti()

max_stars_repo_name: qlete/ANManagement

max_issues_repo_name: qlete/ANManagement

max_forks_repo_name: qlete/ANManagement


text:
# Settings to be used by the scripts in the 'paper' folder

BASE_FOLDERNAME = "C:\\Users\\charl\\Documents\\Current_work\\Learning_controllers\\PkgSRA\\"
FIGURE_FOLDERNAME = BASE_FOLDERNAME*"paper\\figures\\"
EXAMPLE_FOLDERNAME = BASE_FOLDERNAME*"examples\\"
DAT_FOLDERNAME = BASE_FOLDERNAME*"paper\\dat\\"

# COLOR_DICT = Dict("true"=>[0, 0, 0],
#               "data"=>[111, 111, 111],
#               "intrinsic_true"=>[0, 0, 140],
#               "intrinsic"=>[0, 0, 227],
#               "control_true"=>[140, 0, 0],
#               "control_time"=>[227, 0, 0],
#               "control_space"=>[255, 128, 0])
# COLOR_DICT = Dict("true"=>:black, # Should be dashed
#               "data"=>:grey,
#               "model_uncontrolled"=>:teal,
#               "model_controlled"=>:deepskyblue,
#               "prediction"=>:black,
#               "intrinsic_true"=>:mediumblue,
#               "intrinsic"=>:deepskyblue,
#               "residual"=>:blue,
#               "control_true"=>:red,
#               "control_time"=>:darkorange,
#               "control_space"=>:darkorange)

  COLOR_DICT = Dict("true"=>:black, # Should be dashed
                "data"=>:grey,
                "data_uncontrolled"=>:grey,
                "model_uncontrolled"=>:blue,
                "model_partial"=>:blue,
                "model_controlled"=>:purple,
                "prediction"=>:black,
                "intrinsic_true"=>:blue,
                "intrinsic"=>:blue,
                "residual"=>:blue,
                "control_true"=>:black,
                "control_time"=>:red,
                "control_space"=>:red)


text:
using MODELe
using Test

@testset "MODELe.jl" begin
    # Write your own tests here.
end


text:
include("../common.jl")

build_libcurl(ARGS, "LibCURL")


text:
function PreconditionedIterativeLinearSolver(arg0::IterationManager)
    return PreconditionedIterativeLinearSolver((IterationManager,), arg0)
end

function PreconditionedIterativeLinearSolver(arg0::jint)
    return PreconditionedIterativeLinearSolver((jint,), arg0)
end

function solve_in_place(obj::PreconditionedIterativeLinearSolver, arg0::RealLinearOperator, arg1::RealLinearOperator, arg2::RealVector, arg3::RealVector)
    return jcall(obj, "solveInPlace", RealVector, (RealLinearOperator, RealLinearOperator, RealVector, RealVector), arg0, arg1, arg2, arg3)
end


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

"""
    varshape(f::Function)::AbstractValueShape

Get the value shape of the input/argument of an unary function `f`. `f`
should support call syntax

    f(x::T)

with `valshape(x) == varshape(f)` as well all

    f(x::AbstractVector{<:Real})

with `length(eachindex(x)) == vardof(f)` (see [`vardof`](@ref)).
"""
function varshape end
export varshape

"""
    vardof(f::Function)::Integer

Get the number of degrees of freedom of the input/argument of f.

Equivalent to `totalndof(varshape(f))` (see [`varshape`](@ref)).
"""
function vardof end
export vardof

vardof(f::Function) = totalndof(varshape(f))

@inline _maybe_shaped(x::Any, vs::Nothing) = x
Base.@propagate_inbounds _maybe_shaped(x::Any, vs::AbstractValueShape) = stripscalar(vs(x))

@inline _maybe_unshaped(x::Any, vs::Nothing) = x
Base.@propagate_inbounds _maybe_unshaped(x::Any, vs::AbstractValueShape) = unshaped(x, vs)

struct UnshapedFunction{F<:Function,IS<:Union{Nothing,AbstractValueShape},OS<:Union{Nothing,AbstractValueShape}} <: Function
    orig_f::F
    orig_varshape::IS
    orig_valshape::OS
end

"""
    unshaped(
        f::Function,
        orig_varshape::Union{Nothing,AbstractValueShape},
        orig_valshape::Union{Nothing,AbstractValueShape} = nothing
    )

Return a function that
    * Shapes it's input from a flat vector of `Real` using `orig_varshape`, if not `nothing`.
    * Calls `f` with the (optionally) unshaped input.
    * Unshapes the result to a flat vector of `Real` using `orig_valshape`, if not `nothing`.
"""
function unshaped(f::Function, orig_varshape::Union{Nothing,AbstractValueShape}, orig_valshape::Union{Nothing,AbstractValueShape} = nothing)
    UnshapedFunction(f, orig_varshape, orig_valshape)
end

varshape(uf::UnshapedFunction) = uf.varshape

Base.@propagate_inbounds function (uf::UnshapedFunction)(x)
    orig_x =_maybe_shaped(x, uf.orig_varshape)
    orig_y = uf.orig_f(orig_x)
    _maybe_unshaped(orig_y, uf.orig_valshape)
end

# ToDo: Deprecate/remove:
struct FuncWithVarShape{F<:Function,VS<:AbstractValueShape} <: Function
    f::F
    varshape::VS
end

varshape(fws::FuncWithVarShape) = fws.varshape

Base.@propagate_inbounds (fws::FuncWithVarShape)(x::Any) = fws.f(x)

Base.@propagate_inbounds function (fws::FuncWithVarShape)(x::AbstractVector{<:Real})
    # ToDo: Check performance:
    x_shaped = stripscalar(varshape(fws)(x))
    fws.f(x_shaped)
end

# ToDo: Deprecate/remove:
import Base.>>
@inline >>(varshape::AbstractValueShape, f::Function) = FuncWithVarShape(f, varshape)