text:
# Copyright 2017-2019 Eric S. Tellez
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using StatsBase
export accuracy, precision_recall, precision, recall, f1, scores

"""
It computes the recall between the gold dataset and the list of predictions `predict`

It applies the desired weighting scheme for binary and multiclass problems
- `:macro` performs a uniform weigth to each class
- `:weigthed` the weight of each class is proportional to its population in gold
- `:micro` returns the global recall, without distinguishing among classes
"""
function recall(gold, predict; weight=:macro)::Float64
    precision, recall, precision_recall_per_class = precision_recall(gold, predict)
    if weight == :macro
        mean(x -> x.second[2], precision_recall_per_class)
    elseif weight == :weighted
        mean(x -> x.second[2] * x.second[3] / length(gold), precision_recall_per_class)
    elseif :micro
        recall
    else
        throw(Exception("Unknown weighting method $weight"))
    end
end

"""
It computes the precision between the gold dataset and the list of predictions `predict`

It applies the desired weighting scheme for binary and multiclass problems
- `:macro` performs a uniform weigth to each class
- `:weigthed` the weight of each class is proportional to its population in gold
- `:micro` returns the global precision, without distinguishing among classes
"""
function precision(gold, predict; weight=:macro)::Float64
    precision, recall, precision_recall_per_class = precision_recall(gold, predict)
    if weight == :macro
        mean(x -> x.second[1], precision_recall_per_class)
    elseif weight == :weighted
        mean(x -> x.second[1] * x.second[3] / length(gold), precision_recall_per_class)
    elseif weight == :micro
        precision
    else
        throw(Exception("Unknown weighting method $weight"))
    end
end

"""
It computes the F1 score between the gold dataset and the list of predictions `predict`

It applies the desired weighting scheme for binary and multiclass problems
- `:macro` performs a uniform weigth to each class
- `:weigthed` the weight of each class is proportional to its population in gold
- `:micro` returns the global F1, without distinguishing among classes
"""
function f1(gold, predict; weight=:macro)::Float64
    precision, recall, precision_recall_per_class = precision_recall(gold, predict)
    if weight == :macro
        mean(x -> 2 * x.second[1] * x.second[2] / (x.second[1] + x.second[2]), precision_recall_per_class)
    elseif weight == :weighted
        mean(x -> 2 * x.second[1] * x.second[2] / (x.second[1] + x.second[2]) * x.second[3]/length(gold), precision_recall_per_class)
    elseif weight == :micro
        2 * (precision * recall) / (precision + recall)
    else
        throw(Exception("Unknown weighting method $weight"))
    end
end

"""
Computes precision, recall, and f1 scores, for global and per-class granularity
"""
function scores(gold, predicted)
    precision, recall, precision_recall_per_class = precision_recall(gold, predicted)
    m = Dict(
        :micro_f1 => 2 * precision * recall / (precision + recall),
		:precision => precision,
		:recall => recall,
		:class_f1 => Dict(),
		:class_precision => Dict(),
		:class_recall => Dict()
    )

    for (k, v) in precision_recall_per_class
        m[:class_f1][k] = 2 * v[1] * v[2] / (v[1] + v[2])
		m[:class_precision][k] = v[1]
		m[:class_recall][k] = v[2]
    end

    m[:macro_recall] = mean(values(m[:class_recall]))
    m[:macro_f1] = mean(values(m[:class_f1]))
    m[:accuracy] = accuracy(gold, predicted)
    m
end

"""
It computes the global and per-class precision and recall values between the gold standard
and the predicted set
"""
function precision_recall(gold, predicted)
    labels = unique(gold)
    M = Dict{typeof(labels[1]), Tuple}()
    tp_ = 0
    tn_ = 0
    fn_ = 0
    fp_ = 0

    for label in labels
        lgold = label .== gold
        lpred = label .== predicted

        tp = 0
        tn = 0
        fn = 0
        fp = 0
        for i in 1:length(lgold)
            if lgold[i] == lpred[i]
                if lgold[i]
                    tp += 1
                else
                    tn += 1
                end
            else
                if lgold[i]
                    fn += 1
                else
                    fp += 1
                end
            end
        end

        tp_ += tp
        tn_ += tn
        fn_ += fn
        fp_ += fp
        M[label] = (tp / (tp + fp), tp / (tp + fn), sum(lgold) |> Int)  # precision, recall, class-population
    end

    tp_ / (tp_ + fp_), tp_ / (tp_ + fn_), M
end

"""
It computes the accuracy score between the gold and the predicted sets
"""
function accuracy(gold, predicted)
    #  mean(gold .== predicted)
    c = 0
    for i in 1:length(gold)
        c += (gold[i] == predicted[i])
    end

    c / length(gold)
end

######### Regression ########

export pearson, spearman, isqerror
"""
Pearson correlation score
"""
function pearson(X::AbstractVector{F}, Y::AbstractVector{F}) where {F <: AbstractFloat}
    X̄ = mean(X)
    Ȳ = mean(Y)
    n = length(X)
    sumXY = 0.0
    sumX2 = 0.0
    sumY2 = 0.0
    for i in 1:n
        x, y = X[i], Y[i]
        sumXY += x * y
        sumX2 += x * x
        sumY2 += y * y
    end
    num = sumXY - n * X̄ * Ȳ
    den = sqrt(sumX2 - n * X̄^2) * sqrt(sumY2 - n * Ȳ^2)
    num / den
end

"""
Spearman rank correleation score
"""
function spearman(X::AbstractVector{F}, Y::AbstractVector{F}) where {F <: AbstractFloat}
    n = length(X)
    x = invperm(sortperm(X))
    y = invperm(sortperm(Y))
    d = x - y
    1 - 6 * sum(d.^2) / (n * (n^2 - 1))
end

"""
Negative squared error (to be used for maximizing algorithms)
"""
function isqerror(X::AbstractVector{F}, Y::AbstractVector{F}) where {F <: AbstractFloat}
    n = length(X)
    d = 0.0

    @inbounds for i in 1:n
        d += (X[i] - Y[i])^2
    end

    -d
end


text:
using Geodesy
using Geodesy: NAD27, OSGB36, Bounds, center, distance
using Base.Test
using Compat
using FixedSizeArrays

#TODO: Break this up / serious rewrite

#############################################
### Decimal <=> Degrees, Minutes, Seconds ###
#############################################

for (decimal, d, m, s) in [(0.013, 0.0, 0.0, 46.8),
                           (-0.013, -0.0, 0.0, 46.8),
                           (-0.263, -0.0, 15.0, 46.8),
                           (-179.51, -179.0, 30.0, 36.0)]
    @test Geodesy.dms2decimal(d, m, s) === decimal
    d2, m2, s2 = Geodesy.decimal2dms(decimal)
    @test d2 === d
    @test m2 === m
    @test_approx_eq s2 s
end

################################################
### Helpers for testing approximate equality ###
################################################

# TODO: Move this to Compat.jl
if VERSION < v"0.4.0-dev+3616"
    fieldnames = names
end

macro type_approx_eq(a, b)
    quote
        @test fieldnames($(esc(a))) == fieldnames($(esc(b)))
        for n in fieldnames($(esc(a)))
            @test_approx_eq $(esc(a)).(n) $(esc(b)).(n)
        end
    end
end

macro xyz_approx_eq(a, b)
    quote
        @test_approx_eq getX($(esc(a))) getX($(esc(b)))
        @test_approx_eq getY($(esc(a))) getY($(esc(b)))
        @test_approx_eq getZ($(esc(a))) getZ($(esc(b)))
    end
end
macro xy_approx_eq(a, b)
    quote
        @test_approx_eq getX($(esc(a))) getX($(esc(b)))
        @test_approx_eq getY($(esc(a))) getY($(esc(b)))
    end
end

macro xyz_approx_eq_eps(a, b, eps)
    quote
        @test_approx_eq_eps getX($(esc(a))) getX($(esc(b))) $(esc(eps))
        @test_approx_eq_eps getY($(esc(a))) getY($(esc(b))) $(esc(eps))
        @test_approx_eq_eps getZ($(esc(a))) getZ($(esc(b))) $(esc(eps))
    end
end
macro xy_approx_eq_eps(a, b, eps)
    quote
        @test_approx_eq_eps getX($(esc(a))) getX($(esc(b))) $(esc(eps))
        @test_approx_eq_eps getY($(esc(a))) getY($(esc(b))) $(esc(eps))
    end
end
macro z_approx_eq_eps(a, b, eps)
    quote
        @test_approx_eq_eps getZ($(esc(a))) getZ($(esc(b))) $(esc(eps))
    end
end

################################################
### Testing parameter stripping and addition ###
################################################

lla_wgs84 = LLA_WGS84(-71.0960, 42.3673, 0)
lla_null = LLA_NULL(-71.0960, 42.3673, 0)
@xyz_approx_eq convert(LLA_NULL, lla_wgs84) lla_null
@xyz_approx_eq convert(LLA{WGS84}, lla_null) lla_wgs84

###################################
### Testing fixed relationships ###
###################################

lla = LLA_WGS84(-71.0960, 42.3673, 0)
lla_ref = LLA_WGS84(-71.09183, 42.36299, 0)

# LLA -> ECEF
ecef = ECEF_WGS84(1529073.1560519305, -4465040.019013103, 4275835.339260309)
@xyz_approx_eq ECEF(lla) ecef

#LLA -> ENU
enu = ENU(-343.493749083977, 478.764855466788, -0.027242885224325164)
@xyz_approx_eq_eps ENU(lla, lla_ref) enu 1e-8
@xyz_approx_eq_eps ENU{lla_ref}(lla) enu 1e-8

#LLA -> ENU (with height in the reference)
enu_wh = ENU(-343.493749083977, 478.764855466788, -0.027242885224325164 - 10.0)
lla_ref_wh = LLA_WGS84(-71.09183, 42.36299, 10.0)
@xyz_approx_eq_eps ENU(lla, lla_ref_wh) enu_wh 1e-8
@xyz_approx_eq_eps ENU{lla_ref_wh}(lla) enu_wh 1e-8

# ECEF -> ENU
@xyz_approx_eq_eps ENU(ecef, lla_ref) enu 1e-8
@xyz_approx_eq_eps ENU{lla_ref}(ecef) enu 1e-8

# ENU -> LLA
@xyz_approx_eq_eps LLA(enu, lla_ref) lla 1e-8
@xyz_approx_eq_eps LLA(ENU{lla_ref}(enu[1], enu[2], enu[3])) lla 1e-8

# Bounds{LLA} -> Bounds{ENU}
bounds = Bounds{LLA_WGS84}(-71.1, -71.094, 42.365, 42.3695)
bounds_enu_ref = Bounds{ENU}(-247.1268196136449, 247.1268196138187, -249.9308954374605, 249.9353534128848)
@type_approx_eq Bounds{ENU}(bounds) bounds_enu_ref

###################################
### Check CRS / Proj4           ###
###################################

# make sure the above are still valid
lla = LLA_WGS84(-71.0960, 42.3673, 0)
ecef_ref = ECEF_WGS84(1529073.1560519305, -4465040.019013103, 4275835.339260309)

# check projections are as expected
@test Geodesy.proj4_str(SRID(LLA_WGS84)) == "+proj=longlat +datum=WGS84 +no_defs"
@test Geodesy.proj4_str(SRID(ECEF_WGS84)) == "+proj=geocent +datum=WGS84 +units=m +no_defs"

# and make sure they behave as anticipated
lla_crs = CRS{SRID(lla)}(lla[1], lla[2], lla[3])
ecef_crs = CRS{SRID(ecef_ref)}(ecef_ref[1], ecef_ref[2], ecef_ref[3])

# Non SRID -> SRID
@xyz_approx_eq_eps CRS{SRID(ecef_ref)}(lla_crs) ecef_ref 1e-8
@xyz_approx_eq_eps CRS{SRID(lla)}(ecef_crs) lla 1e-8

# SRID -> Non SRID
@xyz_approx_eq_eps ECEF(LLA{WGS84}(lla_crs[1], lla_crs[2], lla_crs[3])) ECEF{WGS84}(lla) 1e-8
@xyz_approx_eq_eps ECEF{WGS84}(ecef_crs[1], ecef_crs[2], ecef_crs[3]) ecef_ref 1e-8

#######################################
### Testing ellipsoid relationships ###
#######################################

ecef = ECEF(5.953150599314804e6, 1.5951418955072558e6, 1.6403589592409942e6)
ecef2wgs = LLA{WGS84}(ecef)
ecef2nad = LLA{NAD27}(ecef)
ecef2osgb = LLA{OSGB36}(ecef)

@test getX(ecef2wgs) == getX(ecef2nad)
@test abs(getY(ecef2wgs) - getY(ecef2nad)) > 1e-4
@test abs(getZ(ecef2wgs) - getZ(ecef2nad)) > 10

@test getX(ecef2wgs) == getX(ecef2osgb)
@test abs(getY(ecef2wgs) - getY(ecef2osgb)) > 1e-4
@test abs(getZ(ecef2wgs) - getZ(ecef2osgb)) > 10

#######################################
### Test UTM zone / bound           ###
#######################################

srid = SRID{:EPSG, 32755}     # WGS 84 / UTM zone 55 S (approx Sydney, Australia)
typealias UTM55S CRS{srid}    # type alias for convenience

# get the cnrs (in LLA)
bounds = Bounds{LLA_WGS84}(144.0, 150.0, -80.0, 0.0)
dlon = bounds.max_x - bounds.min_x
dlat = bounds.max_y - bounds.min_y

# random points inside
randLLA() = bounds.min_x + rand() * dlon, bounds.min_y + rand() * dlat, (rand() - .5) * 18000

for i = 1:1000

    # generate (and store for later)
    lla = LLA_WGS84(randLLA())

    # check the zone and bands are reproted correctly
    (zone, band) = Geodesy.utm_zone(lla)

    # correct band and zone zone?
    @test zone == 55
    @test -10 .<= band .< 0

    # correct SRID? (WGS84 datum only)
    (auth, code) = Geodesy.utm_srid(lla)
    srid_out = SRID{auth, code}
    @test srid_out == srid

end

###############################################
### Test bounding box methods for LLA       ###
###############################################

shift = 180.0 - center(bounds)[1]   # temporarily make lon centered on 180 (to test wrapping)

# add cnr points
lla_vec = [LLA_WGS84(bounds.min_x + shift, bounds.min_y, 0.0), LLA_WGS84(bounds.max_x + shift, bounds.min_y, 0.0), LLA_WGS84(bounds.max_x + shift, bounds.max_y, 0.0), LLA_WGS84(bounds.min_x + shift, bounds.max_y, 0.0)]

# random points inside
randLLA() = bounds.min_x + shift + rand() * dlon, bounds.min_y + rand() * dlat, (rand() - .5) * 18000

# make random points
len = length(lla_vec)
resize!(lla_vec, len+1000)
for i = len+1:len+1000
    lla_vec[i] = randLLA()
end

# and bounds them
bounds_est = Bounds(lla_vec)

# remove the shift
bounds_est.min_x = Geodesy.bound_thetad(bounds_est.min_x - shift)
bounds_est.max_x = Geodesy.bound_thetad(bounds_est.max_x - shift)

@type_approx_eq bounds bounds_est

#############################
### Testing vectorization ###
#############################

#
# Step 1 - Grab some points
#
lla_data = [147.80465755709005 -35.35851810277833 277.427 "EPSG"  4326;
            147.80051130896    -35.35195210962327 273.8   "EPSG"  4326;
            147.80010482270845 -35.36183505613432 267.334 "EPSG"  4326;
            147.81089289501872 -35.36325561300782 287.688 "EPSG"  4326]

# test vectorized version vs non vectorized
Xin = convert(Vector{LLA{WGS84}}, lla_data[:, 1:3]; row=true) # = convert(Vector{CRS{srid}}, raw_data[:, 1:3]; row=true)
ecef = geotransform(ECEF{WGS84}, Xin)
for i = 1:length(Xin)
    @xyz_approx_eq ecef[i] ECEF(Xin[i])
end

# test vectorization with Proj4
auth = symbol(lla_data[1, end-1])
code = lla_data[1, end]
srid = SRID{auth, code}

Xin = convert(Vector{CRS{srid}}, lla_data[:, 1:3]; row=true) # = convert(Vector{CRS{srid}}, raw_data[:, 1:3]; row=true)
ecef = geotransform(ECEF{WGS84}, Xin)
for i = 1:length(Xin)
    @xyz_approx_eq ecef[i] ECEF{WGS84}(Xin[i])
end

#############################
### Testing random errors ###
#############################

# lon lat height
randLLA() = (rand() - .5) * 360, (rand() - .5) * 178, (rand() - .5) * 18000
srand(0)

for _ = 1:10_000

    # a random LLA
    x, y, z = randLLA()

    # TODO: figure out how to make this list and draw randomly from it
    # ellipse = Geodesy.get_datums()

    # get LLA and LL types
    lla = LLA{WGS84}(x, y, z)
    ll = LL{WGS84}(x, y)

    # test the center of the bounds
    lla_bounds = Bounds{LLA{WGS84}}(x - 1, x + 1, y - 1, y + 1)
    ll_bounds = Bounds{LLA{WGS84}}(x - 1, x + 1, y - 1, y + 1)

    # check the center is the center
    @xy_approx_eq center(lla_bounds) lla
    @xy_approx_eq center(ll_bounds) ll

    # geotransform to ecef
    ecefa = ECEF(lla)
    ecef = ECEF(ll)

    # test the round trip accuracy
    @xyz_approx_eq_eps LLA(ecefa) lla 1e-6
    @xy_approx_eq_eps LL(ecefa) ll 1e-6

    # test the LLA height match the Euclidean distance between LL and LLA
    @test_approx_eq_eps distance(ecef, ecefa) abs(getZ(lla)) 1e-8

    # get another LL and LLA point
    y, x, z = randLLA()
    lla2 = LLA{WGS84}(x, y, z)
    ll2 = LL{WGS84}(x, y)

    # null ENU point
    enu000 = ENU(0.0, 0.0, 0.0)

    # ecefa is the same point as lla so...
    @xyz_approx_eq_eps ENU(ecefa, lla) enu000   1e-8
    @xy_approx_eq_eps ENU(ecefa, ll) enu000 1e-8
    @z_approx_eq_eps ENU(ecefa, ll) lla     1e-8

    # get another LL and LLA point to use a reference for ENU transforms
    y, x, z = randLLA()
    lla2 = LLA{WGS84}(x, y, z)
    ll2 = LL{WGS84}(x, y)

    # ECEF version of the point
    ecefa2 = ECEF(lla2)
    ecef2 = ECEF(ll2)

    # test LLA -> ECEF -> ENU vs LLA -> ENU
    enu2 = ENU(ecefa, lla2)
    @xyz_approx_eq enu2 ENU(lla, lla2)
    @xy_approx_eq_eps enu2 ENU(ecefa, ll2) 1e-8
    zdiff = getZ(ENU(ecefa, ll2)) - getZ(enu2)
    @test_approx_eq_eps getZ(lla2) zdiff 1e-8

    # Test the transformation matrix approach as well
    (Rf,tf) = Geodesy.geotransform_params(ENU, ECEF, lla2)
    @xyz_approx_eq_eps ENU(Rf * Vec(ecefa) + tf) enu2 1e-8

    # ECEF => ENU => ECEF w/ little change
    enu2v1 = ENU(ecef2, lla)
    @xyz_approx_eq_eps ECEF(enu2v1, lla) ecef2 1e-8

    # Test the transformation matrix approach as well
    (Rb,tb) = Geodesy.geotransform_params(ECEF, ENU, lla)
    @xyz_approx_eq_eps ECEF(Rb * Vec(enu2v1) + tb) ecef2 1e-8

    # ENU => LL same as ENU => ECEF => LLA
    ecef2v1 = ECEF(enu2v1, lla)
    @xyz_approx_eq LLA(enu2v1, lla) LLA(ecef2v1)
    @xy_approx_eq LL(enu2v1, lla) LL(ecef2v1)

    # test distance functions (IDK)
    dist_ecefa = distance(ecefa, ecefa2)
    dist_enua = distance(ENU(lla, lla), ENU(lla2, lla))
    @test_approx_eq dist_ecefa dist_enua

    dist_ecef = distance(ecef, ecef2)
    dist_enu = distance(ENU(ll, ll), ENU(ll2, ll))
    @test_approx_eq dist_ecef dist_enu

end

# create a bunch of LLA points for testing
function LLATestVec(n::Int=1_000_000)

    # get the cnrs (in LLA)
    bounds = Geodesy.Bounds{LLA_WGS84}(144.0, 150.0, -80.0, 0.0)  # utm zone 55
    dlon = bounds.max_x - bounds.min_x
    dlat = bounds.max_y - bounds.min_y

    lla_vec = Vector{LLA{WGS84}}(n)
    for i = 1:n
        lla_vec[i] = LLA{WGS84}(bounds.min_x + rand() * dlon, bounds.min_y + rand() * dlat, (rand() - .5) * 18000)
    end
    return lla_vec

end

function benchmark()

    lla_vec = LLATestVec(50_000_000)  # random variations have been about ~100ms for this many points

    # make proj4 do it for comparison (so there's no chance for type inference anywhere)
    function P4_conv(lla_vec)

        mat = zeros(length(lla_vec), 3)
        @inbounds for i = 1:length(lla_vec)
            mat[i,1], mat[i,2], mat[i,3] = lla_vec[i][1], lla_vec[i][2], lla_vec[i][3]
        end
        Proj4.geotransform!(Geodesy.get_projection(LLA{WGS84}), Geodesy.get_projection(ECEF{WGS84}), mat)

        ecef_vec = Vector{ECEF{WGS84}}(length(lla_vec))
        @inbounds for i = 1:length(lla_vec)
            ecef_vec[i] = ECEF{WGS84}(mat[i,1], mat[i,2], mat[i,3])
        end
        return ecef_vec

    end

    # make the raw function do it (so there's no chance for type inference anywhere)
    function raw_conv(lla_vec)

        d = Geodesy.ellipsoid(WGS84)

        ecef_vec = Vector{ECEF{WGS84}}(length(lla_vec))
        @inbounds for i = 1:length(lla_vec)

            ϕdeg, λdeg, h = lla_vec[i].lat, lla_vec[i].lon, lla_vec[i].alt

            sinϕ, cosϕ = sind(ϕdeg), cosd(ϕdeg)
            sinλ, cosλ = sind(λdeg), cosd(λdeg)

            N = d.a / sqrt(1 - d.e² * sinϕ^2)  # Radius of curvature (meters)

            x = (N + h) * cosϕ * cosλ
            y = (N + h) * cosϕ * sinλ
            z = (N * (1 - d.e²) + h) * sinϕ

            ecef_vec[i] = ECEF{WGS84}(x,y,z)
        end
        return ecef_vec

    end

    # compile stuff
    raw_conv(lla_vec[1:1])
    geotransform(ECEF, lla_vec[1:1])
    geotransform(ECEF{WGS84}, lla_vec[1:1])
    P4_conv(lla_vec[1:1])

    # and test
    println("Convert to ECEF with no type info in the transform")
    @time raw_conv(lla_vec)

    println("Convert to ECEF")
    @time geotransform(ECEF, lla_vec)

    println("Convert to ECEF{WGS84}")
    @time geotransform(ECEF{WGS84}, lla_vec)

    println("Convert to ECEF{WGS84} via Proj4 (matrix style)")
    @time P4_conv(lla_vec)

end

ext: jl


text:
# oneAPI-specific operations on pointers with address spaces

## adrspace aliases

export AS

module AS

const Private    = 0
const Global     = 1
const Constant   = 2
const Local      = 3
const Generic    = 4
const Input      = 5
const Output     = 6
const Count      = 7

end


text:
export periodicswitching, findsmp, hasrepetition
abstract type AbstractPeriodicSwitching end

# Vector{Int} for unconstrained, Vector{Edge} for constrained
duration(seq::AbstractVector) = length(seq)
#duration(seq::AbstractVector{Tuple{Int, Float64}}) = sum(map(p->p[2], seq))

adaptgrowthrate(g, len::Int) = g^(1/len)
#adaptgrowthrate(g, Δt::Float64) = log(g)/Δt
# Vector{Int} for unconstrained, Vector{Edge} for constrained
adaptgrowthrate(g, period::AbstractVector) = adaptgrowthrate(g, duration(period))
#adaptgrowthrate(g, period::AbstractVector{Tuple{Int,Float64}}) = adaptgrowthrate(g, duration(period))

mutable struct DiscretePeriodicSwitching{S<:AbstractHybridSystem, TT} <: AbstractPeriodicSwitching
    s::S
    period::Vector{TT}
    growthrate::Float64
    hashcomputed::Bool
    hash::UInt
end

function DiscretePeriodicSwitching(s, period, growthrate)
    DiscretePeriodicSwitching(s, period, growthrate, false, zero(UInt))
end

#struct ContinuousPeriodicSwitching <: AbstractPeriodicSwitching
#    s::AbstractContinuousSwitchedSystem
#    period::Vector{Tuple{Int, Float64}}
#    growthrate::Float64
#end

function Base.show(io::IO, s::AbstractPeriodicSwitching)
    #print(io, "Periodic switching of growth rate $(s.growthrate) and modes: $(symbol.(s.s, s.period))")# for the transitions $(s.period)")
    print(io, "PSW($(s.growthrate), $(symbol.(s.s, s.period)))")# for the transitions $(s.period)")
end

function periodicswitching(s::AbstractDiscreteSwitchedSystem, period::Vector, growthrate, args...)
    DiscretePeriodicSwitching(s, period, growthrate)
end

_scale(A, ::Nothing) = A
_scale(A, scaling) = A * scaling
_unscale(A, ::Nothing) = A
_unscale(A, scaling) = A / scaling

function periodicswitching(s::AbstractDiscreteSwitchedSystem, period::Vector, A::AbstractMatrix; scaling = nothing)
    lambda = ρ(A)
    growthrate = _scale(adaptgrowthrate(abs(lambda), period), scaling)
    periodicswitching(s, period, growthrate)
end

function periodicswitching(s::AbstractSwitchedSystem, period::Vector; scaling = nothing)
    A = prod(t -> _unscale(integratorfor(s, t), scaling), reverse(period))
    return periodicswitching(s, period, A, scaling = scaling)
end

# Shortcut
function periodicswitching(s::DiscreteSwitchedLinearSystem, period::Vector{Int})
    periodicswitching(s, HybridSystems.OneStateTransition.(period))
end

isperiodic(sw::HybridSystems.DiscreteSwitchingSequence) = source(sw.s, sw) == target(sw.s, sw)
function periodicswitching(s::HybridSystems.DiscreteSwitchingSequence)
    k = repetition(s.seq)
    if iszero(k)
        periodicswitching(s.s, s.seq, s.A)
    else
        periodicswitching(s.s, s.seq[1:k])
    end
end

#function periodicswitching(s::AbstractContinuousSwitchedSystem, seq, growthrate, dt)
#    # seq is a copy since it has been obtained with seq.seq[i:j]
#    mode, Δt = seq[end]
#    seq[end] = mode, dt
#    ContinuousPeriodicSwitching(s, seq, growthrate)
#end

#hasrepetition(s::AbstractPeriodicSwitching) == !iszero(repetition(s.period))

function Base.:(==)(s1::AbstractPeriodicSwitching, s2::AbstractPeriodicSwitching)
    # Shortcut
    hash(s1) == hash(s2) # assumes norepetition
#    if !(s1.s === s2.s)
#        @assert hash(s1) != hash(s2)
#        false
#    elseif !isapprox(s1.growthrate, s2.growthrate)
#        @assert hash(s1) != hash(s2)
#        false
#    else
#        p1 = s1.period
#        p2 = s2.period
#        k1 = length(p1)
#        k2 = length(p2)
#        if k1 != k2
#            k = lcm(k1, k2)
#            if k != k1
#                p1 = repmat(p1, div(k, k1))
#            end
#            if k != k2
#                p2 = repmat(p2, div(k, k2))
#            end
#        else
#            k = k1
#        end
#        for i in 1:k
#            if p1[1:i] == p2[end-i+1:end] && p1[i+1:end] == p2[1:end-i]
#                if hash(s1) != hash(s2)
#                    @show s1
#                    @show s2
#                end
#                return true
#            end
#        end
#        @assert hash(s1) != hash(s2)
#        false
#    end
end

function Base.hash(sw::DiscretePeriodicSwitching, h::UInt)
    if !sw.hashcomputed
        p = sw.period
        k = length(p)
        hs = map(i -> hash([p[end-i+1:end]; p[1:end-i]]), 1:k)
        sw.hash = hash(sort(hs))
        sw.hashcomputed = true
    end
    hash(sw.hash, h)
end

function isbetter(g1, k1, s2::AbstractPeriodicSwitching)
    g2 = s2.growthrate
    k2 = length(s2.period)
    g1 >= g2 * (1 + eps(g2)) || (g1 >= g2 * (1 - eps(g2)) && k1 < k2)
end
function isbetter(s1::AbstractPeriodicSwitching, s2::AbstractPeriodicSwitching)
    g1 = s1.growthrate
    k1 = length(s1.period)
    isbetter(g1, k1, s2)
end

periodicswitchingtype(s::AbstractDiscreteSwitchedSystem) = DiscretePeriodicSwitching{typeof(s)}
#periodicswitchingtype(s::AbstractContinuousSwitchedSystem) = ContinuousPeriodicSwitching

function bestperiod(s::AbstractDiscreteSwitchedSystem, seq::Vector, I, ::AbstractMatrix, Q::AbstractMatrix)
    periodicswitching(s, seq[I], Q), 1
end

#using Optim
#function bestperiod(s::AbstractContinuousSwitchedSystem, seq::Vector{Tuple{Int,Float64}}, I, P::AbstractMatrix, ::AbstractMatrix)
#    mode, Δt = seq[last(I)]
#    T = duration(@view seq[I]) - Δt
#    function f(dt)
#        -adaptgrowthrate(abs(ρ(integratorfor(s, (mode,dt)) * P)), T+dt)
#    end
#    res = optimize(f, min(1e-5, Δt/2), Δt, iterations = 10)
#    -Optim.minimum(res), Optim.minimizer(res)
#end

function repetition(seq)
    k = length(seq)
    for i in 1:div(k, 2)
        if iszero(k % i)
            ok = true
            for j in 2:div(k, i)
                if (@view seq[i*(j-2) .+ (1:i)]) != (@view seq[i*(j-1) .+ (1:i)])
                    ok = false
                    break
                end
            end
            if ok
                return i
            end
        end
    end
    return 0
end

"""
    findsmp(seq::HybridSystems.DiscreteSwitchingSequence)

Extract the cycle of highest growth rate in the sequence `seq`.
"""
function findsmp(seq::HybridSystems.DiscreteSwitchingSequence)
    s = seq.s
    PS = periodicswitchingtype(s)
    smp = nothing
    for i in 1:seq.len
        startNode = state(s, seq.seq[i], false)
        P = _eyet(s, seq.seq[i])
        k = 0
        for j in i:seq.len
            mode = seq.seq[j]
            k = k + nlabels(s, mode)
            Q = integratorfor(s, mode) * P
            if state(s, mode, true) == startNode
                seqij = seq.seq[i:j]
                if iszero(repetition(symbol.(s, seqij)))
                    newsmp, dt = bestperiod(s, seq.seq, i:j, P, Q)
                    notifyperiodic!(s, newsmp)
                    if smp === nothing || isbetter(newsmp, smp)
                        smp = newsmp
                    end
                end
            end
            P = Q
        end
    end

    if smp !== nothing
        updatesmp!(s, smp)
    end
    smp
end


text:
using FractionalDiffEq

function lorenz(t, x, y, z, k)
    a=40
    b=3
    c=10
    d=15
    if k==1
        return a*(y-x)
    elseif k==2
        return c*x-x*z+d*y
    elseif k==3
        return x*y-b*z
    end
end

α0 = [0.96, 0.96, 0.96]
x0 = [1, 2, 3]
h=0.001
prob=FODESystem(lorenz, α0, x0)
T=20
x, y, z=solve(prob, h, T, GLWithMemory())
using Plots
plot(x, z)

####################################
# Or use the detailed model in FractionalDiffEq.jl

a=40; b=3; c=10; d=15
prob = FractionalLorenz(a, b, c, d, 0.96)
x, y, z = solve(prob, h, T, LorenzADM())

using Plots
plot(x, z)


text:
type PI
	PI()=begin
		[
			:(intTerm = 0),
		],
		[
			"",
		],
		new(
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Controller action",
				:Valid=>["Direct","Reverse"],
				:Default=>"Reverse"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Controller mode",
				:Valid=>["Automatic","Manual"],
				:Default=>"Automatic"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Previous scaled bias",
				:Lower=>0,
				:Upper=>1,
				:Default=>0.5
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Proportional term setPoint change filter",
				:Default=>1
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Controller gain",
				:Lower=>0,
				:Upper=>1,
				:Default=>0.5
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Integral time constant",
				:Unit=>"s"
			)),
			control_signal(Dict{Symbol,Any}(
				:Default=>-1000
			)),
			control_signal(Dict{Symbol,Any}(
				:Default=>1000
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Previous scaled input signal",
				:Default=>0.5,
				:PosX=>0,
				:PosY=>0.5,
				:Protected=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Scaled output signal",
				:Default=>0,
				:PosX=>0.54,
				:PosY=>1,
				:Protected=>true
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Scaled setPoint",
				:Lower=>0,
				:Upper=>1,
				:Default=>0.5
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Proportional term",
				:Default=>0,
				:Protected=>true
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Integral term",
				:Default=>0,
				:Protected=>true
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Scaled input variable",
				:Hidden=>true
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Scaled set point",
				:Hidden=>true
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Controller action: (-1) Direct,(1) Reverse",
				:Default=>-1,
				:Hidden=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Sum of proportional, integral and derivative terms",
				:Hidden=>true
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Error definition for proportional term",
				:Hidden=>true
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Variable outp scaled between -1 and 1",
				:Hidden=>true
			)),
			[
				:(input*(MaxInput-MinInput) = Input-MinInput),
				:(setPoint*(MaxInput-MinInput) = SetPoint-MinInput),
				:(error = beta*setPoint - input),
				:(error = input*(beta-1.0)),
				:(action = -1.0),
				:(action = 1.0),
				:(propTerm=error),
				:(outps=outp),
				:(Output = 1),
				:(Output = -1),
				:(Output = outps),
				:(intTime*diff(intTerm) = setPoint - input),
				:(outp = bias + action*gain*(propTerm + intTerm)),
			],
			[
				"","","Error definition","Error definition","","","Calculate proportional term","Scale outp","","","","Calculate integral term","Sum of proportional, integral and derivative terms",
			],
			[:Action,:Mode,:bias,:beta,:gain,:intTime,:MinInput,:MaxInput,],
			[:Input,:Output,:SetPoint,:propTerm,:intTerm,:input,:setPoint,:action,:outp,:error,:outps,]
		)
	end
	Action::DanaSwitcher
	Mode::DanaSwitcher
	bias::positive
	beta::positive
	gain::positive
	intTime::DanaReal
	MinInput::control_signal
	MaxInput::control_signal
	Input::control_signal
	Output::control_signal
	SetPoint::DanaReal
	propTerm::DanaReal
	intTerm::DanaReal
	input::DanaReal
	setPoint::DanaReal
	action::DanaReal
	outp::control_signal
	error::DanaReal
	outps::DanaReal
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	initials::Array{Expr,1}
	initialNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export PI
function setEquationFlow(in::PI)
	addEquation(1)
	addEquation(2)
	let switch=Mode
		if switch=="Automatic"
			addEquation(3)
		elseif switch=="Manual"
			addEquation(4)
		end
	end
	let switch=Action
		if switch=="Direct"
			addEquation(5)
		elseif switch=="Reverse"
			addEquation(6)
		end
	end
	addEquation(7)
	#outps=2*outp-1;
	addEquation(8)
	if outps > 1
		addEquation(9)
	else
		if outps < -1
			addEquation(10)
		else
			addEquation(11)
		end
	end
	addEquation(12)
	addEquation(13)
end
function initial(in::PI)
	addEquation(1)
end
function atributes(in::PI,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/PI"
	drive!(fields,_)
	return fields
end
PI(_::Dict{Symbol,Any})=begin
	newModel=PI()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
using Omega
using Test
using Statistics: mean

function testreplace2()
  θ = normal(0.0, 1.0)
  x = normal(θ, 1.0)
  θnew = normal(100.0, 1.0)
  xnew = Omega.replace(x, θ => θnew)
  @test isapprox(mean(rand(xnew, 1000)), 100, atol = 1.0)
  θnewnew = normal(200.0, 1.0)
  xnewnew = replace(xnew, θnew => θnewnew)
  @test isapprox(mean(rand(xnewnew, 1000)), 200, atol = 1.0)
  y = xnewnew^2
  ynew = replace(y, xnewnew => 10.0)
  @test rand(ynew) == 10.0^2
end

testreplace2()

function replaceconst()
  μ = uniform(1.0, 2.0)
  σ = uniform(1.0, 2.0)
  x = normal(μ, σ)
  xnew = replace(x, σ => μ)
  rand(xnew)
  xnew = replace(x, σ => 3.432)
  rand(xnew)
  xnew = replace(x, μ => 3.432, σ => 21.2)
  rand(xnew)
end

replaceconst()


text:
module TestPackage
using Memoization: @memoize
@memoize foo(x) = x
@memoize foo(x,y) = nothing
foo(1)
end


text:
"""
The `Traitor` package aims to bring a prototype of traits-based dispatch to
Julia. Users will primarily interact with `Traitor` via the `@traitor` macro.

Traits in `Traitor` follow a simple pattern. Firstly, all traits and groupings
of traits are types, and the traits corresponding to a Julia object can be
determined from its type alone. A "class" of traits is defined by an abstract
type, for example `abstract Size`. Examples of each trait are described by its
subtypes, for example `immutable Big <: Size; end` and
`immutable Small <: Size; end`.

Traits can be grouped in a couple of ways. Traits from the *same* class can be
combined in a `Union`, for example `Union{Big, Small}`. This union represents a
*larger* set of Julia objects than either `Big` or `Small` - they represent all
objects possessing either trait.

Conversely, traits from *different* classes can be combined in a `Tuple`, for
example `Tuple{Big, Smelly}` (where `Smelly <: Odor` is another trait from
the `Odor` class). These combinations represent a *smaller* set of Julia objects
than either `Big` or `Smelly` - they represent only objects that possess both
traits simultaneously.

To tell Julia that a type has a certain trait, you overload the constructor of
the trait class (the abstract supertype, e.g. `Size`). For example, we might
have `Size(::Type{Int}) = Small` and `Size{::Type{BigInt}} == Big`. Note that
types are used as both inputs and outputs. A default trait for all objects
can be given, e.g. `Size(::Any) = UnknownSize` (where `UnknownSize <: Size`).

Trait-based methods are created with the `@traitor` macro. Traits are defined
with an extra (optional) set of `::` in the function signature. For example:

    @traitor howbig(::Integer::Big) = "Huge!"
    @traitor howbig(::Integer::Union{Medium, Small}) = "Smallish"
    @traitor howbig(::Integer::Tuple{Big,Smelly}) = "That integer is really big and quite smelly"
    @traitor howbig(::Integer) = "Dunno"

The function is then called in the standard way, and `Traitor` will automatically
manage dispatch. If your code is statically typed, then this will be done at
compile time with no run-time overhead. For example:

    julia> howbig(1)
    "Smallish"

    julia> howbig(BigInt(1))
    "Huge!"

Finally, it is important to note that trait-dispatch occurs after standard
Julia multiple-dispatch. A function may contain many methods, of which
none, some or all will be split into trait-based submethods. The `@traitor`
macro will automatically take ownership of a single method, overwriting any
previous definition.

(see also @traitor, @betray and supertrait)
"""
module Traitor

using Base: @pure, uncompressed_ast, unwrap_unionall, tuple_type_cons
using Core: SimpleVector, svec, CodeInfo

export @traitor, supertrait, betray!

using Cassette

# A context for doing nothing
Cassette.@context DoNothingCtx

"""
    betray!(f, ::Type{TT}) where {TT <: Tuple}

Et tu?

Take over the methods of `f` and readythem to be compatible with traits-based dispatch,
so that new submethods defined by `@traitor` do not overwrite existing default methods
(and used where no more specific trait match is found).

Betraying functions is not a very kind thing to do, and can cause unintended side effects.
"""
function betray!(f, ::Type{TT}) where {TT<:Tuple}
    whereparams = TT isa UnionAll ? typevars(TT.body) : tuple()
    if !isempty(whereparams)
        throw("Where parameters not currently supported")
    end
    # TODO: make this work with TT of the form (Tuple{T, T} where T)

    m = (last ∘ collect ∘ methods)(f, TT)
    mod = (m.module)
    sig = tuple_type_cons(typeof(f), TT)
    @assert m.sig == sig "the provided signature does not exactly match a method. Given $sig, expected $(m.sig) "
    ci = (last ∘ code_lowered)(f, TT)
    _fname = gensym(Symbol(f))

    N = length((TT).parameters)
    argnames = [gensym(Symbol(:arg, i)) for i ∈ 1:N]
    argstyped = map(1:N) do i
        T = unwrap_unionall(TT).parameters[i]
        :($(argnames[i]) :: $T :: $Tuple{})
    end
    fsym = Symbol(f)

    @eval mod begin
        @generated function $_fname($(argnames...))
            return $ci
        end
        $Traitor.@traitor $fsym($(argstyped...)) = $_fname($(argnames...))
    end
end

"""
    supertrait(trait)

For a simple trait type, returns `supertype(t)`, while for a `Union` of traits
is returns their common supertype (or else throws an error).
"""
@generated function supertrait(::Type{T}) where {T}
    !(T isa Union) && return supertype(T)
    traitclass = supertype(T.a)
    while true
        Tb = T.b
        if traitclass != supertype(Tb)
            error("Unions of traits must be in the same class. Got $T")
        end
        (Tb isa Union) || break
    end
    if traitclass == Any
        error("$T does not have a supertrait.")
    end
    return traitclass
end
supertrait(t) = error("Unknown trait $t")

"""
    extract_arg_trait(ex) -> (normal_argument, trait)

Parse an argument signature to extract both the "normal"
typed-variable declaration and the trait expression.

If no trait is present in the expression, return `Tuple{}` for
`trait`.

# Examples

```julia
extract_arg_trait(:(x::Int::Big)) -> (:(x::Int), :(Tuple{Big}))
```
"""
function extract_arg_trait(ex)
    if isa(ex, Symbol)
        return Expr(:(::), ex, :Any), Tuple{} # argument with no type annotation
    end
    if ex.head == :parameters
        # TODO
        error("Keyword arguments are not supported by @traitor")
    end
    if ex.head == :kw # default values, eg f(x::T::A=1)
        # TODO currently these are ignored!!
        warn("Default values are currently ignored by @traitor...")
        arg, traittype = extract_arg_trait(ex.args[1])
        return Expr(:kw, arg, ex.args[2]), traittype
    end
    @assert ex.head == :(::)
    if length(ex.args) == 1 # of form ::Type, etc. We need an actual symbol name as we refer to it in the generated function
        ex.args = Any[gensym(), ex.args[1]]
    end
    isfourcolon = isa(ex.args[2],Expr) && ex.args[2].head == :(::) # x::::A syntax
    if isa(ex.args[1], Symbol) && !isfourcolon
        return ex, Tuple{} # argument with normal type annotation
    end
    if isfourcolon
        traittype = ex.args[2].args[1]
        if !isa(traittype, Expr) || traittype.head != :curly || traittype.args[1] != :Tuple
            traittype = Expr(:curly, :Tuple, traittype)
        end
        arg = (isa(ex.args[1], Symbol) ? Expr(:(::), ex.args[1], :Any) : ex.args[1])
        return arg, traittype
    end
    normalarg = ex.args[1]
    @assert normalarg.head == :(::)
    traittype = ex.args[2]
    if !isa(traittype, Expr) || traittype.head != :curly || traittype.args[1] != Tuple
        traittype = Expr(:curly, :Tuple, traittype)
    end
    return normalarg, traittype
end

"""
Each standard Julia method of a betrayed function may contain an arbitrary
number of trait-based submethods. This generated function returns (and creates,
if necessary, the first time it is called) the table joining the traits to
the submethods
"""
@generated function get_trait_table(f::Function, ::Type{Signature}) where {Signature <: Tuple}
    # This is a questionable usage of a generated function. The idea is that it's an 'elegant' way
    # to create an independant global dictionary for each method specification which can later be mutated.
    # The problem is that if this generated function gets recompiled for some reason (which the compiler
    # is allowed to do!) then the dictionary will be replaced with a new empty dict, deleting our trait
    # table.
    d = Dict{Any, Function}()
    return quote
        $(Expr(:meta, :inline))
        ($d)
    end
end

"""
    @traitor f(x1::T1::Trait1, x2::T2::Trait, ...) =
    @traitor function f(x1::T1::Trait1, x2::T2::Trait, ...); ...; end

The @traitor takes over a method of function `f` to use trait-based dispatch.
If a prior non-@traitor method is defined, then that method will be deleted.
However, methods of `f` with signatures other than `f(::T1, ::T2, ...)` remain
unaffected.

Dispatch on traits occurs after normal Julian multiple-dispatch for further
specialization to sub-methods. The most specific matching trait-signature is
chosen, or an error is thrown in the case of ambiguities.

NOTE: Currently, default values and keyword arguments are not supported.

(See help for `Traitor` to learn how traits are to be used and defined.)
"""
macro traitor(ex)
    isa(ex, Expr) || error("trait expression must be a function method definition")

    # If it's a compact function definition `f(x) = y`, we convert it over to long form function `f(x); y; end`
    if ex.head == :(=) && isa(ex.args[1], Expr) && ex.args[1].head == :call
        body = ex.args[2]
        if !isa(body, Expr) || body.head != :block # Many one-liners don't have a block...
            body = Expr(:block, body)
        end
        ex = Expr(:function, ex.args[1], ex.args[2])
    end

    ex.head == :function || error("trait expression must be a function method definition")

    # Get the function name, signature and traits
    def = ex.args[1]
    body = ex.args[2]

    (isa(def, Expr) && def.head == :call) || error("trait expression must be a function method definition")
    funcname = def.args[1]

    # Go through all the arguments
    n_args = length(def.args) - 1
    args = Vector{Any}()
    argnames = Vector{Any}()
    quotednames = Vector{Any}()
    argtypes = Vector{Any}()
    traits = Vector{Any}()
    for i = 1:n_args
        arg, trait = extract_arg_trait(def.args[1+i])
        push!(args, arg)
        push!(quotednames, Expr(:quote, arg.args[1]))
        push!(argnames, arg.args[1])
        push!(argtypes, arg.args[2])
        push!(traits, trait)
    end

    argnames = (argnames...,) # I'm conused why $((argnames...)) doesn't work in the quote below
    quotednames = (quotednames...,)

    # Make a new name for this specialized function
    internalname = "_"*string(funcname)*"{"
    for i = 1:length(argtypes)
        internalname = internalname * "::" * string(argtypes[i]) * "::" * string(traits[i])
        if i < length(argtypes)
            internalname = internalname * ","
        end
    end
    internalname = internalname * "}"
    internalname = Symbol(internalname)

    dispatchname = gensym(Symbol(funcname, :_dispatched))

    ex = quote
        $Traitor.@generated function $funcname($(args...))
            dict = Traitor.get_trait_table($funcname, $(Expr(:curly, :Tuple, argtypes...)))
            thunk = () -> Traitor.trait_dispatch(dict, $(Expr(:curly, :Tuple, argnames...)))

            # This cassette overdub pass literally does nothing other than normal evaluation.
            # However, if I don't use it, the backedge attachment later on doesn't appear
            # to have the desired effect.
            $dispatchname = $Cassette.overdub($DoNothingCtx(), thunk)
            #$dispatchname = thunk()
            ex = (Expr(:call, $dispatchname, $(QuoteNode.(argnames)...)))

            # Create a codeinfo
            ci = $expr_to_codeinfo($(__module__), [Symbol("#self#"), $(quotednames...)], [], (), ex)

            # Attached edges from MethodInstrances of the `supertrait` function to to this CodeInfo.
            # This should make it so that adding members to a trait relevant to this function
            # triggers recompilation, fixing the #265 equivalent for trait methods.
            ci.edges = $Core.MethodInstance[]
            for TT in [$(traits...)]
                for T in TT.parameters
                    for mi in $Core.Compiler.method_instances($supertrait, Tuple{Type{T}})
                        push!(ci.edges, mi)
                    end
                end
            end

            return ci
        end
        $internalname($(args...)) = $body
        local d = $Traitor.get_trait_table($funcname, $(Expr(:curly, :Tuple, argtypes...)))
        d[Tuple{$(traits...)}] = $internalname
        $funcname
    end
    esc(ex)
end

"""
    expr_to_codeinfo(m::Module, argnames, spnames, sp, e::Expr)

Take an expr (usually a generated function generator) and convert it into a CodeInfo object
(Julia's internal, linear representation of code).

`m` is the module that the CodeInfo should be generated from (used for name resolution)

`argnames` must be an iterable container of symbols describing the CodeInfo's input arguments.
NOTE: the first argument should be given as `Symbol("#self#")`. So if the function is `f(x) = x + 1`,
then `argnames = [Symbol("#self#"), :x]`

`spnames` should be an iterable container of the names of the static parameters to the CodeInfo body
(e.g.) in `f(x::T) where {T <: Int} = ...`, `T` is a static parameter, so `spnames` should be `[:T]`

`sp` should be an iterable container of the static parameters to the CodeInfo body themselves (as
opposed to their names) (e.g.) in `f(x::T) where {T <: Int} = ...`, `T` is a static parameter,
so `sp` should be `[T]`

`e` is the actual expression to lower to CodeInfo. This must be 'pure' in the same sense as generated
function bodies.
"""
function expr_to_codeinfo(m::Module, argnames, spnames, sp, e::Expr)
    lam = Expr(:lambda, argnames,
               Expr(Symbol("scope-block"),
                    Expr(:block,
                         Expr(:return,
                              Expr(:block,
                                   e,
                                   )))))
    ex = if spnames === nothing || isempty(spnames)
        lam
    else
        Expr(Symbol("with-static-parameters"), lam, spnames...)
    end

    # Get the code-info for the generatorbody in order to use it for generating a dummy
    # code info object.
    ci = ccall(:jl_expand_and_resolve, Any, (Any, Any, Core.SimpleVector), ex, m, Core.svec(sp...))
    @assert ci isa Core.CodeInfo "Failed to create a CodeInfo from the given expression. This might mean it contains a closure or comprehension?\n Offending expression: $e"
    ci
end

"""
This function deterimines the dispatch on the traits. Each instance of
this generated function is specialized on the (standard) signature of the
inputs, so the only task remaining is to find the most specific matching
"trait signature".
"""
function trait_dispatch(trait_dictionary::Dict{Any, Function}, ::Type{Sig}) where {Sig <: Tuple}
    n_args = length(Sig.parameters)
    # First check which (if any) of our trait conditions is satisfied by Sig
    matching_traits = Vector{Any}()
    for traits ∈ keys(trait_dictionary)
        match = true
        # Iterate and search for lack of matches
        for n in 1:n_args
            traits_n = traits.parameters[n]

            # The traits for argument i should be a Tuple type
            for i in 1:length(traits_n.parameters)
                traitclass = supertrait(traits_n.parameters[i])

                if !(traitclass(Sig.parameters[n]) <: traits_n.parameters[i])
                    match = false
                    break
                end
            end
        end
        if match
            push!(matching_traits, traits)
        end
    end

    # Now see if there is exactly one most specific matching trait signature
    if length(matching_traits) == 0
        error("Can't find any matching trait method. Got signature $Sig and traits $(collect(keys(trait_dictionary)))")
    elseif length(matching_traits) == 1
        return trait_dictionary[matching_traits[1]]
    else # length(matching_traits) > 1
        # See if there is any matching trait that is more specific than *all*
        # other matching traits, otherwise throw an ambiguity error
        for i = 1:length(matching_traits)
            most_specific = true
            for j = 1:length(matching_traits)
                if i == j
                    continue
                end
                if is_more_specific_traitsig(matching_traits[j], matching_traits[i])
                    most_specific = false
                    break
                end
            end
            if most_specific
                return trait_dictionary[matching_traits[i]]
            end
        end
        error("Ambiguity. Could not determine most specific matching trait. Options are: $matching_traits")
    end
end

function is_more_specific_traitsig(traits, traits2)
    n = length(traits.parameters)
    more_specific = true
    for i = 1:n
        more_specific = more_specific & is_more_specific_traitarg(traits.parameters[i], traits2.parameters[i])
    end

    return more_specific
end

function is_more_specific_traitarg(traits, traits2)
    # Now arrange these into trait classes
    d = Dict{Any, Any}()
    for trait in traits.parameters
        traitclass = supertrait(trait)
        if haskey(d, traitclass)
            error("Encountered trait class $d more than once in $traits")
        else
            d[traitclass] = trait
        end
    end

    d2 = Dict{Any, Any}()
    for trait in traits2.parameters
        traitclass = supertrait(trait)
        if haskey(d2, traitclass)
            error("Encountered trait class $d more than once in $traits2")
        else
            d2[traitclass] = trait
        end
    end

    # For each trait, calculate what is going on
    more_specific = true
    for traitclass in union(keys(d), keys(d2))
        if haskey(d, traitclass)
            if haskey(d2, traitclass)
                if !(d[traitclass] <: d2[traitclass])
                    more_specific = false
                    break
                end
            end
        elseif haskey(d2, traitclass)
            # only in old trait
            more_specific = false
            break
        end
    end

    return more_specific
end

end # module


text:
using DiffEqMonteCarlo
using Test

@time @testset "Monte Carlo Simulations" begin include("monte.jl") end
@time @testset "Monte Carlo Analysis" begin include("analysis.jl") end


text:
"""
    readbin(fname::String,prec::Type,iosize::Tuple,fldidx=1)

Read in a binary file to an Array.
"""
function readbin(fname::String,prec::Type,iosize::Tuple,fldidx=1)
    n3 = 1
    if length(iosize) == 3
        n1,n2,n3 = iosize
    else
        n1,n2 = iosize
    end

    if prec == Float64
        reclen = 8
    else
        reclen = 4
    end

    field = Array{prec,1}(undef,(n1*n2*n3))
    fid = open(fname)
    if fldidx != 1
        seek(fid,n1*n2*n3*(fldidx-1)*reclen)
    end
    read!(fid,field)
    close(fid)
    field = hton.(field)

    n3>1 ? s=(n1,n2,n3) : s=(n1,n2)
    return reshape(field,s)

end

"""
    readbin(flddata::BinData,tidx=1)

Read in a binary file as an array as specified by BinData
"""
readbin(flddata::BinData,tidx=1) = readbin(flddata.fnames[tidx],flddata.precision,flddata.iosize,flddata.fldidx)

"""
    readncdata(var::NCData,i::Union{Colon,Integer}=:)

Read netcdf file as specified in `NCData` argument. Optional
argument `i` can be used to read a specific records / times.
"""
function readncdata(var::NCData,i::Union{Colon,Integer}=:)
    if var.backend == NCDatasets
        ds = Dataset(var.fname)
        ncfilvar = ds[var.varname]
    else
        ncfilvar = NetCDF.open(var.fname,var.varname)
    end

    numdims = ndims(ncfilvar)
    if numdims == 1
        values = ncfilvar[i]
    elseif numdims == 2
        values = ncfilvar[:,i]
    elseif numdims == 3
        values = ncfilvar[:,:,i]
    elseif numdims == 4
        values = ncfilvar[:,:,:,i]
    else
        println("Variables greater than four dimensions not currently supported")
    end

    if var.backend == NCDatasets
        close(ds)
    else
        if pkg_version("NetCDF") < v"0.10"
            NetCDF.close(var.fname)
        else
            finalize(var.fname)
        end
    end

    return values
end

"""
    readdata(flddata,tidx=1)

Generic wrapper function to read data from a file regardless of file/data type.
"""
function readdata(flddata,tidx=1)
    if isa(flddata,BinData)
        res = readbin(flddata,tidx)
    elseif isa(flddata,NCData)
        res = readncdata(flddata,tidx)
    elseif isa(flddata,TileData)
        if isa(flddata.vals,BinData)
            iosize = flddata.vals.iosize
            grid = flddata.tileinfo["XC"].grid
            if length(iosize) == 3
                f = Array{Array{flddata.vals.precision,2},2}(undef,grid.nFaces,iosize[3])
                exarray = MeshArray(grid,f)
            else
                exarray = flddata.tileinfo["XC"]
            end
            res = read(readbin(flddata.vals,tidx),exarray)
        else
            if isa(flddata.vals,MeshArray) || isa(flddata.vals,MeshArrays.gcmfaces)
                res = flddata.vals
            else
                res = flddata.vals[tidx]
            end
        end
    end
    return res
end

"""
    readncfile(fname,backend::Module=NCDatasets)

Read in a NetCDF file and return variables/dimensions as `NCvar` structs, and
    file attributes as `Dict`. Large variables/dimensions are not loaded into
    memory. This can use either `NCDatasets.jl` or `NetCDF.jl`
"""
function readncfile(fname,backend::Module=NCDatasets)

    if backend==NetCDF
        vars,dims,fileatts = readncfile_NetCDF(fname)
    else
        ds = Dataset(fname)

        dims = Dict{String,NCvar}()
        vars = Dict{String,NCvar}()

        for k in keys(ds)
            if ~haskey(dims,k)
                k_units = get(ds[k].attrib,"units","")
                k_atts = Dict(ds[k].attrib)
                if length(dimnames(ds[k])) == 1 && k == dimnames(ds[k])[1] # this variable is also a dimension
                    k_dims = length(ds[k])
                    if length(size(ds[k])) == 1
                        k_values = ds[k][:]
                    elseif length(size(ds[k])) == 2
                        k_values = ds[k][:,:]
                    else
                        k_values = NCData(fname,k,backend,typeof(ds[k].var[:][1])) # Pointer to this variable in the file
                    end
                    if !isa(k_values,NCData) && !any(isa.(k_values,Ref(Missing)))
                        k_values = typeof(k_values[1]).(k_values)
                    end
                    dims[k] = NCvar(k,k_units,k_dims,k_values,k_atts,backend)
                else # Not a dimension
                    k_dims = Array{NCvar,1}()
                    for d in dimnames(ds[k])
                        #global k_dims,dims
                        if !haskey(dims,d)
                            if haskey(ds,d) # Add the dimension
                                d_units = get(ds[d].attrib,"units","")
                                d_dims = length(ds[d])
                                d_values = NCData(fname,d,backend,typeof(ds[d][1])) # Pointer to this variable in the file
                                d_atts = Dict(ds[d].attrib)
                                dims[d] = NCvar(d,d_units,d_dims,d_values,d_atts,backend)
                            else
                                dims[d] = NCvar(d,"",ds.dim[d],[],Dict(),backend)
                            end

                        end
                        k_dims = cat(k_dims,dims[d],dims=1)
                    end
                    hasTimeDim = hastimedim(k_dims)
                    if  ~hasTimeDim && length(k_dims) == 1# extra variable, not main field, just load the data now
                        k_values = ds[k][:]
                    elseif ~hasTimeDim && length(k_dims) == 2
                        k_values = ds[k][:,:]
                    elseif ~hasTimeDim && length(k_dims) == 3
                        k_values = ds[k][:,:,:]
                    elseif hasTimeDim
                        k_values = NCData(fname,k,backend,eltype(ds[k].var)) # Pointer to this variable in the file
                    end

                    if !hasTimeDim && !any(isa.(k_values,Ref(Missing)))
                        k_values = typeof(k_values[1]).(k_values)
                    end

                    vars[k] = NCvar(k,k_units,k_dims,k_values,k_atts,backend)
                end
            end
        end

        atts = ["_FillValue","missing_value","itile","ntile"]
        fileatts = Dict()
        for a in atts
            if haskey(ds.attrib,a)
                fileatts[a] = ds.attrib[a]
            end
        end
        close(ds)
    end
        return vars,dims,fileatts
end

function readncfile_NetCDF(fname)

    backend = NetCDF
    ncfile = NetCDF.open(fname)
    filvars = ncfile.vars
    fildims = ncfile.dim

    dims = Dict{String,NCvar}()
    vars = Dict{String,NCvar}()

    for d in keys(fildims)
        d_units = get(filvars[d].atts,"units","")
        d_atts = filvars[d].atts
        d_dims = Int(fildims[d].dimlen)
        if haskey(filvars,d)
            d_prec = eltype(filvars[d])
            d_values = NCData(fname,d,backend,d_prec)
            vars[d] = NCvar(d,d_units,d_dims,d_values,d_atts,backend)
        else
            d_values = []
        end
        dims[d] = NCvar(d,d_units,d_dims,d_values,d_atts,backend)
    end

    for v in keys(filvars)
        if ~haskey(vars,v)
            v_units = get(filvars[v].atts,"units","")
            v_atts = filvars[v].atts
            dimnames = [d.name for d in filvars[v].dim]
            v_dims = [dims[d] for d in dimnames]
            v_prec = eltype(filvars[v])
            v_vals = NCData(fname,v,backend,v_prec)
            vars[v] = NCvar(v,v_units,v_dims,v_vals,v_atts,backend)
        end
    end

    fileatts = ncfile.gatts
    # Need to take care of order of attributes, read in as Dict and so does not maintain order from file

    return vars,dims,fileatts
end


text:
# Adders
"""
```julia
function HalfAdder(; name)
```

Takes two bits as input, and outputs the sum and the carry

# States
- `sum(t)`
  The sum of the input bits
- `carry(t)`
  The carry generated by the input bits

# Connectors
- `x1`, `x2`
  The two inputs to add
- `y1`
  Output [`DigitalPin`](@ref) corresponding to the sum
- `y2`
  Output [`DigitalPin`](@ref) corresponding to the carry
"""
function HalfAdder(; name)
    @named x1 = DigitalPin()
    @named x2 = DigitalPin()
    @named y1 = DigitalPin()
    @named y2 = DigitalPin()
    @variables sum(t), carry(t)

    eqs = [
        y1.val ~ _xor(x1.val, x2.val)
        y2.val ~ _and(x1.val, x2.val)
        sum ~ y1.val
        carry ~ y2.val
    ]
    ODESystem(eqs, t, [sum, carry], [], systems=[x1, x2, y1, y2], name=name)
end

"""
```julia
function FullAdder(; name)
```

Takes three bits as input, and outputs the sum and the carry

# States
- `sum(t)`
  The sum of the input bits
- `carry(t)`
  The carry generated by the input bits

# Connectors
- `x1`, `x2`, `x3`
  The three inputs to add
- `y1`
  Output [`DigitalPin`](@ref) corresponding to the sum
- `y2`
  Output [`DigitalPin`](@ref) corresponding to the carry
"""
function FullAdder(; name)
    @named x1 = DigitalPin()
    @named x2 = DigitalPin()
    @named x3 = DigitalPin()
    @named y1 = DigitalPin()
    @named y2 = DigitalPin()
    @variables sum(t), carry(t)

    eqs = [
        y1.val ~ _xor(x1.val, x2.val, x3.val)
        y2.val ~ _or(_and(x3.val, _xor(x1.val, x2.val)), _and(x1.val, x2.val))
        sum ~ y1.val
        carry ~ y2.val
    ]
    ODESystem(eqs, t, [sum, carry], [], systems=[x1, x2, x3, y1, y2], name=name)
end

# Multiplexers

# This selects data from the `N` input ports (`d₀` to `dₙ₋₁`)
# using values of `n` select lines, where `N = 2^n`
"""
```julia
function MUX(; name, N=4)
```

Standard Multiplexer. Selects data from `N` input ports using the values
of `n` select lines, where `N=2ⁿ`. For the `i`th input port to be selected,
the values of the select lines should correspond to the binary representation
of `i`.

# Connectors
- `d1`, `d2`, ...
  The `N` input lines
- `s1`, `s2`, ...
  The `n` select lines
- `y`
  The output, selected from one of the `N` input lines
"""
function MUX(; name, N=4)
    n = log2(N)
    try n = Int(n) catch(e) throw("`N` must be a power of 2") end
    s = map(0:n-1) do i
        DigitalPin(; name=Symbol(:s, i))
    end
    d = map(0:N-1) do i
        DigitalPin(; name=Symbol(:d, i))
    end
    @named y = DigitalPin()

    nodes = Num[]
    for i in 1:N
        bin = digits!(zeros(Int64,n), i-1, base=2)
        statelist = Term{Real, Nothing}[]
        for j in 1:n
            varstate = bin[j] == 0 ? _not(s[j].val) : s[j].val
            push!(statelist, varstate)
        end
        push!(nodes, _and(statelist..., d[i].val))
    end

    eqs = Equation[
        y.val ~ _or(nodes...)
    ]

    ODESystem(eqs, t, [], [], systems=[d..., s..., y], name=name)
end

# This selects one of the `N` output ports (`y₀` to `yₙ₋₁`)
# to transmit data `d` using values of `n` select lines, where `N = 2^n`
"""
```julia
function DEMUX(; name, N=4)
```

Standard Demultiplexer. Performs the reverse operation of a [`MUX`](@ref).
Selects one of the `N` output ports to transmit the input `d` using the
values of `n` select lines, where `N=2ⁿ`. For the `i`th output port to be
selected, the values of the select lines should correspond to the binary
representation of `i`.

# Connectors
- `d`
  The input to be transmitted to one of the output lines
- `s1`, `s2`, ...
  The `n` select lines
- `y1`, `y2`, ...
  The `N` output lines
"""
function DEMUX(; name, N=4)
    n = log2(N)
    try n = Int(n) catch(e) throw("`N` must be a power of 2") end
    @named d = DigitalPin()
    s = map(0:n-1) do i
        DigitalPin(; name=Symbol(:s, i))
    end
    y = map(0:N-1) do i
        DigitalPin(; name=Symbol(:y, i))
    end

    eqs = Equation[]
    for i in 1:N
        bin = digits!(zeros(Int64,n), i-1, base=2)
        statelist = Term{Real, Nothing}[]
        for j in 1:n
            varstate = bin[j] == 0 ? _not(s[j].val) : s[j].val
            push!(statelist, varstate)
        end
        push!(eqs, y[i].val ~ _and(statelist..., d.val))
    end

    ODESystem(eqs, t, [], [], systems=[d, s..., y...], name=name)
end

# Encoder-Decoder

# Encodes `N` inputs to `n` outputs, where `N = 2^n`
"""
```julia
function Encoder(; name, N=4)
```

Encodes `N` inputs to `n` outputs,  where `N=2ⁿ`. Exactly one of the inputs should be `1`.
If the `i`th input is `1`, then the output corresponds to the binary representation of `i`.

# Connectors
- `d1`, `d2`, ...
  The `N` input lines
- `y1`, `y2`, ...
  The `n` output lines
"""
function Encoder(; name, N=4)
    n = log2(N)
    try n = Int(n) catch(e) throw("`N` must be a power of 2") end
    d = map(0:N-1) do i
        DigitalPin(; name=Symbol(:d, i))
    end
    y = map(0:n-1) do i
        DigitalPin(; name=Symbol(:y, i))
    end

    nodes = Vector{Term{Real, Nothing}}[]
    i = 0
    for j in 1:n
        counter = 1
        statelist = Term{Real, Nothing}[]
        while i < N
            while counter <= 2^j
                counter > 2^(j-1) && push!(statelist, d[i+1].val)
                counter += 1
                i = i+1
            end
            counter = 1
        end
        i = 0
        push!(nodes, statelist)
    end

    eqs = Equation[]
    for i in n:-1:1
        push!(eqs, y[i].val ~ _or(nodes[i]...))
    end

    ODESystem(eqs, t, [], [], systems=[d..., y...], name=name)
end

# Decodes `n` inputs to `N` outputs, where `N = 2^n`
"""
```julia
function Decoder(; name, n=2)
```

Performs the reverse operation of an [`Encoder`](@ref). Decodes `n` inputs
to `N` outputs,  where `N=2ⁿ`. The `i`th output is `1` if the values of
the select lines correspond to the binary representation of `1`.

# Connectors
- `d1`, `d2`, ...
  The `n` input lines
- `y1`, `y2`, ...
  The `N` output lines
"""
function Decoder(; name, n=2)
    N = 2^n
    d = map(0:n-1) do i
        DigitalPin(; name=Symbol(:d, i))
    end
    y = map(0:N-1) do i
        DigitalPin(; name=Symbol(:y, i))
    end

    nodes = Vector{Term{Real, Nothing}}[]
    for i in 1:N
        bin = digits!(zeros(Int64, n), i-1, base=2)
        statelist = Term{Real, Nothing}[]
        for j in 1:n
            varst = bin[j] == 0 ? _not(d[j].val) : d[j].val
            push!(statelist, varst)
        end
        push!(nodes, statelist)
    end

    eqs = Equation[]
    for i in N:-1:1
        push!(eqs, y[i].val ~ _and(nodes[i]...))
    end

    ODESystem(eqs, t, [], [], systems=[d..., y...], name=name)
end


text:
MLJFluxModel = Union{NeuralNetworkRegressor,
                     MultitargetNeuralNetworkRegressor,
                     NeuralNetworkClassifier,
                     ImageClassifier}

# # EQUALITY

# to address #124 and #129:
MLJModelInterface.deep_properties(::Type{<:MLJFluxModel}) =
    (:optimiser, :builder)

# # CLEAN METHOD

function MLJModelInterface.clean!(model::MLJFluxModel)
    warning = ""
    if model.lambda < 0
        warning *= "Need `lambda ≥ 0`. Resetting `lambda = 0`. "
        model.lambda = 0
    end
    if model.alpha < 0 || model.alpha > 1
        warning *= "Need alpha in the interval `[0, 1]`. "*
            "Resetting `alpha = 0`. "
        model.alpha = 0
    end
    if model.epochs < 0
        warning *= "Need `epochs ≥ 0`. Resetting `epochs = 10`. "
        model.epochs = 10
    end
    if model.batch_size < 0
        warning *= "Need `batch_size ≥ 0`. Resetting `batch_size = 1`. "
        model.batch_size = 1
    end
    if model.acceleration isa CUDALibs && gpu_isdead()
        warning *= "`acceleration isa CUDALibs` "*
            "but no CUDA device (GPU) currently live. "
    end
    if ! (model.acceleration isa CUDALibs || model.acceleration isa CPU1)
        warning *= "`Undefined acceleration, falling back to CPU`"
        model.acceleration = CPU1()
    end
    return warning
end

# # SUPPORT FOR MLJ ITERATION API

# traits:
MLJModelInterface.supports_training_losses(::Type{<:MLJFluxModel}) =
    true
MLJModelInterface.iteration_parameter(model::Type{<:MLJFluxModel}) = :epochs

# method:
MLJModelInterface.training_losses(::MLJFluxModel, report) =
    report.training_losses[2:end] # exclude pre-training loss


text:
import Mjolnir
using IRTools.Inner: Variable, argtypes, arguments

export reorientapply

"Constant"
struct PIConstant{T}
  value::T
end

# State carried around during inversion
struct PIContext
  ir::IR
  invir::IR
  cfg::CFG
  fwd2inv::VarMap                     # Mapping between variable names in forward and inverse
  fwd2invmerged::Dict{Variable, Variable} # Maps variable to merged variable
  paramarg::Variable                  # Variable of params (todo: make this dynamicctx)
  invinarg::Variable                  # argument for input to the inverse
  vartypes::Dict{Variable, Type}
  fwd2inv_block::Dict{BranchBlock, BlockId} # maps bbr::BranchBlock in forward to inverse block id
end

# Creates a block in inverse for each branch point in forward
# Rather than have one block in inverse for one block in forward we have many
# This is because arriving in inverse direction at different branch points
# means that different variables in the block are known, and potentailly
# different subsets of statements are necessary to invert
function initinverseblocks!(ir, invir, fwd2inv, invinarg)
  fwd2inv_block = Dict{BranchBlock, BlockId}()
  for b in blocks(ir)
    for (i, br) in enumerate(branches(b))
      if IRTools.isreturn(br)
        # @assert false
        # @show IRTools.returnvalue(br)
        entryblock = 1
        add!(fwd2inv, IRTools.returnvalue(br), entryblock, invinarg)
        fwd2inv_block[(branch = i, block = b.id)] = 1 # zt: is this correct
      else
        invb = IRTools.block!(invir)
        for arg in arguments(br)
          invarg = argument!(invb; insert = false)
          isvar(arg) && add!(fwd2inv, arg, invb.id, invarg)
          # if isvar(arg)
          #   add!(fwd2inv, arg, invb.id, invarg)
          # else
          #   # ZT: FIXME if branchpoint has constant then in inverse direction
          #   # we have information about output variables of source/child
          #   nothing
          # end
        end
        fwd2inv_block[(branch = i, block = b.id)] = invb.id
      end
    end
  end
  # display(fwd2inv_block)
  fwd2inv_block
end

# Initialise inverse ir with necessary arguments and analyze fwdir to produce context
function setup!(ir, invir)
  invb = IRTools.block(invir, 1)
  selfarg = IRTools.argument!(invb)     # self
  farg = IRTools.argument!(invb)        # f
  typearg = IRTools.argument!(invb)     # types
  targetarg = IRTools.argument!(invb)   # Target
  knownarg = IRTools.argument!(invb)    # Known
  invinarg = IRTools.argument!(invb)    # input to inverse
  paramarg = IRTools.argument!(invb)    # parameters
  cfg = build_cfg(ir)

  vtypes = vartypes(ir)
  fwd2inv = VarMap()
  fwd2invmerged = Dict{Variable, Variable}()
  fwd2inv_block = initinverseblocks!(ir, invir, fwd2inv, invinarg)
  ctx = PIContext(ir, invir, cfg, fwd2inv, fwd2invmerged, paramarg, invinarg, vtypes, fwd2inv_block)
end

# Add return statement to each fwd input clone
function addreturn!(b::Block, invb::Block, ctx::PIContext)
  invretblocks = filter(bbr -> bbr.block == 1, keys(ctx.fwd2inv_block))
  # display(ctx.fwd2inv)
  if b.id == 1 # zt: this is incomplete -- it's currently not usign invretblocks
    IRTools.return!(invb, invb)
    arginv = [getjoin!(arg, invb, ctx) for arg in fargs(b)]
    tpl = IRTools.xcall(Base, :tuple, arginv...)
    retvar = push!(invb, tpl)
    IRTools.return!(invb, retvar)
  end
end

"invert block `b`, store result in `invb`, assume v ∈ `knownvars` is known"
function invert!(b::Block, invb::Block, ctx::PIContext, knownvars::Dict{Variable, Any})
  reversestatementssimple!(b, invb, ctx, knownvars)
  addbranches!(b, invb, ctx)
  addreturn!(b, invb, ctx)
  invb
end

"Undo each operation statement `%a = f(%x, %y, %z)` in `b`, add to `invb`"
function reversestatements!(b::Block, invb::Block, ctx::PIContext, knownvars::Dict{Variable, Any})
  for lhs in reverse(keys(b))
    knownvars
    stmt = b[lhs]
    stmtvars_ = stmtvars(stmt)
    unknownvars = setdiff(stmtvars_, knownvars)
    axes = [stmtvars_; lhs]
    f = stmt.expr.args[1]

    axesids = [i for (i, axis) in enumerate(axes) if axis in unknownvars]
    want = Places{Tuple{axesids...}}
    atypes = stmtargtypes(stmt, ctx.vartypes)
    # getthething(stmt, axesids)
    inv_stmt = xcall(ParametricInversion, :choose, f, atypes, want, ctx.paramarg)
    union!(knownvars, unknownvars)
    var = push!(invb, inv_stmt)

    # add!(ctx.fwd2inv, lhs, invb.id, )

    ## So ignoring the constants,
    ## Let's assume for the minute that we want all teh parameters
    ## Wanted vars is everything that's not constant
    ##
    # println("\n")
  end
  invb
end

"For statement of the form `v1 = f(v2, v3, ..., vn)` produces [v1, v2, ..., vn]"
function enumplaces(vs)
  coords = [vs.var; vs.stmt.expr.args[2:end]]
  [(i = i + 1, val = v) for (i, v) in enumerate(coords)]
end
# baxes(s) = [s.var; s.stmt.expr.args[2:end]]

function getjoin!(v, b, ctx)
  # display(ctx.fwd2invmerged)
  if v in keys(ctx.fwd2invmerged)
    return ctx.fwd2invmerged[v]
  else
    # display(ctx.fwd2inv)
    invv = ctx.fwd2inv[(v, b.id)]
    # Singleton, don't bother invdupl
    if length(invv) == 1
      v_ = ctx.fwd2invmerged[v] = first(invv)
    else
      mergestmt = xcall(PI, :invdupl, invv...)
      v_ = push!(b, mergestmt)
      ctx.fwd2invmerged[v] = v_
    end
    return v_
  end
end

function reversestatementssimple!(b::Block, invb::Block, ctx::PIContext, knownvars::Dict{Variable, Any})
  i(x) = x.i
  val(x) = x.val
  isretvar(i) = i == 2

  methodid_ = methodid(ctx.ir)

  # In reverse order, for each statement
  for vs in reverse(varstatements(b))
    # axes of all outputs that are variables
    targetaxes = places((i for (i, v) in enumplaces(vs) if isvar(v) && !isretvar(i)))
    # @show targetaxes
    # targetaxes = places Places{Tuple{i.(filter(a -> a.i != 1 && isvar(a.val), saxes(vs)))...}}
    # We know the fwd input and any constants
    knownaxes = places((i for (i, v) in enumplaces(vs) if !isvar(v) || isretvar(i)))
    # knownaxes = Places{Tuple{i.(filter(a -> a.i == 1 || !isvar(a.val), saxes(vs)))...}}

    # @assert !isempty(ctx.fwd2inv[(vs.var, invb.id)])
    @assert haskey(knownvars, vs.var) || !isempty(ctx.fwd2inv[(vs.var, invb.id)])
    if isvar(knownvars[vs.var])
      arg = getjoin!(vs.var, invb, ctx) # zt : does this handle multiple merges, correctly?
    else
      arg = knownvars[vs.var]
    end
    consts = val.(filter(a -> !isvar(a.val), enumplaces(vs)))
    args = [arg; consts] # zt: rename args to inverse

    # Location
    loc = Loc(methodid_, vs.var)

    invstmt = xcall(PI, :choose, ctx.paramarg, loc, head(vs.stmt), stmtargtypes(vs.stmt, ctx.vartypes),
                     targetaxes, knownaxes, args...)
    var = push!(invb, invstmt)

    # detuple
    stmtvars_ = stmtvars(vs.stmt)
    for (i, v) in enumerate(stmtvars_)
      v_ = push!(invb, xcall(Core, :getfield, var, i))
      add!(ctx.fwd2inv, v, invb.id, v_)
      # add new known vars so next iteration
      # knownvars has full set of known variables
      updateknown!(knownvars, v, v)
    end
  end
  invb
end

function choosebranch end

# Add branches to invb from each block to its predecessors (in forward)
function addbranches!(b, invb::Block, ctx)
  branches = incomingbranches(b)
  isempty(branches) && return invb  # No incoming branches, nothing to do!

  # Parametrically choose among possible incoming edges
  parentbids = tuple((ctx.fwd2inv_block[bbr] for bbr in branches)...)
  chosen = push!(invb, xcall(PI, :choosebranch, parentbids, ctx.paramarg)) # zt: specialise case when there's only one parent

  # For each parent make a branch point
  for (i, bbr) in enumerate(branches)
    pa = IRTools.block(b.ir, bbr.block)
    br = IRTools.branches(pa)[bbr.branch] # zt: rename bbr.branch to bbr.brid
    destinvb = ctx.fwd2inv_block[bbr]
    stmt = xcall(Base, :(!=), chosen, destinvb)

    # Leave last branch conditionlessss
    condition = i < length(branches) ? push!(invb, stmt) : nothing

    invargs = []
    # This needs to be fwd to inverse in this block
    # display(ctx.fwd2inv)
    for fwdarg in arguments(b)
      mergedvar = getjoin!(fwdarg, invb, ctx)
      push!(invargs, mergedvar)
    end
    IRTools.branch!(invb, destinvb, invargs...; unless = condition)
  end
end

#dm: we might not want to update the val if one already exists?
#dm: how do we choose? is this possible?
function updateknown!(known, var, val)
  known[var] = val
end

"""
When entering `block` at branchpoint `brid`, which variables are known?
"""
function knownvars(block, brid::BranchId)
  known_ = Dict{Variable, Any}()
  # Go in reverse order from branch upwards
  # Because we know all the conditions above (and including)
  # the branch point
  for i = brid:-1:1
    br = IRTools.branches(block)[i]
    for a in br.args
      isvar(a) && updateknown!(known_, a, a)
    end
    isvar(br.condition) && updateknown!(known_, br.condition, i != brid)
  end
  known_
end

"`ir::IR` that computes inverse inverse of `ir`"
function invert(ir::IR)
  invir = IR()        # Invert IR (has one block already)
  ctx = setup!(ir, invir)

  # For each block in forward create inverse block invir
  for (bbr, invbid) in ctx.fwd2inv_block
    invb = IRTools.block(invir, invbid)
    b = IRTools.block(ir, bbr.block)
    known_ = knownvars(b, bbr.branch)
    invert!(b, invb, ctx, known_)
  end
  return invir
end

function reorientir(f::Type{F}, t::Type{T}, target, known) where {F, T}
  fwdir = Mjolnir.trace(Mjolnir.Defaults(), F, t.parameters...)
  IRTools.explicitbranch!(fwdir)  # IR-transforms assumes no implicit branching
  fwdir |> IRTools.expand! |> IRTools.explicitbranch!
  invir = invert(fwdir)
end

reorientir(f::Function, types::NTuple{N, DataType}) where {N} =
  reorientir(typeof(f), Base.to_tuple_type(types))

function fixir!(ci)
  isempty(ci.linetable) && push!(ci.linetable, Core.LineInfoNode(@__MODULE__, ci.parent, :something, 0, 0))
  Core.Compiler.validate_code(ci)
  ci
end

# TODO, will have to add target to this,
# TODO: What about given.  There's the valus themselves, e.g. the input to
# inverse in case of inversion, as well as the tpe of information
function reorient(f::Type{F}, t::Type{T}, target, known) where {F, T}
  invir = reorientir(f, t, target, known)

  # Finalize
  # zt - I wrote this (I think) but I'm not sure what they do?
  dummy() = return
  argnames_ = [Symbol("#self#"), :f, :t, :target, :known, :arg, :θ]
  ci = code_lowered(dummy, Tuple{})[1]
  ci.slotnames = [argnames_...]
  ci = update!(ci, invir)
  fixir!(ci)
  ci
end

@generated function choose(f, t::Type{T}, target::Type{<:Places}, known::Type{<:Places}, z, θ) where T
  return reorient(f, T, target, known)
end


text:
using Base.Test

include("optunity.jl")

testit(x,y,z) = (x-1)^2 + (y-2)^2 + (z+3)^4

vars, details = minimize(testit, num_evals=10000, solver_name="grid search", x=[-5,5], y=[-5,5], z=[-5,5])

@test_approx_eq_eps vars["x"]  1.0 .2
@test_approx_eq_eps vars["y"]  2.0 .2
@test_approx_eq_eps vars["z"] -3.0 .2

testit2(x,y,z) = x^2 + (y-1)^3 + (z-2)^2

vars, details = minimize(testit2, num_evals=10000, x=[-2,2], y=[-2,2], z=[-3,3])

@test_approx_eq_eps vars["x"]  0.0 .2
@test_approx_eq_eps vars["y"] -2.0 .2
@test_approx_eq_eps vars["z"]  2.0 .2

vars, details = minimize(testit2, num_evals=10000, solver_name="nelder-mead", x=[-2,2], y=[-2,2], z=[-3,3])

@test_approx_eq_eps vars["x"]  0.6 .1
@test_approx_eq_eps vars["y"] -2.0 .1
@test_approx_eq_eps vars["z"]  1.1 .1

vars, details = minimize(testit, x=[-5,5], y=[-5,5], z=[-5,5])

@test_approx_eq_eps vars["x"]  1.0 1.
@test_approx_eq_eps vars["y"]  2.0 1.
@test_approx_eq_eps vars["z"] -3.0 1.

vars, details = minimize((x,y,z) -> (x-1)^2 + (y-2)^2 + (z+3)^4, x=[-5,5], y=[-5,5], z=[-5,5])

@test_approx_eq_eps vars["x"]  1.0 1.
@test_approx_eq_eps vars["y"]  2.0 1.
@test_approx_eq_eps vars["z"] -3.0 1.

max_stars_repo_name: xrounder/optunity

max_issues_repo_name: xrounder/optunity

max_forks_repo_name: xrounder/optunity


text:
# p40.jl - eigenvalues of Orr-Sommerfeld operator (compare p38.jl)

R = 5772; clf();
for N = 40:20:100
    # 2nd- and 4th-order differentiation matrices:
    (D,x) = cheb(N); D2 = D^2; D2 = D2[2:N,2:N];
    S = diagm([0; 1 ./(1-x[2:N].^2); 0]);
    D4 = (diagm(1-x.^2)*D^4 - 8*diagm(x)*D^3 - 12*D^2)*S;
    D4 = D4[2:N,2:N];

    # Orr-Sommerfeld operators A,B and generalized eigenvalues:
    I = eye(N-1);
    A = (D4-2*D2+I)/R - 2im*I - 1im*diagm(1-x[2:N].^2)*(D2-I);
    B = D2-I;
    ee = eigvals(A,B);
    i = Int(N/20-1); subplot(2,2,i);
    plot(real(ee),imag(ee),".",markersize=8);
    grid(true); axis("square"); axis([-.8,.2,-1,0]);
    title("N = $N    \$\\lambda_{max}\$ = $(signif(maximum(real(ee)),7))");
end


text:
using PyPlot
using FileIO

mystep(r::AbstractRange) = step(r)
mystep(r::AbstractArray) = r[2] - r[1]

extendrange(r) = range( (first(r)-mystep(r)/2), last(r)+mystep(r)/2, length=length(r)+1)

function plot_heatmap_tracing_prob_vs_c(results, tracing_probs = 0:0.05:1, Cs=0:0.05:1; cmin=nothing, cmax=nothing, logscale=true, addcbar::Bool=true)
  figure(figsize=(10,5))
  reduction = 1 .- Cs / 1.35  |> collect

  if nothing == cmax
    cmax = maximum(results)
  end

  if nothing == cmin
    cmin = minimum(results)
  end

  if logscale
    pcolor(
      extendrange(reduction),
      extendrange(tracing_probs),
      results,
      norm=matplotlib.colors.LogNorm(vmin=cmin, vmax=cmax),
      cmap="nipy_spectral")
      clim(vmin=cmin)
  else
    pcolor(
      extendrange(reduction),
      extendrange(tracing_probs),
      results,
      norm=matplotlib.colors.Normalize(vmin=cmin, vmax=cmax),
      cmap="nipy_spectral")
  end

  if addcbar
    cbar = colorbar()
  end

  #xlabel("f - Stopień redukcji kontaktów")
  xlabel("f - contact reduction rate")
  xticks(
    0:0.1:1,
    ["$(100*f)%" for f in 0:0.1:1])

  #ylabel("b - skuteczność śledzenia kontaktów")
  ylabel("b - contact tracing efficiency")

  gca().invert_yaxis()
  gca().invert_xaxis()
  #savefig("tracking_heatmap_$title.png", bbox_inches="tight")
end

function plot_heatmap_delay_vs_tracking_prob(
    data_path::AbstractString,
    image_path::AbstractString;
    addcbar::Bool=true,
    relative::Union{Nothing,Real}
  )
#  figure(figsize=(10,5))

  data = load(data_path)

  map::AbstractArray{T,2} where T<:Real = data["map"]'
  if relative!==nothing
    map /= relative
  end

  delays::AbstractVector{T} where T<:Real = data["delays"]
  tracking_probs::AbstractVector{T} where T<:Real = data["tracking_probs"]
  c::Real = data["c"]

  pcolor(
    extendrange(delays),
    extendrange(tracking_probs),
    map,
    norm=matplotlib.colors.LogNorm(vmin=minimum(map), vmax=maximum(map)),
    cmap="nipy_spectral"
  )
  #clim(vmin=10^2)

  if addcbar
    colorbar()
  end
  #title("Łączna liczba zarażonych dla redukcji kontaktów o $(100*(1-c/1.35))%")
  xlabel("d - opóźnienie śledzenia w dniach")
  ylabel("b - skuteczność śledzenia kontaktów")

end

function plot_heatmap_mild_detection_vs_tracking_prob(
    data_path::AbstractString,
    image_path::AbstractString
  )
  data = load(data_path)

  map::AbstractArray{T,2} where T<:Real = data["map"]'

  delay::Real = data["delay"]
  c::Real = data["c"]

  mild_detection_probs::AbstractVector{T} where T<:Real = data["mild_detection_probs"]
  tracking_probs::AbstractVector{T} where T<:Real = data["tracking_probs"]

  pcolor(
    extendrange(mild_detection_probs),
    extendrange(tracking_probs),
    map,
    norm=matplotlib.colors.LogNorm(vmin=minimum(map), vmax=maximum(map)),
    cmap="nipy_spectral"
  )
  clim(vmin=10^2)

  gca().invert_yaxis()
  gca().invert_xaxis()

  colorbar()
  #title("Łączna liczba zarażonych \n dla redukcji kontaktów o $(100*(1-c/1.35))% \n i opóźnienia śledzenia kontaktów o $delay dni")
  xlabel("prawdopodobieństwo wykrycia lekkich przypadków")
  ylabel("b - skuteczność śledzenia kontaktów")

  #xlim(minimum(mild_detection_probs), maximum(mild_detection_probs))
  #ylim(minimum(tracking_probs), maximum(tracking_probs))
  #savefig(image_path, bbox_inches="tight")
end

function plot_heatmap_mild_detection_vs_tracing_prob(
    results,
    mild_detections = 0:0.05:1,
    tracking_probs=0:0.05:1;
    cmin=nothing,
    cmax=nothing,
    logscale=true,
    addcbar=true
  )

  if nothing == cmax
    cmax = maximum(results)
  end

  if nothing == cmin
    cmin = minimum(results)
  end

  if logscale
    im = pcolor(
      extendrange(mild_detections),
      extendrange(tracking_probs),
      results,
      norm=matplotlib.colors.LogNorm(vmin=cmin, vmax=cmax),
      cmap="nipy_spectral")
      clim(vmin=cmin)
  else
    im = pcolor(
      extendrange(mild_detections),
      extendrange(tracking_probs),
      results,
      norm=matplotlib.colors.Normalize(vmin=cmin, vmax=cmax),
      cmap="nipy_spectral")
  end

  clim(vmin=cmin)
  if addcbar
    colorbar()
  end

  gca().invert_yaxis()
  gca().invert_xaxis()
  #xlabel("q' - prawdopodobieństwo wykrycia lekkich przypadków")
  #ylabel("b - skuteczność śledzenia kontaktów")
  xlabel("q' - mild case detection probability")
  ylabel("b - contact tracing efficiency")

  im
end

function plot_heatmap_mild_detection_vs_c(results, mild_detection_probs = 0:0.05:1, Cs=0:0.05:1; cmin=nothing, cmax=nothing, logscale=true, addcbar::Bool=true)
  reduction = 1 .- Cs / 1.35  |> collect

  if nothing == cmax
    cmax = maximum(results)
  end

  if nothing == cmin
    cmin = minimum(results)
  end

  if logscale
    im = pcolor(
      extendrange(reduction),
      extendrange(mild_detection_probs),
      results',
      norm=matplotlib.colors.LogNorm(vmin=cmin, vmax=cmax),
      cmap="nipy_spectral")
      clim(vmin=cmin)
  else
    im = pcolor(
      extendrange(reduction),
      extendrange(mild_detection_probs),
      results',
      norm=matplotlib.colors.Normalize(vmin=cmin, vmax=cmax),
      cmap="nipy_spectral")
  end
  if addcbar
    c = colorbar()
  end

  #xlabel("f - Stopień redukcji kontaktów")
  xlabel("f - contact reduction rate")

  xticks(
    0:0.1:1,
    ["$(100*f)%" for f in 0:0.1:1])

  #ylabel("q' - Skuteczność wykrywania lekkich przypadków")
  ylabel("q' - mild case detection probability")

  gca().invert_yaxis()
  gca().invert_xaxis()

  xticks(
    0:0.1:1,
    ["$(Int(100*f))%" for f in 0:0.1:1],
    rotation=60)

  im
end

function plot_heatmap_mild_detection_vs_c(
    data_path::AbstractString,
    image_path::AbstractString
  )
  data = load(data_path)

  map::AbstractArray{T,2} where T<:Real = data["map"]'

  delay::Real = data["delay"]
  Cs::AbstractRange{T} where T<:Real = data["Cs"]
  mild_detection_probs::AbstractRange{T} where T<:Real = data["mild_detection_probs"]
  tracking_prob::Real = data["tracking_prob"]

  reduction = 1 .- Cs / 1.35  |> collect
  pcolormesh(
    extendrange(reduction),
    extendrange(mild_detection_probs),
#reduction,
#mild_detection_probs,
    map,
    norm=matplotlib.colors.LogNorm(vmin=minimum(map), vmax=maximum(map)),
    cmap="nipy_spectral",
    shading="gouraud"
  )
  clim(vmin=10^2)

  colorbar()
  title("Łączna liczba zarażonych \n dla śledzenia kontaktów z prawdopodobieństwem b=$(tracking_prob*100)% i opóźnieniem $delay dni")

  xlabel("f stopień redukcji kontaktów")
  xticks(
    0.0:0.1:1,
    ["$(100*f)%" for f in 0.0:0.1:1])
  xticks(rotation=30)

  ylabel("q' - prawdopodobieństwo wykrycia lekkich przypadków")
  #xlim(minimum(mild_detection_probs), maximum(mild_detection_probs))
  #ylim(minimum(tracking_probs), maximum(tracking_probs))
  savefig(image_path, bbox_inches="tight")
end

function plot_heatmap_c_vs_phone_tracing_usage(results, Cs=0:0.05:1, phone_tracking_usage = 0:0.05:1; cmin=nothing, cmax=nothing, logscale=true, addcbar::Bool=true)
  reduction = 1 .- Cs / 1.35

  if nothing == cmax
    cmax = maximum(results)
  end

  if nothing == cmin
    cmin = minimum(results)
  end

  if logscale
    im = pcolor(
      extendrange(reduction),
      extendrange(phone_tracking_usage),
      results,
      norm=matplotlib.colors.LogNorm(vmin=cmin, vmax=cmax),
      cmap="nipy_spectral")
      clim(vmin=cmin)
  else
    im = pcolor(
      extendrange(reduction),
      extendrange(phone_tracking_usage),
      results,
      norm=matplotlib.colors.Normalize(vmin=cmin, vmax=cmax),
      cmap="nipy_spectral")
  end

  if addcbar
    cbar = colorbar()
  end

  #xlabel("f - Stopień redukcji kontaktów")
  xlabel("f - contact reduction rate")

  xticks(
    0:0.1:1,
    ["$(Int(100*f))%" for f in 0:0.1:1],
    rotation=60)

  #ylabel("u - część populacji używająca aplikacji \n do śledzenia kontaktów")
  ylabel("u - fraction of population \n using the contact tracing app")

  gca().invert_yaxis()
  gca().invert_xaxis()

  return im
end

function plot_heatmap_phone_tracking_usage_vs_tracking_prob(
  results,
  phone_tracking_usages = 0:0.05:1,
  tracking_probs=0:0.05:1;
  cmin=nothing,
  cmax=nothing,
  logscale=true,
  addcbar=true
)

  if nothing == cmax
    cmax = maximum(results)
  end

  if nothing == cmin
    cmin = minimum(results)
  end

  if logscale
    im = pcolor(
      extendrange(phone_tracking_usages),
      extendrange(tracking_probs),
      results,
      norm=matplotlib.colors.LogNorm(vmin=cmin, vmax=cmax),
      cmap="nipy_spectral")
      clim(vmin=cmin)
  else
    im = pcolor(
      extendrange(phone_tracking_usages),
      extendrange(tracking_probs),
      results,
      norm=matplotlib.colors.Normalize(vmin=cmin, vmax=cmax),
      cmap="nipy_spectral")
  end

  clim(vmin=cmin)

  if addcbar
    colorbar()
  end

  gca().invert_xaxis()
  gca().invert_yaxis()
  xlabel("u - część populacji używająca aplikacji śledzącej")
  ylabel("b - skuteczność śledzenia kontaktów")
  return im
end


text:
# Run tests locally
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
# Pkg.test("Algames")

Pkg.activate(joinpath(@__DIR__, "../test"))
Pkg.activate(joinpath(@__DIR__, "../docs"))

# Pkg.add(Pkg.PackageSpec(;name="RobotDynamics", version="0.3.1"))
# Pkg.add(Pkg.PackageSpec(;name="TrajectoryOptimization", version="0.4.1"))
# Pkg.add(Pkg.PackageSpec(;name="Altro", version="0.3.0"))

# Pkg.add("Plots")
# using Plots
# Pkg.remove("Plots")

# Pkg.add("MeshCat")
# using MeshCat
# Pkg.remove("MeshCat")

# Pkg.add("GeometryBasics")
# using GeometryBasics
# Pkg.remove("GeometryBasics")
# pkgs = Pkg.installed()


text:
module SparseMatrixCSRTests

using Test
using SparseMatricesCSR
using SparseArrays
using LinearAlgebra

function test_csr(Bi,Tv,Ti)
  maxnz=10
  maxrows=5
  maxcols=6

  I = rand(Ti(1):Ti(maxrows),maxnz)
  J = rand(Ti(1):Ti(maxcols),maxnz)
  V = rand(Tv,maxnz)

  #I = Ti[1,2,3]
  #J = Ti[1,1,3]
  #V = Tv[3,4,6]

  CSC = sparse(I,J,V)
  if Bi == 1
    CSR = sparsecsr(I,J,V)
    @test CSR == CSC
    @test copy(CSR) == CSC
  end
  CSR = sparsecsr(Val(Bi),I,J,V)
  show(CSR)
  @test CSR == CSC
  @test copy(CSR) == CSC
  @test eltype(CSR) == Tv
  @test isa(CSR,SparseMatrixCSR{Bi,Tv,Ti})

  for i=1:size(CSR,1)
    for j=1:size(CSR,2)
      if (i,j) in zip(I,J)
         CSR[i,j] = eltype(V)(i+j)
         @test CSR[i,j] ≈ eltype(V)(i+j)
      else
         try
          CSR[i,j] = eltype(V)(i+j)
         catch e
          @test isa(e,ArgumentError)
         end
      end
    end
  end

  CSC = sparse(I,J,V,maxrows,maxcols)
  if Bi == 1
    CSR = sparsecsr(I,J,V,maxrows,maxcols)
    @test CSR == CSC
    @test copy(CSR) == CSC
  end
  CSR = sparsecsr(Val(Bi),I,J,V,maxrows,maxcols)
  @test CSR == CSC
  @test copy(CSR) == CSC
  CSR = sparsecsr(Val(Bi),I,J,V,maxrows,maxcols)
  @test CSR == CSC
  @test copy(CSR) == CSC

  CSR2 = convert(typeof(CSR),CSR)
  @test CSR2 === CSR
  CSR2 = convert(SparseMatrixCSR{0,Float64,Int32},CSR)
  @test CSR2 == CSR
  @test copy(CSR2) == CSR
  CSR2 = convert(SparseMatrixCSR{0,Float64,Int32},CSC)
  @test CSR2 == CSC
  @test copy(CSR2) == CSC
  CSR2 = convert(SparseMatrixCSR{0,Float64,Int32},collect(CSC))
  @test CSR2 == CSC
  @test copy(CSR2) == CSC

  @test size(CSR) == size(CSC)
  @test eltype(CSR) == Tv
  @test isa(CSR,SparseMatrixCSR{Bi,Tv,Ti})
  @test issparse(CSR)
  @test getBi(CSR) == Bi
  @test getoffset(CSR) == 1-Bi
  @test nnz(CSR) == nnz(CSC)
  @test length(nonzeros(CSR)) == nnz(CSR)
  @test nonzeros(CSR) === CSR.nzval
  @test colvals(CSR) === CSR.colval
  i,j,v = findnz(CSR)
  csr = sparsecsr(Val(Bi),i,j,v,maxrows,maxcols)
  @test csr == CSR
  @test count(v->v>0,CSR) == count(v->v>0,nonzeros(CSR))

  x = rand(Tv,maxcols)
  y = Vector{Tv}(undef,maxrows)
  z = Vector{Tv}(undef,maxrows)
  mul!(y,CSR,x)
  mul!(z,CSC,x)
  @test y ≈ z
  @test CSR*x ≈ CSC*x

  mul!(y,CSR,x,1,2)
  mul!(z,CSC,x,1,2)
  @test y ≈ z

  out = LinearAlgebra.fillstored!(CSR,3.33)
  @test out === CSR
  LinearAlgebra.fillstored!(CSC,3.33)
  mul!(y,CSR,x)
  mul!(z,CSC,x)
  @test y ≈ z

  _CSR = copy(CSR)
  out = LinearAlgebra.rmul!(CSR,-1)
  @test out === CSR
  @test _CSR ≈ -1*CSR

end

function test_lu(Bi,I,J,V)
  CSR=sparsecsr(Val(Bi),I,J,V)
  CSC=sparse(I,J,V)
  x=rand(3)
  @test norm(CSR\x-CSC\x) < 1.0e-14
  fact=lu(CSR)
  lu!(fact,CSR)
  y=similar(x)
  ldiv!(y,fact,x)
  @test norm(y-CSC\x) < 1.0e-14
end

for Bi in (0,1)
  for Tv in (Float32,Float64)
    for Ti in (Int32,Int64)
      test_csr(Bi,Tv,Ti)
    end
  end
end

I = [1,1,2,2,2,3,3]
J = [1,2,1,2,3,2,3]
V = [4.0,1.0,-1.0,4.0,1.0,-1.0,4.0]
test_lu(0,I,J,V)
test_lu(1,I,J,V)

end # module


text:
require("ParameterVector"); using ParameterVector
syms = [:left, :middle, :right]
lc = LinearConstraints(Float64, syms, [:(left <= middle <= right), :(5 < middle), :(2*(right-middle) == 7)])
@assert lc.A == [1.0 -1.0 0.0; 0.0 1.0 -1.0]
@assert lc.b == zeros(2)
@assert lc.Aeq == [0.0 -2.0 2.0]
@assert lc.beq == [7.0]
@assert lc.lb == [-Inf, 5.0, -Inf]
@assert lc.ub == [Inf, Inf, Inf]


text:
"""
$SIGNATURES

Convenience function to create pairs (commandname => simple lxdef)
"""
lxd(n, k, d="") = "\\" * n => LxDef("\\" * n, k, subs(d))

const LX_INTERNAL_COMMANDS = [
    # ---------------
    # Hyperreferences (see converter/latex/hyperrefs.jl)
    lxd("eqref",    1), # \eqref{id}
    lxd("cite",     1), # \cite{id}
    lxd("citet",    1), # \citet{id}
    lxd("citep",    1), # \citet{id}
    lxd("label",    1), # \label{id}
    lxd("biblabel", 2), # \biblabel{id}{name}
    lxd("toc",      0), # \toc
    lxd("reflink",  1), # \reflink{id}
    # -------------------
    # inclusion / outputs (see converter/latex/io.jl)
    lxd("input",      2), # \input{what}{rpath}
    lxd("output",     1), # \output{rpath}
    lxd("show",       1), # \show{rpath}
    lxd("textoutput", 1), # \textoutput{rpath}
    lxd("textinput",  1), # \textinput{rpath}
    lxd("figalt",     2), # \figalt{alt}{rpath}
    lxd("tableinput", 2), # \tableinput{header}{rpath}
    lxd("literate",   1), # \literate{rpath}
    # ------------------
    # DERIVED / EXPLICIT
    lxd("fig",             1, "\\figalt{}{#1}"),
    lxd("style",           2, "~~~~~~!#2~~~~~~"),
    lxd("tableofcontents", 0, "\\toc"),
    lxd("codeoutput",      1, "\\output{#1}"), # \codeoutput{rpath}
    # ------------------
    # other
    lxd("nonumber", 1),
    ]

"""
$SIGNATURES

Convenience function to create pairs (envdname => simple envdef)
"""
lxe(n, k, d=Pair("", "")) = n => LxDef(n, k, d)

const LX_INTERNAL_ENVIRONMENTS = [
    lxe("equation",  0, raw"\["           => raw"\]"),
    lxe("equation*", 0, raw"\nonumber{\[" => raw"\]}"),
    lxe("align",     0, raw"\[\begin{aligned}"              => raw"\end{aligned}\]"),
    lxe("align*",    0, raw"\nonumber{\[\begin{aligned}"    => raw"\end{aligned}\]}"),
    lxe("aligned",   0, raw"\[\begin{aligned}"              => raw"\end{aligned}\]"),
    lxe("aligned*",  0, raw"\nonumber{\[\begin{aligned}"    => raw"\end{aligned}\]}"),
    lxe("eqnarray",  0, raw"\[\begin{array}{rcl}"           => raw"\end{array}\]"),
    lxe("eqnarray*", 0, raw"\nonumber{\[\begin{array}{rcl}" => raw"\end{array}\]}")
    ]

"""
    GLOBAL_LXDEFS

List of latex definitions accessible to all pages. This is filled when the
config file is read (via `manager/file_utils.jl:process_config`).
"""
const GLOBAL_LXDEFS = LittleDict{String,LxDef}()

"""
$SIGNATURES

Convenience function to allocate default values for global latex commands and environments
accessible throughout the site. See [`resolve_lxobj`](@ref).
"""
 function def_GLOBAL_LXDEFS!()::Nothing
    empty!(GLOBAL_LXDEFS)
    for store in (LX_INTERNAL_ENVIRONMENTS, LX_INTERNAL_COMMANDS)
        for (name, def) in store
            GLOBAL_LXDEFS[name] = def
        end
    end
    return
end