text:
const MANTISSA_MASK = 0x000fffffffffffff
const EXP_MASK = 0x00000000000007ff

memcpy(d, doff, s, soff, n) = ccall(:memcpy, Cvoid, (Ptr{UInt8}, Ptr{UInt8}, Int), d + doff - 1, s + soff - 1, n)
memmove(d, doff, s, soff, n) = ccall(:memmove, Cvoid, (Ptr{UInt8}, Ptr{UInt8}, Int), d + doff - 1, s + soff - 1, n)

uint(x::Float16) = Core.bitcast(UInt16, x)
uint(x::Float32) = Core.bitcast(UInt32, x)
uint(x::Float64) = Core.bitcast(UInt64, x)

mantissabits(::Type{Float16}) = 10
mantissabits(::Type{Float32}) = 23
mantissabits(::Type{Float64}) = 52

exponentbits(::Type{Float16}) = 5
exponentbits(::Type{Float32}) = 8
exponentbits(::Type{Float64}) = 11

bias(::Type{Float16}) = 15
bias(::Type{Float32}) = 127
bias(::Type{Float64}) = 1023

pow5_bitcount(::Type{Float16}) = 30
pow5_bitcount(::Type{Float32}) = 61
pow5_bitcount(::Type{Float64}) = 121

pow5_inv_bitcount(::Type{Float16}) = 30
pow5_inv_bitcount(::Type{Float32}) = 59
pow5_inv_bitcount(::Type{Float64}) = 122

qinvbound(::Type{Float16}) = 4
qinvbound(::Type{Float32}) = 9
qinvbound(::Type{Float64}) = 21

qbound(::Type{Float16}) = 15
qbound(::Type{Float32}) = 31
qbound(::Type{Float64}) = 63

log10pow2(e) = (e * 78913) >> 18
log10pow5(e) = (e * 732923) >> 20
pow5bits(e) = ((e * 1217359) >> 19) + 1
@inline mulshift(m::UInt64, mula, mulb, j) = ((((UInt128(m) * mula) >> 64) + UInt128(m) * mulb) >> (j - 64)) % UInt64
@inline mulshift(m::UInt32, mul, j) = ((((UInt64(m) * (mul % UInt32)) >> 32) + (UInt64(m) * (mul >> 32))) >> (j - 32)) % UInt32
@inline mulshift(m::UInt16, mul, j) = ((((UInt32(m) * (mul % UInt16)) >> 16) + (UInt32(m) * (mul >> 16))) >> (j - 16))
indexforexp(e) = div(e + 15, 16)
pow10bitsforindex(idx) = 16 * idx + 120
lengthforindex(idx) = div(((Int64(16 * idx) * 1292913986) >> 32) + 1 + 16 + 8, 9)

@inline function pow5(x, p)
    count = 0
    while true
        q = div(x, 5)
        r = x - 5 * q
        r != 0 && return count >= p
        x = q
        count += 1
    end
end

pow2(x, p) = (x & ((Int64(1) << p) - 1)) == 0

@inline function decimallength(v)
    v >= 10000000000000000 && return 17
    v >= 1000000000000000 && return 16
    v >= 100000000000000 && return 15
    v >= 10000000000000 && return 14
    v >= 1000000000000 && return 13
    v >= 100000000000 && return 12
    v >= 10000000000 && return 11
    v >= 1000000000 && return 10
    v >= 100000000 && return 9
    v >= 10000000 && return 8
    v >= 1000000 && return 7
    v >= 100000 && return 6
    v >= 10000 && return 5
    v >= 1000 && return 4
    v >= 100 && return 3
    v >= 10 && return 2
    return 1
end

@inline function decimallength(v::UInt32)
    v >= 100000000 && return 9
    v >= 10000000 && return 8
    v >= 1000000 && return 7
    v >= 100000 && return 6
    v >= 10000 && return 5
    v >= 1000 && return 4
    v >= 100 && return 3
    v >= 10 && return 2
    return 1
end

@inline function decimallength(v::UInt16)
    v >= 10000 && return 5
    v >= 1000 && return 4
    v >= 100 && return 3
    v >= 10 && return 2
    return 1
end

@inline function mulshiftinvsplit(::Type{Float64}, mv, mp, mm, i, j)
    @inbounds mula, mulb = DOUBLE_POW5_INV_SPLIT[i + 1]
    vr = mulshift(mv, mula, mulb, j)
    vp = mulshift(mp, mula, mulb, j)
    vm = mulshift(mm, mula, mulb, j)
    return vr, vp, vm
end

@inline function mulshiftinvsplit(::Type{Float32}, mv, mp, mm, i, j)
    @inbounds mul = FLOAT_POW5_INV_SPLIT[i + 1]
    vr = mulshift(mv, mul, j)
    vp = mulshift(mp, mul, j)
    vm = mulshift(mm, mul, j)
    return vr, vp, vm
end

@inline function mulshiftinvsplit(::Type{Float16}, mv, mp, mm, i, j)
    @inbounds mul = HALF_POW5_INV_SPLIT[i + 1]
    vr = mulshift(mv, mul, j)
    vp = mulshift(mp, mul, j)
    vm = mulshift(mm, mul, j)
    return vr, vp, vm
end

@inline function mulshiftsplit(::Type{Float64}, mv, mp, mm, i, j)
    @inbounds mula, mulb = DOUBLE_POW5_SPLIT[i + 1]
    vr = mulshift(mv, mula, mulb, j)
    vp = mulshift(mp, mula, mulb, j)
    vm = mulshift(mm, mula, mulb, j)
    return vr, vp, vm
end

@inline function mulshiftsplit(::Type{Float32}, mv, mp, mm, i, j)
    @inbounds mul = FLOAT_POW5_SPLIT[i + 1]
    vr = mulshift(mv, mul, j)
    vp = mulshift(mp, mul, j)
    vm = mulshift(mm, mul, j)
    return vr, vp, vm
end

@inline function mulshiftsplit(::Type{Float16}, mv, mp, mm, i, j)
    @inbounds mul = HALF_POW5_SPLIT[i + 1]
    vr = mulshift(mv, mul, j)
    vp = mulshift(mp, mul, j)
    vm = mulshift(mm, mul, j)
    return vr, vp, vm
end

@inline function pow5invsplit(::Type{Float64}, i)
    pow = big(5)^i
    inv = div(big(1) << (bitlength(pow) - 1 + pow5_inv_bitcount(Float64)), pow) + 1
    return (UInt64(inv & ((big(1) << 64) - 1)), UInt64(inv >> 64))
end

@inline function pow5invsplit(::Type{Float32}, i)
    pow = big(5)^i
    inv = div(big(1) << (bitlength(pow) - 1 + pow5_inv_bitcount(Float32)), pow) + 1
    return UInt64(inv)
end

@inline function pow5invsplit(::Type{Float16}, i)
    pow = big(5)^i
    inv = div(big(1) << (bitlength(pow) - 1 + pow5_inv_bitcount(Float16)), pow) + 1
    return UInt32(inv)
end

@inline function pow5split(::Type{Float64}, i)
    pow = big(5)^i
    j = bitlength(pow) - pow5_bitcount(Float64)
    return (UInt64((pow >> j) & ((big(1) << 64) - 1)), UInt64(pow >> (j + 64)))
end

@inline function pow5split(::Type{Float32}, i)
    pow = big(5)^i
    return UInt64(pow >> (bitlength(pow) - pow5_bitcount(Float32)))
end

@inline function pow5split(::Type{Float16}, i)
    pow = big(5)^i
    return UInt32(pow >> (bitlength(pow) - pow5_bitcount(Float16)))
end

const DOUBLE_POW5_INV_SPLIT = map(i->pow5invsplit(Float64, i), 0:291)
const FLOAT_POW5_INV_SPLIT = map(i->pow5invsplit(Float32, i), 0:30)
const HALF_POW5_INV_SPLIT = map(i->pow5invsplit(Float16, i), 0:17)

const DOUBLE_POW5_SPLIT = map(i->pow5split(Float64, i), 0:325)
const FLOAT_POW5_SPLIT = map(i->pow5split(Float32, i), 0:46)
const HALF_POW5_SPLIT = map(i->pow5split(Float16, i), 0:23)

const DIGIT_TABLE = UInt8[
  '0','0','0','1','0','2','0','3','0','4','0','5','0','6','0','7','0','8','0','9',
  '1','0','1','1','1','2','1','3','1','4','1','5','1','6','1','7','1','8','1','9',
  '2','0','2','1','2','2','2','3','2','4','2','5','2','6','2','7','2','8','2','9',
  '3','0','3','1','3','2','3','3','3','4','3','5','3','6','3','7','3','8','3','9',
  '4','0','4','1','4','2','4','3','4','4','4','5','4','6','4','7','4','8','4','9',
  '5','0','5','1','5','2','5','3','5','4','5','5','5','6','5','7','5','8','5','9',
  '6','0','6','1','6','2','6','3','6','4','6','5','6','6','6','7','6','8','6','9',
  '7','0','7','1','7','2','7','3','7','4','7','5','7','6','7','7','7','8','7','9',
  '8','0','8','1','8','2','8','3','8','4','8','5','8','6','8','7','8','8','8','9',
  '9','0','9','1','9','2','9','3','9','4','9','5','9','6','9','7','9','8','9','9'
]

neededdigits(::Type{Float64}) = 309 + 17
neededdigits(::Type{Float32}) = 39 + 9
neededdigits(::Type{Float16}) = 9 + 5

@inline function writeshortest(buf::Vector{UInt8}, pos, x::T,
    plus=false, space=false, hash=true,
    precision=-1, expchar=UInt8('e'), padexp=false, decchar=UInt8('.')) where {T}
    neg = signbit(x)
    # special cases
    if x == 0
        if neg
            buf[pos] = UInt8('-')
            pos += 1
        elseif plus
            buf[pos] = UInt8('+')
            pos += 1
        elseif space
            buf[pos] = UInt8(' ')
            pos += 1
        end
        buf[pos] = UInt8('0')
        pos += 1
        if hash
            buf[pos] = decchar
            pos += 1
        end
        if precision == -1
            buf[pos] = UInt8('0')
            return pos + 1
        end
        while precision > 1
            buf[pos] = UInt8('0')
            pos += 1
            precision -= 1
        end
        return pos
    elseif isnan(x)
        buf[pos] = UInt8('N')
        buf[pos + 1] = UInt8('a')
        buf[pos + 2] = UInt8('N')
        return pos + 3
    elseif !isfinite(x)
        if neg
            buf[pos] = UInt8('-')
        end
        buf[pos + neg] = UInt8('I')
        buf[pos + neg + 1] = UInt8('n')
        buf[pos + neg + 2] = UInt8('f')
        return pos + neg + 3
    end

    bits = uint(x)
    mant = bits & (oftype(bits, 1) << mantissabits(T) - oftype(bits, 1))
    exp = Int((bits >> mantissabits(T)) & ((Int64(1) << exponentbits(T)) - 1))
    m2 = oftype(bits, Int64(1) << mantissabits(T)) | mant
    e2 = exp - bias(T) - mantissabits(T)
    fraction = m2 & ((oftype(bits, 1) << -e2) - 1)
    if e2 > 0 || e2 < -52 || fraction != 0
        if exp == 0
            e2 = 1 - bias(T) - mantissabits(T) - 2
            m2 = mant
        else
            e2 -= 2
        end
        even = (m2 & 1) == 0
        mv = oftype(m2, 4 * m2)
        mp = oftype(m2, mv + 2)
        mmShift = mant != 0 || exp <= 1
        mm = oftype(m2, mv - 1 - mmShift)
        vmIsTrailingZeros = false
        vrIsTrailingZeros = false
        lastRemovedDigit = 0x00
        if e2 >= 0
            q = log10pow2(e2) - (T == Float64 ? (e2 > 3) : 0)
            e10 = q
            k = pow5_inv_bitcount(T) + pow5bits(q) - 1
            i = -e2 + q + k
            vr, vp, vm = mulshiftinvsplit(T, mv, mp, mm, q, i)
            if T == Float32 || T == Float16
                if q != 0 && div(vp - 1, 10) <= div(vm, 10)
                    l = pow5_inv_bitcount(T) + pow5bits(q - 1) - 1
                    mul = T == Float32 ? FLOAT_POW5_INV_SPLIT[q] : HALF_POW5_INV_SPLIT[q]
                    lastRemovedDigit = (mulshift(mv, mul, -e2 + q - 1 + l) % 10) % UInt8
                end
            end
            if q <= qinvbound(T)
                if ((mv % UInt32) - 5 * div(mv, 5)) == 0
                    vrIsTrailingZeros = pow5(mv, q)
                elseif even
                    vmIsTrailingZeros = pow5(mm, q)
                else
                    vp -= pow5(mp, q)
                end
            end
        else
            q = log10pow5(-e2) - (T == Float64 ? (-e2 > 1) : 0)
            e10 = q + e2
            i = -e2 - q
            k = pow5bits(i) - pow5_bitcount(T)
            j = q - k
            vr, vp, vm = mulshiftsplit(T, mv, mp, mm, i, j)
            if T == Float32 || T == Float16
                if q != 0 && div(vp - 1, 10) <= div(vm, 10)
                    j = q - 1 - (pow5bits(i + 1) - pow5_bitcount(T))
                    mul = T == Float32 ? FLOAT_POW5_SPLIT[i + 2] : HALF_POW5_SPLIT[i + 2]
                    lastRemovedDigit = (mulshift(mv, mul, j) % 10) % UInt8
                end
            end
            if q <= 1
                vrIsTrailingZeros = true
                if even
                    vmIsTrailingZeros = mmShift
                else
                    vp -= 1
                end
            elseif q < qbound(T)
                vrIsTrailingZeros = pow2(mv, q - (T != Float64))
            end
        end
        removed = 0
        if vmIsTrailingZeros || vrIsTrailingZeros
            while true
                vpDiv10 = div(vp, 10)
                vmDiv10 = div(vm, 10)
                vpDiv10 <= vmDiv10 && break
                vmMod10 = (vm % UInt32) - UInt32(10) * (vmDiv10 % UInt32)
                vrDiv10 = div(vr, 10)
                vrMod10 = (vr % UInt32) - UInt32(10) * (vrDiv10 % UInt32)
                vmIsTrailingZeros &= vmMod10 == 0
                vrIsTrailingZeros &= lastRemovedDigit == 0
                lastRemovedDigit = vrMod10 % UInt8
                vr = vrDiv10
                vp = vpDiv10
                vm = vmDiv10
                removed += 1
            end
            if vmIsTrailingZeros
                while true
                    vmDiv10 = div(vm, 10)
                    vmMod10 = (vm % UInt32) - UInt32(10) * (vmDiv10 % UInt32)
                    vmMod10 != 0 && break
                    vpDiv10 = div(vp, 10)
                    vrDiv10 = div(vr, 10)
                    vrMod10 = (vr % UInt32) - UInt32(10) * (vrDiv10 % UInt32)
                    vrIsTrailingZeros &= lastRemovedDigit == 0
                    lastRemovedDigit = vrMod10 % UInt8
                    vr = vrDiv10
                    vp = vpDiv10
                    vm = vmDiv10
                    removed += 1
                end
            end
            if vrIsTrailingZeros && lastRemovedDigit == 5 && vr % 2 == 0
                lastRemovedDigit = UInt8(4)
            end
            output = vr + ((vr == vm && (!even || !vmIsTrailingZeros)) || lastRemovedDigit >= 5)
        else
            roundUp = false
            vpDiv100 = div(vp, 100)
            vmDiv100 = div(vm, 100)
            if vpDiv100 > vmDiv100
                vrDiv100 = div(vr, 100)
                vrMod100 = (vr % UInt32) - UInt32(100) * (vrDiv100 % UInt32)
                roundUp = vrMod100 >= 50
                vr = vrDiv100
                vp = vpDiv100
                vm = vmDiv100
                removed += 2
            end
            while true
                vpDiv10 = div(vp, 10)
                vmDiv10 = div(vm, 10)
                vpDiv10 <= vmDiv10 && break
                vrDiv10 = div(vr, 10)
                vrMod10 = (vr % UInt32) - UInt32(10) * (vrDiv10 % UInt32)
                roundUp = vrMod10 >= 5
                vr = vrDiv10
                vp = vpDiv10
                vm = vmDiv10
                removed += 1
            end
            output = vr + (vr == vm || roundUp || lastRemovedDigit >= 5)
        end
        nexp = e10 + removed
    else
        output = m2 >> -e2
        nexp = 0
        while true
            q = div(output, 10)
            r = (output % UInt32) - UInt32(10) * (q % UInt32)
            r != 0 && break
            output = q
            nexp += 1
        end
    end

    if neg
        buf[pos] = UInt8('-')
        pos += 1
    elseif plus
        buf[pos] = UInt8('+')
        pos += 1
    elseif space
        buf[pos] = UInt8(' ')
        pos += 1
    end

    olength = decimallength(output)
    exp_form = true
    pt = nexp + olength
    if -4 < pt <= (precision == -1 ? (T == Float16 ? 5 : 6) : precision) #&& !(pt >= olength && abs(mod(x + 0.05, 10^(pt - olength)) - 0.05) > 0.05)
        exp_form = false
        if pt <= 0
            buf[pos] = UInt8('0')
            pos += 1
            buf[pos] = decchar
            pos += 1
            for _ = 1:abs(pt)
                buf[pos] = UInt8('0')
                pos += 1
            end
        # elseif pt >= olength
            # nothing to do at this point
        # else
            # nothing to do at this point
        end
    else
        pos += 1
    end
    i = 0
    ptr = pointer(buf)
    ptr2 = pointer(DIGIT_TABLE)
    if (output >> 32) != 0
        q = output ÷ 100000000
        output2 = (output % UInt32) - UInt32(100000000) * (q % UInt32)
        output = q

        c = output2 % UInt32(10000)
        output2 = div(output2, UInt32(10000))
        d = output2 % UInt32(10000)
        c0 = (c % 100) << 1
        c1 = (c ÷ 100) << 1
        d0 = (d % 100) << 1
        d1 = (d ÷ 100) << 1
        memcpy(ptr, pos + olength - 2, ptr2, c0 + 1, 2)
        memcpy(ptr, pos + olength - 4, ptr2, c1 + 1, 2)
        memcpy(ptr, pos + olength - 6, ptr2, d0 + 1, 2)
        memcpy(ptr, pos + olength - 8, ptr2, d1 + 1, 2)
        i += 8
    end
    output2 = output % UInt32
    while output2 >= 10000
        c = output2 % UInt32(10000)
        output2 = div(output2, UInt32(10000))
        c0 = (c % 100) << 1
        c1 = (c ÷ 100) << 1
        memcpy(ptr, pos + olength - i - 2, ptr2, c0 + 1, 2)
        memcpy(ptr, pos + olength - i - 4, ptr2, c1 + 1, 2)
        i += 4
    end
    if output2 >= 100
        c = (output2 % UInt32(100)) << 1
        output2 = div(output2, UInt32(100))
        memcpy(ptr, pos + olength - i - 2, ptr2, c + 1, 2)
        i += 2
    end
    if output2 >= 10
        c = output2 << 1
        buf[pos + 1] = DIGIT_TABLE[c + 2]
        buf[pos - exp_form] = DIGIT_TABLE[c + 1]
    else
        buf[pos - exp_form] = UInt8('0') + (output2 % UInt8)
    end

    if !exp_form
        if pt <= 0
            pos += olength
            precision -= olength
            while hash && precision > 0
                buf[pos] = UInt8('0')
                pos += 1
                precision -= 1
            end
        elseif pt >= olength
            pos += olength
            precision -= olength
            for _ = 1:nexp
                buf[pos] = UInt8('0')
                pos += 1
                precision -= 1
            end
            if hash
                buf[pos] = decchar
                pos += 1
                if precision < 0
                    buf[pos] = UInt8('0')
                    pos += 1
                end
                while precision > 0
                    buf[pos] = UInt8('0')
                    pos += 1
                    precision -= 1
                end
            end
        else
            pointoff = olength - abs(nexp)
            memmove(ptr, pos + pointoff + 1, ptr, pos + pointoff, olength - pointoff + 1)
            buf[pos + pointoff] = decchar
            pos += olength + 1
            precision -= olength
            while hash && precision > 0
                buf[pos] = UInt8('0')
                pos += 1
                precision -= 1
            end
        end
    else
        if olength > 1 || hash
            buf[pos] = decchar
            pos += olength
            precision -= olength
        end
        if hash && olength == 1
            buf[pos] = UInt8('0')
            pos += 1
        end
        while hash && precision > 0
            buf[pos] = UInt8('0')
            pos += 1
            precision -= 1
        end

        buf[pos] = expchar
        pos += 1
        exp2 = nexp + olength - 1
        if exp2 < 0
            buf[pos] = UInt8('-')
            pos += 1
            exp2 = -exp2
        elseif padexp
            buf[pos] = UInt8('+')
            pos += 1
        end

        if exp2 >= 100
            c = exp2 % 10
            memcpy(ptr, pos, ptr2, 2 * div(exp2, 10) + 1, 2)
            buf[pos + 2] = UInt8('0') + (c % UInt8)
            pos += 3
        elseif exp2 >= 10
            memcpy(ptr, pos, ptr2, 2 * exp2 + 1, 2)
            pos += 2
        else
            if padexp
                buf[pos] = UInt8('0')
                pos += 1
            end
            buf[pos] = UInt8('0') + (exp2 % UInt8)
            pos += 1
        end
    end

    return pos
end

function writeshortest(x::T,
        plus::Bool=false,
        space::Bool=false,
        hash::Bool=true,
        precision::Integer=-1,
        expchar::UInt8=UInt8('e'),
        padexp::Bool=false,
        decchar::UInt8=UInt8('.')) where {T <: Base.IEEEFloat}
    buf = Base.StringVector(neededdigits(T))
    pos = writeshortest(buf, 1, x)
    @assert pos - 1 <= length(buf)
    return String(resize!(buf, pos - 1))
end


text:
filepath = Base.source_path()

using Pkg
Pkg.activate( joinpath(dirname(filepath),"..") )

import InfrastructureSystems
const IS = InfrastructureSystems

function main(args)
    if length(args) != 2
        println("CLI usage: julia generate_structs.jl INPUT_FILE OUTPUT_DIRECTORY")
        println("REPL usage: IS.generate_structs(INPUT_FILE, OUTPUT_DIRECTORY)")
        return
    end

    IS.generate_structs(args[1], args[2])
end

main(ARGS)

ext: jl


text:
# https://plot.ly/javascript/getting-started

supported_args(::PlotlyBackend) = merge_with_base_supported([
    :annotations,
    :background_color_legend, :background_color_inside, :background_color_outside,
    :foreground_color_legend, :foreground_color_guide,
    # :foreground_color_grid, :foreground_color_axis,
    :foreground_color_text, :foreground_color_border,
    :foreground_color_title,
    :label,
    :seriescolor, :seriesalpha,
    :linecolor, :linestyle, :linewidth, :linealpha,
    :markershape, :markercolor, :markersize, :markeralpha,
    :markerstrokewidth, :markerstrokecolor, :markerstrokealpha, :markerstrokestyle,
    :fillrange, :fillcolor, :fillalpha,
    :bins,
    :title, :title_location, :titlefont,
    :window_title,
    :guide, :lims, :ticks, :scale, :flip, :rotation,
    :tickfont, :guidefont, :legendfont,
    :grid, :legend, :colorbar,
    :marker_z, :levels,
    :ribbon, :quiver,
    :orientation,
    # :overwrite_figure,
    :polar,
    :normalize, :weights,
    # :contours, :aspect_ratio,
    :hover,
    :inset_subplots,
  ])

supported_types(::PlotlyBackend) = [
    :path, :scatter, :bar, :pie, :heatmap,
    :contour, :surface, :path3d, :scatter3d, :shape, :scattergl,
]
supported_styles(::PlotlyBackend) = [:auto, :solid, :dash, :dot, :dashdot]
supported_markers(::PlotlyBackend) = [
    :none, :auto, :circle, :rect, :diamond, :utriangle, :dtriangle,
    :cross, :xcross, :pentagon, :hexagon, :octagon, :vline, :hline
]
supported_scales(::PlotlyBackend) = [:identity, :log10]
is_subplot_supported(::PlotlyBackend) = true
is_string_supported(::PlotlyBackend) = true

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

function _initialize_backend(::PlotlyBackend; kw...)
  @eval begin
    import JSON
    JSON._print(io::IO, state::JSON.State, dt::Union{Date,DateTime}) = print(io, '"', dt, '"')

    _js_path = Pkg.dir("Plots", "deps", "plotly-latest.min.js")
    _js_code = open(@compat(readstring), _js_path, "r")

    # borrowed from https://github.com/plotly/plotly.py/blob/2594076e29584ede2d09f2aa40a8a195b3f3fc66/plotly/offline/offline.py#L64-L71 c/o @spencerlyon2
    _js_script = """

    """

    # if we're in IJulia call setupnotebook to load js and css
    if isijulia()
        display("text/html", _js_script)
    end

    # if isatom()
    #     import Atom
    #     Atom.@msg evaljs(_js_code)
    # end

  end
  # TODO: other initialization
end

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

function plotly_font(font::Font, color = font.color)
    KW(
        :family => font.family,
        :size   => round(Int, font.pointsize*1.4),
        :color  => webcolor(color),
    )
end

function plotly_annotation_dict(x, y, val; xref="paper", yref="paper")
    KW(
        :text => val,
        :xref => xref,
        :x => x,
        :yref => xref,
        :y => y,
        :showarrow => false,
    )
end

function plotly_annotation_dict(x, y, ptxt::PlotText; xref="paper", yref="paper")
    merge(plotly_annotation_dict(x, y, ptxt.str; xref=xref, yref=yref), KW(
        :font => plotly_font(ptxt.font),
        :xanchor => ptxt.font.halign == :hcenter ? :center : ptxt.font.halign,
        :yanchor => ptxt.font.valign == :vcenter ? :middle : ptxt.font.valign,
        :rotation => ptxt.font.rotation,
    ))
end

# function get_annotation_dict_for_arrow(d::KW, xyprev::Tuple, xy::Tuple, a::Arrow)
#     xdiff = xyprev[1] - xy[1]
#     ydiff = xyprev[2] - xy[2]
#     dist = sqrt(xdiff^2 + ydiff^2)
#     KW(
#         :showarrow => true,
#         :x => xy[1],
#         :y => xy[2],
#         # :ax => xyprev[1] - xy[1],
#         # :ay => xy[2] - xyprev[2],
#         # :ax => 0,
#         # :ay => -40,
#         :ax => 10xdiff / dist,
#         :ay => -10ydiff / dist,
#         :arrowcolor => webcolor(d[:linecolor], d[:linealpha]),
#         :xref => "x",
#         :yref => "y",
#         :arrowsize => 10a.headwidth,
#         # :arrowwidth => a.headlength,
#         :arrowwidth => 0.1,
#     )
# end

function plotly_scale(scale::Symbol)
    if scale == :log10
        "log"
    else
        "-"
    end
end

# this method gets the start/end in percentage of the canvas for this axis direction
function plotly_domain(sp::Subplot, letter)
    figw, figh = sp.plt[:size]
    pcts = bbox_to_pcts(sp.plotarea, figw*px, figh*px)
    i1,i2 = (letter == :x ? (1,3) : (2,4))
    [pcts[i1], pcts[i1]+pcts[i2]]
end

function plotly_axis(axis::Axis, sp::Subplot)
    letter = axis[:letter]
    ax = KW(
        :title      => axis[:guide],
        :showgrid   => sp[:grid],
        :zeroline   => false,
    )

    if letter in (:x,:y)
        ax[:domain] = plotly_domain(sp, letter)
        ax[:anchor] = "$(letter==:x ? :y : :x)$(plotly_subplot_index(sp))"
    end

    rot = axis[:rotation]
    if rot != 0
        ax[:tickangle] = rot
    end

    if !(axis[:ticks] in (nothing, :none))
        ax[:titlefont] = plotly_font(axis[:guidefont], axis[:foreground_color_guide])
        ax[:type] = plotly_scale(axis[:scale])
        ax[:tickfont] = plotly_font(axis[:tickfont], axis[:foreground_color_text])
        ax[:tickcolor] = webcolor(axis[:foreground_color_border])
        ax[:linecolor] = webcolor(axis[:foreground_color_border])

        # lims
        lims = axis[:lims]
        if lims != :auto && limsType(lims) == :limits
            ax[:range] = lims
        end

        # flip
        if axis[:flip]
            ax[:autorange] = "reversed"
        end

        # ticks
        ticks = axis[:ticks]
        if ticks != :auto
            ttype = ticksType(ticks)
            if ttype == :ticks
                ax[:tickmode] = "array"
                ax[:tickvals] = ticks
            elseif ttype == :ticks_and_labels
                ax[:tickmode] = "array"
                ax[:tickvals], ax[:ticktext] = ticks
            end
        end
    else
        ax[:showticklabels] = false
        ax[:showgrid] = false
    end

    ax
end

function plotly_layout(plt::Plot)
    d_out = KW()

    w, h = plt[:size]
    d_out[:width], d_out[:height] = w, h
    d_out[:paper_bgcolor] = webcolor(plt[:background_color_outside])
    d_out[:margin] = KW(:l=>0, :b=>0, :r=>0, :t=>20)

    d_out[:annotations] = KW[]

    for sp in plt.subplots
        spidx = plotly_subplot_index(sp)

        # add an annotation for the title... positioned horizontally relative to plotarea,
        # but vertically just below the top of the subplot bounding box
        if sp[:title] != ""
            bb = plotarea(sp)
            tpos = sp[:title_location]
            xmm = if tpos == :left
                left(bb)
            elseif tpos == :right
                right(bb)
            else
                0.5 * (left(bb) + right(bb))
            end
            titlex, titley = xy_mm_to_pcts(xmm, top(bbox(sp)), w*px, h*px)
            titlefont = font(sp[:titlefont], :top, sp[:foreground_color_title])
            push!(d_out[:annotations], plotly_annotation_dict(titlex, titley, text(sp[:title], titlefont)))
        end

        d_out[:plot_bgcolor] = webcolor(sp[:background_color_inside])

        # TODO: x/y axis tick values/labels

        # if any(is3d, seriesargs)
        if is3d(sp)
            d_out[:scene] = KW(
                Symbol("xaxis$spidx") => plotly_axis(sp[:xaxis], sp),
                Symbol("yaxis$spidx") => plotly_axis(sp[:yaxis], sp),
                Symbol("zaxis$spidx") => plotly_axis(sp[:zaxis], sp),
            )
        else
            d_out[Symbol("xaxis$spidx")] = plotly_axis(sp[:xaxis], sp)
            d_out[Symbol("yaxis$spidx")] = plotly_axis(sp[:yaxis], sp)
        end

        # legend
        d_out[:showlegend] = sp[:legend] != :none
        if sp[:legend] != :none
            d_out[:legend] = KW(
                :bgcolor  => webcolor(sp[:background_color_legend]),
                :bordercolor => webcolor(sp[:foreground_color_legend]),
                :font     => plotly_font(sp[:legendfont], sp[:foreground_color_legend]),
            )
        end

        # annotations
        append!(d_out[:annotations], KW[plotly_annotation_dict(ann...; xref = "x$spidx", yref = "y$spidx") for ann in sp[:annotations]])

        # # arrows
        # for sargs in seriesargs
        #     a = sargs[:arrow]
        #     if sargs[:seriestype] in (:path, :line) && typeof(a) <: Arrow
        #         add_arrows(sargs[:x], sargs[:y]) do xyprev, xy
        #             push!(d_out[:annotations], get_annotation_dict_for_arrow(sargs, xyprev, xy, a))
        #         end
        #     end
        # end

        if ispolar(sp)
            d_out[:direction] = "counterclockwise"
        end

        d_out
    end

    # turn off hover if nothing's using it
    if all(series -> series.d[:hover] in (false,:none), plt.series_list)
        d_out[:hovermode] = "none"
    end

    d_out
end

function plotly_layout_json(plt::Plot)
    JSON.json(plotly_layout(plt))
end

function plotly_colorscale(grad::ColorGradient, alpha = nothing)
    [[grad.values[i], webcolor(grad.colors[i], alpha)] for i in 1:length(grad.colors)]
end
plotly_colorscale(c, alpha = nothing) = plotly_colorscale(default_gradient(), alpha)

const _plotly_markers = KW(
    :rect       => "square",
    :xcross     => "x",
    :utriangle  => "triangle-up",
    :dtriangle  => "triangle-down",
    :star5      => "star-triangle-up",
    :vline      => "line-ns",
    :hline      => "line-ew",
)

function plotly_subplot_index(sp::Subplot)
    spidx = sp[:subplot_index]
    spidx == 1 ? "" : spidx
end

# the Shape contructor will automatically close the shape. since we need it closed,
# we split by NaNs and then construct/destruct the shapes to get the closed coords
function plotly_close_shapes(x, y)
    xs, ys = nansplit(x), nansplit(y)
    for i=1:length(xs)
        shape = Shape(xs[i], ys[i])
        xs[i], ys[i] = shape_coords(shape)
    end
    nanvcat(xs), nanvcat(ys)
end

# get a dictionary representing the series params (d is the Plots-dict, d_out is the Plotly-dict)
function plotly_series(plt::Plot, series::Series)
    d = series.d
    sp = d[:subplot]
    d_out = KW()

    # these are the axes that the series should be mapped to
    spidx = plotly_subplot_index(sp)
    d_out[:xaxis] = "x$spidx"
    d_out[:yaxis] = "y$spidx"
    d_out[:showlegend] = should_add_to_legend(series)

    x, y = collect(d[:x]), collect(d[:y])
    d_out[:name] = d[:label]
    st = d[:seriestype]
    isscatter = st in (:scatter, :scatter3d, :scattergl)
    hasmarker = isscatter || d[:markershape] != :none
    # hasline = !isscatter
    hasline = st in (:path, :path3d)

    # set the "type"
    if st in (:path, :scatter, :scattergl)
        d_out[:type] = st==:scattergl ? "scattergl" : "scatter"
        d_out[:mode] = if hasmarker
            hasline ? "lines+markers" : "markers"
        else
            hasline ? "lines" : "none"
        end
        if d[:fillrange] == true || d[:fillrange] == 0
            d_out[:fill] = "tozeroy"
            d_out[:fillcolor] = webcolor(d[:fillcolor], d[:fillalpha])
        elseif !(d[:fillrange] in (false, nothing))
            warn("fillrange ignored... plotly only supports filling to zero. fillrange: $(d[:fillrange])")
        end
        d_out[:x], d_out[:y] = x, y

    elseif st == :shape
        # to draw polygons, we actually draw lines with fill
        d_out[:type] = "scatter"
        d_out[:mode] = "lines"
        d_out[:x], d_out[:y] = plotly_close_shapes(x, y)
        # @show map(length, (x,y,d_out[:x],d_out[:y]))
        # @show d_out[:x] d_out[:y]
        d_out[:fill] = "tozeroy"
        d_out[:fillcolor] = webcolor(d[:markercolor], d[:markeralpha])
        if d[:markerstrokewidth] > 0
            d_out[:line] = KW(
                :color => webcolor(d[:markerstrokecolor], d[:markerstrokealpha]),
                :width => d[:markerstrokewidth],
                :dash => string(d[:markerstrokestyle]),
            )
        end

    elseif st == :bar
        d_out[:type] = "bar"
        d_out[:x], d_out[:y] = x, y
        d_out[:orientation] = isvertical(d) ? "v" : "h"

    # elseif st == :histogram2d
    #     d_out[:type] = "histogram2d"
    #     d_out[:x], d_out[:y] = x, y
    #     if isa(d[:bins], Tuple)
    #         xbins, ybins = d[:bins]
    #     else
    #         xbins = ybins = d[:bins]
    #     end
    #     d_out[:nbinsx] = xbins
    #     d_out[:nbinsy] = ybins
    #     d_out[:colorscale] = plotly_colorscale(d[:fillcolor], d[:fillalpha])

    # elseif st in (:histogram, :density)
    #     d_out[:type] = "histogram"
    #     isvert = isvertical(d)
    #     d_out[isvert ? :x : :y] = y
    #     d_out[isvert ? :nbinsx : :nbinsy] = d[:bins]
    #     if st == :density
    #         d_out[:histogramnorm] = "probability density"
    #     end

    elseif st == :heatmap
        d_out[:type] = "heatmap"
        d_out[:x], d_out[:y], d_out[:z] = d[:x], d[:y], transpose_z(d, d[:z].surf, false)
        d_out[:colorscale] = plotly_colorscale(d[:fillcolor], d[:fillalpha])

    elseif st == :contour
        d_out[:type] = "contour"
        d_out[:x], d_out[:y], d_out[:z] = d[:x], d[:y], transpose_z(d, d[:z].surf, false)
        # d_out[:showscale] = d[:colorbar] != :none
        d_out[:ncontours] = d[:levels]
        d_out[:contours] = KW(:coloring => d[:fillrange] != nothing ? "fill" : "lines")
        d_out[:colorscale] = plotly_colorscale(d[:linecolor], d[:linealpha])

    elseif st in (:surface, :wireframe)
        d_out[:type] = "surface"
        d_out[:x], d_out[:y], d_out[:z] = d[:x], d[:y], transpose_z(d, d[:z].surf, false)
        d_out[:colorscale] = plotly_colorscale(d[:fillcolor], d[:fillalpha])

    elseif st == :pie
        d_out[:type] = "pie"
        d_out[:labels] = pie_labels(sp, series)
        d_out[:values] = y
        d_out[:hoverinfo] = "label+percent+name"

    elseif st in (:path3d, :scatter3d)
        d_out[:type] = "scatter3d"
        d_out[:mode] = if hasmarker
            hasline ? "lines+markers" : "markers"
        else
            hasline ? "lines" : "none"
        end
        d_out[:x], d_out[:y] = x, y
        d_out[:z] = collect(d[:z])

    else
        warn("Plotly: seriestype $st isn't supported.")
        return KW()
    end

    # add "marker"
    if hasmarker
        d_out[:marker] = KW(
            :symbol => get(_plotly_markers, d[:markershape], string(d[:markershape])),
            :opacity => d[:markeralpha],
            :size => 2 * d[:markersize],
            :color => webcolor(d[:markercolor], d[:markeralpha]),
            :line => KW(
                :color => webcolor(d[:markerstrokecolor], d[:markerstrokealpha]),
                :width => d[:markerstrokewidth],
            ),
        )

        # gotta hack this (for now?) since plotly can't handle rgba values inside the gradient
        if d[:marker_z] != nothing
            # d_out[:marker][:color] = d[:marker_z]
            # d_out[:marker][:colorscale] = plotly_colorscale(d[:markercolor], d[:markeralpha])
            # d_out[:showscale] = true
            grad = ColorGradient(d[:markercolor], alpha=d[:markeralpha])
            zmin, zmax = extrema(d[:marker_z])
            d_out[:marker][:color] = [webcolor(getColorZ(grad, (zi - zmin) / (zmax - zmin))) for zi in d[:marker_z]]
        end
    end

    # add "line"
    if hasline
        d_out[:line] = KW(
            :color => webcolor(d[:linecolor], d[:linealpha]),
            :width => d[:linewidth],
            :shape => if st == :steppre
                "vh"
            elseif st == :steppost
                "hv"
            else
                "linear"
            end,
            :dash => string(d[:linestyle]),
            # :dash => "solid",
        )
    end

    # convert polar plots x/y to theta/radius
    if ispolar(d[:subplot])
        d_out[:t] = rad2deg(pop!(d_out, :x))
        d_out[:r] = pop!(d_out, :y)
    end

    # hover text
    hover = d[:hover]
    if hover in (:none, false)
        d_out[:hoverinfo] = "none"
    elseif hover != nothing
        d_out[:hoverinfo] = "text"
        d_out[:text] = hover
    end

    d_out
end

# get a list of dictionaries, each representing the series params
function plotly_series_json(plt::Plot)
    JSON.json(map(series -> plotly_series(plt, series), plt.series_list))
end

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

function html_head(plt::Plot{PlotlyBackend})
    ""
end

function html_body(plt::Plot{PlotlyBackend}, style = nothing)
    if style == nothing
        w, h = plt[:size]
        style = "width:$(w)px;height:$(h)px;"
    end
    uuid = Base.Random.uuid4()
    html = """

    """
    html
end

function js_body(plt::Plot{PlotlyBackend}, uuid)
    js = """
          PLOT = document.getElementById('$(uuid)');
          Plotly.plot(PLOT, $(plotly_series_json(plt)), $(plotly_layout_json(plt)));
    """
end

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

function _writemime(io::IO, ::MIME"image/png", plt::Plot{PlotlyBackend})
    writemime_png_from_html(io, plt)
end

function _writemime(io::IO, ::MIME"image/svg+xml", plt::Plot{PlotlyBackend})
    write(io, html_head(plt) * html_body(plt))
end

function _display(plt::Plot{PlotlyBackend})
    standalone_html_window(plt)
end


text:
using PyPlot
include("MCMC.jl")

function double_well(sparam, θ::Array{Float64,1})
    return [θ[1]^2]
end

function d_double_well(sparam, θ::Array{Float64,1})
    return [θ[1]^2], reshape([2θ[1]],1,1)
end

# function mixture(sparam, θ::Array{Float64,1})
#     if sparam.cond == 1
#         return [0.9 * θ[1]]
#     elseif sparam.cond == 2
#         return [0.1 * θ[1]^3]
#     else
#         error("cond : ", cond, " is not defined")
#     end
# end

struct Setup_Param
    Ny::Int64
    Nθ::Int64
end
s_param = Setup_Param(1, 1)

y = [1.0]
Σ_η = reshape([0.1^2], 1, 1)
μ0 = [2.0]
Σ0 = reshape([5.0^2], 1, 1)
ln_posterior(θ) = log_bayesian_posterior(s_param, θ, double_well,
    y,  Σ_η, μ0, Σ0)
d_ln_posterior(θ) = d_log_bayesian_posterior(s_param, θ, d_double_well,
y,  Σ_η, μ0, Σ0)

# "MALA"  "RWMCMC"  "SVGD"  "emcee"
method = "emcee"
if method == "RWMCMC"
    θ0 = [-0.3]
    step_length, n_ite = 0.2, 10^6
    θs = RWMCMC_Run(ln_posterior, θ0, step_length, n_ite)
elseif method == "MALA"
    θ0 = [-0.3]
    step_length, n_ite = 0.002, 10^6
    θs = MALA_Run(d_ln_posterior, θ0, step_length, n_ite)
elseif method == "SVGD"
    J, n_ite = 100, 1000
    θs0 = randn(J,1)
    θs = SVGD_Run(θs0, d_ln_posterior, n_ite)
elseif method == "emcee"
    J, n_ite = 100, 100
    θs0 = randn(J,1)
    θs = emcee_Run(ln_posterior, θs0, n_ite)
    θs = reshape(θs, (J*n_ite, 1))
end

# Plot
PyPlot.figure()
PyPlot.hist(θs[:, 1], bins=div(size(θs[:, 1], 1),10), density = true, histtype = "step", label = method)
θθ = reshape(Array(LinRange(-2, 2, 1000)), (1000, 1))
ρρ = zeros(size(θθ, 1))
for i =1:size(θθ, 1)
    ρρ[i]= ln_posterior(θθ[i,:])
end
ρρ.= exp.(ρρ)
ρρ /= (sum(ρρ)*(θθ[2,1] - θθ[1,1]))
PyPlot.plot(θθ, ρρ, "--", label="Reference")
PyPlot.legend()


text:
using BenchmarkTools

mutable struct Node
    id::Int
    kids::Vector{Node}
    data::Vector{Float64}
end

mutable struct FlatNode
    id::Int
    data::Vector{Float64}
end

function buildtree!(root, level, mxdepth, n_kids = 4)
    if level < mxdepth

        for i = 1:n_kids
            push!(root.kids, Node(i, Node[], rand(1000)))
        end
        for j = 1:n_kids
            buildtree!(root.kids[j], level+1, mxdepth, n_kids)
        end
    end
end

function count_nodes(mxdepth, n_kids = 4)
    n_nodes = 0
    for l = 1:mxdepth
        n_nodes += n_kids^(l - 1)
    end
    n_nodes
end

function flatnode_array(n_nodes)
    res = Vector{FlatNode}(0)
    for i = 1:n_nodes
        push!(res, FlatNode(i, rand(1000)))
    end
    return res
end

LEVEL = 6
N_KIDS = 6

xroot = Node(1, Node[], rand(1000))
n_nodes = count_nodes(LEVEL, N_KIDS)

@time buildtree!(xroot, 0, LEVEL, N_KIDS)
@time a = flatnode_array(n_nodes)

# With 4 kids, here are the MiB allcoation ratios
2.667 / 10.734
10.654 / 42.969
42.623 / 171.906
170.498 / 687.656

# With 6 kids
2.024 / 12.44
12.133 / 73.508
72.865 / 441.093


text:
#This file has the definitions for the gray Flux-Limited Diffusion Solver
#The code uses units of
# cm for length
# sh for time (1 sh = 10^-8 s)
# GJ (gigajoules) for energy; 1 GJ is also known as a jerk
# keV for temperature (1 keV is 1.1605 * 10^7 K)

using PyCall
using JLD, HDF5

#radiation constant in GJ/keV^4/cm^3
const a = 0.01372
#speed of light in cm/sh; 1 sh is 10^-8 s
const c = 299.98

"This function maps an Nx*Ny 1-D array to a i,j on an Nx by Ny grid"
function vect_to_xy(l,Nx,Ny)
    j = div(l,Nx) + 1
    i = mod(l-j,Nx) + 1
    (i,j)
end

"This function maps i,j on an Nx by Ny grid to an Nx*Ny 1-D array"
function xy_to_vect(i,j,Nx,Ny)
l = Nx*(j-1) + i
l
end

"""
    create_A(D,sigma,Nr,Nz,Lr,Lz; lower_z = "refl",upper_z="refl",upper_r="refl")

Compute the matrix and RHS for an RZ diffusion problem
# Arguments
- `D::Float64(Nr,Nz)`: matrix of diffusion coefficients
- `sigma::Float64(Nr,Nz)`: matrix of absorption coefficients
- `Nr::Integer`: number of cells in r
- `Nz::Integer`: number of cells in z
- `Lr::Float64`: size of domain in r
- `Lz::Float64`: size of domain in z
...
"""
function create_A(Dleft,Dright, Dtop,Dbottom, sigma,Nr,Nz,Lr,Lz; lower_z = "refl",upper_z="refl",upper_r="refl")

    redges = linspace(0,Lr,Nr+1)
    zedges = linspace(0,Lz,Nz+1)

    A = spzeros(Nr*Nz,Nr*Nz)
    b = zeros(Nr*Nz)
    dz = Lz/Nz
    dr = Lr/Nr
    idz = 1/dz
    idz2 = idz*idz
    idr = 1/dr

    #loop over cells
    for i in 1:Nr
        for j in 1:Nz
            here = xy_to_vect(i,j,Nr,Nz)
            Vij = pi*(redges[i+1]^2 - redges[i]^2)*dz
            iVij = 1/Vij
            Siplus = 2*pi*redges[i+1]
            Siminus = 2*pi*redges[i]
            #add in diag
            A[here,here] = sigma[i,j]
            #add in i+1
            if (i< Nr)
                Diplus = 2*Dright[i,j]*Dleft[i+1,j]/(Dright[i,j] + Dleft[i+1,j])
                tmp_zone = xy_to_vect(i+1,j,Nr,Nz)
                A[here,tmp_zone] = -Siplus*iVij*idr*Diplus
                A[here,here] += Siplus*iVij*idr*Diplus
            elseif (upper_r == "vacuum")
                A[here,here] += Siplus*iVij*0.5*c #*idr*D[i,j]
            elseif (isa(upper_r,Float64))
                A[here,here] += Siplus*iVij*0.5*c #*idr*D[i,j]
                b[here] += Siplus*iVij*0.5*upper_r*c
            end
            #add in i-1
            if (i>1)
                Diminus = 2*Dleft[i,j]*Dright[i-1,j]/(Dleft[i,j] + Dright[i-1,j])
                tmp_zone = xy_to_vect(i-1,j,Nr,Nz)
                A[here,tmp_zone] = -Siminus*iVij*idr*Diminus
                A[here,here] += Siminus*iVij*idr*Diminus
            end
            #add in j+1
            if (j < Nz)
                Djplus = 2*Dtop[i,j]*Dbottom[i,j+1]/(Dtop[i,j] + Dbottom[i,j+1])
                tmp_zone = xy_to_vect(i,j+1,Nr,Nz)
                #println("l = $(tmp_zone) i = $(i), j = $(j+1)")
                A[here,tmp_zone] = -idz2*Djplus
                A[here,here] += idz2*Djplus
            elseif (upper_z == "vacuum")
                A[here,here] += idz*0.5*c   #idz2*D[i,j]
            elseif (isa(upper_z,Float64))
                A[here,here] += idz*0.5*c #idz2*D[i,j]
                b[here] += idz*upper_z*0.5*c #idz2*D[i,j]*upper_z
            end
            #add in j-1
            if (j>1)
                Djminus = 2*Dbottom[i,j]*Dtop[i,j-1]/(Dbottom[i,j] + Dtop[i,j-1])
                tmp_zone = xy_to_vect(i,j-1,Nr,Nz)
                A[here,tmp_zone] = -idz2*Djminus
                A[here,here] += idz2*Djminus
            elseif (lower_z == "vacuum")
                A[here,here] += idz*0.5*c #idz2*D[i,j]
            elseif (isa(lower_z,Float64))
                A[here,here] += idz*c #idz2*D[i,j]
                b[here] += idz*lower_z*c #idz2*D[i,j]*lower_z
            end
        end #for j
    end #for i
    A,b
end #function

function time_dep_RT(Tfinal,delta_t,scale_change,D_func,sigma_func,scale,r_apt,sc_length,Q_func,Cv_func,EOS,invEOS,
Nr,Nz,Lr,Lz,theta;lower_z = "refl",upper_z="refl",upper_r="refl", LOUD=0, fname="tmp.jld")
    #sigma,Q,Cv,EOS,invEOS are functions of t,T,Nr,Nz,Lr,Lz
    #D is function of t,T,Nr,Nz,,Lr,Lz,Er,sigma
    done = false
    Lz,Lr,Nr,Nz,T = scale_change(scale,Lr,Lz,Nr,Nz)
    Er = a*T.^4
    time = 0
    dz = Lz/Nz
    dr = Lr/Nr
    times = [0.]
    steps = Int64(ceil(Tfinal/delta_t+1))
    println(steps)
    Er_t = zeros(Nz*Nr,steps+1)
    Er_t[:,1] = Er
    T_t = zeros(Nz*Nr,steps+1)
    T_t[:,1] = T
    step = 1
    while !(done)
        dt = min(Tfinal-time, delta_t)
        t = time+dt
        sigma = sigma_func(t,T,Nr,Nz,Lr,Lz,scale,r_apt,sc_length,theta)
        Cv = Cv_func(t,T,Nr,Nz,Lr,Lz)
        Q = Q_func(t,T,Nr,Nz,Lr,Lz,scale)
        beta = 4*a*c*reshape(T,Nr,Nz).^3./Cv
        f = 1./(1+beta.*sigma*dt)
        sigma_a = sigma.*f
        sigma_star = c*sigma_a + 1.0/(dt)
        Trect = reshape(T,Nr,Nz)

        #define Tleft
        tmp = ones(Nr,Nz)
        tmp[2:Nr,:] = 0.5*(Trect[1:(Nr-1),:] + Trect[2:Nr,:])
        sigma_left = sigma_func(t,tmp,Nr,Nz,Lr,Lz,scale,r_apt,sc_length,theta)
        Dleft = D_func(t,tmp,Nr,Nz,Lr,Lz,Er,sigma_left)
        #define Tright
        tmp = ones(Nr,Nz)
        tmp[1:(Nr-1),:] = 0.5*(Trect[1:(Nr-1),:] + Trect[2:Nr,:])
        sigma_right = sigma_func(t,tmp,Nr,Nz,Lr,Lz,scale,r_apt,sc_length,theta)
        Dright = D_func(t,tmp,Nr,Nz,Lr,Lz,Er,sigma_right)
        #define Ttop
        tmp = ones(Nr,Nz)
        tmp[:,1:(Nz-1)] = 0.5*(Trect[:,1:(Nz-1)] + Trect[:,2:Nz])
        sigma_top = sigma_func(t,tmp,Nr,Nz,Lr,Lz,scale,r_apt,sc_length,theta)
        Dtop = D_func(t,tmp,Nr,Nz,Lr,Lz,Er,sigma_top)
        #define Tbottom
        tmp = ones(Nr,Nz)
        tmp[:,2:Nz] = 0.5*(Trect[:,1:(Nz-1)] + Trect[:,2:Nz])
        sigma_bottom = sigma_func(t,tmp,Nr,Nz,Lr,Lz,scale,r_apt,sc_length,theta)
        Dbottom = D_func(t,tmp,Nr,Nz,Lr,Lz,Er,sigma_bottom)

        A,b = create_A(c*Dleft,c*Dright, c*Dtop,c*Dbottom, sigma_star, Nr, Nz, Lr, Lz, lower_z=lower_z,upper_z=upper_z, upper_r=upper_r);
        b += Er/(dt) + reshape(Q,(Nr*Nz)) + c*a*T.^4.*reshape(sigma_a,Nr*Nz)
        Er = A\b
        dE = EOS(t,T,Nr,Nz,Lr,Lz) + c*dt*reshape(sigma_a,Nr*Nz).*(Er - a*T.^4)
        T = copy(invEOS(t,dE,Nr,Nz,Lr,Lz))

        time += dt
        done = time >= Tfinal
        times = push!(times,time)
        if (LOUD >0)
            println("Step $(step), t = $(time)")
        end
        step += 1
        Er_t[:,step] = copy(Er)
        T_t[:,step] = copy(T)

        save(fname, "Nr", Nr, "Nz", Nz, "dr", dr, "dz", dz, "times", times, "T", T_t, "Er", Er_t)
    end #timestep loop
    if (LOUD == -1)
        println("Step $(step-1), t = $(time)")
    end
    times,Er_t,T_t,Nr,Nz,Lr,Lz
end #function

function meshgrid(x,y)
    Nx = length(x)
    Ny = length(y)
    X = zeros(Nx,Ny)
    Y = zeros(Nx,Ny)
    for i in 1:Nx
        for j in 1:Ny
            X[i,j] = x[i]
            Y[i,j] = y[j]
        end
    end
    X,Y
end

max_stars_repo_name: DrRyanMc/HEDPRadDiffusion

max_issues_repo_name: DrRyanMc/HEDPRadDiffusion

max_forks_repo_name: DrRyanMc/HEDPRadDiffusion


text:
# check that all interace functions are provided

# --- LazySet ---

# support vector
@test check_method_implementation(LazySet, σ,
                                  Function[S -> (Vector{Float64}, S{Float64})])
# dimension
check_method_implementation(LazySet, dim, Function[S -> (S{Float64},)])

# copy
check_method_implementation(LazySet, copy, Function[S -> (S,)])

# --- AbstractPolytope ---

# vertices list
global test_suite_polyhedra
if test_suite_polyhedra
    exclusions = Type[]
else
    exclusions = Type[HPolytope, VPolytope]
end
@test check_method_implementation(AbstractPolytope, vertices_list,
                                  Function[S -> (S{Float64},)],
                                  ignore_types=exclusions)
@test check_method_implementation(AbstractPolytope, constraints_list,
                                  Function[S -> (S{Float64},)])

# --- AbstractCentrallySymmetric ---

# center
@test check_method_implementation(AbstractCentrallySymmetric, center,
                                  Function[S -> (S{Float64},)])
@test check_method_implementation(AbstractCentrallySymmetricPolytope, center,
                                  Function[S -> (S{Float64},)])

# --- AbstractHyperrectangle ---

# radius_hyperrectangle (2x)
@test check_method_implementation(AbstractHyperrectangle, radius_hyperrectangle,
                                  Function[S -> (S{Float64},)])
@test check_method_implementation(AbstractHyperrectangle, radius_hyperrectangle,
                                  Function[S -> (S{Float64}, Int)])

# --- AbstractSingleton ---

# element
@test check_method_implementation(AbstractSingleton, element,
                                  Function[S -> (S{Float64},)])
@test check_method_implementation(AbstractSingleton, element,
                                  Function[S -> (S{Float64}, Int)])

# --- AbstractPolygon ---

# tohrep
@test check_method_implementation(AbstractPolygon, tohrep,
                                  Function[S -> (S{Float64},)])
# tovrep
@test check_method_implementation(AbstractPolygon, tovrep,
                                  Function[S -> (S{Float64},)])

# --- AbstractHPolygon ---

# no functions yet


text:
# Indexing and dimensions

export workitemIdx, workgroupIdx
export threadIdx, blockIdx

@generated function _index(::Val{fname}, ::Val{name}, ::Val{range}) where {fname, name, range}
    T_int32 = LLVM.Int32Type(JuliaContext())

    # create function
    llvm_f, _ = create_function(T_int32)
    mod = LLVM.parent(llvm_f)

    # generate IR
    Builder(JuliaContext()) do builder
        entry = BasicBlock(llvm_f, "entry", JuliaContext())
        position!(builder, entry)

        # call the indexing intrinsic
        intr_typ = LLVM.FunctionType(T_int32)
        intr = LLVM.Function(mod, "llvm.amdgcn.$fname.id.$name", intr_typ)
        idx = call!(builder, intr)

        # attach range metadata
        range_metadata = MDNode([ConstantInt(Int32(range.start), JuliaContext()),
                                 ConstantInt(Int32(range.stop), JuliaContext())],
                                JuliaContext())
        metadata(idx)[LLVM.MD_range] = range_metadata

        ret!(builder, idx)
    end

    call_function(llvm_f, UInt32)
end

# TODO: look these up for the current agent/queue
const max_block_size = (x=1024, y=1024, z=1024)
const max_grid_size  = (x=2147483647, y=65535, z=65535)

for dim in (:x, :y, :z)
    # Workitem index
    fname = Symbol("workitem")
    fn = Symbol("workitemIdx_$dim")
    intr = Symbol("$dim")
    @eval @inline $fn() = Int(_index($(Val(fname)), $(Val(intr)), $(Val(0:max_block_size[dim])))) + 1

    # Workgroup index
    fname = Symbol("workgroup")
    fn = Symbol("workgroupIdx_$dim")
    intr = Symbol("$dim")
    @eval @inline $fn() = Int(_index($(Val(fname)), $(Val(intr)), $(Val(0:max_block_size[dim])))) + 1
end

@inline workitemIdx() = (x=workitemIdx_x(), y=workitemIdx_y(), z=workitemIdx_z())
@inline workgroupIdx() = (x=workgroupIdx_x(), y=workgroupIdx_y(), z=workgroupIdx_z())

# For compat with CUDAnative et. al

@inline threadIdx() = (x=workitemIdx_x(), y=workitemIdx_y(), z=workitemIdx_z())
@inline blockIdx() = (x=workgroupIdx_x(), y=workgroupIdx_y(), z=workgroupIdx_z())


text:
# PeriodicLoss

using LaTeXStrings, LossFunctions, Measures, Latexify

Plots.plot(PeriodicLoss(1.5), leg=:outertopright, c=:olive, linewidth=2, legendtitle= L"L\left( r \right) = 1 - cos(\frac {2r\pi}{c})", size=(850,499))
Plots.plot!(PeriodicLoss(1), c=:orangered4, linewidth=2.5)
Plots.plot!(PeriodicLoss(.5), c=:springgreen4, linewidth=2.5)

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct MetricSpecification <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    displayName::Any # spec type: Union{ Nothing, String } # spec name: displayName
    displayDescription::Any # spec type: Union{ Nothing, String } # spec name: displayDescription
    unit::Any # spec type: Union{ Nothing, String } # spec name: unit
    dimensions::Any # spec type: Union{ Nothing, Vector{Dimension} } # spec name: dimensions
    aggregationType::Any # spec type: Union{ Nothing, String } # spec name: aggregationType
    fillGapWithZero::Any # spec type: Union{ Nothing, Bool } # spec name: fillGapWithZero
    category::Any # spec type: Union{ Nothing, String } # spec name: category
    resourceIdDimensionNameOverride::Any # spec type: Union{ Nothing, String } # spec name: resourceIdDimensionNameOverride

    function MetricSpecification(;name=nothing, displayName=nothing, displayDescription=nothing, unit=nothing, dimensions=nothing, aggregationType=nothing, fillGapWithZero=nothing, category=nothing, resourceIdDimensionNameOverride=nothing)
        o = new()
        validate_property(MetricSpecification, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(MetricSpecification, Symbol("displayName"), displayName)
        setfield!(o, Symbol("displayName"), displayName)
        validate_property(MetricSpecification, Symbol("displayDescription"), displayDescription)
        setfield!(o, Symbol("displayDescription"), displayDescription)
        validate_property(MetricSpecification, Symbol("unit"), unit)
        setfield!(o, Symbol("unit"), unit)
        validate_property(MetricSpecification, Symbol("dimensions"), dimensions)
        setfield!(o, Symbol("dimensions"), dimensions)
        validate_property(MetricSpecification, Symbol("aggregationType"), aggregationType)
        setfield!(o, Symbol("aggregationType"), aggregationType)
        validate_property(MetricSpecification, Symbol("fillGapWithZero"), fillGapWithZero)
        setfield!(o, Symbol("fillGapWithZero"), fillGapWithZero)
        validate_property(MetricSpecification, Symbol("category"), category)
        setfield!(o, Symbol("category"), category)
        validate_property(MetricSpecification, Symbol("resourceIdDimensionNameOverride"), resourceIdDimensionNameOverride)
        setfield!(o, Symbol("resourceIdDimensionNameOverride"), resourceIdDimensionNameOverride)
        o
    end
end # type MetricSpecification

const _property_map_MetricSpecification = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("displayName")=>Symbol("displayName"), Symbol("displayDescription")=>Symbol("displayDescription"), Symbol("unit")=>Symbol("unit"), Symbol("dimensions")=>Symbol("dimensions"), Symbol("aggregationType")=>Symbol("aggregationType"), Symbol("fillGapWithZero")=>Symbol("fillGapWithZero"), Symbol("category")=>Symbol("category"), Symbol("resourceIdDimensionNameOverride")=>Symbol("resourceIdDimensionNameOverride"))
const _property_types_MetricSpecification = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("displayName")=>"String", Symbol("displayDescription")=>"String", Symbol("unit")=>"String", Symbol("dimensions")=>"Vector{Dimension}", Symbol("aggregationType")=>"String", Symbol("fillGapWithZero")=>"Bool", Symbol("category")=>"String", Symbol("resourceIdDimensionNameOverride")=>"String")
Base.propertynames(::Type{ MetricSpecification }) = collect(keys(_property_map_MetricSpecification))
Swagger.property_type(::Type{ MetricSpecification }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_MetricSpecification[name]))}
Swagger.field_name(::Type{ MetricSpecification }, property_name::Symbol) =  _property_map_MetricSpecification[property_name]

function check_required(o::MetricSpecification)
    true
end

function validate_property(::Type{ MetricSpecification }, name::Symbol, val)
end


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

include("../../src/Sparso.jl")
include("../../src/simple-show.jl")
using Sparso

set_options(SA_ENABLE, SA_VERBOSE, SA_USE_SPMP, SA_CONTEXT)

include("./ipm-ref.jl")
include("utils.jl")

# This code is what we expect to generate by context-sensitive optimizations.
# It is used for debugging purpose only
function ipm_with_context_opt(A, b, p) # A: constraint coefficients, b: constraint rhs, p: objective
  t0 = time()
  (m,n) = size(A)

# set initial point, based on largest element in (A,b,p)
  bigM = maximum(A)
  bigM = maximum([norm(b, Inf) norm(p, Inf) bigM])
  x = 100*bigM*ones(n)
  s = copy(x)
  y = zeros(m)

  bc = 1 + maximum([norm(b) norm(p)])

  relResidual = NaN
  iter=1

  blas1_time = 0.
  spgemm_time = 0.
  fact_time = 0.
  trslv_time = 0.
  spmv_time = 0.

  D = speye_int32(n)

  # There are two kinds of matrices we would like to trace: one is constant value,
  # the other constant structure.
  # For a constant-valued matrix, since it is not defined but referenced in the
  # loop, it must already have a reprsentation in a fixed memory. Thus we can
  # safely pass its name into new_matrix_knob(), which will remember its
  # representation.
  # For a constant-structured matrix, however, it is defined in the loop, and in
  # each definition, it might be assigned a representation in a different memory,
  # even though every time the representation assigned has the same structure.
  # Therefore, we cannot pass its name into new_matrix_knob() at this moment, and
  # new_matrix_knob() should not remember its representation, which may not
  # even exist at this moment. Instead, every time is it defined in the loop, its
  # representation should be recorded into its matrix knob.

  # Hoist AT out of loop. This is currently done in pattern replacement.
  # TODO: have a separate loop invariant hoisting phase.
  __AT__ = (Main.ctranspose)(A::Base.SparseMatrix.SparseMatrixCSC{Float64,Int32})

  # Create mknobs for constant-valued matrices. They are of course constant-structured.
  __mknob__AT___8599 = Sparso.new_matrix_knob(:__AT__,true,true,false,false,false,false)
  __mknobA_8601      = Sparso.new_matrix_knob(:A,     true,true,false,false,false,false)

  # Create mknobs for constant-structured matrices in the source code
  # Note R and B are also single-defs: key to enable context info propagation.
  # D is also single-def, but that is not important for this application.
  __mknobR_8596 = Sparso.new_matrix_knob(:R, false,true,false,false,false,true)
  __mknobB_8598 = Sparso.new_matrix_knob(:B, false,true,false,false,false,true)
  __mknobD_8600 = Sparso.new_matrix_knob(:D, false,true,false,false,false,true)

  # Create mknobs for a function call at each call site, representing the
  # call's output, if the output is constant in structure or value. We do not
  # neeed this for cholmod_factor_inverse_divide, since its output is a vector,
  # while we track only matrices: it is OK, because our principle is: if there
  # is context info, make sure it is always correct and we can use it; if there
  # is not, we simply use nothing. Any way, a function call's result must remain
  # the same, no matter it has a fknob with it or not.
  # Any function call's result is a single-def (Defined only the function call).
  __mknob_ADB            = Sparso.new_matrix_knob(:ADB, false,true,false,false,false,true)
  __mknob_cholfact_int32 = Sparso.new_matrix_knob(:cholfact_int32, false,true,false,false,false,true)

  # Create fknobs for a function call at each call site, and add related mknobs
  # For ADB
  __fknob_8602 = Sparso.new_function_knob()
  (Sparso.add_mknob_to_fknob)(__mknob_ADB,__fknob_8602)
  (Sparso.add_mknob_to_fknob)(__mknob__AT___8599,__fknob_8602)
  (Sparso.add_mknob_to_fknob)(__mknobD_8600,__fknob_8602)
  (Sparso.add_mknob_to_fknob)(__mknobA_8601,__fknob_8602)

  (Sparso.set_derivative)(__mknobA_8601, Sparso.DERIVATIVE_TYPE_TRANSPOSE, __mknob__AT___8599)

  # For cholfact_int32
  __fknob_8622 = (Sparso.new_function_knob)()
  (Sparso.add_mknob_to_fknob)(__mknobB_8598,__fknob_8622)

  # For cholmod_factor_inverse_divide
  __fknob_8623 = (Sparso.new_function_knob)()
  (Sparso.add_mknob_to_fknob)(__mknobR_8596,__fknob_8623)

  fknob_spmv1 = (Sparso.new_function_knob)()
  (Sparso.add_mknob_to_fknob)(__mknob__AT___8599, fknob_spmv1)

  fknob_spmv2 = (Sparso.new_function_knob)()
  (Sparso.add_mknob_to_fknob)(__mknobA_8601, fknob_spmv2)

  fknob_spmv3 = (Sparso.new_function_knob)()
  (Sparso.add_mknob_to_fknob)(__mknobA_8601, fknob_spmv3)

  fknob_spmv4 = (Sparso.new_function_knob)()
  (Sparso.add_mknob_to_fknob)(__mknob__AT___8599, fknob_spmv4)

  for iter=1:200

    # compute residuals
    spmv_time -= time()
    #Rd = __AT__*y + s - p
    Rd = Sparso.SpMV(1, __AT__, y, 1, s - p, fknob_spmv1)
    #Rp = A*x - b
    Rp = Sparso.SpMV(1, A, x, -1, b, fknob_spmv2)
    spmv_time += time()

    blas1_time -= time()
    #Rc = x.*s
    Rc = Sparso.element_wise_multiply(x, s)
    #mu = mean(Rc)
    mu = Sparso.sum(Rc)/length(Rc)
    relResidual = norm([Rd; Rp; Rc])/bc
    blas1_time += time()

    if (relResidual <= 1e-7 && mu <= 1e-7) break; end

    blas1_time -= time()
    #Rc = Rc - min(0.1, 100*mu)*mu
    Sparso.WAXPB!(Rc, 1, Rc, -min(0.1, 100*mu)*mu)

    # set up the scaling matrix, and form the coefficient matrix for
    # the linear system
    #d = min(5.e+15, x./s)
    d = Sparso.element_wise_divide(x, s)
    Sparso.min!(d, d, 5.e+15)
    blas1_time += time()

    spgemm_time -= time()
    D.nzval = d
    #B = A*D*A'
    B = Sparso.ADB(A,D,__AT__,__fknob_8602)
    Sparso.propagate_matrix_info(__mknobB_8598, __mknob_ADB)
    spgemm_time += time()

    # use the form of the Cholesky routine "cholinc" that's best
    # suited to interior-point methods
    fact_time -= time()
    #R = cholfact_int32(B)
    R = Sparso.cholfact_int32(B,__fknob_8622)
    fact_time += time()

    # set up the right-hand side
    blas1_time -= time()
    #t1 = x.*Rd - Rc;
    t1 = Sparso.element_wise_multiply(x, Rd)
    Sparso.WAXPBY!(t1, 1, t1, -1, Rc)
    blas1_time += time()

    spmv_time -= time()
    #t2 = -(Rp + A*(t1./s));
    t2 = Sparso.SpMV(-1, A, t1./s, -1, Rp, fknob_spmv3)
    spmv_time += time()

    # solve it and recover the other step components
    trslv_time -= time()
    #dy = R\t2
    dy = Sparso.cholfact_inverse_divide(R,t2,__fknob_8623)
    trslv_time += time()

    spmv_time -= time()
    #temp = A'*dy
    temp = Sparso.SpMV(__AT__, dy, fknob_spmv4)
    spmv_time += time()

    blas1_time -= time()
    #dx = (x.*temp + t1)./s
    Sparso.element_wise_multiply!(temp, x, temp)
    Sparso.WAXPBY!(temp, 1, temp, 1, t1)
    dx = Sparso.element_wise_divide(temp, s)
    #ds = -(s.*dx + Rc)./x
    Sparso.element_wise_multiply!(temp, s, dx)
    Sparso.WAXPBY!(temp, -1, temp, -1, Rc)
    ds = Sparso.element_wise_divide(temp, x)

    tau = max(.9995, 1 - mu)
    #ap = -1/minimum([dx./x; -1])
    Sparso.element_wise_divide!(temp, dx, x)
    ap = -1/min(Sparso.minimum(temp), -1)
    #ad = -1/minimum([ds./s; -1])
    Sparso.element_wise_divide!(temp, ds, s)
    ad = -1/min(Sparso.minimum(temp), -1)
    ap = tau*ap
    ad = tau*ad
    #x = x + ap*dx
    #s = s + ad*ds
    #y = y + ad*dy
    Sparso.WAXPBY!(x, 1, x, ap, dx)
    Sparso.WAXPBY!(s, 1, s, ad, ds)
    Sparso.WAXPBY!(y, 1, y, ad, dy)
    blas1_time += time()
  end

  (Sparso.delete_function_knob)(__fknob_8602)
  (Sparso.delete_function_knob)(__fknob_8622)
  (Sparso.delete_function_knob)(__fknob_8623)
  (Sparso.delete_matrix_knob)(__mknobR_8596)
  (Sparso.delete_matrix_knob)(__mknobB_8598)
  (Sparso.delete_matrix_knob)(__mknob__AT___8599)
  (Sparso.delete_matrix_knob)(__mknobA_8601)
  (Sparso.delete_matrix_knob)(__mknobD_8600) # line 56:
  (Sparso.delete_matrix_knob)(__mknob_ADB)
  (Sparso.delete_matrix_knob)(__mknob_cholfact_int32) # line 56:

  blas1_time -= time()
  f = p'*x
  blas1_time += time()

  x, time() - t0, spgemm_time, fact_time, blas1_time, trslv_time, spmv_time, iter, relResidual, f[1][1]
end

if length(ARGS) == 0
  A, b, p = load_ipm_input("ipm/mps/osa-14")
else
  A, b, p = load_ipm_input(ARGS[1])
end

m = size(A, 1)
n = size(A, 2)
println("Problem size = [$m $n]")
#println("\tsum of A=", sum(A))
#println("\tsum of b=", sum(b))
#println("\tsum of p=", sum(p))

println("Original: ")
x, ref_total_time, spgemm_time, fact_time, blas1_time, trslv_time, spmv_time,
    iter, relResidual, objval = ipm_ref(A, b, p)
x, ref_total_time, spgemm_time, fact_time, blas1_time, trslv_time, spmv_time,
    iter, relResidual, objval = ipm_ref(A, b, p)
@printf "Original iter %2i, resid = %9.2e, objval = %e\n" iter relResidual objval
@printf "\nref_total_time = %f\n" ref_total_time
@printf "spgemm = %f fact = %f blas1 = %f trslv = %f spmv = %f\n" spgemm_time fact_time blas1_time trslv_time spmv_time

println("\n\nWith manual context-sensitive optimization: ")
x, ref_total_time, spgemm_time, fact_time, blas1_time, trslv_time, spmv_time,
    iter, relResidual, objval = ipm_with_context_opt(A, b, p)
Sparso.set_knob_log_level(1)
x, ref_total_time, spgemm_time, fact_time, blas1_time, trslv_time, spmv_time,
    iter, relResidual, objval = ipm_with_context_opt(A, b, p)
@printf "\nopt_total_time = %f\n" ref_total_time
@printf "spgemm = %f fact = %f blas1 = %f trslv = %f spmv = %f\n" spgemm_time fact_time blas1_time trslv_time spmv_time
@printf "Manual_context iter %2i, resid = %9.2e, objval = %e\n" iter relResidual objval

#println("\n\nWith manual context-sensitive optimization: ")
#ipm_ref_simplified_with_context_opt(A, b, p)
#println("\tsum of x=", sum(x))

max_stars_repo_name: IntelLabs/Sparso

max_issues_repo_name: IntelLabs/Sparso

max_forks_repo_name: IntelLabs/Sparso


text:
"""
    findthomsen(p::Vector{VMeasure},
                     sv::Vector{VMeasure},
                     sh::Vector{VMeasure},
                     vp₀::Tuple{Real,Real},
                     vs₀::Tuple{Real,Real},
                     ϵ₀::Tuple{Real,Real},
                     δ₀::Tuple{Real,Real},
                     γ₀::Tuple{Real,Real})

    Estimate Thomsen parameters, wavespeeds and phase angles from measurements of group P, SV ans SH wave velocities, using quasinewton method. A priori values ()₀ are tuples with value and standard deviation (a priori gaussian).

Return (vp0, vs0, ϵ, δ, γ), θp, θsv, θsh, CMpost

"""
function findthomsen(p::Vector{VMeasure},
                     sv::Vector{VMeasure},
                     sh::Vector{VMeasure},
                     vp₀::Tuple{Real,Real},
                     vs₀::Tuple{Real,Real},
                     ϵ₀::Tuple{Real,Real},
                     δ₀::Tuple{Real,Real},
                     γ₀::Tuple{Real,Real})

    dobs, CDi, np, nsv, nsh = inputdobs(p,sv,sh)

    #make vector of model parameters (first guess)
    m₀ = [vp₀[1],vs₀[1],ϵ₀[1],δ₀[1],γ₀[1]]
    σm = [vp₀[2],vs₀[2],ϵ₀[2],δ₀[2],γ₀[2]]
    for x in p
        push!(m₀, phaseangleP(x.angle, vp₀[1]/vs₀[1], ϵ₀[1], δ₀[1]))
        push!(σm, 100.)
    end
    for x in sv
        push!(m₀, phaseangleSV(x.angle, vp₀[1]/vs₀[1], ϵ₀[1], δ₀[1]))
        push!(σm, 100.)
    end
    for x in sh
        push!(m₀, phaseangleSH(x.angle, γ₀[1]))
        push!(σm, 100.)
    end

    CM = Diagonal(σm.^2)
    CMi = inv(CM)

    m, CMpost = quasinewton(dobs, _gt, m₀, CMi, CDi, (np,nsv,nsh))

    return m[1:5], m[6:5+np], m[6+np:5+np+nsv], m[6+np+nsv:5+np+nsv+nsh], CMpost

end

"""
    findthomsen(p::Vector{VMeasure},
                     sv::Vector{VMeasure},
                     sh::Vector{VMeasure},
                     vp₀::Real,
                     vs₀::Real,
                     ϵ₀::Real,
                     δ₀::Real,
                     γ₀::Real)

    Method using just numbers as input for a priori values, assuming std=10.0 for all parameters.

"""
function findthomsen(p::Vector{VMeasure},
                     sv::Vector{VMeasure},
                     sh::Vector{VMeasure},
                     vp₀::Real,
                     vs₀::Real,
                     ϵ₀::Real,
                     δ₀::Real,
                     γ₀::Real)

    return findthomsen(p,sv,sh,(vp₀,10.),(vs₀,10.),(ϵ₀,10.),(δ₀,10.),(γ₀,10.))

end

"""
    findmoduli(p::Vector{VMeasure},
                    sv::Vector{VMeasure},
                    sh::Vector{VMeasure},
                    ρ::AbastractFloat,
                    c11₀::Tuple{Real,Real},
                    c13₀::Tuple{Real,Real},
                    c33₀::Tuple{Real,Real},
                    c55₀::Tuple{Real,Real},
                    c66₀::Tuple{Real,Real}

    Estimate stiffness matrix and phase angles from measurements of group P, SV ans SH wave velocities, using quasinewton method. A priori values ()₀ are tuples with value and standard deviation (a priori gaussian).

Return Cijkl, θp, θsv, θsh, CMpost

"""
function findmoduli(p::Vector{VMeasure},
                    sv::Vector{VMeasure},
                    sh::Vector{VMeasure},
                    ρ::Real,
                    c11₀::Tuple{Real,Real},
                    c13₀::Tuple{Real,Real},
                    c33₀::Tuple{Real,Real},
                    c55₀::Tuple{Real,Real},
                    c66₀::Tuple{Real,Real})

    dobs, CDi, np, nsv, nsh = inputdobs(p,sv,sh)

    #make vector of model parameters (first guess)
    m₀ = [c11₀[1],c13₀[1],c33₀[1],c55₀[1],c66₀[1]]
    σm = [c11₀[2],c13₀[2],c33₀[2],c55₀[2],c66₀[2]]

    C₀ = stiffmatrixTI(m₀...)
    ϵ₀, δ₀, γ₀ = moduli2thomsen(C₀)
    vp,vs = velocities0(C₀,ρ) # here density does not matter

    for x in p
        push!(m₀, phaseangleP(x.angle, vp/vs, ϵ₀, δ₀))
        push!(σm, 100.)
    end
    for x in sv
        push!(m₀, phaseangleSV(x.angle, vp/vs, ϵ₀, δ₀))
        push!(σm, 100.)
    end
    for x in sh
        push!(m₀, phaseangleSH(x.angle, γ₀))
        push!(σm, 100.)
    end

    CM = Diagonal(σm.^2)
    CMi = inv(CM)

    m, CMpost = quasinewton(dobs, _gc, m₀, CMi, CDi, (ρ,np,nsv,nsh))

    return stiffmatrixTI(m[1:5]...), m[6:5+np], m[6+np:5+np+nsv], m[6+np+nsv:5+np+nsv+nsh], CMpost
end

"""
	findmoduli(p::Vector{VMeasure},
                   sv::Vector{VMeasure},
                   sh::Vector{VMeasure},
                   ρ::Real,
                   vp₀::Real,
                   vs₀::Real,
                   errlogC::Real)

Method with isotropic vp,vs as input.
"""
function findmoduli(p::Vector{VMeasure},
                    sv::Vector{VMeasure},
                    sh::Vector{VMeasure},
                    ρ::Real,
                    vp₀::Real,
                    vs₀::Real,
                    errlogC::Real)

    C0 = stiffmatrixTI(vpvsYoungsPoisson(vp₀,vs₀,ρ))

    return findmoduli(p,sv,sh,ρ,
                      (C0[1,1],errlogC),
                      (C0[1,3],errlogC),
                      (C0[3,3],errlogC),
                      (C0[5,5],errlogC),
                      (C0[6,6],errlogC))

end
"""
    findmoduli(p::Vector{VMeasure},
                    sv::Vector{VMeasure},
                    sh::Vector{VMeasure},
                    ρ::Real,
                    c11₀::Real,
                    c13₀::Real,
                    c33₀::Real,
                    c55₀::Real,
                    c66₀::Real,
                    errlogC::Real)

Method using just a unique error on moduli.
"""
function findmoduli(p::Vector{VMeasure},
                    sv::Vector{VMeasure},
                    sh::Vector{VMeasure},
                    ρ::Real,
                    c11₀::Real,
                    c13₀::Real,
                    c33₀::Real,
                    c55₀::Real,
                    c66₀::Real,
                    errlogC::Real)

    return findmoduli(p,sv,sh,ρ,
                      (c11₀,errlogC),
                      (c13₀,errlogC),
                      (c33₀,errlogC),
                      (c55₀,errlogC),
                      (c66₀,errlogC))

end

"""
	inputdobs(p,sv,sh)
"""
function inputdobs(p,sv,sh)
    np = length(p)
    nsv = length(sv)
    nsh = length(sh)

    # make vector of data and data covariance
    dobs = Float64[]
    σd = Float64[]
    for x in vcat(p,sv,sh)
        append!(dobs, [x.angle, log(x.value)])
        append!(σd, [x.err_angle, x.err_value])
    end

    CD = Diagonal(σd.^2)
    CDi = inv(CD)

    return dobs, CDi, np, nsv, nsh
end

function _gc(m,ρ,np,nsv,nsh)
    C₀ = stiffmatrixTI(m[1:5]...)
    ϵ, δ, γ = moduli2thomsen(C₀)
    vp0,vs0 = velocities0(C₀,ρ)
    mt = vcat([vp0,vs0,ϵ,δ,γ], m[6:end])
    return _gt(mt, np, nsv, nsh)
end

function _gt(m,np,nsv,nsh)
    vp0 = m[1]
    vs0 = m[2]
    ϵ = m[3]
    δ = m[4]
    γ = m[5]
    θp = m[6:5+np]
    θsv= m[6+np:5+np+nsv]
    θsh= m[6+np+nsv:5+np+nsv+nsh]

    d = zeros(eltype(m), 2np+2nsv+2nsh)  #NEED eltype HERE to work with ForwardDiff!
    c = 1
    for x in θp
        vp = vphaseP(x,vp0,vs0,ϵ,δ)
        dvvp = dvPdthetav(x,vp0/vs0,ϵ,δ)
        d[c] = groupangle(x,dvvp)
        d[c+1] = log(vgroup(vp, dvvp))
        c+=2
    end
    for x in θsv
        vsv = vphaseSV(x,vp0,vs0,ϵ,δ)
        dvvsv = dvSVdthetav(x,vp0/vs0,ϵ,δ)
        d[c] = groupangle(x,dvvsv)
        d[c+1] = log(vgroup(vsv, dvvsv))
        c+=2
    end
    for x in θsh
        d[c] = groupangleSH(x,γ)
        d[c+1] = log(vgroupSH(d[c], vs0, γ))
        c+=2
    end

    return d
end


text:
using ProgressMeter
using .Threads
using Lasso
using NPZ
using Statistics
using Random
using StatsBase
using LinearAlgebra
import Base: isless

struct TestResult
    i::Int64
    f::Float64
end
isless(a::TestResult, b::TestResult) = isless(a.f, b.f)

struct InvResult{T,S}
    i::Int
    λmin::T
    dn::T
    m::S
end

T = Float64
X1 = npzread("Kernels/G1.npy")
X2 = npzread("Kernels/G2.npy")
X3 = npzread("Kernels/G4.npy")
nfactor = npzread("Testing/nfactor.npy")

ϵ = 1e-2
λ = 1.0
lnorm = -log(ϵ)*(2*size(X1,2))

X = [X3 zeros(T, size(X3));
     zeros(T, size(X3)) X3;
     nfactor*X1 nfactor*X2]
Xn = maximum(sqrt.(sum(X.^2, dims=2)))
X ./= Xn

n = size(X3, 1)
ndas = size(X1, 1)

function Ml(μ, X, ϵ)
    M = ϵ*I + X'*Diagonal(μ)*X
    Mc = cholesky(Hermitian(M))
    l = μ.^2 .*diag(X*(Mc\X'))
    return (M, Mc, l)
end

Ml(μ) = Ml(μ, X, ϵ)

function f(μ, λ, lnorm, n)
    M, Mc, l = Ml(μ)
    l2 = l[1:n] .+ l[(n+1):2*n] # coherence of point sensors only; sum across channels
    return -logdet(M) + λ*maximum(l2) - lnorm
end

f(μ) = f(μ, λ, lnorm, n)

function testf(sta, i)
    u_test = vcat(zeros(64+64), ones(size(X1,1)))
    u_test[sta] .= 1.0
    u_test[sta.+64] .= 1.0
    return TestResult(i, f(u_test))
end

function fit_model(G,dv,i)
    dn = std(dv)
    dv ./= dn
    path = fit(LassoPath, G, dv, λminratio = 1e-5, standardize=false, intercept=false, maxncoef=size(G,2))
    λmin = path.λ[argmin(aicc(path))]
    m = coef(path; select=MinAICc())
    return InvResult(i, λmin, dn, m)
end

T = Float64

G1 = npzread("Kernels/G1.npy")
G2 = npzread("Kernels/G2.npy")
G3 = npzread("Kernels/G4.npy")

d1 = npzread("Testing/das_wvt_data.npy")
d2 = npzread("Testing/nodal_wvt_data_x.npy")
d3 = npzread("Testing/nodal_wvt_data_y.npy")
nfactor = npzread("Testing/nfactor.npy")

G = [nfactor*G1 nfactor*G2;
     G3 zeros(T, size(G3));
     zeros(T, size(G3)) G3]

Gn = std(G)

G ./= Gn

d = [d1; d2; d3]

s_size = [4, 8, 16, 32, 64]
n_runs = 10
n_selecting = 1000

p = Progress(n_runs*(length(s_size)-1)+1, dt=0.5, barglyphs=BarGlyphs("[=> ]"), barlen=50, color=:yellow)

for B in s_size
    pvec = npzread("IncoherenceResults/res_inc_$B.npy")

    for i in 1:(B!=64 ? n_runs : 1)
        rng = MersenneTwister(i+50)

        sensors_vec = [sample(rng, 1:64, aweights(pvec), B, replace=false) for k = 1:n_selecting]
        sensors_res = [testf(sta, i) for (i, sta) in enumerate(sensors_vec)]
        best = minimum(sensors_res)
        sensors = sensors_vec[best.i]

        usensors = sample(rng, 1:64, B, replace=false)
        μ = [(i in sensors ? 1.0 : 0.0) for i in 1:64]
        uμ = [(i in usensors ? 1.0 : 0.0) for i in 1:64]
        μl = vcat(ones(size(G1, 1)), μ, μ)
        uμl = vcat(ones(size(G1, 1)), uμ, uμ)
        Gμ = Diagonal(μl)*G
        Guμ = Diagonal(uμl)*G
        results = Array{Any,1}(undef, size(d,2))
        uresults = Array{Any,1}(undef, size(d,2))

        @threads for j = 1:size(d,2)
            results[j] = fit_model(Gμ, d[:,j], j)
            uresults[j] = fit_model(Guμ, d[:,j], j)
        end

        sres = reduce(hcat, [x.m.*x.dn for x in results])
        usres = reduce(hcat, [x.m.*x.dn for x in uresults])

        npzwrite("Testing/LassoOutputInc/32B_$(B)_run_$i.npy", Array{Float64}(sres))
        npzwrite("Testing/LassoOutputInc/32B_$(B)_urun_$i.npy", Array{Float64}(usres))
        npzwrite("Testing/LassoOutputInc/32B_$(B)_run_$(i)_sensors.npy", sensors)
        npzwrite("Testing/LassoOutputInc/32B_$(B)_urun_$(i)_sensors.npy", usensors)
        next!(p)
    end
end

max_stars_repo_name: jbmuir/DAS-Reconstruction

max_issues_repo_name: jbmuir/DAS-Reconstruction

max_forks_repo_name: jbmuir/DAS-Reconstruction


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

using Behavior.Gherkin.Experimental
using Behavior.Gherkin.Experimental: BadExpectedEOFParseResult, BadUnexpectedEOFParseResult, And
using Behavior.Gherkin: Given, When, Then, Scenario, Feature, ScenarioStep, Background
using Behavior.Gherkin: DataTable, ScenarioOutline
using Test

include("combinators_test.jl")
include("gherkin_combinators_test.jl")
include("scenarios_test.jl")


text:
function solve{uType,tType,isinplace,algType<:OrdinaryDiffEqAlgorithm,F,recompile_flag}(
  prob::AbstractODEProblem{uType,tType,isinplace,F},
  alg::algType,timeseries=[],ts=[],ks=[],recompile::Type{Val{recompile_flag}}=Val{true};
  kwargs...)

  integrator = init(prob,alg,timeseries,ts,ks,recompile;kwargs...)
  solve!(integrator)
  integrator.sol
end

function init{uType,tType,isinplace,algType<:OrdinaryDiffEqAlgorithm,F,recompile_flag}(
  prob::AbstractODEProblem{uType,tType,isinplace,F},
  alg::algType,timeseries_init=uType[],ts_init=tType[],ks_init=[],
  recompile::Type{Val{recompile_flag}}=Val{true};
  save_timeseries = true,
  timeseries_steps = 1,
  dense = save_timeseries && !(typeof(alg) <: Discrete),
  saveat = tType[],tstops = tType[],d_discontinuities= tType[],
  calck = (!isempty(setdiff(saveat,tstops)) || dense),
  dt = typeof(alg) <: Discrete && isempty(tstops) ? tType(1) : tType(0),
  adaptive = isadaptive(alg),
  gamma=9//10,
  abstol=1//10^6,
  reltol=1//10^3,
  qmax=qmax_default(alg),qmin=qmin_default(alg),
  qoldinit=1//10^4, fullnormalize=true,
  beta2=beta2_default(alg),
  beta1=beta1_default(alg,beta2),
  maxiters = 1000000,
  dtmax=tType((prob.tspan[end]-prob.tspan[1])),
  dtmin=tType <: AbstractFloat ? tType(10)*eps(tType) : tType(1//10^(10)),
  internalnorm = ODE_DEFAULT_NORM,
  isoutofdomain = ODE_DEFAULT_ISOUTOFDOMAIN,
  unstable_check = ODE_DEFAULT_UNSTABLE_CHECK,
  verbose = true, force_dtmin = false,
  timeseries_errors = true, dense_errors=false,
  advance_to_tstop = false,stop_at_next_tstop=false,
  progress=false,progress_steps=1000,progress_name="ODE",
  progress_message = ODE_DEFAULT_PROG_MESSAGE,
  userdata=nothing,callback=CallbackSet(),
  allow_extrapolation = alg_extrapolates(alg),
  initialize_integrator=true,kwargs...)

  tspan = prob.tspan
  tdir = sign(tspan[end]-tspan[1])

  t = tspan[1]

  if (!(typeof(alg) <: OrdinaryDiffEqAdaptiveAlgorithm) && !(typeof(alg) <: OrdinaryDiffEqCompositeAlgorithm)) && dt == tType(0) && isempty(tstops)
      error("Fixed timestep methods require a choice of dt or choosing the tstops")
  end

  if tspan[1] == tspan[2]
    error("Timespan is trivial")
  end

  d_discontinuities_col = collect(d_discontinuities)

  if tdir>0
    tstops_internal = binary_minheap(convert(Vector{tType},append!(collect(tstops),d_discontinuities_col)))
  else
    tstops_internal = binary_maxheap(convert(Vector{tType},append!(collect(tstops),d_discontinuities_col)))
  end

  if !isempty(tstops) && tstops[end] != tspan[2]
    push!(tstops_internal,tspan[2])
  elseif isempty(tstops)
    push!(tstops_internal,tspan[2])
  end

  if top(tstops_internal) == tspan[1]
    pop!(tstops_internal)
  end
  f = prob.f
  u0 = prob.u0
  uEltype = eltype(u0)

  # Get the control variables

  (uType<:Array || uType <: Number) ? u = copy(u0) : u = deepcopy(u0)

  ks = Vector{uType}(0)

  order = alg_order(alg)

  uEltypeNoUnits = typeof(recursive_one(u))
  tTypeNoUnits   = typeof(recursive_one(t))

  if dt == zero(dt) && adaptive
    dt = tType(ode_determine_initdt(u0,t,tdir,dtmax,uEltype(abstol),uEltypeNoUnits(reltol),internalnorm,prob,order))
  end

  if sign(dt)!=tdir && dt!=tType(0)
    error("dt has the wrong sign. Exiting")
  end

  if typeof(u) <: AbstractArray
    rate_prototype = similar(u/zero(t),indices(u)) # rate doesn't need type info
  else
    rate_prototype = u/zero(t)
  end
  rateType = typeof(rate_prototype) ## Can be different if united

  saveat_vec =  convert(Vector{tType},collect(saveat))
  if !isempty(saveat_vec) && saveat_vec[end] == tspan[2]
    pop!(saveat_vec)
  end

  if tdir>0
    saveat_internal = binary_minheap(saveat_vec)
  else
    saveat_internal = binary_maxheap(saveat_vec)
  end

  if !isempty(saveat_internal) && top(saveat_internal) == tspan[1]
    pop!(saveat_internal)
  end

  d_discontinuities_vec =  convert(Vector{tType},d_discontinuities_col)

  if tdir>0
    d_discontinuities_internal = binary_minheap(d_discontinuities_vec)
  else
    d_discontinuities_internal = binary_maxheap(d_discontinuities_vec)
  end

  callbacks_internal = CallbackSet(callback)

  ### Algorithm-specific defaults ###
  ksEltype = Vector{rateType}

  # Have to convert incase passed in wrong.
  timeseries = convert(Vector{uType},timeseries_init)
  ts = convert(Vector{tType},ts_init)
  ks = convert(Vector{ksEltype},ks_init)
  alg_choice = Int[]

  copyat_or_push!(ts,1,t)
  copyat_or_push!(timeseries,1,u)
  copyat_or_push!(ks,1,[rate_prototype])

  if typeof(alg) <: Discrete
    abstol_internal = zero(uEltype)
    reltol_internal = zero(uEltypeNoUnits)
  else
    abstol_internal = uEltype(uEltype(1)*abstol)
    reltol_internal = uEltypeNoUnits(reltol)
  end

  opts = DEOptions(Int(maxiters),timeseries_steps,save_timeseries,adaptive,abstol_internal,
    reltol_internal,tTypeNoUnits(gamma),tTypeNoUnits(qmax),tTypeNoUnits(qmin),
    tType(dtmax),tType(dtmin),internalnorm,
    tstops_internal,saveat_internal,d_discontinuities_internal,
    userdata,
    progress,progress_steps,
    progress_name,progress_message,
    timeseries_errors,dense_errors,
    tTypeNoUnits(beta1),tTypeNoUnits(beta2),tTypeNoUnits(qoldinit),dense,
    callbacks_internal,isoutofdomain,unstable_check,verbose,calck,force_dtmin,
    advance_to_tstop,stop_at_next_tstop)

  progress ? (prog = Juno.ProgressBar(name=progress_name)) : prog = nothing

  notsaveat_idxs = Int[1]

  k = ksEltype[]

  if uType <: Array
    uprev = copy(u)
  else
    uprev = deepcopy(u)
  end
  if allow_extrapolation
    if uType <: Array
      uprev2 = copy(u)
    else
      uprev2 = deepcopy(u)
    end
  else
    uprev2 = uprev
  end

  cache = alg_cache(alg,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,Val{isinplace})

  if typeof(alg) <: OrdinaryDiffEqCompositeAlgorithm
    id = CompositeInterpolationData(f,timeseries,ts,ks,alg_choice,notsaveat_idxs,cache)
  else
    id = InterpolationData(f,timeseries,ts,ks,notsaveat_idxs,cache)
  end

  if typeof(alg) <: OrdinaryDiffEqCompositeAlgorithm
    sol = build_solution(prob,alg,ts,timeseries,
                      dense=dense,k=ks,interp=id,
                      alg_choice=alg_choice,
                      calculate_error = false)
  else
    sol = build_solution(prob,alg,ts,timeseries,
                      dense=dense,k=ks,interp=id,
                      calculate_error = false)
  end

  if recompile_flag == true
    FType = typeof(f)
    SolType = typeof(sol)
    cacheType = typeof(cache)
  else
    FType = Function
    SolType = AbstractODESolution
    cacheType =  OrdinaryDiffEqCache
  end

  tprev = t
  dtcache = tType(dt)
  dtpropose = tType(dt)
  iter = 0
  saveiter = 1 # Starts at 1 so first save is at 2
  saveiter_dense = 1
  kshortsize = 1
  reeval_fsal = false
  u_modified = false
  EEst = tTypeNoUnits(1)
  just_hit_tstop = false
  isout = false
  accept_step = false
  dtchangeable = isdtchangeable(alg)
  q11 = tTypeNoUnits(1)

  integrator = ODEIntegrator{algType,uType,tType,
                             tTypeNoUnits,typeof(tdir),eltype(ks),SolType,
                             typeof(rate_prototype),FType,typeof(prog),cacheType,
                             typeof(opts)}(
                             sol,u,k,t,tType(dt),f,uprev,uprev2,tprev,
                             alg,rate_prototype,notsaveat_idxs,dtcache,dtchangeable,
                             dtpropose,tdir,EEst,qoldinit,q11,
                             iter,saveiter,saveiter_dense,prog,cache,
                             kshortsize,just_hit_tstop,accept_step,isout,reeval_fsal,u_modified,opts)
  if initialize_integrator
    initialize!(integrator,integrator.cache)
    initialize!(callbacks_internal,t,u,integrator)
  end
  integrator
end

function solve!(integrator::ODEIntegrator)
  @inbounds while !isempty(integrator.opts.tstops)
    while integrator.tdir*integrator.t < integrator.tdir*top(integrator.opts.tstops)
      loopheader!(integrator)
      @ode_exit_conditions
      perform_step!(integrator,integrator.cache)
      loopfooter!(integrator)
      if isempty(integrator.opts.tstops)
        break
      end
    end
    handle_tstop!(integrator)
  end
  postamble!(integrator)
  if typeof(integrator.sol.prob) <: AbstractODETestProblem
    calculate_solution_errors!(integrator.sol;timeseries_errors=integrator.opts.timeseries_errors,dense_errors=integrator.opts.dense_errors)
  end
  integrator.sol.retcode = :Success
  nothing
end


text:
using ChaosTools
using Test

println("\nTesting type stability...")
@testset "Type stability" begin
for ds ∈ [Systems.towel(), Systems.lorenz()]

    pinteg = parallel_integrator(ds, [get_state(ds), get_state(ds)+1e-9])
    @test_nowarn @inferred lyapunov(pinteg, 1000, 100, 1, 1e-9, 1e-6, 1e-12)
    @test_nowarn @inferred lyapunov(ds, 1000)

    tinteg = tangent_integrator(ds, 3)
    @test_nowarn @inferred lyapunovspectrum(tinteg, 1000, 1, 100, Val{3}())
    @test_nowarn @inferred lyapunovspectrum(ds, 1000)

    tinteg = tangent_integrator(ds, 2)
    @test_nowarn @inferred ChaosTools.gali(tinteg, 1000, 1, 1e-12)
    @test_nowarn @inferred gali(ds, 2, 1000)

end
end


text:
module OopsMorePermutations

using Distributions
using LinearAlgebra: I

# TODO: inverse permutation
# TODO: permutation distance
# TODO: make Permutation extend AbstractArray in order to write
#  v[p] instead of v[p.value]

export Permutation, to_matrix

@doc """
If p is a Permutation instance and v is some vector,
then v[p.value] or permute(v, p) returns the permuted vector whereas
permute!(v, p.value) or permute!(v, p) permutes the
vector in-place.
"""
struct Permutation{T<:Integer}
    value::Vector{T}

    function Permutation{T}(value::Vector{T}) where {T<:Integer}
        if sort(value) == 1:length(value)
            return new(value)
        else
            error("Permutation's value must be a permutation of 1,2,...,N")
        end
    end
end

Permutation(value::Array{T}) where {T<:Integer} = Permutation{T}(value);

permute!(v::Vector, p::Permutation{T}) where {T<:Integer} = permute!(v, p.value)

permute(v::Vector, p::Permutation{T}) where {T<:Integer} = v[p.value]

@doc """
Returns the permutation p1∘p2 obtaining by first applying p2
and then p1.
"""
function Base.:*(p1::Permutation{T}, p2::Permutation{T}) where T <: Integer
    return Permutation(p1.value[p2.value])
end

Base.length(p::Permutation{T}) where {T<:Integer} = length(p.value)

function to_matrix(p::Permutation{T} where T <: Integer)
    N = length(p)
    return [Int8(j == p.value[i]) for i=1:N, j=1:N]
end

function from_matrix(m::Array{T,2} where T <: Real)
    if (all(x in [0,1] for x in m) && all(s == 1 for s in sum(m, dims=1))
                                   && all(s == 1 for s in sum(m, dims=2)))
        p::Vector{Int} = m * Vector(1:size(m)[1])
        return Permutation(p)
    else
        error("Not a permutation matrix")
    end
end

end


text:
import accrue.cryptoerase.runtime.Condition;

class C {
    Condition c;
    final int{L /c H} x = 7;
    void m() {
	this.x = 7;
    }
    public static void main(String[] args) {
	new C().m();
    }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
using Oceananigans.Grids: show_domain

"""
    Cell

A type describing the location at the center of a grid cell.
"""
struct Cell end

"""
	Face

A type describing the location at the face of a grid cell.
"""
struct Face end

"""
    Field{X, Y, Z, A, G} <: AbstractLocatedField{X, Y, Z, A, G}

A field defined at the location (`X`, `Y`, `Z`) which can be either `Cell` or `Face`.
"""
struct Field{X, Y, Z, A, G} <: AbstractLocatedField{X, Y, Z, A, G}
    data :: A
    grid :: G
    function Field{X, Y, Z}(data, grid) where {X, Y, Z}
        return new{X, Y, Z, typeof(data), typeof(grid)}(data, grid)
    end
end

"""
	Field(L::Tuple, data::AbstractArray, grid)

Construct a `Field` on `grid` using the array `data` with location defined by the tuple `L`
of length 3 whose elements are `Cell` or `Face`.
"""
Field(L::Tuple, data::AbstractArray, grid) = Field{L[1], L[2], L[3]}(data, grid)

"""
    Field(L::Tuple, arch::AbstractArchitecture, grid)

Construct a `Field` on architecture `arch` and `grid` at location `L`,
where `L` is a tuple of `Cell` or `Face` types.
"""
Field(L::Tuple, arch::AbstractArchitecture, grid) =
    Field{L[1], L[2], L[3]}(zeros(arch, grid), grid)

"""
    Field(X, Y, Z, arch::AbstractArchitecture, grid)

Construct a `Field` on architecture `arch` and `grid` at location `X`, `Y`, `Z`,
where each of `X, Y, Z` is `Cell` or `Face`.
"""
Field(X, Y, Z, arch::AbstractArchitecture, grid) =  Field((X, Y, Z), arch, grid)

"""
    CellField([T=eltype(grid)], arch, grid)

Return a `Field{Cell, Cell, Cell}` on architecture `arch` and `grid`.
Used for tracers and pressure fields.
"""
CellField(T, arch, grid) = Field{Cell, Cell, Cell}(zeros(T, arch, grid), grid)

"""
    FaceFieldX([T=eltype(grid)], arch, grid)

Return a `Field{Face, Cell, Cell}` on architecture `arch` and `grid`.
Used for the x-velocity field.
"""
FaceFieldX(T, arch, grid) = Field{Face, Cell, Cell}(zeros(T, arch, grid), grid)

"""
    FaceFieldY([T=eltype(grid)], arch, grid)

Return a `Field{Cell, Face, Cell}` on architecture `arch` and `grid`.
Used for the y-velocity field.
"""
FaceFieldY(T, arch, grid) = Field{Cell, Face, Cell}(zeros(T, arch, grid), grid)

"""
    FaceFieldZ([T=eltype(grid)], arch, grid)

Return a `Field{Cell, Cell, Face}` on architecture `arch` and `grid`.
Used for the z-velocity field.
"""
FaceFieldZ(T, arch, grid) = Field{Cell, Cell, Face}(zeros(T, arch, grid), grid)

 CellField(arch, grid) = Field((Cell, Cell, Cell), arch, grid)
FaceFieldX(arch, grid) = Field((Face, Cell, Cell), arch, grid)
FaceFieldY(arch, grid) = Field((Cell, Face, Cell), arch, grid)
FaceFieldZ(arch, grid) = Field((Cell, Cell, Face), arch, grid)

location(a) = nothing
location(::AbstractLocatedField{X, Y, Z}) where {X, Y, Z} = (X, Y, Z)

architecture(f::Field) = architecture(f.data)
architecture(o::OffsetArray) = architecture(o.parent)

@inline size(f::AbstractField) = size(f.grid)
@inline length(f::Field) = length(f.data)

@propagate_inbounds getindex(f::Field, inds...) = @inbounds getindex(f.data, inds...)
@propagate_inbounds setindex!(f::Field, v, inds...) = @inbounds setindex!(f.data, v, inds...)
@inline lastindex(f::Field) = lastindex(f.data)
@inline lastindex(f::Field, dim) = lastindex(f.data, dim)

"Returns `f.data` for `f::Field` or `f` for `f::AbstractArray."
@inline data(a) = a
@inline data(f::Field) = f.data

# Endpoint for recursive `datatuple` function:
@inline datatuple(obj::AbstractField) = data(obj)

"Returns `f.data.parent` for `f::Field`."
@inline Base.parent(f::Field) = f.data.parent

"Returns a view over the interior points of the `field.data`."
@inline interior(f::Field) = view(f.data, 1:f.grid.Nx, 1:f.grid.Ny, 1:f.grid.Nz)

"Returns a reference to the interior points of `field.data.parent.`"
@inline interiorparent(f::Field) = @inbounds f.data.parent[1+f.grid.Hx:f.grid.Nx+f.grid.Hx,
                                                           1+f.grid.Hy:f.grid.Ny+f.grid.Hy,
                                                           1+f.grid.Hz:f.grid.Nz+f.grid.Hz]

iterate(f::Field, state=1) = iterate(f.data, state)

@inline xnode(::Type{Cell}, i, grid) = @inbounds grid.xC[i]
@inline xnode(::Type{Face}, i, grid) = @inbounds grid.xF[i]

@inline ynode(::Type{Cell}, j, grid) = @inbounds grid.yC[j]
@inline ynode(::Type{Face}, j, grid) = @inbounds grid.yF[j]

@inline znode(::Type{Cell}, k, grid) = @inbounds grid.zC[k]
@inline znode(::Type{Face}, k, grid) = @inbounds grid.zF[k]

@inline xnode(i, ϕ::Field{X, Y, Z}) where {X, Y, Z} = xnode(X, i, ϕ.grid)
@inline ynode(j, ϕ::Field{X, Y, Z}) where {X, Y, Z} = ynode(Y, j, ϕ.grid)
@inline znode(k, ϕ::Field{X, Y, Z}) where {X, Y, Z} = znode(Z, k, ϕ.grid)

xnodes(ϕ::AbstractField) = reshape(ϕ.grid.xC, ϕ.grid.Nx, 1, 1)
ynodes(ϕ::AbstractField) = reshape(ϕ.grid.yC, 1, ϕ.grid.Ny, 1)
znodes(ϕ::AbstractField) = reshape(ϕ.grid.zC, 1, 1, ϕ.grid.Nz)

xnodes(ϕ::Field{Face})                    = reshape(ϕ.grid.xF[1:end-1], ϕ.grid.Nx, 1, 1)
ynodes(ϕ::Field{X, Face}) where X         = reshape(ϕ.grid.yF[1:end-1], 1, ϕ.grid.Ny, 1)
znodes(ϕ::Field{X, Y, Face}) where {X, Y} = reshape(ϕ.grid.zF[1:end-1], 1, 1, ϕ.grid.Nz)

nodes(ϕ) = (xnodes(ϕ), ynodes(ϕ), znodes(ϕ))

# Niceties
const AbstractCPUField =
    AbstractField{A, G} where {A<:OffsetArray{T, D, <:Array} where {T, D}, G}

@hascuda const AbstractGPUField =
    AbstractField{A, G} where {A<:OffsetArray{T, D, <:CuArray} where {T, D}, G}

set!(u::Field, v::Number) = @. u.data = v

set!(u::Field{X, Y, Z, A}, v::Field{X, Y, Z, A}) where {X, Y, Z, A} =
    @. u.data.parent = v.data.parent

"Set the CPU field `u` to the array `v`."
function set!(u::AbstractCPUField, v::Array)
    for k in 1:u.grid.Nz, j in 1:u.grid.Ny, i in 1:u.grid.Nx
        u[i, j, k] = v[i, j, k]
    end
    return nothing
end

# Set the GPU field `u` to the array `v`.
@hascuda function set!(u::AbstractGPUField, v::Array)
    v_field = Field(location(u), CPU(), u.grid)
    set!(v_field, v)
    set!(u, v_field)
    return nothing
end

# Set the GPU field `u` to the CuArray `v`.
@hascuda function set!(u::AbstractGPUField, v::CuArray)
    @launch device(GPU()) config=launch_config(u.grid, 3) _set_gpu!(u.data, v, u.grid)
    return nothing
end

function _set_gpu!(u, v, grid)
    @loop for k in (1:grid.Nz; (blockIdx().z - 1) * blockDim().z + threadIdx().z)
        @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
            @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
	            u[i, j, k] = v[i, j, k]
	        end
	    end
    end
    return nothing
end

# Set the GPU field `u` data to the CPU field data of `v`.
@hascuda set!(u::AbstractGPUField, v::AbstractCPUField) = copyto!(u.data.parent, v.data.parent)

# Set the CPU field `u` data to the GPU field data of `v`.
@hascuda set!(u::AbstractCPUField, v::AbstractGPUField) = u.data.parent .= Array(v.data.parent)

"Set the CPU field `u` data to the function `f(x, y, z)`."
set!(u::AbstractCPUField, f::Function) = interior(u) .= f.(nodes(u)...)

# Set the GPU field `u` data to the function `f(x, y, z)`.
@hascuda function set!(u::AbstractGPUField, f::Function)
    u_cpu = Field(location(u), CPU(), u.grid)
    set!(u_cpu, f)
    set!(u, u_cpu)
    return nothing
end

"""
    set!(model; kwargs...)

Set velocity and tracer fields of `model`. The keyword arguments
`kwargs...` take the form `name=data`, where `name` refers to one of the
fields of `model.velocities` or `model.tracers`, and the `data` may be an array,
a function with arguments `(x, y, z)`, or any data type for which a
`set!(ϕ::AbstractField, data)` function exists.

Example
=======
```julia
model = Model(grid=RegularCartesianGrid(size=(32, 32, 32), length=(1, 1, 1))

# Set u to a parabolic function of z, v to random numbers damped
# at top and bottom, and T to some silly array of half zeros,
# half random numbers.

u₀(x, y, z) = z/model.grid.Lz * (1 + z/model.grid.Lz)
v₀(x, y, z) = 1e-3 * rand() * u₀(x, y, z)

T₀ = rand(size(model.grid)...)
T₀[T₀ .< 0.5] .= 0

set!(model, u=u₀, v=v₀, T=T₀)
```
"""
function set!(model::AbstractModel; kwargs...)
    for (fldname, value) in kwargs
        if fldname ∈ propertynames(model.velocities)
            ϕ = getproperty(model.velocities, fldname)
        elseif fldname ∈ propertynames(model.tracers)
            ϕ = getproperty(model.tracers, fldname)
        else
            throw(ArgumentError("name $fldname not found in model.velocities or model.tracers."))
        end
        set!(ϕ, value)
    end
    return nothing
end

function set!(Φ::NamedTuple; kwargs...)
    for (fldname, value) in kwargs
        ϕ = getproperty(Φ, fldname)
        set!(ϕ, value)
    end
    return nothing
end

set_ic!(model; kwargs...) = set!(model; kwargs...) # legacy wrapper

Adapt.adapt_structure(to, field::Field{X, Y, Z}) where {X, Y, Z} =
    Field{X, Y, Z}(adapt(to, data), field.grid)

show_location(X, Y, Z) = string("(", string(typeof(X())), ", ",
                                     string(typeof(Y())), ", ",
                                     string(typeof(Z())), ")")

show_location(field::AbstractLocatedField{X, Y, Z}) where {X, Y, Z} = show_location(X, Y, Z)

short_show(a) = string(typeof(a))
shortname(a::Array) = string(typeof(a).name.wrapper)

show(io::IO, field::Field) =
    print(io,
          short_show(field), '\n',
          "├── data: ", typeof(field.data), '\n',
          "└── grid: ", typeof(field.grid), '\n',
          "    ├── size: ", size(field.grid), '\n',
          "    └── domain: ", show_domain(field.grid), '\n')

short_show(field::AbstractLocatedField) = string("Field at ", show_location(field))