text:
# mri_objects.jl

using MIRT: mri_objects

using MIRT: jim, image_geom, prompt, max_percent_diff
import MIRT: mri_objects_trap
using FFTW: fft, fftshift, ifftshift
using Plots
using Test: @test, @test_throws #, @inferred

# 2D test
#NOTE: Dirac will not plot very well.
#Just make sure it prints a gray image, and that's fine. FFT will also be screwed up.
function mri_objects_test2()
    ig = image_geom(nx = 2^7, ny = 2^7 + 2, dx = 4, offsets = :dsp)

    shift = [0.1, 0.2] .* ig.fovs
    sizes = [.15, .1] .* ig.fovs
    tests = ( # 3 separate tests:
                (:circ2, [shift' sizes[1] 2]),
                (:gauss2, [shift' sizes' 2]),
                (:rect2, [shift' sizes' 2]),
                (:dirac2,[shift' 2])
            )
    ntest = length(tests)
    pl = Matrix{Plots.Plot}(undef, ntest,3)

    for (i,test) = enumerate(tests)
        st = mri_objects(test)
        i2 = st.image(ig.xg, ig.yg, dx=ig.dx, dy=ig.dy)
        s2 = ifftshift(fft(fftshift(i2))) * abs(ig.dx * ig.dy)
        f2 = st.kspace(ig.fg...)
        isinteractive() && (@show max_percent_diff(f2,s2))
    #   @show round.(maximum.([i2, abs.(f2), abs.(s2)]), digits=2)
        pl[i,1] = jim(i2, title = "$(test[1])")
        pl[i,2] = jim(abs.(s2), title = "|fft|")
        pl[i,3] = jim(abs.(f2), title = "|kspace|")
    end
    plot(pl...)
    return plot(pl...)
end

function mri_objects_test_case1()
    ig = image_geom(nx = 2^7, ny = 2^7 + 2, dx = 4, offsets = :dsp)
    xt = mri_objects(:case1,unit = :cm).image(ig.xg, ig.yg)
    jim(xt, title = "case1cm")
end

function mri_objects_test_case4()
    ig = image_geom(nx = 2^5, ny = 2^5 + 2,nz = 2^5, dx = 1,dz = .5, offsets = :dsp)
    xt = mri_objects(:case4,unit = :cm).image(ig.xg, ig.yg, ig.zg)
    jim(xt, title = "case4cm")
end

# 3D test
# NOTE: Dirac will not plot very well.
# Just make sure it prints a gray image, and that's fine. FFT will also be screwed up.
function mri_objects_test3()
    ig = image_geom(nx = 2^7, nz = 2^5, offsets = :dsp, dx = 4, dz = 3)
    shift = [0.1, 0.2, 0.3] .* ig.fovs
    sizes = [0.15, 0.1, 0.2] .* ig.fovs
    tests = ((:cyl3, [shift' sizes[1] sizes[3] 2]),
             (:gauss3, [shift' sizes' 2]),
             (:rect3, [shift' sizes' 2]),
             (:dirac3,[shift' 2]))
    ntest = length(tests)
    pl = Matrix{Plots.Plot}(undef, ntest,3)
    for (i,test) in enumerate(tests)
        st = mri_objects(test)
        i3 = st.image(ig.xg,ig.yg,ig.zg, dx=ig.dx, dy=ig.dy, dz=ig.dz)
        s3 = ifftshift(fft(fftshift(i3)))
        s3 *= abs(ig.dx * ig.dy * ig.dz)
        fg = ig.fg
        f3 = st.kspace(fg...)
        isinteractive() && (@show max_percent_diff(f3,s3))
        pl[i,1] = jim(i3, title = "$(test[1])")
        pl[i,2] = jim(abs.(s3), title = "fft")
        pl[i,3] = jim(abs.(f3), title = "kspace")
        plot(pl[i,:]...)
        prompt()
    end
    return plot(pl...)
end

# plot the trapezoid
function mri_objects_trap_test()
    len = 6
    dz = 2
    z = LinRange(-8,8,201)
    trap = mri_objects_trap(z, dz, len)
    trap0 = mri_objects_trap(z, 0, len)
    plot(z, trap, label="trap")
    plot!(z, trap0, label="rect")
    plot!(xtick = [0 -len/2 len/2-dz/2 len/2 len/2+dz/2])
end

@test_throws String mri_objects(:bad)

@test mri_objects_trap_test() isa Plots.Plot
prompt()

@test mri_objects_test2() isa Plots.Plot
prompt()

@test mri_objects_test_case1() isa Plots.Plot
prompt()

@test mri_objects_test3() isa Plots.Plot

@test mri_objects_test_case4() isa Plots.Plot
prompt()


text:
module IntroCar

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Intro Car" => Ahorn.EntityPlacement(
        Maple.IntroCar
    )
)

barrierSprite = "scenery/car/barrier"
bodySprite = "scenery/car/body"
pavementSprite = "scenery/car/pavement"
wheelsSprite = "scenery/car/wheels"

function Ahorn.selection(entity::Maple.IntroCar)
    x, y = Ahorn.position(entity)

    hasRoadAndBarriers = get(entity.data, "hasRoadAndBarriers", false)

    rectangles = Ahorn.Rectangle[
        Ahorn.getSpriteRectangle(bodySprite, x, y, jx=0.5, jy=1.0),
        Ahorn.getSpriteRectangle(wheelsSprite, x, y, jx=0.5, jy=1.0),
    ]

    if hasRoadAndBarriers
        push!(rectangles, Ahorn.getSpriteRectangle(barrierSprite, x + 32, y, jx=0.0, jy=1.0))
        push!(rectangles, Ahorn.getSpriteRectangle(barrierSprite, x + 41, y, jx=0.0, jy=1.0))
    end

    return Ahorn.coverRectangles(rectangles)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.IntroCar, room::Maple.Room)
    x, y = Ahorn.position(entity)

    hasRoadAndBarriers = get(entity.data, "hasRoadAndBarriers", false)
    rng = Ahorn.getSimpleEntityRng(entity)

    pavementWidth = x - 48
    columns = floor(Int, pavementWidth / 8)

    Ahorn.drawSprite(ctx, wheelsSprite, x, y, jx=0.5, jy=1.0)
    Ahorn.drawSprite(ctx, bodySprite, x, y, jx=0.5, jy=1.0)

    if hasRoadAndBarriers
        Ahorn.drawSprite(ctx, barrierSprite, x + 32, y, jx=0.0, jy=1.0)
        Ahorn.drawSprite(ctx, barrierSprite, x + 41, y, jx=0.0, jy=1.0)

        for col in 0:columns - 1
            choice = col >= columns - 2 ? (col != columns - 2 ? 3 : 2) : rand(rng, 0:2)

            Ahorn.drawImage(ctx, pavementSprite, col * 8, y, choice * 8, 0, 8, 4)
        end
    end
end

end

ext: jl


text:
@testset "CommonRLInterface" begin
@testset "MDPEnv" begin
    struct RLTestMDP <: MDP{Int, Int} end

    POMDPs.actions(m::RLTestMDP) = [-1, 1]
    POMDPs.transition(m::RLTestMDP, s, a) = Deterministic(clamp(s + a, 1, 3))
    POMDPs.initialstate(m::RLTestMDP) = Deterministic(1)
    POMDPs.isterminal(m::RLTestMDP, s) = s == 3
    POMDPs.reward(m::RLTestMDP, s, a, sp) = sp
    POMDPs.states(m::RLTestMDP) = 1:3

    env = convert(RLBase.AbstractEnv, convert(CRL.AbstractEnv, RLTestMDP()))
    RLBase.test_runnable!(env)
end

@testset "POMDPEnv" begin

    struct RLTestPOMDP <: POMDP{Int, Int, Int} end

    POMDPs.actions(m::RLTestPOMDP) = [-1, 1]
    POMDPs.states(m::RLTestPOMDP) = 1:3
    POMDPs.transition(m::RLTestPOMDP, s, a) = Deterministic(clamp(s + a, 1, 3))
    POMDPs.observation(m::RLTestPOMDP, s, a, sp) = Deterministic(sp + 1)
    POMDPs.initialstate(m::RLTestPOMDP) = Deterministic(1)
    POMDPs.initialobs(m::RLTestPOMDP, s) = Deterministic(s + 1)
    POMDPs.isterminal(m::RLTestPOMDP, s) = s == 3
    POMDPs.reward(m::RLTestPOMDP, s, a, sp) = sp
    POMDPs.observations(m::RLTestPOMDP) = 2:4

    env = convert(RLBase.AbstractEnv, convert(CRL.AbstractEnv, RLTestPOMDP()))

    RLBase.test_runnable!(env)
end
end


text:
# ------------------------------------------------------------------------------------------
# # Sequence to sequence model (S2S)
# (c) Deniz Yuret, 2018. Based on ([Sutskever et al.
# 2014](https://papers.nips.cc/paper/5346-sequence-to-sequence-learning-with-neural-
# networks.pdf)).
# S2S models learn to map input sequences to output sequences using an encoder and a decoder
# RNN. Note that this is an instructional example written in low-level Julia/Knet and it is
# slow to train. For a faster and high-level implementation please see `@doc RNN`.
# ------------------------------------------------------------------------------------------

using Pkg; haskey(Pkg.installed(),"Knet") || Pkg.add("Knet")
using Knet

# ------------------------------------------------------------------------------------------
# (image source)
# ------------------------------------------------------------------------------------------

# S2S model definition

function initmodel(H, V; atype=(gpu()>=0 ? KnetArray{Float32} : Array{Float32}))
    init(d...)=atype(xavier(d...))
    model = Dict{Symbol,Any}()
    model[:state0] = [ init(1,H), init(1,H) ]
    model[:embed1] = init(V,H)
    model[:encode] = [ init(2H,4H), init(1,4H) ]
    model[:embed2] = init(V,H)
    model[:decode] = [ init(2H,4H), init(1,4H) ]
    model[:output] = [ init(H,V), init(1,V) ]
    return model
end;

# S2S loss function and its gradient

function s2s(model, inputs, outputs)
    state = initstate(inputs[1], model[:state0])
    for input in inputs
        input = onehotrows(input, model[:embed1])
        input = input * model[:embed1]
        state = lstm(model[:encode], state, input)
    end
    EOS = eosmatrix(outputs[1], model[:embed2])
    input = EOS * model[:embed2]
    sumlogp = 0
    for output in outputs
        state = lstm(model[:decode], state, input)
        ypred = predict(model[:output], state[1])
        ygold = onehotrows(output, model[:embed2])
        sumlogp += sum(ygold .* logp(ypred,dims=2))
        input = ygold * model[:embed2]
    end
    state = lstm(model[:decode], state, input)
    ypred = predict(model[:output], state[1])
    sumlogp += sum(EOS .* logp(ypred,dims=2))
    return -sumlogp
end

s2sgrad = gradloss(s2s);

# ------------------------------------------------------------------------------------------
# (image source)
# ------------------------------------------------------------------------------------------

# A LSTM implementation in Knet

function lstm(param, state, input)
    weight,bias = param
    hidden,cell = state
    h       = size(hidden,2)
    gates   = hcat(input,hidden) * weight .+ bias
    forget  = sigm.(gates[:,1:h])
    ingate  = sigm.(gates[:,1+h:2h])
    outgate = sigm.(gates[:,1+2h:3h])
    change  = tanh.(gates[:,1+3h:4h])
    cell    = cell .* forget + ingate .* change
    hidden  = outgate .* tanh.(cell)
    return (hidden,cell)
end;

# S2S helper functions

function predict(param, input)
    input * param[1] .+ param[2]
end

function initstate(idx, state0)
    h,c = state0
    h = h .+ fill!(similar(value(h), length(idx), length(h)), 0)
    c = c .+ fill!(similar(value(c), length(idx), length(c)), 0)
    return (h,c)
end

function onehotrows(idx, embeddings)
    nrows,ncols = length(idx), size(embeddings,1)
    z = zeros(Float32,nrows,ncols)
    @inbounds for i=1:nrows
        z[i,idx[i]] = 1
    end
    oftype(value(embeddings),z)
end

let EOS=nothing; global eosmatrix
function eosmatrix(idx, embeddings)
    nrows,ncols = length(idx), size(embeddings,1)
    if EOS==nothing || size(EOS) != (nrows,ncols)
        EOS = zeros(Float32,nrows,ncols)
        EOS[:,1] .= 1
        EOS = oftype(value(embeddings), EOS)
    end
    return EOS
end
end;

# Use reversing English words as an example task
# This loads them from /usr/share/dict/words and converts each character to an int.

function readdata(file="words")
    isfile(file) || (file=download("http://people.csail.mit.edu/deniz/models/tutorial/words","words"))
    global strings = map(chomp,readlines(file))
    global tok2int = Dict{Char,Int}()
    global int2tok = Vector{Char}()
    push!(int2tok,'\n'); tok2int['\n']=1 # We use '\n'=>1 as the EOS token
    sequences = Vector{Vector{Int}}()
    for w in strings
        s = Vector{Int}()
        for c in collect(w)
            if !haskey(tok2int,c)
                push!(int2tok,c)
                tok2int[c] = length(int2tok)
            end
            push!(s, tok2int[c])
        end
        push!(sequences, s)
    end
    return sequences
end;

sequences = readdata();
for x in (sequences, strings, int2tok, tok2int); println(summary(x)); end
for x in strings[501:505]; println(x); end

# Minibatch sequences putting equal length sequences together:

function minibatch(sequences, batchsize)
    table = Dict{Int,Vector{Vector{Int}}}()
    data = Any[]
    for s in sequences
        n = length(s)
        nsequences = get!(table, n, Any[])
        push!(nsequences, s)
        if length(nsequences) == batchsize
            push!(data, [[ nsequences[i][j] for i in 1:batchsize] for j in 1:n ])
            empty!(nsequences)
        end
    end
    return data
end

batchsize, statesize, vocabsize = 128, 128, length(int2tok)
data = minibatch(sequences,batchsize)
summary(data)

# Training loop

function train(model, data, opts)
    sumloss = cntloss = 0
    for sequence in data
        grads,loss = s2sgrad(model, sequence, reverse(sequence))
        update!(model, grads, opts)
        sumloss += loss
        cntloss += (1+length(sequence)) * length(sequence[1])
    end
    return sumloss/cntloss
end

file = "rnnreverse113.jld2"; model = opts = nothing; Knet.gc() # clean memory from previous run
if (print("Train from scratch? ");readline()[1]=='y')
    # Initialize model and optimization parameters
    model = initmodel(statesize,vocabsize)
    opts = optimizers(model,Adam)
    @time for epoch=1:10
        @time loss = train(model,data,opts) # ~17 sec/epoch
        println((epoch,loss))
    end
    Knet.save(file,"model",model)
else
    isfile(file) || download("http://people.csail.mit.edu/deniz/models/tutorial/$file",file)
    model = Knet.load(file,"model")
end
summary(model)

# Test on some examples:

function translate(model, str)
    state = model[:state0]
    for c in collect(str)
        input = onehotrows(tok2int[c], model[:embed1])
        input = input * model[:embed1]
        state = lstm(model[:encode], state, input)
    end
    input = eosmatrix(1, model[:embed2]) * model[:embed2]
    output = Char[]
    for i=1:100 #while true
        state = lstm(model[:decode], state, input)
        pred = predict(model[:output], state[1])
        i = argmax(vec(Array(pred)))
        i == 1 && break
        push!(output, int2tok[i])
        input = onehotrows(i, model[:embed2]) * model[:embed2]
    end
    String(output)
end;

translate(model,"capricorn")

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
module RangeEnclosures

using DynamicPolynomials, Requires, Reexport
@reexport using ForwardDiff
@reexport using IntervalArithmetic
const Interval_or_IntervalBox = Union{Interval, IntervalBox}

include("algorithms.jl")
include("intervals.jl")
include("branchandbound.jl")

#================
Optional methods
=================#
function __init__()
    @require SumOfSquares = "4b9e565b-77fc-50a5-a571-1244f986bda1" include("sdp.jl")
    @require TaylorModels = "314ce334-5f6e-57ae-acf6-00b6e903104a" include("taylormodels.jl")
    @require IntervalOptimisation = "c7c68f13-a4a2-5b9a-b424-07d005f8d9d2" include("intervaloptimisation.jl")
end

#================
API
=================#
include("enclose.jl")

export enclose, relative_precision,
  NaturalEnclosure, MooreSkelboeEnclosure, SumOfSquaresEnclosure,
  TaylorModelsEnclosure, BranchAndBoundEnclosure,
  AbstractEnclosureAlgorithm

end # module


text:
import Base.Broadcast: BroadcastStyle, AbstractArrayStyle, Broadcasted, DefaultArrayStyle, materialize!
import Base.Broadcast: _bcs1  # for SOneTo axis information
using Base.Broadcast: _bcsm
# Add a new BroadcastStyle for TupleVectors, derived from AbstractArrayStyle
# A constructor that changes the style parameter N (array dimension) is also required
struct TupleVectorStyle{N} <: AbstractArrayStyle{N} end
TupleVectorStyle{M}(::Val{N}) where {M,N} = TupleVectorStyle{N}()
BroadcastStyle(::Type{<:TupleVector{N,<:Any}}) where N = TupleVectorStyle{N}()
BroadcastStyle(::Type{<:LinearAlgebra.Transpose{<:Any,<:TupleVector{N,<:Any}}}) where N = TupleVectorStyle{N}()
BroadcastStyle(::Type{<:LinearAlgebra.Adjoint{<:Any,<:TupleVector{N,<:Any}}}) where N = TupleVectorStyle{N}()
# Precedence rules
BroadcastStyle(::TupleVectorStyle{M}, ::DefaultArrayStyle{N}) where {M,N} =
    DefaultArrayStyle(Val(max(M, N)))
BroadcastStyle(::TupleVectorStyle{M}, ::DefaultArrayStyle{0}) where {M} =
    TupleVectorStyle{M}()
# copy overload
@inline function Base.copy(B::Broadcasted{TupleVectorStyle{M}}) where M
    flat = Broadcast.flatten(B); as = flat.args; f = flat.f
    argsizes = broadcast_sizes(as...)
    destsize = combine_sizes(argsizes)
    _broadcast(f, destsize, argsizes, as...)
end
# copyto! overloads
@inline Base.copyto!(dest, B::Broadcasted{<:TupleVectorStyle}) = _copyto!(dest, B)
@inline Base.copyto!(dest::AbstractArray, B::Broadcasted{<:TupleVectorStyle}) = _copyto!(dest, B)
@inline function _copyto!(dest, B::Broadcasted{TupleVectorStyle{M}}) where M
    flat = Broadcast.flatten(B); as = flat.args; f = flat.f
    argsizes = broadcast_sizes(as...)
    destsize = combine_sizes((Val(M), argsizes...))
    #=if Val(destsize) === Val{Dynamic()}()
        # destination dimension cannot be determined statically; fall back to generic broadcast!
        return copyto!(dest, convert(Broadcasted{DefaultArrayStyle{M}}, B))
    end=#
    _broadcast!(f, destsize, dest, argsizes, as...)
end

# Resolving priority between dynamic and static axes
_bcs1(a::SOneTo, b::SOneTo) = _bcsm(b, a) ? b : (_bcsm(a, b) ? a : throw(DimensionMismatch("arrays could not be broadcast to a common size")))
_bcs1(a::SOneTo, b::Base.OneTo) = _bcs1(Base.OneTo(a), b)
_bcs1(a::Base.OneTo, b::SOneTo) = _bcs1(a, Base.OneTo(b))

###################################################
## Internal broadcast machinery for TupleVectors ##
###################################################

broadcast_indices(A::TupleVector) = indices(A)

# TODO: just use map(broadcast_size, as)?
@inline broadcast_sizes(a, as...) = (broadcast_size(a), broadcast_sizes(as...)...)
@inline broadcast_sizes() = ()
@inline broadcast_size(a) = Val(0)
@inline broadcast_size(a::AbstractVector) = Val(length(a))
@inline broadcast_size(a::NTuple{N}) where N = Val(N)

function broadcasted_index(oldsize, newindex)
    index = ones(Int, length(oldsize))
    for i = 1:length(oldsize)
        if oldsize[i] != 1
            index[i] = newindex[i]
        end
    end
    return LinearIndices((oldsize,))[index...]
end

# similar to Base.Broadcast.combine_indices:
@generated function combine_sizes(s::Tuple{Vararg{Val}})
    sizes = [sz.parameters[1] for sz ∈ s.parameters]
    ndims = 0
    for i = 1:length(sizes)
        ndims = max(ndims,sizes[i])
    end
    quote
        @_inline_meta
        Val($ndims)
    end
end

scalar_getindex(x) = x
scalar_getindex(x::Ref) = x[]

@generated function _broadcast(f, ::Val{newsize}, s::Tuple{Vararg{Val}}, a...) where newsize
    first_staticarray = a[findfirst(ai -> ai <: Union{TupleVector, LinearAlgebra.Transpose{<:Any, <:TupleVector}, LinearAlgebra.Adjoint{<:Any, <:TupleVector}}, a)]
    if newsize == 0
        # Use inference to get eltype in empty case (see also comments in _map)
        eltys = [:(eltype(a[$i])) for i ∈ 1:length(a)]
        return quote
            @_inline_meta
            T = Core.Compiler.return_type(f, Tuple{$(eltys...)})
            @inbounds return similar_type($first_staticarray, T, Val(newsize))()
        end
    end
    sizes = [sz.parameters[1] for sz ∈ s.parameters]
    indices = CartesianIndices((newsize,))
    exprs = similar(indices, Expr)
    for (j, current_ind) ∈ enumerate(indices)
        exprs_vals = [
            (!(a[i] <: AbstractArray || a[i] <: Tuple) ? :(scalar_getindex(a[$i])) : :(a[$i][$(broadcasted_index(sizes[i], current_ind))]))
            for i = 1:length(sizes)
        ]
        exprs[j] = :(f($(exprs_vals...)))
    end
    return quote
        @_inline_meta
        @inbounds elements = tuple($(exprs...))
        @inbounds return similar_type($first_staticarray, eltype(elements), Val(newsize))(elements)
    end
end

####################################################
## Internal broadcast! machinery for TupleVectors ##
####################################################

@generated function _broadcast!(f, ::Val{newsize}, dest::AbstractArray, s::Tuple{Vararg{Val}}, as...) where {newsize}
    sizes = [sz.parameters[1] for sz ∈ s.parameters]
    sizes = tuple(sizes...)
    indices = CartesianIndices((newsize,))
    exprs = similar(indices, Expr)
    for (j, current_ind) ∈ enumerate(indices)
        exprs_vals = [
            (!(as[i] <: AbstractArray || as[i] <: Tuple) ? :(as[$i][]) : :(as[$i][$(broadcasted_index(sizes[i], current_ind))]))
            for i = 1:length(sizes)
        ]
        exprs[j] = :(dest[$j] = f($(exprs_vals...)))
    end
    return quote
        Base.@_propagate_inbounds_meta
        @inbounds $(Expr(:block, exprs...))
        return dest
    end
end


text:
export PartialInverseOperator

struct PartialInverseOperator{T<:Number,CO<:CachedOperator,BI} <: Operator{T}
    cache::CO
    bandwidths::BI
end

function PartialInverseOperator(CO::CachedOperator{T},bandwidths) where T<:Number
    @assert istriu(CO) # || istril(CO)
    return PartialInverseOperator{T,typeof(CO),typeof(bandwidths)}(CO,bandwidths)
end

PartialInverseOperator(B::Operator, bandwidths) = PartialInverseOperator(cache(B), bandwidths)
PartialInverseOperator(B::Operator) = PartialInverseOperator(B, bandwidths(B))

convert(::Type{Operator{T}},A::PartialInverseOperator) where {T}=PartialInverseOperator(convert(Operator{T},A.cache), A.bandwidths)

domainspace(P::PartialInverseOperator)=rangespace(P.cache)
rangespace(P::PartialInverseOperator)=domainspace(P.cache)
domain(P::PartialInverseOperator)=domain(domainspace(P))
bandwidths(P::PartialInverseOperator) = P.bandwidths

function getindex(P::PartialInverseOperator,k::Integer,j::Integer)
    b = bandwidth(P.cache, 2)
    if k == j
        return inv(P.cache[k,k])
    elseif j > k
        t = zero(T)
        for i = max(k,j-b-1):j-1
            t += ret[k,i]*P.cache[i,j]
        end
        return -t/P.cache[j,j]
    else
        return zero(eltype(P))
    end
end

## These are both hacks that apparently work

function BandedMatrix(S::SubOperator{T,PP,Tuple{UnitRange{Int},UnitRange{Int}}}) where {T,PP<:PartialInverseOperator}
    kr,jr = parentindices(S)
    P = parent(S)
    #ret = BandedMatrix{eltype(S)}(undef, size(S), bandwidths(S))
    ret = BandedMatrix{eltype(S)}(undef, (last(kr),last(jr)), bandwidths(P))
    b = bandwidth(P, 2)
    #@assert first(kr) == first(jr) == 1

    @inbounds for j in 1:last(jr)
        kk = colrange(ret, j)
        if j in kk
            ret[j,j] = inv(P.cache[j,j])
        end
        for k in first(kk):min(last(kk),j-1)
            t = zero(T)
            for i = max(k,j-b-1):j-1
                t += ret[k,i]*P.cache[i,j]
            end
            ret[k,j] = -t/P.cache[j,j]
        end
    end

    ret[kr,jr]
end


text:
using RootFindGF
using Base.Test
using Nemo

function test_all()
    test_Γ()
    test_multi_eval()
    test_interpolation()
    test_resultant()
    test_SRA()
    test_flint()
end

K, x = FiniteField(2, 10, "x")
P, X = PolynomialRing(K, "X")

function test_Γ()
    Γ(map(i -> x^i, 0:9))
end

function test_multi_eval()
    list = map(i -> x^i, 0:400)
    f = X^400 + X^30 + x
    @test RootFindGF.multi_eval(f, list, 10) == [f(p) for p in list]
end

function test_interpolation()
    xs = map(i -> x^i, 0:400)
    f = X^400 + X^30 + x
    ys = map(f, xs)
    @test RootFindGF.interpolate(P, xs, ys) == f
end

function test_resultant()
    PP, Y = PolynomialRing(P, "Y")
    f = X^400 + X^30 + x
    @test RootFindGF.resultant(f, x) == resultant(compose(f, Y), X - Y^2 + x*Y)
end

function test_SRA()
    roots = [x^i for i in 0:300]
    pol = prod([X - r for r in roots])
    @test Set(roots) == Set(SRA(pol))
end

function test_flint()
    roots = [x^i for i in 0:300]
    pol = prod([X - r for r in roots])
    @test Set(roots) == Set(keys(flint_root_find(pol)))
end

test_all()


text:
using Base.Test

allPassed = true
try
	include("sensitivityTest.jl");
catch
	allPassed = false
	warn("including sensitivityTest.jl failed.")
end

try
	runSensitivityTest();
catch
	allPassed = false
	warn("runSensitivityTest() had test errors")
end

@test allPassed == true


text:
module Roots

if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@optlevel"))
    @eval Base.Experimental.@optlevel 1
end

using Printf
import CommonSolve
import CommonSolve: solve, solve!, init
using Setfield

export fzero,
       fzeros,
       secant_method

export find_zero, find_zero!, find_zeros,
       ZeroProblem, solve, solve!, init,
       Order0, Order1, Order2, Order5, Order8, Order16

export Bisection, FalsePosition

## load in files
include("utils.jl")
include("find_zero.jl")
include("bracketing.jl")
include("derivative_free.jl")
include("order0.jl")
include("simple.jl")
include("find_zeros.jl")
include("newton.jl")
include("lith.jl")
include("fzero.jl")

# cf. https://github.com/JuliaDocs/Documenter.jl/pull/1664/files
function _update_module_doc()
    path = joinpath(@__DIR__, "..", "README.md")
    text = read(path, String)
    # The code blocks in the README.md should be julia blocks for the syntax highlighter.
    text = replace(text, "```julia" => "```jldoctest readme")
    @doc text Roots
end
_update_module_doc()

end


text:
@test   TS.write_comment_line( "" ) == "!"
@test   TS.write_comment_line( "comment" ) == "!comment"

@test   TS.write_option_line( Options( 1e9, :ScatteringParameters, :MagnitudeAngle, 50   ) ) == "# GHz S MA R 50.0"
@test   TS.write_option_line( Options( 1e3, :ScatteringParameters, :MagnitudeAngle, 50.0 ) ) == "# kHz S MA R 50.0"
@test   TS.write_option_line( Options( 1e6, :AdmittanceParameters, :MagnitudeAngle, 50.0 ) ) == "# MHz Y MA R 50.0"
@test   TS.write_option_line( Options( 1e9, :ImpedanceParameters,  :DecibelAngle,   50.0 ) ) == "# GHz Z dB R 50.0"
@test   TS.write_option_line( Options( 1.0, :HybridHParameters,    :RealImaginary,   1.0 ) ) == "# Hz H RI R 1.0"
@test   TS.write_option_line( Options( 1.0, :HybridGParameters,    :RealImaginary,  50.0 ) ) == "# Hz G RI R 50.0"

@test TS.write_data(
  DataPoint( 1.0, fill( 2.0 + 3im, 1, 1 ) ), 1,
  Options( 1.0, :ScatteringParameters, :RealImaginary, 1.0 )
  ) == "1.0 2.0 3.0\n"

@test TS.write_data(
  DataPoint( 1e9, fill( cis( π / 2 ), 1, 1 ) ), 1,
  Options( 1e9, :ScatteringParameters, :MagnitudeAngle, 1.0 )
  ) == "1.0 1.0 90.0\n"

@test TS.write_data(
  DataPoint( 1e9, fill( 10cis( π / 2 ), 1, 1 ) ), 1,
  Options( 1e9, :ScatteringParameters, :DecibelAngle, 1.0 )
  ) == """
  1.0 20.0 90.0
  """

@test TS.write_data(
  DataPoint( 1.0, [ 2.0 3.0; 4.0 5.0 ] ), 2,
  Options( 1.0, :ScatteringParameters, :RealImaginary, 1.0 )
  ) == """
  1.0 2.0 0.0 4.0 0.0 3.0 0.0 5.0 0.0
  """

@test TS.write_data(
  DataPoint( 1.0, [ 2.0 3.0 4.0; 5.0 6.0 7.0; 8.0 9.0 10.0 ] ), 3,
  Options( 1.0, :ScatteringParameters, :RealImaginary, 1.0 )
  ) == """
  1.0 2.0 0.0 3.0 0.0 4.0 0.0
  5.0 0.0 6.0 0.0 7.0 0.0
  8.0 0.0 9.0 0.0 10.0 0.0
  """

@test TS.write_data(
  DataPoint( 1.0, [ 2.0 3.0 4.0 5.0; 6.0 7.0 8.0 9.0; 10.0 11.0 12.0 13.0; 14.0 15.0 16.0 17.0 ] ), 4,
  Options( 1.0, :ScatteringParameters, :RealImaginary, 1.0 )
  ) == """
  1.0 2.0 0.0 3.0 0.0 4.0 0.0 5.0 0.0
  6.0 0.0 7.0 0.0 8.0 0.0 9.0 0.0
  10.0 0.0 11.0 0.0 12.0 0.0 13.0 0.0
  14.0 0.0 15.0 0.0 16.0 0.0 17.0 0.0
  """

@test TS.write_data(
  DataPoint( 1.0, [
    4.0+6.0im   12.0+14.0im;
    8.0+10.0im  16.0+18.0im;
  ] ),
  2,
  Options( 1.0, :ImpedanceParameters, :RealImaginary, 2.0 )
  ) == """
  1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0
  """

@test TS.write_data(
  DataPoint( 1.0, [
    4.0+6.0im   12.0+14.0im;
    8.0+10.0im  16.0+18.0im;
  ] ),
  2,
  Options( 1.0, :AdmittanceParameters, :RealImaginary, 0.5 )
  ) == """
  1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0
  """

@test TS.write_data(
  DataPoint( 1.0, [
    1          6.0+7.0im;
    4.0+5.0im  1;
  ] ),
  2,
  Options( 1.0, :HybridGParameters, :RealImaginary, 2 )
  ) == """
  1.0 2.0 0.0 4.0 5.0 6.0 7.0 0.5 0.0
  """

@test TS.write_data(
  DataPoint( 1.0, [
    1          6.0+7.0im;
    4.0+5.0im  1;
  ] ),
  2,
  Options( 1.0, :HybridHParameters, :RealImaginary, 2 )
  ) == """
  1.0 0.5 0.0 4.0 5.0 6.0 7.0 2.0 0.0
  """

@test TS.writeNoiseData( NoiseDataPoint( 1e9, 10, im, 50 ) ) == "1.0 20.0 1.0 90.0 1.0\n"

@test TS.writeNoiseData(
  NoiseDataPoint( 1, 10, im, 1 ),
  Options( 1.0, :ScatteringParameters, :RealImaginary, 1.0 )
  ) == "1.0 20.0 1.0 90.0 1.0\n"

@test write_touchstone_string(
  TouchstoneData(
    [ DataPoint( 100e6, fill( 0.99cis(   -4π / 180 ), 1, 1 ) ) ],
    [
      NoiseDataPoint( 100e6, 10, im, 50 ),
      NoiseDataPoint( 200e6, 10, im, 50 ),
    ],
    Options( 1e6, :ScatteringParameters, :MagnitudeAngle, 50.0 ),
    [],
  )
  ) == """
  # MHz S MA R 50.0
  100.0 0.9899999999999999 -4.0
  100.0 20.0 1.0 90.0 1.0
  200.0 20.0 1.0 90.0 1.0
  """

@test write_touchstone_string( TouchstoneData( DataPoint[] ) ) == """
  # GHz S MA R 50.0
  """

@test write_touchstone_string( TouchstoneData( DataPoint[], NoiseDataPoint[], Options(), [ "Test" ] ) ) == """
  !Test
  # GHz S MA R 50.0
  """

@test write_touchstone_string(
  TouchstoneData(
    [
      DataPoint( 100e6, fill( 0.99cis(   -4π / 180 ), 1, 1 ) ),
      DataPoint( 200e6, fill( 0.80cis(  -22π / 180 ), 1, 1 ) ),
      DataPoint( 300e6, fill( 0.707cis( -45π / 180 ), 1, 1 ) ),
      DataPoint( 400e6, fill( 0.40cis(  -62π / 180 ), 1, 1 ) ),
      DataPoint( 500e6, fill( 0.01cis(  -89π / 180 ), 1, 1 ) ),
    ],
    NoiseDataPoint[],
    Options( 1e6, :ScatteringParameters, :MagnitudeAngle, 75.0 ),
    [ "Test1", "freq magZ11 angZ11" ]
  )
  ) == """
  !Test1
  !freq magZ11 angZ11
  # MHz S MA R 75.0
  100.0 0.9899999999999999 -4.0
  200.0 0.8000000000000002 -21.999999999999996
  300.0 0.707 -44.99999999999999
  400.0 0.4000000000000001 -62.0
  500.0 0.010000000000000002 -89.00000000000001
  """

@test write_touchstone_string(
  TouchstoneData(
    [
      DataPoint( 1e9,  [
        0.3926-0.1211im -0.0003-0.0021im;
        -0.0003-0.0021im 0.3926-0.1211im
      ] ),
      DataPoint( 2e9,  [
        0.3517-0.3054im -0.0096-0.0298im;
        -0.0096-0.0298im 0.3517-0.3054im
      ] ),
      DataPoint( 10e9, [
        0.3419+0.3336im -0.0134+0.0379im;
        -0.0134+0.0379im 0.3419+0.3336im
      ] ),
    ],
    NoiseDataPoint[],
    Options( 1e9, :ScatteringParameters, :RealImaginary, 50.0 ),
    [
      "2-port S-parameter file, three frequency points",
      "freq   ReS11   ImS11   ReS21   ImS21   ReS12   ImS12   ReS22   ImS22",
    ],
  ) ) == """
  !2-port S-parameter file, three frequency points
  !freq   ReS11   ImS11   ReS21   ImS21   ReS12   ImS12   ReS22   ImS22
  # GHz S RI R 50.0
  1.0 0.3926 -0.1211 -0.0003 -0.0021 -0.0003 -0.0021 0.3926 -0.1211
  2.0 0.3517 -0.3054 -0.0096 -0.0298 -0.0096 -0.0298 0.3517 -0.3054
  10.0 0.3419 0.3336 -0.0134 0.0379 -0.0134 0.0379 0.3419 0.3336
  """

@test write_touchstone_string(
  TouchstoneData(
    [
      DataPoint( 1e9,  [
        1.11+1.11im 1.12+1.12im;
        1.21+1.21im 1.22+1.22im
      ] ),
      DataPoint( 2e9,  [
        2.11+2.11im 2.12+2.12im;
        2.21+2.21im 2.22+2.22im
      ] ),
    ],
    NoiseDataPoint[],
    Options( 1e9, :ScatteringParameters, :RealImaginary, 50.0 )
  )
  ) == """
  # GHz S RI R 50.0
  1.0 1.11 1.11 1.21 1.21 1.12 1.12 1.22 1.22
  2.0 2.11 2.11 2.21 2.21 2.12 2.12 2.22 2.22
  """

@test write_touchstone_string(
  TouchstoneData(
    [
      DataPoint( 1e9,  [
        1.11+1.11im 1.12+1.12im 1.13+1.13im;
        1.21+1.21im 1.22+1.22im 1.23+1.23im;
        1.31+1.31im 1.32+1.32im 1.33+1.33im;
      ] ),
      DataPoint( 2e9,  [
        2.11+2.11im 2.12+2.12im 2.13+2.13im;
        2.21+2.21im 2.22+2.22im 2.23+2.23im;
        2.31+2.31im 2.32+2.32im 2.33+2.33im;
      ] ),
    ],
    NoiseDataPoint[],
    Options( 1e9, :ScatteringParameters, :RealImaginary, 50.0 )
  )
  ) == """
  # GHz S RI R 50.0
  1.0 1.11 1.11 1.12 1.12 1.13 1.13
  1.21 1.21 1.22 1.22 1.23 1.23
  1.31 1.31 1.32 1.32 1.33 1.33
  2.0 2.11 2.11 2.12 2.12 2.13 2.13
  2.21 2.21 2.22 2.22 2.23 2.23
  2.31 2.31 2.32 2.32 2.33 2.33
  """

@test write_touchstone_string(
  TouchstoneData(
    [
      DataPoint( 1e9,  [
        1.11+1.11im 1.12+1.12im 1.13+1.13im 1.14+1.14im;
        1.21+1.21im 1.22+1.22im 1.23+1.23im 1.24+1.24im;
        1.31+1.31im 1.32+1.32im 1.33+1.33im 1.34+1.34im;
        1.41+1.41im 1.42+1.42im 1.43+1.43im 1.44+1.44im;
      ] ),
      DataPoint( 2e9,  [
        2.11+2.11im 2.12+2.12im 2.13+2.13im 2.14+2.14im;
        2.21+2.21im 2.22+2.22im 2.23+2.23im 2.24+2.24im;
        2.31+2.31im 2.32+2.32im 2.33+2.33im 2.34+2.34im;
        2.41+2.41im 2.42+2.42im 2.43+2.43im 2.44+2.44im;
      ] ),
    ],
    NoiseDataPoint[],
    Options( 1e9, :ScatteringParameters, :RealImaginary, 50.0 )
  )
  ) == """
  # GHz S RI R 50.0
  1.0 1.11 1.11 1.12 1.12 1.13 1.13 1.14 1.14
  1.21 1.21 1.22 1.22 1.23 1.23 1.24 1.24
  1.31 1.31 1.32 1.32 1.33 1.33 1.34 1.34
  1.41 1.41 1.42 1.42 1.43 1.43 1.44 1.44
  2.0 2.11 2.11 2.12 2.12 2.13 2.13 2.14 2.14
  2.21 2.21 2.22 2.22 2.23 2.23 2.24 2.24
  2.31 2.31 2.32 2.32 2.33 2.33 2.34 2.34
  2.41 2.41 2.42 2.42 2.43 2.43 2.44 2.44
  """

@test write_touchstone_string(
  TouchstoneData(
    [
      DataPoint( 1e9,  [
        1.11+1.11im 1.12+1.12im 1.13+1.13im 1.14+1.14im 1.15+1.15im;
        1.21+1.21im 1.22+1.22im 1.23+1.23im 1.24+1.24im 1.25+1.25im;
        1.31+1.31im 1.32+1.32im 1.33+1.33im 1.34+1.34im 1.35+1.35im;
        1.41+1.41im 1.42+1.42im 1.43+1.43im 1.44+1.44im 1.45+1.45im;
        1.51+1.51im 1.52+1.52im 1.53+1.53im 1.54+1.54im 1.55+1.55im;
      ] ),
      DataPoint( 2e9,  [
        2.11+2.11im 2.12+2.12im 2.13+2.13im 2.14+2.14im 2.15+2.15im;
        2.21+2.21im 2.22+2.22im 2.23+2.23im 2.24+2.24im 2.25+2.25im;
        2.31+2.31im 2.32+2.32im 2.33+2.33im 2.34+2.34im 2.35+2.35im;
        2.41+2.41im 2.42+2.42im 2.43+2.43im 2.44+2.44im 2.45+2.45im;
        2.51+2.51im 2.52+2.52im 2.53+2.53im 2.54+2.54im 2.55+2.55im;
      ] ),
    ],
    NoiseDataPoint[],
    Options( 1e9, :ScatteringParameters, :RealImaginary, 50.0 )
  )
  ) == """
  # GHz S RI R 50.0
  1.0 1.11 1.11 1.12 1.12 1.13 1.13 1.14 1.14
  1.15 1.15
  1.21 1.21 1.22 1.22 1.23 1.23 1.24 1.24
  1.25 1.25
  1.31 1.31 1.32 1.32 1.33 1.33 1.34 1.34
  1.35 1.35
  1.41 1.41 1.42 1.42 1.43 1.43 1.44 1.44
  1.45 1.45
  1.51 1.51 1.52 1.52 1.53 1.53 1.54 1.54
  1.55 1.55
  2.0 2.11 2.11 2.12 2.12 2.13 2.13 2.14 2.14
  2.15 2.15
  2.21 2.21 2.22 2.22 2.23 2.23 2.24 2.24
  2.25 2.25
  2.31 2.31 2.32 2.32 2.33 2.33 2.34 2.34
  2.35 2.35
  2.41 2.41 2.42 2.42 2.43 2.43 2.44 2.44
  2.45 2.45
  2.51 2.51 2.52 2.52 2.53 2.53 2.54 2.54
  2.55 2.55
  """

@test write_touchstone_string(
  TouchstoneData(
    [
      DataPoint( 1e9,  [
        1.11+1.11im 1.12+1.12im 1.13+1.13im 1.14+1.14im 1.15+1.15im 1.16+1.16im 1.17+1.17im 1.18+1.18im 1.19+1.19im 1.110+1.110im;
        1.21+1.21im 1.22+1.22im 1.23+1.23im 1.24+1.24im 1.25+1.25im 1.26+1.26im 1.27+1.27im 1.28+1.28im 1.29+1.29im 1.210+1.210im;
        1.31+1.31im 1.32+1.32im 1.33+1.33im 1.34+1.34im 1.35+1.35im 1.36+1.36im 1.37+1.37im 1.38+1.38im 1.39+1.39im 1.310+1.310im;
        1.41+1.41im 1.42+1.42im 1.43+1.43im 1.44+1.44im 1.45+1.45im 1.46+1.46im 1.47+1.47im 1.48+1.48im 1.49+1.49im 1.410+1.410im;
        1.51+1.51im 1.52+1.52im 1.53+1.53im 1.54+1.54im 1.55+1.55im 1.56+1.56im 1.57+1.57im 1.58+1.58im 1.59+1.59im 1.510+1.510im;
        1.61+1.61im 1.62+1.62im 1.63+1.63im 1.64+1.64im 1.65+1.65im 1.66+1.66im 1.67+1.67im 1.68+1.68im 1.69+1.69im 1.610+1.610im;
        1.71+1.71im 1.72+1.72im 1.73+1.73im 1.74+1.74im 1.75+1.75im 1.76+1.76im 1.77+1.77im 1.78+1.78im 1.79+1.79im 1.710+1.710im;
        1.81+1.81im 1.82+1.82im 1.83+1.83im 1.84+1.84im 1.85+1.85im 1.86+1.86im 1.87+1.87im 1.88+1.88im 1.89+1.89im 1.810+1.810im;
        1.91+1.91im 1.92+1.92im 1.93+1.93im 1.94+1.94im 1.95+1.95im 1.96+1.96im 1.97+1.97im 1.98+1.98im 1.99+1.99im 1.910+1.910im;
        1.101+1.101im 1.102+1.102im 1.103+1.103im 1.104+1.104im 1.105+1.105im 1.106+1.106im 1.107+1.107im 1.108+1.108im 1.109+1.109im 1.1010+1.1010im;
      ] ),
    ],
    NoiseDataPoint[],
    Options( 1e9, :ScatteringParameters, :RealImaginary, 50.0 )
  ) ) == """
  # GHz S RI R 50.0
  1.0 1.11 1.11 1.12 1.12 1.13 1.13 1.14 1.14
  1.15 1.15 1.16 1.16 1.17 1.17 1.18 1.18
  1.19 1.19 1.11 1.11
  1.21 1.21 1.22 1.22 1.23 1.23 1.24 1.24
  1.25 1.25 1.26 1.26 1.27 1.27 1.28 1.28
  1.29 1.29 1.21 1.21
  1.31 1.31 1.32 1.32 1.33 1.33 1.34 1.34
  1.35 1.35 1.36 1.36 1.37 1.37 1.38 1.38
  1.39 1.39 1.31 1.31
  1.41 1.41 1.42 1.42 1.43 1.43 1.44 1.44
  1.45 1.45 1.46 1.46 1.47 1.47 1.48 1.48
  1.49 1.49 1.41 1.41
  1.51 1.51 1.52 1.52 1.53 1.53 1.54 1.54
  1.55 1.55 1.56 1.56 1.57 1.57 1.58 1.58
  1.59 1.59 1.51 1.51
  1.61 1.61 1.62 1.62 1.63 1.63 1.64 1.64
  1.65 1.65 1.66 1.66 1.67 1.67 1.68 1.68
  1.69 1.69 1.61 1.61
  1.71 1.71 1.72 1.72 1.73 1.73 1.74 1.74
  1.75 1.75 1.76 1.76 1.77 1.77 1.78 1.78
  1.79 1.79 1.71 1.71
  1.81 1.81 1.82 1.82 1.83 1.83 1.84 1.84
  1.85 1.85 1.86 1.86 1.87 1.87 1.88 1.88
  1.89 1.89 1.81 1.81
  1.91 1.91 1.92 1.92 1.93 1.93 1.94 1.94
  1.95 1.95 1.96 1.96 1.97 1.97 1.98 1.98
  1.99 1.99 1.91 1.91
  1.101 1.101 1.102 1.102 1.103 1.103 1.104 1.104
  1.105 1.105 1.106 1.106 1.107 1.107 1.108 1.108
  1.109 1.109 1.101 1.101
  """


text:
# Enables deserialization of VariableCost. The default implementation can't figure out the
# variable Union.

function JSON2.read(io::IO, ::Type{VariableCost})
    data = JSON2.read(io)
    if data.cost isa Real
        return VariableCost(Float64(data.cost))
    elseif data.cost[1] isa Array
        variable = Vector{Tuple{Float64, Float64}}()
        for array in data.cost
            push!(variable, Tuple{Float64, Float64}(array))
        end
    else
        @assert data.cost isa Tuple || data.cost isa Array
        variable = Tuple{Float64, Float64}(data.cost)
    end

    return VariableCost(variable)
end

const COMPOSED_COMPONENTS = (Area, Bus, LoadZone)

function JSON2.write(io::IO, component::T) where {T <: Component}
    return JSON2.write(io, encode_components_with_uuids(component, COMPOSED_COMPONENTS))
end

function JSON2.write(component::T) where {T <: Component}
    return JSON2.write(encode_components_with_uuids(component, COMPOSED_COMPONENTS))
end

function encode_components_with_uuids(component::T, types_as_uuids) where {T}
    fields = fieldnames(T)
    vals = []

    for name in fields
        val = getfield(component, name)
        if mapreduce(x -> val isa x, |, types_as_uuids)
            push!(vals, IS.get_uuid(val))
        else
            push!(vals, val)
        end
    end

    return NamedTuple{fields}(vals)
end

# Default JSON conversion for Component to handle composed components stored as UUIDs.
# Keep in sync with COMPOSED_COMPONENTS.

function IS.convert_type(
    ::Type{T},
    data::NamedTuple,
    component_cache::Dict,
) where {T <: Component}
    @debug T data
    values = []
    for (fieldname, fieldtype) in zip(fieldnames(T), fieldtypes(T))
        val = getfield(data, fieldname)
        if !isnothing(val) && (
            fieldtype <: Union{Nothing, Area} ||
            fieldtype <: Union{Nothing, LoadZone} ||
            fieldtype <: Bus
        )
            uuid = Base.UUID(val.value)
            component = component_cache[uuid]
            push!(values, component)
        elseif fieldtype <: Component
            # Recurse.
            push!(values, IS.convert_type(fieldtype, val, component_cache))
        else
            obj = IS.convert_type(fieldtype, val)
            push!(values, obj)
        end
    end

    return T(values...)
end


text:
using Index1024
using Test

include("testfunctions.jl")

@testset "Index1024.jl" begin
    #@test nb(0) == 0
    #@test nb(1) == 16
    @test rw_leafnode()
    @test rw_empty_index()
    @test rw_index1p()
    @test rw_index()
    @test tsearch_exists(8) == ktup(8)
    @test tsearch_exists(8; padding=17) == ktup(8)
    @test tsearch_exists(16) == ktup(16)
    @test tsearch_exists(17) == ktup(17)
    @test tsearch_exists(160) == ktup(160)
    @test tsearch_exists(1_000_000) == ktup(1_000_000)
    @test tsearch_fails(16) === nothing
    @test tsearch_fails(17) === nothing
    @test tsearch_fails(160) === nothing
    @test tget_exists(5) == ktup(5)
    @test tget_exists(16) == ktup(16)
    @test tget_exists(17) == ktup(17)
    @test tget_exists(170) == ktup(170)
    @test tget_fails(16) == 0
    @test tget_fails(17) == 0
    @test tget_fails(170) == 0
    @test tmeta(16)
   # @test sort(collect(keys(node_range(egtree(8), 0x8, 0xa)))) == [0x8, 0x9, 0xa]
  #  @test sort(collect(keys(node_range(egtree(1248), 0x8, 0xa)))) == [0x8, 0x9, 0xa]
  #  @test sort(collect(keys(node_range(egtree(164), 0x9f, 0xffff)))) == [0x9f, 0xa0, 0xa1, 0xa2, 0xa3, 0xa4, 0xa5, 0xa6, 0xa7]
end


text:
### Compute mean recursively given previous mean and new sample

Base.mean!(lastmean::RealVector, k::Integer, x::RealVector) = (lastmean[:] = ((k-1)*lastmean+x)/k)

Base.mean(s::VariableNState{Univariate}) = mean(s.value)

Base.mean(s::VariableNState{Multivariate}, i::Integer) = mean(s.value[i, :])

Base.mean(s::VariableNState{Multivariate}, r::AbstractVector=1:s.size) = eltype(s)[mean(s, i) for i in r]


text:
"""
    shiftswap!(a, inew, iold)

Move an entry in vector `iold`  to `inew`, shifting the other arguments.
"""
function shiftswap!(a::AbstractVector, inew, iold)
    if inew != iold
        v = a[iold]
        deleteat!(a, iold)
        insert!(a, inew, v)
    end
    return a
end

"""
    shiftidx!(idx, fid, inew)

For a vector of unique, consecutive indices `idx`, set the index of entry `fid` to the
new index `inew`. Update the other entries, maintaining relative ordering, such that
`idx` is still a vector of unique, consecutive integers. Runs in O(n).
"""
function shiftidx!(idx::AbstractVector{<:Integer}, fid::Integer, inew::Integer)
    iold = idx[fid]
    ilo = min(iold,inew)
    ihi = max(iold,inew)
    s = -sign(inew - iold)
    for j = 1:length(idx)
        if ilo <= idx[j] <= ihi
            idx[j] += s
        end
    end
    idx[fid] = inew
    return idx
end

"""
    getwidth(fmt::String)

Get the width field from a C-style printf format string.
It looks for a set of numbers followed by a decimal or letter, but not preceded by a decimal
or digit. Returns -1 if no width is found.
"""
function getwidth(fmt::String)::Int
    period = findlast(".", fmt)

    # search for a set of characters followed by a decimal or letter, but not
    #   preceded by a decimal or digit
    regex = r"(?<![0-9|.])[0-9]+(?=[a-zA-Z|.])"
    inds = findfirst(regex, fmt)
    if !isnothing(inds)
        width = parse(Int, fmt[inds])
        if width == 0
            return -1
        else
            return width
        end
    else
        return -1
    end
end


text:
using Random

simulation_end = 1000
num_bikes = 3

bike_positions_x = zeros(num_bikes)
bike_positions_y = zeros(num_bikes)

for bi in 1:num_bikes
    bike_positions_x[bi] = rand(1:100)
    bike_positions_y[bi] = rand(1:100)
end

simulation_steps = zeros(Int64, (simulation_end, 2, num_bikes))
for t in 1:simulation_end
    for bi in 1:num_bikes
        new_x = bike_positions_x[bi]
        new_y = bike_positions_y[bi]
        if bitrand()[1]
            possible_new_x = bike_positions_x[bi] + rand((-1, 1))
            if possible_new_x > -1 && new_x < 100 + 1
                new_x = possible_new_x
            end
        end

        if bitrand()[1]
            possible_new_y = bike_positions_y[bi] + rand((-1, 1))
            if new_y > -1 && new_y < 100 + 1
                new_y = possible_new_y
            end
        end

        move_causes_collision = false
        for other_bi in 1:num_bikes
            if new_x == bike_positions_x[other_bi] && new_y == bike_positions_y[other_bi]
                move_causes_collision = true
                break
            end
        end

        if move_causes_collision
            simulation_steps[t, 1, bi] = bike_positions_x[bi]
            simulation_steps[t, 2, bi] = bike_positions_y[bi]
        else
            simulation_steps[t, 1, bi] = new_x
            simulation_steps[t, 2, bi] = new_y
        end

    end
end

for t in 1:simulation_end
    println("time: $(t)")
    for bi in 1:num_bikes
        println("bike $(bi): x: $(simulation_steps[t, 1, bi]) y: $(simulation_steps[t, 2, bi])")
    end
end


text:
using JuMP

function norm_stress1(N)
    m = Model()
    @variable(m, x[1:N])
    for i in 1:N
        @constraint(m, norm(x[j] for j in 1:i) <= i)
    end
    for i in 1:N
        @constraint(m, norm((i / j) * x[j] for j in 1:i) <= i)
    end
    for i in 1:N
        @constraint(
            m,
            norm((i / j) * x[j] for j in 1:i) + sum(x[j] for j in 1:N) <= i
        )
    end
end

# Warmup
norm_stress1(10)

@time for iter in 1:10
    norm_stress1(100)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function solve()
    a, b = [parse(Int, x) for x in split(readline())]
    (1 - b / a) * 100
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
include("src/weather/load_weather.jl")
include("src/route/domain.jl")
include("src/route/shortest_path.jl")
include("src/performance/polar.jl")

using BenchmarkTools
using Printf
using Dates

"""Run sample optimal routing with no current."""
function run_sample_simulation()
    path = ENV["HOME"]*"/Documents/sail_route.jl/src/data/first40_orgi.csv"
    twa, tws, perf = load_file(path)
    polar = setup_interpolation(tws, twa, perf)
    lon1 = 0.0
    lon2 = 1.0
    lat1 = 5.0
    lat2 = 5.0
    d_an, b_an = haversine(lon1, lat1, lon2, lat2)
    time_an = d_an/perf_interp(polar, min_angle(0.0, 90.0), 10.0)
    d_an_str = @sprintf("%0.2f", d_an)
    b_an_str = @sprintf("%0.2f", b_an)
    t_an_str = @sprintf("%0.2f", time_an)
    println("Analytical distance is $d_an_str nm")
    println("Analytical bearing is $b_an_str deg")
    println("Analytical time is $t_an_str hrs")
    start = Dates.DateTime(2000, 1, 1, 0, 0, 0)
    println("Start time ", start)
    x_nodes = 5
    y_nodes = 5
    y_dist = 100.0
    wisp, widi, cusp, cudi = sample_weather()
    st, ets = route_solve(lon1, lon2, lat1, lat2, x_nodes, y_nodes, y_dist, wisp, widi, cusp, cudi, polar, start)
    t_est_str = @sprintf("%0.2f", st)
    println("Estimated time: ", t_est_str, " hrs")
    if isapprox(st, time_an; rtol=0.000001) == false
        println("failed")
    else
        println("passed")
    end
end

run_sample_simulation()


text:
function _check_columns(db, table; min_nv=nothing, max_nv=nothing)
    println("Checking missing values of table $table")
    query = "SELECT * FROM $table"
    if !isnothing(min_nv) && !isnothing(max_nv)
         query *= " WHERE nb_vertex >= $min_nv AND nb_vertex <= $max_nv"
    end
    results = DBInterface.execute(db, query) |> DataFrame
    for col in names(results)
        values = results[!, col]
        nb_missing = sum(map(ismissing, values))
        nb_total = length(values)
        color = nb_missing == 0 ? Symbol("green") : nb_missing == nb_total ? Symbol("red") : Symbol("yellow")
        printstyled("$col: $nb_missing / $nb_total\n"; color=color)
    end
end

function _check_chordality(row)
    ismissing(row[:graph_path]) && return false
    g = loadgraph(row[:graph_path])
    return ischordal(g)
end

function _check_connectivity(row)
    ismissing(row[:graph_path]) && return false
    g = loadgraph(row[:graph_path])
    return is_connected(g)
end

function _check_self_loops(row)
    ismissing(row[:graph_path]) && return false
    g = loadgraph(row[:graph_path])
    return !has_self_loops(g)
end

function _check_index_clique(row)
    if ismissing(row[:clique_path]) || ismissing(row[:cliquetree_path])
        return false
    end
    clique = read_clique(row[:clique_path])
    cliquetree = read_cliquetree(row[:cliquetree_path])
    for edge in cliquetree
        src = clique[edge[1]]
        dst = clique[edge[2]]
        if isempty(intersect(src, dst))
            return false
        end
    end
    return true
end

function _check_source_graph(db, row)
    if ismissing(row[:origin_name]) || ismissing(row[:origin_scenario])
        return false
    end
    results = DBInterface.execute(db, "SELECT graph_path FROM instances WHERE name = '$(row[:origin_name])' AND scenario = $(row[:origin_scenario])") |> DataFrame
    source_graph = loadgraph(results[1, :graph_path])
    g = loadgraph(row[:graph_path])
    source_edges = edges(source_graph)
    g_edges = edges(g)
    return !isequal(source_graph, g) && all(in(g_edges), source_edges) && vertices(g) == vertices(source_graph)
end

function _check_serialize_graph(row)
    ismissing(row[:graph_path]) && return false
    graph_path = row[:graph_path]
    return isfile(graph_path)
end

function _check_serialize_network(row)
    ismissing(row[:pfn_path]) && return false
    pfn_path = row[:pfn_path]
    return isfile(pfn_path)
end

function _check_basic_feature(row)
    return !ismissing(row[:nbus]) &&
           !ismissing(row[:nbranch]) &&
           !ismissing(row[:nbranch_unique]) &&
           !ismissing(row[:ngen])
end

const check_functions = Dict(
    "chordality" => _check_chordality,
    "connectivity" => _check_connectivity,
    "self_loops" => _check_self_loops,
    "index_clique" => _check_index_clique,
    "source_graph" => _check_source_graph,
    "serialize_graph" => _check_serialize_graph,
    "serialize_network" => _check_serialize_network,
    "basic_feature" => _check_basic_feature
)
const need_db = Dict(
    "chordality" => false,
    "connectivity" => false,
    "self_loops" => false,
    "index_clique" => false,
    "source_graph" => true,
    "serialize_graph" => false,
    "serialize_network" => false,
    "basic_feature" => false
)
const valid_check_instance = Dict(
    "chordality" => true,
    "connectivity" => true,
    "self_loops" => true,
    "index_clique" => false,
    "source_graph" => false,
    "serialize_graph" => true,
    "serialize_network" => true,
    "basic_feature" => true,
)
const valid_check_decomposition = Dict(
    "chordality" => true,
    "connectivity" => true,
    "self_loops" => true,
    "index_clique" => true,
    "source_graph" => true,
    "serialize_graph" => true,
    "serialize_network" => false,
    "basic_feature" => false
)

function _check_sanity(db, rows, checks)
    for c in checks
        check_results = nothing
        if need_db[c]
            check_function(row) = check_functions[c](db, row)
            check_results = check_function.(eachrow(rows))
        else
            check_results = check_functions[c].(eachrow(rows))
        end
        if !all(check_results)
            printstyled("$(c) check failed.\n"; color=:red)
        else
            printstyled("$(c).\n"; color=:green)
        end
    end
end

function check_sanity(db::SQLite.DB, checks; min_nv=nothing, max_nv=nothing, subset_instance=nothing, subset_decomposition=nothing)
    _check_columns(db, "instances"; min_nv=min_nv, max_nv=max_nv)
    _check_columns(db, "decompositions"; min_nv=min_nv, max_nv=max_nv)
    query = "SELECT * FROM instances"
    if !isnothing(min_nv) && !isnothing(max_nv)
         query *= " WHERE nb_vertex >= $min_nv AND nb_vertex <= $max_nv"
    end
    if !isnothing(subset_instance)
        if !isnothing(min_nv) && !isnothing(max_nv)
            query *= " AND id IN ($(join(subset_instance, ',')))"
        else
            query *= " WHERE id IN ($(join(subset_instance, ',')))"
        end
    end
    println(query)
    results = DBInterface.execute(db, query) |> DataFrame
    printstyled("Checking table: instances\n"; bold=true)
    _check_sanity(db, results, filter(c -> valid_check_instance[c], checks))

    query = "SELECT * FROM decompositions"
    if !isnothing(min_nv) && !isnothing(max_nv)
         query *= " WHERE nb_vertex >= $min_nv AND nb_vertex <= $max_nv"
    end
    if !isnothing(subset_decomposition)
        if !isnothing(min_nv) && !isnothing(max_nv)
            query *= " AND id IN ($(join(subset_decomposition, ',')))"
        else
            query *= " WHERE id IN ($(join(subset_decomposition, ',')))"
        end
    end
    println(query)
    results = DBInterface.execute(db, query) |> DataFrame
    printstyled("\nChecking table: decompositions\n"; bold=true)
    _check_sanity(db, results, filter(c -> valid_check_decomposition[c], checks))
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
#
println("Ea) Two-step cascade computations and simulations following the 1s-3p photo-excitation of Si^-: AverageSCA model.")
using JLD
#
setDefaults("method: continuum, asymptotic Coulomb")    ## setDefaults("method: continuum, Galerkin")
setDefaults("method: normalization, pure sine")         ## setDefaults("method: normalization, pure Coulomb")    setDefaults("method: normalization, pure sine")
#
if  true    ## Stepwise decay cascade
setDefaults("print summary: open", "zzz-Cascade-computation.sum")

name = "Cascade after Si- 1s-3p excitation"
grid = Radial.Grid(false)
wa   = Cascade.Computation(Cascade.Computation(); name=name, nuclearModel=Nuclear.Model(14.), grid=grid, approach=Cascade.AverageSCA(),
                           scheme=Cascade.StepwiseDecayScheme([Auger, Radiative], 2, 0, Shell[], Shell[]),
                           initialConfigs=[Configuration("1s^1 2s^2 2p^6 3s^2 3p^4")] )
println(wa)
wb = perform(wa; output=true)
setDefaults("print summary: close", "")

else
setDefaults("print summary: open", "zzz-Cascade-simulation.sum")

## data = [JLD.load("zzz-cascade-decay-computations-2020-01-25T11.jld"), JLD.load("zzz-cascade-decay-computations-2020-01-25T11.jld")]
data = [JLD.load("zzz-cascade-decay-computations-2020-01-25T21.jld")]
name = "Simulation after Si- 1s-3p excitation"

wc   = Cascade.Simulation(Cascade.Simulation(), name=name, properties=Cascade.AbstractSimulationProperty[Cascade.IonDistribution()], ## , Cascade.FinalLevelDistribution()],
                          settings=Cascade.SimulationSettings(0., 0., 0., 0., [(1, 2.0), (2, 1.0), (3, 0.5)]), computationData=data )
println(wc)
wd = perform(wc; output=true)
setDefaults("print summary: close", "")

end


text:
using IbanGen
using Test

# shell> cd /Users/xor/work/julia/Iban
#  pkg> activate .
# test
# or:
# Pkg.test("IbanGen", test_args=["test_iban_from_string", "test_countries", "test_iban", "test_iban_random"])

active_tests = lowercase.(ARGS)

function addtest(fname)
    key = lowercase(splitext(fname)[1])
    if isempty(active_tests) || key in active_tests
        include(fname)
    end
end

addtest("test_countries.jl")
addtest("test_iban_from_string.jl")
addtest("test_iban.jl")
addtest("test_iban_random.jl")


text:
###############################################################################
#
# Update leaf allocation
#
###############################################################################
"""
    leaf_allocation!(
                node::SPACSimple{FT},
                laba::FT
    ) where {FT<:AbstractFloat}
    leaf_allocation!(
                node::SPACSimple{FT},
                isvmax::Bool,
                vmax::FT,
                v2r::FT = FT(0.015)
    ) where {FT<:AbstractFloat}
    leaf_allocation!(
                node::SPACSimple{FT},
                laba::FT,
                vmax::FT
    ) where {FT<:AbstractFloat}

Update leaf area and maximal carboxylation rate, given
- `node` [`SPACSimple`] type struct
- `laba` Given leaf area per basal area
- `vmax` Given Vcmax25
"""
function leaf_allocation!(
            node::SPACSimple{FT},
            laba::FT
) where {FT<:AbstractFloat}
    node.laba         = laba;
    node.hs.leaf.area = node.hs.stem.area * laba;
    node.lai          = laba / node.gaba;

    return nothing
end

function leaf_allocation!(
            node::SPACSimple{FT},
            isvmax::Bool,
            vmax::FT,
            v2r::FT = FT(0.015)
) where {FT<:AbstractFloat}
    node.ps.PSM.v_cmax25_ww = vmax;
    node.ps.PSM.v_cmax25    = vmax;
    node.ps.PSM.j_max25     = vmax * node.vtoj;
    node.ps.PSM.r_d25       = vmax * v2r;

    return nothing
end

function leaf_allocation!(
            node::SPACSimple{FT},
            laba::FT,
            vmax::FT
) where {FT<:AbstractFloat}
    leaf_allocation!(node, laba);
    leaf_allocation!(node, true, vmax);

    return nothing
end


text:
using DelayDiffEq, DiffEqProblemLibrary.DDEProblemLibrary
DDEProblemLibrary.importddeproblems()

using PGFPlotsX

using DelimitedFiles
using LinearAlgebra

function plot_waltman(alg; plotdensity = 10_000)
    # compute solution
    sol = solve(DDEProblemLibrary.prob_dde_RADAR5_waltman_5, MethodOfSteps(alg);
                reltol = 1e-6, abstol = [1e-21, 1e-21, 1e-21, 1e-21, 1e-9, 1e-9])
    ts = range(sol.t[1]; stop = sol.t[end], length = plotdensity)

    # load RADAR5 solution
    radar5 = readdlm(joinpath(@__DIR__, "..", "data", "waltman_radar5.out"))
    radar5_ts = radar5[:, 1]
    radar5_us = radar5[:, 2:7]

    # create log plot of first four components
    @pgf logaxis = SemiLogYAxis()
    for i in 1:4
        @pgf push!(logaxis, PlotInc({ thick }, Table(radar5_ts, radar5_us[:, i])),
                   PlotInc({ thick, dashed }, Table(ts, sol(ts; idxs = i).u)))
    end

    # create standard plot of last two components
    @pgf axis = Axis(
        PlotInc({ thick }, Table(radar5_ts, radar5_us[:, 5])),
        PlotInc({ thick, dashed }, Table(ts, sol(ts; idxs = 5).u)),
        PlotInc({ thick }, Table(radar5_ts, radar5_us[:, 6])),
        PlotInc({ thick, dashed }, Table(ts, sol(ts; idxs = 6).u)),
    )

    # create empty group plot
    @pgf gp = GroupPlot(
        {
            group_style = {
                group_size = "2 by 1",
                xlabels_at = "edge bottom",
                ylabels_at = "edge left",
                horizontal_sep = raw"0.1\textwidth",
            },
            xlabel = raw"$t$", xmin = ts[1], xmax = ts[end],
            width = raw"0.5\textwidth", height = raw"0.7\textheight",
            grid = "major",
            label_style = { font = raw"\footnotesize" },
            tick_label_style = { font = raw"\footnotesize" },
            legend_pos = "north east",
            no_marks
        },
        logaxis,
        axis
    )
end

# specify LinSolveFactorize due to https://github.com/YingboMa/RecursiveFactorization.jl/issues/4
pgfsave(joinpath(@__DIR__, "..", "figures", "generated", "waltman.tex"),
        plot_waltman(Rosenbrock23(linsolve = LinSolveFactorize(lu)));
        include_preamble = false)


text:
using StochasticOptimalTransport
using Distributions
using Documenter
using Random
using Test

Random.seed!(1234)

const SOT = StochasticOptimalTransport

@testset "StochasticOptimalTransport.jl" begin
    @testset "Utilities" begin
        include("utils.jl")
    end
    @testset "Discrete OT" begin
        include("discrete.jl")
    end
    @testset "Semi-discrete OT" begin
        include("semidiscrete.jl")
    end
    @testset "doctests" begin
        DocMeta.setdocmeta!(
            StochasticOptimalTransport,
            :DocTestSetup,
            :(using StochasticOptimalTransport);
            recursive=true,
        )
        doctest(StochasticOptimalTransport)
    end
end


text:
"""
DataPackage module
https://github.com/frictionlessdata/DataPackage.jl
"""
module DataPackage

export Package
export Resource
# export Profile

export PackageError

# using Base.Iterators: filter
# using Base.Iterators: Repeated, repeated

using JSON

using TableSchema
import TableSchema: Schema
import TableSchema: read, infer, validate, is_empty

import HTTP: request

include("exceptions.jl")
include("resource.jl")
include("package.jl")

export import_table

end # module

max_issues_repo_name: frictionlessdata/datapackage-jl

max_forks_repo_name: frictionlessdata/datapackage-jl


text:
"""
    one_fold(x,y,path,folds,fold) -> Vector

For a regularization path given by the `Int` vector `path`,
this function performs penalized linear regression on `x` and `y` and computes an out-of-sample error based on the indices given in `folds`.

Arguments:

- `x` is the `n` x `p` design matrix.
- `y` is the `n`-vector of responses.
- `path` is an `Int` vector that specifies which model sizes to include in the path, e.g. `path = collect(k0:increment:k_end)`.
- `folds` is an `Int` vector indicating which component of `y` goes to which fold, e.g. `folds = IHT.cv_get_folds(n,nfolds)`
- `fold` is the current fold to compute.

Optional Arguments:

- `tol` is the convergence tolerance to pass to the path computations. Defaults to `1e-4`.
- `max_iter` caps the number of permissible iterations in the IHT algorithm. Defaults to `1000`.
- `max_step` caps the number of permissible backtracking steps. Defaults to `50`.
- `quiet` is a Boolean to activate output. Defaults to `true` (no output).

Output:

- `errors` is a vector of out-of-sample errors (MSEs) for the current fold.
"""
function one_fold{T <: Float}(
    x        :: DenseMatrix{T},
    y        :: DenseVector{T},
    path     :: DenseVector{Int},
    folds    :: DenseVector{Int},
    fold     :: Int;
    tol      :: T    = convert(T, 1e-4),
    max_iter :: Int  = 100,
    max_step :: Int  = 50,
    quiet    :: Bool = true,
)

    # make vector of indices for folds
    test_idx  = folds .== fold
    test_size = sum(test_idx)

    # train_idx is the vector that indexes the TRAINING set
    train_idx = .!test_idx

    # allocate the arrays for the training set
    x_train = x[train_idx,:]
    y_train = y[train_idx]

    # compute the regularization path on the training set
    betas = iht_path(x_train, y_train, path, tol=tol, max_iter=max_iter, quiet=quiet, max_step=max_step)

    # compute the mean out-of-sample error for the TEST set
    Xb = view(x, test_idx, :) * betas
#    r  = broadcast(-, view(y, test_idx, 1), Xb)
    r  = broadcast(-, y[test_idx], Xb)
    er = vec(sum(abs2, r, 1)) ./ (2*test_size)

    return er :: Vector{T}
end

"""
    pfold(x,y,path,folds,q [, pids=procs()]) -> Vector

This function is the parallel execution kernel in `cv_iht()`. It is not meant to be called outside of `cv_iht()`.
For floating point data `x` and `y` and an integer vector `path`, `pfold` will distribute `q` crossvalidation folds across the processes supplied by the optional argument `pids`.
It calls `one_fold()` for each fold, then collects the vectors of MSEs each process, applys a reduction, and finally returns the average MSE vector across all folds.
"""
function pfold{T <: Float}(
    x        :: DenseMatrix{T},
    y        :: DenseVector{T},
    path     :: DenseVector{Int},
    folds    :: DenseVector{Int},
    q        :: Int;
    pids     :: Vector{Int} = procs(),
    tol      :: T    = convert(T, 1e-4),
    max_iter :: Int  = 100,
    max_step :: Int  = 50,
    quiet    :: Bool = true,
)
    # how many CPU processes can pfold use?
    np = length(pids)

    # report on CPU processes
    quiet || println("pfold: np = ", np)
    quiet || println("pids = ", pids)

    # set up function to share state (indices of folds)
    i = 1
    nextidx() = (idx=i; i+=1; idx)

    # preallocate cell array for results
    results = SharedArray{T}((length(path),q), pids=pids) :: SharedMatrix{T}

    # master process will distribute tasks to workers
    # master synchronizes results at end before returning
    @sync begin

        # loop over all workers
        for worker in pids

            # exclude process that launched pfold, unless only one process is available
            if worker != myid() || np == 1

                # asynchronously distribute tasks
                @async begin
                    while true

                        # grab next fold
                        current_fold = nextidx()

                        # if current fold exceeds total number of folds then exit loop
                        current_fold > q && break

                        # report distribution of fold to worker and device
                        quiet || print_with_color(:blue, "Computing fold $current_fold on worker $worker.\n\n")

                        # launch job on worker
                        # worker loads data from file paths and then computes the errors in one fold
                        r = remotecall_fetch(worker) do
                                one_fold(x, y, path, folds, current_fold, tol=tol, max_iter=max_iter, max_step=max_step, quiet=quiet)
                        end # end remotecall_fetch()
                        setindex!(results, r, :, current_fold)
                    end # end while
                end # end @async
            end # end if
        end # end for
    end # end @sync

    # return reduction (row-wise sum) over results
    return (vec(sum(results, 2) ./ q)) :: Vector{T}
end

"""
    cv_iht(x,y) -> Vector

This function will perform `q`-fold cross validation for the ideal model size in IHT least squares regression using the `n` x `p` design matrix `x` and the response vector `y`.
Each path is asynchronously spawned using any available processor.
For each path, one fold is held out of the analysis for testing, while the rest of the data are used for training.
The function to compute each path, `one_fold()`, will return a vector of out-of-sample errors (MSEs).

Arguments:

- `x` is the `n` x `p` design matrix.
- `y` is the `n`-vector of responses.

Optional Arguments:

- `q` is the number of folds to compute. Defaults to `max(3, min(CPU_CORES, 5))`, where `CPU_CORES`is the Julia variable to query the number of available CPU cores.
- `path` is an `Int` vector that specifies which model sizes to include in the path. Defaults to `path = collect(1:min(p,20))`.
- `folds` is the partition of the data. Defaults to `IHT.cv_get_folds(n,q)`.
- `pids`, a vector of process IDs. Defaults to `procs()`, which recruits all available processes.
- `tol` is the convergence tolerance to pass to the path computations. Defaults to `1e-4`.
- `max_iter` caps the number of permissible iterations in the IHT algorithm. Defaults to `100`.
- `max_step` caps the number of permissible backtracking steps. Defaults to `50`.
- `quiet` is a Boolean to activate output. Defaults to `true` (no output).
   *NOTA BENE*: each processor outputs feed to the console without regard to the others,
   so setting `quiet = false` can yield very messy output!

Output:

An `IHTCrossvalidationResults` object with the following fields:

- `mses` is the averaged MSE over all folds.
- `k` is the best crossvalidated model size.
- `path` is the regularization path used in the crossvalidation.
- `b`, a vector of `k_star` floats
- `bidx`, a vector of `k_star` indices indicating the support of the best model.
"""
function cv_iht{T <: Float}(
    x        :: DenseMatrix{T},
    y        :: DenseVector{T};
    q        :: Int   = cv_get_num_folds(3, 5),
    path     :: DenseVector{Int} = collect(1:min(size(x,2),20)),
    folds    :: DenseVector{Int} = cv_get_folds(sdata(y),q),
    pids     :: Vector{Int}      = procs(),
    tol      :: T    = convert(T, 1e-4),
    max_iter :: Int  = 100,
    max_step :: Int  = 50,
    quiet    :: Bool = true,
)
    # do not allow crossvalidation with fewer than 3 folds
    q > 2 || throw(ArgumentError("Number of folds q = $q must be at least 3."))

    # problem dimensions?
    n,p = size(x)

    # check for conformable arrays
    n == length(y) || throw(DimensionMismatch("Row dimension of x ($n) must match number of rows in y ($(length(y)))"))

    # how many elements are in the path?
    nmodels = length(path)

    # want to compute a path for each fold
    # the folds are computed asynchronously over processes enumerated by pids
    # master process then reduces errors across folds and returns MSEs
    mses = pfold(x, y, path, folds, q, pids=pids, tol=tol, max_iter=max_iter, max_step=max_step, quiet=quiet)

    # what is the best model size?
    # if there are multiple model sizes of EXACTLY the same MSE,
    # then this chooses the smaller of the two
    k = path[indmin(mses)] :: Int

    # print results
    !quiet && print_cv_results(mses, path, k)

    # refit the best model
#    b, bidx = refit_iht(sdata(x), sdata(y), k, tol=tol, max_iter=max_iter, max_step=max_step, quiet=quiet)
    output = L0_reg(sdata(x), sdata(y), k, max_iter=max_iter, max_step = max_step, quiet=quiet, tol=tol)
    bidx = find(output.beta)
    b = output.beta[bidx]

    return IHTCrossvalidationResults(mses, sdata(path), b, bidx, k)
end


text:
import Base.<=

export HPolygonOpt

"""
    HPolygonOpt{N, VN<:AbstractVector{N}} <: AbstractHPolygon{N}

Type that represents a convex polygon in constraint representation whose edges
are sorted in counter-clockwise fashion with respect to their normal directions.
This is a refined version of [`HPolygon`](@ref).

### Fields

- `constraints`       -- list of linear constraints, sorted by the normal
                         direction in counter-clockwise fashion
- `ind`               -- index in the list of constraints to begin the search
                         to evaluate the support function
- `sort_constraints`  -- (optional, default: `true`) flag for sorting the
                         constraints (sortedness is a running assumption of this
                         type)
- `check_boundedness` -- (optional, default: `false`) flag for checking if the
                         constraints make the polygon bounded; (boundedness is a
                         running assumption of this type)
- `prune`             -- (optional, default: `true`) flag for removing redundant
                         constraints

### Notes

This structure is optimized to evaluate the support function/vector with a large
sequence of directions that are close to each other. The strategy is to have an
index that can be used to warm-start the search for optimal values in the
support vector computation.

The option `sort_constraints` can be used to deactivate automatic sorting of
constraints in counter-clockwise fashion, which is an invariant of this type.
Alternatively, one can construct an `HPolygonOpt` with empty constraints list,
which can then be filled iteratively using `addconstraint!`.

Similarly, the option `prune` can be used to deactivate automatic pruning of
redundant constraints.

Another type assumption is that the polygon is bounded.
The option `check_boundedness` can be used to assert this.
This option is deactivated by default because we explicitly want to allow the
iterative addition of the constraints, and hence one has to initially construct
an empty list of constraints (which represents an unbounded set).
The user has to make sure that the `HPolygonOpt` is not used before the
constraints actually describe a bounded set.
The function `isbounded` can be used to manually assert boundedness.
"""
mutable struct HPolygonOpt{N, VN<:AbstractVector{N}} <: AbstractHPolygon{N}
    constraints::Vector{LinearConstraint{N, VN}}
    ind::Int

    # default constructor that applies sorting of the given constraints
    function HPolygonOpt(constraints::Vector{LinearConstraint{N, VN}},
                         ind::Int=1;
                         sort_constraints::Bool=true,
                         check_boundedness::Bool=false,
                         prune::Bool=true) where {N, VN<:AbstractVector{N}}
        if sort_constraints
            sorted_constraints = Vector{LinearConstraint{N, VN}}()
            sizehint!(sorted_constraints, length(constraints))
            for ci in constraints
                addconstraint!(sorted_constraints, ci; prune=prune)
            end
            P = new{N, VN}(sorted_constraints, ind)
        else
            P = new{N, VN}(constraints, ind)
        end
        @assert (!check_boundedness ||
                 isbounded(P, false)) "the polygon is not bounded"
        return P
    end
end

isoperationtype(::Type{<:HPolygonOpt}) = false
isconvextype(::Type{<:HPolygonOpt}) = true

# constructor for an HPolygon with no constraints
function HPolygonOpt{N, VN}() where {N, VN<:AbstractVector{N}}
    HPolygonOpt(Vector{LinearConstraint{N, VN}}())
end

# constructor for an HPolygonOpt with no constraints and given numeric type
function HPolygonOpt{N}() where {N}
    HPolygonOpt(Vector{LinearConstraint{N, Vector{N}}}())
end

# constructor for an HPolygonOpt without explicit numeric type, defaults to Float64
function HPolygonOpt()
    HPolygonOpt{Float64}()
end

# constructor from a simple H-representation
HPolygonOpt(A::AbstractMatrix,
            b::AbstractVector;
            sort_constraints::Bool=true,
            check_boundedness::Bool=false,
            prune::Bool=true) =
    HPolygonOpt(constraints_list(A, b); sort_constraints=sort_constraints,
                check_boundedness=check_boundedness, prune=prune)

# --- LazySet interface functions ---

"""
    σ(d::AbstractVector, P::HPolygonOpt;
      [linear_search]::Bool=(length(P.constraints) < BINARY_SEARCH_THRESHOLD))

Return the support vector of an optimized polygon in a given direction.

### Input

- `d`             -- direction
- `P`             -- optimized polygon in constraint representation
- `linear_search` -- (optional, default: see below) flag for controlling whether
                     to perform a linear search or a binary search

### Output

The support vector in the given direction.
The result is always one of the vertices; in particular, if the direction has
norm zero, any vertex is returned.

### Algorithm

Comparison of directions is performed using polar angles; see the overload of
`<=` for two-dimensional vectors.

For polygons with `BINARY_SEARCH_THRESHOLD = 10` or more constraints we use a
binary search by default.
"""
function σ(d::AbstractVector, P::HPolygonOpt;
           linear_search::Bool=(length(P.constraints) < BINARY_SEARCH_THRESHOLD))
    n = length(P.constraints)
    @assert n > 0 "the polygon has no constraints"
    if linear_search
        # linear search
        if (d <= P.constraints[P.ind].a)
            # search backward
            k = P.ind-1
            while (k >= 1 && d <= P.constraints[k].a)
                k -= 1
            end
            if (k == 0)
                P.ind = n
                # corner case: wrap-around in constraints list
                return element(intersection(Line2D(P.constraints[n]),
                                            Line2D(P.constraints[1])))
            else
                P.ind = k
            end
        else
            # search forward
            k = P.ind+1
            while (k <= n && P.constraints[k].a <= d)
                k += 1
            end
            if (k == n+1)
                P.ind = n
                # corner case: wrap-around in constraints list
                return element(intersection(Line2D(P.constraints[n]),
                                            Line2D(P.constraints[1])))
            else
                P.ind = k-1
            end
        end
        return element(intersection(Line2D(P.constraints[P.ind]),
                                    Line2D(P.constraints[P.ind + 1])))
    else
        # binary search
        k = binary_search_constraints(d, P.constraints, n, P.ind)
        if k == 1 || k == n+1
            P.ind = 1
            # corner cases: wrap-around in constraints list
            return element(intersection(Line2D(P.constraints[n]),
                                        Line2D(P.constraints[1])))
        else
            P.ind = k
            return element(intersection(Line2D(P.constraints[k-1]),
                                        Line2D(P.constraints[k])))
        end
    end
end

"""
    translate(P::HPolygonOpt, v::AbstractVector; share::Bool=false)

Translate (i.e., shift) an optimized polygon in constraint representation by a
given vector.

### Input

- `P`     -- optimized polygon in constraint representation
- `v`     -- translation vector
- `share` -- (optional, default: `false`) flag for sharing unmodified parts of
             the original set representation

### Output

A translated optimized polygon in constraint representation.

### Notes

The normal vectors of the constraints (vector `a` in `a⋅x ≤ b`) are shared with
the original constraints if `share == true`.

### Algorithm

We translate every constraint.
"""
function translate(P::HPolygonOpt, v::AbstractVector; share::Bool=false)
    @assert length(v) == dim(P) "cannot translate a $(dim(P))-dimensional " *
                                "set by a $(length(v))-dimensional vector"
    constraints = [translate(c, v; share=share) for c in constraints_list(P)]
    return HPolygonOpt(constraints, P.ind;
                       sort_constraints=false, check_boundedness=false,
                       prune=false)
end


text:
using NumericalMethodsforEngineers, Test

a = [0. 0. 16.; 0. 4. 5.; 8. -4. 22.]
b = Float64[4, 2, 5]

c = a \ b
c |> display

println("\nCholesky factorization using banded storage\n")
println("Banded Coefficient Matrix:")
a |> display

println("\nRight Hand side:")
b |> display

cholin!(a)
println("\nL in Band Form:")
a |> display

chobac!(a, b)
println("\nSolution Vector:")
b |> display


text:
export odn_to_grid,grid_to_odn,fwi_wavelet,index_block

function index_block(idx,block_size::I) where {I<:Integer}
    m = div(length(idx),block_size)
    out_idx = [(i+1:i+block_size) for i in (0:m-1)*block_size]
    if length(idx) > m*block_size
        push!(out_idx, m*block_size+1:length(idx))
    end
    return [idx[out_idx[i]] for i in 1:length(out_idx)]
end

function fwi_wavelet(freqs::AbstractArray{F,1},t0::F,f0::F) where {F<:Real}
    w = exp.(-2*π*im*t0*freqs)
    if f0 > 0
        @. w = freqs^2 * exp(-(freqs./f0)^2)*w
    end
    return w
end

# Convert odn coordinates to their full grid coordinates
function odn_to_grid(o::AbstractArray{F,1},d::AbstractArray{F,1},n::AbstractArray{I,1}) where {F<:AbstractFloat, I<:Integer}
    (length(o)==length(d) && length(o)==length(n)) || throw(Exception("o,d,n must have the same length"))
    x = Array{StepRangeLen,1}(length(o))
    for i in 1:length(o)
        x[i] = o[i] + (0:(n[i]-1))*d[i]
    end
    return tuple(x...)
end

# Convert grid coordinates to odn coordinates
function grid_to_odn(x::AbstractArray{StepRangeLen{F},1}) where {F<:AbstractFloat}
    o = Array{F,1}()
    d = Array{F,1}()
    n = Array{Int64,1}()
    for i in 1:length(o)
        o[i] = x[i][1]
        d[i] = x[i][2] - x[i][1]
        n[i] = length(x[i])
    end
    return (o,d,n)
end


text:
mutable struct Fractal{D, T, N, M}
    segments::Deque{Segment{D, T}}
    step_in::MVector{N, Segment{D, T}}
    step_out::MVector{M, Segment{D, T}}
    step!::Function

    state::Int
end

function Fractal{D, T}(input::Vector{SVector{D, T}}, M::Int, step!::Function)
    @assert length(input) > 1

    segments = Deque{Segment{D, T}}()
    for (v0, v1) in zip(input[1:end-1], input[2:end])
        push!(segments, Segment(v0, v1))
    end

    N = length(input) - 1
    step_in = MVector((Segment(SVector(zeros(T, D)...), SVector(zeros(T, D)...)) for i in 1:N)...)

    step_out = MVector((Segment(SVector(zeros(T, D)...), SVector(zeros(T, D)...)) for i in 1:M)...)

    Fractal(segments, step_in, step_out, step!, length(segments))
end

start{D, T, N, M}(f::Fractal{D, T, N, M}) = f.state

function next{D, T, N, M}(f::Fractal{D, T, N, M}, state::Int)

    taken, state = state, 0
    while taken > 0
        taken -= N

        for i in 1:N
            f.step_in[i] = shift!(f.segments)
        end
        f.step!(f.step_in, f.step_out)

        for j in 1:M
            push!(f.segments, f.step_out[j])
        end
        state += M
    end
    @assert length(f.segments) == state
    f.state = state
    state
end

done{D, T, N, M}(f::Fractal{D, T, N, M}) = true

function refine!{D, T, N, M}(f::Fractal{D, T, N, M}, count::Int)
    @assert count >= 1
    next(f, f.state)
    (count == 1) ? f : refine!(f, count-1)
end

function koch_t!{T}(input::MVector{1, Segment{2, T}}, output::MVector{4, Segment{2, T}})
    v0, v1 = input[1].A, input[1].B

    v13 = v0*2/3 + v1*1/3
    v23 = v0*1/3 + v1*2/3

    t = (v1 - v0)

    shift = SVector(t[2], -t[1])

    vmid = 0.5*(v0 + v1) + shift*sqrt(3)/6

    # New points
    output[1] = Segment(v0, v13)
    output[2] = Segment(v13, vmid)
    output[3] = Segment(vmid, v23)
    output[4] = Segment(v23, v1)
end

function KochTFractal{T}(A::SVector{2, T}=SVector(0., 0.), B::SVector{2, T}=SVector(1., 0.))
    Fractal([A, B], 4, koch_t!)
end

# --------------------------------------------------------------------

function koch_q1!{T}(input::MVector{1, Segment{2, T}}, output::MVector{5, Segment{2, T}})
    v0, v1 = input[1].A, input[1].B

    v13 = v0*2/3 + v1*1/3
    v23 = v0*1/3 + v1*2/3

    t = (v1 - v0)
    shift = SVector{2, Float64}(t[2], -t[1])
    shift /= 3

    # New points
    #[v13, v13+shift, v23+shift, v23]
    shift_v13 = v13 + shift
    shift_v23 = v23 + shift

    output[1] = Segment(v0, v13)
    output[2] = Segment(v13, shift_v13)
    output[3] = Segment(shift_v13, shift_v23)
    output[4] = Segment(shift_v23, v23)
    output[5] = Segment(v23, v1)
end

function KochQ1Fractal{T}(A::SVector{2, T}=SVector(0., 0.), B::SVector{2, T}=SVector(1., 0.))
    Fractal([A, B], 5, koch_q1!)
end

# --------------------------------------------------------------------

function koch_q2!{T}(input::MVector{1, Segment{2, T}}, output::MVector{8, Segment{2, T}})
    v0, v1 = input[1].A, input[1].B

    v14 = v0*3/4 + v1*1/4
    v34 = v0*1/4 + v1*3/4
    vmid = v0/2 + v1/2

    t = (v1 - v0)

    shift = SVector{2, Float64}(t[2], -t[1])
    shift /= 4

    shift_v14 = v14 + shift
    up_vmid = vmid + shift
    down_vmid = vmid - shift
    shift_v34 = v34 - shift

    # New points
    # [v14, v14+shift, vmid+shift, vmid, vmid-shift, v34-shift, v34]
    output[1] = Segment(v0, v14)
    output[2] = Segment(v14, shift_v14)
    output[3] = Segment(shift_v14, up_vmid)
    output[4] = Segment(up_vmid, vmid)
    output[5] = Segment(vmid, down_vmid)
    output[6] = Segment(down_vmid, shift_v34)
    output[7] = Segment(shift_v34, v34)
    output[8] = Segment(v34, down_v1)

end

function KochQ2Fractal{T}(A::SVector{2, T}=SVector(0., 0.), B::SVector{2, T}=SVector(1., 0.))
    Fractal([A, B], 8, koch_q2!)
end

# --------------------------------------------------------------------

function levy!{T}(input::MVector{1, Segment{2, T}}, output::MVector{2, Segment{2, T}})
    v0, v1 = input[1].A, input[1].B

    t = (v1 - v0)

    shift = SVector(t[2], -t[1])

    vmid = 0.5*(v0 + v1) + shift/2

    # New points
    output[1] = Segment(v0, vmid)
    output[2] = Segment(vmid, v1)
end

function LevyFractal{T}(A::SVector{2, T}=SVector(0., 0.), B::SVector{2, T}=SVector(1., 0.))
    Fractal([A, B], 2, levy!)
end

# --------------------------------------------------------------------

function dragon!{T}(input::MVector{2, Segment{2, T}}, output::MVector{4, Segment{2, T}})
    v0, v1, v2 = input[1].A, input[1].B, input[2].B

    v01 = v0/2 + v1/2
    v12 = v1/2 + v2/2

    t1 = (v1 - v0)
    v01 += SVector(t1[2], -t1[1])/2

    t2 = (v2 - v1)
    v12 -= SVector(t2[2], -t2[1])/2

    output[1] = Segment(v0, v01)
    output[2] = Segment(v01, v1)
    output[3] = Segment(v1, v12)
    output[4] = Segment(v12, v2)
end

function DragonFractal{T}(A::SVector{2, T}=SVector(0., 0.),
                          B::SVector{2, T}=SVector(0.5, -0.5),
                          C::SVector{2, T}=SVector(1., 0.))
    Fractal([A, B, C], 4, dragon!)
end

# --------------------------------------------------------------------

function sierpinski_arrow_head!{T}(input::MVector{3, Segment{2, T}}, output::MVector{9, Segment{2, T}})
    v0, v1, v2, v3 = input[1].A, input[2].A, input[3].A, input[3].B

    vmid = (v0 + v3)/2

    t = (v3 - v0)

    v01 = (v0 + vmid)/2
    v10 = (v1 + vmid)/2

    v23 = (v2 + vmid)/2
    v32 = (v3 + vmid)/2

    vtop = vmid + 2*((v1+v2)/2 - vmid)

    v12 = (v1 + vtop)/2
    v21 = (v2 + vtop)/2

    output[1] = Segment(v0, v01)
    output[2] = Segment(v01, v10)
    output[3] = Segment(v10, v1)
    output[4] = Segment(v1, v12)
    output[5] = Segment(v12, v21)
    output[6] = Segment(v21, v2)
    output[7] = Segment(v2, v23)
    output[8] = Segment(v23, v32)
    output[9] = Segment(v32, v3)
end

function SierpinskiArrowHeadFractal{T}(A::SVector{2, T}=SVector(0., 0.),
                                       B::SVector{2, T}=SVector(0.25, sqrt(3)/4.),
                                       C::SVector{2, T}=SVector(0.75, sqrt(3)/4.),
                                       D::SVector{2, T}=SVector(1., 0.))
    Fractal([A, B, C, D], 9, sierpinski_arrow_head!)
end

# --------------------------------------------------------------------

function fractal_curve{D, T, N, M}(f::Fractal{D, T, N, M})
    segments = Vector{Segment{D, T}}(length(f.segments))
    for (i, seg) in enumerate(f.segments)
        segments[i] = seg
    end

    Curve(segments)
end

Curve{D, T, N, M}(f::Fractal{D, T, N, M}) = Curve(collect(f.segments))


text:
# constant and time-varying bounds and bound pairs

# TODO:
# - add bounds constructors with maxt, to extend/shorten given bounds

# single bounds

abstract type AbstractBound end
getdt(b::AbstractBound) = b.dt

struct ConstBound <: AbstractBound
    b::Float64
    dt::Float64

    function ConstBound(b::Real, dt::Real)
        b > zero(b) || error("bound needs to be positive")
        dt > zero(dt) || error("dt needs to be positive")
        return new(float(b), float(dt))
    end
end
getbound(b::ConstBound, n::Int) = b.b
getbound(b::ConstBound) = b.b
getboundgrad(b::ConstBound, n::Int) = 0.0
getmaxn(b::ConstBound) = typemax(Int)

struct VarBound <: AbstractBound
    b::Vector{Float64}
    bg::Vector{Float64}
    dt::Float64

    function VarBound(b::AbstractVector{T1}, bg::AbstractVector{T2},
                      dt::Real) where {T1 <: Real, T2 <: Real}
        dt > zero(dt) || error("dt needs to be positive")
        length(b) > 0 || error("bounds need to be of non-zero length")
        length(b) == length(bg) || error("b and bg need to be of same length")
        b[1] > 0.0 || error("b[1] needs to be positive")
        return new(b, bg, float(dt))
    end
    VarBound(b::AbstractVector{T}, dt::Real) where T <: Real = VarBound(b, fdgrad(b, dt), dt)
end
getbound(b::VarBound, n::Int) = b.b[n]
getboundgrad(b::VarBound, n::Int) = b.bg[n]
getmaxn(b::VarBound) = length(b.b)

struct LinearBound <: AbstractBound
    b0::Float64
    bslope::Float64
    dtbslope::Float64
    dt::Float64
    maxn::Int

    function LinearBound(b0::Real, bslope::Real, dt::Real, maxt::Real=Inf)
        dt > zero(dt) || error("dt needs to be positive")
        b0 > zero(dt) || error("b0 needs to be positive")
        maxt > zero(maxt) || error("maxt needs to be positive")
        maxn = maxt / dt
        # subtracting dt * bslope ensures that getbound(b, 1) = b0
        return new(b0 - dt * bslope, bslope, dt * bslope, dt,
            isfinite(maxn) ? ceil(Int, maxn) : typemax(Int))
    end
end
getbound(b::LinearBound, n::Int) = b.b0 + n * b.dtbslope
getboundgrad(b::LinearBound, n::Int) = b.bslope
getmaxn(b::LinearBound) = b.maxn

# bound pairs

abstract type AbstractBounds end

struct SymBounds{T <: AbstractBound} <: AbstractBounds
    b::T

    SymBounds{T}(b::T) where T = new(b)
end
getdt(b::SymBounds) = getdt(b.b)
getubound(b::SymBounds, n::Int) = getbound(b.b, n)
getlbound(b::SymBounds, n::Int) = -getbound(b.b, n)
getuboundgrad(b::SymBounds, n::Int) = getboundgrad(b.b, n)
getlboundgrad(b::SymBounds, n::Int) = -getboundgrad(b.b, n)
getmaxn(b::SymBounds) = getmaxn(b.b)

const VarSymBounds = SymBounds{VarBound}

const ConstSymBounds = SymBounds{ConstBound}
ConstSymBounds(b::Real, dt::Real) = ConstSymBounds(ConstBound(b, dt))
getbound(b::ConstSymBounds) = b.b.b

const LinearSymBounds = SymBounds{LinearBound}
LinearSymBounds(b0::Real, bslope::Real, dt::Real) = LinearSymBounds(LinearBound(b0, bslope, dt))

struct AsymBounds{T1 <: AbstractBound, T2 <: AbstractBound} <: AbstractBounds
    upper::T1
    lower::T2

    function AsymBounds{T1,T2}(upper::T1, lower::T2) where {T1,T2}
        getdt(upper) == getdt(lower) || error("Bounds need to have equal dt")
        return new(upper, lower)
    end
end
getdt(b::AsymBounds) = getdt(b.upper)
getubound(b::AsymBounds, n::Int) = getbound(b.upper, n)
getlbound(b::AsymBounds, n::Int) = -getbound(b.lower, n)
getuboundgrad(b::AsymBounds, n::Int) = getboundgrad(b.upper, n)
getlboundgrad(b::AsymBounds, n::Int) = -getboundgrad(b.lower, n)
getmaxn(b::AsymBounds) = min(getmaxn(b.upper), getmaxn(b.lower))

const VarAsymBounds = AsymBounds{VarBound, VarBound}

const ConstAsymBounds = AsymBounds{ConstBound, ConstBound}
ConstAsymBounds(upper::Real, lower::Real, dt::Real) =
    ConstAsymBounds(ConstBound(upper, dt), ConstBound(-lower, dt))
getubound(b::ConstAsymBounds) = getbound(b.upper)
getlbound(b::ConstAsymBounds) = -getbound(b.lower)

const ConstBounds = Union{ConstSymBounds, ConstAsymBounds}


text:
# Copyright (c) Microsoft Corporation.
# Licensed under the MIT license.

#Pkg.add("ForwardDiff")
#Pkg.update()
#Pkg.checkout("ForwardDiff")
#Pkg.status()
using ForwardDiff

include("common_io.jl")

typealias iVec Vector{Int}
typealias Vec Vector{Float64}
typealias Mat Matrix{Float64}
typealias Tensor3 Array{Float64,3}

type HandModel
  bone_names::Vector{String}
  parents::iVec
  base_relatives::Tensor3
  inverse_base_absolutes::Tensor3
  base_positions::Mat
  weights::Mat
  is_mirrored::Bool
end

type HandData
  model::HandModel
  correspondences::iVec
  points::Mat
end

#################### IO ############################
function read_hand_instance(model_dir,fn_in)
  bones_fn = joinpath(model_dir,"bones.txt")

  A = readdlm(bones_fn, ':')

  n_bones = size(A,1)
  bone_names = A[:,1]
  parents = A[:,2] + 1 #julia indexing

  transforms = A[:,3:18]
  transforms = permutedims(reshape(transforms,(n_bones,4,4)), (1,3,2))

  inverse_absolute_transforms  = A[:,19:34]
  inverse_absolute_transforms = permutedims(reshape(
    inverse_absolute_transforms,(n_bones,4,4)),(1,3,2));

  vertices_fn = joinpath(model_dir,"vertices.txt")
  A = readdlm(vertices_fn,':');

  base_positions = A[:,1:3]';
  n_vertices = size(A,1);

  weights = zeros(Float64,n_bones,n_vertices);
  for i_vert = 1:n_vertices
      for i=0:A[i_vert,9]-1
        i_bone = A[i_vert,9 + i*2 + 1] + 1; #julia indexing
        weights[i_bone, i_vert] = A[i_vert,9 + i*2 + 2];
      end
  end

  fid = open(fn_in)
  lines = readlines(fid)
  close(fid)
  line=split(lines[1]," ")
  n_corrs = parse(Int,line[1])
  n_params = parse(Int,line[2])
  off = 1

  corrs = zeros(Int,n_corrs)
  pts = zeros(Float64,3,n_corrs)
  for i in 1:n_corrs
    line=split(lines[i+off]," ")
    corrs[i] = parse(Int,line[1]) + 1 #julia indexing
    for j in 1:3
      pts[j,i] = parse(Float64,line[j+1])
    end
  end
  off += n_corrs

  params = zeros(Float64,n_params)
  for i in 1:n_params
    params[i] = parse(Float64,lines[i+off])
  end

  model = HandModel(bone_names,parents,transforms,
                    inverse_absolute_transforms,
                    [base_positions; ones(Float64,1,size(base_positions,2))],
                    weights,false)

  data = HandData(model, corrs, pts)

  (params, data)
end

##################### objective #############################

function euler_angles_to_rotation_matrix(xyz)
  tx = xyz[1]
  ty = xyz[2]
  tz = xyz[3]
  Rx = [1 0 0; 0  cos(tx) -sin(tx); 0 sin(tx) cos(tx)]
  Ry = [cos(ty) 0 sin(ty); 0 1 0; -sin(ty) 0 cos(ty)]
  Rz = [cos(tz) -sin(tz) 0; sin(tz) cos(tz) 0; 0 0 1]
  return Rz*Ry*Rx;
end

function get_posed_relatives(model, pose_params)
# default parametrization xzy # Flexion, Abduction, Twist
  order = [1, 3, 2]
  offset = 3
  n_bones = size(model.bone_names,1)
  relatives = zeros(eltype(pose_params),n_bones,4,4)

  for i_bone = 1:n_bones
    T = eye(eltype(pose_params),4)
    T[1:3,1:3] = euler_angles_to_rotation_matrix(pose_params[order,i_bone+offset])
    relatives[i_bone,:,:] = squeeze(model.base_relatives[i_bone,:,:],1) * T
  end
  relatives
end

function relatives_to_absolutes(relatives, parents)
  absolutes = zeros(eltype(relatives),size(relatives));
  for i=1:length(parents)
    if parents[i] == 0
      absolutes[i,:,:] = relatives[i,:,:]
    else
      absolutes[i,:,:] = squeeze(absolutes[parents[i],:,:],1) *
               squeeze(relatives[i,:,:],1)
    end
  end
  absolutes
end

function angle_axis_to_rotation_matrix(angle_axis)
  n = sqrt(sum(angle_axis.^2));
  if n < .0001
      return eye(eltype(angle_axis),3)
  end

  x = angle_axis[1] / n
  y = angle_axis[2] / n
  z = angle_axis[3] / n

  s = sin(n)
  c = cos(n)

  R = [x*x+(1-x*x)*c x*y*(1-c)-z*s x*z*(1-c)+y*s;
    x*y*(1-c)+z*s y*y+(1-y*y)*c y*z*(1-c)-x*s;
    x*z*(1-c)-y*s z*y*(1-c)+x*s z*z+(1-z*z)*c]
  return R
end

function apply_global_transform(pose_params, positions)
  T = eye(eltype(pose_params),3,4)
  T[:,1:3] = angle_axis_to_rotation_matrix(pose_params[:,1]);
  T[:,1:3] = T[1:3,1:3] .* pose_params[:,2]'
  T[:,4] = pose_params[:,3]

  return T * [positions; ones(eltype(positions),1,size(positions,2))]
end

function get_skinned_vertex_positions(model, pose_params)
  relatives = get_posed_relatives(model, pose_params)

  absolutes = relatives_to_absolutes(relatives, model.parents)

  transforms = zeros(eltype(pose_params),size(absolutes))
  for i=1:size(transforms,1)
    transforms[i,:,:] = squeeze(absolutes[i,:,:],1) *
        squeeze(model.inverse_base_absolutes[i,:,:],1)
  end

  n_verts = size(model.base_positions,2)
  positions = zeros(eltype(transforms),3,n_verts)
  for i=1:size(transforms,1)
    positions = positions +
        (squeeze(transforms[i,1:3,:],1) * model.base_positions) .* model.weights[i,:]
  end

  if model.is_mirrored
    positions[1,:] = -positions[1,:]
  end

  apply_global = true
  if apply_global
    positions = apply_global_transform(pose_params, positions)
  end
  positions
end

function to_pose_params(theta,n_bones)
# to_pose_params !!!!!!!!!!!!!!! fixed order pose_params !!!!!
#       1) global_rotation 2) scale 3) global_translation
#       4) wrist
#       5) thumb1, 6)thumb2, 7) thumb3, 8) thumb4
#       similarly: index, middle, ring, pinky
#       end) forearm

  n = 3 + n_bones
  pose_params = zeros(eltype(theta),3,n)

  pose_params[:,1] = theta[1:3]
  pose_params[:,2] = 1
  pose_params[:,3] = theta[4:6]

  i_theta = 7
  i_pose_params = 6
  n_fingers = 5
  for finger = 1:n_fingers
    for i=2:4
      pose_params[1,i_pose_params] = theta[i_theta]
      i_theta = i_theta + 1
      if i==2
        pose_params[2,i_pose_params] = theta[i_theta]
        i_theta = i_theta + 1
      end
      i_pose_params = i_pose_params+1
    end
    i_pose_params = i_pose_params+1
  end
  pose_params
end

function hand_objective(params, data)
  pose_params = to_pose_params(params, length(data.model.bone_names))

  vertex_positions = get_skinned_vertex_positions(data.model, pose_params)

  n_corr = length(data.correspondences)
  err = zeros(eltype(params),3, n_corr)
  for i=1:n_corr
      err[:,i] = data.points[:,i] - vertex_positions[:,data.correspondences[i]]
  end
  err
end

##################### run it ###########################

# Read instance
dir_in = ARGS[1]
dir_out = ARGS[2]
fn = ARGS[3]
nruns_f = parse(Int,ARGS[4])
nruns_J = parse(Int,ARGS[5])

model_dir = joinpath(dir_in, "model")
fn_in = joinpath(dir_in, fn)
fn_out = string(dir_out, fn)

params,data = read_hand_instance(model_dir,string(fn_in,".txt"))

precompile(euler_angles_to_rotation_matrix,(Vec,))
precompile(get_posed_relatives,(HandModel, Mat))
precompile(relatives_to_absolutes,(Tensor3, iVec))
precompile(angle_axis_to_rotation_matrix,(Vec,))
precompile(apply_global_transform,(Mat, Mat))
precompile(get_skinned_vertex_positions,(HandModel, Mat))
precompile(to_pose_params,(Vec, Int))

err = zeros(Float64,3,size(data.points,2))
tic()
for i in 1:nruns_f
  err = hand_objective(params, data)
end
tf = toq()/nruns_f;
@printf "tf: %g\n" tf
#println(err)

# Gradient helper
function wrapper_hand_objective(params)
  err=hand_objective(params,data)
  err[:]
end

# Gradient
g = jacobian(wrapper_hand_objective)
precompile(g,(Vec,))
J = zeros(3*size(data.points,2),length(params))
tic()
for i in 1:nruns_J
  J = g(params)
end
tJ = toq()/nruns_J;
@printf "tJ: %g\n" tJ
#println("J:")
#println(J)

name = "Julia_F"

write_J(string(fn_out,"_J_",name,".txt"),J)
write_times(string(fn_out,"_times_",name,".txt"),tf,tJ)

max_stars_repo_name: dsyme/ADBench

max_issues_repo_name: dsyme/ADBench

max_forks_repo_name: dsyme/ADBench


text:
function loadColors!(colData, colorMat, shapeStride)
	if size(colorMat, 1) == 1
		colData[1:4:end] = colorMat[1].r
		colData[2:4:end] = colorMat[1].g
		colData[3:4:end] = colorMat[1].b
		colData[4:4:end] = 1.0
	else
		@inbounds @simd for c = 1:size(colorMat, 1)
			colData[(c-1)*shapeStride+1:4:c*shapeStride] = colorMat[c].r
			colData[(c-1)*shapeStride+2:4:c*shapeStride] = colorMat[c].g
			colData[(c-1)*shapeStride+3:4:c*shapeStride] = colorMat[c].b
			colData[(c-1)*shapeStride+4:4:c*shapeStride] = 1.0
		end
	end
end