text:
function Base.show(io::IO, ::MIME"text/plain", node::TreeNode)
    if !get(io, :compact, false)
        node.isroot ? println(io, "Node $(node.label) (root)") : println(io, "Node $(node.label): ")
        node.isroot ? println(io, "Ancestor : `nothing` (root)") : println(io, "Ancestor: $(node.anc.label), branch length = $(node.tau)")
        print(io, "$(length(node.child)) children: $([x.label for x in node.child])")
    end
    return nothing
end
function Base.show(io::IO, node::TreeNode)
	print(io, "$(typeof(node)): $(node.label)")
	node.isroot && print(io, " (root)")
	return nothing
end

function Base.show(io::IO, t::Tree{T}) where T
	nn = length(nodes(t))
	nl = length(leaves(t))
	long = begin
		base = "Tree{$T}: "
		base *= nn > 1 ? "$nn nodes, " : "$nn node, "
		base *= nl > 1 ? "$nl leaves" : "$nl leaf"
		base
	end
	if length(long) < 0.8*displaysize(io)[2]
		print(io, long)
		return nothing
	end

	short = begin
		base = "Tree w. "
		base *= nl > 1 ? "$nl leaves" : "$nl leaf"
		base
	end
	print(io, short)
	return nothing
end
function Base.show(io::IO, ::MIME"text/plain", t::Tree; maxnodes=40)
    if length(nodes(t)) < maxnodes
        print_tree_ascii(io, t)
    else
    	show(io, t)
    end
end

function print_tree_(io, node, cdepth; vindent=2, hindent=5, hoffset=0, maxdepth=5)
    hspace = ""
    for i in 1:hindent
        hspace *= "-"
    end
    offset = ""
    for i in 1:hoffset
        offset *= " "
    end
    if cdepth < maxdepth
    	println(io, "$offset $hspace $(node.label):$(node.tau)")
    elseif cdepth == maxdepth
    	if node.isleaf
    		println(io, "$offset $hspace $(node.label):$(node.tau)")
    	else
    		println(io, "$offset $hspace $(node.label):$(node.tau) ...")
    	end
    end
        #
    if cdepth <= maxdepth
        if !node.isleaf
            for c in node.child
                for i in 1:vindent
                    cdepth < maxdepth && println(io, "$offset $(" "^hindent)|")
                end
                print_tree_(io, c, cdepth + 1, vindent=vindent, hindent=hindent, hoffset=hoffset+hindent, maxdepth=maxdepth)
            end
        end
        #
    end
end
function print_tree(io, node::TreeNode; vindent=2, hindent=5, maxdepth=5)
    print_tree_(io, node, 1, vindent=vindent, hindent=hindent, hoffset=0, maxdepth=maxdepth)
end

print_tree(io, t::Tree; vindent=2, hindent=5, maxdepth=5) = print_tree(io, t.root; vindent=2, hindent=5, maxdepth=maxdepth)

"""
    print_tree_ascii(io, t::Tree)

Julia implementation of Bio.Phylo.draw_ascii function:
https://github.com/biopython/biopython/blob/master/Bio/Phylo/_utils.py
"""
function print_tree_ascii(io, t::Tree)
    column_width = 80
    taxa = [ node.label for node in POTleaves(t)] #need leaves in Post-Order Traversal
    max_label_width = maximum([length(taxon) for taxon in taxa])
    drawing_width = column_width - max_label_width - 1
    drawing_height = 2 * length(taxa) - 1

    function get_col_positions(t::Tree)
        depths = [divtime(node, t.root) for node in values(t.lnodes)]
        # If there are no branch lengths, assume unit branch lengths
        if ismissing(maximum(depths))
            println(io, "\n not all branch lengths known, assuming identical branch lengths")
            depths = [node_depth(node) for node in values(t.lnodes)]
        end
        # Potential drawing overflow due to rounding -- 1 char per tree layer
        fudge_margin = ceil(Int64, log2(length(taxa)))
        if maximum(depths)==0
            cols_per_branch_unit = (drawing_width - fudge_margin)
        else
            cols_per_branch_unit = (drawing_width - fudge_margin) / maximum(depths)
        end
        return Dict(zip(keys(t.lnodes), round.(Int64,depths*cols_per_branch_unit .+2.0)))
    end

    function get_row_positions(t::Tree)
        positions = Dict{Any, Int64}(zip(taxa, 2 *(1:length(taxa)) ) )
        function calc_row(clade::TreeNode)
            for subclade in clade.child
                if !haskey(positions, subclade.label)
                    calc_row(subclade)
                end
            end
            if !haskey(positions, clade.label)
                positions[clade.label] = floor(Int64, (positions[clade.child[1].label] + positions[clade.child[end].label])/2)
            end
        end
        calc_row(t.root)
        return positions
    end

    col_positions = get_col_positions(t)
    row_positions = get_row_positions(t)
    char_matrix = [[" " for x in 1:(drawing_width+1)] for y in 1:(drawing_height+1)]

    function draw_clade(clade::TreeNode, startcol::Int64)
        thiscol = col_positions[clade.label]
        thisrow = row_positions[clade.label]
        # Draw a horizontal line
        for col in startcol:thiscol
            char_matrix[thisrow][col] = "_"
        end
        if !isempty(clade.child)
            # Draw a vertical line
            toprow = row_positions[clade.child[1].label]
            botrow = row_positions[clade.child[end].label]
            for row in (toprow+1):botrow
                char_matrix[row][thiscol] = "|"
            end
            # Short terminal branches need something to stop rstrip()
            if (col_positions[clade.child[1].label] - thiscol) < 2
                char_matrix[toprow][thiscol] = ","
            end
            # Draw descendents
            for child in clade.child
                draw_clade(child, thiscol + 1)
            end
        end
    end
    draw_clade(t.root, 1)
    # Print the complete drawing
    for i in 1:length(char_matrix)
        line = rstrip(join(char_matrix[i]))
        # Add labels for terminal taxa in the right margin
        if i % 2 == 0
            line = line * " " * strip(taxa[round(Int64, i/2)]) #remove white space from labels to make more tidy
        end
        println(io, line)
    end
end

"""
    check_tree(t::Tree; strict=true)
- Every non-leaf node should have at least one child (two if `strict`)
- Every non-root node should have exactly one ancestor
- If n.child[...] == c, c.anc == n is true
- Tree has only one root
"""
function check_tree(tree::Tree; strict=true)
    labellist = Dict{String, Int64}()
    nroot = 0
    flag = true
    for n in values(tree.lnodes)
        if !n.isleaf && length(n.child)==0
        	(flag = false) || (@warn "Node $(n.label) is non-leaf and has no child.")
        elseif !n.isroot && n.anc == nothing
        	(flag = false) || (@warn "Node $(n.label) is non-root and has no ancestor.")
        elseif strict && length(n.child) == 1
            if !(n.isroot && n.child[1].isleaf)
        	   (flag = false) || (@warn "Node $(n.label) has only one child.")
            end
        elseif length(n.child) == 0 && !haskey(tree.lleaves, n.label)
            (flag = false) || (@warn "Node $(n.label) has no child but is not in `tree.lleaves`")
        end
        for c in n.child
            if c.anc != n
                (flag = false) || (@warn "One child of $(n.label) does not satisfy `c.anc == n`.")
            end
        end
        if get(labellist, n.label, 0) == 0
            labellist[n.label] = 1
        else
            labellist[n.label] += 1
            (flag = false) || (@warn "Label $(n.label) already exists!")
        end
        if n.isroot
            nroot += 1
        end
    end
    if nroot > 1
        (flag = false) || (@warn "Tree has multiple roots")
    elseif nroot ==0
        (flag = false) || (@warn "Tree has no root")
    end
    return flag
end

make_random_label(base="NODE") = make_random_label(base, 8)
make_random_label(base, i) = "$(base)_$(randstring(i))"

"""
    create_label(t::Tree, base="NODE")
Create new node label in tree `t` with format `\$(base)_i` with `i::Int`.
"""
function create_label(t::Tree, base="NODE")
    label_init = 1
    pattern = Regex(base)
    for n in values(t.lnodes)
        if match(pattern, n.label)!=nothing && parse(Int64, n.label[length(base)+2:end]) >= label_init
            label_init = parse(Int64, n.label[length(base)+2:end]) + 1
        end
    end
    return "$(base)_$(label_init)"
end

function set_unique_label!(node::TreeNode, t::Tree; delim = "__")
	id = randstring(8)
	node.label *= delim * id
	while haskey(t.lnodes, node.label)
		node.label = node.label[1:end-5] * randstring(8)
	end
end

"""
    map_dict_to_tree!(t::Tree{MiscData}, dat::Dict; symbol=false, key = nothing)
Map data in `dat` to nodes of `t`. All node labels of `t` should be keys of `dat`. Entries of `dat` should be dictionaries, or iterable similarly, and are added to `n.data.dat`.
If `!isnothing(key)`, only a specific key of `dat` is added. It's checked for by `k == key || Symbol(k) == key` for all keys `k` of `dat`.
If `symbol`, data is added to nodes of `t` with symbols as keys.
"""
function map_dict_to_tree!(t::Tree{MiscData}, dat::Dict; symbol=false, key = nothing)
    for (name, n) in t.lnodes
        for (k,v) in dat[name]
            if !isnothing(key) && (k == key || Symbol(k) == key)
                n.data.dat[key] = v
            elseif isnothing(key)
                n.data.dat[symbol ? Symbol(k) : k] = v
            end
        end
    end
    nothing
end
"""
    map_dict_to_tree!(t::Tree{MiscData}, dat::Dict, key)
Map data in `dat` to nodes of `t`. All node labels of `t` should be keys of `dat`. Entries of `dat` corresponding to `k` are added to `t.lnodes[k].data.dat[key]`.
"""
function map_dict_to_tree!(t::Tree{MiscData}, dat::Dict, key)
    for (name, n) in t.lnodes
        n.data.dat[key] = dat[name]
    end
end

"""
    rand_times!(t, p)
Add random branch lengths to tree.
"""
function rand_times!(t, p)
    for n in values(t.lnodes)
        if !n.isroot
            n.tau = rand(p)
        end
    end
end


text:
# Autogenerated wrapper script for ssht_jll for i686-linux-gnu
export libssht

using FFTW_jll
JLLWrappers.@generate_wrapper_header("ssht")
JLLWrappers.@declare_library_product(libssht, "libssht.so")
function __init__()
    JLLWrappers.@generate_init_header(FFTW_jll)
    JLLWrappers.@init_library_product(
        libssht,
        "lib/libssht.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
using HTTP

include("utils.jl")

function request(; url::String,
                 query::Union{Nothing, String}=nothing,
                 filter::Union{Nothing, Dict{String, String}}=nothing,
                 offset::Union{Nothing, Int}=nothing,
                 limit::Union{Nothing, Int}=nothing,
                 sample::Union{Nothing, Int}=nothing,
                 sort::Union{Nothing, String}=nothing,
                 order::Union{Nothing, String}=nothing,
                 facet::Union{Nothing, Dict{String, Int}}=nothing,
                 select::Union{Nothing, String}=nothing,
                 cursor::Union{Nothing, String}=nothing)

        if !isnothing(offset)
                if offset > 10000 # max offset
                        error("Offset; Max: 10000")
                end
        end

        if !isnothing(sample)
                if sample > 100 # max sample size
                        error("Sample size; Max: 100")
                end
        end

        headers = []

        try
                push!(headers, ("User-Agent", string("mailto:", ENV["MAILTO"])))
                push!(headers, ("X-USER-AGENT", string(ENV["MAILTO"])))
        catch e
                #println(e)
        end

        try
                push!(headers, ("Crossref-Plus-API-Token", string(ENV["CROSSREF_PLUS_API_TOKEN"])))
        catch e
                #println(e)
        end

        payload = Dict("query" => query,
                       "filter" => filter_handler(filter),
                       "offset" => offset,
                       "rows" => limit,
                       "sample" => sample,
                       "sort" => sort,
                       "order" => order,
                       "facet" => filter_handler(facet),
                       "select" => select,
                       "cursor" => cursor)

        for (k, v) in payload
                if isnothing(v)
                        delete!(payload, k)
                end
        end

        return HTTP.get(url, headers; query=payload)
end

max_stars_repo_name: naustica/Pitaya

max_issues_repo_name: naustica/Pitaya

max_forks_repo_name: naustica/Pitaya


text:
using VisualEFM
using Statistics
using Test

@testset "VisualEFM.jl" begin

    f = "test_figures/".*["rpm000.jpg", "rpm150.jpg", "rpm175.jpg", "rpm225.jpg"]
    im_original = read_image.(f, roi=[500:2050,:])
    mask = binarymask("test_figures/teste_mascara.png")
    im = read_image.(f, roi=[500:2050,:], mask=mask)

    @testset "Common functions" begin
        @test_throws MethodError read_image(f)
        @test length(im_original) == 4
        @test size(im_original[1]) == (1551, 1456)
        @test typeof(mask) == BitArray{2}
        @test mask[1:10,1:10] == fill(false, (10,10))
        @test size(im[1]) == (1551, 1456)
        @test length(im) == 4
    end

    _,mo = erosion_one(im_original[3], im_original[1], showPlot=false)
    _,m = erosion_one(im[3], im[1], showPlot=false)
    a,_ = erosion_colormap(im_original, [150,175,225], showPlot=false)
    b,_ = erosion_colormap(im, [150,175,225], showPlot=false)

    @testset "Sand erosion" begin
        @test count(mo) > 231400 & count(mo) < 231800
        @test count(m) > 102620 & count(m) < 102820
        @test_throws DimensionMismatch erosion_anim("test.gif",im_original[2:4], im_original[1:2])
        @test_throws DimensionMismatch erosion_colormap(im_original, [150,175])
        @test_throws DimensionMismatch erosion_colormap(im, [150,175])
        @test maximum(x->isnan(x) ? -Inf : x, a) == 225
        @test maximum(x->isnan(x) ? -Inf : x, b) == 225
        @test minimum(x->isnan(x) ? Inf : x, a) == 150
        @test minimum(x->isnan(x) ? Inf : x, b) == 150
    end

    f_smoke = "test_figures/".*["frame1.png", "frame2.png", "frame3.png"]
    im_smoke = read_image.(f_smoke)
    s_all,s_m = smoke_statsmap(mean, im_smoke, nothing, showPlot=false)
    _,s_sd = smoke_statsmap(std, im_smoke, nothing, showPlot=false)

    @testset "Smoke" begin
        @test_throws MethodError smoke_anim("as.gif", im_smoke, im_smoke)
        @test_throws MethodError smoke_statsmap(mean, im_smoke, im_smoke)
        @test s_m[250,300] ≈ (s_all[1][250,300]+s_all[2][250,300]+s_all[3][250,300])/3
        @test s_m[250,300] ≈ mean([s_all[1][250,300],s_all[2][250,300],s_all[3][250,300]])
        @test s_sd[250,300] ≈ std([s_all[1][250,300],s_all[2][250,300],s_all[3][250,300]])
        @test s_m[250,500] ≈ mean([s_all[1][250,500],s_all[2][250,500],s_all[3][250,500]])
        @test s_sd[250,500] ≈ std([s_all[1][250,500],s_all[2][250,500],s_all[3][250,500]])
    end

end


text:
# Inspired from SOSTools
export randpsd, randsos

cfld(x::NTuple{2,Int}, n) = (cld(x[1], n), fld(x[2], n))

function sub_extdegree(X::AbstractVector{<:MP.AbstractMonomial}, vars)
    if isempty(X)
        return (0, 0)
    else
        return extrema(map(mono -> sum(var -> MP.degree(mono, var), vars), X))
    end
end

# Cheap approximation of the convex hull as the approximation of:
#
# z such that mindegree < sum(z) < maxdegree
# |\
# |#\ <-------- sum(z) = maxdegree
# |##\
# |\<-\-------- sum(z) = mindegree
# | \##\
# +---------
#
# and:
#
# z such that minmultideg < z < maxmultideg
# | +----+ <--- maxmultidegree
# | |####|
# | |####|
# | +----+
# | ^---------- minmultidegree
# +---------

function _sub_half_newton_polytope(X::AbstractVector{<:MP.AbstractMonomial},
                                   extdeg, exp, vars)
    mindeg, maxdeg = cfld(extdeg, 2)
    n = length(vars)
    minmultideg, maxmultideg = Vector{Int}(undef, n), Vector{Int}(undef, n)
    for i in 1:n
        exponent_i(mono) = exp(mono, i)
        minmultideg[i] = cld(mapreduce(exponent_i, min, X), 2)
        maxmultideg[i] = fld(mapreduce(exponent_i, max, X), 2)
    end
    MP.monomials(vars, mindeg:maxdeg,
                 mono -> all(i -> minmultideg[i] <= exp(mono, i) <= maxmultideg[i],
                             1:n))
end

function sub_half_newton_polytope(X::AbstractVector{<:MP.AbstractMonomial},
                                  vars)
    _sub_half_newton_polytope(X, sub_extdegree(X, vars),
                              (mono, i) -> MP.degree(mono, vars[i]), vars)
end

# Multipartite
# TODO we might do this recursively : do 2 parts, merge them, merge with next
#      one and so on so that the filter at the end prunes more.
function half_newton_polytope(X::AbstractVector, parts::Tuple)
    if !all(i -> all(j -> i == j || isempty(parts[i] ∩ parts[j]),
                     1:length(parts)),
            1:length(parts))
        throw(ArgumentError("Parts are not disjoint in multipartite Newton polytope estimation: $parts"))
    end
    # Some variables might be in no part...
    missing = setdiff(MP.variables(X), reduce(union, parts))
    if isempty(missing)
        all_parts = parts
    else
        # in that case, create a part with the missing ones
        all_parts = (parts..., missing)
    end
    if length(all_parts) == 1
        # all variables on same part, fallback to shortcut
        return half_newton_polytope(X, tuple())
    end
    monovecs = map(vars -> sub_half_newton_polytope(X, vars), all_parts)
    # Cartesian product of the newton polytopes of the different parts
    product = [prod(monos) for monos in Iterators.product(monovecs...)]
    mindeg, maxdeg = cfld(MP.extdegree(X), 2)
    # We know that the degree inequalities are satisfied variable-wise and
    # part-wise but for all variables together so we filter with that
    return MP.monovec(filter(mono -> mindeg <= MP.degree(mono) <= maxdeg, product))
end
# Shortcut for more efficient `extdeg` and `exp` function in case all the
# variables are in the same part
function half_newton_polytope(X::AbstractVector, parts::Tuple{})
    return _sub_half_newton_polytope(X, MP.extdegree(X),
                                     (mono, i) -> MP.exponents(mono)[i],
                                     MP.variables(X))
end

function monomials_half_newton_polytope(X::AbstractVector, parts)
    half_newton_polytope(MP.monovec(X), parts)
end

function randpsd(n; r=n, eps=0.1)
    Q = randn(n,n)
    d = zeros(Float64, n)
    d[1:r] = eps .+ abs.(randn(r))
    Q' * Diagonal(d) * Q
end

function _randsos(X::AbstractVector{<:MP.AbstractMonomial}; r=-1, monotype=:Classic, eps=0.1)
    if monotype == :Classic
        x = monomials_half_newton_polytope(X, tuple())
    elseif monotype == :Gram
        x = X
    else
        throw(ArgumentError("Monotype $monotype not known"))
    end
    n = length(x)
    if r < 0
        r = n
    end
    GramMatrix(randpsd(n, r=r, eps=eps), x)
end

randsos(X::AbstractVector; kws...) = _randsos(MP.monovec(X); kws...)


text:
import Base.isless
import Base.isequal
import Base.==
import Base.hash

import Base.show

"""
SAS has fields species::String and pos::Int
it represents an atom slot in a species
"""
struct SAS
    species::String
    pos::Int64
end

function Base.isless(x::SAS,y::SAS)
    if( x.species < y.species ); return true; end
    if( x.species == y.species)
        return x.pos < y.pos
    else
        return false
    end
end
function ==(x::SAS,y::SAS)
    return isequal(x.species,y.species) && x.pos==y.pos
end
function Base.hash(x::SAS)
    return hash(x.species)+hash("xxx$(x.pos)")
end
function Base.isequal(x::SAS,y::SAS)
    return isequal(x.species,y.species) && x.pos==y.pos
end

"""
RAS has fields flux::String, sas::SAS and rpos::Int
it represents an atom slot in a reactand molecule of a specific flux
"""
struct RAS
    flux::String
    sas::SAS
    rpos::Int64
end
# sort: 1. flux 2. species, 3. (!!) rpos 4. sas.pos
function Base.isless(x::RAS,y::RAS)
    if( x.flux < y.flux ); return true; end
    if( x.flux == y.flux )
        if( x.sas.species < y.sas.species ); return true; end
        if( x.species == y.species)
            if(x.rpos < y.rpos)
                return true
            end
            if(x.rpos == y.rpos)
                return x.sas.pos < y.sas.pos
            end
        end
        return false
    end
    return false
end
function ==(x::RAS,y::RAS)
    return  x.flux==y.flux && isequal(x.sas,y.sas) && x.rpos==y.rpos
end
function Base.hash(x::RAS)
    return hash("$(x.flux)ab$(x.sas.species)xx$(x.rpos)yy$(x.sas.pos)")
end
function Base.isequal(x::RAS,y::RAS)
    return x.flux==y.flux && Base.isequal(x.sas,y.sas) && x.rpos==y.rpos
end
function Base.show(io::IO, x::RAS)
  print(io, "RAS:($(x.flux),$(x.sas.species),$(x.sas.pos));rpos=$(x.rpos))")
end

"""
PAS has fields flux::String, sas::SAS and rpos::Int
it represents an atom slot in a product molecule of a specific flux
"""
struct PAS
    flux::String
    sas::SAS
    ppos::Int64
end
# sort: 1. flux, 2. species, 3. (!!) rpos 4. sas.pos
function Base.isless(x::PAS,y::PAS)
    if(x.flux <  y.flux); return true; end
    if(x.flux == y.flux)
        if( x.sas.species < y.sas.species ); return true; end
        if( x.sas.species == y.sas.species)
            if(x.ppos < y.ppos)
                return true
            end
            if(x.ppos == y.ppos)
                return x.sas.pos < y.sas.pos
            end
        end
        return false
    end
    return false
end

function ==(x::PAS,y::PAS)
    return x.flux==y.flux && Base.isequal(x.sas,y.sas) && x.ppos==y.ppos
end
function Base.hash(x::PAS)
    return return hash("$(x.flux)ab$(x.sas.species)xx$(x.ppos)yy$(x.sas.pos)")
end
function Base.isequal(x::PAS,y::PAS)
    return x.flux==y.flux && Base.isequal(x.sas,y.sas) && x.ppos==y.ppos
end
function Base.show(io::IO, x::PAS)
  print(io, "PAS:($(x.flux),$(x.sas.species),$(x.sas.pos));ppos=$(x.ppos))")
end

"""
MassDistributionMeasurement has fields sass::Set{SAS} and data::Dict{Int64,Array{Float64,1}}
It represents a mass distribution measurement for the set of SAS given by .sass
.data contains measurement data for the different isotopic weights, i.e. the
key is the weight, and the value is the data.
Format of .data[m] is usually either:
  .data[m][1] : mean of measurement
  .data[m][2] : std. dev. of measurement
or (in case of multiple data points):
  .data[m][i] : i'th measurement
"""
struct MassDistributionMeasurement
    sass::Set{SAS}
    data::Dict{Int64,Array{Float64,1}}
end

function MassDistributionMeasurement( sass::Set{SAS} , data::Vector{Tuple{Int64,Array{Float64}}} )
    dict_data = Dict{Int64,Array{Float64,1}}()
    for zi in 1:length(data)
        dict_data[data[zi][1]] = data[zi][2]
    end
    for zm in 0:length(sass)
        if( ! contains(==,collect(keys(dict_data)),zm) )
            dict_data[zm] = [0.0;0.0;0.0]
        end
    end
    return MassDistributionMeasurement(sass,dict_data)
end

#struct ConfiguredSAS
#    sass::Set{SAS}
"""
ConfiguredSEMU contains fields .semu::Set{SAS} and .config::Vector{Int64}
.config specifies unlabeled / labeled atoms (0/1)
"""
struct ConfiguredSEMU
    semu::Set{SAS}
    config::Vector{Int64}
end
function ==(x::ConfiguredSEMU,y::ConfiguredSEMU)
    return Base.isequal(x.semu,y.semu) && x.config==y.config
end
function Base.hash(x::ConfiguredSEMU)
    return hash(x.semu)+hash(x.config)
end
function Base.isequal(x::ConfiguredSEMU,y::ConfiguredSEMU)
    return Base.isequal(x.semu,y.semu) && x.config==y.config
end
function Base.show(io::IO, x::SAS)
    print(io, "SAS:($(x.species),$(x.pos))")
end

ConfiguredSEMUMixture = Dict{ConfiguredSEMU,Float64}

SpeciesEMU = Set{SAS}
ReactandEMU = Set{RAS}
ProductEMU = Set{PAS}

function constructSpeciesEMU(species::String,pos::Vector{Integer})
    return constructSpeciesEMU(species,convert(Array{Int64},pos))
end

function constructSpeciesEMU(species::String,pos::Vector{Int64})
    semu1 = SpeciesEMU()
    for zi=1:length(pos)
        push!(semu1,SAS(species,pos[zi]))
    end
    return semu1
end

function Base.isless(a::SpeciesEMU,b::SpeciesEMU)
    if(length(a)<length(b)); return true; end
    if(length(a)==length(b))
        as = sort(collect(a))
        bs = sort(collect(b))
        for i in 1:length(a)
            if(as[i]<bs[i])
                return true
            else
                if(bs[i]<as[i])
                    return false
                end
            end
        end
        return false
    end
    return false
end

function Base.show(io::IO, x::SpeciesEMU)
    xs = sort(collect(x))
    if( all( map( x1 -> x1.species == xs[1].species , xs ) ) )
        print(io,"SEMU::$(xs[1].species)::(")
        str_slots = ""
        for xi in xs; str_slots = string(str_slots,"$(xi.pos),") end
        print(io,str_slots[1:length(str_slots)-1])
        print(io,")")
    else
        print(io,"SEMU::")
        for xi in xs
            print(io,"($(xi.species),$(xi.pos))")
        end
        #print(io,"")
    end
end

function get_positions_of_semu(s::SpeciesEMU)
    if(!isinsinglespecies(s)); error("getpositionofsemu requires single species EMU!"); end
    return map( si->si.pos  ,s)
end

function remu2semu(re::ReactandEMU)
    return SpeciesEMU( map( x -> x.sas , re ) )
end
function pemu2semu(pe::ProductEMU)
    return SpeciesEMU( map( x -> x.sas , pe ) )
end

function isinsinglespecies(se::SpeciesEMU)
    if(isempty(se)); return true; end
    e_species = collect(se)[1].species
    return all( map( x -> x.species==e_species , collect(se) ) )
end

"""""
matching_mass_combinations( target_mass::Int64 , max_individual_masses::Vector{Int64} )
return Vector{Vector{Int64}}

returns all combinations of masses which (i) sum up to target_mass, and where
mass i has value inbetween 0 and max_individual_masses[i]

Example:
matching_mass_combinations( 2 , [1,3,2] )
returns: [0,0,2] , [0,1,1] , [0,2,0] , [1,0,1] , [1,1,0]

matching_mass_combinations( 3 , [1,4] )
return [0,4] , [1,3]
etc..
"""
function matching_mass_combinations( target_mass::Int64 , max_individual_masses::Vector{Int64} )
    # compute recursively: take the first mim value, and loop over all possibilities..
    combis = Vector{Vector{Int64}}()

    if(length(max_individual_masses)==1)
        # end of recursion..
        if(max_individual_masses[1]>=target_mass)
            return [[target_mass]]
        else
            return Array{Array{Int64,1},1}(0)
        end
    end

    mim_rec = max_individual_masses[2:end]
    for mi_w in 0:max_individual_masses[1]
        if(target_mass-mi_w<0); continue; end
        mmc_i = matching_mass_combinations( target_mass-mi_w , mim_rec )
        mmc_i_2 = map( xi -> [mi_w;xi] , mmc_i )
        append!( combis , mmc_i_2 )
    end
    return combis
end

function benchmarks_01()

    f_randi = (max_i) -> trunc(Int,rand()*max_i*0.99999)+1

    flux_names    = map( x -> string("v_",randstring()) ,1:100)
    species_names = map( x -> string("s_",randstring()) ,1:100)

    sas1::Vector{SAS} = Vector{SAS}()
    for zi=1:200
        push!( sas1 , SAS(species_names[f_randi(100)] , f_randi(12) ) )
    end

    semus1::Vector{SpeciesEMU} = Vector{SpeciesEMU}()
    for zi=1:200
        i_spc   = species_names[f_randi(100)]

        i_slots_pre = collect( 1:f_randi(12) )
        i_slots     = shuffle(i_slots_pre)
        i_sl        = i_slots[1:f_randi(length(i_slots))]
        sas_i::Vector{SAS} = Vector{SAS}()
        for zj=1:length(i_sl); push!(sas_i,SAS(i_spc,i_sl[zj])); end
        semui       = SpeciesEMU(Set(sas_i))
        push!( semus1 , semui )
    end

    # test dicts:
    tic()
    d_semu::Dict{SpeciesEMU , Float64} = Dict{SpeciesEMU , Float64}()
    for zi=1:1e4
        d_semu = Dict{SpeciesEMU , Float64}()
        for zj=1:100
            d_semu[semus1[zj]] = rand()
        end
    end
    t_dicts_1 = toc()

    # test dicts 2:
    semus_to_choose = collect( keys( d_semu ) )

    tic()
    sum_1::Float64 = 0.0
    for zi=1:1000000
        sum_1 += 0.0001 * d_semu[ semus_to_choose[ f_randi(length(semus_to_choose)) ] ]
    end
    t_dicts_2 = toc()
    println("\ndict acccess speed: $(1000000/t_dicts_2) / s  \n")
end


text:
# Memory-level Operations
# =======================

module Mem

function copy(dst, src, n::Integer)
    return copy(pointer(dst), pointer(src), n)
end

function copy(dst::Ptr, src::Ptr, n::Integer)
    return ccall(:memcpy, Ptr{Void}, (Ptr{Void}, Ptr{Void}, Csize_t), dst, src, n)
end

end


text:
import HTTP2
using HTTP2.Frame
using Test

## Run `nghttpd --verbose --no-tls 9000` to make this test pass
(headers, body) = HTTP2.request(ip"127.0.0.1", 9000, "/")

println()
println("Results of the request")
println("======================")
println()
println("Headers")
println("======================")
for header in headers
    print(ascii(header[1]))
    print(": ")
    print(ascii(header[2]))
    print("\n")
end
println()
println("Body")
println("======================")
println(ascii(body))


text:
@testset "function counter" begin
    prob = MultivariateProblems.UnconstrainedProblems.examples["Rosenbrock"]

    let
        global fcount = 0
        global fcounter
        function fcounter(reset::Bool = false)
            if reset
                fcount = 0
            else
                fcount += 1
            end
            fcount
        end
        global gcount = 0
        global gcounter
        function gcounter(reset::Bool = false)
            if reset
                gcount = 0
            else
                gcount += 1
            end
            gcount
        end
        global hcount = 0
        global hcounter
        function hcounter(reset::Bool = false)
            if reset
                hcount = 0
            else
                hcount += 1
            end
            hcount
        end
    end

    f(x) = begin
        fcounter()
        MVP.objective(prob)(x)
    end
    g!(out, x) = begin
        gcounter()
        MVP.gradient(prob)(out, x)
    end
    h!(out, x) = begin
        hcounter()
        MVP.hessian(prob)(out, x)
    end

    ls = LineSearches.HagerZhang()

    for solver in (AcceleratedGradientDescent, BFGS, ConjugateGradient,
                   GradientDescent, LBFGS, MomentumGradientDescent,
                   NGMRES, OACCEL)
        fcounter(true); gcounter(true)
        res = Optim.optimize(f, g!, prob.initial_x,
                             solver(linesearch = ls))
        @test fcount == Optim.f_calls(res)
        @test gcount == Optim.g_calls(res)
    end

    for solver in (Newton(linesearch = ls), NewtonTrustRegion())
        fcounter(true); gcounter(true); hcounter(true)
        res = Optim.optimize(f,g!, h!, prob.initial_x,
                             solver)
        @test fcount == Optim.f_calls(res)
        @test gcount == Optim.g_calls(res)
        @test hcount == Optim.h_calls(res)
    end

    # Need to define fg! and hv! for KrylovTrustRegion
    fg!(out,x) = begin
        g!(out,x)
        f(x)
    end
    hv!(out, x, v) = begin
        n = length(x)
        H = Matrix{Float64}(undef, n, n)
        h!(H, x)
        out .= H * v
    end
    begin
        solver = Optim.KrylovTrustRegion()
        fcounter(true); gcounter(true); hcounter(true)
        df = Optim.TwiceDifferentiableHV(f, fg!, hv!, prob.initial_x)
        res = Optim.optimize(df, prob.initial_x, solver)
        @test fcount == Optim.f_calls(res)
        @test gcount == Optim.g_calls(res)
        @test hcount == Optim.h_calls(res)
    end
end


text:
using Gtk, Colors, Graphics, Dates

const radius = 300
const win = GtkWindow("Clock", radius, radius)
const can = GtkCanvas()
push!(win, can)

global drawcontext = []

function drawline(ctx, l, color)
    isempty(l) && return
    p = first(l)
    move_to(ctx, p.x, p.y)
    set_source(ctx, color)
    for i = 2:length(l)
        p = l[i]
        line_to(ctx, p.x, p.y)
    end
    stroke(ctx)
end

function clockbody(ctx)
    set_coordinates(ctx, BoundingBox(0, 100, 0, 100))
    rectangle(ctx, 0, 0, 100, 100)
    set_source(ctx, colorant"yellow")
    fill(ctx)
    set_source(ctx, colorant"blue")
    arc(ctx, 50, 50, 45, 45, 360)
    stroke(ctx)
    for hr in 1:12
        radians = hr * pi / 6.0
        drawline(ctx, [Point(50 + 0.95 * 45 * sin(radians),
            50 - 0.95 * 45 * cos(radians)),
            Point(50 + 1.0 * 45 * sin(radians),
            50 - 1.0 * 45 * cos(radians))], colorant"blue")
    end
end

Gtk.draw(can) do widget
    ctx = getgc(can)
    if length(drawcontext) < 1
        push!(drawcontext, ctx)
    else
        drawcontext[1] = ctx
    end
    clockbody(ctx)
end

function update(can)
    dtim = now()
    hr = hour(dtim)
    mi = minute(dtim)
    sec = second(dtim)
    if length(drawcontext) < 1
        return
    end
    ctx = drawcontext[1]
    clockbody(ctx)
    rad = (hr % 12) * pi / 6.0 + mi * pi / 360.0
    drawline(ctx, [Point(50, 50),
        Point(50 + 45 * 0.5 * sin(rad), 50 - 45 * 0.5 * cos(rad))], colorant"black")
    stroke(ctx)
    rad = mi * pi / 30.0  + sec * pi / 1800.0
    drawline(ctx, [Point(50, 50),
        Point(50 + 0.7 * 45 * sin(rad), 50 - 0.7 * 45 * cos(rad))], colorant"darkgreen")
    stroke(ctx)
    rad = sec * pi / 30.0
    drawline(ctx, [Point(50, 50),
        Point(50 + 0.9 * 45 * sin(rad), 50 - 0.9 * 45 * cos(rad))], colorant"red")
    stroke(ctx)
    reveal(can)
end

Gtk.showall(win)
sloc = Base.Threads.SpinLock()
lock(sloc)
signal_connect(win, :destroy) do widget
    unlock(sloc)
end
while !trylock(sloc)
    update(win)
    sleep(1.0)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
export Geometry, Circle, Rectangle, Triangle
export checkCollisionYZCircle
abstract type Geometry end
@mustimplement name(geometry::Geometry)

const geometryMinimum = 1.0Unitful.mm

struct Circle <: Geometry
  diameter::typeof(1.0Unitful.mm)
  name::String

  function Circle(;diameter::Unitful.Length, name::String)
    if diameter < geometryMinimum
      error("Circle is too small")
    else
      new(diameter, name)
    end
  end
end

name(geometry::Circle) = geometry.name
Circle(dict::Dict) = params_from_dict(Circle, dict)

struct Rectangle <: Geometry
  width::typeof(1.0Unitful.mm)
  height::typeof(1.0Unitful.mm)
  name::String

  function Rectangle(;width::Unitful.Length, height::Unitful.Length, name::String)
    if width < geometryMinimum || height < geometryMinimum
      error("Rectangle is too small")
    else
      new(width, height, name)
    end
  end
end
name(geometry::Rectangle) = geometry.name
Rectangle(dict::Dict) = params_from_dict(Rectangle, dict)

struct Triangle <: Geometry
  width::typeof(1.0Unitful.mm)
  height::typeof(1.0Unitful.mm)
  name::String

  function Triangle(;width::Unitful.Length, height::Unitful.Length, name::String)
    if width < geometryMinimum || height < geometryMinimum
      error("Triangle does not fit in scanner...")
    else
      new(width, height, name)
    end
  end
end
name(geometry::Triangle) = geometry.name
Triangle(dict::Dict) = params_from_dict(Triangle, dict)

Geometry(dict::Dict) = params_from_dict(eval(Symbol(pop!(dict,"type"))), dict)
convert(::Type{Geometry}, dict::Dict) = Geometry(dict)

"Returns true, 0mm, 0mm if given geometry does not collide with the scanner radius based on  given position and clearance"
@mustimplement checkCollisionYZCircle(geo::Geometry, scannerRad::Unitful.Length, posY::Unitful.Length, posZ::Unitful.Length, clearance::Unitful.Length)

function checkCollisionYZCircle(geo::Circle, scannerRad::Unitful.Length, posY::Unitful.Length, posZ::Unitful.Length, clearance::Unitful.Length)
  r = sqrt(posY^2 + posZ^2)

  if scannerRad - clearance >= r + geo.diameter / 2
    return true, zero(0.0u"mm"), zero(0.0u"mm")
  else
    delta_r = r - (scannerRad - clearance)
    delta_y1 = (delta_r+geo.diameter/2) * posY/r
    delta_z1 = (delta_r+geo.diameter/2) * posZ/r
    return false, 1.0001*delta_y1, 1.0001*delta_z1
  end
end

function checkCollisionYZCircle(geo::Rectangle, scannerRad::Unitful.Length, posY::Unitful.Length, posZ::Unitful.Length, clearance::Unitful.Length)
  corners = [ [posY - geo.width / 2 posZ - geo.height / 2]
              [posY - geo.width / 2 posZ + geo.height / 2]
              [posY + geo.width / 2 posZ + geo.height / 2]
              [posY + geo.width / 2 posZ - geo.height / 2] ];

  dist = norm.(eachrow(corners))
  if all(dist .< (scannerRad - clearance))
    return true, zero(0.0u"mm"), zero(0.0u"mm")
  else
    maxdist, maxidx = findmax(dist)
    delta = maxdist - (scannerRad - clearance)
    if abs(posY) < geo.width/2 + 1u"mm"
      return false, 0.0u"mm", (corners[maxidx,2] - sign(corners[maxidx,2])*sqrt((scannerRad-clearance)^2-(1.0001*corners[maxidx,1])^2))
    end
    if abs(posZ) < geo.height/2 + 1u"mm"
      return false, (corners[maxidx,1] - sign(corners[maxidx,1])*sqrt((scannerRad-clearance)^2-(1.0001*corners[maxidx,2])^2)), 0.0u"mm"
    end

    return false, 1.0001*delta * corners[maxidx,1] / maxdist, 1.0001*delta * corners[maxidx,2] / maxdist
  end
end

function checkCollisionYZCircle(geo::Triangle, scannerRad::Unitful.Length, posY::Unitful.Length, posZ::Unitful.Length, clearance::Unitful.Length)
  # following code only for symmetric triangles / center of gravity = center of plug adapter

  corners = [ [posY   posZ + 2/3 * geo.height]
              [posY-geo.width/2   posZ - 1/3 * geo.height]
              [posY+geo.width/2   posZ - 1/3 * geo.height] ];

  dist = norm.(eachrow(corners))
  if all(dist .< (scannerRad - clearance))
    return true, zero(0.0u"mm"), zero(0.0u"mm")
  else
    maxdist, maxidx = findmax(dist)
    delta = maxdist - (scannerRad - clearance)

    return false, 1.0001*delta * corners[maxidx,1] / maxdist, 1.0001*delta * corners[maxidx,2] / maxdist
  end
end


text:
function get_min_idx( Nmax::Int64, y )
    k = Nmax
    while abs(y[k-1]) < abs(y[k])
        k =  k - 1
        if k == 1 # last index
            break
        end
    end
    k = k - 1
    return k
end


text:
# This assumes that ps is normalized.
function rand(ps::Vector{Float64})
    acc = 0.
    r = rand()
    n = length(ps)
    for i in 1:n-1
        @inbounds acc += ps[i]
        isnan(acc) && nanerror()
        r < acc && return i
    end
    n == 1 && isnan(ps[1]) && nanerror()
    n
end

nanerror() = error("unexpected nan in rand")


text:
using PyPlot, Revise, HiQGA.transD_GP, Random, Statistics, BenchmarkTools
profile_within_juno  = false
## model
zfixed   = [-1e5,   0,    20,   50]
rho      = [1e12,   10,   1,   100]
nmax = 100
##  geometry and modeling parameters
rRx = 17.
zRx = -37.0
zTx = -35.0
rTx = 12.607
lowpassfcs = [300000, 450000.0]
ntimesperdecade = 10
nfreqsperdecade = 5
freqlow = 1e-3
# Note that the receiver depth needs to be in same model layer as transmitter.
## LM times and ramp
LM_times = [
1.539e-5  1.939e-5  2.439e-5  3.139e-5  3.939e-5  4.939e-5  6.239e-5  7.839e-5  9.939e-5  0.00012539  0.00015739  0.00019939  0.00025039  0.00031539  0.00039739  0.00050039  0.00063039  0.00079339
 1.9e-5    2.4e-5    3.1e-5    3.9e-5    4.9e-5    6.2e-5    7.8e-5    9.9e-5    0.000125  0.000157    0.000199    0.00025     0.000315    0.000397    0.0005      0.00063     0.000793    0.000999]'
LM_times = exp.(mean(log.(LM_times),dims=2))[:]
LM_ramp =[
-0.001  -0.0009146  -0.0007879  -0.0005964  0.0  4.629e-7  8.751e-7  1.354e-6  2.54e-6  3.972e-6  5.404e-6  5.721e-6  6.113e-6  6.663e-6   8.068e-6  0.00125
 0.0     0.6264      0.9132      0.9905     1.0  0.9891    0.9426    0.8545    0.6053   0.303     0.04077   0.01632   0.004419  0.0006323  0.0       0.0]'
## HM times and ramp
HM_times = [
7.539e-5  9.639e-5  0.00012239  0.00015439  0.00019639  0.00024739  0.00031239  0.00039439  0.00049739  0.00062739  0.00079039  0.00099639  0.00125539  0.00158139  0.00199139  0.00250839  0.00315839  0.00397739  0.00500839  0.00630639  0.00793939
9.6e-5    0.000122  0.000154    0.000196    0.000247    0.000312    0.000394    0.000497    0.000627    0.00079     0.000996    0.001255    0.001581    0.001991    0.002508    0.003158    0.003977    0.005008    0.006306    0.007939    0.009739]'
HM_times = exp.(mean(log.(HM_times),dims=2))[:]
HM_ramp = [
-0.01  -0.008386  -0.00638  -0.003783  0.0  3.96e-7  7.782e-7  1.212e-6  3.44e-6  1.981e-5  3.619e-5  3.664e-5  3.719e-5   3.798e-5  3.997e-5  0.01
  0.0    0.4568     0.7526    0.9204    1.0  0.9984   0.9914    0.9799    0.9175   0.4587    0.007675  0.003072  0.0008319  0.000119  0.0       0.0]'
## LM operator
Flm = transD_GP.AEM_VMD_HMD.HFieldDHT(
                      ntimesperdecade = ntimesperdecade,
                      nfreqsperdecade = nfreqsperdecade,
                      lowpassfcs = lowpassfcs,
                      times  = LM_times,
                      ramp   = LM_ramp,
                      nmax   = nmax,
                      zTx    = zTx,
                      rRx    = rRx,
                      rTx    = rTx,
                      zRx    = zRx,
		      freqlow = freqlow
		      )
## HM operator
Fhm = transD_GP.AEM_VMD_HMD.HFieldDHT(
                      ntimesperdecade = ntimesperdecade,
                      nfreqsperdecade = nfreqsperdecade,
                      lowpassfcs = lowpassfcs,
                      times  = HM_times,
                      ramp   = HM_ramp,
                      nmax   = nmax,
                      zTx    = zTx,
                      rRx    = rRx,
                      rTx    = rTx,
                      zRx    = zRx,
		      freqlow = freqlow)
## get the fields in place
μ = transD_GP.AEM_VMD_HMD.μ
@time transD_GP.AEM_VMD_HMD.getfieldTD!(Flm, zfixed, rho)
@time transD_GP.AEM_VMD_HMD.getfieldTD!(Fhm, zfixed, rho)
##
figure()
loglog(Flm.times,μ*abs.(Flm.dBzdt), label="lm")
loglog(Fhm.times,μ*abs.(Fhm.dBzdt), label="hm")
grid()
## timing FD
ntimes = 1000
t = time()
for i = 1:ntimes
    transD_GP.AEM_VMD_HMD.getfieldFD!(Flm, zfixed, rho)
    transD_GP.AEM_VMD_HMD.getfieldFD!(Fhm, zfixed, rho)
end
t = time() - t
@info "FD timing is $(t/ntimes) s"
## timing TD
ntimes = 1000
t = time()
for i = 1:ntimes
    transD_GP.AEM_VMD_HMD.getfieldTD!(Flm, zfixed, rho)
    transD_GP.AEM_VMD_HMD.getfieldTD!(Fhm, zfixed, rho)
end
t = time() - t
@info "TD timing is $(t/ntimes) s for a $(length(zfixed)) layer model"
## time many random layers
Random.seed!(435)
nlayers = 50
z = [-1e5, 0, 20, 50 .+ cumsum(15*rand(nlayers-3))...]
ρ = [1e13, 10, 1, 10 .^(-0.5 .+ 1.5*rand(nlayers-3))...]
figure()
step(log10.(ρ[2:end]), z[2:end])
xlim(-0.5,1)
grid(); gca().invert_yaxis()
ntimes = 100
## solely timing
function mytest(Fl, Fh, zz, rr)
   transD_GP.AEM_VMD_HMD.getfieldTD!(Fl, zz, rr)
   transD_GP.AEM_VMD_HMD.getfieldTD!(Fh, zz, rr)
end
@info "timing for a $(length(z)) layer model is:"
@btime mytest($Flm, $Fhm, $z, $ρ)
## using Profile
if profile_within_juno
    using Profile
    function doprofile(ntimes::Int)
        for i = 1:ntimes
            transD_GP.AEM_VMD_HMD.getfieldTD!(Flm, z, ρ)
        end
    end
    @profiler doprofile(ntimes) combine=true
end

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA

max_forks_repo_name: GeoscienceAustralia/HiQGA


text:
# FusionTreeIterator:
# iterate over fusion trees for fixed coupled and uncoupled sector labels
#==============================================================================#
function fusiontrees(uncoupled::NTuple{N, I}, coupled::I = one(I),
                        isdual::NTuple{N, Bool} = ntuple(n->false, Val(N))) where
                        {N, I<:Sector}
    FusionTreeIterator{I, N}(uncoupled, coupled, isdual)
end

struct FusionTreeIterator{I<:Sector, N}
    uncoupled::NTuple{N, I}
    coupled::I
    isdual::NTuple{N, Bool}
end

Base.IteratorSize(::FusionTreeIterator) = Base.HasLength()
Base.IteratorEltype(::FusionTreeIterator) = Base.HasEltype()
Base.eltype(T::Type{FusionTreeIterator{I, N}}) where {I<:Sector, N} =
    fusiontreetype(I, N)

Base.length(iter::FusionTreeIterator) = _fusiondim(iter.uncoupled, iter.coupled)
_fusiondim(u::Tuple{}, c::I) where {I<:Sector} = Int(one(c) == c)
_fusiondim(u::Tuple{I}, c::I) where {I<:Sector} = Int(u[1] == c)
_fusiondim((a, b)::Tuple{I, I}, c::I) where {I<:Sector} = Int(Nsymbol(a, b, c))
function _fusiondim(u::Tuple{I, I, Vararg{I}}, c::I) where {I<:Sector}
    a = u[1]
    b = u[2]
    d = 0
    for c′ in a ⊗ b
        d += Nsymbol(a, b, c′)*_fusiondim((c′, TupleTools.tail2(u)...), c)
    end
    return d
end

# * Iterator methods:
#   Start with special cases:
function Base.iterate(it::FusionTreeIterator{I, 0},
                        state = (it.coupled != one(I))) where {I<:Sector}
    state && return nothing
    T = vertex_labeltype(I)
    tree = FusionTree{I, 0, 0, 0, T}((), one(I), (), (), ())
    return tree, true
end

function Base.iterate(it::FusionTreeIterator{I, 1},
                        state = (it.uncoupled[1] != it.coupled)) where {I<:Sector}
    state && return nothing
    T = vertex_labeltype(I)
    tree = FusionTree{I, 1, 0, 0, T}(it.uncoupled, it.coupled, it.isdual, (), ())
    return tree, true
end

#   General case:
function Base.iterate(it::FusionTreeIterator{I, N} where {N}) where {I<:Sector}
    next = _iterate(it.uncoupled, it.coupled)
    next === nothing && return nothing
    lines, vertices, states = next
    vertexlabels = labelvertices(it.uncoupled, it.coupled, lines, vertices)
    f = FusionTree(it.uncoupled, it.coupled, it.isdual, lines, vertexlabels)
    return f, (lines, vertices, states)
end
function Base.iterate(it::FusionTreeIterator{I, N} where {N}, state) where {I<:Sector}
    next = _iterate(it.uncoupled, it.coupled, state...)
    next === nothing && return nothing
    lines, vertices, states = next
    vertexlabels = labelvertices(it.uncoupled, it.coupled, lines, vertices)
    f = FusionTree(it.uncoupled, it.coupled, it.isdual, lines, vertexlabels)
    return f, (lines, vertices, states)
end

labelvertices(uncoupled::NTuple{2, I}, coupled::I, lines::Tuple{},
                vertices::Tuple{Int}) where {I<:Sector} =
    (vertex_ind2label(vertices[1], uncoupled..., coupled),)

function labelvertices(uncoupled::NTuple{N, I}, coupled::I, lines,
                        vertices) where {I<:Sector, N}
    c = lines[1]
    resttree = tuple(c, TupleTools.tail2(uncoupled)...)
    rest = labelvertices(resttree, coupled, tail(lines), tail(vertices))
    l = vertex_ind2label(vertices[1], uncoupled[1], uncoupled[2], c)
    return (l, rest...)
end

# Actual implementation
@inline function _iterate(uncoupled::NTuple{2, I}, coupled::I, lines = (),
                            vertices = (0,), states = ()) where {I<:Sector}
    a, b = uncoupled
    n = vertices[1] + 1
    n > Nsymbol(a, b, coupled) && return nothing
    return (), (n,), ()
end

function _iterate(uncoupled::NTuple{N, I}, coupled::I) where {N, I<:Sector}
    a, b, = uncoupled
    it = a ⊗ b
    next = iterate(it)
    next === nothing && return nothing
    # this should not happen: there should always be at least one fusion output
    c, s = next
    resttree = tuple(c, TupleTools.tail2(uncoupled)...)
    rest = _iterate(resttree, coupled)
    while rest === nothing
        next = iterate(it, s)
        next === nothing && return nothing
        c, s = next
        resttree = tuple(c, TupleTools.tail2(uncoupled)...)
        rest = _iterate(resttree, coupled)
    end
    n = 1
    restlines, restvertices, reststates = rest
    lines = (c, restlines...)
    vertices = (n, restvertices...)
    states = (s, reststates...)
    return lines, vertices, states
end

function _iterate(uncoupled::NTuple{N, I}, coupled::I, lines, vertices,
                            states) where {N, I<:Sector}
    a, b, = uncoupled
    it = a ⊗ b
    c = lines[1]
    n = vertices[1]
    s = states[1]
    restlines = tail(lines)
    restvertices = tail(vertices)
    reststates = tail(states)
    if n < Nsymbol(a, b, c)
        n += 1
        return lines, (n, restvertices...), states
    end
    n = 1
    resttree = tuple(c, TupleTools.tail2(uncoupled)...)
    rest = _iterate(resttree, coupled, restlines, restvertices, reststates)
    while rest === nothing
        next = iterate(it, s)
        next === nothing && return nothing
        c, s = next
        resttree = tuple(c, TupleTools.tail2(uncoupled)...)
        rest = _iterate(resttree, coupled)
    end
    restlines, restvertices, reststate = rest
    lines = (c, restlines...)
    vertices = (n, restvertices...)
    states = (s, reststate...)
    return lines, vertices, states
end


text:
using PrettyTables
export writetdf, readtdf

function writetdf(filename, df)
    out = DataFrame()
    align = Symbol[]
    @assert nrow(df) > 1
    for name in Symbol.(names(df))
        col = deepcopy(df[:, name])
        t0 = eltype(col)
        tt = nonmissingtype(t0)
        ii = findall(.!ismissing.(col))
        if tt <: AbstractString
            col[ii] .= String.(strip.(df[ii, name]))
            push!(align, :l)
        elseif (tt == Bool)  ||
            (tt == UnitRange{Int64})  ||
            (tt == StepRange{Int64, Int64})
            ;
            push!(align, :c)
        elseif tt <: Integer
            col[ii] .= Int64.(df[ii, name])
            push!(align, :r)
        elseif tt <: AbstractFloat
            col[ii] .= Float64.(df[ii, name])
            push!(align, :r)
        else
            error("Unsupported data type: $(t0)")
        end
        out[!, name] = col
    end

    f = open(filename, "w")
    pretty_table(f, out, tf=tf_ascii_rounded, alignment=align)
    close(f)
end

function readtdf(filename)
    lines = readlines(filename)
    names = Symbol.(strip.(split(lines[2], '|'))); names = names[2:end-1]
    types = String.(strip.(split(lines[3], '|'))); types = types[2:end-1]
    out = DataFrame()
    for icol in 1:length(names)
        tt = types[icol]
        if tt[end] == '?'
            out[!, names[icol]] = eval(Meta.parse("Union{Missing, " * tt[1:end-1] * "}[]"))
        else
            out[!, names[icol]] = eval(Meta.parse(tt * "[]"))
        end
    end
    for i in 5:length(lines)-1
        row = []
        line = String.(strip.(split(lines[i], '|')));
        line = line[2:end-1]

        for icol in 1:length(names)
            t0 = types[icol]
            tt = types[icol]
            hasmissing = false
            if tt[end] == '?'
                hasmissing = true
                tt = tt[1:end-1]
            end
            vv = line[icol]

            if hasmissing  &&  ((vv == "missing")  ||  (vv == ""))
                push!(row, missing)
            elseif tt == "String"
                push!(row, vv)
            elseif tt == "Bool"
                if (vv == "true")  ||  ((length(vv) != 0)  &&  (vv != "false"))
                    push!(row, true)
                else
                    push!(row, false)
                end
            elseif (tt == "UnitRange{Int64}")  ||
                (tt == "StepRange{Int64, Int64}")
                push!(row, eval(Meta.parse(vv)))
            elseif tt == "Int64"
                push!(row, Meta.parse(vv))
            elseif tt == "Float64"
                if vv == "NaN"
                    push!(row, NaN)
                else
                    push!(row, Meta.parse(vv))
                end
            else
                error("Unsupported data type: $(t0)")
            end
        end
        push!(out, row)
    end
    return out
end


text:
using CUDA

include("LineOfResponse.jl")
include("Util.jl")

function ray_tracing(event::Event, t, DIMX, DIMY, DIMZ)
    x1 = event.lor.P1.x
    x2 = event.lor.P2.x
    y1 = event.lor.P1.y
    y2 = event.lor.P2.y
    z1 = event.lor.P1.z
    z2 = event.lor.P2.z

    k_y = (y2 - y1) / (x2 - x1)
    k_z = (z2 - z1) / (x2 - x1)

    i = t - DIMX/2 - 1

    if k_y >= 0
        y_min = y1 + k_y*(i-x1) + DIMY/2
        y_plus = y1 + k_y*(i+1-x1) + DIMY/2
    else
        y_min = y1 + k_y*(i+1-x1) + DIMY/2
        y_plus = y1 + k_y*(i-x1) + DIMY/2
    end

    Y_min = floor(Int32, y_min) # Zero based indices
    Y_plus = floor(Int32, y_plus) # Zero based indices

    if k_z >= 0
        z_min = z1 + k_z*(i-x1) + DIMZ/2
        z_plus = z1 + k_z*(i+1-x1) + DIMZ/2
    else
        z_min = z1 + k_z*(i+1-x1) + DIMZ/2
        z_plus = z1 + k_z*(i-x1) + DIMZ/2
    end

    Z_min = floor(Int32, z_min) # Zero based indices
    Z_plus = floor(Int32, z_plus) # Zero based indices

    l = sqrt(1 + (y_plus - y_min)^2 + (z_plus - z_min)^2)

    l_min_min = 0.0
    l_min_plus = 0.0
    l_plus_min = 0.0
    l_plus_plus = 0.0

    if Y_min == Y_plus
        if Y_min >= 0 && Y_min < DIMY
            if Z_min == Z_plus
                # 1-voxel case
                if Z_min >= 0 && Z_min < DIMZ
                    l_min_min = l
                end
            else
                # 2-voxel case in z-direction
                if Z_min >= -1 && Z_min < DIMZ
                    r_z = (Z_plus - z_min)/(z_plus - z_min)

                    if Z_min > -1
                        l_min_min = l * r_z
                    end
                    if Z_plus < DIMZ
                        l_min_plus = l * (1 - r_z)
                    end
                end
            end
        end
    else
        if Y_min >= -1 && Y_min < DIMY
            if Z_min == Z_plus
                # 2-voxel case in y-direction
                if Z_min >= 0 && Z_min < DIMZ
                    r_y = (Y_plus - y_min)/(y_plus - y_min)

                    if Y_min > -1
                        l_min_min = l * r_y
                    end
                    if Y_plus < DIMY
                        l_plus_min = l * (1 - r_y)
                    end
                end
            else
                # 3-voxel cases
                if Z_min >= -1 && Z_min < DIMZ
                    r_y = (Y_plus - y_min)/(y_plus - y_min)
                    r_z = (Z_plus - z_min)/(z_plus - z_min)

                    if r_y > r_z
                        if Y_min > -1 && Z_min > -1
                            l_min_min = l * r_z
                        end
                        if Y_min > -1 && Z_plus < DIMZ
                            l_min_plus = l * (r_y - r_z)
                        end
                        if Y_plus < DIMY && Z_plus < DIMZ
                            l_plus_plus = l * (1 - r_y)
                        end
                    else
                        if Y_min > -1 && Z_min > -1
                            l_min_min = l * r_y
                        end
                        if Y_plus < DIMY && Z_min > -1
                            l_plus_min = l * (r_z - r_y)
                        end
                        if Y_plus < DIMY && Z_plus < DIMZ
                            l_plus_plus = l * (1 - r_z)
                        end
                    end
                end
            end
        end
    end

    return Y_min, Y_plus, Z_min, Z_plus, l_min_min, l_min_plus, l_plus_min, l_plus_plus
end

function perm(main_plane::Plane, x, y, z)
	if main_plane == x_plane
		return x, y, z
	elseif main_plane == y_plane
		return y, x, z
	else
		return z, y, x
	end
end

function gpu_kernel(events::CuDeviceArray{Event}, image::CuDeviceArray{T,3}, corr::CuDeviceArray{T,3}, tmp_total_values::CuDeviceArray{T,1}, DIMX, DIMY, DIMZ) where {T}
    t = threadIdx().x
    index = blockIdx().x

    event = events[index]

    # Ray tracing
    Y_min, Y_plus, Z_min, Z_plus, l_min_min, l_min_plus, l_plus_min, l_plus_plus = ray_tracing(event, t, DIMX, DIMY, DIMZ)

    if Y_min >= 0 && Z_min >=0 && Y_plus < DIMY && Z_plus < DIMZ
        value = 0.0
		value += image[perm(event.mainPlane, t, Y_min+1, Z_min+1)...] * l_min_min
		value += image[perm(event.mainPlane, t, Y_min+1, Z_plus+1)...] * l_min_plus
		value += image[perm(event.mainPlane, t, Y_plus+1, Z_min+1)...] * l_plus_min
		value += image[perm(event.mainPlane, t, Y_plus+1, Z_plus+1)...] * l_plus_plus

        # Forward project
        CUDA.atomic_add!(pointer(tmp_total_values, index), Float32(value))
        CUDA.sync_threads()
        total_value = tmp_total_values[index]

        if total_value > 0
            # Compare and back project
			corr[perm(event.mainPlane, t, Y_min+1, Z_min+1)...] += l_min_min / total_value
            corr[perm(event.mainPlane, t, Y_min+1, Z_plus+1)...] += l_min_plus / total_value
            corr[perm(event.mainPlane, t, Y_plus+1, Z_min+1)...] += l_plus_min / total_value
            corr[perm(event.mainPlane, t, Y_plus+1, Z_plus+1)...] += l_plus_plus / total_value
        end
    end

    return nothing
end

function reconstruct3D(events, DIMX, DIMY, DIMZ, recon_iters)
  # Transfer the events from CPU memory to GPU memory
  c_events = CUDA.CuArray(events);

  # Initialise the image array on GPU
  c_image = CUDA.ones(Float32, DIMX, DIMY, DIMZ);

  for k = 1:recon_iters
    # Initialise correction matrix memory on GPU
    c_corr = CUDA.zeros(Float32, DIMX, DIMY, DIMZ);
    # Initialise array to store forward projection values on GPU
    c_tmp_total_values = CUDA.zeros(Float32, length(c_events))

    # Schedule a kernel that performs:
    #   - Forward projection
    #   - Compare
    #   - Back projection
    @cuda blocks=length(c_events) threads=DIMX gpu_kernel(c_events, c_image, c_corr, c_tmp_total_values, DIMX, DIMY, DIMZ)

    # Perform the update step
    c_image = c_image .* c_corr
  end

  # Transfer the image estimate from GPU to CPU memory
  image = Array(c_image)

  return image
end

max_stars_repo_name: michielvangendt/PET-Julia

max_issues_repo_name: michielvangendt/PET-Julia

max_forks_repo_name: michielvangendt/PET-Julia


text:
################################################################################
### A polar domain described by a single variable function
################################################################################
export PolarDomain
struct PolarDomain{T} <: EuclideanDomain{2,T}
    charFun
    box     ::  EuclideanDomain{2,T}
end
export polardomain
function polardomain(char, domain)
    T = subeltype(domain)
    # F = Fun(char, Fourier(100, -T(pi), T(pi)), Interval(-T(pi), T(pi)))
    # PolarDomain(F, boundingbox(domain))
    PolarDomain(char, boundingbox(domain))
end

indomain(x, c::PolarDomain) = sqrt(x[1]^2+x[2]^2) < real(c.charFun(atan(x[2],x[1])))

boundingbox(c::PolarDomain) = c.box

show(io::IO, c::PolarDomain) = print(io, "a domain in polar coordinates")

distance(x, c::PolarDomain) = real(c.charFun(atan(x[2],x[1])))-norm(x)

function normal(x, c::PolarDomain)
    phi = atan(x[2], x[1])
    y=c.charFun'(phi)*sin(phi)+c.charFun(phi)*cos(phi)
    x=c.charFun'(phi)*cos(phi)-c.charFun(phi)*sin(phi)
    return [real(y/(norm([x;y]))),-real(x/(norm([x;y])))]
end


text:
module ParallelRun

import FITSIO
import JLD

import ..Log
using ..Model
import ..SDSSIO
import ..Infer
import ..SDSSIO: RunCamcolField
import ..PSF

import ..DeterministicVI: infer_source

include("joint_infer.jl")

#set this to false to use source-division parallelism
const SKY_DIVISION_PARALLELISM=true

const MIN_FLUX = 2.0

# In production mode, rather the development mode, always catch exceptions
const is_production_run = haskey(ENV, "CELESTE_PROD") &&
                          ENV["CELESTE_PROD"] != ""

# Use distributed parallelism (with Dtree)
if haskey(ENV, "USE_DTREE") && ENV["USE_DTREE"] != ""
    const Distributed = true
    using Dtree
else
    const Distributed = false
    const dt_nodeid = 1
    const dt_nnodes = 1
    DtreeScheduler(n, f) = ()
    initwork(dt) = 0, (1, 0)
    getwork(dt) = 0, (1, 0)
    runtree(dt) = 0
    cpu_pause() = ()
end

# Use threads (on the loop over sources)
const Threaded = true
if Threaded && VERSION > v"0.5.0-dev"
    using Base.Threads
else
    # Pre-Julia 0.5 there are no threads
    nthreads() = 1
    threadid() = 1
    macro threads(x)
        x
    end
    SpinLock() = 1
    lock(l) = ()
    unlock(l) = ()
end

# A workitem is of this ra / dec size
const wira = 0.025
const widec = 0.025

"""
Timing information.
"""
type InferTiming
    query_fids::Float64
    num_infers::Int64
    read_photoobj::Float64
    read_img::Float64
    fit_psf::Float64
    opt_srcs::Float64
    num_srcs::Int64
    write_results::Float64
    wait_done::Float64

    InferTiming() = new(0.0, 0, 0.0, 0.0, 0.0, 0.0, 0, 0.0, 0.0)
end

function add_timing!(i::InferTiming, j::InferTiming)
    i.query_fids = i.query_fids + j.query_fids
    i.num_infers = i.num_infers + j.num_infers
    i.read_photoobj = i.read_photoobj + j.read_photoobj
    i.read_img = i.read_img + j.read_img
    i.fit_psf = i.fit_psf + j.fit_psf
    i.opt_srcs = i.opt_srcs + j.opt_srcs
    i.num_srcs = i.num_srcs + j.num_srcs
    i.write_results = i.write_results + j.write_results
    i.wait_done = i.wait_done + j.wait_done
end

immutable BoundingBox
    ramin::Float64
    ramax::Float64
    decmin::Float64
    decmax::Float64
end

function BoundingBox(ramin::String, ramax::String, decmin::String, decmax::String)
    BoundingBox(parse(Float64, ramin),
                parse(Float64, ramax),
                parse(Float64, decmin),
                parse(Float64, decmax))
end

@inline nputs(nid, s) = ccall(:puts, Cint, (Ptr{Int8},), string("[$nid] ", s))
@inline phalse(b) = b[] = false

include("source_division_inference.jl")

"""
Given a BoundingBox that is to be divided into `nra` x `ndec` subareas,
return the `i`th subarea. `i` is a linear index between 1 and
`nra * ndec`.

This function assumes a cartesian (rather than spherical) coordinate system!
"""
function divide_skyarea(box, nra, ndec, i)
    global wira, widec
    ix, iy = ind2sub((nra, ndec), i)

    return (box.ramin + (ix - 1) * wira,
            min(box.ramin + ix * wira, box.ramax),
            box.decmin + (iy - 1) * widec,
            min(box.decmin + iy * widec, box.decmax))
end

function time_puts(elapsedtime, bytes, gctime, allocs)
    s = @sprintf("%10.6f seconds", elapsedtime/1e9)
    if bytes != 0 || allocs != 0
        bytes, mb = Base.prettyprint_getunits(bytes, length(Base._mem_units),
                            Int64(1024))
        allocs, ma = Base.prettyprint_getunits(allocs, length(Base._cnt_units),
                            Int64(1000))
        if ma == 1
            s = string(s, @sprintf(" (%d%s allocation%s: ", allocs,
                                   Base._cnt_units[ma], allocs==1 ? "" : "s"))
        else
            s = string(s, @sprintf(" (%.2f%s allocations: ", allocs,
                                   Base._cnt_units[ma]))
        end
        if mb == 1
            s = string(s, @sprintf("%d %s%s", bytes,
                                   Base._mem_units[mb], bytes==1 ? "" : "s"))
        else
            s = string(s, @sprintf("%.3f %s", bytes, Base._mem_units[mb]))
        end
        if gctime > 0
            s = string(s, @sprintf(", %.2f%% gc time", 100*gctime/elapsedtime))
        end
        s = string(s, ")")
    elseif gctime > 0
        s = string(s, @sprintf(", %.2f%% gc time", 100*gctime/elapsedtime))
    end
    nputs(dt_nodeid, s)
end

"""
Divide `N` into `np` parts as evenly as possible, returning `my` part as
a (first, last) tuple.
"""
function divparts(N, np, my)
    len, rem = divrem(N, np)
    if len == 0
        if my > rem
            return 1, 0
        end
        len, rem = 1, 0
    end
    # compute my part
    f = 1 + ((my-1) * len)
    l = f + len - 1
    # distribute remaining evenly
    if rem > 0
        if my <= rem
            f = f + (my-1)
            l = l + my
        else
            f = f + rem
            l = l + rem
        end
    end
    return f, l
end

"""
Divide the given ra, dec range into sky areas of `wira`x`widec` and
use Dtree to distribute these sky areas to nodes. Within each node
use `one_node_infer()` to fit the Celeste model to sources in each sky area.
"""
function divide_sky_and_infer(
                box::BoundingBox,
                stagedir::String;
                timing=InferTiming(),
                outdir=".")
    if dt_nodeid == 1
        nputs(dt_nodeid, "running on $dt_nnodes nodes")
    end

    # how many `wira` X `widec` sky areas (work items)?
    global wira, widec
    nra = ceil(Int64, (box.ramax - box.ramin) / wira)
    ndec = ceil(Int64, (box.decmax - box.decmin) / widec)

    num_work_items = nra * ndec
    each = ceil(Int64, num_work_items / dt_nnodes)

    if dt_nodeid == 1
        nputs(dt_nodeid, "work item dimensions: $wira X $widec")
        nputs(dt_nodeid, "$num_work_items work items, ~$each per node")
    end

    # create Dtree and get the initial allocation
    dt, is_parent = DtreeScheduler(num_work_items, 0.4)
    ni, (ci, li) = initwork(dt)
    rundt = Ref(runtree(dt))
    @inline function rundtree(again)
        if again[]
            again[] = runtree(dt)
            cpu_pause()
        end
        again[]
    end

    # work item processing loop
    nputs(dt_nodeid, "initially $ni work items ($ci to $li)")
    itimes = InferTiming()
    while ni > 0
        li == 0 && break
        if ci > li
            nputs(dt_nodeid, "consumed allocation (last was $li)")
            ni, (ci, li) = getwork(dt)
            nputs(dt_nodeid, "got $ni work items ($ci to $li)")
            continue
        end
        item = ci
        ci = ci + 1

        # map item to subarea
        iramin, iramax, idecmin, idecmax = divide_skyarea(box, nra, ndec, item)

        # Get vector of (run, camcol, field) triplets overlapping this patch
        tic()
        box = BoundingBox(iramin, iramax, idecmin, idecmax)
        rcfs = get_overlapping_fields(box, stagedir)
        itimes.query_fids = toq()

        # run inference for this subarea
        results = one_node_infer(rcfs, stagedir;
                                 box=BoundingBox(iramin, iramax, idecmin, idecmax),
                                 reserve_thread=rundt,
                                 thread_fun=rundtree,
                                 timing=itimes)
        tic()
        save_results(outdir, iramin, iramax, idecmin, idecmax, results)
        itimes.write_results = toq()

        timing.num_infers = timing.num_infers+1
        add_timing!(timing, itimes)
        rundtree(rundt)
    end
    nputs(dt_nodeid, "out of work")
    tic()
    while rundt[]
        rundtree(rundt)
    end
    finalize(dt)
    timing.wait_done = toq()
end

function infer_init(rcfs::Vector{RunCamcolField},
                    stagedir::String;
                    objid="",
                    box=BoundingBox(-1000., 1000., -1000., 1000.),
                    target_rcfs=RunCamcolField[],
                    primary_initialization=true)
    # Read all primary objects in these fields.
    duplicate_policy = primary_initialization ? :primary : :first
    catalog = SDSSIO.read_photoobj_files(rcfs, stagedir,
                        duplicate_policy=duplicate_policy)
    Log.info("$(length(catalog)) primary sources")

    # Filter out low-flux objects in the catalog.
    catalog = filter(entry->(maximum(entry.star_fluxes) >= MIN_FLUX), catalog)
    Log.info("$(length(catalog)) primary sources after MIN_FLUX cut")

    # Get indicies of entries in the RA/Dec range of interest.
    entry_in_range = entry->((box.ramin < entry.pos[1] < box.ramax) &&
                             (box.decmin < entry.pos[2] < box.decmax))
    target_sources = find(entry_in_range, catalog)

    Log.info(string("Found $(length(target_sources)) target sources in ",
          "$(box.ramin), $(box.ramax), $(box.decmin), $(box.decmax)"))

    # For infer-box.jl, target sources are everything in the box.
    # For infer-rcf.jl, target sources are primary detections in the target rcf.
    if !isempty(target_rcfs)
        target_entries = SDSSIO.read_photoobj_files(target_rcfs,
                                                    stagedir,
                                                    duplicate_policy=:primary)
        target_ids = Set([entry.objid for entry in target_entries])
        target_sources = filter(ts->(catalog[ts].objid in target_ids), target_sources)
        Log.info("$(length(target_sources)) target light sources after target_rcf cut")
    end

    # Filter any object not specified, if an objid is specified
    if objid != ""
        target_sources = filter(ts->(catalog[ts].objid == objid), target_sources)
        Log.info("$(length(target_sources)) target light sources after objid cut")
    end

    return catalog, target_sources
end

"""
Use mulitple threads on one node to fit the Celeste model to sources in a given
bounding box.
"""
function one_node_infer(rcfs::Vector{RunCamcolField},
                        stagedir::String;
                        infer_callback=one_node_single_infer,
                        objid="",
                        box=BoundingBox(-1000., 1000., -1000., 1000.),
                        target_rcfs=RunCamcolField[],
                        primary_initialization=true)
    catalog, target_sources = infer_init(rcfs,
                                 stagedir;
                                 objid=objid,
                                 box=box,
                                 target_rcfs=target_rcfs,
                                 primary_initialization=primary_initialization)

    # If there are no objects of interest, return early.
    if length(target_sources) == 0
        images = Image[]
        neighbor_map = Vector{Int}[]
    else
        # Read in images for all (run, camcol, field).
        images = SDSSIO.load_field_images(rcfs, stagedir)

        tic()
        neighbor_map = Infer.find_neighbors(target_sources, catalog, images)
        Log.info("neighbors found in $(toq()) seconds")
    end

    Log.info("Running with $(nthreads()) threads")

    # NB: All I/O happens above. The methods below don't touch disk.
    infer_callback(catalog, target_sources, neighbor_map, images)
end

immutable OptimizedSource
    thingid::Int64
    objid::String
    init_ra::Float64
    init_dec::Float64
    vs::Vector{Float64}
end

function one_node_single_infer(catalog::Vector{CatalogEntry},
                               target_sources::Vector{Int},
                               neighbor_map::Vector{Vector{Int}},
                               images::Vector{Image};
                               infer_source_callback=infer_source,
                               reserve_thread=Ref(false),
                               thread_fun=phalse,
                               timing=InferTiming())
    curr_source = 1
    last_source = length(target_sources)
    sources_lock = SpinLock()
    results = OptimizedSource[]
    results_lock = SpinLock()

    # iterate over sources
    function process_sources()
        tid = threadid()

        if reserve_thread[] && tid == 1
            while reserve_thread[]
                thread_fun(reserve_thread)
                cpu_pause()
            end
        else
            while true
                lock(sources_lock)
                ts = curr_source
                curr_source += 1
                unlock(sources_lock)
                if ts > last_source
                    break
                end

                try
                    s = target_sources[ts]
                    entry = catalog[s]

                    # could subset images to images_local here too.
                    neighbors = catalog[neighbor_map[ts]]

                    t0 = time()
                    vs_opt = infer_source_callback(images, neighbors, entry)
                    runtime = time() - t0

                    result = OptimizedSource(entry.thing_id,
                                             entry.objid,
                                             entry.pos[1],
                                             entry.pos[2],
                                             vs_opt)
                    lock(results_lock)
                    push!(results, result)
                    unlock(results_lock)

                    rt1 = round(runtime, 1)
                    Log.debug("source $s [objid $(entry.objid)] took $rt1 seconds")
                catch ex
                    if is_production_run || nthreads() > 1
                        Log.error(string(ex))
                    else
                        rethrow(ex)
                    end
                end
            end
        end
    end

    tic()
    if nthreads() == 1
        process_sources()
    else
        ccall(:jl_threading_run, Void, (Any,), Core.svec(process_sources))
        ccall(:jl_threading_profile, Void, ())
    end
    timing.opt_srcs = toq()
    timing.num_srcs = length(target_sources)

    results
end

"""
Query the SDSS database for all fields that overlap the given RA, Dec range.
"""
function get_overlapping_field_extents(query::BoundingBox, stagedir::String)
    f = FITSIO.FITS("$stagedir/field_extents.fits")

    hdu = f[2]::FITSIO.TableHDU

    # read in the entire table.
    all_run = read(hdu, "run")::Vector{Int16}
    all_camcol = read(hdu, "camcol")::Vector{UInt8}
    all_field = read(hdu, "field")::Vector{Int16}
    all_ramin = read(hdu, "ramin")::Vector{Float64}
    all_ramax = read(hdu, "ramax")::Vector{Float64}
    all_decmin = read(hdu, "decmin")::Vector{Float64}
    all_decmax = read(hdu, "decmax")::Vector{Float64}

    close(f)

    ret = Tuple{RunCamcolField, BoundingBox}[]

    # The ramin, ramax, etc is a bit unintuitive because we're looking
    # for any overlap.
    for i in eachindex(all_ramin)
        if (all_ramax[i] > query.ramin && all_ramin[i] < query.ramax &&
                all_decmax[i] > query.decmin && all_decmin[i] < query.decmax)
            cur_box = BoundingBox(all_ramin[i], all_ramax[i],
                                  all_decmin[i], all_decmax[i])
            cur_fe = RunCamcolField(all_run[i], all_camcol[i], all_field[i])
            push!(ret, (cur_fe, cur_box))
        end
    end

    return ret
end

"""
Like `get_overlapping_fields`, but return a Vector of
(run, camcol, field) triplets.
"""
function get_overlapping_fields(query::BoundingBox, stagedir::String)
    fes = get_overlapping_field_extents(query, stagedir)
    [fe[1] for fe in fes]
end

immutable UnknownRCF <: Exception
    rcf::RunCamcolField
end

"""
This function is called directly from the `bin` directory.
It optimizes all the primary detections in the specified run-camcol-field,
using essentially all relevant images. Typically relevant image data includes
run-camcol-field's in addition to the specified one, that overlap
with the specified run-camcol-field. I say "essentially" because
some large light sources may contribute photons to images that do
not overlap with the rcf containig the center of the light source's
primary detection, and these images are excluded. This is more of a
feature than a bug because we need to limit the amount of computation
per light souce.
"""
function infer_rcf(rcf::RunCamcolField,
                   stagedir::String,
                   outdir::String;
                   objid="")
    whole_sky = BoundingBox(-999, 999, -999, 999)
    all_fes = get_overlapping_field_extents(whole_sky, stagedir)
    this_fe = filter(fe->(fe[1] == rcf), all_fes)

    @assert(length(this_fe) <= 1)
    if isempty(this_fe)
        throw(UnknownRCF(rcf))
    end

    rcf_bounds = this_fe[1][2]
    overlapping_rcfs = get_overlapping_fields(rcf_bounds, stagedir)
    @assert(rcf in overlapping_rcfs, "$rcf doesn't overlap with itself?")

    tic()
    println(rcf_bounds)
    results = one_node_infer(overlapping_rcfs,
                             stagedir;
                             objid=objid,  # just for making unit tests run fast
                             box=rcf_bounds,  # could exclude this argument
                             target_rcfs=[rcf,])
    Log.info("Inferred $rcf in $(toq()) seconds")

    fname = @sprintf "%s/celeste-%06d-%d-%04d.jld" outdir rcf.run rcf.camcol rcf.field
    JLD.save(fname, "results", results)
end

"""
Save provided results to a JLD file.
"""
function save_results(outdir, ramin, ramax, decmin, decmax, results)
    fname = @sprintf("%s/celeste-%.4f-%.4f-%.4f-%.4f.jld",
                     outdir, ramin, ramax, decmin, decmax)
    JLD.save(fname, "results", results)
end

save_results(outdir, box, results) =
    save_results(outdir, box.ramin, box.ramax, box.decmin, box.decmax, results)

"""
called from main entry point.
"""
function infer_box(box::BoundingBox, stagedir::String, outdir::String)
    # Base.@time hack for distributed environment
    gc_stats = ()
    gc_diff_stats = ()
    elapsed_time = 0.0
    gc_stats = Base.gc_num()
    elapsed_time = time_ns()

    times = InferTiming()
    if !SKY_DIVISION_PARALLELISM
        Log.debug("source division parallelism")
        divide_sources_and_infer(box, stagedir; timing=times, outdir=outdir)
    elseif dt_nnodes > 1
        Log.debug("sky division parallelism")
        divide_sky_and_infer(box, stagedir; timing=times, outdir=outdir)
    else
        Log.debug("multithreaded parallelism only")
        tic()
        # Get vector of (run, camcol, field) triplets overlapping this patch
        rcfs = get_overlapping_fields(box, stagedir)
        times.query_fids = toq()

        results = one_node_infer(rcfs, stagedir; box=box, timing=times)

        tic()
        save_results(outdir, box, results)
        times.write_results = toq()
    end

    # Base.@time hack for distributed environment
    elapsed_time = time_ns() - elapsed_time
    gc_diff_stats = Base.GC_Diff(Base.gc_num(), gc_stats)
    time_puts(elapsed_time, gc_diff_stats.allocd, gc_diff_stats.total_time,
              Base.gc_alloc_count(gc_diff_stats))

    times.num_srcs = max(1, times.num_srcs)
    nputs(dt_nodeid, "timing: query_fids=$(times.query_fids)")
    nputs(dt_nodeid, "timing: num_infers=$(times.num_infers)")
    nputs(dt_nodeid, "timing: read_photoobj=$(times.read_photoobj)")
    nputs(dt_nodeid, "timing: read_img=$(times.read_img)")
    nputs(dt_nodeid, "timing: fit_psf=$(times.fit_psf)")
    nputs(dt_nodeid, "timing: opt_srcs=$(times.opt_srcs)")
    nputs(dt_nodeid, "timing: num_srcs=$(times.num_srcs)")
    nputs(dt_nodeid, "timing: average opt_srcs=$(times.opt_srcs/times.num_srcs)")
    nputs(dt_nodeid, "timing: write_results=$(times.write_results)")
    nputs(dt_nodeid, "timing: wait_done=$(times.wait_done)")
end

end


text:
#!/usr/bin/env julia

using Plots
x = 1:10; y = rand(10); # These are the plotting data
plot(x, y)

max_stars_repo_name: jdurbin/sandbox

max_issues_repo_name: jdurbin/sandbox

max_forks_repo_name: jdurbin/sandbox


text:
@safetestset "Peer info" begin
  @safetestset "Peer info is disabled by default" begin
    using Genie, Genie.Requests
    using HTTP

    port = nothing
    port = rand(8500:8900)

    route("/") do
      "$(peer().ip)-$(peer().port)"
    end

    server = up(port)

    response = try
      HTTP.request("GET", "http://127.0.0.1:$port")
    catch ex
      ex.response
    end

    @test Genie.config.features_peerinfo == false
    @test response.status == 200
    @test String(response.body) == "-"

    down()
    sleep(1)
    server = nothing
  end;

  @safetestset "Peer info can be activated" begin
    using Genie, Genie.Requests
    using HTTP

    port = nothing
    port = rand(8500:8900)
    Genie.config.features_peerinfo = true

    route("/") do
      "$(peer().ip)-$(peer().port)"
    end

    server = up(port)

    response = try
      HTTP.request("GET", "http://127.0.0.1:$port")
    catch ex
      ex.response
    end

    @test Genie.config.features_peerinfo == true
    @test response.status == 200
    @test_broken String(response.body) == "127.0.0.1-$port"

    down()
    sleep(1)
    server = nothing
    port = nothing
  end;

end


text:
# Test default model and model with other solver
m = FlexibilityModel(solver = PavitoSolver(mip_solver = GLPKSolverMIP(),
                     cont_solver = IpoptSolver(print_level = 0), log_level = 0, mip_solver_drives = false))
m2 = FlexibilityModel()
true


text:
function figures_reg(df::DataFrame, treated::String, control::Vector{String}, brexit_dates::Vector{Date}, interval::Int64, products::Vector{String})

    # create directory for specification
    folder = ""
    for i in [treated; control] folder = folder * "_" * i end
    folder = dir_dropbox * "results/regression/" * folder

    if isdir(folder)
        println(" ✓ The specification with treated: $treated and controls: $control are overwritten.")
    else
        mkdir(folder)
        println(" ✓ The specification with treated: $treated and controls: $control are added.")
    end

    # timeframe
    timeframe = brexit_dates[1]-Month(interval+1):Month(1):brexit_dates[end]+Month(interval)

    # subset data
    df = transform(df, :PERIOD => ByRow(x -> Date(string(x), DateFormat("yyyymm"))) => :DATE)
    # subset!(df, :PARTNER_ISO => ByRow(x -> x in ["WORLD"; treated; control]), :PRODUCT_NC => ByRow(x -> x in products), :DATE => ByRow(x -> x in timeframe))
    subset!(df, :PARTNER_ISO => ByRow(x -> x in ["WORLD"; treated; control]), :PRODUCT_NC => ByRow(x -> x in products))

    # create figures

    return df, folder
end

function figure_1(df::DataFrame, folder::String, treated::String)

    # aggregate over products
    cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PERIOD", "DATE"]
    gdf = groupby(df, cols_grouping)
    df = combine(gdf, :VALUE_IN_EUROS => sum, renamecols=false)
    sort!(df)

    # HP filter
    λ = 20
    cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW"]
    gdf = groupby(df, cols_grouping)
    df = transform(gdf, :VALUE_IN_EUROS => (x -> HP(x, λ)) => :VALUE_IN_EUROS_HP)

    # create 3MMA as observations are very volatile
    cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW"]
    gdf = groupby(df, cols_grouping)
    df = transform(gdf, :VALUE_IN_EUROS => (x -> movingaverage(x,3)) => :VALUE_3MMA)

    # rolling standard deviation
    cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW"]
    gdf = groupby(df, cols_grouping)
    df = transform(gdf, [:VALUE_IN_EUROS, :VALUE_3MMA] .=> (x -> rolling_std(x, 6)) .=> [:STD_VALUE, :STD_VALUE_3MMA])

    for flow in ["exports", "imports"]

        # values
        p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PARTNER_ISO => ByRow(x -> x != "WORLD")) plot(:DATE, :VALUE_IN_EUROS/1e9,
        group=:PARTNER_ISO, lw=2, legend=:topleft, ylabel="euros (billion)", title="Flemish "* flow*": values")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
        savefig(p, folder * "/" * "fig1_" * flow * "_values" * ".png") # export image dropbox

        # HP
        p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PARTNER_ISO => ByRow(x -> x != "WORLD")) plot(:DATE, :VALUE_IN_EUROS_HP/1e9,
                group=:PARTNER_ISO, lw=2, legend=:topleft, ylabel="euros (billion)", title="Flemish "*flow*": values (HP, λ=$λ)")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
        savefig(p, folder * "/" * "fig1_" * flow * "_values" * "_HP" * ".png") # export image dropbox

        # 3MMA
        p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PARTNER_ISO => ByRow(x -> x != "WORLD")) plot(:DATE, :VALUE_3MMA/1e9,
                group=:PARTNER_ISO, lw=2, legend=:topleft, ylabel="euros (billion)", title="Flemish "*flow*": values (3MMA)")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
        savefig(p, folder * "/" * "fig1_" * flow * "_values" * "_3MMA" * ".png") # export image dropbox

        # STD
        p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PARTNER_ISO => ByRow(x -> x != "WORLD")) plot(:DATE, :STD_VALUE/1e9,
                group=:PARTNER_ISO, lw=2, legend=:topleft, ylabel="euros (billion)", title="Flemish "*flow*": 6 months STD \n (values, rolling window)")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
        savefig(p, folder * "/" * "fig4_" * flow * "_values" * "_STD" * ".png") # export image dropbox

    end

    # STD
    p = @df subset(df, :PARTNER_ISO => ByRow(x -> x == treated)) plot(:DATE, :STD_VALUE/1e9,
            group=:FLOW, lw=2, legend=:topleft, ylabel="euros (billion)", title="Flemish vs $treated: 6 months STD \n (values, rolling window)")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, folder * "/" * "fig4_" * "_values" * "_STD" * ".png") # export image dropbox

end

function figure_3(df::DataFrame, folder::String, treated::String)

    # aggregate over products
    cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PERIOD", "DATE"]
    gdf = groupby(df, cols_grouping)
    df = combine(gdf, :VALUE_IN_EUROS => sum, renamecols=false)

    # create total
    df_WORLD = subset(df, :PARTNER_ISO => ByRow(x -> x == "WORLD"))
    cols_grouping = ["DECLARANT_ISO", "FLOW", "PERIOD", "DATE"]
    gdf = groupby(df_WORLD, cols_grouping)
    df_WORLD = combine(gdf, :VALUE_IN_EUROS => sum => :VALUE_WORLD)

    # join to create shares
    df = leftjoin(subset(df, :PARTNER_ISO => ByRow(x -> x != "WORLD")), df_WORLD, on=[:DECLARANT_ISO, :FLOW, :PERIOD, :DATE])
    transform!(df, [:VALUE_IN_EUROS, :VALUE_WORLD] => ByRow((x,s) -> x/s) => :SHARE_VALUE)
    sort!(df)

    # HP filter
    λ = 20
    cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW"]
    gdf = groupby(df, cols_grouping)
    df = transform(gdf, :SHARE_VALUE => (x -> HP(x, λ)) => :SHARE_VALUE_HP)

    for flow in ["imports", "exports"]

        # values
        p = @df subset(df, :FLOW => ByRow(x -> x == flow)) plot(:DATE, :SHARE_VALUE*100,
                group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="percentages", title="Flemish "* flow*" share (values)")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
        savefig(p, folder * "/" * "fig3_" * flow * "_value_share" * ".png") # export image dropbox

        # HP
        p = @df subset(df, :FLOW => ByRow(x -> x == flow)) plot(:DATE, :SHARE_VALUE_HP*100,
            group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="percentages", title="Flemish "* flow*" share (values, λ=$λ)")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
        savefig(p, folder * "/" * "fig3_" * flow * "_value_share" * "_HP" * ".png") # export image dropbox

    end

    # VLA vis-a-vis treated
    p = @df subset(df, :PARTNER_ISO => ByRow(x -> x == treated)) plot(:DATE, :SHARE_VALUE_HP*100,
            group=:FLOW, lw=2, legend=:bottomleft, ylabel="percentages", title="Total trade share (values, λ=$λ): \n Flanders vs $treated")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, folder * "/" * "fig3_" * "VLA_treated" * "_value_share" * "_HP" * ".png") # export image

end

function figure_5(df::DataFrame, folder::String)

    # aggregate
    cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PERIOD", "DATE"]
    gdf = groupby(df, cols_grouping)
    df = combine(gdf, nrow => :COUNT)

    for flow in ["imports", "exports"]

        p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PARTNER_ISO => ByRow(x -> x != "WORLD")) plot(:DATE, :COUNT,
                group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="number of products", title="Flemish "* flow*": number of products")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
        savefig(p, folder * "/" * "fig5_" * flow * "_product_count" * ".png") # export image dropbox

    end

end

function figure_7(df::DataFrame, folder::String)

    # aggregate
    cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PERIOD", "DATE"]
    gdf = groupby(df, cols_grouping)
    df = combine(gdf, :VALUE_IN_EUROS => sum, renamecols=false)

    # compute YOY monthly percentage change
    cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW"]
    gdf = groupby(df, cols_grouping)
    df = transform(gdf, [:PERIOD, :VALUE_IN_EUROS] => yoy_change => :YOY_VALUE)
    subset!(df, :YOY_VALUE => ByRow(x -> !ismissing(x))) # need to drop missing for MAD/HP

    # # MAD adjustment
    # #   - need to drop missing
    # subset!(df, :YOY_VALUE => ByRow(x -> !ismissing(x)))
    # outlier_cutoff = 4
    # gdf = groupby(df, "PRODUCT_NC_digits")
    # df = transform(gdf, :YOY_VALUE => MAD_method => :MAD)
    # subset!(df, :MAD => ByRow(x -> x < outlier_cutoff))
    # df = df[:, Not(:MAD)]

    # HP filter
    λ = 20
    cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW"]
    gdf = groupby(df, cols_grouping)
    df = transform(gdf, :YOY_VALUE => (x -> HP(x, λ)) => :YOY_VALUE_HP)

    # std
    cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW"]
    gdf = groupby(df, cols_grouping)
    df = transform(gdf, [:YOY_VALUE, :YOY_VALUE_HP] .=> (x -> rolling_std(x, 6)) .=> [:STD_YOY, :STD_YOY_HP])

    for flow in ["imports", "exports"]

        # YOY HP
        p = @df subset(df, :FLOW => ByRow(x -> x == flow)) plot(:DATE, :STD_YOY_HP,
                group=:PARTNER_ISO, lw=2, legend=:topleft, ylabel="percentages", title="Flemish "*flow*" : 6 months STD \n (YOY change, rolling window, λ=$λ)")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
        savefig(p, folder * "/" * "fig7_" * flow * "_YOY" * "_STD_HP" * ".png") # export image dropbox

    end

end

function table_3(df::DataFrame, folder::String, brexit_dates::Vector{Date})

    # aggregate
    cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PERIOD", "DATE"]
    gdf = groupby(df, cols_grouping)
    df = combine(gdf, :VALUE_IN_EUROS => sum, renamecols=false)

    # compute YOY monthly percentage change
    cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW"]
    gdf = groupby(df, cols_grouping)
    df = transform(gdf, [:PERIOD, :VALUE_IN_EUROS] => yoy_change => :YOY_VALUE)
    subset!(df, :YOY_VALUE => ByRow(x -> !ismissing(x))) # need to drop missing for MAD/HP

    tab_std = transform(df, :DATE => ByRow(x -> ifelse(x in brexit_dates[1]-Month(1)-Year(1):Month(1):brexit_dates[1]-Month(1), "pre",
        ifelse(x in brexit_dates[1]:Month(1):brexit_dates[2]-Month(1), "brexit_1", ifelse(x in brexit_dates[2]:Month(1):brexit_dates[3]-Month(1), "brexit_2",
            ifelse(x in brexit_dates[3]:Month(1):brexit_dates[3]+Month(10), "brexit_3", "outside"))))) => :BREXIT)

    cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "BREXIT"]
    gdf = groupby(tab_std, cols_grouping)
    tab_std = combine(gdf, :YOY_VALUE => std => :STD_VALUE)

    tab_std_wide = unstack(tab_std, :BREXIT, :STD_VALUE)
    tab_std_wide = tab_std_wide[:, Not(:outside)]
    sort!(tab_std_wide, :FLOW)

    XLSX.writetable(folder * "/" * "table3_std" * ".xlsx", tab_std_wide, overwrite=true)

end


text:
using DirectDetections, Distributions, Plots
using MCMCChains

# Create a function mapping (age_Myr, mass_Mjup) -> temp_K
const cooling_tracks = sonora_cooling_interpolator()

# Create functions mapping (temp_K, mass_Mjup) -> absolute magnitude
const sonora_temp_mass_Z = sonora_photometry_interpolator(:MKO_Z)
const sonora_temp_mass_J = sonora_photometry_interpolator(:MKO_J)
const sonora_temp_mass_L = sonora_photometry_interpolator(:Keck_L′)

##
@named b = DirectDetections.Planet(
    Derived(
        Z = (sys, pl) -> sonora_temp_mass_Z(cooling_tracks(sys.age, pl.mass), pl.mass),
        J = (sys, pl) -> sonora_temp_mass_J(cooling_tracks(sys.age, pl.mass), pl.mass),
        L = (sys, pl) -> sonora_temp_mass_L(cooling_tracks(sys.age, pl.mass), pl.mass),
    ),
    Priors(
        # mass = Uniform(0.53, 95)
        # mass = TruncatedNormal(7, 7, 1, 80)
        # mass = NNor(1, 50)
        mass = Uniform(2, 60)
    ),
    # Photometry
    Photometry(
        (band = :Z, phot=16.0, σ_phot=3.),
        (band = :J, phot=17.0, σ_phot=2.),
        (band = :L, phot=12.0, σ_phot=3.0)
    )
)

@named HD12345 = System(
    Priors(
        age = TruncatedNormal(40, 15, 5, Inf)
    ),
    b,
)

##
chain = DirectDetections.hmc(
    HD12345, 0.85,
    adaptation =   1_000,
    iterations =  50_000,
)

##
gelmandiag(chain)
##
using PairPlots
PairPlots.corner(chain)

##
table = (
    mass=chain["b[mass]"],
    age=chain["age"],
)
PairPlots.corner(table)

##
p = plot()
for planet in chain.info.model.planets
    phot = planet.photometry
    N= 1500
    ii = rand(1:size(chain,1)*size(chain,3),N)
    Zs = chain["b[Z]"][ii]
    Js = chain["b[J]"][ii]
    Ls = chain["b[L]"][ii]
    plot!(vcat(Zs', Js', Ls'), color=:black, lw=1,  alpha=0.05, label="")
    scatter!(vcat(Zs', Js', Ls'), color=:black, lw=1, ms=4, alpha=0.01, label="")
    scatter!(
        collect(phot.phot),
        yerr=collect(phot.σ_phot),
        xticks=(collect(eachindex(phot.phot)), string.(phot.band)),
        color=:red,
        markerstrokecolor=:red,
        label="Measured"
    )
end
p