text:
##
## Interface functions
##

function sample_init(
    rng::Union{AbstractRNG,AbstractVector{<:AbstractRNG}},
    h::Hamiltonian,
    θ::AbstractVecOrMat{<:AbstractFloat}
)
    # Ensure h.metric has the same dim as θ.
    h = update(h, θ)
    # Initial transition
    t = Transition(phasepoint(rng, θ, h), NamedTuple())
    return h, t
end

# A step is a momentum refreshment plus a transition
function step(
    rng::Union{AbstractRNG,AbstractVector{<:AbstractRNG}},
    h::Hamiltonian,
    τ::AbstractProposal,
    z::PhasePoint
)
    # Refresh momentum
    z = refresh(rng, z, h)
    # Make transition
    return transition(rng, τ, h, z)
end

adapt!(
    h::Hamiltonian,
    τ::AbstractProposal,
    adaptor::Adaptation.NoAdaptation,
    i::Int,
    n_adapts::Int,
    θ::AbstractVecOrMat{<:AbstractFloat},
    α::AbstractScalarOrVec{<:AbstractFloat}
) = h, τ, false

function adapt!(
    h::Hamiltonian,
    τ::AbstractProposal,
    adaptor::Adaptation.AbstractAdaptor,
    i::Int,
    n_adapts::Int,
    θ::AbstractVecOrMat{<:AbstractFloat},
    α::AbstractScalarOrVec{<:AbstractFloat}
)
    isadapted = false
    if i <= n_adapts
        i == 1 && Adaptation.initialize!(adaptor, n_adapts)
        adapt!(adaptor, θ, α)
        i == n_adapts && finalize!(adaptor)
        h, τ = update(h, τ, adaptor)
        isadapted = true
    end
    return h, τ, isadapted
end

"""
Progress meter update with all trajectory stats, iteration number and metric shown.
"""
function pm_next!(pm, stat::NamedTuple, i::Int, metric::AbstractMetric)
    # Add current iteration and mass matrix
    stat = (iterations=i, stat..., mass_matrix=metric)
    ProgressMeter.next!(pm; showvalues=[tuple(s...) for s in pairs(stat)])
end

"""
Simple progress meter update without any show values.
"""
simple_pm_next!(pm, stat::NamedTuple, ::Int, ::AbstractMetric) = ProgressMeter.next!(pm)

##
## Sampling functions
##

sample(
    h::Hamiltonian,
    τ::AbstractProposal,
    θ::AbstractVecOrMat{<:AbstractFloat},
    n_samples::Int,
    adaptor::Adaptation.AbstractAdaptor=Adaptation.NoAdaptation(),
    n_adapts::Int=min(div(n_samples, 10), 1_000);
    drop_warmup=false,
    verbose::Bool=true,
    progress::Bool=false,
    (pm_next!)::Function=pm_next!
) = sample(
    GLOBAL_RNG,
    h,
    τ,
    θ,
    n_samples,
    adaptor,
    n_adapts;
    drop_warmup=drop_warmup,
    verbose=verbose,
    progress=progress,
    (pm_next!)=pm_next!,
)

"""
    sample(
        rng::AbstractRNG,
        h::Hamiltonian,
        τ::AbstractProposal,
        θ::AbstractVecOrMat{T},
        n_samples::Int,
        adaptor::Adaptation.AbstractAdaptor=Adaptation.NoAdaptation(),
        n_adapts::Int=min(div(n_samples, 10), 1_000);
        drop_warmup::Bool=false,
        verbose::Bool=true,
        progress::Bool=false
    )

Sample `n_samples` samples using the proposal `τ` under Hamiltonian `h`.
- The randomness is controlled by `rng`.
    - If `rng` is not provided, `GLOBAL_RNG` will be used.
- The initial point is given by `θ`.
- The adaptor is set by `adaptor`, for which the default is no adaptation.
    - It will perform `n_adapts` steps of adaptation, for which the default is the minimum of `1_000` and 10% of `n_samples`
- `drop_warmup` controls to drop the samples during adaptation phase or not
- `verbose` controls the verbosity
- `progress` controls whether to show the progress meter or not
"""
function sample(
    rng::Union{AbstractRNG,AbstractVector{<:AbstractRNG}},
    h::Hamiltonian,
    τ::AbstractProposal,
    θ::T,
    n_samples::Int,
    adaptor::Adaptation.AbstractAdaptor=Adaptation.NoAdaptation(),
    n_adapts::Int=min(div(n_samples, 10), 1_000);
    drop_warmup=false,
    verbose::Bool=true,
    progress::Bool=false,
    (pm_next!)::Function=pm_next!
) where {T<:AbstractVecOrMat{<:AbstractFloat}}
    @assert !(drop_warmup && (adaptor isa Adaptation.NoAdaptation)) "Cannot drop warmup samples if there is no adaptation phase."
    # Prepare containers to store sampling results
    n_keep = n_samples - (drop_warmup ? n_adapts : 0)
    θs, stats = Vector{T}(undef, n_keep), Vector{NamedTuple}(undef, n_keep)
    # Initial sampling
    h, t = sample_init(rng, h, θ)
    # Progress meter
    pm = progress ? ProgressMeter.Progress(n_samples, desc="Sampling", barlen=31) : nothing
    time = @elapsed for i = 1:n_samples
        # Make a step
        t = step(rng, h, τ, t.z)
        # Adapt h and τ; what mutable is the adaptor
        tstat = stat(t)
        h, τ, isadapted = adapt!(h, τ, adaptor, i, n_adapts, t.z.θ, tstat.acceptance_rate)
        tstat = merge(tstat, (is_adapt=isadapted,))
        # Update progress meter
        if progress
            pm_next!(pm, tstat, i, h.metric)
        # Report finish of adapation
        elseif verbose && isadapted && i == n_adapts
            @info "Finished $n_adapts adapation steps" adaptor τ.integrator h.metric
        end
        # Store sample
        if !drop_warmup || i > n_adapts
            j = i - drop_warmup * n_adapts
            θs[j], stats[j] = t.z.θ, tstat
        end
    end
    # Report end of sampling
    if verbose
        EBFMI_est = EBFMI(map(s -> s.hamiltonian_energy, stats))
        average_acceptance_rate = mean(map(s -> s.acceptance_rate, stats))
        if θ isa AbstractVector
            n_chains = 1
        else
            n_chains = size(θ, 2)
            # TODO: see if there is other trick to make
            EBFMI_est = "[" * join(EBFMI_est, ", ") * "]"
            average_acceptance_rate = "[" * join(average_acceptance_rate, ", ") * "]"
        end
        @info "Finished $n_samples sampling steps for $n_chains chains in $time (s)" h τ EBFMI_est average_acceptance_rate
    end
    return θs, stats
end


text:
# TODO: define VertexNamedDimDataGraph, a graph with only data on the
# vertices, and EdgeNamedDimDataGraph, a graph with only data on the edges.
struct NamedDimDataGraph{VD,ED,V,E,G<:AbstractGraph} <: AbstractDataGraph{VD,ED,V,E}
  underlying_graph::G
  vertex_data::MultiDimDictionary{V,VD}
  edge_data::Dictionary{E,ED}
end
underlying_graph(graph::NamedDimDataGraph) = graph.underlying_graph
vertex_data(graph::NamedDimDataGraph) = graph.vertex_data
edge_data(graph::NamedDimDataGraph) = graph.edge_data

function copy(graph::NamedDimDataGraph)
  return NamedDimDataGraph(
    copy(underlying_graph(graph)), copy(vertex_data(graph)), copy(edge_data(graph))
  )
end

function NamedDimDataGraph{VD,ED}(
  underlying_graph::AbstractGraph,
  vertex_data::MultiDimDictionary{V,VD},
  edge_data::Dictionary{E,ED},
) where {VD,ED,V,E}
  G = typeof(underlying_graph)
  return NamedDimDataGraph{VD,ED,V,E,G}(underlying_graph, vertex_data, edge_data)
end

function NamedDimDataGraph{VD,ED}(underlying_graph::AbstractGraph) where {VD,ED}
  V = eltype(underlying_graph)
  E = edgetype(underlying_graph)
  vertex_data = MultiDimDictionary{V,VD}()
  edge_data = Dictionary{E,ED}()
  return NamedDimDataGraph{VD,ED}(underlying_graph, vertex_data, edge_data)
end

function NamedDimDataGraph(underlying_graph::AbstractGraph)
  return NamedDimDataGraph{Any,Any}(underlying_graph)
end
function NamedDimDataGraph{VD}(underlying_graph::AbstractGraph) where {VD}
  return NamedDimDataGraph{VD,Any}(underlying_graph)
end

## function default_data(
##   index_type::Type,
##   data_type::Type,
##   indices_function::Function,
##   underlying_graph::AbstractGraph,
##   ::Nothing
## )
##   return similar(Indices(indices_function(underlying_graph)), data_type)
## end
##
## # XXX: Assumes `is_directed(underlying_graph)`.
## function default_data(
##   index_type::Type{<:AbstractEdge},
##   data_type::Type,
##   indices_function::Function,
##   underlying_graph::AbstractGraph,
##   ::Nothing
## )
##   out_indices = indices_function(underlying_graph)
##   in_indices = reverse.(out_indices)
##   # Interleave the indices.
##   indices = collect(Iterators.flatten(zip(out_indices, in_indices)))
##   return similar(Indices(indices), data_type)
## end

## function NamedDimDataGraph{VD,ED}(
##   underlying_graph::G,
##   vertex_data,
##   edge_data
## ) where {VD,ED,G<:AbstractGraph}
##   V = eltype(underlying_graph)
##   E = edgetype(underlying_graph)
##   vertex_data = MultiDimDictionary{V,VD}(vertices(underlying_graph), vertex_data)
##   edge_data = MultiDimDictionary{E,ED}(edges(underlying_graph), edge_data)
##
##   @show vertex_data
##   @show edge_data
##
##   return NamedDimDataGraph{VD,ED,V,E,G}(underlying_graph, vertex_data, edge_data)
## end
##
## copy(graph::NamedDimDataGraph) = NamedDimDataGraph(copy(underlying_graph(graph)), copy(vertex_data(graph)), copy(edge_data(graph)))
##
## function NamedDimDataGraph{VD}(
##   underlying_graph::AbstractGraph,
##   vertex_data,
##   edge_data
## ) where {VD}
##   ED = data_type(edge_data)
##   return NamedDimDataGraph{VD,ED}(underlying_graph, vertex_data, edge_data)
## end
##
## function (NamedDimDataGraph{VD,ED} where {VD})(
##   underlying_graph::AbstractGraph,
##   vertex_data,
##   edge_data
## ) where {ED}
##   VD = data_type(vertex_data)
##   return NamedDimDataGraph{VD,ED}(underlying_graph, vertex_data, edge_data)
## end
##
## function NamedDimDataGraph(
##   underlying_graph::AbstractGraph,
##   vertex_data,
##   edge_data
## )
##   VD = data_type(vertex_data)
##   ED = data_type(edge_data)
##   return NamedDimDataGraph{VD,ED}(underlying_graph, vertex_data, edge_data)
## end
##
## #
# kwarg versions call arg versions
#

## function NamedDimDataGraph{VD,ED}(
##   underlying_graph::AbstractGraph;
##   vertex_data=nothing,
##   edge_data=nothing
## ) where {VD,ED}
##   return NamedDimDataGraph{VD,ED}(underlying_graph, vertex_data, edge_data)
## end
##
## function NamedDimDataGraph{VD}(
##   underlying_graph::AbstractGraph;
##   vertex_data=nothing,
##   edge_data=nothing
## ) where {VD}
##   return NamedDimDataGraph{VD}(underlying_graph, vertex_data, edge_data)
## end
##
## function (NamedDimDataGraph{VD,ED} where {VD})(
##   underlying_graph::AbstractGraph;
##   vertex_data=nothing,
##   edge_data=nothing
## ) where {ED}
##   return (NamedDimDataGraph{VD,ED} where {VD})(underlying_graph, vertex_data, edge_data)
## end
##
## function NamedDimDataGraph(
##   underlying_graph::AbstractGraph;
##   vertex_data=nothing,
##   edge_data=nothing
## )
##   return NamedDimDataGraph(underlying_graph, vertex_data, edge_data)
## end

ext: jl


text:
using ISAtmosphere
using Test

@testset "ISAtmosphere.jl" begin
    @test T_K(9000.0, -1.5) ≈ 228.15 atol = 0.01
    @test T_K(13000.0, 0.0) ≈ 216.65 atol = 0.01
    @test T_K(13000) ≈ 216.65 atol = 0.01
    @test p_Pa(9000.0, -1.5) ≈ 30742.5 atol = 0.1
    @test p_Pa(13000.0, 0.0) ≈ 16510.4 atol = 0.1
    @test p_Pa(13000.0, 0) ≈ 16510.4 atol = 0.1
    @test ρ_kg_m³(30742.5, 228.15) ≈ 0.469414 atol = 0.0001
    @test ρ_kg_m³(16510.4, 216.65) ≈ 0.265483 atol = 0.0001
    @test ρ_kg_m³(16510, 217) ≈ 0.265049 atol = 0.0001
    @test a_m_s(228.15) ≈ 302.8 atol = 0.1
    @test a_m_s(216.65) ≈ 295.07 atol = 0.1
    @test a_m_s(217) ≈ 295.308 atol = 0.01

    @test Vcas2Vtas(128.611, p_Pa(9000.0), T_K(9000.0)) ≈ 201.01 atol = 0.1
    @test Vtas2Vcas(201.01, p_Pa(9000.0), T_K(9000.0)) ≈ 128.611 atol = 0.1
    @test Vtas2Vcas(201, p_Pa(9000.0), T_K(9000.0)) ≈ 128.611 atol = 0.1
    @test Vtas2M(201.01, T_K(9000.0)) ≈ 0.661667 atol = 0.001
    @test Vtas2M(201, T_K(9000.0)) ≈ 0.661667 atol = 0.001
    @test M2Vtas(0.661667, T_K(9000.0)) ≈ 201.01 atol = 0.1
    @test Vcas2M(128.611, p_Pa(9000.0), T_K(9000.0)) ≈ 0.661667 atol = 0.001
    @test M2Vcas(0.661667, p_Pa(9000.0), T_K(9000.0)) ≈ 128.611 atol = 0.1

    @test Hp_trans_m(133.756, 0.8) ≈ 11260.287 atol = 0.1
    @test Hp_trans_m(133.756, 0.6) ≈ 6952.80 atol = 0.1

    @test ISAtmosphere.θ(288) ≈ 0.999479 atol = 0.001
    @test ISAtmosphere.δ(101325) ≈ 1.000 atol = 0.001
    @test ISAtmosphere.σ(1) ≈ 0.816327 atol = 0.001

    cond = ISAtmosphere.conditions(9000.0, -1.5)
    @test cond.Hp_m ≈ 9000.0 atol = 0.1
    @test cond.T_K ≈ 228.15 atol = 0.01
    @test cond.p_Pa ≈ 30742.5 atol = 0.1
    @test cond.ρ_kg_m³ ≈ 0.469414 atol = 0.0001
    @test cond.a_m_s ≈ 302.8 atol = 0.1

    cond2 = conditions(9000.0, 0.0)
    @test cond2.Hp_m ≈ 9000.0 atol = 0.1
    @test cond2.T_K ≈ 229.65 atol = 0.01
    @test cond2.p_Pa ≈ 30742.5 atol = 0.1
    @test cond2.ρ_kg_m³ ≈ 0.466348 atol = 0.0001
    @test cond2.a_m_s ≈ 303.793 atol = 0.1

    @test Vcas2Vtas(128.611, cond2) ≈ 201.01 atol = 0.1
    @test Vtas2Vcas(201.01, cond2) ≈ 128.611 atol = 0.1
    @test Vtas2Vcas(201, cond2) ≈ 128.611 atol = 0.1
    @test Vtas2M(201.01, cond2) ≈ 0.661667 atol = 0.001
    @test Vtas2M(201, cond2) ≈ 0.661667 atol = 0.001
    @test M2Vtas(0.661667, cond2) ≈ 201.01 atol = 0.1
    @test Vcas2M(128.611, cond2) ≈ 0.661667 atol = 0.001
    @test M2Vcas(0.661667, cond2) ≈ 128.611 atol = 0.1
end


text:
function reconstr_x(u, x, xmin, xmax)

    ntau  = size(u)[1]
    nx    = size(x)[1]
    L     = xmax - xmin
    UA    = zeros(ComplexF64, (ntau,nx))
    v     = fft(u) / nx

    for j in 1:nx÷2
        UA .+= v[:, j] .* exp.(1im * 2π / L * (j-1) * (x' .- xmin))
    end

    for j in nx÷2+1:nx
        UA .+= v[:, j] .* exp.(1im * 2π / L * (j-1-nx) * (x' .- xmin))
    end

    UA

end


text:
## reductions ##

###### higher level reduction functions ######

# Note that getting type-stable results from reduction functions,
# or at least having type-stable loops, is nontrivial (#6069).

# reduce

function reduce(op::Callable, itr) # this is a left fold
    if is(op,+)
        return sum(itr)
    elseif is(op,*)
        return prod(itr)
    elseif is(op,|)
        return any(itr)
    elseif is(op,&)
        return all(itr)
    end
    s = start(itr)
    if done(itr, s)
        return op()  # empty collection
    end
    (v, s) = next(itr, s)
    if done(itr, s)
        return v
    else # specialize for length > 1 to have a hopefully type-stable loop
        (x, s) = next(itr, s)
        result = op(v, x)
        while !done(itr, s)
            (x, s) = next(itr, s)
            result = op(result, x)
        end
        return result
    end
end

# pairwise reduction, requires n > 1 (to allow type-stable loop)
function r_pairwise(op::Callable, A::AbstractArray, i1,n)
    if n < 128
        @inbounds v = op(A[i1], A[i1+1])
        for i = i1+2:i1+n-1
            @inbounds v = op(v,A[i])
        end
        return v
    else
        n2 = div(n,2)
        return op(r_pairwise(op,A, i1,n2), r_pairwise(op,A, i1+n2,n-n2))
    end
end

function reduce(op::Callable, A::AbstractArray)
    n = length(A)
    n == 0 ? op() : n == 1 ? A[1] : r_pairwise(op,A, 1,n)
end

function reduce(op::Callable, v0, A::AbstractArray)
    n = length(A)
    n == 0 ? v0 : n == 1 ? op(v0, A[1]) : op(v0, r_pairwise(op,A, 1,n))
end

function reduce(op::Callable, v0, itr)
    v = v0
    if is(op,+)
        for x in itr
            v = v+x
        end
    elseif is(op,*)
        for x in itr
            v = v*x
        end
    else
        u = v0
        for x in itr
            u = op(u,x)
        end
        return u
    end
    return v
end

# foldl & foldr

function foldl(op::Callable, v0, itr, i=start(itr))
    v = v0
    while !done(itr,i)
        x, i = next(itr,i)
        v = op(v,x)
    end
    return v
end
function foldl(op::Callable, itr)
    v0, i = next(itr,start(itr))
    foldl(op, v0, itr, i)
end

function foldr(op::Callable, v0, itr, i=endof(itr))
    v = v0
    while i > 0
        x = itr[i]
        v = op(x,v)
        i = i - 1
    end
    return v
end
function foldr(op::Callable, itr)
    i = endof(itr)
    foldr(op, itr[i], itr, i-1)
end

###### Specific reduction functions ######

## in & contains

function in(x, itr)
    for y in itr
        if y==x
            return true
        end
    end
    return false
end

function contains(itr, x)
    depwarn("contains(collection, item) is deprecated, use in(item, collection) instead", :contains)
    in(x, itr)
end

function contains(eq::Function, itr, x)
    for y in itr
        if eq(y,x)
            return true
        end
    end
    return false
end

## countnz & count

function countnz(itr)
    n = 0
    for x in itr
        n += (x != 0)
    end
    return n
end

function countnz{T}(a::AbstractArray{T})
    n = 0
    z = zero(T)
    for i = 1:length(a)
        @inbounds n += (a[i] != z)
    end
    return n
end

function countnz(a::AbstractArray{Bool})
    n = 0
    for x in a
        if x; n += 1; end
    end
    return n
end

function count(pred::Function, itr)
    s = 0
    for x in itr
        if pred(x)
            s+=1
        end
    end
    s
end

## sum

function sum(itr)
    s = start(itr)
    if done(itr, s)
        if method_exists(eltype, (typeof(itr),))
            T = eltype(itr)
            return zero(T) + zero(T)
        else
            throw(ArgumentError("sum(itr) is undefined for empty collections; instead, do isempty(itr) ? z : sum(itr), where z is the correct type of zero for your sum"))
        end
    end
    (v, s) = next(itr, s)
    done(itr, s) && return v + zero(v) # adding zero for type stability
    # specialize for length > 1 to have type-stable loop
    (x, s) = next(itr, s)
    result = v + x
    while !done(itr, s)
        (x, s) = next(itr, s)
        result += x
    end
    return result
end

sum(A::AbstractArray{Bool}) = countnz(A)

# a fast implementation of sum in sequential order (from left to right).
# to allow type-stable loops, requires length > 1
function sum_seq{T}(a::AbstractArray{T}, ifirst::Int, ilast::Int)

    @inbounds if ifirst + 6 >= ilast  # length(a) < 8
        i = ifirst
        s = a[i] + a[i+1]
        i = i+1
        while i < ilast
            s += a[i+=1]
        end
        return s

    else # length(a) >= 8

        # more effective utilization of the instruction
        # pipeline through manually unrolling the sum
        # into four-way accumulation. Benchmark shows
        # that this results in about 2x speed-up.

        s1 = a[ifirst] + a[ifirst + 4]
        s2 = a[ifirst + 1] + a[ifirst + 5]
        s3 = a[ifirst + 2] + a[ifirst + 6]
        s4 = a[ifirst + 3] + a[ifirst + 7]

        i = ifirst + 8
        il = ilast - 3
        while i <= il
            s1 += a[i]
            s2 += a[i+1]
            s3 += a[i+2]
            s4 += a[i+3]
            i += 4
        end

        while i <= ilast
            s1 += a[i]
            i += 1
        end

        return s1 + s2 + s3 + s4
    end
end

# Pairwise (cascade) summation of A[i1:i1+n-1], which has O(log n) error growth
# [vs O(n) for a simple loop] with negligible performance cost if
# the base case is large enough.  See, e.g.:
#        http://en.wikipedia.org/wiki/Pairwise_summation
#        Higham, Nicholas J. (1993), "The accuracy of floating point
#        summation", SIAM Journal on Scientific Computing 14 (4): 783–799.
# In fact, the root-mean-square error growth, assuming random roundoff
# errors, is only O(sqrt(log n)), which is nearly indistinguishable from O(1)
# in practice.  See:
#        Manfred Tasche and Hansmartin Zeuner, Handbook of
#        Analytic-Computational Methods in Applied Mathematics (2000).
#

# Note: sum_seq uses four accumulators, so each accumulator gets at most 256 numbers
const PAIRWISE_SUM_BLOCKSIZE = 1024

# note: requires length > 1, due to sum_seq
function sum_pairwise(a::AbstractArray, ifirst::Int, ilast::Int)
    # bsiz: maximum block size

    if ifirst + PAIRWISE_SUM_BLOCKSIZE >= ilast
        sum_seq(a, ifirst, ilast)
    else
        imid = (ifirst + ilast) >>> 1
        sum_pairwise(a, ifirst, imid) + sum_pairwise(a, imid+1, ilast)
    end
end

function sum{T}(a::AbstractArray{T})
    n = length(a)
    n == 0 && return zero(T) + zero(T)
    n == 1 && return a[1] + zero(a[1])
    sum_pairwise(a, 1, length(a))
end

function sum{T<:Integer}(a::AbstractArray{T})
    n = length(a)
    n == 0 && return zero(T) + zero(T)
    n == 1 && return a[1] + zero(a[1])
    sum_seq(a, 1, length(a))
end

# Kahan (compensated) summation: O(1) error growth, at the expense
# of a considerable increase in computational expense.
function sum_kbn{T<:FloatingPoint}(A::AbstractArray{T})
    n = length(A)
    if (n == 0)
        return zero(T)+zero(T)
    end
    s = A[1]+zero(T)
    c = zero(T)+zero(T)
    for i in 2:n
        Ai = A[i]
        t = s + Ai
        if abs(s) >= abs(Ai)
            c += ((s-t) + Ai)
        else
            c += ((Ai-t) + s)
        end
        s = t
    end

    s + c
end

## prod

function prod(itr)
    s = start(itr)
    if done(itr, s)
        if method_exists(eltype, (typeof(itr),))
            T = eltype(itr)
            return one(T) * one(T)
        else
            throw(ArgumentError("prod(itr) is undefined for empty collections; instead, do isempty(itr) ? o : prod(itr), where o is the correct type of identity for your product"))
        end
    end
    (v, s) = next(itr, s)
    done(itr, s) && return v * one(v) # multiplying by one for type stability
    # specialize for length > 1 to have type-stable loop
    (x, s) = next(itr, s)
    result = v * x
    while !done(itr, s)
        (x, s) = next(itr, s)
        result *= x
    end
    return result
end

prod(A::AbstractArray{Bool}) =
    error("use all() instead of prod() for boolean arrays")

function prod_rgn{T}(A::AbstractArray{T}, first::Int, last::Int)
    if first > last
        return one(T) * one(T)
    end
    i = first
    @inbounds v = A[i]
    i == last && return v * one(v)
    @inbounds result = v * A[i+=1]
    while i < last
        @inbounds result *= A[i+=1]
    end
    return result
end
prod{T}(A::AbstractArray{T}) = prod_rgn(A, 1, length(A))

## maximum & minimum

function maximum(itr)
    s = start(itr)
    if done(itr, s)
        error("argument is empty")
    end
    (v, s) = next(itr, s)
    while !done(itr, s)
        (x, s) = next(itr, s)
        v = scalarmax(v,x)
    end
    return v
end

function minimum(itr)
    s = start(itr)
    if done(itr, s)
        error("argument is empty")
    end
    (v, s) = next(itr, s)
    while !done(itr, s)
        (x, s) = next(itr, s)
        v = scalarmin(v,x)
    end
    return v
end

function maximum_rgn{T<:Real}(A::AbstractArray{T}, first::Int, last::Int)
    if first > last; error("argument range must not be empty"); end

    # locate the first non NaN number
    v = A[first]
    i = first + 1
    while v != v && i <= last
        @inbounds v = A[i]
        i += 1
    end

    while i <= last
        @inbounds x = A[i]
        if x > v
            v = x
        end
        i += 1
    end
    v
end

function minimum_rgn{T<:Real}(A::AbstractArray{T}, first::Int, last::Int)
    if first > last; error("argument range must not be empty"); end

    # locate the first non NaN number
    v = A[first]
    i = first + 1
    while v != v && i <= last
        @inbounds v = A[i]
        i += 1
    end

    while i <= last
        @inbounds x = A[i]
        if x < v
            v = x
        end
        i += 1
    end
    v
end

maximum{T<:Real}(A::AbstractArray{T}) = maximum_rgn(A, 1, length(A))
minimum{T<:Real}(A::AbstractArray{T}) = minimum_rgn(A, 1, length(A))

## extrema

function extrema(itr)
    s = start(itr)
    if done(itr, s)
        error("argument is empty")
    end
    (v, s) = next(itr, s)
    vmin = v
    vmax = v
    while !done(itr, s)
        (x, s) = next(itr, s)
        if x == x
            if x > vmax
                vmax = x
            elseif x < vmin
                vmin = x
            end
        end
    end
    return (vmin, vmax)
end

function extrema{T<:Real}(A::AbstractArray{T})
    if isempty(A); error("argument must not be empty"); end
    n = length(A)

    # locate the first non NaN number
    v = A[1]
    i = 2
    while v != v && i <= n
        @inbounds v = A[i]
        i += 1
    end

    vmin = v
    vmax = v
    while i <= n
        @inbounds v = A[i]
        if v > vmax
            vmax = v
        elseif v < vmin
            vmin = v
        end
        i += 1
    end

    return (vmin, vmax)
end

## all & any

function all(itr)
    for x in itr
        if !x
            return false
        end
    end
    return true
end

function any(itr)
    for x in itr
        if x
            return true
        end
    end
    return false
end

###### mapreduce ######

function mapreduce(f::Callable, op::Callable, itr)
    s = start(itr)
    if done(itr, s)
        return op()  # empty collection
    end
    (x, s) = next(itr, s)
    v = f(x)
    if done(itr, s)
        return v
    else # specialize for length > 1 to have a hopefully type-stable loop
        (x, s) = next(itr, s)
        result = op(v, f(x))
        while !done(itr, s)
            (x, s) = next(itr, s)
            result = op(result, f(x))
        end
        return result
    end
end

function mapreduce(f::Callable, op::Callable, v0, itr)
    v = v0
    for x in itr
        v = op(v,f(x))
    end
    return v
end

# pairwise reduction, requires n > 1 (to allow type-stable loop)
function mr_pairwise(f::Callable, op::Callable, A::AbstractArray, i1,n)
    if n < 128
        @inbounds v = op(f(A[i1]), f(A[i1+1]))
        for i = i1+2:i1+n-1
            @inbounds v = op(v,f(A[i]))
        end
        return v
    else
        n2 = div(n,2)
        return op(mr_pairwise(f,op,A, i1,n2), mr_pairwise(f,op,A, i1+n2,n-n2))
    end
end
function mapreduce(f::Callable, op::Callable, A::AbstractArray)
    n = length(A)
    n == 0 ? op() : n == 1 ? f(A[1]) : mr_pairwise(f,op,A, 1,n)
end
function mapreduce(f::Callable, op::Callable, v0, A::AbstractArray)
    n = length(A)
    n == 0 ? v0 : n == 1 ? op(v0, f(A[1])) : op(v0, mr_pairwise(f,op,A, 1,n))
end

# specific mapreduce functions

maximum(f::Function, itr) = mapreduce(f, scalarmax, itr)
minimum(f::Function, itr) = mapreduce(f, scalarmin, itr)
sum(f::Function, itr)     = mapreduce(f, +        , itr)
prod(f::Function, itr)    = mapreduce(f, *        , itr)

function any(pred::Function, itr)
    for x in itr
        if pred(x)
            return true
        end
    end
    return false
end

function all(pred::Function, itr)
    for x in itr
        if !pred(x)
            return false
        end
    end
    return true
end

max_stars_repo_name: alyst/julia

max_issues_repo_name: alyst/julia

max_forks_repo_name: alyst/julia


text:
function rotate(pos, instr)
    alefbet = "abcdefghijklmnopqrstuvwxyz"
    res = ""
    for c in instr
        isup = isuppercase(c)
        orgpos = findfirst(lowercase(c), alefbet)
        if orgpos == nothing
            res = res * c
        else
            newpos = (orgpos + pos) % 26
            if newpos == 0
                newpos = 26
            end
            newletter = !isup ? alefbet[newpos] : uppercase(alefbet[newpos])
            res = res * newletter
        end
    end
    if instr isa String
        return res
    else
        return res[1]
    end
end

max_stars_repo_name: stepinski/julia-in

max_issues_repo_name: stepinski/julia-init

max_forks_repo_name: stepinski/julia-init


text:
const M = 10
function generate_f(N, d, k; X=rand(N, d))
    f = rand(MvNormal(zeros(N), kernelmatrix(k, X; obsdim=1) + 1e-5I))
    return X, f
end
const VIs = ["AnalyticVI", "QuadratureVI", "MCIntegrationVI"]
const VIcodes = ["AVI", "QVI", "MCVI"]

function tests(model1, model2, X, f, y, problem)
    model1, state1 = train!(model1, X, y, 1)
    L = AGP.objective(model1, X, y)
    train!(model1, X, y, 5)
    @test testconv(model1, problem, X, f, y)
    @test all(proba_y(model1, X)[2] .> 0)
    train!(model2, X, y, 6)
    @test testconv(model2, problem, X, f, y)
    @test all(proba_y(model2, X)[2] .> 0)
end

function tests(model1::OnlineSVGP, model2, X, f, y, problem)
    state1 = nothing
    for (X_, y_) in eachbatch((X, y); obsdim=1, size=10)
        model1, state1 = train!(model1, X_, y_, state1; iterations=5)
        L = AGP.ELBO(model1, state1, X_, y_)
    end
    @test testconv(model1, problem, X, f, y)
    @test all(proba_y(model1, X)[2] .> 0)
    state2 = nothing
    for (X_, y_) in eachbatch((X, y); obsdim=1, size=10)
        model2, state2 = train!(model2, X_, y_, state2; iterations=5)
    end
    @test testconv(model2, problem, X, f, y)
    @test all(proba_y(model2, X)[2] .> 0)
end

function tests(
    model1::AbstractGPModel{T,<:BernoulliLikelihood{<:AGP.SVMLink}},
    model2,
    X,
    f,
    y,
    problem,
) where {T}
    model1, state1 = train!(model1, X, y, 1)
    @test_nowarn AGP.objective(model1, X, y)
    train!(model1, X, y, 5)
    @test testconv(model1, problem, X, f, y)
    @test all(proba_y(model1, X)[2] .> 0)
    model2, state2 = train!(model2, X, y, 6)
    @test testconv(model2, problem, X, f, y)
    @test all(proba_y(model2, X)[2] .> 0)
end

function tests_likelihood(
    l, # likelihood
    ltype, # likelihood types
    dict::Dict,
    floattypes,
    problem,
    n_latent,
    X,
    f,
    y,
    k,
)
    Z = inducingpoints(KmeansAlg(M), X)
    @testset "VGP" begin
        for floattype in floattypes
            dictvgp = dict["VGP"]
            for (name, code, inference) in
                zip(VIs, VIcodes, [AnalyticVI(), QuadratureVI(), MCIntegrationVI()])
                @testset "$name" begin
                    test_inference_VGP(
                        X,
                        y,
                        f,
                        k,
                        l,
                        ltype,
                        floattype,
                        n_latent,
                        problem,
                        inference;
                        valid=dictvgp[code],
                    )
                end
            end
        end # Loop on float types
    end # VGP
    @testset "OSVGP" begin
        dictosvgp = dict["OSVGP"]
        for floattype in floattypes
            @testset "AnalyticVI" begin
                if dictosvgp["AVI"]
                    model = OnlineSVGP(
                        k,
                        l,
                        AnalyticVI(),
                        AGP.InducingPoints.OIPS();
                        optimiser=false,
                        verbose=0,
                    )
                    @test eltype(model) == floattype
                    @test AGP.likelihood(model) isa ltype
                    @test AGP.inference(model) isa AnalyticVI
                    @test AGP.getf(model) isa NTuple{n_latent,AGP.OnlineVarLatent}
                    @test AGP.n_latent(model) == n_latent
                    model_opt = OnlineSVGP(
                        k,
                        l,
                        AnalyticVI(),
                        AGP.InducingPoints.OIPS();
                        optimiser=true,
                        verbose=0,
                    )
                    tests(model, model_opt, X, f, y, problem)
                else
                    @test_throws ErrorException OnlineSVGP(
                        k, l, AnalyticVI(), UniGrid(10); optimiser=false, verbose=0
                    )
                end
            end  # Analytic VI
            @testset "NumericalVI" begin
                if dictosvgp["QVI"]
                    model = OnlineSVGP(k, l, QuadratureVI(), UniGrid(10); optimiser=false)
                    @test eltype(model) == floattype
                    @test AGP.likelihood(model) isa ltype
                    @test AGP.inference(model) isa QuadratureVI
                    @test AGP.getf(model) isa NTuple{n_latent,AGP.OnlineVarLatent}
                    @test AGP.n_latent(model) == n_latent
                    model_opt = OnlineSVGP(
                        k, l, QuadratureVI(), UniGrid(10); optimiser=true, verbose=0
                    )
                    tests(model, model_opt, X, f, y, problem)
                else
                    @test_throws ErrorException OnlineSVGP(
                        k, l, QuadratureVI(), UniGrid(10); optimiser=false, verbose=0
                    )
                end
                if dictosvgp["MCVI"]
                    model = OnlineSVGP(
                        k, l, MCIntegrationVI(), UniGrid(10); optimiser=false
                    )
                    @test eltype(model) == floattype
                    @test AGP.likelihood(model) isa ltype
                    @test AGP.inference(model) isa MCIntegrationVI
                    @test AGP.getf(model) isa NTuple{n_latent,AGP.SparseVarLatent}
                    @test AGP.n_latent(model) == n_latent
                    model_opt = OnlineSVGP(
                        k, l, MCIntegrationVI(), UniGrid(10); optimiser=true, verbose=0
                    )
                    tests(model, model_opt, X, f, y, problem)
                else
                    @test_throws ErrorException OnlineSVGP(
                        k, l, MCIntegrationVI(), UniGrid(10); optimiser=false, verbose=0
                    )
                end
            end # Loop on Numerical VI
        end # Loop on float types
    end #Loop on Online SVGP
    @testset "SVGP" begin
        dictsvgp = dict["SVGP"]
        for floattype in floattypes
            for (name, code, inference, stoch_inference) in zip(
                VIs,
                VIcodes,
                [AnalyticVI(), QuadratureVI(), MCIntegrationVI()],
                [AnalyticSVI(10), QuadratureSVI(10), MCIntegrationSVI(10)],
            )
                @testset "$name" begin
                    test_inference_SVGP(
                        X,
                        y,
                        f,
                        k,
                        l,
                        ltype,
                        floattype,
                        Z,
                        n_latent,
                        problem,
                        inference,
                        stoch_inference;
                        valid=dictsvgp[code],
                    )
                end #
            end
        end # Loop on float types
    end # SVGP
    @testset "MCGP" begin
        dictmcgp = dict["MCGP"]
        for floattype in floattypes
            @testset "Gibbs Sampling" begin
                if dictmcgp["Gibbs"]
                    model = MCGP(X, y, k, l, GibbsSampling())
                    @test eltype(model) == floattype
                    @test AGP.likelihood(model) isa ltype
                    @test AGP.inference(model) isa GibbsSampling
                    @test AGP.getf(model) isa NTuple{n_latent,AGP.SampledLatent}
                    @test AGP.n_latent(model) == n_latent
                    samples = sample(model, 100; progress=false)
                    @test_broken samples2 = sample(model, 100; cat=true, progress=false)
                else
                    @test_throws ErrorException MCGP(X, y, k, l, GibbsSampling())
                end
                if dictmcgp["HMC"]
                    model = MCGP(X, y, k, l, HMCSampling())
                    @test eltype(model) == floattype
                    @test AGP.likelihood(model) isa ltype
                    @test AGP.inference(model) isa HMCSampling
                    @test AGP.getf(model) isa NTuple{n_latent,AGP.SampledLatent}
                    @test AGP.n_latent(model) == n_latent
                    samples = sample(model, 20; progress=false)
                    @test_broken samples2 = sample(model, 20; cat=true)
                else
                    @test_throws ErrorException MCGP(X, y, k, l, HMCSampling())
                end
            end # Gibbs Sampling
        end # Loop floattypes
    end # MCGP
end

function testconv(
    model::AbstractGPModel,
    problem_type::String,
    X::AbstractArray,
    f::AbstractArray,
    y::AbstractArray,
)
    μ, Σ = predict_f(model, X; cov=true, diag=false)
    μ, diagΣ = predict_f(model, X; cov=true, diag=true)
    y_pred = predict_y(model, X)
    py_pred = proba_y(model, X)
    if problem_type == "Regression"
        err = mean(abs.(y_pred - f))
        # @info "Regression Error" err
        return err < 1.5
    elseif problem_type == "Classification"
        err = mean(y_pred .!= y)
        # @info "Classification Error" err
        return err < 0.5
    elseif problem_type == "MultiClass"
        err = mean(y_pred .!= y)
        # @info "Multiclass Error" err
        return err < 0.9
    elseif problem_type == "Poisson" || problem_type == "NegBinomial"
        err = mean(abs.(y_pred - y))
        # @info "Event Error" err
        return err < 20.0
    else
        throw(ErrorException("Wrong problem description $problem_type"))
    end
end

function test_inference_VGP(
    X, y, f, k, l, ltype, floattype, n_latent, problem, inference; valid=false
)
    if valid
        model = VGP(X, y, k, l, inference; optimiser=false, verbose=0)
        @test eltype(model) == floattype
        @test AGP.likelihood(model) isa ltype
        @test AGP.inference(model) isa typeof(inference)
        @test AGP.getf(model) isa NTuple{n_latent,AGP.VarLatent}
        @test AGP.n_latent(model) == n_latent
        model_opt = VGP(X, y, k, l, inference; optimiser=true, verbose=0)
        tests(model, model_opt, X, f, y, problem)
    else
        @test_throws ErrorException VGP(X, y, k, l, inference; optimiser=false, verbose=0)
    end
end

function test_inference_SVGP(
    X,
    y,
    f,
    k,
    l,
    ltype,
    floattype,
    Z,
    n_latent,
    problem,
    inference,
    stoch_inference;
    valid=false,
)
    if valid
        model = SVGP(k, l, inference, Z; optimiser=false, verbose=0)
        @test eltype(model) == floattype
        @test AGP.likelihood(model) isa ltype
        @test AGP.inference(model) isa typeof(inference)
        @test AGP.n_latent(model) == n_latent
        model_opt = SVGP(k, l, inference, Z; optimiser=true, Zoptimiser=true, verbose=0)
        tests(model, model_opt, X, f, y, problem)
        @test AGP.getf(model) isa NTuple{n_latent,AGP.SparseVarLatent}

        model_svi = SVGP(k, l, stoch_inference, Z; optimiser=false, verbose=0)
        tests(model_svi, model, X, f, y, problem)
    else
        @test_throws ErrorException SVGP(k, l, inference, Z; optimiser=false, verbose=0)
    end
end


text:
include("run2.jl")
import Gadfly
Gadfly.set_default_plot_size(20Gadfly.cm, 12Gadfly.cm)

intrvl = Any[
    (1,1,600),
    (153,0.01,154),
    (272,0.01,273),
    (368,0.01,369)
]

# generate data for the plots on given intervals and save it to disk
function data2(intrvl)

    @sync @parallel for i in 1:length(intrvl)
        eigval, Δe = q2(32, [intrvl[i][1]:intrvl[i][2]:intrvl[i][3]]);
        writedlm("p2eigvals-$(intrvl[i][1])-to-$(intrvl[i][3]).dat",eigval)
    end

end

# create plots and save them as images
function plot2(intrvl)

    plots = Array(Gadfly.Plot,0)

    @sync @parallel for i in 1:length(intrvl)
        eigval = readdlm("p2eigvals-$(intrvl[i][1])-to-$(intrvl[i][3]).dat")
        p = Gadfly.plot(
            x=[intrvl[i][1]:intrvl[i][2]:intrvl[i][3]],
            y=eigval,
            Gadfly.Geom.line, Gadfly.Geom.point,
            Gadfly.Guide.xlabel("Eigenvalue Estimate"),
            Gadfly.Guide.ylabel("Magnitude of ψ*ψ"),
            Gadfly.Guide.title("Magnitude of wavefunction squared vs. Estimated Eigenvalue"),
            Gadfly.Guide.xticks(ticks=[intrvl[i][1]:10intrvl[i][2]:intrvl[i][3]])
        )
        # save as SVG
        Gadfly.draw(Gadfly.SVG("p2eigvals-$(intrvl[i][1])-to-$(intrvl[i][3]).svg", 20Gadfly.cm, 12Gadfly.cm), p)
        push!(plots,p)

    end

    return plots

end


text:
#=
   Binomial Trial Count.

   From infer.net test/Tests/BlogTests.cs
   """
   // Example of inferring the size of a population.
   // Reference: A. Raftery, "Inference for the binomial N parameter: A hierarchical Bayes approach", Biometrika 1988
   // http://pluto.huji.ac.il/~galelidan/52558/Material/Raftery.pdf

   ....
   ExpectationPropagation
   impala
   theta = Beta(80.24,91.17)[mean=0.4681]
   N mode = 44, median = 45
   waterbuck
   theta = Beta(238,233.8)[mean=0.5044]
   N mode = 124, median = 125
   impala
   N mode = 37, median = 65
   waterbuck
   N mode = 122, median = 208

   VariationalMessagePassing
   impala
   theta = Beta(106,370.9)[mean=0.2223]
   N mode = 94, median = 94
   waterbuck
   theta = Beta(316,393)[mean=0.4457]
   N mode = 141, median = 141
   impala
   N mode = 37, median = 65
   waterbuck
   N mode = 122, median = 208
   """

   Model 1: Impala
   Summary Statistics
   parameters      mean       std   naive_se      mcse        ess      rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64    Float64   Float64       Float64

     beta[1]   43.9012   26.9240     0.2692    1.6380   232.9437    1.0011      110.3999
     beta[2]   59.2836   25.3251     0.2533    1.2813   328.1459    1.0000      155.5194
       theta    0.4091    0.1625     0.0016    0.0090   272.1842    1.0010      128.9972
           n   73.4049   87.6808     0.8768    4.6531   335.0041    1.0000      158.7697

   Model 2: Waterbuck
   Summary Statistics
   parameters       mean        std   naive_se      mcse        ess      rhat   ess_per_sec
      Symbol    Float64    Float64    Float64   Float64    Float64   Float64       Float64

     beta[1]    42.1176    25.6791     0.2568    1.7376   152.0795    1.0034       90.3085
     beta[2]    63.5964    24.1648     0.2416    1.7421   106.0723    1.0363       62.9883
       theta     0.3738     0.1557     0.0016    0.0104   164.0217    1.0003       97.4000
           n   215.9076   134.5262     1.3453    9.2679   136.5959    1.0000       81.1139

   Cf ~/webppl/binomial_trial_count.wppl

=#

using Turing
using Plots, StatsPlots
include("jl_utils.jl")

@model function binomial_trial_count(data)
    maxN = 1000
    beta ~ filldist(Uniform(0.1,100),2)
    theta ~ Beta(beta[1],beta[2])

    n ~ Categorical(simplex([1/i for i in 1:maxN]))
    data ~ filldist(Binomial(n, theta), length(data))

end

impala = [15,20,21,23,26]
waterbuck = [53,57,66,67,72]
println("Model 1: Impala")
model = binomial_trial_count(impala)

num_chains = 4
# chs = sample(model, Prior(), 1000)
# chs = sample(model, MH(), 10_000)
# chs = sample(model, PG(5), 10_000)
# chs = sample(model, IS(), 10_000)
chs = sample(model, SMC(), 10_000)
# chs = sample(model, SMC(), MCMCThreads(), 10_000, num_chains)

# chs = sample(model, SGLD(), 10_000)
# chs = sample(model,NUTS(), 10_000)
# chs = sample(model,HMC(0.01,5), 10_000)

display(chs)
# display(plot(chs))

println("\n\nModel 2: Waterbuck")
model = binomial_trial_count(waterbuck)

num_chains = 4
# chs = sample(model, Prior(), 1000)
# chs = sample(model, MH(), 10_000)
# chs = sample(model, PG(5), 10_000)
# chs = sample(model, IS(), 10_000)
chs = sample(model, SMC(), 10_000)
# chs = sample(model, SMC(), MCMCThreads(), 10_000, num_chains)

# chs = sample(model, SGLD(), 10_000)
# chs = sample(model,NUTS(), 10_000)
# chs = sample(model,HMC(0.01,5), 10_000)

display(chs)
# display(plot(chs))

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
function linkbuildingheights(env,project,activetiles)
    # Find active tile data paths
    activetilespaths(pth,tileReg,tile) = joinpath(pth,tileReg.tileRefDirs[tileReg.tileRefs.==tile][1])
    tilepaths = activetilespaths.((string(env.paths[:OS]),),(project.tileRegister,),activetiles)

    function loadbha(tilepath,tile)
        # Identify subdir
        dirs = readdir(tilepath)
        dir = dirs[occursin.("building_heights",dirs)][1]
        # Identify lower subdir
        pth = joinpath(tilepath,dir,lowercase(tile[1:2]))
        fname = readdir(pth)
        # Identify full path
        fullpath = joinpath(pth,fname[1])
        # Load
        bha = CSV.read(fullpath,header=false,DataFrame)
    end
    bha = loadbha.(tilepaths,activetiles)
    bha = vcat(bha...)

    # Parse OSGB as int
    parsereplace(str,str_rm,T) = parse(T,replace(str,str_rm=>""))
    bha[!,:osgb] = parsereplace.(bha.Column1,"osgb",Int64)
    # Find indices in BHA data, for all master data (based on OSGB)
    findfirst_arr(arr,val) = findfirst(arr.==val)
    i = findfirst_arr.((bha.osgb,),project.dat["master"].osgb)

    # For non-missing entries in master stock data, fill BHA data
    i_linked = findall(isnothing.(i).==false)
    L = nrow(project.dat["master"])
    project.dat["master"][!,:AbsHMin] = fill!(Array{Union{Missing,Float32}}(undef,L),missing)
    project.dat["master"][!,:AbsH2] = fill!(Array{Union{Missing,Float32}}(undef,L),missing)
    project.dat["master"][!,:AbsHMax] = fill!(Array{Union{Missing,Float32}}(undef,L),missing)
    project.dat["master"][!,:RelH2] = fill!(Array{Union{Missing,Float32}}(undef,L),missing)
    project.dat["master"][!,:RelHMax] = fill!(Array{Union{Missing,Float32}}(undef,L),missing)

    project.dat["master"].AbsHMin[i_linked] .= bha.Column5[i[i_linked]]
    project.dat["master"].AbsH2[i_linked] .= bha.Column6[i[i_linked]]
    project.dat["master"].AbsHMax[i_linked] .= bha.Column7[i[i_linked]]
    project.dat["master"].RelH2[i_linked] .= bha.Column8[i[i_linked]]
    project.dat["master"].RelHMax[i_linked] .= bha.Column9[i[i_linked]]

    return project
end


text:
# irflpest.m
# Luca Brugnolini 2014
# Based on Kilian code
# This function estimates reduced form IRFs by local projections

function irflpest(ys, yt, x)

# Input : dependent variables ys, lagged RHS variables (yt and x)
# Output: rf IRF for each horizon (K by K)
# Impulse Response by Local Projection

Mx = eye(size(x,1))-x/(x'*x)*x';         # annhilation matrix
B1 = (yt'*Mx*yt)\(yt'*Mx*ys);            # rf IRF by Local Projection

return B1

end


text:
using Copulas
using Base.Test

tests = [
         "utils_test.jl",
         "paircop_test.jl",
         "condTree_test.jl",
         "graphviz_test.jl",
         "vine_test.jl",
         "rvine_matrix_test.jl",
         ## "pcc_test.jl",
         ## "vineConstruction_test.jl",
         "ParamPC_tests.jl",
         "ParamPC_Cpp_tests.jl",
         "pitmatrix_test.jl"
         ]

println("Running copula tests:")

for t in tests
    println("")
    println("----------------------")
    println(" * $(t)")
    println("----------------------")
    println("")
    include(string(Pkg.dir("Copulas"), "/test/", t))
end


text:
export Geometry

function Geometry(con::H5SeisCon)
  if isnothing(con.seis)
    @error "Seis is Nothing"
    return
  end

  if isnothing(con.pkey) || isnothing(con.xkey) || isnothing(con.ykey) || isnothing(con.zkey)
    @error "pkey/xkey/ykey/zkey is Nothing"
    return
  end

  if isnothing(con.pkeyvals)
    @error "pkeyvals is Nothing"
    return
  end

  dt = Float32(abs(con.seis.getSampRate("ms")))
  nt = Integer(con.seis.getNSamp())
  t = Float32((nt-1)*dt)

  nsrc = length(con.pkeyvals)
  xCell = Vector{Vector{Float32}}(undef, nsrc)
  yCell = Vector{Vector{Float32}}(undef, nsrc)
  zCell = Vector{Vector{Float32}}(undef, nsrc)
  dtCell = Vector{Float32}(undef, nsrc)
  ntCell = Vector{Integer}(undef, nsrc)
  tCell = Vector{Float32}(undef, nsrc)
  keylist = [con.pkey]
  for block in 1:nsrc
    minlist = [con.pkeyvals[block]]
    maxlist = [con.pkeyvals[block]]
    _, _, ind = con.seis.getSortedData(keylist, minlist, maxlist, 0, 0)
    xy = Float32.(con.seis.getXYTraceHeaders([con.xkey, con.ykey], ind, "m", false))  # MUST BE TRUE

    xCell[block] = xy[:,1]
    yCell[block] = xy[:,2]
    zCell[block] = Float32.(con.seis.getTraceHeader([con.zkey], ind, [con.seis.getLengthUnits()], ["m"])[:])
    dtCell[block] = dt
    ntCell[block] = nt
    tCell[block] = t
  end

  return GeometryIC{Float32}(xCell, yCell, zCell, dtCell, ntCell, tCell)
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

"""A Cluster is a repeating pattern of indices defined in a lattice called the element lattice. The cluster of elements has its own lattice, which by definition is different from the element lattice unless each cluster consists of a single element. For example this defines a Cluster of three hexagonal elements:

```
julia> function hex3cluster()
    clusterelts = SVector((0,0),(-1,0),(-1,1))
    eltlattice = HexBasis1()
    clusterbasis = LatticeBasis(( -1,2),(2,-1))
    return LatticeCluster(clusterbasis,eltlattice,clusterelts)
end

julia> lattice = hex3cluster()
LatticeCluster{3, 2, Float64, LatticeBasis{2, Int64}, HexBasis1{2, Float64}}(LatticeBasis{2, Int64}([-1 2; 2 -1]), HexBasis1{2, Float64}(), [(0, 0), (-1, 0), (-1, 1)])

julia> lattice[1,1]  #returns the locations of the 3 elements in the cluster at the cluster offset of 1,1
2×3 SMatrix{2, 3, Float64, 6} with indices SOneTo(2)×SOneTo(3):
 1.0  -0.5        1.0
 1.0   0.133975  -0.732051

```
"""
struct LatticeCluster{N1, N, T<:Real, B1<:Basis{N,Int},B2<:Basis{N,T}}
    clusterbasis::B1 #this basis defines the offsets of the clusters
    elementbasis::B2 #this basis defines the underlying lattice

    clusterelements::SVector{N1,NTuple{N,Int}} #vector containing the lattice indices of the elements in the cluster. Each column is one lattice coordinate. These indices are assumed to be offsets from the origin of the lattice.

    LatticeCluster(clusterbasis::B1,eltlattice::B2,clusterelements::SVector{N1,NTuple{N,Int}}) where{N1,N,T<:Real,B1<:Basis{N,Int},B2<:Basis{N,T}} = new{N1,N,T,B1,B2}(clusterbasis,eltlattice,clusterelements)
end
export LatticeCluster

clusterelements(a::LatticeCluster) = a.clusterelements
export clusterelements
elementbasis(a::LatticeCluster) = a.elementbasis
clustersize(a::LatticeCluster) = length(a.clusterelements)

"""returns the positions of every element in a cluster given the cluster indices"""
function Base.getindex(A::LatticeCluster{N1,N,T,B1,B2}, indices::Vararg{Int, N}) where{N1,N,T,B1<:Basis{N,Int},B2<:Basis{N,T}}
    clusteroffset = A.clusterbasis[indices...]
    temp = MMatrix{N,N1,T}(undef)
    for i in 1:N1
        temp[:,i] = A.elementbasis[A.clusterelements[i]...] + clusteroffset
    end
    return SMatrix{N,N1,T}(temp) #positions of the cluster elements, not the lattice coordinates.

end

"""Can access any cluster in a cluster lattice so the range of indices is unlimited"""
function Base.size(::LatticeCluster{N1,N}) where{N1,N}
    return ntuple((i)->Base.IsInfinite(),N)
end

Base.setindex!(A::Basis{N}, v, I::Vararg{Int, N}) where{N} = nothing #can't set lattice points. Might want to throw an exception instead.

""" returns the lattice indices of the elements in the cluster. These are generally not the positions of the elements"""
function clustercoordinates(a::LatticeCluster{N1,N},indices::Vararg{Int,N}) where{N1,N}
    clusteroffset = a.clusterbasis[indices...]
    temp = MMatrix{N,N1,Int}(undef)
    for i in 1:N1
        temp[:,i] = SVector{N,Int}(a.clusterelements[i]) + clusteroffset
    end
    return SMatrix{N,N1,Int}(temp) #positions of the cluster elements, not the lattice coordinates.
end
export clustercoordinates

"""
May want to have many properties associated with the elements in a cluster, which is why properties is represented as a DataFrame. The DataFrame in the properties field should have as many rows as there are elements in a cluster. At a minimum it must have a :Color and a :Name column.

Example:
lenslets = DataFrame()
"""
struct ClusterWithProperties{N1,N,T}
    cluster::LatticeCluster{N1,N,T}
    properties::DataFrame
end
export ClusterWithProperties

Base.getindex(A::ClusterWithProperties{N1,N,T}, indices::Vararg{Int,N}) where{N1,N,T} = cluster(A)[indices...]
Base.setindex!(A::ClusterWithProperties, v, I::Vararg{Int, N}) where{N} = nothing #can't set lattice points. Might want to throw an exception instead.

properties(a::ClusterWithProperties) = a.properties
export properties
cluster(a::ClusterWithProperties) = a.cluster
export cluster
clustercoordinates(a::ClusterWithProperties,indices...) = clustercoordinates(cluster(a),indices...)
elementbasis(a::ClusterWithProperties) = elementbasis(cluster(a))
export elementbasis
clustersize(a::ClusterWithProperties) = clustersize(a.cluster)
export clustersize


text:
# O3DFT.jl
#  - Orbital 3 Density Functional Theory

push!(LOAD_PATH,"../src/") # load module from local directory
using O3DFT

function main()
    # Does anyone know what units we should be in?
    for n in 1E25:1E25:1E26 # density /m^3
        @printf("n: %g \t T(n): %g\t T_fnderiv(n): %g\t E_xc(n): %g\n",n,ThomasFermi_T(n),ThomasFermi_T_fnderiv(n),ThomasFermi_Exc(n))
    end

    for n=1:5:20
        AtomicTest(n,verbose=false) #Oh uh huh make it magnificent
    end
end

main() # C-party like the mid-90s


text:
# SegmentTree map operations

const SegmentMap{K, V} = SegmentTree{K, SegmentValue{K, V}}

function Base.setindex!(t::SegmentBTree{K, V, B}, value,
                        key::Tuple{Any, Any}) where {K, V, B}
    push!(t, V(key[1], key[2], value), true)
end

function Base.setindex!(t::SegmentBTree{K, V, B}, value,
                        key::Segment{K}) where {K, V, B}
    push!(t, V(key.first, key.last, value), true)
end

function Base.setindex!(t::SegmentBTree{K, V, B}, value,
                        first, last) where {K, V, B}
    push!(t, V(first, last, value), true)
end

function Base.getindex(t::SegmentBTree{K, V, B}, key::AbstractSegment{K}) where {K, V, B}
    return _getindex(t.root, key)
end

function _getindex(t::InternalNode{K, V, B}, key::AbstractSegment{K}) where {K, V, B}
    i = findidx(t, key)
    if 1 <= length(t) - 1 && key >= t.keys[i]
        return _getindex(t.children[i+1], key)
    else
        return _getindex(t.children[i], key)
    end
end

function _getindex(t::LeafNode{K, V, B}, key::AbstractSegment{K}) where {K, V, B}
    i = findidx(t, key)
    if 1 <= i <= length(t) && first(t.entries[i]) == first(key) &&
        last(t.entries[i]) == last(key)
        return t.entries[i]
    else
        error(KeyError((first(key), last(key))))
    end
end

function Base.get(t::SegmentBTree{K, V, B}, key::Tuple{Any, Any}, default) where {K, V, B}
    return _get(t.root, Segment{K}(key[1], key[2]), default)
end

function Base.get(t::SegmentBTree{K, V, B}, key::AbstractSegment{K}, default) where {K, V, B}
    return _get(t.root, key, default)
end

function _get(t::InternalNode{K, V, B}, key::AbstractSegment{K}, default) where {K, V, B}
    i = findidx(t, key)
    if 1 <= length(t) - 1 && key >= t.keys[i]
        return _get(t.children[i+1], key, default)
    else
        return _get(t.children[i], key, default)
    end
end

function _get(t::LeafNode{K, V, B}, key::AbstractSegment{K}, default) where {K, V, B}
    i = findidx(t, key)
    if 1 <= i <= length(t) &&
        first(t.entries[i]) == first(key) && last(t.entries[i]) == last(key)
        return t.entries[i]
    else
        return default
    end
end

function Base.get!(t::SegmentBTree{K, V, B}, key::Tuple{Any, Any}, default) where {K, V, B}
    return push!(t, V(key[1], key[2], default), true, false)
end

function Base.get!(t::SegmentBTree{K, V, B}, key::AbstractSegment{K}, default) where {K, V, B}
    return push!(t, V(first(key), last(key), default), true, false)
end

function Base.delete!(t::SegmentBTree{K, V, B}, first, last) where {K, V, B}
    return deletefirst!(t, first, last)
end

function Base.delete!(t::SegmentBTree{K, V, B}, key::Tuple{Any, Any}) where {K, V, B}
    return deletefirst!(t, key)
end

function Base.delete!(t::SegmentBTree{K, V, B}, key::Segment{K}) where {K, V, B}
    return deletefirst!(t, key)
end

struct SegmentKeyIterator{K, V, B}
    t::SegmentBTree{K, V, B}
end

Base.eltype(::Type{SegmentKeyIterator{K,V,B}}) where {K,V,B} = Segment{K}
Base.IteratorSize(::Type{SegmentKeyIterator{K,V,B}}) where {K,V,B} = Base.SizeUnknown()

function Base.keys(t::SegmentBTree)
    return SegmentKeyIterator(t)
end

function Base.iterate(
        it::SegmentKeyIterator{K, V, B},
        state::SegmentBTreeIteratorState{K, V, B}=iterinitstate(it)) where {K, V, B}
    if state.leaf === nothing || isempty(notnothing(state.leaf))
        return nothing
    end
    leaf = state.leaf
    key = Segment{K}(first(leaf.entries[state.i]), last(leaf.entries[state.i]))
    if state.i < length(leaf)
        state = SegmentBTreeIteratorState{K, V, B}(leaf, state.i + 1)
    else
        state = SegmentBTreeIteratorState{K, V, B}(leaf.right, 1)
    end
    return key, state
end

function iterinitstate(it::SegmentKeyIterator{K, V, B}) where {K, V, B}
    # traverse to the first leaf node
    node = it.t.root
    while !isa(node, LeafNode{K, V, B})
        node = node.children[1]
    end
    return SegmentBTreeIteratorState(node, 1)
end

struct SegmentValueIterator{K, V <: SegmentValue, B}
    t::SegmentBTree{K, V, B}
end

Base.eltype(::Type{SegmentValueIterator{K,V,B}}) where {K,V,B} = valtype(V)
Base.IteratorSize(::Type{SegmentValueIterator{K,V,B}}) where {K,V,B} = Base.SizeUnknown()

function Base.values(t::SegmentBTree{K, V, B}) where {K, V<:SegmentValue, B}
    return SegmentValueIterator(t)
end

function Base.iterate(
        it::SegmentValueIterator{K, V, B},
        state::SegmentBTreeIteratorState{K, V, B}=iterinitstate(it)) where {K, V, B}
    if state.leaf === nothing || isempty(notnothing(state.leaf))
        return nothing
    end
    leaf = state.leaf
    value = leaf.entries[state.i].value
    if state.i < length(leaf)
        state = SegmentBTreeIteratorState{K, V, B}(leaf, state.i + 1)
    else
        state = SegmentBTreeIteratorState{K, V, B}(leaf.right, 1)
    end
    return value, state
end

function iterinitstate(it::SegmentValueIterator{K, V, B}) where {K, V, B}
    # traverse to the first leaf node
    node = it.t.root
    while !isa(node, LeafNode{K, V, B})
        node = node.children[1]
    end
    return SegmentBTreeIteratorState(node, 1)
end


text:
# Copied from base for use on 0.6 (only in base as of 0.7).
# Copyright (c) 2009-2016: Jeff Bezanson, Stefan Karpinski, Viral B. Shah, and other contributors:
using Base.LibGit2: GitHash, SignatureStruct, AbstractGitObject, REFCOUNT

"""
    LibGit2.BlameOptions

Matches the [`git_blame_options`](https://libgit2.github.com/libgit2/#HEAD/type/git_blame_options) struct.
"""
@Base.LibGit2.kwdef struct BlameOptions
    version::Cuint                    = 1
    flags::UInt32                     = 0
    min_match_characters::UInt16      = 20
    newest_commit::GitHash
    oldest_commit::GitHash
    min_line::Csize_t                 = 1
    max_line::Csize_t                 = 0
end

for (typ, owntyp, sup, cname) in [
    (:GitBlame,          :GitRepo,              :AbstractGitObject, :git_blame),
    ]

    if owntyp === nothing
        @eval mutable struct $typ <: $sup
            ptr::Ptr{Void}
            function $typ(ptr::Ptr{Void}, fin::Bool=true)
                # fin=false should only be used when the pointer should not be free'd
                # e.g. from within callback functions which are passed a pointer
                @assert ptr != C_NULL
                obj = new(ptr)
                if fin
                    Threads.atomic_add!(REFCOUNT, UInt(1))
                    finalizer(obj, Base.close)
                end
                return obj
            end
        end
    else
        @eval mutable struct $typ <: $sup
            owner::$owntyp
            ptr::Ptr{Void}
            function $typ(owner::$owntyp, ptr::Ptr{Void}, fin::Bool=true)
                @assert ptr != C_NULL
                obj = new(owner, ptr)
                if fin
                    Threads.atomic_add!(REFCOUNT, UInt(1))
                    finalizer(obj, Base.close)
                end
                return obj
            end
        end
        if isa(owntyp, Expr) && owntyp.args[1] == :Nullable
            @eval begin
                $typ(ptr::Ptr{Void}, fin::Bool=true) = $typ($owntyp(), ptr, fin)
                $typ(owner::$(owntyp.args[2]), ptr::Ptr{Void}, fin::Bool=true) =
                    $typ($owntyp(owner), ptr, fin)
            end
        end
    end
    @eval function Base.close(obj::$typ)
        if obj.ptr != C_NULL
            ccall(($(string(cname, :_free)), :libgit2), Void, (Ptr{Void},), obj.ptr)
            obj.ptr = C_NULL
            if Threads.atomic_sub!(REFCOUNT, UInt(1)) == 1
                # will the last finalizer please turn out the lights?
                ccall((:git_libgit2_shutdown, :libgit2), Cint, ())
            end
        end
    end
end

"""
    LibGit2.BlameHunk

Matches the [`git_blame_hunk`](https://libgit2.github.com/libgit2/#HEAD/type/git_blame_hunk) struct.
The fields represent:
    * `lines_in_hunk`: the number of lines in this hunk of the blame.
    * `final_commit_id`: the [`GitHash`](@ref) of the commit where this section was last changed.
    * `final_start_line_number`: the *one based* line number in the file where the
       hunk starts, in the *final* version of the file.
    * `final_signature`: the signature of the person who last modified this hunk. You will
       need to pass this to [`Signature`](@ref) to access its fields.
    * `orig_commit_id`: the [`GitHash`](@ref) of the commit where this hunk was first found.
    * `orig_path`: the path to the file where the hunk originated. This may be different
       than the current/final path, for instance if the file has been moved.
    * `orig_start_line_number`: the *one based* line number in the file where the
       hunk starts, in the *original* version of the file at `orig_path`.
    * `orig_signature`: the signature of the person who introduced this hunk. You will
       need to pass this to [`Signature`](@ref) to access its fields.
    * `boundary`: `'1'` if the original commit is a "boundary" commit (for instance, if it's
       equal to an oldest commit set in `options`).
"""
@Base.LibGit2.kwdef struct BlameHunk
    lines_in_hunk::Csize_t

    final_commit_id::GitHash
    final_start_line_number::Csize_t
    final_signature::Ptr{SignatureStruct}

    orig_commit_id::GitHash
    orig_path::Cstring
    orig_start_line_number::Csize_t
    orig_signature::Ptr{SignatureStruct}

    boundary::Char
end

function GitBlame(repo::GitRepo, path::AbstractString; options::BlameOptions=BlameOptions())
    blame_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @Base.LibGit2.check ccall((:git_blame_file, :libgit2), Cint,
                  (Ptr{Ptr{Void}}, Ptr{Void}, Cstring, Ptr{BlameOptions}),
                   blame_ptr_ptr, repo.ptr, path, Ref(options))
    return GitBlame(repo, blame_ptr_ptr[])
end

function counthunks(blame::GitBlame)
    return ccall((:git_blame_get_hunk_count, :libgit2), Int32, (Ptr{Void},), blame.ptr)
end

function Base.getindex(blame::GitBlame, i::Integer)
    if !(1 <= i <= counthunks(blame))
        throw(BoundsError(blame, (i,)))
    end
    hunk_ptr = ccall((:git_blame_get_hunk_byindex, :libgit2),
                      Ptr{BlameHunk},
                      (Ptr{Void}, Csize_t), blame.ptr, i-1)
    return unsafe_load(hunk_ptr)
end

function Base.show(io::IO, blame_hunk::BlameHunk)
    println(io, "GitBlameHunk:")
    println(io, "Original path: ", unsafe_string(blame_hunk.orig_path))
    println(io, "Lines in hunk: ", blame_hunk.lines_in_hunk)
    println(io, "Final commit oid: ", blame_hunk.final_commit_id)
    print(io, "Final signature: ")
    show(io, Signature(blame_hunk.final_signature))
    println(io)
    println(io, "Original commit oid: ", blame_hunk.orig_commit_id)
    print(io, "Original signature: ")
    show(io, Signature(blame_hunk.orig_signature))
end


text:
#####################################################################
##                      Data Type Declarations                     ##
#####################################################################

struct StateSpace{T, MT<:AbstractMatrix{T}} <: LTISystem
    A::MT
    B::MT
    C::MT
    D::MT
    Ts::Float64
    nx::Int
    nu::Int
    ny::Int
    function StateSpace{T, MT}(A::MT, B::MT,
            C::MT, D::MT, Ts::Float64) where {T, MT <: AbstractMatrix{T}}
        nx = size(A, 1)
        nu = size(B, 2)
        ny = size(C, 1)

        if size(A, 2) != nx && nx != 0
            error("A must be square")
        elseif size(B, 1) != nx
            error("B must have the same row size as A")
        elseif size(C, 2) != nx
            error("C must have the same column size as A")
        elseif nu != size(D, 2)
            error("D must have the same column size as B")
        elseif ny != size(D, 1)
            error("D must have the same row size as C")
        end

        # Validate sampling time
        if Ts < 0 && Ts != -1
            error("Ts must be either a positive number, 0
                   (continuous system), or -1 (unspecified)")
        end
        new{T, MT}(A, B, C, D, Ts, nx, nu, ny)
    end
end

function StateSpace{T,MT}(A::AbstractArray, B::AbstractArray, C::AbstractArray, D::AbstractArray, Ts::Real) where {T, MT <: AbstractMatrix{T}}
        return StateSpace{T,Matrix{T}}(MT(to_matrix(T, A)), MT(to_matrix(T, B)), MT(to_matrix(T, C)),
            MT(to_matrix(T, D)), Float64(Ts))
end

function StateSpace(A::AbstractArray, B::AbstractArray, C::AbstractArray, D::AbstractArray, Ts::Real)
        # TODO: change back in 0.7 T = promote_type(eltype(A),eltype(B),eltype(C),eltype(D))
        T = promote_type(promote_type(eltype(A),eltype(B)), promote_type(eltype(C),eltype(D)))
        @assert (typeof(to_matrix(T, A)) == typeof(to_matrix(T, B)) == typeof(to_matrix(T, C)) == typeof(to_matrix(T, D)))
        return StateSpace{T,Matrix{T}}(to_matrix(T, A), to_matrix(T, B), to_matrix(T, C),
            to_matrix(T, D), Float64(Ts))
end

# Getter functions
get_A(sys::StateSpace) = sys.A
get_B(sys::StateSpace) = sys.B
get_C(sys::StateSpace) = sys.C
get_D(sys::StateSpace) = sys.D

get_Ts(sys::StateSpace) = sys.Ts

ssdata(sys::StateSpace) = get_A(sys), get_B(sys), get_C(sys), get_D(sys)

# Funtions for number of intputs, outputs and states
ninputs(sys::StateSpace) = size(get_D(sys), 2)
noutputs(sys::StateSpace) = size(get_D(sys), 1)
nstates(sys::StateSpace) = size(get_A(sys), 1)

#####################################################################
##                         Math Operators                          ##
#####################################################################

## EQUALITY ##
function ==(sys1::StateSpace, sys2::StateSpace)
    return all(getfield(sys1, f) == getfield(sys2, f) for f in fieldnames(StateSpace))
end

## Approximate ##
function isapprox(sys1::StateSpace, sys2::StateSpace)
    return all(getfield(sys1, f) ≈ getfield(sys2, f) for f in fieldnames(StateSpace))
end

## ADDITION ##
function +(s1::StateSpace{T,MT}, s2::StateSpace{T,MT}) where {T, MT}
    #Ensure systems have same dimensions
    if size(s1) != size(s2)
        error("Systems have different shapes.")
    elseif s1.Ts != s2.Ts
        error("Sampling time mismatch")
    end

    A = [s1.A                   fill(zero(T), nstates(s1), nstates(s2));
         fill(zero(T), nstates(s2), nstates(s1))        s2.A]
    B = [s1.B ; s2.B]
    C = [s1.C s2.C;]
    D = [s1.D + s2.D;]

    return StateSpace(A, B, C, D, s1.Ts)
end

+(sys::StateSpace, n::Number) = StateSpace(sys.A, sys.B, sys.C, sys.D .+ n, sys.Ts)
+(n::Number, sys::StateSpace) = +(sys, n)

## SUBTRACTION ##
-(sys1::StateSpace, sys2::StateSpace) = +(sys1, -sys2)
-(sys::StateSpace, n::Number) = +(sys, -n)
-(n::Number, sys::StateSpace) = +(-sys, n)

## NEGATION ##
-(sys::StateSpace) = StateSpace(sys.A, sys.B, -sys.C, -sys.D, sys.Ts)

## MULTIPLICATION ##
function *(sys1::StateSpace{T,MT}, sys2::StateSpace{T,MT}) where {T, MT}
    #Check dimension alignment
    #Note: sys1*sys2 = y <- sys1 <- sys2 <- u
    if sys1.nu != sys2.ny
        error("sys1*sys2: sys1 must have same number of inputs as sys2 has outputs")
    elseif sys1.Ts != sys2.Ts
        error("Sampling time mismatch")
    end

    A = [sys1.A    sys1.B*sys2.C;
         fill(zero(T), sys2.nx, sys1.nx)  sys2.A]
    B = [sys1.B*sys2.D ; sys2.B]
    C = [sys1.C   sys1.D*sys2.C;]
    D = [sys1.D*sys2.D;]

    return StateSpace{T,MT}(A, B, C, D, sys2.Ts)
end

*(sys::StateSpace, n::Number) = StateSpace(sys.A, sys.B, sys.C*n, sys.D*n, sys.Ts)
*(n::Number, sys::StateSpace) = *(sys, n)

## DIVISION ##
/(sys1::StateSpace, sys2::StateSpace) = sys1*inv(sys2)

function /(n::Number, sys::StateSpace)
    # Ensure s.D is invertible
    A, B, C, D = ssdata(sys)
    Dinv = try
        inv(D)
    catch
        error("D isn't invertible")
    end
    return StateSpace(A - B*Dinv*C, B*Dinv, -n*Dinv*C, n*Dinv, get_Ts(sys))
end

Base.inv(sys::StateSpace) = 1/sys
/(sys::StateSpace, n::Number) = StateSpace(sys.A, sys.B, sys.C/n, sys.D/n, sys.Ts)

Base.:^(sys::StateSpace, p::Integer) = Base.power_by_squaring(sys, p)

#####################################################################
##                       Indexing Functions                        ##
#####################################################################
Base.ndims(::StateSpace) = 2 # NOTE: Also for SISO systems?
Base.size(sys::StateSpace) = (noutputs(sys), ninputs(sys)) # NOTE: or just size(get_D(sys))
Base.size(sys::StateSpace, d) = d <= 2 ? size(sys)[d] : 1
Base.eltype(::Type{S}) where {S<:StateSpace} = S

function Base.getindex(sys::StateSpace, inds...)
    if size(inds, 1) != 2
        error("Must specify 2 indices to index statespace model")
    end
    rows, cols = index2range(inds...) # FIXME: ControlSystems.index2range(inds...)
    return StateSpace(copy(sys.A), sys.B[:, cols], sys.C[rows, :], sys.D[rows, cols], sys.Ts)
end

#####################################################################
##                        Display Functions                        ##
#####################################################################

function _string_mat_with_headers(X::Matrix)
    #mat = [[""] reshape(cols,1,length(cols));
    #       rows X]
    p = (io, m) -> Base.print_matrix(io, m)
    return replace(sprint(p, X), "\"" => " ")
end

Base.print(io::IO, sys::StateSpace) = show(io, sys)

function Base.show(io::IO, sys::StateSpace)
    # Compose the name vectors
    #inputs = format_names(s.inputnames, "u", "?")
    #outputs = format_names(s.outputnames, "y", "?")
    #println(io, "StateSpace:")
    println(io, typeof(sys))
    if nstates(sys) > 0
        #states = format_names(s.statenames, "x", "?")
        println(io, "A = \n", _string_mat_with_headers(sys.A))
        println(io, "B = \n", _string_mat_with_headers(sys.B))
        println(io, "C = \n", _string_mat_with_headers(sys.C))
    end
    println(io, "D = \n", _string_mat_with_headers(sys.D), "\n")
    # Print sample time
    if sys.Ts > 0
        println(io, "Sample Time: ", sys.Ts, " (seconds)")
    elseif sys.Ts == -1
        println(io, "Sample Time: unspecified")
    end
    # Print model type
    if nstates(sys) == 0
        print(io, "Static gain") # NOTE: Not quite...still has a time type
    elseif iscontinuous(sys)
        print(io, "Continuous-time state-space model")
    else
        print(io, "Discrete-time state-space model")
    end
end

"""
`minsys = minreal(s::StateSpace, tol=sqrt(eps()))` is implemented via `baltrunc` and returns a system on diagonal form.
"""
function minreal(s::StateSpace, tol=sqrt(eps()))
    s = baltrunc(s, atol=tol, rtol = 0)[1]
    try
        return diagonalize(s)
    catch
        error("Minreal only implemented for diagonalizable systems.")
    end
end

"""
`dsys = diagonalize(s::StateSpace, digits=12)` Diagonalizes the system such that the A-matrix is diagonal.
"""
function diagonalize(s::StateSpace, digits = 12)
    r = x -> round(x, digits=digits)
    S,V = eigen(s.A)
    try
        A = diagm(0 => S)
        B = V\s.B
        C = s.C*V
        D = s.D
        return ss(A,B,C,D)
    catch e
        error("System not diagonalizable", e)
    end
end


text:
entropy(s) = -sum(x -> x * log(2, x), count(x -> x == c, s) / length(s) for c in unique(s))
@show entropy("1223334444")
@show entropy([1, 2, 3, 1, 2, 1, 2, 3, 1, 2, 3, 4, 5])

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# # bard1.mod	QQR2-MN-8-5
# # Original AMPL coding by Sven Leyffer
#
# # An MPEC from J.F. Bard, Convex two-level optimization,
# # Mathematical Programming 40(1), 15-27, 1988.
#
# # Number of variables:   2 + 3 slack + 3 multipliers
# # Number of constraints: 4
#
# var x >= 0;
# var y >= 0;
#
# # ... multipliers
# var l{1..3};
#
# minimize f:(x - 5)^2 + (2*y + 1)^2;
#
# subject to
#
#    KKT:    2*(y-1) - 1.5*x + l[1] - l[2]*0.5 + l[3] = 0;
#
#    lin_1:  0 <= 3*x - y - 3        complements l[1] >= 0;
#    lin_2:  0 <= - x + 0.5*y + 4    complements l[2] >= 0;
#    lin_3:  0 <= - x - y + 7        complements l[3] >= 0;

using JuMP, Complementarity
using Ipopt
using Test

@testset "bard1_test.jl" begin

    m = Model(Ipopt.Optimizer)

    @variable(m, x>=0)
    @variable(m, y>=0)
    @variable(m, l[1:3])

    @NLobjective(m, Min, (x - 5)^2 + (2*y + 1)^2)

    # @NLconstraint(m, 2*(y-1) - 1.5*x + l[1] - l[2]*0.5 + l[3] == 0)
    @NLconstraint(m, 2*(y-1) - 1.5*x + l[1] + l[2]*0.5 - l[3] == 0)

    # Original Form
    # @complements(m, 0 <= 3*x - y - 3,        l[1] >= 0)
    # @complements(m, 0 <= - x + 0.5*y + 4,    l[2] >= 0)
    # @complements(m, 0 <= - x - y + 7,        l[3] >= 0)
    @complements(m, 0 >= - 3*x + y + 3,        l[1] >= 0, smooth)
    @complements(m, 0 <= - x + 0.5*y + 4,      l[2] <= 0, simple)
    @complements(m, 0 >= + x + y - 7,          l[3] <= 0, smooth)

    # just testing
    @macroexpand(@complements(m, 3*x - y - 3,    0 <=  l[1] <= 5))
    @macroexpand(@complements(m, 3*x - y - 3,    5 >=  l[1] >= 0))

    @test_throws LoadError @macroexpand(@complements(m,  0 <= 3*x - y - 3 <=10,  l[1]))
    @test_throws LoadError @macroexpand(@complements(m, 10 >= 3*x - y - 3 >= 5,  l[1]))

    JuMP.optimize!(m)
    @show JuMP.value.(x)
    @show JuMP.value.(y)
    @show JuMP.value.(l)
    @show JuMP.objective_value(m)
    @test isapprox(JuMP.objective_value(m), 17.0000, atol=1e-4)
end


text:
# Flux vNext
# ==========

# Optimisation Basics
# -------------------

using Flux
using Flux: step!

# Using Flux is very simple. You write a program, and we'll make tweaks to that
# program so that it gets gradually better.

# What does it mean to make a program "better"? That's up to you – your program
# returns a score, called the "loss", which determines how well the program
# is doing. Flux's job is to minimise that score. For example, let's take the
# simplest possible program, one that simply returns a constant.

w = 1
loss = () -> w

loss()

# This program doesn't look very interesting, but we can still do something
# interesting with it. The core function that optimises programs for us is
# `step!`. We have to pass `step!` an *optimiser* called `Descent`; this basically
# tells Flux how agressive to be, but we'll talk more about that later.

opt = Descent(1)
step!(loss, opt)

# `step!` returns the same loss as before, `1`. But something more interesting
# has happened; try running `loss()` again.

loss()

# It went down! And if we keep calling `step!` in a loop, it'll keep going down.

for i = 1:10
  @show step!(loss, opt)
end

loss()

# Of course, this case is pretty easy: we can always improve the parameter by
# making `w` less.

w

# Here's something harder; now our `loss` is always positive, so it can't keep
# improving indefinitely. Things will stop improving once we hit the *minimum*
# of this function (which we happen to know is at $w = 0$, where $loss = 0$.)

w = 1
loss = () -> w^2

opt = Descent(0.2)
for i = 1:10
  @show step!(loss, opt)
end

w

# You can see that our loss gradually tends towards $0$, and so does $w$. Note,
# however, that Flux will never say: "Ok, we're done here, here's the best value
# for $w$." Though there are tools that can do this in simple cases, Flux is
# designed to scale to extremely complex problems where this is no longer
# possible. So we only make tweaks and it's up to you when to finish.

# Let's put these ideas towards something a little more interesting. Say we want
# to solve $5x = 10$, to find an $x$ that makes this true. What's our program?
# Well, to start with we want to take $f(x) = 5x$. Then our loss should be something like
# $f(x) - 10$, so that it measures how far the $f(x)$ is from where we want it to
# be. This doesn't quite work, however, since the loss will be low (negative) if
# $f(x)$ is `-Inf`! So we can use our squaring trick again here, to make
# sure that $f(x) - 10$ tends to zero.

x = 1 # Our initial guess
f = x -> 5x

opt = Descent(0.01)

for i = 1:10
  l = step!(opt) do
    (f(x) - 10)^2
  end
  @show l
end

# Our loss ended up being pretty low. How's our function looking?

5x

# That looks pretty good. So we're beginning to be able to use Flux to solve
# problems where we know what the *output* should look like, but we're not
# sure what the *input* should be to get there.

# You now arguably understand everything you need to do productive ML. But let's
# look over a few more examples to see how it looks in practice.

# Optimising Colours
# ------------------

# Just like Julia more generally, Flux has good support for custom types.
# This means we can carry out optimisation on things like colours!

# This example uses the excellent Colors.jl. Colors contains, among other
# things, a `colordiff` function which uses fancy colour theory algorithms to
# estimate the *perceptual* difference between two colours. We can use this
# directly in our loss function.

using Colors

target = RGB(1, 0, 0)
colour = RGB(1, 1, 1)
[target, colour]
#-
opt = Descent(0.01)

for i = 1:10
  step!(opt, target) do y
    colordiff(colour, y)
  end
end

[colour, target]

# `colour` started out white and is now red. That makes sense, as we've
# minimised the distance between the two colours. But we can also *maximize*
# with a simple minus sign.

colour1 = RGB(1, 1, 1)

for i = 1:10
  step!(opt, target) do y
    -colordiff(colour1, y)
  end
end

[colour1, target]

# Now we have green, a colour that's arguably very different from red. However,
# there's a subtlety here; notice what happens if we use a different colour as
# our starting point.

colour2 = RGB(0, 0, 1)

for i = 1:10
  step!(opt, target) do y
    -colordiff(colour2, y)
  end
end

[colour2, target]

# Now we have a dark blue! If we look directly at `colourdiff` we'll see that
# green is better.

colordiff(target, colour1), colordiff(target, colour2)

# So why do we get blue here? This is another case where it's important that
# Flux optimises programs through a series of small tweaks. In this case, even
# though green is better overall, making our colour slightly more green actually
# makes our score worse temporarily.

colordiff(target, RGB(0, 0, 0.4)), colordiff(target, RGB(0, 0.1, 0.4))

# This is known as a *local optimimum*. It's important to understand how Flux
# optimises programs and what this means for you, so we'll cover this in more
# detail in future.

# Keras in 5 lines
# -----------------

# [Working on making this an MNIST demo, but here's the gist of it.]

# Dummy data.

x = rand(10)
y = [1, 0]

# Logistic regresion.

using Flux: crossentropy

W = randn(2, 10)
b = zeros(2)

predict = x -> softmax(W * x .+ b)

opt = Descent(0.1)

loss = (x, y) -> crossentropy(predict(x), y)

step!(loss, opt, x, y)

# Multi-layer perceptron.

function dense(in, out, σ = identity)
  W = randn(out, in)
  b = zeros(out)
  x -> σ.(W * x .+ b)
end

chain(fs...) = x -> foldl((x, m) -> m(x), fs, init = x)

model = chain(dense(10, 5, relu), dense(5, 2), softmax)

# Doesn't quite work yet.

## step!(opt, x, y) do x, y
##   crossentropy(model(x), y)
## end

max_stars_repo_name: jstrube/model-zoo

max_issues_repo_name: jstrube/model-zoo

max_forks_repo_name: jstrube/model-zoo