text:
##################################################################################
#
#  LIFTING TRANSFORMS
#  Periodic boundaries, dyadic length (powers of 2)
#
##################################################################################

for (Xwt) in (:dwt!, :wpt!)
@eval begin
    # pseudo "out of place" by copying
    function $Xwt{T<:FloatingPoint}(y::AbstractArray{T}, x::AbstractArray{T}, scheme::GLS, L::Union(Integer, BitVector), fw::Bool)
        copy!(y, x)
        $Xwt(y, scheme, L, fw)
        return y
    end
end
end

# return scheme parameters adjusted for direction and type
function makescheme{T<:FloatingPoint}(::Type{T}, scheme::GLS, fw::Bool)
    if fw
        return scheme.step, convert(T, scheme.norm1), convert(T, scheme.norm2)
    else
        return reverse(scheme.step), convert(T, 1/scheme.norm1), convert(T, 1/scheme.norm2)
    end
end

# 1-D
# inplace transform of y, no vector allocation
# tmp: size at least n>>2
function dwt!{T<:FloatingPoint}(y::AbstractVector{T}, scheme::GLS, L::Integer, fw::Bool, tmp::Vector{T}=Array(T,length(y)>>2))

    n = length(y)
    J = nscales(n)
    @assert isdyadic(y)
    @assert 0 <= L <= J
    @assert length(tmp) >= n>>2
    L == 0 && return y          # do nothing

    if fw
        jrange = (J-1):-1:(J-L)
        ns = n
        half = ns>>1
    else
        jrange = (J-L):(J-1)
        ns = 2^(jrange[1]+1)
        half = ns>>1
    end
    s = y
    stepseq, norm1, norm2 = makescheme(T, scheme, fw)

    for j in jrange
        if fw
            split!(s, ns, tmp)
            for step in stepseq
                stepcoef = convert(Array{T}, step.coef)
                if step.stept == 'p'
                    predictfw!(s, half, stepcoef, step.shift)
                elseif step.stept == 'u'
                    updatefw!(s, half, stepcoef, step.shift)
                end
            end
            normalize!(s, half, ns, norm1, norm2)
            ns = ns>>1
            half = half>>1
        else
            normalize!(s, half, ns, norm1, norm2)
            for step in stepseq
                stepcoef = convert(Array{T}, step.coef)
                if step.stept == 'p'
                    predictbw!(s, half, stepcoef, step.shift)
                elseif step.stept == 'u'
                    updatebw!(s, half, stepcoef, step.shift)
                end
            end
            merge!(s, ns, tmp)        # inverse split
            ns = ns<<1
            half = half<<1
        end
    end
    return y
end
# 1-D, unsafe simple 1 level transform
# inplace transform of y, no vector allocation
# tmp: size at least n>>2
# oopc: use oop computation, if false iy and incy are assumed to be 1
# oopv: the out of place location
function unsafe_dwt1level!{T<:FloatingPoint}(y::AbstractArray{T}, iy::Integer, incy::Integer, oopc::Bool, oopv::Vector{T}, scheme::GLS, fw::Bool, stepseq::Vector, norm1::T, norm2::T, tmp::Vector{T})
    if !oopc
        oopv = y
    end
    ns = length(oopv)
    half = ns>>1

    if fw
        if oopc
            split!(oopv, y, iy, incy, ns)
        else
            split!(oopv, ns, tmp)
        end
        for step in stepseq
            stepcoef = convert(Vector{T}, step.coef)
            if step.stept == 'p'
                predictfw!(oopv, half, stepcoef, step.shift)
            else # 'u'
                updatefw!(oopv, half, stepcoef, step.shift)
            end
        end
        if oopc
            normalize!(y, iy, incy, oopv, half, ns, norm1, norm2)
        else
            normalize!(oopv, half, ns, norm1, norm2)
        end
    else
        if oopc
            normalize!(oopv, y, iy, incy, half, ns, norm1, norm2)
        else
            normalize!(oopv, half, ns, norm1, norm2)
        end
        for step in stepseq
            stepcoef = convert(Vector{T}, step.coef)
            if step.stept == 'p'
                predictbw!(oopv, half, stepcoef, step.shift)
            else # 'u'
                updatebw!(oopv, half, stepcoef, step.shift)
            end
        end
        if oopc
            merge!(y, iy, incy, oopv, ns)
        else
            merge!(oopv, ns, tmp)
        end
    end

    return y
end

# 2-D
# inplace transform of y, no vector allocation
# tmp: size at least n>>2
# tmpvec: size at least n
function dwt!{T<:FloatingPoint}(y::Matrix{T}, scheme::GLS, L::Integer, fw::Bool, tmp::Vector{T}=Array(T,size(y,1)>>2), tmpvec::Vector{T}=Array(T,size(y,1)))

    n = size(y,1)
    J = nscales(n)
    @assert iscube(y)
    @assert isdyadic(y)
    @assert 0 <= L <= J
    @assert length(tmp) >= n>>2
    @assert length(tmpvec) >= n
    L == 0 && return y          # do nothing

    if fw
        jrange = (J-1):-1:(J-L)
        nsub = n
    else
        jrange = (J-L):(J-1)
        nsub = int(2^(J-L+1))
    end
    stepseq, norm1, norm2 = makescheme(T, scheme, fw)

    xm = 0
    xs = n
    for j in jrange
        tmpsub = unsafe_vectorslice(tmpvec, 1, nsub)
        if fw
            # rows
            for i=1:nsub
                xi = i
                # out of place in a dense array for speed
                unsafe_dwt1level!(y, xi, xs, true, tmpsub, scheme, fw,
                                    stepseq, norm1, norm2, tmp)
            end
            # columns
            for i=1:nsub
                xi = 1+(i-1)*n
                ya = unsafe_vectorslice(y, xi, nsub)
                unsafe_dwt1level!(ya, 1, 1, false, tmpsub, scheme, fw,
                                    stepseq, norm1, norm2, tmp)
            end
        else
            # columns
            for i=1:nsub
                xi = 1+(i-1)*n
                ya = unsafe_vectorslice(y, xi, nsub)
                unsafe_dwt1level!(ya, 1, 1, false, tmpsub, scheme, fw,
                                    stepseq, norm1, norm2, tmp)
            end
            # rows
            for i=1:nsub
                xi = i
                # out of place in a dense array for speed
                unsafe_dwt1level!(y, xi, xs, true, tmpsub, scheme, fw,
                                    stepseq, norm1, norm2, tmp)
            end
        end

        fw  && (nsub = nsub>>1)
        !fw && (nsub = nsub<<1)
    end

    return y
end

# WPT
# 1-D
function wpt!{T<:FloatingPoint}(y::AbstractVector{T}, scheme::GLS, L::Integer, fw::Bool, tmp::Vector{T}=Array(T,length(y)>>2))
    wpt!(y, scheme, maketree(length(y), L, :full), fw, tmp)
end
function wpt!{T<:FloatingPoint}(y::AbstractVector{T}, scheme::GLS, tree::BitVector, fw::Bool, tmp::Vector{T}=Array(T,length(y)>>2))

    n = length(y)
    J = nscales(n)
    @assert isdyadic(y)
    @assert isvalidtree(y, tree)
    @assert length(tmp) >= n>>2
    tree[1] || return y          # do nothing

    stepseq, norm1, norm2 = makescheme(T, scheme, fw)

    L = J
    while L > 0
        ix = 1
        k = 1
        fw  && (Lfw = J-L)
        !fw && (Lfw = L-1)
        nj = detailn(tl2level(n, Lfw))
        treeind = 2^(Lfw)-1

        while ix <= n
            if tree[treeind+k]
                dy = unsafe_vectorslice(y, ix, nj)
                unsafe_dwt1level!(dy, 1, 1, false, tmp, scheme, fw, stepseq, norm1, norm2, tmp)
            end
            ix += nj
            k += 1
        end
        L -= 1
    end

    return y
end

function normalize!{T<:FloatingPoint}(x::AbstractVector{T}, half::Int, ns::Int, n1::T, n2::T)
    for i = 1:half
        @inbounds x[i] *= n1
    end
    for i = half+1:ns
        @inbounds x[i] *= n2
    end
    return x
end
# out of place normalize from x to y
function normalize!{T<:FloatingPoint}(y::AbstractVector{T}, x::AbstractVector{T}, half::Int, ns::Int, n1::T, n2::T)
    for i = 1:half
        @inbounds y[i] = n1*x[i]
    end
    for i = half+1:ns
        @inbounds y[i] = n2*x[i]
    end
    return y
end
function normalize!{T<:FloatingPoint}(y::AbstractArray{T}, iy::Int, incy::Int, x::AbstractVector{T}, half::Int, ns::Int, n1::T, n2::T)
    for i = 1:half
        @inbounds y[iy + (i-1)*incy] = n1*x[i]
    end
    for i = half+1:ns
        @inbounds y[iy + (i-1)*incy] = n2*x[i]
    end
    return y
end
function normalize!{T<:FloatingPoint}(y::AbstractVector{T}, x::AbstractArray{T}, ix::Int, incx::Int, half::Int, ns::Int, n1::T, n2::T)
    for i = 1:half
        @inbounds y[i] = n1*x[ix + (i-1)*incx]
    end
    for i = half+1:ns
        @inbounds y[i] = n2*x[ix + (i-1)*incx]
    end
    return y
end

# predict and update lifting steps inplace on x, forward and backward
# half: half of the length under consideration, shift: shift to left, c: lift coefs
# For predict: writes to range 1:half, reads from 1:2*half
# For update : writes to range half+1:2*half, reads from 1:2*half
for (fname,lift_inb,lift_bound,pred) in (
                        (:predictfw!,:lift_inboundsfw!, :liftp_perboundaryfw!, true),
                        (:predictbw!,:lift_inboundsbw!, :liftp_perboundarybw!, true),
                        (:updatefw!, :lift_inboundsfw!, :liftu_perboundaryfw!, false),
                        (:updatebw!, :lift_inboundsbw!, :liftu_perboundarybw!, false) )
@eval begin
function ($fname){T<:FloatingPoint}(x::AbstractVector{T}, half::Int, c::Vector{T}, shift::Int)
    lhsr, irange, rhsr, rhsis = getliftranges(half, length(c), shift, $pred)

    # left boundary
    ($lift_bound)(x, half, c, lhsr, rhsis)
    # main loop
    ($lift_inb)(x, c, irange, rhsis)
    # right boundary
    ($lift_bound)(x, half, c, rhsr, rhsis)
    return x
end
end # eval begin
end # for

function getliftranges(half::Int, nc::Int, shift::Int, pred::Bool)
    # define index shift rhsis
    if pred
    	rhsis = -shift+half
    else
        rhsis = -shift-half
    end
    # conditions for every element i in irange to be in bounds
    # 1 <= i <= half
    # 1 <= i+1-1-shift <= half
    # 1 <= i+nc-1-shift <= half
    irmin = max(shift+1, 1-nc+shift)
    irmax = min(half+1+shift-nc, half+shift)
    if irmin > half || irmax < 1
        irange = 1:0  # empty
    else
        irmin = max(irmin,1)
        irmax = min(irmax,half)
        irange = irmin:irmax
    end
    # periodic boundary
    if length(irange)==0
        lhsr = 1:half
        rhsr = 1:0
    else
        lhsr = 1:irmin-1
        rhsr = irmax+1:half
    end
    if !(pred)  # shift ranges for update
    	irange += half
    	lhsr += half
    	rhsr += half
    end
    return (lhsr, irange, rhsr, rhsis)
end

# periodic boundary
for (fname,op,puxind) in (  (:liftp_perboundaryfw!,:-,:(mod1(i+k-1+rhsis-half,half)+half)),
                            (:liftp_perboundarybw!,:+,:(mod1(i+k-1+rhsis-half,half)+half)),
                            (:liftu_perboundaryfw!, :-,:(mod1(i+k-1+rhsis,half))),
                            (:liftu_perboundarybw!, :+,:(mod1(i+k-1+rhsis,half)))
                            )
@eval begin
function ($fname){T<:FloatingPoint}(x::AbstractVector{T}, half::Int, c::Vector{T}, irange::Range, rhsis::Int)
    nc = length(c)
    for i in irange
        for k = 1:nc
            @inbounds x[i] = ($op)(x[i], c[k]*x[$puxind] )
        end
    end
    return x
end
end # eval begin
end # for

# main lift loop
for (fname,op) in ( (:lift_inboundsfw!,:-), (:lift_inboundsbw!,:+,) )
@eval begin
function ($fname){T<:FloatingPoint}(x::AbstractVector{T}, c::Vector{T}, irange::Range, rhsis::Int)
    nc = length(c)
    if nc == 1  # hard code the most common cases (1, 2, 3) for speed
        c1 = c[1]
        for i in irange
            @inbounds x[i] = ($op)(x[i], c1*x[i+rhsis] )
        end
    elseif nc == 2
        c1,c2 = c[1],c[2]
        rhsisp1 = rhsis+1
        for i in irange
            @inbounds x[i] = ($op)(x[i], c1*x[i+rhsis] + c2*x[i+rhsisp1] )
        end
    elseif nc == 3
        c1,c2,c3 = c[1],c[2],c[3]
        rhsisp1 = rhsis+1
        rhsisp2 = rhsis+2
        for i = irange
            @inbounds x[i] = ($op)(x[i], c1*x[i+rhsis] + c2*x[i+rhsisp1] + c3*x[i+rhsisp2] )
        end
    else
        for i in irange
            for k = 0:nc-1
                @inbounds x[i] = ($op)(x[i], c[k]*x[i+k+rhsis] )
            end
        end
    end
    return x
end
end # eval begin
end # for


text:
# # CUDA acceleration of Quantum Simulation

# For example, a SWAP gate can be instructed on a GPU register like

# ```julia
# function instruct!(reg::GPUReg, ::Val{:SWAP}, locs::Tuple{Int,Int})
#     b1, b2 = locs
#     state = statevec(reg)
#     mask1 = bmask(b1)
#     mask2 = bmask(b2)

#     function kf(state, mask1, mask2)
#         inds = ((blockIdx().x-1) * blockDim().x + threadIdx().x,
#                        (blockIdx().y-1) * blockDim().y + threadIdx().y)
#         b = inds[1]-1
#         c = inds[2]
#         c <= size(state, 2) || return nothing
#         if b&mask1==0 && b&mask2==mask2
#             i = b+1
#             i_ = b  (mask1|mask2) + 1
#             temp = state[i, c]
#             state[i, c] = state[i_, c]
#             state[i_, c] = temp
#         end
#         nothing
#     end
#     X, Y = cudiv(size(state)...)
#     @cuda threads=X blocks=Y kf(state, mask1, mask2)
#     state
# end
# ```

# Here, we devide the threads and blocks into a two dimensional grid with a same shape as the input GPUReg storage (i.e. ``2^a\times 2^rB``).
# Only if two qubits at `locs` are $0$ and $1$ respectively, they are exchanged, otherwise do nothing.
# Although $3/4$ of threads are idle and plenty room for optimization, from this example, we see how easy CUDA programming is with
# [CUDAnative](https://github.com/JuliaGPU/CUDAnative.jl).

ext: jl

max_stars_repo_name: logankilpatrick/tutorials

max_issues_repo_name: logankilpatrick/tutorials

max_forks_repo_name: logankilpatrick/tutorials


text:
using Price-Index
using Test

@testset "Price-Index.jl" begin

end


text:
"""
originally from:
https://gist.github.com/RobBlackwell/10a1aeabeb85bbf1a17cc334e5e60acf
"""

using PyCall
pickle = pyimport("pickle")

"""
Save pickle from python pickle file.
"""
function save_pickle end

"""
Load pickle from python pickle file.
"""
function load_pickle end

function save_pickle(filename, obj)
    out = open(filename, "w")
    pickle.dump(obj, out)
    close(out)
end

function load_pickle(filename)
    r = nothing
    @pywith pybuiltin("open")(filename, "rb") as f begin
        r = pickle.load(f)
    end
    return r
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

let
    function foo_error(c, n)
        c[1] += 1
        if c[1] <= n
            error("foo")
        end
        return 7
    end

    # Success on first attempt
    c = [0]
    @test retry(foo_error)(c,0) == 7
    @test c[1] == 1

    # Success on second attempt
    c = [0]
    @test retry(foo_error;n=1)(c,1) == 7
    @test c[1] == 2

    # 2 failed retry attempts, so exception is raised
    c = [0]
    ex = try retry(foo_error;n=2)(c,3) catch e; e end
    @test ex.msg == "foo"
    @test c[1] == 3

    c = [0]
    ex = try retry(foo_error, ErrorException)(c,2) catch e; e end
    @test typeof(ex) == ErrorException
    @test ex.msg == "foo"
    @test c[1] == 2

    c = [0]
    ex = try retry(foo_error, e->e.msg == "foo")(c,2) catch e; e end
    @test typeof(ex) == ErrorException
    @test ex.msg == "foo"
    @test c[1] == 2

    # No retry if condition does not match
    c = [0]
    ex = try retry(foo_error, e->e.msg == "bar"; n=3)(c,2) catch e; e end
    @test typeof(ex) == ErrorException
    @test ex.msg == "foo"
    @test c[1] == 1

    c = [0]
    ex = try retry(foo_error, e->e.http_status_code == "503")(c,2) catch e; e end
    @test typeof(ex) == ErrorException
    @test ex.msg == "foo"
    @test c[1] == 1

    c = [0]
    ex = try retry(foo_error, SystemError)(c,2) catch e; e end
    @test typeof(ex) == ErrorException
    @test ex.msg == "foo"
    @test c[1] == 1
end

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
# WARNING: This file was generated from bigbed-parser.rl using ragel. Do not edit!
type BigBedData <: IntervalStream{BEDMetadata}
stream::BufferedInputStream
header::BigBedHeader
zoom_headers::Vector{BigBedZoomHeader}
autosql::AbstractString
summary::BigBedTotalSummary
btree_header::BigBedBTreeHeader
rtree_header::BigBedRTreeHeader
data_count::UInt32

# preallocated space for reading and searching the B-tree
btree_internal_nodes::Vector{BigBedBTreeInternalNode}
btree_leaf_nodes::Vector{BigBedBTreeLeafNode}
key::Vector{UInt8}
node_keys::Vector{Vector{UInt8}}
uncompressed_data::Vector{UInt8}
end

	# Parser for data blocks in a BigBed file. This is very similar
	# to the BED parser in bed.rl, with the following exceptions:
	#
	#    * BigBed has binary chrom_index, start, and end, insteado of ASCII
#      chromosome name, start, end.
#    * BigBed entries are null ('\0') terminated, rather than newline separated.
#
const _bigbedparser_start  = 39
const _bigbedparser_first_final  = 39
const _bigbedparser_error  = 0
const _bigbedparser_en_main  = 39
const __bigbedparser_nfa_targs = Int8[ 0, 0 ,  ]
const __bigbedparser_nfa_offsets = Int8[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ,  ]
const __bigbedparser_nfa_push_actions = Int8[ 0, 0 ,  ]
const __bigbedparser_nfa_pop_trans = Int8[ 0, 0 ,  ]
typealias StringFieldVector Vector{StringField}
typealias NullableStringFieldVector Nullable{StringFieldVector}
typealias NullableStringField Nullable{StringField}

type BigBedDataParser <: AbstractParser
state::Ragel.State

# intermediate values used during parsing
chrom_id::UInt32
red::Float32
green::Float32
blue::Float32
block_size_idx::Int
block_first_idx::Int
seq_names::Nullable{StringFieldVector}
assumed_seqname::Nullable{StringField}

function BigBedDataParser(input::BufferedInputStream;
seq_names::NullableStringFieldVector=NullableStringFieldVector(),
assumed_seqname::NullableStringField=NullableStringField())

begin
cs = convert( Int , _bigbedparser_start );

end
return new(Ragel.State(cs, input), 0, 0.0, 0.0, 0.0, 1, 1, seq_names, assumed_seqname)
end
end

Ragel.@generate_read!_function("_bigbedparser", BigBedDataParser, BEDInterval,
begin
begin
if ( p == pe  )
@goto _test_eof

end
if ( cs  == 39 )
@goto st_case_39
elseif ( cs  == 1 )
@goto st_case_1
elseif ( cs  == 2 )
@goto st_case_2
elseif ( cs  == 3 )
@goto st_case_3
elseif ( cs  == 4 )
@goto st_case_4
elseif ( cs  == 5 )
@goto st_case_5
elseif ( cs  == 6 )
@goto st_case_6
elseif ( cs  == 7 )
@goto st_case_7
elseif ( cs  == 8 )
@goto st_case_8
elseif ( cs  == 9 )
@goto st_case_9
elseif ( cs  == 10 )
@goto st_case_10
elseif ( cs  == 11 )
@goto st_case_11
elseif ( cs  == 12 )
@goto st_case_12
elseif ( cs  == 40 )
@goto st_case_40
elseif ( cs  == 0 )
@goto st_case_0
elseif ( cs  == 13 )
@goto st_case_13
elseif ( cs  == 14 )
@goto st_case_14
elseif ( cs  == 15 )
@goto st_case_15
elseif ( cs  == 16 )
@goto st_case_16
elseif ( cs  == 17 )
@goto st_case_17
elseif ( cs  == 18 )
@goto st_case_18
elseif ( cs  == 19 )
@goto st_case_19
elseif ( cs  == 20 )
@goto st_case_20
elseif ( cs  == 21 )
@goto st_case_21
elseif ( cs  == 22 )
@goto st_case_22
elseif ( cs  == 23 )
@goto st_case_23
elseif ( cs  == 24 )
@goto st_case_24
elseif ( cs  == 25 )
@goto st_case_25
elseif ( cs  == 26 )
@goto st_case_26
elseif ( cs  == 27 )
@goto st_case_27
elseif ( cs  == 28 )
@goto st_case_28
elseif ( cs  == 29 )
@goto st_case_29
elseif ( cs  == 30 )
@goto st_case_30
elseif ( cs  == 31 )
@goto st_case_31
elseif ( cs  == 32 )
@goto st_case_32
elseif ( cs  == 33 )
@goto st_case_33
elseif ( cs  == 34 )
@goto st_case_34
elseif ( cs  == 35 )
@goto st_case_35
elseif ( cs  == 36 )
@goto st_case_36
elseif ( cs  == 37 )
@goto st_case_37
elseif ( cs  == 38 )
@goto st_case_38
end
@goto st_out
@label st_case_39
begin
@goto ctr64

end
@label ctr64
begin
Ragel.@anchor!
end
@goto st1
@label ctr65
begin
input.block_size_idx = 1
input.block_first_idx = 1

# // fbreak causes will cause the pushmark action for the next seqname
# // to be skipped, so we do it here
Ragel.@anchor!
Ragel.@yield 1
end
begin
Ragel.@anchor!
end
@goto st1
@label st1
p+= 1;
if ( p == pe  )
@goto _test_eof1

end
@label st_case_1
begin
@goto st2

end
@label st2
p+= 1;
if ( p == pe  )
@goto _test_eof2

end
@label st_case_2
begin
@goto st3

end
@label st3
p+= 1;
if ( p == pe  )
@goto _test_eof3

end
@label st_case_3
begin
@goto st4

end
@label st4
p+= 1;
if ( p == pe  )
@goto _test_eof4

end
@label st_case_4
begin
@goto ctr3

end
@label ctr3
begin
input.chrom_id = Ragel.@load_from_anchor!(UInt32)
if !isnull(input.seq_names)
output.seqname = copy(get(input.seq_names)[input.chrom_id + 1])
else
output.seqname = copy(get(input.assumed_seqname))
end
output.metadata.used_fields = 0
end
begin
Ragel.@anchor!
end
@goto st5
@label st5
p+= 1;
if ( p == pe  )
@goto _test_eof5

end
@label st_case_5
begin
@goto st6

end
@label st6
p+= 1;
if ( p == pe  )
@goto _test_eof6

end
@label st_case_6
begin
@goto st7

end
@label st7
p+= 1;
if ( p == pe  )
@goto _test_eof7

end
@label st_case_7
begin
@goto st8

end
@label st8
p+= 1;
if ( p == pe  )
@goto _test_eof8

end
@label st_case_8
begin
@goto ctr7

end
@label ctr7
begin
output.first = 1 + Ragel.@load_from_anchor!(UInt32)
end
begin
Ragel.@anchor!
end
@goto st9
@label st9
p+= 1;
if ( p == pe  )
@goto _test_eof9

end
@label st_case_9
begin
@goto st10

end
@label st10
p+= 1;
if ( p == pe  )
@goto _test_eof10

end
@label st_case_10
begin
@goto st11

end
@label st11
p+= 1;
if ( p == pe  )
@goto _test_eof11

end
@label st_case_11
begin
@goto st12

end
@label st12
p+= 1;
if ( p == pe  )
@goto _test_eof12

end
@label st_case_12
if ( (data[1+(p )])== 0  )
begin
@goto ctr11

end

end
if ( 32 <= (data[1+(p )])&& (data[1+(p )])<= 126  )
begin
@goto ctr12

end

end
begin
@goto st0

end
@label ctr11
begin
output.last = Ragel.@load_from_anchor!(UInt32)
end
@goto st40
@label ctr14
begin
Ragel.@copy_from_anchor!(output.metadata.name)
end
begin
output.metadata.used_fields += 1
end
@goto st40
@label ctr18
begin
output.metadata.score = Ragel.@int64_from_anchor!
end
begin
output.metadata.used_fields += 1
end
@goto st40
@label ctr22
begin
output.metadata.used_fields += 1
end
@goto st40
@label ctr25
begin
output.metadata.thick_first = 1 + Ragel.@int64_from_anchor!
end
begin
output.metadata.used_fields += 1
end
@goto st40
@label ctr29
begin
output.metadata.thick_last = Ragel.@int64_from_anchor!
end
begin
output.metadata.used_fields += 1
end
@goto st40
@label ctr33
begin
input.red = input.green = input.blue = (Ragel.@int64_from_anchor!) / 255.0
end
begin
output.metadata.item_rgb = RGB{Float32}(input.red, input.green, input.blue)
end
begin
output.metadata.used_fields += 1
end
@goto st40
@label ctr48
begin
input.blue = (Ragel.@int64_from_anchor!) / 255.0
end
begin
output.metadata.item_rgb = RGB{Float32}(input.red, input.green, input.blue)
end
begin
output.metadata.used_fields += 1
end
@goto st40
@label ctr51
begin
output.metadata.block_count = Ragel.@int64_from_anchor!

if (output.metadata.block_count > length(output.metadata.block_sizes))
resize!(output.metadata.block_sizes, output.metadata.block_count)
end

if (output.metadata.block_count > length(output.metadata.block_firsts))
resize!(output.metadata.block_firsts, output.metadata.block_count)
end
end
begin
output.metadata.used_fields += 1
end
@goto st40
@label ctr55
begin
if input.block_size_idx > length(output.metadata.block_sizes)
error("More size blocks encountered than BED block count field suggested.")
end
output.metadata.block_sizes[input.block_size_idx] = Ragel.@int64_from_anchor!
input.block_size_idx += 1
end
begin
output.metadata.used_fields += 1
end
@goto st40
@label ctr60
begin
if input.block_first_idx > length(output.metadata.block_firsts)
error("More start blocks encountered than BED block count field suggested.")
end
output.metadata.block_firsts[input.block_first_idx] = 1 + Ragel.@int64_from_anchor!
input.block_first_idx += 1
end
begin
output.metadata.used_fields += 1
end
@goto st40
@label st40
p+= 1;
if ( p == pe  )
@goto _test_eof40

end
@label st_case_40
begin
@goto ctr65

end
@label st_case_0
@label st0
cs = 0;
@goto _out
@label ctr12
begin
output.last = Ragel.@load_from_anchor!(UInt32)
end
begin
Ragel.@anchor!
end
@goto st13
@label st13
p+= 1;
if ( p == pe  )
@goto _test_eof13

end
@label st_case_13
if ( (data[1+(p )]) == 0 )
begin
@goto ctr14

end
elseif ( (data[1+(p )]) == 9 )
begin
@goto ctr15

end
end
if ( 32 <= (data[1+(p )])&& (data[1+(p )])<= 126  )
begin
@goto st13

end

end
begin
@goto st0

end
@label ctr15
begin
Ragel.@copy_from_anchor!(output.metadata.name)
end
begin
output.metadata.used_fields += 1
end
@goto st14
@label st14
p+= 1;
if ( p == pe  )
@goto _test_eof14

end
@label st_case_14
if ( 48 <= (data[1+(p )])&& (data[1+(p )])<= 57  )
begin
@goto ctr17

end

end
begin
@goto st0

end
@label ctr17
begin
Ragel.@anchor!
end
@goto st15
@label st15
p+= 1;
if ( p == pe  )
@goto _test_eof15

end
@label st_case_15
if ( (data[1+(p )]) == 0 )
begin
@goto ctr18

end
elseif ( (data[1+(p )]) == 9 )
begin
@goto ctr19

end
end
if ( 48 <= (data[1+(p )])&& (data[1+(p )])<= 57  )
begin
@goto st15

end

end
begin
@goto st0

end
@label ctr19
begin
output.metadata.score = Ragel.@int64_from_anchor!
end
begin
output.metadata.used_fields += 1
end
@goto st16
@label st16
p+= 1;
if ( p == pe  )
@goto _test_eof16

end
@label st_case_16
if ( (data[1+(p )]) == 43 )
begin
@goto ctr21

end
elseif ( (data[1+(p )]) == 63 )
begin
@goto ctr21

end
end
if ( 45 <= (data[1+(p )])&& (data[1+(p )])<= 46  )
begin
@goto ctr21

end

end
begin
@goto st0

end
@label ctr21
begin
output.strand = convert(Strand, (Ragel.@char))
end
@goto st17
@label st17
p+= 1;
if ( p == pe  )
@goto _test_eof17

end
@label st_case_17
if ( (data[1+(p )]) == 0 )
begin
@goto ctr22

end
elseif ( (data[1+(p )]) == 9 )
begin
@goto ctr23

end
end
begin
@goto st0

end
@label ctr23
begin
output.metadata.used_fields += 1
end
@goto st18
@label st18
p+= 1;
if ( p == pe  )
@goto _test_eof18

end
@label st_case_18
if ( 48 <= (data[1+(p )])&& (data[1+(p )])<= 57  )
begin
@goto ctr24

end

end
begin
@goto st0

end
@label ctr24
begin
Ragel.@anchor!
end
@goto st19
@label st19
p+= 1;
if ( p == pe  )
@goto _test_eof19

end
@label st_case_19
if ( (data[1+(p )]) == 0 )
begin
@goto ctr25

end
elseif ( (data[1+(p )]) == 9 )
begin
@goto ctr26

end
end
if ( 48 <= (data[1+(p )])&& (data[1+(p )])<= 57  )
begin
@goto st19

end

end
begin
@goto st0

end
@label ctr26
begin
output.metadata.thick_first = 1 + Ragel.@int64_from_anchor!
end
begin
output.metadata.used_fields += 1
end
@goto st20
@label st20
p+= 1;
if ( p == pe  )
@goto _test_eof20

end
@label st_case_20
if ( 48 <= (data[1+(p )])&& (data[1+(p )])<= 57  )
begin
@goto ctr28

end

end
begin
@goto st0

end
@label ctr28
begin
Ragel.@anchor!
end
@goto st21
@label st21
p+= 1;
if ( p == pe  )
@goto _test_eof21

end
@label st_case_21
if ( (data[1+(p )]) == 0 )
begin
@goto ctr29

end
elseif ( (data[1+(p )]) == 9 )
begin
@goto ctr30

end
end
if ( 48 <= (data[1+(p )])&& (data[1+(p )])<= 57  )
begin
@goto st21

end

end
begin
@goto st0

end
@label ctr30
begin
output.metadata.thick_last = Ragel.@int64_from_anchor!
end
begin
output.metadata.used_fields += 1
end
@goto st22
@label st22
p+= 1;
if ( p == pe  )
@goto _test_eof22

end
@label st_case_22
if ( 48 <= (data[1+(p )])&& (data[1+(p )])<= 57  )
begin
@goto ctr32

end

end
begin
@goto st0

end
@label ctr32
begin
Ragel.@anchor!
end
@goto st23
@label st23
p+= 1;
if ( p == pe  )
@goto _test_eof23

end
@label st_case_23
if ( (data[1+(p )]) == 0 )
begin
@goto ctr33

end
elseif ( (data[1+(p )]) == 9 )
begin
@goto ctr34

end
elseif ( (data[1+(p )]) == 11 )
begin
@goto ctr35

end
elseif ( (data[1+(p )]) == 32 )
begin
@goto ctr35

end
elseif ( (data[1+(p )]) == 44 )
begin
@goto ctr36

end
end
if ( 48 <= (data[1+(p )])&& (data[1+(p )])<= 57  )
begin
@goto st23

end

end
begin
@goto st0

end
@label ctr34
begin
input.red = input.green = input.blue = (Ragel.@int64_from_anchor!) / 255.0
end
begin
output.metadata.item_rgb = RGB{Float32}(input.red, input.green, input.blue)
end
begin
output.metadata.used_fields += 1
end
@goto st24
@label st24
p+= 1;
if ( p == pe  )
@goto _test_eof24

end
@label st_case_24
if ( (data[1+(p )]) == 9 )
begin
@goto st25

end
elseif ( (data[1+(p )]) == 11 )
begin
@goto st25

end
elseif ( (data[1+(p )]) == 32 )
begin
@goto st25

end
elseif ( (data[1+(p )]) == 44 )
begin
@goto st26

end
end
if ( 48 <= (data[1+(p )])&& (data[1+(p )])<= 57  )
begin
@goto ctr40

end

end
begin
@goto st0

end
@label ctr35
begin
input.red = input.green = input.blue = (Ragel.@int64_from_anchor!) / 255.0
end
@goto st25
@label st25
p+= 1;
if ( p == pe  )
@goto _test_eof25

end
@label st_case_25
if ( (data[1+(p )]) == 9 )
begin
@goto st25

end
elseif ( (data[1+(p )]) == 11 )
begin
@goto st25

end
elseif ( (data[1+(p )]) == 32 )
begin
@goto st25

end
elseif ( (data[1+(p )]) == 44 )
begin
@goto st26

end
end
begin
@goto st0

end
@label ctr36
begin
input.red = input.green = input.blue = (Ragel.@int64_from_anchor!) / 255.0
end
@goto st26
@label st26
p+= 1;
if ( p == pe  )
@goto _test_eof26

end
@label st_case_26
if ( (data[1+(p )]) == 9 )
begin
@goto st26

end
elseif ( (data[1+(p )]) == 11 )
begin
@goto st26

end
elseif ( (data[1+(p )]) == 32 )
begin
@goto st26

end
end
if ( 48 <= (data[1+(p )])&& (data[1+(p )])<= 57  )
begin
@goto ctr41

end

end
begin
@goto st0

end
@label ctr41
begin
Ragel.@anchor!
end
@goto st27
@label st27
p+= 1;
if ( p == pe  )
@goto _test_eof27

end
@label st_case_27
if ( (data[1+(p )]) == 9 )
begin
@goto ctr42

end
elseif ( (data[1+(p )]) == 11 )
begin
@goto ctr42

end
elseif ( (data[1+(p )]) == 32 )
begin
@goto ctr42

end
elseif ( (data[1+(p )]) == 44 )
begin
@goto ctr43

end
end
if ( 48 <= (data[1+(p )])&& (data[1+(p )])<= 57  )
begin
@goto st27

end

end
begin
@goto st0

end
@label ctr42
begin
input.green = (Ragel.@int64_from_anchor!) / 255.0
end
@goto st28
@label st28
p+= 1;
if ( p == pe  )
@goto _test_eof28

end
@label st_case_28
if ( (data[1+(p )]) == 9 )
begin
@goto st28

end
elseif ( (data[1+(p )]) == 11 )
begin
@goto st28

end
elseif ( (data[1+(p )]) == 32 )
begin
@goto st28

end
elseif ( (data[1+(p )]) == 44 )
begin
@goto st29

end
end
begin
@goto st0

end
@label ctr43
begin
input.green = (Ragel.@int64_from_anchor!) / 255.0
end
@goto st29
@label st29
p+= 1;
if ( p == pe  )
@goto _test_eof29

end
@label st_case_29
if ( (data[1+(p )]) == 9 )
begin
@goto st29

end
elseif ( (data[1+(p )]) == 11 )
begin
@goto st29

end
elseif ( (data[1+(p )]) == 32 )
begin
@goto st29

end
end
if ( 48 <= (data[1+(p )])&& (data[1+(p )])<= 57  )
begin
@goto ctr47

end

end
begin
@goto st0

end
@label ctr47
begin
Ragel.@anchor!
end
@goto st30
@label st30
p+= 1;
if ( p == pe  )
@goto _test_eof30

end
@label st_case_30
if ( (data[1+(p )]) == 0 )
begin
@goto ctr48

end
elseif ( (data[1+(p )]) == 9 )
begin
@goto ctr49

end
end
if ( 48 <= (data[1+(p )])&& (data[1+(p )])<= 57  )
begin
@goto st30

end

end
begin
@goto st0

end
@label ctr49
begin
input.blue = (Ragel.@int64_from_anchor!) / 255.0
end
begin
output.metadata.item_rgb = RGB{Float32}(input.red, input.green, input.blue)
end
begin
output.metadata.used_fields += 1
end
@goto st31
@label st31
p+= 1;
if ( p == pe  )
@goto _test_eof31

end
@label st_case_31
if ( 48 <= (data[1+(p )])&& (data[1+(p )])<= 57  )
begin
@goto ctr40

end

end
begin
@goto st0

end
@label ctr40
begin
Ragel.@anchor!
end
@goto st32
@label st32
p+= 1;
if ( p == pe  )
@goto _test_eof32

end
@label st_case_32
if ( (data[1+(p )]) == 0 )
begin
@goto ctr51

end
elseif ( (data[1+(p )]) == 9 )
begin
@goto ctr52

end
end
if ( 48 <= (data[1+(p )])&& (data[1+(p )])<= 57  )
begin
@goto st32

end

end
begin
@goto st0

end
@label ctr52
begin
output.metadata.block_count = Ragel.@int64_from_anchor!

if (output.metadata.block_count > length(output.metadata.block_sizes))
resize!(output.metadata.block_sizes, output.metadata.block_count)
end

if (output.metadata.block_count > length(output.metadata.block_firsts))
resize!(output.metadata.block_firsts, output.metadata.block_count)
end
end
begin
output.metadata.used_fields += 1
end
@goto st33
@label st33
p+= 1;
if ( p == pe  )
@goto _test_eof33

end
@label st_case_33
if ( 48 <= (data[1+(p )])&& (data[1+(p )])<= 57  )
begin
@goto ctr54

end

end
begin
@goto st0

end
@label ctr54
begin
Ragel.@anchor!
end
@goto st34
@label st34
p+= 1;
if ( p == pe  )
@goto _test_eof34

end
@label st_case_34
if ( (data[1+(p )]) == 0 )
begin
@goto ctr55

end
elseif ( (data[1+(p )]) == 9 )
begin
@goto ctr56

end
elseif ( (data[1+(p )]) == 44 )
begin
@goto ctr57

end
end
if ( 48 <= (data[1+(p )])&& (data[1+(p )])<= 57  )
begin
@goto st34

end

end
begin
@goto st0

end
@label ctr63
begin
output.metadata.used_fields += 1
end
@goto st35
@label ctr56
begin
if input.block_size_idx > length(output.metadata.block_sizes)
error("More size blocks encountered than BED block count field suggested.")
end
output.metadata.block_sizes[input.block_size_idx] = Ragel.@int64_from_anchor!
input.block_size_idx += 1
end
begin
output.metadata.used_fields += 1
end
@goto st35
@label st35
p+= 1;
if ( p == pe  )
@goto _test_eof35

end
@label st_case_35
if ( 48 <= (data[1+(p )])&& (data[1+(p )])<= 57  )
begin
@goto ctr59

end

end
begin
@goto st0

end
@label ctr59
begin
Ragel.@anchor!
end
@goto st36
@label st36
p+= 1;
if ( p == pe  )
@goto _test_eof36

end
@label st_case_36
if ( (data[1+(p )]) == 0 )
begin
@goto ctr60

end
elseif ( (data[1+(p )]) == 44 )
begin
@goto ctr61

end
end
if ( 48 <= (data[1+(p )])&& (data[1+(p )])<= 57  )
begin
@goto st36

end

end
begin
@goto st0

end
@label ctr61
begin
if input.block_first_idx > length(output.metadata.block_firsts)
error("More start blocks encountered than BED block count field suggested.")
end
output.metadata.block_firsts[input.block_first_idx] = 1 + Ragel.@int64_from_anchor!
input.block_first_idx += 1
end
@goto st37
@label st37
p+= 1;
if ( p == pe  )
@goto _test_eof37

end
@label st_case_37
if ( (data[1+(p )])== 0  )
begin
@goto ctr22

end

end
if ( 48 <= (data[1+(p )])&& (data[1+(p )])<= 57  )
begin
@goto ctr59

end

end
begin
@goto st0

end
@label ctr57
begin
if input.block_size_idx > length(output.metadata.block_sizes)
error("More size blocks encountered than BED block count field suggested.")
end
output.metadata.block_sizes[input.block_size_idx] = Ragel.@int64_from_anchor!
input.block_size_idx += 1
end
@goto st38
@label st38
p+= 1;
if ( p == pe  )
@goto _test_eof38

end
@label st_case_38
if ( (data[1+(p )]) == 0 )
begin
@goto ctr22

end
elseif ( (data[1+(p )]) == 9 )
begin
@goto ctr63

end
end
if ( 48 <= (data[1+(p )])&& (data[1+(p )])<= 57  )
begin
@goto ctr54

end

end
begin
@goto st0

end
@label st_out
@label _test_eof1
cs = 1;
@goto _test_eof
@label _test_eof2
cs = 2;
@goto _test_eof
@label _test_eof3
cs = 3;
@goto _test_eof
@label _test_eof4
cs = 4;
@goto _test_eof
@label _test_eof5
cs = 5;
@goto _test_eof
@label _test_eof6
cs = 6;
@goto _test_eof
@label _test_eof7
cs = 7;
@goto _test_eof
@label _test_eof8
cs = 8;
@goto _test_eof
@label _test_eof9
cs = 9;
@goto _test_eof
@label _test_eof10
cs = 10;
@goto _test_eof
@label _test_eof11
cs = 11;
@goto _test_eof
@label _test_eof12
cs = 12;
@goto _test_eof
@label _test_eof40
cs = 40;
@goto _test_eof
@label _test_eof13
cs = 13;
@goto _test_eof
@label _test_eof14
cs = 14;
@goto _test_eof
@label _test_eof15
cs = 15;
@goto _test_eof
@label _test_eof16
cs = 16;
@goto _test_eof
@label _test_eof17
cs = 17;
@goto _test_eof
@label _test_eof18
cs = 18;
@goto _test_eof
@label _test_eof19
cs = 19;
@goto _test_eof
@label _test_eof20
cs = 20;
@goto _test_eof
@label _test_eof21
cs = 21;
@goto _test_eof
@label _test_eof22
cs = 22;
@goto _test_eof
@label _test_eof23
cs = 23;
@goto _test_eof
@label _test_eof24
cs = 24;
@goto _test_eof
@label _test_eof25
cs = 25;
@goto _test_eof
@label _test_eof26
cs = 26;
@goto _test_eof
@label _test_eof27
cs = 27;
@goto _test_eof
@label _test_eof28
cs = 28;
@goto _test_eof
@label _test_eof29
cs = 29;
@goto _test_eof
@label _test_eof30
cs = 30;
@goto _test_eof
@label _test_eof31
cs = 31;
@goto _test_eof
@label _test_eof32
cs = 32;
@goto _test_eof
@label _test_eof33
cs = 33;
@goto _test_eof
@label _test_eof34
cs = 34;
@goto _test_eof
@label _test_eof35
cs = 35;
@goto _test_eof
@label _test_eof36
cs = 36;
@goto _test_eof
@label _test_eof37
cs = 37;
@goto _test_eof
@label _test_eof38
cs = 38;
@goto _test_eof
@label _test_eof
begin

end
if ( p == eof  )
begin
if ( cs  == 40 )
begin
input.block_size_idx = 1
input.block_first_idx = 1

# // fbreak causes will cause the pushmark action for the next seqname
# // to be skipped, so we do it here
Ragel.@anchor!
Ragel.@yield 0
end

break;
end

end

end
@label _out
begin

end

end
end)


text:
export @leaf

# from https://thautwarm.github.io/MLStyle.jl/latest/tutorials/capture.html
function capture(template, ex, action)
    let template = Expr(:quote, template)
        quote
            @match $ex begin
                $template => $action
                _         => nothing
            end
        end
    end
end

macro capture(template, ex, action)
    capture(template, ex, action) |> esc
end

"""
    _leaf(expr)

This function is called by the `@leaf` macro.
The function turns:

```julia
@leaf RNorm(m=1,s=0) RealInterval(0.0, 1.0)
```

into

```julia
struct RNorm{P} <: AdvancedProbabilisticCircuits.AbstractLeaf
    scope::Int
    params::P
end
function RNorm(scope::Int, nt::NamedTuple)
    P = typeof(nt)
    return RNorm{P}(scope, nt)
end
RNorm(scope, m, s) = RNorm(scope, (; m, s))
RNorm(scope::Int; m = 1, s = 1) = RNorm(scope, m, s)
(printnodetype(::Type{T}) where T <: RNorm) = RNorm
support(n::RNorm) = n.scope => RealInterval(0.0, 1.0)
````

"""
function _leaf(expr1, expr2)
    @capture $n($(p...)) expr1 begin

        err1str = "Field \'"
        err2str = "\' is missing a default value."

        kwargs_ = Any[]
        args_ = Any[]
        for k in filter(pi -> pi isa Expr, p)
            push!(args_, k.args[1])
            push!(kwargs_, k)
        end

        for k in filter(pi -> pi isa Symbol, p)
            push!(args_, k)
            push!(kwargs_, Expr(:kw, k, :(error($err1str * $k * $err2str))))
        end

        # sanity check
        if !(eval(expr2) isa RealInterval) && !(eval(expr2) isa Vector{Int})
            throw(ErrorException("Unknown support specification!"))
        end

        q = quote
            struct $n{P} <: AdvancedProbabilisticCircuits.AbstractLeaf
                scope::Int
                params::P
            end

            function $n(scope::Int, nt::NamedTuple)
                P = typeof(nt)
                return $n{P}(scope, nt)
            end
            $n(scope, $(args_...)) = $n(scope, (;$(args_...)))
            $n(scope::Int; $(kwargs_...)) = $n(scope, $(args_...))

            printnodetype(::Type{T}) where {T<:$n} = $n

            support(n::$n) = (n.scope => $expr2)
        end
        return q
    end
end

"""
    _leaf(expr)

This function is called by the `@leaf` macro.
The function turns:

```julia
@leaf RNorm(m=1,s=0)
```

into

```julia
struct RNorm{P} <: AdvancedProbabilisticCircuits.AbstractLeaf
    scope::Int
    params::P
end
function RNorm(scope::Int, nt::NamedTuple)
    P = typeof(nt)
    return RNorm{P}(scope, nt)
end
RNorm(scope, m, s) = RNorm(scope, (; m, s))
RNorm(scope::Int; m = 1, s = 1) = RNorm(scope, m, s)
(printnodetype(::Type{T}) where T <: RNorm) = RNorm
````

"""
function _leaf(expr1)
    @capture $n($(p...)) expr1 begin

        err1str = "Field \'"
        err2str = "\' is missing a default value."

        kwargs_ = Any[]
        args_ = Any[]
        for k in filter(pi -> pi isa Expr, p)
            push!(args_, k.args[1])
            push!(kwargs_, k)
        end

        for k in filter(pi -> pi isa Symbol, p)
            push!(args_, k)
            push!(kwargs_, Expr(:kw, k, :(error($err1str * $k * $err2str))))
        end

        q = quote
            struct $n{P} <: AdvancedProbabilisticCircuits.AbstractLeaf
                scope::Int
                params::P
            end

            function $n(scope::Int, nt::NamedTuple)
                P = typeof(nt)
                return $n{P}(scope, nt)
            end
            $n(scope, $(args_...)) = $n(scope, (;$(args_...)))
            $n(scope::Int; $(kwargs_...)) = $n(scope, $(args_...))

            printnodetype(::Type{T}) where {T<:$n} = $n
        end
        return q
    end
end

"""
Macro which generates a leaf node.
Example usage to create a Gaussian leaf with support on the whole real line:

```julia
@leaf Normal(μ = 0.0, σ = 1.0) RealInterval(-Inf, Inf)

# define the logpdf
function logpdf(n::Normal{P}, x::Real) where {P<:NamedTuple{(:μ, :σ)}}
    return -log(n.params.σ) - (x-n.params.μ)^2 / (2 * n.params.σ^2)
end
```

Example usage to create a Bernoulli leaf with support on [0,1]:

```julia
@leaf Bernoulli(p = 0.5) [0, 1]

# define the logpdf
function logpdf(n::Bernoulli{P}, x::Real) where {P<:NamedTuple{(:p)}}
    !(isone(x) | iszero(x)) && return -Inf
    return isone(x) ? log(n.params.p) : log(1-n.params.p)
end
```

If necessary, we can also define the support manually.
Example usage to create an Indicator leaf:

```julia
@leaf Indicator(v = 1)

# define the logpdf
function logpdf(n::Indicator{P}, x::Real) where {P<:NamedTuple{(:v)}}
    return x == n.params.v ? 0 : -Inf
end

# define support manually
support(n::Indicator{P}) where {P<:NamedTuple{(:v,)}} = n.scope => [n.params.v]
```

"""
macro leaf(expr1)
    esc(_leaf(expr1))
end

macro leaf(expr1, expr2)
    esc(_leaf(expr1, expr2))
end


text:
mutable struct Entry <: AbstractEntry
  name::AbstractString
  score::Int
end


text:
module problem116

description = """
Red, green or blue tiles
A row of five black square tiles is to have a number of its tiles replaced with coloured oblong tiles chosen from red (length two), green (length three), or blue (length four).
If red tiles are chosen there are exactly seven ways this can be done.
If green tiles are chosen there are three ways.
And if blue tiles are chosen there are two ways.
Assuming that colours cannot be mixed there are 7 + 3 + 2 = 12 ways of replacing the black tiles in a row measuring five units in length.

How many different ways can the black tiles in a row measuring fifty units in length be replaced if colours cannot be mixed and at least one coloured tile must be used?
"""

# Fill function, based on problems 114, 115
function F(m::Integer, n :: Integer, memo = zeros(Int, n))
  # Base case: no space remaining in the row
  if n < m return 0 end

  # Recursive case, count 1 for each placement of block at positions 1:n-m+1, and all ways of
  # filling the remaining space at the end of the row after each placement, memoizing results.
  if memo[n] == 0
    total = 0
    for blockpos = 1:(n-m+1)
      total += 1 + F(m, n - blockpos - m + 1, memo)
    end
    memo[n] = total
  end
  return memo[n]
end

# total ways to fill a row with only red,blue or green tiles.
function homogenous_replacements(n :: Integer)
  return F(2, n) + F(3, n) + F(4, n)
end

using Base.Test

@test F(2, 5) == 7
@test F(3, 5) == 3
@test F(4, 5) == 2
@test homogenous_replacements(5) == (7+3+2)

end

max_stars_repo_name: mbuhot/mbuhot-euler-solutions

max_issues_repo_name: mbuhot/mbuhot-euler-solutions

max_forks_repo_name: mbuhot/mbuhot-euler-solutions


text:
using DrWatson
@quickactivate "UltrasoundDesignGallery"
include("../src/UltrasoundDesignGallery.jl")

import FileIO
import Colors
import CUDA

using TestImages
using ImageView
using ImageCore

f   = ImageCore.scaleminmax(0.0, 1.0)

img = FileIO.load(DrWatson.datadir("image", "forearm.png"))
#img = FileIO.load(DrWatson.datadir("phantom", "cyst_phantom.png"))
img = Float32.(Colors.Gray.(img))

M = size(img, 1)
N = size(img, 2)

M_pad = 2^4 * ceil(Int, M / 2^4)
N_pad = 2^4 * ceil(Int, N / 2^4)

img_pad = zeros(Float32, M_pad, N_pad)
img_pad[1:M, 1:N] = img
img = img_pad

factor = 2
device = CUDA.CuDevice(0)
#device = :cpu

pyr    = Images.gaussian_pyramid(img, 4, factor, Float32(16.0))
pyr    = laplacian_pyramid(pyr, factor)

niter  = 8
Δt     = 0.7
ρ      = 0.001
σ      = 0.1
α      = 0.8
β      = 1.5
x      = srad(CUDA.CuArray(pyr[1]), Δt, ρ, niter; device=device)
#pyr[1]  = diffusion(pyr[1], 0.1, 2.0, niter)
pyr[1] = x#homo_contrast(x, σ, α, β; device=device)

x      = srad(CUDA.CuArray(pyr[2]), Δt, ρ, niter; device=device)
#pyr[2] = diffusion(pyr[2], 0.1, 1.5, niter)
pyr[2] = x#homo_contrast(x, σ, α, β; device=device)

x      = srad(CUDA.CuArray(pyr[3]), Δt, ρ, niter; device=device)
#pyr[3] = diffusion(pyr[3], 0.1, 1.5, niter)
#pyr[3] = homo_contrast(x, σ, α, β; device=device)

x      = srad(CUDA.CuArray(pyr[4]), Δt, ρ, niter; device=device)
#pyr[4] = diffusion(pyr[4], 0.1, 1.5, niter)
pyr[4] = x#homo_contrast(x, σ, α, β; device=device)

img  = synthesize_pyramid(pyr, factor)

img = f.(img)
imshow(img)
FileIO.save("pmad.png", img)


text:
@testset "Basic 2D Orientation Tests" begin
	o = (0f0,0f0); x = (1,0); y = (0,1); p = (2,2); q = (3,4); r = (1f0,1f0);
	@test ccw(o,x,p) == ccw(o,x,toPoint2H(p))
	@test ccw(o,x,p) == ccw(o,toPoint2H(x),(p))
	@test ccw(o,x,p) == ccw(toPoint2H(o),x,(p))
	@test ccw(o,x,p) == ccw(x,p,toPoint2H(o))
	@test ccw(o,x,p) == ccw(toPoint2H(o),toPoint2H(x),toPoint2H(p))
	@test ccw(o,x,p) == ccw(toPoint2H(o),toPoint2H(x),(p))
	@test ccw(o,x,p) == ccw(toPoint2H(o),(x),toPoint2H(p))
	@test ccw(o,x,p) == ccw((o),toPoint2H(x),toPoint2H(p))
	@test 0 == ccw(o,r,toPoint2H(p))
	@test 0 == ccw(r,p,toPoint2H(o))
	@test 0 == ccw((o),toPoint2H(r),toPoint2H(p))
	@test 0 == ccw(toPoint2H(o),toPoint2H(r),(p))
	@test 0 == ccw(toPoint2H(o),toPoint2H(r),toPoint2H(p))

	@test leftturn(o,x,p)
	@test !rightAhead(o,x,p)
	@test rightturn(o,y,q)
	@test !leftAhead(o,y,q)
	@test !leftturn(x,p,q)
	@test leftAhead(x,p,q) # colinear
	@test rightAhead(x,p,q) # colinear
	@test !leftAhead(o,p,r) # colinear
	@test leftOnAhead(y,o,q)
	@test leftOnAhead(o,p,r) # colinear
	@test rightOnAhead(o,p,r) # colinear
	@test rightOnAhead(x,p,q) # colinear
	@test !rightOnAhead(o,x,p) # colinear
	@test leftofseg(p,y,x)  # query point with segment
	@test rightofseg(p,x,y)
	@test leftward(p,r, o,y)
	@test !rightward(p,r, o,y)
	@test !leftward(p,r, o,r) # parallel
	@test !rightward(p,r, o,r) # parallel
	@test line(o,y) == (0,-1,0)
	@test line(y,r) == (-1,0,1)
	@test leftturn(line(o,x), line(x,p), line(p,o))
	@test rightof(p,line(r,q))
	@test leftof(r,line(p,q))
	@test !leftof(o, line(p,r))
	@test !rightof(o, line(p,r))
	@test rightof(toPoint2H(p),line(r,q))
	@test leftof(toPoint2H(r),line(p,q))
	@test !leftof(toPoint2H(o), line(p,r))
	@test !rightof(toPoint2H(o), line(p,r))
end;

max_stars_repo_name: snoeyink/CGinJulia

max_issues_repo_name: snoeyink/CGinJulia

max_forks_repo_name: snoeyink/CGinJulia


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct ProbabilisticInternal <: ForecastInternal
        label::String
        resolution::Dates.Period
        initial_time::Dates.DateTime
        percentiles::Vector{Float64}
        time_series_uuid::UUIDs.UUID
        horizon::Int
        internal::InfrastructureSystemsInternal
    end

A Probabilistic forecast for a particular data field in a Component.

# Arguments
- `label::String`: label of component parameter forecasted
- `resolution::Dates.Period`
- `initial_time::Dates.DateTime`: forecast availability time
- `percentiles::Vector{Float64}`: Percentiles for the probabilistic forecast
- `time_series_uuid::UUIDs.UUID`: reference to time series data; timestamp - scalingfactor
- `horizon::Int`: length of this forecast
- `internal::InfrastructureSystemsInternal`
"""
mutable struct ProbabilisticInternal <: ForecastInternal
    "label of component parameter forecasted"
    label::String
    resolution::Dates.Period
    "forecast availability time"
    initial_time::Dates.DateTime
    "Percentiles for the probabilistic forecast"
    percentiles::Vector{Float64}
    "reference to time series data; timestamp - scalingfactor"
    time_series_uuid::UUIDs.UUID
    "length of this forecast"
    horizon::Int
    internal::InfrastructureSystemsInternal
end

function ProbabilisticInternal(label, resolution, initial_time, percentiles, time_series_uuid, horizon, )
    ProbabilisticInternal(label, resolution, initial_time, percentiles, time_series_uuid, horizon, InfrastructureSystemsInternal(), )
end

function ProbabilisticInternal(; label, resolution, initial_time, percentiles, time_series_uuid, horizon, )
    ProbabilisticInternal(label, resolution, initial_time, percentiles, time_series_uuid, horizon, )
end

"""Get ProbabilisticInternal label."""
get_label(value::ProbabilisticInternal) = value.label
"""Get ProbabilisticInternal resolution."""
get_resolution(value::ProbabilisticInternal) = value.resolution
"""Get ProbabilisticInternal initial_time."""
get_initial_time(value::ProbabilisticInternal) = value.initial_time
"""Get ProbabilisticInternal percentiles."""
get_percentiles(value::ProbabilisticInternal) = value.percentiles
"""Get ProbabilisticInternal time_series_uuid."""
get_time_series_uuid(value::ProbabilisticInternal) = value.time_series_uuid
"""Get ProbabilisticInternal horizon."""
get_horizon(value::ProbabilisticInternal) = value.horizon
"""Get ProbabilisticInternal internal."""
get_internal(value::ProbabilisticInternal) = value.internal


text:
const compilecache = Dict{UInt,Any}()
const operation_annotation_color = :light_black

"""
    Routine

Abstract type for general YaoLang programs.
"""
abstract type Routine end

"""
    Operation

Abstract type operations on registers. `Operation`s
are objects that can be applied to a register. This
can be a non-parameterized gate, or parameterized gate
with known parameters, or pulses.
"""
abstract type Operation end

struct GenericRoutine{name} <: Routine end

abstract type IntrinsicRoutine <: Routine end

Base.parent(x::Routine) = x

function print_routine(io::IO, x::IntrinsicRoutine)
    print(io, routine_name(x))
    fnames = fieldnames(typeof(x))
    if !isempty(fnames)
        print(io, "(")
        for k in 1:length(fnames)
            print(io, getfield(x, fnames[k]))
            if k != length(fnames)
                print(io, ", ")
            end
        end
        print(io, ")")
    end
end

function Base.show(io::IO, x::IntrinsicRoutine)
    print_routine(io, x)
end

function Base.show(io::IO, ::MIME"text/plain", x::IntrinsicRoutine)
    print_routine(io, x)
    printstyled(io, " (intrinsic operation)"; color = operation_annotation_color)
end

function Base.show(io::IO, fn::GenericRoutine{name}) where {name}
    print(io, name)
end

function Base.show(io::IO, ::MIME"text/plain", fn::GenericRoutine{name}) where {name}
    print(io, name)
    printstyled(
        io,
        " (generic routine with ",
        length(methods(fn).ms),
        " methods)";
        color = operation_annotation_color,
    )
end

# NOTE: kwargs is not supported for now
struct RoutineSpec{P,Vars} <: Operation
    parent::P
    variables::Vars

    function RoutineSpec(parent, vars...)
        new{typeof(parent),typeof(vars)}(parent, vars)
    end
end

Base.parent(x::RoutineSpec) = x.parent
Base.:(==)(::GenericRoutine{A}, ::GenericRoutine{A}) where {A} = true

# NOTE: the reason we don't use a gensym
# here for stub function is
"stub for routine CodeInfo"
struct RoutineStub end

const routine_stub = RoutineStub()

struct Adjoint{P<:Operation} <: Operation
    parent::P
end

Base.adjoint(x::Operation) = Adjoint(x)
Base.adjoint(x::Adjoint) = x.parent

routine_name(::Type) = nothing
routine_name(x) = routine_name(typeof(x))
routine_name(::Type{<:GenericRoutine{name}}) where {name} = name
routine_name(::Type{T}) where {T<:IntrinsicRoutine} = nameof(T)
routine_name(::Type{<:RoutineSpec{P}}) where {P} = routine_name(P)
routine_name(::Type{<:Adjoint{P}}) where {P} = Symbol(routine_name(P), "_dag")

struct DeviceError <: Exception
    msg::String
end

"""
Semantic extension for Julia.
"""
module Semantic
using ExprTools
using ..YaoCompiler: DeviceError, Operation, Locations, CtrlLocations, IntrinsicRoutine

# semantic stubs
const SEMANTIC_STUBS = Symbol[]
const PRESERVED_MACROS = Symbol[]

macro semantic_stub(ex::Expr)
    ex.head === :call || throw(ParseError("expect function call"))
    ex.args[1] isa Symbol || throw(ParseError("stub must be a function"))

    def = Dict{Symbol,Any}()
    name = ex.args[1]
    macroname = Symbol("@$name")
    def[:name] = name
    def[:args] = ex.args[2:end]
    err = DeviceError("@$name is not executed as a quantum program")
    def[:body] = quote
        throw($err)
        # nothing
    end

    quote
        push!(SEMANTIC_STUBS, $(QuoteNode(name)))
        push!(PRESERVED_MACROS, $(QuoteNode(macroname)))
        @noinline $(combinedef(def))
    end |> esc
end

@semantic_stub main(gate::Union{Operation,IntrinsicRoutine})
@semantic_stub ctrl(gate::Union{Operation,IntrinsicRoutine}, loc::Locations, ctrl::CtrlLocations)
@semantic_stub gate(gate::Union{Operation,IntrinsicRoutine}, loc::Locations)
# NOTE: other measurement options are just syntax sugars
# TODO: check in typeinf
@semantic_stub measure(locs::Locations)
@semantic_stub barrier(locs::Locations)

end # Semantic

function gate_sugar(ex::Expr)
    is_gate_location(ex) || throw(ParseError("invalid syntax for gate: $ex"))
    return ex.args[3], ex.args[2]
end

function gate_m(ex::Expr)
    gate, locs = gate_sugar(ex)
    return :($(GlobalRef(Semantic, :gate))($gate, $Locations($locs)))
end

function ctrl_m(ctrl, ex::Expr)
    gate, locs = gate_sugar(ex)
    return :($(GlobalRef(Semantic, :ctrl))($gate, $Locations($locs), $CtrlLocations($ctrl)))
end

macro ctrl(ctrl, ex::Expr)
    return esc(ctrl_m(ctrl, ex))
end

macro gate(ex::Expr)
    return esc(gate_m(ex))
end

macro measure(locs)
    return Expr(:call, GlobalRef(Semantic, :measure), :($Locations($locs)))
end

macro barrier(locs)
    return esc(Expr(:call, GlobalRef(Semantic, :barrier), locs))
end

macro device(ex::Expr)
    def = splitdef(ex; throw = false)

    if isnothing(def)
        ex.head === :call || throw(ParseError("invalid syntax of @device: $ex"))
        return esc(device_call(ex))
    end
    return esc(device_def(def))
end

macro device(ex...)
    call = ex[end]
    kwargs = ex[1:end-1]
    return esc(device_call(call, kwargs))
end

# frontend, only handles some syntax sugars
function device_def(def::Dict)
    args = get(def, :args, Any[])
    stub_def = copy(def)
    device_def = copy(def)
    code = Expr(:block)
    isfunction = !(haskey(def, :name) && !(def[:name] isa Symbol))

    if haskey(def, :name) && !(def[:name] isa Symbol)
        self = rm_annotations(def[:name])
        self_annotation = annotations(def[:name])
        name = nothing
    else
        self = gensym(:self)
        name = get(def, :name, gensym(:device))
        self_annotation = GenericRoutine{name}
    end

    device_def[:name] = :($self::$self_annotation)
    device_def[:body] = quote
        $RoutineSpec($self, $(rm_annotations.(args)...))
    end

    stub_def[:name] = :(::$RoutineStub)
    stub_def[:args] = [:($self::$self_annotation), args...]
    stub_def[:body] = preprocess_device_gate_syntax(def[:body])

    push!(code.args, combinedef(device_def))
    push!(code.args, combinedef(stub_def))

    if !isnothing(name)
        push!(code.args, :(const $name = $self_annotation()))
    end

    push!(code.args, name)
    return code
end

function preprocess_device_gate_syntax(ex)
    ex isa Expr || return ex
    if is_gate_location(ex)
        location = ex.args[2]
        gate = ex.args[3]
        return Expr(:call, GlobalRef(Semantic, :gate), gate, :($Locations($location)))
    end

    #= make @inbounds etc. work =#
    # exclude YaoCompiler macros
    if ex.head in [:block, :if, :for] || (ex.head === :macrocall && !is_preserved_macro(ex))
        return Expr(ex.head, map(preprocess_device_gate_syntax, ex.args)...)
    end
    return ex
end

function device_call(ex::Expr, kwargs) end

function is_gate_location(ex)
    ex isa Expr || return false
    return (ex.head === :call) && (ex.args[1] == :(=>))
end

function is_preserved_macro(ex::Expr)
    ex.head === :macrocall || return false
    return ex.args[1] in Semantic.PRESERVED_MACROS
end

@inline function (spec::RoutineSpec)(r::AbstractRegister)
    return execute(Semantic.main, r, spec)
end

@inline function (spec::RoutineSpec)(r::AbstractRegister, locs::Locations)
    return execute(Semantic.gate, r, spec, locs)
end

@inline function (spec::RoutineSpec)(r::AbstractRegister, locs::Locations, ctrl::CtrlLocations)
    return execute(Semantic.ctrl, r, spec, locs, ctrl)
end


text:
# NOTE: run with $ julia --depwarn=yes for these tests to pass
seed!(RNG, 1337)

@testset "histogram" begin
  x = randn(RNG, 10000)
  @test_logs (:warn, r"`symb`.+deprecated") histogram(x, symb = "#", closed = :right)
  @test_logs (:warn, r"`bins`.+deprecated") histogram(x, bins = 5, closed = :right)
  @test_logs (:warn, r"nbins::Int.+deprecated") histogram(x, 5, closed = :right)
end

@testset "barplot" begin
    @test_logs (:warn, r"`symb`.+deprecated") barplot(["Paris", "Moskau"], [2.244, 11.92]; symb = "#")
end

@testset "plot" begin
    p = @inferred Plot([1, 2], [-1, 1])
    @test_logs (:warn, r"`annotate!`.+renamed") annotate!(p, :l, "text 1")
    @test_logs (:warn, r"`annotate!`.+renamed") annotate!(p, :l, "text 2", :red)
    @test_logs (:warn, r"`annotate!`.+renamed") annotate!(p, :l, "text 3", color=:blue)

    for i in 1:2
      @test_logs (:warn, r"`annotate!`.+renamed") annotate!(p, :l, i, "$i", color=:green)
      @test_logs (:warn, r"`annotate!`.+renamed") annotate!(p, :l, i, "$i", :yellow)
    end

    @test_logs (:warn, r"`UnicodePlots.savefig`.+text") savefig(p, tempname() * ".png")
end

@testset "spy" begin
    @test_logs (:warn, r"`color = :automatic`.+deprecated") spy(_stable_sprand(RNG, 10, 10, .15), color=:automatic)
end


text:
######## FARIMA Process #########
"""
Fractional Auto-regression Integrated Moving Average (FARIMA) process.
"""
struct FARIMA <: DiscreteTimeStochasticProcess
    d::Real
    ar::Vector{AbstractFloat}
    ma::Vector{AbstractFloat}

    FARIMA(d::Real, ar::Vector{<:AbstractFloat}, ma::Vector{<:AbstractFloat}) = new(d, ar, ma)

    # function FARIMA(d::Real, ar::Vector{AbstractFloat}, ma::Vector{AbstractFloat})
    #     # (typeof(P)==Int64 && P==length(ar)) || error("Inconsistent parameters for auto-regression.")
    #     # (typeof(Q)==Int64 && Q==length(ma)) || error("Inconsistent parameters for moving average.")
    #     # abs(d) < 0.5 || error("Order of fractional differential must be in the range (-0.5, 0.5).")
    #     # new(d, all(ar.==0) ? Float64[] : ar, all(ma.==0) ? AbstractFloat[] : ma)
    #     new(d, ar, ma)
    # end
end

FARIMA(d::Real) = FARIMA(d, Float64[], Float64[])

# parameters(X::FARIMA{P,Q}) where {P,Q} = Dict('p'=>P, 'd'=>X.d, 'q'=>Q)
# parameters(X::FARIMA{P,Q}) where {P,Q} = (P, X.d, Q)

# FARIMA(d::Real, ar::Vector{AbstractFloat}, ma::Vector{AbstractFloat}) = FARIMA{length(ar), length(ma)}(d, ar, ma)

######## Fractional Integrated Process ########
"""
Fractional integrated process.

This is a stationary process defined by ∇^d X(t) = ε(t), where ∇^d, d in (-1/2, 1/2) is the fractional differential operator, and ε(t) are i.i.d. standard Gaussian variables.
"""
struct FractionalIntegrated <: DiscreteTimeStationaryProcess
    d::Real

    function FractionalIntegrated(d::Real)
        abs(d) < 0.5 || error("Order of fractional differential must be in the range (-0.5, 0.5).")
        # new(d, all(ar.==0) ? Float64[] : ar, all(ma.==0) ? Float64[] : ma)
        new(d)
    end
end

ss_exponent(X::FractionalIntegrated) = X.d + 1/2

partcorr(X::FractionalIntegrated, k::DiscreteTime) = X.d/(k-X.d)

function autocov(X::FractionalIntegrated, n::DiscreteTime)
    return n > 0 ? (n-1+X.d) / (n-X.d) * autocov(X, n-1) : gamma(1-2*X.d) / gamma(1-X.d)^2
end

"""
Note: The covariance of a fractional integrated process is computed recursively, and `autocov!()` is reloaded for reason of efficiency.
"""
function autocov!(C::Vector{<:AbstractFloat}, X::FractionalIntegrated, G::AbstractVector{<:DiscreteTime})
    # check dimension
    # @assert length(C) == length(G)
    # @assert isregulargrid(G) && eltype(G)<:Int

    V = zeros(G[end])
    V[1] = gamma(1-2*X.d) / gamma(1-X.d)^2  # cov(0)
    for n = 1:length(V)-1
        V[n+1] = V[n] * (n-1+X.d) / (n-X.d)
    end
    for n = 1:length(C)
        C[n] = V[G[n]]
    end
    return C
end

max_stars_repo_name: yanncalec/FracFin

max_issues_repo_name: yanncalec/FracFin

max_forks_repo_name: yanncalec/FracFin


text:
function eigen(x::Matrix{WElem})
    @info "use wolframscript"
    u, v = WS"Eigensystem"(x)
    Eigen(u, transpose(v))
end

function svd(x::Matrix{WElem})
    @info "use wolframscript"
    u, v, w = WS"SingularValueDecomposition"(x)
    SVD(u, diag(v), w')
end

function qr(x::Matrix{WElem})
    @info "use wolframscript"
    q, r = WS"QRDecomposition"(x)
    (Q = q', R = r)
end

function lq(x::Matrix{WElem})
    @info "use wolframscript"
    q, r = WS"QRDecomposition"(permutedims(x, (2,1)))
    (L = transpose(r), Q = conj(q))
end

function schur(x::Matrix{WElem})
    @info "use wolframscript"
    z, t = WS"SchurDecomposition"(x)
    Schur(t, z, diag(t))
end

inv(x::Matrix{WElem}) = WS"Inverse"(x, true)
pinv(x::Matrix{WElem}) = WS"PseudoInverse"(x, true)
nullspace(x::Matrix{WElem}) = permutedims(WS"NullSpace"(x, true), (2,1))


text:
#using Revise

# Add PATH
push!(LOAD_PATH, joinpath(@__DIR__, "..", "src"))

#using PyPlot
using Plots
using BenchmarkTools
using DifferentialEquations

using LSM.PhysCon
using LSM.WaterMod
#using LSM.LeafPhotosynthesisMod
using LSM.LeafEnergyWaterBalanceMod

output_dir = joinpath(@__DIR__,"..","output")
mkpath(output_dir)

# T=270:0.1:315
# muT = zeros(size(T));

# for i=1:length(T)
#     muT[i] = μ_l(T[i]);
# end
# plot(T-273.15*ones(size(T)),muT/μ_l(273.15+20),label="mu/mu(20C)")

# function lorenz!(du,u,p,t)
#  du[1] = 10.0*(u[2]-u[1]);
#  du[2] = u[1]*(28.0-u[3]) - u[2];
#  du[3] = u[1]*u[2] - (8/3)*u[3];
# end
# u0 = [1.0;0.0;0.0];
# tspan = (0.0,100.0);
# prob = ODEProblem(lorenz!,u0,tspan);
# sol = solve(prob);

# plot(sol,vars=(1,2,3))

# Create a leaf structure
l = leaf_params{Float32}();
# Create a Flux structure
f = fluxes{Float32}();
# Create a meteo structure
met = meteo{Float32}();

#met = meteo();

# initialize some reasonable values
f.Je   = 100; f.gbc  = 100; f.gbv  = 100; f.ceair= 1500; f.eair = 1500; f.APAR = 500; f.H=0;f.LE=0; # leaf should not have eair
l.Kn = 2.44; l.α=0.2; l.ε=0.98; l.LMA=100e-3; l.RWC=80/100;l.psi_l=-1e6;l.psi_l50 = -2e6;l.ck=3;met.zscreen = 2.0;
l.height   = 1.0; met.zscreen  = 2.0;
met.stab_type_stable = 2;

# A diurnal cycle for radiation and Tair
Deltat  = 60.0;
Samp    = 500.0; # W/m2 amplitude
DeltaSamp = 100.0;
Tmean   = 273.15+20;
DeltaT  = 3.0;
t       = range(0, stop=24*3600, step=Deltat); # diurnal cycle in seconds
#print(t)
Sdown_t = Samp*ones(size(t));
Tair_t  = Tmean*ones(size(t));
zeros_t = zeros(size(t));
for i = 1:length(Sdown_t)
    if(t[i]>3*3600)
         Sdown_t[i] = Samp+DeltaSamp;
    end
end
#print(typeof(phase[1]));
#print(typeof(Sdown_t[1]));
#print(Sdown_t);
#clf();
#fig = figure("Diurnal cycle", figsize=(10,5));
# plot(t/3600,Sdown_t)
#plot!(t/3600,10*(Tair_t-273.15*ones(size(Tair_t))))
#    title = "Shortwave Incoming radiation (W/m^2)");#,
#    xlabel = "Hours",
#    ylabel = "S_{down} (W/m^2)")
#plot(t/3600,Tair_t-273.15*ones(size(Tair_t)),
#    title = "Temperature (C)",
#    xlabel = "Hours",
#    ylabel = "Temperature (C)")

plot(t/3600,Sdown_t,label="Sdown (W/m2)")

plot(t/3600,Tair_t-273.15*ones(size(Tair_t)),label="Tair (C)")

psi_s      = -0.3e6 ; # soil water potential (Pa)
U          =  1.0;
RH         =  65/100;
eps_air    =  0.75;
tspan      =  (0.0,Deltat);
N          =  length(Sdown_t);
mutable struct parameters_ode
    l::leaf_params;
    met::meteo;
    f::fluxes;
    psi_s;
end

# small time stepping
dt     = 0.1*60; # in s
T_t    = zeros(size(Sdown_t));
psil_t = zeros(size(Sdown_t));
Cc_t   = zeros(size(Sdown_t));
Rn_t   = zeros(size(Sdown_t));
GPP_t  = zeros(size(Sdown_t));
GPP_diffusion_t  = zeros(size(Sdown_t));
LUE_t  = zeros(size(Sdown_t));
H_t    = zeros(size(Sdown_t));
LE_t   = zeros(size(Sdown_t));
rs_t   = zeros(size(Sdown_t));
ra_t   = zeros(size(Sdown_t));

function f_ode!(du,u,p,t) # p are parameters
    du .= LeafEnergyWaterBalance(u[1], u[2], u[3], p.met, p.l, p.f, p.psi_s);
    #println("du_inside = $(du), u_inside = $(u)")
end

let
    # initial conditions
    met.T_air  = Tair_t[1];
    l.T        = met.T_air;
    l.psi_l    = psi_s;
    l.Cc       = 0.6*met.Ca;
    #println("Tair1=",met.T_air," Tleaf1=",l.T," psi_leaf1=",l.psi_l)

    for i=1:N
        met.S_down = Sdown_t[i];
        met.L_down = eps_air*physcon.σ*(Tair_t[i])^4;
        met.T_air  = Tair_t[i];
        met.e_air  = RH*SatVap(Tair_t[i])[1];
        met.PAR    = 45/100*physcon.Wtoμmole_s*Sdown_t[i];
        f.APAR     = met.PAR;
        met.U      = U;
        met.Ca     = 400.0; #ppm
        f.APAR     = met.PAR;
        #println("Tair=",met.T_air," Tleaf=",l.T," psi_leaf=",l.psi_l)
        for j=1:trunc(Deltat/dt)
            u    = [l.T;l.psi_l;l.Cc];
            p    = parameters_ode(l,met,f,psi_s);
            #(p.met, p.l,  p.psi_s, p.U)    = [l;met;psi_s;U];
            #prob = ODEProblem(f_ode!,u0,tspan,p);
            #du   = zeros(size(u));
            #f_ode!(du,u,p,t,dummy);
            #println("du_outside = $(du), u_outside = $(u)")
            # Rn_t[i] = dummy[3]; H_t[i] = dummy[4]; LE_t[i] = dummy[5];

#             du   = zeros(size(u));
#             f_ode!(du,u,p,t);
#             (l.T,l.psi_l,l.Cc) = du*dt+u;

            u0   = [l.T;l.psi_l;l.Cc];
            prob = ODEProblem(f_ode!,u0,tspan,p);
            sol  = solve(prob);
            # save values
            met = p.met;
            f   = p.f;
            l   = p.l;
            (l.T,l.psi_l,l.Cc) = sol[1:3,end];

            #println("Cc_out=",l.Cc)
            Rn_t[i] = p.f.Rn; H_t[i] = p.f.H; LE_t[i] = p.f.LE;
            rs_t[i] = 1.0/(p.l.gs/p.f.g_m_s_to_mol_m2_s);
            ra_t[i] = p.f.ra;
            GPP_t[i]= p.f.An_biochemistry;
            GPP_diffusion_t[i]= p.f.An_diffusion;
            LUE_t[i]= p.f.An_biochemistry/f.APAR;

            #T_old   = l.T;
            #(l.T,l.psi_l) = du*dt+u;
            T_t[i]    = l.T ; #  = T_old;
            psil_t[i] = l.psi_l;
            Cc_t[i]   = l.Cc;

            if(abs(H_t[i])>500)
                println("index ($i) ($j)")
            end

#             (dumb,l.psi_l) = du*dt+u;
#             T_t[i]  = dumb ; #  = T_old;

        end
    end

#     l.T        = Tair_t[1]; # initialize temperature of the leaf
#     l.psi_l    = psi_s;

#     for i=1:N
#         met.S_down = Sdown_t[i];
#         met.L_down = eps_air*physcon.σ*(Tair_t[i])^4;
#         met.T_air  = Tair_t[i];
#         met.ea_air = RH*SatVap(Tair_t[i])[1];
#         u0   = [l.T;l.psi_l];
#         p    = parameters_ode(l,met,psi_s);
#         #(p.met, p.l,  p.psi_s, p.U)    = [l;met;psi_s;U];
#         prob = ODEProblem(f_ode!,u0,tspan,p);
#         sol  = solve(prob);
#         (l.T,l.psi_l) = sol[:,end];
#         T_t[i] = l.T; psil_t[i]=l.psi_l;
#         #print((l.T,l.psi_l) )
#         #dT_dt,dH2Ol_dt = LeafEnergyWaterBalance(met, l, psi_s);
#         #l.T        = l.T + Deltat*dT_dt;
#     end

end

plot(t/3600, T_t-273.15*ones(size(T_t)),xlabel = "t (hours)",ylabel = "T (C)",label="T (C)",ylim=0:100)
#plot!(t/3600,-psil_t/100000,xlabel = "t (s)",ylabel = "|psi_l (1000 hPa)|",label="psi_l",ylim=0:100)

plot(t-3*3600*ones(size(t)), T_t-273.15*ones(size(T_t)),xlabel = "t (s)",ylabel = "T (C)",label="T (C)",xlims=(0*3600,0.1*3600))

savefig(joinpath(output_dir, "T_step_t.png"))

plot(t-3*3600*ones(size(t)), Cc_t,xlabel = "t (s)",ylabel = "Cc (ppm)",label="Cc",xlims=(0*3600,0.1*3600))

plot(t-3*3600*ones(size(t)),  Rn_t,xlabel = "t (s)",ylabel = "Rn (W/m2)",label="Rn",xlims=(0*3600,0.1*3600))
plot!(t-3*3600*ones(size(t)), H_t,xlabel  = "t (s)",ylabel = "H (W/m2)" ,label="H",xlims=(0*3600,0.1*3600))
plot!(t-3*3600*ones(size(t)), LE_t,xlabel = "t (s)",ylabel = "LE (W/m2)",label="LE",xlims=(0*3600,0.1*3600))

savefig(joinpath(output_dir, "Fluxes_step_t.png"))

plot(t-3*3600*ones(size(t)),  rs_t,xlabel = "t (s)",ylabel = "rs (s/m)",label="rs",xlims=(0*3600,0.1*3600))
plot!(t-3*3600*ones(size(t)), ra_t,xlabel  = "t (s)",ylabel = "ra (s/m)" ,label="ra",xlims=(0*3600,0.1*3600))

savefig(joinpath(output_dir, "resistances_step_t.png"))

plot(t/3600, GPP_t,xlabel = "t (hr)",ylabel = "GPP (micomles/m2/s)",label="GPP",ylim=0:100)
plot!(t/3600, GPP_diffusion_t,xlabel = "t (hr)",ylabel = "GPP  (micomles/m2/s)",label="GPP diffusion",ylim=0:100)

plot(t-3*3600*ones(size(t)), GPP_t-GPP_diffusion_t,xlabel = "t (s)",ylabel = "GPP - GPP diffusion (micomles/m2/s)",label="GPP biochemistry",xlims=(0*3600,0.5*3600),ylim=0:100)

savefig(joinpath(output_dir, "GPP_step_t.png"))

plot(t/3600, LUE_t,xlabel = "t (hr)",ylabel = "LUE (-)",label="LUE",ylim=0:100)

savefig(joinpath(output_dir, "LUE_step_t.png"))

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
module DarcyCubedSphereTestsMPI
  using PartitionedArrays
  using Test
  using FillArrays
  using Gridap
  using GridapPETSc
  using GridapGeosciences

  include("../DarcyCubedSphereTests.jl")
  include("../ConvergenceAnalysisTools.jl")

  function petsc_gamg_options()
    """
      -ksp_type gmres -ksp_rtol 1.0e-06 -ksp_atol 0.0
      -ksp_monitor -pc_type gamg -pc_gamg_type agg
      -mg_levels_esteig_ksp_type gmres -mg_coarse_sub_pc_type lu
      -mg_coarse_sub_pc_factor_mat_ordering_type nd -pc_gamg_process_eq_limit 50
      -pc_gamg_square_graph 9 pc_gamg_agg_nsmooths 1
    """
  end
  function petsc_mumps_options()
    """
      -ksp_type preonly -ksp_error_if_not_converged true
      -pc_type lu -pc_factor_mat_solver_type mumps
    """
  end
  function main(parts)
    GridapPETSc.with(args=split(petsc_gamg_options())) do
       num_refs=[2,3,4,5]
       hs=[2.0/2^n for n in num_refs]
       model_args_series=zip(Fill(parts,length(num_refs)),num_refs)
       t = PartitionedArrays.PTimer(parts,verbose=true)
       PartitionedArrays.tic!(t)
       hs1,k1errors,s1=convergence_study(solve_darcy,hs,model_args_series,1,4,PETScLinearSolver())
       PartitionedArrays.toc!(t,"DarcyCubedSphereConvergenceStudy")
       display(t)
       println(hs1)
       println(k1errors)
       println(s1)
    end
  end
  prun(main,mpi,4)
end #module

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
@with_kw struct DiffusionModel
    drift::Array{Float64,1}
    sigma::Array{Float64,1}
    bound_hi::Array{Float64, 1}
    bound_lo::Array{Float64, 1}
    ndt::NTuple{2, Float64}
    Δt::Float64 = 0.01
    function DiffusionModel(drift::Array{Float64,1},
                            sigma::Array{Float64,1},
                            bound_hi::Array{Float64,1},
                            bound_lo::Array{Float64},
                            ndt::NTuple{2, Float64},
                            Δt::Float64)
            new(drift, sigma, bound_hi, bound_lo, ndt, Δt)
    end
end

function DiffusionModel(drift::Real, sigma::Real,
                bound_hi::Real, bound_lo::Real,
                ndt::NonDecisionTime;
                Δt::Real = 0.01)
    @assert bound_hi > 0
    @assert bound_lo < 0

    drift = convert(Array{Float64, 1}, [drift])
    sigma = convert(Array{Float64, 1}, [sigma])
    bound_hi = convert(Array{Float64, 1}, [bound_hi])
    bound_lo = convert(Array{Float64, 1}, [bound_lo])
    ndt = convert(NTuple{2, Float64}, (ndt.lower, ndt.upper))
    Δt = convert(Float64, Δt)

    out = DiffusionModel(drift, sigma,
            bound_hi, bound_lo, ndt, Δt)
    return out
end

function DiffusionModel(drift::AbstractDrift, sigma::AbstractSigma,
                bounds::AbstractBounds, ndt::NonDecisionTime;
                Δt::Real = 0.01)
    hi, lo = get(bounds)
    drift = convert(Array{Float64, 1}, get(drift))
    sigma = convert(Array{Float64, 1}, get(sigma))
    bound_hi = convert(Array{Float64, 1}, hi)
    bound_lo = convert(Array{Float64, 1}, lo)
    ndt = convert(NTuple{2, Float64}, (ndt.lower, ndt.upper))
    Δt = convert(Float64, Δt)

    out = DiffusionModel(drift, sigma,
            bound_hi, bound_lo, ndt, Δt)
    return out
end


text:
module IntervalSets

using Base: @pure
import Base: eltype, convert, show, in, length, isempty, isequal, issubset, ==, hash,
             union, intersect, minimum, maximum, extrema, range, clamp, float, ⊇, ⊊, ⊋

using Statistics
import Statistics: mean

using Dates

export AbstractInterval, Interval, OpenInterval, ClosedInterval,
            ⊇, .., ±, ordered, width, duration, leftendpoint, rightendpoint, endpoints,
            isopenset, isclosedset, isleftclosed, isrightclosed,
            isleftopen, isrightopen, closedendpoints,
            infimum, supremum

"""
A subtype of `Domain{T}` represents a subset of type `T`, that provides `in`.
"""
abstract type Domain{T} end

Base.IteratorSize(::Type{<:Domain}) = Base.SizeUnknown()

"""
A subtype of `AbstractInterval{T}` represents an interval subset of type `T`, that provides
`endpoints`, `closedendpoints`.
"""
abstract type AbstractInterval{T} <: Domain{T} end

"A tuple containing the left and right endpoints of the interval."
endpoints(d::AI) where AI<:AbstractInterval = error("Override endpoints(::$(AI))")

"The left endpoint of the interval."
leftendpoint(d::AbstractInterval) = endpoints(d)[1]

"The right endpoint of the interval."
rightendpoint(d::AbstractInterval) = endpoints(d)[2]

"A tuple of `Bool`'s encoding whether the left/right endpoints are closed."
closedendpoints(d::AI) where AI<:AbstractInterval = error("Override closedendpoints(::$(AI))")

"Is the interval closed at the left endpoint?"
isleftclosed(d::AbstractInterval) = closedendpoints(d)[1]

"Is the interval closed at the right endpoint?"
isrightclosed(d::AbstractInterval) = closedendpoints(d)[2]

# open_left and open_right are implemented in terms of closed_* above, so those
# are the only ones that should be implemented for specific intervals
"Is the interval open at the left endpoint?"
isleftopen(d::AbstractInterval) = !isleftclosed(d)

"Is the interval open at the right endpoint?"
isrightopen(d::AbstractInterval) = !isrightclosed(d)

# Only closed if closed at both endpoints, and similar for open
isclosedset(d::AbstractInterval) = isleftclosed(d) && isrightclosed(d)

"Is the interval open?"
isopenset(d::AbstractInterval) = isleftopen(d) && isrightopen(d)

eltype(::Type{AbstractInterval{T}}) where {T} = T
@pure eltype(::Type{I}) where {I<:AbstractInterval} = eltype(supertype(I))

convert(::Type{AbstractInterval}, i::AbstractInterval) = i
convert(::Type{AbstractInterval{T}}, i::AbstractInterval{T}) where T = i

ordered(a::T, b::T) where {T} = ifelse(a < b, (a, b), (b, a))
ordered(a, b) = ordered(promote(a, b)...)

checked_conversion(::Type{T}, a, b) where {T} = _checked_conversion(T, convert(T, a), convert(T, b))
_checked_conversion(::Type{T}, a::T, b::T) where {T} = a, b
_checked_conversion(::Type{Any}, a, b) = throw(ArgumentError("$a and $b promoted to type Any"))
_checked_conversion(::Type{T}, a, b) where {T} = throw(ArgumentError("$a and $b are not both of type $T"))

function infimum(d::AbstractInterval{T}) where T
    a = leftendpoint(d)
    b = rightendpoint(d)
    a > b && throw(ArgumentError("Infimum not defined for empty intervals"))
    a
end

function supremum(d::AbstractInterval{T}) where T
    a = leftendpoint(d)
    b = rightendpoint(d)
    a > b && throw(ArgumentError("Supremum not defined for empty intervals"))
    b
end

mean(d::AbstractInterval) = (leftendpoint(d) + rightendpoint(d))/2

"""
    w = width(iv)

Calculate the width (max-min) of interval `iv`. Note that for integers
`l` and `r`, `width(l..r) = length(l:r) - 1`.
"""
function width(A::AbstractInterval)
    _width = rightendpoint(A) - leftendpoint(A)
    max(zero(_width), _width)   # this works when T is a Date
end

"""
A subtype of `TypedEndpointsInterval{L,R,T}` where `L` and `R` are `:open` or `:closed`,
that represents an interval subset of type `T`, and provides `endpoints`.
"""
abstract type TypedEndpointsInterval{L,R,T} <: AbstractInterval{T} end

closedendpoints(d::TypedEndpointsInterval{:closed,:closed}) = (true,true)
closedendpoints(d::TypedEndpointsInterval{:closed,:open}) = (true,false)
closedendpoints(d::TypedEndpointsInterval{:open,:closed}) = (false,true)
closedendpoints(d::TypedEndpointsInterval{:open,:open}) = (false,false)

in(v, I::TypedEndpointsInterval{:closed,:closed}) = leftendpoint(I) ≤ v ≤ rightendpoint(I)
in(v, I::TypedEndpointsInterval{:open,:open}) = leftendpoint(I) < v < rightendpoint(I)
in(v, I::TypedEndpointsInterval{:closed,:open}) = leftendpoint(I) ≤ v < rightendpoint(I)
in(v, I::TypedEndpointsInterval{:open,:closed}) = leftendpoint(I) < v ≤ rightendpoint(I)

in(v::Complex, I::TypedEndpointsInterval{:closed,:closed}) = isreal(v) && in(real(v), I)
in(v::Complex, I::TypedEndpointsInterval{:open,:open}) = isreal(v) && in(real(v), I)
in(v::Complex, I::TypedEndpointsInterval{:closed,:open}) = isreal(v) && in(real(v), I)
in(v::Complex, I::TypedEndpointsInterval{:open,:closed}) = isreal(v) && in(real(v), I)

in(::Missing, I::TypedEndpointsInterval{:closed,:closed}) = !isempty(I) && missing
in(::Missing, I::TypedEndpointsInterval{:open,:open}) = !isempty(I) && missing
in(::Missing, I::TypedEndpointsInterval{:closed,:open}) = !isempty(I) && missing
in(::Missing, I::TypedEndpointsInterval{:open,:closed}) = !isempty(I) && missing

isempty(A::TypedEndpointsInterval{:closed,:closed}) = leftendpoint(A) > rightendpoint(A)
isempty(A::TypedEndpointsInterval) = leftendpoint(A) ≥ rightendpoint(A)

isequal(A::TypedEndpointsInterval{L,R}, B::TypedEndpointsInterval{L,R}) where {L,R} = (isequal(leftendpoint(A), leftendpoint(B)) & isequal(rightendpoint(A), rightendpoint(B))) | (isempty(A) & isempty(B))
isequal(A::TypedEndpointsInterval, B::TypedEndpointsInterval) = isempty(A) & isempty(B)

==(A::TypedEndpointsInterval{L,R}, B::TypedEndpointsInterval{L,R}) where {L,R} = (leftendpoint(A) == leftendpoint(B) && rightendpoint(A) == rightendpoint(B)) || (isempty(A) && isempty(B))
==(A::TypedEndpointsInterval, B::TypedEndpointsInterval) = isempty(A) && isempty(B)

function issubset(A::TypedEndpointsInterval, B::TypedEndpointsInterval)
    Al, Ar = endpoints(A)
    Bl, Br = endpoints(B)
    return isempty(A) | ( (Bl ≤ Al) & (Ar ≤ Br) )
end
function issubset(A::TypedEndpointsInterval{:closed,R1} where R1, B::TypedEndpointsInterval{:open,R2} where R2)
    Al, Ar = endpoints(A)
    Bl, Br = endpoints(B)
    return isempty(A) | ( (Bl < Al) & (Ar ≤ Br) )
end
function issubset(A::TypedEndpointsInterval{L1,:closed} where L1, B::TypedEndpointsInterval{L2,:open} where L2)
    Al, Ar = endpoints(A)
    Bl, Br = endpoints(B)
    return isempty(A) | ( (Bl ≤ Al) & (Ar < Br) )
end
function issubset(A::TypedEndpointsInterval{:closed,:closed}, B::TypedEndpointsInterval{:open,:open})
    Al, Ar = endpoints(A)
    Bl, Br = endpoints(B)
    return isempty(A) | ( (Bl < Al) & (Ar < Br) )
end

⊇(A::AbstractInterval, B::AbstractInterval) = issubset(B, A)
⊊(A::AbstractInterval, B::AbstractInterval) = (A ≠ B) & (A ⊆ B)
⊋(A::AbstractInterval, B::AbstractInterval) = (A ≠ B) & (A ⊇ B)

const _interval_hash = UInt == UInt64 ? 0x1588c274e0a33ad4 : 0x1e3f7252

hash(I::TypedEndpointsInterval, h::UInt) = hash(leftendpoint(I), hash(rightendpoint(I), hash(_interval_hash, h)))

minimum(d::TypedEndpointsInterval{:closed}) = infimum(d)
minimum(d::TypedEndpointsInterval{:open}) = throw(ArgumentError("$d is open on the left. Use infimum."))
maximum(d::TypedEndpointsInterval{L,:closed}) where L = supremum(d)
maximum(d::TypedEndpointsInterval{L,:open}) where L = throw(ArgumentError("$d is open on the right. Use supremum."))

extrema(I::TypedEndpointsInterval) = (infimum(I), supremum(I))

# Open and closed at endpoints
isleftclosed(d::TypedEndpointsInterval{:closed}) = true
isleftclosed(d::TypedEndpointsInterval{:open}) = false
isrightclosed(d::TypedEndpointsInterval{L,:closed}) where {L} = true
isrightclosed(d::TypedEndpointsInterval{L,:open}) where {L} = false

# UnitRange construction
# The third is the one we want, but the first two are needed to resolve ambiguities
Base.Slice{T}(i::TypedEndpointsInterval{:closed,:closed,I}) where {T<:AbstractUnitRange,I<:Integer} =
    Base.Slice{T}(minimum(i):maximum(i))
Base.Slice(i::TypedEndpointsInterval{:closed,:closed,I}) where I<:Integer =
    Base.Slice(minimum(i):maximum(i))
function Base.OneTo{T}(i::TypedEndpointsInterval{:closed,:closed,I}) where {T<:Integer,I<:Integer}
    @noinline throwstart(i) = throw(ArgumentError("smallest element must be 1, got $(minimum(i))"))
    minimum(i) == 1 || throwstart(i)
    Base.OneTo{T}(maximum(i))
end
Base.OneTo(i::TypedEndpointsInterval{:closed,:closed,I}) where {I<:Integer} =
    Base.OneTo{I}(i)
UnitRange{T}(i::TypedEndpointsInterval{:closed,:closed,I}) where {T<:Integer,I<:Integer} = UnitRange{T}(minimum(i), maximum(i))
UnitRange(i::TypedEndpointsInterval{:closed,:closed,I}) where {I<:Integer} = UnitRange{I}(i)
range(i::TypedEndpointsInterval{:closed,:closed,I}) where {I<:Integer} = UnitRange{I}(i)

"""
    range(i::ClosedInterval; step, length)
    range(i::ClosedInterval, len::Integer)

Constructs a range of a specified step or length.
"""
range(i::TypedEndpointsInterval{:closed,:closed}; step=nothing, length=nothing) =
    range(leftendpoint(i); stop=rightendpoint(i), step=step, length=length)
range(i::TypedEndpointsInterval{:closed,:closed}, len::Integer) = range(i; length=len)

"""
    range(i::Interval{:closed,:open}; length)
    range(i::Interval{:closed,:open}, len::Integer)

Constructs a range of a specified length with `step=width(i)/length`.
"""
range(i::TypedEndpointsInterval{:closed,:open}; length::Integer) =
    range(leftendpoint(i); step=width(i)/length, length=length)
range(i::TypedEndpointsInterval{:closed,:open}, len::Integer) = range(i; length=len)

"""
    clamp(t, i::ClosedInterval)

Clamp the scalar `t` such that the result is in the interval `i`.
"""
clamp(t, i::TypedEndpointsInterval{:closed,:closed}) =
    clamp(t, leftendpoint(i), rightendpoint(i))

"""
   duration(iv)

calculates the the total number of integers or dates of an integer or date
valued interval. For example, `duration(0..1)` is 2, while `width(0..1)` is 1.
"""
function duration(A::TypedEndpointsInterval{:closed,:closed,T}) where {T<:Integer}
    # TODO: Remove this method in the next breaking release
    Base.depwarn("`duration` will be removed in the next breaking release because it is ill-defined concept. " *
    "If you need 3 from 2..4, replace duration(2..4) with width(2..4)+1", :duration)
    max(0, Int(A.right - A.left) + 1)
end
function duration(A::TypedEndpointsInterval{:closed,:closed,Date})
    # TODO: Remove this method in the next breaking release
    Base.depwarn("`duration` will be removed in the next breaking release because it is ill-defined concept. " *
    "If you need 961 from A=Date(2018,08,08)..Date(2021,03,25), replace duration(A) with Dates.days(width(A))+1", :duration)
    max(0, Dates.days(A.right - A.left) + 1)
end

include("interval.jl")

# convert numbers to intervals
# TODO: These conversions will be removed in the next breaking release (#97)
function convert(::Type{AbstractInterval}, x::Number)
    Base.depwarn("`The conversion number to interval will be removed; construct an interval explicitly, e.g., `x..x`.", :convert)
    x..x
end
function convert(::Type{AbstractInterval{T}}, x::Number) where T
    Base.depwarn("`The conversion number to interval will be removed; construct an interval explicitly, e.g., `x..x`.", :convert)
    convert(AbstractInterval{T}, convert(AbstractInterval, x))
end
function convert(::Type{TypedEndpointsInterval{:closed,:closed}}, x::Number)
    Base.depwarn("`The conversion number to interval will be removed; construct an interval explicitly, e.g., `x..x`.", :convert)
    x..x
end
function convert(::Type{TypedEndpointsInterval{:closed,:closed,T}}, x::Number) where T
    Base.depwarn("`The conversion number to interval will be removed; construct an interval explicitly, e.g., `x..x`.", :convert)
    convert(AbstractInterval{T}, convert(AbstractInterval, x))
end
function convert(::Type{ClosedInterval}, x::Number)
    Base.depwarn("`The conversion number to interval will be removed; construct an interval explicitly, e.g., `x..x`.", :convert)
    x..x
end
function convert(::Type{ClosedInterval{T}}, x::Number) where T
    Base.depwarn("`The conversion number to interval will be removed; construct an interval explicitly, e.g., `x..x`.", :convert)
    convert(AbstractInterval{T}, convert(AbstractInterval, x))
end

include("findall.jl")

end # module


text:
"""
    TimeType

Abstract type for time.
"""
abstract type TimeType end

"""
    ContinuousTime <: TimeType

Abstract type for processes that run in continuous time and need to be discretized.
"""
abstract type ContinuousTime <: TimeType end

"""
    DiscreteTime <: TimeType

Abstract type for processes that run in discrete time without a notion of physical units (i.e. time steps).
"""
abstract type DiscreteTime <: TimeType end

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
module SiriusB

using CSV
using DataFrames
using Dates
using DelimitedFiles
using Dierckx
using FITSIO
import ImageTransformations: center
using Printf
using PyCall
using SkyCoords

export rootdir, datadir, srcdir, notebookdir, paperdir, figuredir, load_or_produce,
       parallactic_angles, gaussian_fit, gaussian_fit_offset, center, table_interpolate, ATMO2020,
       contrast_to_dmag, distance_modulus, SonoraSolar, SonoraMetalRich

# constants associated with Sirius B
const parallax = 374.48958852876103e-3 # arcsecond
const distance = inv(parallax) # pc
const pxscale = 0.01 # arcseconds/px
const auscale = pxscale / parallax # AU/px
const appmag = 9.1 # Lp band from Bonnet-Biduad 2008

# path configuration
rootdir(args...) = joinpath(@__DIR__(), "..", args...)
datadir(args...) = rootdir("data", args...)
srcdir(args...) = rootdir("src", args...)
notebookdir(args...) = rootdir("notebooks", args...)
paperdir(args...) = rootdir("paper", args...)
figuredir(args...) = paperdir("figures", args...)

# simple calculations
contrast_to_dmag(contrast) = -2.5log10(contrast)
distance_modulus(distance) = 5log10(distance) - 5

"""
    load_or_produce(f, filename; force=false)

Load the data from `filename` if it exists, otherwise (or if `force` is true) produce the data using the function `f`. This is a simplified version of the system DrWatson.jl uses for organization and reproducibility. To facilitate loading and saving, only ".fits" and ".csv" files are supported.

# Examples

```julia
using ADI
load_or_produce("flat_res_pca-2.csv") do
    PCA(2)(cube, angles)
end
```
"""
function load_or_produce(func, filename; force=false)
    filename = normpath(filename)
    ext = last(splitext(filename))

    if !isfile(filename) || force
        @info "$(normpath(filename)) not found, producing..."
        res = func()
        SAVER[ext](filename, res)
        @info "file saved to $filename"
        return res
    else
        @info "file loaded from $filename"
        return LOADER[ext](filename)
    end
end

function load_fits(filename; kwargs...)
    PYFITS = pyimport("astropy.io.fits")
    Array(PYFITS.getdata(filename; kwargs...))
end
load_csv(filename; kwargs...) = DataFrame(CSV.File(filename; kwargs...))
const LOADER = Dict(
    ".csv" => load_csv,
    ".fits" => load_fits
)

function save_fits(filename, data; kwargs...)
    mkpath(dirname(filename))
    PYFITS = pyimport("astropy.io.fits")
    PYFITS.writeto(filename, Array(data); overwrite=true, kwargs...)
end
function save_csv(filename, data; kwargs...)
    mkpath(dirname(filename))
    CSV.write(filename, data; kwargs...)
end
const SAVER = Dict(
    ".csv" => save_csv,
    ".fits" => save_fits
)

include("angles.jl")
include("fitting.jl")
include("badpix.jl")
include("tables.jl")

end

max_stars_repo_name: mileslucas/sirius-b

max_issues_repo_name: mileslucas/sirius-b

max_forks_repo_name: mileslucas/sirius-b


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule PCRE2_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("PCRE2")
JLLWrappers.@generate_main_file("PCRE2", UUID("efcefdf7-47ab-520b-bdef-62a2eaa19f15"))
end  # module PCRE2_jll


text:
# function top_down(point::Vector{SVector{3,Float64}}, vec_tri_tet::Vector{SVector{N,Int64}}) where {N}
#     (N == 3) || (N == 4) || error("N is $N. N should be 3 for triangle mesh or 4 for tetrahedral mesh")
#
#     sv_aabb = [calc_obb(point[k]) for k = vec_tri_tet]
#     bb_aabb = [bin_BB_Tree(k, sv_aabb[k]) for k = 1:length(sv_aabb)]
#     return recursive_top_down(bb_aabb)
# end

function recursive_top_down(bb_aabb::Vector{bin_BB_Tree{OBB}})
    n = length(bb_aabb)
    (n == 1) && (return bb_aabb[1])
    (n == 2) && (return bin_BB_Tree(bb_aabb[1], bb_aabb[2]))

    final_AABB = bb_aabb[1].box
    for box_k = bb_aabb
        final_AABB = OBB(final_AABB, box_k.box)
    end
    _, mi_ = findmax(final_AABB.e)

    all_c = [bb_aabb[k].box.c[mi_] for k = 1:n]
    s_perm = sortperm(all_c)
    n_mid = Int64(ceil(n / 2))
    i_pos = s_perm[1:(n_mid - 1)]
    i_neg = s_perm[n_mid:n]

    bb_aabb_1 = bb_aabb[i_pos]
    bb_aabb_2 = bb_aabb[i_neg]
    tree_1 = recursive_top_down(bb_aabb_1)
    tree_2 = recursive_top_down(bb_aabb_2)
    return bin_BB_Tree(tree_1, tree_2)
end


text:
struct CoherentSubspaces{M,S<:Union{Matrix{ComplexF64},NTuple{M,Matrix{ComplexF64}}},T,D} <: ObjectiveFunction
    dims::Dims{D}
    δt::Float64 # timestep (fixed)
    Ut::S # target unitary
    s::BitSet # coherent subspaces
    Ne::Int # effective target Operator dimension (rank of target operator)
    Hd::Matrix{T} # drift Hamiltonian
    Hc::Vector{Matrix{T}} # control Hamiltonian(s)
    u_last::Vector{Float64} # last control amplitudes, ordered by control, and then time
    # Propagator and eigendecompotion storage
    U::Vector{Matrix{ComplexF64}} # individual propagators
    X::Vector{Matrix{ComplexF64}} # forward cumulative propagators
    P::Vector{S} # backwards cumulative propagators
    D::Vector{Vector{Float64}} # eigenvalues associated to each propagator
    V::Vector{Matrix{T}} # eigenvector matrix associated to each propagator
    # Workspace variables
    H::Hermitian{T,Matrix{T}} # temporary storage for full Hamiltonian
    A::Matrix{ComplexF64} # temporary
    Jkj::Matrix{ComplexF64} # temporary
    cisDj::Vector{ComplexF64} # temporary
end

CoherentSubspaces(dims,δt,Ut::NTuple{M,<:Matrix},s,Ne,Hd,Hc,u_last,U,X,P,D,V,H,A,Jkj,cisDj) where M = CoherentSubspaces{M,typeof(Ut),eltype(H),length(dims)}(dims,δt,Ut,s,Ne,Hd,Hc,u_last,U,X,P,D,V,H,A,Jkj,cisDj)
CoherentSubspaces(dims,δt,Ut::Matrix,s,Ne,Hd,Hc,u_last,U,X,P,D,V,H,A,Jkj,cisDj) = CoherentSubspaces{1,typeof(Ut),eltype(H),length(dims)}(dims,δt,Ut,s,Ne,Hd,Hc,u_last,U,X,P,D,V,H,A,Jkj,cisDj)

function CoherentSubspaces(Ut::NTuple{M,<:Operator},s::IntCol,Hd::Operator,Hc::Vector{<:Operator},t::Real,n::Integer,Ne=prod(dims(Hd))) where M
    # TODO: Prob should do some error checking on inputs
    N = prod(dims(Hd))
    m = length(Hc)
    # Make sure we pass dense operators
    Ut_d = complex.(Array.(Ut))
    Hd_d = Array(Hd)
    Hc_d = Array.(Hc)
    # Storage for last control ampitudes, NaN for first run
    u_last = fill(NaN64,m*n)
    # Generate cache for various objects
    H = Hermitian(zeros(promote_eltype(typeof(Hd_d),eltype(Hc_d)),N,N))
    A = Matrix{ComplexF64}(undef,N,N)
    U = [Matrix{ComplexF64}(undef,N,N) for i=1:n]
    X = [Matrix{ComplexF64}(undef,N,N) for i=1:n]
    P = [similar.(Ut_d) for i=1:n]
    for i=1:M; copyto!(P[end][i],Ut_d[i]); end
    # For the exact derivative we need to store eigenvectors and eigenvalues
    D = [Vector{Float64}(undef,N) for i=1:n] # eigenvalues are always real
    V = [similar(H.data) for i=1:n]
    # More temporary storage
    Jkj = Matrix{ComplexF64}(undef,N,N)
    cisDj = Vector{ComplexF64}(undef,N)
    return CoherentSubspaces(dims(Hd),t/n,Ut_d,s,Ne,Hd_d,Hc_d,u_last,U,X,P,D,V,H,A,Jkj,cisDj)
end

function CoherentSubspaces(Ut::Operator,s::IntCol,Hd::Operator,Hc::Vector{<:Operator},t::Real,n::Integer,Ne=prod(dims(Hd)))
    # TODO: Prob should do some error checking on inputs
    N = prod(dims(Hd))
    m = length(Hc)
    # Make sure we pass dense operators
    Ut_d = complex(Array(Ut))
    Hd_d = Array(Hd)
    Hc_d = Array.(Hc)
    # Storage for last control ampitudes, NaN for first run
    u_last = fill(NaN64,m*n)
    # Generate cache for various objects
    H = Hermitian(zeros(promote_eltype(typeof(Hd_d),eltype(Hc_d)),N,N))
    A = Matrix{ComplexF64}(undef,N,N)
    U = [Matrix{ComplexF64}(undef,N,N) for i=1:n]
    X = [Matrix{ComplexF64}(undef,N,N) for i=1:n]
    P = [Matrix{ComplexF64}(undef,N,N) for i=1:n]
    copyto!(P[end],Ut_d)
    # For the exact derivative we need to store eigenvectors and eigenvalues
    D = [Vector{Float64}(undef,N) for i=1:n] # eigenvalues are always real
    V = [similar(H.data) for i=1:n]
    # More temporary storage
    Jkj = Matrix{ComplexF64}(undef,N,N)
    cisDj = Vector{ComplexF64}(undef,N)
    return CoherentSubspaces(dims(Hd),t/n,Ut_d,s,Ne,Hd_d,Hc_d,u_last,U,X,P,D,V,H,A,Jkj,cisDj)
end

function objective(O::CoherentSubspaces,u)
    Uf = O.X[end]
    # Calculate forward propagators
    calc_fprops!(O,u)
    # Calculate coherent subspaces norm infidelity:
    # fₑ = 1 - |⟨Ut,Uf⟩|cs/N
    return 1 - inner_cs(O.Ut,Uf,O.s)/O.Ne
end

function gradient!(O::CoherentSubspaces,fp,u)
    Uf = O.X[end]
    n = length(O.U); m = length(O.Hc)
    # Calculate forward and backward propagators
    calc_fprops!(O,u)
    calc_bprops!(O)
    # Calculate exact derivative of coherent subspaces norm infidelity
    x,y = _inner_cs_1(Uf,O.Ut,O.s)
    for j = 1:n
        O.cisDj .= cis.(O.D[j])
        for k = 1:m
            Jmat!(O,k,j)
            j==1 ? copyto!(O.A,O.Jkj) : mul!(O.A,O.Jkj,O.X[j-1])
            fp[(k-1)*n+j] = -inner_cs_grad(O.P[j],O.A,x,y,O.s)/O.Ne
        end
    end
    return fp
end


text:
# map_ranges.jl

range(m::AbstractMap) = m * domain(m)

"""
Decide whether a given point `y` is in the range of the map `y=f(x)`, possibly
up to a given tolerance if a tolerance is supplied.
"""
in_range(m::AbstractMap, y) = in(y, range(m))

in_range(m::AbstractMap, y, tolerance) = approx_in(y, range(m), tolerance)

# We can't include domain and ranges of maps in the maps code, because it is
# interpreted before the domain code. Therefore, for a few specific maps, we
# define the domain and range here.

domain(m::IdentityMap{T}) where {T} = FullSpace{T}()
range(m::IdentityMap) = domain(m)

domain(m::ConstantMap{T,S}) where {T,S} = FullSpace{S}()
range(m::ConstantMap) = Point(constant(m))

domain(m::ProductMap) = ProductDomain(map(domain, elements(pm))...)
range(m::ProductMap) = ProductDomain(map(range, elements(pm))...)