text:
function genRcode(v,av,g,n::Array{ADnode,1},mess,fname=[])
# generate source code for reverse Autodiff pass
# (c) David Barber, University College London 2015
   if isempty(fname)
       iostream=STDOUT
       str="function ("
   else
       iostream=open(fname*"!.jl","w")
       str="function "*fname*"!("
   end
    N=length(n)

    str=str*"v,av,g)"
    println(iostream,str)

    # reverse pass:
    reverseorder=N:-1:1
    println(iostream,"g[$N]=1.0")
    for i in reverseorder[2:N]
        if any(n[i].takederivative)
            count=0
            for c in n[i].children
                str="g[$i]="
                count+=1
                if count==1
                    preval=""
                else
                    preval="g[$i]+"
                end
                str=str*preval*string(mess[i,c])*"("
                for j in n[c].parents
                    tstr="::"*string(typeof(v[j]))
                    str=str*"v[$j]"*tstr*","
                end
                tstr="::"*string(typeof(v[c]))
                str=str*"v[$c]"*tstr*","
                tstr="::"*string(typeof(av[c]))
                str=str*"av[$c]"*tstr*","
                tstr="::"*string(typeof(g[c]))
                str=str*"g[$c]"*tstr*")"
                println(iostream,str)
            end
        end
    end
    println(iostream,"end")
    if !isempty(fname)
        close(iostream)
    end
end

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
module Thermo

using Unitful
using Unitful: K, °C, Pa, hPa
using Unitful: NoUnits, Temperature, Pressure, Amount
#using Unitful: Length, Mass
#using Unitful: Units
#=
# extend Unitful.unit function to be able to find the units of arrays
Unitful.unit(a::AbstractArray)=Unitful.unit(a[1])
Unitful.dimension(a::AbstractArray)=Unitful.dimension(a[1])
=#

# Practical extensions to Unitful:
# Define a nondimensional number to include nondimensional unitless builtin Number types but not dimensional Quantities.
# Define a union of NondimensionalNumber and DimensionlessQuantity types.
# Includes dimensionless ratios of Pressure, Mass, Temperature quantities.
# It's easier to dispatch on Unions than typeintersections.
typealias NondimensionalNumber Union{Real, Complex} # not including Quantities
typealias NondimensionalQuantity Union{NondimensionalNumber, DimensionlessQuantity}
NdN = NondimensionalNumber # for short
NdQ = NondimensionalQuantity

# Extend uconvert to facilitate consistent interfaces for functions of Unitful and/or NondimensionalNumber inputs
# by converting to expected units, dimensional or nondimensional, for all inputs as they come in. E.g:
# most specific method calculates:
#   theta(T,::Temperature, p::Pressure, p0::Pressure, kappa::NondimensionalNumber) = K(T)/Exner(p, p0, kappa)
# one catch-all method handles all unit conversions:
#   theta(T, p, p0, kappa) = theta(uconvert(K,T), uconvert(hPa,p), uconvert(hPa,p0), kappa)
"""
ucon(args...)=Unitful.uconvert(args...)
ucon(a::Units, x::NondimensionalNumber) = a*x
extends and shadows Unitful.uconvert to add units to NondimensionalNumbers
x::NondimensionalNumber assumed to be in a::Units.
examples:
ucon(K, x*°C)      # Unitful.uconvert standard method
ucon(K, x+273.15)  # extended to NondimensionalNumbers
"""
ucon(args...)=Unitful.uconvert(args...)
ucon(a::Unitful.Units, x::NondimensionalNumber) = a*x
ucon{T<:NondimensionalNumber}(a::Unitful.Units, x::AbstractArray{T}) = ucon.(a,x)
#ucon(a::Unitful.Units, x::NondimensionalNumber, sourceunits::Unitful.Units) = ucon(a, x*sourceunits) # superfluous

"""
Nondimensional unit convert
undconvert(a::Unitful.Units, x::Unitful.Quantity) = NoUnits(ucon(a,x)/a)
undconvert(a::Unitful.Units, x::NondimensionalNumber) = x
outputs nondimensional number assumed in the requested units. Assumes nondimensional x already is in units of a.
"""
undconvert(a::Unitful.Units, x::Unitful.Quantity) = NoUnits(ucon(a,x)/a)
undconvert(a::Unitful.Units, x::NondimensionalNumber) = x # == NoUnits(ucon(a,a*x)/a)

export theta, ev,qv,es,qs, Tlcl, Lv, dqsdT, moistad

include("ThermoConstants.jl")
using .constants_unitful
using .constants_unitful: pref, Rd, Rv, Cp, Cpv, Cw
using .constants_unitful: EarthGravity

"T/Θ = Exner(p, p0=1000hPa, kappa=Rd/Cp)  Exner function"
Exner(p::NdN  , p0::NdN, kappa=NoUnits(Rd/Cp)) = (p/p0)^kappa # both unitless
Exner(p::Pressure, p0::Pressure, kappa=NoUnits(Rd/Cp)) = NoUnits(p/p0)^kappa # both pressure, but may be of different units
#Exner(p::Pressure, p0::Pressure=hPa(pref), kappa=NoUnits(Rd/Cp)) = NoUnits(p/p0)^kappa # define intersecting specific case to avoid method ambiguity
# weird cases where one of p,p0 has units but the other doesn't probably should be errors, but here we convert them to have same units as the unitful one:
Exner(p::NdN, p0::Pressure=hPa(pref), kappa=NoUnits(Rd/Cp)) = (p/NoUnits(p0/unit(p0)))^kappa
Exner(p::Pressure, p0::NdN=NoUnits(pref/hPa), kappa=NoUnits(Rd/Cp)) = (NoUnits(p/unit(p))/p0)^kappa
# extend to broadcasting arrays, iterables with .
Exner(p::AbstractArray, p0, kappa)=Exner.(p, p0, kappa)

"Exner exponent kappa(qv[kg/kg])=(Rd/Cp)*(1-0.28*qv) varies with the specific humidity."
kappa(qv::NdQ)=NoUnits(Rd/Cp*(1-0.28*qv))

"""
theta(T,p, p0=1000hPa, kappa=Rd/Cp)
calculates potential temperature from
temperature T (default K), pressure p (default hPa).
Quantities with other compatible units may be specified
by the Unitful module.
"""
# unitless input assumes T in K and returns unitless T values corresponding to K
theta(T::NdN,      p, p0=hPa(pref), kappa=NoUnits(Rd/Cp)) =   T /Exner(p, p0,kappa)
theta(T::Temperature, p, p0=hPa(pref), kappa=NoUnits(Rd/Cp)) = K(T)/Exner(p, p0,kappa)
# output unitful quantities with units K, automatically converts Celsius to K
theta(T::AbstractArray,p,p0=hPa(pref),kappa=NoUnits(Rd/Cp))=theta.(T,p,p0,kappa) # extend to arrays
# theta(T,p::AbstractArray,p0=hPa(pref),kappa=NoUnits(Rd/Cp))=theta.(T,p,p0,kappa)
# multiple dispatch of p is handed by Exner

"""
qv(p/ev)
qv(p::Pressure,ev::Pressure)
qv(p::Amount,ev::Amount)
specific humidity [kg/kg]
"""
#qv(p::Pressure,ev::Pressure) = NoUnits(Rd/(Rv*(p/ev + (Cp/Rv-1.))))
#qv(p::Amount,ev::Amount) = NoUnits(Rd/(Rv*(p/ev + (Cp/Rv-1.))))
qv(poev::Number) = ucon(u"kg/kg",Rd/(Rv*(poev + (Cp/Rv-1.))))
qv(p::Number,ev::Number) = qv(p/ev)
#qv(p::Union{AbstractArray,Number},ev::Union{AbstractArray,Number})=qv.(p,ev) # handle abstractarrays

"""
es(T,p) [hPa] computes saturation vapor pressure based on Wexler's formula,
with enhancement factor for moist air rather than water vapor.
The enhancement factor requires a pressure.
T [degrees C], p [hPa] (note the reversed input order), es [hPa]
From A. L. Buck 1981: JAM, 20, 1527-1532.
SPdeS 7 July 2004
"""
# procedure: convert to expected units, nondimensionalize, calculate
# methods "diagonally dispatch" nondimensionalizing one argument at a time
es(T::NdN,p::NdN=NoUnits(pref/hPa)) = 6.1121hPa*(1.0007 + 3.46e-8*p).*exp((17.502*T)./(240.97 + T))
es(T::Temperature,p::Pressure=pref) = es(NoUnits(°C(T)/°C),NoUnits(p/hPa)) # resolves ambiguity below
# weird corner cases probably should be errors, but are allowed for now:
es(T::Temperature,p::NdN) = es(NoUnits(°C(T)/°C),p)
es(T::NdN,   p::Pressure) = es(T,NoUnits(p/hPa))
# NoUnits form expects T in °C, p in hPa
# uses Buck 1981 numerical constants
#es(T::AbstractArray,p               )=es.(T,p) # broadcasts
#es(T::AbstractArray,p::AbstractArray)=es.(T,p)
#es(T               ,p::AbstractArray)=es.(T,p)

"""
qs(T[°C],p[hPa]) [kg/kg]
Saturation specific humidity based on Wexler's formula for es
with enhancement factor. See es(T,p).
From A. L. Buck 1981: JAM, 20, 1527-1532.
SPdeS 7 July 2004
"""
qs(T::Temperature,p::Pressure=pref)=qv(p,es(T,p))
qs(T::NdN,p::NdN=NoUnits(pref/hPa))=qv(p,NoUnits(es(T,p)/hPa))
#qs(T::AbstractArray,p::AbstractArray)=qs.(T,p)  # broadcast
#qs(T::AbstractArray,p)=qs.(T,p)
#qs(T,p::AbstractArray)=qs.(T,p)

"ev(p,qv) vapor pressure in units of p"
ev(p::Pressure,qv::NdQ) = p/NoUnits(Rd/(Rv*qv) +1)
ev(p::NdN,qv::NdQ) = p/NoUnits(Rd/(Rv*qv) +1)
#ev(p,qv)=ev.(p,qv)
#  ev = p.*qv./(Rd/Rv + qv)
#  Tl = 2840./(3.5*log(T) - log(.01*ev) - 4.805) + 55

"""
moistad(T[°C], p[hPa])
Moist adaibatic lapse rate of _potential_ temperature [K/m]
from Wood and Bretherton 2006
(also Rogers and Yau, but with mixing ratio instead of specific humidity).
(c) Simon de Szoeke, 2011
"""
moistad(T::Temperature,p=pref)=moistad(NoUnits(°C(T)/°C), p)
moistad(T::NdN,p::Pressure=pref)=moistad(T, NoUnits(p/hPa))
function moistad(T::NdN,p::NdN=NoUnits(pref/hPa))
    L=Lv(T)
    qsat=qs(T,p)
    (EarthGravity/Cp) * (1 - (1+L*qsat/(Rd*(T+273.15))) / (1+((L/(T+273.15))^2 *qsat/(Cp*Rv))) )
end

"Lv(T) [J/kg] Latent heat of vaporization of water as a function of temperature [°C]."
Lv(T::NdN=0.0) = (2.501e6 + (Cpv-Cw)*T/unit(Cpv))*u"J/kg"
Lv(T::Temperature=0.0°C) = (2.501e6u"J/kg" + (Cpv-Cw)*°C(T)) |> u"J/kg" # more specific method called first
Lv{Y<:Union{Range,AbstractArray}}(T::Y)=Lv.(T)
# Lv(T::AbstractArray)=Lv.(T)

"dqsdT(T[K],p[hPa]) = dqs/dT = qs(p,T)*Lv(T)/(Rv*T^2) [1/K] using Wexler and Clausius-Clapeyron formulas"
dqsdT(T::Temperature,p::Pressure=pref) = qs(T,p)*Lv(T)/(Rv*(K(T)^2))
dqsdT(T::Temperature,p::NdN) = dqsdT(T,p*hPa)
dqsdT(T::NdN,p::Pressure=pref) = dqsdT(T*K,p)
dqsdT(T::NdN,p::NdN) = dqsdT(T*K,p*hPa)
dqsdT(T::AbstractArray,p=pref) = dqsdT.(T,p) # broadcasts
dqsdT(T,p::AbstractArray) = dqsdT.(T,p)

"""
Tlcl(T[K], p[hPa], qv[kg/kg])
Tlcl(T[K], p[hPa], ev[Pa])
Temperature at the LCL [K]. From Bolton, 1980, MWR, 108, 1046-1053.
"""
Tlcl(T::Temperature,p,qv) = Tlcl(NoUnits(K(T)/K),p,vapor)
Tlcl(T::NdN,p::Pressure,qv) = Tlcl(T,NoUnits(p/Pa),vapor)
Tlcl(T::NdN,p::NdN,ev::Pressure) = (2840./(3.5*log(T) - log(NoUnits(ev/hPa)           ) - 4.805) + 55.)*K
Tlcl(T::NdN,p::NdN,qv::NdQ) =   (2840./(3.5*log(T) - log(p/(Rd/(Rv*NoUnits(qv)) +1)) - 4.805) + 55.)*K
Tlcl(T,p,q)=Tlcl.(T,p,q)

"theta_e(T[K],p[hPa],qv[kg/kg]) from eqn. 43 of Bolton, 1980, MWR, 108, 1046-1053."
theta_e(T::Temperature,p::Pressure,qv::NdQ) = T/Exner(p,p0,kappa(qv)) * exp((3376./Tlcl(T,p,qv) - 2.54)*qv*(1 + 0.81*qv))
theta_e(T,p,qv) = theta_e(ucon(K,T),ucon(hPa,p),NoUnits(ucon(u"kg/kg",qv)))

include("theta_w.jl")

"Blackbody radiation bb(T)"
bb(t::Temperature)=ucon(u"W/m^2", Unitful.σ*K(t)^4) # specific method calculates
bb(t)=bb(ucon(K, t)) # wrapper converts arguments' units

end

ext: jl


text:
abstract type Scheduler end

mutable struct Work
    f::Any
    scheduler::Union{Scheduler,Nothing}
    isenqueued::Bool
    @atomic isdone::Bool
    result::Any
    iserror::Bool
    Work(@nospecialize(f)) = new(f, nothing, true, false, nothing, false)
end

# TODO: Better way to handle this.... parametric type?
scheduler(work::Work) = work.scheduler::Union{WorkStealingScheduler,DepthFirstScheduler}

function run!(work::Work)
    work.result = try
        Base.invokelatest(work.f)
    catch err
        work.iserror = true
        err
    end
    @atomic work.isdone = true
end

mutable struct WorkerState
    @atomic state::UInt
end

struct Worker
    state::WorkerState
    task::Task
end

mutable struct SchedulerState
    @atomic counter::Int
end

struct WorkStealingScheduler <: Scheduler
    queues::Vector{BoundedWorkStealingDeque{Work}}
    workers::Vector{Worker}
    state::SchedulerState
end

function WorkStealingScheduler(; dequesize = 2^7)
    sch = WorkStealingScheduler(
        [BoundedWorkStealingDeque{Work}(dequesize) for _ in 1:Threads.nthreads()],
        Worker[],
        SchedulerState(Threads.nthreads()),
    )
    for tid in 1:Threads.nthreads()
        function workerloop()
            try
                workerloop!(sch)
            catch err
                @error "Worker loop exited" exception = (err, catch_backtrace())
                rethrow()
            end
        end
        task = Task(workerloop)
        @assert task.sticky == true
        ccall(:jl_set_task_tid, Cvoid, (Any, Cint), task, tid - 1)
        push!(sch.workers, Worker(WorkerState(WKR_WORKING), task))
    end
    for wkr in sch.workers
        schedule(wkr.task)
    end
    return sch
end

function trypush!(sch::WorkStealingScheduler, work::Work)
    @assert work.scheduler === nothing
    work.scheduler = sch
    if !trypush!(sch.queues[Threads.threadid()], work)
        work.isenqueued = false
        run!(work)
        @trace(label = :push_failed, threadid = Threads.threadid())
        return false
    end
    @trace(
        label = :push_success,
        threadid = Threads.threadid(),
        length = length(sch.queues[Threads.threadid()]),
    )
    trywakeupall!(sch)

    return true
end

function trywakeupall!(sch)
    state = sch.state
    if (@atomic :monotonic state.counter) < Threads.nthreads() - 1
        for wkr in sch.workers
            if trywakeup!(wkr)
                @atomic state.counter += 1
            end
        end
    end
end

const WKR_WAITING = UInt(0)
const WKR_WORKING = UInt(1)
const WKR_NOTIFIED = UInt(2)

function workerloop!(sch::WorkStealingScheduler)
    nspins = 1000
    wkr::WorkerState = sch.workers[Threads.threadid()].state
    schstate::SchedulerState = sch.state
    while true
        for _ in 1:nspins
            i = Threads.threadid() + 1
            for _ in 1:Threads.nthreads()-1
                if i == Threads.threadid()
                    i += 1
                end
                if i > Threads.nthreads()
                    i = 1
                end
                while true
                    item = maybepopfirst!(sch.queues[i])
                    if item === nothing
                        break
                    else
                        run!(something(item))
                        GC.safepoint()
                    end
                end
            end
            GC.safepoint()
            ccall(:jl_cpu_pause, Cvoid, ())
        end

        @trace(label = :worker_start_sleep, threadid = Threads.threadid())
        @atomic wkr.state = WKR_WAITING
        @atomic schstate.counter -= 1
        wait()
        @assert (@atomic wkr.state) == WKR_NOTIFIED
        @atomic wkr.state = WKR_WORKING
        @trace(label = :worker_done_sleep, threadid = Threads.threadid())
    end
end

function trywakeup!(wkr::Worker)
    wkrstate = wkr.state
    _, issuccess = @atomicreplace wkrstate.state WKR_WAITING => WKR_NOTIFIED
    if issuccess
        schedule(wkr.task)
        @trace(label = :wakeup_scheduled, threadid = Threads.threadid())
        return true
    else
        return false
    end
end

function helpself!(sch::WorkStealingScheduler)
    item = maybepop!(sch.queues[Threads.threadid()])
    if item === nothing
        @trace(label = :helpself_fail, threadid = Threads.threadid())
        return false
    else
        run!(something(item))
        @trace(label = :helpself_success, threadid = Threads.threadid())
        GC.safepoint()
        return true
    end
end

function helpself_until!(f, sch::WorkStealingScheduler)
    while helpself!(sch)
        f() && return true
    end
    return false
end

function helpothers_until!(f, sch::WorkStealingScheduler)
    Threads.nthreads() == 1 && return
    i = Threads.threadid()
    while true
        i += 1
        if i == Threads.threadid()
            i += 1
        end
        if i > Threads.nthreads()
            i = 1
            GC.safepoint()
        end

        while true
            item = maybepopfirst!(sch.queues[i])
            if item === nothing
                f() && return
                @trace(label = :helpothers_empty, threadid = Threads.threadid())
                break
            else
                run!(something(item))
                @trace(label = :helpothers_success, threadid = Threads.threadid())
                f() && return
            end
        end
    end
end

function wait_nothrow(work::Work)
    @inline isdone() = @atomic work.isdone
    if !work.isenqueued
        @assert isdone()
        return
    end
    sch = scheduler(work)
    helpself_until!(isdone, sch) && return
    # TODO: use yieldto to avoid stackoverflow
    helpothers_until!(isdone, sch)
    # Note: assuming `work` is in `sch`, `helpothers_until!` will finish
    # eventually.
end

function Base.wait(work::Work)
    wait_nothrow(work)
    if work.iserror
        throw(work.result)
    end
end

function Base.fetch(work::Work)
    wait(work)
    return work.result
end

function spawn!(@nospecialize(f), sch::Scheduler)
    work = Work(f)
    trypush!(sch, work)
    return work
end

spawn(@nospecialize(f)) = spawn!f, (WORK_STEALING_SCHEDULER[])

const WORK_STEALING_SCHEDULER = Ref{WorkStealingScheduler}()

function init_bounded_scheduler()
    WORK_STEALING_SCHEDULER[] = WorkStealingScheduler()
end


text:
#  Copyright 2019, Jerome Le Ny

using PrivateDynamicData

using Test  # Note: test-specific dependency

using Random
using Distributions

using ControlSystems
using LinearAlgebra

@testset "Complete Test Set" begin
  include("basics_tests.jl")
  if PrivateDynamicData.dpkf_ok
      include("dpkf_tests.jl")
  end
end


text:
export
  AngularMotion,
  set!, interpolate!, update!

using ..Math
using ..Ranger

# ----------------------------------------------------------------------
# Angular motion
# ----------------------------------------------------------------------
mutable struct AngularMotion{T <: AbstractFloat} <: Ranger.AbstractMotion
    # Range is [from:to]
    from::T
    to::T
    rate::T
    time_scale::T  # typically 1000.0ms = 1sec

    # It is enabled by default.
    auto_wrap::Bool

    function AngularMotion{T}() where {T <: AbstractFloat}
        new(0.0, 0.0, 0.0, 1000.0,  true)
    end
end

# Inferred from:
# https://gameprogrammingpatterns.com/game-loop.html

# interpolate is called during rendering (aka visits) which happens more often than updates
# For example, if rendering is 60fps and updates are 30ups then interpolate is
# called twice or more for updates.
function interpolate!(angMo::AngularMotion{T}, t::T) where {T <: AbstractFloat}
    angle = Math.lerp(angMo.from, angMo.to, t)

    # Adjust angle if exited the range
    if angMo.auto_wrap
        if angMo.rate > 0.0
            if angle >= 360.0
                # Wrap range back around
                angMo.from = 0.0
                angMo.to = 361.0 - angle
                # Calc new angle from the adjusted range
                angle = lerp(angMo.from, angMo.to, t)
            end
        else
            if angle <= 0.0
                angMo.from = 360.0
                angMo.to = 358.0 - angle
                angle = lerp(angMo.from, angMo.to, t)
            end
        end
    end

    angle
end

# ----------------------------------------------------------------------
# Linear motion
# ----------------------------------------------------------------------
mutable struct LinearMotion{T <: AbstractFloat} <: Ranger.AbstractMotion
    # Range is [from:to]
    from::T
    to::T
    rate::T
    time_scale::T  # typically 1000.0ms = 1sec

    function LinearMotion{T}() where {T <: AbstractFloat}
        new(0.0, 0.0, 0.0, 1000.0)
    end
end

function interpolate!(angMo::LinearMotion{T}, t::T) where {T <: AbstractFloat}
    Math.lerp(angMo.from, angMo.to, t)
end

# ----------------------------------------------------------------------
# Linear 2D motion
# ----------------------------------------------------------------------
mutable struct Linear2DMotion{T <: AbstractFloat} <: Ranger.AbstractMotion
    # Range is [from:to]
    from::Math.Vector2D{T}
    to::Math.Vector2D{T}
    p::Math.Vector2D{T}
    rate::T
    time_scale::T  # typically 1000.0ms = 1sec

    function Linear2DMotion{T}() where {T <: AbstractFloat}
        new(Math.Vector2D{Float64}(),
            Math.Vector2D{Float64}(),
            Math.Vector2D{Float64}(),
            0.0,
            1000.0)
    end
end

function interpolate!(angMo::Linear2DMotion{T}, t::T) where {T <: AbstractFloat}
    Math.lerp(angMo.from, angMo.to, t, angMo.p)
    angMo.p
end

# ----------------------------------------------------------------------
# Abstract motion
# ----------------------------------------------------------------------
function set!(angMo::Ranger.AbstractMotion, from::T, to::T) where {T <: AbstractFloat}
    angMo.from = from
    angMo.to = to;
end

function set!(angMo::Ranger.AbstractMotion, from::T, to::T, rate::T) where {T <: AbstractFloat}
    set!(angMo, from, to)
    angMo.rate = rate;
end

# dt = milliseconds
# Each update sets a new time window that rendering passes will interpolate
# between.
# Only used for non-2D motion
function update!(angMo::Ranger.AbstractMotion, dt::T) where {T <: AbstractFloat}
    # During each frame the "from" becomes the current "to"
    angMo.from = angMo.to

    # "to" is now moved to the next value
    # divide by time_scale
    angMo.to += angMo.rate * (dt / angMo.time_scale);
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
module test_lupin_partition

using Test
using Lupin

V = collect(1:3)
A = partition(V, cols=5, fill=0)
@test A == [ 1  2  3  0  0]

V = collect(1:13)
A = partition(V, cols=5, fill=0)
@test A == [ 1   2   3  4   5
             6   7   8  9  10
            11  12  13  0   0]

V = split("a b c d e f g h i j")
A = partition(V, cols=5, fill="")
@test A == ["a"  "b"  "c"  "d"  "e"
            "f"  "g"  "h"  "i"  "j"]

@test partition(1:5, cols=2, fill=0) == [1 2
                                         3 4
                                         5 0]

end # module test_lupin_partition


text:
macro dshow( args... )
        local d = esc( :debug_print )
        return :( if $d ; @show( debug_print ) ;  @show( $(args...) ); end )
end

macro dprint( args... )
        local d = esc( :debug_print )
        return :( if $d ; print( $(args...) ); end )
end

macro dd( args... )
        if isdefined( :implement_debug_print ) && implement_debug_print
                local d = esc( :debug_print )
                return :( if $d ; $(args...) ; end )
        else
                return :()
        end
end


text:
@testset "Dual" begin
  QZ = Hecke.QmodnZ()
  for T in [fmpz, Int, fmpq, Rational{BigInt}]
    @test @inferred QZ(0) == T(0)
    @test @inferred T(0) == QZ(0)
  end
  @test @inferred iszero(QZ(0))
  @test @inferred !iszero(QZ(1//2))
  for T in [fmpq, Rational{BigInt}]
    @test @inferred QZ(1//2) != T(1)
    @test @inferred T(1) != QZ(1//2)
  end
end


text:
#=
This file is auto-generated. Do not edit.
=#

#! format: off

"""
    mutable struct EXAC2 <: AVR
        Tr::Float64
        Tb::Float64
        Tc::Float64
        Ka::Float64
        Ta::Float64
        Va_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        Kb::Float64
        Vr_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        Te::Float64
        Kl::Float64
        Kh::Float64
        Kf::Float64
        Tf::Float64
        Kc::Float64
        Kd::Float64
        Ke::Float64
        V_lr::Float64
        E_sat::Tuple{Float64, Float64}
        Se::Tuple{Float64, Float64}
        V_ref::Float64
        saturation_coeffs::Tuple{Float64, Float64}
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
        states_types::Vector{StateTypes}
        internal::InfrastructureSystemsInternal
    end

Modified AC2. This excitation systems consists of an alternator main exciter feeding its output via non-controlled rectifiers.
The exciter does not employ self-excitation, and the voltage regulator power is taken from a source that is not affected by external transients.
Parameters of IEEE Std 421.5 Type AC2A Excitacion System. The alternator main exciter is used, feeding its output via non-controlled rectifiers. The Type AC2C model is similar to that of Type AC1C except for the inclusion of exciter time constant compensation and exciter field current limiting elements. EXAC2 in PSSE and PSLF.

# Arguments
- `Tr::Float64`: Regulator input filter time constant in s, validation range: `(0, 0.5)`, action if invalid: `warn`
- `Tb::Float64`: Regulator denominator (lag) time constant in s, validation range: `(0, 20)`, action if invalid: `warn`
- `Tc::Float64`: Regulator numerator (lead) time constant in s, validation range: `(0, 20)`, action if invalid: `warn`
- `Ka::Float64`: Regulator output gain, validation range: `(0, 1000)`
- `Ta::Float64`: Regulator output time constant in s, validation range: `(0, 10)`, action if invalid: `warn`
- `Va_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Limits for regulator output `(Va_min, Va_max)`
- `Kb::Float64`: Second Stage regulator gain, validation range: `(eps(), 500)`, action if invalid: `error`
- `Vr_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Limits for exciter field voltage `(Vr_min, Vr_max)`
- `Te::Float64`: Exciter field time constant, validation range: `(eps(), 2)`, action if invalid: `error`
- `Kl::Float64`: Exciter field current limiter gain, validation range: `(0, 1.1)`, action if invalid: `warn`
- `Kh::Float64`: Exciter field current regulator feedback gain, validation range: `(0, 1.1)`, action if invalid: `warn`
- `Kf::Float64`: Rate feedback excitation system stabilizer gain, validation range: `(0, 0.3)`, action if invalid: `warn`
- `Tf::Float64`: Rate feedback time constant, validation range: `(eps(), nothing)`, action if invalid: `error`
- `Kc::Float64`: Rectifier loading factor proportional to commutating reactance, validation range: `(0, 1)`, action if invalid: `warn`
- `Kd::Float64`: Demagnetizing factor, function of exciter alternator reactances, validation range: `(0, 1)`, action if invalid: `warn`
- `Ke::Float64`: Exciter field proportional constant, validation range: `(0, 1)`, action if invalid: `warn`
- `V_lr::Float64`: Maximum exciter field current, validation range: `(0, nothing)`, action if invalid: `warn`
- `E_sat::Tuple{Float64, Float64}`: Exciter output voltage for saturation factor: (E1, E2)
- `Se::Tuple{Float64, Float64}`: Exciter saturation factor at exciter output voltage: (Se(E1), Se(E2))
- `V_ref::Float64`: Reference Voltage Set-point, validation range: `(0, nothing)`
- `saturation_coeffs::Tuple{Float64, Float64}`: Coefficients (A,B) of the function: Se(V) = B(V - A)^2/V
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: The states are:
	Vm: Sensed terminal voltage,
	Vr1: Lead-lag state,
	Vr2: Regulator output state,
	Ve: Integrator output state,
	Vr3: Feedback output state
- `n_states::Int`: EXAC2 has 5 states
- `states_types::Vector{StateTypes}`: EXAC2 has 5 states
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct EXAC2 <: AVR
    "Regulator input filter time constant in s"
    Tr::Float64
    "Regulator denominator (lag) time constant in s"
    Tb::Float64
    "Regulator numerator (lead) time constant in s"
    Tc::Float64
    "Regulator output gain"
    Ka::Float64
    "Regulator output time constant in s"
    Ta::Float64
    "Limits for regulator output `(Va_min, Va_max)`"
    Va_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Second Stage regulator gain"
    Kb::Float64
    "Limits for exciter field voltage `(Vr_min, Vr_max)`"
    Vr_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Exciter field time constant"
    Te::Float64
    "Exciter field current limiter gain"
    Kl::Float64
    "Exciter field current regulator feedback gain"
    Kh::Float64
    "Rate feedback excitation system stabilizer gain"
    Kf::Float64
    "Rate feedback time constant"
    Tf::Float64
    "Rectifier loading factor proportional to commutating reactance"
    Kc::Float64
    "Demagnetizing factor, function of exciter alternator reactances"
    Kd::Float64
    "Exciter field proportional constant"
    Ke::Float64
    "Maximum exciter field current"
    V_lr::Float64
    "Exciter output voltage for saturation factor: (E1, E2)"
    E_sat::Tuple{Float64, Float64}
    "Exciter saturation factor at exciter output voltage: (Se(E1), Se(E2))"
    Se::Tuple{Float64, Float64}
    "Reference Voltage Set-point"
    V_ref::Float64
    "Coefficients (A,B) of the function: Se(V) = B(V - A)^2/V"
    saturation_coeffs::Tuple{Float64, Float64}
    ext::Dict{String, Any}
    "The states are:
	Vm: Sensed terminal voltage,
	Vr1: Lead-lag state,
	Vr2: Regulator output state,
	Ve: Integrator output state,
	Vr3: Feedback output state"
    states::Vector{Symbol}
    "EXAC2 has 5 states"
    n_states::Int
    "EXAC2 has 5 states"
    states_types::Vector{StateTypes}
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function EXAC2(Tr, Tb, Tc, Ka, Ta, Va_lim, Kb, Vr_lim, Te, Kl, Kh, Kf, Tf, Kc, Kd, Ke, V_lr, E_sat, Se, V_ref=1.0, saturation_coeffs=PowerSystems.get_avr_saturation(E_sat, Se), ext=Dict{String, Any}(), )
    EXAC2(Tr, Tb, Tc, Ka, Ta, Va_lim, Kb, Vr_lim, Te, Kl, Kh, Kf, Tf, Kc, Kd, Ke, V_lr, E_sat, Se, V_ref, saturation_coeffs, ext, [:Vm, :Vr1, :Vr2, :Ve, :Vr3], 5, [StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Differential, StateTypes.Differential], InfrastructureSystemsInternal(), )
end

function EXAC2(; Tr, Tb, Tc, Ka, Ta, Va_lim, Kb, Vr_lim, Te, Kl, Kh, Kf, Tf, Kc, Kd, Ke, V_lr, E_sat, Se, V_ref=1.0, saturation_coeffs=PowerSystems.get_avr_saturation(E_sat, Se), ext=Dict{String, Any}(), states=[:Vm, :Vr1, :Vr2, :Ve, :Vr3], n_states=5, states_types=[StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Differential, StateTypes.Differential], internal=InfrastructureSystemsInternal(), )
    EXAC2(Tr, Tb, Tc, Ka, Ta, Va_lim, Kb, Vr_lim, Te, Kl, Kh, Kf, Tf, Kc, Kd, Ke, V_lr, E_sat, Se, V_ref, saturation_coeffs, ext, states, n_states, states_types, internal, )
end

# Constructor for demo purposes; non-functional.
function EXAC2(::Nothing)
    EXAC2(;
        Tr=0,
        Tb=0,
        Tc=0,
        Ka=0,
        Ta=0,
        Va_lim=(min=0.0, max=0.0),
        Kb=0,
        Vr_lim=(min=0.0, max=0.0),
        Te=0,
        Kl=0,
        Kh=0,
        Kf=0,
        Tf=0,
        Kc=0,
        Kd=0,
        Ke=0,
        V_lr=0,
        E_sat=(0.0, 0.0),
        Se=(0.0, 0.0),
        V_ref=0,
        saturation_coeffs=(0.0, 0.0),
        ext=Dict{String, Any}(),
    )
end

"""Get [`EXAC2`](@ref) `Tr`."""
get_Tr(value::EXAC2) = value.Tr
"""Get [`EXAC2`](@ref) `Tb`."""
get_Tb(value::EXAC2) = value.Tb
"""Get [`EXAC2`](@ref) `Tc`."""
get_Tc(value::EXAC2) = value.Tc
"""Get [`EXAC2`](@ref) `Ka`."""
get_Ka(value::EXAC2) = value.Ka
"""Get [`EXAC2`](@ref) `Ta`."""
get_Ta(value::EXAC2) = value.Ta
"""Get [`EXAC2`](@ref) `Va_lim`."""
get_Va_lim(value::EXAC2) = value.Va_lim
"""Get [`EXAC2`](@ref) `Kb`."""
get_Kb(value::EXAC2) = value.Kb
"""Get [`EXAC2`](@ref) `Vr_lim`."""
get_Vr_lim(value::EXAC2) = value.Vr_lim
"""Get [`EXAC2`](@ref) `Te`."""
get_Te(value::EXAC2) = value.Te
"""Get [`EXAC2`](@ref) `Kl`."""
get_Kl(value::EXAC2) = value.Kl
"""Get [`EXAC2`](@ref) `Kh`."""
get_Kh(value::EXAC2) = value.Kh
"""Get [`EXAC2`](@ref) `Kf`."""
get_Kf(value::EXAC2) = value.Kf
"""Get [`EXAC2`](@ref) `Tf`."""
get_Tf(value::EXAC2) = value.Tf
"""Get [`EXAC2`](@ref) `Kc`."""
get_Kc(value::EXAC2) = value.Kc
"""Get [`EXAC2`](@ref) `Kd`."""
get_Kd(value::EXAC2) = value.Kd
"""Get [`EXAC2`](@ref) `Ke`."""
get_Ke(value::EXAC2) = value.Ke
"""Get [`EXAC2`](@ref) `V_lr`."""
get_V_lr(value::EXAC2) = value.V_lr
"""Get [`EXAC2`](@ref) `E_sat`."""
get_E_sat(value::EXAC2) = value.E_sat
"""Get [`EXAC2`](@ref) `Se`."""
get_Se(value::EXAC2) = value.Se
"""Get [`EXAC2`](@ref) `V_ref`."""
get_V_ref(value::EXAC2) = value.V_ref
"""Get [`EXAC2`](@ref) `saturation_coeffs`."""
get_saturation_coeffs(value::EXAC2) = value.saturation_coeffs
"""Get [`EXAC2`](@ref) `ext`."""
get_ext(value::EXAC2) = value.ext
"""Get [`EXAC2`](@ref) `states`."""
get_states(value::EXAC2) = value.states
"""Get [`EXAC2`](@ref) `n_states`."""
get_n_states(value::EXAC2) = value.n_states
"""Get [`EXAC2`](@ref) `states_types`."""
get_states_types(value::EXAC2) = value.states_types
"""Get [`EXAC2`](@ref) `internal`."""
get_internal(value::EXAC2) = value.internal

"""Set [`EXAC2`](@ref) `Tr`."""
set_Tr!(value::EXAC2, val) = value.Tr = val
"""Set [`EXAC2`](@ref) `Tb`."""
set_Tb!(value::EXAC2, val) = value.Tb = val
"""Set [`EXAC2`](@ref) `Tc`."""
set_Tc!(value::EXAC2, val) = value.Tc = val
"""Set [`EXAC2`](@ref) `Ka`."""
set_Ka!(value::EXAC2, val) = value.Ka = val
"""Set [`EXAC2`](@ref) `Ta`."""
set_Ta!(value::EXAC2, val) = value.Ta = val
"""Set [`EXAC2`](@ref) `Va_lim`."""
set_Va_lim!(value::EXAC2, val) = value.Va_lim = val
"""Set [`EXAC2`](@ref) `Kb`."""
set_Kb!(value::EXAC2, val) = value.Kb = val
"""Set [`EXAC2`](@ref) `Vr_lim`."""
set_Vr_lim!(value::EXAC2, val) = value.Vr_lim = val
"""Set [`EXAC2`](@ref) `Te`."""
set_Te!(value::EXAC2, val) = value.Te = val
"""Set [`EXAC2`](@ref) `Kl`."""
set_Kl!(value::EXAC2, val) = value.Kl = val
"""Set [`EXAC2`](@ref) `Kh`."""
set_Kh!(value::EXAC2, val) = value.Kh = val
"""Set [`EXAC2`](@ref) `Kf`."""
set_Kf!(value::EXAC2, val) = value.Kf = val
"""Set [`EXAC2`](@ref) `Tf`."""
set_Tf!(value::EXAC2, val) = value.Tf = val
"""Set [`EXAC2`](@ref) `Kc`."""
set_Kc!(value::EXAC2, val) = value.Kc = val
"""Set [`EXAC2`](@ref) `Kd`."""
set_Kd!(value::EXAC2, val) = value.Kd = val
"""Set [`EXAC2`](@ref) `Ke`."""
set_Ke!(value::EXAC2, val) = value.Ke = val
"""Set [`EXAC2`](@ref) `V_lr`."""
set_V_lr!(value::EXAC2, val) = value.V_lr = val
"""Set [`EXAC2`](@ref) `E_sat`."""
set_E_sat!(value::EXAC2, val) = value.E_sat = val
"""Set [`EXAC2`](@ref) `Se`."""
set_Se!(value::EXAC2, val) = value.Se = val
"""Set [`EXAC2`](@ref) `V_ref`."""
set_V_ref!(value::EXAC2, val) = value.V_ref = val
"""Set [`EXAC2`](@ref) `saturation_coeffs`."""
set_saturation_coeffs!(value::EXAC2, val) = value.saturation_coeffs = val
"""Set [`EXAC2`](@ref) `ext`."""
set_ext!(value::EXAC2, val) = value.ext = val
"""Set [`EXAC2`](@ref) `states_types`."""
set_states_types!(value::EXAC2, val) = value.states_types = val


text:
using ElectrostaticPIC1D, Random, Test; Random.seed!(0)

using ForwardDiff, QuadGK

@testset "BasisFunctions" begin

numtests = 10

function isnormalisedtest(Shape)
  outcome = true
  for _ ∈ 1:numtests
    w = rand()
    s = Shape(w)
    x0 = randn()
    integral = quadgk(x->s(x, x0), x0 - 10*w, x0 + 10*w, order=27, rtol=eps())[1]
    outcome = integral ≈ 1
    outcome || break
  end
  @test outcome
end
function supporttest(Shape)
  for _ ∈ 1:10
    w = rand()
    x0 = randn()
    b = BasisFunction(Shape(w), x0)
    outcome = true
    for _ in 1:1000
      x = 3 * randn() * w + x0
      if x ∈ b
        outcome = b(x) > 0
        outcome || @show b, x, b(x)
      else
        outcome = b(x) <= b(x0) * eps()
        outcome || @show b, x, b(x)
      end
      outcome || break
    end
    @test outcome
  end
end
function energytest(Shape)
  outcome = true
  for _ ∈ 1:numtests
    w = rand()
    x0 = randn()
    b = BasisFunction(Shape(w), x0)
    expected = quadgk(x->b(x)^2, x0 - 10*w, x0 + 10*w, order=27, rtol=eps())[1] * weight(b)^2 / 2
    result = integral(b, b, PeriodicGridBC(x0-10*w, x0+10*w)) * weight(b)^2 / 2
    outcome = expected ≈ result
    outcome || break
  end
  @test outcome
end

for S ∈ (BSpline{0}, BSpline{1}, BSpline{2}, GaussianShape)
 @testset "$S" begin
   isnormalisedtest(S)
   supporttest(S)
   energytest(S)
 end
end

@testset "Pairwise Integrals" begin
for U ∈ (DeltaFunctionShape, BSpline{0}, BSpline{1}, BSpline{2}, GaussianShape)
for V ∈ (BSpline{0}, BSpline{1}, BSpline{2}, GaussianShape)
@testset "$U $V" begin
  noverlaps = 0
  while true
    v = BasisFunction(V(5*rand()), randn())
    u = if U <: DeltaFunctionShape
      BasisFunction(U(), lower(v) + rand() * width(v))
    else
      BasisFunction(U(5*rand()), randn())
    end
    expected = if U <: DeltaFunctionShape
      v(centre(u))
    else
      quadgk(x-> u(x) * v(x), lower(u, v), upper(u, v);
             order=51, atol=0, rtol=eps())[1]
    end
    expected < eps() && continue
    try
      @inferred integral(u, v)
      @test true
    catch
      @test false
    end
    try
      @inferred integral(u, v, PeriodicGridBC(1.0))
      @test true
    catch
      @warn "integral($U, $V) cannot be inferred"
      @test false
    end
    result = integral(u, v)
    @test result ≈ expected atol=100eps()
    break
  end
end # testset U, V
end # V
end # U
end # testset integrals

end


text:
using Colors

include("helpers/page.jl")
include("helpers/doc.jl")

intro = md"""
$(title(2, "Higher Order Layouts"))

$(vskip(1em))
# Functions
$(vskip(1em))

$(
    showdocs([ tabs, pages, menu, submenu, dropdownmenu, icon,
                iconbutton, toolbar])
)
"""

function main(window)
    docpage(intro)
end


text:
struct TsptwTrajectoryState <: GraphTrajectoryState
    fg::FeaturedGraph
    variableIdx::Int
    possibleValuesIdx::AbstractVector{Int}
end

Flux.@functor TsptwTrajectoryState


text:
export EFTfitterModel
export to_correlation_matrix

export get_observables, get_parameters, get_measurements, get_correlations
export get_nuisance_correlations, get_total_covariance, get_covariances, get_measurement_distributions

"""
    struct EFTfitterModel

This is the central type for using EFTfitter.
It comprises all information necessary for performing an analysis.
Only active `Measurement` and `Correlation` objects will be considered.

Fields:
* `parameters::BAT.NamedTupleDist`
* `measurements::NamedTuple{<:Any, <:Tuple{Vararg{Measurement}}}`
* `measurementdistributions::NamedTuple{<:Any, <:Tuple{Vararg{MeasurementDistribution}}}`
* `correlations::NamedTuple{<:Any, <:Tuple{Vararg{Correlation}}}`
* `nuisances::Union{NamedTuple{<:Any, <:Tuple{Vararg{NuisanceCorrelation}}}, Nothing}`

Constructors:

```julia
EFTfitterModel(
    parameters::BAT.NamedTupleDist,
    measurements::NamedTuple{<:Any, <:Tuple{Vararg{AbstractMeasurement}}},
    correlations::NamedTuple{<:Any, <:Tuple{Vararg{AbstractCorrelation}}},
    nuisances::Union{NamedTuple{<:Any, <:Tuple{Vararg{NuisanceCorrelation}}}, Nothing} = nothing
)
```

Examples:

```julia
model = EFTfitterModel(parameters, measurements, correlations) # no nuisance correlations

model = EFTfitterModel(parameters, measurements, correlations, nuisances) # with nuisance correlations
)
```
"""
struct EFTfitterModel
    parameters::BAT.NamedTupleDist
    measurements::NamedTuple{<:Any, <:Tuple{Vararg{Measurement}}}
    measurementdistributions::NamedTuple{<:Any, <:Tuple{Vararg{MeasurementDistribution}}}
    correlations::NamedTuple{<:Any, <:Tuple{Vararg{Correlation}}}
    nuisances::Union{NamedTuple{<:Any, <:Tuple{Vararg{NuisanceCorrelation}}}, Nothing}
end

function EFTfitterModel(
    parameters::BAT.NamedTupleDist,
    measurements::NamedTuple{<:Any, <:Tuple{Vararg{AbstractMeasurement}}},
    correlations::NamedTuple{<:Any, <:Tuple{Vararg{AbstractCorrelation}}},
    nuisances::Union{NamedTuple{<:Any, <:Tuple{Vararg{NuisanceCorrelation}}}, Nothing} = nothing
)
    measurement_vec, measurement_keys = unpack(measurements)
    correlation_vec, uncertainty_keys = unpack(correlations)

    # convert elements of MeasurementDistribution to Measurement for each bin
    binned_measurements, binned_measurement_keys = convert_to_bins(measurement_vec, measurement_keys)
    # use only active measurements/bins
    active_measurements, active_measurement_keys, corrs = only_active_measurements(binned_measurements, binned_measurement_keys, correlation_vec)
    # use only active uncertainties and correlations
    active_measurements, active_correlations, uncertainty_keys = only_active_uncertainties(active_measurements, corrs, uncertainty_keys)

    correlation_nt = namedtuple(uncertainty_keys, active_correlations)
    measurement_nt = namedtuple(active_measurement_keys, active_measurements)
    meas_dists_nt  = create_distributions(measurements, uncertainty_keys)
    nuisances_nt   = only_active_nuisances(nuisances, active_measurement_keys, uncertainty_keys)

    params = add_nuisance_parameters(parameters, nuisances_nt)

    return EFTfitterModel(params, measurement_nt, meas_dists_nt, correlation_nt, nuisances_nt)
end
#=============================================================#
"""
    get_parameters(m::EFTfitterModel)

Returns model parameters.
"""
get_parameters(m::EFTfitterModel) = m.parameters

"""
    get_measurements(m::EFTfitterModel)

Returns a `NamedTuple` with the `Measurement`s in the `EFTfitterModel`.
"""
get_measurements(m::EFTfitterModel) = m.measurements

"""
    get_measurement_distributions(m::EFTfitterModel)

Returns a `NamedTuple` with the `MeasurementDistribution`s in the `EFTfitterModel`.
"""
get_measurement_distributions(m::EFTfitterModel) = m.measurementdistributions

"""
    get_correlations(m::EFTfitterModel)

Returns a `NamedTuple` with the `Correlations`s in the `EFTfitterModel`.
"""
get_correlations(m::EFTfitterModel) = m.correlations

"""
    get_nuisance_correlations(m::EFTfitterModel)

Returns a `NamedTuple` with the `NuisanceCorrelation`s in the `EFTfitterModel`.
"""
get_nuisance_correlations(m::EFTfitterModel) = m.nuisances

"""
    get_observables(m::EFTfitterModel)

Returns a `NamedTuple` with the `Observable`s in the `EFTfitterModel`.
Note: The upper and lower limits are ignored and for each unique `Function`s only one `Observable` is returned.
"""
function  get_observables(model::EFTfitterModel)
    meas = get_measurements(model)
    obs = unique(Observable.([m.observable.func for m in values(meas)]))
    obs_names = [string(o.func) for o in obs]
    observables_nt = namedtuple(obs_names, obs)
end

#========================================================#

function unpack(nt::NamedTuple{N,T}) where {N,T}
     return collect(nt), collect(keys(nt))
end

function convert_to_bins(
    measurements::Array{<:AbstractMeasurement, 1},
    measurement_keys::Array{Symbol, 1}
)
    binned_measurements = reduce(vcat, convert_to_bins.(measurements))
    binned_meas_keys = reduce(vcat, keys_of_bins.(measurements, measurement_keys))

    if isa(binned_measurements, Measurement) # when only one measurement is active
        return [binned_measurements], [binned_meas_keys]
    end

    return binned_measurements, binned_meas_keys
end

function convert_to_bins(m::Measurement)
    return m
end

function convert_to_bins(md::MeasurementDistribution)
    nbins = length(md.value)
    uncertainties = [[u[i] for u in md.uncertainties] for i in 1:nbins]

    unc_nt = namedtuple.(Ref(keys(md.uncertainties)), uncertainties)
    return Measurement.(md.observable, md.value, unc_nt, md.active)
end

function only_active_measurements(
    measurements::Array{Measurement, 1},
    measurement_keys::Array{Symbol, 1},
    correlations::Array{<:AbstractCorrelation, 1}
)
    active_idxs = filter(x->measurements[x].active, 1:length(measurements))
    measurements = measurements[active_idxs]
    measurement_keys = measurement_keys[active_idxs]

    corrs = get_correlations.(correlations, Ref(active_idxs))

    return measurements, measurement_keys, corrs
end

#TODO: rename?
function get_correlations(
    correlation::Correlation,
    idxs::Array{<:Integer, 1}
)
    n=length(idxs)
    # ignore size of given correlation matrix if active==false
    if correlation.active == false
        return Correlation(Symmetric(zeros(n, n)), false)
    end

    return Correlation(Symmetric(correlation.matrix[idxs, idxs]), correlation.active)
end

function get_correlations(
    correlation::NoCorrelation,
    idxs::Array{<:Integer, 1}
)
    n = length(idxs)
    return Correlation(Symmetric(Matrix{Float64}(I, n, n)), correlation.active)
end

#--------------------
function create_distributions(
    m::NamedTuple{<:Any, <:Tuple{Vararg{AbstractMeasurement}}},
    uncertainty_keys::Vector{Symbol}
)
    active_idxs = [i for i in 1:length(m) if (isa(m[i], MeasurementDistribution) && any(m[i].active))]
    dists = [only_active_bins(m[i], uncertainty_keys) for i in active_idxs]
    dists_keys = [keys(m)[i] for i in active_idxs]

    if isempty(dists)
        return NamedTuple()
    else
        return namedtuple(dists_keys, dists)
    end
end

function only_active_bins(md::MeasurementDistribution, uncertainty_keys::Vector{Symbol})
    obs = md.observable[md.active]
    vals = md.value[md.active]

    unc = [md.uncertainties[u][md.active] for u in uncertainty_keys]
    uncertainties = namedtuple(uncertainty_keys, unc)
    bin_names = md.bin_names[md.active]

    return MeasurementDistribution(obs, vals, uncertainties=uncertainties, bin_names=bin_names)
end

function keys_of_bins(m::Measurement, key::Symbol)
    return key
end

function keys_of_bins(md::MeasurementDistribution, key::Symbol)
    return [Symbol(String(key)*"_"*String(b)) for b in md.bin_names]
end

function only_active_uncertainties(
    measurements::Array{Measurement, 1},
    correlations::Array{<:AbstractCorrelation, 1},
    uncertainty_keys::Array{Symbol, 1}
)
    active_idxs = filter(x->correlations[x].active, 1:length(correlations))
    active_correlations = correlations[active_idxs]
    uncertainty_keys = uncertainty_keys[active_idxs]

    measurements = [@set m.uncertainties = select(m.uncertainties, uncertainty_keys) for m in measurements]

    return measurements, active_correlations, uncertainty_keys
end

function only_active_nuisances(
    nuisances::Nothing,
    active_meas_keys::Array{Symbol, 1},
    unc_keys::Array{Symbol, 1}
)
    return nothing
end

function only_active_nuisances(
    nuisances::NamedTuple{<:Any, <:Tuple{Vararg{NuisanceCorrelation}}},
    active_meas_keys::Array{Symbol, 1},
    unc_keys::Array{Symbol, 1}
)
    nuisance_vec, nuisance_keys = unpack(nuisances)

    active_nuisances = NuisanceCorrelation[]
    active_nuisance_keys = Symbol[]

    for (nui, nui_k) in zip(nuisance_vec, nuisance_keys)
        if check_nuisance(nui, active_meas_keys, unc_keys)
            push!(active_nuisances, nui)
            push!(active_nuisance_keys, nui_k)
        end
    end

    if length(active_nuisances) > 0
        return namedtuple(active_nuisance_keys, active_nuisances)
    else
        return nothing
    end
end

function check_nuisance(
    nuisance::NuisanceCorrelation,
    active_meas_keys::Array{Symbol, 1},
    unc_keys::Array{Symbol, 1}
)
    if nuisance.meas1 == nuisance.meas2
        return false
    elseif !(nuisance.unc_key in unc_keys)
        return false
    elseif !(nuisance.meas1 in active_meas_keys)
        return false
    elseif !(nuisance.meas2 in active_meas_keys)
        return false
    end

    return true
end

function get_total_covariance(m::EFTfitterModel)
    covs = get_covariances(m)
    total_cov = Symmetric(sum(covs))

    return total_cov
end

function get_covariances(m::EFTfitterModel)
    unc_values = [[Float64(meas.uncertainties[u]) for meas in m.measurements] for u in keys(m.correlations)]
    corrs = [c.matrix for c in m.correlations]

    covs = [Symmetric(σ*ρ*σ) for (σ, ρ) in zip(diagm.(unc_values), corrs)]
end

function Nbins(measurements::NamedTuple{<:Any, <:Tuple{Vararg{AbstractMeasurement}}})
    sum([length(m.value) for m in measurements])
end

function get_indices(measurements, symb)
    meas_vec, meas_keys = unpack(measurements)
    names = keys_of_bins(meas_vec, meas_keys)
    ind = [findfirst(x -> x==symb, names)]

    if ind[1] != nothing
        return ind
    else
        ind = [findfirst(x -> x==Symbol(String(symb)*"_bin$i"), names) for i in 1:length(measurements[symb].value)]
    end
end

"""
    to_correlation_matrix(
        measurements::NamedTuple{<:Any, <:Tuple{Vararg{AbstractMeasurement}}},
        correlations::Tuple{Symbol,Symbol, Union{<:Real, Array{<:Real, 2}}}...
    )

    Returns a `Matrix{Float64}` as the correlation matrix for the measurements with all diagonal-elements being unity.
    The correlations are specified by passing tuples of two `Symbols` (the keys of the `Measurements`) with a value or matrix for the correlation coefficients.
    If the matrix is not positive-definite, a warning is shown but the matrix is still returned.

    Example:
    ```julia
    dist_corr = [1.0 0.5 0.0;
                 0.5 1.0 0.0;
                 0.0 0.0 1.0]

    another_corr_matrix = to_correlation_matrix(
        measurements,
        (:Meas1, :Meas2, 0.4),
        (:Meas1, :MeasDist, 0.1),
        (:MeasDist, :MeasDist, dist_corr),
        (:MeasDist_bin2, :MeasDist_bin3, 0.3),
    )
    ```
"""
function to_correlation_matrix(
    measurements::NamedTuple{<:Any, <:Tuple{Vararg{AbstractMeasurement}}},
    correlations::Any...#Tuple{Symbol,Symbol, Union{<:Real, Array{<:Real, 2}}}... #TODO: typing
)
    nmeas = Nbins(measurements)
    corr = Matrix{Float64}(I, nmeas, nmeas)

    for c in correlations #TODO: split into functions
        if isa(c, Tuple{Symbol,Symbol, Union{<:Real, Array{<:Real, 2}}})
            i = get_indices(measurements, c[1])
            j = get_indices(measurements, c[2])

            corr[i, j] .= c[3]
            corr[j, i] .= c[3]

        elseif isa(c, Tuple{Any, Union{<:Real, Array{<:Real, 2}}}) #TODO: add testing
            i = reduce(vcat, get_indices.(Ref(measurements), c[1]))
            corr[i, i] .= c[2]
        end
    end

    corr[diagind(corr)] .= 1

    if !isposdef(corr)
        @warn "The correlation matrix $corr is not positive semidefinite"
    end

    return corr
end

function has_nuisance_correlations(m::EFTfitterModel)
    m.nuisances == nothing ? (return false) : (return true)
end

function add_nuisance_parameters(
    parameters::NamedTupleDist,
    nuisances_nt::NamedTuple{<:Any, <:Tuple{Vararg{NuisanceCorrelation}}}
)
    param_keys = [keys(parameters)...,  keys(nuisances_nt)...]

    nuisances = [nui.prior for nui in nuisances_nt]
    param_dists = [values(parameters)...,  nuisances...]
    nt = namedtuple(param_keys, param_dists)

    return BAT.NamedTupleDist(nt)
end

function add_nuisance_parameters(
    parameters::NamedTupleDist,
    nuisances_nt::Nothing
)
    return parameters
end

# Todo: rename?
function keys_of_bins(
    measurements::Array{T, 1},
    measurement_keys::Array{Symbol, 1}
) where T <: AbstractMeasurement
    keys = reduce(vcat, keys_of_bins.(measurements, measurement_keys))
    return keys
end


text:
using GLMakie
#input format: scattering cross section magnitude (R), emiss cross section ((x,y) point_list)
function main()
    granularity = 20

    function convert_to_spherical(point_list)
        spherical_list = []
        for row in 1:size(point_list,1)
            point_x = point_list[row, 1]
            point_y = point_list[row, 2]
            point_z = point_list[row, 3]
            r = point_x^2 + point_y^2 + point_z^2

            if point_y == 0
                θ = 0
            else
                θ = atan(point_x/point_y)
            end

            if point_y == 0
                ϕ = 0
            else
                ϕ = acos(point_z/r)
            end

            if row == 1
                spherical_list = [r θ ϕ]
            else
                spherical_list = vcat(spherical_list, [r θ ϕ])
            end
        end
        return(spherical_list)
    end

    function draw_mesh(points, faces)
        scene = mesh(points, faces, shading = false)
    end

    function make_mesh(granularity::Int, point_list)
        num_points = size(point_list,1)
        new_point_list = []
        for row in 1:num_points
            for arc in 0:granularity-1
                theta = arc*pi/granularity
                point_x = cos(theta)*point_list[row,1]
                point_y = point_list[row,2]
                point_z = sin(theta)*point_list[row,1]
                if new_point_list == []
                    new_point_list = [point_x point_y point_z]
                else
                    new_point_list = vcat(new_point_list, [point_x point_y point_z])
                end
            end
        end

        face_list = []
        for v in 0:(num_points-1)
            for u in 0:(granularity-1)
                p1 = num_points*v+u+1
                p2 = mod(num_points*v + mod(u + 1, granularity), num_points*granularity)+1
                p3 = mod(num_points*(v + 1) + u, num_points*granularity)+1
                p4 = mod(num_points*(v + 1) + mod(u + 1, granularity), num_points*granularity)+1
                if face_list == []
                    face_list =
                        [
                            p1 p2 p3 p4
                        ]

                else
                    face_list = vcat(
                        face_list,
                        [
                            p1 p2 p3 p4
                        ]
                    )
                end
            end
        end

        return (new_point_list, face_list)
    end

    function calculate_T(
        face_list::Vector{Any}, #make this better typed
        point_list::Vector{Vector{Float64}}
    )

        spherical_list = convert_to_spherical(point_list) #format is [r, θ, ϕ], use coordinatetransformations.jl to convert

        #TODO: bundle these into a new struct and make that not break
        r_array, θ_array, ϕ_array = ([element[1] for element in spherical_list], [element[2] for element in spherical_list], [element[3] for element in spherical_list])

        # calculate r and n̂ for the geometry

        n̂_array = point_list #copying so that size is same

        for element in n̂_array
            element = [0,0,0]
        end

        for face in face_list
            for (vertex, i) in enumerate(face)
                n_array[vertex] += vertex × point_list[face[(i+1)%size(face)]] #TODO this should be edges not vertices
            end
        end

        for element in n̂_array
            element = normalize(element) #TODO make sure this works, use normalize!(element) later if at all
        end

        T = sum(n̂_array)
        return T
    end

    for i in 1:1

        point_list = [1 1 1; 1 8 1]
        for i in 1:2
            point_list = vcat(point_list,([rand((2:8)) rand((0:8)) rand((2:8))]))
        end

        mesh_point_list, mesh_face_list = make_mesh(granularity, point_list)
        scene = mesh(mesh_point_list, mesh_face_list, shading = false)

        println(mesh_face_list)
        println(mesh_point_list)

    end
end

main()


text:
@testset "240.search-a-2d-matrix-ii.jl" begin
    arr = [
        1 4 7 11 15
        2 5 8 12 19
        3 6 9 16 22
        10 13 14 17 24
        18 21 23 26 30
    ]

    @test search_matrix(arr, 5) == true
    @test search_matrix(arr, 20) == false
end


text:
type Liver # lumped portal vein, sinusoid, hepatic veins
    C::Vector{Float64}
    R::Vector{Float64}
    L::Float64
    V0::Vector{Float64}
    P::Any
    V::Any
    Q::Any

    function Liver()
        this = new()
        this.R = [0.004*mmHgToPa/cm3Tom3];
        append!(this.R,[0.005,0.005,0.004]*mmHgToPa/cm3Tom3);
        this.C = [3*cm3Tom3/mmHgToPa];
        append!(this.C,[10,15,45]*cm3Tom3/mmHgToPa);
        this.L = 5e-5*mmHgToPa/cm3Tom3;
        this.V0 = [50*cm3Tom3];
        append!(this.V0,[30,53,75]*cm3Tom3);
        this.P = Array{Float64,2}[];
        this.V = Array{Float64,2}[];
        this.Q = Array{Float64,2}[];
        return this
    end
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
# Blocked axes

function append_ranges(ranges::AbstractVector{<:AbstractUnitRange})
  blockedrange(map(length,ranges))
end

function append_ranges(ranges::AbstractVector{<:BlockedUnitRange})
  MultiLevelBlockedUnitRange(ranges)
end

"""
Implementation of a multilevel blocked range
It is a generalization of BlockedUnitRange that stores more information
about the local ranges.
"""
struct MultiLevelBlockedUnitRange{L} <: AbstractUnitRange{Int}
  local_ranges::Vector{L}
  global_range::BlockedUnitRange{Vector{Int}}
  function MultiLevelBlockedUnitRange(local_ranges::Vector{L}) where L<:AbstractUnitRange
    local_lengths = map(length,local_ranges)
    global_range = blockedrange(local_lengths)
    new{L}(local_ranges,global_range)
  end
end

Base.first(a::MultiLevelBlockedUnitRange) = first(a.global_range)
Base.last(a::MultiLevelBlockedUnitRange) = last(a.global_range)
BlockArrays.blockaxes(a::MultiLevelBlockedUnitRange) = blockaxes(a.global_range)
BlockArrays.blocklasts(a::MultiLevelBlockedUnitRange) = blocklasts(a.global_range)
BlockArrays.findblock(a::MultiLevelBlockedUnitRange,k::Integer) = findblock(a.global_range,k)
Base.getindex(a::MultiLevelBlockedUnitRange,i::Integer) = a.global_range[i]
Base.getindex(a::MultiLevelBlockedUnitRange,i::Block{1}) = a.global_range[i]
Base.getindex(a::MultiLevelBlockedUnitRange, i::BlockRange{1}) = a.global_range[i]
Base.axes(a::MultiLevelBlockedUnitRange) = axes(a.global_range)
Base.Broadcast.axistype(a::T, b::T) where T<:MultiLevelBlockedUnitRange =  Base.Broadcast.axistype(a.global_range,b.global_range)
Base.Broadcast.axistype(a::MultiLevelBlockedUnitRange, b::MultiLevelBlockedUnitRange) = Base.Broadcast.axistype(a.global_range,b.global_range)
Base.Broadcast.axistype(a::MultiLevelBlockedUnitRange, b) = Base.Broadcast.axistype(a.global_range,b)
Base.Broadcast.axistype(a, b::MultiLevelBlockedUnitRange) =Base.Broadcast.axistype(a,b.global_range)
#Base.print_matrix_row(
#  io::IO,
#  X::MultiLevelBlockedUnitRange,
#  A::Vector,
#  i::Integer,
#  cols::AbstractVector,
#  sep::AbstractString) = print_matrix_row(io,X.global_range,A,i,cols,sep)

local_range(a::AbstractUnitRange,k::Integer) = Base.OneTo(length(a[Block(k)]))
local_range(a::MultiLevelBlockedUnitRange,k::Integer) = a.local_ranges[k]
function Base.similar(
  t::Type{Array{T,N}},shape::NTuple{N,MultiLevelBlockedUnitRange}) where {T,N}
  sizes = map(length,shape)
  similar(t,sizes)
end

function similar_range(r::Base.OneTo,n::Integer)
  Base.OneTo(Int(n))
end

function similar_range(r::BlockedUnitRange,n::Integer)
  blockedrange([n])
end

function similar_range(r::MultiLevelBlockedUnitRange,n::Integer)
  r = similar_range(first(r.local_ranges),n)
  append_ranges([r])
end

"""
Check if the full multi-level block structure is the same.
This is in contrast to BlockArrays.blockisequal that only checks
the top level block structure.
"""
@inline blocks_equal(a::AbstractUnitRange,b::AbstractUnitRange) = BlockArrays.blockisequal(a,b)
blocks_equal(a::AbstractUnitRange,b::MultiLevelBlockedUnitRange) = false
blocks_equal(a::MultiLevelBlockedUnitRange,b::AbstractUnitRange) = false
blocks_equal(a::Base.OneTo,b::Base.OneTo) = a==b
function blocks_equal(a::MultiLevelBlockedUnitRange,b::MultiLevelBlockedUnitRange)
  if a === b
    return true
  end
  r = blocks_equal(a.global_range,b.global_range)
  if r == false
    return false
  end
  la = length(a.local_ranges)
  lb = length(b.local_ranges)
  if la!=lb
    return false
  else
    for i in 1:la
      @inbounds ra = a.local_ranges[i]
      @inbounds rb = b.local_ranges[i]
      r = r && blocks_equal(ra,rb)
    end
    return r
  end
end

"""
Check if the number ob blocks in all levels is the same.
"""
@inline function num_blocks_equal(a::AbstractUnitRange,b::AbstractUnitRange)
  sa = BlockArrays.blocksize(a)
  sb = BlockArrays.blocksize(b)
  sa == sb
end
num_blocks_equal(a::AbstractUnitRange,b::MultiLevelBlockedUnitRange) = false
num_blocks_equal(a::MultiLevelBlockedUnitRange,b::AbstractUnitRange) = false
function num_blocks_equal(a::MultiLevelBlockedUnitRange,b::MultiLevelBlockedUnitRange)
  if a === b
    return true
  end
  if length(a.local_ranges) != length(b.local_ranges)
    return false
  end
  for i in 1:length(a.local_ranges)
    if ! num_blocks_equal(a.local_ranges[i],b.local_ranges[i])
      return false
    end
  end
  true
end

blocks_equal(a::Tuple,b::Tuple) = all(map(blocks_equal,a,b))
num_blocks_equal(a::Tuple,b::Tuple) = all(map(num_blocks_equal,a,b))

"""
Block array in which some of the blocks are known to be zero
"""
struct BlockArrayCoo{T,N,A,X} <: AbstractBlockArray{T,N}
  axes::X
  blockids::Vector{NTuple{N,Int}}
  blocks::Vector{A}
  ptrs::Array{Int,N}
  zero_blocks::Vector{A}
  function BlockArrayCoo(
    axes::NTuple{N},
    blockids::Vector{NTuple{N,Int}},
    blocks::Vector{A},
    ptrs::Array{Int,N},
    zero_blocks::Vector{A}) where {T,N,A<:AbstractArray{T,N}}

    @assert length(blockids) == length(blocks)
    @check length(unique(blockids)) == length(blockids) "We cannot built a BlockArrayCoo from repeated blocks"
    @check _valid_block_sizes(axes,ptrs,blocks) "The given blocks do not match with the given axes"

    X = typeof(axes)
    new{T,N,A,X}(axes,blockids,blocks,ptrs,zero_blocks)
  end
end

const BlockMatrixCoo = BlockArrayCoo{T,2} where T
const BlockVectorCoo = BlockArrayCoo{T,1} where T

function _valid_block_sizes(_axes,ptrs,blocks)
  s = map(i->first(blocksize(i)),_axes)
  cis = CartesianIndices(s)
  for ci in cis
    p = ptrs[ci]
    if p>0
      I = Tuple(ci)
      laxs = map( local_range, _axes, I)
      if !blocks_equal(axes(blocks[p]),laxs)
        return false
      end
    end
  end
  true
end

# Minimal constructor

function BlockArrayCoo(
  axes::NTuple{N},
  blockids::Vector{NTuple{N,Int}},
  blocks::Vector{A}) where {T,N,A<:AbstractArray{T,N}}

  ptrs = _compute_ptrs(blockids,axes)
  zero_blocks = _compute_zero_blocks(A,ptrs,axes)
  BlockArrayCoo(axes,blockids,blocks,ptrs,zero_blocks)
end

function _compute_ptrs(blockids,axes)
  s = map(i->first(blocksize(i)),axes)
  ptrs = zeros(Int,s)
  for (i,c) in enumerate(blockids)
    ptrs[c...] = i
  end
  m = 1
  for (k,p) in enumerate(ptrs)
    if p ==0
      ptrs[k] = -m
      m += 1
    end
  end
  ptrs
end

function _compute_zero_blocks(::Type{A},ptrs,axs) where A
  cis = CartesianIndices(ptrs)
  zero_blocks = A[]
  for ci in cis
    p = ptrs[ci]
    if p<0
      I = Tuple(ci)
      laxs = map( local_range, axs, I)
      block = _zero_block(A,laxs)
      push!(zero_blocks,block)
    end
  end
  zero_blocks
end

function _zero_block(::Type{A},axs::Tuple) where A <: AbstractArray
  a = similar(A,axs)
  fill!(a,zero(eltype(a)))
  a
end

function _zero_block(::Type{<:Transpose{T,A}},axs::Tuple) where {T,A}
  a = similar(A,(axs[2],))
  fill!(a,zero(eltype(a)))
  Transpose(a)
end

function _zero_block(::Type{<:BlockArrayCoo{T,N,A}},axs::Tuple) where {T,N,A}
  blocks = A[]
  blockids = NTuple{N,Int}[]
  BlockArrayCoo(axs,blockids,blocks)
end

# Minimal constructor (for lazy_map)

struct BlockArrayCooMap{N} <: Map
  blocksize::NTuple{N,Int}
  blockids::Vector{NTuple{N,Int}}
  ptrs::Array{Int,N}
  function BlockArrayCooMap(blocksize::NTuple{N,Int}, blockids::Vector{NTuple{N,Int}}) where N
    ptrs = fill(-1,blocksize)
    for (p,I) in enumerate(blockids)
      ptrs[I...] = p
      for i in 1:N
        @check blocksize[i] >= I[i]
      end
    end
    new{N}(blocksize,blockids,ptrs)
  end
end

@inline function lazy_map(k::BlockArrayCooMap,T::Type,f::AbstractArray...)
  s = _common_size(f...)
  N = length(s)
  LazyArray(T,Val(N),Fill(k,s),f...)
end

#function BlockArrayCoo(
#  axes::NTuple{N},
#  blockids::Vector{NTuple{N,Int}},
#  blocks::A...) where {T,N,A<:AbstractArray{T,N}}
#
#  BlockArrayCoo(axes,blockids,collect(blocks))
#end

function return_cache(
  k::BlockArrayCooMap,
  axes::NTuple{N},
  blocks::A...) where {T,N,A<:AbstractArray{T,N}}

  @check map(i->first(blocksize(i)),axes) == k.blocksize "The given axes are not compatible with the given BlockArrayCooMap"

  if _valid_block_sizes(axes,_compute_ptrs(k.blockids,axes),blocks)
    goodblocks = collect(blocks)
  else
    goodblocks = A[]
    for (i,I) in enumerate(k.blockids)
      laxs = map( local_range, axes, I)
      block = similar(typeof(blocks[i]),laxs)
      push!(goodblocks,block)
    end
  end

  r = BlockArrayCoo(axes,k.blockids,goodblocks)
  CachedArray(r)
end

@inline function evaluate!(
  cache,
  ::BlockArrayCooMap,
  axes::NTuple{N},
  blocks::A...) where {T,N,A<:AbstractArray{T,N}}

  setaxes!(cache,axes)
  r = cache.array
  copyto!(r.blocks,blocks)
  r
end

BlockArrays.blocksize(a::BlockArrayCooMap) = a.blocksize
is_zero_block(a::BlockArrayCooMap,i::Integer) = a.ptrs[i] < 0
is_zero_block(a::BlockArrayCooMap{N},i::Vararg{Integer,N}) where N = a.ptrs[i...] < 0
is_zero_block(a::BlockArrayCooMap{N},i::Vararg{Block,N}) where N = is_zero_block(a,map(Int,i)...)
is_zero_block(a::BlockArrayCooMap,i::Block) = is_zero_block(a,convert(Tuple,i)...)
is_zero_block(a::BlockArrayCooMap,i::CartesianIndex) = is_zero_block(a,Tuple(i)...)

# Specific API

@inline is_nonzero_block(args...) = ! is_zero_block(args...)

function is_zero_block(a::AbstractArray,b::Block)
  i = convert(Tuple,b)
  is_zero_block(a,i...)
end

function is_zero_block(a::AbstractArray{T,N},b::Vararg{Block,N}) where {T,N}
  i = map(Int,b)
  is_zero_block(a,i...)
end

function is_zero_block(a::AbstractArray,b::CartesianIndex)
  i = Tuple(b)
  is_zero_block(a,i...)
end

function is_zero_block(a::BlockArrayCoo{T,N},i::Vararg{Integer,N}) where {T,N}
  #@boundscheck BlockArrays.blockcheckbounds(a, i...)
  p = a.ptrs[i...]
  @check p != 0
  p < 0
end

function is_zero_block(a::AbstractArray{T,N},i::Vararg{Integer,N}) where {T,N}
  @boundscheck BlockArrays.blockcheckbounds(a, i...)
  false
end

function eachblockid(a)
  cis = CartesianIndices(blocksize(a))
  map(ci->Block(Tuple(ci)),cis)
end

BlockArrays.eachblock(a::BlockArrayCoo) = ( a[I]  for I in eachblockid(a) )
enumerateblocks(a) = zip(eachblockid(a),eachblock(a))

# AbstractBlockArray

# Avoid the internal copy in BlockArrayCoo 0.13
@propagate_inbounds Base.getindex( block_arr::BlockArrayCoo{T,N}, block::Block{N}) where {T,N} =  getblock(block_arr, block.n...)
@propagate_inbounds Base.getindex( block_arr::BlockVectorCoo, block::Block{1}) =  getblock(block_arr, block.n[1])
@inline Base.getindex(block_arr::BlockArrayCoo{T,N}, block::Vararg{Block{1}, N}) where {T,N} =  getblock(block_arr, (Block(block).n)...)

@inline function BlockArrays.getblock(a::BlockArrayCoo{T,N}, block::Vararg{Integer, N}) where {T,N}
  #@boundscheck BlockArrays.blockcheckbounds(a, block...)
  p = a.ptrs[block...]
  if p>0
    a.blocks[p]
  else
    a.zero_blocks[-p]
  end
end

function BlockArrays.getblock!(c,a::BlockArrayCoo{T,N}, block::Vararg{Integer, N}) where {T,N}
  b = a[Block(block...)]
  copy!(c,b)
  c
end

@inline function BlockArrays.setblock!(a::BlockArrayCoo{T, N}, v, block::Vararg{Integer, N}) where {T,N}
  #@boundscheck BlockArrays.blockcheckbounds(a, block...)
  p = a.ptrs[block...]
  if p>0
    a.blocks[p] = v
  else
    @unreachable "It is not possible to set a zero block in BlockArrayCoo"
  end
  a
end

# AbstractArray

Base.size(a::BlockArrayCoo) = map(length,Base.axes(a))
Base.axes(a::BlockArrayCoo) = a.axes
Base.IndexStyle(::Type{<:BlockArrayCoo}) = IndexCartesian()

function Base.getindex(a::BlockArrayCoo{T,N},i::Vararg{Integer,N}) where {T,N}
  s = map(findblockindex,a.axes,i)
  ai = a[s...]
  ai
end

function Base.setindex!(a::BlockArrayCoo{T,N},v,i::Vararg{Integer,N}) where {T,N}
  s = map(findblockindex,a.axes,i)
  I = Block(map(i->i.I[1],s)...)
  α = CartesianIndex(map(BlockArrays.blockindex,s))
  a[I][α] = v
end

# Fast copy

function Base.copy(a::BlockArrayCoo)
  blocks = copy.(a.blocks)
  blockids = copy(a.blockids)
  axes = copy.(a.axes)
  ptrs = copy(a.ptrs)
  zero_blocks = copy.(a.zero_blocks)
  BlockArrayCoo(axes,blockids,blocks,ptrs,zero_blocks)
end

function Base.copy!(a::BlockArrayCoo,b::BlockArrayCoo)
  copyto!(a,b)
end

function Base.copyto!(a::BlockArrayCoo,b::BlockArrayCoo)
  if a.ptrs === b.ptrs || a.ptrs == b.ptrs
    for p in 1:length(a.blocks)
      copyto!(a.blocks[p],b.blocks[p])
    end
  else
    for I in eachblockid(a)
      if is_nonzero_block(a,I)
        copyto!(a[I],b[I])
      end
    end
  end
 a
end

# Similar
# We customize similar only for matching axis type
# Otherwise we do not have enough info to create a BlockArrayCoo

function Base.similar(a::BlockArrayCoo)
  similar(a,eltype(a),axes(a))
end

function Base.similar(a::BlockArrayCoo,::Type{T}) where T
  similar(a,T,axes(a))
end

function Base.similar(a::BlockArrayCoo{S,N,A,X},::Type{T}, _axes::X) where {S,A,T,N,X}
  _similar_block_array_coo(a,T,_axes)
end

function Base.similar(
  a::BlockArrayCoo{S,N,A,X},::Type{T}, _axes::X) where {S,A,T,N,X<:Tuple{Vararg{BlockedUnitRange}}}
  _similar_block_array_coo(a,T,_axes)
end

function _similar_block_array_coo(a,::Type{T},_axes) where T
  if num_blocks_equal(axes(a),_axes)
    _similar_block_array_coo_preserving(a,T,_axes)
  else
    _similar_block_array_coo_non_preserving(a,T,_axes)
  end
end

# similar preserving zero block structure
function _similar_block_array_coo_preserving(a,::Type{T},_axes) where T
  ai = first(eachblock(a))
  A = typeof(similar(ai,T))
  blocks = A[]
  for i in 1:length(a.blocks)
    ai = a.blocks[i]
    I = a.blockids[i]
    laxs = map( local_range, _axes, I)
    block = similar(ai,T,laxs)
    push!(blocks,block)
  end
  blockids = copy(a.blockids)
  BlockArrayCoo(_axes,blockids,blocks)
end

# similar without zero block structure (all blocks are assumed to be non-zero)
function _similar_block_array_coo_non_preserving(a::BlockArrayCoo{S,N,A},::Type{T},_axes) where {S,N,A,T}
  s = map(i->first(blocksize(i)),_axes)
  cis = CartesianIndices(s)
  blocks = A[]
  blockids = NTuple{N,Int}[]
  for ci in cis
    I = Tuple(ci)
    laxs = map( local_range, _axes, I)
    block = similar(similar(A,laxs),T,laxs)
    push!(blocks,block)
    push!(blockids,I)
  end
  BlockArrayCoo(_axes,blockids,blocks)
end

# In this case we can not preserve the zero block structure neither (all blocks are assumed to be non-zero).
function Base.similar(::Type{<:BlockArrayCoo{S,N,A,X}},_axes::X) where {S,N,A,X}
  s = map(i->first(blocksize(i)),_axes)
  cis = CartesianIndices(s)
  blocks = A[]
  blockids = NTuple{N,Int}[]
  for ci in cis
    I = Tuple(ci)
    laxs = map( local_range, _axes, I)
    block = similar(A,laxs)
    push!(blocks,block)
    push!(blockids,I)
  end
  BlockArrayCoo(_axes,blockids,blocks)
end

# More API

function Base.fill!(a::BlockArrayCoo,v)
  @notimplementedif v != zero(v)
  for b in a.blocks
    fill!(b,v)
  end
  a
end

function fill_entries!(a::BlockArrayCoo,v)
  for b in a.blocks
    fill_entries!(b,v)
  end
  a
end

function scale_entries!(c::BlockArrayCoo,β)
  for block in c.blocks
    scale_entries!(block,β)
  end
  c
end

Base.zero(a::BlockArrayCoo) = _zero_block(typeof(a),axes(a))

# + and -

for op in (:+,:-)
  @eval begin
    function Base.$op(a::BlockArrayCoo{Ta,N},b::BlockArrayCoo{Tb,N}) where {Ta,Tb,N}
      @check blocks_equal(axes(a),axes(b))
      I1 = first(eachblockid(a))
      A = typeof($op(a[I1],b[I1]))
      blocks = A[]
      blockids = NTuple{N,Int}[]
      for (I,aI) in enumerateblocks(a)
        bI = b[I]
        if is_nonzero_block(a,I) || is_nonzero_block(b,I)
          block = $op(aI,bI)
          push!(blocks,block)
          push!(blockids,I.n)
        end
      end
      BlockArrayCoo(a.axes,blockids,blocks)
    end
  end
end

Base.:+(a::BlockArrayCoo) = a

function Base.:-(a::BlockArrayCoo)
  b = copy(a)
  scale_entries!(b,-1)
  b
end

# Products

function Base.:*(a::BlockArrayCoo,b::Number)
  *(b,a)
end

function Base.:*(a::Number,b::BlockArrayCoo)
  f(block)= a*block
  blocks = f.(b.blocks)
  BlockArrayCoo(b.axes,b.blockids,blocks,b.ptrs,f.(b.zero_blocks))
end

function Base.:*(a::BlockMatrixCoo,b::BlockVectorCoo)
  @assert blocksize(a,2) == blocksize(b,1)
  A = typeof(a[Block(1,1)]*b[Block(1)])
  blocks = A[]
  blockids = Tuple{Int}[]
  for i in 1:blocksize(a,1)
    block = zero(a[Block(i,1)]*b[Block(1)])
    for j in 1:blocksize(a,2)
      if is_nonzero_block(a,Block(i,j)) && is_nonzero_block(b,Block(j))
        block += a[Block(i,j)]*b[Block(j)]
      end
    end
    push!(blocks,block)
    push!(blockids,(i,))
  end
  axs = (axes(a)[1],)
  BlockArrayCoo(axs,blockids,blocks)
end

function Base.:*(a::BlockMatrixCoo,b::BlockMatrixCoo)
  @assert blocksize(a,2) == blocksize(b,1)
  A = typeof(a[Block(1,1)]*b[Block(1,1)])
  blocks = A[]
  blockids = Tuple{Int,Int}[]
  for i in 1:blocksize(a,1)
    for j in 1:blocksize(b,2)
      block = zero(a[Block(i,1)]*b[Block(1,j)])
      for k in 1:blocksize(a,2)
        if is_nonzero_block(a,Block(i,k)) && is_nonzero_block(b,Block(k,j))
          block += a[Block(i,k)]*b[Block(k,j)]
        end
      end
      push!(blocks,block)
      push!(blockids,(i,j))
    end
  end
  axs = (axes(a)[1],axes(b)[2])
  BlockArrayCoo(axs,blockids,blocks)
end

function LinearAlgebra.mul!(c::BlockVectorCoo,a::BlockMatrixCoo,b::BlockVectorCoo)
  fill_entries!(c,zero(eltype(c)))
  mul!(c,a,b,1,0)
end

function LinearAlgebra.mul!(c::BlockMatrixCoo,a::BlockMatrixCoo,b::BlockMatrixCoo)
  fill_entries!(c,zero(eltype(c)))
  mul!(c,a,b,1,0)
end

function LinearAlgebra.mul!(c::BlockVectorCoo,a::BlockMatrixCoo,b::BlockVectorCoo,α::Number,β::Number)
  for I in 1:blocksize(a,1)
    cI = c[Block(I)]
    if is_nonzero_block(c,Block(I))
      scale_entries!(cI,β)
    end
    for J in 1:blocksize(a,2)
      if is_nonzero_block(a,Block(I,J)) && is_nonzero_block(b,Block(J))
        @check is_nonzero_block(c,Block(I))
        aIJ = a[Block(I,J)]
        bJ = b[Block(J)]
        mul!(cI,aIJ,bJ,α,1)
      end
    end
  end
  c
end

function LinearAlgebra.mul!(c::BlockMatrixCoo,a::BlockMatrixCoo,b::BlockMatrixCoo,α::Number,β::Number)
  for I in 1:blocksize(a,1)
    for J in 1:blocksize(b,2)
      cIJ = c[Block(I,J)]
      if is_nonzero_block(c,Block(I,J))
        scale_entries!(cIJ,β)
      end
      for K in 1:blocksize(a,2)
        if is_nonzero_block(a,Block(I,K)) && is_nonzero_block(b,Block(K,J))
          @check is_nonzero_block(c,Block(I,J))
          cIJ = c[Block(I,J)]
          aIK = a[Block(I,K)]
          bKJ = b[Block(K,J)]
          mymul!(cIJ,aIK,bKJ,α) # Hack to avoid allocations
          #mul!(cIJ,aIK,bKJ,α,1) # Why this leads to memory allocations?? Specially with Ints
        end
      end
    end
  end
  c
end

# Hack to avoid allocations
@inline function mymul!(cIJ,aIK,bKJ,α)
  @boundscheck begin
    @assert size(cIJ,1) == size(aIK,1)
    @assert size(cIJ,2) == size(bKJ,2)
    @assert size(aIK,2) == size(bKJ,1)
  end
  for i in 1:size(cIJ,1)
    for j in 1:size(bKJ,2)
      for k in 1:size(bKJ,1)
        @inbounds cIJ[i,j] += α*aIK[i,k]*bKJ[k,j]
      end
    end
  end
end

# Hack to avoid allocations
@inline function mymul!(cIJ::BlockMatrixCoo,aIK::BlockMatrixCoo,bKJ::BlockMatrixCoo,α)
  mul!(cIJ,aIK,bKJ,α,1)
end

function LinearAlgebra.transpose(a::BlockMatrixCoo)
  blocks = transpose.(a.blocks)
  zero_blocks = transpose.(a.zero_blocks)
  blockids = [ (j,i) for (i,j) in a.blockids ]
  ax,ay = axes(a)
  axs = (ay,ax)
  ptrs = collect(transpose(a.ptrs))
  BlockArrayCoo(axs,blockids,blocks,ptrs,zero_blocks)
end