text:
export iauAtoc13
"""
Observed place at a groundbased site to to ICRS astrometric RA,Dec.
The caller supplies UTC, site coordinates, ambient air conditions
and observing wavelength.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   type   char[]   type of coordinates - "R", "H" or "A" (Notes 1,2)
   ob1    double   observed Az, HA or RA (radians; Az is N=0,E=90)
   ob2    double   observed ZD or Dec (radians)
   utc1   double   UTC as a 2-part...
   utc2   double   ...quasi Julian Date (Notes 3,4)
   dut1   double   UT1-UTC (seconds, Note 5)
   elong  double   longitude (radians, east +ve, Note 6)
   phi    double   geodetic latitude (radians, Note 6)
   hm     double   height above ellipsoid (m, geodetic Notes 6,8)
   xp,yp  double   polar motion coordinates (radians, Note 7)
   phpa   double   pressure at the observer (hPa = mB, Note 8)
   tc     double   ambient temperature at the observer (deg C)
   rh     double   relative humidity at the observer (range 0-1)
   wl     double   wavelength (micrometers, Note 9)

Returned:
   rc,dc  double   ICRS astrometric RA,Dec (radians)

Returned (function value):
         int      status: +1 = dubious year (Note 4)
                           0 = OK
                           -1 = unacceptable date

Notes:

   1.  "Observed" Az,ZD means the position that would be seen by a
   perfect geodetically aligned theodolite.  (Zenith distance is
   used rather than altitude in order to reflect the fact that no
   allowance is made for depression of the horizon.)  This is
   related to the observed HA,Dec via the standard rotation, using
   the geodetic latitude (corrected for polar motion), while the
   observed HA and RA are related simply through the Earth rotation
   angle and the site longitude.  "Observed" RA,Dec or HA,Dec thus
   means the position that would be seen by a perfect equatorial
   with its polar axis aligned to the Earth's axis of rotation.

   2.  Only the first character of the type argument is significant.
   "R" or "r" indicates that ob1 and ob2 are the observed right
   ascension and declination;  "H" or "h" indicates that they are
   hour angle (west +ve) and declination;  anything else ("A" or
   "a" is recommended) indicates that ob1 and ob2 are azimuth
   (north zero, east 90 deg) and zenith distance.

   3.  utc1+utc2 is quasi Julian Date (see Note 2), apportioned in any
   convenient way between the two arguments, for example where utc1
   is the Julian Day Number and utc2 is the fraction of a day.

   However, JD cannot unambiguously represent UTC during a leap
   second unless special measures are taken.  The convention in the
   present function is that the JD day represents UTC days whether
   the length is 86399, 86400 or 86401 SI seconds.

   Applications should use the function iauDtf2d to convert from
   calendar date and time of day into 2-part quasi Julian Date, as
   it implements the leap-second-ambiguity convention just
   described.

   4.  The warning status "dubious year" flags UTCs that predate the
   introduction of the time scale or that are too far in the
   future to be trusted.  See iauDat for further details.

   5.  UT1-UTC is tabulated in IERS bulletins.  It increases by exactly
   one second at the end of each positive UTC leap second,
   introduced in order to keep UT1-UTC within +/- 0.9s.  n.b. This
   practice is under review, and in the future UT1-UTC may grow
   essentially without limit.

   6.  The geographical coordinates are with respect to the WGS84
   reference ellipsoid.  TAKE CARE WITH THE LONGITUDE SIGN:  the
   longitude required by the present function is east-positive
   (i.e. right-handed), in accordance with geographical convention.

   7.  The polar motion xp,yp can be obtained from IERS bulletins.  The
   values are the coordinates (in radians) of the Celestial
   Intermediate Pole with respect to the International Terrestrial
   Reference System (see IERS Conventions 2003), measured along the
   meridians 0 and 90 deg west respectively.  For many
   applications, xp and yp can be set to zero.

   8.  If hm, the height above the ellipsoid of the observing station
   in meters, is not known but phpa, the pressure in hPa (=mB), is
   available, an adequate estimate of hm can be obtained from the
   expression

         hm = -29.3 * tsl * log ( phpa / 1013.25 );

   where tsl is the approximate sea-level air temperature in K
   (See Astrophysical Quantities, C.W.Allen, 3rd edition, section
   52).  Similarly, if the pressure phpa is not known, it can be
   estimated from the height of the observing station, hm, as
   follows:

         phpa = 1013.25 * exp ( -hm / ( 29.3 * tsl ) );

   Note, however, that the refraction is nearly proportional to
   the pressure and that an accurate phpa value is important for
   precise work.

   9.  The argument wl specifies the observing wavelength in
   micrometers.  The transition from optical to radio is assumed to
   occur at 100 micrometers (about 3000 GHz).

   10. The accuracy of the result is limited by the corrections for
   refraction, which use a simple A*tan(z) + B*tan^3(z) model.
   Providing the meteorological parameters are known accurately and
   there are no gross local effects, the predicted astrometric
   coordinates should be within 0.05 arcsec (optical) or 1 arcsec
   (radio) for a zenith distance of less than 70 degrees, better
   than 30 arcsec (optical or radio) at 85 degrees and better
   than 20 arcmin (optical) or 30 arcmin (radio) at the horizon.

   Without refraction, the complementary functions iauAtco13 and
   iauAtoc13 are self-consistent to better than 1 microarcsecond
   all over the celestial sphere.  With refraction included,
   consistency falls off at high zenith distances, but is still
   better than 0.05 arcsec at 85 degrees.

   11. It is advisable to take great care with units, as even unlikely
   values of the input parameters are accepted and processed in
   accordance with the models used.

Called:
   iauApco13    astrometry parameters, ICRS-observed
   iauAtoiq     quick observed to CIRS
   iauAticq     quick CIRS to ICRS

This revision:   2013 October 9

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""
# int iauAtoc13(const char *type, double ob1, double ob2,
#               double utc1, double utc2, double dut1,
#               double elong, double phi, double hm, double xp, double yp,
#               double phpa, double tc, double rh, double wl,
#               double *rc, double *dc)

function iauAtoc13(type_::Char, ob1::Real, ob2::Real,
               utc1::Real, utc2::Real, dut1::Real,
               elong::Real, phi::Real, hm::Real, xp::Real, yp::Real,
               phpa::Real, tc::Real, rh::Real, wl::Real)

   # Allocate return values
   ref_type_  = Ref{UInt8}(convert(UInt8, type_))
   ref_rc     = Ref{Float64}(0.0)
   ref_dc     = Ref{Float64}(0.0)

   status = ccall((:iauAtoc13, libsofa_c), Cint,
       (Ref{UInt8}, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble,
       Cdouble, Cdouble, Cdouble, Cdouble, Cdouble,
       Cdouble, Cdouble, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
       ref_type_, convert(Float64, ob1), convert(Float64, ob2),
       convert(Float64, utc1), convert(Float64, utc2),
       convert(Float64, dut1), convert(Float64, elong),
       convert(Float64, phi), convert(Float64, hm),
       convert(Float64, xp), convert(Float64, yp),
       convert(Float64, phpa), convert(Float64, tc),
       convert(Float64, rh), convert(Float64, wl),
       ref_rc, ref_dc)

   return status, ref_rc[], ref_dc[]
end


text:
function Base.show(io::IO, ms::ModelingStruct{T,V,K}) where {T,V,K}
    key = K[1]
    print(io, "($key = $(ms[key])")
    for idx in 2:length(K)
        key = K[idx]
        print(io, ", $key = $(ms[key])")
    end
    print(io, ")")
    return nothing
end
function Base.show(io::IO, ::MIME"text/plain", ms::ModelingStruct{T}) where T
    print(io, "ModelingStruct{" , T, "}")
    show(io, ms)
    return nothing
end
function Base.show(io::IO, a::AbstractVector{<:T}) where T<:ModelingStruct
    elem = a[1]
    print(io, "[$elem")
    for idx in 2:length(a)
        print(io, ", $(a[idx])")
    end
    print(io, "]")
    return nothing
end

ext: jl

max_stars_repo_name: jonniedie/NamedViewVectors


text:
#=
Popis
  Funkcia vypocita pre kazdy vrchol hodnotu jeho degree centrality.

Parametre
  g : datova struktura grafu :: GenericGraph

Navratova hodnota
  degree : slovnik obsahujuci mapovanie vrcholov a k nim prisluchajucich hodnot
           degree centrality (index_vrchola => hodnota_degree_centrality) :: Dict{Int64, Float64}()

Poznamky
  Hodnota centrality je normalizovana (norm_coefficient).
=#

function centrality_degree(g::GenericGraph)

  degree = Dict{Int64, Float64}()

  norm_coefficient = 1.0 / (num_vertices(g) - 1.0)

  for vertex in vertices(g)

    degree[vertex.index] = out_degree(vertex, g) * norm_coefficient

  end

  return degree

end

max_stars_repo_name: mpuk/GraphAnalysis

max_issues_repo_name: mpuk/GraphAnalysis

max_forks_repo_name: mpuk/GraphAnalysis


text:
"""
    is_foggy()

Return `true` if it's foggy.
"""
is_foggy() = rand(Bool)

"""
    is_dog_distracted()

Return `true` if Annalyn's dog has been distracted.
"""
is_dog_distracted() = rand() < 0.25

"""
    loot()

Return the number of coins Annalyn finds in the camp.
"""
loot() = rand(3:13)

"""
    loot(crate)

Return which out of the given items Annalyn finds in the chest.
"""
loot(crate) = rand(crate)


text:
"""
PoolQueues facilitate sharing pools of items between producer Tasks and consumer
Tasks.  See `PoolQueue` for more information.
"""
module PoolQueues

export PoolQueue
export acquire!
export produce!
export consume!
export recycle!
export produce_on_command

"""
A PoolQueue facilitates sharing a pool of items between a producer Task and a
consumer Task.  The four main operations on a PoolQueue are: `acquire!`,
`produce!`, `consume!`, and `recycle!`.  The first two are called by producer
Tasks; the latter two by consumer Tasks.  Typically the PoolQueue's pool is
prepopluated with preallocated items (such as Arrays or user defined structs).
PoolQueue constructor methods exist to facilitate this.  By recycling these
items in the PoolQueue, memory allocations (and garbage collection) can be
minimized.

The general flow for the producer Task is:

```julia
while true
    # Acquire an available item from the PoolQueue's pool
    item = acquire!(poolqueue)

    # Preprare item for consumer task (application specific)

    # Produce the item to the PoolQueue's queue
    produce!(poolqueue, item)
end
```

The general flow for the consumer Task is:

```julia
while true
    # Consume an item from the PoolQueue's queue
    item = consume!(poolqueue)

    # Process the item (application specific)

    # Recycle the item back to the PoolQueue's pool
    produce!(poolqueue, item)
end
```

# Example use case

One usage scenario is for the producer task to read a portion of a data file,
send that to the consumer task for processing, and then read the next portion of
the data file.  The producer Task's reading of the next portion of data happens
in parallel with the consumer Task's processing the previous data.  Instead of:

    main task: read0 process0 read1 process1 read2 process2 ... [time -->]

using a PoolQueue with two (or more) items allows:

    producer task: read0 read1    read2    ... [time -->]
    consumer task:       process0 process1 ... [time -->]
"""
struct PoolQueue{C}
    pool::C
    queue::C
    PoolQueue{C}(p::C, q::C) where {T, C<:AbstractChannel{T}} = new(p, q)
end

"""
    PoolQueue(p::C, q::C) where {T, C<:AbstractChannel{T}}

Construct a PoolQueue from two `AbstractChannel{T}` instances.
"""
function PoolQueue(p::C, q::C) where {T, C<:AbstractChannel{T}}
    PoolQueue{C}(p, q)
end

"""
    PoolQueue{T}(np::Integer, nq::Integer=np) where {T}

Construct a PoolQueue using `Channel{T}` channels.  The pool channel will hold
up to `np` items of type `T` and the queue channel will hold up to `np` items of
type `T`.
"""
function PoolQueue{T}(np::Integer, nq::Integer=np) where {T}
    np > 0 || throw(ArgumentError("pool size must be positive"))
    nq > 0 || throw(ArgumentError("queue size must be positive"))
    PoolQueue{Channel{T}}(Channel{T}(np), Channel{T}(nq))
end

"""
    PoolQueue{T}(f::Function, np::Integer, nq::Integer=np, fargs...; fkwargs...) where {T}

Construct a PoolQueue using `Channel{T}` channels.  The pool channel will hold
up to `np` items of type `T` and the queue channel will hold up to `np` items of
type `T`.  The function `f`, which should return a single item of type `T`, will
be called `np` times as `f(fargs...; fkwargs...)` to prepopulate the PoolQueue's
pool.  If `fargs` is used, `nq` must be passed explicitly.
"""
function PoolQueue{T}(f::Function, np::Integer, nq::Integer=np, fargs...; fkwargs...) where {T}
    pq = PoolQueue{T}(np, nq)
    for _ in 1:np
        recycle!(pq, f(fargs...; fkwargs...))
    end
    pq
end

"""
    PoolQueue(::Type{T}, np::Integer, nq::Integer=np, fargs...; fkwargs...) where {T}

Construct a PoolQueue using `Channel{T}` channels.  The pool channel will hold
up to `np` items of type `T` and the queue channel will hold up to `np` items of
type `T`.  The constructor of `T` will be called `np` times as `T(fargs...;
fkwargs...)` to prepopulate the PoolQueue's pool.  If `fargs` is used, `nq` must
be passed explicitly.
"""
function PoolQueue(::Type{T}, np::Integer, nq::Integer=np, fargs...; fkwargs...) where {T}
    PoolQueue{T}((a...; k...)->T(a...; k...), np, nq, fargs...; fkwargs...)
end

"""
Close the `pool` and `queue` channels associated with `pq`.
"""
function Base.close(pq::PoolQueue)
    close(pq.queue)
    close(pq.pool)
end

"""
    acquire!(pq::PoolQueue{C})::T where {T, C<:AbstractChannel{T}}

Acquire an available item from `pq.pool`.
"""
function acquire!(pq::PoolQueue{C})::T where {T, C<:AbstractChannel{T}}
    take!(pq.pool)
end

"""
    produce!(pq::PoolQueue{C}, item::T)::T where {T, C<:AbstractChannel{T}}

Produce `item` to `pq.queue`.
"""
function produce!(pq::PoolQueue{C}, item::T)::T where {T, C<:AbstractChannel{T}}
    put!(pq.queue, item)
end

"""
    produce!(f::Function, pq::PoolQueue{C}, fargs...)::T where {T, C<:AbstractChannel{T}}

Produce an item by acquiring an available item from `pq.pool`, call `f(item,
fargs...)`, and `produce!` the value returned by `f` unless it is `nothing`.  If
`f` returns `nothing`, the item is recycled without being produced.  The value
returned by `f`, which is of type `T` or `nothing`, is returned from this
function.  returned by this call.
"""
function produce!(f::Function, pq::PoolQueue{C}, fargs...)::Union{T, Nothing} where {T, C<:AbstractChannel{T}}
    poolitem = take!(pq.pool)
    produceitem = f(poolitem, fargs...)
    produceitem === nothing ? put!(pq.pool, poolitem) : put!(pq.queue, produceitem)
    produceitem
end

"""
    consume!(pq::PoolQueue{C})::T where {T, C<:AbstractChannel{T}}

Consume an item from `pq.queue`.
"""
function consume!(pq::PoolQueue{C})::T where {T, C<:AbstractChannel{T}}
    take!(pq.queue)
end

"""
    consume!(f::Function, pq::PoolQueue{C}, fargs....)::T where {T, C<:AbstractChannel{T}}

Consume an item from `pq.queue`, call `f(item, fargs...)`, and `recycle!` the
value returned by `f`.
"""
function consume!(f::Function, pq::PoolQueue{C}, fargs...)::T where {T, C<:AbstractChannel{T}}
    take!(pq.queue) |> item->f(item, fargs...) |> item->put!(pq.pool, item)
end

"""
    recycle!(pq::PoolQueue{C}, item::T)::T where {T, C<:AbstractChannel{T}}

Recycle `item` back to `pq.pool`.
"""
function recycle!(pq::PoolQueue{C}, item::T)::T where {T, C<:AbstractChannel{T}}
    put!(pq.pool, item)
end

"""
    produce_on_command(produce, cmd::AbstractChannel{T}, pq::PoolQueue{C};
                       autoclose=true)

For each `command` taken from `cmd`, call `produce(command, pq)`.  This function
will typically be called within a task:

    producer_task = @task produce_on_command(myproducerfunc, cmd, pq)

If an exception is thrown while taking from `cmd`, an `@info` message is logged
and the funtion will return.  If an exception is thrown within `produce`, a
`@warn` message is logged at the function will return.  If `autoclose` is
`true` (the default), then `cmd` and `pq` will be closed before returning.
"""
function produce_on_command(produce, cmd::AbstractChannel{T}, pq::PoolQueue{C};
                            autoclose=true
                           ) where {T<:AbstractString, C<:AbstractChannel}
    # Command loop
    while true
        command = try
            # Take rawstem from cmd channel
            take!(cmd)
        catch
            @info "got exception from command channel [done]"
            break
        end

        @debug "calling $produce with command $command"
        try
            produce(command, pq)
        catch
            @warn "got exception from produce function [done]"
            break
        end
    end

    if autoclose
        close(pq.queue)
        close(cmd)
    end
    nothing
end

end # module PoolQueues


text:
using .Tapir: ConcreteMaybe

"""
    WorkStealingTaskGroup()

A custom taskgroup implementation for `Base.Experimental.Tapir` using
work-stealing scheduler.

NOTE: This scheduler does not support any possibly blocking synchronization
between Tapir tasks between *any* syncregions (including the syncregions placed
in different concurrent `Task`s).

# Examples
```julia
Tapir.@sync WorkStealingTaskGroup() begin
    ...
end
```
"""
WorkStealingTaskGroup

function _WorkStealingTaskGroup end

refresh() = @eval invalidator() = nothing
invalidator() = nothing

struct WorkStealingTaskGroup
    global _WorkStealingTaskGroup() = new()
end

function WorkStealingTaskGroup()
    invalidator()
    return _WorkStealingTaskGroup()
end

# TODO: implement
function Tapir.spawn!(tg::WorkStealingTaskGroup, @nospecialize(f))
    error("Tapir.spawn!(::WorkStealingTaskGroup, _) not implemented yet")
    push!(tg, spawn(f))
end

#=
function Tapir.sync!(tg::WorkStealingTaskGroup)
    foreach(wait_nothrow, tg)
    ref = Ref{Union{Nothing,CompositeException}}(nothing)
    foreach(tg) do work
        if work.iserror
            local ex = ref[]
            if ex === nothing
                ex = ref[] = CompositeException()
            end
            push!(ex, work.result)
        end
    end
    ex = ref[]
    if ex !== nothing
        throw(ex)
    end
end
=#

Tapir.spawn(::Type{WorkStealingTaskGroup}, @nospecialize(f)) =
    spawn!(f, WORK_STEALING_SCHEDULER[])

function Tapir.synctasks(args::ConcreteMaybe{Work}...)
    ex = nothing
    for t in args
        if t !== nothing
            work::Work = something(t)
            wait_nothrow(work)
            if work.iserror
                if ex === nothing
                    ex = CompositeException()
                end
                push!(ex::CompositeException, work.result)
            end
        end
    end
    if ex !== nothing
        throw(ex)
    end
end

macro sync_ws(block)
    esc(:($Tapir.@sync $WorkStealingTaskGroup() $block))
end

"""
    DepthFirstTaskGroup()

A custom taskgroup implementation for `Base.Experimental.Tapir` using a relaxed
depth-first scheduler.

**NOTE**: Currently, it must used via `@sync_df`.
"""
DepthFirstTaskGroup

function _DepthFirstTaskGroup end

struct DepthFirstTaskGroup
    global _DepthFirstTaskGroup() = new()
end

function DepthFirstTaskGroup()
    invalidator()
    return _DepthFirstTaskGroup()
end

# TODO: implement
function Tapir.spawn!(tg::DepthFirstTaskGroup, @nospecialize(f))
    error("Tapir.spawn!(::DepthFirstTaskGroup, _) not implemented yet")
    push!(tg, spawn!(f, DEPTH_FIRST_SCHEDULER[]))
end

Tapir.spawn(::Type{DepthFirstTaskGroup}, @nospecialize(f)) =
    spawn!(f, DEPTH_FIRST_SCHEDULER[])

"""
    @sync_df block

Run Tapir tasks inside a `DepthFirstTaskGroup`.
"""
macro sync_df(block)
    @gensym pr
    if Meta.isexpr(block, :block)
        block = Expr(:block, __source__, block)
    end
    quote
        # TODO: make it possible to define this via Tapir entry points
        $pr = $get_current_priority_range()
        try
            $Tapir.@sync $DepthFirstTaskGroup() $block
        finally
            $set_current_priority_range($pr)
        end
    end |> esc
end

"""
    ConstantPriorityTaskGroup()

Almost equivalent to `DepthFirstTaskGroup` but the priority is constant.
"""
ConstantPriorityTaskGroup

function _ConstantPriorityTaskGroup end

struct ConstantPriorityTaskGroup
    global _ConstantPriorityTaskGroup() = new()
end

function ConstantPriorityTaskGroup()
    invalidator()
    return _ConstantPriorityTaskGroup()
end

# TODO: implement
function Tapir.spawn!(tg::ConstantPriorityTaskGroup, @nospecialize(f))
    error("Tapir.spawn!(::ConstantPriorityTaskGroup, _) not implemented yet")
    push!(tg, spawn!(f, CONSTANT_PRIORITY_SCHEDULER[]))
end

Tapir.spawn(::Type{ConstantPriorityTaskGroup}, @nospecialize(f)) =
    spawn!(f, CONSTANT_PRIORITY_SCHEDULER[])

"""
    @sync_cp block

Run Tapir tasks inside a `ConstantPriorityTaskGroup`.
"""
macro sync_cp(block)
    @gensym pr
    if Meta.isexpr(block, :block)
        block = Expr(:block, __source__, block)
    end
    quote
        # TODO: make it possible to define this via Tapir entry points
        $pr = $get_current_priority_range()
        try
            $Tapir.@sync $ConstantPriorityTaskGroup() $block
        finally
            $set_current_priority_range($pr)
        end
    end |> esc
end

"""
    RandomPriorityTaskGroup()

Almost equivalent to `DepthFirstTaskGroup` but the priority is randomized.
"""
RandomPriorityTaskGroup

function _RandomPriorityTaskGroup end

struct RandomPriorityTaskGroup
    global _RandomPriorityTaskGroup() = new()
end

function RandomPriorityTaskGroup()
    invalidator()
    return _RandomPriorityTaskGroup()
end

# TODO: implement
function Tapir.spawn!(tg::RandomPriorityTaskGroup, @nospecialize(f))
    error("Tapir.spawn!(::RandomPriorityTaskGroup, _) not implemented yet")
    push!(tg, spawn!(f, RANDOM_PRIORITY_SCHEDULER[]))
end

Tapir.spawn(::Type{RandomPriorityTaskGroup}, @nospecialize(f)) =
    spawn!(f, RANDOM_PRIORITY_SCHEDULER[])

"""
    @sync_rp block

Run Tapir tasks inside a `RandomPriorityTaskGroup`.
"""
macro sync_rp(block)
    @gensym pr
    if Meta.isexpr(block, :block)
        block = Expr(:block, __source__, block)
    end
    quote
        # TODO: make it possible to define this via Tapir entry points
        $pr = $get_current_priority_range()
        try
            $Tapir.@sync $RandomPriorityTaskGroup() $block
        finally
            $set_current_priority_range($pr)
        end
    end |> esc
end


text:
"""
    GuidedProposal{T,Ttarget,Taux,TL,Txobs0,TxobsT,Tnshapes,TmT,F} <: ContinuousTimeProcess{T}

struct that contains target, auxiliary process for each shape, time grid, observation at time 0, observations
at time T, number of shapes, and momenta in final state used for constructing the auxiliary processes
`guidrec` is a vector of GuidRecursions, which contains the results from the backward recursions and gpupdate step at time zero

## Arguments
```julia
target::Ttarget                 # target diffusion P
aux::Vector{Taux}               # auxiliary diffusion for each shape (Ptilde for each shape)
tt::Vector{Float64}             # grid of time points on single segment (S,T]
xobs0::Txobs0                   # observation at time 0
xobsT::Vector{TelxobsT}         # vector of observations at time T
guidrec::Vector{TL}             # guided recursions on grid tt
nshapes::Int64                  # number of shapes
mT::Vector{TmT}                 # vector of artificial momenta used for constructing auxiliary process (one for each shape)
endpoint::F
```

## Example (for landmarksmatching)
```julia
    n = 5
    P = MarslandShardlow(1.0, 2.0, 0.5, 0.0, n)
    xT = State(rand(PointF,n),rand(PointF,n))
    Paux = [BridgeLandmarks.auxiliary(P,xT)]
    tt = collect(0:0.05:1.0)
    xobs0 = rand(PointF,n)
    xobsT = [rand(PointF,n)]
    σobs = 0.01
    Σobs = fill([σobs^2 * one(UncF) for i in 1:n],2)
    obsinfo = BridgeLandmarks.set_obsinfo(xobs0,xobsT,Σobs, true,false)
    gr = [BridgeLandmarks.init_guidrec(tt,obsinfo)]
    nshapes = 1
    mT = [rand(PointF,n)]
    Q = BridgeLandmarks.GuidedProposal(P,Paux,tt,xobs0,xobsT,gr,nshapes,mT)
```
"""
struct GuidedProposal{T,Ttarget,Taux,Txobs0,TelxobsT,TL,TmT,F} <: ContinuousTimeProcess{T}
    target::Ttarget                 # target diffusion P
    aux::Vector{Taux}               # auxiliary diffusion for each shape (Ptilde for each shape)
    tt::Vector{Float64}             # grid of time points on single segment (S,T]
    xobs0::Txobs0                   # observation at time 0
    xobsT::Vector{TelxobsT}         # vector of observations at time T
    guidrec::Vector{TL}             # guided recursions on grid tt
    nshapes::Int64                  # number of shapes
    mT::Vector{TmT}                 # vector of artificial momenta used for constructing auxiliary process (one for each shape)
    endpoint::F

    function GuidedProposal(target, aux, tt_,xobs0,xobsT, guidrec,nshapes,mT, endpoint=Bridge.endpoint)
        tt = collect(tt_)
        new{Bridge.valtype(target),typeof(target),eltype(aux),typeof(xobs0), eltype(xobsT),eltype(guidrec),eltype(mT),typeof(endpoint)}(
                    target, aux, tt,xobs0,xobsT, guidrec, nshapes,mT, endpoint)
    end
end

"""
    set_guidrec(Q::GuidedProposal, gr) = GuidedProposal(Q.target, Q.aux, X.tt, Q.xobs, Q.xobsT, gr,Q.nshapes, Q.endpoint)

Update `guidrec` field of Q to `gr`
"""
set_guidrec(Q::GuidedProposal, gr) = GuidedProposal(Q.target, Q.aux, Q.tt, Q.xobs0, Q.xobsT, gr, Q.nshapes,Q.mT, Q.endpoint)

"""
    update_mT(Q, mTv, obsinfo)

Update State vector of auxiliary process for each shape.
For the k-th shape, the momentum gets replaced with `mTv[k]`

## Example
See documentation for GuidedProposal to contruct an example instance, say `Q`

```julia
    mTv = [zeros(PointF,5)]

    BridgeLandmarks.update_mT!(Q, mTv, obsinfo)
"""
function update_mT!(Q, mTv, obsinfo)
    for k ∈ 1:Q.nshapes
        Q.aux[k] = auxiliary(Q.target,State(Q.xobsT[k],mTv[k]))  # auxiliary process for each shape
    end
    Q = update_guidrec!(Q, obsinfo)
    Q
end

"""
    construct_gp_xobsT(Q, xobsTᵒ)

Update xobsTᵒ into auxiliary process of Q, following by recomputing the backwards ODEs
"""
function construct_gp_xobsT(Q, xobsTᵒ)
    aux = [auxiliary(Q.target,State(xobsTᵒ[k],Q.mT[k])) for k ∈ 1:Q.nshapes]
    GuidedProposal(Q.target, aux, Q.tt, Q.xobs0, xobsTᵒ, Q.guidrec, Q.nshapes, Q.mT)
end

"""
    _b!((i,t), x::State, out::State, Q::GuidedProposal,k)

Evaluate drift bᵒ of guided proposal at (t,x), write into out
"""
function _b!((i,t), x::State, out::State, Q::GuidedProposal,k)
    Bridge.b!(t, x, out, Q.target)
    out .+= amul(t,x,Q.guidrec[k].Lt[i]' *
        (Q.guidrec[k].Mt[i] *(Q.xobsT[k]-Q.guidrec[k].μt[i]-Q.guidrec[k].Lt[i]*vec(x))),Q.target)
    out
end

"""
    σ!(t, x, dw, out, Q::GuidedProposal)

Evaluate σ(t,x) dw and write into out
"""
σ!(t, x, dw, out, Q::GuidedProposal) = σ!(t, x, dw, out, Q.target)

"""
    _r!((i,t), x::State, out::State, Q::GuidedProposal,k)

Evaluate tilde_r (appearing in guiding term of guided proposal) at (t,x), write into out
"""
function _r!((i,t), x::State, out::State, Q::GuidedProposal,k)
    out .= vecofpoints2state(Q.guidrec[k].Lt[i]' * (Q.guidrec[k].Mt[i] * (Q.xobsT[k]-Q.guidrec[k].μt[i]-Q.guidrec[k].Lt[i]*vec(x))))
    out
end

"""
    lρtilde(x0, Q,k)

Compute log ρ̃(0,x_0,k), where k indexes shape
"""
function lρtilde(x0, Q,k)
  y = deepvec([Q.xobs0; Q.xobsT[k]] - Q.guidrec[k].μt0 - Q.guidrec[k].Lt0*vec(x0))
  M⁺0deep = deepmat(Q.guidrec[k].Mt⁺0)
  -0.5*logdet(M⁺0deep) -0.5*dot(y, M⁺0deep\y)
end

"""
    gp!(::LeftRule,  Xᵒ, x0, W, Q::GuidedProposal,k; skip = sk, ll0 = true)

Simulate guided proposal as specified in `Q` and compute loglikelihood for one shape,
starting from `x0`, using Wiener increments `W`

## Write into
Xᵒ

## Returns
logliklihood.
"""
function gp!(::LeftRule,  Xᵒ, x0, W, Q::GuidedProposal, k, At, xT=auxiliary(Q,k).xT; skip = sk, ll0 = true)
    Pnt = eltype(x0)
    tt =  Xᵒ.tt
    Xᵒ.yy[1] .= deepvalue(x0)
    som::deepeltype(x0)  = 0.
    # initialise objects to write into srout and strout are vectors of Points
    dwiener = dimwiener(Q.target)
    srout = zeros(Pnt, dwiener)
    strout = zeros(Pnt, dwiener)
    x = copy(x0)
    rout = copy(x0)
    bout = copy(x0)
    btout = copy(x0)
    wout = copy(x0)
    if !constdiff(Q)
        #At = Bridge.a((1,0), x0, auxiliary(Q,k), xT)  # auxtimehomogeneous switch
        Bridge.a!((1,0), x0, At, auxiliary(Q,k), xT)  # auxtimehomogeneous switch
        A = zeros(Unc{deepeltype(x0)}, 2Q.target.n,2Q.target.n)
    end
    for i ∈ 1:length(tt)-1
        dt = tt[i+1]-tt[i]
        b!(tt[i], x, bout, target(Q)) # b(t,x)
        _r!((i,tt[i]), x, rout, Q, k) # tilder(t,x)
        σt!(tt[i], x, rout, srout, target(Q))      #  σ(t,x)' * tilder(t,x) for target(Q)
        Bridge.σ!(tt[i], x, srout*dt + W.yy[i+1] - W.yy[i], wout, target(Q)) # σ(t,x) (σ(t,x)' * tilder(t,x) + dW(t))
        # likelihood terms
        if i<=length(tt)-1-skip
            b!(tt[i], x, btout, auxiliary(Q,k),xT)
            som += dot(bout-btout, rout) * dt
            if !constdiff(Q)
                σt!(tt[i], x, rout, strout, auxiliary(Q,k),xT)  #  tildeσ(t,x)' * tilder(t,x) for auxiliary(Q)
                som += 0.5*Bridge.inner(srout) * dt    # |σ(t,x)' * tilder(t,x)|^2
                som -= 0.5*Bridge.inner(strout) * dt   # |tildeσ(t,x)' * tilder(t,x)|^2
                Bridge.a!((i,tt[i]), x, A, target(Q))
                som += 0.5*(dot(At,Q.guidrec[k].Ht[i]) - dot(A,Q.guidrec[k].Ht[i])) * dt
            end
        end
        x .= x + dt * bout + wout
        Xᵒ.yy[i+1] .= deepvalue(x)
    end
    if ll0
        logρ0 = lρtilde(x0,Q,k)
    else
        logρ0 = 0.0 # don't compute
    end
    copyto!(Xᵒ.yy[end], Bridge.endpoint(Xᵒ.yy[end],Q))
    som + logρ0, Xᵒ
end

"""
    gp!(::LeftRule,  X, x0, W, Q::GuidedProposal; skip = sk, ll0 = true)

Simulate guided proposal and compute loglikelihood (vector version, multiple shapes)
`x0` is assumed to be of type `State`
"""
function gp!(::LeftRule,  X::Vector, x0, W, Q::GuidedProposal, At; skip = sk, ll0 = true)
     logliks  = zeros(deepeltype(x0), Q.nshapes)
     for k ∈ 1:Q.nshapes
         logliks[k], X[k] = gp!(LeftRule(), X[k],x0,W[k],Q, k, At ;skip=skip,ll0=ll0)
     end
     logliks, X
end

"""
    gp!(::LeftRule,  X::Vector, q, p , W, Q::GuidedProposal; skip = sk, ll0 = true)

Simulate guided proposal and compute loglikelihood (vector version, multiple shapes)
`q` and `p` make up a state (could in fact be turned into the intial state by `merge_state(q,p)`)
"""
function gp!(::LeftRule,  X::Vector, q, p , W, Q::GuidedProposal, At; skip = sk, ll0 = true)
    T = typeof(q[1] + p[1])
    x0 = NState(reinterpret(Point{T}, T.(q)),reinterpret(Point{T}, T.(p)))
    logliks  = zeros(deepeltype(x0), Q.nshapes)
    for k ∈ 1:Q.nshapes
        logliks[k], X[k] = gp!(LeftRule(), X[k],x0,W[k],Q, k, At ;skip=skip,ll0=ll0)
    end
    logliks, X
end

"""
    getpars(Q::GuidedProposal)

Extract parameters from GuidedProposal `Q`, that is, `(a,c,γ)``
"""
function getpars(Q::GuidedProposal)
    P = Q.target
    [P.a, getγ(P)]
end

"""
    adjust_to_newpars(Q::GuidedProposal,(aᵒ,cᵒ,γᵒ),obsinfo)

Provide new parameter values for GuidedProposal `Q`, these are written into fields `target` and `aux`.
Returns a new instance of `GuidedProposal`, adjusted to the new set of parameters.
"""
function adjust_to_newpars(Q::GuidedProposal,θᵒ, obsinfo, gramT_container)
    (aᵒ,γᵒ) = θᵒ
    if isa(Q.target,MarslandShardlow)
        target = MarslandShardlow(aᵒ, Q.target.c, γᵒ, Q.target.λ, Q.target.n)
    end
    if isa(Q.target,Landmarks)
        if γconstant
            nfs = Q.target.nfs
        else
            nfs = construct_nfs(Q.target.db, Q.target.nfstd, γᵒ)
        end
        target = Landmarks(aᵒ,Q.target.c, Q.target.n, Q.target.db, Q.target.nfstd, nfs)
    end
    for k ∈ eachindex(gramT_container)
         gram_matrix!(gramT_container[k], Q.xobsT[k], Q.aux[k])
    end
    aux = [auxiliary(target,State(Q.xobsT[k],Q.mT[k]), gramT_container[k]) for k ∈ 1:Q.nshapes]
    Qnew = GuidedProposal(target, aux, Q.tt, Q.xobs0, Q.xobsT, Q.guidrec, Q.nshapes, Q.mT)
    Qnew = update_guidrec!(Qnew, obsinfo)
    Qnew
end

"""
    target(Q::GuidedProposal) = Q.target
"""
    target(Q::GuidedProposal) = Q.target

"""
    auxiliary(Q::GuidedProposal,k::Int64) = Q.aux[k]

Extract auxiliary process of k-th shape.
"""
auxiliary(Q::GuidedProposal,k::Int64) = Q.aux[k] # auxiliary process of k-th shape

"""
    constdiff(Q::GuidedProposal)

If true, both the target and auxiliary process have constant diffusion coefficient.
"""
constdiff(Q::GuidedProposal) = constdiff(target(Q)) && constdiff(auxiliary(Q,1))


text:
module TestUtils
    using StatsBase
    using Base.Test
    using NullableArrays

    @testset "describe" begin
        io = IOBuffer()
        describe(io, NullableArray(1:10))
        @test String(take!(io)) == """
                                   Summary Stats:
                                   Mean:           5.500000
                                   Minimum:        1.000000
                                   1st Quartile:   3.250000
                                   Median:         5.500000
                                   3rd Quartile:   7.750000
                                   Maximum:        10.000000
                                   Length:         10
                                   Type:           $Int
                                   Number Missing: 0
                                   % Missing:      0.000000
                                   """
        describe(io, NullableArray([1, Nullable()]))
        @test String(take!(io)) == """
                                   Summary Stats:
                                   Mean:           1.000000
                                   Minimum:        1.000000
                                   1st Quartile:   1.000000
                                   Median:         1.000000
                                   3rd Quartile:   1.000000
                                   Maximum:        1.000000
                                   Length:         1
                                   Type:           $Int
                                   Number Missing: 1
                                   % Missing:      50.000000
                                   """
        describe(io, NullableArray(["s"]))
        @test String(take!(io)) == """
                                   Summary Stats:
                                   Length:         1
                                   Type:           Nullable{String}
                                   Number Unique:  1
                                   Number Missing: 0
                                   % Missing:      0.000000
                                   """
        describe(io, NullableArray(["s", Nullable()]))
        @test String(take!(io)) == """
                                   Summary Stats:
                                   Length:         2
                                   Type:           Nullable{String}
                                   Number Unique:  2
                                   Number Missing: 1
                                   % Missing:      50.000000
                                   """
        describe(io, ["s", Nullable()])
        @test String(take!(io)) == """
                                   Summary Stats:
                                   Length:         2
                                   Type:           Nullable{String}
                                   Number Unique:  2
                                   Number Missing: 1
                                   % Missing:      50.000000
                                   """
        describe(io, [1, Nullable()])
        @test String(take!(io)) == """
                                   Summary Stats:
                                   Mean:           1.000000
                                   Minimum:        1.000000
                                   1st Quartile:   1.000000
                                   Median:         1.000000
                                   3rd Quartile:   1.000000
                                   Maximum:        1.000000
                                   Length:         1
                                   Type:           $Int
                                   Number Missing: 1
                                   % Missing:      50.000000
                                   """
        describe(io, NullableArray{Any}(5))
        @test String(take!(io)) == """
                                   Summary Stats:
                                   Length:         5
                                   Type:           Nullable{Any}
                                   Number Unique:  1
                                   Number Missing: 5
                                   % Missing:      100.000000
                                   """
        describe(io, NullableArray{Float64}(5))
        @test String(take!(io)) == """
                                   Summary Stats:
                                   Type:           Nullable{Float64}
                                   Number Missing: 5
                                   % Missing:      100.000000
                                   """
        describe(io, fill(Nullable{String}(), 5))
        @test String(take!(io)) == """
                                   Summary Stats:
                                   Length:         5
                                   Type:           Nullable{String}
                                   Number Unique:  1
                                   Number Missing: 5
                                   % Missing:      100.000000
                                   """
        describe(io, fill(Nullable{Float64}(), 5))
        @test String(take!(io)) == """
                                   Summary Stats:
                                   Type:           Nullable{Float64}
                                   Number Missing: 5
                                   % Missing:      100.000000
                                   """
    end
end


text:
#=

Code to simulate the model in Globalization and Synchronization of Innovation
Cycles by Kiminori Matsuyama, Laura Gardini and Iryna Sushko

See <http://www.centreformacroeconomics.ac.uk/Discussion-Papers/2015/CFMDP2015-27-Paper.pdf>

=#

using PyPlot

"""
If we expand the implicit function for h_j(n_k) then we find that
it is a quadratic. We know that h_j(n_k) > 0 so we can get its
value by using the quadratic form
"""

function h_j(j::Int64, nk, s1::Float64, s2::Float64, theta::Float64, delta::Float64, rho::Float64)
    # Find out who's h we are evaluating
    if j == 1
        sj = s1
        sk = s2
    else
        sj = s2
        sk = s1
    end

    # Coefficients on the quadratic a x^2 + b x + c = 0
    a = 1.0
    b = ((rho + 1/rho)*nk - sj - sk)
    c = (nk*nk - (sj*nk)/rho - sk*rho*nk)

    # Positive solution of quadratic form
    root = (-b + sqrt(b*b - 4*a*c))/(2*a)

    return root
end

"Determine whether (n1, n2) is in the set DLL"
function DLL(n1::Float64, n2::Float64, s1_rho::Float64, s2_rho::Float64,
        s1::Float64, s2::Float64, theta::Float64, delta::Float64, rho::Float64)
    return (n1 <= s1_rho) && (n2 <= s2_rho)
end

"Determine whether (n1, n2) is in the set DHH"
function DHH(n1::Float64, n2::Float64, s1_rho::Float64, s2_rho::Float64,
        s1::Float64, s2::Float64, theta::Float64, delta::Float64, rho::Float64)
    return (n1 >= h_j(1, n2, s1, s2, theta, delta, rho)) && (n2 >= h_j(2, n1, s1, s2, theta, delta, rho))
end

"Determine whether (n1, n2) is in the set DHL"
function DHL(n1::Float64, n2::Float64, s1_rho::Float64, s2_rho::Float64,
        s1::Float64, s2::Float64, theta::Float64, delta::Float64, rho::Float64)
    return (n1 >= s1_rho) && (n2 <= h_j(2, n1, s1, s2, theta, delta, rho))
end

"Determine whether (n1, n2) is in the set DLH"
function DLH(n1::Float64, n2::Float64, s1_rho::Float64, s2_rho::Float64,
        s1::Float64, s2::Float64, theta::Float64, delta::Float64, rho::Float64)
    return (n1 <= h_j(1, n2, s1, s2, theta, delta, rho)) && (n2 >= s2_rho)
end

"""
Takes a current value for (n_{1, t}, n_{2, t}) and returns the
values (n_{1, t+1}, n_{2, t+1}) according to the law of motion.
"""
function one_step(n1::Float64, n2::Float64, s1_rho::Float64, s2_rho::Float64,
        s1::Float64, s2::Float64, theta::Float64, delta::Float64, rho::Float64)
    # Depending on where we are, evaluate the right branch
    if DLL(n1, n2, s1_rho, s2_rho, s1, s2, theta, delta, rho)
        n1_tp1 = delta*(theta*s1_rho + (1-theta)*n1)
        n2_tp1 = delta*(theta*s2_rho + (1-theta)*n2)
    elseif DHH(n1, n2, s1_rho, s2_rho, s1, s2, theta, delta, rho)
        n1_tp1 = delta*n1
        n2_tp1 = delta*n2
    elseif DHL(n1, n2, s1_rho, s2_rho, s1, s2, theta, delta, rho)
        n1_tp1 = delta*n1
        n2_tp1 = delta*(theta*h_j(2, n1, s1, s2, theta, delta, rho) + (1-theta)*n2)
    elseif DLH(n1, n2, s1_rho, s2_rho, s1, s2, theta, delta, rho)
        n1_tp1 = delta*(theta*h_j(1, n2, s1, s2, theta, delta, rho) + (1-theta)*n1)
        n2_tp1 = delta*n2
    end

    return n1_tp1, n2_tp1
end

"""
Given an initial condition, continues to yield new values of
n1 and n2
"""
function new_n1n2(n1_0::Float64, n2_0::Float64, s1_rho::Float64, s2_rho::Float64,
        s1::Float64, s2::Float64, theta::Float64, delta::Float64, rho::Float64)
    n1_tp1, n2_tp1 = one_step(n1_0, n2_0, s1_rho, s2_rho, s1, s2, theta, delta, rho)
    return n1_tp1, n2_tp1
end

"""
Takes initial values and iterates forward to see whether
the histories eventually end up in sync.

If countries are symmetric then as soon as the two countries have the
same measure of firms then they will by synchronized -- However, if
they are not symmetric then it is possible they have the same measure
of firms but are not yet synchronized. To address this, we check whether
firms stay synchronized for `npers` periods with Euclidean norm

Parameters
----------
n1_0 : scalar(Float)
    Initial normalized measure of firms in country one
n2_0 : scalar(Float)
    Initial normalized measure of firms in country two
maxiter : scalar(Int)
    Maximum number of periods to simulate
npers : scalar(Int)
    Number of periods we would like the countries to have the
    same measure for

Returns
-------
synchronized : scalar(Bool)
    Did they two economies end up synchronized
pers_2_sync : scalar(Int)
    The number of periods required until they synchronized
"""
function pers_till_sync(n1_0::Float64, n2_0::Float64, s1_rho::Float64, s2_rho::Float64,
        s1::Float64, s2::Float64, theta::Float64, delta::Float64, rho::Float64,maxiter::Int64, npers)

    # Initialize the status of synchronization
    synchronized = false
    pers_2_sync = maxiter
    iters = 0

    nsync = 0

    while (~synchronized) && (iters < maxiter)
        # Increment the number of iterations and get next values
        iters += 1

        n1_t, n2_t = new_n1n2(n1_0, n2_0, s1_rho, s2_rho, s1, s2, theta, delta, rho)

        # Check whether same in this period
        if abs(n1_t - n2_t) < 1e-8
            nsync += 1
        # If not, then reset the nsync counter
        else
            nsync = 0
        end

        # If we have been in sync for npers then stop and countries
        # became synchronized nsync periods ago
        if nsync > npers
            synchronized = true
            pers_2_sync = iters - nsync
        end
        n1_0, n2_0 = n1_t, n2_t
    end
    return synchronized, pers_2_sync
end

function create_attraction_basis(s1_rho, s2_rho, s1, s2, theta, delta, rho,
        maxiter, npers, npts)
    # Create unit range with npts
    synchronized, pers_2_sync = false, 0
    unit_range = linspace(0.0, 1.0, npts)

    # Allocate space to store time to sync
    time_2_sync = Array{Float64}(npts, npts)
    # Iterate over initial conditions
    for (i, n1_0) in enumerate(unit_range)
        for (j, n2_0) in enumerate(unit_range)
            #println((i,j))
            synchronized, pers_2_sync = pers_till_sync(n1_0, n2_0, s1_rho, s2_rho,
                                                        s1, s2, theta, delta, rho,
                                                        maxiter, npers)
            time_2_sync[i, j] = pers_2_sync
        end
    end

    return time_2_sync
end

# == Now we define a class for the model == #

"""
The paper "Globalization and Synchronization of Innovation Cycles" presents
a two country model with endogenous innovation cycles. Combines elements
from Deneckere Judd (1985) and Helpman Krugman (1985) to allow for a
model with trade that has firms who can introduce new varieties into
the economy.

We focus on being able to determine whether two countries eventually
synchronize their innovation cycles. To do this, we only need a few
of the many parameters. In particular, we need the parameters listed
below

Parameters
----------
s1 : scalar(Float)
    Amount of total labor in country 1 relative to total worldwide labor
theta : scalar(Float)
    A measure of how mcuh more of the competitive variety is used in
    production of final goods
delta : scalar(Float)
    Percentage of firms that are not exogenously destroyed every period
rho : scalar(Float)
    Measure of how expensive it is to trade between countries
"""
type MSGSync
    s1::Float64
    s2::Float64
    s1_rho::Float64
    s2_rho::Float64
    theta::Float64
    delta::Float64
    rho::Float64
end

function MSGSync(s1::Float64=0.5, theta::Float64=2.5,
                 delta::Float64=0.7, rho::Float64=0.2)
    # Store other cutoffs and parameters we use
    s2 = 1 - s1
    s1_rho = min((s1 - rho*s2) / (1 - rho),1)
    s2_rho = 1 - s1_rho

    model=MSGSync(s1,s2,s1_rho,s2_rho,theta,delta,rho)

    return model
end

function unpack_params(model::MSGSync)
    return model.s1, model.s2, model.theta, model.delta, model.rho, model.s1_rho, model.s2_rho
end

"""
Simulates the values of (n1, n2) for T periods

Parameters
----------
n1_0 : scalar(Float)
Initial normalized measure of firms in country one
n2_0 : scalar(Float)
Initial normalized measure of firms in country two
T : scalar(Int)Number of periods to simulate

Returns
-------
n1 : Array{Float64}(ndim=1)
A history of normalized measures of firms in country one
n2 : Array{Float64}(ndim=1)
A history of normalized measures of firms in country two
"""
function simulate_n(model::MSGSync, n1_0::Float64, n2_0::Float64, T::Int64)
    # Unpack parameters
    s1, s2, theta, delta, rho, s1_rho, s2_rho = unpack_params(model)

    # Allocate space
    n1 = Array{Float64}(T)
    n2 = Array{Float64}(T)

    # Simulate for T periods
    for t in 1:T
        # Get next values
        n1[t], n2[t] = n1_0, n2_0
        n1_0, n2_0 = new_n1n2(n1_0, n2_0, s1_rho, s2_rho, s1, s2, theta, delta, rho)
    end

    return n1, n2
end

"""
Takes initial values and iterates forward to see whether
the histories eventually end up in sync.

If countries are symmetric then as soon as the two countries have the
same measure of firms then they will by synchronized -- However, if
they are not symmetric then it is possible they have the same measure
of firms but are not yet synchronized. To address this, we check whether
firms stay synchronized for `npers` periods with Euclidean norm

Parameters
----------
n1_0 : scalar(Float)
    Initial normalized measure of firms in country one
n2_0 : scalar(Float)
    Initial normalized measure of firms in country two
maxiter : scalar(Int)
    Maximum number of periods to simulate
npers : scalar(Int)
    Number of periods we would like the countries to have the
    same measure for

Returns
-------
synchronized : scalar(Bool)
    Did they two economies end up synchronized
pers_2_sync : scalar(Int)
    The number of periods required until they synchronized
"""
function pers_till_sync(model::MSGSync, n1_0::Float64, n2_0::Float64,
                        maxiter=500, npers=3)
    # Unpack parameters
    s1, s2, theta, delta, rho, s1_rho, s2_rho = unpack_params(model)

    return pers_till_sync(n1_0, n2_0, s1_rho, s2_rho, s1, s2, theta, delta, rho, maxiter, npers)
end

"""
Creates an attraction basis for values of n on [0, 1] X [0, 1] with npts in each dimension
"""
function create_attraction_basis(model::MSGSync; maxiter=250, npers=3, npts=50)
    # Unpack parameters
    s1, s2, theta, delta, rho, s1_rho, s2_rho = unpack_params(model)

    ab = create_attraction_basis(s1_rho, s2_rho, s1, s2, theta, delta,
                                  rho, maxiter, npers, npts)
    return ab
end


text:
@horizontal_average(
    AtmosLESConfigType,
    rho,
    "kg m^-3",
    "air density",
    "air_density",
) do (atmos::AtmosModel, states::States, curr_time, cache)
    states.prognostic.ρ
end

@horizontal_average(
    AtmosLESConfigType,
    u,
    "m s^-1",
    "x-velocity",
    "",
    rho,
) do (atmos::AtmosModel, states::States, curr_time, cache)
    states.prognostic.ρu[1]
end

@horizontal_average(
    AtmosLESConfigType,
    v,
    "m s^-1",
    "y-velocity",
    "",
    rho,
) do (atmos::AtmosModel, states::States, curr_time, cache)
    states.prognostic.ρu[2]
end

@horizontal_average(
    AtmosLESConfigType,
    w,
    "m s^-1",
    "z-velocity",
    "",
    rho,
) do (atmos::AtmosModel, states::States, curr_time, cache)
    states.prognostic.ρu[3]
end

@horizontal_average(
    AtmosLESConfigType,
    uu,
    rho,
) do (atmos::AtmosModel, states::States, curr_time, cache)
    states.prognostic.ρu[1]^2 / states.prognostic.ρ
end

@horizontal_average(
    AtmosLESConfigType,
    vv,
    rho,
) do (atmos::AtmosModel, states::States, curr_time, cache)
    states.prognostic.ρu[2]^2 / states.prognostic.ρ
end

@horizontal_average(
    AtmosLESConfigType,
    ww,
    rho,
) do (atmos::AtmosModel, states::States, curr_time, cache)
    states.prognostic.ρu[3]^2 / states.prognostic.ρ
end

@horizontal_average(
    AtmosLESConfigType,
    www,
    rho,
) do (atmos::AtmosModel, states::States, curr_time, cache)
    states.prognostic.ρu[3]^3 / states.prognostic.ρ^2
end

@horizontal_average(
    AtmosLESConfigType,
    wu,
    rho,
) do (atmos::AtmosModel, states::States, curr_time, cache)
    states.prognostic.ρu[3] * states.prognostic.ρu[1] / states.prognostic.ρ
end

@horizontal_average(
    AtmosLESConfigType,
    wv,
    rho,
) do (atmos::AtmosModel, states::States, curr_time, cache)
    states.prognostic.ρu[3] * states.prognostic.ρu[2] / states.prognostic.ρ
end

@horizontal_average(
    AtmosLESConfigType,
    wrho,
    rho,
) do (atmos::AtmosModel, states::States, curr_time, cache)
    states.prognostic.ρu[3] * states.prognostic.ρ
end

@horizontal_average(
    AtmosLESConfigType,
    temp,
    "K",
    "air temperature",
    "air_temperature",
    rho,
) do (atmos::AtmosModel, states::States, curr_time, cache)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    air_temperature(ts) * states.prognostic.ρ
end

@horizontal_average(
    AtmosLESConfigType,
    pres,
    "Pa",
    "air pressure",
    "air_pressure",
    rho,
) do (atmos::AtmosModel, states::States, curr_time, cache)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    air_pressure(ts) * states.prognostic.ρ
end

@horizontal_average(
    AtmosLESConfigType,
    thd,
    "K",
    "dry potential temperature",
    "air_potential_temperature",
    rho,
) do (atmos::AtmosModel, states::States, curr_time, cache)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    dry_pottemp(ts) * states.prognostic.ρ
end

@horizontal_average(
    AtmosLESConfigType,
    et,
    "J kg^-1",
    "total specific energy",
    "specific_dry_energy_of_air",
    rho,
) do (
    energy::TotalEnergyModel,
    atmos::AtmosModel,
    states::States,
    curr_time,
    cache,
)
    states.prognostic.energy.ρe
end

@horizontal_average(
    AtmosLESConfigType,
    ei,
    "J kg^-1",
    "specific internal energy",
    "internal_energy",
    rho,
) do (atmos::AtmosModel, states::States, curr_time, cache)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    internal_energy(ts) * states.prognostic.ρ
end

@horizontal_average(
    AtmosLESConfigType,
    ht,
    "J kg^-1",
    "specific enthalpy based on total energy",
    "",
    rho,
) do (
    energy::TotalEnergyModel,
    atmos::AtmosModel,
    states::States,
    curr_time,
    cache,
)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    e_tot = states.prognostic.energy.ρe / states.prognostic.ρ
    total_specific_enthalpy(ts, e_tot) * states.prognostic.ρ
end

@horizontal_average(
    AtmosLESConfigType,
    hi,
    "J kg^-1",
    "specific enthalpy based on internal energy",
    "atmosphere_enthalpy_content",
    rho,
) do (atmos::AtmosModel, states::States, curr_time, cache)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    specific_enthalpy(ts) * states.prognostic.ρ
end

@horizontal_average(
    AtmosLESConfigType,
    w_ht_sgs,
    "kg kg^-1 m s^-1",
    "vertical sgs flux of total specific enthalpy",
    "",
    rho,
) do (
    energy::TotalEnergyModel,
    atmos::AtmosModel,
    states::States,
    curr_time,
    cache,
)
    D_t = get!(cache, :D_t) do
        _, D_t, _ = turbulence_tensors(
            atmos,
            states.prognostic,
            states.gradient_flux,
            states.auxiliary,
            curr_time,
        )
        D_t
    end
    d_h_tot = -D_t .* states.gradient_flux.energy.∇h_tot
    d_h_tot[end] * states.prognostic.ρ
end

@horizontal_average(
    AtmosLESConfigType,
    eiei,
    rho,
) do (atmos::AtmosModel, states::States, curr_time, cache)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    internal_energy(ts)^2 * states.prognostic.ρ
end

@horizontal_average(
    AtmosLESConfigType,
    wthd,
    rho,
) do (atmos::AtmosModel, states::States, curr_time, cache)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    states.prognostic.ρu[3] * dry_pottemp(ts)
end

@horizontal_average(
    AtmosLESConfigType,
    wei,
    rho,
) do (atmos::AtmosModel, states::States, curr_time, cache)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    states.prognostic.ρu[3] * internal_energy(ts)
end

@horizontal_average(
    AtmosLESConfigType,
    qt,
    "kg kg^-1",
    "mass fraction of total water in air (qv+ql+qi)",
    "mass_fraction_of_water_in_air",
    rho,
) do (
    moisture::Union{EquilMoist, NonEquilMoist},
    atmos::AtmosModel,
    states::States,
    curr_time::AbstractFloat,
    cache::Dict{Symbol, Any},
)
    states.prognostic.moisture.ρq_tot
end

@horizontal_average(
    AtmosLESConfigType,
    ql,
    "kg kg^-1",
    "mass fraction of liquid water in air",
    "mass_fraction_of_cloud_liquid_water_in_air",
    rho,
) do (
    moisture::Union{EquilMoist, NonEquilMoist},
    atmos::AtmosModel,
    states::States,
    curr_time::AbstractFloat,
    cache::Dict{Symbol, Any},
)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    liquid_specific_humidity(ts) * states.prognostic.ρ
end

@horizontal_average(
    AtmosLESConfigType,
    qi,
    "kg kg^-1",
    "mass fraction of ice in air",
    "mass_fraction_of_cloud_ice_in_air",
    rho,
) do (
    moisture::Union{EquilMoist, NonEquilMoist},
    atmos::AtmosModel,
    states::States,
    curr_time::AbstractFloat,
    cache::Dict{Symbol, Any},
)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    ice_specific_humidity(ts) * states.prognostic.ρ
end

@horizontal_average(
    AtmosLESConfigType,
    qv,
    "kg kg^-1",
    "mass fraction of water vapor in air",
    "specific_humidity",
    rho,
) do (
    moisture::Union{EquilMoist, NonEquilMoist},
    atmos::AtmosModel,
    states::States,
    curr_time::AbstractFloat,
    cache::Dict{Symbol, Any},
)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    vapor_specific_humidity(ts) * states.prognostic.ρ
end

@horizontal_average(
    AtmosLESConfigType,
    thv,
    "K",
    "virtual potential temperature",
    "virtual_potential_temperature",
    rho,
) do (
    moisture::Union{EquilMoist, NonEquilMoist},
    atmos::AtmosModel,
    states::States,
    curr_time::AbstractFloat,
    cache::Dict{Symbol, Any},
)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    virtual_pottemp(ts) * states.prognostic.ρ
end

@horizontal_average(
    AtmosLESConfigType,
    thl,
    "K",
    "liquid-ice potential temperature",
    "",
    rho,
) do (
    moisture::Union{EquilMoist, NonEquilMoist},
    atmos::AtmosModel,
    states::States,
    curr_time::AbstractFloat,
    cache::Dict{Symbol, Any},
)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    liquid_ice_pottemp(ts) * states.prognostic.ρ
end

@horizontal_average(
    AtmosLESConfigType,
    w_qt_sgs,
    "kg kg^-1 m s^-1",
    "vertical sgs flux of total specific humidity",
    "",
    rho,
) do (
    moisture::Union{EquilMoist, NonEquilMoist},
    atmos::AtmosModel,
    states::States,
    curr_time::AbstractFloat,
    cache::Dict{Symbol, Any},
)
    D_t = get!(cache, :D_t) do
        _, D_t, _ = turbulence_tensors(
            atmos,
            states.prognostic,
            states.gradient_flux,
            states.auxiliary,
            curr_time,
        )
        D_t
    end
    d_q_tot = -D_t .* states.gradient_flux.moisture.∇q_tot
    d_q_tot[end] * states.prognostic.ρ
end

@horizontal_average(
    AtmosLESConfigType,
    qtqt,
    rho,
) do (
    moisture::Union{EquilMoist, NonEquilMoist},
    atmos::AtmosModel,
    states::States,
    curr_time::AbstractFloat,
    cache::Dict{Symbol, Any},
)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    states.prognostic.moisture.ρq_tot^2 / states.prognostic.ρ
end

@horizontal_average(
    AtmosLESConfigType,
    thlthl,
    rho,
) do (
    moisture::Union{EquilMoist, NonEquilMoist},
    atmos::AtmosModel,
    states::States,
    curr_time::AbstractFloat,
    cache::Dict{Symbol, Any},
)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    liquid_ice_pottemp(ts)^2 * states.prognostic.ρ
end

@horizontal_average(
    AtmosLESConfigType,
    wqt,
    rho,
) do (
    moisture::Union{EquilMoist, NonEquilMoist},
    atmos::AtmosModel,
    states::States,
    curr_time::AbstractFloat,
    cache::Dict{Symbol, Any},
)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    q_tot = states.prognostic.moisture.ρq_tot / states.prognostic.ρ
    states.prognostic.ρu[3] * q_tot
end

@horizontal_average(
    AtmosLESConfigType,
    wql,
    rho,
) do (
    moisture::Union{EquilMoist, NonEquilMoist},
    atmos::AtmosModel,
    states::States,
    curr_time::AbstractFloat,
    cache::Dict{Symbol, Any},
)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    states.prognostic.ρu[3] * liquid_specific_humidity(ts)
end

@horizontal_average(
    AtmosLESConfigType,
    wqi,
    rho,
) do (
    moisture::Union{EquilMoist, NonEquilMoist},
    atmos::AtmosModel,
    states::States,
    curr_time::AbstractFloat,
    cache::Dict{Symbol, Any},
)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    states.prognostic.ρu[3] * ice_specific_humidity(ts)
end

@horizontal_average(
    AtmosLESConfigType,
    wqv,
    rho,
) do (
    moisture::Union{EquilMoist, NonEquilMoist},
    atmos::AtmosModel,
    states::States,
    curr_time::AbstractFloat,
    cache::Dict{Symbol, Any},
)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    states.prognostic.ρu[3] * vapor_specific_humidity(ts)
end

@horizontal_average(
    AtmosLESConfigType,
    wthv,
    rho,
) do (
    moisture::Union{EquilMoist, NonEquilMoist},
    atmos::AtmosModel,
    states::States,
    curr_time::AbstractFloat,
    cache::Dict{Symbol, Any},
)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    states.prognostic.ρu[3] * virtual_pottemp(ts)
end

@horizontal_average(
    AtmosLESConfigType,
    wthl,
    rho,
) do (
    moisture::Union{EquilMoist, NonEquilMoist},
    atmos::AtmosModel,
    states::States,
    curr_time::AbstractFloat,
    cache::Dict{Symbol, Any},
)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    states.prognostic.ρu[3] * liquid_ice_pottemp(ts)
end

@horizontal_average(
    AtmosLESConfigType,
    qtthl,
    rho,
) do (
    moisture::Union{EquilMoist, NonEquilMoist},
    atmos::AtmosModel,
    states::States,
    curr_time::AbstractFloat,
    cache::Dict{Symbol, Any},
)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    states.prognostic.moisture.ρq_tot * liquid_ice_pottemp(ts)
end

@horizontal_average(
    AtmosLESConfigType,
    qtei,
) do (
    moisture::Union{EquilMoist, NonEquilMoist},
    atmos::AtmosModel,
    states::States,
    curr_time::AbstractFloat,
    cache::Dict{Symbol, Any},
)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    states.prognostic.moisture.ρq_tot * internal_energy(ts)
end

#= TODO
    Variables["cld_frac"] = DiagnosticVariable(
        "cld_frac",
        diagnostic_var_attrib(
            "",
            "cloud fraction",
            "cloud_area_fraction_in_atmosphere_layer",
        ),
    )
    Variables["cld_cover"] = DiagnosticVariable(
        "cld_cover",
        diagnostic_var_attrib("", "cloud cover", "cloud_area_fraction"),
    )
    Variables["cld_top"] = DiagnosticVariable(
        "cld_top",
        diagnostic_var_attrib("m", "cloud top", "cloud_top_altitude"),
    )
    Variables["cld_base"] = DiagnosticVariable(
        "cld_base",
        diagnostic_var_attrib("m", "cloud base", "cloud_base_altitude"),
    )
    Variables["lwp"] = DiagnosticVariable(
        "lwp",
        diagnostic_var_attrib(
            "kg m^-2",
            "liquid water path",
            "atmosphere_mass_content_of_cloud_condensed_water",
        ),
    )
=#


text: calc_I_M(reactor, R_0, B_0) = K_I(reactor) / ( B_0 * sqrt(R_0) )


text:
Base.utf16(glypharray::Array{GLGlyph{Uint16}}) = utf16(Uint16[c.glyph for c in glypharray])
Base.utf8(glypharray::Array{GLGlyph{Uint16}})  = utf8(Uint8[uint8(c.glyph) for c in glypharray])

function escape_regex(x::String)
  result = ""
  for elem in x
      if elem in regex_literals
          result *= string('\\')
      end
      result *= string(elem)
  end
  result
end
regreduce(arr, prefix="(", suffix=")") = Regex(reduce((v0, x) -> v0*"|"*prefix*escape_regex(x)*suffix, prefix*escape_regex(arr[1])*suffix, arr[2:end]))

function update_groups!{T}(textGPU::Texture{GLGlyph{T}, 4, 2}, regexs::Dict{T, Regex}, start=1, stop=length(text_array))
  textRam = textGPU.data[start:stop]
  text    = utf8(map(textRam) do x
    char(x.glyph)
  end)
  for (group,regex) in regexs
      for match in matchall(regex, text)
        startorigin = match.offset+1
        stoporigin  = match.offset+match.endof
        setindex1D!(textGPU, group, startorigin:stoporigin, 4) # Set group
      end
  end
end

function update_glyphpositions!{T}(text_array::AbstractArray{GLGlyph{T}}, start=1, stop=length(text_array))
  line = text_array[start].line
  row  = text_array[start].row
  for i=1:stop
    glyph = text_array[i].glyph
    setindex1D!(text_array, T[line, row], i, 2:3)
    if glyph == '\n'
      row = zero(T)
      line += one(T)
    else
      row += one(T)
    end
  end
end
function update_glyphpositions!{T}(text_array::Texture{GLGlyph{T}, 4, 2}, start=1, stop=length(text_array))
  textarray = data(text_array)
  line = textarray[start].line
  row  = textarray[start].row
  for i=1:stop-1
    glyph = textarray[i].glyph
    setindex1D!(textarray, T[line, row], i, 2:3)
    if glyph == '\n'
      row = 0
      line += 1
    else
      row += 1
    end
  end
  text_array[1:end, 1:end] = textarray
end
function makedisplayable(text::String, tab=3)
  result = map(collect(text)) do x
    str = string(x)
    if !is_valid_utf8(str)
      return utf8([one(Uint8)]) # replace with something that yields a missing symbol
    elseif str == "\r"
      return "\n"
    else
      return str == "\t" ? utf8(" "^tab) : utf8(str) # also replace tabs
    end
  end
  join(result)
end

function toglypharray(text::String, tab=3)
  #@assert is_valid_utf16(text) # future support for utf16
  text = makedisplayable(text,tab)
  #Allocate some more memory, to reduce growing the texture residing on VRAM
  texturesize = div(length(text),     1024)+1 # a texture size of 1024 should be supported on every GPU
  text_array  = Array(GLGlyph{Uint16}, 1024, texturesize)
  setindex1D!(text_array, 1, 1, 2) # set first line
  setindex1D!(text_array, 0, 1, 3) # set first row
  #Set text
  for (i, elem) in enumerate(text)
    setindex1D!(text_array, uint16(char(elem)), i, 1)
    setindex1D!(text_array, 0, i, 4)
  end
  update_glyphpositions!(text_array) # calculate glyph positions
  text_array
end

operators = [":", ";","=", "+", "-", "!", "¬", "~", "<", ">","=", "/", "&", "|", "\$", "*"]
brackets  = ["(", ")", "[", "]", "{", "}"]
keywords  = ["for", "end", "while", "if", "elseif", "using", "return", "in", "function", "local", "global", "let", "quote", "begin", "const", "do", "false", "true"]
regex_literals = ['|', '[', ']', '*', '.', '?', '\\', '(', ')', '{', '}', '+', '-', '$']

julia_groups = @compat Dict(
  1 => regreduce(operators),
  2 => regreduce(brackets),
  3 => regreduce(keywords, "((?<![[:alpha:]])", "(?![[:alpha:]]))"),
  4 => r"(#=.*=#)|(#.*[\n\r])", #Comments
  5 => r"(\".*\")|('.*')|((?<!:):[[:alpha:]][[:alpha:]_]*)", #String alike
  6 => r"(?<![[:alpha:]])[[:graph:]]*\(" # functions
)


text:
import MaximumLikelihoodProblems

import Documenter
import Literate
import Random

_this_filename = @__FILE__
_docs = dirname(_this_filename)
root = dirname(_docs)
docs = joinpath(root, "docs")
docs_src = joinpath(root, "docs", "src")
docs_src_examples = joinpath(root, "docs", "src", "examples")
examples_directory = joinpath(root, "examples")

examples_list = ["linear_regression",
                 "logistic_regression",
                 "multinomial_logistic_regression"]

for example in examples_list
    example_filename = joinpath(examples_directory,
                                "$(example).jl")
    Literate.markdown(example_filename,
                      docs_src_examples)
end

Random.seed!(123)

generated_examples = [joinpath("examples", "$(example).md") for example in examples_list]

Documenter.makedocs(;
    root = docs,
    modules = [MaximumLikelihoodProblems],
    format = Documenter.HTML(),
    pages = vcat(["Home" => "index.md",],
                  generated_examples,
                  ["api_public.md",]),
    repo = "https://github.com/bcbi/MaximumLikelihoodProblems.jl/blob/{commit}{path}#L{line}",
    sitename = "MaximumLikelihoodProblems.jl",
    authors = "Dilum P. Aluthge",
)

Documenter.deploydocs(;
    repo="github.com/bcbi/MaximumLikelihoodProblems.jl",
)


text:
# This file was generated, do not modify it. # hide
unimportant_features = filter(:importance => <(0.005), feature_importance_table).feature

pipe2 = ContinuousEncoder() |>
    FeatureSelector(features=unimportant_features, ignore=true) |> booster

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
using Test
import MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOU = MOI.Utilities

@testset "Start value attributes and nothing" begin
    model = MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()))
    MOI.Test.start_values_unset_test(model, MOI.Test.TestConfig())
end

@testset "Failed copy" begin
    include("../dummy.jl")
    model = DummyModelWithAdd()
    MOIT.failcopytestc(model)
    MOIT.failcopytestia(model)
    MOIT.failcopytestva(model)
    MOIT.failcopytestca(model)
end


text:
module LoadCSV

include("../Types.jl")
using Dates, CSV

export load_csv

timescalingfactor = 1000000 # Convert ns to ms.

function load_csv(output:: Channel)
    for row in CSV.File("~/_dev/pedometer/data/indoor_100steps.csv")
        if row[1] === "T"
            continue
        end
        T = row[1] / timescalingfactor
        x = row[2]
        y = row[3]
        z = row[4]
        magnitude = x^2 + y^2 + z^2
        dp = Types.Datapoint(T, magnitude)
        put!(output, dp)
    end
end

end

max_stars_repo_name: eliabiraschi/pedometer

max_issues_repo_name: eliabiraschi/pedometer

max_forks_repo_name: eliabiraschi/pedometer


text:
include("SortSparse.jl")
include("KoLesky.jl")
include("Cholesky.jl")
using LinearAlgebra
using RecursiveArrayTools
using Distributions
using Distances
using Random
using Profile
using JLD
using MAT
using BlockArrays
Random.seed!(123)

NTrain = 10000
NTest = 5
N = NTest + NTrain
d =  2
y = randn(N)
x = randn(d, N)

setSplit = randperm(N)
trainInds = setSplit[1:NTrain]
testInds = setSplit[(NTrain + 1) : end]

xTrain = x[ :, trainInds ]
yTrain = y[ trainInds ]
xTest = x[:, testInds ]
yTest = y[ testInds ]

n = 10
h = 1 / (n - 1)

xTrain = Matrix(VectorOfArray([[i * h; j * h] for i in 1:n, j in 1:n][:]))
NTrain = size(xTrain,2)
yTrain = randn(NTrain)
N = NTrain + NTest

l = 1.0
σ = 1e-15
cov = r -> exp(-r/l) # + Float64(r <= 1e-10)

rhovals = [2.0, 3.0, 4.0, 5.0]

ρ = 3.0
λ = 1.3

# reordering points and forming skeletons
@time colptr, rowval, P, revP, distances = sortSparseRev(xTrain,ρ, 1)
xTrain = xTrain[:, P]
yTrain = yTrain[P]
x = hcat(xTest, xTrain)
skeletons = construct_skeletons(colptr, rowval, distances, λ)

# Creating
Θ = cov.(pairwise(Euclidean(), x, dims=2))
L_Θ = cholesky(Matrix(Θ))
y = L_Θ.L * randn(N)
Θ = Θ + σ * I

Θ = PseudoBlockArray(Θ, [NTest, NTrain], [NTest, NTrain])
y = PseudoBlockArray(y, [NTest, NTrain])
y[Block(2)] += √(σ) * randn(NTrain)

L = assembleL(skeletons, xTrain, cov)
U = sparse(triu(L * L') + I / max(σ, 1e-30))
icholU_high_level!(U)

μTrue = Θ[Block(1,2)] * (Θ[Block(2,2)] \ y[Block(2)])

# Direct application of Vecchia
# μ = Θ[Block(1,2)] * (L * (L' * y[Block(2)]))

# Version with Cholesky factor
μ = Θ[Block(1,2)] * (σ \ (  U \ (U' \ (L * (L' * y[Block(2)])))))

display(μTrue)
display(μ)
display(y[Block(1)])

@show norm(μTrue - μ) / norm(μTrue)


text:
# Autogenerated wrapper script for gocloc_jll for x86_64-apple-darwin
export gocloc

JLLWrappers.@generate_wrapper_header("gocloc")
JLLWrappers.@declare_executable_product(gocloc)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        gocloc,
        "bin/gocloc",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Tests for Revenue Hölder DEA Models
@testset "RevenueHölderDEAModel" begin

    # Test using Book data
    X = [1; 1; 1; 1; 1; 1; 1; 1];
    Y = [7 7; 4 8; 8 4; 3 5; 3 3; 8 2; 6 4; 1.5 5];
    P = [1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1];

    # Revenue Hölder L1
    revenueholderl1 = dearevenueholder(X, Y, P, l = 1)

    @test typeof(revenueholderl1) == RevenueHolderDEAModel

    @test nobs(revenueholderl1) == 8
    @test ninputs(revenueholderl1) == 1
    @test noutputs(revenueholderl1) == 2
    @test ismonetary(revenueholderl1) == false

    @test efficiency(revenueholderl1) == efficiency(revenueholderl1, :Economic)
    @test efficiency(revenueholderl1, :Economic)   ≈ [0.0; 2.0; 2.0; 6.0; 8.0; 4.0; 4.0; 7.5]
    @test efficiency(revenueholderl1, :Technical)  ≈ [0.0; 0.0; 0.0; 3.0; 5.0; 0.0; 2.0; 3.0]
    @test efficiency(revenueholderl1, :Allocative) ≈ [0.0; 2.0; 2.0; 3.0; 3.0; 4.0; 2.0; 4.5]
    @test normfactor(revenueholderl1) ≈ [1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0]

    # With monetary = true
    revenueholderl1mom = dearevenueholder(X, Y, P, l = 1, monetary = true)

    @test ismonetary(revenueholderl1mom) == true

    @test efficiency(revenueholderl1mom) ≈ efficiency(revenueholderl1) .* normfactor(revenueholderl1)
    @test efficiency(revenueholderl1mom, :Technical) ≈ efficiency(revenueholderl1, :Technical) .* normfactor(revenueholderl1)
    @test efficiency(revenueholderl1mom, :Allocative) ≈ efficiency(revenueholderl1, :Allocative) .* normfactor(revenueholderl1)

    # With CRS (Same result in this example)
    revenueholderl1crs = dearevenueholder(X, Y, P, l = 1, rts = :CRS)

    @test efficiency(revenueholderl1crs, :Economic)   ≈ [0.0; 2.0; 2.0; 6.0; 8.0; 4.0; 4.0; 7.5]
    @test efficiency(revenueholderl1crs, :Technical)  ≈ [0.0; 0.0; 0.0; 3.0; 5.0; 0.0; 2.0; 3.0]
    @test efficiency(revenueholderl1crs, :Allocative) ≈ [0.0; 2.0; 2.0; 3.0; 3.0; 4.0; 2.0; 4.5]

    # Revenue Hölder L1 Weighted (weakly)
    revenueholderl1weight = dearevenueholder(X, Y, P, l = 1, weight = true)

    @test efficiency(revenueholderl1weight, :Economic)   ≈ [0.0; 0.25; 0.25; 1.2; 8/3; 0.5; 2/3; 1.5]
    @test efficiency(revenueholderl1weight, :Technical)  ≈ [0.0; 0.0; 0.0; 0.6; 5/3; 0.0; 1/3; 0.6]
    @test efficiency(revenueholderl1weight, :Allocative) ≈ [0.0; 0.25; 0.25; 0.6; 1.0; 0.5; 1/3; 0.9]

    # Revenue Hölder L2
    logs, value = Test.collect_test_logs() do
        revenueholderl2 = dearevenueholder(X, Y, P, l = 2, optimizer = DEAOptimizer(:NLP))

        @test normfactor(revenueholderl2) ≈ sqrt(2) .* ones(8)
        show(IOBuffer(), revenueholderl2)
    end

    # Revenue Hölder L2 Weighted (weakly)
    logs, value = Test.collect_test_logs() do
        revenueholderl2weight = dearevenueholder(X, Y, P, l = 2, weight = true, optimizer = DEAOptimizer(:NLP))

        @test normfactor(revenueholderl2weight) ≈ sqrt.(sum((P .* Y).^2, dims = 2))
        show(IOBuffer(), revenueholderl2weight)
    end

    # Revenue Hölder LInf
    revenueholderlInf = dearevenueholder(X, Y, P, l = Inf)

    @test efficiency(revenueholderlInf, :Economic)   ≈ [0.0; 1.0; 1.0; 3.0; 4.0; 2.0; 2.0; 3.75]
    @test efficiency(revenueholderlInf, :Technical)  ≈ [0.0; 0.0; 0.0; 2.5; 4.0; 0.0; 1.5; 2.875]
    @test efficiency(revenueholderlInf, :Allocative) ≈ [0.0; 1.0; 1.0; 0.5; 0.0; 2.0; 0.5; 0.875]
    @test normfactor(revenueholderlInf) ≈ [2.0; 2.0; 2.0; 2.0; 2.0; 2.0; 2.0; 2.0]

    # Revenue Hölder LInf Weighted (weakly)
    revenueholderlInfweight = dearevenueholder(X, Y, P, l = Inf, weight = true)

    @test efficiency(revenueholderlInfweight, :Economic)   ≈ [0.0; 1/6; 1/6; 0.75; 4/3; 0.4; 0.4; 15/13]
    @test efficiency(revenueholderlInfweight, :Technical)  ≈ [0.0; 0.0; 0.0; 5/9; 4/3; 0.0; 3/11; 0.6]
    @test efficiency(revenueholderlInfweight, :Allocative) ≈ [0.0; 1/6; 1/6; 7/36; 0.0; 0.4; 7/55; 36/65]
    @test normfactor(revenueholderlInfweight) ≈ [14.0; 12.0; 12.0; 8.0; 6.0; 10.0; 10.0; 6.5]

    # Print
    show(IOBuffer(), revenueholderl1)
    show(IOBuffer(), revenueholderl1weight)
    show(IOBuffer(), revenueholderlInf)
    show(IOBuffer(), revenueholderlInfweight)

    # Test errors
    @test_throws DimensionMismatch dearevenueholder([1; 2 ; 3], [4 ; 5], [1; 1; 1], l = 1) #  Different number of observations
    @test_throws DimensionMismatch dearevenueholder([1; 2; 3], [4; 5; 6], [1; 2; 3; 4], l = 1) # Different number of observation in prices
    @test_throws DimensionMismatch dearevenueholder([1 1; 2 2; 3 3 ], [4; 5; 6], [1 1 1; 2 2 2; 3 3 3], l = 1) # Different number of input prices and inputs
    @test_throws ArgumentError dearevenueholder([1; 2; 3], [4; 5; 6], [1; 2; 3], l = 5) # Invalid l
    @test_throws ArgumentError dearevenueholder([1; 2; 3], [4; 5; 6], [1; 2; 3], l = 1, rts = :Error) # Invalid returns to scale

end


text:
#========== cost "model" ==========#

const BLOCK = Ref(2^18)
# matmul: crossover about 70x70 on my laptop, 70^3 = 343_000, log2(70^3) = 18.3, but only 30% effect at 100^3=10^6
# batchmul: crossover between 20 & 30, log2(20^4) == 17.3, log2(30^4) == 19.6
# contract01: 1500 * 100, length 15_000, doesn't want threading
# cosine01: block 65_536, not sure if it wants
# log: vector crossover about length 10_000

"""
    COSTS = Dict(:* => 0, :log =>10, ...)

Initial cost is `1`, and every other function call adds the value from this dictionary.
Then `n = BLOCK[] ÷ cost` is the number of iterations at which the macro thinks it
worthwhile to turn on threading; you can override this with keyword `threads=n`.
"""
const COSTS = Dict(:+ => 0, :- => 0, :* => 0,
    :conj => 0, :adjoint => 0, :abs =>0, abs2 => 0,
    :getindex => 0, :getproperty => 0, :getfield => 0,
    :^ => 2, :/ => 2, :div =>2, :rem =>2, :mod =>2,
    :log => 10, :exp => 10) # and all others 10, plus 1 initially

callcost(sy, store) = store.cost += get(COSTS, sy, 10)

#========== runtime functions ==========#

"""
    threader(f!,T, Z, (A,B), (1:5,1:6), (1:7), +, block=100, keep=nothing)

Calling `f!(T, Z,A,B, 1:5,1:6, 1:7, nothing)` should do the work.
But if there are enough elements (meaning `5*6*7 > 100`)
then this will call `f!` many times in different threads.
(`block=nothing` turns the whole thing off.)

The first tuple of ranges are supposed to be safe to thread over,
probably the axes of the output `Z`.
It will subdivide the longest until either there are too few elements,
or it has spent its spawning budget, `nthreads()`.

For a scalar reduction the first tuple will be empty, and `length(Z)==1`.
Then it divides up the other axes, each accumulating in its own copy of `Z`.

`keep=nothing` means that it overwrites the array, anything else (`keep=true`) adds on.
"""
@inline function threader(fun!::F, ::Type{T}, Z::AbstractArray, As::Tuple, I0s::Tuple, J0s::Tuple, redfun, block, keep=nothing) where {F <: Function, T}
    if isnothing(block) # then threading is disabled
        fun!(T, Z, As..., I0s..., J0s..., keep)
        return nothing
    elseif !all(r -> r isa AbstractUnitRange, I0s) || !all(r -> r isa AbstractUnitRange, J0s)
        # don't thread ranges like 10:-1:1, and disable @avx too
        fun!(Array, Z, As..., I0s..., J0s..., keep)
        return nothing
    end

    Is = map(UnitRange, I0s)
    Js = map(UnitRange, J0s)
    Ielements = productlength(Is)
    Jelements = productlength(Js)
    threads = min(Threads.nthreads(), cld(Ielements * Jelements, block), Ielements)

    if length(Is) >= 1 && threads>1
        thread_halves(fun!, T, (Z, As...), Is, Js, threads, keep)
    else
        tile_halves(fun!, T, (Z, As...), Is, Js, keep)
    end
    nothing
end

"""
    ∇threader(f!,T, (dA,dB,dZ,A,B), (1:5), (1:6,1:7), block)

Again, calling `f!(T, dA,dB,dZ,A,B, 1:5,1:6, 1:7)` should do the work.

The first tuple of ranges should be safe to thread over, e.g. those in common
to all output arrays.

If there are none, then it should to take a second strategy
of dividing up the other ranges into tiles disjoint in every index,
and giving those to different threads. But this was only right for 2 indices,
and is now disabled.
"""
function ∇threader(fun!::F, ::Type{T}, As::Tuple, I0s::Tuple, J0s::Tuple, block) where {F <: Function, T}
    if isnothing(block) # then threading is disabled
        fun!(T, As..., I0s..., J0s...)
        return nothing
    elseif !all(r -> r isa AbstractUnitRange, I0s) || !all(r -> r isa AbstractUnitRange, J0s)
        # don't thread ranges like 10:-1:1, and disable @avx too
        fun!(Array, As..., I0s..., J0s...)
        return nothing
    end

    Is = map(UnitRange, I0s)
    Js = map(UnitRange, J0s)
    Ielements = productlength(Is)
    Jelements = productlength(Js)
    threads = min(Threads.nthreads(), cld(Ielements * Jelements, block), Ielements)

    if threads > 1
        thread_halves(fun!, T, As, Is, Js, threads)
    else
        tile_halves(fun!, T, As, Is, Js)
    end
    nothing
end

function thread_halves(fun!::F, ::Type{T}, As::Tuple, Is::Tuple, Js::Tuple, threads::Int, keep=nothing) where {F <: Function, T}
    if threads > 2 && rem(threads,3) == 0 # not always halves!
        I1s, I2s, I3s = trisect(Is)
        task1 = Threads.@spawn begin
            thread_halves(fun!, T, As, I1s, Js, threads÷3, keep)
        end
        task2 = Threads.@spawn begin
            thread_halves(fun!, T, As, I2s, Js, threads÷3, keep)
        end
        thread_halves(fun!, T, As, I3s, Js, threads÷3, keep)
        wait(task1)
        wait(task2)
    elseif threads > 1
        I1s, I2s = cleave(Is, maybe32divsize(T))
        task = Threads.@spawn begin
            thread_halves(fun!, T, As, I1s, Js, threads÷2, keep)
        end
        thread_halves(fun!, T, As, I2s, Js, threads÷2, keep)
        wait(task)
    else
        tile_halves(fun!, T, As, Is, Js, keep)
    end
    nothing
end

function tile_halves(fun!::F, ::Type{T}, As::Tuple, Is::Tuple, Js::Tuple, keep=nothing, final=true) where {F <: Function, T}
    # keep == nothing || keep == true || error("illegal value for keep")
    # final == nothing || final == true || error("illegal value for final")
    maxI, maxJ = maximumlength(Is), maximumlength(Js)
    maxL = tile_maxiter(T)
    if maxI < maxL && maxJ < maxL
        fun!(T, As..., Is..., Js..., keep, final)
    elseif maxI > maxJ
        I1s, I2s = cleave(Is)
        tile_halves(fun!, T, As, I1s, Js, keep, final)
        tile_halves(fun!, T, As, I2s, Js, keep, final)
    else
        J1s, J2s = cleave(Js)
        tile_halves(fun!, T, As, Is, J1s, keep, nothing)
        tile_halves(fun!, T, As, Is, J2s, true, final)
    end
    nothing
end

"""
    TILE[] = $(TILE[])
    tile_maxiter(Array{Float64}) == 64?

This now sets the maximum length of iteration of any index,
before it gets broken in half to make smaller tiles.
`TILE[]` is in bytes.
"""
const TILE = Ref(512) # this is now a length, in bytes!

function tile_maxiter(::Type{<:AbstractArray{T}}) where {T}
    isbitstype(T) || return TILE[] ÷ 8
    max(TILE[] ÷ sizeof(T), 4)
end
tile_maxiter(::Type{AT}) where {AT} = TILE[] ÷ 8 # treat anything unkown like Float64

#=

using Tullio
Z = zeros(Int, 11,9);
cnt = 0
f!(::Type, Z, i, j, ♻, 💀) = begin
    global cnt
    Z[i,j] .= (global cnt+=1)
end
Tullio.tile_halves(f!, Array, (Z,), UnitRange.(axes(Z)), (), 4, nothing, true)
Z

  1   1   3   3   5   5   7   7   7
  1   1   3   3   5   5   7   7   7
  2   2   4   4   6   6   8   8   8
  2   2   4   4   6   6   8   8   8
  9   9  10  10  13  13  14  14  14
  9   9  10  10  13  13  14  14  14
  9   9  10  10  13  13  14  14  14
 11  11  11  11  15  15  16  16  16
 11  11  11  11  15  15  16  16  16
 12  12  12  12  15  15  16  16  16
 12  12  12  12  15  15  16  16  16

using TiledIteration
function colour!(A, n=1)
    for (i,t) in enumerate(TileIterator(axes(A), ntuple(_->n, ndims(A))))
        A[t...] .= i
    end
    A
end;
colour!(zeros(Int, 11,9), 2)

 1  1   7   7  13  13  19  19  25
 1  1   7   7  13  13  19  19  25
 2  2   8   8  14  14  20  20  26
 2  2   8   8  14  14  20  20  26
 3  3   9   9  15  15  21  21  27
 3  3   9   9  15  15  21  21  27
 4  4  10  10  16  16  22  22  28
 4  4  10  10  16  16  22  22  28
 5  5  11  11  17  17  23  23  29
 5  5  11  11  17  17  23  23  29
 6  6  12  12  18  18  24  24  30

=#

#========== scalar case ==========#

"""
    thread_scalar(f,T, z, (A,B), (1:5,1:6), +, block=100, keep=nothing)

Just like `threader`, but doesn't take any safe indices `Is`.
And `f` doesn't actually mutate anything, it returns the value...
and `z` is now a scalar, the `init` value for the reduction.
"""
@inline function thread_scalar(fun!::F, ::Type{T}, z::RT, As::Tuple, J0s::Tuple, redfun, block, keep=nothing)::RT where {F <: Function, T, RT}
    if isnothing(block) # then threading is disabled
        return fun!(T, z, As..., J0s..., keep)
    elseif !all(r -> r isa AbstractUnitRange, J0s)
        # don't thread ranges like 10:-1:1, and disable @avx too
        return fun!(Array, z, As..., J0s..., keep)
    end

    Js = map(UnitRange, J0s)
    Jelements = productlength(Js)
    threads = min(Threads.nthreads(), cld(Jelements, block), Jelements)

    if threads < 2
        return fun!(T, z, As..., Js..., keep)
    else
        return scalar_halves(fun!, T, z, As, Js, redfun, threads, keep)
    end
end

function scalar_halves(fun!::F, ::Type{T}, z::RT, As::Tuple, Js::Tuple, redfun, threads, keep=nothing)::RT where {F <: Function, T, RT}
    if threads < 1
        return fun!(T, z, As..., Js..., keep)
    else
        J1s, J2s = cleave(Js)
        S1 = z # scope
        task = Threads.@spawn begin
            S1 = scalar_halves(fun!, T, z, As, J1s, redfun, threads÷2, nothing)
        end
        S2 = scalar_halves(fun!, T, z, As, J2s, redfun, threads÷2, keep)
        wait(task)
        return redfun(S1, S2)
    end
end

#========== tuple functions ==========#

@inline productlength(Is::Tuple) = prod(length.(Is))
@inline productlength(Is::Tuple, Js::Tuple) = productlength(Is) * productlength(Js)

@inline maximumlength(Is::Tuple) = max(length.(Is)...)
@inline maximumlength(::Tuple{}) = 0

@inline maybe32divsize(::Type{<:AbstractArray{T}}) where T<:Number = max(1, 32 ÷ sizeof(T))
@inline maybe32divsize(::Type) = 4

"""
    cleave((1:10, 1:20, 5:15)) -> lo, hi
Picks the longest of a tuple of ranges, and divides that one in half.
"""
@inline cleave(::Tuple{}, n::Int=4) = (), ()
@inline function cleave(ranges::Tuple{UnitRange}, step::Int=4)
    r1 = first(ranges)
    cleft = findcleft(r1, step)
    tuple(first(r1):cleft), tuple(cleft+1:last(r1))
end
@inline function cleave(ranges::Tuple{UnitRange,UnitRange}, step::Int=4)
    r1, r2 = ranges
    if length(r1) > length(r2)
        cleft = findcleft(r1, step)
        return tuple(first(r1):cleft, r2), tuple(cleft+1:last(r1), r2)
    else
        cleft = findcleft(r2, step)
        return tuple(r1, first(r2):cleft), tuple(r1, cleft+1:last(r2))
    end
end
@inline @generated function cleave(ranges::Tuple{Vararg{<:UnitRange,N}}, step::Int=4) where {N}
    ex_finds = [quote
        li = length(ranges[$i])
        if li>l
            c = $i
            l = li
        end
    end for i in 1:N]
    ex_alpas = [:($i==c ? (first(ranges[$i]):cleft) : (ranges[$i])) for i in 1:N]
    ex_betas = [:($i==c ? (cleft+1:last(ranges[$i])) : (ranges[$i])) for i in 1:N]
    quote
        c, l = 0, 0
        $(ex_finds...)
        cleft = findcleft(ranges[c], step)
        tuple($(ex_alpas...)), tuple($(ex_betas...))
    end
end

@inline function findcleft(r::UnitRange, step::Int)
    if length(r) >= 2*step
        minimum(r) - 1 + step * div(length(r), step * 2)
    else
        # minimum(r) - 1 + div(length(r), 2, RoundNearest) # not in Julia 1.3
        minimum(r) - 1 + round(Int, length(r)/2)
    end
end

#=
@btime Tullio.cleave(z[],4)  setup=(z=Ref((1:200, 1:500, 1:300)))
@btime Tullio.cleave(z[],4)  setup=(z=Ref((1:200, 1:50)))
@btime Tullio.cleave(z[],4)  setup=(z=Ref((5:55,)))
=#

"""
    trisect((1:10, 1:20, 5:15)) -> lo, mid, hi

Just like `cleave`, but makes 3 pieces, for 6-core machines.
"""
@inline trisect(::Tuple{}) = (), (), ()
@inline trisect(ranges::Tuple{UnitRange}) = map(tuple, findthree(first(ranges)))
@inline function trisect(ranges::Tuple{UnitRange,UnitRange})
    r1, r2 = ranges
    if length(r1) > length(r2)
        a,b,c = findthree(r1)
        return (a,r2), (b,r2), (c,r2)
    else
        a,b,c = findthree(r2)
        return (r1,a), (r1,b), (r1,c)
    end
end
@inline @generated function trisect(ranges::Tuple{Vararg{<:UnitRange,N}}) where {N}
    ex_finds = [quote
        li = length(ranges[$i])
        if li>l
            c = $i
            l = li
        end
    end for i in 1:N]
    ex_alpas = [:($i==c ? (lo) : (ranges[$i])) for i in 1:N]
    ex_betas = [:($i==c ? (mid) : (ranges[$i])) for i in 1:N]
    ex_gammas = [:($i==c ? (hi) : (ranges[$i])) for i in 1:N]
    quote
        c, l = 0, 0
        $(ex_finds...)
        lo,mid,hi = findthree(ranges[c])
        tuple($(ex_alpas...)), tuple($(ex_betas...)), tuple($(ex_gammas...))
    end
end

@inline function findthree(r::UnitRange)
    d = div(length(r), 3)
    i0 = first(r)
    (i0 : i0+d-1), (i0+d : i0+2d-1), (i0+2d : i0+length(r)-1)
end

#========== the end ==========#


text:
using Pkg
using BinaryBuilder, BinaryBuilder.BinaryBuilderBase
using BinaryBuilder.BinaryBuilderBase: preferred_runner, platform_dlext, platform_exeext
using Base.BinaryPlatforms
using Random, LibGit2, Test, ObjectFile, SHA
import Libdl

# The platform we're running on
const platform = HostPlatform()
const build_tests_dir = joinpath(@__DIR__, "build_tests")

# Helper function to run fortran code with the path to libgfortran/libquadmath
# embedded in the appropriate environment variables (JLL packages we love you so)
csl_path = dirname(first(filter(x -> occursin("libgfortran", x), Libdl.dllist())))
LIBPATH_var, envsep = if Sys.iswindows()
    ("PATH", ";")
elseif Sys.isapple()
    ("DYLD_LIBRARY_PATH", ":")
else
    ("LD_LIBRARY_PATH", ":")
end
function with_libgfortran(f::Function)
    libpath_list = [csl_path split(get(ENV, LIBPATH_var, ""), envsep)]
    libpath = join(filter(x -> !isempty(x), libpath_list), envsep)
    withenv(f, LIBPATH_var => libpath)
end

## Tests involing building packages and whatnot
libfoo_products = [
    LibraryProduct("libfoo", :libfoo),
    ExecutableProduct("fooifier", :fooifier),
]

libfoo_make_script = raw"""
cd ${WORKSPACE}/srcdir/libfoo
make install
install_license ${WORKSPACE}/srcdir/libfoo/LICENSE.md
"""

libfoo_cmake_script = raw"""
mkdir ${WORKSPACE}/srcdir/libfoo/build && cd ${WORKSPACE}/srcdir/libfoo/build
cmake -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} ..
make install
install_license ${WORKSPACE}/srcdir/libfoo/LICENSE.md
"""

libfoo_meson_script = raw"""
mkdir ${WORKSPACE}/srcdir/libfoo/build && cd ${WORKSPACE}/srcdir/libfoo/build
meson .. -Dprefix=${prefix} --cross-file="${MESON_TARGET_TOOLCHAIN}"
meson install

# grumble grumble meson!  Why do you go to all the trouble to build it properly
# in `build`, then screw it up when you `install` it?!  Silly willy.
if [[ ${target} == *apple* ]]; then
    install_name_tool ${prefix}/bin/fooifier -change ${prefix}/lib/libfoo.0.dylib @rpath/libfoo.0.dylib
fi
install_license ${WORKSPACE}/srcdir/libfoo/LICENSE.md
"""

libfoo_autotools_script = raw"""
cd ${WORKSPACE}/srcdir/libfoo
autoreconf -fiv
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --disable-static
make install
install_license ${WORKSPACE}/srcdir/libfoo/LICENSE.md
"""

# Run all our tests
include("basic.jl")
include("building.jl")
include("auditing.jl")
include("jll.jl")
include("wizard.jl")
include("declarative.jl")


text:
module Ion

using Comonicon
using PkgTemplates
using LibGit2
using Pkg

include("project.jl")
include("registry.jl")
include("utils.jl")

@command_main name = "ion" version = "0.1.0" doc = """
      The Ion manager.
      """

# include("cmd.jl")

include("precompile.jl")
_precompile_()

end # module