text:
using HTTP
using SimpleMock

@testset "save input" begin
    mktempdir() do path
        cd(path) do
            # without environment variable
            @test_throws KeyError Advent.save_input(2020, 10)

            response = HTTP.Response(500, "some error")
            mock(HTTP.get => Mock(response)) do mock_response
                @test_throws ErrorException Advent.save_input(2020, 10, "session")
                @test called_once(mock_response)
            end

            response = HTTP.Response(200, "the response")
            mock(HTTP.get => Mock(response)) do mock_response
                @test "day10.in" === Advent.save_input(2020, 10, "session")
                @test readlines("day10.in") == ["the response"]

                # HTTP.get should only be called the first time since we cache input
                @test "day10.in" === Advent.save_input(2020, 10, "session")
                @test called_once(mock_response)
            end
        end
    end
end


text:
#=
Copyright (c) 2016, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

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

module DomainPass

import ParallelAccelerator
import ParallelAccelerator.DomainIR

using CompilerTools
import CompilerTools.DebugMsg
DebugMsg.init()

using CompilerTools.LambdaHandling
using CompilerTools.Helper
#using Debug

import HPAT
import HPAT.CaptureAPI
import HPAT.CaptureAPI.getColName
import HPAT.CaptureAPI.revColName

import CompilerTools.AstWalker
import CompilerTools.AliasAnalysis

mk_alloc(typ, s) = Expr(:alloc, typ, s)
mk_call(fun,args) = Expr(:call, fun, args...)

const generatedFuncs = [:__hpat_data_source_HDF5_open,
                        :__hpat_data_source_HDF5_size,
                        :__hpat_get_H5_dim_size,
                        :__hpat_data_source_HDF5_read,
                        :__hpat_data_source_HDF5_close,
                        :__hpat_data_source_TXT_open,
                        :__hpat_data_source_TXT_size,
                        :__hpat_get_TXT_dim_size,
                        :__hpat_data_source_TXT_read,
                        :__hpat_data_source_TXT_close,
                        :__hpat_filter,
                        :__hpat_join,
                        :__hpat_aggregate,
                        ]

const generatedExprHeads = [:alloc,
                            :filter,
                            :join,
                            :aggregate]

# ENTRY to DomainPass
function from_root(function_name, ast)

    linfo, body = CompilerTools.LambdaHandling.lambdaToLambdaVarInfo(ast)
    @dprintln(1,"Starting main DomainPass.from_root.  function = ", function_name, " ast = ", linfo, body)

    tableCols, tableTypes, tableIds = get_table_meta(body)
    @dprintln(3,"HPAT tables: ", tableCols,tableTypes)
    state::DomainState = DomainState(linfo, tableCols, tableTypes, tableIds, 0)

    # transform body
    body.args = from_toplevel_body(body.args, state)
    @dprintln(1,"DomainPass.from_root returns function = ", function_name, " body = ", body)
    #println("DomainPass.from_root returns function = ", function_name, " body = ", body)
    return LambdaVarInfoToLambda(state.linfo, body.args, ParallelAccelerator.DomainIR.AstWalk)
end

# information about AST gathered and used in DomainPass
type DomainState
    linfo  :: LambdaVarInfo
    tableCols::Dict{Symbol,Vector{Symbol}}
    tableTypes::Dict{Symbol,Vector{Symbol}}
    tableIds::Dict{Int,Symbol}
    # a unique id for domain operations (data sources/sinks, table operations)
    unique_id::Int
end

"""
    get a unique id for renaming
"""
function get_unique_id(state)
  state.unique_id += 1
  return state.unique_id
end

function get_table_meta(body)
    first_arg = body.args[1]
    if isa(first_arg, Expr) && first_arg.head==:meta
        for meta in first_arg.args
            if meta.head==:hpat_tables
                @dprintln(3, "hpat tables found: ", meta)
                return meta.args[1], meta.args[2], meta.args[3]
            end
        end
    end
    return Dict{Symbol,Vector{Symbol}}(),Dict{Symbol,Vector{Symbol}}()
end

# nodes are :body of AST
function from_toplevel_body(nodes::Array{Any,1}, state::DomainState)
    res::Array{Any,1} = []
    nodes = translate_table_oprs(nodes,state)
    @dprintln(3,"body after table translation: ", Expr(:body,nodes...))

    for node in nodes
        new_exprs = from_expr(node, state)
        append!(res, new_exprs)
    end
    return res
end

function from_expr(node::Expr, state::DomainState)
    head = node.head
    if head==:(=)
        return from_assignment(node, state)
    elseif getHPATcall(node)==:data_sink_HDF5
        return translate_data_sink_HDF5(node.args[2], node.args[3], node.args[4], state)
    elseif head==:return && getHPATcall(node.args[1])==:data_sink_HDF5
        snode = node.args[1]
        s_expr = translate_data_sink_HDF5(snode.args[2], snode.args[3], snode.args[4], state)
        return Any[s_expr;Expr(:return,nothing)]
    else
        return [node]
    end
end

function from_expr(node::Any, state::DomainState)
    return [node]
end

function translate_table_oprs(nodes::Array{Any,1}, state::DomainState)
    # array of new nodes after translation
    new_nodes = []
    # number of nodes to skip in the loop
    # translating operations requires skipping some nodes after current node
    skip = 0
    for i in 1:length(nodes)
        if skip!=0
            skip-=1
            continue
        end
        if isa(nodes[i],Expr) && nodes[i].head==:(=) && (isCall(nodes[i].args[2]) || isInvoke(nodes[i].args[2]))
            local func_call = getCallFunction(nodes[i].args[2])
            local args = getCallArguments(nodes[i].args[2])
            # convert :invoke to :call to be consistent
            if isInvoke(nodes[i].args[2])
                nodes[i].args[2].head = :call
                nodes[i].args[2].args = [func_call; args]
            end
            if func_call==GlobalRef(HPAT.API, :join)
                @dprintln(3,"new nodes before join: ", Expr(:body,new_nodes...))
                remove_before,remove_after,ast = translate_join(nodes[i], nodes, i,state)
                skip += remove_after
                s_start = (length(new_nodes)-remove_before)+1
                s_end = length(new_nodes)
                # replace ast nodes with new node
                splice!(new_nodes, s_start:s_end, ast)
                @dprintln(3,"new nodes after join: ", Expr(:body,new_nodes...))
                continue
            elseif func_call==GlobalRef(HPAT.API, :aggregate)
                @dprintln(3,"new nodes before aggregate: ", Expr(:body,new_nodes...))
                remove_before,remove_after,ast = translate_aggregate(nodes, i, nodes[i],state)
                skip += remove_after
                s_start = (length(new_nodes)-remove_before)+1
                s_end = length(new_nodes)
                # replace ast nodes with new node
                splice!(new_nodes, s_start:s_end, ast)
                @dprintln(3,"new nodes after aggregate: ", Expr(:body,new_nodes...))
                continue
            elseif func_call==GlobalRef(HPAT.API, :table_filter!)
                @dprintln(3,"new nodes before filter: ", Expr(:body,new_nodes...))
                # returns: new ast :filter node
                # number of junk nodes to remove AFTER the filter call
                # number of junk nodes to remove BEFORE the filter call
                remove_before,remove_after,ast = translate_filter(nodes,i,nodes[i],state)
                skip += remove_after
                s_start = (length(new_nodes)-remove_before)+1
                s_end = length(new_nodes)
                # replace ast nodes with new node
                splice!(new_nodes, s_start:s_end, ast)
                @dprintln(3,"new nodes after filter: ", Expr(:body,new_nodes...))
                continue
            end
        # TODO: any recursive case?
        # elseif isa(nodes[i],Expr) && nodes[i].head==:block
        elseif isCall(nodes[i]) || isInvoke(nodes[i])
            func_call = getCallFunction(nodes[i])
            args = getCallArguments(nodes[i])
            # convert :invoke to :call to be consistent
            if isInvoke(nodes[i])
                nodes[i].head = :call
                nodes[i].args = [func_call; args]
            end
        end
        push!(new_nodes, nodes[i])
    end
    return new_nodes
end

""" Translate table_filter to Expr(:filter, cond_arr, t1 ,col_arrs...) and remove array of array garbage

    returns: number of junk nodes to remove before the filter call
             number of junk nodes to remove after the filter call
             new ast :filter node

    example:
    _6 = (_4 .> 2)::BitArray{1}
    _7 = (Core.ccall)(:jl_alloc_array_1d,(Core.apply_type)(Core.Array,Array{T,1},1)::Type{Array{Array{T,1},1}},(Core.svec)(Core.Any,Core.Int)::SimpleVector,(Core.apply_type)(Core.Array,Array{T,1},1)::Type{Array{Array{T,1},1}},0,2,0)::Array{Array{T,1},1}
    _8 = _4
    (ParallelAccelerator.API.setindex!)(_7,_8,1)::Array{Array{T,1},1}
    _8
    _9 = _5
    (ParallelAccelerator.API.setindex!)(_7,_9,2)::Array{Array{T,1},1}
    _9
    _10 = (HPAT.API.table_filter!)(1,3,_6,_7)
    SSAValue(0) = (ParallelAccelerator.API.getindex)(_10,1)
    _11 = (Base.convert)(Array{Int64,1},SSAValue(0))::Array{Int64,1}
    SSAValue(1) = (ParallelAccelerator.API.getindex)(_10,2)
    _12 = (Core.typeassert)((Base.convert)(Array{Float64,1},SSAValue(1)),Array{Float64,1})::Array{Float64,1}
"""
function translate_filter(nodes::Array{Any,1},curr_pos,filter_node::Expr,state)
    @dprintln(3,"translating filter: ",filter_node)

    opr_num = get_unique_id(state)

    cond_arr = toLHSVar(filter_node.args[2].args[4])
    in_arr_of_arrs = toLHSVar(filter_node.args[2].args[5])
    out_arr_of_arrs = toLHSVar(filter_node.args[1])
    # TODO: remove arr_of_arrs from Lambda
    # convert _filter_t1 to t1
    in_table_name = state.tableIds[filter_node.args[2].args[2]]
    out_table_name = state.tableIds[filter_node.args[2].args[3]]

    in_cols = state.tableCols[in_table_name]
    in_col_arrs = []  #map(x->getColName(in_table_name, x), in_cols)
    in_num_cols = length(in_cols)

    i = curr_pos-1
    while !isAllocAssignment(nodes[i])
        if isArraySet(nodes[i])
            set_in = nodes[i].args[3]
            # previous node is array_col = set_in
            @assert nodes[i-1].head==:(=) && nodes[i-1].args[1]==set_in
            array_col = nodes[i-1].args[2]
            # remove variable?
            #CompilerTools.LambdaHandling.
            push!(in_col_arrs, array_col)
        end
        i -= 1
    end
    in_col_arrs = in_col_arrs[end:-1:1]
    remove_before = curr_pos-i

    out_cols = state.tableCols[out_table_name]
    out_col_arrs = []  #map(x->getColName(out_table_name, x), out_cols)
    out_num_cols = length(out_cols)
    for k in curr_pos+2:2:curr_pos+2*out_num_cols
        #@assert nodes[i].head==:(=) && nodes[i].args[2].args[1]==GlobalRef(Core,:typeassert)
        #println(nodes[k])
        push!(out_col_arrs, nodes[k].args[1])
    end

    # remove temp array assignment and setindex!() for each column, remove array of array allocation

    # remove type convert calls after filter()
    remove_after = 2*out_num_cols
    new_filter_node = Expr(:filter, cond_arr, out_table_name, in_table_name,
                                in_cols, out_col_arrs, in_col_arrs, opr_num)
    @dprintln(3,"filter remove_before: ",remove_before," remove_after: ",remove_after," filter_node: ",filter_node)
    return remove_before, remove_after, Any[new_filter_node]
end

""" Translate join to Expr(:join, t3,t1,t2,out_cols, in1_cols, in2_cols) and remove array of array garbage

    returns: number of junk nodes to remove before the join call
             number of junk nodes to remove after the join call
             new ast :join node
    example:
    _8 = (Core.ccall)(:jl_alloc_array_1d,(Core.apply_type)(Core.Array,Array{T,1},1)::Type{Array{Array{T,1},1}},(Core.svec)(Core.Any,Core.Int)::SimpleVector,(Core.apply_type)(Core.Array,Array{T,1},1)::Type{Array{Array{T,1},1}},0,2,0)::Array{Array{T,1},1}
    _9 = _4
    (ParallelAccelerator.API.setindex!)(_8,_9,1)::Array{Array{T,1},1}
    _9
    _10 = _5
    (ParallelAccelerator.API.setindex!)(_8,_10,2)::Array{Array{T,1},1}
    _10
    _11 = (Core.ccall)(:jl_alloc_array_1d,(Core.apply_type)(Core.Array,Array{T,1},1)::Type{Array{Array{T,1},1}},(Core.svec)(Core.Any,Core.Int)::SimpleVector,(Core.apply_type)(Core.Array,Array{T,1},1)::Type{Array{Array{T,1},1}},0,2,0)::Array{Array{T,1},1}
    _12 = _6
    (ParallelAccelerator.API.setindex!)(_11,_12,1)::Array{Array{T,1},1}
    _12
    _13 = _7
    (ParallelAccelerator.API.setindex!)(_11,_13,2)::Array{Array{T,1},1}
    _13
    _14 = (HPAT.API.join)(3,1,2,_8,_11)
    SSAValue(0) = (ParallelAccelerator.API.getindex)(_14,1)
    _15 = (Base.convert)(Array{Int64,1},SSAValue(0))::Array{Int64,1}
    SSAValue(1) = (ParallelAccelerator.API.getindex)(_14,2)
    _16 = (Core.typeassert)((Base.convert)(Array{Float64,1},SSAValue(1)),Array{Float64,1})::Array{Float64,1}
    SSAValue(2) = (ParallelAccelerator.API.getindex)(_14,3)
    _17 = (Core.typeassert)((Base.convert)(Array{Float64,1},SSAValue(2)),Array{Float64,1})::Array{Float64,1}
"""
function translate_join(join_node, nodes, curr_pos, state)
    @dprintln(3, "translating join: ", join_node)

    opr_num = get_unique_id(state)

    out_arr = toLHSVar(join_node.args[1])
    local t3_id::Int = join_node.args[2].args[2]
    local t1_id::Int = join_node.args[2].args[3]
    local t2_id::Int = join_node.args[2].args[4]
    in1_arr = toLHSVar(join_node.args[2].args[5])
    in2_arr = toLHSVar(join_node.args[2].args[6])

    # convert _join_out_t3 to t3
    t3 = state.tableIds[t3_id]
    # convert _join_t1 to t1
    t1 = state.tableIds[t1_id]
    t2 = state.tableIds[t2_id]

    t3_cols = state.tableCols[t3]
    t1_cols = state.tableCols[t1]
    t2_cols = state.tableCols[t2]
    t3_num_cols = length(t3_cols)
    t1_num_cols = length(t1_cols)
    t2_num_cols = length(t2_cols)

    t2_end = 2*t2_num_cols

    t1_start = t2_end + 3
    t1_end = 2*(t1_num_cols-1) + t1_start
    # Extract inputs columns from nodes above join node
    t2_cols_sorted = []
    t1_cols_sorted = []

    # assuming Julia doesn't reorder nodes, produces same AST!
    # read columns until array of array allocation in AST
    # read t2 columns
    i = curr_pos-1
    while !isAllocAssignment(nodes[i])
        if isArraySet(nodes[i])
            set_in = nodes[i].args[3]
            # previous node is array_col = set_in
            @assert nodes[i-1].head==:(=) && nodes[i-1].args[1]==set_in
            array_col = nodes[i-1].args[2]
            # remove variable?
            #CompilerTools.LambdaHandling.
            push!(t2_cols_sorted, array_col)
        end
        i -= 1
    end
    # read t1 columns
    i -= 1
    while !isAllocAssignment(nodes[i])
        if isArraySet(nodes[i])
            set_in = nodes[i].args[3]
            # previous node is array_col = set_in
            @assert nodes[i-1].head==:(=) && nodes[i-1].args[1]==set_in "array column set expected"
            array_col = nodes[i-1].args[2]
            # remove variable?
            #CompilerTools.LambdaHandling.
            push!(t1_cols_sorted, array_col)
        end
        i -= 1
    end

    # As we read columns from bottom up we need to reverse them
    t2_cols_sorted = t2_cols_sorted[end:-1:1]
    t1_cols_sorted = t1_cols_sorted[end:-1:1]
    t3_arrs = []
    for k in curr_pos+2:2:curr_pos+2*t3_num_cols
        #@assert nodes[i].head==:(=) && nodes[i].args[2].args[1]==GlobalRef(Core,:typeassert)
        #println(nodes[k])
        push!(t3_arrs, nodes[k].args[1])
    end
    remove_before = curr_pos-i
    remove_after =  2*t3_num_cols
    new_join_node = Expr(:join, t3, t1, t2, t3_cols,
        map(x->get_col_name_from_arr(x,state.linfo), t1_cols_sorted),
        map(x->get_col_name_from_arr(x,state.linfo), t2_cols_sorted),
        t3_arrs, t1_cols_sorted, t2_cols_sorted, opr_num)
    return remove_before, remove_after, Any[new_join_node]
end

function isAllocAssignment(node::Expr)
    if node.head==:(=) && isa(node.args[2],Expr) && node.args[2].head==:call &&
          node.args[2].args[1]==GlobalRef(Core,:ccall)
        return true
    end
    return false
end

isAllocAssignment(node::ANY) = false

function isArraySet(node::Expr)
    if node.head==:call && node.args[1].name==:setindex!
        return true
    end
    return false
end

isArraySet(node::ANY) = false

# Extract col name from array var (e.g. _16 -> :col)
function get_col_name_from_arr(arr, linfo)
    col = CompilerTools.LambdaHandling.lookupVariableName(arr, linfo)
    arr = split(string(col),"@")
    return Symbol(arr[3])
end

"""   Example:
        _5 = _4
        _6 = \$(Expr(:invoke, LambdaInfo for .>(::Array{Float64,1}, ::Float64), :(ParallelAccelerator.API..>), :(_4), 2.0))
        _7 = \$(Expr(:invoke, LambdaInfo for .==(::Array{Float64,1}, ::Float64), :(ParallelAccelerator.API..==), :(_4), 1.0))
        \$(Expr(:invoke, LambdaInfo for sum(::Array{Float64,1}), :(ParallelAccelerator.API.sum), :(_5)))
        _8 = Float64
        _9 = Int64
        _10 = Int64
        SSAValue(4) = (Core.tuple)((Core.tuple)(_5,Main.sum)::Tuple{Array{Float64,1},Base.#sum},(Core.tuple)(_6,Main.length)::Tuple{BitArray{1},Base.#length},(Core.tuple)(_7,Main.length)::Tuple{BitArray{1},Base.#length})::Tuple{Tuple{Array{Float64,1},Base.#sum},Tuple{BitArray{1},Base.#length},Tuple{BitArray{1},Base.#length}}
        # meta: location abstractarray.jl vect 24
        _16 = \$(QuoteNode(Tuple{DenseArray{T,1},Function})) # line 27:
        SSAValue(5) = (Base.nfields)(SSAValue(4))::Int64
        # meta: pop location
        _11 = (HPAT.API.aggregate)(1,2,_3,\$(Expr(:invoke, LambdaInfo for copy!(::Array{Tuple{DenseArray{T,1},Function},1}, ::Tuple{Tuple{Array{Float64,1},Base.#sum},Tuple{BitArray{1},Base.#length},Tuple{BitArray{1},Base.#length}}), :(Base.copy!), :((Core.ccall)(:jl_alloc_array_1d,(Core.apply_type)(Core.Array,Tuple{DenseArray{T,1},Function},1)::Type{Array{Tuple{DenseArray{T,1},Function},1}},(Core.svec)(Core.Any,Core.Int)::SimpleVector,(Core.apply_type)(Core.Array,Tuple{DenseArray{T,1},Function},1)::Type{Array{Tuple{DenseArray{T,1},Function},1}},0,SSAValue(5),0)::Array{Tuple{DenseArray{T,1},Function},1}), SSAValue(4))))
        SSAValue(0) = (ParallelAccelerator.API.getindex)(_11,1)
        _12 = (Base.convert)(Array{Int64,1},SSAValue(0))::Array{Int64,1}
        SSAValue(1) = (ParallelAccelerator.API.getindex)(_11,2)
        _13 = (Core.typeassert)((Base.convert)(Array{Float64,1},SSAValue(1)),Array{Float64,1})::Array{Float64,1}
        SSAValue(2) = (ParallelAccelerator.API.getindex)(_11,3)
        _14 = (Base.convert)(Array{Int64,1},SSAValue(2))::Array{Int64,1}
        SSAValue(3) = (ParallelAccelerator.API.getindex)(_11,4)
        _15 = (Base.convert)(Array{Int64,1},SSAValue(3))::Array{Int64,1} # line 15:
"""
function translate_aggregate(nodes, curr_pos, aggregate_node, state)
    @dprintln(3,"translating aggregate: ",aggregate_node)

    opr_num = get_unique_id(state)

    out_arr = toLHSVar(aggregate_node.args[1])

    t1 = state.tableIds[aggregate_node.args[2].args[2]]
    t2 = state.tableIds[aggregate_node.args[2].args[3]]

    t1_cols = state.tableCols[t1]
    t2_cols = state.tableCols[t2]
    t1_num_cols = length(t1_cols)
    t2_num_cols = length(t2_cols)

    # one typeof() call for each output column except key
    remove_before = t2_num_cols-1
    # extra assignments
    remove_after =  2*t2_num_cols

    # t1's key to aggregate on
    key_arr = toLHSVar(aggregate_node.args[2].args[4])

    agg_list = []

    #   Julia sometimes inlines the aggregate expresson
    aggregate_expr = aggregate_node.args[2].args[5]
    @assert aggregate_expr.head==:invoke || aggregate_expr.head==:call "invalid aggregate AST format"
    local aggregate_expr_call = getCallFunction(aggregate_expr)
    local aggregate_expr_args = getCallArguments(aggregate_expr)

    if aggregate_expr_call.name==:vect
        agg_list = aggregate_expr_args
    else
        @assert aggregate_expr_call.name==:copy! "invalid aggregate AST format"
        i = curr_pos-1
        while !isTupleAssignment(nodes[i])
            i -= 1
        end
        remove_before += curr_pos-i
        @assert nodes[i].args[2].args[1]==GlobalRef(Core,:tuple) "expected aggregate tuple assignment"
        # example: :((Core.tuple)((Core.tuple)(_5,Main.sum)::Tuple{Array{Float64,1},Base.#sum},
        #    (Core.tuple)(_6,Main.length)...
        agg_list = nodes[i].args[2].args[2:end]
    end
    in_e_arr_list = map(x->toLHSVar(x.args[2]), agg_list)
    in_func_list = map(x->x.args[3], agg_list)
    replace_length_unique(in_func_list)

    out_col_arrs = []
    for k in curr_pos+2:2:curr_pos+2*t2_num_cols
        #@assert nodes[i].head==:(=) && nodes[i].args[2].args[1]==GlobalRef(Core,:typeassert)
        #println(nodes[k])
        push!(out_col_arrs, nodes[k].args[1])
    end
    new_aggregate_node = Expr(:aggregate, t2, t1, key_arr, in_e_arr_list, in_func_list, out_col_arrs, opr_num)
    return remove_before, remove_after, Any[new_aggregate_node]
end

# dummy function so liveness in parallel-ir can know the output type (typeOfOpr)
function length_unique(a::Array)
    r::Int = length(a)
    return r
end

function replace_length_unique(func_list)
    for (i,f) in enumerate(func_list)
        if f==GlobalRef(Main,:length_unique)
            func_list[i] = GlobalRef(HPAT.DomainPass, :length_unique)
        end
    end
end

function isTupleAssignment(node::Expr)
    if node.head==:(=) && isa(node.args[2],Expr) &&
         node.args[2].args[1]==GlobalRef(Core,:tuple)
        return true
    end
    return false
end

isTupleAssignment(node::ANY) = false

# :(=) assignment (:(=), lhs, rhs)
function from_assignment(node::Expr, state)

    # pattern match distributed calls that need domain translation
    hpat_call::Symbol = getHPATcall(node.args[2])
    if hpat_call==:null
        return [node]
    end
    return translate_hpat_dist_calls(node.args[1], node.args[2], hpat_call, state)
end

function translate_hpat_dist_calls(lhs::LHSVar, rhs::Expr, hpat_call::Symbol, state)
    if hpat_call==:data_source_HDF5
        return translate_data_source_HDF5(lhs, rhs, state)
    elseif hpat_call==:data_source_TXT
        return translate_data_source_TXT(lhs, rhs, state)
    elseif hpat_call in [:Kmeans,:LinearRegression,:NaiveBayes]
        # enable OpenMP for DAAL calls
        #HPAT.enableOMP()
        # set type for output since type inference is broken!
        in_typ = CompilerTools.LambdaHandling.getType(rhs.args[2], state.linfo)
        CompilerTools.LambdaHandling.setType(lhs, in_typ, state.linfo)
        # no change
        #return [Expr(:(=),lhs,rhs)]
        # instead of assignment, return a call with lhs as 1st argument
        return [Expr(:call,GlobalRef(HPAT.API,hpat_call), lhs,rhs.args[2:end]...)]
    else # not handled yet
        return [Expr(:(=),lhs,rhs)]
    end
end

function translate_hpat_dist_calls(lhs::ANY, rhs::ANY, hpat_call::Symbol, state)
    return Any[]
end

function getHPATcall(call::Expr)
    if isCall(call) || isInvoke(call)
        local func = getCallFunction(call)
        local args = getCallArguments(call)
        # convert :invoke to :call to be consistent
        if isInvoke(call)
            call.head = :call
            call.args = [func; args]
        end
        # TODO: hack to get around Julia 0.4 function resolution issue (q26)
        # remove in 0.5
        if call.args[1]==GlobalRef(Main,:Kmeans)
            call.args[1]=GlobalRef(HPAT.API,:Kmeans)
        end
        if call.args[1]==GlobalRef(Main,:LinearRegression)
            call.args[1]=GlobalRef(HPAT.API,:LinearRegression)
        end
        if call.args[1]==GlobalRef(Main,:NaiveBayes)
            call.args[1]=GlobalRef(HPAT.API,:NaiveBayes)
        end
        return getHPATcall_inner(call.args[1])
    end
    return :null
end

function getHPATcall(call::ANY)
    return :null
end

function getHPATcall_inner(func::GlobalRef)
    if func.mod==HPAT.API
        return func.name
    end
    return :null
end

function getHPATcall_inner(func::ANY)
    return :null
end

function translate_data_source_HDF5(lhs::LHSVar, rhs::Expr, state)
    res = Any[]
    dprintln(3,"HPAT data source found ", rhs)
    hdf5_var = rhs.args[3]
    hdf5_file = rhs.args[4]
    # update counter and get data source number
    dsrc_num = get_unique_id(state)
    dsrc_id_var = addTempVariable(Int64, state.linfo)
    push!(res, TypedExpr(Int64, :(=), dsrc_id_var, dsrc_num))
    # get array type
    arr_typ = getType(lhs, state.linfo)
    dims = ndims(arr_typ)
    elem_typ = eltype(arr_typ)
    # generate open call
    # lhs is dummy argument so ParallelIR wouldn't reorder
    open_call = mk_call(GlobalRef(HPAT.API,:__hpat_data_source_HDF5_open), [dsrc_id_var, hdf5_var, hdf5_file, lhs])
    push!(res, open_call)
    # generate array size call
    # arr_size_var = addTempVariable(Tuple, state.linfo)
    # assume 1D for now
    arr_size_var = addTempVariable(ParallelAccelerator.H5SizeArr_t, state.linfo)
    size_call = mk_call(GlobalRef(HPAT.API,:__hpat_data_source_HDF5_size), [dsrc_id_var, lhs])
    push!(res, TypedExpr(arr_size_var, :(=), arr_size_var, size_call))
    # generate array allocation
    size_expr = Any[]
    for i in dims:-1:1
        size_i_name = Symbol("__hpat_h5_dim_size_"*string(dsrc_num)*"_"*string(i))
        size_i = toLHSVar(CompilerTools.LambdaHandling.addLocalVariable(
            size_i_name, Int64, ISASSIGNEDONCE | ISASSIGNED, state.linfo))
        # size_i = addTempVariable(Int64, state.linfo)
        size_i_call = mk_call(GlobalRef(HPAT.API,:__hpat_get_H5_dim_size), [arr_size_var, i])
        push!(res, TypedExpr(Int64, :(=), size_i, size_i_call))
        push!(size_expr, size_i)
    end
    arrdef = TypedExpr(arr_typ, :alloc, elem_typ, size_expr)
    push!(res, TypedExpr(arr_typ, :(=), lhs, arrdef))
    # generate read call
    read_call = mk_call(GlobalRef(HPAT.API,:__hpat_data_source_HDF5_read), [dsrc_id_var, lhs])
    push!(res, read_call)
    close_call = mk_call(GlobalRef(HPAT.API,:__hpat_data_source_HDF5_close), [dsrc_id_var])
    push!(res, close_call)
    return res
end

function translate_data_source_TXT(lhs::LHSVar, rhs::Expr, state)
    dprintln(3,"TXT data source found ", rhs)
    res = Any[]
    txt_file = rhs.args[3]
    # update counter and get data source number
    dsrc_num = get_unique_id(state)
    dsrc_id_var = addTempVariable(Int64, state.linfo)
    push!(res, TypedExpr(Int64, :(=), dsrc_id_var, dsrc_num))
    # get array type
    arr_typ = getType(lhs, state.linfo)
    dims = ndims(arr_typ)
    elem_typ = eltype(arr_typ)
    # generate open call
    # lhs is dummy argument so ParallelIR wouldn't reorder
    open_call = mk_call(:__hpat_data_source_TXT_open, [dsrc_id_var, txt_file, lhs])
    push!(res, open_call)
    # generate array size call
    # arr_size_var = addTempVariable(Tuple, state.linfo)
    arr_size_var = addTempVariable(ParallelAccelerator.SizeArr_t, state.linfo)
    size_call = mk_call(:__hpat_data_source_TXT_size, [dsrc_id_var, lhs])
    push!(res, TypedExpr(arr_size_var, :(=), arr_size_var, size_call))
    # generate array allocation
    size_expr = Any[]
    for i in dims:-1:1
        size_i_name = Symbol("__hpat_txt_dim_size_"*string(dsrc_num)*"_"*string(i))
        size_i = toLHSVar(CompilerTools.LambdaHandling.addLocalVariable(
            size_i_name, Int64, ISASSIGNEDONCE | ISASSIGNED, state.linfo))
        #size_i = addTempVariable(Int64, state.linfo)
        size_i_call = mk_call(:__hpat_get_TXT_dim_size, [arr_size_var, i])
        push!(res, TypedExpr(Int64, :(=), size_i, size_i_call))
        push!(size_expr, size_i)
    end
    arrdef = TypedExpr(arr_typ, :alloc, elem_typ, size_expr)
    push!(res, TypedExpr(arr_typ, :(=), lhs, arrdef))
    # generate read call
    read_call = mk_call(:__hpat_data_source_TXT_read, [dsrc_id_var, lhs])
    push!(res, read_call)
    close_call = mk_call(:__hpat_data_source_TXT_close, [dsrc_id_var])
    push!(res, close_call)
    return res
end

            #=
            elseif isa(inner_call.args[1],GlobalRef) && inner_call.args[1].name==:__hpat_Kmeans
                dprintln(3,"kmeans found ", inner_call)
                HPAT.enableOMP()
                lib_call = mk_call(:__hpat_Kmeans, [lhs,inner_call.args[2], inner_call.args[3]])
                return [lib_call]
            elseif isa(inner_call.args[1],GlobalRef) && inner_call.args[1].name==:__hpat_LinearRegression
                dprintln(3,"LinearRegression found ", inner_call)
                HPAT.enableOMP()
                lib_call = mk_call(:__hpat_LinearRegression, [lhs,inner_call.args[2], inner_call.args[3]])
                return [lib_call]
            elseif isa(inner_call.args[1],GlobalRef) && inner_call.args[1].name==:__hpat_NaiveBayes
                dprintln(3,"NaiveBayes found ", inner_call)
                HPAT.enableOMP()
                lib_call = mk_call(:__hpat_NaiveBayes, [lhs,inner_call.args[2], inner_call.args[3], inner_call.args[4]])
                return [lib_call]
            end
        end
    end

    return Any[]
end

function pattern_match_hpat_dist_calls(lhs::Any, rhs::Any, state)
    return Any[]
end
=#
function translate_data_sink_HDF5(y, hdf5_var, hdf5_file, state)
    res = Any[]
    dprintln(3,"HPAT data sink found ", y)
    # update counter and get data source number
    dsrc_num = get_unique_id(state)
    dsrc_id_var = addTempVariable(Int64, state.linfo)
    push!(res, TypedExpr(Int64, :(=), dsrc_id_var, dsrc_num))
    # get array type
    arr_typ = getType(y, state.linfo)
    dims = ndims(arr_typ)
    elem_typ = eltype(arr_typ)
    # generate open call
    # y is dummy argument so ParallelIR wouldn't reorder
    open_call = mk_call(GlobalRef(HPAT.API,:__hpat_data_sink_HDF5_open), [dsrc_id_var, hdf5_var, hdf5_file, y])
    push!(res, open_call)
    #=
    # generate array size call
    # arr_size_var = addTempVariable(Tuple, state.linfo)
    # assume 1D for now
    arr_size_var = addTempVariable(ParallelAccelerator.H5SizeArr_t, state.linfo)
    size_call = mk_call(GlobalRef(HPAT.API,:__hpat_data_source_HDF5_size), [dsrc_id_var, y])
    push!(res, TypedExpr(arr_size_var, :(=), arr_size_var, size_call))
    # generate array allocation
    size_expr = Any[]
    for i in dims:-1:1
        size_i = symbol("__hpat_h5_dim_size_"*string(dsrc_num)*"_"*string(i))
        CompilerTools.LambdaHandling.addLocalVariable(size_i, Int64, ISASSIGNEDONCE | ISASSIGNED, state.linfo)
        # size_i = addTempVariable(Int64, state.linfo)
        size_i_call = mk_call(GlobalRef(HPAT.API,:__hpat_get_H5_dim_size), [arr_size_var, i])
        push!(res, TypedExpr(Int64, :(=), size_i, size_i_call))
        push!(size_expr, size_i)
    end
    arrdef = TypedExpr(arr_typ, :alloc, elem_typ, size_expr)
    push!(res, TypedExpr(arr_typ, :(=), y, arrdef))
    =#
    # generate read call
    write_call = mk_call(GlobalRef(HPAT.API,:__hpat_data_sink_HDF5_write), [dsrc_id_var, hdf5_var, y])
    push!(res, write_call)
    close_call = mk_call(GlobalRef(HPAT.API,:__hpat_data_sink_HDF5_close), [dsrc_id_var])
    push!(res, close_call)
    return res
end

function AstWalkCallback(node::Expr,dw)

    if node.head==:filter
        # Structure: cond, output table, input table, cols, output table cols, input table cols, id
        cond_arr = node.args[1]
        node.args[1] = AstWalker.AstWalk(cond_arr, ParallelAccelerator.DomainIR.AstWalkCallback, dw)
        #out_t = node.args[2]
        in_t = node.args[3]
        node.args[3] = AstWalker.AstWalk(node.args[3], ParallelAccelerator.DomainIR.AstWalkCallback, dw)
        node.args[3]!=in_t && @dprintln(3,"Mapping from :",in_col," to ",in_col_arrs[i])
        col_arrs = node.args[4]
        #out_col_arrs = node.args[5]
        in_col_arrs = node.args[6]
        for i in 1:length(col_arrs)
            #out_col = out_col_arrs[i]
            in_col = in_col_arrs[i]
            #out_col_arrs[i] = AstWalker.AstWalk(out_col_arrs[i], ParallelAccelerator.DomainIR.AstWalkCallback, dw)
            in_col_arrs[i] = AstWalker.AstWalk(in_col_arrs[i], ParallelAccelerator.DomainIR.AstWalkCallback, dw)
            #out_col_arrs!=out_col && @dprintln(3,"Mapping from :",out_col," to ",out_col_arrs[i])
            in_col_arrs!=in_col && @dprintln(3,"Mapping from :",in_col," to ",in_col_arrs[i])
            #@assert col_arrs[i]==old_arr "Table column name shouldn't change for now"
        end
        return node
    elseif node.head==:join
        # Structure: output table, input table 1, input table 2, out cols, in 1 cols, in 2 cols , output table cols, input table 1 cols , input table 2 cols, id
        in_t1 = node.args[2]
        in_t2 = node.args[3]
        node.args[2] = AstWalker.AstWalk(node.args[2], ParallelAccelerator.DomainIR.AstWalkCallback, dw)
        node.args[2]!=in_t1 && @dprintln(3,"Mapping from :",in_t1," to ",node.args[2])
        node.args[3] = AstWalker.AstWalk(node.args[3], ParallelAccelerator.DomainIR.AstWalkCallback, dw)
        node.args[3]!=in_t2 && @dprintln(3,"Mapping from :",in_t2," to ",node.args[3])

        t3_arrs = node.args[7]
        t1_arrs = node.args[8]
        t2_arrs = node.args[9]
        for i in 1:length(t3_arrs)
            old_arr = t3_arrs[i]
            t3_arrs[i] = AstWalker.AstWalk(t3_arrs[i], ParallelAccelerator.DomainIR.AstWalkCallback, dw)
            t3_arrs[i]!=old_arr && @dprintln(3,"Mapping from :",old_arr," to ",t3_arrs[i])
        end
        for i in 1:length(t1_arrs)
            old_arr = t1_arrs[i]
            t1_arrs[i] = AstWalker.AstWalk(t1_arrs[i], ParallelAccelerator.DomainIR.AstWalkCallback, dw)
            t1_arrs[i]!=old_arr && @dprintln(3,"Mapping from :",old_arr," to ",t1_arrs[i])
        end
        for i in 1:length(t2_arrs)
            old_arr = t2_arrs[i]
            t2_arrs[i] = AstWalker.AstWalk(t2_arrs[i], ParallelAccelerator.DomainIR.AstWalkCallback, dw)
            t2_arrs[i]!=old_arr && @dprintln(3,"Mapping from :",old_arr," to ",t2_arrs[i])
        end
        return node
    elseif node.head==:aggregate
        # Structure: output table, input table, groupby key, output expr list, expr rhs list, func list, output list, id
        in_t = node.args[2]
        key_arr = node.args[3]
        out_e_arrs = node.args[4]
        func_arrs = node.args[5]
        out_col_arrs = node.args[6]

        node.args[2] = AstWalker.AstWalk(in_t, ParallelAccelerator.DomainIR.AstWalkCallback, dw)
        node.args[3] = AstWalker.AstWalk(key_arr, ParallelAccelerator.DomainIR.AstWalkCallback, dw)
        for i in 1:length(out_col_arrs)
            old_arr = out_col_arrs[i]
            out_col_arrs[i] = AstWalker.AstWalk(out_col_arrs[i], ParallelAccelerator.DomainIR.AstWalkCallback, dw)
        end
        return node
    end
    return CompilerTools.AstWalker.ASTWALK_RECURSE
end

function live_cb(node::Expr)

    if node.head==:filter
        @dprintln(3,"DomainPass filter live CB on: ",node)
        # read condition array and column arrays
        # write all column arrays
        cond_arr = node.args[1]
        t = node.args[2]
        cols = node.args[3]
        col_arrs = node.args[4]
        # artificial assignments to signify write to LivenessAnalysis
        assign_exprs = map(x->Expr(:(=),x,1), col_arrs)
        exprs_to_process = Any[cond_arr;col_arrs;assign_exprs]
        @dprintln(3,"DomainPass filter live CB returns: ", exprs_to_process)
        return exprs_to_process
    elseif node.head==:join
        @dprintln(3,"DomainPass join live CB on: ",node)
        t3_arrs = node.args[7]
        t1_arrs = node.args[8]
        t2_arrs = node.args[9]
        # only t3 is written
        assign_exprs = map(x->Expr(:(=),x,1), t3_arrs)
        exprs_to_process = Any[t1_arrs;t2_arrs;assign_exprs]
        @dprintln(3,"DomainPass join live CB returns: ", exprs_to_process)
        return exprs_to_process
    elseif node.head==:aggregate
        @dprintln(3,"DomainPass aggregate live CB on: ",node)
        key_arr = node.args[3]
        in_e_arrs = node.args[4]
        out_col_arrs = node.args[5]
        assign_exprs = map(x->Expr(:(=),x,1), out_col_arrs)
        exprs_to_process = [key_arr;in_e_arrs;assign_exprs]
        @dprintln(3,"DomainPass aggregate live CB returns: ", exprs_to_process)
        return exprs_to_process
    end
    return nothing
end

function alias_cb(node::Expr)
    if node.head==:filter || node.head==:join || node.head==:aggregate
        return CompilerTools.AliasAnalysis.NotArray
    end
    return nothing
end

end # module

ext: jl


text:
module MimiPAGE2020

using Mimi
using Random
using Distributions
using DataFrames

import Random

export getpage

include("utils/load_parameters.jl")
include("utils/mctools.jl")

include("models/mcs_def.jl")
include("compute_scc.jl")

include("components/RCPSSPScenario.jl")
include("components/CO2emissions.jl")
include("components/CO2cycle.jl")
include("components/CO2forcing.jl")
include("components/CH4emissions.jl")
include("components/CH4cycle.jl")
include("components/CH4forcing.jl")
include("components/N2Oemissions.jl")
include("components/N2Ocycle.jl")
include("components/N2Oforcing.jl")
include("components/LGemissions.jl")
include("components/LGcycle.jl")
include("components/LGforcing.jl")
include("components/SulphateForcing.jl")
include("components/TotalForcing.jl")
include("components/ClimateTemperature.jl")
include("components/SeaLevelRise.jl")
include("components/PAGE09SeaLevelRise.jl")
include("components/GDP.jl")
include("components/MarketDamages.jl")
include("components/MarketDamagesBurke.jl")
include("components/NonMarketDamages.jl")
include("components/Discontinuity.jl")
include("components/PAGE09Discontinuity.jl")
include("components/AdaptationCosts.jl")
include("components/SLRDamages.jl")
include("components/AbatementCostParameters.jl")
include("components/AbatementCosts.jl")
include("components/TotalAbatementCosts.jl")
include("components/TotalAdaptationCosts.jl")
include("components/Population.jl")
include("components/TotalCosts.jl")
include("components/EquityWeighting.jl")
include("components/PermafrostSiBCASA.jl")
include("components/PermafrostJULES.jl")
include("components/PermafrostTotal.jl")

include("models/main_model_def.jl")
include("models/climate_model_def.jl")

end


text:
using SerialCOBS
using Test

# %%
portname = "/dev/tty.usbmodem14201"
holybro_baudrate = 57600

include("$(@__DIR__)/test.jl")


text:
### BasicGibbsJob

mutable struct BasicGibbsJob <: GibbsJob
  model::GenericModel
  dpindex::IntegerVector # Indices of dependent variables (parameters and transformations) in model.vertices
  dependent::Vector{Union{Parameter, Transformation}} # Points to model.vertices[dpindex] for faster access
  dpjob::Vector{Union{BasicMCJob, Void}} # BasicMCJobs for parameters that will be sampled via Monte Carlo methods
  range::BasicMCRange
  vstate::VariableStateVector # Vector of variable states ordered according to variables in model.vertices
  dpstate::VariableStateVector # Points to vstate[dpindex] for faster access
  outopts::Vector # Options related to output
  output::Vector{Union{VariableNState, VariableIOStream, Void}} # Output of model's dependent variables
  ndp::Integer # Number of dependent variables, i.e. length(dependent)
  hasdpjob::Bool
  hasoutput::Bool
  hasiostream::Bool
  count::Integer # Current number of post-burnin iterations
  verbose::Bool

  function BasicGibbsJob(
    model::GenericModel,
    dpindex::IntegerVector,
    dpjob::Vector,
    range::BasicMCRange,
    vstate::VariableStateVector,
    outopts::Vector,
    verbose::Bool,
    check::Bool
  )
    instance = new()

    instance.model = model
    instance.dpindex = dpindex
    instance.dpjob = isa(dpjob, Vector{Union{BasicMCJob, Void}}) ? dpjob : convert(Vector{Union{BasicMCJob, Void}}, dpjob)
    instance.vstate = vstate

    if check
      checkin(instance)
    end

    instance.ndp = length(dpindex)

    instance.hasdpjob = any(j -> j != nothing, instance.dpjob)

    instance.outopts = isa(outopts, Vector{Dict{Symbol, Any}}) ? outopts : convert(Vector{Dict{Symbol, Any}}, outopts)

    instance.range = range
    instance.verbose = verbose

    instance.dependent = instance.model.vertices[instance.dpindex]
    for i in 1:instance.ndp
      instance.dependent[i].states = instance.vstate
    end

    instance.dpstate = instance.vstate[instance.dpindex]

    instance.output = Array{Union{VariableNState, VariableIOStream, Void}}(instance.ndp)
    for i in 1:instance.ndp
      if isa(instance.dependent[i], Parameter)
        augment_parameter_outopts!(instance.outopts[i])
      else
        augment_variable_outopts!(instance.outopts[i])
      end
      instance.output[i] = initialize_output(instance.dpstate[i], range.npoststeps, instance.outopts[i])
    end

    instance.hasoutput = any(o -> o != nothing, instance.output)

    instance.hasiostream = any(o -> o == VariableIOStream, instance.output)

    instance.count = 0

    instance
  end
end

BasicGibbsJob(
  model::GenericModel,
  dpjob::Vector,
  range::BasicMCRange,
  v0::VariableStateVector;
  dpindex::IntegerVector=find(v::Variable -> isa(v, Parameter) || isa(v, Transformation), model.vertices),
  outopts::Vector=[Dict(:destination=>:nstate, :monitor=>[:value]) for i in 1:length(dpindex)],
  verbose::Bool=false,
  check::Bool=false
) =
  BasicGibbsJob(model, dpindex, dpjob, range, v0, outopts, verbose, check)

function BasicGibbsJob(
  model::GenericModel,
  dpjob::Dict,
  range::BasicMCRange,
  v0::Dict{Symbol, S};
  dpindex::IntegerVector=find(v::Variable -> isa(v, Parameter) || isa(v, Transformation), model.vertices),
  outopts::Dict=Dict([(k, Dict(:destination=>:nstate, :monitor=>[:value])) for k in keys(model.vertices[dpindex])]),
  verbose::Bool=false,
  check::Bool=false
) where S<:VariableState
  ndpindex = length(dpindex)

  jobs = Array{Union{BasicMCJob, Void}}(ndpindex)
  opts = Array{Dict}(ndpindex)

  for i in 1:ndpindex
    key = vertex_key(model.vertices[dpindex[i]])
    jobs[i] = haskey(dpjob, key) ? dpjob[key] : nothing
    opts[i] = haskey(outopts, key) ? outopts[key] : Dict(:destination=>:nstate, :monitor=>[:value])
  end

  vstate = Array{S}(length(v0))
  for (k, v) in v0
    vstate[model.ofkey[k]] = v
  end

  BasicGibbsJob(model, dpindex, jobs, range, vstate, opts, verbose, check)
end

function BasicGibbsJob(
  model::GenericModel,
  dpjob::Vector,
  range::BasicMCRange,
  v0::Vector;
  dpindex::IntegerVector=find(v::Variable -> isa(v, Parameter) || isa(v, Transformation), model.vertices),
  outopts::Vector=[Dict(:destination=>:nstate, :monitor=>[:value]) for i in 1:length(dpindex)],
  verbose::Bool=false,
  check::Bool=false
)
  vstate = default_state(model.vertices, v0, outopts, dpindex)
  BasicGibbsJob(model, dpindex, dpjob, range, vstate, outopts, verbose, check)
end

function BasicGibbsJob(
  model::GenericModel,
  dpjob::Dict,
  range::BasicMCRange,
  v0::Dict;
  dpindex::IntegerVector=find(v::Variable -> isa(v, Parameter) || isa(v, Transformation), model.vertices),
  outopts::Dict=Dict([(k, Dict(:destination=>:nstate, :monitor=>[:value])) for k in keys(model.vertices[dpindex])]),
  verbose::Bool=false,
  check::Bool=false
)
  vstate = Dict{Symbol, VariableState}()
  for (k, v) in v0
    vstate[k] = default_state(model[k], v, get!(outopts, k, Dict()))
  end

  BasicGibbsJob(model, dpjob, range, vstate, dpindex=dpindex, outopts=outopts, verbose=verbose, check=check)
end

function close(job::BasicGibbsJob)
  for i in 1:job.ndp
    if isa(job.output[i], VariableIOStream)
      close(job.output[i])
    end
  end
end

function reset(job::BasicGibbsJob)
  for i in 1:job.ndp
    if isa(job.dpjob[i], BasicMCJob)
      if job.dpjob[i].resetpstate
        reset(job.dpjob[i], rand(job.dependent[i].prior))
      else
        reset(job.dpjob[i])
      end
    end
  end
end

function save(job::BasicGibbsJob, i::Integer)
  for j in 1:job.ndp
    if isa(job.output[j], VariableNState)
      copy!(job.output[j], job.dpstate[j], i)
    elseif isa(job.output[j], VariableIOStream)
      write(job.output[j], job.dpstate[j])
      if job.outopts[j][:flush]
        flush(job.outopts[j])
      end
    else
      error("To save output, :destination must be set to :nstate or :iostream, got $(job.outopts[j][:destination])")
    end
  end
end

function iterate!(job::BasicGibbsJob)
  for i in 1:job.ndp
    if isa(job.dependent[i], Parameter)
      if isa(job.dpjob[i], BasicMCJob)
        run(job.dpjob[i])
        job.dpstate[i] = job.dpjob[i].pstate
      else
        setpdf!(job.dependent[i], job.dpstate[i])
        job.dpstate[i].value = rand(job.dependent[i].pdf)
      end
    else
      job.dependent[i].transform!(job.dpstate[i])
    end
  end
end

function run(job::BasicGibbsJob)
  fmt_iter::Union{Function, Void} = job.verbose ? format_iteration(ndigits(job.range.nsteps)) : nothing

  if job.verbose
    println("Iteration ", fmt_iter(i), " of ", job.range.nsteps)
  end

  for i = 1:job.range.nsteps
    iterate!(job)

    if in(i, job.range.postrange)
      job.count += 1

      if job.hasoutput
        save(job, job.count)
      end
    end

    if job.hasdpjob
      reset(job)
    end
  end

  if isa(job.output, VariableIOStream)
    close(job.output)
  end

  if job.hasiostream
    close(job)
  end
end

function checkin(job::BasicGibbsJob)
  dpindex = find(v::Variable -> isa(v, Parameter) || isa(v, Transformation), job.model.vertices)
  ndp = length(dpindex)

  if ndp <= 0
    error("The model has neither parameters nor transformations, but at least one of them is required in a BasicGibbsJob")
  end

  ndpindex = length(job.dpindex)

  @assert ndpindex > 0 "Length of dpindex must be positive, got $ndpindex"

  if ndp == ndpindex
    if dpindex != job.dpindex
      error("Indices of parameters and transformations in model.vertices not same as dpindex")
    end
  elseif ndp < ndpindex
    error("Number of parameters and transformations ($ndp) in model.vertices less than length of dpindex ($ndpindex)")
  else # ndp > ndpindex
    warn("Number of parameters and transformations ($ndp) in model.vertices greater than length of dpindex ($ndpindex)")
    if in(job.dpindex, dpindex)
      error("Indices of parameters and transformations in model.vertices do not contain dpindex")
    end
  end

  ndpjob = length(job.dpjob)
  if ndpindex != ndpjob
    error("Length of dpindex ($ndpindex) not equal to length of length of dpjob ($ndpjob)")
  else
    for i in 1:ndpindex
      if isa(job.dpjob[i], BasicMCJob) && !isa(job.model.vertices[job.dpindex[i]], Parameter)
        error("BasicMCJob specified for model.vertices[$(job.dpindex[i])] variable, which is not a parameter")
      end
    end
  end

  nv = num_vertices(job.model)
  nvstate = length(job.vstate)

  if nv != nvstate
    warn("Number of variables ($nv) not equal to number of variable states ($nvstate)")
  end

  for i in ndpindex
    if isa(job.model.vertices[i], Parameter)
      if !isa(job.vstate[i], ParameterState)
        error("The parameter's state must be saved in a ParameterState subtype, got $(typeof(job.vstate[i])) state type")
      else
        check_support(job.model.vertices[i], job.vstate[i])
      end
    end
  end
end

num_dp(job::BasicGibbsJob) = job.ndp
num_dptransforms(job::BasicGibbsJob) = count(v -> isa(v, Transformation), job.dependent)
num_randdp(job::BasicGibbsJob) = count(v -> isa(v, Parameter), job.dependent)
num_randdp_viamcmc(job::BasicGibbsJob) = count(j -> j != nothing, job.dpjob)
num_radndp_viadistribution(job::BasicGibbsJob) = num_dp(job)-num_dptransforms(job)-num_randdp_viamcmc(job)

dpkeys(job::BasicGibbsJob) = Symbol[dp.key for dp in job.dependent]

output(job::BasicGibbsJob) = job.output

Dict(job::BasicGibbsJob, field::Symbol=:output) = Dict(zip(dpkeys(job), getfield(job, field)))

function show(io::IO, job::BasicGibbsJob)
  ndptransforms = num_dptransforms(job)
  ndpviamcmc = num_randdp_viamcmc(job)
  ndpviadistribution = job.ndp-ndptransforms-ndpviamcmc

  indentation = "  "

  println(io, "BasicGibbsJob:")
  print(io, string(
    indentation,
    "$(num_dp(job)) dependent variables: ",
    "$ndptransforms transformations, ",
    "$ndpviadistribution sampled from their distribution, ",
    "$ndpviamcmc sampled via MCMC-within-Gibbs"
  ))
  print(io, "\n"*indentation)
  show(io, job.model)
  print(io, "\n"*indentation)
  show(io, job.range)
end

function job2dot(stream::IOStream, job::BasicGibbsJob)
  graphkeyword, edgesign = is_directed(job.model) ? ("digraph", "->") : ("graph", "--")
  dotindentation, dotspacing = "  ", " "

  write(stream, "$graphkeyword BasicGibbsJob {\n")

  invdpindex = Dict{Symbol, Integer}(zip(keys(job.model.vertices[job.dpindex]), 1:job.ndp))

  for v in vertices(job.model)
    vstring = string(dotindentation, v.key, dotspacing, "[shape=", dotshape(v))

    if isa(v, Parameter) || isa(v, Transformation)
      vstring *= ","*dotspacing*"peripheries=2"

      i = invdpindex[v.key]

      if job.outopts[i][:destination] != :none
        vstring *= ","*dotspacing*"label=<<u>$(v.key)</u>>"
      end

      if isa(v, Parameter) && (job.dpjob[i] != nothing)
        vstring *= ","*dotspacing*"style=diagonals"
      end
    end

    write(stream, vstring*"]\n")
  end

  for d in edges(job.model)
    write(stream, string(dotindentation, d.source.key, dotspacing, edgesign, dotspacing, d.target.key, "\n"))
  end

  write(stream, "}\n")
end

function job2dot(filename::AbstractString, job::BasicGibbsJob, mode::AbstractString="w")
  stream = open(filename, mode)
  job2dot(stream, job)
  close(stream)
end


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

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

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

T = Float64

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

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

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

Gn = std(G)

G ./= Gn

d = [d1; d2; d3]

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

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

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

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

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

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

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

max_stars_repo_name: jbmuir/DAS-Reconstruction

max_issues_repo_name: jbmuir/DAS-Reconstruction

max_forks_repo_name: jbmuir/DAS-Reconstruction


text:
using Test, Hijack

@test true

@testset "runtests" begin
    @test true
end

include("included.jl")


text:
module PhotonTransferCurves

using Distributions
using Random
using Statistics
using StatsBase

export testdata, ShotNoise, ReadNoise, noise

img_to_curve(cube::AbstractArray; dims) = map(img_to_curve, eachslice(cube; dims))
img_to_curve(img::AbstractMatrix) = mean_and_std(img)
img_to_curve(pair::Tuple) = pair

function signal_and_noise(cube::AbstractArray; dims)
    data = img_to_curve(cube; dims)
    S = map(first, data)
    σ = map(last, data)
    return S, σ
end

function signal_and_noise(images::Vararg{<:AbstractMatrix})
    data = map(img_to_curve, images)
    S = map(first, data)
    σ = map(last, data)
    return S, σ
end

include("data.jl")
include("sources.jl")

end


text:
# Copyright (c) 2020-2021, Bank of Canada
# All rights reserved.

using OrderedCollections

# -------------------------------------------------------------------------------
# MVTSeries -- multivariate TSeries
# -------------------------------------------------------------------------------

mutable struct MVTSeries{
    F<:Frequency,
    T<:Number,
    C<:AbstractMatrix{T}
} <: AbstractMatrix{T}

    firstdate::MIT{F}
    columns::OrderedDict{Symbol,TSeries{F,T}}
    values::C

    # inner constructor enforces constraints
    function MVTSeries(firstdate::MIT{F},
        names::NTuple{N,Symbol},
        values::AbstractMatrix
    ) where {F<:Frequency,N}
        if N != size(values, 2)
            ArgumentError("Number of names and columns don't match:" *
                          " $N ≠ $(size(values, 2)).") |> throw
        end
        columns = OrderedDict(nm => TSeries(firstdate, view(values, :, ind))
                              for (nm, ind) in zip(names, axes(values, 2)))
        new{F,eltype(values),typeof(values)}(firstdate, columns, values)
    end

end

@inline _names_as_tuple(names::Symbol) = (names,)
@inline _names_as_tuple(names::AbstractString) = (Symbol(names),)
@inline _names_as_tuple(names) = tuple(Symbol.(names)...)

# standard constructor with default empty values
MVTSeries(fd::MIT = 1U, names = ()) = (names = _names_as_tuple(names); MVTSeries(fd, names, zeros(0, length(names))))
MVTSeries(fd::MIT, names, data::AbstractMatrix) = (names = _names_as_tuple(names); MVTSeries(fd, names, data))
# see more constructors below

@inline _vals(x::MVTSeries) = getfield(x, :values)
@inline _cols(x::MVTSeries) = getfield(x, :columns)

columns(x::MVTSeries) = getfield(x, :columns)

@inline colnames(x::MVTSeries) = keys(_cols(x))
@inline rawdata(x::MVTSeries) = _vals(x)
@inline Base.pairs(x::MVTSeries) = pairs(_cols(x))

@inline firstdate(x::MVTSeries) = getfield(x, :firstdate)
@inline lastdate(x::MVTSeries) = firstdate(x) + size(_vals(x), 1) - one(firstdate(x))
@inline frequencyof(::Type{<:MVTSeries{F}}) where {F<:Frequency} = F
@inline rangeof(x::MVTSeries) = firstdate(x):lastdate(x)

# -------------------------------------------------------------------------------
# Make MVTSeries work properly as an AbstractArray

@inline Base.size(x::MVTSeries) = size(_vals(x))
@inline Base.axes(x::MVTSeries) = (rangeof(x), tuple(colnames(x)...))
@inline Base.axes1(x::MVTSeries) = rangeof(x)

# the following are needed for copy() and copyto!() (and a bunch of Julia internals that use them)
@inline Base.IndexStyle(x::MVTSeries) = IndexStyle(_vals(x))
@inline Base.dataids(x::MVTSeries) = Base.dataids(_vals(x))

# normally only the first of the following is sufficient.
# we add few other versions of similar below
Base.similar(x::MVTSeries) = MVTSeries(firstdate(x), colnames(x), similar(_vals(x)))

# -------------------------------------------------------------------------------
# Indexing with integers and booleans - same as matrices

# Indexing with integers falls back to AbstractArray
const _FallbackType = Union{Integer,Colon,AbstractUnitRange{<:Integer},AbstractArray{<:Integer},CartesianIndex}
Base.getindex(sd::MVTSeries, i1::_FallbackType...) = getindex(_vals(sd), i1...)
Base.setindex!(sd::MVTSeries, val, i1::_FallbackType...) = setindex!(_vals(sd), val, i1...)

# -------------------------------------------------------------
# Some other constructors
# -------------------------------------------------------------

# Empty (0 variables) from range
@inline MVTSeries(T::Type{<:Number}, rng::AbstractUnitRange{<:MIT}) = MVTSeries(first(rng), (), Matrix{T}(undef, length(rng), 0))
@inline MVTSeries(rng::UnitRange{<:MIT}) = MVTSeries(Float64, rng)

# Empty from a list of variables and of specified type (first date must also be given, Frequency is not enough)
# @inline MVTSeries(fd::MIT, vars) = MVTSeries(Float64, fd, vars)
MVTSeries(T::Type{<:Number}, fd::MIT, vars) = MVTSeries(fd, vars, Matrix{T}(undef, 0, length(vars)))

# Uninitialized from a range and list of variables
@inline MVTSeries(rng::UnitRange{<:MIT}, vars) = MVTSeries(Float64, rng, vars, undef)
@inline MVTSeries(rng::UnitRange{<:MIT}, vars, ::UndefInitializer) = MVTSeries(Float64, rng, vars, undef)
@inline MVTSeries(T::Type{<:Number}, rng::UnitRange{<:MIT}, vars) = MVTSeries(T, rng, vars, undef)
MVTSeries(T::Type{<:Number}, rng::UnitRange{<:MIT}, vars, ::UndefInitializer) =
    MVTSeries(first(rng), vars, Matrix{T}(undef, length(rng), length(vars)))
MVTSeries(T::Type{<:Number}, rng::UnitRange{<:MIT}, vars::Symbol, ::UndefInitializer) =
    MVTSeries(first(rng), (vars,), Matrix{T}(undef, length(rng), 1))

# initialize with a function like zeros, ones, rand.
MVTSeries(rng::UnitRange{<:MIT}, vars, init::Function) = MVTSeries(first(rng), vars, init(length(rng), length(vars)))
# no type-explicit version because the type is determined by the output of init()

#initialize with a constant
MVTSeries(rng::UnitRange{<:MIT}, vars, v::Number) = MVTSeries(first(rng), vars, fill(v, length(rng), length(vars)))

# construct with a given range (rather than only the first date). We must check the range length matches the data size 1
function MVTSeries(rng::UnitRange{<:MIT}, vars, vals::AbstractMatrix{<:Number})
    lrng = length(rng)
    lvrs = length(vars)
    nrow, ncol = size(vals)
    if lrng != nrow || lvrs != ncol
        throw(ArgumentError("Number of periods and variables do not match size of data." *
                            " ($(lrng)×$(lvrs)) ≠ ($(nrow)×$(ncol))"
        ))
    end
    return MVTSeries(first(rng), vars, vals)
end

# construct if data is given as a vector (it must be exactly 1 variable)
MVTSeries(fd::MIT, vars, data::AbstractVector) = MVTSeries(fd, vars, reshape(data, :, 1))
MVTSeries(fd::MIT, vars::Symbol, data::AbstractVector) = MVTSeries(fd, (vars,), reshape(data, :, 1))
MVTSeries(rng::UnitRange{<:MIT}, vars, data::AbstractVector) = MVTSeries(rng, vars, reshape(data, :, 1))
MVTSeries(rng::UnitRange{<:MIT}, vars::Symbol, data::AbstractVector) = MVTSeries(rng, (vars,), reshape(data, :, 1))

# construct uninitialized by way of calling similar
Base.similar(::Type{<:AbstractArray}, T::Type{<:Number}, shape::Tuple{UnitRange{<:MIT},NTuple{N,Symbol}}) where {N} = MVTSeries(T, shape[1], shape[2])
Base.similar(::Type{<:AbstractArray{T}}, shape::Tuple{UnitRange{<:MIT},NTuple{N,Symbol}}) where {T<:Number,N} = MVTSeries(T, shape[1], shape[2])
Base.similar(::AbstractArray, T::Type{<:Number}, shape::Tuple{UnitRange{<:MIT},NTuple{N,Symbol}}) where {N} = MVTSeries(T, shape[1], shape[2])
Base.similar(::AbstractArray{T}, shape::Tuple{UnitRange{<:MIT},NTuple{N,Symbol}}) where {T<:Number,N} = MVTSeries(T, shape[1], shape[2])

# construct from range and fill with the given constant or array
Base.fill(v::Number, rng::UnitRange{<:MIT}, vars::NTuple{N,Symbol}) where {N} = MVTSeries(first(rng), vars, fill(v, length(rng), length(vars)))

# -------------------------------------------------------------------------------
# Dot access to columns

Base.propertynames(x::MVTSeries) = tuple(colnames(x)...)

function Base.getproperty(x::MVTSeries, col::Symbol)
    col ∈ fieldnames(typeof(x)) && return getfield(x, col)
    columns = _cols(x)
    haskey(columns, col) && return columns[col]
    Base.throw_boundserror(x, [col,])
end

function Base.setproperty!(x::MVTSeries, name::Symbol, val)
    name ∈ fieldnames(typeof(x)) && return setfield!(x, name, val)
    columns = _cols(x)
    if !haskey(columns, name)
        error("Cannot append new column this way. " *
              "Use hcat(x; $name = value) or push!(x; $name = value).")
    end
    col = columns[name]
    ####  Do we need this mightalias check here? I think it's done in copyto!() so we should be safe without it.
    # if Base.mightalias(col, val)
    #     val = copy(val)
    # end
    if val isa TSeries
        rng = intersect(rangeof(col), rangeof(val))
        return copyto!(col, rng, val)
    elseif val isa Number
        return fill!(col.values, val)
    else
        return copyto!(col.values, val)
    end
end

# -------------------------------------------------------------------------------
# Indexing other than integers

# some check bounds that plug MVTSeries into the Julia infrastructure for AbstractArrays
@inline Base.checkbounds(::Type{Bool}, x::MVTSeries, p::MIT) = checkindex(Bool, rangeof(x), p)
@inline Base.checkbounds(::Type{Bool}, x::MVTSeries, p::UnitRange{<:MIT}) = checkindex(Bool, rangeof(x), p)
@inline Base.checkbounds(::Type{Bool}, x::MVTSeries, c::Symbol) = haskey(_cols(x), c)
@inline function Base.checkbounds(::Type{Bool}, x::MVTSeries, INDS::Union{Vector{Symbol},NTuple{N,Symbol}}) where {N}
    cols = _cols(x)
    for c in INDS
        haskey(cols, c) || return false
    end
    return true
end

@inline function Base.checkbounds(::Type{Bool}, x::MVTSeries, p::Union{MIT,UnitRange{<:MIT}}, c::Union{Symbol,Vector{Symbol},NTuple{N,Symbol}}) where {N}
    return checkbounds(Bool, x, p) && checkbounds(Bool, x, c)
end

# ---- single argument access

# single argument - MIT point - return the row as a vector (slice of .values)
Base.getindex(x::MVTSeries, p::MIT) = mixed_freq_error(x, p)
@inline function Base.getindex(x::MVTSeries{F}, p::MIT{F}) where {F<:Frequency}
    @boundscheck checkbounds(x, p)
    fi = firstindex(_vals(x), 1)
    getindex(x.values, fi + oftype(fi, p - firstdate(x)), :)
end

Base.setindex!(x::MVTSeries, val, p::MIT) = mixed_freq_error(x, p)
@inline function Base.setindex!(x::MVTSeries{F}, val, p::MIT{F}) where {F<:Frequency}
    @boundscheck checkbounds(x, p)
    fi = firstindex(_vals(x), 1)
    setindex!(_vals(x), val, fi + oftype(fi, p - firstdate(x)), :)
end

# single argument - MIT range
Base.getindex(x::MVTSeries, rng::UnitRange{MIT}) = mixed_freq_error(x, rng)
@inline function Base.getindex(x::MVTSeries{F}, rng::UnitRange{MIT{F}}) where {F<:Frequency}
    start, stop = _ind_range_check(x, rng)
    return MVTSeries(first(rng), axes(x, 2), getindex(_vals(x), start:stop, :))
end

Base.setindex!(x::MVTSeries, val, rng::UnitRange{MIT}) = mixed_freq_error(x, rng)
@inline function Base.setindex!(x::MVTSeries{F}, val, rng::UnitRange{MIT{F}}) where {F<:Frequency}
    start, stop = _ind_range_check(x, rng)
    setindex!(_vals(x), val, start:stop, :)
end

# single argument - variable - return a TSeries of the column
@inline Base.getindex(x::MVTSeries, col::AbstractString) = getindex(x, Symbol(col))
@inline function Base.getindex(x::MVTSeries, col::Symbol)
    @boundscheck checkbounds(x, col)
    getproperty(x, col)
end

@inline Base.setindex!(x::MVTSeries, val, col::AbstractString) = setindex!(x, val, Symbol(col))
@inline function Base.setindex!(x::MVTSeries, val, col::Symbol)
    setproperty!(x, col, val)
end

# single argument - list/tuple of variables - return a TSeries of the column
@inline function Base.getindex(x::MVTSeries, col::Union{Vector{Symbol},NTuple{N,Symbol}}) where {N}
    @boundscheck checkbounds(x, col)
    names = [colnames(x)...]
    inds = indexin(col, names)
    return MVTSeries(firstdate(x), tuple(names[inds]...), getindex(_vals(x), :, inds))
end

@inline function Base.setindex!(x::MVTSeries, val, col::Union{Vector{Symbol},NTuple{N,Symbol}}) where {N}
    @boundscheck checkbounds(x, col)
    names = [colnames(x)...]
    inds = indexin(col, names)
    setindex!(x.values, val, :, inds)
end

# ---- two arguments indexing

const _SymbolOneOrCollection = Union{Symbol,Vector{Symbol},NTuple{N,Symbol}} where {N}
const _MITOneOrRange = Union{MIT,UnitRange{<:MIT}}

@inline Base.getindex(x::MVTSeries, p::_MITOneOrRange, c::_SymbolOneOrCollection) = mixed_freq_error(x, p)
@inline Base.setindex!(x::MVTSeries, val, p::_MITOneOrRange, c::_SymbolOneOrCollection) = mixed_freq_error(x, p)

# if one argument is Colon, fall back on single argument indexing
@inline Base.getindex(x::MVTSeries, p::_MITOneOrRange, ::Colon) = getindex(x, p)
@inline Base.getindex(x::MVTSeries, ::Colon, c::_SymbolOneOrCollection) = getindex(x, c)

@inline Base.setindex!(x::MVTSeries, val, p::_MITOneOrRange, ::Colon) = setindex!(x, val, p, axes(x, 2))
@inline Base.setindex!(x::MVTSeries, val, ::Colon, c::_SymbolOneOrCollection) = setindex!(x, val, axes(x, 1), c)

#

_colind(x, c::Symbol) = indexin([Symbol(c)], [colnames(x)...])[1]
_colind(x, c::Tuple) = indexin(Symbol.(c), [colnames(x)...])
_colind(x, c::Vector) = indexin(Symbol.(c), [colnames(x)...])

# with a single MIT and single Symbol we return a number
# with a single MIT and multiple Symbol-s we return a Vector
# the appropriate dispatch is done in getindex on the values, so we wrap both cases in a single function
@inline function Base.getindex(x::MVTSeries{F}, p::MIT{F}, c::_SymbolOneOrCollection) where {F<:Frequency}
    @boundscheck checkbounds(x, c)
    @boundscheck checkbounds(x, p)
    fi = firstindex(_vals(x), 1)
    i1 = oftype(fi, fi + (p - firstdate(x)))
    i2 = _colind(x, c)
    getindex(x.values, i1, i2)
end

# with an MIT range and a Symbol (single column) we return a TSeries
@inline function Base.getindex(x::MVTSeries{F}, p::UnitRange{MIT{F}}, c::Symbol) where {F<:Frequency}
    @boundscheck checkbounds(x, c)
    @boundscheck checkbounds(x, p)
    start, stop = _ind_range_check(x, p)
    i1 = start:stop
    i2 = _colind(x, c)
    return TSeries(first(p), getindex(_vals(x), i1, i2))
end

# with an MIT range and a sequence of Symbol-s we return an MVTSeries
@inline function Base.getindex(x::MVTSeries{F}, p::UnitRange{MIT{F}}, c::Union{NTuple{N,Symbol},Vector{Symbol}}) where {F<:Frequency,N}
    @boundscheck checkbounds(x, c)
    @boundscheck checkbounds(x, p)
    start, stop = _ind_range_check(x, p)
    i1 = start:stop
    i2 = _colind(x, c)
    return MVTSeries(first(p), axes(x, 2)[i2], getindex(_vals(x), i1, i2))
end

# assignments

# with a single MIT we assign a number or a row-Vector
@inline function Base.setindex!(x::MVTSeries{F}, val, p::MIT{F}, c::_SymbolOneOrCollection) where {F<:Frequency}
    @boundscheck checkbounds(x, p)
    @boundscheck checkbounds(x, c)
    fi = firstindex(_vals(x), 1)
    i1 = oftype(fi, fi + (p - firstdate(x)))
    i2 = _colind(x, c)
    setindex!(x.values, val, i1, i2)
end

# with a range of MIT and a single column - we fall back on TSeries assignment
@inline function Base.setindex!(x::MVTSeries{F}, val, r::UnitRange{MIT{F}}, c::Symbol) where {F<:Frequency}
    col = getindex(x, c)
    setindex!(col, val, r)
end

@inline function Base.setindex!(x::MVTSeries{F}, val, r::UnitRange{MIT{F}}, c::Union{Vector{Symbol},NTuple{N,Symbol}}) where {F<:Frequency,N}
    @boundscheck checkbounds(x, r)
    @boundscheck checkbounds(x, c)
    start, stop = _ind_range_check(x, r)
    i1 = start:stop
    i2 = _colind(x, c)
    setindex!(_vals(x), val, i1, i2)
end

@inline function Base.setindex!(x::MVTSeries{F}, val::MVTSeries{F}, r::UnitRange{MIT{F}}, c::Union{Vector{Symbol},NTuple{N,Symbol}}) where {F<:Frequency,N}
    @boundscheck checkbounds(x, r)
    @boundscheck checkbounds(x, c)
    @boundscheck checkbounds(val, r)
    @boundscheck checkbounds(val, c)
    xcols = _cols(x)
    vcols = _cols(val)
    for col in c
        setindex!(xcols[col], vcols[col], r)
    end
    return x
end

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

Base.copyto!(dest::MVTSeries, src::AbstractArray) = (copyto!(dest.values, src); dest)
Base.copyto!(dest::MVTSeries, src::MVTSeries) = (copyto!(dest.values, src.values); dest)

# -------------------------------------------------------------------------------
# ways add new columns (variables)

function Base.hcat(x::MVTSeries; KW...)
    y = MVTSeries(rangeof(x), tuple(colnames(x)..., keys(KW)...))
    # copyto!(y, x)
    for (k, v) in pairs(x)
        setproperty!(y, k, v)
    end
    for (k, v) in KW
        setproperty!(y, k, v)
    end
    return y
end

####   Views

Base.fill!(x::MVTSeries, val) = fill!(_vals(x), val)

# @inline Base.view(x::MVTSeries, I::Union{MIT,AbstractRange{<:MIT}}, J) = mixed_freq_error(x, eltype(I))

@inline Base.view(x::MVTSeries, I, J) = view(_vals(x), I, J)

@inline Base.view(x::MVTSeries, ::Colon, J::_SymbolOneOrCollection) = view(x, axes(x, 1), J)
@inline Base.view(x::MVTSeries, I::_MITOneOrRange, ::Colon) = view(x, I, axes(x, 2))
@inline Base.view(x::MVTSeries, ::Colon, ::Colon) = view(x, axes(x, 1), axes(x, 2))
function Base.view(x::MVTSeries, I::_MITOneOrRange, J::_SymbolOneOrCollection) where {F<:Frequency}
    @boundscheck checkbounds(x, I)
    @boundscheck checkbounds(x, J)
    start, stop = _ind_range_check(x, I)
    i1 = start:stop
    i2 = _colind(x, J)
    return MVTSeries(first(I), axes(x, 2)[i2], view(_vals(x), i1, i2))
end

####

include("mvtseries/mvts_broadcast.jl")
include("mvtseries/mvts_show.jl")

####  arraymath

@inline Base.promote_shape(x::MVTSeries, y::MVTSeries) =
    (intersect(rangeof(x),rangeof(y)),
    sym_common_axes(axes(x,2), axes(y,2)))

@inline Base.promote_shape(x::MVTSeries, y::AbstractMatrix) =
    promote_shape(_vals(x), y)

@inline Base.LinearIndices(x::MVTSeries) = LinearIndices(_vals(x))


text:
import Term: remove_markup, remove_ansi, truncate, reshape_text
import Term.style: apply_style

@testset "\e[31mTU - remove" begin
    # ------------------------------- remove markup ------------------------------ #
    @test remove_markup("[red]text[/red]") == "text"
    @test remove_markup("[red]sffsdf[green on_blue italic]") == "sffsdf"
    @test remove_markup("[red]aa[green on_blue italic]bb[/green on_blue italic]cc[/]") ==
          "aabbcc"
    @test remove_markup("[(255, 1, 3)]aa[#ff0022 italic]bb[/ff0022]") == "aabb"

    # test that parentheses are ignored
    strings = [
        "1111",
        "[222a",
        "333asd[adsada",
        "a444asdas]asda]asda",
        "[a5555adsad[[sada",
        "a6666[[asdasd]]",
    ]
    for str in strings
        @test remove_markup(str) == str
    end

    # -------------------------------- remove ansi ------------------------------- #
    strings = [
        ("[red]asdasdad[/red]", "asdasdad"),
        ("aaa[bold on_green]bb[/]", "aaabb"),
        ("[red]aa[bold]bb[#ffffff]cc[/#ffffff]dd", "aabbccdd"),
        ("[bold italic underline]aa[/bold italic underline]", "aa"),
    ]
    for (str, clean) in strings
        @test remove_ansi(apply_style(str)) == clean
        @test remove_markup(str) == clean
    end
end

@testset "\e[31mTU - reshape text" begin
    nlines(x) = length(split(x, "\n"))
    lw(x) = length(split(x, "\n")[1])

    # --------------------------------- truncate --------------------------------- #
    strings = [("a"^20, 6), ("asd"^33, 12), ("c"^3, 22)]
    for (str, w) in strings
        @test length(truncate(str, w)) <= w
    end

    # ---------------------------------- reshape --------------------------------- #
    s1 = "."^500

    sizes = [(50, 10), (25, 20), (5, 100), (100, 5)]

    for (w, h) in sizes
        rs = reshape_text(s1, w)
        @test lw(rs) == w
        @test nlines(rs) == h
    end

    s1 = "[red]aa[green]bb[/green]c[/red]"^100

    sizes = [(50, 10), (25, 20), (5, 100), (100, 5)]

    for (w, h) in sizes
        rs = remove_markup(reshape_text(s1, w))
        @test lw(rs) == w
        @test nlines(rs) == h
    end

end


text:
function parseinfix2rpn(s)
    outputq = []
    opstack = []
    infix = split(s)
    for tok in infix
        if all(isnumber, tok)
            push!(outputq, tok)
        elseif tok == "("
            push!(opstack, tok)
        elseif tok == ")"
            while !isempty(opstack) && (op = pop!(opstack)) != "("
               push!(outputq, op)
            end
        else # operator
            while !isempty(opstack)
                op = pop!(opstack)
                if Base.operator_precedence(Symbol(op)) > Base.operator_precedence(Symbol(tok)) ||
                   (Base.operator_precedence(Symbol(op)) ==
                     Base.operator_precedence(Symbol(tok)) && op != "^")
                    push!(outputq, op)
                else
                    push!(opstack, op)  # undo peek
                    break
                end
            end
            push!(opstack, tok)
        end
        println("The working output stack is $outputq")
    end
    while !isempty(opstack)
        if (op = pop!(opstack)) == "("
            throw("mismatched parentheses")
        else
            push!(outputq, op)
        end
    end
    outputq
end

teststring = "3 + 4 * 2 / ( 1 - 5 ) ^ 2 ^ 3"
println("\nResult: $teststring becomes $(join(parseinfix2rpn(teststring), ' '))")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function LidDrivenCavity(;
                         tfinal = Inf,
                         Lx = 1, Ly = 1, CFL = 0.5, Re = 100,
                         Nx = 65, Ny = 65,
                         u_n = 1, u_s = 0, v_w = 0, v_e = 0,
                         printfreq = 10)
    t0 = time()
    println("------------------Ny = $(Ny), Nx = $(Nx) ---------------")
    Δy  = Ly/(Ny-1)
    Δx  = Lx/(Nx-1)
    x = 0:Δx:Lx
    y = 0:Δy:Ly
    Δt = CFL * Δx

    # Construct matrix for Poisson equation
    A_poisson = BuildPoissonMatrix(Ny,Nx,Δx,Δy) # for coNxgrad
    # Construct matrix for advection-diffusion equation
    ap,an,as,ae,aw = BuildAdvectionDiffusionCoefficients(Re,Δt,Δx,Δy)
    # allocate empty matrices for Gauss-Siedel solver
    Rp = zeros(Ny,Nx); res = zeros(Ny,Nx)

    # initialize ω and ψ
    ω = zeros(Ny,Nx)
    ψ = zeros(Ny,Nx)

    # keep track of changes with Δ's
    ω_old = zeros(Ny,Nx)
    ψ_old = zeros(Ny,Nx)
    ω_hist = []
    ψ_hist = []
    residual = 1

    ###########################
    k0,t = 0,0
    while t < tfinal && maximum(residual) > 1e-8
        t += Δt
        k0 += 1
        # Solve Poisson equation for ψ:
        LinearSolve!(A_poisson,ψ,-ω)

        VorticityBoundaryConditions!(ω,ψ,Δx,Δy,u_n,u_s,v_e,v_w)
        # Modify the explicit part of advection-diffusion equation
        BuildAdvectionDiffusionRHS!(Rp,ω,ψ,Δt,Δx,Δy,Ny,Nx,Re)

        # Solve advection-diffusion equation for ω:
        GaussSiedel!(ω,ap,an,as,ae,aw,Rp,res)

        residual = RecordHistory!.([ω,ψ],[ω_old,ψ_old],[ω_hist,ψ_hist])

        if (k0 % printfreq == 0)
            println("Step: $k0 \t Time: $(round(t,digits=3))\t",
                    "|Δω|: $(round((residual[1]),digits=8)) \t",
                    "|Δψ|: $(round((residual[2]),digits=8)) \t")
        end
    end
    tt = round(time() - t0, digits=3)
    println("This took $tt seconds.")
    println("--------------------------------------------------------")
    Results(ψ,ω,hcat(ω_hist,ψ_hist),x,y,t,k0,Re)
end


text:
struct TreeView{T, N, Tnode <: AbstractNode{<: Any, N}} <: AbstractArray{T, N}
    node::Tnode
end

function TreeView(rootnode::AbstractNode{<: Any, N}) where {N}
    T = leafeltype(rootnode)
    TreeView{T, N, typeof(rootnode)}(rootnode)
end

rootnode(x::TreeView) = x.node
Base.size(x::TreeView) = totalsize(rootnode(x))
leaftype(x::TreeView) = leaftype(rootnode(x))
leafeltype(x::TreeView) = leafeltype(rootnode(x))
nleafblocks(x::TreeView) = nleafblocks(rootnode(x))
leafblockunit(x::TreeView) = leafblockunit(rootnode(x))

@inline TreeLinearIndex(x::TreeView{<: Any, N}, I::Vararg{Integer, N}) where {N} = TreeLinearIndex(rootnode(x), I...)
@inline TreeCartesianIndex(x::TreeView{<: Any, N}, I::Vararg{Integer, N}) where {N} = TreeCartesianIndex(rootnode(x), I...)

@inline function Base.getindex(x::TreeView{<: Any, N}, I::Vararg{Int, N}) where {N}
    @boundscheck checkbounds(x, I...)
    index = TreeLinearIndex(x, I...)
    @inbounds x[index]
end

@inline function Base.setindex!(x::TreeView{<: Any, N}, v, I::Vararg{Int, N}) where {N}
    @boundscheck checkbounds(x, I...)
    index = TreeLinearIndex(x, I...)
    @inbounds x[index] = v
end

for f in (:isactive, :allocate!, :isallocated)
    @eval @inline function $f(x::TreeView{<: Any, N}, I::Vararg{Int, N}) where {N}
        @boundscheck checkbounds(x, I...)
        index = TreeLinearIndex(x, I...)
        @inbounds $f(x, index)
    end
end

wrap_treeview(x) = x
wrap_treeview(x::AbstractNode) = TreeView(x)
@generated function Base.getindex(x::TreeView, I::TreeIndex{depth}) where {depth}
    ex = :(rootnode(x))
    for i in 1:depth
        ex = :($ex[I[$i]])
    end
    quote
        @_inline_meta
        @_propagate_inbounds_meta
        wrap_treeview($ex)
    end
end

@inline function Base.setindex!(x::TreeView, v, I::TreeIndex)
    @_propagate_inbounds_meta
    _setindex!_getleaf(x, v, I)
    x
end

@generated function Base.setindex!(x::TreeView, ::Nothing, I::TreeIndex{depth}) where {depth}
    ex = :(rootnode(x))
    for i in 1:depth
        sym = Symbol(:node, i)
        ex = quote
            $sym = $ex
            isactive($sym, I[$i]) || return x
            $(i == depth ? :($sym[I[$i]] = nothing) : :(unsafe_getindex($sym, I[$i])))
        end
    end
    quote
        @_inline_meta
        @_propagate_inbounds_meta
        $ex
    end
end

@generated function isactive(x::TreeView, I::TreeIndex{depth}) where {depth}
    ex = :(rootnode(x))
    for i in 1:depth
        sym = Symbol(:node, i)
        ex = quote
            $sym = $ex
            isactive($sym, I[$i]) || return false
            $(i == depth ? true : :(unsafe_getindex($sym, I[$i])))
        end
    end
    quote
        @_inline_meta
        @_propagate_inbounds_meta
        $ex
    end
end

# used in ContinuousView
@inline function _setindex!_getleaf(x::TreeView{<: Any, N}, v, I::Vararg{Int, N}) where {N}
    @boundscheck checkbounds(x, I...)
    index = TreeLinearIndex(x, I...)
    @inbounds _setindex!_getleaf(x, v, index)
end

@generated function _setindex!_getleaf(x::TreeView, v, I::TreeIndex{depth}) where {depth}
    ex = :(rootnode(x))
    for i in 1:depth
        if i == depth
            ex = :(setindex!($ex, v, I[$i]))
        else
            ex = :(allocate!($ex, I[$i]))
        end
    end
    quote
        @_inline_meta
        @_propagate_inbounds_meta
        $ex
    end
end

@generated function allocate!(x::TreeView, I::TreeIndex{depth}) where {depth}
    ex = :(rootnode(x))
    for i in 1:depth
        ex = :(allocate!($ex, I[$i]))
    end
    quote
        @_inline_meta
        @_propagate_inbounds_meta
        $ex
    end
end

@generated function isallocated(x::TreeView, I::TreeIndex{depth}) where {depth}
    exps = Expr[:(node = rootnode(x))]
    for i in 1:depth
        ex = quote
            isallocated(node, I[$i]) || return false
            node = unsafe_getindex(node, I[$i])
        end
        push!(exps, ex)
    end
    quote
        @_inline_meta
        @_propagate_inbounds_meta
        $(exps...)
        true
    end
end

cleanup!(x::TreeView) = (cleanup!(rootnode(x)); x)

function Base.fill!(x::TreeView, ::Nothing)
    node = rootnode(x)
    isnull(node) && return x
    if length(x) > THREADS_THRESHOLD
        fill!(getmask(node), false)
        Threads.@threads for i in eachindex(node)
            @inbounds deactivate!(node[i])
        end
    else
        deactivate!(node)
    end
    x
end

function Base.fill!(A::SubArray{<: Any, <: Any, <: TreeView}, ::Nothing)
    for I in eachindex(A)
        @inbounds A[I] = nothing
    end
    A
end

function nleaves(x::TreeView)
    node = rootnode(x)
    if length(x) > THREADS_THRESHOLD
        isnull(node) && return 0
        counts = zeros(Int, Threads.nthreads())
        @inbounds Threads.@threads for i in eachindex(node)
            if isactive(node, i)
                child = unsafe_getindex(node, i)
                counts[Threads.threadid()] += nleaves(child)
            end
        end
        sum(counts)
    else
        nleaves(node)
    end
end


text:
# container for nongenetic covariates
immutable CovariateMatrix{T <: Float} <: AbstractArray{T, 2}
    x  :: SharedMatrix{T}
    p  :: Int
    xt :: SharedMatrix{T}
    h  :: Vector{String} # header

    #CovariateMatrix(x::SharedMatrix{T}, p::Int, xt::SharedMatrix{T}, h::Vector{String}) = new(x,p,xt,h)
    CovariateMatrix{T}(x, p, xt, h) where {T <: Float} = new{T}(x,p,xt,h)
end

### 22 Sep 2016: not needed in Julia v0.5?
#function CovariateMatrix{T <: Float}(
#    x  :: SharedMatrix{T},
#    p  :: Int,
#    xt :: SharedMatrix{T},
#    h  :: Vector{String}
#)
#    CovariateMatrix{eltype(x)}(x, p, xt, h) :: CovariateMatrix{eltype(x)}
#end

function CovariateMatrix(
    T        :: Type,
    filename :: String;
    pids     :: Vector{Int} = procs(),
    header   :: Bool = false
)
    # first load data
    if header
        x_temp, h = readdlm(filename, T, header=true) :: Matrix{T}, Matrix{AbstractString}
        h = convert(Vector{String}, h) :: Vector{String}
    else
        x_temp = readdlm(filename, T, header=false) :: Matrix{T}
        h = ["V" * "$i" for i = 1:size(x_temp,2)] :: Vector{String}
    end

    # make SharedArray container for x_temp
    x = SharedArray{T}(size(x_temp), pids=pids) :: SharedMatrix{eltype(x_temp)}
    copy!(x, x_temp)

    # do same for x'
    xt = SharedArray{T}(reverse(size(x_temp)), pids=pids) :: SharedMatrix{eltype(x_temp)}
    transpose!(xt, x_temp)

    # save p for convenience
    p = size(x,2)

    return CovariateMatrix{eltype(x)}(x, p, xt, h) :: CovariateMatrix{eltype(x)}
end

# default to Float64 constructor
function CovariateMatrix(
    filename :: String;
    pids     :: Vector{Int} = procs(),
    header   :: Bool = false
)
    CovariateMatrix(Float64, filename, pids=pids, header=header) :: CovariateMatrix{Float64}
end

# subroutines to define AbstractArray
Base.size(x::CovariateMatrix)         = size(x.x)
Base.size(x::CovariateMatrix, d::Int) = size(x.x, d)

Base.length(x::CovariateMatrix) = prod(size(x.x))

Base.ndims(x::CovariateMatrix) = ndims(x.x)

Base.endof(x::CovariateMatrix) = length(x.x)

Base.eltype(x::CovariateMatrix) = eltype(x.x)

Base.IndexStyle(::Type{CovariateMatrix}) = Base.IndexLinear()

# assume that SharedArray processes are from principal covariate matrix
Base.procs(x::CovariateMatrix) = procs(x.x)

# here we annotate return value with type of x
# this hack is **VERY** important for efficient indexing x from a BEDFile!
# experience is that compiler does not know return value of getindex(x::CovariateMatrix{T},...)!
# consequently, return value of getindex(::BEDFile,...) is **Any**
# annotating type here changes return value of getindex(::BEDFile,...) to T <: Float !!!
Base.getindex{T <: Float}(x::CovariateMatrix{T}, i::Int) = getindex(x.x, i)::T
Base.getindex{T <: Float}(x::CovariateMatrix{T}, i::Int, j::Int) = getindex(x.x, i, j)::T

# fix this to change index for both x, xt?
#function Base.setindex!(x::CovariateMatrix, i::Int)
#  setindex!(x.x, i), setindex!(x.xt, i)
#end

function Base.setindex!(x::CovariateMatrix, i::Int, j::Int)
    setindex!(x.x, i, j)
    setindex!(x.xt, j, i)
end

function Base.similar(x::CovariateMatrix, T::Real, dims::Dims)
    CovariateMatrix(zeros(T, dims), dims[2], zeros(T, reverse(dims)), ["V" * "$i" for i = 1:dims[2]])
end # function Base.similar

function Base.similar(x::CovariateMatrix)
    n,p = size(x)
    T   = eltype(x)
    CovariateMatrix(SharedArray{T}((n,p), pids=procs(x)), p, SharedArray{T}((p,n), pids=procs(x)), ["V" * "$i" for i = 1:p])
end # function Base.similar

display(x::CovariateMatrix) = display(x.x)

max_stars_repo_name: klkeys/PLINK

max_issues_repo_name: klkeys/PLINK


text:
#
#
"""
A piecewise_linear_generator. Used for generating budget constraints used in
conventional microeconomics. See Stark and Duncan 1992 in the biblio.
"""

struct Point2DG{T<:Real}
      x::T
      y::T
end

const Point2D = Point2DG{Float64}

struct Line2DG{T<:Real}
      a::T
      b::T
end

const Line2D = Line2DG{Float64}

const BudgetConstraint = Array{Point2DG,1}
const PointsSet = Set{Point2DG}

const VERTICAL   = 9999999999.9999;
const TOLERANCE  = 0.0001;
const INCREMENT  = 0.0001;
const MAX_DEPTH  = 50;
const MAX_INCOME = 20000.0;
const MIN_INCOME = 0.0;

const ROUND_OUTPUT = false;

struct BCSettings
    mingross :: Float64
    maxgross :: Float64
    increment :: Float64
    tolerance :: Float64
    round_output :: Bool
    maxdepth :: Integer
end

const DEFAULT_SETTINGS = BCSettings( MIN_INCOME, MAX_INCOME, INCREMENT, TOLERANCE, true, MAX_DEPTH )

function makeline( point_1 :: Point2D, point_2 :: Point2D )::Line2D
    a :: Float64 = 0.0
    b :: Float64 = 0.0
    if point_1.x == point_2.x # ?? never taken??
        b = point_1.x;
        a = VERTICAL;
    else
        b = (point_1.y - point_2.y)/(point_1.x - point_2.x );
        b = min( b, VERTICAL );
        a = ( point_1.y - point_1.x*b );
    end
    return Line2D( a, b );
end

function findintersection( line_1::Line2D, line_2 :: Line2D ) :: Point2D
    x :: Float64 = 0.0
    y :: Float64 = 0.0
    if !( line_1.b ≈ line_2.b )
        x = (line_2.a - line_1.a) / (line_1.b - line_2.b);
        y = line_1.a + (x * line_1.b );
    else
        x = 0.0;
        y = line_1.a;
    end
    return Point2D( x, y );
end

function comparepoints( point_1 :: Point2DG, point_2 :: Point2DG ) :: Integer
    if( point_1.x > point_2.x )
        return 1;
    end
    if( point_1.x < point_2.x )
        return -1;
    end
    if( point_1.y > point_2.y )
        return 1;
    end
    if( point_1.y < point_2.y )
        return -1;
    end
    return 0;
end

# sorting for points
import Base.isless

function isless( point_1 :: Point2DG, point_2 :: Point2DG ) :: Bool
    return comparepoints( point_1, point_2 ) < 0
end

function marginalrate( point_1 :: Point2D, point_2 :: Point2D )::Float64
    mr :: Float64 = 0.0
    if !( point_2.x ≈ point_1.x )
        mr = 100.0 * (1-(point_2.y-point_1.y) / (point_2.x - point_1.x));
    else
        mr = VERTICAL;
        if( point_2.y < point_1.y )
            mr *= -1;
        end
    end
    return mr;
end

# approximate equality for points and lines
import Base.≈

function ≈(left :: Point2D, right::Point2D )::Bool
   (left.x ≈ right.x) && ( left.y ≈ right.y )
end

function ≈(left :: Line2D, right::Line2D )::Bool
   # return (left.a ≈ right.a) && ( left.b ≈ right.b )
   #  FIXME this uses the global tolerance even
   # if the user overrides it in settings
   return ((( abs(left.a-right.a)) <= TOLERANCE ) && (( abs(left.b-right.b)) <= TOLERANCE ));
end

# round a float to 2dps
function round2pl( x::Float64 )::Float64
    x *= 100.0
    i = trunc(x)
    x = Float64(i/100)
end

function round2pl!( bc :: BudgetConstraint )
    nbc = size( bc )[1]
    for i in 1:nbc
        p = Point2D( round2pl( bc[i].x ), round2pl( bc[i].y))
        bc[i] = p
    end
end

function toarray( ps :: PointsSet ) :: BudgetConstraint
    bc = BudgetConstraint();
    for p in ps
        push!( bc, p )
    end
    sort!( bc ) # , lt=isless
    bc
end

function censor( ps :: PointsSet, round :: Bool=true ) :: BudgetConstraint
    # println( "censor; ps=$ps")
    bc = toarray( ps )
    sort!( bc )
    nbc = size( bc )[1]
    first = bc[1]
    last = bc[nbc]
    if( nbc < 3 )
        return bc
    end
    i = 2
    while i < (nbc - 1)
        p1 = bc[i-1]
        p2 = bc[i]
        p3 = bc[i+1]
        l1 = makeline( p1, p2 )
        l2 = makeline( p2, p3 )
        if l1 ≈ l2
            deleteat!( bc, i )
            nbc -= 1
        else
            i += 1
        end
    end # while
    i = 1
    if round
        round2pl!( bc )
    end
    # 1 liner which sorts and removes dups
    bc = BudgetConstraint( toarray( PointsSet( bc )))
    if bc[1] != first
        # bc = vcat(first,bc)
    end
    if bc[nbc] != last
        # bc = vcat(bc,last)
    end
    return bc
end

function generate!(
    bc       :: PointsSet,
    getnet,  # TODO give this a signature ( Float64 ) :: Float64
    depth    :: Integer,
    startpos :: Float64,
    endpos   :: Float64,
    settings :: BCSettings ) :: Integer
    diff = abs( startpos - endpos )
    tolerance = settings.tolerance
    if( diff < settings.tolerance )
        return depth
    end
    if depth > settings.maxdepth
        throw( "max depth exceeded $depth"  )
    end
    p1 = Point2D( startpos, getnet(startpos) )
    p2 = Point2D( startpos+settings.increment, getnet(startpos+settings.increment))
    p4 = Point2D( endpos, getnet(endpos) )
    p3 = Point2D( endpos-settings.increment, getnet(endpos-settings.increment) )

    line1 = makeline( p1, p2 )
    line2 = makeline( p3, p4 )

    if line1 ≈ line2
        push!( bc, p1 )
        push!( bc, p4 )

        return depth
    end
    p5 = findintersection( line1, line2 )
    if( p5.x <= startpos ) || ( p5.x >= endpos )
        anchor = startpos + (( endpos - startpos)/2.0)
    else
        anchor = p5.x
    end
    depth += 1
    #
    # expore to the left
    #
    depth = generate!( bc, getnet, depth, startpos, anchor, settings )
    #
    # then the right
    #
    depth = generate!( bc, getnet, depth, anchor, endpos, settings )
    return depth - 1
end

"""
Make a budget constraint using function `getnet` to extract net incomes and `settings` (see above on this struct).
getnet should be a function of the form `net=f(gross)`. See the testcase for an example.
"""
function makebc( getnet, settings :: BCSettings = DEFAULT_SETTINGS ) :: BudgetConstraint
    bc = BudgetConstraint()
    ps = PointsSet()
    depth = 0
    try
        depth = generate!( ps, getnet, depth, settings.mingross, settings.maxgross, settings )
        bc = censor( ps, settings.round_output )
    catch e
        ## FIXME print a fuller stack trace here
        println( "failed! $e")
        println(stacktrace())
    end
    bc;
end

function pointstoarray( bc :: BudgetConstraint ) :: Array{Float64,2}
    sz = size( bc, 1 )
    pts = zeros( Float64, sz, 2 )
    for b in 1:sz
        pts[b,1] = bc[b].x
        pts[b,2] = bc[b].y
    end
    pts
end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

find a polynomial f such that f² >= Σᵢ gᵢ² or f >= Σᵢ |gᵢ|
where gᵢ are arbitrary polynomials, and the volume under f is minimized
=#

struct PolyNormJuMP{T <: Real} <: ExampleInstanceJuMP{T}
    n::Int
    rand_halfdeg::Int # maximum degree of randomly generated polynomials
    epi_halfdeg::Int # half of max degree of cone polynomials
    num_polys::Int
    use_l1::Bool # use epigraph of one norm, otherwise Euclidean norm
    use_norm_cone::Bool # use Euclidean / L1 norm WSOS cones
    use_wsos_scalar::Bool # use WSOS scalar cone
end

function build(inst::PolyNormJuMP{T}) where {T <: Float64}
    (n, num_polys, epi_halfdeg, rand_halfdeg) =
        (inst.n, inst.num_polys, inst.epi_halfdeg, inst.rand_halfdeg)
    @assert epi_halfdeg >= rand_halfdeg

    dom = PolyUtils.BoxDomain{T}(-ones(T, n), ones(T, n))
    rand_U = binomial(n + 2 * rand_halfdeg, n)
    rand_coeffs = rand(-9:9, rand_U, num_polys)
    (U, pts, Ps, V, w) = PolyUtils.interpolate(dom, epi_halfdeg,
        calc_V = true, get_quadr = true)
    polys = V[:, 1:rand_U] * rand_coeffs

    model = JuMP.Model()
    JuMP.@variable(model, f[1:U])
    JuMP.@objective(model, Min, dot(w, f))

    if inst.use_norm_cone
        cone = (inst.use_l1 ? Hypatia.WSOSInterpEpiNormOneCone :
            Hypatia.WSOSInterpEpiNormEuclCone)
        JuMP.@constraint(model, vcat(f, vec(polys)) in
            cone{T}(num_polys + 1, U, Ps))
    elseif inst.use_l1
        lhs = one(T) * f
        for i in 1:inst.num_polys
            p_plus_i = JuMP.@variable(model, [1:U])
            p_minus_i = JuMP.@variable(model, [1:U])
            lhs .-= p_plus_i
            lhs .-= p_minus_i
            JuMP.@constraints(model, begin
                polys[:, i] .== p_plus_i - p_minus_i
                p_plus_i in Hypatia.WSOSInterpNonnegativeCone{T, T}(U, Ps)
                p_minus_i in Hypatia.WSOSInterpNonnegativeCone{T, T}(U, Ps)
            end)
        end
        JuMP.@constraint(model, lhs in
            Hypatia.WSOSInterpNonnegativeCone{T, T}(U, Ps))
    else
        R = num_polys + 1
        svec_dim = Cones.svec_length(R)
        polyvec = zeros(JuMP.AffExpr, svec_dim * U)
        polyvec[1:U] .= f

        if !inst.use_wsos_scalar
            idx = 2
            rt2 = sqrt(T(2))
            for i in 1:inst.num_polys
                polyvec[Cones.block_idxs(U, idx)] = rt2 * polys[:, i]
                idx += i
                polyvec[Cones.block_idxs(U, idx)] = f
                idx += 1
            end
            cone = Hypatia.WSOSInterpPosSemidefTriCone{T}(R, U, Ps)
        else
            ypts = zeros(svec_dim, R)
            ypts[1, 1] = 1
            idx = 2
            for j in 2:R
                polyvec[Cones.block_idxs(U, idx)] = 2 * (polys[:, j - 1] + f)
                ypts[idx, 1] = ypts[idx, j] = 1
                idx += 1
                for i in 2:(j - 1)
                    polyvec[Cones.block_idxs(U, idx)] = 2 * f
                    ypts[idx, i] = ypts[idx, j] = 1
                    idx += 1
                end
                polyvec[Cones.block_idxs(U, idx)] = f
                ypts[idx, j] = 1
                idx += 1
            end
            new_Ps = Matrix{T}[]
            for P in Ps
                push!(new_Ps, kron(ypts, P))
            end
            cone = Hypatia.WSOSInterpNonnegativeCone{T, T}(U * svec_dim, new_Ps)
        end
        JuMP.@constraint(model, polyvec in cone)
    end

    return model
end


text:
export  E_EXPECTED_EOF, E_UNEXPECTED_EOF, E_UNEXPECTED_CHAR, E_BAD_KEY,
        E_BAD_ESCAPE, E_BAD_CONTROL, E_LEADING_ZERO, E_BAD_NUMBER, E_BAD_HEADER,
        E_BAD_TRAILER, E_BAD_TYPE, E_NOT_IMPLEMENTED,
        E_INVALID_OBJECT, E_INVALID_PAGE_NUMBER, E_INVALID_PAGE_LABEL,
        E_NOT_TAGGED_PDF, E_INVALID_PASSWORD

# The following errors may be thrown by the reader
const E_EXPECTED_EOF    = "Expected end of input"
const E_UNEXPECTED_EOF  = "Unexpected end of input"
const E_UNEXPECTED_CHAR = "Unexpected character"
const E_BAD_KEY         = "Invalid object key"
const E_BAD_ESCAPE      = "Invalid escape sequence"
const E_BAD_CONTROL     = "ASCII control character in string"
const E_BAD_NUMBER      = "Invalid number"
const E_BAD_HEADER      = "Invalid file header"
const E_BAD_TRAILER     = "Invalid file trailer"
const E_BAD_TYPE        = "Invalid object type"

const E_FAILED_COMPRESSION  = "Data not compressed properly."
const E_DECRYPT_DOCUMENT    = "Unable to decrypt document"
const E_INVALID_CATALOG   = "Invalid document catalog"
const E_INVALID_DATE        = "Invalid date format in input"
const E_INVALID_DELIMITER = "Invalid delimiter character"
const E_INVALID_OBJECT      = "Invalid Object Found"
const E_INVALID_PAGE_NUMBER = "Page number is invalid"
const E_INVALID_PAGE_LABEL  = "Page label is invalid or page label definitions are not present in the document (Table 28 PDF Specification 1.7)"
const E_INVALID_PASSWORD    = "The password supplied to open the document is invalid"
const E_INVALID_CRYPT       = "The crypt handler is not supported"
const E_NOT_TAGGED_PDF      = "PDF file is not tagged"
const E_NOT_IMPLEMENTED     = "Not Implemented"


text:
abstract type Component end
using MLStyle: @λ

mutable struct Eqn
    sym
    rhs
    val
end
mutable struct Par
    sym
    val
end
mutable struct Coupling
    name
    ex
    flag
    targets
end

eqnhandle = @λ begin
    Expr(:tuple, Expr(:(=), sym, rhs), ic) -> Eqn(Symbol(String(sym)[2:end]),rhs,ic)
    Expr(:(=), sym, Expr(:tuple,rhs,ic)) -> Eqn(Symbol(String(sym)[2:end]),rhs,ic)
    _ -> nothing
end
couplinghandle = @λ begin
    Expr(:(=), name, Expr(:(->), flag, ex)) -> Coupling(name, ex.args[2], flag, Symbol[])
    _ -> nothing
end
parhandle = @λ begin
    Expr(:(=), sym, val) -> Par(sym, eval(val))
    _ -> nothing
end
macro Component(name, eqnex, parex)
    equations = filter(x -> x != nothing, eqnhandle.(eqnex.args))
    parameters = filter(x -> x != nothing, parhandle.(parex.args))
    couplings = filter(x -> x != nothing, couplinghandle.(eqnex.args))
    return quote
        struct $name <: Component
            name
            eqns::Vector{Eqn}
            pars::Vector{Par}
            couplings::Vector{Coupling}
        end
        function $(esc(name))(id; kwargs...)
            eqns = deepcopy($equations)
            pars = deepcopy($parameters)
            coups = deepcopy($couplings)
            # update default initial conditions
            for eqn in eqns
                if Symbol(eqn.sym,0) in keys(kwargs)
                    eqn.val = kwargs[Symbol(eqn.sym,0)]
                end
            end
            for par in pars
                if par.sym in keys(kwargs)
                    par.val = kwargs[par.sym]
                end
            end
            for coup in coups
                if coup.name in keys(kwargs)
                    val = kwargs[coup.name]
                    val isa Array ? coup.targets = val : coup.targets = [val]
                end
            end
            return $(esc(name))(id, eqns, pars, coups)
        end
    end
end


text:
include("line.jl")
include("rectangle.jl")
include("sphere.jl")
include("graph.jl")
include("partition.jl")

max_stars_repo_name: CliMA/GCMMesh

max_issues_repo_name: CliMA/GCMMesh

max_forks_repo_name: CliMA/GCMMesh