text:
# Copy the specified kind of state's variables into an `MArray`.
# This is GPU-safe.
function extract_state(bl, state, ijk, e, st::AbstractStateType)
    FT = eltype(state)
    num_state = number_states(bl, st)
    local_state = MArray{Tuple{num_state}, FT}(undef)
    for s in 1:num_state
        local_state[s] = state[ijk, s, e]
    end
    return Vars{vars_state(bl, st, FT)}(local_state)
end

# Return `true` if the specified symbol is a type name that is a subtype
# of `BalanceLaw` and `false` otherwise.
function isa_bl(sym::Symbol)
    symstr = String(sym)
    bls = map(bl -> String(Symbol(bl)), subtypes(BalanceLaw))
    any(bl -> (isequal(bl, symstr) || endswith(bl, "." * symstr)), bls)
end
isa_bl(ex) = false

uppers_in(s) = foldl((f, c) -> isuppercase(c) ? f * c : f, s, init = "")

# Return a name for the array generated for storing the values of the
# diagnostic variables of `dvtype`.
function dvar_array_name(::InterpolationType, dvtype)
    dvt_short = uppers_in(split(String(Symbol(dvtype)), ".")[end])
    Symbol("vars_", dvt_short, "_array")
end
function dvar_array_name(::InterpolateBeforeCollection, dvtype)
    dvt_short = uppers_in(split(String(Symbol(dvtype)), ".")[end])
    Symbol("acc_vars_", dvt_short, "_array")
end

# Generate the common definitions used in many places.
function generate_common_defs()
    quote
        mpicomm = DiagnosticsMachine.Settings.mpicomm
        dg = DiagnosticsMachine.Settings.dg
        Q = DiagnosticsMachine.Settings.Q
        mpirank = MPI.Comm_rank(mpicomm)
        bl = dg.balance_law
        grid = dg.grid
        grid_info = basic_grid_info(dg)
        topl_info = basic_topology_info(grid.topology)
        Nq1 = grid_info.Nq[1]
        Nq2 = grid_info.Nq[2]
        Nqh = grid_info.Nqh
        Nqk = grid_info.Nqk
        npoints = prod(grid_info.Nq)
        nrealelem = topl_info.nrealelem
        nvertelem = topl_info.nvertelem
        nhorzelem = topl_info.nhorzrealelem
        FT = eltype(Q)
        interpol = dgngrp.interpol
        params = dgngrp.params
    end
end

# Generate the `dims` dictionary for `Writers.init_data`.
function generate_init_dims(::NoInterpolation, cfg_type, dvtype_dvars_map)
    dimslst = Any[]
    for dvtype in keys(dvtype_dvars_map)
        dimnames = DiagnosticsMachine.dv_dg_dimnames(cfg_type, dvtype)
        dimranges = DiagnosticsMachine.dv_dg_dimranges(cfg_type, dvtype)
        for (dimname, dimrange) in zip(dimnames, dimranges)
            lhs = :($dimname)
            rhs = :(collect($dimrange), Dict())
            push!(dimslst, :($lhs => $rhs))
        end
    end

    quote
        OrderedDict($(Expr(:tuple, dimslst...))...)
    end
end
function generate_init_dims(::InterpolationType, cfg_type, dvtype_dvars_map)
    quote
        dims = dimensions(interpol)
        if interpol isa InterpolationCubedSphere
            # Adjust `level` on the sphere.
            level_val = dims["level"]
            dims["level"] = (
                level_val[1] .-
                FT(planet_radius(DiagnosticsMachine.Settings.param_set)),
                level_val[2],
            )
        end
        dims
    end
end

get_dimnames(::NoInterpolation, cfg_type, dvtype) =
    DiagnosticsMachine.dv_dg_dimnames(cfg_type, dvtype)
get_dimnames(::InterpolationType, cfg_type, dvtype) =
    DiagnosticsMachine.dv_i_dimnames(cfg_type, dvtype)

# Generate the `vars` dictionary for `Writers.init_data`.
function generate_init_vars(intrp, cfg_type, dvtype_dvars_map)
    varslst = Any[]
    for (dvtype, dvlst) in dvtype_dvars_map
        for dvar in dvlst
            lhs = :($(DiagnosticsMachine.dv_name(cfg_type, dvar)))
            dimnames = get_dimnames(intrp, cfg_type, dvtype)
            rhs = :(
                $dimnames,
                FT,
                $(DiagnosticsMachine.dv_attrib(cfg_type, dvar)),
            )
            push!(varslst, :($lhs => $rhs))
        end
    end

    quote
        # TODO: add code to filter this based on what's actually in `bl`.
        OrderedDict($(Expr(:tuple, varslst...))...)
    end
end

# Generate `Diagnostics.$(name)_init(...)` which will initialize the
# `DiagnosticsGroup` when called.
function generate_init_fun(
    name,
    config_type,
    params_type,
    init_fun,
    interpolate,
    dvtype_dvars_map,
)
    init_name = Symbol(name, "_init")
    cfg_type_name = getfield(ConfigTypes, config_type)
    intrp = getfield(@__MODULE__, interpolate)
    quote
        function $init_name(dgngrp, curr_time)
            $(generate_common_defs())

            $(init_fun)(dgngrp, curr_time)

            # TODO: uncomment when old diagnostics groups are removed
            #if dgngrp.onetime
            DiagnosticsMachine.collect_onetime(mpicomm, dg, Q)
            #end

            if mpirank == 0
                dims =
                    $(generate_init_dims(
                        intrp(),
                        cfg_type_name(),
                        dvtype_dvars_map,
                    ))
                vars =
                    $(generate_init_vars(
                        intrp(),
                        cfg_type_name(),
                        dvtype_dvars_map,
                    ))

                # create the output file
                dprefix = @sprintf("%s_%s", dgngrp.out_prefix, dgngrp.name)
                dfilename =
                    joinpath(DiagnosticsMachine.Settings.output_dir, dprefix)
                noov = DiagnosticsMachine.Settings.no_overwrite
                init_data(dgngrp.writer, dfilename, noov, dims, vars)
            end

            return nothing
        end
    end
end

# Generate code snippet for copying arrays to the CPU if needed. Ideally,
# this will be removed when diagnostics are made to run on GPU.
function generate_array_copies()
    quote
        # get needed arrays onto the CPU
        if array_device(Q) isa CPU
            prognostic_data = Q.realdata
            gradient_flux_data = dg.state_gradient_flux.realdata
            auxiliary_data = dg.state_auxiliary.realdata
            vgeo = grid.vgeo
        else
            prognostic_data = Array(Q.realdata)
            gradient_flux_data = Array(dg.state_gradient_flux.realdata)
            auxiliary_data = Array(dg.state_auxiliary.realdata)
            vgeo = Array(grid.vgeo)
        end
    end
end

# Generate code to create the necessary arrays for the diagnostics
# variables.
function generate_create_vars_arrays(
    ::InterpolateBeforeCollection,
    cfg_type,
    dvtype_dvars_map,
)
    quote end
end
function generate_create_vars_arrays(
    intrp::InterpolationType,
    cfg_type,
    dvtype_dvars_map,
)
    cva_exs = []
    for (dvtype, dvlst) in dvtype_dvars_map
        arr_name = dvar_array_name(intrp, dvtype)
        npoints = DiagnosticsMachine.dv_dg_points_length(cfg_type, dvtype)
        nvars = length(dvlst)
        nelems = DiagnosticsMachine.dv_dg_elems_length(cfg_type, dvtype)
        cva_ex = quote
            $arr_name = Array{FT}(undef, $npoints, $nvars, $nelems)
            fill!($arr_name, 0)
        end
        push!(cva_exs, cva_ex)
    end
    return Expr(:block, (cva_exs...))
end

# Generate the LHS of the assignment expression into the diagnostic
# variables array.
function dvar_array_assign_expr(
    ::InterpolateBeforeCollection,
    cfg_type,
    dvtype,
    dvidx,
    arr_name,
)
    return :($(arr_name)[x, y, z, $dvidx])
end
function dvar_array_assign_expr(
    ::InterpolationType,
    cfg_type,
    dvtype,
    dvidx,
    arr_name,
)
    pt = DiagnosticsMachine.dv_dg_points_index(cfg_type, dvtype)
    elem = DiagnosticsMachine.dv_dg_elems_index(cfg_type, dvtype)
    :($(arr_name)[$pt, $dvidx, $elem])
end

# Generate calls to the implementations for the `DiagnosticVar`s in this
# group and store the results.
function generate_collect_calls(
    intrp::InterpolationType,
    cfg_type,
    dvtype_dvars_map,
)
    cc_exs = []
    for (dvtype, dvlst) in dvtype_dvars_map
        arr_name = dvar_array_name(intrp, dvtype)
        dvt = split(String(Symbol(dvtype)), ".")[end]
        var_impl = Symbol("dv_", dvt)

        for (v, dvar) in enumerate(dvlst)
            lhs_ex =
                dvar_array_assign_expr(intrp, cfg_type, dvtype, v, arr_name)
            impl_args = DiagnosticsMachine.dv_args(cfg_type, dvar)
            AT1 = impl_args[1][2] # the type of the first argument
            if isa_bl(AT1)
                impl_extra_params = ()
            else
                AT2 = impl_args[2][2] # the type of the second argument
                @assert isa_bl(AT2)
                AT1 = impl_args[1][2] # the type of the first argument
                impl_extra_params = (:(BalanceLaws.sub_model(bl, $AT1)),)
            end
            cc_ex = DiagnosticsMachine.dv_op(
                cfg_type,
                dvtype,
                lhs_ex,
                :(DiagnosticsMachine.$(var_impl)(
                    $cfg_type,
                    $dvar,
                    $(impl_extra_params...),
                    bl,
                    states,
                    curr_time,
                    cache,
                )),
            )
            push!(cc_exs, cc_ex)
        end
    end

    return Expr(:block, (cc_exs...))
end

# Generate the nested loops to traverse the DG grid within which we extract
# the various states and then generate the individual collection calls.
function generate_dg_collections(
    ::InterpolateBeforeCollection,
    cfg_type,
    dvtype_dvars_map,
)
    quote end
end
function generate_dg_collections(
    intrp::InterpolationType,
    cfg_type,
    dvtype_dvars_map,
)
    quote
        cache = Dict{Symbol, Any}()
        for eh in 1:nhorzelem, ev in 1:nvertelem
            e = ev + (eh - 1) * nvertelem
            for k in 1:Nqk, j in 1:Nq2, i in 1:Nq1
                ijk = i + Nq1 * ((j - 1) + Nq2 * (k - 1))
                MH = vgeo[ijk, grid.MHid, e]
                prognostic = DiagnosticsMachine.extract_state(
                    bl,
                    prognostic_data,
                    ijk,
                    e,
                    Prognostic(),
                )
                gradient_flux = DiagnosticsMachine.extract_state(
                    bl,
                    gradient_flux_data,
                    ijk,
                    e,
                    GradientFlux(),
                )
                auxiliary = DiagnosticsMachine.extract_state(
                    bl,
                    auxiliary_data,
                    ijk,
                    e,
                    Auxiliary(),
                )
                states = States(prognostic, gradient_flux, auxiliary)
                $(generate_collect_calls(intrp, cfg_type, dvtype_dvars_map))
                empty!(cache)
            end
        end
    end
end

# Generate any reductions needed for the data collected thus far.
function generate_dg_reductions(
    intrp::NoInterpolation,
    cfg_type,
    dvtype_dvars_map,
)
    red_exs = []
    for (dvtype, dvlst) in dvtype_dvars_map
        arr_name = dvar_array_name(intrp, dvtype)
        red_ex = DiagnosticsMachine.dv_reduce(cfg_type, dvtype, arr_name)
        push!(red_exs, red_ex)
    end

    return Expr(:block, (red_exs...))
end
function generate_dg_reductions(::InterpolationType, cfg_type, dvtype_dvars_map)
    quote end
end

# Generate code to perform density averaging, as needed.
function generate_density_averaging(
    intrp::InterpolationType,
    cfg_type,
    dvtype_dvars_map,
)
    da_exs = []
    for (dvtype, dvlst) in dvtype_dvars_map
        arr_name = dvar_array_name(intrp, dvtype)

        for (v, dvar) in enumerate(dvlst)
            scale_dvar = DiagnosticsMachine.dv_scale(cfg_type, dvar)
            if isnothing(scale_dvar)
                continue
            end
            sdv_idx = findfirst(dvar -> scale_dvar == dvar, dvlst)
            @assert !isnothing(sdv_idx)
            da_ex = quote
                $(arr_name)[:, $v, :] ./= $(arr_name)[:, $sdv_idx, :]
            end
            push!(da_exs, da_ex)
        end
    end

    return Expr(:block, (da_exs...))
end

# Generate interpolation calls as needed. None for `NoInterpolation`.
function generate_interpolations(::NoInterpolation, cfg_type, dvtype_dvars_map)
    quote end
end
# Interpolate only the diagnostic variables arrays.
function generate_interpolations(
    intrp::InterpolateAfterCollection,
    cfg_type,
    dvtype_dvars_map,
)
    ic_exs = []
    for (dvtype, dvlst) in dvtype_dvars_map
        nvars = length(dvlst)
        arr_name = dvar_array_name(intrp, dvtype)
        iarr_name = Symbol("i", arr_name)
        acc_arr_name = Symbol("acc_", arr_name)
        i_pr = findall(
            dvar -> DiagnosticsMachine.dv_project(cfg_type, dvar),
            dvlst,
        )
        ic_ex = quote
            # TODO: `interpolate_local!` expects arrays of the same type
            # as `grid`. This should be solved when we make this GPU-capable.
            if !(array_device(Q) isa CPU)
                ArrayType = arraytype(grid)
                $arr_name = ArrayType($arr_name)
            end
            $iarr_name = similar($arr_name, interpol.Npl, $nvars)
            interpolate_local!(interpol, $arr_name, $iarr_name)

            if interpol isa InterpolationCubedSphere
                project_cubed_sphere!(
                    interpol,
                    $iarr_name,
                    tuple(collect($i_pr)...),
                )
            end

            $acc_arr_name =
                accumulate_interpolated_data(mpicomm, interpol, $iarr_name)
            if !(array_device(Q) isa CPU)
                $acc_arr_name = Array($acc_arr_name)
            end
        end
        push!(ic_exs, ic_ex)
    end

    return Expr(:block, (ic_exs...))
end
# Interpolate all the arrays needed for `States`.
function generate_interpolations(
    ::InterpolateBeforeCollection,
    cfg_type,
    dvtype_dvars_map,
)
    quote
        iprognostic_array =
            similar(Q.realdata, interpol.Npl, number_states(bl, Prognostic()))
        interpolate_local!(interpol, Q.realdata, iprognostic_array)
        igradient_flux_array =
            similar(Q.realdata, interpol.Npl, number_states(bl, GradientFlux()))
        interpolate_local!(
            interpol,
            dg.state_gradient_flux.realdata,
            igradient_flux_array,
        )
        iauxiliary_array =
            similar(Q.realdata, interpol.Npl, number_states(bl, Auxiliary()))
        interpolate_local!(
            interpol,
            dg.state_auxiliary.realdata,
            iauxiliary_array,
        )

        if interpol isa InterpolationCubedSphere
            i_ρu = varsindex(vars_state(bl, Prognostic(), FT), :ρu)
            project_cubed_sphere!(
                interpol,
                iprognostic_array,
                tuple(collect(i_ρu)...),
            )
        end

        # FIXME: accumulating to rank 0 is not scalable
        all_prognostic_data =
            accumulate_interpolated_data(mpicomm, interpol, iprognostic_array)
        all_gradient_flux_data = accumulate_interpolated_data(
            mpicomm,
            interpol,
            igradient_flux_array,
        )
        all_auxiliary_data =
            accumulate_interpolated_data(mpicomm, interpol, iauxiliary_array)
    end
end

# Generate code to create the necessary arrays to collect the diagnostics
# variables on the interpolated grid.
function generate_create_i_vars_arrays(
    intrp::InterpolateBeforeCollection,
    cfg_type,
    dvtype_dvars_map,
)
    cva_exs = []
    for (dvtype, dvlst) in dvtype_dvars_map
        arr_name = dvar_array_name(intrp, dvtype)
        nvars = length(dvlst)
        cva_ex = quote
            $arr_name = Array{FT}(
                undef,
                tuple(map(d -> length(d[1]), values(dims))..., $nvars),
            )
            fill!($arr_name, 0)
        end
        push!(cva_exs, cva_ex)
    end

    return Expr(:block, (cva_exs...))
end
function generate_create_i_vars_arrays(
    ::InterpolationType,
    cfg_type,
    dvtype_dvars_map,
)
    quote end
end

# Generate the nested loops to traverse the interpolated grid within
# which we extract the various (interpolated) states and then generate
# the individual collection calls.
function generate_i_collections(
    intrp::InterpolateBeforeCollection,
    cfg_type,
    dvtype_dvars_map,
)
    quote
        cache = Dict{Symbol, Any}()
        (x1, x2, x3) = map(d -> length(d[1]), values(dims))
        for x in 1:x1, y in 1:x2, z in 1:x3
            iprognostic = Vars{vars_state(bl, Prognostic(), FT)}(view(
                all_prognostic_data,
                x,
                y,
                z,
                :,
            ))
            igradient_flux = Vars{vars_state(bl, GradientFlux(), FT)}(view(
                all_gradient_flux_data,
                x,
                y,
                z,
                :,
            ))
            iauxiliary = Vars{vars_state(bl, Auxiliary(), FT)}(view(
                all_auxiliary_data,
                x,
                y,
                z,
                :,
            ))
            states = States(iprognostic, igradient_flux, iauxiliary)
            $(generate_collect_calls(intrp, cfg_type, dvtype_dvars_map))
            empty!(cache)
        end
    end
end
function generate_i_collections(::InterpolationType, cfg_type, dvtype_dvars_map)
    quote end
end

# Generate assignments into `varvals` for writing.
function generate_varvals(intrp::NoInterpolation, cfg_type, dvtype_dvars_map)
    vv_exs = []
    for (dvtype, dvlst) in dvtype_dvars_map
        arr_name = dvar_array_name(intrp, dvtype)
        for (v, dvar) in enumerate(dvlst)
            rhs = :(view($(arr_name), :, $v, :))
            if length(DiagnosticsMachine.dv_dg_dimnames(cfg_type, dvtype)) == 1
                rhs = :(reshape($(rhs), :))
            end
            vv_ex = quote
                varvals[$(DiagnosticsMachine.dv_name(cfg_type, dvar))] = $rhs
            end
            push!(vv_exs, vv_ex)
        end
    end

    return Expr(:block, (vv_exs...))
end
function generate_varvals(::InterpolationType, cfg_type, dvtype_dvars_map)
    vv_exs = []
    for (dvtype, dvlst) in dvtype_dvars_map
        dvt_short = uppers_in(split(String(Symbol(dvtype)), ".")[end])
        acc_arr_name = Symbol("acc_vars_", dvt_short, "_array")
        for (v, dvar) in enumerate(dvlst)
            vv_ex = quote
                varvals[$(DiagnosticsMachine.dv_name(cfg_type, dvar))] =
                    $(acc_arr_name)[:, :, :, $v]
            end
            push!(vv_exs, vv_ex)
        end
    end

    return Expr(:block, (vv_exs...))
end

# Generate `Diagnostics.$(name)_collect(...)` which when called,
# performs a collection of all the diagnostic variables in the group
# and writes them out.
function generate_collect_fun(
    name,
    config_type,
    params_type,
    init_fun,
    interpolate,
    dvtype_dvars_map,
)
    collect_name = Symbol(name, "_collect")
    cfg_type_name = getfield(ConfigTypes, config_type)
    intrp = getfield(@__MODULE__, interpolate)
    quote
        function $collect_name(dgngrp, curr_time)
            $(generate_common_defs())
            $(generate_array_copies())
            $(generate_create_vars_arrays(
                intrp(),
                cfg_type_name(),
                dvtype_dvars_map,
            ))

            # Traverse the DG grid and collect diagnostics as needed.
            $(generate_dg_collections(
                intrp(),
                cfg_type_name(),
                dvtype_dvars_map,
            ))

            # Perform any reductions necessary.
            $(generate_dg_reductions(
                intrp(),
                cfg_type_name(),
                dvtype_dvars_map,
            ))

            # Perform density averaging if needed.
            $(generate_density_averaging(
                intrp(),
                cfg_type_name(),
                dvtype_dvars_map,
            ))

            # Interpolate and accumulate if needed.
            $(generate_interpolations(
                intrp(),
                cfg_type_name(),
                dvtype_dvars_map,
            ))

            if mpirank == 0
                dims = dimensions(interpol)

                $(generate_create_i_vars_arrays(
                    intrp(),
                    cfg_type_name(),
                    dvtype_dvars_map,
                ))

                # Traverse the interpolated grid and collect diagnostics if needed.
                $(generate_i_collections(
                    intrp(),
                    cfg_type_name(),
                    dvtype_dvars_map,
                ))

                # Assemble the diagnostic variables and write them.
                varvals = OrderedDict()
                $(generate_varvals(intrp(), cfg_type_name(), dvtype_dvars_map))
                append_data(dgngrp.writer, varvals, curr_time)
            end

            MPI.Barrier(mpicomm)
            return nothing
        end
    end
end

# Generate `Diagnostics.$(name)_fini(...)`, which does nothing right now.
function generate_fini_fun(name, args...)
    fini_name = Symbol(name, "_fini")
    quote
        function $fini_name(dgngrp, curr_time) end
    end
end

# Generate `$(name)(...)` which will create the `DiagnosticsGroup` for
# $name when called.
function generate_setup_fun(
    name,
    config_type,
    params_type,
    init_fun,
    interpolate,
    dvtype_dvars_map,
)
    init_name = Symbol(name, "_init")
    collect_name = Symbol(name, "_collect")
    fini_name = Symbol(name, "_fini")

    setup_name = Symbol(name)
    intrp = getfield(@__MODULE__, interpolate)

    no_intrp_err = quote end
    some_intrp_err = quote end
    if intrp() isa NoInterpolation
        no_intrp_err = quote
            @warn "$($name) does not specify interpolation, but an " *
                  "`InterpolationTopology` has been provided; ignoring."
            interpol = nothing
        end
    else
        some_intrp_err = quote
            throw(ArgumentError(
                "$($name) specifies interpolation, but no " *
                "`InterpolationTopology` has been provided.",
            ))
        end
    end
    quote
        function $setup_name(
            interval::String,
            out_prefix::String,
            params::$params_type = nothing;
            writer = NetCDFWriter(),
            interpol = nothing,
        ) where {
            $config_type <: ClimateMachineConfigType,
            $params_type <: Union{Nothing, DiagnosticsGroupParams},
        }
            if isnothing(interpol)
                $(some_intrp_err)
            else
                $(no_intrp_err)
            end

            return DiagnosticsMachine.DiagnosticsGroup(
                $name,
                $init_name,
                $fini_name,
                $collect_name,
                interval,
                out_prefix,
                writer,
                interpol,
                # TODO: uncomment when old diagnostics groups are removed
                #$(intrp() isa NoInterpolation),
                params,
            )
        end
    end
end


text:
module Xsum

export xsum, XAccumulator

using xsum_jll

include("xsum_h.jl")

"""
`acc = XAccumulator()` creates an object `acc` that can be used to accumulate an
exactly rounded double-precision sum.   The sum is accumulated by repeatedly
calling `accumulate!(acc, x)` on summands `x`, and then the final exactly-rounded
result is obtained by calling `float(acc)`.
"""
struct XAccumulator <: Real
    acc::Base.RefValue{xsum_small_accumulator}
    XAccumulator(acc::Base.RefValue{xsum_small_accumulator}) = new(acc)
    function XAccumulator()
        acc = Base.RefValue{xsum_small_accumulator}()
        xsum_small_init(acc)
        return new(acc)
    end
end
XAccumulator(s::XAccumulator) = XAccumulator(Ref(s.acc[]))
Base.float(s::XAccumulator) = xsum_small_round(s.acc)
XAccumulator(init::Real) = accumulate!(XAccumulator(), init)
Base.show(io::IO, s::XAccumulator) = print(io, "XAccumulator(", float(s), ")")

"""
    accumulate!(s::XAccumulator, x::Real)

Add `Float64(x)` to the accumulated sum `s`, with the sum performed in (effectively) infinite precision.
"""
Base.accumulate!(s::XAccumulator, x::Real) = begin; xsum_small_add1(s.acc, x); s; end

function Base.accumulate!(s::XAccumulator, s′::XAccumulator)
    xsum_small_add_accumulator(s.acc, s′.acc)
    return s
end

function Base.sum(a::AbstractVector{XAccumulator})
    s = XAccumulator()
    for s′ in a
        accumulate!(s, s′)
    end
    return s
end

# slower out-of-place variants for completeness
Base.:+(s::XAccumulator, s′::XAccumulator) = accumulate!(XAccumulator(s), s′)
Base.:+(s::XAccumulator, x::Real) = accumulate!(XAccumulator(s), x)
Base.:+(x::Real, s::XAccumulator) = s + x

function Base.accumulate!(s::XAccumulator, a::StridedVector{Float64})
    if stride(a,1) != 1
        for x in a
            accumulate!(s, x)
        end
    else
        n = length(a)
        if n < 256 # empirical threshold for small accumulator to be faster
            xsum_small_addv(s.acc, a)
        else
            acc = Base.RefValue{xsum_small_accumulator}()
            xsum_large_to_small_accumulator(acc, _xsum_large(a))
            xsum_small_add_accumulator(s.acc, acc)
        end
    end
    return s
end

Base.accumulate!(s::Complex{XAccumulator}, x::Complex) = Complex(accumulate!(real(s), real(x)), accumulate!(imag(s), imag(x)))
Base.accumulate!(s::Complex{XAccumulator}, x::Real) = Complex(accumulate!(real(s), x), imag(s))
Base.accumulate!(s::XAccumulator, x::Complex) = Complex(accumulate!(s, real(x)), XAccumulator(imag(x)))
XAccumulator(init::Complex) = Complex(XAccumulator(real(init)), XAccumulator(imag(init)))

"""
    negate!(s::XAccumulator)

Change `s` in-place to negate the sign of the sum it represents.
"""
negate!(s::XAccumulator) = begin xsum_small_negate(s.acc); s; end
Base.:-(s::XAccumulator) = negate!(XAccumulator(s))
Base.:-(s::XAccumulator, s′::XAccumulator) = accumulate!(-s′, s)
Base.:-(s::XAccumulator, x::Real) = accumulate!(XAccumulator(s), -x)
Base.:-(x::Real, s::XAccumulator) = accumulate!(-s, x)

"""
    xsum([f,] itr)

Compute the exactly rounded double-precision sum of calling `f` (which defaults to `identity`) on
the elements of the iterable collection `itr`. That is, return the sum as if (1) the elements of `f.(itr)`
were converted to double (`Float64`) precision, (2) summed in *infinite* precision, and
(3) rounded to the closest double-precision value in the final result.    Both real and complex
sums are supported.
"""
function xsum(f, itr)
    i = iterate(itr)
    if i === nothing
        zero = Base.mapreduce_empty_iter(f, +, itr, Base.IteratorEltype(itr))
        return zero isa Real ? Float64(zero) : zero isa Complex ? ComplexF64(zero) : zero
    end
    v, state = i
    s = XAccumulator(f(v))
    while true
        i = iterate(itr, state)
        i === nothing && return float(s)
        v, state = i
        s = accumulate!(s, f(v))
    end
end

xsum(itr) = xsum(identity, itr)

function _xsum_large(a::StridedVector{Float64})
    @assert stride(a,1) == 1
    acc = Base.RefValue{xsum_large_accumulator}()
    xsum_large_init(acc)
    xsum_large_addv(acc, a)
    return acc
end

function xsum(a::StridedVector{Float64})
    stride(a,1) != 1 && return xsum(identity, a)
    n = length(a)
    if n < 256 # empirical threshold for small accumulator to be faster
        let acc = Base.RefValue{xsum_small_accumulator}()
            xsum_small_init(acc)
            xsum_small_addv(acc, a)
            return xsum_small_round(acc)
        end
    else
        return xsum_large_round(_xsum_large(a))
    end
end

end # module

ext: jl


text:
using Distributions, PDMats, Primes, Random, LinearAlgebra

function qsimvnv(Σ::AbstractMatrix, a, b;m=nothing)
    #= rev 1.13

    This function uses an algorithm given in the paper
    "Numerical Computation of Multivariate Normal Probabilities", in
     J. of Computational and Graphical Stat., 1(1992), pp. 141-149, by
    Alan Genz, WSU Math, PO Box 643113, Pullman, WA 99164-3113
    Email : alangenz@wsu.edu
    The primary references for the numerical integration are
    "On a Number-Theoretical Integration Method"
    H. Niederreiter, Aequationes Mathematicae, 8(1972), pp. 304-11, and
    "Randomization of Number Theoretic Methods for Multiple Integration"
    R. Cranley and T.N.L. Patterson, SIAM J Numer Anal, 13(1976), pp. 904-14.

    Re-coded in Julia from the MATLAB function qsimvnv(m,r,a,b)

    Alan Genz is the author the MATLAB qsimvnv() function.
    Alan Genz software website: http://archive.is/jdeRh
    Source code to MATLAB qsimvnv() function: http://archive.is/h5L37
    % QSIMVNV(m,r,a,b) and _chlrdr(r,a,b)
    %
    % Copyright (C) 2013, Alan Genz,  All rights reserved.
    %
    % Redistribution and use in source and binary forms, with or without
    % modification, are permitted provided the following conditions are met:
    %   1. Redistributions of source code must retain the above copyright
    %      notice, this list of conditions and the following disclaimer.
    %   2. Redistributions in binary form must reproduce the above copyright
    %      notice, this list of conditions and the following disclaimer in
    %      the documentation and/or other materials provided with the
    %      distribution.
    %   3. The contributor name(s) may not be used to endorse or promote
    %      products derived from this software without specific prior
    %      written permission.
    % THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
    % "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
    % LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
    % FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
    % COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
    % INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
    % BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
    % OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
    % ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
    % TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF USE
    % OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
    %

    Julia dependencies
    Distributions
    PDMats
    Primes
    Random
    LinearAlgebra =#

    if isnothing(m)
        m = 1000 * size(Σ, 1)  # default is 1000 * dimension
    end

    # check for proper dimensions
    n = size(Σ, 1)
    nc = size(Σ, 2) 	# assume square Cov matrix nxn
    # check dimension > 1
    n >= 2   || throw(ErrorException("dimension of Σ must be 2 or greater. Σ dimension: $(size(Σ))"))
    n == nc  || throw(DimensionMismatch("Σ matrix must be square. Σ dimension: $(size(Σ))"))

    # check dimensions of lower vector, upper vector, and cov matrix match
    (n == size(a, 1) == size(b, 1)) || throw(DimensionMismatch("inconsistent argument dimensions. Sizes: Σ $(size(Σ))  a $(size(a))  b $(size(b))"))

    # check that a and b are column vectors; if row vectors, fix it
    if size(a, 1) < size(a, 2)
        a = transpose(a)
    end
    if size(b, 1) < size(b, 2)
        b = transpose(b)
    end

    # check that lower integration limit a < upper integration limit b for all elements
    all(a .<= b) || throw(ArgumentError("lower integration limit a must be <= upper integration limit b"))

    # check that Σ is positive definate; if not, print warning
    isa(Σ, AbstractPDMat) || isposdef(Σ) || @warn "covariance matrix Σ fails positive definite check"

    # check if Σ, a, or b contains NaNs
    if any(isnan.(Σ)) || any(isnan.(a)) || any(isnan.(b))
        p = NaN
        e = NaN
        return (p, e)
    end

    # check if a==b
    if a == b
        p = 0.0
        e = 0.0
        return (p, e)
    end

    # check if a = -Inf & b = +Inf
    if all(a .== -Inf) && all(b .== Inf)
        p = 1.0
        e = 0.0
        return (p, e)
    end

    # check input Σ, a, b are floats; otherwise, convert them
    if eltype(Σ) <: Signed
        Σ = float(Σ)
    end

    if eltype(a) <: Signed
        a = float(a)
    end

    if eltype(b) <: Signed
        b = float(b)
    end

    ##################################################################
    #
    # Special cases: positive Orthant probabilities for 2- and
    # 3-dimesional Σ have exact solutions. Integration range [0,∞]
    #
    ##################################################################

    if all(a .== zero(eltype(a))) && all(b .== Inf) && n <= 3
        Σstd = sqrt.(diag(Σ))
        Rcorr = cov2cor(Σ, Σstd)

        if n == 2
            p = 1 / 4 + asin(Rcorr[1,2]) / (2π)
            e = eps()
        elseif n == 3
            p = 1 / 8 + (asin(Rcorr[1,2]) + asin(Rcorr[2,3]) + asin(Rcorr[1,3])) / (4π)
            e = eps()
        end

        return (p, e)

    end

    ##################################################################
    #
    # get lower cholesky matrix and (potentially) re-ordered integration vectors
    #
    ##################################################################

    (ch, as, bs) = _chlrdr(Σ, a, b) # ch =lower cholesky; as=lower vec; bs=upper vec

    ##################################################################
    #
    # quasi-Monte Carlo integration of MVN integral
    #
    ##################################################################

    ### setup initial values
    ai = as[1]
    bi = bs[1]
    ct = ch[1,1]

    unitnorm = Normal() # unit normal distribution
    rng = RandomDevice()

    # if ai is -infinity, explicity set c=0
    # implicitly, the algorith classifies anythign > 9 std. deviations as infinity
    if ai > -9 * ct
        if ai < 9 * ct
            c1 = cdf.(unitnorm, ai / ct)
        else
            c1 = 1.0
        end
    else
        c1 = 0.0
    end

    # if bi is +infinity, explicity set d=0
    if bi > -9 * ct
        if bi < 9 * ct
            d1 = cdf(unitnorm, bi / ct)
        else
            d1 = 1.0
        end
    else
        d1 = 0.0
    end

    # n=size(Σ,1) 	# assume square Cov matrix nxn
    cxi = c1			# initial cxi; genz uses ci but it conflicts with Lin. Alg. ci variable
    dci = d1 - cxi		# initial dcxi
    p = 0.0			# probablity = 0
    e = 0.0			# error = 0

    # Richtmyer generators
    ps = sqrt.(primes(Int(floor(5 * n * log(n + 1) / 4)))) # Richtmyer generators
    q = ps[1:n - 1,1]
    ns = 12
    nv = Int(max(floor(m / ns), 1))

    Jnv    = ones(1, nv)
    cfill  = transpose(fill(cxi, nv)) 	# evaulate at nv quasirandom points row vec
    dpfill = transpose(fill(dci, nv))
    y      = zeros(n - 1, nv)				# n-1 rows, nv columns, preset to zero

    #= Randomization loop for ns samples
     j is the number of samples to integrate over,
         but each with a vector nv in length
     i is the number of dimensions, or integrals to comptue =#

    for j in 1:ns					# loop for ns samples
        c  = copy(cfill)
        dc = copy(dpfill)
        pv = copy(dpfill)
        for i in 2:n
            x = transpose(abs.(2.0 .* mod.((1:nv) .* q[i - 1] .+ rand(rng), 1) .- 1))	 # periodizing transformation
                # note: the rand() is not broadcast -- it's a single random uniform value added to all elements
            y[i - 1,:] = quantile.(unitnorm, c .+ x .* dc)
            s = transpose(ch[i,1:i - 1]) * y[1:i - 1,:]
            ct = ch[i,i]										# ch is cholesky matrix
            ai = as[i] .- s
            bi = bs[i] .- s
            c = copy(Jnv)										# preset to 1 (>9 sd, +∞)
            d = copy(Jnv)										# preset to 1 (>9 sd, +∞)

            c[findall(x -> isless(x, -9 * ct), ai)] .= 0.0		# If < -9 sd (-∞), set to zero
            d[findall(x -> isless(x, -9 * ct), bi)] .= 0.0		# if < -9 sd (-∞), set to zero
            tstl = findall(x -> isless(abs(x), 9 * ct), ai)		# find cols inbetween -9 and +9 sd (-∞ to +∞)
            c[tstl] .= cdf.(unitnorm, ai[tstl] / ct)			# for those, compute Normal CDF
            tstl = (findall(x -> isless(abs(x), 9 * ct), bi))	# find cols inbetween -9 and +9 sd (-∞ to +∞)
            d[tstl] .= cdf.(unitnorm, bi[tstl] / ct)
            @. dc = d - c
            @. pv = pv * dc
        end # for i=
        d = (mean(pv) - p) / j
        p += d
        e = (j - 2) * e / j + d^2
    end # for j=

    e = 3 * sqrt(e) 	# error estimate is 3 times standard error with ns samples

    return (p, e)  	# return probability value and error estimate
end # function qsimvnv

function _chlrdr(Σ, a, b)

# Rev 1.13

# define constants
# 64 bit machien error 1.0842021724855e-19 ???
# 32 bit machine error 2.220446049250313e-16 ???
    ep = 1e-10 # singularity tolerance
    if Sys.WORD_SIZE == 64
        fpsize = Float64
        ϵ = eps(0.0) # 64-bit machine error
    else
        fpsize = Float32
        ϵ = eps(0.0f0) # 32-bit machine error
    end

    if !@isdefined sqrt2π
        sqrt2π = √(2π)
    end

# unit normal distribution
    unitnorm = Normal()

    n = size(Σ, 1) # covariance matrix n x n square

    ckk = 0.0
    dem = 0.0
    am = 0.0
    bm = 0.0
    ik = 0.0

    if eltype(Σ) <: Signed
        c = copy(float(Σ))
    else
        c = copy(Σ)
    end

    if eltype(a) <: Signed
        ap = copy(float(a))
    else
        ap = copy(a)
    end

    if eltype(b) <: Signed
        bp = copy(float(b))
    else
        bp = copy(b)
    end

    d = sqrt.(diag(c))
    for i in 1:n
        if d[i] > 0.0
            c[:,i] /= d[i]
            c[i,:] /= d[i]
            ap[i] = ap[i] / d[i]     # ap n x 1 vector
            bp[i] = bp[i] / d[i]     # bp n x 1 vector
        end
    end

    y = zeros(fpsize, n) # n x 1 zero vector to start

    for k in 1:n
        ik = k
        ckk = 0.0
        dem = 1.0
        s = 0.0
    # pprinta(c)
        for i in k:n
            if c[i,i] > ϵ  # machine error
                cii = sqrt(max(c[i,i], 0))

                if i > 1 && k > 1
                    s = (c[i,1:(k - 1)] .* y[1:(k - 1)])[1]
                end

                ai = (ap[i] - s) / cii
                bi = (bp[i] - s) / cii
                de = cdf(unitnorm, bi) - cdf(unitnorm, ai)

                if de <= dem
                    ckk = cii
                    dem = de
                    am = ai
                    bm = bi
                    ik = i
                end
            end # if c[i,i]> ϵ
        end # for i=
        i = n

        if ik > k
            ap[ik], ap[k] = ap[k], ap[ik]
            bp[ik], bp[k] = bp[k], bp[ik]

            c[ik,ik] = c[k,k]

            if k > 1
                c[ik,1:(k - 1)], c[k,1:(k - 1)] = c[k,1:(k - 1)], c[ik,1:(k - 1)]
            end

            if ik < n
                c[(ik + 1):n,ik], c[(ik + 1):n,k] = c[(ik + 1):n,k], c[(ik + 1):n,ik]
            end

            if k <= (n - 1) && ik <= n
                c[(k + 1):(ik - 1),k], c[ik,(k + 1):(ik - 1)] = transpose(c[ik,(k + 1):(ik - 1)]), transpose(c[(k + 1):(ik - 1),k])
            end
        end # if ik>k

        if ckk > k * ep
            c[k,k] = ckk
            if k < n
                c[k:k,(k + 1):n] .= 0.0
            end

            for i in (k + 1):n
                c[i,k] /= ckk
                c[i:i,(k + 1):i] -= c[i,k] * transpose(c[(k + 1):i,k])
            end

            if abs(dem) > ep
                y[k] = (exp(-am^2 / 2) - exp(-bm^2 / 2)) / (sqrt2π * dem)
            else
                if am < -10
                    y[k] = bm
                elseif bm > 10
                y[k] = am
            else
                y[k] = (am + bm) / 2
                end
            end # if abs
        else
            c[k:n,k] .== 0.0
            y[k] = 0.0
        end # if ckk>ep*k
    end # for k=

    return (c, ap, bp)
end # function _chlrdr

function testmvn(;m=nothing)

# Typical Usage/Example Code
# Example multivariate Normal CDF for various vectors
    println()

# from MATLAB documentation
    r = [4 3 2 1;3 5 -1 1;2 -1 4 2;1 1 2 5]
 	a = [-Inf; -Inf; -Inf; -Inf] # -inf for each
	b = [1; 2; 3; 4 ]
	m = 5000
	# m=4000 # rule of thumb: 1000*(number of variables)
	(myp, mye) = qsimvnv(r, a, b;m=m)
	println("Answer should about 0.6044 to 0.6062")
	println(myp)
	println("The Error should be <= 0.001 - 0.0014");
	println(mye)

	r = [1 0.6 0.333333;
	   0.6 1 0.733333;
	   0.333333 0.733333 1]
	r  = [  1   3 / 5   1 / 3;
		  3 / 5    1    11 / 15;
		  1 / 3  11 / 15    1]
	a = [-Inf;-Inf;-Inf]
	b = [1;4;2]
	# m=3000;
	(myp, mye) = qsimvnv(r, a, b;m=4000)
	println()
	println("Answer shoudl be about 0.82798")
	# answer from Genz paper 0.827984897456834
	println(myp)
	println("The Error should be <= 2.5-05")
	println(mye)

	r = [1 0.25 0.2; 0.25 1 0.333333333; 0.2 0.333333333 1]
	a = [-1;-4;-2]
	b = [1;4;2]
	# m=3000;
	(myp, mye) = qsimvnv(r, a, b;m=4000);
	println()
	println("Answer should be about 0.6537")
	println(myp)
	println("The Error should be <= 2.5-05")
	println(mye)

	# Genz problem 1.5 p. 4-5  & p. 63
	# correct answer F(a,b) = 0.82798
    r = [1 / 3 3 / 5 1 / 3;
         3 / 5 1.0 11 / 15;
         1 / 3 11 / 15 1.0]
    a = [-Inf; -Inf; -Inf]
    b = [1; 4; 2]
	(myp, mye) = qsimvnv(r, a, b;m=4000)
	println()
	# println("Answer shoudl be about 0.82798")
	println("Answer should be 0.9432")
	println(myp)
	println("The error should be < 6e-05")
	println(mye)

	# Genz page 63 uses a different r Matrix
	r = [1 0 0;
		3 / 5 1 0;
		1 / 3 11 / 15 1]
		a = [-Inf; -Inf; -Inf]
		b = [1; 4; 2]
	(myp, mye) = qsimvnv(r, a, b;m=4000)
	println()
	println("Answer shoudl be about 0.82798")
	println(myp)
	println("The error should be < 6e-05")
	println(mye)
	# mystery solved - he used the wrong sigma Matrix
	# when computing the results on page 6

	# singular cov Example
	r = [1 1 1; 1 1 1; 1 1 1]
	a = [-Inf, -Inf, -Inf]
	b = [1, 1, 1]
	(myp, mye) = qsimvnv(r, a, b)
	println()
	println("Answer should be 0.841344746068543")
	println(myp)
	println("The error should be 0.0")
	println(mye)
	println("Cov matrix is singular")
	println("Problem reduces to a univariate problem with")
	println("p = cdf.(Normal(),1) = ", cdf.(Normal(), 1))

	# 5 dimensional Example
	# c = LinearAlgebra.tri(5)
	r = [1 1 1 1 1;
		 1 2 2 2 2;
		 1 2 3 3 3;
		 1 2 3 4 4;
		 1 2 3 4 5]
	a = [-1,-2,-3,-4,-5]
	b = [2,3,4,5,6]
	(myp, mye) = qsimvnv(r, a, b)
	println()
	println("Answer should be ~ 0.7613")
	# genz gives 0.4741284  p. 5 of book
	# Julia, MATLAB, and R all give ~0.7613 !
	println(myp)
	println("The error should be < 0.001")
	println(mye)

	# genz reversed the integration limits when computing
	# see p. 63
	a = sort!(a)
	b = 1 .- a
	(myp, mye) = qsimvnv(r, a, b)
	println()
	println("Answer should be ~ 0.4741284")
	# genz gives 0.4741284  p. 5 of book
	# now matches with reversed integration limits
	println(myp)
	println("The error should be < 0.001")
	println(mye)

	# positive orthant of above
	a = [0,0,0,0,0]
	(myp, mye) = qsimvnv(r, a, b)
	println()
	println("Answer should be ~  0.11353418")
	# genz gives 0.11353418   p. 6 of book
	println(myp)
	println("The error should be < 0.001")
	println(mye)

	# now with a = -inf
	a = [-Inf,-Inf,-Inf,-Inf,-Inf]
	(myp, mye) = qsimvnv(r, a, b)
	println()
	println("Answer should be ~ 0.81031455")
	# genz gives 0.81031455  p. 6 of book
	println(myp)
	println("The error should be < 0.001")
	println(mye)

	# eight dimensional Example
	r = [1 1 1 1 1 1 1 1;
		 1 2 2 2 2 2 2 2;
		 1 2 3 3 3 3 3 3;
		 1 2 3 4 4 4 4 4;
		 1 2 3 4 5 5 5 5;
		 1 2 3 4 5 6 6 6;
		 1 2 3 4 5 6 7 7;
		 1 2 3 4 5 6 7 8]
	a = -1 * [1,2,3,4,5,6,7,8]
	b = [2,3,4,5,6,7,8,9]
	(myp, mye) = qsimvnv(r, a, b)
	println()
	println("Answer should be ~ 0.7594")
	# genz gives 0.32395   p. 6 of book
	# MATLAB gives 0.7594362
	println(myp)
	println("The error should be < 0.001")
	println(mye)

	# eight dim orthant
	a = [0,0,0,0,0,0,0,0]
	b = [Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf]
	(myp, mye) = qsimvnv(r, a, b)
	println()
	println("Answer should be ~ 0.196396")
	# genz gives 0.076586  p. 6 of book
	# MATLB gives 0.196383
	println(myp)
	println("The error should be < 0.001")
	println(mye)

	# eight dim with a = -inf
	a = -Inf * [1,1,1,1,1,1,1,1]
	b = [2,3,4,5,6,7,8,9]
	# b = [0,0,0,0,0,0,0,0]
	(myp, mye) = qsimvnv(r, a, b)
	println()
	println("Answer should be ~ 0.9587")
	# genz gives 0.69675    p. 6 of book
	# MATLAB gives 0.9587
	println(myp)
	println("The error should be < 0.001")
	println(mye)
end # testmvn


text:
@testset "Pressure camera" begin

@testset "Parallel ray generator" begin
    ray = GRAYTR.ray_parallel(2.0, 0.0, 0.0)
    @test ray.origin ≈ Point3(0.0, 0.0, 1.0)
    @test ray.direction ≈ Vector3(0.0, 0.0, -1.0)
    @test norm(ray.direction) ≈ 1.0

    # Need to explicitly loosen test tolerance here because of float instability
    # with trinometry. There are some vector components that should be zero but get
    # values on the order of 1e-16.

    ray = GRAYTR.ray_parallel(2.0, 1.0, 0.0)
    @test ray.direction ≈ Vector3(0.0, 0.0, -1.0)
    @test isapprox(ray.origin, Point3(2.0, 0.0, 1.0) ; atol=1e-15)

    ray = GRAYTR.ray_parallel(2.0, 0.0, 1.0)
    @test ray.direction ≈ Vector3(0.0, 0.0, -1.0)
    @test ray.origin ≈ Point3(0.0, 0.0, 1.0)

    ray = GRAYTR.ray_parallel(2.0, 1.0, 1.0)
    @test ray.direction ≈ Vector3(0.0, 0.0, -1.0)
    @test isapprox(ray.origin, Point3(2.0, 0.0, 1.0) ; atol=1e-15)

    ray = GRAYTR.ray_parallel(2.0, 1.0, 0.5)
    @test isapprox(ray.origin, Point3(-2.0, 0.0, 1.0) ; atol=1e-15)

end

@testset "Point-source ray generator" begin
    ray = GRAYTR.ray_from_point(2.0, 4.0, 1.0, 0.0)
    @test ray.origin ≈ Point3(0.0, 0.0, 4.0)
    @test ray.direction ≈ Vector3(0.0, 0.0, -1.0)
    @test norm(ray.direction) ≈ 1.0

    # Need to explicitly loosen test tolerance here because of float instability
    # with trinometry. There are some vector components that should be zero but get
    # values on the order of 1e-16.

    ray = GRAYTR.ray_from_point(2.0, 4.0, 1.0, 1.0)
    @test ray.origin ≈ Point3(0.0, 0.0, 4.0)
    @test ray.direction ≈ Vector3(0.0, 0.0, -1.0)
    @test norm(ray.direction) ≈ 1.0

    ct = sqrt(1.0 - (2.0 / 4.0)^2)

    ray = GRAYTR.ray_from_point(2.0, 4.0, 0.0, 0.0)
    @test ray.origin ≈ Point3(0.0, 0.0, 4.0)
    @test isapprox(ray.direction, Vector3(sqrt(1-ct^2), 0.0, -ct), atol=1e-15)
    @test norm(ray.direction) ≈ 1.0

    ray = GRAYTR.ray_from_point(2.0, 4.0, 0.0, 0.5)
    @test ray.origin ≈ Point3(0.0, 0.0, 4.0)
    @test isapprox(ray.direction, Vector3(-sqrt(1-ct^2), 0.0, -ct), atol=1e-15)
    @test norm(ray.direction) ≈ 1.0

    ray = GRAYTR.ray_from_point(2.0, 4.0, 0.0, 0.25)
    @test ray.origin ≈ Point3(0.0, 0.0, 4.0)
    @test isapprox(ray.direction, Vector3(0.0, sqrt(1-ct^2), -ct), atol=1e-15)
    @test norm(ray.direction) ≈ 1.0

    ray = GRAYTR.ray_from_point(2.0, 4.0, 0.0, 0.75)
    @test ray.origin ≈ Point3(0.0, 0.0, 4.0)
    @test isapprox(ray.direction, Vector3(0.0, -sqrt(1-ct^2), -ct), atol=1e-15)
    @test norm(ray.direction) ≈ 1.0

end

end

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
# For testing the inner merging functionality of `addentry!`
struct TestMerger
    val::Int
end
Counters.datamerge(a::TestMerger, b::TestMerger) = TestMerger(a.val + b.val)

@testset "Testing Database" begin
    # Test ntsorter
    nt = (c = 1, b = 2, a = 3)
    @test Counters.ntsort(nt) == (a = 3, b = 2, c = 1)

    # ismatch
    a = (a = 1, b = 2)
    @test Counters.ismatch(a, a)

    # Make sure "data" is ignored
    a = (a = 1, b = 2, data = 10)
    b = (a = 1, b = 2, data = 5)
    @test Counters.ismatch(a, b)

    # Test missing
    a = (a = nothing, b = 2)
    b = (a = nothing, b = 2)
    @test Counters.ismatch(a, b)

    b = (a = 1, b = 2)
    @test Counters.ismatch(a, b) == false

    # Test subset
    b = (a = nothing,)
    @test Counters.ismatch(a, b)

    # Test expansion.
    @test Counters.expand((a = 1, c = 2), (:b, :c)) == (a = 1, b = nothing, c = 2)

    #####
    ##### The database
    #####
    # Now try some schenanigans with a DataBase.
    db = Counters.database()
    @test length(db) == 0

    # Make sure showing an empty database doesn't error out.
    println(db)

    # Start adding entries to make sure that everything works more or less as we want.
    db = Counters.addentry!(db, (a = 1, b = 2), Dict(:data => TestMerger(10)))
    @test length(db) == 1
    println(db)

    # Add data that matches
    db = Counters.addentry!(db, (a = 1, b = 2), Dict(:data => TestMerger(100), :a => 1))
    # Ensure inner `datamerge` was called correctly
    @test db[1].data[:data] == TestMerger(110)
    @test db[1].data[:a] == 1

    # Make sure information propogation works.
    db = Counters.addentry!(db, (b = 3, c = 4))
    @test length(db) == 2
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
abstract type AbstractCoefficients{T} end

Base.isequal(coeff1::AbstractCoefficients{T1}, coeff2::AbstractCoefficients{T2}) where {T1,T2} = (coeff1 == coeff2 && T1 == T2 && typeof(coeff1) == typeof(coeff2))


text:
"""
Package: Forecast

Store results from the function `summarize`

# Arguments
`quantiles::DataFrame`        DataFrame with the data quantiles for each column
`moment::DataFrame`           DataFrame with the first four moments for each column
`format::DataFrame`           DataFrame with number of types per column
"""
mutable struct SUMMARIZE
    quantiles::DataFrame
    moments::DataFrame
    format::DataFrame
end

function Base.show(io::IO, s::SUMMARIZE)
    pretty_table(s.quantiles, nosubheader = true, show_row_number=false, crop = :none)
    pretty_table(s.moments, nosubheader = true, show_row_number=false, crop = :none)
    pretty_table(s.format, nosubheader = true, show_row_number=false, crop = :none)
end


text:
export
    occult,
    ordc,
    ordd,
    ordi,
    orderc,
    orderd,
    orderi,
    oscelt,
    oscltx

"""
    occult(targ1, shape1, frame1, targ2, shape2, frame2, abcorr, obsrvr, et)

Determines the occultation condition (not occulted, partially, etc.) of one target relative to
another target as seen by an observer at a given time.

The surfaces of the target bodies may be represented by triaxial ellipsoids or by topographic data
provided by DSK files.

### Arguments ###

- `targ1`: Name or ID of first target.
- `shape1`: Type of shape model used for first target.
- `frame1`: Body-fixed, body-centered frame for first body.
- `targ2`: Name or ID of second target.
- `shape2`: Type of shape model used for second target.
- `frame2`: Body-fixed, body-centered frame for second body.
- `abcorr`: Aberration correction flag.
- `obsrvr`: Name or ID of the observer.
- `et`: Time of the observation (seconds past J2000).

### Output ###

Returns the occultation identification code.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/occult_c.html)
"""
function occult(targ1, shape1, frame1, targ2, shape2, frame2, abcorr, obsrvr, et)
    frame1 = isempty(frame1) ? " " : frame1
    frame2 = isempty(frame2) ? " " : frame2
    ocltid = Ref{SpiceInt}()
    ccall((:occult_c, libcspice), Cvoid,
          (Cstring, Cstring, Cstring, Cstring, Cstring, Cstring, Cstring, Cstring, SpiceDouble,
           Ref{SpiceInt}),
          targ1, shape1, frame1, targ2, shape2, frame2, abcorr, obsrvr, et, ocltid)
    handleerror()
    Int(ocltid[])
end

"""
    ordc(set, item)

The function returns the ordinal position of any given item in a character set.

### Arguments ###

- `set`: A set to search for a given item
- `item`: An item to locate within a set

### Output ###

Returns the ordinal position or `nothing` if the items does not appear in the set.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ordc_c.html)
"""
function ordc(set, item)
    res = ccall((:ordc_c, libcspice), SpiceInt,
                (Cstring, Ref{Cell{SpiceChar}}),
                item, set.cell)
    handleerror()
    res == -1 ? nothing : res + 1
end

"""
    ordd(set, item)

The function returns the ordinal position of any given item in a character set.

### Arguments ###

- `set`: A set to search for a given item
- `item`: An item to locate within a set

### Output ###

Returns the ordinal position or `nothing` if the items does not appear in the set.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ordd_c.html)
"""
function ordd(set, item)
    res = ccall((:ordd_c, libcspice), SpiceInt,
                (SpiceDouble, Ref{Cell{SpiceDouble}}),
                item, set.cell)
    handleerror()
    res == -1 ? nothing : res + 1
end

"""
    ordi(set, item)

The function returns the ordinal position of any given item in a character set.

### Arguments ###

- `set`: A set to search for a given item
- `item`: An item to locate within a set

### Output ###

Returns the ordinal position or `nothing` if the items does not appear in the set.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ordi_c.html)
"""
function ordi(set, item)
    res = ccall((:ordi_c, libcspice), SpiceInt,
                (SpiceInt, Ref{Cell{SpiceInt}}),
                item, set.cell)
    handleerror()
    res == -1 ? nothing : res + 1
end

function _orderc(array)
    array_, ndim, lenvals = chararray(array)
    iorder = Array{SpiceInt}(undef, ndim)
    ccall((:orderc_c, libcspice), Cvoid,
          (SpiceInt, Ref{SpiceChar}, SpiceInt, Ref{SpiceInt}),
          lenvals, array_, ndim, iorder)
    handleerror()
    Int.(iorder) .+ 1
end

@deprecate orderc sortperm

"""
    orderc(array)

!!! warning "Deprecated"
    Use `sortperm` instead.
"""
orderc

function _orderd(array)
    ndim = length(array)
    iorder = Array{SpiceInt}(undef, ndim)
    ccall((:orderd_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, SpiceInt, Ref{SpiceInt}),
          array, ndim, iorder)
    handleerror()
    Int.(iorder) .+ 1
end

@deprecate orderd sortperm

"""
    orderd(array)

!!! warning "Deprecated"
    Use `sortperm` instead.
"""
orderd

function _orderi(array)
    array_ = SpiceInt.(array)
    ndim = length(array)
    iorder = Array{SpiceInt}(undef, ndim)
    ccall((:orderi_c, libcspice), Cvoid,
          (Ref{SpiceInt}, SpiceInt, Ref{SpiceInt}),
          array_, ndim, iorder)
    handleerror()
    Int.(iorder) .+ 1
end

@deprecate orderi sortperm

"""
    orderi(array)

!!! warning "Deprecated"
    Use `sortperm` instead.
"""
orderi

"""
    oscelt(state, et, mu)

Determine the set of osculating conic orbital elements that corresponds to
the state (position, velocity) of a body at some epoch.

### Arguments ###

- `state`: State of body at epoch of elements
- `et`: Epoch of elements
- `mu`: Gravitational parameter (GM) of primary body

### Output ###

Returns the equivalent conic elements:

- `rp`: Perifocal distance
- `ecc`: Eccentricity
- `inc`: Inclination
- `lnode`: Longitude of the ascending node
- `argp`: Argument of periapsis
- `m0`: Mean anomaly at epoch
- `t0`: Epoch
- `mu`: Gravitational parameter

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/oscelt_c.html)
"""
function oscelt(state, et, mu)
    @checkdims 6 state
    elts = Array{SpiceDouble}(undef, 8)
    ccall((:oscelt_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, SpiceDouble, SpiceDouble, Ref{SpiceDouble}),
          state, et, mu, elts)
    handleerror()
    elts
end

"""
    oscltx(state, et, mu)

Determine the set of osculating conic orbital elements that corresponds to the state
(position, velocity) of a body at some epoch. In addition to the classical elements,
return the true anomaly, semi-major axis, and period, if applicable.

### Arguments ###

- `state`: State of body at epoch of elements
- `et`: Epoch of elements
- `mu`: Gravitational parameter (GM) of primary body

### Output ###

Returns the extended set of classical conic elements:

- `rp`: Perifocal distance.
- `ecc`: Eccentricity.
- `inc`: Inclination.
- `lnode`: Longitude of the ascending node.
- `argp`: Argument of periapsis.
- `m0`: Mean anomaly at epoch.
- `t0`: Epoch.
- `mu`: Gravitational parameter.
- `nu`: True anomaly at epoch.
- `a`: Semi-major axis. A is set to zero if it is not computable.
- `tau`: Orbital period. Applicable only for elliptical orbits. Set to zero otherwise.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/oscltx_c.html)
"""
function oscltx(state, et, mu)
    @checkdims 6 state
    elts = Array{SpiceDouble}(undef, 11)
    ccall((:oscltx_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, SpiceDouble, SpiceDouble, Ref{SpiceDouble}),
          state, et, mu, elts)
    handleerror()
    elts
end


text:
## writer

type Writer
    _cache::Dict
    _partialCache::Dict
    _loadPartial ## Function or nothing
end

Writer() = Writer(Dict(), Dict(), nothing)

function clearCache(w::Writer)
    w._cache=Dict()
    w._partialCache=Dict()
end

function compile(io::IO, w::Writer, template, tags)
#    if haskey(w._cache, template)
#        return(w._cache[template])
#    end

##    tokens = parse(template, tags)
    tokens = template
    compileTokens(io, w, tokens.tokens, template)
#    w._cache[template] = compileTokens(io, w, tokens.tokens, template)

#    return(w._cache[template])
end

function compilePartial(w::Writer, name, template, tags)
    fn = compile(w, template, tags)
    w._partialCache[name] = fn
    fn
end

function getPartial(w::Writer, name)
## didn't do loadPartial, as not sure where template is
#    if !haskey(w._partialCache, name) && is(w._loadPartial, Function)
#        compilePartial(w,

    w._partialCache[name]
end

function compileTokens(io, w::Writer, tokens, template)
    ## return a function
    function f(w::Writer, view) #  no partials
       renderTokens(io, tokens, w, Context(view), template) # io in closure
    end
    return(f)
end

function render(io::IO, w::Writer, template, view)
    f = compile(io, w, template, ["{{", "}}"])
    f(w, view)
end


text:
include(joinpath(@__DIR__,"..","src","open_file.jl"))

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "--diff"
        help = "print to STDOUT the output of `diff`"
        action = :store_true
    "--two"
        help = "open and display both files"
        action = :store_true
    "--bw"
        help = "generate, save and display a B&W difference image for PNG and SVG files.  requires Rsvg, Cairo, and Images"
        action = :store_true
    "filter"
        help = "a regular expression describing the filenames to compare"
        default=""
end

args = parse_args(s)

if args["bw"]
    using Rsvg, Cairo, Images

    function svg2img(filename)
        r = Rsvg.handle_new_from_file(filename)
        d = Rsvg.handle_get_dimensions(r)
        cs = Cairo.CairoImageSurface(d.width,d.height,Cairo.FORMAT_ARGB32)
        c = Cairo.CairoContext(cs)
        Rsvg.handle_render_cairo(c,r)
        fout = tempname()
        Cairo.write_to_png(cs,fout)
        png = load(fout)
        rm(fout)
        png
    end
end

# delete diffedoutput/
repo = LibGit2.GitRepo(dirname(@__DIR__))
options = LibGit2.StatusOptions(flags=LibGit2.Consts.STATUS_OPT_INCLUDE_IGNORED |
                                      LibGit2.Consts.STATUS_OPT_RECURSE_IGNORED_DIRS)
status = LibGit2.GitStatus(repo, status_opts=options)
for i in 1:length(status)
    entry = status[i]
    index_to_workdir = unsafe_load(entry.index_to_workdir)
    if index_to_workdir.status == Int(LibGit2.Consts.DELTA_IGNORED)
        filepath = unsafe_string(index_to_workdir.new_file.path)
        startswith(filepath,joinpath("test/diffedoutput")) || continue
        rm(filepath)
    end
end

# Compare with cached output
cachedout = joinpath((@__DIR__), "cachedoutput")
gennedout = joinpath((@__DIR__), "gennedoutput")
diffedout = joinpath((@__DIR__), "diffedoutput")
ndifferentfiles = 0
const creator_producer = r"(Creator|Producer)"
filter_mkdir_git(x) = !mapreduce(y->x==y,|,[".mkdir","git.log","git.status"])
filter_regex(x) = ismatch(Regex(args["filter"]), x)
cached_files = filter(x->filter_mkdir_git(x) && filter_regex(x), readdir(cachedout))
genned_files = filter(x->filter_mkdir_git(x) && filter_regex(x), readdir(gennedout))
cached_notin_genned = setdiff(cached_files, genned_files)
isempty(cached_notin_genned) ||
      warn("files in cachedoutput/ but not in gennedoutput/: ", join(cached_notin_genned,", "))
genned_notin_cached = setdiff(genned_files, cached_files)
isempty(genned_notin_cached) ||
      warn("files in gennedoutput/ but not in cachedoutput/: ", join(genned_notin_cached,", "))
for file in intersect(cached_files,genned_files)
    print("Comparing ", file, " ... ")
    cached = open(readlines, joinpath(cachedout, file))
    genned = open(readlines, joinpath(gennedout, file))
    same = length(cached) == length(genned)
    if same
        lsame = Bool[cached[i] == genned[i] for i = 1:length(cached)]
        if !all(lsame)
            for idx in find(lsame.==false)
                # Don't worry about lines that are due to
                # Creator/Producer (e.g., Cairo versions)
                if !isempty(search(cached[idx], creator_producer))
                    lsame[idx] = true
                end
            end
        end
        same = same & all(lsame)
    end
    if same
        println("same!")
    else
        ndifferentfiles +=1
        println("different :(")
        if args["diff"]
            diffcmd = `diff $(joinpath(cachedout, file)) $(joinpath(gennedout, file))`
            run(ignorestatus(diffcmd))
        end
        if args["two"]
            open_file("$(joinpath(cachedout,file))")
            open_file("$(joinpath(gennedout,file))")
        end
        if args["bw"] && (endswith(file,".svg") || endswith(file,".png"))
            wait_for_user = false
            if endswith(file,".svg")
                gimg = svg2img(joinpath(gennedout,file));
                cimg = svg2img(joinpath(cachedout,file));
            elseif endswith(file,".png")
                gimg = load(joinpath(gennedout,file));
                cimg = load(joinpath(cachedout,file));
            end
            if size(gimg)==size(cimg)
                dimg = convert(Matrix{Gray}, gimg.==cimg)
                if any(dimg.==0)
                    fout = joinpath(diffedout,file*".png")
                    Images.save(fout, dimg)
                    wait_for_user = true
                    open_file("$fout")
                else
                    println("files are different but PNGs are the same")
                end
            else
                wait_for_user = true
                println("PNGs are different sizes :(")
            end
        end
        args["diff"] || args["two"] || (args["bw"] && wait_for_user) || continue
        println("Press ENTER to continue, CTRL-C to quit")
        readline()
    end
end

infoorwarn = ndifferentfiles==0 ? info : warn
infoorwarn("# different images = ",ndifferentfiles)


text:
compile(path::String, bf) = open(io -> write(io, compile(bf)), path, "w")

bfrun(t::Tape, x) = interpret(t, compile(x))
bfrun(x) = bfrun(Tape(), x)

macro run(ex)
  :(bfrun(@bf($(esc(ex)))))
end

function stack(t::Tape)
  stk = similar(t.tape, 0)
  interpret(t, "[<<]")
  while get(t.tape, t.pos + 2, 0) ≠ 0
    push!(stk, t.tape[t.pos + 1])
    t.pos += 2
  end
  return stk
end

for T in [Native, Flip, Word, Quote]
  @eval a::$T == b::$T = a.code == b.code
end

max_stars_repo_name: JuliaTagBot/BrainForth

max_issues_repo_name: JuliaTagBot/BrainForth

max_forks_repo_name: JuliaTagBot/BrainForth


text:
import GaussianMixtureTest
using Distributions
using KernelEstimator
using RCall

mu_true = [log(1/0.779 - 1)]
wi_true = [1.0]
sigmas_true = [1.2]
n = 282
m = MixtureModel(map((u, v) -> Normal(u, v), mu_true, sigmas_true), wi_true)

T1 = zeros(100)
P = zeros(100)
for b in 1:100
    srand(b)
    x = rand(m, n)
    T1[b], P[b] = GaussianMixtureTest.kstest(x, 1)
    print(b,"->", T1[b], " | ")
end

xs = linspace(0.01, 14, 500)
den2 = kerneldensity(T1, xeval=xs, kernel=gammakernel, lb=0.)

@rput xs den2 T1 P;
rprint("""
hist(T1, breaks=15, freq=F, ylim=c(0, .4))
rug(T1)
lines(xs, dchisq(xs, 2))
lines(xs, den2, lwd=2, col="blue")
""")

import GaussianMixtureTest
using Distributions
using KernelEstimator
using RCall

mu_true = [-2.0858,-1.4879]
wi_true = [0.0828,0.9172]
sigmas_true = [0.6735,0.2931]
n = 282
m = MixtureModel(map((u, v) -> Normal(u, v), mu_true, sigmas_true), wi_true)
srand(35);x = rand(m, n);GaussianMixtureTest.kstest(x, 2)

T1 = zeros(100)
P = zeros(100)
for b in 1:100
    srand(b)
    x = rand(m, n)
    T1[b], P[b] = GaussianMixtureTest.kstest(x, 2)
    print(b,"->", T1[b], " | ")
end

x = rand(m, n)
wi,mu,sigmas= GaussianMixtureTest.gmm(x, 2)
Ttrue = GaussianMixtureTest.asymptoticdistribution(x, wi, mu, sigmas)

xs = linspace(0.01, 14, 500)
den1=kerneldensity(Ttrue, xeval=xs, kernel=gammakernel, lb=0.)
den2 = kerneldensity(T1, xeval=xs, kernel=gammakernel, lb=0.)

@rput xs den1 den2 T1 P;
rprint("""
hist(T1, breaks=15, freq=F, ylim=c(0, .4))
rug(T1)
lines(xs, den1, lwd=2)
lines(xs, den2, lwd=2, col="blue")
""")

rprint("""
z = qnorm(1-P)
hist(z, freq=F, xlim=c(-4,4))
lines(density(z))
curve(dnorm, -4, 4, col="red",add=T)
NULL
""")

### The components

import GaussianMixtureTest
using Distributions
using KernelEstimator
using RCall

mu_true = [log(1/0.779 - 1)/3 - 4.0, log(1/0.779 - 1)/3 + 1.0, log(1/0.779 - 1)/3 + 4.0;]
wi_true = [.3, .4, .3]
sigmas_true = [1.2, .8, .9]
n = 282
m = MixtureModel(map((u, v) -> Normal(u, v), mu_true, sigmas_true), wi_true)
#srand(35);x = rand(m, n);GaussianMixtureTest.kstest(x, 2)

T1 = zeros(100)
P = zeros(100)
for b in 1:100
    srand(b)
    x = rand(m, n)
    T1[b], P[b] = GaussianMixtureTest.kstest(x, 3)
    print(b,"->", T1[b], " | ")
end

x = rand(m, n)
wi,mu,sigmas= GaussianMixtureTest.gmm(x, 3)
Ttrue = GaussianMixtureTest.asymptoticdistribution(x, wi, mu, sigmas);

xs = linspace(0.01, 14, 500);
den1=kerneldensity(Ttrue, xeval=xs, kernel=gammakernel, lb=0.);
den2 = kerneldensity(T1, xeval=xs, kernel=gammakernel, lb=0.);

@rput xs den1 den2 T1 P;
rprint("""
hist(T1,freq=F, xlim=c(0, 15), ylim=c(0, .22))
rug(T1)
lines(xs, den1, lwd=2)
lines(xs, den2, lwd=2, col="blue")
""")

## The fitting

import GaussianMixtureTest
using Distributions
using KernelEstimator
using RCall

mu_true = [log(1/0.779 - 1)/3 - 4.0, log(1/0.779 - 1)/3 + 1.0, log(1/0.779 - 1)/3 + 4.0;]
wi_true = [.3, .4, .3]
sigmas_true = [1.2, .8, .9]
m = MixtureModel(map((u, v) -> Normal(u, v), mu_true, sigmas_true), wi_true)

x = rand(m, 500)
wi, mu, sigmas = GaussianMixtureTest.gmm(x, 3)
mhat = MixtureModel(map((u, v) -> Normal(u, v), mu, sigmas), wi)

xs = linspace(-6, 5, 500)
dentrue = pdf(m, xs)
denhat = pdf(mhat, xs)

@rput xs dentrue denhat
rprint("""
plot(xs, dentrue, lwd=2, type="l")
lines(xs, denhat, lwd=2, col="blue")
""")

Ttrue = GaussianMixtureTest.asymptoticdistribution(x, wi_true, mu_true, sigmas_true)
xs = linspace(0.01, 14, 500)
den1=kerneldensity(Ttrue, xeval=xs, kernel=gammakernel, lb=0.)

@rput xs den1
rprint("""
plot(xs, den1, type="l")
""")


text:
# included in runtests.jl
# Very similar to client_serverWS_test.jl
# Test sending / receiving messages correctly,
# closing from within websocket handlers,
# symmetry of client and server side websockets,
# stress tests opening and closing a sequence of servers.
# At this time, we unfortunately get irritating messages
# 'Workqueue inconsistency detected:...'

# @info "External server http request"
# @test 200 == HTTP.request("GET", EXTERNALHTTP).status

@info "Listen: Open, http response, close. Repeat three times. Takes a while."
for i = 1:3
    let
        server = startserver(url=SURL,port=PORT)
        status = HTTP.request("GET", "http://$SURL:$PORT").status
        println("Status($(i)): $(status)")
        @test 200 == status
        close(server)
    end
end

@info "Listen: Client side initates message exchange."
let
    server = startserver(url=SURL,port=PORT)
    WebSockets.open(initiatingws, "ws://$SURL:$PORT")
    close(server)
end

@info "Listen: Server side initates message exchange."
let
    server = startserver(url=SURL,port=PORT)
    WebSockets.open(echows, "ws://$SURL:$PORT", subprotocol = SUBPROTOCOL)
    close(server)
end

@info "Listen: Server side initates message exchange. Close from within server side handler."
let
    server = startserver(url=SURL,port=PORT)
    WebSockets.open(echows, "ws://$SURL:$PORT", subprotocol = SUBPROTOCOL_CLOSE)
    close(server)
end
nothing


text:
# modified from https://github.com/mohamed82008/LinearElasticity.jl
using Einsum: @einsum
using LinearAlgebra: I

function get_Kσs(problem::TrussProblem{xdim, TT}, u_dofs, cellvalues) where {xdim, TT}
    Es = getE(problem)
    νs = getν(problem)
    As = getA(problem)
    dh = problem.ch.dh

    # usually ndof_pc = xdim * n_basefuncs
    # ? number of nodes per cell == n_basefuncs per cell
    ndof_pc = ndofs_per_cell(dh)
    n_basefuncs = getnbasefunctions(cellvalues)
    global_dofs = zeros(Int, ndof_pc)
    Kσs = [zeros(TT, ndof_pc, ndof_pc) for i in 1:getncells(dh.grid)]
    Kσ_e = zeros(TT, ndof_pc, ndof_pc)
    # block-diagonal - block σ_e = σ_ij, i,j in xdim
    ψ_e = zeros(TT, xdim*xdim, xdim*xdim)
    G = zeros(TT, xdim*xdim, xdim*n_basefuncs)
    δ = Matrix(TT(1.0)I, xdim, xdim)
    ϵ = zeros(TT, xdim, xdim)
    σ = zeros(TT, xdim, xdim)
    # u_i,j: partial derivative
    u_p = zeros(TT, xdim, xdim)

    for (cellidx, cell) in enumerate(CellIterator(dh))
        Kσ_e .= 0
        truss_reinit!(cellvalues, cell, As[cellidx])
        # get cell's dof's global dof indices, i.e. CC_a^e
        celldofs!(global_dofs, dh, cellidx)
        E = Es[cellidx]
        ν = νs[cellidx]
        for q_point in 1:getnquadpoints(cellvalues)
            dΩ = getdetJdV(cellvalues, q_point)
            for d in 1:xdim
                ψ_e[(d-1)*xdim+1:d*xdim, (d-1)*xdim+1:d*xdim] .= 0
            end
            for a in 1:n_basefuncs
                ∇ϕ = shape_gradient(cellvalues, q_point, a)
                # given displacement values of the cell nodes
                _u = @view u_dofs[(@view global_dofs[xdim*(a-1) .+ (1:xdim)])]
                # u_i,j, i for spatial xdim, j for partial derivative
                @einsum u_p[i,j] = _u[i]*∇ϕ[j]
                # effect of the quadratic term in the strain formula have on the stress field is ignored
                @einsum ϵ[i,j] = 1/2*(u_p[i,j] + u_p[j,i])
                # isotropic solid
                # @einsum σ[i,j] = E*ν/(1-ν^2)*δ[i,j]*ϵ[k,k] + E*ν*(1+ν)*ϵ[i,j]
                # ! truss element special treatment here
                σ = E .* ϵ

                for d in 1:xdim
                    ψ_e[(d-1)*xdim+1:d*xdim, (d-1)*xdim+1:d*xdim] .+= σ
                    G[(xdim*(d-1)+1):(xdim*d), (a-1)*xdim+d] .= ∇ϕ
                end
            end
            # @show G
            Kσ_e .+= G'*ψ_e*G*dΩ
        end
        # @show Kσ_e
        Kσs[cellidx] .= Kσ_e
    end
    return Kσs
end

function buckling(problem::TrussProblem{xdim, T}, ginfo, einfo; u=undef) where {xdim, T}
    dh = problem.ch.dh

    if u === undef
        u = ginfo.K \ ginfo.f
    end
    Kσs = get_Kσs(problem, u, einfo.cellvalues)
    Kσ = deepcopy(ginfo.K)

    # @show Kσs

    if Kσ isa Symmetric
        Kσ.data.nzval .= 0
        assembler = Ferrite.AssemblerSparsityPattern(Kσ.data, T[], Int[], Int[])
    else
        Kσ.nzval .= 0
        assembler = Ferrite.AssemblerSparsityPattern(Kσ, T[], Int[], Int[])
    end

    # * assemble global geometric stiffness matrix
    global_dofs = zeros(Int, ndofs_per_cell(dh))
    Kσ_e = zeros(T, size(Kσs[1]))
    celliteratortype = CellIterator{typeof(dh).parameters...}
    _celliterator::celliteratortype = CellIterator(dh)
    TK = eltype(Kσs)
    for (i,cell) in enumerate(_celliterator)
        celldofs!(global_dofs, dh, i)
        if TK <: Symmetric
            Ferrite.assemble!(assembler, global_dofs, Kσs[i].data)
        else
            Ferrite.assemble!(assembler, global_dofs, Kσs[i])
        end
    end

    # f = copy(ginfo.f)
    # apply!(Kσ, f, problem.ch)

    return ginfo.K, Kσ
end


text:
module TestKernels
using GaussianProcesses, Calculus
using Test, LinearAlgebra, Statistics, Random

Random.seed!(1)

@testset "Kernels" begin
    d, n, n2 = 3, 10, 5
    X = randn(d, n)
    X2 = randn(d, n2)
    y = randn(n)
    ll = rand(d)

    # Random columns
    i, j = rand(1:n), rand(1:n2)
    Xi = view(X, :, i)
    Xj = view(X, :, j) # works since n2 < n
    X2j = view(X2, :, j)

    # Preallocate some matrices
    cK = zeros(n, n)
    cK2 = zeros(n, n2)

    kernels = [# Isotropic kernels
               SEIso(1.0, 1.0), Mat12Iso(1.0,1.0), Mat32Iso(1.0,1.0), Mat52Iso(1.0,1.0),
               RQIso(1.0, 1.0, 1.0), Periodic(1.0, 1.0, 2π),
               # Non-isotropic
               Lin(1.0), Poly(0.0, 0.0, 2), Noise(1.0),
               # Constant kernel
               Const(1.0),
               # ARD kernels
               SEArd(ll, 1.0), Mat12Ard(ll, 1.0), Mat32Ard(ll, 1.0), Mat52Ard(ll, 1.0),
               RQArd(ll, 0.0, 2.0), LinArd(ll),
               # Composite kernels
               SEIso(1.0, 1.0) + Mat12Iso(1.0, 1.0),
               (SEIso(1.0, 1.0) + Mat12Iso(1.0, 1.0)) + Lin(1.0),
               SEIso(1.0, 1.0) * Mat12Iso(1.0, 1.0),
               (SEIso(1.0, 1.0) * Mat12Iso(1.0, 1.0)) * Lin(1.0),
               # Fixed Kernel
               fix(RQIso(1.0, 1.0, 1.0), :lσ), fix(RQIso(1.0, 1.0, 1.0)),
               # Sum and Product and Fix
               SEIso(1.0, 1.0) * Mat12Iso(1.0, 1.0) +
               Lin(1.0) * fix(RQIso(1.0, 1.0, 1.0), :lσ)]

    @testset "Kernel $(typeof(kernel))" for kernel in kernels
        println("\tTesting ", nameof(typeof(kernel)), "...")
        @test length(GaussianProcesses.get_param_names(kernel)) ==
            length(GaussianProcesses.get_params(kernel)) ==
            GaussianProcesses.num_params(kernel)

        @testset for masked in (true, false)
            if masked
                # This is a bit of a hack, to remove one parameter
                # from an ARD kernel to make it have the right number
                # of parameters when masked
                if isa(kernel, LinArd) || isa(kernel, GaussianProcesses.StationaryARD)
                    par = GaussianProcesses.get_params(kernel)
                    k_masked = typeof(kernel)([par[1]], par[d+1:end]...)
                    kern = Masked(k_masked, [1])
                else
                    kern = Masked(kernel, [1])
                end
            else
                kern = kernel
            end

            @testset "Variance" begin
                spec = cov(kern, X)
                @test spec ≈ invoke(cov, Tuple{Kernel, Matrix{Float64}}, kern, X)
                @test spec[i, j] ≈ cov(kern, Xi, Xj)

                key = GaussianProcesses.kernel_data_key(kern, X)
                @test typeof(key) == String
                # check we've overwritten the default if necessary
                kdata = GaussianProcesses.KernelData(kern, X)
                if typeof(kdata) != GaussianProcesses.EmptyData
                    @test key != "EmptyData"
                end
            end

            @testset "Covariance" begin
                spec = cov(kern, X, X2)
                @test spec[i,j] ≈ cov(kern, Xi, X2j)
            end

            @testset "Gradient" begin
                nparams = GaussianProcesses.num_params(kern)
                init_params = GaussianProcesses.get_params(kern)
                data = GaussianProcesses.KernelData(kern, X)
                stack1 = Array{Float64}(undef, n, n, nparams)
                stack2 = Array{Float64}(undef, n, n, nparams)

                GaussianProcesses.grad_stack!(stack1, kern, X, data)
                invoke(GaussianProcesses.grad_stack!,
                       Tuple{AbstractArray, Kernel, Matrix{Float64},
                             GaussianProcesses.EmptyData},
                       stack2, kern, X, GaussianProcesses.EmptyData())
                @test stack1 ≈ stack2

                theor_grad = vec(sum(stack1; dims=[1,2]))
                numer_grad = Calculus.gradient(init_params) do params
                    set_params!(kern, params)
                    sum(cov(kern, X))
                end
                @test theor_grad ≈ numer_grad rtol=1e-1 atol=1e-2
            end

            @testset "dtarget" begin
                gp = GPE(X, y, MeanConst(0.0), kern, -3.0)
                init_params = GaussianProcesses.get_params(gp)
                GaussianProcesses.update_target_and_dtarget!(gp)
                theor_grad = copy(gp.dtarget)
                numer_grad = Calculus.gradient(init_params) do params
                    set_params!(gp, params)
                    update_target!(gp)
                    gp.target
                end
                @test theor_grad ≈ numer_grad rtol=1e-3 atol=1e-3
            end
        end
    end
end
end


text:
using NBody
using ArgParse

function parse_commandline()
  s = ArgParseSettings()

  @add_arg_table s begin
    "--seed", "-s"
      help = "Position only"
      arg_type = Int128
      default = abs(rand(Int128))

    "--number", "-n"
      help = "Number of particles"
      arg_type = Int
      default = 3
      range_tester = x-> x>0
  end

  return parse_args(s)
end

function main(args)
  if args["seed"] > 0
    srand(args["seed"])
  end
  n = args["number"]
  nb = mk_plummer(n)

  if n > 0
    NBody.shift_to_center_of_mass!(nb)
  end
  if n > 1
    NBody.adjust_units!(nb)
  end

  NBody.write_snapshot(nb)

  println(STDERR, "Seed: $(args["seed"])")
end

function frand(low, high)
  low + rand() * (high - low)
end

function sperical(r)
  vec = Array{Float64}(3)

  theta = acos(frand(-1, 1))
  phi = frand(0, 2*pi)

  vec[1] = r * sin(theta) * cos(phi)
  vec[2] = r * sin(theta) * sin(phi)
  vec[3] = r * cos(theta)
  vec
end

function mk_plummer(n)
  res = NBody.NBodySystem(n)

  res.bodies = map(1:n) do i
    b = plummer_sample()
    b.mass = 1 / n
    b.id = i
    b
  end
  res
end

function plummer_sample()
  # scale to viral radius
  scalefactor = 16 / (3 * pi)

  radius = 1 / sqrt( rand() ^ (-2//3) - 1)
  pos = sperical(radius) / scalefactor

  x = 0.0
  y = 0.1
  while y > x*x*(1 - x*x) ^ 3.5
    x = frand(0, 1)
    y = frand(0, 0.1)
  end

  velocity = x * sqrt(2) * (1 + radius*radius) ^ (-1//4)
  vel = sperical(velocity) * sqrt(scalefactor)

  NBody.Body(
    pos = pos,
    vel = vel,
  )
end

main(parse_commandline())


text:
using TaylorModels
using TaylorModels: validated_integ

# parameters of the model
const g = 9.81           # gravity constant in m/s^2
const R = 0.1            # radius of center mass in m
const l = 0.5            # distance of motors to center mass in m
const Mrotor = 0.1       # motor mass in kg
const M = 1.0            # center mass in kg
const m = M + 4*Mrotor   # total mass in kg
const mg = m*g

# moments of inertia
const Jx = (2/5)*M*R^2 + 2*l^2*Mrotor
const Jy = Jx
const Jz = (2/5)*M*R^2 + 4*l^2*Mrotor
const Cyzx = (Jy - Jz)/Jx
const Czxy = (Jz - Jx)/Jy
const Cxyz = 0.0 #(Jx - Jy)/Jz

# considering the control parameters as *parameters*
const u₁ = 1.0
const u₂ = 0.0
const u₃ = 0.0

@inline function quad_property(t, x)
    b1 = (x[3] < 1.4)
    b2 = t ≥ 1.0 ? (x[3] > 0.9) : true
    b3 = t ≥ 5.0 ? (0.98 ≤ x[3] ≤ 1.02) : true
    return b1 && b2 && b3
end

@taylorize function quadrotor!(t, x, dx)
    # unwrap the variables and the controllers; the last three are the controllers
    # x₁, x₂, x₃, x₄, x₅, x₆, x₇, x₈, x₉, x₁₀, x₁₁, x₁₂, u₁, u₂, u₃ = x
    x₁  = x[1]
    x₂  = x[2]
    x₃  = x[3]
    x₄  = x[4]
    x₅  = x[5]
    x₆  = x[6]
    x₇  = x[7]
    x₈  = x[8]
    x₉  = x[9]
    x₁₀ = x[10]
    x₁₁ = x[11]
    x₁₂ = x[12]

    # equations of the controllers
    F = (mg - 10*(x₃ - u₁)) + 3*x₆  # height control
    τϕ = -(x₇ - u₂) - x₁₀            # roll control
    τθ = -(x₈ - u₃) - x₁₁            # pitch control
    local τψ = 0.0                   # heading is uncontrolled
    #
    Tx = τϕ/Jx
    Ty = τθ/Jy
    Tz = τψ/Jz
    F_m = F/m

    # Some abbreviations
    sx7 = sin(x₇)
    cx7 = cos(x₇)
    sx8 = sin(x₈)
    cx8 = cos(x₈)
    sx9 = sin(x₉)
    cx9 = cos(x₉)
    #
    sx7sx9 = sx7*sx9
    sx7cx9 = sx7*cx9
    cx7sx9 = cx7*sx9
    cx7cx9 = cx7*cx9
    sx7cx8 = sx7*cx8
    cx7cx8 = cx7*cx8
    sx7_cx8 = sx7/cx8
    cx7_cx8 = cx7/cx8
    #
    x4cx8 = cx8*x₄
    #
    p11 = sx7_cx8*x₁₁
    p12 = cx7_cx8*x₁₂
    xdot9 = p11 + p12

    # differential equations for the quadrotor
    #
    dx[1] = (cx9*x4cx8 + (sx7cx9*sx8 - cx7sx9)*x₅) + (cx7cx9*sx8 + sx7sx9)*x₆
    dx[2] = (sx9*x4cx8 + (sx7sx9*sx8 + cx7cx9)*x₅) + (cx7sx9*sx8 - sx7cx9)*x₆
    dx[3] = (sx8*x₄ - sx7cx8*x₅) - cx7cx8*x₆
    dx[4] = (x₁₂*x₅ - x₁₁*x₆) - g*sx8
    dx[5] = (x₁₀*x₆ - x₁₂*x₄) + g*sx7cx8
    dx[6] = (x₁₁*x₄ - x₁₀*x₅) + (g*cx7cx8 - F_m)
    dx[7] = x₁₀ + sx8*xdot9
    dx[8] = cx7*x₁₁ - sx7*x₁₂
    dx[9] = xdot9
    dx[10] = Cyzx * (x₁₁ * x₁₂) + Tx
    dx[11] = Czxy * (x₁₀ * x₁₂) + Ty
    dx[12] = Cxyz * (x₁₀ * x₁₁) + Tz
     #
    return dx
end

"""
    quad_TMJets(; t0=0.0, T=5.0, abs_tol=1e-7, orderT=5, orderQ=2, maxsteps=500,
                  property=quad_property)

Build and run the quadrotor model.

### Input

- `t0`       -- (optional, default: `0.0`) initial time
- `T`        -- (optional, default: `5.0`) time horizon
- `abs_tol`  -- (optional, default: `1e-7`) absolute tolerance used for time step
- `orderT`   -- (optional, default: `5`) order of the Taylor model in t
- `orderQ`   -- (optional, default: `2`) order of the Taylor model for Jet transport
                variables
- `maxsteps` -- (optional, default: `1000`) use this maximum number of steps in
                the validated integration
- `property` -- (optional, default: see `quad_property`) safe states property;
                see below

### Notes

The task is to change the height from 0 [m] to 1 [m] within 5 [s]. A goal region [0.98, 1.02] of the
height x3 has to be reached within 5 [s] and the height has to stay below 1.4 for all times. After
1 [s] the height should stay above 0.9 [m]. The initial position of the quadrotor is uncertain in
all directions within [−0.4, 0.4] [m] and also the velocity is uncertain within [−0.4, 0.4] [m/s] for
all directions. All other values are initialized as 0.
"""
function quad_TMJets(; t0=0.0, T=5.0, abs_tol=1e-7,
                       orderT=5, orderQ=2, maxsteps=1000,
                       property=quad_property)

    # initial conditions, deviations and initial box
    Wpos = 0.4
    Wvel = 0.4

    q0 = zeros(12)
    δq0 = IntervalBox(-Wpos..Wpos, -Wpos..Wpos, -Wpos..Wpos, -Wvel..Wvel, -Wvel..Wvel, -Wvel..Wvel,
                      0..0, 0..0, 0..0, 0..0, 0..0, 0..0 )

    # set variables
    set_variables("δ", numvars=length(q0), order=orderQ)

    tTM, xTM = validated_integ(quadrotor!, q0, δq0, t0, T, orderQ, orderT,
                               abs_tol, maxsteps=maxsteps,
                               check_property=property)
    return tTM, xTM
end

"""
    add_time_and_project(tTM, xTM, i)

### Input

- `tTM`  -- vector with the time slices
- `xTM`  -- vector with the interval boxes
- `i`    -- index of the dimension to be plotted, e.g. if `i = 3` then return the
            reachsets for x3 as a function of time

### Notes

For plotting, to be used as

```julia
julia> include("models/ARCH/NLN/Quadrotor/quad_TMJets.jl")

julia> tTM, xTM = quad_TMJets();

julia> x3 = add_time_and_project(tTM, xTM, 3);

julia> using Plots; pyplot()

julia> begin
            plot(x3[:], legend=:none, color=:blue, lw=0, alpha=.5, xlab="t", ylab="x3")
            ylims!(-0.5, 1.5)
            plot!([0, 5], [1.4, 1.4], color=:red, legend=:none)
            plot!([0, 5], [0.98, 0.98], color=:red, legend=:none)
            plot!([0, 5], [1.02, 1.02], color=:red, legend=:none)
       end
```
"""
function add_time_and_project(tTM, xTM, i)
    N = length(tTM) # number of reachsets
    t_x = Vector{IntervalBox{2,Float64}}(undef, N-1)
    for ind in 1:N-1
        Δt = Interval(tTM[ind], tTM[ind+1])
        Δxi = union(xTM[ind][i], xTM[ind+1][i])
        t_x[ind] = IntervalBox(Δt, Δxi)
    end
    return t_x
end

max_stars_repo_name: Aadesha/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


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

JLLWrappers.@generate_main_file_header("CoolProp")
JLLWrappers.@generate_main_file("CoolProp", UUID("3351c21f-4feb-5f29-afb9-f4fcb0e27549"))
end  # module CoolProp_jll


text:
"""
    function plots_TCRDD(sys::System; kwargs...)

Runs loops of OPFs and plots a series of graphs to facilitate the analisys of the
Transmission Constrained Residual Demand Derivative (TCRDD) approach to maximise the slack
generator profit vs the actual behaviour of the system under different ammounts of bids.

Detailed description: The function runs a series of Optimal Power Flows (OPFs) for a series
of bids which range from the minimum (Pmin) power of the slack generator to the maximum
power of the slack generator (Pmax) for a set number of samples. Then for each sample the
values of the Profit, Cost, LMP, TCRDD, Residual demand and Binding Lines are calculated and
plotted. The TCRDD algorithm is taken from [1]. Conceptually with the TCRDD is possible to
create an approximation of the profit function. The Plot of the approximated profit function
is designed to take the initial solved OPF point, and create the approximated profit
function. This approximation will not be recalculated unless the value of the TCRDD changes.
The parameter "gap" allows you to force a reaproximation of the profit function for every
"gap" ammount of iterations. This helps the user to see how good or bad is the approximation
for a specific OPF point.

[1] L. Xu, R. Baldick and Y. Sutjandra, "Bidding Into Electricity Markets:
    A Transmission-Constrained Residual Demand Derivative Approach," in IEEE
    Transactions on Power Systems, vol. 26, no. 3, pp. 1380-1388, Aug. 2011,
    doi: 10.1109/TPWRS.2010.2083702.

#Arguments
- `sys::System`                             Power system in p.u. (from PowerSystems.jl)

#Keywords
- `div::Int64 = 1000`                       Number of divisions (samples) from Pmin to Pmax
- `network::DataType = StandardPTDFModel`   Network Representation of the Power System
- `solver= optimizer_with_attributes
    (Ipopt.Optimizer)`                      Solver to be used for OPF
- `tcrdd_tol::Float64 = 1e-3`               Tolerance to identify a change in the tcrdd
- `gap::Int64 = 50`                         Gap to force a recalculation of the tcrdd
- `dual_lines_tol::Float64 = 1e-1`          Value to identify any binding lines
- `dual_gen_tol::Float64 = 1e-1`            Value to identify any binding generators
- `examplecase4::Bool = false`              If running Case 4 example, provide an extra plot
"""
function plots_TCRDD(
    sys::System;
    div::Int64 = 1000,
    network::DataType = StandardPTDFModel,
    solver = optimizer_with_attributes(Ipopt.Optimizer),
    tcrdd_tol::Float64 = 1e-3,
    gap::Int64 = 50,
    dual_lines_tol::Float64 = 1e-1,
    dual_gen_tol::Float64 = 1e-1,
    examplecase4::Bool = false
    )

    # Getting Slack Generator and Bus
    # Get Slack Generator component, ID and name
    gen_thermal_slack, gen_thermal_slack_id, gen_thermal_slack_name = get_thermal_slack(sys)
    bus_slack = gen_thermal_slack.bus
    bus_slack_name = bus_slack.name
    # Get Slack Generator Costs gen_cost = αPg² + βPg + γ
    α, β = get_cost(get_variable(get_operation_cost(gen_thermal_slack)))
    γ = get_fixed(get_operation_cost(gen_thermal_slack))
    # Get System Power Base
    BaseMVA = sys.units_settings.base_value

    # Defining Plot limits
    # Get Active Power Limits
    Pmin_orig, Pmax_orig = get_active_power_limits(gen_thermal_slack)
    # Define divisions and get step
    Bidmin = Pmin_orig
    Bidmax = Pmax_orig
    step = (Bidmax - Bidmin)/div

    # Dimentionalize Arrays for Plots
    bids = Array{Float64}(undef, div+1)
    Pg_slack = Array{Float64}(undef, div+1)
    residualD_slack_Pg = Array{Float64}(undef, div+1)
    gen_profit = Array{Float64}(undef, div+1)
    gen_profit_tcrdd = Array{Float64}(undef, div+1)
    gen_profit_noCost = Array{Float64}(undef, div+1)
    gen_profit_tcrdd_noCost = Array{Float64}(undef, div+1)
    gen_cost_slack = Array{Float64}(undef, div+1)
    lmp_slack = Array{Float64}(undef, div+1)
    lmp_slack_apx_tcrdd = Array{Float64}(undef, div+1)
    tcrdd = Array{Float64}(undef, div+1)
    nbind_lines = Array{Int64}(undef, div+1)
    # TCRDD plot settings
    same_tcrdd = false
    gap_count = 0
    tcrdd0 = 0
    lmp_slack0 = 0
    Pg_slack0 = 0

    # Only for Case 4 example (needs to be declared even if empty to avoid conflicts)
    branch1_3= Array{Float64}(undef, div+1)

    # Initial Conditions
    iter_opf = 0
    # Run Loops of OPF
    for bid in Bidmin:step:Bidmax #minimum bid to maximum bid
        # Increment iterations
        iter_opf = iter_opf + 1
        # Save bid
        bids[iter_opf] = bid

        # Set the bid as the maximum power for OPF
        set_active_power_limits!(gen_thermal_slack, (min = Bidmin, max = bid))
        # Solve PTDF OPF
        lmp, res = opf_PTDF(sys; network, solver)
        # Calculate TCRDD
        tcrdd[iter_opf] = f_TCRDD(sys, res; dual_lines_tol, dual_gen_tol)

        # Getting optimised generation Outputs
        all_PGenThermal = get_variables(res)[:P__ThermalStandard]
        Pg_slack[iter_opf] = all_PGenThermal[1, gen_thermal_slack_name]
        # Calculate Residual demand
        residualD_slack_Pg[iter_opf] = residual_demand(sys, res)
        # lmp slackbus
        lmp_slack[iter_opf] = lmp[1, bus_slack_name]

        # Calculate cost of slack generator [$/hr] gen_cost = αPg² + βPg + γ
        gen_cost_slack[iter_opf] =
            α*(Pg_slack[iter_opf]*BaseMVA)^2 + β*(Pg_slack[iter_opf]*BaseMVA) + γ

        # Clear price
        clearprice = lmp_slack[iter_opf] * Pg_slack[iter_opf] * BaseMVA

        # Calculating Gen Profit of Slack bus
        gen_profit[iter_opf] = clearprice - gen_cost_slack[iter_opf]
        gen_profit_noCost[iter_opf] = clearprice

        # Number of Binding Lines
        bind_lines = get_binding_lines(sys, res; dual_lines_tol)
        nbind_lines[iter_opf] = length(bind_lines[:, 1])

        # Only for Case 4 example
        if examplecase4
            # Saving Loadability of branch 1-3
            branch1_3[iter_opf] = res.variable_values[:Fp__Line][1, "Line1"]
        end

        # Gen Profit TCRDD
        # Check if TCRDD is has changed within a tolerance
        if iter_opf > 1
            A_tcrdd = abs(tcrdd[iter_opf] - tcrdd[iter_opf - 1])
            if A_tcrdd < tcrdd_tol
                same_tcrdd = true
            else
                same_tcrdd = false
            end
        end

        if iter_opf == 1
            # Calculate the values for the approximation
            tcrdd0 = tcrdd[iter_opf]/(BaseMVA^2)
            lmp_slack0 = lmp_slack[iter_opf]
            Pg_slack0 = Pg_slack[iter_opf]
            # Evaluate the approximation for the current Pg_slack
            lmp_slack_apx_tcrdd[iter_opf] =
                ( (1/tcrdd0) * (Pg_slack[iter_opf] - Pg_slack0) ) + lmp_slack0
            clearprice_apx_tcrdd =
                lmp_slack_apx_tcrdd[iter_opf] * Pg_slack[iter_opf] * BaseMVA
            gen_profit_tcrdd[iter_opf] = clearprice_apx_tcrdd - gen_cost_slack[iter_opf]
            gen_profit_tcrdd_noCost[iter_opf] = clearprice_apx_tcrdd
            gap_count = 1

        elseif iter_opf > 1 && same_tcrdd
            # Don't recalculte the approximation, just evaluate it for the current Pg_slack
            gap_count = gap_count +1;
            if gap_count < gap
                lmp_slack_apx_tcrdd[iter_opf] =
                    ( (1/tcrdd0) * (Pg_slack[iter_opf] - Pg_slack0) ) + lmp_slack0
                clearprice_apx_tcrdd =
                    lmp_slack_apx_tcrdd[iter_opf] * Pg_slack[iter_opf] * BaseMVA
                gen_profit_tcrdd[iter_opf] = clearprice_apx_tcrdd - gen_cost_slack[iter_opf]
                gen_profit_tcrdd_noCost[iter_opf] = clearprice_apx_tcrdd
            else
                # The TCRDD is the same, but the gap has been reached,
                # So force a re-calculation of the values for the approximation
                tcrdd0 = tcrdd[iter_opf]/(BaseMVA^2);
                lmp_slack0 = lmp_slack[iter_opf];
                Pg_slack0 = Pg_slack[iter_opf];
                # Evaluate the approximation for the current Pg_slack
                lmp_slack_apx_tcrdd[iter_opf] =
                    ( (1/tcrdd0) * (Pg_slack[iter_opf] - Pg_slack0) ) + lmp_slack0
                clearprice_apx_tcrdd =
                    lmp_slack_apx_tcrdd[iter_opf] * Pg_slack[iter_opf] * BaseMVA
                gen_profit_tcrdd[iter_opf] = clearprice_apx_tcrdd - gen_cost_slack[iter_opf]
                gen_profit_tcrdd_noCost[iter_opf] = clearprice_apx_tcrdd
                # Restart the gap count
                gap_count = 0
            end

        elseif iter_opf > 1 && !same_tcrdd
            # The TCRDD changed, so re-calculate the values for the approximation
            tcrdd0 = tcrdd[iter_opf]/(BaseMVA^2)
            lmp_slack0 = lmp_slack[iter_opf]
            Pg_slack0 = Pg_slack[iter_opf]
            # Evaluate the approximation for the current Pg_slack
            lmp_slack_apx_tcrdd[iter_opf] =
                ( (1/tcrdd0) * (Pg_slack[iter_opf] - Pg_slack0) ) + lmp_slack0
            clearprice_apx_tcrdd =
                lmp_slack_apx_tcrdd[iter_opf] * Pg_slack[iter_opf] * BaseMVA
            gen_profit_tcrdd[iter_opf] = clearprice_apx_tcrdd - gen_cost_slack[iter_opf]
            gen_profit_tcrdd_noCost[iter_opf] = clearprice_apx_tcrdd
            # Restart the gap count
            gap_count = 0
        end
    end

    # Display Plots
    disp_plots(
        bids,
        gen_profit,
        gen_profit_tcrdd,
        gen_profit_noCost,
        gen_profit_tcrdd_noCost,
        gen_cost_slack,
        Pg_slack,
        lmp_slack,
        lmp_slack_apx_tcrdd,
        residualD_slack_Pg,
        tcrdd,
        nbind_lines,
        div;
        examplecase4,
        branch1_3,
        )

    return (nothing)
end

max_stars_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_issues_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_forks_repo_name: invenia/MaximumGeneratorProfitTCRDD


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 07293e64-8d8f-11eb-22d5-177de857c830
using Revise

# ╔═╡ 113d66b4-8d8f-11eb-0f3f-91dccc68148b
using OmegaCore

# ╔═╡ 15a29128-8d8f-11eb-2482-1fd2c6a81e5b
using Distributions

# ╔═╡ 71e48cca-8d91-11eb-06be-83e5843a8052
using Plots

# ╔═╡ 93e8a4d6-8d97-11eb-3bc3-57b4d317b063
md"# Latent Confounders"

# ╔═╡ 6863b6da-8daf-11eb-1c79-69ff3008f7a4
md"""A common query that we want to evaluate is th"""

# ╔═╡ 6e9c4090-8d90-11eb-14e0-516763c50379

# ╔═╡ 75a9f7f8-8d91-11eb-394a-d34b0788f105
gr()

# ╔═╡ 01f47360-8d90-11eb-0c84-d1ec8bb940d2
md"People who consistently smoked an average of less than one cigarette per day over their lifetimes had a 64 percent higher risk of earlier death than never smokers, and those who smoked between one and 10 cigarettes a day had an 87"

# ╔═╡ 3aa8b382-8d8f-11eb-2080-3fe9076830c4
basecancer = 2/1000

# ╔═╡ 19429ae4-8d8f-11eb-250d-4d8b713bca11
nsmoke = 1 ~ Poisson(10)

# ╔═╡ 7eda4c30-8d91-11eb-0d49-ff1f78b66095
histogram(randsample(nsmoke, 1000))

# ╔═╡ fe6fd072-8d92-11eb-0a1c-a5eee1623839
linear(x, m, c) = m*x + c

# ╔═╡ 11bcd9fe-8d93-11eb-3cfa-ed9ba46c74e0
c = 10 ~ Uniform(0.0, 64.0)

# ╔═╡ 16128454-8d93-11eb-3954-5dea0be6c6d7
m = 2 ~ Uniform((87.0 - 64.0)/10.0, 87.0 - 80.0)

# ╔═╡ 208be9e6-8d90-11eb-0fb0-bdc23ba88e5e
imultiplier(ω) = nsmoke(ω) * m(ω) + c(ω)

# ╔═╡ 5523ba5e-8d94-11eb-1ad0-f15982a3766a
histogram(randsample(imultiplier, 1000);
		  title = "Multiplicative factor: increase in cancer rate due to smoking",
           label = false)

# ╔═╡ 980054c2-8daf-11eb-216c-e7e0bd422fff
plot([x -> x*m_ + c_ for (m_, c_) in randsample((m, c), 100)], 0, 100;
	 label = false,
	 color=:blues,
	 title="Function samples",
	 xlabel = "Number of cigarettes smoked over lifetime",
     ylabel = "Multiplicative increase in cancer rate")

# ╔═╡ ee40c62a-8db0-11eb-1363-fd2bdbafdbfc
md"The probability of acquiring cancer is a random variable: the base rate multiplied by the multiplicative factor"

# ╔═╡ 7b8daf52-8d93-11eb-14da-3999f86e9ffd
p(ω) = min(1.0, imultiplier(ω) * basecancer)

# ╔═╡ 1de871d2-8d95-11eb-2d74-e3146bdb32ed
histogram(randsample(p, 1000);
		  title = "Probability of cancer",
          label = false)

# ╔═╡ bda18ac4-8d90-11eb-3560-8b9d65d24aa7
cancer = 3 ~ Bernoulli(p)

# ╔═╡ 64864aac-8d8f-11eb-0cc1-3f12e9d3bf22
prob(x; n = 1000) = mean(randsample(x, n))

# ╔═╡ 6497eeb8-8db1-11eb-2b5c-c314d28ebc53
md"The *treatment effect* is"

# ╔═╡ 88b5094e-8d90-11eb-0ebd-ff49cf28fa0c
treatment_effect = cancer -ₚ intervene(cancer, nsmoke => nsmoke *ₚ 0.5);

# ╔═╡ b5b67b82-8db1-11eb-29e0-0d62a2f307c0
md"The average treatment effect is simply the expectation: $E(t)$"

# ╔═╡ ab056dd6-8d95-11eb-34ce-2b9540520c55
ate = prob(treatment_effect)

# ╔═╡ b61b3768-8d96-11eb-2971-7b6bd7f468b1
lprob(x) = ω -> prob(x(ω))

# ╔═╡ e8292952-8db1-11eb-1208-21ed90d499e3
md"""We can use the random conditional distribution (or in this case, equivalently but more efficiently the random interventional distribution) to compute the distribution over the average treatment effect (which is a distribution over probabilities)"""

# ╔═╡ 9a5887c4-8d96-11eb-0d06-ed7b4b6bbc69
ate_high = lprob(rid(treatment_effect, m))

# ╔═╡ cec6ff0e-8d96-11eb-3683-1d894ce6524f
histogram(randsample(ate_high, 1000))

# ╔═╡ 299258b4-8db2-11eb-1fb2-594026226ef3
ate_high_c = lprob(rid(treatment_effect, c))

# ╔═╡ 41fa4df8-8db2-11eb-13f3-bd0d4c27e269
histogram(randsample(ate_high_c, 1000))

# ╔═╡ Cell order:
# ╟─93e8a4d6-8d97-11eb-3bc3-57b4d317b063
# ╠═6863b6da-8daf-11eb-1c79-69ff3008f7a4
# ╠═07293e64-8d8f-11eb-22d5-177de857c830
# ╠═113d66b4-8d8f-11eb-0f3f-91dccc68148b
# ╠═15a29128-8d8f-11eb-2482-1fd2c6a81e5b
# ╠═6e9c4090-8d90-11eb-14e0-516763c50379
# ╠═71e48cca-8d91-11eb-06be-83e5843a8052
# ╠═75a9f7f8-8d91-11eb-394a-d34b0788f105
# ╟─01f47360-8d90-11eb-0c84-d1ec8bb940d2
# ╠═3aa8b382-8d8f-11eb-2080-3fe9076830c4
# ╠═19429ae4-8d8f-11eb-250d-4d8b713bca11
# ╠═7eda4c30-8d91-11eb-0d49-ff1f78b66095
# ╠═fe6fd072-8d92-11eb-0a1c-a5eee1623839
# ╠═11bcd9fe-8d93-11eb-3cfa-ed9ba46c74e0
# ╠═16128454-8d93-11eb-3954-5dea0be6c6d7
# ╠═208be9e6-8d90-11eb-0fb0-bdc23ba88e5e
# ╠═5523ba5e-8d94-11eb-1ad0-f15982a3766a
# ╠═980054c2-8daf-11eb-216c-e7e0bd422fff
# ╠═ee40c62a-8db0-11eb-1363-fd2bdbafdbfc
# ╠═7b8daf52-8d93-11eb-14da-3999f86e9ffd
# ╠═1de871d2-8d95-11eb-2d74-e3146bdb32ed
# ╠═bda18ac4-8d90-11eb-3560-8b9d65d24aa7
# ╠═64864aac-8d8f-11eb-0cc1-3f12e9d3bf22
# ╠═6497eeb8-8db1-11eb-2b5c-c314d28ebc53
# ╠═88b5094e-8d90-11eb-0ebd-ff49cf28fa0c
# ╠═b5b67b82-8db1-11eb-29e0-0d62a2f307c0
# ╠═ab056dd6-8d95-11eb-34ce-2b9540520c55
# ╠═b61b3768-8d96-11eb-2971-7b6bd7f468b1
# ╠═e8292952-8db1-11eb-1208-21ed90d499e3
# ╠═9a5887c4-8d96-11eb-0d06-ed7b4b6bbc69
# ╠═cec6ff0e-8d96-11eb-3683-1d894ce6524f
# ╠═299258b4-8db2-11eb-1fb2-594026226ef3
# ╠═41fa4df8-8db2-11eb-13f3-bd0d4c27e269

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect