text:
abstract type AbstractUpdateMode end
struct UpdateMode <: AbstractUpdateMode end
struct RegenerateMode <: AbstractUpdateMode end

const retdiff = gensym("retdiff")
const discard = gensym("discard")

const calldiff_prefix = gensym("calldiff")
calldiff_var(node::GenerativeFunctionCallNode) = Symbol("$(calldiff_prefix)_$(node.addr)")

const choice_discard_prefix = gensym("choice_discard")
choice_discard_var(node::RandomChoiceNode) = Symbol("$(choice_discard_prefix)_$(node.addr)")

const call_discard_prefix = gensym("call_discard")
call_discard_var(node::GenerativeFunctionCallNode) = Symbol("$(call_discard_prefix)_$(node.addr)")

########################
# forward marking pass #
########################

struct ForwardPassState
    input_changed::Set{Union{RandomChoiceNode,GenerativeFunctionCallNode}}
    value_changed::Set{StaticIRNode}
    constrained_or_selected_choices::Set{RandomChoiceNode}
    constrained_or_selected_calls::Set{GenerativeFunctionCallNode}
    discard_calls::Set{GenerativeFunctionCallNode}
end

function ForwardPassState()
    input_changed = Set{Union{RandomChoiceNode,GenerativeFunctionCallNode}}()
    value_changed = Set{StaticIRNode}()
    constrained_or_selected_choices = Set{RandomChoiceNode}()
    constrained_or_selected_calls = Set{GenerativeFunctionCallNode}()
    discard_calls = Set{GenerativeFunctionCallNode}()
    ForwardPassState(input_changed, value_changed, constrained_or_selected_choices,
        constrained_or_selected_calls, discard_calls)
end

function forward_pass_argdiff!(state::ForwardPassState,
                               arg_nodes::Vector{ArgumentNode},
                               argdiffs_type::Type)
    for (node, diff_type) in zip(arg_nodes, argdiffs_type.parameters)
        if diff_type != NoChange
            push!(state.value_changed, node)
        end
    end
end

function process_forward!(::AddressSchema, ::ForwardPassState, ::TrainableParameterNode) end

function process_forward!(::AddressSchema, ::ForwardPassState, node::ArgumentNode) end

function process_forward!(::AddressSchema, state::ForwardPassState, node::JuliaNode)
    if any(input_node in state.value_changed for input_node in node.inputs)
        push!(state.value_changed, node)
    end
end

function process_forward!(schema::AddressSchema, state::ForwardPassState,
                          node::RandomChoiceNode)
    @assert isa(schema, StaticAddressSchema) || isa(schema, EmptyAddressSchema) || isa(schema, AllAddressSchema)
    if isa(schema, AllAddressSchema) || (isa(schema, StaticAddressSchema) && (node.addr in keys(schema)))
        push!(state.constrained_or_selected_choices, node)
        push!(state.value_changed, node)
    end
    if any(input_node in state.value_changed for input_node in node.inputs)
        push!(state.input_changed, node)
    end
end

function process_forward!(schema::AddressSchema, state::ForwardPassState,
                          node::GenerativeFunctionCallNode)
    @assert isa(schema, StaticAddressSchema) || isa(schema, EmptyAddressSchema) || isa(schema, AllAddressSchema)
    if isa(schema, AllAddressSchema) || (isa(schema, StaticAddressSchema) && (node.addr in keys(schema)))
        push!(state.constrained_or_selected_calls, node)
        push!(state.value_changed, node)
        push!(state.discard_calls, node)
    end
    if any(input_node in state.value_changed for input_node in node.inputs)
        push!(state.input_changed, node)
        push!(state.value_changed, node) # TODO can check whether the node is satically absorbing
        push!(state.discard_calls, node)
    end
end

#########################
# backward marking pass #
#########################

# this pass is used to determine which JuliaNodes need to be re-run (their
# return value is not currently cached in the trace)

struct BackwardPassState
    marked::Set{StaticIRNode}
end

function BackwardPassState()
    BackwardPassState(Set{StaticIRNode}())
end

function process_backward!(::ForwardPassState, ::BackwardPassState, ::TrainableParameterNode, options) end

function process_backward!(::ForwardPassState, ::BackwardPassState, ::ArgumentNode, options) end

function process_backward!(fwd::ForwardPassState, back::BackwardPassState,
                           node::JuliaNode, options)
    # if the node needs to be re-run which then we need the value of its inputs
    if ((options.cache_julia_nodes && node in fwd.value_changed) ||
        (!options.cache_julia_nodes && node in back.marked))
        for input_node in node.inputs
            push!(back.marked, input_node)
        end
    end
end

function process_backward!(fwd::ForwardPassState, back::BackwardPassState,
                           node::RandomChoiceNode, options)
    if node in fwd.input_changed || node in fwd.constrained_or_selected_choices
        for input_node in node.inputs
            push!(back.marked, input_node)
        end
    end
end

function process_backward!(fwd::ForwardPassState, back::BackwardPassState,
                           node::GenerativeFunctionCallNode, options)
    if node in fwd.input_changed || node in fwd.constrained_or_selected_calls
        for input_node in node.inputs
            push!(back.marked, input_node)
        end
    end
end

########################
# code generation pass #
########################

function arg_values_and_diffs_from_tracked_diffs(input_nodes)
    arg_values = map((node) -> Expr(:call, qn_strip_diff, node.name), input_nodes)
    arg_diffs = map((node) -> Expr(:call, qn_get_diff, node.name), input_nodes)
    (arg_values, arg_diffs)
end

function process_codegen!(stmts, ::ForwardPassState, back::BackwardPassState,
                          node::TrainableParameterNode, ::AbstractUpdateMode, options)
    if node in back.marked
        push!(stmts, :($(node.name) = $(QuoteNode(get_param))($(QuoteNode(get_gen_fn))(trace), $(QuoteNode(node.name)))))
    end
end

function process_codegen!(stmts, ::ForwardPassState, ::BackwardPassState,
                          node::ArgumentNode, ::AbstractUpdateMode, options)
    if options.track_diffs
        push!(stmts, :($(get_value_fieldname(node)) = $qn_strip_diff($(node.name))))
    else
        push!(stmts, :($(get_value_fieldname(node)) = $(node.name)))
    end
end

function process_codegen!(stmts, fwd::ForwardPassState, back::BackwardPassState,
                         node::JuliaNode, ::AbstractUpdateMode, options)
    run_it = ((options.cache_julia_nodes && node in fwd.value_changed) ||
              (!options.cache_julia_nodes && node in back.marked))
    if options.track_diffs

        # track diffs
        if run_it
            arg_values, arg_diffs = arg_values_and_diffs_from_tracked_diffs(node.inputs)
            args = map((v, d) -> Expr(:call, qn_Diffed, v, d), arg_values, arg_diffs)
            push!(stmts, :($(node.name) = $(QuoteNode(node.fn))($(args...))))
        elseif options.cache_julia_nodes
            push!(stmts, :($(node.name) = $qn_Diffed(trace.$(get_value_fieldname(node)), $qn_no_change)))
        end
        if options.cache_julia_nodes
            push!(stmts, :($(get_value_fieldname(node)) = $qn_strip_diff($(node.name))))
        end
    else

        # no track diffs
        if run_it
            arg_values = map((n) -> n.name, node.inputs)
            push!(stmts, :($(node.name) = $(QuoteNode(node.fn))($(arg_values...))))
        elseif options.cache_julia_nodes
            push!(stmts, :($(node.name) = trace.$(get_value_fieldname(node))))
        end
        if options.cache_julia_nodes
            push!(stmts, :($(get_value_fieldname(node)) = $(node.name)))
        end
    end
end

function process_codegen!(stmts, fwd::ForwardPassState, back::BackwardPassState,
                          node::RandomChoiceNode, ::UpdateMode,
                          options)
    if options.track_diffs

        # track diffs
        arg_values, _ = arg_values_and_diffs_from_tracked_diffs(node.inputs)
        new_logpdf = gensym("new_logpdf")
        addr = QuoteNode(node.addr)
        dist = QuoteNode(node.dist)
        if node in fwd.constrained_or_selected_choices || node in fwd.input_changed
            if node in fwd.constrained_or_selected_choices
                push!(stmts, :($(node.name) = $qn_Diffed($qn_static_get_value(constraints, Val($addr)), $qn_unknown_change)))
                push!(stmts, :($(choice_discard_var(node)) = trace.$(get_value_fieldname(node))))
            else
                push!(stmts, :($(node.name) = $qn_Diffed(trace.$(get_value_fieldname(node)), NoChange())))
            end
            push!(stmts, :($new_logpdf = $qn_logpdf($dist, $qn_strip_diff($(node.name)), $(arg_values...))))
            push!(stmts, :($weight += $new_logpdf - trace.$(get_score_fieldname(node))))
            push!(stmts, :($total_score_fieldname += $new_logpdf - trace.$(get_score_fieldname(node))))
            push!(stmts, :($(get_score_fieldname(node)) = $new_logpdf))
        else
            push!(stmts, :($(node.name) = $qn_Diffed(trace.$(get_value_fieldname(node)), $qn_no_change)))
            push!(stmts, :($(get_score_fieldname(node)) = trace.$(get_score_fieldname(node))))
        end
        push!(stmts, :($(get_value_fieldname(node)) = $(node.name)))

    else

        # no track diffs
        arg_values = map((n) -> n.name, node.inputs)
        new_logpdf = gensym("new_logpdf")
        addr = QuoteNode(node.addr)
        dist = QuoteNode(node.dist)
        if node in fwd.constrained_or_selected_choices || node in fwd.input_changed
            if node in fwd.constrained_or_selected_choices
                push!(stmts, :($(node.name) = $qn_static_get_value(constraints, Val($addr))))
                push!(stmts, :($(choice_discard_var(node)) = trace.$(get_value_fieldname(node))))
            else
                push!(stmts, :($(node.name) = trace.$(get_value_fieldname(node))))
            end
            push!(stmts, :($new_logpdf = $qn_logpdf($dist, $(node.name), $(arg_values...))))
            push!(stmts, :($weight += $new_logpdf - trace.$(get_score_fieldname(node))))
            push!(stmts, :($total_score_fieldname += $new_logpdf - trace.$(get_score_fieldname(node))))
            push!(stmts, :($(get_score_fieldname(node)) = $new_logpdf))
        else
            push!(stmts, :($(node.name) = trace.$(get_value_fieldname(node))))
            push!(stmts, :($(get_score_fieldname(node)) = trace.$(get_score_fieldname(node))))
        end
        push!(stmts, :($(get_value_fieldname(node)) = $(node.name)))
    end
end

function process_codegen!(stmts, fwd::ForwardPassState, back::BackwardPassState,
                          node::RandomChoiceNode, ::RegenerateMode,
                          options)
    if options.track_diffs

        # track diffs
        arg_values, _ = arg_values_and_diffs_from_tracked_diffs(node.inputs)
        new_logpdf = gensym("new_logpdf")
        addr = QuoteNode(node.addr)
        dist = QuoteNode(node.dist)
        if node in fwd.constrained_or_selected_choices || node in fwd.input_changed
            output_value = Expr(:call, qn_strip_diff, node.name)
            if node in fwd.constrained_or_selected_choices
                # the choice was selected, it does not contribute to the weight
                push!(stmts, :($(node.name) = $qn_Diffed($qn_random($dist, $(arg_values...)), UnknownChange())))
                push!(stmts, :($new_logpdf = $qn_logpdf($dist, $output_value, $(arg_values...))))
            else
                # the choice was not selected, and the input to the choice changed
                # it does contribute to the weight
                push!(stmts, :($(node.name) = $qn_Diffed(trace.$(get_value_fieldname(node)), NoChange())))
                push!(stmts, :($new_logpdf = $qn_logpdf($dist, $output_value, $(arg_values...))))
                push!(stmts, :($weight += $new_logpdf - trace.$(get_score_fieldname(node))))
            end
            push!(stmts, :($total_score_fieldname += $new_logpdf - trace.$(get_score_fieldname(node))))
            push!(stmts, :($(get_score_fieldname(node)) = $new_logpdf))
        else
            push!(stmts, :($(node.name) = $qn_Diffed(trace.$(get_value_fieldname(node)), NoChange())))
            push!(stmts, :($(get_score_fieldname(node)) = trace.$(get_score_fieldname(node))))
        end
        push!(stmts, :($(get_value_fieldname(node)) = $(node.name)))
    else

        # no track diffs
        arg_values = map((n) -> n.name, node.inputs)
        new_logpdf = gensym("new_logpdf")
        addr = QuoteNode(node.addr)
        dist = QuoteNode(node.dist)
        if node in fwd.constrained_or_selected_choices || node in fwd.input_changed
            if node in fwd.constrained_or_selected_choices
                # the choice was selected, it does not contribute to the weight
                push!(stmts, :($(node.name) = $qn_random($dist, $(arg_values...))))
                push!(stmts, :($new_logpdf = $qn_logpdf($dist, $(node.name), $(arg_values...))))
            else
                # the choice was not selected, and the input to the choice changed
                # it does contribute to the weight
                push!(stmts, :($(node.name) = trace.$(get_value_fieldname(node))))
                push!(stmts, :($new_logpdf = $qn_logpdf($dist, $(node.name), $(arg_values...))))
                push!(stmts, :($weight += $new_logpdf - trace.$(get_score_fieldname(node))))
            end
            push!(stmts, :($total_score_fieldname += $new_logpdf - trace.$(get_score_fieldname(node))))
            push!(stmts, :($(get_score_fieldname(node)) = $new_logpdf))
        else
            push!(stmts, :($(node.name) = trace.$(get_value_fieldname(node))))
            push!(stmts, :($(get_score_fieldname(node)) = trace.$(get_score_fieldname(node))))
        end
        push!(stmts, :($(get_value_fieldname(node)) = $(node.name)))
    end
end

function process_codegen!(stmts, fwd::ForwardPassState, back::BackwardPassState,
                          node::GenerativeFunctionCallNode, ::UpdateMode,
                          options)
    if options.track_diffs
        arg_values, arg_diffs = arg_values_and_diffs_from_tracked_diffs(node.inputs)
    else
        arg_values = map((n) -> n.name, node.inputs)
        arg_diffs = map((n) -> QuoteNode(n in fwd.value_changed ? UnknownChange() : NoChange()), node.inputs)
    end
    addr = QuoteNode(node.addr)
    gen_fn = QuoteNode(node.generative_function)
    subtrace = get_subtrace_fieldname(node)
    prev_subtrace = :(trace.$subtrace)
    call_weight = gensym("call_weight")
    call_constraints = gensym("call_constraints")
    if node in fwd.constrained_or_selected_calls || node in fwd.input_changed
        if node in fwd.constrained_or_selected_calls
            push!(stmts, :($call_constraints = $qn_static_get_submap(constraints, Val($addr))))
        else
            push!(stmts, :($call_constraints = $qn_empty_choice_map))
        end
        push!(stmts, :(($subtrace, $call_weight, $(calldiff_var(node)), $(call_discard_var(node))) =
            $qn_update($prev_subtrace, $(Expr(:tuple, arg_values...)), $(Expr(:tuple, arg_diffs...)), $call_constraints)))
        push!(stmts, :($weight += $call_weight))
        push!(stmts, :($total_score_fieldname += $qn_get_score($subtrace) - $qn_get_score($prev_subtrace)))
        push!(stmts, :($total_noise_fieldname += $qn_project($subtrace, $qn_empty_selection) - $qn_project($prev_subtrace, $qn_empty_selection)))
        push!(stmts, :(if !$qn_isempty($qn_get_choices($subtrace)) && $qn_isempty($qn_get_choices($prev_subtrace))
                            $num_nonempty_fieldname += 1 end))
        push!(stmts, :(if $qn_isempty($qn_get_choices($subtrace)) && !$qn_isempty($qn_get_choices($prev_subtrace))
                            $num_nonempty_fieldname -= 1 end))
        if options.track_diffs
            push!(stmts, :($(node.name) = $qn_Diffed($qn_get_retval($subtrace), $(calldiff_var(node)))))
        else
            push!(stmts, :($(node.name) = $qn_get_retval($subtrace)))
        end
    else
        push!(stmts, :($subtrace = $prev_subtrace))
        if options.track_diffs
            push!(stmts, :($(node.name) = $qn_Diffed($qn_get_retval($subtrace), $(QuoteNode(NoChange())))))
        else
            push!(stmts, :($(node.name) = $qn_get_retval($subtrace)))
        end
    end
end

function process_codegen!(stmts, fwd::ForwardPassState, back::BackwardPassState,
                          node::GenerativeFunctionCallNode, ::RegenerateMode,
                          options)
    if options.track_diffs
        arg_values, arg_diffs = arg_values_and_diffs_from_tracked_diffs(node.inputs)
    else
        arg_values = map((n) -> n.name, node.inputs)
        arg_diffs = map((n) -> QuoteNode(n in fwd.value_changed ? UnknownChange() : NoChange()), node.inputs)
    end
    addr = QuoteNode(node.addr)
    gen_fn = QuoteNode(node.generative_function)
    subtrace = get_subtrace_fieldname(node)
    prev_subtrace = :(trace.$subtrace)
    call_weight = gensym("call_weight")
    call_subselection = gensym("call_subselection")
    if node in fwd.constrained_or_selected_calls || node in fwd.input_changed
        if node in fwd.constrained_or_selected_calls
            push!(stmts, :($call_subselection = $qn_static_getindex(selection, Val($addr))))
        else
            push!(stmts, :($call_subselection = $qn_empty_selection))
        end
        push!(stmts, :(($subtrace, $call_weight, $(calldiff_var(node))) =
            $qn_regenerate($prev_subtrace, $(Expr(:tuple, arg_values...)), $(Expr(:tuple, arg_diffs...)), $call_subselection)))
        push!(stmts, :($weight += $call_weight))
        push!(stmts, :($total_score_fieldname += $qn_get_score($subtrace) - $qn_get_score($prev_subtrace)))
        push!(stmts, :($total_noise_fieldname += $qn_project($subtrace, $qn_empty_selection) - $qn_project($prev_subtrace, $qn_empty_selection)))
        push!(stmts, :(if !$qn_isempty($qn_get_choices($subtrace)) && !$qn_isempty($qn_get_choices($prev_subtrace))
                            $num_nonempty_fieldname += 1 end))
        push!(stmts, :(if $qn_isempty($qn_get_choices($subtrace)) && !$qn_isempty($qn_get_choices($prev_subtrace))
                            $num_nonempty_fieldname -= 1 end))
        if options.track_diffs
            push!(stmts, :($(node.name) = $qn_Diffed($qn_get_retval($subtrace), $(calldiff_var(node)))))
        else
            push!(stmts, :($(node.name) = $qn_get_retval($subtrace)))
        end
    else
        push!(stmts, :($subtrace = $prev_subtrace))
        if options.track_diffs
            push!(stmts, :($(node.name) = $qn_Diffed($qn_get_retval($subtrace), $qn_no_change)))
        else
            push!(stmts, :($(node.name) = $qn_get_retval($subtrace)))
        end
    end
end

function initialize_score_weight_num_nonempty!(stmts::Vector{Expr})
    push!(stmts, :($total_score_fieldname = trace.$total_score_fieldname))
    push!(stmts, :($total_noise_fieldname = trace.$total_noise_fieldname))
    push!(stmts, :($weight = 0.))
    push!(stmts, :($num_nonempty_fieldname = trace.$num_nonempty_fieldname))
end

function unpack_arguments!(stmts::Vector{Expr}, arg_nodes::Vector{ArgumentNode}, options)
    if options.track_diffs
        arg_names = Symbol[arg_node.name for arg_node in arg_nodes]
        push!(stmts, :($(Expr(:tuple, arg_names...)) = $(QuoteNode(map))($qn_Diffed, args, argdiffs)))
    else
        arg_names = Symbol[arg_node.name for arg_node in arg_nodes]
        push!(stmts, :($(Expr(:tuple, arg_names...)) = args))
    end
end

function generate_return_value!(stmts::Vector{Expr}, fwd::ForwardPassState, return_node::StaticIRNode, options)
    if options.track_diffs
        push!(stmts, :($return_value_fieldname = $qn_strip_diff($(return_node.name))))
        push!(stmts, :($retdiff = $qn_get_diff($(return_node.name))))
    else
        push!(stmts, :($return_value_fieldname = $(return_node.name)))
        push!(stmts, :($retdiff = $(QuoteNode(return_node in fwd.value_changed ? UnknownChange() : NoChange()))))
    end
end

function generate_new_trace!(stmts::Vector{Expr}, trace_type::Type, options)
    if options.track_diffs
        # note that the generative function is the last field
        constructor_args = map((name) -> Expr(:call, QuoteNode(strip_diff), name),
                            fieldnames(trace_type)[1:end-1])
        push!(stmts, :($trace = $(QuoteNode(trace_type))($(constructor_args...),
                        $(Expr(:(.), :trace, QuoteNode(static_ir_gen_fn_ref))))))
    else
        push!(stmts, :($static_ir_gen_fn_ref = $(Expr(:(.), :trace, QuoteNode(static_ir_gen_fn_ref)))))
        push!(stmts, :($trace = $(QuoteNode(trace_type))($(fieldnames(trace_type)...))))
    end
end

function generate_discard!(stmts::Vector{Expr},
                           constrained_choices::Set{RandomChoiceNode},
                           discard_calls::Set{GenerativeFunctionCallNode})
    discard_leaf_nodes = Dict{Symbol,Symbol}()
    for node in constrained_choices
        discard_leaf_nodes[node.addr] = choice_discard_var(node)
    end
    discard_internal_nodes = Dict{Symbol,Symbol}()
    for node in discard_calls
        discard_internal_nodes[node.addr] = call_discard_var(node)
    end
    if length(discard_leaf_nodes) > 0
        (leaf_keys, leaf_nodes) = collect(zip(discard_leaf_nodes...))
    else
        (leaf_keys, leaf_nodes) = ((), ())
    end
    if length(discard_internal_nodes) > 0
        (internal_keys, internal_nodes) = collect(zip(discard_internal_nodes...))
    else
        (internal_keys, internal_nodes) = ((), ())
    end
    leaf_keys = map((key::Symbol) -> QuoteNode(key), leaf_keys)
    internal_keys = map((key::Symbol) -> QuoteNode(key), internal_keys)
    expr = :($(QuoteNode(StaticChoiceMap))(
            $(QuoteNode(NamedTuple)){($(leaf_keys...),)}(($(leaf_nodes...),)),
            $(QuoteNode(NamedTuple)){($(internal_keys...),)}(($(internal_nodes...),))))
    push!(stmts, :($discard = $expr))
end

#######################
# generated functions #
#######################

function codegen_update(trace_type::Type{T}, args_type::Type, argdiffs_type::Type,
                        constraints_type::Type) where {T<:StaticIRTrace}
    gen_fn_type = get_gen_fn_type(trace_type)
    schema = get_address_schema(constraints_type)

    # convert the constraints to a static assignment if it is not already one
    if !(isa(schema, StaticAddressSchema) || isa(schema, EmptyAddressSchema))
        return quote $qn_update(trace, args, argdiffs, $(QuoteNode(StaticChoiceMap))(constraints)) end
    end

    ir = get_ir(gen_fn_type)
    options = get_options(gen_fn_type)

    # forward marking pass
    fwd_state = ForwardPassState()
    forward_pass_argdiff!(fwd_state, ir.arg_nodes, argdiffs_type)
    for node in ir.nodes
        process_forward!(schema, fwd_state, node)
    end

    # backward marking pass
    # TODO computes which values do we need to recompute - if we cache all Julia values, then we don't need to recompute anything.
    # TODO: we still need to extract only:
    # - trainable parameters nodes that we need
    # - julia node values that we need
    bwd_state = BackwardPassState()
    push!(bwd_state.marked, ir.return_node)
    for node in reverse(ir.nodes)
        process_backward!(fwd_state, bwd_state, node, options)
    end

    # forward code generation pass
    stmts = Expr[]
    initialize_score_weight_num_nonempty!(stmts)
    unpack_arguments!(stmts, ir.arg_nodes, options)
    for node in ir.nodes
        process_codegen!(stmts, fwd_state, bwd_state, node, UpdateMode(), options)
    end
    generate_return_value!(stmts, fwd_state, ir.return_node, options)
    generate_new_trace!(stmts, trace_type, options)
    generate_discard!(stmts, fwd_state.constrained_or_selected_choices, fwd_state.discard_calls)

    # return trace and weight and discard and retdiff
    push!(stmts, :(return ($trace, $weight, $retdiff, $discard)))

    Expr(:block, stmts...)
end

function codegen_regenerate(trace_type::Type{T}, args_type::Type, argdiffs_type::Type,
                            selection_type::Type) where {T<:StaticIRTrace}
    gen_fn_type = get_gen_fn_type(trace_type)
    schema = get_address_schema(selection_type)

    # convert a hierarchical selection to a static selection if it is not alreay one
    if !(isa(schema, StaticAddressSchema) || isa(schema, EmptyAddressSchema) || isa(schema, AllAddressSchema))
        return quote $qn_regenerate(trace, args, argdiffs, $(QuoteNode(StaticSelection))(selection)) end
    end

    ir = get_ir(gen_fn_type)
    options = get_options(gen_fn_type)

    # forward marking pass
    fwd_state = ForwardPassState()
    forward_pass_argdiff!(fwd_state, ir.arg_nodes, argdiffs_type)
    for node in ir.nodes
        process_forward!(schema, fwd_state, node)
    end

    # backward marking pass
    bwd_state = BackwardPassState()
    push!(bwd_state.marked, ir.return_node)
    for node in reverse(ir.nodes)
        process_backward!(fwd_state, bwd_state, node, options)
    end

    # forward code generation pass
    stmts = Expr[]
    initialize_score_weight_num_nonempty!(stmts)
    unpack_arguments!(stmts, ir.arg_nodes, options)
    for node in ir.nodes
        process_codegen!(stmts, fwd_state, bwd_state, node, RegenerateMode(), options)
    end
    generate_return_value!(stmts, fwd_state ,ir.return_node, options)
    generate_new_trace!(stmts, trace_type, options)

    # return trace and weight and retdiff
    push!(stmts, :(return ($trace, $weight, $retdiff)))

    Expr(:block, stmts...)
end

let T = gensym()
    push!(generated_functions, quote
    @generated function $(Expr(:(.), Gen, QuoteNode(:update)))(trace::$T, args::Tuple, argdiffs::Tuple,
                                   constraints::$(QuoteNode(ChoiceMap))) where {$T<:$(QuoteNode(StaticIRTrace))}
        $(QuoteNode(codegen_update))(trace, args, argdiffs, constraints)
    end
    end)

    push!(generated_functions, quote
    @generated function $(Expr(:(.), Gen, QuoteNode(:regenerate)))(trace::$T, args::Tuple, argdiffs::Tuple,
                                       selection::$(QuoteNode(Selection))) where {$T<:$(QuoteNode(StaticIRTrace))}
        $(QuoteNode(codegen_regenerate))(trace, args, argdiffs, selection)
    end
    end)
end

max_stars_repo_name: ali-ramadhan/Gen

max_issues_repo_name: ali-ramadhan/Gen

max_forks_repo_name: ali-ramadhan/Gen


text:
module test_goto

using Test

function gotojump(x)
    @goto jump
    x += 2
    @label jump
    x
end

@test 1 == gotojump(1)

function gototop(x)
    @label top
    x += 2
    x < 2 && @goto top
    x
end

@test 3 == gototop(1)

end

ext: jl


text:
"""
    is_spatial(x::StaticSymbol) -> StaticBool

Returns `static(true)` if `x` refers to a spatial dimension. By default, any dimension that
is not one of the following is considered spatial:

* time, Time
* channels, Channels
* observations, Observations, obs
"""
is_spatial(x::Symbol) = is_spatial(static(x))
is_spatial(x) = is_spatial(typeof(x))
is_spatial(::Type{StaticSymbol{sym}}) where {sym} = static(true)
const NotSpatial = Union{StaticSymbol{:time},StaticSymbol{:Time}, StaticSymbol{:Channels},
    StaticSymbol{:channels},StaticSymbol{:color}, StaticSymbol{:Color},
    StaticSymbol{:observations},StaticSymbol{:Observations},StaticSymbol{:obs}}
is_spatial(::Type{<:NotSpatial}) = static(false)

@inline function ArrayInterface.to_dims(::Type{T}, ::typeof(is_spatial)) where {T}
    _to_sdims(static(0), dimnames(T), Static.nstatic(Val(ndims(T))))
end

_to_sdims(::StaticInt, ::Tuple{}, ::Tuple{}) = ()
@inline function _to_sdims(
    ::StaticInt{N},
    dn::Tuple{StaticSymbol{name},Vararg{Any}},
    dims::Tuple{StaticInt{D},Vararg{Any}}
) where {N,name,D}

    if known(is_spatial(static(name)))
        if N === 2
            return (static(D),)
        else
            return (static(D), _to_sdims(static(N + 1), tail(dn), tail(dims))...)
        end
    else
        return _to_sdims(static(N + 1), tail(dn), tail(dims))
    end
end

"""
    spatialdims(x) -> Tuple{Vararg{StaticInt}}

Return a tuple of the spatial dimensions of `x`.
"""
spatialdims(@nospecialize(x)) = spatialdims(typeof(x))
@inline spatialdims(::Type{T}) where {T} = to_dims(T, is_spatial)

@inline spatial_axes(x) = _spatial_axes(x, spatialdims(x))
_spatial_axes(x, ::Tuple{}) = ()
@inline function _spatial_axes(x, d::Tuple{StaticInt{N},Vararg{Any}}) where {N}
    (ArrayInterface.axes(x, static(N)), _spatial_axes(x, tail(d))...)
end

"""
    spatial_order(x) -> Tuple

Returns the `dimnames` of `x` that correspond to spatial dimensions, as determined by
[`is_spatial`](@ref). At most three dimensions will be considered spatial.
"""
spatial_order(@nospecialize(x)) = spatial_order(typeof(x))
@inline spatial_order(::Type{T}) where {T} = Static.permute(dimnames(T), spatialdims(T))

"""
    sdims(x) -> Int

Return the number of spatial dimensions that `x` has.

See also: [`spatialdims`](@ref), [`is_spatial`](@ref)
"""
@inline sdims(x) = length(spatialdims(x))

"""
    spatial_indices(x)

Returns the keys along each spatial dimension.
"""
@inline spatial_indices(x) = map(keys, spatial_axes(x))
spatial_indices(x::MetaArray) = spatial_indices(parent(x))
# first get the spatial_indices of the parent array and then index into the ones that aren't
# dropped
function spatial_indices(x::SubArray)
    map(getindex,
        Static.permute(spatial_indices(parent(x)), to_spatialdims(typeof(x))),
        map(subdim -> getfield(x, :indices)[subdim], map(d -> to_parent_dims(x, d), spatialdims(x)))
    )
end
@inline function spatial_indices(x::Union{PermutedDimsArray,Adjoint,Transpose})
    Static.permute(spatial_indices(parent(x)), to_spatialdims(typeof(x)))
end

"""
    pixel_spacing(x)

Return a tuple representing the separation between adjacent pixels along each axis
of the image. Derived from the step size of each element of `spatial_indices`.
"""
@inline pixel_spacing(x) = map(_axis_pixel_spacing, spatial_indices(x))
@inline function pixel_spacing(x::SubArray)
    Static.permute(pixel_spacing(parent(x)), to_spatialdims(typeof(x)))
end
_axis_pixel_spacing(x::AbstractRange) = step(x)
_axis_pixel_spacing(x) = 1

"""
    origin(x) -> Tuple

Returns the spatial origin of `x`.
"""
origin(x) = spatial_first(x)
@inline origin(x::MetaArray) = getmeta(spatial_first, x, :origin)
@inline function origin(x::Union{PermutedDimsArray,SubArray,Adjoint,Transpose})
    Static.permute(origin(parent(x)), to_spatialdims(typeof(x)))
end

"""
    spatial_directions(img) -> (axis1, axis2, ...)

Return a tuple-of-tuples, each `axis[i]` representing the displacement
vector between adjacent pixels along spatial axis `i` of the image
array, relative to some external coordinate system ("physical
coordinates").

By default this is computed from `pixel_spacing`, but you can set this
manually using ImageMeta.
"""
@inline spatial_directions(x::X) where {X} = _spatial_directions(pixel_spacing(x))
function _spatial_directions(ps::NTuple{N,Any}) where {N}
    return ntuple(Val(N)) do i
        ntuple(Val(N)) do d
            if d === i
                getfield(ps, d)
            else
                zero(getfield(ps, d))
            end
        end
    end
end
@inline function spatial_directions(x::MetaArray)
    if has_metadata(x, :spatial_directions)
        return metadata(x, :spatial_directions)
    else
        return spatial_directions(parent(x))
    end
end
@inline function spatial_direction(x::Union{PermutedDimsArray,SubArray,Adjoint,Transpose})
    Static.permute(spatial_directions(parent(x)), to_spatialdims(typeof(x)))
end

"""
    spatial_first(x)

The first position along each spatial indices (i.e., where each spatial axis starts).
"""
spatial_first(x) = map(first, spatial_indices(x))

"""
    spatial_last(x)

The last position along each spatial indices (i.e., where each spatial axis stops).
"""
spatial_last(x) = map(last, spatial_indices(x))

function to_spatialdims(::Type{T}) where {T}
    _to_spatialdims(spatialdims(parent_type(T)), to_parent_dims(T))
end
@generated function _to_spatialdims(::SD, ::TD) where {SD,TD}
    out = Expr(:tuple)
    for p in TD.parameters
        i = findfirst(==(known(p)), known(SD))
        if i !== nothing
            push!(out.args, :(static($i)))
        end
    end
    Expr(:block, Expr(:meta, :inline), out)
end

"""
    width(x)

Returns the size of the dimension corresponding to width.
"""
width(x) = _width(spatialdims(x), x)
_width(::Union{Tuple{},Tuple{Any}}, x) = 1
_width(dims::Union{Tuple{Any,Any},Tuple{Any,Any,Any}}, x) = size(x, getfield(dims, 2))

"""
    height(x)

Returns the size of the dimension corresponding to height.
"""
height(x) = size(x, getfield(spatialdims(x), 1))

"""
    depth(x)

Returns the size of the dimension corresponding to depth.
"""
depth(x) = _depth(spatialdims(x), x)
_depth(::Union{Tuple{},Tuple{Any},Tuple{Any,Any}}, x) = 1
_depth(dims::Tuple{Any,Any,Any}, x) = size(x, getfield(dims, 3))


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct Rename <: Edit
        from::Symbol
        to::Symbol
    end

Change column name `from` -> `to`.

# Arguments
- `from::Symbol`: original column name
- `to::Symbol`: new column name
"""
mutable struct Rename <: Edit
    "original column name"
    from::Symbol
    "new column name"
    to::Symbol
end

function Rename(; from, to, )
    Rename(from, to, )
end

"""Get [`Rename`](@ref) `from`."""
get_from(value::Rename) = value.from
"""Get [`Rename`](@ref) `to`."""
get_to(value::Rename) = value.to

"""Set [`Rename`](@ref) `from`."""
set_from!(value::Rename, val) = value.from = val
"""Set [`Rename`](@ref) `to`."""
set_to!(value::Rename, val) = value.to = val

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
using Images, SIUnits
using Base.Test

const savedir = joinpath(tempdir(), "Images")
const writedir = joinpath(savedir, "write")

if !isdir(savedir)
    mkdir(savedir)
end
if !isdir(writedir)
    mkdir(writedir)
end

# Gray, raw
img = imread(joinpath(Pkg.dir(), "Images", "test", "io", "small.nrrd"))
@assert colorspace(img) == "Gray"
@assert ndims(img) == 3
@assert colordim(img) == 0
@assert eltype(img) == Float32
outname = joinpath(writedir, "small.nrrd")
imwrite(img, outname)
imgc = imread(outname)
@assert img.data == imgc.data

img = imread(joinpath(Pkg.dir(), "Images", "test", "io", "units.nhdr"))
ps = pixelspacing(img)
@test_approx_eq ps[1]/(0.1*Milli*Meter) 1
@test_approx_eq ps[2]/(0.2*Milli*Meter) 1
@test_approx_eq ps[3]/(1*Milli*Meter) 1

# Gray, compressed (gzip)
img = imread(joinpath(Pkg.dir(), "Images", "test", "io", "smallgz.nrrd"))
@assert colorspace(img) == "Gray"
@assert ndims(img) == 3
@assert colordim(img) == 0
@assert eltype(img) == Float32
outname = joinpath(writedir, "smallgz.nrrd")
imwrite(img, outname)
imgc = imread(outname)
@assert img.data == imgc.data


text:
# matrix of zeros

export joZeros
"""
    julia> op = joZeros(m[,n];[DDT=...,][RDT=...,][name=...])

Operator equivalent to matrix of ones

# Signature

    joZeros(m::Integer,n::Integer=m;
        DDT::DataType=joFloat,RDT::DataType=DDT,name::String="joZeros")

# Arguments

- `m`: size
- optional
    - `n`: 2nd dimension if not square
- keywords
    - `DDT`: domain data type
    - `RDT`: range data type
    - `name`: custom name

# Examples

    joZeros(m)
    joZeros(m,n)

examples with DDT/RDT

    joZeros(m; DDT=Float32)
    joZeros(m; DDT=Float32,RDT=Float64)

"""
joZeros(m::Integer,n::Integer=m;DDT::DataType=joFloat,RDT::DataType=DDT,name::String="joZeros") =
    joMatrix{DDT,RDT}(name,m,n,
        v1->(size(v1,2)>1 ? zeros(RDT,m,size(v1,2)) : zeros(RDT,m)),
        v2->(size(v2,2)>1 ? zeros(DDT,n,size(v2,2)) : zeros(DDT,n)),
        v3->(size(v3,2)>1 ? zeros(DDT,n,size(v3,2)) : zeros(DDT,n)),
        v4->(size(v4,2)>1 ? zeros(RDT,m,size(v4,2)) : zeros(RDT,m)),
        @joNF, @joNF, @joNF, @joNF
        )


text:
function build_nabla2_x(LF)
  D2jl = LF.LFx.D2jl
  Nx = LF.Nx
  Ny = LF.Ny
  Nz = LF.Nz
  Npoints = Nx*Ny*Nz

  # initialize with
  rowval = zeros(Int64,Nx*Npoints)
  colptr = zeros(Int64,Npoints+1)
  colptr[1] = 1
  nzval = zeros(Float64,Nx*Npoints)

  idx = 1
  for colLoc = 1:Nx
    colGbl_start = (colLoc-1)*(Ny*Nz)+1
    colGbl_stop = colLoc*Ny*Nz
    for colGbl = colGbl_start:colGbl_stop
      for rowLoc = 1:Nx
        rowGbl = (rowLoc-1)*Ny*Nz + colGbl - (colLoc-1)*Ny*Nz
        rowval[idx] = rowGbl
        nzval[idx] = D2jl[rowLoc,colLoc]
        idx = idx + 1
      end
      colptr[colGbl+1] = colptr[colGbl] + Nx
    end
  end # colLoc

  return SparseMatrixCSC( Npoints, Npoints, colptr, rowval, nzval )
end


text:
@testset "VariationalDenseBE" begin
    ensemble_size = 4
    samples_per_model = 4
    input_dim = 5
    output_dim = 5
    rank = 1
    inputs = rand(Float32, input_dim, samples_per_model) |> gpu
    layer = VariationalDenseBE(
        input_dim,
        output_dim,
        rank,
        ensemble_size;
        device = gpu
    ) |> gpu
    batch_inputs = repeat(inputs, 1, ensemble_size) |> gpu
    batch_outputs = layer(batch_inputs) |> gpu

    # Test gradients
    grads = gradient(params(layer)) do
        ŷ = layer(batch_inputs)
        return Flux.logitcrossentropy(ŷ, batch_outputs)
    end
    for param in params(layer)
        @test size(param) == size(grads[param])
        @test grads[param] isa CuArray
    end
end

@testset "VariationalConvBE" begin
    ensemble_size = 4
    samples_per_model = 4
    input_dim = 5
    output_dim = 10
    rank = 1
    inputs = rand(Float32, 10, 10, input_dim, samples_per_model)
    layer = VariationalConvBE(
        (5, 5),
        5 => 10,
        rank,
        ensemble_size,
        relu;
        device = gpu
    ) |> gpu
    batch_inputs = repeat(inputs, 1, 1, 1, ensemble_size) |> gpu
    batch_outputs = layer(batch_inputs)

    # Test gradients
    grads = gradient(params(layer)) do
        ŷ = layer(batch_inputs)
        return Flux.logitcrossentropy(ŷ, batch_outputs)
    end
    for param in params(layer)
        @test size(param) == size(grads[param])
        @test grads[param] isa CuArray
    end
end


text:
function isdf_pp1( filename )  #=

perform the preprocessing stage 1 which includes

- error checking
- expansion of the use sections

and returns the result as an EzXML document.

Author: Taku Noda
Started on Nov. 29, 2020
Update on June 29, 2021

Copyright 2021 Central Research Institute of Electric Power Industry

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

=#

if length( filename ) <= 4
  omegasim_error( "The file name is too short." )
end
ext = filename[ ( end - 3 ):end ]
if ext != ".xml"
  omegasim_error( "The file extension must be \".xml\"." )
end
filename_pp1 = ( filename[ 1:( end - 4 ) ] )*".pp1.xml"

sys = read_isdf( filename )
sys_pp1 = expand_uses( sys )

io = open( filename_pp1, "w" )
prettyprint( io, sys_pp1 )
close( io )

return sys_pp1

end

max_stars_repo_name: takunoda/omegasim

max_issues_repo_name: takunoda/omegasim

max_forks_repo_name: takunoda/omegasim


text:
struct ParenthesizeNegative <: FormatRequirement
    symloc::Symbol  # :inside, :outside, or :unspecified
end
ParenthesizeNegative() = ParenthesizeNegative(:unspecified)

struct DigitSeparator <: FormatRequirement
    sep::String        # e.g. ",", ".", "'", " "; "\0" for unspecified
    rule::Tuple{Int, Int}  # first, rest, e.g. (3, 3); (0, 0) for unspecified
end
DigitSeparator(c::AbstractString) = DigitSeparator(c, (3, 3))

struct DecimalSeparator <: FormatRequirement
    sep::String  # e.g. ",", "."; "\0" for unspecified
end

struct CurrencySymbol <: FormatRequirement
    symtype::Symbol   # :short, :long, :iso4217, or :unspecified
    location::Symbol  # :before, :after, :none, :dependent, or :unspecified
    spacing::Symbol   # :space, :none, :dependent, or :unspecified
    glued::Symbol     # :require, :disallow, :unspecified
    compose::Vector{Function}  # apply each function in turn to result
end
function CurrencySymbol(;
    symtype=:unspecified, location=:unspecified,
    spacing=:unspecified, glued=:unspecified,
    compose=Function[])

    CurrencySymbol(symtype, location, spacing, glued, compose)
end

struct RenderAs <: FormatRequirement
    sym::Symbol                            # name of symbol to require
    options::Dict{String, Int}  # Allowable symbols => priority
end

reconcile(p::ParenthesizeNegative, p′::ParenthesizeNegative) =
    [ParenthesizeNegative(takenonzero(p.symloc, p′.symloc, :unspecified))]
function reconcile(d::DigitSeparator, d′::DigitSeparator)
    [DigitSeparator(
        takenonzero(d.sep, d′.sep, "\0"),
        takenonzero(d.rule, d′.rule, (0, 0)))]
end
reconcile(d::DecimalSeparator, d′::DecimalSeparator) =
    [DecimalSeparator(takenonzero(d.sep, d′.sep, "\0"))]
reconcile(s::CurrencySymbol, s′::CurrencySymbol) =
    [CurrencySymbol(
        takenonzero(s.symtype, s′.symtype, :unspecified),
        takenonzero(s.location, s′.location, :unspecified),
        takenonzero(s.spacing, s′.spacing, :unspecified),
        takenonzero(s.glued, s′.glued, :unspecified),
        s.compose ∪ s′.compose)]
function reconcile(s::RenderAs, s′::RenderAs)
    @assert s.sym ≡ s′.sym
    allowed = allowable(s.options, s′.options)
    if isempty(allowed)
        throw(IncompatibleFormatException("Cannot agree on symbol: $(s.sym)."))
    end
    [RenderAs(s.sym, allowed)]
end

# require reconciliation if same type
conflict(x::RenderAs, y::RenderAs) = x.sym ≡ y.sym

function getsymboltable(spec::FormatSpecification)
    default = Dict(
        :minus_sign => "−",
        :zero_dash => "—",
        :thin_space => " ")
    table = Dict{Symbol, Dict{String, Int}}()
    for req in spec.reqs
        if isa(req, RenderAs)
            table[req.sym] = req.options
        end
    end
    for (k, v) in table
        default[k] = best(v)
    end
    default
end


text:
module Mesh

# using Triangle

export define_line_element, create_line_element_property_vector, extrude_nodes, extrude_elements, open_cross_section_tessellation, surface,
   calculate_weights

"""
    define_line_element(member_definitions)

Accepts `member_definitions` tuple and returns element discretization `dz`, line coordinates 'z', and nodal properties at each node in the line `node_props`.
"""

function define_line_element(member_definitions)

   #mesh along the length
   for i in eachindex(member_definitions)

      L = member_definitions[i][1]
      dL = member_definitions[i][2]
      num_segments = round(Int64, L/dL)

      if i == 1
         dz = ones(num_segments)*dL  #member discretization
      else
         dz = [dz; ones(num_segments)*dL]
      end

   end
   dz = dz
   z = [0; cumsum(dz)]

   #define what properties to apply at each node
   node_props = assign_line_element_nodal_properties(member_definitions)

   return dz, z, node_props

end

"""
   assign_line_element_nodal_properties(member_definitions)

Accepts `member_definitions` tuple and returns a properties assignment for each node line segment `dm`.

There are some tricks to defining the nodal properties. If the number of `member_definitions` is even, the first and last segments are added as `num_nodes - 1`.  If the number of `member_definitions` is odd, middle segment is added as 'num_nodes - 1'.

"""

function assign_line_element_nodal_properties(member_definitions)

   #define member properties to use at each node

   if length(member_definitions) == 1

      L = member_definitions[1][1]
      dL = member_definitions[1][2]
      num_segments = round(Int64, L/dL)
      num_nodes = num_segments + 1   #number of nodes in a segment

      dm = ones(Int8, num_nodes)

   elseif iseven(length(member_definitions))

      for i in eachindex(member_definitions)

         L = member_definitions[i][1]
         dL = member_definitions[i][2]
         num_segments = floor(Int64, L/dL)
         num_nodes = num_segments+1

         if i == 1
            dm = ones(Int8, num_nodes-1)*i
         elseif i == length(member_definitions)
            dm = [dm; ones(Int8, num_nodes)*i]
         else
            dm = [dm; ones(Int8, num_nodes-1)*i]
         end

      end

   elseif isodd(length(member_definitions))

      middle_segment = (length(member_definitions) - 1)/2 + 1

      for i in eachindex(member_definitions)

         L = member_definitions[i][1]
         dL = member_definitions[i][2]
         num_segments = floor(Int64, L/dL)
         num_nodes = num_segments + 1

         if i == 1
            dm = ones(Int8, num_nodes-1)*i
         elseif i < middle_segment
            dm = [dm; ones(Int8, num_nodes-1)*i]
         elseif i == middle_segment
            dm = [dm; ones(Int8, num_nodes)*i]
         elseif i > middle_segment
            dm = [dm; ones(Int8, num_nodes-1)*i]
         end

      end

   end

   dm = dm

   return dm

end

"""
   create_line_element_property_vector(member_definitions, dm, dz, property, property_order, property_type)

Accepts the `member_definitions` tuple, property array 'dm' assigning a property to each node, discretization along the line segment `dz`, the `property` magnitude, `property_order` which defines how the properties are ordered in `member_definitions` , and `property_type` which is the flavor of the property (e.g., `Ix` or `Iy` or `J` or ...) .  The function returns `property_array` with the magnitude of a property assigned node by node.

"""

function create_line_element_property_array(member_definitions, dm, dz, property, property_order, property_type)

   z = [0.0; cumsum(dz)]

   num_nodes = length(dm)

   property_array = zeros(num_nodes)

   for i=1:num_nodes
      property_index = member_definitions[dm[i]][property_order]  #maps properties to each segment along a line element
      property_array[i] = property[property_index][property_type]
   end

   return property_array

end

"""
   extrude_nodes(x, y, z)

Accepts the `x` and `y` cross-section coordinates and `z` coordinates along the extrusion. Returns the 3D `coordinates` in the form `[x y z]` from the extrusion.
"""

function extrude_nodes(x, y, z)

   num_nodes_section = length(x)
   coordinates = zeros(num_nodes_section, 3)   #initialize coordinates array

   for i in 1:length(z)
      z_section = z[i] .* ones(length(x))
      if i==1
         coordinates .= [x y z_section]   #start with first cross-section
      else
         coordinates = vcat(coordinates, [x y z_section])  #add following cross-sections along z
      end
   end

   return coordinates

end

"""
   extrude_elements(num_cross_section_nodes, num_cross_section_elements, num_extruded_layers)

Returns a triangular surface mesh `connectivity` matrix as `[node_i node_j node_k]`.   If `num_cross_section_nodes = num_cross_section_elements`, then the cross-section is closed or multi-branched.

"""

function extrude_elements(num_cross_section_nodes, num_cross_section_elements, num_extruded_layers)

   num_elements = num_cross_section_elements * (num_extruded_layers - 1)

   uppoly = zeros(Int, num_elements, 3)
   downpoly = zeros(Int, num_elements, 3)

   for j = 1:num_extruded_layers - 1

      node_offset =  num_cross_section_nodes * (j - 1)
      elem_offset =  (j - 1) * num_cross_section_elements

      for i = 1:num_cross_section_elements

         element_count = i + elem_offset

         #triangular element
                     #c
                  #  #
               #     #
            #        #
         #           #
        #a###########b

         a_up = i
         b_up = i + 1
         c_up = i + 1 + num_cross_section_nodes

         #address closed or multi-branch cross-sections
         if i == num_cross_section_nodes
            b_up = 1 + 1
            c_up = 1 + 1 + num_cross_section_nodes
         end

         uppoly[element_count, :] = [a_up b_up c_up] .+ node_offset

         #triangular element
         #b#########c
         #        #
         #      #
         #    #
         #  #
         #a

         a_down = i
         b_down = i + num_cross_section_nodes
         c_down = i + num_cross_section_nodes + 1

         #address closed or multi-branch cross-sections
         if i == num_cross_section_nodes
            c_down = 1 + num_cross_section_nodes + 1
         end

         downpoly[element_count, :] = [a_down b_down c_down] .+ node_offset

      end

   end

   connectivity=vcat(uppoly,downpoly)

   return connectivity

end

# """
#    open_cross_section_tessellation(x, y)

# Accepts `x` and `y` coordinates for an open cross-section and determines a through-thickness triangular mesh `shape_mesh`.  The Array `cross_section_edges = [node(i) node(j)]` is also output, providing the nodal connectivity that describes the cross-section boundary.

# """

# function open_cross_section_tessellation(x, y)

#     #define cross-section edges
#     cross_section_edges = zeros(Int, length(x), 2)
#     for i = 1:length(x) - 1
#         cross_section_edges[i, 1] = i
#         cross_section_edges[i, 2] = i + 1
#     end
#     cross_section_edges[end, :] = [length(x) 1]  #for open cross-section only

#     points = [x y]
#     points_map = [i for i=1:length(x)]
#     shape_mesh = Triangle.constrained_triangulation(points, points_map, cross_section_edges)

#     return shape_mesh, cross_section_edges

# end

"""
   surface(xcoords, ycoords, zcoords)

Accepts 3D surface coordinates of a structural member and defines a triangular surface mesh with `coordinates` and `connectivity`.

"""

function surface(xcoords, ycoords, zcoords)

   #generate triangular mesh at the ends
   section_mesh, cross_section_edges = open_cross_section_tessellation(xcoords, ycoords)

   #get member 3D coordinates via extrusion
   coordinates = extrude_nodes(xcoords, ycoords, zcoords)

   #define triangular mesh along extrusion
   num_cross_section_nodes = length(xcoords)
   num_extruded_layers = length(zcoords)
   num_cross_section_elements = size(cross_section_edges)[1]
   extrusion_connectivity = extrude_elements(num_cross_section_nodes, num_cross_section_elements, num_extruded_layers)

   # mesh beginning and end faces of extrusion
   #beginning face
   begin_connectivity = zeros(Int, length(section_mesh), 3)
   for i=1:length(section_mesh)
       begin_connectivity[i, :] = section_mesh[i]
   end

   #end face
   end_connectivity = begin_connectivity .+ num_cross_section_nodes * (num_extruded_layers - 1)

   #combine to define the full surface with end caps
   connectivity = [extrusion_connectivity; begin_connectivity; end_connectivity]

   return coordinates, connectivity

end

#############################################################
# Fornberg algorithm

# This implements the Fornberg (1988) algorithm (https://doi.org/10.1090/S0025-5718-1988-0935077-0)
# to obtain Finite Difference weights over arbitrary points to arbitrary order.

function calculate_weights(order::Int, x0::T, x::AbstractVector) where T<:Real
   #=
       order: The derivative order for which we need the coefficients
       x0   : The point in the array 'x' for which we need the coefficients
       x    : A dummy array with relative coordinates, e.g., central differences
              need coordinates centred at 0 while those at boundaries need
              coordinates starting from 0 to the end point
       The approximation order of the stencil is automatically determined from
       the number of requested stencil points.
   =#
   N = length(x)
   @assert order < N "Not enough points for the requested order."
   M = order
   c1 = one(T)
   c4 = x[1] - x0
   C = zeros(T, N, M+1)
   C[1,1] = 1
   @inbounds for i in 1 : N-1
       i1 = i + 1
       mn = min(i, M)
       c2 = one(T)
       c5 = c4
       c4 = x[i1] - x0
       for j in 0 : i-1
           j1 = j + 1
           c3 = x[i1] - x[j1]
           c2 *= c3
           if j == i-1
               for s in mn : -1 : 1
                   s1 = s + 1
                   C[i1,s1] = c1*(s*C[i,s] - c5*C[i,s1]) / c2
               end
               C[i1,1] = -c1*c5*C[i,1] / c2
          end
           for s in mn : -1 : 1
               s1 = s + 1
               C[j1,s1] = (c4*C[j1,s1] - s*C[j1,s]) / c3
           end
           C[j1,1] = c4 * C[j1,1] / c3
       end
       c1 = c2
   end
   #=
       This is to fix the problem of numerical instability which occurs when the sum of the stencil_coefficients is not
       exactly 0.
       https://scicomp.stackexchange.com/questions/11249/numerical-derivative-and-finite-difference-coefficients-any-update-of-the-fornb
       Stack Overflow answer on this issue.
       http://epubs.siam.org/doi/pdf/10.1137/S0036144596322507 - Modified Fornberg Algorithm
   =#
   _C = C[:,end]
   _C[div(N,2)+1] -= sum(_C)
   return _C
end

end #module


text:
# indicator of second-order cones

"""
  IndSOC()

Returns the indicator of the second-order cone (ice-cream cone) of R^n.
"""

immutable IndSOC <: IndicatorConvexCone end

function (f::IndSOC){T <: Real}(x::AbstractArray{T,1})
  # the tolerance in the following line should be customizable
  if norm(x[2:end]) - x[1] <= 1e-14
    return 0.0
  end
  return +Inf
end

function prox!{T <: Real}(f::IndSOC, x::AbstractArray{T,1}, y::AbstractArray{T,1}, gamma::Real=1.0)
  nx = norm(x[2:end])
  t = x[1]
  if t <= -nx
    y[:] = 0.0
  elseif t >= nx
    y[:] = x
  else
    r = 0.5 * (1 + t / nx)
    y[1] = r * nx
    y[2:end] = r * x[2:end]
  end
  return 0.0
end

fun_name(f::IndSOC) = "indicator of the second-order cone"
fun_dom(f::IndSOC) = "AbstractArray{Real,1}"
fun_expr(f::IndSOC) = "x ↦ 0 if x[1] >= ||x[2:end]||, +∞ otherwise"
fun_params(f::IndSOC) = "none"

function prox_naive{T <: Real}(f::IndSOC, x::AbstractArray{T,1}, gamma::Real=1.0)
  nx = norm(x[2:end])
  t = x[1]
  if t <= -nx
    y = zeros(x)
  elseif t >= nx
    y = x
  else
    y = zeros(x)
    r = 0.5 * (1 + t / nx)
    y[1] = r * nx
    y[2:end] = r * x[2:end]
  end
  return y, 0.0
end

# ########################
# ROTATED SOC
# ########################

"""
  IndRotatedSOC()

Returns the indicator of the *rotated* second-order cone of R^n, that is {(p,q,x) : norm(x)^2 ⩽ 2pq, p ⩾ 0, q ⩾ 0}
"""

immutable IndRotatedSOC <: IndicatorConvexCone end

function (f::IndRotatedSOC){T <: Real}(x::AbstractArray{T,1})
  if x[1] >= -1e-14 && x[2] >= -1e-14 && norm(x[3:end])^2 - 2*x[1]*x[2] <= 1e-14
    return 0.0
  end
  return +Inf
end

function prox!{T <: Real}(f::IndRotatedSOC, x::AbstractArray{T,1}, y::AbstractArray{T,1}, gamma::Real=1.0)
  # sin(pi/4) = cos(pi/4) = 0.7071067811865475
  # rotate x ccw by pi/4
  x1 = 0.7071067811865475*x[1] + 0.7071067811865475*x[2]
  x2 = 0.7071067811865475*x[1] - 0.7071067811865475*x[2]
  # project rotated x onto SOC
  nx = sqrt(x2^2+norm(x[3:end])^2)
  t = x1
  if t <= -nx
    y[:] = 0.0
  elseif t >= nx
    y[1] = x1
    y[2] = x2
    y[3:end] = x[3:end]
  else
    r = 0.5 * (1 + t / nx)
    y[1] = r * nx
    y[2] = r * x2
    y[3:end] = r * x[3:end]
  end
  # rotate back y cw by pi/4
  y1 = 0.7071067811865475*y[1] + 0.7071067811865475*y[2]
  y2 = 0.7071067811865475*y[1] - 0.7071067811865475*y[2]
  y[1] = y1
  y[2] = y2
  return 0.0
end

fun_name(f::IndRotatedSOC) = "indicator of the rotated second-order cone"
fun_dom(f::IndRotatedSOC) = "AbstractArray{Real,1}"
fun_expr(f::IndRotatedSOC) = "x ↦ 0 if x[1] ⩾ 0, x[2] ⩾ 0, norm(x[3:end])² ⩽ 2*x[1]*x[2], +∞ otherwise"
fun_params(f::IndRotatedSOC) = "none"

function prox_naive{T <: Real}(f::IndRotatedSOC, x::AbstractArray{T,1}, gamma::Real=1.0)
  g = IndSOC()
  z = copy(x)
  z[1] = 0.7071067811865475*x[1] + 0.7071067811865475*x[2]
  z[2] = 0.7071067811865475*x[1] - 0.7071067811865475*x[2]
  y, = prox_naive(g, z, gamma)
  y1 = 0.7071067811865475*y[1] + 0.7071067811865475*y[2]
  y2 = 0.7071067811865475*y[1] - 0.7071067811865475*y[2]
  y[1] = y1
  y[2] = y2
  return y, 0.0
end


text:
using MOCNeutronTransport
@testset "Quadrilateral8_3D" begin
    for T in [Float32, Float64]
        @testset "Constructors" begin
            p₁ = Point_3D(T, 0)
            p₂ = Point_3D(T, 1)
            p₃ = Point_3D(T, 1, 1)
            p₄ = Point_3D(T, 0, 1)
            p₅ = Point_3D(T, 1//2,    0)
            p₆ = Point_3D(T,    1, 1//2)
            p₇ = Point_3D(T, 1//2,    1)
            p₈ = Point_3D(T,    0, 1//2)
            quad8 = Quadrilateral8_3D((p₁, p₂, p₃, p₄, p₅, p₆, p₇, p₈))
            @test quad8.points == (p₁, p₂, p₃, p₄, p₅, p₆, p₇, p₈)

            # single constructor
            quad8 = Quadrilateral8_3D((p₁, p₂, p₃, p₄, p₅, p₆, p₇, p₈))
            @test quad8.points == (p₁, p₂, p₃, p₄, p₅, p₆, p₇, p₈)
        end

        @testset "Methods" begin
            p₁ = Point_3D(T, 0)
            p₂ = Point_3D(T, 1)
            p₃ = Point_3D(T, 1, 1)
            p₄ = Point_3D(T, 0, 1)
            p₅ = Point_3D(T, 1//2,    0)
            p₆ = Point_3D(T,    1, 1//2)
            p₇ = Point_3D(T, 1//2,    1)
            p₈ = Point_3D(T,    0, 1//2)
            quad8 = Quadrilateral8_3D((p₁, p₂, p₃, p₄, p₅, p₆, p₇, p₈))

            # interpolation
            @test quad8(0, 0) ≈ p₁
            @test quad8(1, 0) ≈ p₂
            @test quad8(1, 1) ≈ p₃
            @test quad8(0, 1) ≈ p₄
            @test quad8(1//2,    0) ≈ p₅
            @test quad8(   1, 1//2) ≈ p₆
            @test quad8(1//2,    1) ≈ p₇
            @test quad8(   0, 1//2) ≈ p₈
            @test quad8(1//2, 1//2) ≈ Point_3D(T, 1//2, 1//2)

            # area
            p₁ = Point_3D(T, 0)
            p₂ = Point_3D(T, 2)
            p₃ = Point_3D(T, 2, 3)
            p₄ = Point_3D(T, 0, 3)
            p₅ = Point_3D(T, 3//2, 1//2)
            p₆ = Point_3D(T, 5//2, 3//2)
            p₇ = Point_3D(T, 3//2, 5//2)
            p₈ = Point_3D(T, 0,    3//2)
            quad8 = Quadrilateral8_3D((p₁, p₂, p₃, p₄, p₅, p₆, p₇, p₈))

            p₈ = Point_3D(T, 0, 3//2, 2)
            quad8 = Quadrilateral8_3D((p₁, p₂, p₃, p₄, p₅, p₆, p₇, p₈))
            @test isapprox(area(quad8), 10.09431, atol=1.0e-4)

            # intersect
            p₁ = Point_3D(T, 0)
            p₂ = Point_3D(T, 2)
            p₃ = Point_3D(T, 2, 3)
            p₄ = Point_3D(T, 0, 3)
            p₅ = Point_3D(T, 3//2, 1//2)
            p₆ = Point_3D(T, 5//2, 3//2)
            p₇ = Point_3D(T, 3//2, 5//2)
            p₈ = Point_3D(T, 0,    3//2)
            quad8 = Quadrilateral8_3D((p₁, p₂, p₃, p₄, p₅, p₆, p₇, p₈))
            l = LineSegment_3D(Point_3D(T, 1, 3//2, -2),
                               Point_3D(T, 1, 3//2,  2))
            intersection = l ∩ quad8

            # 1 intersection
            @test intersection[1]
            @test intersection[2] == 1
            @test intersection[3] ≈ Point_3D(T, 1, 3//2, 0)

            # 0 intersection
            l = LineSegment_3D(Point_3D(T, 1, 0, -2),
                               Point_3D(T, 1, 0,  2))
            intersection = l ∩ quad8
            @test !intersection[1]

            # 2 intersections
            p₈ = Point_3D(T, 0, 3//2, 2)
            quad8 = Quadrilateral8_3D((p₁, p₂, p₃, p₄, p₅, p₆, p₇, p₈))
            l = LineSegment_3D(Point_3D(T, 1, 0, 1//2),
                               Point_3D(T, 1, 3, 1//2))
            intersection = l ∩ quad8
            @test intersection[1]
            @test intersection[2] == 2
            @test norm(intersection[3] - Point_3D(T, 1, 0.6099032181207209, 1//2)) < 1.0e-3
            @test norm(intersection[4] - Point_3D(T, 1, 2.390096778327549, 1//2)) < 1.0e-3
        end
    end
end

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
"""
Plots the index of refraction ath certain wavelength (usually λ0) of the multilayer structure computed.
using PyPlot, PlotUtils, LaTeXStrings
using Printf: @sprintf
See examples folder for more details.
"""

using PyPlot, PlotUtils, LaTeXStrings
using Printf: @sprintf

function nplot(λ::AbstractArray{T1,N1}, θ::AbstractArray{T2,N2}, λ0::T3, d::AbstractArray{T3,N4}, τ::AbstractArray{T4,N5}, emf::AbstractArray{T3,N6}, nλ0::AbstractArray{T5,N7}, n) where {T1<:Number, N1, T2<:Number, N2, T3<:Float64, N4, T4<:Float64, N5, N6, T5<:ComplexF64, N7}

    # check input
    nn = lastindex(n)
    nd = lastindex(d)
    nd == lastindex(nλ0)-2 || throw(DimensionMismatch("lastindex(nλ0)-2 not equal to lastindex(d)"))

    # define the offset for the incident and substrate medium
    doffset = 200. # nm

    # Maps each certain to a certain color
    aux2 = unique(n)
    naux2 = lastindex(aux2)
    cm = cgrad(:Spectral)
    β = [cm[i] for i in LinRange(0,1,naux2)] # range(0,stop=1,length=naux2) should work also
    aux3 = [getfield(β[j],i) for i=1:4, j=1:lastindex(β)]

    βassigned = Array{Float64,2}(undef, (4,nn))
    for idx2 = 1 : naux2, idx1 = 1 : nn
        if aux2[idx2] == n[idx1]
            βassigned[:,idx1] = aux3[:,idx2]
        end
    end

    nλ0 = real.(nλ0)
    figure()
    ax_profile = gca()
    # draw incident medium
    ax_profile[:fill_between]([-doffset, 0], 0, nλ0[1], facecolor=βassigned[:,1], interpolate="true", alpha=0.8, edgecolor = "black")

    new_d = cumsum([0; d], dims=1)
    for idx = 1 : nd
        # draw successive layers
        ax_profile[:fill_between]([new_d[idx], new_d[idx+1]], 0, nλ0[idx+1], facecolor=βassigned[:,idx+1], interpolate="true", alpha=.6, edgecolor = "black")
    end
    # draw substrate medium
    ax_profile[:fill_between]([new_d[end], new_d[end]+doffset], 0, nλ0[end], facecolor=βassigned[:,end], interpolate="true", alpha=.6, edgecolor = "black")

    ax_profile[:tick_params](which="both", direction="in", pad=10, labelsize=22) # ticks offset
    ax_profile[:yaxis][:set_ticks_position]("both")
    ax_profile[:xaxis][:set_ticks_position]("both")
    # axis("tight")
    xlabel("Thickness profile [nm]")
    ylabel(L"\mathrm{Refractive index at }\,\, \lambda_0")

    if size(emf,1) > 1
        if (length(λ) == 1) | (length(θ) == 1)

            figure()
            ax_profile = gca()
            # draw incident medium
            ax_profile[:fill_between]([-doffset, 0], 0, nλ0[1], facecolor=βassigned[:,1], interpolate="true", alpha=0.8, edgecolor = "black")

            new_d = cumsum([0; d], dims=1)
            for idx = 1 : nd
                # draw successive layers
                ax_profile[:fill_between]([new_d[idx], new_d[idx+1]], 0, nλ0[idx+1], facecolor=βassigned[:,idx+1], interpolate="true", alpha=.6, edgecolor = "black")
            end
            # draw substrate medium
            ax_profile[:fill_between]([new_d[end], new_d[end]+doffset], 0, nλ0[end], facecolor=βassigned[:,end], interpolate="true", alpha=.6, edgecolor = "black")

            ax_profile[:tick_params](which="both", direction="in", pad=10, labelsize=22) # ticks offset
            ax_profile[:yaxis][:set_ticks_position]("both")
            ax_profile[:xaxis][:set_ticks_position]("both")
            # axis("tight")
            xlabel("Thickness profile [nm]")
            ylabel(L"\mathrm{Refractive index at }\,\, \lambda_0")

            # plot the emf for the reference wavelength
            aux1 = findmin(abs.(λ .- λ0))[2][1]
            plot(τ, emf[aux1, 1, :])
            title((@sprintf "EMF-resonance at %0.0f [nm]" λ[aux1]), y=1.02)

        end # if (length(λ) == 1) | (length(θ) == 1)
    end # size(emf,1) > 1

end # function nplot(...)


text:
# ADAPTING, MERGING AND FILTERING PRIOR TO PCA
#####################################################################

function merge_beds(parsed_args)
    merge_plink(parsed_args["input"]; des = parsed_args["output"])
end

issnp(x::String) = length(x) == 1

all_batches_ok(row, batchcols) = all(x == 1 for x in row[batchcols])

bounds(pos, lb, ub) =  (min(lb, pos - 10^4), max(ub, pos + 10^4))

function notin_ldblocks(row, ldblocks)
    if (chr=row.chromosome,) ∉ keys(ldblocks)
        return true
    else
        group = ldblocks[(chr=row.chromosome,)]
        inblock = any(b.lower_bound < row.position < b.upper_bound for b in eachrow(group))
        return ~inblock
    end
end

"""
    filter_chromosome(parsed_args)

The purpose of this method is to filter SNPS before applying a
dimensionality reduction technique such as PCA that is later used
to extract population stratification compomemts.
We filter SNPs using quality control metrics from the following resource:
    - https://biobank.ndph.ox.ac.uk/showcase/refer.cgi?id=1955
"""
function filter_chromosome(parsed_args)

    qc_df = CSV.File(parsed_args["qcfile"]) |> DataFrame

    snp_data = SnpData(parsed_args["input"])
    # Load and redefine LD bounds
    ld_blocks = CSV.File(
        parsed_args["ld-blocks"];
        header=["rsid","chr","pos","LDblock_lower","LDblock_upper","LDblock_length","lower_bound","upper_bound"]) |> DataFrame
    ld_blocks.chr = string.(ld_blocks.chr)
    transform!(ld_blocks,
        :,
        [:pos, :lower_bound, :upper_bound] => ByRow(bounds) => [:lower_bound, :upper_bound],
        )
    ld_blocks = groupby(ld_blocks, :chr)

    # Remove SNP's with MAF < maf-threshold
    maf_threshold = parsed_args["maf-threshold"]
    snp_data.snp_info[!, "MAF"] = SnpArrays.maf(snp_data.snparray)
    mafpassed = filter(:MAF => >=(maf_threshold), snp_data.snp_info)

    # Remove LD regions specified by ld_blocks
    ld_pruned = filter(x -> notin_ldblocks(x, ld_blocks), mafpassed)

    # The QC file contains information on fully genotyped SNPS
    # We only keep those
    fully_genotyped_snps = innerjoin(
        ld_pruned,
        qc_df,
        on = :snpid => :rs_id,
        makeunique = true
        )

    # If an RSID appears multiple times, it is because it has
    # more than 2 possible alleles: we remove them
    # (why? maybe because the PCA then cannot tackle them)
    duplicate_rsids = Set(fully_genotyped_snps.snpid[nonunique(fully_genotyped_snps, ["snpid"])])
    biallelic = filter(:snpid=>∉(duplicate_rsids), fully_genotyped_snps)

    # Keep only actual SNPs and not other kinds of variants
    actual_snps = subset(biallelic, :allele1 => ByRow(issnp), :allele2 => ByRow(issnp))

    # All batches pass QC
    batch_cols = [x for x in names(actual_snps) if occursin("Batch", x)]
    batches_ok = filter(row -> all_batches_ok(row, batch_cols), actual_snps)

    # Assayed in both genotyping arrays
    final = filter(:array => ==(2), batches_ok)

    rsids = Set(final.snpid)
    SnpArrays.filter(parsed_args["input"]; des=parsed_args["output"], f_snp = x -> x[:snpid] ∈ rsids)
end

function adapt_flashpca(parsed_args)
    # I think FID and IID are just duplicates
    pcs = CSV.File(parsed_args["input"], drop=["IID"]) |> DataFrame
    rename!(pcs, :FID => :SAMPLE_ID)
    CSV.write(parsed_args["output"], pcs)
end


text:
mutable struct TestClassifier <: AbstractClassifier
    dummy_loss::Float64
    classes::Vector
end

Lighthouse.classes(c::TestClassifier) = c.classes

function Lighthouse.train!(c::TestClassifier, dummy_batches, logger)
    for (dummy_input_batch, loss_delta) in dummy_batches
        c.dummy_loss += loss_delta
        Lighthouse.log_value!(logger, "train/loss_per_batch", c.dummy_loss)
    end
    return c.dummy_loss
end

const RNG_LOSS = StableRNG(22)

function Lighthouse.loss_and_prediction(c::TestClassifier, dummy_input_batch)
    dummy_soft_label_batch = rand(RNG_LOSS, length(c.classes), size(dummy_input_batch)[end])
    # Fake a softmax
    dummy_soft_label_batch .= dummy_soft_label_batch ./ sum(dummy_soft_label_batch; dims=1)
    return c.dummy_loss, dummy_soft_label_batch
end

@testset "Multi-class learn!(::TestModel, ...)" begin
    mktempdir() do tmpdir
        model = TestClassifier(1000000.0, ["class_$i" for i in 1:5])
        k, n = length(model.classes), 3
        rng = StableRNG(22)
        train_batches = [(rand(rng, 4 * k, n), -rand(rng)) for _ in 1:100]
        test_batches = [((rand(rng, 4 * k, n),), (n * i - n + 1):(n * i)) for i in 1:10]
        possible_vote_labels = collect(0:k)
        votes = [rand(rng, possible_vote_labels) for sample in 1:(n * 10), voter in 1:7]
        votes[:, [1, 2, 3]] .= votes[:, 4] # Voter 1-3 voted identically to voter 4 (force non-zero agreement)
        logger = LearnLogger(joinpath(tmpdir, "logs"), "test_run")
        limit = 5
        let counted = 0
            upon_loss_decrease = Lighthouse.upon(logger,
                                                 "test_set_prediction/mean_loss_per_epoch";
                                                 condition=<, initial=Inf)
            callback = n -> begin
                upon_loss_decrease() do _
                    counted += n
                    @info counted n
                end
            end
            elected = majority.((rng,), eachrow(votes), (1:length(Lighthouse.classes(model)),))
            Lighthouse.learn!(model, logger, () -> train_batches, () -> test_batches, votes, elected;
                              epoch_limit=limit, post_epoch_callback=callback)
            @test counted == sum(1:limit)
        end
        @test length(logger.logged["train/loss_per_batch"]) == length(train_batches) * limit
        for key in ["test_set_prediction/loss_per_batch",
             "test_set_prediction/time_in_seconds_per_batch",
             "test_set_prediction/gc_time_in_seconds_per_batch",
             "test_set_prediction/allocations_per_batch",
             "test_set_prediction/memory_in_mb_per_batch"]
            @test length(logger.logged[key]) == length(test_batches) * limit
        end
        for key in ["test_set_prediction/mean_loss_per_epoch",
             "test_set_evaluation/time_in_seconds_per_epoch",
             "test_set_evaluation/gc_time_in_seconds_per_epoch",
             "test_set_evaluation/allocations_per_epoch",
             "test_set_evaluation/memory_in_mb_per_epoch"]
            @test length(logger.logged[key]) == limit
        end
        @test length(logger.logged["test_set_evaluation/metrics_per_epoch"]) == limit

        # Test multiclass optimal_threshold param invalid
        @test_throws ArgumentError Lighthouse.learn!(model, logger, () -> train_batches,
                                                     () -> test_batches, votes;
                                                     epoch_limit=limit,
                                                     optimal_threshold_class=1)

        # Test `predict!`
        num_samples = sum(b -> size(b[1][1], 2), test_batches)
        predicted_soft = zeros(num_samples, length(model.classes))
        predict!(model, predicted_soft, test_batches, logger; logger_prefix="halloooo")
        @test !all(predicted_soft .== 0)
        @test length(logger.logged["halloooo/mean_loss_per_epoch"]) == 1
        @test length(logger.logged["halloooo/loss_per_batch"]) == length(test_batches)
        @test length(logger.logged["halloooo/time_in_seconds_per_batch"]) ==
              length(test_batches)

        # Test `evaluate!`
        n_examples = 40
        num_voters = 20
        predicted_soft = rand(rng, Float32, n_examples, length(model.classes))
        predicted_hard = map(label -> Lighthouse.onecold(model, label),
                             eachrow(predicted_soft))
        votes = [rand(rng, possible_vote_labels)
                 for sample in 1:n_examples, voter in 1:num_voters]
        votes[:, 3] .= votes[:, 4] # Voter 4 voted identically to voter 3 (force non-zero agreement)
        elected_hard = map(row -> majority(rng, row, 1:length(model.classes)),
                           eachrow(votes))
        evaluate!(predicted_hard, predicted_soft, elected_hard, model.classes, logger;
                  logger_prefix="wheeeeeee", logger_suffix="_for_all_time", votes=votes)
        @test length(logger.logged["wheeeeeee/time_in_seconds_for_all_time"]) == 1
        @test length(logger.logged["wheeeeeee/metrics_for_all_time"]) == 1

        # Test startified eval
        strata = [Set("group $(j % Int(ceil(sqrt(j))))" for j in 1:(i - 1))
                  for i in 1:size(votes, 1)]
        plot, plot_data = evaluation_metrics_plot(predicted_hard, predicted_soft,
                                                  elected_hard, model.classes, 0.0:0.01:1.0;
                                                  votes=votes, strata=strata)
        @test haskey(plot_data, "stratified_kappas")

        # Test plotting
        plot_data = last(logger.logged["test_set_evaluation/metrics_per_epoch"])
        @test isa(plot_data["thresholds"], AbstractVector)

        pr = plot_pr_curves(plot_data["per_class_pr_curves"], plot_data["class_labels"])
        @testplot pr

        roc = plot_roc_curves(plot_data["per_class_roc_curves"],
                              plot_data["per_class_roc_aucs"], plot_data["class_labels"])
        @testplot roc

        # Kappa no IRA
        kappas_no_ira = plot_kappas(vcat(plot_data["multiclass_kappa"], plot_data["per_class_kappas"]),
                        vcat("Multiclass", plot_data["class_labels"]))
        @testplot kappas_no_ira

        # Kappa with IRA
        kappas_ira = plot_kappas(vcat(plot_data["multiclass_kappa"], plot_data["per_class_kappas"]),
                        vcat("Multiclass", plot_data["class_labels"]),
                        vcat(plot_data["multiclass_IRA_kappas"],
                             plot_data["per_class_IRA_kappas"]))
        @testplot kappas_ira

        reliability_calibration = plot_reliability_calibration_curves(plot_data["per_class_reliability_calibration_curves"],
                                                                      plot_data["per_class_reliability_calibration_scores"],
                                                                      plot_data["class_labels"])
        @testplot reliability_calibration

        confusion_row = plot_confusion_matrix(plot_data["confusion_matrix"],
                                              plot_data["class_labels"], :Row)
        @testplot confusion_row

        confusion_col = plot_confusion_matrix(plot_data["confusion_matrix"],
                                              plot_data["class_labels"], :Column)
        @testplot confusion_col

        all_together_2 = evaluation_metrics_plot(plot_data)
        @testplot all_together_2

        #savefig(all_together_2, "/tmp/multiclass.png")
    end
end

@testset "2-class `learn!(::TestModel, ...)`" begin
    mktempdir() do tmpdir
        model = TestClassifier(1000000.0, ["class_$i" for i in 1:2])
        k, n = length(model.classes), 3
        rng = StableRNG(23)
        train_batches = [(rand(rng, 4 * k, n), -rand(rng)) for _ in 1:100]
        test_batches = [((rand(rng, 4 * k, n),), (n * i - n + 1):(n * i)) for i in 1:10]
        possible_vote_labels = collect(0:k)
        votes = [rand(rng, possible_vote_labels) for sample in 1:(n * 10), voter in 1:7]
        votes[:, [1, 2, 3]] .= votes[:, 4] # Voter 1-3 voted identically to voter 4 (force non-zero agreement)
        logger = LearnLogger(joinpath(tmpdir, "logs"), "test_run")
        limit = 5
        let counted = 0
            upon_loss_decrease = Lighthouse.upon(logger,
                                                 "test_set_prediction/mean_loss_per_epoch";
                                                 condition=<, initial=Inf)
            callback = n -> begin
                upon_loss_decrease() do _
                    counted += n
                    @info counted n
                end
            end
            elected = majority.((rng,), eachrow(votes), (1:length(Lighthouse.classes(model)),))
            Lighthouse.learn!(model, logger, () -> train_batches, () -> test_batches, votes, elected;
                              epoch_limit=limit, post_epoch_callback=callback)
            @test counted == sum(1:limit)
        end
        # Binary classification logs some additional metrics
        @test length(logger.logged["test_set_evaluation/spearman_correlation_per_epoch"]) ==
              limit
        plot_data = last(logger.logged["test_set_evaluation/metrics_per_epoch"])
        @test haskey(plot_data, "spearman_correlation")

        # No `optimal_threshold_class` during learning...
        @test !haskey(plot_data, "optimal_threshold")
        @test !haskey(plot_data, "optimal_threshold_class")

        # And now, `optimal_threshold_class` during learning
        elected = majority.((rng,), eachrow(votes), (1:length(Lighthouse.classes(model)),))
        Lighthouse.learn!(model, logger, () -> train_batches, () -> test_batches, votes, elected;
                          epoch_limit=limit, optimal_threshold_class=2,
                          test_set_logger_prefix="validation_set")
        plot_data = last(logger.logged["validation_set_evaluation/metrics_per_epoch"])
        @test haskey(plot_data, "optimal_threshold")
        @test haskey(plot_data, "optimal_threshold_class")
        @test plot_data["optimal_threshold_class"] == 2

        # `optimal_threshold_class` param invalid
        @test_throws ArgumentError Lighthouse.learn!(model, logger, () -> train_batches,
                                                     () -> test_batches, votes;
                                                     epoch_limit=limit,
                                                     optimal_threshold_class=3)

        # Test `evaluate!` for votes, no votes
        n_examples = 40
        num_voters = 20
        predicted_soft = rand(rng, Float32, n_examples, length(model.classes))
        predicted_soft .= predicted_soft ./ sum(predicted_soft; dims=2)
        predicted_hard = map(label -> Lighthouse.onecold(model, label),
                             eachrow(predicted_soft))
        votes = [rand(rng, possible_vote_labels)
                 for sample in 1:n_examples, voter in 1:num_voters]
        votes[:, 3] .= votes[:, 4] # Voter 4 voted identically to voter 3 (force non-zero agreement)
        elected_hard = map(row -> majority(rng, row, 1:length(model.classes)),
                           eachrow(votes))

        evaluate!(predicted_hard, predicted_soft, elected_hard, model.classes, logger;
                  logger_prefix="wheeeeeee", logger_suffix="_for_all_time", votes=nothing)
        plot_data = last(logger.logged["wheeeeeee/metrics_for_all_time"])
        @test !haskey(plot_data, "per_class_IRA_kappas")
        @test !haskey(plot_data, "multiclass_IRA_kappas")

        evaluate!(predicted_hard, predicted_soft, elected_hard, model.classes, logger;
                  logger_prefix="wheeeeeee", logger_suffix="_for_all_time", votes=votes)
        plot_data = last(logger.logged["wheeeeeee/metrics_for_all_time"])
        @test haskey(plot_data, "per_class_IRA_kappas")
        @test haskey(plot_data, "multiclass_IRA_kappas")

        # Test `evaluate` for different optimal_threshold classes
        evaluate!(predicted_hard, predicted_soft, elected_hard, model.classes, logger;
                  logger_prefix="wheeeeeee", logger_suffix="_for_all_time", votes=votes,
                  optimal_threshold_class=1)
        plot_data_1 = last(logger.logged["wheeeeeee/metrics_for_all_time"])
        evaluate!(predicted_hard, predicted_soft, elected_hard, model.classes, logger;
                  logger_prefix="wheeeeeee", logger_suffix="_for_all_time", votes=votes,
                  optimal_threshold_class=2)
        plot_data_2 = last(logger.logged["wheeeeeee/metrics_for_all_time"])

        # The thresholds should not be identical (since they are *inclusive* when applied:
        # values greater than _or equal to_ the threshold are given the class value)
        @test plot_data_1["optimal_threshold"] != plot_data_2["optimal_threshold"]

        # The two threshold options yield different results
        thresh_from_roc = Lighthouse._get_optimal_threshold_from_ROC(plot_data_2["per_class_roc_curves"];
                                                                     thresholds=plot_data_2["thresholds"],
                                                                     class_of_interest_index=plot_data_2["optimal_threshold_class"])
        thresh_from_calibration = Lighthouse._calculate_optimal_threshold_from_discrimination_calibration(predicted_soft,
                                                                                votes;
                                                                                thresholds=plot_data_2["thresholds"],
                                                                                class_of_interest_index=plot_data_2["optimal_threshold_class"]).threshold
        @test !isequal(thresh_from_roc, thresh_from_calibration)
        @test isequal(thresh_from_calibration, plot_data_2["optimal_threshold"])

        # Test that plotting succeeds (no specialization relative to the multi-class tests)
        plot_data = last(logger.logged["validation_set_evaluation/metrics_per_epoch"])
        all_together = evaluation_metrics_plot(plot_data)
        #savefig(all_together, "/tmp/binary.png")
        @testplot all_together
    end
end

@testset "`_calculate_ira_kappas`" begin
    classes = ["roy", "gee", "biv"]
    @test isnothing(Lighthouse._calculate_ira_kappas([1; 1; 1; 1], classes))  # Only one voter...
    @test isnothing(Lighthouse._calculate_ira_kappas([1 0; 1 0; 0 1], classes))  # No observations in common...
end

@testset "Calculate `_spearman_corr`" begin
    # Test failures for...
    # ...nonbinary input
    @test_throws ArgumentError Lighthouse._calculate_spearman_correlation([0.9 0.1], [1],
                                                                          ["oh" "em" "gee"])
    # ...non-softmax input
    @test_throws ArgumentError Lighthouse._calculate_spearman_correlation([0.9 0.9], [1],
                                                                          ["oh" "em"])

    # Test class order doesn't matter
    predicted = [0.1 0.9; 0.3 0.7; 1 0]
    votes = [0 1 1 2; 2 2 0 0; 0 2 2 1]
    predicted_flipped = [0.9 0.1; 0.7 0.3; 0 1]
    votes_flipped = [0 2 2 1; 1 1 0 0; 0 1 1 2]
    @test Lighthouse._calculate_spearman_correlation(predicted, votes, ["oh" "em"]) ==
          Lighthouse._calculate_spearman_correlation(predicted_flipped, votes_flipped,
                                                     ["em" "oh"])

    # Test complete agreement
    predicted_soft = [0.0 1.0; 1.0 0.0; 1.0 0.0]
    votes = [2; 1; 1]
    sp = Lighthouse._calculate_spearman_correlation(predicted_soft, votes, ["oh" "em"])
    @test sp.ρ == 1

    # Test complete disagreement
    votes = [1; 2; 2]
    sp = Lighthouse._calculate_spearman_correlation(predicted_soft, votes, ["oh" "em"])
    @test sp.ρ == -1

    # Test NaN spearman due to unranked input
    votes = [1; 2; 2]
    predicted_soft = [
        0.3 0.7
        0.3 0.7
        0.3 0.7
    ]
    sp = Lighthouse._calculate_spearman_correlation(predicted_soft, votes, ["oh" "em"])
    @test isnan(sp.ρ)

    # Test CI
    sp = Lighthouse._spearman_corr([0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.2],
                                   [0.1, 0.3, 0.2, 0.1, 0.3, 0.1, 0.6])
    @test -1.0 < sp.ci_lower < sp.ρ < sp.ci_upper < 1.0
end

@testset "Discrimination calibration" begin
    # Test single voter has same calibration as "whole" by definition
    votes = [1; 1; 1; 2; 2; 2]
    @test length(Lighthouse._elected_probabilities(votes, 1)) == 6
    @test Lighthouse._elected_probabilities(votes, 1) == [1; 1; 1; 0; 0; 0]

    # Test single voter discrimination calibration
    single_voter_calibration = Lighthouse._calculate_voter_discrimination_calibration(votes;
                                                                                class_of_interest_index=1)
    @test length(single_voter_calibration.mse) == 1

    # Test multi-voter voter discrimination calibration
    votes = [
        0 1 1 1
        1 2 0 0
        2 1 2 2
    ] # Note: voters 3 and 4 have voted identically
    voter_calibration = Lighthouse._calculate_voter_discrimination_calibration(votes;
                                                                               class_of_interest_index=1)
    @test length(voter_calibration.mse) == size(votes, 2)
    @test length(voter_calibration.plot_curve_data) == size(votes, 2)
    @test voter_calibration.mse[1] > voter_calibration.mse[3]
    @test !isequal(voter_calibration.plot_curve_data[1],
                   voter_calibration.plot_curve_data[2])
    @test isequal(voter_calibration.plot_curve_data[3],
                  voter_calibration.plot_curve_data[4])

    # Test predicted voter discrimination calibration
    predicted_soft_labels = [1.0 0.0; 0.0 1.0; 0.0 1.0]
    cal = Lighthouse._calculate_optimal_threshold_from_discrimination_calibration(predicted_soft_labels,
                                                                                votes;
                                                                                thresholds=0.0:0.01:1.0,
                                                                                class_of_interest_index=1)
    @test all(cal.mse .<= voter_calibration.mse)
    cal2 = Lighthouse._calculate_optimal_threshold_from_discrimination_calibration(predicted_soft_labels,
                                                                                votes;
                                                                                thresholds=0.0:0.01:1.0,
                                                                                class_of_interest_index=2)
    @test cal.mse == cal2.mse
    @test cal.plot_curve_data[2] != cal2.plot_curve_data[2]
end

@testset "2-class per_class_confusion_statistics" begin
    predicted_soft_labels = [
        0.51 0.49
        0.49 0.51
        0.1 0.9
        0.9 0.1
        0.0 1.0
    ]
    elected_hard_labels = [1, 2, 2, 2, 1]
    thresholds = [0.25, 0.5, 0.75]
    class_1, class_2 = Lighthouse.per_class_confusion_statistics(predicted_soft_labels,
                                                                 elected_hard_labels,
                                                                 thresholds)
    stats_1, stats_2 = class_1[1], class_2[1] # threshold == 0.25
    @test stats_1.actual_negatives == stats_2.actual_positives == 3
    @test stats_1.actual_positives == stats_2.actual_negatives == 2
    @test stats_1.predicted_positives == 3
    @test stats_2.predicted_positives == 4
    @test stats_1.predicted_negatives == 2
    @test stats_2.predicted_negatives == 1
    @test stats_1.true_positives == 1
    @test stats_2.true_positives == 2
    @test stats_1.true_negatives == 1
    @test stats_2.true_negatives == 0
    @test stats_1.false_positives == 2
    @test stats_2.false_positives == 2
    @test stats_1.false_negatives == 1
    @test stats_2.false_negatives == 1
    @test stats_1.true_positive_rate == 0.5
    @test stats_2.true_positive_rate == 2 / 3
    @test stats_1.true_negative_rate == 1 / 3
    @test stats_2.true_negative_rate == 0.0
    @test stats_1.false_positive_rate == 2 / 3
    @test stats_2.false_positive_rate == 1.0
    @test stats_1.false_negative_rate == 0.5
    @test stats_2.false_negative_rate == 1 / 3
    @test stats_1.precision == 1 / 3
    @test stats_2.precision == 0.5

    stats_1, stats_2 = class_1[2], class_2[2] # threshold == 0.5
    @test stats_1.actual_negatives == stats_2.actual_positives == 3
    @test stats_1.actual_positives == stats_2.actual_negatives == 2
    @test stats_1.predicted_negatives == stats_2.predicted_positives == 3
    @test stats_1.predicted_positives == stats_2.predicted_negatives == 2
    @test stats_1.true_negatives == stats_2.true_positives == 2
    @test stats_1.true_positives == stats_2.true_negatives == 1
    @test stats_1.false_positives == stats_2.false_negatives == 1
    @test stats_1.false_negatives == stats_2.false_positives == 1
    @test stats_1.true_positive_rate == stats_2.true_negative_rate == 0.5
    @test stats_1.true_negative_rate == stats_2.true_positive_rate == 2 / 3
    @test stats_1.false_positive_rate == stats_2.false_negative_rate == 1 / 3
    @test stats_1.false_negative_rate == stats_2.false_positive_rate == 0.5
    @test stats_1.precision == 0.5 && stats_2.precision == 2 / 3

    stats_1, stats_2 = class_1[3], class_2[3] # threshold == 0.75
    @test stats_1.actual_negatives == stats_2.actual_positives == 3
    @test stats_1.actual_positives == stats_2.actual_negatives == 2
    @test stats_1.predicted_positives == 1
    @test stats_2.predicted_positives == 2
    @test stats_1.predicted_negatives == 4
    @test stats_2.predicted_negatives == 3
    @test stats_1.true_positives == 0
    @test stats_2.true_positives == 1
    @test stats_1.true_negatives == 2
    @test stats_2.true_negatives == 1
    @test stats_1.false_positives == 1
    @test stats_2.false_positives == 1
    @test stats_1.false_negatives == 2
    @test stats_2.false_negatives == 2
    @test stats_1.true_positive_rate == 0.0
    @test stats_2.true_positive_rate == 1 / 3
    @test stats_1.true_negative_rate == 2 / 3
    @test stats_2.true_negative_rate == 0.5
    @test stats_1.false_positive_rate == 1 / 3
    @test stats_2.false_positive_rate == 0.5
    @test stats_1.false_negative_rate == 1.0
    @test stats_2.false_negative_rate == 2 / 3
    @test stats_1.precision == 0.0
    @test stats_2.precision == 0.5
end

@testset "3-class per_class_confusion_statistics" begin
    predicted_soft_labels = [
        1/3 1/3 1/3
        0.1 0.7 0.2
        0.25 0.25 0.5
        0.4 0.5 0.1
        0.0 0.0 1.0
        0.2 0.5 0.3
        0.5 0.4 0.1
    ]
    elected_hard_labels = [1, 2, 2, 1, 3, 3, 1]
    # TODO would be more robust to have multiple thresholds, but our naive tests
    # here will have to be refactored to avoid becoming a nightmare if we do that
    thresholds = [0.5]
    class_1, class_2, class_3 = Lighthouse.per_class_confusion_statistics(predicted_soft_labels,
                                                                          elected_hard_labels,
                                                                          thresholds)
    stats_1, stats_2, stats_3 = class_1[], class_2[], class_3[] # threshold == 0.5
    @test stats_1.predicted_positives == 1
    @test stats_2.predicted_positives == 3
    @test stats_3.predicted_positives == 2
    @test stats_1.predicted_negatives == 6
    @test stats_2.predicted_negatives == 4
    @test stats_3.predicted_negatives == 5
    @test stats_1.actual_positives == 3
    @test stats_2.actual_positives == 2
    @test stats_3.actual_positives == 2
    @test stats_1.actual_negatives == 4
    @test stats_2.actual_negatives == 5
    @test stats_3.actual_negatives == 5
    @test stats_1.true_positives == 1
    @test stats_2.true_positives == 1
    @test stats_3.true_positives == 1
    @test stats_1.true_negatives == 4
    @test stats_2.true_negatives == 3
    @test stats_3.true_negatives == 4
    @test stats_1.false_positives == 0
    @test stats_2.false_positives == 2
    @test stats_3.false_positives == 1
    @test stats_1.false_negatives == 2
    @test stats_2.false_negatives == 1
    @test stats_3.false_negatives == 1
    @test stats_1.true_positive_rate == 1 / 3
    @test stats_2.true_positive_rate == 0.5
    @test stats_3.true_positive_rate == 0.5
    @test stats_1.true_negative_rate == 1.0
    @test stats_2.true_negative_rate == 0.6
    @test stats_3.true_negative_rate == 0.8
    @test stats_1.false_positive_rate == 0.0
    @test stats_2.false_positive_rate == 0.4
    @test stats_3.false_positive_rate == 0.2
    @test stats_1.false_negative_rate == 2 / 3
    @test stats_2.false_negative_rate == 0.5
    @test stats_3.false_negative_rate == 0.5
    @test stats_1.precision == 1.0
    @test stats_2.precision == 1 / 3
    @test stats_3.precision == 0.5
end


text:
#InspectDR:
#-------------------------------------------------------------------------------
#=
TAGS:
	#WANTCONST, HIDEWARN_0.7
=#

const _devmode = false
if _devmode
@warn("Development mode - precompile deactivated.")
end
__precompile__(!_devmode)

module InspectDR

using Colors
using Graphics
using NumericIO
using Pkg
import Printf: @sprintf
import Cairo
import Cairo: CairoContext

const GtkAvailable = true
import Gtk
const _Gtk = Gtk.ShortNames
#=
REMINDER:
Conditional inclusion causes issues with __precompile__() (source: tkelman).
Not worth risking precompile issues just for current troubles with JuliaBox/Gtk.

try
	import Gtk
	eval(:(const _Gtk = Gtk.ShortNames))
	GtkAvailable = true
catch
	@warn("InspectDR: Error loading Gtk.  GUI-based features are unavailable.")
end
=#

#==Constants: Initial default values
===============================================================================#
const DTPPOINTS_PER_INCH = 72 #Typography (desktop publishing) "points per inch"

#Default font:
const DEFAULT_FONTNAME = (@static Sys.iswindows() ? "Cambria" : "Serif")
#Cairo "built-in": Serif, Sans, Serif, Fantasy, Monospace
#NOTE: "Serif" does not work well in Windows.

#==Aliases
===============================================================================#
NullOr{T} = Union{Nothing, T}

include("codegen.jl")
include("styles.jl")
include("math.jl")
include("math_graphics.jl")
include("math_coordinates.jl")
include("numericfmt.jl")
include("events.jl")
include("datasets.jl")
include("grids.jl")
include("base.jl")
include("stylesheets.jl")
include("defaults.jl")
include("graph2d.jl")
include("glyphs.jl")
include("cairo_ext.jl")
include("cairo_base.jl")
include("cairo_annotation.jl")
include("cairo_smithpolar.jl")
include("cairo_legends.jl")
include("cairo_top.jl")
include("cairo_io.jl")
include("templates.jl")
include("compat.jl")

if GtkAvailable
include("gtk_base.jl")
include("gtk_events.jl")
include("gtk_zoom.jl")
include("gtk_markers.jl")
include("gtk_traces.jl")
include("gtk_input.jl")
include("gtk_top.jl")
end

include("show.jl")

#==Comments
================================================================================
-SVG MIME output (using show) does not show up properly in Notebook.  There
 appears to be an issue with determining the extents of the image (bounding box).
-Data area referred to as "graph" area.  This might not be correct.

DReal vs Float64:
	Float64: Used in API dealing with Gtk/Cairo coordinates.
	DReal: Used in API dealing with plot data itself
		(Attempt to reduce code bloat/different code paths by using concrete
		data type in internal data structures/algorithms... despite potential
		inefficiencies).

Transformations between coordinate systems:
         nlxfrm      lxfrm
	data/ -----> axis -----> device
	world                    /view

	where:
		lxfrm is a linear transformation
		nlxfrm is potentially a nonlinear transformation

	- Data coordinates hold whichever units is associated with the data (D).
	- Axis units are transformations of D (ex: D, log(D), dB20(D), ...)
	- Device units can typically be thought of in pixels.

Auxiliary transformations:
	     nlxfrm
	axis -----> readable
                (read)
	Readable coordinates are basically axis coordinates in a readable form.
	Ex (for log10-X scale): X{axis}=2.5 --> X{read}=10^2.5

PExtents2D vs BoundingBox:
	PExtents2D: Extents of data (data coordinate system).
	BoundingBox: Extents of graphical element ("device" coordinate system).

About extents:
	Extents are typically stored using user-level coordinates (USER).
	When drawing the plot, extents are manipulated in a transformed
	coordinate system (XFRM).  For example, XFRM coordinates of a log-x
	plot are given by XFRM=log10(USER).

	ext_data::PExtents2D #Maximum extents of data
	ext_full::PExtents2D #Defines "full" zoom when combined with ext_data
		#Set fields to NaN to keep ext_data values
	ext::PExtents2D #Current/active extents

	setextents(::Plot2D, ::PExtents2D): Also invalidates rendered plot.
	getextents(...): Complementary accessor.

About bounds:
	plotbounds(...): Bounding box of entire plot.
	databounds(...): Bounding box of entire data area (multiple graph strips).
	graphbounds(...): Bounding box of individual graph strip.
=#

#==Exported interface
===============================================================================#
add = _add #Danger: high risk of collision (common name)
export add, line, glyph
export vmarker, hmarker, atext

#==Unexported interface
================================================================================
Supported plot types:
	InspectDR.Multiplot() #Figure supporting multiple plots
	InspectDR.Plot2D <: Plot #2D plot object

Convenience functions:
	write_png(path, ::Plot)
	write_svg(path, ::Plot)
	write_eps(path, ::Plot)
	write_pdf(path, ::Plot)

Other:
	refresh(::PlotWidget)
	refresh(::GtkPlot)
	clearsubplots(::GtkPlot)
=#

#Create new plot window (Is this a good idea?):
#	figure(args...; kwargs...) = GtkPlot(args...; kwargs...)

#Accessors:
if GtkAvailable
	Plot(pwidget::PlotWidget) = pwidget.src
	Plot(gplot::GtkPlot, i::Int) = gplot.subplots[i]
end

#==Already exported
================================================================================
#Displaying a plot object:
	Base.display(d::InspectDR.GtkDisplay, p::Plot) #Gtk mutltiplot window

#Writing plot to IO stream:
	Base.show(::IO, ::MIME, ::Plot2D)
for the following MIMEs:
		MIME"image/png"
		MIME"image/svg+xml"
		MIME"image/eps"
		MIME"application/pdf"
=#

function __init__()
	global defaults
	_initialize(defaults)

	if GtkAvailable
		initialize_cursors()
		_setdefaults(keybindings)
	end

	checkcompat_plots()
end

#include("precompile.jl")

end #module