text:
## ExoplanetsSysSim/examples/dr25_gaia_fgk/abc_setup.jl
## (c) 2019 Eric B. Ford & Danley C. Hsu
# Collection of functions which specific ABC simulation parameters

module EvalSysSimModel
export setup, get_param_vector, get_ss_obs
export gen_data, calc_summary_stats, calc_distance, is_valid_uniform, is_valid_beta, is_valid_dirichlet, normalize_dirch
using ExoplanetsSysSim
using ApproximateBayesianComputing
const ABC = ApproximateBayesianComputing
include(joinpath(abspath(joinpath(dirname(Base.find_package("ExoplanetsSysSim")),"..")),"examples","dr25_gaia_m", "dr25_binrates_func.jl"))

sim_param_closure = SimParam()
summary_stat_ref_closure =  CatalogSummaryStatistics()

function is_valid_uniform(param_vector::Vector{Float64})
    global sim_param_closure
    update_sim_param_from_vector!(param_vector,sim_param_closure)
    local rate_tab::Array{Float64,2} = get_any(sim_param_closure, "obs_par", Array{Float64,2})
    limitP::Array{Float64,1} = get_any(sim_param_closure, "p_lim_arr", Array{Float64,1})
    #const lambda = sum_kbn(rate_tab)
    if any(x -> x < 0., rate_tab) || any([floor(3*log(limitP[i+1]/limitP[i])/log(2)) for i in 1:length(limitP)-1] .< sum(rate_tab, dims=1)')
        return false
    end
    return true
end

function is_valid_dirichlet(param_vector::Vector{Float64})
    global sim_param_closure
    update_sim_param_from_vector!(param_vector,sim_param_closure)
    local rate_tab::Array{Float64,2} = get_any(sim_param_closure, "obs_par", Array{Float64,2})
    limitP::Array{Float64,1} = get_any(EvalSysSimModel.sim_param_closure, "p_lim_arr", Array{Float64,1})
    #const lambda = sum_kbn(rate_tab)
    if any(x -> x < 0., rate_tab) || any([floor(3*log(limitP[i+1]/limitP[i])/log(2)) for i in 1:length(limitP)-1] .< rate_tab[1,:])
        return false
    end
    return true
end

function normalize_dirch(param_vector::Vector{Float64})
    global sim_param_closure
    local p_dim = length(get_any(sim_param_closure, "p_lim_arr", Array{Float64,1}))-1
    local r_dim = length(get_any(sim_param_closure, "r_lim_arr", Array{Float64,1}))-1

    for i in 1:p_dim
        param_vector[((i-1)*(r_dim+1)+2):((i-1)*(r_dim+1)+(r_dim+1))] ./= sum(param_vector[((i-1)*(r_dim+1)+2):((i-1)*(r_dim+1)+(r_dim+1))])
    end

    update_sim_param_from_vector!(param_vector,sim_param_closure)
    return param_vector
end

function is_valid_mfgk_ratio(param_vector::Vector{Float64})
    global sim_param_closure
    update_sim_param_from_vector!(param_vector,sim_param_closure)
    limitP::Array{Float64,1} = get_any(sim_param_closure, "p_lim_arr", Array{Float64,1})
    const rate_tab::Array{Float64,2} = get_any(sim_param_closure, "obs_par", Array{Float64,2})[:,1:length(limitP)]
    const mfgk_ratio::Float64 = get_real(sim_param, "mfgk_ratio")
    #const lambda = sum(rate_tab)
    if any(x -> x < 0., rate_tab) || any([floor(3*log(limitP[i+1]/limitP[i])/log(2)) for i in 1:length(limitP)-1] .< mfgk_ratio*sum(rate_tab, dims=1)')
        return false
    end
    return true
end

function gen_data(param_vector::Vector{Float64})
    global sim_param_closure
    update_sim_param_from_vector!(param_vector,sim_param_closure)
    cat_phys = generate_kepler_physical_catalog(sim_param_closure)
    #cat_phys_cut = ExoplanetsSysSim.generate_obs_targets(cat_phys, sim_param_closure)
    #cat_obs = ExoplanetsSysSim.observe_kepler_targets_single_obs(cat_phys_cut, sim_param_closure)
    cat_obs = ExoplanetsSysSim.observe_kepler_targets_sky_avg(cat_phys, sim_param_closure)
    return cat_obs
end

# TODO OPT: Eventually, could adapt ABC.jl to use distance from first pass to decide if should compute additional summary statistics
function calc_summary_stats(cat::KeplerObsCatalog)
    global sim_param_closure
    sum_stat = calc_summary_stats_obs_binned_rates(cat, sim_param_closure, obs_skyavg = true)
    return sum_stat
end

function calc_distance(sum_stat_obs::CatalogSummaryStatistics,sum_stat_sim::CatalogSummaryStatistics, n::Integer = 0)
    global sim_param_closure
    dist1 = calc_distance_vector_binned(sum_stat_obs,sum_stat_sim, 1, sim_param_closure)
    num_available = length(dist1)
    num_to_use = n>0 ? min(n,num_available) : num_available
    return calc_scalar_distance(dist1[1:num_to_use])
end

function setup(prior_choice::String, bin_size_factor::Float64)
    global sim_param_closure = setup_sim_param_dr25binrates()
    add_param_fixed(sim_param_closure,"bin_size_factor",bin_size_factor)
    if prior_choice == "dirichlet"
        sim_param_closure = set_test_param_total(sim_param_closure)
        add_param_fixed(sim_param_closure,"generate_num_planets",generate_num_planets_binrates_dirichlet)
        if (length(get_any(sim_param_closure, "r_lim_arr", Array{Float64,1}))-1) > 1
            add_param_fixed(sim_param_closure,"generate_period_and_sizes", generate_period_and_sizes_binrates_dirichlet)
        end
    elseif prior_choice == "beta"
        sim_param_closure = set_test_param(sim_param_closure)
        add_param_fixed(sim_param_closure,"generate_num_planets",generate_num_planets_binrates_beta)
        add_param_fixed(sim_param_closure,"generate_period_and_sizes", generate_period_and_sizes_binrates_beta)
    elseif prior_choice == "uniform"
        sim_param_closure = set_test_param(sim_param_closure)
    else
        println("# Invalid prior given!")
        quit()
    end

    ### Use simulated planet candidate catalog data
    # df_star = setup_star_table_christiansen(sim_param_closure)
    # println("# Finished reading in stellar data")
    # add_param_fixed(sim_param_closure,"num_kepler_targets",1000000)  # For "observed" catalog
    # cat_obs = simulated_read_kepler_observations(sim_param_closure)
    # println("# Finished setting up simulated true catalog")
    ###

    ### Use real planet candidate catalog data
    df_koi,usable_koi = read_koi_catalog(sim_param_closure)
    println("# Finished reading in KOI data")
    df_star = setup_star_table_dr25(sim_param_closure)
    println("# Finished reading in stellar data")
    cat_obs = setup_actual_pc_catalog_kepler(df_star, df_koi, usable_koi, sim_param_closure)
    println("# Finished setting up true catalog")
    ###

    global summary_stat_ref_closure = calc_summary_stats_obs_binned_rates(cat_obs,sim_param_closure, trueobs_cat = true)
end

get_param_vector() = make_vector_of_sim_param(sim_param_closure)
get_ss_obs() = summary_stat_ref_closure

function set_simparam_ss(sim_param::ExoplanetsSysSim.SimParam, ss_true::ExoplanetsSysSim.CatalogSummaryStatistics)
    global sim_param_closure = sim_param
    global summary_stat_ref_closure = ss_true
end

end  # module EvalSysSimModel

#include(joinpath(Pkg.dir("ABC"),"src/composite.jl"))

module SysSimABC
export setup_abc, run_abc, run_abc_largegen, setup_abc_p2
using Distributions, Random, Distributed
using ApproximateBayesianComputing
const ABC = ApproximateBayesianComputing
import ApproximateBayesianComputing.CompositeDistributions.CompositeDist
import ApproximateBayesianComputing.TransformedBetaDistributions.LinearTransformedBeta

#using Compat
import ExoplanetsSysSim
import ..EvalSysSimModel
include(joinpath(abspath(joinpath(dirname(Base.find_package("ExoplanetsSysSim")),"..")),"examples","dr25_gaia_m", "dr25_binrates_func.jl"))
include(joinpath(abspath(joinpath(dirname(Base.find_package("ExoplanetsSysSim")),"..")),"examples","dr25_gaia_m", "beta_proposal.jl"))

function setup_abc(num_dist::Integer = 0; prior_choice::String = "uniform", bin_size_factor::Float64 = 1.5)
    EvalSysSimModel.setup(prior_choice, bin_size_factor)

    limitP::Array{Float64,1} = get_any(EvalSysSimModel.sim_param_closure, "p_lim_arr", Array{Float64,1})
    limitR::Array{Float64,1} = get_any(EvalSysSimModel.sim_param_closure, "r_lim_arr", Array{Float64,1})
    limitR_full::Array{Float64,1} = get_any(EvalSysSimModel.sim_param_closure, "r_lim_full", Array{Float64,1})
    local r_dim = length(limitR)-1

    prior_arr = ContinuousDistribution[]
    ss_obs_table = EvalSysSimModel.get_ss_obs().stat["planets table"]

    if prior_choice == "ratio"
        prior_arr = vcat(prior_arr, Uniform(0.0, 10.0))
    elseif prior_choice == "dirichlet"
        weights_arr = [log(limitR[j+1]/limitR[j]) for j in 1:r_dim]/minimum([log(limitR_full[k+1]/limitR_full[k]) for k in 1:(length(limitR_full)-1)])
        for i in 1:(length(limitP)-1)
            max_in_col = 3*log(limitP[i+1]/limitP[i])/log(2)
            lambda_col = Uniform(0.0, max_in_col)
            prior_arr = vcat(prior_arr, lambda_col)
            if r_dim > 1
                dirch_dist = Dirichlet(weights_arr)
                prior_arr = vcat(prior_arr, dirch_dist)
            end
        end
    else
        for i in 1:(length(limitP)-1)
            max_in_col = bin_size_factor*log(limitP[i+1]/limitP[i])/log(2)
            for j in 1:r_dim
                uniform_dist = Uniform(0.0, max_in_col*log(limitR[j+1]/limitR[j])/log(2))
                prior_arr = vcat(prior_arr, uniform_dist)
            end
        end
    end

    param_prior = CompositeDist(prior_arr)
    in_parallel = nworkers() > 1 ? true : false

    calc_distance_ltd(sum_stat_obs::ExoplanetsSysSim.CatalogSummaryStatistics,sum_stat_sim::ExoplanetsSysSim.CatalogSummaryStatistics) = EvalSysSimModel.calc_distance(sum_stat_obs,sum_stat_sim,num_dist)

    global abc_plan = ABC.abc_pmc_plan_type(EvalSysSimModel.gen_data,EvalSysSimModel.calc_summary_stats, calc_distance_ltd, param_prior, make_proposal_dist=make_proposal_dist_multidim_beta, is_valid=EvalSysSimModel.is_valid_uniform, num_part=500, num_max_attempt=200, num_max_times=200, epsilon_init=9.9e99, target_epsilon=1.0e-100, in_parallel=in_parallel, adaptive_quantiles = false, epsilon_reduction_factor=0.9, tau_factor=2.0);

    if prior_choice == "ratio"
        abc_plan.make_proposal_dist = make_proposal_dist_multidim_beta_ratio
        abc_plan.is_valid = EvalSysSimModel.is_valid_mfgk_ratio
    elseif prior_choice == "dirichlet" && r_dim > 1
        abc_plan.make_proposal_dist = make_proposal_dist_multidim_beta_dirichlet
        abc_plan.is_valid = EvalSysSimModel.is_valid_dirichlet
        abc_plan.normalize = EvalSysSimModel.normalize_dirch
    end

    return abc_plan
end

# function setup_abc_p2(abc_plan::ABC.abc_pmc_plan_type)
#     abc_plan.tau_factor = 2.0
#     abc_plan.num_max_times = 200
#     ExoplanetsSysSim.add_param_fixed(EvalSysSimModel.sim_param_closure,"num_targets_sim_pass_one",10000)
#     return abc_plan
# end

function run_abc_largegen(abc_plan::ABC.abc_pmc_plan_type, pop::ABC.abc_population_type, ss_true::ExoplanetsSysSim.CatalogSummaryStatistics, epshist_targ::Float64; npart::Integer = 1000, num_dist::Integer = 0)

    abc_plan.num_max_times = 1

    println("# run_abc_largegen: ",EvalSysSimModel.sim_param_closure)
    sampler_largegen = abc_plan.make_proposal_dist(pop, abc_plan.tau_factor)
    theta_largegen = Array{Float64}(undef, size(pop.theta, 1), npart)
    weight_largegen = Array{Float64}(undef, npart)
    for i in 1:npart
        theta_val, dist_largegen, attempts_largegen = ABC.generate_theta(abc_plan, sampler_largegen, ss_true, epshist_targ)
        theta_largegen[:,i] = theta_val
        prior_logpdf = Distributions.logpdf(abc_plan.prior,theta_val)
        sampler_logpdf = Distributions.logpdf(sampler_largegen, theta_val)
        weight_largegen[i] = exp(prior_logpdf-sampler_logpdf)
    end
    return theta_largegen, weight_largegen
end

function run_abc(abc_plan::ABC.abc_pmc_plan_type)
    #global sim_param_closure
    println("# run_abc: ",EvalSysSimModel.sim_param_closure)
    ss_true = EvalSysSimModel.get_ss_obs()
    #println("True catalog SS: ", ss_true)
    pop_out = ABC.run_abc(abc_plan,ss_true;verbose=true)
end

function run_abc(abc_plan::ABC.abc_pmc_plan_type, pop::ABC.abc_population_type)
    #global sim_param_closure
    dist_threshold = maximum(pop.dist)
    EvalSysSimModel.add_param_fixed(EvalSysSimModel.sim_param_closure,"minimum ABC dist skip pass 2",dist_threshold)
    println("# run_abc: ",EvalSysSimModel.sim_param_closure)
    ss_true = EvalSysSimModel.get_ss_obs()
    pop_out = ABC.run_abc(abc_plan,ss_true,pop;verbose=true)
end

end # module SysSimABC


text:
"""
    ParticleFlowMethod

Abstract type for a method used to flow particles according to a prescribed change of the density.
"""
abstract type ParticleFlowMethod end

"""
    flow!(ensemble::UnweightedParticleEnsemble, method::ParticleFlowMethod) --> eq

Performs a particle flow of `ensemble` according to particle flow method `method`.
"""
function flow!(ensemble::UnweightedParticleRepresentation, method::ParticleFlowMethod) end

ext: jl

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
# Tokenize
# ========

# Token type
# ----------

# token kind
const STRING  = 0b0000
const INTEGER = 0b0001
const FLOAT   = 0b0010
const BOOL    = 0b0100
const QSTRING = 0b1000  # string with quotation marks
const MISSING = 0b1111  # missing can be any data type

const MAX_TOKEN_START = 2^36
const MAX_TOKEN_LENGTH = 2^24

struct Token
    # From most significant
    #    4bit: kind (+ missing)
    #   36bit: start position (64 GiB)
    #   24bit: length (16 MiB)
    value::UInt64

    function Token(kind::UInt8, start::Int, len::Int)
        @assert start < MAX_TOKEN_START
        @assert len < MAX_TOKEN_LENGTH
        return new((UInt64(kind) << 60) | (UInt64(start) << 24) | UInt64(len))
    end
end

const TOKEN_NULL = Token(0x00, 0, 0)

function kind(token::Token)
    return (token.value >> 60) % UInt8
end

function ismissing(token::Token)
    return (token.value & (UInt64(1) << 63)) != 0
end

function location(token::Token)
    x = token.value & (~UInt64(0) >> 4)
    return (x >> 24) % Int, (x & (~UInt64(0) >> 40)) % Int
end

# Line scanner
# ------------

struct ReadError <: Exception
    msg::String
end

macro state(name, ex)
    @assert name isa Symbol
    @assert ex isa Expr && ex.head == :block
    quote
        @label $(name)
        #println($(QuoteNode(name)))
        #@show quoted
        pos += 1
        #if pos > pos_end
        #    @goto END
        #end
        @inbounds c = mem[pos]
        #@show Char(c)
        #println()
        $(ex)
        # This area is unreachable because transition must be exhaustive.
        @assert false
    end |> esc
end

macro multibytestring()
    quote
        # multibyte UTF8 character
        if 0b110_00000 ≤ c ≤ 0b110_11111
            if pos + 1 ≤ pos_end && (mem[pos+1] >> 6) ≤ 0b10  # same as: 0b10_000000 ≤ mem[pos+1] ≤ 0b10_111111
                pos += 1
                @goto STRING
            end
        elseif 0b1110_0000 ≤ c ≤ 0b1110_1111
            if pos + 2 ≤ pos_end && (max(mem[pos+1], mem[pos+2]) >> 6) ≤ 0b10
                pos += 2
                @goto STRING
            end
        elseif 0b11110_000 ≤ c ≤ 0b11110_111
            if pos + 3 ≤ pos_end && (max(mem[pos+1], mem[pos+2], mem[pos+3]) >> 6) ≤ 0b10
                pos += 3
                @goto STRING
            end
        else
            @goto ERROR
        end
    end |> esc
end

macro follows(s)
    @assert s isa String && isascii(s)
    i = 0
    foldl(s, init = :(pos + $(sizeof(s)) ≤ pos_end)) do ex, c
        up = UInt8(uppercase(c))
        lo = UInt8(lowercase(c))
        i += 1
        :($(ex) && (mem[pos+$(i)] == $(up) || mem[pos+$(i)] == $(lo)))
    end |> esc
end

macro begintoken()
    esc(:(start = pos))
end

macro recordtoken(kind)
    quote
        token = Token($(kind), start, pos - start)
    end |> esc
end

macro endheadertoken()
    quote
        push!(tokens, token)
        quoted = false
        qstring = false
    end |> esc
end

function scanheader(mem::Memory, pos::Int, nl::Int, delim::UInt8, quot::UInt8, trim::Bool)
    tokens = Token[]
    token = TOKEN_NULL
    quoted = false
    qstring = false
    start = 0
    pos_end = nl

    @state BEGIN begin
        @begintoken
        if c == quot
            if quoted
                @recordtoken STRING
                @goto QUOTE_END
            end
            quoted = true
            @goto BEGIN
        elseif c == delim
            if quoted
                @goto STRING
            end
            @recordtoken STRING
            @endheadertoken
            @goto BEGIN
        elseif c == SP
            if trim && !quoted
                @goto BEGIN
            end
            @goto STRING
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == LF
            @recordtoken STRING
            @goto LF
        elseif c == CR
            @recordtoken STRING
            @goto CR
        else
            @multibytestring
        end
    end

    @state STRING begin
        if quoted && c == quot
            if qstring
                @recordtoken QSTRING
            else
                @recordtoken STRING
            end
            @goto QUOTE_END
        elseif c == delim
            if quoted
                @goto STRING
            end
            if qstring
                @recordtoken QSTRING
            else
                @recordtoken STRING
            end
            @endheadertoken
            @goto BEGIN
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == SP
            if trim && !quoted
                if qstring
                    @recordtoken QSTRING
                else
                    @recordtoken STRING
                end
                @goto STRING_SPACE
            end
            @goto STRING
        elseif c == LF
            if qstring
                @recordtoken QSTRING
            else
                @recordtoken STRING
            end
            @recordtoken STRING
            @goto LF
        elseif c == CR
            @recordtoken STRING
            @goto CR
        else
            @multibytestring
        end
    end

    @state STRING_SPACE begin
        if c == SP
            @goto STRING_SPACE
        elseif c == delim
            @endheadertoken
            @goto BEGIN
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == LF
            @goto LF
        elseif c == CR
            @goto CR
        else
            @multibytestring
        end
    end

    @state QUOTE_END begin
        if c == delim
            @endheadertoken
            @goto BEGIN
        elseif c == quot
            qstring = true
            @goto STRING
        elseif c == SP
            if trim
                @goto QUOTE_END_SPACE
            end
            @goto ERROR
        elseif c == LF
            quoted = false
            @goto LF
        elseif c == CR
            quoted = false
            @goto CR
        else
            @goto ERROR
        end
    end

    @state QUOTE_END_SPACE begin
        if c == delim
            @endheadertoken
            @goto BEGIN
        elseif c == SP
            @goto QUOTE_END_SPACE
        elseif c == LF
            quoted = false
            @goto LF
        elseif c == CR
            quoted = false
            @goto CR
        else
            @goto ERROR
        end
    end

    @label ERROR
    throw(ReadError("invalid file header format"))

    @label CR
    if quoted
        throw(ReadError("quoted multiline string is not allowed in the header"))
    elseif pos + 1 ≤ pos_end && mem[pos + 1] == LF
        pos += 1
    end

    @label LF
    if quoted
        throw(ReadError("quoted multiline string is not allowed in the header"))
    end
    @endheadertoken

    @label END
    return pos, tokens
end

macro endtoken()
    quote
        if i > ncols
            msg = "unexpected number of columns at line $(line)"
            @goto ERROR
        end
        @inbounds tokens[i,row] = token
        i += 1
        quoted = false
        qstring = false
    end |> esc
end

# Scan a line in mem; mem must include one or more lines.
function scanline!(
        # output info
        tokens::Matrix{Token}, row::Int,
        # input info
        mem::Memory, pos::Int, lastnl::Int, line::Int,
        # parser parameters
        delim::UInt8, quot::UInt8, trim::Bool
    )

    # Check parameters.
    @assert delim != quot
    @assert !trim || delim != SP
    @assert !trim || quot != SP
    @assert mem[lastnl] == CR || mem[lastnl] == LF

    # Initialize variables.
    pos_end = lastnl
    ncols = size(tokens, 1)
    quoted = false
    qstring = false
    token = TOKEN_NULL
    msg = ""  # error message
    start = 0  # the starting position of a token
    i = 1  # the column of a token

    @state BEGIN begin
        @begintoken
        if c == quot
            if quoted
                @recordtoken MISSING
                @goto QUOTE_END
            else
                quoted = true
                @goto BEGIN
            end
        elseif c == delim
            if quoted
                @goto STRING
            end
            @recordtoken MISSING
            @endtoken
            @goto BEGIN
        elseif c == UInt8('-') || c == UInt8('+')
            @goto SIGN
        elseif UInt8('0') ≤ c ≤ UInt8('9')
            @goto INTEGER
        elseif c == SP
            if trim && !quoted
                @goto BEGIN
            end
            @goto STRING
        elseif c == UInt8('.')
            @goto DOT
        elseif c == UInt8('I') || c == UInt8('i')
            # Infinity
            if @follows("NFINITY")
                pos += 7
                @goto SPECIAL_FLOAT
            elseif @follows("NF")
                pos += 2
                @goto SPECIAL_FLOAT
            end
            @goto STRING
        elseif (c == UInt8('N') || c == UInt8('n')) && @follows("AN")  # case-insensitive
            # NaN
            pos += 2  # for 'A' and 'N'
            @goto SPECIAL_FLOAT
        elseif c == UInt8('F') || c == UInt8('f')
            if @follows("alse")
                pos += 4
                @goto BOOL
            end
            @goto BOOL
        elseif c == UInt8('T') || c == UInt8('t')
            if @follows("rue")
                pos += 3
                @goto BOOL
            end
            @goto BOOL
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == LF
            if i == ncols
                @recordtoken MISSING
            end
            @goto LF
        elseif c == CR
            if i == ncols
                @recordtoken MISSING
            end
            @goto CR
        else
            @multibytestring
        end
    end

    @state SIGN begin
        if quoted && c == quot
            @recordtoken STRING
            @goto QUOTE_END
        elseif c == delim
            if quoted
                @goto STRING
            end
            @recordtoken STRING
            @endtoken
            @goto BEGIN
        elseif UInt8('0') ≤ c ≤ UInt8('9')
            @goto INTEGER
        elseif c == UInt8('.')
            @goto DOT
        elseif c == SP
            if trim && !quoted
                @recordtoken STRING
                @goto STRING_SPACE
            end
            @goto STRING
        elseif (c == UInt8('N') || c == UInt8('n')) && @follows("AN")
            pos += 2
            @goto SPECIAL_FLOAT
        elseif c == UInt8('I') || c == UInt8('i')
            # Infinity
            if @follows("NFINITY")
                pos += 7
                @goto SPECIAL_FLOAT
            elseif @follows("NF")
                pos += 2
                @goto SPECIAL_FLOAT
            end
            @goto STRING
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == LF
            @recordtoken STRING
            @goto LF
        elseif c == CR
            @recordtoken STRING
            @goto CR
        else
            @multibytestring
        end
    end

    @state INTEGER begin
        if quoted && c == quot
            @recordtoken INTEGER|FLOAT
            @goto QUOTE_END
        elseif c == delim
            if quoted
                @goto STRING
            end
            @recordtoken INTEGER|FLOAT
            @endtoken
            @goto BEGIN
        elseif UInt8('0') ≤ c ≤ UInt8('9')
            @goto INTEGER
        elseif c == UInt8('.')
            @goto POINT_FLOAT
        elseif c == SP
            if trim && !quoted
                @recordtoken INTEGER|FLOAT
                @goto INTEGER_SPACE
            end
            @goto STRING
        elseif c == UInt8('e') || c == UInt8('E')
            @goto EXPONENT
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == LF
            @recordtoken INTEGER|FLOAT
            @goto LF
        elseif c == CR
            @recordtoken INTEGER|FLOAT
            @goto CR
        else
            @multibytestring
        end
    end

    @state INTEGER_SPACE begin
        if c == SP
            @goto INTEGER_SPACE
        elseif c == delim
            @endtoken
            @goto BEGIN
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == LF
            @goto LF
        elseif c == CR
            @goto CR
        else
            @multibytestring
        end
    end

    @state DOT begin
        if quoted && c == quot
            @recordtoken STRING
            @goto QUOTE_END
        elseif c == delim
            if quoted
                @goto STRING
            end
            @recordtoken STRING
            @endtoken
            @goto BEGIN
        elseif UInt8('0') ≤ c ≤ UInt8('9')
            @goto POINT_FLOAT
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == SP
            if trim && !quoted
                @recordtoken STRING
                @goto STRING_SPACE
            end
            @goto STRING
        elseif c == LF
            @recordtoken STRING
            @goto LF
        elseif c == CR
            @recordtoken STRING
            @goto CR
        else
            @multibytestring
        end
    end

    @state POINT_FLOAT begin
        if quoted && c == quot
            @recordtoken FLOAT
            @goto QUOTE_END
        elseif c == delim
            if quoted
                @goto STRING
            end
            @recordtoken FLOAT
            @endtoken
            @goto BEGIN
        elseif UInt8('0') ≤ c ≤ UInt8('9')
            @goto POINT_FLOAT
        elseif c == UInt8('e') || c == UInt8('E')
            @goto EXPONENT
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == SP
            if trim && !quoted
                @recordtoken FLOAT
                @goto FLOAT_SPACE
            end
            @goto STRING
        elseif c == LF
            @recordtoken FLOAT
            @goto LF
        elseif c == CR
            @recordtoken FLOAT
            @goto CR
        else
            @multibytestring
        end
    end

    @state EXPONENT begin
        if quoted && c == quot
            @recordtoken STRING
            @goto QUOTE_END
        elseif c == delim
            if quoted
                @goto STRING
            end
            @recordtoken STRING
            @endtoken
            @goto BEGIN
        elseif UInt8('0') ≤ c ≤ UInt8('9')
            @goto EXPONENT_FLOAT
        elseif c == UInt8('-') || c == UInt8('+')
            @goto EXPONENT_SIGN
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == SP
            if trim && !quoted
                @recordtoken STRING
                @goto STRING_SPACE
            end
            @goto STRING
        elseif c == LF
            @recordtoken STRING
            @goto LF
        elseif c == CR
            @recordtoken STRING
            @goto CR
        else
            @multibytestring
        end
    end

    @state EXPONENT_SIGN begin
        if quoted && c == quot
            @recordtoken STRING
            @goto QUOTE_END
        elseif c == delim
            if quoted
                @goto STRING
            end
            @recordtoken STRING
            @endtoken
            @goto BEGIN
        elseif UInt8('0') ≤ c ≤ UInt8('9')
            @goto EXPONENT_FLOAT
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == SP
            if trim && !quoted
                @recordtoken STRING
                @goto STRING_SPACE
            end
            @goto STRING
        elseif c == LF
            @recordtoken STRING
            @goto LF
        elseif c == CR
            @recordtoken STRING
            @goto CR
        else
            @multibytestring
        end
    end

    @state EXPONENT_FLOAT begin
        if quoted && c == quot
            @recordtoken FLOAT
            @goto QUOTE_END
        elseif c == delim
            if quoted
                @goto STRING
            end
            @recordtoken FLOAT
            @endtoken
            @goto BEGIN
        elseif UInt8('0') ≤ c ≤ UInt8('9')
            @goto EXPONENT_FLOAT
        elseif c == SP
            if trim && !quoted
                @recordtoken FLOAT
                @goto FLOAT_SPACE
            end
            @goto STRING
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == LF
            @recordtoken FLOAT
            @goto LF
        elseif c == CR
            @recordtoken FLOAT
            @goto CR
        else
            @multibytestring
        end
    end

    @state SPECIAL_FLOAT begin
        if quoted && c == quot
            @recordtoken FLOAT|STRING
            @goto QUOTE_END
        elseif c == delim
            if quoted
                @goto STRING
            end
            @recordtoken FLOAT|STRING
            @endtoken
            @goto BEGIN
        elseif c == SP
            if trim && !quoted
                @recordtoken FLOAT|STRING
                @goto FLOAT_SPACE
            end
            @goto STRING
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == LF
            @recordtoken FLOAT|STRING
            @goto LF
        elseif c == CR
            @recordtoken FLOAT|STRING
            @goto CR
        else
            @multibytestring
        end
    end

    @state FLOAT_SPACE begin
        if c == SP
            @goto FLOAT_SPACE
        elseif c == delim
            @endtoken
            @goto BEGIN
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == LF
            @goto LF
        elseif c == CR
            @goto CR
        else
            @multibytestring
        end
    end

    @state BOOL begin
        if quoted && c == quot
            @recordtoken BOOL
            @goto QUOTE_END
        elseif c == delim
            if quoted
                @goto STRING
            end
            @recordtoken BOOL
            @endtoken
            @goto BEGIN
        elseif c == SP
            if trim && !quoted
                @recordtoken BOOL
                @goto BOOL_SPACE
            end
            @goto STRING
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == LF
            @recordtoken BOOL
            @goto LF
        elseif c == CR
            @recordtoken BOOL
            @goto CR
        else
            @multibytestring
        end
    end

    @state BOOL_SPACE begin
        if c == SP
            @goto BOOL_SPACE
        elseif c == delim
            @endtoken
            @goto BEGIN
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == LF
            @recordtoken BOOL
            @goto LF
        elseif c == CR
            @recordtoken BOOL
            @goto CR
        else
            @multibytestring
        end
    end

    @state STRING begin
        if quoted && c == quot
            if qstring
                @recordtoken QSTRING
            else
                @recordtoken STRING
            end
            @goto QUOTE_END
        elseif c == delim
            if quoted
                @goto STRING
            end
            if qstring
                @recordtoken QSTRING
            else
                @recordtoken STRING
            end
            @endtoken
            @goto BEGIN
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == SP
            if trim && !quoted
                if qstring
                    @recordtoken QSTRING
                else
                    @recordtoken STRING
                end
                @goto STRING_SPACE
            end
            @goto STRING
        elseif c == LF
            if qstring
                @recordtoken QSTRING
            else
                @recordtoken STRING
            end
            @goto LF
        elseif c == CR
            if qstring
                @recordtoken QSTRING
            else
                @recordtoken STRING
            end
            @goto CR
        else
            @multibytestring
        end
    end

    @state STRING_SPACE begin
        if c == SP
            @goto STRING_SPACE
        elseif c == delim
            @endtoken
            @goto BEGIN
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == LF
            @goto LF
        elseif c == CR
            @goto CR
        else
            @multibytestring
        end
    end

    @state QUOTE_END begin
        if c == delim
            @endtoken
            @goto BEGIN
        elseif c == quot  # e.g. xxx," foo ""bar""",xxx
            qstring = true
            @goto STRING
        elseif c == SP
            if trim
                @goto QUOTE_END_SPACE
            end
            @goto ERROR
        elseif c == LF
            quoted = false
            @goto LF
        elseif c == CR
            quoted = false
            @goto CR
        else
            @goto ERROR
        end
    end

    @state QUOTE_END_SPACE begin
        if c == delim
            @endtoken
            @goto BEGIN
        elseif c == SP
            @goto QUOTE_END_SPACE
        elseif c == LF
            quoted = false
            @goto LF
        elseif c == CR
            quoted = false
            @goto CR
        else
            @goto ERROR
        end
    end

    @label ERROR
    if isempty(msg)
        # default error message
        msg = "invalid file format at line $(line), byte $(repr(c))"
    end
    throw(ReadError(msg))

    @label CR  # carriage return
    if pos + 1 ≤ pos_end && mem[pos + 1] == LF
        pos += 1
        # fall through
    elseif quoted
        if pos == lastnl
            # need more data
            return 0, 0
        end
        @goto STRING
    else
        @endtoken
        @goto END
    end

    @label LF  # line feed
    if quoted
        if pos == lastnl
            # need more data
            return 0, 0
        end
        @goto STRING
    else
        @endtoken
        # fall through
    end

    @label END
    return pos, i - 1
end


text:
#---------------------------------------------------------------------------------------------------
# Kraus operator
#---------------------------------------------------------------------------------------------------
import LinearAlgebra: mul!, adjoint
import Base: *, eltype, Array
export kraus
struct Kraus{TL, TU}
    KL::Array{TL, 3}
    KU::Array{TU, 3}
    dir::Symbol
end
eltype(::Kraus{TL, TU}) where TL where TU = promote_type(TL, TU)
adjoint(k::Kraus) = k.dir == :r ? Kraus(k.KL, k.KU, :l) : Kraus(k.KL, k.KU, :r)
kraus(KL, KR; dir::Symbol=:r) = Kraus(Array(KL), Array(KR), dir)
kraus(KL; dir::Symbol=:r) = Kraus(Array(KL), conj.(KL), dir)
function Array(k::Kraus)
    KL, KU, dir = k.KL, k.KU, k.dir
    α, β = size(KL, 1), size(KU, 1)
    dim = α * β
    mat = Matrix{eltype(k)}(undef, dim, dim)
    if dir == :r
        mat_r = reshape(mat, α, β, α, β)
        @tensor mat_r[:] = KL[-1,1,-3] * KU[-2,1,-4]
    else
        mat_r = reshape(mat, β, α, β, α)
        @tensor mat_r[:] = KU[-3,1,-1] * KL[-4,1,-2]
    end
    mat
end
#---------------------------------------------------------------------------------------------------
function mul!(ρ::Matrix, k::Kraus, ρ0::AbstractMatrix)
    KL, KU, dir = k.KL, k.KU, k.dir
    if dir == :r
        @tensor ρ[:] = KL[-1, 3, 1] * ρ0[1, 2] * KU[-2, 3, 2]
    elseif dir == :l
        @tensor ρ[:] = KU[1, 3, -1] * ρ0[1, 2] * KL[2, 3, -2]
    else
        error("Illegal direction: $dir.")
    end
end

function mul!(ρ::Vector, k::Kraus, ρ0::AbstractVector)
    α = size(k.KL, 1)
    ρ_r = reshape(ρ, α, α)
    ρ0_r = reshape(ρ0, α, α)
    mul!(ρ_r, k, ρ0_r)
end
#---------------------------------------------------------------------------------------------------
function *(k::Kraus, ρ0::AbstractVecOrMat)
    ctype = promote_type(eltype.((k, ρ0))...)
    ρ = Array{ctype}(undef, size(ρ0))
    mul!(ρ, k, ρ0)
    ρ
end

(k::Kraus)(v::AbstractVector) = k * v

#---------------------------------------------------------------------------------------------------
# Eigen system using power iteration

# Find dominent eigensystem by iterative multiplication.
# Krylov method ensures Hermicity and semi-positivity.
#---------------------------------------------------------------------------------------------------
function complex_iter!(K, ρ1, ρ2)
    mul!(ρ2, K, ρ1)
    mul!(ρ1, K, ρ2)
    normalize!(ρ1)
end
#---------------------------------------------------------------------------------------------------
function real_iter!(K, ρ1, ρ2)
    mul!(ρ2, K, ρ1)
    mul!(ρ1, K, ρ2)
    ρ1 .+= ρ2
    normalize!(ρ1)
end
#---------------------------------------------------------------------------------------------------
function power_iteration!(K, ρ1::Matrix, itr::Integer; method::Symbol=:r)
    ρ2 = similar(ρ1)
    if method == :r
        for i = 1:itr
            real_iter!(K, ρ1, ρ2)
        end
    elseif method == :c
        for i = 1:itr
            complex_iter!(K, ρ1, ρ2)
        end
    else
        error("Illegal method: $method.")
    end
    ρ1
end

#---------------------------------------------------------------------------------------------------
# Eigen system using Arnoldi algorithm

# Ard/noldi method ensures Hermicity but NOT semi-positivity.
#---------------------------------------------------------------------------------------------------
function krylov_arnoldi(K, ρ::Matrix)
    α = size(ρ, 1)
    v0 = reshape(ρ, :)
    e,v = eigsolve(K, v0, 1, :LR)
    reshape(v[1], α, α)
end
#---------------------------------------------------------------------------------------------------
# steady state from identity mat
function steady_mat(
    K::Array{<:Number, 3},
    itr::Integer=100;
    dir::Symbol=:r
)
    α = size(K, 3)
    Kc = conj(K)
    kraus = Kraus(K, Kc, dir)
    ρ = Array{eltype(K)}(I(α))
    mat = if α < 20
        power_iteration!(kraus, ρ, itr) |> Hermitian
    else
        m = krylov_arnoldi(kraus, ρ)
        if real(tr(m)) < 0.0
            m .= -1
        end
        m
    end
    Hermitian(mat)
end
#---------------------------------------------------------------------------------------------------
# Random fixed-point matrix.
function rand_fixed_mat(
    K::Array{<:Number, 3},
    itr::Integer=100;
    dir::Symbol=:r
)
    α = size(K, 3)
    Kc = conj(K)
    kraus = Kraus(K, Kc, dir)
    ρ = rand(ComplexF64, α, α) |> Hermitian |> Array
    mat = if α < 20
        power_iteration!(kraus, ρ, itr) |> Hermitian
    else
        krylov_arnoldi(kraus, ρ)
    end
    Hermitian(mat)
end


text:
"""
    versioninfo()

Print information about the libjpeg-turbo in use.
"""
function versioninfo()
    println("libjpeg version: ", LibJpeg.JPEG_LIB_VERSION)
    println("libjpeg-turbo version: ", LibJpeg.LIBJPEG_TURBO_VERSION)
    println("bit mode: ", LibJpeg.BITS_IN_JSAMPLE)
    println("SIMD: ", LibJpeg.WITH_SIMD == 1 ? "enabled" : "disabled")
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export daq_darkthemeprovider

"""
    daq_darkthemeprovider(;kwargs...)
    daq_darkthemeprovider(children::Any;kwargs...)
    daq_darkthemeprovider(children_maker::Function;kwargs...)

A DarkThemeProvider component.
DarkThemeProvider is a component that is placed at the root of
the component tree to make all components match the dark theme
Keyword arguments:
- `children` (Array of a list of or a singular dash component, string or numbers | a list of or a singular dash component, string or number; optional): The children of this component
- `theme` (optional): Theme object to override with a custom theme. theme has the following type: lists containing elements 'primary', 'secondary', 'detail', 'dark'.
Those elements have the following types:
  - `primary` (String; optional): Highlight color
  - `secondary` (String; optional): Supporting color
  - `detail` (String; optional): Color used for UI details, like borders
  - `dark` (Bool; optional): True for Dark mode, false for Light
"""
function daq_darkthemeprovider(; kwargs...)
        available_props = Symbol[:children, :theme]
        wild_props = Symbol[]
        return Component("daq_darkthemeprovider", "DarkThemeProvider", "dash_daq", available_props, wild_props; kwargs...)
end

daq_darkthemeprovider(children::Any; kwargs...) = daq_darkthemeprovider(;kwargs..., children = children)
daq_darkthemeprovider(children_maker::Function; kwargs...) = daq_darkthemeprovider(children_maker(); kwargs...)


text:
include("compute_edge_face_ring.jl")

function compute_boundary(face)
    """
        compute_boundary - compute the vertices on the boundary of a 3D mesh

          boundary=compute_boundary(face);

          Copyright (c) 2007 Gabriel Peyre
    """

    # compute edges (i,j) that are adjacent to only 1 face
    A = compute_edge_faces_ring(face)
    (i,j,v) = findnz(A)

    i = i[v .== -1]
    j = j[v .== -1]

    # build the boundary by traversing the edges
    boundary = [i[1]]
    i = i[2:end]
    j = j[2:end]
    while !(isempty(i))
        b = boundary[end]
        I = find(i.==b)
        if isempty(I)
            I = find(j.==b)
            if isempty(I)
                warn("Problem with boundary")
            end
            append!(boundary, i[I])
        else
            append!(boundary, j[I])
        end
        deleteat!(i, I)
        deleteat!(j, I)
    end

    #cyclic boundary
    #boundary = boundary + [boundary[0]]

    return boundary
end


text:
import Base: values
using LinearAlgebra

"""
An instance of the m-set problem.

Values are not restricted.

It can be formalised as follows:

``\\max \\sum_i \\mathrm{values}_i x_i``
``\\mathrm{s.t.} \\sum_i x_i \\leq m, \\quad x \\in \\{0, 1\\}^d``
"""
struct MSetInstance
  values::Vector{Float64}
  m::Int

  function MSetInstance(values::Vector{Float64}, m::Int)
    # Error checking.
    if m < 0
      error("m is less than zero: there is no solution.")
    end

    if m == 0
      error("m is zero: the only solution is to take no items.")
    end

    # Return a new instance.
    new(values, m)
  end
end

values(i::MSetInstance) = i.values
m(i::MSetInstance) = i.m

value(i::MSetInstance, o::Int) = values(i)[o]
values(i::MSetInstance, o) = values(i)[o]
dimension(i::MSetInstance) = length(values(i))

struct MSetSolution
  instance::MSetInstance
  items::Vector{Int} # Indices to the chosen items.
end

function value(s::MSetSolution)
  return sum(s.instance.values[i] for i in s.items)
end

function msets_greedy(instance::MSetInstance)
  # Algorithm: sort the weights, take the m largest ones, this is the optimum m-set solution.
  # Implementation: no need for sorting, partialsortperm returns the largest items.
  items = collect(partialsortperm(instance.values, 1:instance.m, rev=true))
  return MSetSolution(instance, items)
end

function msets_dp(i::MSetInstance)
  # V[µ, δ]: µ items, up to δ.
  V = Matrix{Float64}(undef, m(i), dimension(i))
  S = Dict{Tuple{Int, Int}, Vector{Int}}()

  # Initialise: µ == 1, just take the best element among [δ+1, d]; δ == d, take nothing.
  V[1, dimension(i)] = 0.0
  S[1, dimension(i)] = Int[]

  for δ in (dimension(i) - 1):-1:1
     v, x = findmax(values(i, (δ + 1):dimension(i)))
     V[1, δ] = v
     S[1, δ] = [x]
  end

  for µ in 2:m(i)
    V[µ, dimension(i)] = 0.0
    S[µ, dimension(i)] = Int[]
  end

  # Dynamic part.
  for µ in 2:m(i)
    for δ in (dimension(i) - 1):-1:0
      obj_idx = δ + 1
      take_δ = value(i, obj_idx) + V[µ - 1, δ + 1]
      dont_take_δ = V[µ, δ + 1]

      if take_δ > dont_take_δ
        V[µ, δ + 1] = take_δ
        S[µ, δ + 1] = vcat(δ, S[µ - 1, δ + 1])
      else
        V[µ, δ + 1] = dont_take_δ
        S[µ, δ + 1] = S[µ, δ + 1]
      end
    end
  end

  return MSetSolution(i, S[m(i), 1] .+ 1)
end

function msets_lp(i::MSetInstance; solver=nothing)
  model = Model(solver)
  @variable(model, x[1:length(values(i))], Bin)
  @objective(model, Max, dot(x, values(i)))
  @constraint(model, sum(x) <= m(i))

  set_silent(model)
  optimize!(model)

  return MSetSolution(i, findall(JuMP.value.(x) .>= 0.5))
end


text:
"""
A basic implementation of the order-based algorithm for queues with unobservable service.
"""
module UnobservableQueue

using SimJulia, ResumableFunctions
using Distributions, DataFrames, StatsBase

export
    Paraminf,
    infer,
    inferLCFS

include("queue.jl")
include("estimators.jl")
include("inference.jl")
include("lcfs.jl")

end


text:
function log_cost_to_go_CIR_keep_fixed_number(δ, γ, σ, λ, data, fixed_number::Int64; silence = false)
    # println("filter_WF_mem2")

    function prune_keeping_fixed_number(Λ_of_t, wms_of_t)
        return keep_fixed_number_of_weights(Λ_of_t, wms_of_t, fixed_number)
    end

    Λ_tilde_prime_of_t, logwms_tilde_of_t, θ_tilde_prime_of_t = compute_all_log_cost_to_go_functions_CIR_pruning(δ, γ, σ, λ, data, prune_keeping_fixed_number; silence = silence)

    return Λ_tilde_prime_of_t, logwms_tilde_of_t, θ_tilde_prime_of_t
end

function smooth_CIR_keep_fixed_number(δ, γ, σ, λ, data, fixed_number::Int64; silence = false)
    # println("filter_WF_mem2")
    α = δ/2
    β = γ/σ^2

    function prune_keeping_fixed_number(Λ_of_t, wms_of_t)
        Λ_of_t_kept, wms_of_t_kept = keep_fixed_number_of_weights(Λ_of_t, wms_of_t, fixed_number)
        return Λ_of_t_kept, normalise(wms_of_t_kept)
    end

    Λ_of_t, logwms_of_t, θ_of_t = filter_CIR_pruning_logweights(δ, γ, σ, λ, data, prune_keeping_fixed_number; silence = silence)
    Λ_tilde_prime_of_t, logwms_tilde_of_t, θ_tilde_prime_of_t = log_cost_to_go_CIR_keep_fixed_number(δ, γ, σ, λ, data, fixed_number::Int64; silence = silence)

    Λ_of_t_smooth, wms_of_t_smooth, θ_of_t_smooth = merge_filtering_and_cost_to_go_logweights_CIR(Λ_of_t, logwms_of_t, θ_of_t, Λ_tilde_prime_of_t, logwms_tilde_of_t, θ_tilde_prime_of_t, data, α, β)

    return Λ_of_t_smooth, wms_of_t_smooth, θ_of_t_smooth
end

function log_cost_to_go_CIR_keep_above_threshold(δ, γ, σ, λ, data, logε::Real; silence = false)
    # println("filter_WF_mem2")

    function prune_keeping_above_logthreshold(Λ_of_t, logwms_of_t)
        Λ_of_t_kept, logwms_of_t_kept = keep_above_threshold(Λ_of_t, logwms_of_t, logε)
        return Λ_of_t_kept, logwms_of_t_kept
    end

    Λ_tilde_prime_of_t, logwms_tilde_of_t, θ_tilde_prime_of_t = compute_all_log_cost_to_go_functions_CIR_pruning(δ, γ, σ, λ, data, prune_keeping_above_logthreshold; silence = silence)

    return Λ_tilde_prime_of_t, logwms_tilde_of_t, θ_tilde_prime_of_t
end

# function smooth_CIR_keep_above_threshold(δ, γ, σ, λ, data, ε_filter::Real, logε_cost_to_go::Real; silence = false)
#     # println("filter_WF_mem2")
#     α = δ/2
#     β = γ/σ^2
#
#     # function prune_keeping_above_logthreshold(Λ_of_t, logwms_of_t)
#     #     Λ_of_t_kept, logwms_of_t_kept = keep_above_threshold(Λ_of_t, logwms_of_t, log(ε_filter))
#     #     return Λ_of_t_kept, lognormalise(logwms_of_t_kept)
#     # end
#     # This is so because there is no pruning function with log internals at the moment.
#     function prune_keeping_above_threshold(Λ_of_t, wms_of_t)
#         Λ_of_t_kept, wms_of_t_kept = keep_above_threshold(Λ_of_t, wms_of_t, ε_filter)
#         return Λ_of_t_kept, normalise(wms_of_t_kept)
#     end
#
#     Λ_of_t, logwms_of_t, θ_of_t = filter_CIR_pruning_logweights(δ, γ, σ, λ, data, prune_keeping_above_threshold; silence = silence)
#
#     Λ_tilde_prime_of_t, logwms_tilde_of_t, θ_tilde_prime_of_t = log_cost_to_go_CIR_keep_above_threshold(δ, γ, σ, λ, data, logε_cost_to_go; silence = silence)
#
#     Λ_of_t_smooth, wms_of_t_smooth, θ_of_t_smooth = merge_filtering_and_cost_to_go_logweights_CIR(Λ_of_t, logwms_of_t, θ_of_t, Λ_tilde_prime_of_t, logwms_tilde_of_t, θ_tilde_prime_of_t, data, α, β)
#
#     return Λ_of_t_smooth, wms_of_t_smooth, θ_of_t_smooth
# end

function log_cost_to_go_CIR_keep_fixed_fraction(δ, γ, σ, λ, data, fraction::Float64; silence = false)
    # println("filter_WF_mem2")

    function prune_keeping_fixed_fraction(Λ_of_t, logwms_of_t)
        keep_fixed_fraction_logw(Λ_of_t, logwms_of_t, fraction, logtotal = logsumexp(logwms_of_t))
    end

    Λ_tilde_prime_of_t, logwms_tilde_of_t, θ_tilde_prime_of_t = compute_all_log_cost_to_go_functions_CIR_pruning(δ, γ, σ, λ, data, prune_keeping_fixed_fraction; silence = silence)

    return Λ_tilde_prime_of_t, logwms_tilde_of_t, θ_tilde_prime_of_t
end

function smooth_CIR_keep_fixed_fraction(δ, γ, σ, λ, data, fraction::Float64; silence = false)
    # println("filter_WF_mem2")
    α = δ/2
    β = γ/σ^2

    # function prune_keeping_above_logthreshold(Λ_of_t, logwms_of_t)
    #     Λ_of_t_kept, logwms_of_t_kept = keep_above_threshold(Λ_of_t, logwms_of_t, log(ε_filter))
    #     return Λ_of_t_kept, lognormalise(logwms_of_t_kept)
    # end
    # This is so because there is no pruning function with log internals at the moment.
    function prune_keeping_fixed_fraction(Λ_of_t, wms_of_t)
        Λ_of_t_kept, wms_of_t_kept = keep_fixed_fraction(Λ_of_t, wms_of_t, fraction)
        return Λ_of_t_kept, normalise(wms_of_t_kept)
    end

    Λ_of_t, logwms_of_t, θ_of_t = filter_CIR_pruning_logweights(δ, γ, σ, λ, data, prune_keeping_fixed_fraction; silence = silence)

    Λ_tilde_prime_of_t, logwms_tilde_of_t, θ_tilde_prime_of_t = log_cost_to_go_CIR_keep_fixed_fraction(δ, γ, σ, λ, data, fraction::Float64; silence = silence)

    Λ_of_t_smooth, wms_of_t_smooth, θ_of_t_smooth = merge_filtering_and_cost_to_go_logweights_CIR(Λ_of_t, logwms_of_t, θ_of_t, Λ_tilde_prime_of_t, logwms_tilde_of_t, θ_tilde_prime_of_t, data, α, β)

    return Λ_of_t_smooth, wms_of_t_smooth, θ_of_t_smooth
end

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


text:
"""Calculate the number of grains on square `square`."""
function on_square(square)

end

"""Calculate the total number of grains after square `square`."""
function total_after(square)

end


text:
#!/bin/env julia

# install the POMDPs.jl interface
Pkg.clone("https://github.com/sisl/POMDPs.jl.git")

# install the SARSOP wrapper
Pkg.clone("https://github.com/sisl/SARSOP.jl")
# build SARSOP, it builds from source, so this may take some time
Pkg.build("SARSOP")

# install the QMDP solver
Pkg.clone("https://github.com/sisl/QMDP.jl")

# install two helper modules
Pkg.clone("https://github.com/sisl/POMDPToolbox.jl") # this provides implementations of discrete belief updating
Pkg.clone("https://github.com/sisl/POMDPDistributions.jl") # helps with sampling

# install discrete value iteration
Pkg.clone("https://github.com/sisl/DiscreteValueIteration.jl")

# install MCTS
Pkg.clone("https://github.com/sisl/MCTS.jl")

# update
Pkg.update()

max_stars_repo_name: abrahambotros/MarsHabBot

max_issues_repo_name: abrahambotros/MarsHabBot

max_forks_repo_name: abrahambotros/MarsHabBot


text:
module Datetime

importall Base

export Calendar, ISOCalendar, Offsets, TimeZone, Offset, CALENDAR, OFFSET, Period,
	Date, DateTime, DateRange, DateRange1, DateTimeRange, DateTimeRange1,
    Year, Month, Week, Day, Hour, Minute, Second,
    year, month, week, day, hour, minute, second, millisecond,
    years,months,weeks,days,hours,minutes,seconds, milliseconds,
    addwrap, subwrap, date, datetime, unix2datetime, totaldays,
    isleap, lastdayofmonth, dayofweek, dayofyear, isdate, isdatetime,
    dayofweekinmonth, daysofweekinmonth, firstdayofweek, lastdayofweek,
    recur, now, today, calendar, timezone, offset, setcalendar, settimezone,
    Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, Sunday,
    Mon,Tue,Wed,Thu,Fri,Sat,Sun,
    January, February, March, April, May, June, July,
    August, September, October, November, December,
	Jan,Feb,Mar,Apr,Jun,Jul,Aug,Sep,Oct,Nov,Dec

abstract AbstractTime
abstract Calendar <: AbstractTime
abstract ISOCalendar <: Calendar
#Set the default calendar to use; overriding this will affect module-wide defaults
CALENDAR = ISOCalendar
setcalendar{C<:Calendar}(cal::Type{C}) = (global CALENDAR = cal)

abstract Offsets <: AbstractTime
abstract Offset{n} <: Offsets
abstract TimeZone <: Offsets
include("Timezone.jl")
#Set the default timezone to use; overriding this will affect module-wide defaults
OFFSET = UTC
settimezone{T<:Offsets}(tz::Type{T}) = (global OFFSET = tz)

abstract TimeType <: AbstractTime
bitstype 64 Date{C <: Calendar}						<: TimeType
bitstype 64 DateTime{C <: Calendar, T <: Offsets} 	<: TimeType

abstract Period 	<: AbstractTime
abstract DatePeriod <: Period
abstract TimePeriod <: Period
typealias PeriodMath   Union(Real,Period)

bitstype 32 Year{C<:Calendar}   <: DatePeriod
bitstype 32 Month{C<:Calendar}  <: DatePeriod
bitstype 32 Week{C<:Calendar}   <: DatePeriod
bitstype 32 Day{C<:Calendar}    <: DatePeriod
bitstype 32 Hour{C<:Calendar}   <: TimePeriod
bitstype 32 Minute{C<:Calendar} <: TimePeriod
bitstype 32 Second{C<:Calendar} <: TimePeriod

#Conversion/Promotion
convert{T<:TimeType}(::Type{T},x::Int64) = Base.box(T,Base.unbox(Int64,x))
convert{T<:TimeType}(::Type{Int64},x::T) = Base.box(Int64,Base.unbox(T,x))
promote_rule{T<:TimeType,R<:Real}(::Type{T},::Type{R}) = R
convert{T<:TimeType}(::Type{T},x::Real) = convert(T,int64(x))
convert{R<:Real}(::Type{R}, x::TimeType) = convert(R,int64(x))
promote_rule{C<:Calendar,T<:Offsets}(::Type{Date{C}},::Type{DateTime{C,T}}) = DateTime{C,T}
convert{T<:Offsets}(::Type{DateTime{ISOCalendar,T}},x::Date) = convert(DateTime{ISOCalendar,T},86400000*int64(x))
convert(::Type{Date{ISOCalendar}},x::DateTime) = convert(Date{ISOCalendar},_days(x))
date(x::DateTime) = convert(Date{CALENDAR},x)
datetime(x::Date) = convert(DateTime{CALENDAR,OFFSET},x)
#Date type safety
promote_rule{C<:Calendar,CC<:Calendar}(::Type{Date{C}},::Type{Date{CC}}) = Date{ISOCalendar}
convert(::Type{Date{ISOCalendar}},x::Date) = convert(Date{ISOCalendar},int64(x))
convert(::Type{Date},x::Int64) = convert(Date{ISOCalendar},x)
#DateTime type safety
promote_rule{C<:Calendar,CC<:Calendar,T<:Offsets,TT<:Offsets}(::Type{DateTime{C,T}},::Type{DateTime{CC,TT}}) = DateTime{ISOCalendar,UTC}
convert(::Type{DateTime{ISOCalendar,UTC}},x::DateTime) = convert(DateTime{ISOCalendar,UTC},int64(x))
convert(::Type{DateTime},x::Int64) = convert(DateTime{ISOCalendar,UTC},x)
hash(x::TimeType) = hash(int64(x))
length(::TimeType) = 1
isless{T<:TimeType}(x::T,y::T) = isless(int64(x),int64(y))
isless(x::TimeType,y::Real) = isless(int64(x),int64(y))
isless(x::Real,y::TimeType) = isless(int64(x),int64(y))
isless(x::TimeType,y::TimeType) = isless(promote(x,y)...)
isequal{T<:TimeType}(x::T,y::T) = isequal(int64(x),int64(y))
isequal(x::TimeType,y::Real) = isequal(int64(x),int64(y))
isequal(x::Real,y::TimeType) = isequal(int64(x),int64(y))
isequal(x::TimeType,y::TimeType) = isequal(promote(x,y)...)
isfinite(x::TimeType) = true

#Date algorithm implementations
#Source: Peter Baum - http://mysite.verizon.net/aesir_research/date/date0.htm
#Convert y,m,d to # of Rata Die days
yeardays(y::Int64) = 365y + fld(y,4) - fld(y,100) + fld(y,400)
const monthdays = [306, 337, 0, 31, 61, 92, 122, 153, 184, 214, 245, 275]
totaldays(y::Int64,m::Int64,d::Int64) = d + monthdays[m] + yeardays(m<3 ? y-1 : y) - 306
#Convert # of Rata Die days to proleptic Gregorian calendar y,m,d,w
function _year(dt::Int64)
	z = dt + 306; h = 100z - 25;	a = fld(h,3652425)
	b = a - fld(a,4); y = fld(100b+h,36525); c = b + z - 365y - fld(y,4)
	m = div(5c+456,153); return m > 12 ? y+1 : y
end
function _month(dt::Int64)
	z = dt + 306; h = 100z - 25;	a = fld(h,3652425)
	b = a - fld(a,4); y = fld(100b+h,36525); c = b + z - 365y - fld(y,4)
	m = div(5c+456,153); return m > 12 ? m-12 : m
end
function _day(dt::Int64)
	z = dt + 306; h = 100z - 25;	a = fld(h,3652425)
	b = a - fld(a,4); y = fld(100b+h,36525); c = b + z - 365y - fld(y,4)
	m = div(5c+456,153); d = c - div(153m-457,5); return d
end
function _day2date(dt::Int64)
	z = dt + 306; h = 100z - 25;	a = fld(h,3652425)
	b = a - fld(a,4); y = fld(100b+h,36525); c = b + z - 365y - fld(y,4)
	m = div(5c+456,153); d = c - div(153m-457,5); return m > 12 ? (y+1,m-12,d) : (y,m,d)
end
#https://en.wikipedia.org/wiki/Talk:ISO_week_date#Algorithms
function _week(dt::Int64)
	w = fld(abs(dt-1),7) % 20871
	c = fld((w + (w >= 10435)),5218)
	w = (w + (w >= 10435)) % 5218
	w = (w*28+(15,23,3,11)[c+1]) % 1461
	return fld(w,28) + 1
end
const DAYSINMONTH = [31,28,31,30,31,30,31,31,30,31,30,31]
isleap(y::Int64) = ((y % 4 == 0) && (y % 100 != 0)) || (y % 400 == 0)
lastdayofmonth(y::Int64,m::Int64) = DAYSINMONTH[m] + (m == 2 && isleap(y))

#TimeType constructors
setoffset(tz::Type{UTC},secs,y,s) = 0 - (y < 1972 ? 0 : s == 60 ? leaps1(secs) : leaps(secs))
getoffset(tz::Type{UTC},secs) = 0 - leaps(secs)
getoffset_secs(tz::Type{UTC},secs) = 0 - leaps1(secs)
getabr(tz::Type{UTC},secs,y) = "UTC"
#DateTime constructor with timezone and/or leap seconds
ff = open(FILEPATH*"/tzdata/"*"_leaps")
const _leaps = deserialize(ff)
ff = open(FILEPATH*"/tzdata/"*"_leaps1")
const _leaps1 = deserialize(ff)
close(ff)
# const _leaps  = [62214479999000,62230377599000,62261913599000,62293449599000,62324985599000,62356607999000,62388143999000,62419679999000,62451215999000,62498476799000,62530012799000,62561548799000,62624707199000,62703676799000,62766835199000,62798371199000,62845631999000,62877167999000,62908703999000,62956137599000,63003398399000,63050831999000,63271756799000,63366451199000,63476783999000,9223372036854775807]
# const _leaps1 = [62214480000000,62230377601000,62261913602000,62293449603000,62324985604000,62356608005000,62388144006000,62419680007000,62451216008000,62498476809000,62530012810000,62561548811000,62624707212000,62703676813000,62766835214000,62798371215000,62845632016000,62877168017000,62908704018000,62956137619000,63003398420000,63050832021000,63271756822000,63366451223000,63476784024000,9223372036854775807]
leaps(secs::Union(DateTime,Int64))  = (i = 1; while true; @inbounds (_leaps[i]  >= secs && break); i+=1 end; return 1000*(i-1))
leaps1(secs::Union(DateTime,Int64)) = (i = 1; while true; @inbounds (_leaps1[i] >= secs && break); i+=1 end; return 1000*(i-1))

date{C<:Calendar}(y::Int64,m::Int64=1,d::Int64=1,cal::Type{C}=CALENDAR) = convert(Date{cal},totaldays(y,m,d))
date{C<:Calendar}(y::PeriodMath,m::PeriodMath=1,d::PeriodMath=1,cal::Type{C}=CALENDAR) = date(int64(y),int64(m),int64(d),cal)
datetime{C<:Calendar,T<:Offsets}(y::Int64,m::Int64=1,d::Int64=1,h::Int64=0,mi::Int64=0,s::Int64=0,milli::Int64=0,tz::Type{T}=OFFSET,cal::Type{C}=CALENDAR) =
	(secs = milli + 1000*(s + 60mi + 3600h + 86400*totaldays(y,m,d)); return convert(DateTime{cal,tz}, secs - setoffset(tz,secs,y,s)))
datetime{C<:Calendar,T<:Offsets}(y::PeriodMath,m::PeriodMath=1,d::PeriodMath=1,h::PeriodMath=0,mi::PeriodMath=0,s::PeriodMath=0,cal::Type{C}=CALENDAR,tz::Type{T}=OFFSET) =
	datetime(int64(y),int64(m),int64(d),int64(h),int64(mi),int64(s),tz,cal)
datetime(y::Int64,m::Int64,d::Int64,h::Int64,mi::Int64,s::Int64,milli::Int64,tz::String) = datetime(y,m,d,h,mi,s,milli,timezone(tz),CALENDAR)
function date(s::String)
	if ismatch(r"[\/|\-|\.|,|\s]",s)
		m = match(r"[\/|\-|\.|,|\s]",s)
		a,b,c = split(s,m.match)
		y = length(a) == 4 ? int64(a) : length(c) == 4 ? int64(c) : 0
		a,b,c = int64(a),int64(b),int64(c)
		y == 0 && (y = c > 49 ? c + 1900 : c + 2000)
		m,d = y == a ? (b,c) : (a,b)
		return m > 12 ? date(y,d,m) : date(y,m,d)
	else
		error("Can't parse Date, please use date(format,datestring)")
	end
end
type DateTimeFormat
	year::Range1
	month::Range1
	day::Range1
	hour::Range1
	minute::Range1
	second::Range1
	fraction::Range1
	tz::Range1
	sep::String
	#ampm::Bool
end
function DateTimeFormat(f::String)
	y = first(search(f,"y")):last(rsearch(f,"y"))
	mon = first(search(f,"M")):last(rsearch(f,"M"))
	d = first(search(f,"d")):last(rsearch(f,"d"))
	h = first(search(f,"H")):last(rsearch(f,"H"))
	# if h == 0:-1
	# 	h = first(search(f,"H")):last(rsearch(f,"H"))
	# 	ampm = false
	# else
	# 	ampm = true
	# end
	min = first(search(f,"m")):last(rsearch(f,"m"))
	s = first(search(f,"s")):last(rsearch(f,"s"))
	frac = first(search(f,"S")):last(rsearch(f,"S"))
	sep = ismatch(r"[\/|\-|\.|,|\s]",f) ? match(r"[\/|\-|\.|,|\s]",f).match : ""
	tz = first(search(f,"z")):last(rsearch(f,"z"))
	return DateTimeFormat(y,mon,d,h,min,s,frac,tz,sep) #,ampm)
end
function dt2string(format::DateTimeFormat,dt::String)
	y = format.year == 0:-1 ? 1 : int(dt[format.year])
	mon = format.month == 0:-1 ? 1 : int(dt[format.month])
	d = format.day == 0:-1 ? 1 : int(dt[format.day])
	h = format.hour == 0:-1 ? 0 : int(dt[format.hour])
	min = format.minute == 0:-1 ? 0 : int(dt[format.minute])
	s = format.second == 0:-1 ? 0 : int(dt[format.second])
	frac = format.fraction == 0:-1 ? 0 : int(dt[format.fraction])*100
	tz = format.tz == 0:-1 ? UTC : dt[format.tz]
	return datetime(y,mon,d,h,min,s,frac,tz)
end
datetime(f::String,dt::String) = dt2string(DateTimeFormat(f),dt)
function datetime{T<:String}(f::T,t::Array{T})
	ff = DateTimeFormat(f)
	ret = DateTime{CALENDAR,OFFSET}[]
	for i in t
		push!(ret,dt2string(ff,i))
	end
	return ret
end

#Accessor/trait functions
typealias DateTimeDate Union(DateTime,Date)
_days(dt::Date) = int64(dt)
_days{C<:Calendar,T<:Offsets}(dt::DateTime{C,T}) = fld(int64(dt)+getoffset(T,dt),86400000)
year(dt::DateTimeDate) = _year(_days(dt))
month(dt::DateTimeDate) = _month(_days(dt))
week(dt::DateTimeDate) = _week(_days(dt))
day(dt::DateTimeDate) = _day(_days(dt))
hour{C<:Calendar,T<:Offsets}(dt::DateTime{C,T})     = fld(int64(dt)+getoffset(T,dt),3600000) % 24
minute{C<:Calendar,T<:Offsets}(dt::DateTime{C,T})   = fld(int64(dt)+getoffset(T,dt),60000) % 60
second{C<:Calendar,T<:Offsets}(dt::DateTime{C,T})   = (s = fld(int64(dt)+getoffset_secs(T,dt),1000) % 60; return s != 0 ? s : dt in _leaps1 ? 60 : 0)
milliseconds = (millisecond(dt::DateTime) = int64(dt) % 1000)
calendar{C<:Calendar}(dt::Date{C}) = C
typemax{D<:Date}(::Type{D}) = date(252522163911149,12,31)
typemin{D<:Date}(::Type{D}) = date(-252522163911150,1,1)
isdate(n) = typeof(n) <: Date
calendar{C<:Calendar,T<:Offsets}(dt::DateTime{C,T}) = C
timezone{C<:Calendar,T<:Offsets}(dt::DateTime{C,T}) = T
typemax{D<:DateTime}(::Type{D}) = datetime(292277024,12,31,23,59,59)
typemin{D<:DateTime}(::Type{D}) = datetime(-292277023,1,1,0,0,0)
isdatetime(x) = typeof(x) <: DateTime
#Functions to work with timezones
convert{C<:Calendar,T<:Offsets,TT<:Offsets}(::Type{DateTime{C,T}},x::DateTime{C,TT}) = convert(DateTime{C,TT},int64(x))
timezone{C<:Calendar,T<:Offsets,TT<:Offsets}(dt::DateTime{C,T},tz::Type{TT}) = convert(DateTime{C,TT},int64(dt))
offset{C<:Calendar,T<:Offsets,TT<:Offsets}(dt::DateTime{C,T},tz::Type{TT}) = convert(DateTime{C,TT},int64(dt))
timezone(x::String) = get(TIMEZONES,x,Zone0)
timezone(dt::DateTime,tz::String) = timezone(dt,timezone(tz))

#String/print/show
_s(x::Int64) = @sprintf("%02d",x)
function string(dt::Date)
	y,m,d = _day2date(_days(dt))
	y = y < 0 ? @sprintf("%05d",y) : @sprintf("%04d",y)
	return string(y,"-",_s(m),"-",_s(d))
end
function string{C<:Calendar,T<:Offsets}(dt::DateTime{C,T})
    y,m,d = _day2date(_days(dt))
    h,mi,s = hour(dt),minute(dt),second(dt)
    return string(y < 0 ? @sprintf("%05d",y) : @sprintf("%04d",y),"-",
    	_s(m),"-",_s(d),"T",_s(h),":",_s(mi),":",_s(s)," ",getabr(T,dt,y))
end
print(io::IO,x::TimeType) = print(io,string(x))
show(io::IO,x::TimeType) = print(io,string(x))

#Generic date functions
isleap(dt::DateTimeDate) = isleap(year(dt))
lastdayofmonth(dt::DateTimeDate) = lastdayofmonth(year(dt),month(dt))
dayofweek(dt::DateTimeDate) = (_days(dt) % 7)
dayofweekinmonth(dt::DateTimeDate) = (d = day(dt); return d < 8 ? 1 : d < 15 ? 2 : d < 22 ? 3 : d < 29 ? 4 : 5)
function daysofweekinmonth(dt::DateTimeDate)
	d,ld = day(dt),lastdayofmonth(dt)
	return ld == 28 ? 4 : ld == 29 ? ((d in (1,8,15,22,29)) ? 5 : 4) :
		   ld == 30 ? ((d in (1,2,8,9,15,16,22,23,29,30)) ? 5 : 4) :
		   (d in (1,2,3,8,9,10,15,16,17,22,23,24,29,30,31)) ? 5 : 4
end
firstdayofweek{C<:Calendar}(dt::Date{C}) = (d = dayofweek(dt); return convert(Date{C},int64(dt) - d + 1))
firstdayofweek{C<:Calendar,T<:Offsets}(dt::DateTime{C,T}) = (d = dayofweek(dt); return convert(DateTime{C,T},int64(dt) - 86400000*(d - 1)))
lastdayofweek{C<:Calendar}(dt::Date{C}) = (d = dayofweek(dt); return convert(Date{C},int64(dt) + (7-d)))
lastdayofweek{C<:Calendar,T<:Offsets}(dt::DateTime{C,T}) = (d = dayofweek(dt); return convert(DateTime{C,T},int64(dt) + 86400000*(7-d)))
dayofyear(dt::DateTimeDate) = _days(dt) - totaldays(year(dt),1,1) + 1
@vectorize_1arg DateTimeDate isleap
@vectorize_1arg DateTimeDate lastdayofmonth
@vectorize_1arg DateTimeDate dayofweek
@vectorize_1arg DateTimeDate dayofweekinmonth
@vectorize_1arg DateTimeDate daysofweekinmonth
@vectorize_1arg DateTimeDate firstdayofweek
@vectorize_1arg DateTimeDate lastdayofweek
@vectorize_1arg DateTimeDate dayofyear

#TimeType-specific functions (now, unix2datetime, today)
const UNIXEPOCH = 62135683200000 #Rata Die milliseconds for 1970-01-01T00:00:00 UTC
unix2datetime{T<:Offsets}(x::Int64,tz::Type{T}) = convert(DateTime{CALENDAR,tz},UNIXEPOCH + x + leaps(UNIXEPOCH + x))
now() = unix2datetime(int64(1000*time()),OFFSET)
now{T<:Offsets}(tz::Type{T}) = unix2datetime(int64(1000*time()),tz)
today() = date(now())

#TimeType arithmetic
(+)(x::TimeType) = x
(-)(x::Date) = date(-int64(year(x)),month(x),day(x))
(-){C<:Calendar,T<:Offsets}(x::DateTime{C,T}) = datetime(-year(x),month(x),day(x),hour(x),minute(x),second(x))
for op in (:+,:*,:/)
    @eval ($op)(x::DateTimeDate,y::DateTimeDate) = Base.no_op_err($op,"Date/DateTime")
end
(-){C<:Calendar}(x::Date{C},y::Date{C}) = days(-(int64(x),int64(y)))
(-){C<:Calendar,T<:Offsets}(x::DateTime{C,T},y::DateTime{C,T}) = int64(x)-int64(y)
(-)(x::DateTimeDate,y::DateTimeDate) = (-)(promote(x,y)...)
#years
function (+){C<:Calendar}(dt::Date{C},y::Year{C})
	oy,m,d = _day2date(_days(dt)); ny = oy+int64(y); ld = lastdayofmonth(ny,m)
	return date(ny,m,d <= ld ? d : ld,C)
end
function (-){C<:Calendar}(dt::Date{C},y::Year{C})
	oy,m,d = _day2date(_days(dt)); ny = oy-int64(y); ld = lastdayofmonth(ny,m)
	return date(ny,m,d <= ld ? d : ld,C)
end
function (+){C<:Calendar,T<:Offsets}(dt::DateTime{C,T},y::Year{C})
	oy,m,d = _day2date(_days(dt)); ny = oy+int64(y); ld = lastdayofmonth(ny,m)
	return datetime(ny,m,d <= ld ? d : ld,hour(dt),minute(dt),second(dt))
end
function (-){C<:Calendar,T<:Offsets}(dt::DateTime{C,T},y::Year{C})
	oy,m,d = _day2date(_days(dt)); ny = oy-int64(y); ld = lastdayofmonth(ny,m)
	return datetime(ny,m,d <= ld ? d : ld,hour(dt),minute(dt),second(dt))
end
#months
function (+){C<:Calendar}(dt::Date{C},z::Month{C})
	y,m,d = _day2date(_days(dt))
	ny = addwrap(y,m,int64(z))
	mm = addwrap(m,int64(z)); ld = lastdayofmonth(ny,mm)
	return date(ny,mm,d <= ld ? d : ld)
end
function (-){C<:Calendar}(dt::Date{C},z::Month{C})
	y,m,d = _day2date(_days(dt))
	ny = subwrap(y,m,int64(z))
	mm = subwrap(m,int64(z)); ld = lastdayofmonth(ny,mm)
	return date(ny,mm,d <= ld ? d : ld)
end
function (+){C<:Calendar,T<:Offsets}(dt::DateTime{C,T},z::Month{C})
	y,m,d = _day2date(_days(dt))
	ny = addwrap(y,m,int64(z))
	mm = addwrap(m,int64(z)); ld = lastdayofmonth(ny,mm)
	return datetime(ny,mm,d <= ld ? d : ld,hour(dt),minute(dt),second(dt))
end
function (-){C<:Calendar,T<:Offsets}(dt::DateTime{C,T},z::Month{C})
	y,m,d = _day2date(_days(dt))
	ny = subwrap(y,m,int64(z))
	mm = subwrap(m,int64(z)); ld = lastdayofmonth(ny,mm)
	return datetime(ny,mm,d <= ld ? d : ld,hour(dt),minute(dt),second(dt))
end
#weeks,days,hours,minutes,seconds
(+){C<:Calendar}(x::Date{C},y::Week{C}) = convert(Date{C},(int64(x) + 7*int64(y)))
(-){C<:Calendar}(x::Date{C},y::Week{C}) = convert(Date{C},(int64(x) - 7*int64(y)))
(+){C<:Calendar}(x::Date{C},y::Day{C})  = convert(Date{C},(int64(x) + int64(y)))
(-){C<:Calendar}(x::Date{C},y::Day{C})  = convert(Date{C},(int64(x) - int64(y)))
(+){C<:Calendar,T<:Offsets}(x::DateTime{C,T},y::Week{C})   = (x = int64(x)+604800000*int64(y) - leaps(int64(x)); convert(DateTime{C,T},x+leaps(x)))
(-){C<:Calendar,T<:Offsets}(x::DateTime{C,T},y::Week{C})   = (x = int64(x)-604800000*int64(y) - leaps(int64(x)); convert(DateTime{C,T},x+leaps(x)))
(+){C<:Calendar,T<:Offsets}(x::DateTime{C,T},y::Day{C})    = (x = int64(x)+86400000*int64(y) - leaps(int64(x)); convert(DateTime{C,T},x+leaps(x)))
(-){C<:Calendar,T<:Offsets}(x::DateTime{C,T},y::Day{C})    = (x = int64(x)-86400000*int64(y) - leaps(int64(x)); convert(DateTime{C,T},x+leaps(x)))
(+){C<:Calendar,T<:Offsets}(x::DateTime{C,T},y::Hour{C})   = (x = int64(x)+3600000*int64(y) - leaps(int64(x)); convert(DateTime{C,T},x+leaps(x)))
(-){C<:Calendar,T<:Offsets}(x::DateTime{C,T},y::Hour{C})   = (x = int64(x)-3600000*int64(y) - leaps(int64(x)); convert(DateTime{C,T},x+leaps(x)))
(+){C<:Calendar,T<:Offsets}(x::DateTime{C,T},y::Minute{C}) = (x = int64(x)+60000*int64(y) - leaps(int64(x)); convert(DateTime{C,T},x+leaps(x)))
(-){C<:Calendar,T<:Offsets}(x::DateTime{C,T},y::Minute{C}) = (x = int64(x)-60000*int64(y) - leaps(int64(x)); convert(DateTime{C,T},x+leaps(x)))
(+){C<:Calendar,T<:Offsets}(x::DateTime{C,T},y::Second{C}) = convert(DateTime{C,T},int64(x)+1000*int64(y))
(-){C<:Calendar,T<:Offsets}(x::DateTime{C,T},y::Second{C}) = convert(DateTime{C,T},int64(x)-1000*int64(y))
(+)(y::Period,x::DateTimeDate) = x + y
(-)(y::Period,x::DateTimeDate) = x - y
typealias DateTimePeriod Union(DateTimeDate,Period)
(+){T<:DateTimePeriod}(x::DateTimePeriod, y::AbstractArray{T}) = reshape([x + y[i] for i in 1:length(y)], size(y))
(+){T<:DateTimePeriod}(x::AbstractArray{T}, y::DateTimePeriod) = reshape([x[i] + y for i in 1:length(x)], size(x))
(-){T<:DateTimePeriod}(x::DateTimePeriod, y::AbstractArray{T}) = reshape([x - y[i] for i in 1:length(y)], size(y))
(-){T<:DateTimePeriod}(x::AbstractArray{T}, y::DateTimePeriod) = reshape([x[i] - fyor i in 1:length(x)], size(x))

#DateRange: for creating fixed period frequencies
immutable DateRange{C<:Calendar} <: Ranges{Date{C}}
	start::Date{C}
	step::DatePeriod
	len::Int
end
#DateRange1: for creating recurring events/dates, takes a function as step which return true for dates to be included
immutable DateRange1{C<:Calendar} <: Ranges{Date{C}}
	start::Date{C}
	step::Function
	len::Int
end
typealias DateRanges Union(DateRange,DateRange1)
size(r::DateRanges) = (r.len,)
length(r::DateRanges) = r.len
step(r::DateRanges)  = r.step
start(r::DateRanges) = 0
first(r::DateRanges) = r.start
done(r::DateRanges, i) = length(r) <= i
last(r::DateRange) = r.start + convert(typeof(r.step),(r.len-1)*r.step)
next(r::DateRange, i) = (r.start + convert(typeof(r.step),i*r.step), i+1)
function show(io::IO,r::DateRanges)
	print(io, r.start, ':', r.step, ':', last(r))
end
colon{C<:Calendar}(t1::Date{C}, y::Year{C}, t2::Date{C}) = DateRange{C}(t1, y, fld(year(t2) - year(t1),y) + 1)
colon{C<:Calendar}(t1::Date{C}, m::Month{C}, t2::Date{C}) = DateRange{C}(t1, m, fld((year(t2)-year(t1))*12+month(t2)-month(t1),m) + 1)
colon{C<:Calendar}(t1::Date{C}, w::Week{C}, t2::Date{C}) = DateRange{C}(t1, w, fld(int64(t2)-int64(t1),7w) + 1)
colon{C<:Calendar}(t1::Date{C}, d::Day{C}, t2::Date{C}) = DateRange{C}(t1, d, fld(int64(t2)-int64(t1),d) + 1)
colon{C<:Calendar}(t1::Date{C}, t2::Date{C}) = DateRange{C}(t1, day(1), fld(int64(t2)-int64(t1),day(1)) + 1)
(+){C<:Calendar}(r::DateRange{C},p::DatePeriod) = DateRange{C}(r.start+p,r.step,r.len)
(-){C<:Calendar}(r::DateRange{C},p::DatePeriod) = DateRange{C}(r.start-p,r.step,r.len)
function last(r::DateRange1)
	t = r.start
	len = 0
	while true
		r.step(t) && (len += 1)
		len == r.len && break
		t += day(1)
	end
	return t
end
function next(r::DateRange1, i)
	t = r.start
	len = 0
	while true
		r.step(t) && (len += 1)
		len == i+1 && break
		t += day(1)
	end
	return (t,i+1)
end
function colon{C<:Calendar}(t1::Date{C},fun::Function,t2::Date{C})
	len = 0
	t = d1 = t1
	while t < t2
		if fun(t)
			len += 1
			len == 1 && (d1 = t)
		end
		t += day(1)
	end
	return DateRange1{C}(d1, fun, len)
end
(+){C<:Calendar}(r::DateRange1{C},p::DatePeriod) = DateRange1{C}(r.start+p,r.step,r.len)
(-){C<:Calendar}(r::DateRange1{C},p::DatePeriod) = DateRange1{C}(r.start-p,r.step,r.len)
recur{C<:Calendar}(fun::Function,di::DateRange{C}) = colon(di.start,fun,last(di))
const Monday,Tuesday,Wednesday,Thursday,Friday,Saturday,Sunday = 1,2,3,4,5,6,7
const January,February,March,April,May,June,July,August,September,October,November,December = 1,2,3,4,5,6,7,8,9,10,11,12
const Mon,Tue,Wed,Thu,Fri,Sat,Sun = 1,2,3,4,5,6,7
const Jan,Feb,Mar,Apr,Jun,Jul,Aug,Sep,Oct,Nov,Dec = 1,2,3,4,5,6,7,8,9,10,11,12

immutable DateTimeRange{C<:Calendar,T<:Offsets} <: Ranges{DateTime{C,T}}
	start::DateTime{C,T}
	step::Period
	len::Int
end
size(r::DateTimeRange) = (r.len,)
length(r::DateTimeRange) = r.len
step(r::DateTimeRange)  = r.step
start(r::DateTimeRange) = 0
first(r::DateTimeRange) = r.start
last(r::DateTimeRange) = r.start + convert(typeof(r.step),(r.len-1)*r.step)
next(r::DateTimeRange, i) = (r.start + convert(typeof(r.step),i*r.step), i+1)
done(r::DateTimeRange, i) = length(r) <= i
function show{C<:Calendar}(io::IO,r::DateTimeRange{C})
	print(io, r.start, ':', r.step, ':', last(r))
end
colon{C<:Calendar,T<:Offsets}(t1::DateTime{C,T}, y::Year{C}, t2::DateTime{C,T}) = DateTimeRange{C,T}(t1, y, fld(year(t2) - year(t1),y) + 1)
colon{C<:Calendar,T<:Offsets}(t1::DateTime{C,T}, m::Month{C}, t2::DateTime{C,T}) = DateTimeRange{C,T}(t1, m, fld((year(t2)-year(t1))*12+month(t2)-month(t1),m) + 1)
colon{C<:Calendar,T<:Offsets}(t1::DateTime{C,T}, w::Week{C}, t2::DateTime{C,T}) = DateTimeRange{C,T}(t1, w, fld(_days(t2)-_days(t1),7w)+1)
colon{C<:Calendar,T<:Offsets}(t1::DateTime{C,T}, d::Day{C}, t2::DateTime{C,T}) = DateTimeRange{C,T}(t1, d, fld(_days(t2)-_days(t1),d)+1)
colon{C<:Calendar,T<:Offsets}(t1::DateTime{C,T}, h::Hour{C}, t2::DateTime{C,T}) = DateTimeRange{C,T}(t1, h, fld(fld(int64(t2),3600000)-fld(int64(t1),3600000),h)+1)
colon{C<:Calendar,T<:Offsets}(t1::DateTime{C,T}, mi::Minute{C}, t2::DateTime{C,T}) = DateTimeRange{C,T}(t1, mi, fld(fld(int64(t2),60000)-fld(int64(t1),60000),mi)+1)
colon{C<:Calendar,T<:Offsets}(t1::DateTime{C,T}, s::Second{C}, t2::DateTime{C,T}) = DateTimeRange{C,T}(t1, s, fld(fld(int64(t2),1000)-fld(int64(t1),1000),s)+1)
(+){C<:Calendar,T<:Offsets}(r::DateTimeRange{C,T},p::DatePeriod) = DateTimeRange{C,T}(r.start+p,r.step,r.len)
(-){C<:Calendar,T<:Offsets}(r::DateTimeRange{C,T},p::DatePeriod) = DateTimeRange{C,T}(r.start-p,r.step,r.len)

#Period types
convert{P<:Period}(::Type{P},x::Int32) = Base.box(P,Base.unbox(Int32,x))
convert{P<:Period}(::Type{Int32},x::P) = Base.box(Int32,Base.unbox(P,x))
convert{P<:Period}(::Type{P},x::Real) = convert(P,int32(x))
convert{R<:Real}(::Type{R},x::Period) = convert(R,int32(x))
for period in (Year,Month,Week,Day,Hour,Minute,Second)
	@eval convert(::Type{$period},x::Int32) = convert($period{CALENDAR},int32(x))
	@eval promote_rule{P<:$period}(::Type{P},::Type{P}) = $period{CALENDAR}
	@eval promote_rule{P<:$period,PP<:$period}(::Type{P},::Type{PP}) = $period{CALENDAR}
	@eval convert{P<:$period}(::Type{$period{ISOCalendar}},x::P) = convert($period{ISOCalendar},int32(x))
end
promote_rule{P<:Period,R<:Real}(::Type{P},::Type{R}) = R
promote_rule{A<:Period,B<:Period}(::Type{A},::Type{B}) =
	A == Second{CALENDAR} || B == Second{CALENDAR} ? Second{CALENDAR} :
    A == Minute{CALENDAR} || B == Minute{CALENDAR} ? Minute{CALENDAR} :
    A == Hour{CALENDAR}   || B == Hour{CALENDAR}   ? Hour{CALENDAR}   :
	A == Day{CALENDAR}    || B == Day{CALENDAR}    ? Day{CALENDAR}    :
    A == Week{CALENDAR}   || B == Week{CALENDAR}   ? Week{CALENDAR}   : Month{CALENDAR}
#conversion rules between periods; ISO/Greg/Julian-specific, but we define as the catchall for all Calendars
convert{C<:Calendar}(::Type{Month{C}},x::Year{C})    = convert(Month{C},  int32(12x))
convert{C<:Calendar}(::Type{Week{C}},x::Year{C})     = convert(Week{C},   int32(52x))
convert{C<:Calendar}(::Type{Day{C}},x::Year{C})      = convert(Day{C},    int32(365x))
convert{C<:Calendar}(::Type{Day{C}},x::Week{C})      = convert(Day{C},    int32(7x))
convert{C<:Calendar}(::Type{Hour{C}},x::Year{C})     = convert(Hour{C},   int32(8760x))
convert{C<:Calendar}(::Type{Hour{C}},x::Week{C})     = convert(Hour{C},   int32(168x))
convert{C<:Calendar}(::Type{Hour{C}},x::Day{C})      = convert(Hour{C},   int32(24x))
convert{C<:Calendar}(::Type{Minute{C}},x::Year{C})   = convert(Minute{C}, int32(525600x)) #Rent anybody?
convert{C<:Calendar}(::Type{Minute{C}},x::Week{C})   = convert(Minute{C}, int32(10080x))
convert{C<:Calendar}(::Type{Minute{C}},x::Day{C})    = convert(Minute{C}, int32(1440x))
convert{C<:Calendar}(::Type{Minute{C}},x::Hour{C})   = convert(Minute{C}, int32(60x))
convert{C<:Calendar}(::Type{Second{C}},x::Year{C})   = convert(Second{C}, int32(31536000x))
convert{C<:Calendar}(::Type{Second{C}},x::Week{C})   = convert(Second{C}, int32(604800x))
convert{C<:Calendar}(::Type{Second{C}},x::Day{C})    = convert(Second{C}, int32(86400x))
convert{C<:Calendar}(::Type{Second{C}},x::Hour{C})   = convert(Second{C}, int32(3600x))
convert{C<:Calendar}(::Type{Second{C}},x::Minute{C}) = convert(Second{C}, int32(60x))
#Constructors; default to CALENDAR; with (s) too cutesy? seems natural/expected though
years   = (year(x::PeriodMath)    = convert(Year{CALENDAR},x))
months  = (month(x::PeriodMath)   = convert(Month{CALENDAR},x))
weeks   = (week(x::PeriodMath)    = convert(Week{CALENDAR},x))
days    = (day(x::PeriodMath) 	  = convert(Day{CALENDAR},x))
hours   = (hour(x::PeriodMath)    = convert(Hour{CALENDAR},x))
minutes = (minute(x::PeriodMath)  = convert(Minute{CALENDAR},x))
seconds = (second(x::PeriodMath)  = convert(Second{CALENDAR},x))
#Print/show/traits
_units(x::Period) = " " * lowercase(string(typeof(x).name)) * (x == 1 ? "" : "s")
print(io::IO,x::Period) = print(io,int32(x),_units(x))
show(io::IO,x::Period) = print(io,x)
string(x::Period) = string(int32(x),_units(x))
typemin{P<:Period}(::Type{P}) = convert(P,typemin(Int32))
typemax{P<:Period}(::Type{P}) = convert(P,typemax(Int32))
zero{P<:Period}(::Union(Type{P},P)) = convert(P,int32(0))
one{P<:Period}(::Union(Type{P},P)) = convert(P,int32(1))
offset(x::Period) = return Offset{int(minutes(x))}
offset(x::Period...) = return Offset{int(minutes(sum(x...)))}
#Period Arithmetic:
isless{P<:Period}(x::P,y::P) = isless(int32(x),int32(y))
isless(x::PeriodMath,y::PeriodMath) = isless(promote(x,y)...)
isequal{P<:Period}(x::P,y::P) = isequal(int32(x),int32(y))
isequal(x::Period,y::Period) = isequal(promote(x,y)...)
isequal(x::Period,y::Real) = isequal(promote(x,y)...)
isequal(x::Real,y::Period) = isequal(promote(x,y)...)
(-){P<:Period}(x::P) = convert(P,-int32(x))
for op in (:+,:-,:%,:fld)
	@eval begin
    ($op != /) && (($op){P<:Period}(x::P,y::P) = convert(P,($op)(int32(x),int32(y))))
	($op){P<:Period}(x::P,y::Real) = ($op)(promote(x,y)...)
	($op){P<:Period}(x::Real,y::P) = ($op)(promote(x,y)...)
	!($op in (/,%,fld)) && ( ($op)(x::Period,y::Period) = ($op)(promote(x,y)...) )
    if !($op in (/,%,fld))
        ($op){P<:Period}(x::Period, y::AbstractArray{P}) = reshape([($op)(x, y[i]) for i in 1:length(y)], size(y))
        ($op){P<:Period}(x::AbstractArray{P}, y::Period) = reshape([($op)(x[i], y) for i in 1:length(x)], size(x))
	end
    end
end

for op in (:*, :/)
    @eval begin
    ($op){P<:Period}(x::P,y::P) = convert(P,($op)(int32(x),int32(y)))
	($op){P<:Period}(x::P,y::Real) = convert(P, ($op)(promote(x,y)...))
	($op){P<:Period}(x::Real,y::P) = convert(P, ($op)(promote(x,y)...))
	($op)(x::Period,y::Period) = ($op)(promote(x,y)...)
    ($op){P<:Period}(x::Period, y::AbstractArray{P}) = reshape([($op)(x, y[i]) for i in 1:length(y)], size(y))
    ($op){P<:Period}(x::AbstractArray{P}, y::Period) = reshape([($op)(x[i], y) for i in 1:length(x)], size(x))
    end
end

(/){P<:Period}(x::P,y::P) = int32(x) / int32(y)

#wrapping arithmetic
addwrap(x::Int64,y::Int64) = (v = (x + y)      % 12; return v == 0 ? 12 : v)
subwrap(x::Int64,y::Int64) = (v = (x - y + 12) % 12; return v == 0 ? 12 : v)
addwrap(x::Int64,y::Int64,z::Int64) = x + (y + z > 12 ? max(1, fld(z,12)) : 0)
subwrap(x::Int64,y::Int64,z::Int64) = x - (y - z < 1  ? max(1, fld(z,12)) : 0)

#DateRange: for creating fixed period frequencies
immutable PeriodRange{P<:Period} <: Ranges{Period}
	start::P
	step::P
	len::Int
end
size(r::PeriodRange) = (r.len,)
length(r::PeriodRange) = r.len
step(r::PeriodRange)  = r.step
start(r::PeriodRange) = 0
first(r::PeriodRange) = r.start
done(r::PeriodRange, i) = length(r) <= i
last(r::PeriodRange) = r.start + convert(typeof(r.step),(r.len-1)*r.step)
next(r::PeriodRange, i) = (r.start + convert(typeof(r.step),i*r.step), i+1)
function show(io::IO,r::PeriodRange)
	print(io, r.start, ':', r.step, ':', last(r))
end
colon{P<:Period}(t1::P, s::P, t2::P) = PeriodRange{P}(t1, s, fld(t2-t1,s) + int32(1))
colon{P<:Period}(t1::P, t2::P) = PeriodRange{P}(t1, one(P), int32(t2)-int32(t1) + int32(1))
(+){P<:Period}(r::PeriodRange{P},p::P) = PeriodRange{P}(r.start+p,r.step,r.len)
(-){P<:Period}(r::PeriodRange{P},p::P) = PeriodRange{P}(r.start-p,r.step,r.len)

end #module