text:
using Distances: pairwise, Euclidean, euclidean
using Statistics: mean

# TODO: inplace
function cal_novelty(X::AM{<:AF})
    D = pairwise(Euclidean(), X, dims=2)
    return mean(D; dims=2)
end

function cal_novelty(x::AV{<:AF}, X::AM{<:AF})
    @assert length(x) == size(X, 1)
    return mean(y->euclidean(x, y), eachcol(X))
end

function debi!(de::StdMoDE, evaluator::Function, ngen; show_progress::Bool=true)::Vector{Vector{Int}}
    n = nindividuals(de)
    d = ndims(de)
    ss = senses(de)
    # initial evaluation
    de.fitness[1, :] .= evaluator.(eachcol(de.pop))
    de.fitness[2, :] .= vec(cal_novelty(de.pop))
    # de.fitness[2, :] .= 0
    # in DEMO, we will have an intermediate population of size [n, 2n]
    im_pop = zeros(eltype(de.pop), d, 2*n)
    im_fitness = zeros(eltype(de.fitness), nobjectives(de), 2*n)
    im_indices = collect(1:2*n)
    U = de.pop[:, 1] # copy to avoid allocation
    V = de.pop[:, 1]
    F = de.options.F
    Cr = de.options.Cr
    if show_progress
        pg = Progress(ngen, 0.1)
        ProgressMeter.ijulia_behavior(:clear)
    end
    fit = zeros(2)
    archive = de.pop
    for g = 1:ngen
        im_pop[:, 1:n] .= de.pop
        im_fitness[:, 1:n] .= de.fitness
        im_indices .= 1: 2*n
        cursor = n # the actual number of individuals in `im_pop`
        mean_obj = mean(@view im_fitness[1, :])
        for i = 1:n  # produce n offspring
            mutate_rand_1!(V, i, im_pop, F, @view im_indices[1:cursor]) # donor
            X = view(im_pop, :, i)  # target (base)
            bounce_back!(V, X, de.bounds)
            crossover_binomial!(U, X, V, Cr)  # trial
            fit[1] = evaluator(U)
            if fit[1] > mean_obj
                # continue
            end

            fit[2] = cal_novelty(U, view(im_pop, :, 1:cursor))
            # fit[2] = 0
            parent_fit = @view im_fitness[:, i]
            if dominate(fit, parent_fit, ss) # better, replace the parent
                im_pop[:, i] .= U
                im_fitness[:, i] .= fit
            elseif !dominate(parent_fit, fit, ss)  # equally good, add to the population
                cursor += 1
                im_pop[:, cursor] .= U
                im_fitness[:, cursor] .= fit
            end
        end
        # now we have a population of size `cursor` ∈ [n, 2n]
        # follow NSGA-II to reduce the population to size n
        # println("$g", ": $cursor")
        if cursor > n
            fronts = nondominated_sort(@view(im_fitness[:, 1:cursor]), ss)
            i_ind = 0
            i_front = 1
            # copy individuals in preferential fronts to de
            while i_ind + length(fronts[i_front]) <= n
                front = fronts[i_front]
                l = length(front)
                de.pop[:, i_ind+1:i_ind+l] .= @view im_pop[:, front]
                de.fitness[:, i_ind+1:i_ind+l] .= @view im_fitness[:, front]
                i_ind += l
                i_front += 1
            end
            # if there are still vacancies
            # we need only a portion of the next front accroding to crowding distance
            if i_ind < n
                front = fronts[i_front]
                l = length(front)
                im_indices[1:l] .= 1:l
                cds_i_front = assign_crowding_distance(view(im_fitness, :, front), [@view im_indices[1:l]])
                # arg-sort of the front in decending order w.r.t crowding distance
                ix = @view im_indices[1:l] # just a buffer
                sortperm!(ix, cds_i_front; rev=true)
                # pick the ones with larger crowding distance
                for i = 1:n-i_ind
                    picked = front[ix[i]]
                    de.pop[:, i_ind+i] .= @view im_pop[:, picked]
                    de.fitness[:, i_ind+i] .= @view im_fitness[:, picked]
                end
            end
        else # we have coincidentally n individuals
            front = Vector{Int}[]
            de.pop .= @view im_pop[:, 1:n]
            de.fitness .= @view im_fitness[:, 1:n]
        end

        if show_progress
            ProgressMeter.next!(pg; showvalues=[(:gen, "$g/$ngen")])
        end
    end
    # finally, perform nondominated sorting again to obtain fronts
    return nondominated_sort(de.fitness, ss)
end


text:
module date_naming

import Base
import Dates

function rename(cmd_message)
    @info "Starting date_naming"
    starttime     = convert_string_to_datetime(cmd_message[2])
    endtime       = convert_string_to_datetime(cmd_message[3])
    file_prefix   = cmd_message[4]
    file_postfix  = cmd_message[5]
    indir         = cmd_message[6]
    leadtimechars = parse(Int, cmd_message[7])
    maxleadtime   = parse(Int, cmd_message[8])
    outputdir     = cmd_message[9]

    make_output_dir(outputdir)

    rename_vfld_files(file_prefix, file_postfix, indir, starttime, endtime, leadtimechars, maxleadtime, outputdir)

    @info "Finished"
end

function convert_string_to_datetime(time)
    df = Dates.DateFormat("y-m-d-H")
    return Dates.DateTime(time, df)
end

function make_output_dir(dir::String)
    Base.mkpath(dir)
end

function rename_vfld_files(file_prefix::String, file_postfix::String, indir::String, starttime::Dates.DateTime, endtime::Dates.DateTime,
                         leadtimechars::Integer, maxleadtime::Integer, outdir::String)

    files = readdir(indir, join=false)
    vfld_files = [x for x in files if (startswith(x, file_prefix) && endswith(x, file_postfix))]
    vfld_files_within_range = []

    for f in vfld_files
        datestring = f[lastindex(f)-11:lastindex(f)-leadtimechars]
        leadtime = parse(Int, f[lastindex(f)-leadtimechars+1:lastindex(f)])
        dl = Dates.DateTime(datestring, "yyyymmddHH")
        if dl >= starttime && dl<endtime && leadtime <= maxleadtime
           dl_new = dl + Dates.Hour(leadtime)
           ff_new = file_prefix * Dates.format(dl_new, "yyyymmddHH")
           Base.cp(indir*f, outdir*ff_new)
        end
    end
end

end

ext: jl


text:
"""
    ModelCheckingSolver

A probabilistic model checker for MDPs and POMDPs with LTL specification.
The solver takes as input an LTL formula and the underlying MDP/POMDP planning algorithm used to perform the model checking.
It supports any solver from POMDPs.jl.
Internally, this solver requires a discrete state and discrete action model.

# Fields
- `property::SpotFormula`
- `solver::Solver` any MDP/POMDP solver
- `tolerance::Float64 = 1e-3`
- `verbose::Bool = true`
"""
@with_kw mutable struct ModelCheckingSolver <: Solver
    property::SpotFormula = ltl"true"
    solver::Solver = ValueIterationSolver() # can use any solver that returns a value function :o
    tolerance::Float64 = 1e-3
    verbose::Bool = false
end

mutable struct ModelCheckingPolicy{P <: Policy, M <:Union{ProductMDP, ProductPOMDP}, Q} <: Policy
    policy::P
    problem::M
    memory::Q
end

# NOTE: the returned value function will be 0. at the accepting states instead of 1, this is overriden by the implementation
# of POMDPs.value below but only for MDPs. A solution would be to add a sink state
function POMDPs.solve(solver::ModelCheckingSolver, problem::M) where M<:Union{MDP,POMDP}
    verbose = solver.verbose
    # parse formula first
    automata = DeterministicRabinAutomata(solver.property)
    pmdp = nothing
    sink_state = ProductState(first(states(problem)), -1)
    # build product with automata
    if isa(problem, POMDP)
        pmdp = ProductPOMDP(problem, automata, sink_state, 1.0 - solver.tolerance)
    else
        pmdp = ProductMDP(problem, automata, sink_state, 1.0 - solver.tolerance)
    end
    if isempty(pmdp.accepting_states)
        mec_time = @elapsed begin
            accepting_states!(pmdp, verbose=verbose) # compute the maximal end components via a graph analysis
            verbose ? println("Found ", length(pmdp.accepting_states), " accepting states") : nothing
        end
    end
    println("MEC comp time: ", mec_time)
    policy = solve(solver.solver, pmdp) # solve using your favorite method
    return ModelCheckingPolicy(policy, pmdp, get_init_state_number(automata))
end

# For MDPs

function POMDPs.action(policy::ModelCheckingPolicy{P, M}, s) where {P <: Policy, M <: ProductMDP}
    a = action(policy.policy, ProductState(s, policy.memory))
    update_memory!(policy, s, a)
    return a
end

function POMDPs.action(policy::ModelCheckingPolicy{P, M}, s::ProductState{S, Q}) where {P <: Policy, M <: ProductMDP, S, Q}
    return action(policy.policy, s)
end

function POMDPs.value(policy::ModelCheckingPolicy{P, M}, s) where {P <: Policy, M <: ProductMDP}
    return value(policy, ProductState(s, policy.memory))
end

function POMDPs.value(policy::ModelCheckingPolicy{P, M}, s::ProductState{S, Q}) where {P <: Policy, M <: ProductMDP, S, Q}
    if s ∈ policy.problem.accepting_states  # see comment in POMDPs.solve
        return 1.0
    else
        return value(policy.policy, s)
    end
end

function POMDPPolicies.actionvalues(policy::ModelCheckingPolicy{P, M}, s) where {P <: Policy, M <: ProductMDP}
    return actionvalues(policy, ProductState(s, policy.memory))
end

function POMDPPolicies.actionvalues(policy::ModelCheckingPolicy{P, M}, s::ProductState{S, Q}) where {P <: Policy, M <: ProductMDP, S, Q}
    if s ∈ policy.problem.accepting_states
        return ones(length(actions(policy.problem)))
    else
        return actionvalues(policy.policy, s)
    end
end

function reset_memory!(policy::ModelCheckingPolicy)
    policy.memory = get_init_state_number(policy.problem.automata)
end

function update_memory!(policy::ModelCheckingPolicy, s, a)
    l = labels(policy.problem.problem, s, a)
    qp = nextstate(policy.problem.automata, policy.memory, l)
    # println("current automata state: " , policy.memory, "  next state " , qp)
    if qp != nothing
        policy.memory = qp
    end
end

# For POMDPs

POMDPs.action(policy::ModelCheckingPolicy{P, M}, b) where {P <: Policy, M <: ProductPOMDP} = action(policy.policy, b)
POMDPs.value(policy::ModelCheckingPolicy{P, M}, b) where {P <: Policy, M <: ProductPOMDP} = value(policy.policy, b)
POMDPPolicies.actionvalues(policy::ModelCheckingPolicy{P, M}, b) where {P <: Policy, M <: ProductPOMDP} = actionvalues(policy.policy, b)


text:
function HarrisPriester(arg0::PVCoordinatesProvider, arg1::OneAxisEllipsoid)
    return HarrisPriester((PVCoordinatesProvider, OneAxisEllipsoid), arg0, arg1)
end

function HarrisPriester(arg0::PVCoordinatesProvider, arg1::OneAxisEllipsoid, arg2::Vector{Vector{jdouble}})
    return HarrisPriester((PVCoordinatesProvider, OneAxisEllipsoid, Vector{Vector{jdouble}}), arg0, arg1, arg2)
end

function HarrisPriester(arg0::PVCoordinatesProvider, arg1::OneAxisEllipsoid, arg2::Vector{Vector{jdouble}}, arg3::jdouble)
    return HarrisPriester((PVCoordinatesProvider, OneAxisEllipsoid, Vector{Vector{jdouble}}, jdouble), arg0, arg1, arg2, arg3)
end

function HarrisPriester(arg0::PVCoordinatesProvider, arg1::OneAxisEllipsoid, arg2::jdouble)
    return HarrisPriester((PVCoordinatesProvider, OneAxisEllipsoid, jdouble), arg0, arg1, arg2)
end

function get_density(obj::HarrisPriester, arg0::AbsoluteDate, arg1::Vector3D, arg2::Frame)
    return jcall(obj, "getDensity", jdouble, (AbsoluteDate, Vector3D, Frame), arg0, arg1, arg2)
end

function get_density(obj::HarrisPriester, arg0::FieldAbsoluteDate, arg1::FieldVector3D, arg2::Frame)
    return jcall(obj, "getDensity", RealFieldElement, (FieldAbsoluteDate, FieldVector3D, Frame), arg0, arg1, arg2)
end

function get_density(obj::HarrisPriester, arg0::Vector3D, arg1::FieldVector3D)
    return jcall(obj, "getDensity", RealFieldElement, (Vector3D, FieldVector3D), arg0, arg1)
end

function get_density(obj::HarrisPriester, arg0::Vector3D, arg1::Vector3D)
    return jcall(obj, "getDensity", jdouble, (Vector3D, Vector3D), arg0, arg1)
end

function get_frame(obj::HarrisPriester)
    return jcall(obj, "getFrame", Frame, ())
end

function get_max_alt(obj::HarrisPriester)
    return jcall(obj, "getMaxAlt", jdouble, ())
end

function get_min_alt(obj::HarrisPriester)
    return jcall(obj, "getMinAlt", jdouble, ())
end

function get_tab_density(obj::HarrisPriester)
    return jcall(obj, "getTabDensity", Vector{Vector{jdouble}}, ())
end


text:
"""
    $SIGNATURES

Construct (new) subgraph and draw the subgraph associated with clique `frontalSym::Symbol`.

Notes
- See `drawGraphCliq`/`writeGraphPdf` for details on keyword options.

Related

drawGraphCliq, spyCliqMat, drawTree, buildCliqSubgraphUp, buildSubgraphFromLabels!
"""
function drawCliqSubgraphUpMocking(fgl::G,
                                   treel::AbstractBayesTree,
                                   frontalSym::Symbol;
                                   show::Bool=true,
                                   filepath::String="/tmp/caesar/random/cliq_sfg.dot",
                                   engine::AS1="sfdp",
                                   viewerapp::AS2="xdot"  ) where {G <: AbstractDFG, AS1 <: AbstractString, AS2 <: AbstractString}
  #
  sfg = buildCliqSubgraphUp(fgl, treel, frontalSym)
  drawGraph(sfg, show=show, viewerapp=viewerapp, engine=engine, filepath=filepath)
  nothing
end

"""
    $SIGNATURES

Draw and show the factor graph `<:AbstractDFG` via system graphviz and xdot app.

Notes
- Requires system install on Linux of `sudo apt-get install xdot`
- Should not be calling outside programs.
- Need long term solution
- DFG's `toDotFile` a better solution -- view with `xdot` application.
- also try `engine={"sfdp","fdp","dot","twopi","circo","neato"}`

Notes:
- Calls external system application `xdot` to read the `.dot` file format
  - ```toDot(fg,file=...); @async run(`xdot file.dot`)```

Related

drawGraphCliq, [`drawTree`](@ref), printCliqSummary, spyCliqMat
"""
function drawGraph( fgl::AbstractDFG;
                    viewerapp::AbstractString="xdot",
                    filepath::AbstractString="/tmp/caesar/random/fg.dot",
                    engine::AbstractString="neato", #sfdp
                    show::Bool=true )
  #
  mkpath(dirname(filepath))
  #   mkpath(joinpath( "/", (split(filepath, '/')[1:(end-1)])...) )

  @debug "Writing factor graph file"
  fext = split(filepath, '.')[end]
  fpwoext = filepath[1:(end-length(fext)-1)] # split(filepath, '.')[end-1]
  dotfile = fpwoext*".dot"

  # create the dot file
  DFG.toDotFile(fgl, dotfile)

  try
    # run(`$(engine) $(dotfile) -T$(fext) -o $(filepath)`)
    show ? (@async run(`$(viewerapp) $(dotfile)`)) : nothing
  catch e
    @warn "not able to show $(filepath) with viewerapp=$(viewerapp). Exception e=$(e)"
  end
  nothing
end

"""
    $SIGNATURES

Draw the factor graph from a clique state machine history at a particular step as pdf and show.

Related

drawCliqSubgraphUpMocking, drawGraph, drawTree
"""
function drawGraphCliq( hists::Dict{Int, <: Tuple},
                        step::Int,
                        tree::AbstractBayesTree,
                        frontal::Symbol;
                        show::Bool=true  )
  #
  cid = getId(getClique(tree, frontal))
  cfg = hists[cid][step][4].cliqSubFg
  drawGraph(cfg, show=show)
end

"""
    $SIGNATURES

Print basic statistics about a clique variables and factors.

Related

printCliqHistorySummary
"""
function printCliqSummary(dfg::G,
                          cliq::TreeClique,
                          logger=ConsoleLogger() ) where G <: AbstractDFG
  #
  frtl = getCliqFrontalVarIds(cliq)
  seps = getCliqSeparatorVarIds(cliq)
  fcts = getCliqFactorIdsAll(cliq)

  isinit = map(x->isInitialized(dfg,x), [frtl;seps])
  # infdim = map(x->getVariableInferredDim(dfg, x), [frtl;seps])

  with_logger(logger) do
    @info "Clique $(getId(cliq)) summary:"
    @info "  num frontals:    $(length(frtl))"
    @info "  num separators:  $(length(seps))"
    @info "  num factors:     $(length(fcts))"
    @info "  num initialized: $(sum(isinit)) of $(length(isinit))"
    @info ""
    @info "  frontals:  $(frtl)"
    @info "  separator: $(seps)"
    @info "  factors:   $(fcts)"
    @info "  init'ed:   $(Int.(isinit))"
    # @info "  infr'dims: $(infdim)"
  end
  nothing
end

#


text:
println(3+4)
☎ = 2
println(☎)
println(rationalize(1.20))

max_stars_repo_name: dan-seol/Julia-playground

max_issues_repo_name: dan-seol/Julia-playground

max_forks_repo_name: dan-seol/Julia-playground


text:
category(ex) =
    ismethod(ex) ? :method :
    ismacro(ex)  ? :macro  :
    istype(ex)   ? :type   :
    isglobal(ex) ? :global :
    issymbol(ex) ? :symbol :
    error("@doc: cannot document object:\n$(ex)")

ismethod(ex) = isexpr(ex, [:function, :(=)]) && isexpr(ex.args[1], :call)
isglobal(ex) = isexpr(ex, [:global, :const, :(=)]) && !isexpr(ex.args[1], :call)
istype(ex)   = isexpr(ex, [:type, :abstract, :typealias])
ismacro(ex)  = isexpr(ex, :macro)

# handle module/function as symbols at later stage of execution
issymbol(s::Symbol) = true
issymbol(ex) = false

guess(f::Function) = :function
guess(m::Module)   = :module
guess(unknown)     = error("@doc: cannot document a $(unknown)")

function lateguess(curmod, symb)
    isdefined(curmod, symb) || error("@doc: undefined object: $(symb)")
    guess(getfield(curmod, symb))
end

function lastmethod(fn)
    res = nothing
    for f in fn.env; res = f; end
    res
end

name(ex::Expr) = name(isa(ex.args[1], Bool) ? ex.args[2] : ex.args[1])
name(s::Symbol) = s

function separate(expr)
    data, obj = expr.args
    (data,), obj
end

function separate(docs, expr)
    meta, obj = expr.args
    (docs, meta), obj
end

const METADATA = :__METADATA__

## macros –––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––

macro tex_mstr(text)
    triplequoted(text)
end

@docref () -> REF_DOCSTRINGS_MACRO
macro docstrings()
    esc(:($METADATA = Docile.Documentation(current_module())))
end

@docref () -> REF_DOC_MACRO
macro doc(args...)
    isexpr(last(args), :(->)) || error("@doc: use `->` to separate docs/object:\n$(args)")

    # Separate out the documentation and metadata from the object.
    data, obj = separate(args...)

    # Find the category and name of an object. Build corresponding quoted expressions
    # for use in the `quote` returned. Macros names are prefixed by `@` here.
    c, n   = category(obj.args[2]), name(obj.args[2])
    qc, qn = Expr(:quote, c), Expr(:quote, c == :macro ? symbol("@$(n)") : n)

    # Capture the line and file.
    loc    = obj.args[1].args
    source = (loc[1], string(loc[2]))

    # Prebuilt expressions to avoid packing lines into the destination module.
    cat = :($qc == :symbol ? Docile.lateguess(current_module(), $qn) : $qc)
    var = :($qc == :method ? Docile.lastmethod($n) : $qc in (:type, :symbol) ? $n : $qn)

    esc(:($obj; push!($METADATA, $var, Docile.Entry{$cat}($source, $(data...)))))
end


text:
using Pkg
Pkg.activate(".")
using Distributed, ClusterManagers

# Start 15 workers on SLURM - uses the allocated ones
n_tasks = 15
ClusterManagers.addprocs_slurm(
    n_tasks;
    exeflags=["--project=."]
)

# Activate the current Julia environment
@everywhere using Pkg;
@everywhere Pkg.activate(".")
@everywhere using BenchmarkTools;

# Gives all workers access to the source code (walk1D function)
@everywhere include("parallel_code_in_julia.jl")

N = 4096
T = 10^6
# Example of getting the array of positions:
@btime positions = pmap(walk1D, T for _ in 1:N);

# Example of a reduction
mean_position = @distributed (+) for _ = 1:N
    walk1D(T)
end
mean_position /= N
println("Mean Position is $mean_position")

max_stars_repo_name: JamieMair/intro-to-julia-talk

max_issues_repo_name: JamieMair/intro-to-julia-talk

max_forks_repo_name: JamieMair/intro-to-julia-talk


text:
function update!(::Val{nstates}, ::Val{Np}, rhs::CuArray, Q, elems, rka, rkb,
                 dt) where {nstates, Np}
  nelem = length(elems)
  @cuda(threads=Np, blocks=nelem,
        knl_update!(Val(nstates), Val(Np), rhs, Q, nelem, rka, rkb, dt))
end

function knl_update!(::Val{nstates}, ::Val{Np}, rhs, Q, nelem, rka, rkb,
                     dt) where {nstates, Np}
  n = threadIdx().x
  e = blockIdx().x

  @inbounds if n <= Np && e <= nelem
    for s = 1:nstates
      Q[n, s, e] += rkb * dt * rhs[n, s, e]
      rhs[n, s, e] *= rka
    end
  end
  nothing
end

max_stars_repo_name: christophernhill-dev/CLIMA

max_issues_repo_name: christophernhill-dev/CLIMA

max_forks_repo_name: christophernhill-dev/CLIMA


text:
include("tests.jl")

runtests(slow_tests = !("fast" in ARGS))


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct AvailableProvidersListState <: SwaggerModel
    stateName::Any # spec type: Union{ Nothing, String } # spec name: stateName
    providers::Any # spec type: Union{ Nothing, Vector{String} } # spec name: providers
    cities::Any # spec type: Union{ Nothing, Vector{AvailableProvidersListCity} } # spec name: cities

    function AvailableProvidersListState(;stateName=nothing, providers=nothing, cities=nothing)
        o = new()
        validate_property(AvailableProvidersListState, Symbol("stateName"), stateName)
        setfield!(o, Symbol("stateName"), stateName)
        validate_property(AvailableProvidersListState, Symbol("providers"), providers)
        setfield!(o, Symbol("providers"), providers)
        validate_property(AvailableProvidersListState, Symbol("cities"), cities)
        setfield!(o, Symbol("cities"), cities)
        o
    end
end # type AvailableProvidersListState

const _property_map_AvailableProvidersListState = Dict{Symbol,Symbol}(Symbol("stateName")=>Symbol("stateName"), Symbol("providers")=>Symbol("providers"), Symbol("cities")=>Symbol("cities"))
const _property_types_AvailableProvidersListState = Dict{Symbol,String}(Symbol("stateName")=>"String", Symbol("providers")=>"Vector{String}", Symbol("cities")=>"Vector{AvailableProvidersListCity}")
Base.propertynames(::Type{ AvailableProvidersListState }) = collect(keys(_property_map_AvailableProvidersListState))
Swagger.property_type(::Type{ AvailableProvidersListState }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_AvailableProvidersListState[name]))}
Swagger.field_name(::Type{ AvailableProvidersListState }, property_name::Symbol) =  _property_map_AvailableProvidersListState[property_name]

function check_required(o::AvailableProvidersListState)
    true
end

function validate_property(::Type{ AvailableProvidersListState }, name::Symbol, val)
end


text:
# Julia wrapper for header: /opt/local/include/gdal_frmts.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function GDALRegister_GDB()
    ccall((:GDALRegister_GDB,libgdal),Void,())
end

function GDALRegister_GTiff()
    ccall((:GDALRegister_GTiff,libgdal),Void,())
end

function GDALRegister_GXF()
    ccall((:GDALRegister_GXF,libgdal),Void,())
end

function GDALRegister_OGDI()
    ccall((:GDALRegister_OGDI,libgdal),Void,())
end

function GDALRegister_HFA()
    ccall((:GDALRegister_HFA,libgdal),Void,())
end

function GDALRegister_AAIGrid()
    ccall((:GDALRegister_AAIGrid,libgdal),Void,())
end

function GDALRegister_GRASSASCIIGrid()
    ccall((:GDALRegister_GRASSASCIIGrid,libgdal),Void,())
end

function GDALRegister_AIGrid()
    ccall((:GDALRegister_AIGrid,libgdal),Void,())
end

function GDALRegister_AIGrid2()
    ccall((:GDALRegister_AIGrid2,libgdal),Void,())
end

function GDALRegister_CEOS()
    ccall((:GDALRegister_CEOS,libgdal),Void,())
end

function GDALRegister_SAR_CEOS()
    ccall((:GDALRegister_SAR_CEOS,libgdal),Void,())
end

function GDALRegister_SDTS()
    ccall((:GDALRegister_SDTS,libgdal),Void,())
end

function GDALRegister_ELAS()
    ccall((:GDALRegister_ELAS,libgdal),Void,())
end

function GDALRegister_EHdr()
    ccall((:GDALRegister_EHdr,libgdal),Void,())
end

function GDALRegister_GenBin()
    ccall((:GDALRegister_GenBin,libgdal),Void,())
end

function GDALRegister_PAux()
    ccall((:GDALRegister_PAux,libgdal),Void,())
end

function GDALRegister_ENVI()
    ccall((:GDALRegister_ENVI,libgdal),Void,())
end

function GDALRegister_DOQ1()
    ccall((:GDALRegister_DOQ1,libgdal),Void,())
end

function GDALRegister_DOQ2()
    ccall((:GDALRegister_DOQ2,libgdal),Void,())
end

function GDALRegister_DTED()
    ccall((:GDALRegister_DTED,libgdal),Void,())
end

function GDALRegister_MFF()
    ccall((:GDALRegister_MFF,libgdal),Void,())
end

function GDALRegister_HKV()
    ccall((:GDALRegister_HKV,libgdal),Void,())
end

function GDALRegister_PNG()
    ccall((:GDALRegister_PNG,libgdal),Void,())
end

function GDALRegister_DDS()
    ccall((:GDALRegister_DDS,libgdal),Void,())
end

function GDALRegister_GTA()
    ccall((:GDALRegister_GTA,libgdal),Void,())
end

function GDALRegister_JPEG()
    ccall((:GDALRegister_JPEG,libgdal),Void,())
end

function GDALRegister_JPEG2000()
    ccall((:GDALRegister_JPEG2000,libgdal),Void,())
end

function GDALRegister_JP2KAK()
    ccall((:GDALRegister_JP2KAK,libgdal),Void,())
end

function GDALRegister_JPIPKAK()
    ccall((:GDALRegister_JPIPKAK,libgdal),Void,())
end

function GDALRegister_MEM()
    ccall((:GDALRegister_MEM,libgdal),Void,())
end

function GDALRegister_JDEM()
    ccall((:GDALRegister_JDEM,libgdal),Void,())
end

function GDALRegister_RASDAMAN()
    ccall((:GDALRegister_RASDAMAN,libgdal),Void,())
end

function GDALRegister_GRASS()
    ccall((:GDALRegister_GRASS,libgdal),Void,())
end

function GDALRegister_PNM()
    ccall((:GDALRegister_PNM,libgdal),Void,())
end

function GDALRegister_GIF()
    ccall((:GDALRegister_GIF,libgdal),Void,())
end

function GDALRegister_BIGGIF()
    ccall((:GDALRegister_BIGGIF,libgdal),Void,())
end

function GDALRegister_Envisat()
    ccall((:GDALRegister_Envisat,libgdal),Void,())
end

function GDALRegister_FITS()
    ccall((:GDALRegister_FITS,libgdal),Void,())
end

function GDALRegister_ECW()
    ccall((:GDALRegister_ECW,libgdal),Void,())
end

function GDALRegister_JP2ECW()
    ccall((:GDALRegister_JP2ECW,libgdal),Void,())
end

function GDALRegister_ECW_JP2ECW()
    ccall((:GDALRegister_ECW_JP2ECW,libgdal),Void,())
end

function GDALRegister_FujiBAS()
    ccall((:GDALRegister_FujiBAS,libgdal),Void,())
end

function GDALRegister_FIT()
    ccall((:GDALRegister_FIT,libgdal),Void,())
end

function GDALRegister_VRT()
    ccall((:GDALRegister_VRT,libgdal),Void,())
end

function GDALRegister_USGSDEM()
    ccall((:GDALRegister_USGSDEM,libgdal),Void,())
end

function GDALRegister_FAST()
    ccall((:GDALRegister_FAST,libgdal),Void,())
end

function GDALRegister_HDF4()
    ccall((:GDALRegister_HDF4,libgdal),Void,())
end

function GDALRegister_HDF4Image()
    ccall((:GDALRegister_HDF4Image,libgdal),Void,())
end

function GDALRegister_L1B()
    ccall((:GDALRegister_L1B,libgdal),Void,())
end

function GDALRegister_LDF()
    ccall((:GDALRegister_LDF,libgdal),Void,())
end

function GDALRegister_BSB()
    ccall((:GDALRegister_BSB,libgdal),Void,())
end

function GDALRegister_XPM()
    ccall((:GDALRegister_XPM,libgdal),Void,())
end

function GDALRegister_BMP()
    ccall((:GDALRegister_BMP,libgdal),Void,())
end

function GDALRegister_GSC()
    ccall((:GDALRegister_GSC,libgdal),Void,())
end

function GDALRegister_NITF()
    ccall((:GDALRegister_NITF,libgdal),Void,())
end

function GDALRegister_RPFTOC()
    ccall((:GDALRegister_RPFTOC,libgdal),Void,())
end

function GDALRegister_ECRGTOC()
    ccall((:GDALRegister_ECRGTOC,libgdal),Void,())
end

function GDALRegister_MrSID()
    ccall((:GDALRegister_MrSID,libgdal),Void,())
end

function GDALRegister_MG4Lidar()
    ccall((:GDALRegister_MG4Lidar,libgdal),Void,())
end

function GDALRegister_PCIDSK()
    ccall((:GDALRegister_PCIDSK,libgdal),Void,())
end

function GDALRegister_BT()
    ccall((:GDALRegister_BT,libgdal),Void,())
end

function GDALRegister_DODS()
    ccall((:GDALRegister_DODS,libgdal),Void,())
end

function GDALRegister_GMT()
    ccall((:GDALRegister_GMT,libgdal),Void,())
end

function GDALRegister_netCDF()
    ccall((:GDALRegister_netCDF,libgdal),Void,())
end

function GDALRegister_LAN()
    ccall((:GDALRegister_LAN,libgdal),Void,())
end

function GDALRegister_CPG()
    ccall((:GDALRegister_CPG,libgdal),Void,())
end

function GDALRegister_AirSAR()
    ccall((:GDALRegister_AirSAR,libgdal),Void,())
end

function GDALRegister_RS2()
    ccall((:GDALRegister_RS2,libgdal),Void,())
end

function GDALRegister_ILWIS()
    ccall((:GDALRegister_ILWIS,libgdal),Void,())
end

function GDALRegister_PCRaster()
    ccall((:GDALRegister_PCRaster,libgdal),Void,())
end

function GDALRegister_IDA()
    ccall((:GDALRegister_IDA,libgdal),Void,())
end

function GDALRegister_NDF()
    ccall((:GDALRegister_NDF,libgdal),Void,())
end

function GDALRegister_RMF()
    ccall((:GDALRegister_RMF,libgdal),Void,())
end

function GDALRegister_BAG()
    ccall((:GDALRegister_BAG,libgdal),Void,())
end

function GDALRegister_HDF5()
    ccall((:GDALRegister_HDF5,libgdal),Void,())
end

function GDALRegister_HDF5Image()
    ccall((:GDALRegister_HDF5Image,libgdal),Void,())
end

function GDALRegister_MSGN()
    ccall((:GDALRegister_MSGN,libgdal),Void,())
end

function GDALRegister_MSG()
    ccall((:GDALRegister_MSG,libgdal),Void,())
end

function GDALRegister_RIK()
    ccall((:GDALRegister_RIK,libgdal),Void,())
end

function GDALRegister_Leveller()
    ccall((:GDALRegister_Leveller,libgdal),Void,())
end

function GDALRegister_SGI()
    ccall((:GDALRegister_SGI,libgdal),Void,())
end

function GDALRegister_SRTMHGT()
    ccall((:GDALRegister_SRTMHGT,libgdal),Void,())
end

function GDALRegister_DIPEx()
    ccall((:GDALRegister_DIPEx,libgdal),Void,())
end

function GDALRegister_ISIS3()
    ccall((:GDALRegister_ISIS3,libgdal),Void,())
end

function GDALRegister_ISIS2()
    ccall((:GDALRegister_ISIS2,libgdal),Void,())
end

function GDALRegister_PDS()
    ccall((:GDALRegister_PDS,libgdal),Void,())
end

function GDALRegister_IDRISI()
    ccall((:GDALRegister_IDRISI,libgdal),Void,())
end

function GDALRegister_Terragen()
    ccall((:GDALRegister_Terragen,libgdal),Void,())
end

function GDALRegister_WCS()
    ccall((:GDALRegister_WCS,libgdal),Void,())
end

function GDALRegister_WMS()
    ccall((:GDALRegister_WMS,libgdal),Void,())
end

function GDALRegister_HTTP()
    ccall((:GDALRegister_HTTP,libgdal),Void,())
end

function GDALRegister_SDE()
    ccall((:GDALRegister_SDE,libgdal),Void,())
end

function GDALRegister_GSAG()
    ccall((:GDALRegister_GSAG,libgdal),Void,())
end

function GDALRegister_GSBG()
    ccall((:GDALRegister_GSBG,libgdal),Void,())
end

function GDALRegister_GS7BG()
    ccall((:GDALRegister_GS7BG,libgdal),Void,())
end

function GDALRegister_GRIB()
    ccall((:GDALRegister_GRIB,libgdal),Void,())
end

function GDALRegister_INGR()
    ccall((:GDALRegister_INGR,libgdal),Void,())
end

function GDALRegister_ERS()
    ccall((:GDALRegister_ERS,libgdal),Void,())
end

function GDALRegister_PALSARJaxa()
    ccall((:GDALRegister_PALSARJaxa,libgdal),Void,())
end

function GDALRegister_DIMAP()
    ccall((:GDALRegister_DIMAP,libgdal),Void,())
end

function GDALRegister_GFF()
    ccall((:GDALRegister_GFF,libgdal),Void,())
end

function GDALRegister_COSAR()
    ccall((:GDALRegister_COSAR,libgdal),Void,())
end

function GDALRegister_TSX()
    ccall((:GDALRegister_TSX,libgdal),Void,())
end

function GDALRegister_ADRG()
    ccall((:GDALRegister_ADRG,libgdal),Void,())
end

function GDALRegister_SRP()
    ccall((:GDALRegister_SRP,libgdal),Void,())
end

function GDALRegister_COASP()
    ccall((:GDALRegister_COASP,libgdal),Void,())
end

function GDALRegister_BLX()
    ccall((:GDALRegister_BLX,libgdal),Void,())
end

function GDALRegister_LCP()
    ccall((:GDALRegister_LCP,libgdal),Void,())
end

function GDALRegister_PGCHIP()
    ccall((:GDALRegister_PGCHIP,libgdal),Void,())
end

function GDALRegister_TMS()
    ccall((:GDALRegister_TMS,libgdal),Void,())
end

function GDALRegister_EIR()
    ccall((:GDALRegister_EIR,libgdal),Void,())
end

function GDALRegister_GEOR()
    ccall((:GDALRegister_GEOR,libgdal),Void,())
end

function GDALRegister_TIL()
    ccall((:GDALRegister_TIL,libgdal),Void,())
end

function GDALRegister_R()
    ccall((:GDALRegister_R,libgdal),Void,())
end

function GDALRegister_Rasterlite()
    ccall((:GDALRegister_Rasterlite,libgdal),Void,())
end

function GDALRegister_EPSILON()
    ccall((:GDALRegister_EPSILON,libgdal),Void,())
end

function GDALRegister_PostGISRaster()
    ccall((:GDALRegister_PostGISRaster,libgdal),Void,())
end

function GDALRegister_NWT_GRD()
    ccall((:GDALRegister_NWT_GRD,libgdal),Void,())
end

function GDALRegister_NWT_GRC()
    ccall((:GDALRegister_NWT_GRC,libgdal),Void,())
end

function GDALRegister_SAGA()
    ccall((:GDALRegister_SAGA,libgdal),Void,())
end

function GDALRegister_KMLSUPEROVERLAY()
    ccall((:GDALRegister_KMLSUPEROVERLAY,libgdal),Void,())
end

function GDALRegister_GTX()
    ccall((:GDALRegister_GTX,libgdal),Void,())
end

function GDALRegister_LOSLAS()
    ccall((:GDALRegister_LOSLAS,libgdal),Void,())
end

function GDALRegister_Istar()
    ccall((:GDALRegister_Istar,libgdal),Void,())
end

function GDALRegister_NTv2()
    ccall((:GDALRegister_NTv2,libgdal),Void,())
end

function GDALRegister_CTable2()
    ccall((:GDALRegister_CTable2,libgdal),Void,())
end

function GDALRegister_JP2OpenJPEG()
    ccall((:GDALRegister_JP2OpenJPEG,libgdal),Void,())
end

function GDALRegister_XYZ()
    ccall((:GDALRegister_XYZ,libgdal),Void,())
end

function GDALRegister_HF2()
    ccall((:GDALRegister_HF2,libgdal),Void,())
end

function GDALRegister_PDF()
    ccall((:GDALRegister_PDF,libgdal),Void,())
end

function GDALRegister_JPEGLS()
    ccall((:GDALRegister_JPEGLS,libgdal),Void,())
end

function GDALRegister_MAP()
    ccall((:GDALRegister_MAP,libgdal),Void,())
end

function GDALRegister_OZI()
    ccall((:GDALRegister_OZI,libgdal),Void,())
end

function GDALRegister_ACE2()
    ccall((:GDALRegister_ACE2,libgdal),Void,())
end

function GDALRegister_CTG()
    ccall((:GDALRegister_CTG,libgdal),Void,())
end

function GDALRegister_E00GRID()
    ccall((:GDALRegister_E00GRID,libgdal),Void,())
end

function GDALRegister_SNODAS()
    ccall((:GDALRegister_SNODAS,libgdal),Void,())
end

function GDALRegister_WEBP()
    ccall((:GDALRegister_WEBP,libgdal),Void,())
end

function GDALRegister_ZMap()
    ccall((:GDALRegister_ZMap,libgdal),Void,())
end

function GDALRegister_NGSGEOID()
    ccall((:GDALRegister_NGSGEOID,libgdal),Void,())
end

function GDALRegister_MBTiles()
    ccall((:GDALRegister_MBTiles,libgdal),Void,())
end

function GDALRegister_ARG()
    ccall((:GDALRegister_ARG,libgdal),Void,())
end

function GDALRegister_IRIS()
    ccall((:GDALRegister_IRIS,libgdal),Void,())
end

function GDALRegister_KRO()
    ccall((:GDALRegister_KRO,libgdal),Void,())
end


text:
"""
    download(mission::Mission, master::DataFrames.DataFrame, obsid::String; overwrite=false)

Finds the FTP server-side path via `_ftp_dir`, downloads folder using `lftp`, currently excludes the `uf` files
assuming calibrations are up to date. Saves download folder to local, dot-free, path
"""
function download(mission::Mission, obs_row::DataFrames.DataFrameRow{DataFrames.DataFrame}; overwrite=false)
    dir_down = _obs_path_server(mission, obs_row)
    dir_dest = _obs_path_local(mission, obs_row, kind=:download)

    @info "heasarc.gsfc.nasa.gov:$dir_down --> $dir_dest"

    download_command = `lftp heasarc.gsfc.nasa.gov -e "mirror \"$dir_down\" \"$dir_dest\" --parallel=10 --only-newer --exclude-glob *ufa.evt.gz --exclude-glob *ufa.evt --exclude-glob *uf.evt.gz && exit"`

    if _log_query(mission, obs_row, "meta", :downloaded) && !overwrite
        @warn "Skipping download, dir exists"
    else
        mkpath(dir_dest)

        println(download_command)
        try
            run(download_command)
        catch err
            if typeof(err) != InterruptException
                throw(JAXTAMError("Error downloading heasarc.gsfc.nasa.gov:$dir_down", :download, err))
            else
                rethrow(err)
            end
        end
        _log_add(mission, obs_row, Dict("meta"=>Dict(:downloaded=>true)))
    end
end

function download(mission::Mission, obs_rows::DataFrames.DataFrame; overwrite=false)
    for obs_row in DataFrames.eachrow(obs_rows)
        try
            download(mission, obs_row; overwrite=overwrite)
        catch err
            if typeof(err) == JAXTAMError
                _log_add(mission, obs_row, Dict{String,Any}("errors"=>Dict(err.step=>err)))
                @warn err
                continue
            else
                rethrow(err)
            end
        end
    end
end

function download(mission::Mission, obsid::String; overwrite=false)
    obs_rows = master_query(mission, :obsid, obsid)

    download(mission, obs_rows; overwrite=overwrite)
end


text:
using DataFrames, TableOperations, Tables, Random, TreeParzen
using Recommenders:
    Movielens1M,
    load_dataset,
    ratio_split,
    Randomwalk,
    evaluate_u2i,
    MeanPrecision,
    MeanRecall,
    MeanNDCG,
    fit!,
    predict_u2i,
    make_u2i_dataset,
    name

function main()
    ml1M = Movielens1M()
    download(ml1M)
    rating, _, _ = load_dataset(ml1M)

    rating = rating |> TableOperations.transform(Dict(:rating => x -> 1.0))

    Random.seed!(1234)
    train_valid_table, test_table = ratio_split(rating, 0.8)
    train_table, valid_table = ratio_split(train_valid_table, 0.8)

    prec10 = MeanPrecision(10)
    recall10 = MeanRecall(10)
    ndcg10 = MeanNDCG(10)
    metrics = [prec10, recall10, ndcg10]

    space = Dict(
        :terminate_prob => HP.QuantUniform(:terminate_prob, 0.1, 0.9, 0.1),
        :total_walk_length => HP.LogUniform(:total_walk_length, log(1e3), log(1e6)),
        :min_high_visited_candidates =>
            HP.Choice(:min_high_visited_candidates, [10, 50, 100, 250, 500, 1000]),
        :high_visited_count_threshold =>
            HP.Choice(:high_visited_count_threshold, [2, 4, 8, 16, 32, 64, Inf]),
        :pixie_walk_length_scaling =>
            HP.Choice(:pixie_walk_length_scaling, [true, false]),
        :pixie_multi_hit_boosting =>
            HP.Choice(:pixie_multi_hit_boosting, [true, false]),
    )

    model = Randomwalk()
    fit!(model, train_table, col_item = :movieid)
    userids_valid, gts_valid = make_u2i_dataset(valid_table, col_item = :movieid)

    function invert_output(params)
        terminate_prob = params[:terminate_prob]
        total_walk_length = round(Int, params[:total_walk_length])
        min_high_visited_candidates = convert(Int, params[:min_high_visited_candidates])
        high_visited_count_threshold = params[:high_visited_count_threshold]
        if high_visited_count_threshold < Inf
            high_visited_count_threshold = convert(Int, high_visited_count_threshold)
        end
        pixie_walk_length_scaling = params[:pixie_walk_length_scaling]
        pixie_multi_hit_boosting = params[:pixie_multi_hit_boosting]

        recoms = predict_u2i(
            model,
            userids_valid,
            10,
            col_user = :userid,
            col_item = :movieid,
            col_rating = :rating,
            drop_history = true,
            terminate_prob = terminate_prob,
            total_walk_length = total_walk_length,
            min_high_visited_candidates = min_high_visited_candidates,
            high_visited_count_threshold = high_visited_count_threshold,
            pixie_walk_length_scaling = pixie_walk_length_scaling,
            pixie_multi_hit_boosting = pixie_multi_hit_boosting,
        )
        result = Dict()
        for metric in metrics
            result[Symbol(name(metric))] = metric(recoms, gts_valid)
        end
        result = NamedTuple(result)
        @info params, result
        return -result[:ndcg10]
    end

    @info "Tuning start."
    best = fmin(invert_output, space, 100, logging_interval = -1)
    @info best

    @info "Evaluate best model."

    terminate_prob = best[:terminate_prob]
    total_walk_length = round(Int, best[:total_walk_length])
    min_high_visited_candidates = convert(Int, best[:min_high_visited_candidates])
    high_visited_count_threshold = best[:high_visited_count_threshold]
    if high_visited_count_threshold < Inf
        high_visited_count_threshold = convert(Int, high_visited_count_threshold)
    end
    pixie_walk_length_scaling = best[:pixie_walk_length_scaling]
    pixie_multi_hit_boosting = best[:pixie_multi_hit_boosting]

    model = Randomwalk()
    result = evaluate_u2i(
        model,
        train_valid_table,
        test_table,
        metrics,
        10,
        col_user = :userid,
        col_item = :movieid,
        col_rating = :rating,
        drop_history = true,
        terminate_prob = terminate_prob,
        total_walk_length = total_walk_length,
        min_high_visited_candidates = min_high_visited_candidates,
        high_visited_count_threshold = high_visited_count_threshold,
        pixie_walk_length_scaling = pixie_walk_length_scaling,
        pixie_multi_hit_boosting = pixie_multi_hit_boosting,
    )
    @info result

end

main()


text:
## Main loop to extract in e2e the DBSCAN parameters.
## The list of votable is selected directly from the directory
##
## That version is to run in parallel the optimization

## OPTIONS
# -s : first centile to consider in the vot files (0-100)
# -e : last centile to consider in the vot files (0-100)
# -m : metafile for all split

using DataFrames, Query
using CSV, Glob, Dates
using Statistics, Random
using Printf

rootdir =  ENV["GAIA_ROOT"]

push!(LOAD_PATH,"$rootdir/run/src")
using GaiaClustering

## directory
wdir    = "$rootdir/e2e_products"
votdir  = "$wdir/votable"
plotdir = "$wdir/plotSelect"
ocdir   = "$wdir/oc"

## load a liist of votable and update the file if done
## add results
##

function readlist_votable(filelist::String)
    vot = DataFrames.copy(CSV.read(filelist))
    return(vot)
end

function getdata(filevot)
    voname = filevot

    data       = read_votable(voname)
    df         = filter_data(data, [0,2100])
    dfcart     = add_cartesian(df)
    blck       = [[1,2,3],[4,5], [6,7,8]]
    wghtblck   = [4.0,5.0,1.0]
    norm       = "identity"

    dfcartnorm , scale8 = normalization_PerBlock(dfcart, blck, wghtblck , norm, false)
    return(df, dfcart , dfcartnorm)
end

## to check if done and record
## check and updt if votname analyzed. If not done return false
function updt_votcompleted(fileres, votname , cycletot=1, flag= 0 , onlycheck=true)
    let
        if  onlycheck
            if !isfile(fileres)
                return(0, false)
            else
                res = CSV.read(fileres, delim=";" , DataFrame)
                if votname in res.votname
                    x = @from i in res begin
                        @where i.votname == votname
                        @select i
                        @collect DataFrame
                    end
                    return(x, true)
                else
                    return(0, false)
                end
            end
            ## UPDTE
        else
            if !isfile(fileres)
                res = DataFrame(votname=votname, cycle=cycletot, flag=flag)
                CSV.write(fileres,res,delim=';')
                println("## $fileres created...")
                return(res,true)
            else
                res = CSV.read(fileres, delim=";", DataFrame)
                newrow = DataFrame(votname=votname,cycle=cycletot, flag=flag)
                append!(res,newrow)
                CSV.write(fileres,res,delim=';')
                return(res,true)
            end
        end
    end
end
## check for blacklist
##
function read_blacklist(blackname)
    if isfile(blackname)
        df= DataFrames.copy(CSV.read(blackname, delim=";"))
        blacklist= df.votname
    else
        blacklist= [""]
    end

    return(blacklist)
end

#######################
## Main loop
##
function main(filelist, metafile, prefile)
    let
        wd= pwd() ; nfile= size(filelist)[1] ; totalTime= 0.
        println("## Starting main loop using optimization with cycles...")
        println("## It can be very long but will be resumed to the last reduced file.")
        println("## Working directory: $wd")
        println("## $nfile files to analyze...")

        m= read_params(metafile, false)
        m.plotdir= plotdir
        m.ocdir= ocdir
        m.prefile= prefile
        fileres= "$(m.prefile).done.csv"

        # read a possible votname blacklist
        blackname= "blacklist.csv"
        blacklist= read_blacklist(blackname)
        println("## blacklist read...")

        for i in 1:nfile
            votname = filelist[i]
            res, votfound= updt_votcompleted(fileres, votname, 0, 0, true)

            ## test blacklist
            if votname in blacklist
                println("##")
                println("## $votname in blacklist...")
                println("##")
                votfound= true
            end

            if !votfound
                tstart= now()
                println("###########################")
                println("## Starting with $votname")
                println("## Starting at $tstart")

                df , dfcart , dfcartnorm = getdata(votdir*"/"*votname)
                m.votname= votname
                cycle, flag= cycle_extraction(df, dfcart, m)
                res,  votfound= updt_votcompleted(fileres, votname , cycle, flag, false)

                tend= now()
                println("## Ending at $tend")
                println("## number of cycle: $cycle , flag:$flag ")
                println("##########################")
                println("##")

                duration= Dates.value(tend-tstart) / (1000*3600)
                totalTime += duration
                meanTime= totalTime / i
                ETA= meanTime * (nfile-i) / 24
                nleft= nfile-i
                ETAstr= @sprintf("%3.3f", ETA) ; durationstr= @sprintf("%3.3f", duration)
                @printf("## %s \n",specialstr("Duration: $durationstr hours","YELLOW"))
                @printf("## %s \n",specialstr("ETA: $ETAstr days","YELLOW"))
                @printf("## %s \n",specialstr("Votable done: $votname","YELLOW"))
                @printf("## %s \n",specialstr("Files analyzed: $i","YELLOW"))
                @printf("## %s \n",specialstr("Files to go: $nleft","YELLOW"))
                println("##\n##")
            end
        end
    end
    print("## Main loop done.")
end
###############################################################################
let
    cd(votdir)
    votlist= glob("*.vot")
    cd(wdir)

    istart= 0 ; iend= 100
    metafile= "configAll.ext"

    for i in 1:length(ARGS)
        if ARGS[i] == "-s"
            istart= parse(Int, ARGS[i+1])
        end
        if ARGS[i] == "-e"
            iend= parse(Int, ARGS[i+1])
        end
        if ARGS[i] == "-m"
            metafile= ARGS[i+1]
        end
    end

    if istart < 0 istart= 0 end
    if iend > 100 iend= 100 end

###############
    header_extract()
    m= read_params(metafile, false)
    prefile= "$(m.prefile)-$istart-$iend"

    nfile= length(votlist)
    i1= convert(Int,floor(nfile*istart/100))
    if i1== 0 i1= 1 end

    if iend < 100
        i2= convert(Int,floor(nfile*iend/100)-1)
    else
        i2= convert(Int,floor(nfile*iend/100))
    end

    println("## DBSCAN optimization starting")
    println("## split files: $prefile")
    println("##")

    votsublist= votlist[i1:i2]
    main(votsublist, metafile, prefile)
    println("##\n## END of split optimization for $prefile \n##")
end

max_stars_repo_name: bosscha/gaia-shock

max_issues_repo_name: bosscha/gaia-shock

max_forks_repo_name: bosscha/gaia-shock


text:
# SPDX-License-Identifier: MIT

module BayesianModels

using CUDA
import ExpFamilyDistributions
const EFD = ExpFamilyDistributions
using LinearAlgebra

include("bmobj.jl")

export getparams, isbayesianparam, todict, fromdict
include("params/params.jl")
include("params/bayesparam.jl")
include("params/constparam.jl")

export gpu!, cpu!, init_gpu
include("utils.jl")

export loglikelihood, posterior, predict
include("models/models.jl")

export Mixture, Normal, NormalDiag
include("models/mixture.jl")
include("models/normal.jl")

export elbo, ∇elbo, gradstep
include("elbo.jl")

end # module

max_stars_repo_name: lucasondel/PPCA

max_issues_repo_name: lucasondel/PPCA

max_forks_repo_name: lucasondel/PPCA


text:
using ReHorde
using Test

@testset "ReHorde.jl" begin
    # Write your tests here.
end


text:
module LegacyMeansTestedBenefits
#
# This module implements the pre-universal credit means-tested benefit system.
#
using ArgCheck
using Base: Bool
using Dates
using ScottishTaxBenefitModel
using .Definitions

using .STBIncomes

using .Utils: to_md_table

using .ModelHousehold:
    BenefitUnit,
    Household,
    Person,
    between_ages,
    count,
    empl_status_in,
    ge_age,
    get_benefit_units,
    has_children,
    is_lone_parent,
    is_severe_disability,
    is_single,
    le_age,
    num_adults,
    num_carers,
    pers_is_disabled,
    search

using .STBParameters:
    AgeLimits,
    ChildTaxCredit,
    HoursLimits,
    HousingBenefits,
    HousingRestrictions,
    IncomeRules,
    IncomeRules,
    LegacyMeansTestedBenefitSystem,
    MinimumWage,
    NonMeansTestedSys,
    PersonalAllowances,
    Premia,
    SavingsCredit,
    WorkingTaxCredit,
    reached_state_pension_age,
    state_pension_age

using .GeneralTaxComponents:
    RateBands,
    TaxResult,
    calctaxdue

using .Results:
    BenefitUnitResult,
    HouseholdResult,
    IndividualResult,
    LMTIncomes,
    LMTResults,
    LMTCanApplyFor,
    aggregate_tax,
    gross_pension_contributions,
    has_any,
    to_string,
    total

using .Intermediate:
    MTIntermediate,
    HHIntermed,
    apply_2_child_policy,
    born_before,
    is_working_hours,
    make_recipient,
    num_born_before,
    working_disabled

using .LocalLevelCalculations:
    apply_rent_restrictions

export
    calc_legacy_means_tested_benefits,
    calc_allowances,
    calc_incomes,
    calc_NDDs,
    calc_premia,
    calculateHB_CTR!,
    calcWTC_CTC!,
    make_lmt_benefit_applicability,
    num_qualifying_for_severe_disability,
    tariff_income

"""
Incomes for old style mt benefits
The CPAG guide ch 21/21 has over 100 pages on this stuff
this can no more than catch the gist.
"""
function calc_incomes(
    which_ben :: LMTBenefitType, # esa hb is jsa pc wtc ctr pc sc
    bu        :: BenefitUnit,
    bur       :: BenefitUnitResult,
    intermed  :: MTIntermediate,
    incrules  :: IncomeRules,
    hours     :: HoursLimits ) :: LMTIncomes
    T = typeof( incrules.permitted_work )
    mntr = bur.legacy_mtbens # shortcut
    inc = LMTIncomes{T}()
    gross_earn = zero(T)
    net_earn = zero(T)
    other = zero(T)

    if which_ben in [hb,ctr]
        inclist = incrules.hb_incomes
    elseif which_ben == pc
        inclist = incrules.pc_incomes
    elseif which_ben == sc
        inclist = incrules.sc_incomes
    else
        inclist = incrules.incomes
    end
    # children's income doesn't count see cpag p421, so:
    for pid in bu.adults
        pers = bu.people[pid]
        pres = bur.pers[pid]
        gpc = gross_pension_contributions(pres)
        gross =
            pres.income[WAGES] +
            pres.income[SELF_EMPLOYMENT_INCOME] # this includes losses
        if which_ben in [pc,is,jsa,esa,hb,ctr]
            net =
                gross - ## FIXME parameterise this so we can use gross/net
                pres.it.non_savings_tax - ## FIXME?? income[INCOME_TAX] ??
                pres.income[NATIONAL_INSURANCE] -
                0.5*gpc # .income[PENSION_CONTRIBUTIONS_EMPLOYEE]
            # println( "net=$net; gross=$gross pres.it.non_savings_tax = $(pres.it.non_savings_tax) gpc=$gpc ni=$(pres.income[NATIONAL_INSURANCE])")
        else
            # wtc,ctr all pension contributions but not IT/NI
            gross = max( 0.0, gross - gross_pension_contributions(pres)) # .income[PENSION_CONTRIBUTIONS_EMPLOYEE])
            net = gross
        end
        gross_earn += gross
        net_earn += max( 0.0, net )
        other += isum( pres.income, inclist )
        #    data=pers.income,
        #    calculated=pres.incomes,
        #    included=inclist )
    end
    # disregards
    # if which_ben in [hb,jsa,is,]
    # FIXME this is not quite right for ESA
    disreg = intermed.is_sing ?  incrules.low_single : incrules.low_couple

    if( which_ben == esa )
        if ! search( bu, is_working_hours, hours.lower )
            disreg = incrules.high
            # and some others ... see CPAG
        end
    elseif which_ben in [hb,ctr,jsa,is,pc]
        if intermed.is_sparent
            disreg = which_ben == hb ? incrules.lone_parent_hb : incrules.high
        elseif ! isdisjoint( mntr.premia, [
                carer_single,
                carer_couple,
                disability_couple,
                disability_single,
                severe_disability_couple,
                severe_disability_single,
                enhanced_disability_couple,
                enhanced_disability_single] )
            disreg = incrules.high
        end
    end

    if( which_ben in [hb,ctr] )
        # fixme do this above
        if has_any( bur, CONTRIB_EMPLOYMENT_AND_SUPPORT_ALLOWANCE, NON_CONTRIB_EMPLOYMENT_AND_SUPPORT_ALLOWANCE )
            disreg = incrules.high
        end
        # HB disregard CPAG p432 this, too, is very approximate
        # work 30+ hours - should really check premia if haskeys( mtr.premia )
        extra = 0.0
        # println( "hours.higher $(hours.higher)")
        if search( bu, is_working_hours, hours.higher )
            extra = incrules.hb_additional
            # println( "extra=$extra")
        elseif search(  bu, is_working_hours, hours.lower )
            if intermed.is_sparent || (intermed.num_children > 0) || intermed.is_disabled
                extra = incrules.hb_additional
            end
        end
        disreg += extra
        # childcare in HB - costs are assigned in frs to the children
        if ( intermed.num_children > 0 )
            maxcc = intermed.num_children == 1 ? incrules.childcare_max_1 : incrules.childcare_max_2
            cost_of_childcare = 0.0
            for pid in bu.children
                cost_of_childcare += bu.people[pid].cost_of_childcare
            end
            inc.childcare = min(cost_of_childcare, maxcc )
        end
    end

    """
    not even remotely right ... cpag 21
    """

    cap = 0.0
    for pid in bu.adults
        if bu.people[pid].over_20_k_saving
            cap = 100_000 # some huge number bigger than any reasonable capmax
        else
            for (at,val) in bu.people[pid].assets
                cap += val
            end
        end
    end
    inc.other_income = other
    inc.capital = cap
    inc.gross_earnings = gross_earn
    inc.net_earnings = max(0.0, net_earn - disreg ) - inc.childcare

    capmin = incrules.capital_min
    capmax = incrules.capital_max
    tariff = incrules.capital_tariff
    if intermed.someone_pension_age
        tariff = incrules.pensioner_tariff
        if which_ben == pc
            capmin = incrules.pc_capital_min
            capmax = incrules.pc_capital_max
        else
            capmin = incrules.pensioner_capital_min
            capmax = incrules.pensioner_capital_max
        end
    end

    inc.tariff_income = tariff_income(cap, capmin, tariff )
    inc.disqualified_on_capital = cap > capmax
    inc.total_income = inc.net_earnings + inc.other_income + inc.tariff_income
    inc.disregard = disreg
    return inc
end

"""
The strategy here is to include *all* benefits the BU is entitled to
and then decide later on which ones to route to. Source: CPAG chs 9-15
'Who can get XX' sections.
"""
function make_lmt_benefit_applicability(
    mt_ben_sys   :: LegacyMeansTestedBenefitSystem,
    intermed :: MTIntermediate,
    hrs      :: HoursLimits ) :: LMTCanApplyFor
    whichb = LMTCanApplyFor()
    if intermed.someone_pension_age #!! fixme both pension age for new claims 2019=>
        whichb.pc = true
    end
    if intermed.someone_pension_age_2016
        whichb.sc = true
    end
    # println( "intermed: intermed.num_working_16_or_less=$(intermed.num_working_16_or_less) intermed.economically_active=$(intermed.economically_active) intermed.num_working_16_or_less=$(intermed.num_working_16_or_less)")
    if (! intermed.all_pension_age ) && (intermed.num_working_16_or_less >= 1)
        if ((intermed.num_adults == 1) ||
           ((intermed.num_adults == 2) && (intermed.num_working_24_plus == 0)))
            if intermed.limited_capacity_for_work
                whichb.esa = true
            # CPAG 21/2 p236/7
            elseif( intermed.num_adults == 1) && (intermed.num_allowed_children>0) &&
                ((intermed.age_youngest_child <= 5) || (intermed.age_oldest_adult < 18))
                    whichb.is = true
            elseif intermed.economically_active
                whichb.jsa = true
            else
                whichb.is = true
            end
        end
    end # all pens age

    #
    # Tax credits
    # CTC - easy - but route pensioners to pension credit even
    # if they have children - this is how Age Concern do it, anyway.
    #
    if intermed.has_children && ( ! whichb.pc )
        whichb.ctc = true
    end
    #
    # WTC - not quite so easy
    #
    if intermed.someone_working_ft_and_25_plus
        whichb.wtc = true
    elseif (intermed.total_hours_worked >= hrs.med) && (intermed.num_working_pt>0) && intermed.has_children
        # ie. 24 hrs worked total and one person  >= 16 hrs and has children
        whichb.wtc = true
    elseif (intermed.num_working_pt>0) && intermed.someone_pension_age
        whichb.wtc = true
    elseif (intermed.num_working_pt>0) && intermed.is_sparent
        whichb.wtc = true
    elseif intermed.working_disabled
        whichb.wtc = true
    end
    # FIXME not really true
    if intermed.benefit_unit_number == 1
        whichb.hb = true
        whichb.ctr = true
    end
    # check!!! is this correct? is, etc. or wtc over pc/sc
    if( whichb.esa || whichb.jsa || whichb.is ) ## || whichb.wtc)
        whichb.pc = false
        whichb.sc = false
    end
    # FIXME we have abolished code twice!
    # override these if something abolished
    if mt_ben_sys.isa_jsa_esa_abolished
        whichb.jsa = false
        whichb.is = false
        whichb.esa = false
        # FIXME problem with pension credit here
        whichb.pc = false
    end
    whichb.ctc &= (! mt_ben_sys.child_tax_credit.abolished)
    whichb.wtc &= (! mt_ben_sys.working_tax_credit.abolished)
    whichb.hb &= (! mt_ben_sys.hb.abolished)
    whichb.ctr &= (! mt_ben_sys.ctr.abolished)
    whichb.sc &= (! mt_ben_sys.savings_credit.abolished)
    return whichb
end # make_lmt_benefit_applicability

# CTC disabled child definition:
#
# Getting Disability Living Allowance (DLA), Personal Independence Payment (PIP) or Armed Forces Independence Payment (AFIP)
# In hospital but otherwise would get DLA, PIP or AFIP
# Certified blind (or were until less than 28 weeks before you made the claim)

#
#

# Disabled adult for wtc
# Answer yes if you work at least 16 hours a week and any of the following apply to you:
#
#    you get a disability benefit, eg Disability Living Allowance, Attendance Allowance or Personal Independence Payment
#    you have recently been getting disability benefit, eg Incapacity Benefit or Employment and Support Allowance
#    you have a disability that makes it hard for you to get a job, eg you’re deaf or blind
#

"""
tariff income from capital.
See CPAG p488 £1 pw for every £250, or part of £250 above 6,000
"""
function tariff_income( cap :: Real, capital_min::Real, tariff :: Real )::Real
    return ceil( max(0.0, cap-capital_min)/tariff)
end

"""
WRONG AND NOT USED!
CPAG 2019/20 p347. I *thik* this is what it means ...
Moved to own function since it's convoluted.
FIXME move this as it's being used in the UC module
"""
function num_qualifying_for_severe_disability(
    bu    :: BenefitUnit,
    bures :: BenefitUnitResult,
    num_bus :: Int ) :: Int
    if num_bus > 1
        return 0;
    end
    n = 0
    for pid in bu.adults
        pers = bu.people[pid]
        if (pers.dla_self_care_type in [high,mid] ||
              pers.attendance_allowance_type != missing_lmh ||
              pers.registered_blind ||
              pers.pip_daily_living_type == enhanced_pip )
              n += 1
        end
        if bures.pers[pid].income[CARERS_ALLOWANCE] > 0
            n -= 1
        end
    end
    return max( 0, n )
end

"""
See: https://www.gov.uk/disability-premiums/eligibility
"""
function qualifies_for_disability_premium(
    pers :: Person,
    pres :: IndividualResult,
    prem_sys :: Premia ) :: Bool
    if pers.registered_blind
        return true
    end
    if any_positive( pres.income, prem_sys.disability_premium_qualifying_benefits)
        return true
    end
    return false
end

function qualifies_for_enhanced_disability(
    pers     :: Person,
    pres     :: IndividualResult,
    prem_sys :: Premia,
    nmt      :: NonMeansTestedSys,
    age_limits :: AgeLimits ) :: Bool
    if reached_state_pension_age(
        age_limits,
        pers.age,
        pers.sex )
        return false
    end
    if any_positive( pres.income, [NON_CONTRIB_EMPLOYMENT_AND_SUPPORT_ALLOWANCE,CONTRIB_EMPLOYMENT_AND_SUPPORT_ALLOWANCE])
        return true
    end
    for p1 in [
        PERSONAL_INDEPENDENCE_PAYMENT_DAILY_LIVING,
        SCOTTISH_DISABILITY_ASSISTANCE_WORKING_AGE_DAILY_LIVING]
        if p1 in prem_sys.enhanced_disability_premium_qualifying_benefits
            if pres.income[p1] >= nmt.pip.dl_enhanced
                return true
            end
        end
    end
    for p1 in [DLA_SELF_CARE, SCOTTISH_DISABILITY_ASSISTANCE_CHILDREN_DAILY_LIVING]
        if p1 in prem_sys.enhanced_disability_premium_qualifying_benefits
            if pres.income[p1] >= nmt.dla.care_high
                return true
            end
        end
    end
    ## FIXME support group for ESA
    return false
end

"""
FIXME this should be over the HOUSEHOLD not the benefit unit,
I think
"""
function num_adults_qualifiying_for_disability_premium(
    bu    :: BenefitUnit,
    bures :: BenefitUnitResult,
    prem_sys :: Premia ) :: Integer
    n = 0
    for pid in bu.adults
        if qualifies_for_disability_premium( bu.people[pid], bures.pers[pid], prem_sys )
            n += 1
        end
    end
    return n
end

function qualifies_for_severe_disability_premium(
    pers     :: Person,
    pres     :: IndividualResult,
    prem_sys :: Premia,
    nmt      :: NonMeansTestedSys ) :: Bool
    return qualifies_for_disability_premium( pers, pres, prem_sys ) &&
        any_positive( pres.income,
            [PERSONAL_INDEPENDENCE_PAYMENT_DAILY_LIVING,
            SCOTTISH_DISABILITY_ASSISTANCE_WORKING_AGE_DAILY_LIVING,
            DLA_SELF_CARE,
            SCOTTISH_DISABILITY_ASSISTANCE_CHILDREN_DAILY_LIVING,
            ATTENDANCE_ALLOWANCE,
            SCOTTISH_DISABILITY_ASSISTANCE_OLDER_PEOPLE] )
end

function calc_premia(
    which_ben :: LMTBenefitType,
    bu        :: BenefitUnit,
    bures     :: BenefitUnitResult,
    intermed  :: MTIntermediate,
    prem_sys  :: Premia,
    nmt       :: NonMeansTestedSys,
    age_limits:: AgeLimits ) :: Tuple
    premium = 0.0
    premset = LMTPremiaSet()
    # disabled child premium
    num_ads = num_adults( bu )
    if which_ben in [hb,ctr]
        if intermed.num_disabled_children > 0
            premium += intermed.num_disabled_children*prem_sys.disabled_child
            union!( premset, [disabled_child])
        end
    end
    # carer's premium - all benefits, I think.
    if intermed.num_carers == 1
        premium += prem_sys.carer_single
        union!( premset,[carer_single] )
    elseif intermed.num_carers == 2
        # FIXME what if 1 is a child?
        premium += prem_sys.carer_couple
        union!( premset, [carer_couple] )
    end
    # disability premiums
    ndis = num_adults_qualifiying_for_disability_premium(
        bu, bures, prem_sys )
    if which_ben in [is, jsa, hb, ctr] # FIXME not quite right but perhaps near enough - should be reached pension age? CPAG19 p344
        if ndis == 1
            premium += prem_sys.disability_single
            union!( premset,[disability_single])
        elseif ndis == 2
            premium += prem_sys.disability_couple
            union!( premset,[disability_couple])
        end
    end
    # enhanced disability - also for ESA
    if which_ben in [is, jsa, hb, ctr, esa]
        nenh = 0
        for pid in bu.adults
            nenh += qualifies_for_enhanced_disability(
                    bu.people[pid],
                    bures.pers[pid],
                    prem_sys,
                    nmt,
                    age_limits )
        end
        if nenh == 1
            premium += prem_sys.enhanced_disability_single
            union!( premset, [enhanced_disability_single] )
        elseif nenh == 2
            premium += prem_sys.enhanced_disability_couple
            union!( premset,[enhanced_disability_couple])
        end
    end
    # severe disability premium
    if ndis == num_ads # all adults disabled - check for severe dis
        nsev = 0
        for pid in bu.adults
            nsev += qualifies_for_severe_disability_premium(
                bu.people[pid],
                bures.pers[pid],
                prem_sys,
                nmt
            )
        end
        if nsev == 1
            premium += prem_sys.severe_disability_single
            union!( premset, [severe_disability_single])
        elseif nsev == 2
            premium += prem_sys.severe_disability_couple
            union!( premset, [severe_disability_couple])
        end
    end
    # FIXME this is not the right test for sev disabled child
    if which_ben in [hb,ctr,is,jsa,esa,pc]
        premium += intermed.num_severely_disabled_children*prem_sys.enhanced_disability_child
        if intermed.num_severely_disabled_children > 0
            union!( premset,[enhanced_disability_child] )
        end
    end
    #=
    if which_ben in [ hb, ctr, is, jsa, esa, pc, sc ]
        # CPAG 19/20 - I *think* this is what it means..
        nsd = num_qualifying_for_severe_disability( bu, bures, intermed.num_benefit_units ) > 0
        if nsd == 1
            premium += prem_sys.severe_disability_single
            union!( premset, [severe_disability_single])
        elseif nsd == 2
            premium += prem_sys.severe_disability_couple
            union!( premset, [severe_disability_couple])
        end
    end
    =#
    if which_ben in [ is, jsa, esa ] # this should almost never happen given our routing; cpag p345
        if intermed.someone_pension_age
            premium += prem_sys.pensioner_is
            union!( premset, [pensioner_is] )
        end
    end
    # end
    #
    # we're ignoring support components (p355-) for now.
    #
    return (premium, premset)
end

function calc_allowances(
    which_ben :: LMTBenefitType,
    intermed :: MTIntermediate,
    pas :: PersonalAllowances,
    ages :: AgeLimits,
    oo_housing_costs :: Real  ) :: Real
    pers_allow = 0.0
    @assert intermed.num_adults in [1,2]
    if which_ben == pc
        if intermed.num_adults == 1
            pers_allow = pas.pc_mig_single
        elseif intermed.num_adults == 2
            pers_allow = pas.pc_mig_couple
        end
        # children - cpag p 272 says there's an allowance for children
        # if not claiming CTC, but so far as I can see there's no
        # upper limit on CTC
        # .. anyway, the AC calculator suggests this:
        pers_allow += pas.pc_child * intermed.num_allowed_children
        ## FIXME also an allowance for 1st child born before 6th april 2018
    else
        if intermed.age_oldest_adult < 18
            # argh .. not there's a conditional on ESA that we don't cover here
            # seems to be it - no change for sps, marriage..
            # but some change
            if intermed.num_adults == 2
                pers_allow = pas.couple_both_under_18
            else
                if which_ben != esa
                    pers_allow = pas.age_18_24
                else # FIXME rename this allowance !!
                    pers_allow = pas.age_25_and_over
                end
            end
            # should be somethinh like ...
            #if which_ben in [is,jsa,esa]
            #
            #elseif which_ben == hb
            #    pers_allow = pas.age_18_24
            #end
        elseif intermed.age_youngest_adult < 18 && intermed.num_adults == 2
            pers_allow = pas.couple_one_over_18_high
            ## FIXME some cases lower than this see p 335 CPAG
        else # all over 17
            if intermed.num_adults == 1
                if intermed.num_children > 0 # single parent
                    if intermed.someone_pension_age
                        pers_allow = pas.lone_parent_over_pension_age
                    else
                        pers_allow = pas.lone_parent
                    end
                else
                    if intermed.someone_pension_age
                        pers_allow = pas.over_pension_age
                    elseif intermed.age_oldest_adult < 25
                        pers_allow = pas.age_18_24
                    else
                        pers_allow = pas.age_25_and_over
                    end
                end
            else # 2 adults, both at least 18
                if intermed.someone_pension_age
                    pers_allow = pas.couple_over_pension_age
                else
                    pers_allow = pas.couple_both_over_18
                end
            end # 2 adults
        end # no 16-17 yos
        if which_ben in [hb,ctr] #
            pers_allow += pas.child * intermed.num_allowed_children
        end
    end
    pers_allow += oo_housing_costs
    @assert pers_allow > 0
    return pers_allow
end

function calc_full_ctc(
    bu :: BenefitUnit,
    intermed :: MTIntermediate,
    ctc :: ChildTaxCredit )
    ctc_elements = 0.0
    # p282
    if num_born_before( bu, Date( 2017, 4, 6 )) > 0
        ctc_elements += ctc.family
    end
    ctc_elements += (intermed.num_allowed_children*ctc.child)
    ctc_elements += intermed.num_disabled_children*ctc.disability
    ctc_elements += intermed.num_severely_disabled_children*ctc.severe_disability
    return ctc_elements
end

function calcWTC_CTC!(
    benefit_unit_result :: BenefitUnitResult,
    benefit_unit :: BenefitUnit,
    intermed :: MTIntermediate,
    wtc :: WorkingTaxCredit,
    ctc :: ChildTaxCredit )

    bu = benefit_unit # aliases
    bures = benefit_unit_result
    bu_lmt = bures.legacy_mtbens
    can_apply_for = bu_lmt.can_apply_for
    it, ni = aggregate_tax( bures )

    other_income = it.savings_income + it.dividends_income
    # FIXME does this tread pensions correctly - they're disgregarded
    if other_income < wtc.non_earnings_minima
        other_income = 0.0
    end
    nsi = it.non_savings_income - (total( bures, PENSION_CONTRIBUTIONS_EMPLOYEE )+it.pension_relief_at_source )
    income = other_income + nsi
    wtc_elements = 0.0
    ctc_elements = 0.0
    threshold = wtc.threshold
    cost_of_childcare = 0.0
    if can_apply_for.wtc && ( ! wtc.abolished )
        wtc_elements = wtc.basic
        if intermed.is_sparent
            wtc_elements += wtc.lone_parent
        elseif intermed.num_adults > 1
            wtc_elements += wtc.couple
        end
        if intermed.someone_working_ft > 0
            wtc_elements += wtc.hours_ge_30
        end
        if intermed.working_disabled
            wtc_elements += wtc.disability
        end
        # FIXME do servere_working disabled
        cost_of_childcare = 0.0
        for pid in bu.children
            cost_of_childcare += bu.people[pid].cost_of_childcare
        end

        cost_of_childcare *= wtc.childcare_proportion

        if intermed.num_children > 1
            cost_of_childcare = min( wtc.childcare_max_2_plus_children, cost_of_childcare )
        else
            cost_of_childcare = min( wtc.childcare_max_1_child, cost_of_childcare )
        end
        wtc_elements += cost_of_childcare
    end
    if can_apply_for.ctc && ( ! ctc.abolished )
        ctc_elements = calc_full_ctc( bu, intermed, ctc)
        if ! can_apply_for.wtc
            threshold = ctc.threshold
        end
    end

    elements = ctc_elements + wtc_elements
    excess = wtc.taper * max( 0.0, income - threshold )
    wtc_ctc = max( 0.0, elements - excess )
    # allocate
    ctc_amt = min( wtc_ctc, ctc_elements )
    wtc_amt = wtc_ctc - ctc_amt
    ## assign to an individual

    bu_lmt.wtc_ctc_tapered_excess = excess
    bu_lmt.wtc_income = income
    bu_lmt.ctc_elements = ctc_elements
    bu_lmt.wtc_elements = wtc_elements
    bu_lmt.cost_of_childcare = cost_of_childcare
    bu_lmt.wtc_ctc_threshold = threshold

    # to spouse if has one
    recipient = make_recipient( bu, WORKING_TAX_CREDIT )
    bures.legacy_mtbens.wtc_recipient = recipient
    bures.pers[recipient].income[WORKING_TAX_CREDIT] = wtc_amt

    recipient = make_recipient( bu, CHILD_TAX_CREDIT )
    bures.pers[recipient].income[CHILD_TAX_CREDIT] = ctc_amt
end

function calc_NDDs(
    bu       :: BenefitUnit,
    bur      :: BenefitUnitResult,
    intermed :: MTIntermediate,
    incomes  :: LMTIncomes,
    hb       :: HousingBenefits )::Real
    ndd = 0.0
    # income based on the couples income if it is a couple
    any_pay_ndds = false
    for pid in bu.adults
        pays_ndd = true
        pers = bu.people[pid]
        persr = bur.pers[pid]
        # fixme check this list & reference
        if any_positive( persr.income,
            [ATTENDANCE_ALLOWANCE,
             DLA_SELF_CARE,
             PERSONAL_INDEPENDENCE_PAYMENT_DAILY_LIVING,
             SCOTTISH_DISABILITY_ASSISTANCE_WORKING_AGE_DAILY_LIVING,
             PENSION_CREDIT] )
            pays_ndd = false
        elseif pers.registered_blind || pers.registered_partially_sighted
            pays_ndd = false
        elseif pers.age < 18
            pays_ndd = false
        elseif pers.age < 25
            if any_positive(
                    persr.income,
                    [INCOME_SUPPORT,
                    NON_CONTRIB_JOBSEEKERS_ALLOWANCE,
                    NON_CONTRIB_EMPLOYMENT_AND_SUPPORT_ALLOWANCE] )
                pays_ndd = false
            end
        end
        any_pay_ndds = any_pay_ndds || pays_ndd
    end # pids loop
    if any_pay_ndds
        if intermed.someone_working_ft
            n = size(hb.ndd_incomes)[1]
            w = n
            for i in 1:n
                if hb.ndd_incomes[i] > incomes.gross_earnings
                    w = i
                    break
                end
            end
            ndd = hb.ndd_deductions[w]
        else
            ndd = hb.ndd_deductions[1]
        end
    end # any ndd payable
    return ndd
end

function calculateHB_CTR!(
    household_result :: HouseholdResult,
    which_ben        :: LMTBenefitType,
    hh               :: Household,
    intermed         :: HHIntermed,
    lmt_ben_sys      :: LegacyMeansTestedBenefitSystem,
    age_limits       :: AgeLimits,
    nmt_bens         :: NonMeansTestedSys )
    @argcheck which_ben in [ctr, hb]

    eligible_amount = which_ben == ctr ?
        total(household_result, LOCAL_TAXES ) :
        household_result.housing.allowed_rent

    bus = get_benefit_units( hh )
    nbus = size(bus)[1]
    ndds = 0.0
    for bn in nbus:-1:1 # fixme bn=>buno for consistency
        bu = bus[bn]
        bures = household_result.bus[bn]
        incomes = calc_incomes(
            hb,
            bu,
            bures,
            intermed.buint[bn],
            lmt_ben_sys.income_rules,
            lmt_ben_sys.hours_limits )
        if bn == 1
            benefit = max(0.0, eligible_amount - ndds ) # ndds deducted from eligible rent
            premium = 0.0
            allowances = 0.0
            passported = false
            # FIXME we're doing passpored twice
            if has_any( bures, lmt_ben_sys.hb.passported_bens... )
                # no need to do anything
                passported = true
            else
                ## FIXME pass this in
                premium, premset = calc_premia(
                    hb,
                    bu,
                    bures,
                    intermed.buint[bn],
                    lmt_ben_sys.premia,
                    nmt_bens,
                    age_limits )
                union!(bures.legacy_mtbens.premia, premset)
                allowances = calc_allowances(
                    hb,
                    intermed.buint[bn],
                    lmt_ben_sys.allowances,
                    age_limits,
                    0.0  )
                excess = max( 0.0, incomes.total_income - (premium+allowances))
                if excess > 0
                    taper = which_ben == ctr ? lmt_ben_sys.ctr.taper : lmt_ben_sys.hb.taper
                    benefit = max( 0.0, benefit - taper*excess )
                end

            end

            recipient = make_recipient( bu, HOUSING_BENEFIT )
            bures.legacy_mtbens.hb_recipient = recipient
            bures.legacy_mtbens.ctr_recipient = recipient

            if which_ben == hb
                bures.pers[recipient].income[HOUSING_BENEFIT] = benefit
                bures.legacy_mtbens.hb_passported = passported
                bures.legacy_mtbens.hb_eligible_rent = eligible_amount
                if ! passported
                    bures.legacy_mtbens.hb_premia = premium
                    bures.legacy_mtbens.hb_allowances = allowances
                    bures.legacy_mtbens.hb_incomes = incomes
                end
            elseif which_ben == ctr
                bures.pers[recipient].income[COUNCIL_TAX_BENEFIT] = benefit
                bures.legacy_mtbens.ctr = benefit
                bures.legacy_mtbens.ctr_passported = passported
                bures.legacy_mtbens.ctr_eligible_amount = eligible_amount
                if ! passported
                    bures.legacy_mtbens.ctr_premia = premium
                    bures.legacy_mtbens.ctr_allowances = allowances
                    bures.legacy_mtbens.ctr_incomes = incomes
                end
            end
        else # ndds for hb, not ctr
            if which_ben == hb
                ndds += calc_NDDs(
                    bu,
                    bures,
                    intermed.buint[bn],
                    incomes,
                    lmt_ben_sys.hb )
            end
        end
    end
end # calculateHB_CTR

function calc_legacy_means_tested_benefits!(
    benefit_unit_result :: BenefitUnitResult,
    benefit_unit :: BenefitUnit,
    intermed     :: MTIntermediate,
    mt_ben_sys   :: LegacyMeansTestedBenefitSystem,
    age_limits   :: AgeLimits,
    hours        :: HoursLimits,
    nmt_bens     :: NonMeansTestedSys,
    hh :: Union{Nothing,Household} ) # passed on 1st BU only
    # aliases
    bures = benefit_unit_result
    bu = benefit_unit
    premium = 0.0
    bures.legacy_mtbens.can_apply_for = make_lmt_benefit_applicability(
        mt_ben_sys,
        intermed,
        mt_ben_sys.hours_limits
    )
    # FIXME we need to check the fuck out of this charge
    oo_housing_costs = 0.0
    if hh !== nothing
        if is_owner_occupier( hh.tenure )
            oo_housing_costs = hh.other_housing_charges
        end
    end
    # alias
    can_apply_for = bures.legacy_mtbens.can_apply_for
    # FIXME MIG not really the right name for is,jsa,esa
    which_mig = nothing
    if can_apply_for.esa
        which_mig = esa
    elseif can_apply_for.jsa
        which_mig = jsa
    elseif can_apply_for.is
        which_mig = is
    end
    mig_entitlement = 0.0 # we need to record this for wtc/ctc routing below
    if which_mig !== nothing
        premium,premset = calc_premia(
            which_mig,
            bu,
            bures,
            intermed,
            mt_ben_sys.premia,
            nmt_bens,
            age_limits )
        union!(bures.legacy_mtbens.premia, premset)
        incomes = calc_incomes(
            which_mig,
            bu,
            bures,
            intermed,
            mt_ben_sys.income_rules,
            hours )
        allowances = calc_allowances(
            which_mig,
            intermed,
            mt_ben_sys.allowances,
            age_limits,
            oo_housing_costs )
        bures.legacy_mtbens.mig = premium+allowances
        bures.legacy_mtbens.mig_incomes = incomes
        bures.legacy_mtbens.mig_allowances = allowances
        bures.legacy_mtbens.mig_premia = premium
        mig_entitlement = max( 0.0, bures.legacy_mtbens.mig - incomes.total_income );
        if (! incomes.disqualified_on_capital) && (mig_entitlement > 0)
            # FIXME we just allocate payment to the head of
            # the BU - make this a function or make can_apply_for
            # specify the payee as well as whether the BU qualifies
            if can_apply_for.ctc
                recipient = make_recipient( bu, CHILD_TAX_CREDIT )
                bures.pers[recipient].income[CHILD_TAX_CREDIT] =
                    calc_full_ctc( bu, intermed, mt_ben_sys.child_tax_credit )
            end
            if can_apply_for.esa
                recipient = make_recipient( bu, NON_CONTRIB_EMPLOYMENT_AND_SUPPORT_ALLOWANCE )
                bures.pers[recipient].income[NON_CONTRIB_EMPLOYMENT_AND_SUPPORT_ALLOWANCE] = mig_entitlement
                # bures.legacy_mtbens.esa = entitlement
            elseif can_apply_for.jsa
                recipient = make_recipient( bu, NON_CONTRIB_JOBSEEKERS_ALLOWANCE )
                bures.pers[recipient].income[NON_CONTRIB_JOBSEEKERS_ALLOWANCE] = mig_entitlement
                # bures.legacy_mtbens.jsa = entitlement
            elseif can_apply_for.is
                recipient = make_recipient( bu, INCOME_SUPPORT )
                bures.pers[recipient].income[INCOME_SUPPORT] = mig_entitlement
                # bures.legacy_mtbens.is = entitlement
            end
            bures.legacy_mtbens.mig_recipient = recipient
            can_apply_for.wtc = false # no overlapping
        end
    end

    if can_apply_for.pc && ( ! mt_ben_sys.pen_credit_abolished)
        premium, premset = calc_premia(
            pc,
            bu,
            bures,
            intermed,
            mt_ben_sys.premia,
            nmt_bens,
            age_limits )
        union!( bures.legacy_mtbens.premia, premset )
        incomes = calc_incomes(
            pc,
            bu,
            bures,
            intermed,
            mt_ben_sys.income_rules,
            hours )
        allowances = calc_allowances(
            pc,
            intermed,
            mt_ben_sys.allowances,
            age_limits,
            oo_housing_costs
        )
        # NOTE - there can be 'qualifying housing costs' in the MIG;
        # see: CPAG 2122
        bures.legacy_mtbens.mig = premium+allowances # = max( 0.0, entitlement - incomes.total_income );
        bures.legacy_mtbens.pc_incomes = incomes
        bures.legacy_mtbens.pc_allowances = allowances
        bures.legacy_mtbens.pc_premia = premium
        pc_entitlement = max( 0.0,
            bures.legacy_mtbens.mig -
            bures.legacy_mtbens.pc_incomes.total_income )

        recipient = make_recipient( bu, SAVINGS_CREDIT )

        if can_apply_for.sc && ( ! incomes.disqualified_on_capital ) && ( ! mt_ben_sys.savings_credit.abolished )
            scsys = mt_ben_sys.savings_credit #  shortcut
            sc_incomes = calc_incomes(
                sc,
                bu,
                bures,
                intermed,
                mt_ben_sys.income_rules,
                hours
                 )
            thresh = scsys.threshold_single
            maxpay = scsys.max_single
            if intermed.num_adults > 1
                thresh = scsys.threshold_couple
                maxpay = scsys.max_couple
            end
            # see: CPAG 2011/12 edition ch 19
            scent = 0.0
            income_over_thresh = max(0.0, sc_incomes.total_income - thresh )
            if income_over_thresh > 0
                sc_maximum = min( maxpay, scsys.withdrawal_rate * income_over_thresh)
                income_over_mig = incomes.total_income - bures.legacy_mtbens.mig
                if income_over_mig <= 0.0
                    scent = sc_maximum
                else
                    scent = max( 0.0, sc_maximum - (1-scsys.withdrawal_rate)*income_over_mig)
                end
                bures.pers[recipient].income[SAVINGS_CREDIT] = scent
                bures.legacy_mtbens.sc_incomes = sc_incomes
            end
        end

        bures.pers[recipient].income[PENSION_CREDIT] = pc_entitlement
    end

    #
    # Do WTC iff you can apply for it or if you qualify for CTC
    # but not the automatic full amount
    #
    if can_apply_for.wtc || (can_apply_for.ctc && ( mig_entitlement == 0.0 ))
        calcWTC_CTC!(
                bures,
                bu,
                intermed,
                mt_ben_sys.working_tax_credit,
                mt_ben_sys.child_tax_credit )

    end

    if has_any( bures,
            PENSION_CREDIT,
            NON_CONTRIB_JOBSEEKERS_ALLOWANCE,
            NON_CONTRIB_EMPLOYMENT_AND_SUPPORT_ALLOWANCE,
            INCOME_SUPPORT )
        bures.legacy_mtbens.hb_passported = true
        bures.legacy_mtbens.ctr_passported = true

    end
end

"""
 TODO FIXME housing costs for oo if qualifies for PC,JSA etc.
    see CPAG ch 18
"""
function calc_legacy_means_tested_benefits!(
            household_result :: HouseholdResult,
            household        :: Household,
            intermed         :: HHIntermed,
            lmt_ben_sys      :: LegacyMeansTestedBenefitSystem,
            age_limits       :: AgeLimits,
            hours            :: HoursLimits,
            nmt_bens         :: NonMeansTestedSys,
            hr               :: HousingRestrictions )
    # fixme not just for renters? fixme do this earlier
    household_result.housing = apply_rent_restrictions(
        household, intermed.hhint, hr )
    bus = get_benefit_units(household)
    nbus = size( bus )[1]
    for buno in 1:nbus
        calc_legacy_means_tested_benefits!(
            household_result.bus[buno],
            bus[buno],
            intermed.buint[buno],
            lmt_ben_sys,
            age_limits,
            hours,
            nmt_bens,
            buno == 1 ? household : nothing )
    end
    # hb using the whole hhls but assigned to 1st bu
    if ! lmt_ben_sys.hb.abolished
        calculateHB_CTR!(
            household_result,
            hb,
            household,
            intermed,
            lmt_ben_sys,
            age_limits,
            nmt_bens )
    end

    if ! lmt_ben_sys.ctr.abolished
        calculateHB_CTR!(
            household_result,
            ctr,
            household,
            intermed,
            lmt_ben_sys,
            age_limits,
            nmt_bens )
    end
    #
end

end # module LegacyMeansTestedBenefits

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
"""
    scaled_LHC_sampling_plan(search_range::Array{Tuple{Float64,Float64},1},num_samples,sampling_plan_opt_gens;trace::Symbol=:silent)

Create a Latin Hypercube sampling plan scaled to the `search_range` where
`size(plan) = (num_dimensions,num_samples)`.
"""
function scaled_LHC_sampling_plan(search_range::Array{Tuple{Float64,Float64},1},num_samples,sampling_plan_opt_gens;trace::Symbol=:silent)

    #Printing of creation
    (trace == :verbose) && print_lhc_create()

    unscaled_plan, audze = LHCoptim(num_samples,length(search_range),sampling_plan_opt_gens)

    unscaled_plan = unscaled_plan'

    #Verbose printing of resulting plan
    (trace == :verbose) && print_lhc_plan(audze[end],num_samples,length(search_range),sampling_plan_opt_gens)

    plan = Array{eltype(eltype(search_range))}(undef,size(unscaled_plan))
    for (i,sr) in enumerate(search_range)
        plan[i,:] = scale(unscaled_plan[i,:],sr[1],sr[2])
    end
    return plan
end

function print_lhc_create()
    println("Creating optimized Latin Hypercube Sampling Plan")
end

function print_lhc_plan(audze,num_samples,dims,sampling_plan_opt_gens)
    println("    Optimised LHC plan summary")
    println(@sprintf("    %-15s %-40s %-10s %-0s", "Iterations","Audze–Eglais obj. (higher is better)","Samples","Dimensions"))
    println(@sprintf("    %-15.7g %-40.7g %-10.7g %-0.7g \n", sampling_plan_opt_gens,audze,num_samples,dims))
end


text:
@testset "Storage data misspecification" begin
    # See https://discourse.julialang.org/t/how-to-use-test-warn/15557/5 about testing for warning throwing
    warn_message = "The data doesn't include devices of type GenericBattery, consider changing the device models"
    model = DeviceModel(PSY.GenericBattery, PSI.BookKeeping)
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5)
    @test_logs (:warn, warn_message) construct_device!(op_problem, :Storage, model)
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys14)
    @test_logs (:warn, warn_message) construct_device!(op_problem, :Storage, model)
end

@testset "Storage Basic Storage With DC - PF" begin
    model = DeviceModel(PSY.GenericBattery, PSI.BookKeeping)
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_bat)
    construct_device!(op_problem, :Storage, model)
    moi_tests(op_problem, false, 72, 72, 0, 0, 24, false)
    psi_checkobjfun_test(op_problem, GAEVF)
end

@testset "Storage Basic Storage With AC - PF" begin
    model = DeviceModel(PSY.GenericBattery, PSI.BookKeeping)
    op_problem = OperationsProblem(TestOpProblem, ACPPowerModel, c_sys5_bat)
    construct_device!(op_problem, :Storage, model)
    moi_tests(op_problem, false, 96, 96, 0, 0, 24, false)
    psi_checkobjfun_test(op_problem, GAEVF)
end

@testset "Storage with Reservation DC - PF" begin
    model = DeviceModel(PSY.GenericBattery, PSI.BookKeepingwReservation)
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_bat)
    construct_device!(op_problem, :Storage, model)
    moi_tests(op_problem, false, 96, 24, 48, 48, 24, true)
    psi_checkobjfun_test(op_problem, GAEVF)
end

@testset "Storage with Reservation With AC - PF" begin
    model = DeviceModel(PSY.GenericBattery, PSI.BookKeepingwReservation)
    op_problem = OperationsProblem(TestOpProblem, ACPPowerModel, c_sys5_bat)
    construct_device!(op_problem, :Storage, model)
    moi_tests(op_problem, false, 120, 48, 48, 48, 24, true)
    psi_checkobjfun_test(op_problem, GAEVF)
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/AsterReader.jl/blob/master/LICENSE

using Test, AsterReader

@testset "AsterReader.jl" begin
    include("test_read_aster_mesh.jl")
    include("test_read_aster_results.jl")
    include("test_read_gmsh_med.jl")
    include("test_read_file_not_found.jl")
end


text:
# This script  will be called from runtests.jl as
#   ./argparse.jl 33 9
using ArgParse

parser = ArgParseSettings()
@add_arg_table! parser begin
    "--num1"
        help = "First number"
        arg_type = Int
        default = 33
    "--num2"
        help = "Second number"
        arg_type = Int
        default = 42
    "--ans"
        help = "Answer to check"
        arg_type = Int
        default = 75
end

@testset "ArgParseTest" begin
    parsed_args = parse_args(parser)
    ans = parsed_args["num1"] + parsed_args["num2"]

    @test ans == parsed_args["ans"]
end


text:
using HackerNewsSlackApp
using Test

@testset "HackerNewsSlackApp.jl" begin
    # Write your own tests here.
end


text:
using JuliaDB

using Dates

function get_genotype(filtered_vcf_table)

    start_time = now()

    genotypes = []

    for row in filtered_vcf_table

        genotype = []

        allele_ref = row[:REF]

        allele_alt = row[:ALT]

        numeric_genotype = split(split(row[:GERM], ":")[1], "/")

        if occursin("|", numeric_genotype[1])

            numeric_genotype = split(numeric_genotype[1], "|")

        end

        counter = 1

            while counter < 3

                if "0" == numeric_genotype[counter]

                    push!(genotype, allele_ref)

                elseif "1" == numeric_genotype[counter]

                    push!(genotype, split(allele_alt, ",")[1])

                elseif "2" == numeric_genotype[counter]

                    push!(genotype, split(allele_alt, ",")[2])

                end

            counter += 1

            end

        push!(genotypes, genotype)

    end

    return genotypes

    end_time = now()

    println("\nTook $(canonicalize(Dates.CompoundPeriod(end_time - start_time))).\n")

end

export get_genotype


text:
##################
# Load Libraries #
##################
push!(LOAD_PATH,pwd())
using Printf
using JuliaPBF

##################
# Initialization #
##################
INPUT_DIR   = "./runtests.xml" # Directory xml Inputs
ANALYDATA   = JuliaPBF.IO_dev.Parsing_xml_file(INPUT_DIR) # Parsing Inputs
SIMULDATA   = JuliaPBF.Solver.PreProcessing(ANALYDATA) # FillInitialBox, GenInitGrid

#################
# Update Solver #
#################
JuliaPBF.Solver.Update(SIMULDATA, ANALYDATA)

max_stars_repo_name: Hipgineer/JuliaPBF

max_issues_repo_name: Hipgineer/JuliaPBF

max_forks_repo_name: Hipgineer/JuliaPBF


text:
# Function for loading sample data

"""
    get_data(dataset, poollength, datalength, horizon)

Return features and labels from one of the sample datasets in the repository. `dataset`
can be one of `:solar`, `:traffic`, `:exchange_rate` or `:electricity`. `poollength` gives
the number of timesteps to pool for the model, `datalength` determines the number of time steps
included into the output, and `horizon` determines the number of time steps that should be
forecasted by the model.

See also: [`prepare_data`](@ref), [`load_data`](@ref)
"""
get_data(dataset, poollength, datalength, horizon; normalise=true) =
    prepare_data(load_data(dataset), poollength, datalength, horizon; normalise=normalise)

"""
    prepare_data(data, poollength, datalength, horizon)
    prepare_data(data, poollength, datalength, horizon; normalise=true)

Cast 2D time series data into the format used by FluxArchitectures. `data` is a matrix or
Tables.jl compatible datasource containing data in the form `timesteps x features` (i.e.
each column contains the time series for one feature). `poollength` defines the number of
timesteps to pool when preparing a single frame of data to be fed to the model. `datalength`
determines the number of time steps included into the output, and `horizon` determines the
number of time steps that should be forecasted by the model. Outputs features and labels.
"""
function prepare_data(data::AbstractMatrix, poollength, datalength, horizon; normalise=true)
    extendedlength = datalength + poollength - 1
    extendedlength > size(data, 1) && throw(ArgumentError("datalength $(datalength) larger than available data $(size(data, 1) - poollength + 1)"))
    (normalise == true) && (data = Flux.normalise(data, dims=1))
    features = similar(data, size(data, 2), poollength, 1, datalength)
    for i in 1:datalength
      features[:,:,:,i] .= permutedims(data[i:(i + poollength - 1) ,:])
    end
    labels = circshift(data[1:datalength,1], -horizon)
    return features, labels
end

prepare_data(data, poollength, datalength, horizon; normalise=true) = prepare_data(Tables.matrix(data), poollength, datalength, horizon; normalise=true)

"""
    load_data(dataset)

Load the raw data from one of the available datasets.
The following example data from https://github.com/laiguokun/multivariate-time-series-data
is included:

* `:solar`: The raw data is coming from http://www.nrel.gov/grid/solar-power-data.html:
  It contains the solar power production records in the year of 2006, which is sampled every
  10 minutes from 137 PV plants in Alabama State.

* `:traffic`: The raw data is coming from http://pems.dot.ca.gov. The data in this repo is a
  collection of 48 months (2015-2016) hourly data from the California Department of
  Transportation. The data describes the road occupancy rates (between 0 and 1) measured by
  different sensors on San Francisco Bay area freeways.

* `:electricity`: The raw dataset is from https://archive.ics.uci.edu/ml/datasets/ElectricityLoadDiagrams20112014.
  It is the electricity consumption in kWh was recorded every 15 minutes from 2011 to 2014
  for 321 clients. The data has been cleaned and converted to hourly consumption.

* `:exchange_rate`: The collection of daily exchange rates of eight foreign countries
  including Australia, Great Britain, Canada, Switzerland, China, Japan, New Zealand and
  Singapore ranging from 1990 to 2016.
"""
function load_data(dataset)
    admissible = [:solar, :traffic, :exchange_rate, :electricity]
    dataset in admissible || throw(ArgumentError("Sample data $dataset not found"))

    datadir = joinpath(artifact"sample_data", "data")
    (dataset == :solar) && (BSON.@load joinpath(datadir, "solar_AL.bson") inp_raw)
    (dataset == :traffic) && (BSON.@load joinpath(datadir, "traffic.bson") inp_raw)
    (dataset == :exchange_rate) && (BSON.@load joinpath(datadir, "exchange_rate.bson") inp_raw)
    (dataset == :electricity) && (BSON.@load joinpath(datadir, "electricity.bson") inp_raw)
    return inp_raw
end

max_stars_repo_name: sdobber/FluxArchitectures

max_issues_repo_name: sdobber/FluxArchitectures

max_forks_repo_name: sdobber/FluxArchitectures


text:
using KernelAbstractions, Test
using CUDA

if has_cuda_gpu()
    CUDA.allowscalar(false)
end

@kernel function kernel_print()
  I = @index(Global)
  @print("Hello from thread ", I, "!\n")
end

function test_print(backend)
  kernel = kernel_print(backend, 4)
  kernel(ndrange=(4,))
end

@testset "print test" begin
    if has_cuda_gpu()
        wait(test_print(CUDADevice()))
        @test true
    end

    wait(test_print(CPU()))
    @test true

    @print("Why this should work")
    @test true
end