text:
module AxisArrayConversion
using Requires
export to

const PUBLIC = "This function is part of the public API and subject to SemVer guarantees."
const PRIVATE = "This function is considered private and may silently change or be removed."

"""
    namedtuple(arr)::NamedTuple{(:axes, :values)}

Create a `NamedTuple` from an array `arr` with axis data.

```jldoctest
julia> using AxisArrayConversion: namedtuple

julia> namedtuple([1 2; 3 4])
(axes = (x1 = Base.OneTo(2), x2 = Base.OneTo(2)), values = [1 2; 3 4])
```

$PUBLIC
"""
function namedtuple end

"""
    from_namedtuple(::Type{T}, nt::NamedTuple)

Construct an instance of `T` from a `NamedTuple` `nt` with propertynames `(:axes, :values)`.

$PUBLIC
"""
function from_namedtuple(::Type{T}, nt::NamedTuple)::T where {T <:NamedTuple}
    return nt::T
end

"""
    to(SomeArray, ::OtherArray)::SomeArray

Construct an instance of `SomeArray` from an instance of `OtherArray`.

# Examples
```jldoctest
julia> using AxisArrayConversion: SimpleAxisArray, to

julia> to(SimpleAxisArray, (axes=(a=1:3,), values=[10,20,30]))
SimpleAxisArray(axes = (a = 1:3,), values = [10, 20, 30])

julia> arr = to(SimpleAxisArray, (axes=(a=1:3,), values=[10,20,30]))
SimpleAxisArray(axes = (a = 1:3,), values = [10, 20, 30])

julia> to(NamedTuple, arr)
(axes = (a = 1:3,), values = [10, 20, 30])

julia> to(NamedTuple, [1 2; 3 4])
(axes = (x1 = Base.OneTo(2), x2 = Base.OneTo(2)), values = [1 2; 3 4])
```

This function should not be overloaded. Instead [`namedtuple`](@ref) and [`from_namedtuple`](@ref) should be overloaded.

$PUBLIC
"""
function to(::Type{T}, o) where {T}
    nt = namedtuple(o)
    return from_namedtuple(T,nt)::T
end

"""
    roottype(::Type{T})

Construct a more general variant of an AxisArray like type `T`. Typically by stripping all type parameters.
```jldoctest
julia> using AxisArrayConversion: roottype

julia> roottype(Vector{Int})
Array

julia> roottype(typeof((a=1, b=2)))
NamedTuple
```
$PUBLIC
"""
roottype(::Type{T}) where {T <: NamedTuple} = NamedTuple

@generated function roottype(::Type{T}) where {T <: AbstractArray}
    getfield(parentmodule(T), nameof(T))
end

ntfieldnames(::Type{<:NamedTuple{fnames}}) where {fnames} = fnames

function check_consistency(arr::AbstractArray)
end

"""
    check_consistency(axis_array_like)

Check that a `NamedTuple` represents an axis array like object.
Also accepts `AbstractArray` inputs, which means a noop.

# Examples
```jldoctest
julia> using AxisArrayConversion: check_consistency

julia> check_consistency((axs=(1:2,), values=[1,2]))
ERROR: ArgumentError: Unexpected properties.
Expected exactly two properties `axes`, `values`
Got: `propertynames(nt) = (:axs, :values)
[...]

julia> check_consistency((axes=(1:2,), values=[1 2; 3 4]))
ERROR: ArgumentError: Inconsistent number of dimensions.
Expected: length(nt.axes) == ndims(nt.values)
Got: length(nt.axes)  = 1
ndims(nt.values) = 2
[...]

julia> check_consistency((axes=(1:2,), values=[1,2,3]))
ERROR: ArgumentError: Inconsistent size.
Expected: map(length, Tuple(nt.axes)) == size(nt.values)
Got: map(length, Tuple(nt.axes)) = (2,)
Got: size(nt.values) = (3,)
[...]

julia> check_consistency((axes=(1:2,), values=[1,2])) # works

julia> check_consistency([1,2,3]) # noop for AbstractArray
```
$PUBLIC
"""
function check_consistency(nt::NamedTuple)
    pnames = ntfieldnames(typeof(nt))
    if (length(nt) !== 2) || !(:axes in pnames) || !(:values in pnames)
        msg = """
        Unexpected properties.
        Expected exactly two properties `axes`, `values`
        Got: `propertynames(nt) = $pnames"""
        throw(ArgumentError(msg))
    end
    ax_ndims = length(nt.axes)
    vals_ndims = ndims(nt.values)
    if ax_ndims != vals_ndims
        msg = """
        Inconsistent number of dimensions.
        Expected: length(nt.axes) == ndims(nt.values)
        Got: length(nt.axes)  = $(ax_ndims)
        ndims(nt.values) = $(vals_ndims)"""
        throw(ArgumentError(msg))
    end
    ax_size = map(length, Tuple(nt.axes))
    vals_size = size(nt.values)
    if ax_size != vals_size
        msg = """
        Inconsistent size.
        Expected: map(length, Tuple(nt.axes)) == size(nt.values)
        Got: map(length, Tuple(nt.axes)) = $ax_size
        Got: size(nt.values) = $vals_size"""
        throw(ArgumentError(msg))
    end
end
################################################################################
##### Helpers
################################################################################

@inline function axissymbol_from_Int(i::Int)
    return Symbol("x$i")
end

function axissymbols(arr::AbstractArray{T, dim}) where {T,dim}
    make_axissymbols(Val(dim))
end

"""

    make_axissymbols(::Val{dim})::NTuple{dim, Symbol}

Create default names for `dim` many axes.

$PRIVATE
"""
@generated function make_axissymbols(::Val{dim}) where {dim}
    ntuple(axissymbol_from_Int, dim)
end

"""
    name_axes(axes)::NamedTuple

Add names to `axes` or return them as is, if they already have names.

# Examples
```jldoctest
julia> using AxisArrayConversion: name_axes

julia> name_axes((1:3, ["a", "b"]))
(x1 = 1:3, x2 = ["a", "b"])

julia> name_axes((a=1:3, b=["a", "b"]))
(a = 1:3, b = ["a", "b"])
```

$PUBLIC
"""
function name_axes(axes::NamedTuple)
    return axes
end

function name_axes(axes::NTuple{N, Any}) where {N}
    names = make_axissymbols(Val(N))
    return NamedTuple{names}(axes)
end

function axisvalues(arr::AbstractArray)
    Base.axes(arr)
end

function values(arr::AbstractArray)
    arr
end

function namedaxes(o)
    return NamedTuple{axissymbols(o)}(axisvalues(o))
end
function namedaxes(o::NamedTuple)
    return name_axes(o.axes)
end
function values(o::NamedTuple)
    return o.values
end
function namedtuple(o)
    return (axes=namedaxes(o), values=values(o))
end

function __init__()
    @require AxisArrays = "39de3d68-74b9-583c-8d2d-e117c070f3a9" include("AxisArrays.jl")
    @require AxisKeys = "94b1ba4f-4ee9-5380-92f1-94cde586c3c5"   include("AxisKeys.jl")
    @require DimensionalData = "0703355e-b756-11e9-17c0-8b28908087d0" include("DimensionalData.jl")

end

################################################################################
##### SimpleAxisArray
################################################################################
"""
    SimpleAxisArray

Simple `AxisArray` like type. Meant for demonstrations, testing and as an implementing examples.

$PUBLIC
"""
struct SimpleAxisArray{T, N, NT<:NamedTuple, V <: AbstractArray} <: AbstractArray{T,N}
    axes::NT
    values::V
    function SimpleAxisArray(axes::NamedTuple, values::AbstractArray)
        check_consistency((axes=axes, values=values))
        NT = typeof(axes)
        V = typeof(values)
        T = eltype(values)
        N = ndims(values)
        return new{T,N,NT,V}(axes, values)
    end
end

SimpleAxisArray(;axes, values) = SimpleAxisArray(axes, values)
_show(io::IO, sa::SimpleAxisArray) = print(io, "SimpleAxisArray", to(NamedTuple, sa))
Base.show(io::IO, sa::SimpleAxisArray) = _show(io, sa)
Base.show(io::IO, ::MIME"text/plain", sa::SimpleAxisArray) = _show(io, sa)

Base.size(arr::SimpleAxisArray) = size(arr.values)
Base.getindex(arr::SimpleAxisArray, I...) = arr.values[I...]

function namedtuple(o::SimpleAxisArray)
    return (axes=o.axes, values=o.values)
end

function from_namedtuple(::Type{T}, nt::NamedTuple) where {T <: SimpleAxisArray}
    SimpleAxisArray(name_axes(nt.axes), nt.values)
end

include("Base.jl")

end


text:
using CSV, DataFrames, Rematch, MATLAB, StaticArrays, CoordinateTransformations, Dierckx, Unitful

abstract type AbstractExperiment end

abstract type PresentOrgNest <: AbstractExperiment end
abstract type AbsentOrgNest <: AbstractExperiment end

######################### Experiments ##############

struct ClosedNest <: AbsentOrgNest
    path::String
    name::String
end

struct TransferNest <: PresentOrgNest
    path::String
    name::String
end

struct Transfer <: AbsentOrgNest
    path::String
    name::String
end

struct TransferNestBelen <: AbsentOrgNest
    path::String
    name::String
end

struct DawaySandpaper <: PresentOrgNest
    path::String
    name::String
end

struct Daway <: PresentOrgNest
    path::String
    name::String
end

function AbstractExperiment(path::String)
    _, name = splitdir(path)
    @match name begin
        "Closed nest Therese" || "Closed nest_Belén"                                                => ClosedNest(path, name)
        "Transfer nest Therese"                                                                     => TransferNest(path, name)
        "Transfer Therese" || "Transfer pellet moved Therese"  || "Transfer pellet removed Therese" => Transfer(path, name)
        "Transfer_Belén"                                                                            => TransferNestBelen(path, name)
        "Daway_Belén" || "Awayd_closed_nest_Belén" || "Dright_Belén" || "Dleft_closednest_Belén" || "Dzerovector_Rot_Elin" || "Daway_Rot_Elin" || "Daway_S_Elin" || "Dforward_S_Elin" || "Dleft_S_Elin" || "Dright_S_Elin" => DawaySandpaper(path, name)
        "Daway Therese" || "Dleft Therese" || "Dtowards Therese" || "Dright Therese"                => Daway(path, name)
    end
end

# filter!(x -> x ∉ bad, runs)

const Coordinate = SVector{2, Float64}

struct Point
    xy::Coordinate
    i::Int
end

struct Final
    name::String
    track::Vector{Coordinate}
    time::StepRangeLen
    angle::Vector{Float64}
    length::Vector{Float64}
    speed::Vector{Float64}
    pellet::Union{Missing, Point}
    turn_point::Point
    nest::Coordinate
    feeder::Coordinate
    org_nest::Coordinate
end

struct Experiment{E <: AbstractExperiment}
    exp_type::E
    runs::Vector{Final}
end

struct Basic{E<:AbstractExperiment}
    exp_type::E
    name::String
    path::String
    Basic{E}(exp_type::E, name::String) where {E<:AbstractExperiment} = new(exp_type, name, joinpath(exp_type.path, name))
end
Basic(exp_type::E, name::String) where {E<:AbstractExperiment} = Basic{E}(exp_type, name)

#=function get_experiment(path::String)
    name, run = splitdir(path)
    _, _experiment = splitdir(name)
    match_experiment(_experiment, run, path)
end=#

# get_experiment_type(d::Final{E, <:FinalPellet}) where {E <: Experiment} = E

function Experiment(path::String)
    exp_type = AbstractExperiment(path)
    runs = Final[]
    for r in readdir(path)
        if r"^[^.]"(r)
            a = Basic(exp_type, r)
            if isdir(a.path)
                if CHECK
                    msgs = check(a)
                    if !isempty(msgs)
                        print_with_color.(:magenta, "Errors in $(a.exp_type.name)/$(a.name):\n")
                        print_with_color.(:red, msgs, '\n')
                        continue
                    end
                end
                push!(runs, Final(a))
            end
        end
    end
    Experiment(exp_type, runs)
end

abstract type Calibration end

struct Calib <: Calibration
    file::String
end
struct JPG <: Calibration
    file::String
    checker::Float64 # in cm
end
struct MATcalib <: Calibration
    file::String
end

const factor_filename = "factors.csv"

function get_calibration(path, calibfile)
    name, ext = splitext(calibfile)
    return @match ext begin
        ".calib" => Calib(joinpath(path, calibfile))
        ".mat" => MATcalib(joinpath(path, calibfile))
        ".jpg" || ".JPG" where r"calib"(name) => begin
            factors = CSV.read(joinpath(path, factor_filename), header=["Factor", "Level"], types=[String, Union{String, Missing}], transforms = Dict(1 => strip))
            d, u = split(first(factors[factors[:Factor] .== "check", :Level]))
            checker = ustrip(uconvert(Unitful.cm, parse(d)*getfield(Unitful, Symbol(u))))
            JPG(joinpath(path, calibfile), checker)
        end
    end
end
get_calibrated_data(path::String, calib::Calib, mtlbs::Vector{Dict{String, Any}}) = mxcall(:calib_calibrate, 2, joinpath(path, calib.file), mtlbs)
get_calibrated_data(path::String, calib::JPG, mtlbs::Vector{Dict{String, Any}}) = mxcall(:jpg_calibrate, 2, joinpath(path, calib.file), calib.checker, mtlbs)
get_calibrated_data(path::String, calib::MATcalib, mtlbs::Vector{Dict{String, Any}}) = mxcall(:mat_calibrate, 2, joinpath(path, calib.file), mtlbs)

const poi_filename = "columns.csv"

function get_data(a::Basic{<:AbstractExperiment})
    fname = joinpath(a.path, poi_filename)
    pois = CSV.read(fname, header=["resfile", "columnnumber", "poi"], types=[String, Int, String], transforms = Dict(1 => strip, 3 => strip))
    # pellet = "pellet" ∈ pois[:poi] ? BasicPresent(pois) : BasicAbsent(pois)
    # Basic(experiment, pellet)
    fname = joinpath(a.path, "calibrations.csv")
    calibs = CSV.read(fname, header=["resfile", "calibfile"], types=[String, String], transforms = Dict(1 => strip, 2 => strip))
    todo = join(calibs, pois, on = :resfile)
    data = Dict{Symbol,Tuple{Array{Float64,2},Array{Float64,1}}}()
    for ucalib in groupby(todo, :calibfile)
        mtlbs = Dict{String, Any}[]
        for ures in groupby(ucalib, :resfile)
            mtlb = Dict{String, Any}()
            mtlb["resfile"] = joinpath(a.path, first(ures[:resfile]))
            mtlb["pois"] = ures[:columnnumber]
            push!(mtlbs, mtlb)
        end
        calib = get_calibration(a.path, first(ucalib[:calibfile]))
        xys, ts = get_calibrated_data(a.path, calib, mtlbs)
        for (poi, xy, t) in zip(ucalib[:poi], xys, ts)
            _t = t isa Float64 ? [t] : t
            data[Symbol(poi)] = (xy, _t)
        end
    end
    data
end

const Coordinate3D = SVector{3, Float64}

# abstract type DataPellet end

# struct DataAbsent <: DataPellet
# end

struct DataPoint
    xy::Coordinate3D
    t::Float64
end

struct Data{E<:AbstractExperiment}
    exp_type::E
    name::String
    pellet::Union{DataPoint, Missing}
    track::Vector{Coordinate3D}
    time::Vector{Float64}
    data::Dict{Symbol,Tuple{Matrix{Float64},Vector{Float64}}}
end

function retrievetrack(_track::Matrix{Float64}, _time::Vector{Float64})
    track = [Coordinate3D(_track[i,1], _track[i,2], 0.0) for i in 1:size(_track, 1)]
    time = vec(_time)
    track, time
end

function retrievepoint(data::Dict{Symbol,Tuple{Matrix{Float64},Vector{Float64}}}, f::Symbol)
    x, t = data[f]
    Coordinate3D(x[1], x[2], 0.0), t[1]
end

function Data(a::Basic{<:AbstractExperiment})
    data = get_data(a)
    track, time = retrievetrack(data[:track]...)
    pellet = haskey(data, :pellet) ? DataPoint(retrievepoint(data, :pellet)...) : missing
    Data(a.exp_type, a.name, pellet, track, time, data)
end

struct Common
    name::String
    pellet::Union{DataPoint, Missing}
    track::Vector{Coordinate3D}
    time::Vector{Float64}
    feeder::Coordinate3D
    nest::Coordinate3D
    org_nest::Coordinate3D
end

######################### Common methods ###########

function Common(b::Data{ClosedNest})
    feeder, _ = retrievepoint(b.data, :feeder)
    nest, _ = retrievepoint(b.data, :nest)
    Common(b.name, b.pellet, b.track, b.time, feeder, nest, nest)
end

function Common(b::Data{TransferNestBelen})
    southbefore, _ = retrievepoint(b.data, :southbefore)
    northbefore, _ = retrievepoint(b.data, :northbefore)
    feederbefore, _ = retrievepoint(b.data, :feederbefore)
    nestbefore, _ = retrievepoint(b.data, :nestbefore)
    v = northbefore - southbefore
    u = nestbefore - feederbefore
    azimuth = atan2(v[2], v[1]) - atan2(u[2], u[1])
    nest2feeder = norm(nestbefore - feederbefore)

    south, _ = retrievepoint(b.data, :south)
    north, _ = retrievepoint(b.data, :north)
    v = north - south
    α = atan2(v[2], v[1]) + azimuth
    u = Coordinate3D(cos(α), sin(α), 0)
    feeder, _  = retrievepoint(b.data, :feeder)
    nest = feeder + u*nest2feeder
    Common(b.name, b.pellet, b.track, b.time, feeder, nest, Coordinate3D(NaN, NaN, NaN))
end

function Common(b::Data{Transfer})
    south, _ = retrievepoint(b.data, :south)
    north, _ = retrievepoint(b.data, :north)
    factors = CSV.read(joinpath(b.exp_type.path, b.name, factor_filename), header=["Factor", "Level"], types=[String, Union{String, Missing}], transforms = Dict(1 => strip))
    d, u = split(first(factors[factors[:Factor] .== "nest2feeder", :Level]))
    nest2feeder = ustrip(uconvert(Unitful.cm, parse(d)*getfield(Unitful, Symbol(u))))
    d, u = split(first(factors[factors[:Factor] .== "azimuth", :Level]))
    azimuth = ustrip(uconvert(Unitful.rad, parse(d)*getfield(Unitful, Symbol(u))))
    v = north - south
    α = atan2(v[2], v[1]) + azimuth - π
    u = Coordinate3D(cos(α), sin(α), 0)
    feeder, _ = retrievepoint(b.data, :feeder)
    nest = feeder + u*nest2feeder
    Common(b.name, b.pellet, b.track, b.time, feeder, nest, Coordinate3D(NaN, NaN, NaN))
end

function Common(b::Data{TransferNest})
    btransfer = Data(Transfer(b.exp_type.path, b.exp_type.name), b.name, b.pellet, b.track, b.time, b.data)
    ctransfer = Common(btransfer)
    org_nest, _ = retrievepoint(b.data, :originalnest)
    Common(b.name, ctransfer.pellet, ctransfer.track, ctransfer.time, ctransfer.feeder, ctransfer.nest, Coordinate3D(org_nest))
end

function Common(b::Data{DawaySandpaper})
    org_nest, _ = retrievepoint(b.data, :nest)
    initial = mean(first(retrievepoint(b.data, k)) for k in [:rightdowninitial, :leftdowninitial, :rightupinitial, :leftupinitial])
    final = mean(first(retrievepoint(b.data, k)) for k in [:rightdownfinal, :leftdownfinal, :rightupfinal, :leftupfinal])
    v = final - initial
    nest = org_nest + v
    _feeder, _ = retrievepoint(b.data, :feeder)
    feeder = _feeder + v
    Common(b.name, b.pellet, b.track, b.time, feeder, nest, org_nest)
end

function Common(b::Data{Daway})
    org_nest, _ = retrievepoint(b.data, :nest)
    initialfeeder, _ = retrievepoint(b.data, :initialfeeder)
    feeder, _ = retrievepoint(b.data, :feeder)
    v = feeder - initialfeeder
    nest = org_nest + v
    Common(b.name, b.pellet, b.track, b.time, feeder, nest, org_nest)
end

##################################################

function get_standardization(nest::Coordinate3D, feeder::Coordinate3D)
    trans = Translation(-vec(nest))
    v = feeder - nest
    rot = LinearMap(RotZ(-atan2(v[2], v[1]) - π/2))
    rot∘trans
    # f(x::Coordinate3D) = Coordinate(composed(x)[1:2])
    # return f
end

standardize(f::T, pellet::DataPoint, t1) where {T <: AbstractAffineMap} = DataPoint(f(pellet.xy), pellet.t - t1)
standardize(f::T, pellet::Missing, t1) where {T <: AbstractAffineMap} = pellet
const atol = 1e-11

function standardize(f::T, c::Common) where {T<:AbstractAffineMap}
    track = f.(c.track)
    t1 = c.time[1]
    time = c.time - t1
    pellet = standardize(f, c.pellet, t1)
    feeder = f(c.feeder)
    nest = f(c.nest)
    org_nest = f(c.org_nest)
    @assert isapprox(feeder[1], 0, atol = atol) "feeder x value: $feederX > $atol"
    return Common(c.name, pellet, track, time, feeder, nest, org_nest)
end

const smoothing_factor = 10
function smooth(track::Vector{Coordinate3D}, time::Vector{Float64})
    X = Array{Float64, 2}(2, length(track))
    for i in eachindex(track)
        X[:,i] = track[i][1:2]
    end
    Δt = mean(diff(time))
    ParametricSpline(time, X, s = 0.4smoothing_factor/Δt, k=1)
end
function outlier_fences(data, α = 0.2, q = 6)#3
    q1, q3 = quantile(data,[α, 1 - α])
    iqr = q3-q1
    major = q*iqr#10-14
    q3+major
end
function find_outlier(dang, i1, i)
    old = shift!(i)
    for new in i
        M = outlier_fences(dang[i1[1]:old[1]])
        dang[new] > M && return new
        old = new
    end
    return old
end
function findlocalmaxima(x::Vector{T}) where T <: Real
    inds = Int[]
    for i in 2:length(x) - 1
        if x[i - 1] < x[i] > x[i + 1]
            push!(inds, i)
        end
    end
    return inds
end
function detect_turningpoint(pathlength, angles)
    i1 = findfirst(pathlength .> 20)
    if i1 == 0
        i1 = 1
    end
    i2 = pathlength[end] ≤ 400 ? length(pathlength) : findfirst(pathlength .> 400)
    angs = angles[1:i2]
    dang = abs.(diff(angs))
    i = findlocalmaxima(dang)
    filter!(x -> x > i1, i)
    if isempty(i)
        push!(i, length(dang))
    end
    j = find_outlier(dang, i1, i)
    μ = atan2(sum(sin.(angles[i1:j])), sum(cos.(angles[i1:j])))
    δ = abs.(angles[i1:end] - μ)
    turningpoint_i = findfirst(δ .> π/4) + i1 - 1
    if turningpoint_i == 0
        turningpoint_i = 1
    end
    return turningpoint_i
end

function common(track, time)
    spl = smooth(track, time)
    n = 2length(time)
    times = linspace(time[1], time[end], n)
    tracks = [Coordinate(evaluate(spl, t)...) for t in times]
    angles = Vector{Float64}(n)
    pathlengths = Vector{Float64}(n)
    for i in eachindex(angles)
        dx, dy = derivative(spl, times[i])
        angles[i] = atan2(dy, dx)
        pathlengths[i] = sqrt(dx^2 + dy^2)
    end
    pathlengths .= cumsum(pathlengths*step(times))
    turningpoint_i = detect_turningpoint(pathlengths, angles)
    speeds = pathlengths./times
    speeds[1] = speeds[2]
    # range = 1:turningpoint_i
    # feeder2turning = Coordinate.(tracks[range], times[range], angles[range], pathlengths[range], speeds[range])
    tracks, times, angles, pathlengths, speeds, turningpoint_i
end

# abstract type FinalPellet end

# struct FinalAbsent <: FinalPellet
# end

# struct FinalPresent <: FinalPellet
# xy::Coordinate
# i::Int
# end

function find_pellet_i(pellet::DataPoint, time, track::Vector{Coordinate})
    m = Inf
    ind = 0
    for i in eachindex(time)
        if time[i] > pellet.t - 2
            _m = norm(track[i] - pellet.xy[1:2])
            if _m < m
                m = _m
                ind = i
            end
        elseif time[i] > pellet.t + 2
            break
        end
    end
    return Point(Coordinate(pellet.xy[1:2]), ind)
end
find_pellet_i(pellet::Missing, time, track::Vector{Coordinate}) = pellet

larger(i::Int, p::Point) = i ≥ p.i
larger(i::Int, p::Missing) = false
function Final(c::Common)
    f = get_standardization(c.nest, c.feeder)
    d = standardize(f, c)
    track, time, angle, pathlength, speed, turningpoint_i = common(d.track, d.time)
    pellet = find_pellet_i(d.pellet, time, track)
    turningpoint_i = larger(turningpoint_i, pellet) ? pellet.i - 1 : turningpoint_i
    Final(c.name, track, time, angle, pathlength, speed, pellet, Point(track[turningpoint_i], turningpoint_i), Coordinate(d.nest[1:2]), Coordinate(d.feeder[1:2]), Coordinate(d.org_nest[1:2]))
end

Final(a::Basic) = Final(Common(Data(a)))

ext: jl


text:
function trunc_fixed(x)
if x < 0
[ 30 ]
return zero(x)
else
return x
end
end

max_stars_repo_path: Module 2/Chapter03/trunc_fixed.jl

max_issues_repo_path: Module 2/Chapter03/trunc_fixed.jl

max_forks_repo_path: Module 2/Chapter03/trunc_fixed.jl


text:
@rule GammaMixture((:b, k), Marginalisation) (q_out::Any, q_switch::Any, q_a::Any) = begin
    π_k = probvec(q_switch)[k]
    return GammaShapeRate(1 + π_k * mean(q_a), π_k * mean(q_out))
end


text:
"""
    PInterval <: SimpleProbabilitySet

Credal set of categorical distributions defined by lower and upper probabilities.
"""
struct PInterval <: SimpleProbabilitySet
    plower::Vector{Float64}
    pupper::Vector{Float64}
end

"""
    PInterval(p::Vector{Float64}, uncsize::Float64)

Probability intervals centered on 'p' of width ± 'uncsize' (with saturated bounds)
"""
function PInterval(p::Vector{Float64}, uncsize::Float64)
    pl = zeros(length(p))
    pu = ones(length(p))
    for (i, pi) in enumerate(p)
        pl[i] = max(0.0, pi - uncsize)
        pu[i] = min(1.0, pi + uncsize)
    end
    PInterval(pl, pu)
end

"""
    PInterval(p::Vector{Float64})

Degenerate probability interval of only pmf 'p'
"""
PInterval(p::Vector{Float64}) = PInterval(p,p) # degenerate PInterval

"""
    PInterval()

Default probability interval of p1 ∈ [0.2, 0.3], p2 ∈ [0.7, 0.8]
"""
PInterval() = PInterval([0.2, 0.7], [0.3, 0.8]) # default PInterval

"""
    psample(plower::Vector{Float64}, pupper::Vector{Float64})

Generate a probability distribution, uniformly at random, from the ambiguity set defined by the `pupper` and `plower` probabilities.
"""
function psample(plower::Vector{Float64}, pupper::Vector{Float64})
    n = length(plower)
    A = zeros(Float64, 2n+1, n)
    for i = 1:n
        A[2i - 1,i] = 1.0
        A[2i, i] = -1.0
    end
    A[end,:] = ones(n)
    bc = Array{Float64}(2n+1)
    for i = 1:n
        # due to limitations in hitandrun, set upper limit to 0+eps if upper limit is actually 0.0
        bc[2i - 1] = pupper[i] == 0.0 ? pupper[i] + 0.0000001 : pupper[i]
        bc[2i] = -1.0 * plower[i]
    end
    bc[end] = 1.0
    d = vcat(fill("<=", 2n),"=")
    nsample = 1
    rseed = rand(1:10_000)
    rout = R"""
    library(hitandrun)
    set.seed($rseed)
    constr <- list(constr = $A, rhs = $bc, dir = $d)
    samples <- hitandrun(constr, n.samples = $nsample, thin = ($n) ^ 3)
    """
    p = reshape(rcopy(rout), n)
end

psample(pint::PInterval) = psample(pint.plower, pint.pupper)


text:
using ParameterJuMP
using MathOptInterface
const MOI = MathOptInterface

export BendersBridge

function parametrized_function(var::StructuredVariableRef,
                               structured_model::StructuredModel,
                               model::JuMP.Model,
                               variable_map::Dict{Int, JuMP.VariableRef},
                               parameter_map::Dict{Int, ParameterJuMP.ParameterRef})
    if var.model === structured_model
        return variable_map[var.idx]
    elseif var.model === structured_model.parent
        if !haskey(parameter_map, var.idx)
            parameter_map[var.idx] = ParameterJuMP.add_parameter(model)
        end
        return parameter_map[var.idx]
    else
        error("A StructuredModel cannot have expression using variables of a StructuredModel different from itself or its parent")
    end
end
function parametrized_function(aff::JuMP.GenericAffExpr{C, StructuredVariableRef},
                               structured_model::StructuredModel,
                               model::JuMP.Model,
                               variable_map::Dict{Int, JuMP.VariableRef},
                               parameter_map::Dict{Int, ParameterJuMP.ParameterRef}) where C
    param_aff = ParameterJuMP.PAE{C}(JuMP.GenericAffExpr{C, JuMP.VariableRef}(aff.constant),
                                     JuMP.GenericAffExpr{C, ParameterJuMP.ParameterRef}(zero(C)))
    for (coef, var) in JuMP.linear_terms(aff)
        JuMP.add_to_expression!(param_aff,
                                coef,
                                parametrized_function(var, structured_model, model, variable_map, parameter_map))
    end
    # The inferred return value is an union of two concrete types but
    # small unions are efficiently handled in Julia v0.7
    if JuMP.iszero(param_aff.p)
        return param_aff.v
    else
        return param_aff
    end
end
function parametrized_function(quad::JuMP.GenericQuadExpr{C, StructuredVariableRef},
                               structured_model::StructuredModel,
                               model::JuMP.Model,
                               variable_map::Dict{Int, JuMP.VariableRef},
                               parameter_map::Dict{Int, ParameterJuMP.ParameterRef}) where C
    param_aff = parametrized_function(quad.aff, structured_model, model, variable_map, parameter_map)
    if param_aff isa ParameterJuMP.PAE
        error("parametrized quadratic functions are not supported yet")
        #quadv = JuMP.GenericQuadExpr{C, JuMP.VariableRef}(param_aff.v)
        #quadp = JuMP.GenericQuadExpr{C, ParameterJuMP.ParameterRef}(param_aff.p)
    else
        quadv = JuMP.GenericQuadExpr{C, JuMP.VariableRef}(param_aff)
        quadp = JuMP.GenericQuadExpr{C, ParameterJuMP.ParameterRef}(zero(C))
    end
    for (coef, var1, var2) in JuMP.quadterms(quad)
        if var1.model === structured_model && var2.model == structured_model
            JuMP.add_to_expression!(quadv, coef, variable_map[var1.idx], variable_map[var2.idx])
        else
            error("parametrized quadratic functions are not supported yet")
        end
    end
    if JuMP.iszero(quadp)
        return quadv
    else
        error("parametrized quadratic functions are not supported yet")
    end
end
parametrized_function(funcs::Vector, args...) = map(f -> parametrized_function(f, args...), funcs)

struct ParametrizedModel
    structured_model::StructuredModel
    model::JuMP.Model
    # Map between index of structured variable in `structured_model` and
    # the corresponding variable in `model`.
    variable_map::Dict{Int, JuMP.VariableRef}
    # Map between index of structured variable in `structured_model.parent`
    # and the corresponding parameter in `model`.
    parameter_map::Dict{Int, ParameterJuMP.ParameterRef}
    # Cost of children
    θ::Dict{Int, JuMP.VariableRef}
end
function ParametrizedModel(structured_model::StructuredModel, args...; kwargs...)
    if structured_model.parent === nothing
        model = Model(args...; kwargs...)
    else
        model = ModelWithParams(args...; kwargs...)
    end
    variable_map = Dict{Int, JuMP.VariableRef}()
    for (index, var) in structured_model.variables
        name = structured_model.varnames[index]
        variable_map[index] = JuMP.add_variable(model, var, name)
    end
    parameter_map = Dict{Int, ParameterJuMP.ParameterRef}()
    for (index, con) in structured_model.constraints
        name = structured_model.connames[index]
        param_fun = parametrized_function(con.func, structured_model, model,
                                          variable_map, parameter_map)
        param_con = JuMP.build_constraint(error, param_fun, con.set)
        JuMP.add_constraint(model, param_con, name)
    end
    objective_function = parametrized_function(structured_model.objective_function,
                                               structured_model, model,
                                               variable_map, parameter_map)
    if structured_model.objective_sense == MOI.MAX_SENSE
        objective_function = -objective_function
    end
    θ = Dict{Int, VariableRef}()
    for (id, proba) in getprobability(structured_model)
        θid = @variable(model, lower_bound = 0.0)
        θ[id] = θid
        JuMP.add_to_expression!(objective_function, proba, θid)
    end
    JuMP.set_objective(model, structured_model.objective_sense, objective_function)
    ParametrizedModel(structured_model, model, variable_map, parameter_map, θ)
end

include("Benders_pmap.jl")

# The optimizers are function returning an empty optimizer,
# this will be replaced by factories once this is implemented in JuMP
function BendersBridge(structured_model::StructuredModel,
                       master_optimizer::JuMP.OptimizerFactory,
                       sub_optimizer::JuMP.OptimizerFactory)
    master_model = ParametrizedModel(structured_model, master_optimizer)
    children = getchildren(structured_model)
    sub_models = Dict{Int, ParametrizedModel}()
    for (id, child) in getchildren(structured_model)
        sub_models[id] = ParametrizedModel(child, sub_optimizer)
    end
    return Benders_pmap(master_model, sub_models)
end


text:
## Spectral clustering on artists
using DelimitedFiles
using SparseArrays
names = readlines("../2-unit-1-demos/artistsim.names")
data = readdlm("../2-unit-1-demos/artistsim.smat")
xy = readdlm("artistsim.xy")
A = sparse(Int.(data[2:end,1]).+1,
           Int.(data[2:end,2]).+1,
           Int.(data[2:end,3]),
           Int(data[1,1]),Int(data[1,1]))
n = size(A,1)
##
function simplepagerank(A,α,v)
  @assert(0 ≤ α < 1, "needs probably α")
  @assert(all(vi -> vi ≥ 0, v), "needs non-negative v")
  v = v ./ sum(v) # we can normalize for them.
  d = vec(sum(A,dims=2)) # compute the degrees
  x = copy(v) # start of with v
  nsteps = 2*ceil(Int,log(eps(1.0))/log(α)) # upper bound on steps
  for i=1:nsteps
    x = α*(A'*(x./d)) .+ (1-α).*v
  end
  return x./sum(x) # renormalize to probability
end
##
using MatrixNetworks
pr = MatrixNetworks.pagerank(A,0.85)
##
x = simplepagerank(A,0.85,ones(n)./n)
##
using LinearAlgebra
@show norm(x-pr)

##
@time pr = MatrixNetworks.pagerank(A,0.85)
@time x = simplepagerank(A,0.85,ones(n)./n)
@show norm(x-pr)
## Let's watch PageRank converge
using Plots
function pagerank_converence_animation(A,xy,α,nsteps)
  n = size(A,1)
  v = ones(n)/n
  @assert(0 ≤ α < 1, "needs probably α")
  @assert(all(vi -> vi ≥ 0, v), "needs non-negative v")
  v = v ./ sum(v) # we can normalize for them.
  d = vec(sum(A,dims=2)) # compute the degrees
  x = copy(v) # start of with v

  @assert(nsteps <= size(A,1))
  anim = @animate for i=1:nsteps
    x = α*(A'*(x./d)) .+ (1-α).*v
    f2 = copy(x)
    f2 = f2/sum(f2)
    p = sortperm(f2) # this gets things to draw in the right z-order
    scatter(xy[p,1],xy[p,2], marker_z=-abs.(log10.(f2[p])).^0.5,colorbar=false,alpha=0.5,
      color=:magma,
      size=(1200,1200),
      markerstrokecolor=nothing,markerstrokewidth=0,label="",framestyle=:none,background=:black)
  end
  return anim
end
anim = pagerank_converence_animation(A,xy,0.85,25)
gif(anim,"pagerank-convergence-artistsim.gif",fps=5)

## Change between iterates
function pagerank_converence_plot(A,xy,α,nsteps)
  n = size(A,1)
  v = ones(n)/n
  @assert(0 ≤ α < 1, "needs probably α")
  @assert(all(vi -> vi ≥ 0, v), "needs non-negative v")
  v = v ./ sum(v) # we can normalize for them.
  d = vec(sum(A,dims=2)) # compute the degrees
  x = copy(v) # start of with v
  y = copy(v)
  hist = zeros(0)
  for i=1:nsteps
    x = α*(A'*(x./d)) .+ (1-α).*v
    push!(hist, norm(y - x/sum(x),1))
    y = x/sum(x)
  end
  return hist
end
hist = pagerank_converence_plot(A,xy,0.85,100)
##
pyplot()
theme(:dark)
plot(hist, yscale=:log10, legend=false, size=(300,300))
xlabel!("Iteration")
ylabel!("1-norm change in normalized vector")
savefig("pagerank-convergence-artistsim-norm.pdf")
##
x85 = simplepagerank(A,0.85,ones(n)./n)
x99 = simplepagerank(A,0.99,ones(n)./n)
x5 = simplepagerank(A,0.5,ones(n)./n)
##

## Show a few different vectors
gr()
p=sortperm(x99)
scatter(xy[p,1],xy[p,2],marker_z=-abs.(log10.(x99[p])).^0.5,markerstrokewidth=0,
color=:magma)
##
gr()
p=sortperm(x5)
scatter(xy[p,1],xy[p,2],marker_z=-abs.(log10.(x5[p])).^0.5,markerstrokewidth=0,
color=:magma)

## Now show the entire graph
using GraphRecipes
using LinearAlgebra
using Plots

ei, ej = findnz(A)
graphplot(ei, ej, x =xy[:,1], y=xy[:,2],
  markercolor=:black, markerstrokecolor=:white,
  size=(1200,1200),dpi=300,
  markersize=0, linecolor=1, linealpha=0.01, linewidth=0.5,
  markeralpha=0.2,colorbar=false,
  axis_buffer=0.02, background=nothing)
p=sortperm(x5)
scatter!(xy[p,1],xy[p,2],marker_z=-abs.(log10.(x5[p])).^0.5,markerstrokewidth=0)
savefig("artistsim-pagerank-50.png")

##
ei, ej = findnz(A)
graphplot(ei, ej, x =xy[:,1], y=xy[:,2],
  markercolor=:black, markerstrokecolor=:white,
  size=(1200,1200),dpi=300,
  markersize=0, linecolor=1, linealpha=0.01, linewidth=0.5,
  markeralpha=0.2,colorbar=false,
  axis_buffer=0.02, background=nothing)
p=sortperm(x99)
scatter!(xy[p,1],xy[p,2],marker_z=-abs.(log10.(x99[p])).^0.5,markerstrokewidth=0)
savefig("artistsim-pagerank-99.png")


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct Bus
        zone::Int
        name::String
        status::Int
        area::Int
        vm_ub::Float64
        vm_lb::Float64
        index::Int
        va::Float64
        vm::Float64
        base_kv::Float64
        type::Int
    end

# Arguments
- `zone::Int`: assigned control zone
- `name::String`: a flexible name for components, non required to be unique
- `status::Int`: a 0/1 value indicating if the component should be omitted or not (0 => omitted)
- `area::Int`: assigned control area
- `vm_ub::Float64`: an upper limit on voltage magnitude
- `vm_lb::Float64`: a lower limit on voltage magnitude
- `index::Int`: a unique identifier for components
- `va::Float64`: voltage angle
- `vm::Float64`: voltage magnitude
- `base_kv::Float64`: base voltage
- `type::Int`: bus type for power flow (1=PQ, 2=PV, 3=Slack)
"""
mutable struct Bus
    "assigned control zone"
    zone::Int
    "a flexible name for components, non required to be unique"
    name::String
    "a 0/1 value indicating if the component should be omitted or not (0 => omitted)"
    status::Int
    "assigned control area"
    area::Int
    "an upper limit on voltage magnitude"
    vm_ub::Float64
    "a lower limit on voltage magnitude"
    vm_lb::Float64
    "a unique identifier for components"
    index::Int
    "voltage angle"
    va::Float64
    "voltage magnitude"
    vm::Float64
    "base voltage"
    base_kv::Float64
    "bus type for power flow (1=PQ, 2=PV, 3=Slack)"
    type::Int
end

function Bus(; zone, name, status, area, vm_ub, vm_lb, index, va, vm, base_kv, type, )
    Bus(zone, name, status, area, vm_ub, vm_lb, index, va, vm, base_kv, type, )
end

# Constructor for demo purposes; non-functional.
function Bus(::Nothing)
    Bus(;
        zone=1,
        name="1",
        status=1,
        area=1,
        vm_ub=1.2,
        vm_lb=1.2,
        index=1,
        va=1.3,
        vm=1.3,
        base_kv=1.2,
        type=1,
    )
end

"""Get Bus zone."""
get_zone(value::Bus) = value.zone
"""Get Bus name."""
get_name(value::Bus) = value.name
"""Get Bus status."""
get_status(value::Bus) = value.status
"""Get Bus area."""
get_area(value::Bus) = value.area
"""Get Bus vm_ub."""
get_vm_ub(value::Bus) = value.vm_ub
"""Get Bus vm_lb."""
get_vm_lb(value::Bus) = value.vm_lb
"""Get Bus index."""
get_index(value::Bus) = value.index
"""Get Bus va."""
get_va(value::Bus) = value.va
"""Get Bus vm."""
get_vm(value::Bus) = value.vm
"""Get Bus base_kv."""
get_base_kv(value::Bus) = value.base_kv
"""Get Bus type."""
get_type(value::Bus) = value.type

"""Set Bus zone."""
set_zone!(value::Bus, val::Int) = value.zone = val
"""Set Bus name."""
set_name!(value::Bus, val::String) = value.name = val
"""Set Bus status."""
set_status!(value::Bus, val::Int) = value.status = val
"""Set Bus area."""
set_area!(value::Bus, val::Int) = value.area = val
"""Set Bus vm_ub."""
set_vm_ub!(value::Bus, val::Float64) = value.vm_ub = val
"""Set Bus vm_lb."""
set_vm_lb!(value::Bus, val::Float64) = value.vm_lb = val
"""Set Bus index."""
set_index!(value::Bus, val::Int) = value.index = val
"""Set Bus va."""
set_va!(value::Bus, val::Float64) = value.va = val
"""Set Bus vm."""
set_vm!(value::Bus, val::Float64) = value.vm = val
"""Set Bus base_kv."""
set_base_kv!(value::Bus, val::Float64) = value.base_kv = val
"""Set Bus type."""
set_type!(value::Bus, val::Int) = value.type = val

max_stars_repo_name: NREL-SIIP/CommonTransmissionModel

max_issues_repo_name: NREL-SIIP/CommonTransmissionModel

max_forks_repo_name: NREL-SIIP/CommonTransmissionModel


text:
plot(sin, 0, π, label="sine", lw=0.1)
scatter!(cos, 0, π, length=25, label="cosine")
xlim(0, π)
legend(fontsize=12, font="texcmtt", nobox=true, position="bottom-left")


text:
ss_cu_model(t) = sample(t, "a", Normal())
ss_du_model(t) = sample(t, "a", Poisson(3.0))
ss_cm_model(t) = sample(t, "a", MvNormal(3, 1.0))

function ms_model(t)
    log_rate = sample(t, "log_rate", Normal(2.0, 1.0))
    sample(t, "data", Poisson(exp(log_rate)))
end

@testset "single site" begin
    @testset "unrestricted continuous kernel" begin
        t0 = trace()
        ss_cu_model(t0)
        k = make_kernel(ss_cu_model, t0["a"])

        t1 = trace()
        k(t0, t1)
        @info "Original trace: $t0"
        @info "New trace: $t1"
        @test mean(t1["a"].dist) == t0["a"].value
    end

    @testset "discrete univariate kernel" begin
        t0 = trace()
        ss_du_model(t0)
        k = make_kernel(ss_du_model, t0["a"])

        t1 = trace()
        k(t0, t1)
        @info "Original trace: $t0"
        @info "New trace: $t1"
        @test length(t1["a"].dist.support) == 3

        cond = condition(ss_du_model, Dict("a" => 0))
        t2 = trace()
        cond(t2)
        t3 = trace()
        k(t2, t3)
        @info "Original trace: $t2"
        @info "New trace: $t3"
        @test t3["a"].dist.support == [0, 1]
    end

    @testset "continuous multivariate kernel" begin
        t0 = trace()
        ss_cm_model(t0)
        k = make_kernel(ss_cm_model, t0["a"])

        t1 = trace()
        k(t0, t1)
        @info "Original trace: $t0"
        @info "New trace: $t1"
        @test mean(t1["a"].dist) == t0["a"].value
    end

    @testset "generate kernels" begin
        kernels = make_kernels(ms_model)
        t0 = trace()
        ms_model(t0)
        t1 = trace()
        map(k -> k(t0, t1), kernels)
        @info "Original trace: $t0"
        @info "New trace: $t1"
    end

end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
;ls D:\\texlive\\2021\\bin\\win32\\lualatex.exe

# %%
using PGFPlotsX

p1 = range(0, 360; length = 100)' # parameter 1
p2 = range(0, 360; length = 70)   # parameter 2
w1 = @. sind(3*(p1 + 2 * p2)) + 1.25 # wave 1
w2 = @. 6 + w1 * cosd(p1)            # wave 2
x = vec(@. w2 * cosd(p2))
y = vec(@. w2 * sind(p2))
z = vec(@. w1 * sind(p1))
logo_colors = [(77,100,174), (57,151,79), (255,255,255), (146,89,163), (202,60,50)]

@pgf Axis(
    {
        axis_equal,
        axis_lines = "none"
    },
    Plot3(
        {
            surf,
            z_buffer = "sort",
            colormap = "{Julia}{$(join(map(c -> "rgb255=$c", logo_colors), ", "))}",
            "mesh/rows" = length(p1)
        },
        Table(x, y, z)))

# %%

max_stars_repo_path: 0015/PGFPlotsX example.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0015/PGFPlotsX example.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0015/PGFPlotsX example.jl

max_forks_repo_name: genkuroki/public


text:
#CircuitAnalysis: Base tools

#==Main data structures
===============================================================================#

#Object used to provide type information to something:
struct DataTag{ID, T}
	v::T
end
#=TODO: Figure out how to broadcast operations on DataTag (ex: make 5 .* ycap
work)... but not ycap .* ycap.
=#

#==Type aliases
===============================================================================#

#Tags:
const TImpedance = DataTag{:Z}
const TAdmittance = DataTag{:Y}
const TInductance = DataTag{:L}
const TCapacitance = DataTag{:C}

#Base.size(x::DataTag) = size(x.v)
#Base.length(x::DataTag) = length(x.v)
#Base.eltype(x::DataTag{ID, T}) where {ID, T} = T

#==Useful validations/assertions
===============================================================================#

#==Constructor interfaces
===============================================================================#

(::Type{DataTag{ID}})(v::T) where {ID, T} = DataTag{ID, T}(v)

#==Helper functions
===============================================================================#
#multMD: disambiguate multiplication to use
multMD(a, b) = a*b #DataMD values don't require .*
multMD(a::AbstractArray, b::AbstractArray) = a .* b
divMD(a, b) = a/b #DataMD values don't require .*
divMD(a::AbstractArray, b::AbstractArray) = a ./ b

#==Operations
===============================================================================#
#Mostly for unit testing:
Base.:(≈)(a::DataTag, b::DataTag; kwargs...) = false
Base.:(≈)(a::DataTag{ID}, b::DataTag{ID}; kwargs...) where ID =
	Base.:(≈)(a.v, b.v; kwargs...)

Base.:(==)(a::DataTag, b::DataTag) = false
Base.:(==)(a::DataTag{ID}, b::DataTag{ID}) where ID = (a.v == b.v)

#_dotop(x)=Symbol(".$x")
_operators1 = [:+, :-] #Unary operator on a tagged object
_operators2 = [:+, :-, :*, :/]

for op in _operators1; @eval begin #CODEGEN-------------------------------------

Base.$op(d1::DataTag{ID}) where ID = DataTag{ID}(Base.$op(d1.v))

end; end #CODEGEN---------------------------------------------------------------

#=

for op in _operators2; @eval begin #CODEGEN--------------------------------

Base.$op(d1::DataTag{ID}, d2::Number) where ID =
	DataTag{ID}(Base.$op(d1.v, d2))

Base.$op(d1::Number, d2::DataTag{ID}) where ID =
	DataTag{ID}(Base.$op(d1, d2.v))

Base.$op(d1::DataTag{ID}, d2::DataTag{ID}) where ID =
	DataTag{ID}(Base.$op(d1.v, d2.v))

end; end #CODEGEN---------------------------------------------------------------

=#
#==Main functions
===============================================================================#
impedance(v) = TImpedance(v)
admittance(v) = TAdmittance(v)
capacitance(v) = TCapacitance(v)
inductance(v) = TInductance(v)

impedance(y::TAdmittance) = impedance(divMD(1, y.v))
admittance(z::TImpedance) = admittance(divMD(1, z.v))

#Capacitor values:
_Y(c::TCapacitance, f) = admittance(multMD((2*pi*c.v)im, f))
_Y(c::TCapacitance, ::Nothing) = throw(ArgumentError("Missing kwarg :f"))
admittance(c::TCapacitance; f = nothing) = _Y(c, f)
impedance(c::TCapacitance, args...; kwargs...) = impedance(admittance(c, args...; kwargs...))
#TODO: create admittance(:L, value, f=x)??

#Inductor values:
_Z(l::TInductance, f) = admittance(multMD((2*pi*c.v)im, f))
_Z(l::TInductance, ::Nothing) = throw(ArgumentError("Missing kwarg :f"))
impedance(l::TInductance; f = nothing) = _Z(l, f)
admittance(l::TInductance, args...; kwargs...) = admittance(impedance(l, args...; kwargs...))

#==Helper functions
===============================================================================#

#==Traits (utility functions for instances and types)
===============================================================================#

#Last line


text:
export ARS111, ARS121, ARS232, ARS343

struct ARSAlgorithm{name,L} <: DistributedODEAlgorithm
    linsolve::L
end
ARSAlgorithm{name}(;linsolve) where {name} = ARSAlgorithm{name,typeof(linsolve)}(linsolve)

struct ARSTableau{RT}
    a::Matrix{RT}
    ahat::Matrix{RT}
    c::Vector{RT}
    chat::Vector{RT}
    Q0::Vector{RT}
    Q::Matrix{RT}
    Qhat::Matrix{RT}
    γ::RT
end
function ARSTableau(a::AbstractMatrix{RT}, ahat::AbstractMatrix{RT}) where {RT}
    #=
    utilde[i] = u0 +
    dt * sum(j -> ahat[i+1,j] * g(U[j-1]), 1:i-1) +
    dt * sum(j -> a[i,j]      * f(U[j]), 1:i-1)

    Uhat[i] = utilde[i] + dt * ahat[i+1,i] g(U[i-1])
    U[i] = Uhat[i] + dt * a[i,i] f(U[i])
        ≈ Uhat[i] + dt * a[i,i] (I - dt * a[i,i] * J) \ f(U[i])

    we don't store g or f, or utilde

    g(U[i-1]) = (Uhat[i] - utilde[i]) / (dt * ahat[i+1,i])
    f(U[i]) = (U[i] - Uhat[i]) / (dt * a[i,i])

    utilde[i] = u0 +
    sum(j -> (ahat[i+1,j]/ahat[j+1,j]) * (Uhat[j] - utilde[j]), 1:i-1) +
    sum(j -> (a[i,j]/a[j,j]) (U[j] - Uhat[j]), 1:i-1)

    write this as

    utilde[i] = Q0[i] +
    sum(k -> Qhat[i,k] * Uhat[k], 1:i-1) +
    sum(k -> Q[i,k] * U[k], 1:i-1)

    =#

    N = size(a,2)

    c = vec(sum(a, dims=2))
    chat = vec(sum(ahat, dims=2))

    Q = zeros(RT, N+1,N)
    Q0 = zeros(RT,N+1)
    Qhat = zeros(RT, N+1,N)

    for i = 1:N+1
      Q0[i] = 1 - sum(j -> ahat[i+1,j] / ahat[j+1,j] * Q0[j], 1:i-1;init=0.0)
      for k = 1:i-1
        Q[i, k] = a[i,k]/a[k,k] - sum(j -> ahat[i+1,j] / ahat[j+1,j] * Q[j,k], 1:i-1)
        Qhat[i,k] = ahat[i+1,k]/ahat[k+1,k] - a[i,k]/a[k,k] - sum(j -> ahat[i+1,j] / ahat[j+1,j] * Qhat[j,k], 1:i-1)
      end
    end

    for i = 1:N+1
      @assert Q0[i] + sum(Q[i,:]) + sum(Qhat[i,:]) ≈ 1
    end

    γ = a[1,1]
    for i = 2:N
      @assert a[i,i] == γ
    end
    ARSTableau(a, ahat, c, chat, Q0, Q, Qhat, γ)
end

"""
    ARS111

The Forward-Backward (1,1,1) implicit-explicit (IMEX) Runge-Kutta scheme of
[ARS1997](@cite), section 2.1.

This is equivalent to the `OrdinaryDiffEq.IMEXEuler` algorithm.
"""
const ARS111 = ARSAlgorithm{:ARS111}
function tableau(::ARS111, RT)
    # https://github.com/SciML/OrdinaryDiffEq.jl/issues/1590
    a = RT[1;
           1;;]
    ahat = RT[0 0;
              1 0;
              1 0]
    return ARSTableau(a, ahat)
end

"""
    ARS121

The Forward-Backward (1,2,1) implicit-explicit (IMEX) Runge-Kutta scheme of
[ARS1997](@cite), section 2.2.

This is equivalent to the `OrdinaryDiffEq.IMEXEulerARK` algorithm.
"""
const ARS121 = ARSAlgorithm{:ARS121}
function tableau(::ARS121, RT)
    a = RT[1;
           1;;]
    ahat = RT[0 0;
              1 0;
              0 1]
    return ARSTableau(a, ahat)
end

"""
    ARS232

The Forward-Backward (2,3,2) implicit-explicit (IMEX) Runge-Kutta scheme of
[ARS1997](@cite), section 2.5.
"""
const ARS232 = ARSAlgorithm{:ARS232}
function tableau(::ARS232, RT)
    γ = (2 - sqrt(2))/2
    δ = -2*sqrt(2)/3
    # implicit
    a = RT[γ 0;
        1-γ γ;
        1-γ γ];
    # explicit
    ahat = RT[ 0 0 0;
            γ 0 0;
            δ 1-δ 0;
            0 1-γ γ];
    return ARSTableau(a, ahat)
end

"""
    ARS343

The L-stable, third-order (3,4,3) implicit-explicit (IMEX) Runge-Kutta scheme of
[ARS1997](@cite), section 2.7.
"""
const ARS343 = ARSAlgorithm{:ARS343}
function tableau(::ARS343, RT)
    N = 3
    γ = 0.4358665215084590
    b1 = -3/2 * γ^2 + 4 * γ - 1/4
    b2 =  3/2 * γ^2 - 5 * γ + 5/4
    # implicit tableau
    a = RT[γ 0 0;
        (1-γ)/2 γ 0;
        b1 b2 γ;
        b1 b2 γ]

    dA42 = 0.5529291480359398;
    dA43 = 0.5529291480359398;
    dA31 =(
        (1.0 - 4.5 * γ + 1.5 * γ * γ) * dA42
          + (2.75 - 10.5 * γ + 3.75 * γ * γ) * dA43
          - 3.5 + 13 * γ - 4.5 * γ * γ)
    dA32 = (
        (-1.0 + 4.5 * γ - 1.5 * γ * γ) * dA42
          + (-2.75 + 10.5 * γ - 3.75 * γ * γ) * dA43
          + 4.0 - 12.5 * γ + 4.5 * γ * γ)
    dA41 = 1.0 - dA42 - dA43;
    # explicit tableau
    ahat = RT[
      0 0 0 0;
      γ 0 0 0;
      dA31 dA32 0 0;
      dA41 dA42 dA43 0;
      0 b1 b2 γ
    ]
    ARSTableau(a, ahat)
end

struct ARSCache{Nstages, RT, A}
    tableau::ARSTableau{RT}
    U::NTuple{Nstages, A}
    Uhat::NTuple{Nstages, A}
    idu::A
    W
    linsolve!
end

function cache(
    prob::DiffEqBase.AbstractODEProblem,
    alg::ARSAlgorithm; kwargs...)

    tab = tableau(alg, eltype(prob.u0))
    Nstages = length(tab.c) - 1
    U = ntuple(i -> similar(prob.u0), Nstages)
    Uhat = ntuple(i -> similar(prob.u0), Nstages)
    idu = similar(prob.u0)
    W = prob.f.f1.jac_prototype
    linsolve! = alg.linsolve(Val{:init}, W, prob.u0; kwargs...)

    return ARSCache(tab, U, Uhat, idu, W, linsolve!)
end

function step_u!(int, cache::ARSCache{Nstages}) where {Nstages}

    f = int.prob.f
    f1! = f.f1
    f2! = f.f2
    Wfact! = f1!.Wfact

    u = int.u
    p = int.p
    t = int.t
    dt = int.dt

    tab = cache.tableau
    W = cache.W
    U = cache.U
    Uhat = cache.Uhat
    idu = cache.idu
    linsolve! = cache.linsolve!

    # Update W
    # Wfact!(W, u, p, dt*tab.γ, t)

    # implicit eqn:
    #   ux = u + dt * f(ux, p, t)
    # Newton iteration:
    #   ux <- ux + (I - dt J) \ (u + dt f(ux, p, t) - ux)
    # initial iteration
    #   ux <- u + dt (I - dt J) \ f(u, p, t)

    function implicit_step!(ux, u, p, t, dt)
        Wfact!(W, u, p, dt, t)
        # currently this does just a single Newton iteration
        # and assumes we compute f(u, p, t) directly
        # need to figure out how we would do multiple iterations
        # du = f(u, p, t)
        f1!(ux, u, p, t)
        # @show ux
        # solve for idu = (I - dt J) \ du
        linsolve!(idu, W, ux)
        # @show W
        # @show idu
        @. ux = u + dt * idu
        # @show ux
    end

    #### stage 1
    # explicit
    Uhat[1] .= u # utilde[i],  Q0[1] == 1
    f2!(Uhat[1], u, p, t+dt*tab.chat[1], dt*tab.ahat[2,1])

    # implicit
    implicit_step!(U[1], Uhat[1], p, t+dt*tab.c[1], dt*tab.a[1,1])
    if Nstages == 1
        u .= tab.Q0[2] .* u .+
            tab.Qhat[2,1] .* Uhat[1] .+ tab.Q[2,1] .* U[1] # utilde[2]
        f2!(u, U[1], p, t+dt*tab.chat[2], dt*tab.ahat[3,2])
        return
    end

    #### stage 2
    Uhat[2] .= tab.Q0[2] .* u .+
            tab.Qhat[2,1] .* Uhat[1] .+ tab.Q[2,1] .* U[1] # utilde[2]
    f2!(Uhat[2], U[1], p, t+dt*tab.chat[2], dt*tab.ahat[3,2])

    implicit_step!(U[2], Uhat[2], p, t+dt*tab.c[2], dt*tab.a[2,2])

    if Nstages == 2
        u .= tab.Q0[3] .* u .+
            tab.Qhat[3,1] .* Uhat[1] .+ tab.Q[3,1] .* U[1] .+
            tab.Qhat[3,2] .* Uhat[2] .+ tab.Q[3,2] .* U[2] # utilde[3]
        f2!(u, U[2], p, t+dt*tab.chat[3], dt*tab.ahat[4,3])
        return
    end

    #### stage 3
    Uhat[3] .= tab.Q0[3] .* u .+
            tab.Qhat[3,1] .* Uhat[1] .+ tab.Q[3,1] .* U[1] .+
            tab.Qhat[3,2] .* Uhat[2] .+ tab.Q[3,2] .* U[2] # utilde[3]
    f2!(Uhat[3], U[2], p, t+dt*tab.chat[3], dt*tab.ahat[4,3])
    # @show Uhat[3] t+dt*tab.chat[3]

    implicit_step!(U[3], Uhat[3], p, t+dt*tab.c[3], dt*tab.a[3,3])
    # @show U[3] t+dt*tab.c[3]

    ### final update
    u .= tab.Q0[4] .* u .+
    tab.Qhat[4,1] .* Uhat[1] .+ tab.Q[4,1] .* U[1] .+
    tab.Qhat[4,2] .* Uhat[2] .+ tab.Q[4,2] .* U[2] .+
    tab.Qhat[4,3] .* Uhat[3] .+ tab.Q[4,3] .* U[3]

    # @show u
    f2!(u, U[3], p, t+dt*tab.chat[4], dt*tab.ahat[5,4])
    # @show u t+dt*tab.chat[4]
    return
end


text:
module TestExamplesP4estMesh2D

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "p4est_2d_dgsem")

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
isdir(outdir) && rm(outdir, recursive=true)

@testset "P4estMesh2D" begin
  @trixi_testset "elixir_advection_basic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"),
      # Expected errors are exactly the same as with StructuredMesh!
      l2   = [8.311947673061856e-6],
      linf = [6.627000273229378e-5])
  end

  @trixi_testset "elixir_advection_nonconforming.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_nonconforming.jl"),
      l2   = [2.7905480848832338e-5],
      linf = [0.00022847020768290704])
  end

  @trixi_testset "elixir_advection_nonconforming_unstructured.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_nonconforming_unstructured.jl"),
      l2   = [0.0026958321660563362],
      linf = [0.04122573088346193])
  end

  @trixi_testset "elixir_advection_amr_solution_independent.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_solution_independent.jl"),
      # Expected errors are exactly the same as with StructuredMesh!
      l2   = [4.949660644033807e-5],
      linf = [0.0004867846262313763])
  end

  @trixi_testset "elixir_advection_amr_unstructured_flag.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_unstructured_flag.jl"),
      l2   = [0.0012766060609964525],
      linf = [0.01750280631586159])
  end

  @trixi_testset "elixir_advection_restart.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"),
      l2   = [4.219208035626337e-6],
      linf = [3.4384344042126536e-5])
  end

  @trixi_testset "elixir_euler_source_terms_nonperiodic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_nonperiodic.jl"),
      l2   = [2.2594405120861626e-6, 2.3188881560096627e-6, 2.318888156064081e-6, 6.332786324236605e-6],
      linf = [1.4987382633613322e-5, 1.9182011925522602e-5, 1.9182011924634423e-5, 6.0526717124531615e-5])
  end

  @trixi_testset "elixir_euler_source_terms_nonperiodic_unstructured.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_nonperiodic_unstructured.jl"),
      l2   = [0.005238881525717585, 0.0043246899191607775, 0.0043246899191606925, 0.00986166157942196],
      linf = [0.052183952487556695, 0.05393708345945791, 0.05393708345946191, 0.09119199890635965])
  end

  @trixi_testset "elixir_euler_free_stream.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream.jl"),
      l2   = [2.063350241405049e-15, 1.8571016296925367e-14, 3.1769447886391905e-14, 1.4104095258528071e-14],
      linf = [1.9539925233402755e-14, 2e-12, 4.8e-12, 4e-12],
      atol = 2.0e-12, # required to make CI tests pass on macOS
    )
  end

  @trixi_testset "elixir_eulergravity_convergence.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulergravity_convergence.jl"),
      l2   = [0.00024871265138964204, 0.0003370077102132591, 0.0003370077102131964, 0.0007231525513793697],
      linf = [0.0015813032944647087, 0.0020494288423820173, 0.0020494288423824614, 0.004793821195083758],
      tspan = (0.0, 0.1))
  end
end

# Clean up afterwards: delete Trixi output directory
@test_nowarn rm(outdir, recursive=true)

end # module


text:
module RData

using DataFrames, CategoricalArrays, FileIO, TimeZones, Unicode

export
    sexp2julia,
    DictoVec,
    RDataException,
    UnsupportedROBJ,
    CodecMissingError,
    load # reexport FileIO.load()

include("config.jl")
include("sxtypes.jl")

"""
Abstract RDA format IO stream wrapper.
"""
abstract type RDAIO end

include("errors.jl")

include("io/XDRIO.jl")
include("io/ASCIIIO.jl")
include("io/NativeIO.jl")
include("io/utils.jl")

include("DictoVec.jl")
include("identifier.jl")
include("convert.jl")

include("context.jl")
include("readers.jl")

include("decompression.jl")

function __init__()
    define_optional_decompressor_streams()
end

##############################################################################
##
## FileIO integration.
## supported `kwoptions`:
## convert::Bool (true by default) for converting R objects into Julia equivalents,
##               otherwise load() returns R internal representation (ROBJ-derived objects)
## TODO option for disabling names checking (e.g. column names)
##
##############################################################################

function fileio_load(s::Stream{format"RData"}; kwoptions...)
    io = stream(s)
    @assert FileIO.detect_rdata(io)
    ctx = RDAContext(rdaio(io, chomp(readline(io))); kwoptions...)
    @assert ctx.fmtver == 2 || ctx.fmtver == 3  # supported format versions
    @debug "written by R version $(ctx.Rver)"
    @debug "minimal R version: $(ctx.Rmin)"

    convert2julia = get(ctx.kwdict, :convert, true)

    # top level read -- must be a paired list of objects
    # we read it here to be able to convert to julia objects inplace
    fl = readuint32(ctx.io)
    sxtype(fl) == LISTSXP || error("Top level R object is not a paired list")
    !hasattr(fl) || error("Top level R paired list should have no attributes")

    res = Dict{RString,Any}()
    while sxtype(fl) != NILVALUE_SXP
        hastag(fl) || error("Top level list element has no name")
        tag = readitem(ctx)
        obj_name = convert(RString, isa(tag, RSymbol) ? tag.displayname : "\0")
        obj = readitem(ctx)
        setindex!(res, (convert2julia ? sexp2julia(obj) : obj), obj_name)
        fl = readuint32(ctx.io)
        readattrs(ctx, fl)
    end

    return res
end

function fileio_load(s::Stream{format"RDataSingle"}; kwoptions...)
    io = stream(s)
    @assert FileIO.detect_rdata_single(io)
    ctx = RDAContext(rdaio(io, chomp(readline(io))); kwoptions...)
    @assert ctx.fmtver == 2 || ctx.fmtver == 3  # supported format versions
    convert2julia = get(ctx.kwdict, :convert, true)
    return convert2julia ? sexp2julia(readitem(ctx)) : readitem(ctx)
end

function fileio_load(f::Union{File{format"RData"}, File{format"RDataSingle"}};
                     kwoptions...)
    io = open(filename(f), "r")
    try
        io = decompress(io)
        return fileio_load(Stream(f, io); kwoptions...)
    catch
        rethrow()
    finally
        close(io)
    end
end

end # module


text:
### A Pluto.jl notebook ###
# v0.18.0

using Markdown
using InteractiveUtils

# ╔═╡ b7772def-6186-44a2-b0d8-bbf55855948a
import Pkg

# ╔═╡ a194ae49-7def-4020-b047-995f15ad6516
Pkg.activate(Base.current_project())

# ╔═╡ 4e023a0d-8860-4d3e-ad20-7f6f2026efad
using Omega, Distributions, UnicodePlots

# ╔═╡ 1d48fee0-8bc6-11ec-3c86-3f1be0c6595c
md" # Basic Tutorial

In this tutorial we will run through the basics of creating a model and conditioning it.

First load Omega and Distributions:"

# ╔═╡ 6ef3eaa6-87e0-4206-b1b4-17e9fbc820a0
md"If you tossed a coin and observed the sequqnce `HHHHH`, you would be a little suspicious, `HHHHHHHH` would make you very suspicious.
Elementary probability theory tells us that for a fair coin, `HHHHHHHH` is just a likely outcome as `HHTTHHTH`.  What gives?
We will use Omega to model this behaviour, and see how that belief about a coin changes after observing a number of tosses.

Model the coin as a bernoulli distribution.  The weight of a bernoulli determines the probability it comes up true (which represents heads). Use a [beta distribution](https://en.wikipedia.org/wiki/Beta_distribution) to represent our prior belief weight of the coin."

# ╔═╡ eb7d27e9-2784-48be-b96c-4810edbfd07c
weight = @~ Beta(2.0, 2.0)

# ╔═╡ 4c985b30-307e-444c-bebf-94bba60500fe
beta_samples = randsample(weight, 10000)

# ╔═╡ e6405695-7577-4654-8971-496aad5f7e82
md"Let's see what this distribution looks like using UnicodePlots.  If you don't have it installed already install with:"

# ╔═╡ eb64defe-bce0-46d1-aa4b-83ea6f7927ab
UnicodePlots.histogram(beta_samples)

# ╔═╡ 8c796c06-0413-4430-8bc1-652471a4bc49
md"""
The distribution is symmetric around 0.5 and has support over the the interval [0, 1].

So far we have not done anything we couldn't do with `Distributions.jl`.
A primary distinction between a package like `Distribution.jl`, is that `Omega.jl` allows you to __condition__ probability distributions.

Create a model representing four flips of the coin.
Since a coin can be heads or tales, the appropriate distribution is the [bernouli distribution](https://en.wikipedia.org/wiki/Bernoulli_distribution):"""

# ╔═╡ e92d4bdb-b0f5-40f9-ac72-ba6512e40e7a
nflips = 4

# ╔═╡ 4c80763b-36ef-4dd7-9fcc-77865d410dd5
typeof(~ Bernoulli.(weight))

# ╔═╡ 0f6b489e-37c9-4728-ac1f-e9fa8db41a45
randsample(@~ (Bernoulli.(weight)))

# ╔═╡ 5ca61b39-a4bc-42d9-8bad-8a234f223b70
coinflips = Mv(Bernoulli.(weight), 1:nflips)

# ╔═╡ b8c6c52c-1bed-4b3f-a970-a6c4b30a6457
md"""

Take note that `weight` is the random variable defined previously.
`bernoulli` takes a type as its secoond argument; `Bool` indicates the result will be a `Bool` rather than an `Int`.

`coinflips` is a normal Julia array of Random Variables (`RandVar`s).
For reasons we will elaborate in later sections, it will be useful to have an `Array`-valued `RandVar` (instead of an `Array` of `RandVar`).

One way to do this (there are several ways discuseed later), is to use the function `randarray`
"""

# ╔═╡ 6351d968-668b-4df6-9fb8-d2f2ce97de1d
md"`coinflips` is a random variable and hence we can sample from it with `randsample`"

# ╔═╡ b6d44330-6b40-4a42-bea6-d915e9e35945
randsample(coinflips)

# ╔═╡ 3d76d33f-4149-404c-b929-143562919d2c
md"""Now we can condition the model.
We want to find the conditional distribution over the weight of the coin given some observations.
"""

# ╔═╡ ffb7ae8b-3d3f-4098-a79a-52bd4c541790
observations = [true, true, true, false]

# ╔═╡ d1134d33-2dfa-47b0-90cc-28fa20f6bbdb
md"Create a predicate that tests whether simulating from the model matches the observed data:"

# ╔═╡ 3e1f5056-0c16-4cbd-af84-10f835bf656d
condition = coinflips .== observations

# ╔═╡ 2f6c2ed4-614a-4d6d-83bb-32d2a2232a8e
md"""
`condition` is a random variable; we can sample from it.  The function `==ᵣ` (and more generally functions subscripted with ᵣ) should be read as "a realization of coinflips == observations"

We can use `rand` to sample from the model conditioned on `condition` being true:"""

# ╔═╡ 370da9a6-6e6f-42dd-83af-ba02d629d6be
weight_samples = randsample(weight |ᶜ condition, 10; alg = RejectionSample)

# ╔═╡ 728d77db-64ed-45d1-8c69-85d83ff2d3b6
md"""
`weight_samples` is a set of `10` samples from the conditional (sometimes called posterior) distribution of `weight` condition on the fact that coinflips == observations.

In this case, `rand` takes
- A random variable we want to sample from
- A predicate (type `RandVar` which evaluates to a `Bool`) that we want to condition on, i.e. assert that it is true
- An inference algorithm.  Here we use rejection sampling.

Plot a histogram of the weights like before:
"""

# ╔═╡ 8c7684b3-9a1f-4e09-b8ba-038600150f27
UnicodePlots.histogram(weight_samples)

# ╔═╡ 671aae80-cb2d-4ad4-9964-a07f5f9cd0ef
md"""Observe that our belief about the weight has now changed.
We are more convinced the coin is biased towards heads (`true`)."""

# ╔═╡ Cell order:
# ╠═b7772def-6186-44a2-b0d8-bbf55855948a
# ╠═a194ae49-7def-4020-b047-995f15ad6516
# ╟─1d48fee0-8bc6-11ec-3c86-3f1be0c6595c
# ╠═4e023a0d-8860-4d3e-ad20-7f6f2026efad
# ╟─6ef3eaa6-87e0-4206-b1b4-17e9fbc820a0
# ╠═eb7d27e9-2784-48be-b96c-4810edbfd07c
# ╠═4c985b30-307e-444c-bebf-94bba60500fe
# ╟─e6405695-7577-4654-8971-496aad5f7e82
# ╠═eb64defe-bce0-46d1-aa4b-83ea6f7927ab
# ╟─8c796c06-0413-4430-8bc1-652471a4bc49
# ╠═e92d4bdb-b0f5-40f9-ac72-ba6512e40e7a
# ╠═4c80763b-36ef-4dd7-9fcc-77865d410dd5
# ╠═0f6b489e-37c9-4728-ac1f-e9fa8db41a45
# ╠═5ca61b39-a4bc-42d9-8bad-8a234f223b70
# ╟─b8c6c52c-1bed-4b3f-a970-a6c4b30a6457
# ╟─6351d968-668b-4df6-9fb8-d2f2ce97de1d
# ╠═b6d44330-6b40-4a42-bea6-d915e9e35945
# ╟─3d76d33f-4149-404c-b929-143562919d2c
# ╠═ffb7ae8b-3d3f-4098-a79a-52bd4c541790
# ╟─d1134d33-2dfa-47b0-90cc-28fa20f6bbdb
# ╠═3e1f5056-0c16-4cbd-af84-10f835bf656d
# ╟─2f6c2ed4-614a-4d6d-83bb-32d2a2232a8e
# ╠═370da9a6-6e6f-42dd-83af-ba02d629d6be
# ╟─728d77db-64ed-45d1-8c69-85d83ff2d3b6
# ╠═8c7684b3-9a1f-4e09-b8ba-038600150f27
# ╟─671aae80-cb2d-4ad4-9964-a07f5f9cd0ef

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
include("../src/interfaces/highs_julia_api.jl")

cc = [1.0, -2.0]
cl = [0.0, 0.0]
cu = [10.0, 10.0]
ru = [2.0, 1.0]
rl = [0.0, 0.0]
astart = [0, 2, 4]
aindex = [0, 1, 0, 1]
avalue = [1.0, 2.0, 1.0, 3.0]

model = HighsModel(cc, cl, cu, rl, ru, astart, aindex, avalue)

# solve model in one go
status, solution, basis, modelstatus = Highs_call(model)

# create highs object and work with it
highs = Highs_create()
Highs_loadModel(highs, model)
Highs_run(highs)
Highs_getSolution(highs)
Highs_getBasis(highs)
Highs_getObjectiveValue(highs)
Highs_getIterationCount(highs)
Highs_destroy(highs)

# read model from file
highs = Highs_create()
Highs_readFromFile(highs, "../check/instances/adlittle.mps")
Highs_run(highs)
Highs_getSolution(highs)
Highs_getBasis(highs)
Highs_getObjectiveValue(highs)
Highs_getIterationCount(highs)
Highs_destroy(highs)

# create model from scratch
highs = Highs_create()
Highs_addCols(highs, cc, cl, cu, [], [], [])
Highs_addRows(highs, rl, ru, astart, aindex, avalue)
Highs_run(highs)
Highs_getSolution(highs)
Highs_getBasis(highs)
Highs_getObjectiveValue(highs)
Highs_getIterationCount(highs)
Highs_destroy(highs)

max_stars_repo_name: ycollet/HiGHS

max_issues_repo_name: ycollet/HiGHS

max_forks_repo_name: ycollet/HiGHS


text:
using JLD2, JSON
include("../../src/SFFM.jl")
@load pwd()*"/dev/CMEParamsData/CMEParams1.jld2" tempDict
CMEParams = tempDict

erlangDParams = Dict()

function integrateD(evals,order)
    # evals is an integer specifying how many points to eval the function at
    # order is the order of the erlang distribution

    erlang = SFFM.MakeErlang(order)
    λ = -erlang.S[1,1]

    edges = range(0,1+20/order,length=evals+1) # points at which to evaluate the fn
    h = (4/order)/(evals)

    orbit_LHS = erlang.a'
    orbit_RHS = zeros(order)
    v_RHS = zeros(order)
    v_RHS[1] = 1
    v_LHS = ones(order)
    D = zeros(order,order)
    fct = factorial.(big.(0:order-1))
    poly(x) = convert.(Float64,((λ*x).^(0:order-1))./fct)
    for t in edges[2:end]
        vec = poly(t)
        orbit_RHS = vec./sum(vec)
        v_RHS = exp(-λ*t).*(cumsum(vec)[end:-1:1])
        orbit = (orbit_LHS+orbit_RHS)./2

        v = v_LHS - v_RHS

        Dᵢ = v*orbit'
        D += Dᵢ

        orbit_LHS = copy(orbit_RHS)
        v_LHS = copy(v_RHS)
    end
    return D
end

# numerical approximation of D
k = 100_000_000
T=[]
N=[]
for n in keys(CMEParams)
    evals = Int(ceil(k/(n/3)))# number of function evals
    display(n)
    display(evals)
    erlangDParams[n] = integrateD(evals,n)
    push!(N,n)
end

# tempDict = Dict()
# open("dev/erlangParamsData/erlangDParams.json", "r") do f
#     global tempDict
#     tempDict=JSON.parse(f)  # parse and transform data
# end

# erlangDParams = Dict()
# for k in keys(tempDict)
#     Dtemp = zeros(parse(Int,k),parse(Int,k))
#     for i in 1:parse(Int,k)
#         try
#             Dtemp[:,i] = Float64.(tempDict[k][i])
#         catch

#         end
#     end
#     erlangDParams[k] = Dtemp
# end

open("dev/erlangParamsData/erlangDParams.json","w") do f
    JSON.print(f, erlangDParams)
end

@save "dev/erlangParamsData/erlangDParams.jld2" erlangDParams

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
module Rasch
    using Optim
    using LineSearches

    import ElementarySymmetricFunctions: poisbin_sum_taub!, poisbin_sum_taub_dervs_2!

    function logsumexp(x::AbstractVector{T}) where T <: Real
        u, idx = findmax(x)
        u + log1p(sum(exp.(x[1:end .!= idx] .- u)))
    end
    function logsumexpfast(x::AbstractVector{T}) where T <: Real
        u = maximum(x)
        u + log(sum(exp.(x .- u)))
    end

    include("haberman.jl")
    export haberman

end


text:
@testset "Benchmark vs power flow" begin

    # set measurement path
    msr_path = joinpath(mktempdir(),"temp.csv")

    @testset "ACP-rwlav" begin
        # set model
        crit = "rwlav"
        model = _PMD.ACPUPowerModel

        # load data
        data = _PMD.parse_file(_PMDSE.get_enwl_dss_path(ntw, fdr))
        if rm_transfo _PMDSE.rm_enwl_transformer!(data) end
        if rd_lines   _PMDSE.reduce_enwl_lines_eng!(data) end

        # insert the load profiles
        _PMDSE.insert_profiles!(data, season, elm, pfs, t = time_step)

        # transform data model
        data = _PMD.transform_data_model(data);

        # solve the power flow
        pf_result = _PMD.solve_mc_pf(data, model, ipopt_solver)

        # write measurements based on power flow
        _PMDSE.write_measurements!(model, data, pf_result, msr_path, exclude = ["vr","vi"])

        # read-in measurement data and set initial values
        _PMDSE.add_measurements!(data, msr_path, actual_meas = true)
        _PMDSE.assign_start_to_variables!(data)
        _PMDSE.update_all_bounds!(data; v_min = 0.8, v_max = 1.2, pg_min=-1.0, pg_max = 1.0, qg_min=-1.0, qg_max=1.0, pd_min=-1.0, pd_max=1.0, qd_min=-1.0, qd_max=1.0 )

        # set se settings
        data["se_settings"] = Dict{String,Any}("criterion" => crit,
                                           "rescaler" => 100)

        # solve the state estimation
        se_result = _PMDSE.solve_mc_se(data, model, ipopt_solver)

        # tests
        delta, max, avg = _PMDSE.calculate_voltage_magnitude_error(se_result, pf_result)
        @test isapprox(max, 0.0; atol = 1e-6)
        @test isapprox(avg, 0.0; atol = 1e-8)
    end
    @testset "ACP-WLS" begin
        # set model
        crit = "wls"
        model = _PMD.ACPUPowerModel

        # load data
        data = _PMD.parse_file(_PMDSE.get_enwl_dss_path(ntw, fdr))
        if rm_transfo _PMDSE.rm_enwl_transformer!(data) end
        if rd_lines   _PMDSE.reduce_enwl_lines_eng!(data) end

        # insert the load profiles
        _PMDSE.insert_profiles!(data, season, elm, pfs, t = time_step)

        # transform data model
        data = _PMD.transform_data_model(data);

        # solve the power flow
        pf_result = _PMD.solve_mc_pf(data, model, ipopt_solver)

        # write measurements based on power flow
        _PMDSE.write_measurements!(model, data, pf_result, msr_path)

        # read-in measurement data and set initial values
        _PMDSE.add_measurements!(data, msr_path, actual_meas = true)
        _PMDSE.assign_start_to_variables!(data)
        _PMDSE.update_all_bounds!(data; v_min = 0.8, v_max = 1.2, pg_min=-1.0, pg_max = 1.0, qg_min=-1.0, qg_max=1.0, pd_min=-1.0, pd_max=1.0, qd_min=-1.0, qd_max=1.0 )

        # set se settings
        data["se_settings"] = Dict{String,Any}("criterion" => crit,
                                           "rescaler" => 100)

        # solve the state estimation
        se_result = _PMDSE.solve_mc_se(data, model, ipopt_solver)

        # tests
        delta, max, avg = _PMDSE.calculate_voltage_magnitude_error(se_result, pf_result)
        @test isapprox(max, 0.0; atol = 1e-6)
        @test isapprox(avg, 0.0; atol = 1e-8)
    end
    @testset "rACP-rwlav" begin
        # set model
        crit = "rwlav"
        model = _PMDSE.ReducedACPUPowerModel

        # load data
        data = _PMD.parse_file(_PMDSE.get_enwl_dss_path(ntw, fdr))
        if rm_transfo _PMDSE.rm_enwl_transformer!(data) end
        if rd_lines   _PMDSE.reduce_enwl_lines_eng!(data) end

        # insert the load profiles
        _PMDSE.insert_profiles!(data, season, elm, pfs, t = time_step)

        # transform data model
        data = _PMD.transform_data_model(data);

        # solve the power flow
        pf_result = _PMD.solve_mc_pf(data, model, ipopt_solver)

        # write measurements based on power flow
        _PMDSE.write_measurements!(model, data, pf_result, msr_path, exclude = ["vr","vi"])

        # read-in measurement data and set initial values
        _PMDSE.add_measurements!(data, msr_path, actual_meas = true)
        _PMDSE.assign_start_to_variables!(data)
        _PMDSE.update_all_bounds!(data; v_min = 0.8, v_max = 1.2, pg_min=-1.0, pg_max = 1.0, qg_min=-1.0, qg_max=1.0, pd_min=-1.0, pd_max=1.0, qd_min=-1.0, qd_max=1.0 )

        # set se settings
        data["se_settings"] = Dict{String,Any}("criterion" => crit,
                                           "rescaler" => 100)

        # solve the state estimation
        se_result = _PMDSE.solve_mc_se(data, model, ipopt_solver)

        # tests
        delta, max, avg = _PMDSE.calculate_voltage_magnitude_error(se_result, pf_result)
        @test isapprox(max, 0.0; atol = 1e-6)
        @test isapprox(avg, 0.0; atol = 1e-8)
    end
    @testset "rACP-WLS" begin
        # set model
        crit = "wls"
        model = _PMDSE.ReducedACPUPowerModel

        # load data
        data = _PMD.parse_file(_PMDSE.get_enwl_dss_path(ntw, fdr))
        if rm_transfo _PMDSE.rm_enwl_transformer!(data) end
        if rd_lines   _PMDSE.reduce_enwl_lines_eng!(data) end

        # insert the load profiles
        _PMDSE.insert_profiles!(data, season, elm, pfs, t = time_step)

        # transform data model
        data = _PMD.transform_data_model(data);

        # solve the power flow
        pf_result = _PMD.solve_mc_pf(data, model, ipopt_solver)

        # write measurements based on power flow
        _PMDSE.write_measurements!(model, data, pf_result, msr_path)

        # read-in measurement data and set initial values
        _PMDSE.add_measurements!(data, msr_path, actual_meas = true)
        _PMDSE.assign_start_to_variables!(data)
        _PMDSE.update_all_bounds!(data; v_min = 0.8, v_max = 1.2, pg_min=-1.0, pg_max = 1.0, qg_min=-1.0, qg_max=1.0, pd_min=-1.0, pd_max=1.0, qd_min=-1.0, qd_max=1.0 )

        # set se settings
        data["se_settings"] = Dict{String,Any}("criterion" => crit,
                                           "rescaler" => 100)

        # solve the state estimation
        se_result = _PMDSE.solve_mc_se(data, model, ipopt_solver)

        # tests
        delta, max, avg = _PMDSE.calculate_voltage_magnitude_error(se_result, pf_result)
        @test isapprox(max, 0.0; atol = 1e-6)
        @test isapprox(avg, 0.0; atol = 1e-8)
    end
    @testset "ACR-rwlav" begin
        # set model
        crit = "rwlav"
        model = _PMD.ACRUPowerModel

        # load data
        data = _PMD.parse_file(_PMDSE.get_enwl_dss_path(ntw, fdr))
        if rm_transfo _PMDSE.rm_enwl_transformer!(data) end
        if rd_lines   _PMDSE.reduce_enwl_lines_eng!(data) end

        # insert the load profiles
        _PMDSE.insert_profiles!(data, season, elm, pfs, t = time_step)

        # transform data model
        data = _PMD.transform_data_model(data);

        # solve the power flow
        pf_result = _PMD.solve_mc_pf(data, model, ipopt_solver)

        # write measurements based on power flow
        _PMDSE.write_measurements!(model, data, pf_result, msr_path)

        # read-in measurement data and set initial values
        _PMDSE.add_measurements!(data, msr_path, actual_meas = true)
        _PMDSE.assign_start_to_variables!(data)
        _PMDSE.update_all_bounds!(data; v_min = 0.8, v_max = 1.2, pg_min=-1.0, pg_max = 1.0, qg_min=-1.0, qg_max=1.0, pd_min=-1.0, pd_max=1.0, qd_min=-1.0, qd_max=1.0 )

        # set se settings
        data["se_settings"] = Dict{String,Any}("criterion" => crit,
                                           "rescaler" => 100)

        # solve the state estimation
        se_result = _PMDSE.solve_mc_se(data, model, ipopt_solver)

        # tests
        delta, max, avg = _PMDSE.calculate_voltage_magnitude_error(se_result, pf_result)
        @test isapprox(max, 0.0; atol = 1e-6)
        @test isapprox(avg, 0.0; atol = 1e-8)
    end
    @testset "ACR-WLS" begin
        # set model
        crit = "wls"
        model = _PMD.ACRUPowerModel

        # solve the feeders
        data = _PMD.parse_file(_PMDSE.get_enwl_dss_path(ntw, fdr))
        if rm_transfo _PMDSE.rm_enwl_transformer!(data) end
        if rd_lines   _PMDSE.reduce_enwl_lines_eng!(data) end

        # insert the load profiles
        _PMDSE.insert_profiles!(data, season, elm, pfs, t = time_step)

        # transform data model
        data = _PMD.transform_data_model(data);

        # solve the power flow
        pf_result = _PMD.solve_mc_pf(data, model, ipopt_solver)

        # write measurements based on power flow
        _PMDSE.write_measurements!(model, data, pf_result, msr_path)

        # read-in measurement data and set initial values
        _PMDSE.add_measurements!(data, msr_path, actual_meas = true)
        _PMDSE.assign_start_to_variables!(data)
        _PMDSE.update_all_bounds!(data; v_min = 0.8, v_max = 1.2, pg_min=-1.0, pg_max = 1.0, qg_min=-1.0, qg_max=1.0, pd_min=-1.0, pd_max=1.0, qd_min=-1.0, qd_max=1.0 )

        # set se settings
        data["se_settings"] = Dict{String,Any}("criterion" => crit,
                                           "rescaler" => 100)

        # solve the state estimation
        se_result = _PMDSE.solve_mc_se(data, model, ipopt_solver)

        # tests
        delta, max, avg = _PMDSE.calculate_voltage_magnitude_error(se_result, pf_result)
        @test isapprox(max, 0.0; atol = 1e-6)
        @test isapprox(avg, 0.0; atol = 1e-8)
    end
    @testset "rACR-rwlav" begin
        # set model
        crit = "rwlav"
        model = _PMDSE.ReducedACRUPowerModel

        # load data
        data = _PMD.parse_file(_PMDSE.get_enwl_dss_path(ntw, fdr))
        if rm_transfo _PMDSE.rm_enwl_transformer!(data) end
        if rd_lines   _PMDSE.reduce_enwl_lines_eng!(data) end

        # insert the load profiles
        _PMDSE.insert_profiles!(data, season, elm, pfs, t = time_step)

        # transform data model
        data = _PMD.transform_data_model(data);

        # solve the power flow
        pf_result = _PMD.solve_mc_pf(data, _PMD.ACRUPowerModel, ipopt_solver)

        # write measurements based on power flow
        _PMDSE.write_measurements!(_PMD.ACRUPowerModel, data, pf_result, msr_path)

        # read-in measurement data and set initial values
        _PMDSE.add_measurements!(data, msr_path, actual_meas = true)
        _PMDSE.assign_start_to_variables!(data)
        _PMDSE.update_all_bounds!(data; v_min = 0.8, v_max = 1.2, pg_min=-1.0, pg_max = 1.0, qg_min=-1.0, qg_max=1.0, pd_min=-1.0, pd_max=1.0, qd_min=-1.0, qd_max=1.0 )

        # set se settings
        data["se_settings"] = Dict{String,Any}("criterion" => crit,
                                           "rescaler" => 100)

        # solve the state estimation
        se_result = _PMDSE.solve_mc_se(data, model, ipopt_solver)

        # tests
        delta, max, avg = _PMDSE.calculate_voltage_magnitude_error(se_result, pf_result)
        @test isapprox(max, 0.0; atol = 1e-6)
        @test isapprox(avg, 0.0; atol = 1e-8)
    end
    @testset "rACR-WLS" begin
        # set model
        crit = "wls"
        model = _PMDSE.ReducedACRUPowerModel

        # load data
        data = _PMD.parse_file(_PMDSE.get_enwl_dss_path(ntw, fdr))
        if rm_transfo _PMDSE.rm_enwl_transformer!(data) end
        if rd_lines   _PMDSE.reduce_enwl_lines_eng!(data) end

        # insert the load profiles
        _PMDSE.insert_profiles!(data, season, elm, pfs, t = time_step)

        # transform data model
        data = _PMD.transform_data_model(data);

        # solve the power flow
        pf_result = _PMD.solve_mc_pf(data, _PMD.ACRUPowerModel, ipopt_solver)

        # write measurements based on power flow
        _PMDSE.write_measurements!(_PMD.ACRUPowerModel, data, pf_result, msr_path)

        # read-in measurement data and set initial values
        _PMDSE.add_measurements!(data, msr_path, actual_meas = true)
        _PMDSE.assign_start_to_variables!(data)
        _PMDSE.update_all_bounds!(data; v_min = 0.8, v_max = 1.2, pg_min=-1.0, pg_max = 1.0, qg_min=-1.0, qg_max=1.0, pd_min=-1.0, pd_max=1.0, qd_min=-1.0, qd_max=1.0 )

        # set se settings
        data["se_settings"] = Dict{String,Any}("criterion" => crit,
                                           "rescaler" => 100)

        # solve the state estimation
        se_result = _PMDSE.solve_mc_se(data, model, ipopt_solver)

        # tests
        delta, max, avg = _PMDSE.calculate_voltage_magnitude_error(se_result, pf_result)
        @test isapprox(max, 0.0; atol = 1e-6)
        @test isapprox(avg, 0.0; atol = 1e-8)
    end
    @testset "IVR-rwlav" begin
        # set model
        crit = "rwlav"
        model = _PMD.IVRUPowerModel

        # load data
        data = _PMD.parse_file(_PMDSE.get_enwl_dss_path(ntw, fdr))
        if rm_transfo _PMDSE.rm_enwl_transformer!(data) end
        if rd_lines   _PMDSE.reduce_enwl_lines_eng!(data) end

        # insert the load profiles
        _PMDSE.insert_profiles!(data, season, elm, pfs, t = time_step)

        # transform data model
        data = _PMD.transform_data_model(data);

        # solve the power flow
        pf_result = _PMD.solve_mc_pf(data, model, ipopt_solver)

        # write measurements based on power flow
        _PMDSE.write_measurements!(model, data, pf_result, msr_path)

        # read-in measurement data and set initial values
        _PMDSE.add_measurements!(data, msr_path, actual_meas = true)
        _PMDSE.assign_start_to_variables!(data)
        _PMDSE.update_all_bounds!(data; v_min = 0.8, v_max = 1.2, pg_min=-1.0, pg_max = 1.0, qg_min=-1.0, qg_max=1.0, pd_min=-1.0, pd_max=1.0, qd_min=-1.0, qd_max=1.0 )

        # set se settings
        data["se_settings"] = Dict{String,Any}("criterion" => crit,
                                           "rescaler" => 100)

        # solve the state estimation
        se_result = _PMDSE.solve_mc_se(data, model, ipopt_solver)

        # tests
        delta, max, avg = _PMDSE.calculate_voltage_magnitude_error(se_result, pf_result)
        @test isapprox(max, 0.0; atol = 1e-6)
        @test isapprox(avg, 0.0; atol = 1e-8)
    end
    @testset "IVR-WLS" begin
        # set model
        crit = "wls"
        model = _PMD.IVRUPowerModel

        # load data
        data = _PMD.parse_file(_PMDSE.get_enwl_dss_path(ntw, fdr))
        if rm_transfo _PMDSE.rm_enwl_transformer!(data) end
        if rd_lines   _PMDSE.reduce_enwl_lines_eng!(data) end

        # insert the load profiles
        _PMDSE.insert_profiles!(data, season, elm, pfs, t = time_step)

        # transform data model
        data = _PMD.transform_data_model(data);

        # solve the power flow
        pf_result = _PMD.solve_mc_pf(data, model, ipopt_solver)

        # write measurements based on power flow
        _PMDSE.write_measurements!(model, data, pf_result, msr_path)

        # read-in measurement data and set initial values
        _PMDSE.add_measurements!(data, msr_path, actual_meas = true)
        _PMDSE.assign_start_to_variables!(data)
        _PMDSE.update_all_bounds!(data; v_min = 0.8, v_max = 1.2, pg_min=-1.0, pg_max = 1.0, qg_min=-1.0, qg_max=1.0, pd_min=-1.0, pd_max=1.0, qd_min=-1.0, qd_max=1.0 )

        # set se settings
        data["se_settings"] = Dict{String,Any}("criterion" => crit,
                                           "rescaler" => 100)

        # solve the state estimation
        se_result = _PMDSE.solve_mc_se(data, model, ipopt_solver)

        # tests
        delta, max, avg = _PMDSE.calculate_voltage_magnitude_error(se_result, pf_result)
        @test isapprox(max, 0.0; atol = 1e-6)
        @test isapprox(avg, 0.0; atol = 1e-8)
    end
    @testset "rIVR-rwlav" begin
        # set model
        crit = "rwlav"
        model = _PMDSE.ReducedIVRUPowerModel

        # load data
        data = _PMD.parse_file(_PMDSE.get_enwl_dss_path(ntw, fdr))
        if rm_transfo _PMDSE.rm_enwl_transformer!(data) end
        if rd_lines   _PMDSE.reduce_enwl_lines_eng!(data) end

        # insert the load profiles
        _PMDSE.insert_profiles!(data, season, elm, pfs, t = time_step)

        # transform data model
        data = _PMD.transform_data_model(data);

        # solve the power flow
        pf_result = _PMD.solve_mc_pf(data, _PMD.IVRUPowerModel, ipopt_solver)

        # write measurements based on power flow
        _PMDSE.write_measurements!(_PMD.IVRUPowerModel, data, pf_result, msr_path)

        # read-in measurement data and set initial values
        _PMDSE.add_measurements!(data, msr_path, actual_meas = true)
        _PMDSE.assign_start_to_variables!(data)
        _PMDSE.update_all_bounds!(data; v_min = 0.8, v_max = 1.2, pg_min=-1.0, pg_max = 1.0, qg_min=-1.0, qg_max=1.0, pd_min=-1.0, pd_max=1.0, qd_min=-1.0, qd_max=1.0 )

        # set se settings
        data["se_settings"] = Dict{String,Any}("criterion" => crit,
                                           "rescaler" => 100)

        # solve the state estimation
        se_result = _PMDSE.solve_mc_se(data, model, ipopt_solver)

        # tests
        delta, max, avg = _PMDSE.calculate_voltage_magnitude_error(se_result, pf_result)
        @test isapprox(max, 0.0; atol = 1e-6)
        @test isapprox(avg, 0.0; atol = 1e-8)
    end
    @testset "rIVR-WLS" begin
        # set model
        crit = "wls"
        model = _PMDSE.ReducedIVRUPowerModel

        # load data
        data = _PMD.parse_file(_PMDSE.get_enwl_dss_path(ntw, fdr))
        if rm_transfo _PMDSE.rm_enwl_transformer!(data) end
        if rd_lines   _PMDSE.reduce_enwl_lines_eng!(data) end

        # insert the load profiles
        _PMDSE.insert_profiles!(data, season, elm, pfs, t = time_step)

        # transform data model
        data = _PMD.transform_data_model(data);

        # solve the power flow
        pf_result = _PMD.solve_mc_pf(data, _PMD.IVRUPowerModel, ipopt_solver)

        # write measurements based on power flow
        _PMDSE.write_measurements!(_PMD.IVRUPowerModel, data, pf_result, msr_path)

        # read-in measurement data and set initial values
        _PMDSE.add_measurements!(data, msr_path, actual_meas = true)
        _PMDSE.assign_start_to_variables!(data)
        _PMDSE.update_all_bounds!(data; v_min = 0.8, v_max = 1.2, pg_min=-1.0, pg_max = 1.0, qg_min=-1.0, qg_max=1.0, pd_min=-1.0, pd_max=1.0, qd_min=-1.0, qd_max=1.0 )

        # set se settings
        data["se_settings"] = Dict{String,Any}("criterion" => crit,
                                           "rescaler" => 100)

        # solve the state estimation
        se_result = _PMDSE.solve_mc_se(data, model, ipopt_solver)

        # tests
        delta, max, avg = _PMDSE.calculate_voltage_magnitude_error(se_result, pf_result)
        @test isapprox(max, 0.0; atol = 1e-6)
        @test isapprox(avg, 0.0; atol = 1e-8)
    end
end


text:
abstract type Transform end

function otherargs(P)
    return [(k,getfield(P,k)) for k in fieldnames(typeof(P))]
end

struct Shift2D <: Transform
    x::Real
    y::Real
end
Shift(x,y) = Shift2D(x,y)
*(S::Shift2D,V::Vec2) = Vec2(V.x + S.x, V.y + S.y)

struct Scale2D <: Transform
    x::Real
    y::Real
end
Scale(x,y) = Scale2D(x,y)
Scale(x) = Scale(x,x)
*(S::Scale2D,V::Vec2) = Vec2(V.x * S.x, V.y * S.y)
*(S::Scale2D,C::Circle2D) =
        Circle(S*C.center,S.x*C.radius;otherargs(C)[3:end]...)

*(S::Transform,P::Point2D) = Point(S*P.P;otherargs(P)[2:end]...)
*(S::Transform,P::Path2D) =
        Path([S*v for v in P.points];otherargs(P)[2:end]...)
*(S::Transform,P::Polygon2D) =
        Polygon([S*v for v in P.points];otherargs(P)[2:end]...)
*(S::Transform,C::Circle2D) =
        Circle(S*C.center,C.radius;otherargs(C)[3:end]...)
*(S::Transform,L::Label2D) =
        Label(L.s,S*L.location;otherargs(L)[3:end]...)

*(S::Transform,P::Plot2D) = Plot2D([S*e for e in P.elements],P.options)

struct TransformComposition <: Transform
    transforms::Vector{<:Transform}
end

*(S::TransformComposition, T::TransformComposition) =
    TransformComposition([S.transforms;T.transforms])

*(S::Transform,T::TransformComposition) =
    TransformComposition([S;T.transforms])

*(S::TransformComposition, T::Transform) =
    TransformComposition([S.transforms;T])

*(S::Transform,T::Transform) = TransformComposition([S,T])

*(T::TransformComposition, P::Plot2D) =
    foldr(*,[T.transforms;P])

*(T::TransformComposition, P::GraphicElement2D) =
    foldr(*,[T.transforms;P])

function fitsquare(P::Plot2D,ignoreaspect)
    bbox = AsyPlots.boundingbox(P)
    w,h = bbox.xmax-bbox.xmin, bbox.ymax-bbox.ymin
    S = ignoreaspect ? Scale(1/w,1/h) : Scale(1/w)
    S*(Shift(-bbox.xmin,-bbox.ymin)*P)
end

fitsquare(P::Plot2D;ignoreaspect=false) = fitsquare(P,ignoreaspect)

function layout(V::Vector{Plot2D};cols=isqrt(length(V)),
                                  ignoreaspect=false,
                                  margin=0.05,
                                  hmargin=margin,
                                  vmargin=margin)
    U = fitsquare.(V,false)
    push!(U[end].options,(:width,
            floor(Integer,(1+hmargin)*cols*Statistics.mean(
                       merge(_DEFAULT_PLOT2D_KWARGS,
                             Dict(P.options))[:width] for P in U))))
    sum([Shift(mod(i-1,cols)*(1+hmargin),
                    -(i-1)÷cols*(1+vmargin))*P for (i,P) in enumerate(U)])
end