text:
function schedule(employees, env)
    # Wrap the scheduling in windowing!
    start_index, end_index, time_delta = get_window(env["coverage"])

    # Apply the window
    env["time_between_coverage"] += time_delta
    if "cycle_length" in keys(env)
        env["cycle_length"] -= time_delta
    end
    if "no_shifts_after" in keys(env)
        env["no_shifts_after"] -= time_delta
    end
    env["coverage"] = apply_window(env["coverage"], start_index, end_index)
    # todo - do we have a "time since last shift" at start of week?
    for e in keys(employees)
        employees[e]["availability"] = apply_window(employees[e]["availability"], start_index, end_index)
    end
    Logging.info( "Windowing removed $time_delta hours per day (start index $start_index / end index $end_index)")

    ok, shifts = run_schedule(employees, env)
    if !ok
        Logging.info( "Scheduling not ok at windower")
        return ok, {}
    end

    # Exit the window and return
    return ok, remove_window(shifts, start_index)
end

function run_schedule(employees, env)
    if length(employees) == 0
        Logging.info( "Switching to unassigned shifts mode")
        ok, shifts = run_schedule_unassigned(env)
    else
        ok, shifts = run_schedule_standard(employees, env)
    end
    return ok, shifts
end

function run_schedule_standard(employees, env)
    attempt = 1
    employees = meet_base_coverage(env, employees)

    @label schedule_start

    ok, employees, env = build_week(employees, env)
    if !ok
        Logging.info( "Build failed")
        return false, {}
    end

    ok, message = validate_input(employees, env)
    if !ok
        Logging.info( message)
        return false, {}
    end

    # Check days per week is large enough for consec
    if size(env["coverage"], 1) > 2
        # Schedule consec first, and do non-consec as a last resort
        consecutive = true
        ok, weekly_schedule_consecutive = schedule_week(employees, env, consecutive)
        if ok
            Logging.info( "Consecutive schedule passed.")
            return true, weekly_schedule_consecutive
        end
        Logging.info( "Consecutive schedule validation failed. Trying inconsecutive.")
    end

    consecutive = false
    ok, weekly_schedule = schedule_week(employees, env, consecutive)
    if ok
        return true, weekly_schedule
    end

    # Otherwise - make an unassigned shift and restart
    name, val = generate_unassigned_shift(env)
    Logging.info( "Unable to meet coverage. Adding an unassigned shift and repeating. (end of attempt $attempt)")
    employees[name] = val
    attempt += 1
    @goto schedule_start
end

function run_schedule_unassigned(env)
    total_hours = sum(sum(env["coverage"]))
    Logging.info( "Coverage total hours: $total_hours")
    if total_hours > BIFURCATE_THRESHOLD
        # Generate subproblems
        Logging.info( "Split env into sub problems")
        ok, shifts = bifurcate_unassigned_scheduler(env)
        if !ok
            Logging.err("Bifurcatin failed in unassigned shift")
        end
        return ok, shifts
    end

    # Schedule day by day
    week_length = size(env["coverage"], 1)
    week_employees = Dict()
    shifts = Dict()
    for day in 1:week_length
        Logging.info( "Starting day scheduler day $day of $week_length")
        # Build a day env
        day_env = deepcopy(env)
        day_env["coverage"] = env["coverage"][day]

        # Generate lots of shifts
        day_employees = meet_unassigned_base_coverage(env, day)

        ok = false
        attempt = 1
        # Define day_shifts up here due to scoping
        day_shifts = Dict()
        while !ok
            # Change week employees to day employees
            day_employees_filtered = Dict()
            for e in keys(day_employees)
                d = deepcopy(day_employees[e])
                d["availability"] = day_employees[e]["availability"][day]
                d["min"] = d["shift_time_min"]
                d["max"] = d["shift_time_max"]
                day_employees_filtered[e] = d
            end

            # Run day schedule
            ok, day_shifts = schedule_day(day_employees_filtered, day_env)

            if !ok
                Logging.info( "Unable to meet coverage. Adding an unassigned shift and repeating. (end of attempt $attempt on day $day)")
                name, val = generate_unassigned_shift(env)
                day_employees[name] = val
                attempt += 1
            end
        end

        # push for testing later
        for k in keys(day_employees)
            week_employees[k] = day_employees[k]
        end

        # Now we  have to do post-processing to adjust the day correctly
        for name in keys(day_shifts)
            shift =  day_shifts[name]
            # Being really careful
            shift["day"] = day
            if !(name in keys(shifts))
                shifts[name] = Any[]
            end
            push!(shifts[name], shift)
        end
    end

    ok, message= validate_schedule(week_employees, env, shifts)
    if !ok
        Logging.err("Unassigned shift schedule invalid after post-processing - $message")
    end
    return ok, shifts
end

function schedule_week(employees, env, consecutive)

    # Setup
    valid_schedules = Any[]
    prior_lift = false
    lift_ceiling = false
    weekly_schedule = Dict{Any,Any}
    start_day = 0
    start_time = time()
    perfect_optimality = week_sum_coverage(env)

    # Start the loop, captain!
    while (prior_lift == false || (prior_lift * LYFT_THROTTLE > MIN_LYFT))
        if (time() > start_time + ITERATION_TIMEOUT) && (length(valid_schedules) == 0)
            Logging.info("Iteration timeout hit")
            break
        end

        if (time() > start_time + SEARCH_TIMEOUT) && length(valid_schedules) > 0
            Logging.info("Search timeout hit")
            break
        end

        if prior_lift != false
            lift_ceiling = prior_lift * LYFT_THROTTLE
        end

        ok, employees, prior_lift = assign_employees_to_days(employees, env, consecutive, lift_ceiling)
        if !ok
            Logging.info( "consecutive $consecutive day assignment failed")
            break
        end

        # Set up parallel calculations
        calculations = Any[]

        # Only run one scheduler if it's a single day.
        if size(env["coverage"], 1) > 1
            methods = ["balanced", "greedy"]
        else
            methods = ["balanced"]
        end

        for start_day=1:size(env["coverage"], 1)
            for method in methods
                push!(calculations, {
                    "start_day" => start_day,
                    "method" => method,
                    "employees" => employees,
                    "env" => env,
                })
            end
        end

        results = pmap(schedule_by_day, calculations)
        lift_schedules = map(return_schedules, filter(valid_schedules_net, results))
        num = length(lift_schedules)
        Logging.info( "Found $num schedules this lift")
        for valid_schedule=lift_schedules
            hours = week_hours_scheduled(valid_schedule)
            if hours == perfect_optimality
                return true, valid_schedule
            end
            push!(valid_schedules, valid_schedule)
        end
    end

    if length(valid_schedules) == 0
        return false, {}
    end

    best_hours = 0
    best_schedule = None
    for test_schedule=valid_schedules
        test_hours = week_hours_scheduled(test_schedule)
        if (test_hours < best_hours) || (best_schedule == None)
            best_hours = test_hours
            best_schedule = test_schedule
        end
    end

    if best_schedule == None
        num = length(valid_schedules)
        Logging.info( "Found $num schedules")
        return false, {}
    end
    elapsed = (time() - start_time) / 60
    num_schedules = length(valid_schedules)
    hours = week_hours_scheduled(best_schedule)
    overage = (hours - perfect_optimality) / perfect_optimality
    Logging.info( "Full time took $elapsed minutes and ended at lift $prior_lift. Found $num_schedules valid schedules. Best overage $overage")
    return true, best_schedule

end

function build_week(employees, env)
    # First check for some required shit in the environment
    if !("shift_time_min" in keys(env))
        # turn employees into a message variable
        return false, "shift_time_min not in environment", {}
    end

    if !("shift_time_max" in keys(env))
        Logging.info( "shift_time_max not in environment")
        return false, {}, {}
    end

    days_per_week = size(env["coverage"], 1)
    # Pick number of shifts per employee
    for e in keys(employees)

        # This has to be run now due to later function dependencies
        if (size(employees[e]["availability"], 1) != days_per_week)
            Logging.info( string("Employee ", e, " availability wrong number days per week"))
            return false, {}, {}
        end

        # set employee to unavailable when coverage is 0
        for day in 1:days_per_week
            for t in size(env["coverage"][day])
                if env["coverage"][day][t] == 0
                    employees[e]["availability"][day][t] = 0
                end
            end
        end

        # shift time
        if !("shift_time_min" in keys(employees[e]))
            # inherit from env if not defined
            employees[e]["shift_time_min"] = env["shift_time_min"]
        end

        if !("shift_time_max" in keys(employees[e]))
            employees[e]["shift_time_max"] = env["shift_time_max"]
        end

        if (
            "shift_count_min" in keys(employees[e])
            && "shift_count_max" in keys(employees[e])
        )
            continue
        end

        # calculate based on min/max
        if "shift_count" in keys(employees[e])
            employees[e]["shift_count_min"] = employees[e]["shift_count"]
            employees[e]["shift_count_max"] = employees[e]["shift_count"]
            continue
        end

        if !("shift_count_min" in keys(employees[e]))
            # calculate from limits
            employees[e]["shift_count_min"] = int(ceil(employees[e]["hours_min"] / employees[e]["shift_time_max"]))
        end

        if !("shift_count_max" in keys(employees[e]))
            max = int(floor(employees[e]["hours_max"] / employees[e]["shift_time_min"]))
            # can't exceed days per week
            if max > days_per_week
                max = days_per_week
            end
            employees[e]["shift_count_max"] = max
        end

        # Set ceiling of shift_count_max at availability
        days_available = days_available_per_week(employees[e])
        if employees[e]["shift_count_max"] > days_available
            employees[e]["shift_count_max"] = days_available
        end

    end

    # Anneal availability - get rid of availability that isn't feasible
    for e in keys(employees)
        employees[e]["availability"] = anneal_availability(employees[e]["availability"], employees[e]["shift_time_min"])
    end

    # day availability - num hours using longest_shift
    for e in keys(employees)
        longest_per_day = Int[]
        for day in 1:days_per_week
            push!(longest_per_day, get_longest_availability(employees[e], day))
        end
        employees[e]["longest_availability"] = longest_per_day
    end

    return true, employees, env
end

function validate_input(employees, env)
    # first build an array that is the length of each day of the week
    day_length = Int[]
    days_per_week = size(env["coverage"], 1)

    for day in 1:days_per_week
        push!(day_length, length(env["coverage"][day]))
    end

    # Days must be the same length
    overage = 0
    for day in 2:days_per_week
        if day_length[day] != day_length[1]
            return false, "ERROR: excessive day length mismatch"
        end
    end

    # Build the coverage count
    coverage_count_per_hour = Any[]
    for day in 1:size(day_length, 1)
        push!(coverage_count_per_hour, zeros(Int, day_length[day]))
    end

    # sum this from employee longest_availability variable
    coverage_count_per_day = zeros(Int, days_per_week)

    # Loop through employees
    # do shit, including adding to coverage
    for e in keys(employees)
        for t in 1:days_per_week
            if size(employees[e]["availability"][t], 1) != day_length[t]
                return false, string("Employee ", e, " availability for day ", t, " wrong number of hours")
            end
            # 1) Push hourly availability to per-hour coverage
            for h in day_length[t]
                coverage_count_per_hour[t][h] += employees[e]["availability"][t][h]
            end
            # Check their max availability per day
            coverage_count_per_day[t] += employees[e]["longest_availability"][t]
        end
        # Make copy of longest_availability then sort it; largest first
        avail = sort(employees[e]["longest_availability"], rev=true)
        # check that the sum of largest "longest avail" >= hours_min
        avail_sum_max = 0
        for i in 1:employees[e]["shift_count_max"]
            # bounds check - threw an error before
            if i <= length(avail)
                avail_sum_max += avail[i]
            end
        end

        if avail_sum_max < employees[e]["hours_min"]
            # check that sum of longest_availability > hours_min
            return false, string("Employee ", e, " not available for min_hours")
        end
    end

    for e in keys(employees)
        # Shift count validations
        shift_count_min = employees[e]["shift_count_min"]
        shift_count_max = employees[e]["shift_count_max"]
        shift_time_min = employees[e]["shift_time_min"] # inherited from env
        shift_time_max = employees[e]["shift_time_max"]
        hours_min = employees[e]["hours_min"]
        hours_max = employees[e]["hours_max"]
        days_available = days_available_per_week(employees[e])

        if shift_count_min > days_available
            return false, "Employee $e has fewer days available than minimum shift count"
        end

        if shift_count_max > days_per_week
            return false, "shift_count_max exceeds days per week"
        end

        if shift_count_min > shift_count_max
            return false, "shift count min is larger than shift count max!"
        end

        if shift_count_min * shift_time_max < hours_min
            return false, "Shift count min is too low for employee $e"
        end

        if shift_count_max * shift_time_min > hours_max
            return false, "Shift count max is too high for employee $e"
        end
    end

    return true, ""
end

function assign_employees_to_days(employees, env, consecutive=false, lift_ceiling=false)
    #=
    Consecutive days off works by looking at days_per_week, and if an employee
    is scheduled for <= (days_per_week-2) shifts, then at least two of their
    days off must be consecutive.
    =#
    days_per_week = size(env["coverage"], 1)
    number_of_employees = length(employees)
    decision_index_max = days_per_week * number_of_employees
    day_length = Int[]

    for day in 1:days_per_week
        push!(day_length, length(env["coverage"][day]))
    end

    if GUROBI
        m = Model(
            solver=GurobiSolver(
                LogToConsole=0,
                TimeLimit=CALCULATION_TIMEOUT,
                OutputFlag=0,
            )
        )
    else
        m = Model(
            solver=CbcSolver(
                threads=Base.CPU_CORES,
            )
        )
    end

    # Decision Variables
    @defVar(m, decision_variable[1:decision_index_max], Bin)

    if consecutive
        @defVar(m, decision_consecutive[1:decision_index_max], Bin)
    end

    #=
    Lyft tries to evenly allocate extra resources between days. This way,
    if you have excess capacity, it is spread evenly over days and thus
    maximizes likelihood of feasibility.
    Ignoring lift may mean that one day gets all the excess capacity. This sucks.
    =#
    @defVar(m, lift >= 1)

    decision_index_to_employee = String[]
    decision_index_to_day = Int[]
    decision_index_hours = Int[]
    for e in keys(employees)
        for t in 1:days_per_week
            push!(decision_index_to_employee, e)
            push!(decision_index_to_day, t)
            push!(decision_index_hours, employees[e]["longest_availability"][t])
        end
    end

    # Now we do per-employee constraint
    # Note: decision var 1st e.g. use x*3 instead of 3*x
    for e in keys(employees)
        # days assigned = # shifts
        # Break into to sections due to Gurobi range issue in Julia
        @addConstraint(m,
            employees[e]["shift_count_min"] <= sum{decision_variable[i], i=1:decision_index_max; decision_index_to_employee[i] == e}
        )

        @addConstraint(m,
            sum{decision_variable[i], i=1:decision_index_max; decision_index_to_employee[i] == e} <= employees[e]["shift_count_max"]
        )

        # assigned days * longest availability > min week hours
        @addConstraint(m,
            sum{decision_index_hours[i] * decision_variable[i], i=1:decision_index_max; decision_index_to_employee[i] == e} >= employees[e]["hours_min"]
        )

        if consecutive
            if employees[e]["shift_count_max"] < (days_per_week - 1)
                @addConstraint(m,
                    sum{decision_consecutive[i], i=1:decision_index_max; decision_index_to_employee[i] == e } >= 1
                )
            end
        end
    end

    for i in 1:decision_index_max
        # If an employee has no availability - don't schedule them!
        if decision_index_hours[i] == 0
            @addConstraint(m, decision_variable[i] == 0)
        end

        if consecutive
            if decision_index_to_day[i] == 1
                # See if employee worked day before this week
                if ("worked_day_preceding_week" in keys(employees[decision_index_to_employee[i]])) && (employees[decision_index_to_employee[i]]["worked_day_preceding_week"] == false)
                    @addConstraint(m, decision_consecutive[i] + decision_variable[i] == 1)
                else
                    # throw away the first index - it defaults to zero
                    @addConstraint(m, decision_consecutive[i] == 0)
                end
            else # day > 1
                # Goal: decision_consecutive[i] = 1 IF AND ONLY IF
                # decision_variable[i]=0 and decision_variable[i-1]=0
                @addConstraint(m, (decision_variable[i] + decision_variable[i-1]) >= 1 - decision_consecutive[i])
                @addConstraint(m, (2 - decision_variable[i] - decision_variable[i-1]) >= 2*decision_consecutive[i])
            end

        end
    end

    for d in 1:days_per_week
        # sum of day coverage is sufficient based on longest avail.
        @addConstraint(m,
            sum{decision_index_hours[i] * decision_variable[i], i=1:decision_index_max; decision_index_to_day[i] == d} >= sum(env["coverage"][d])*lift
        )

        # ensure hourly coverage
        for t in 1:day_length[d]
            @addConstraint(m,
                sum{employees[decision_index_to_employee[i]]["availability"][d][t] * decision_variable[i], i=1:decision_index_max; decision_index_to_day[i] == d} >= env["coverage"][d][t]
            )
        end
    end

    # For looping if infeasible models
    if lift_ceiling != false
        @addConstraint(m, lift <= lift_ceiling)
    end

    # Objective is to maximize that people get assigned to days
    # when they are free
    @setObjective(m, Max, lift)

    status = solve(m)
    if status == :Optimal
        i = 0
        for e in keys(employees)
            assigned = Int[]
            for t in 1:days_per_week
                i += 1
                push!(assigned, round(getValue(decision_variable[i])))
            end
            employees[e]["days_assigned"] = assigned
        end
        lift_value = getValue(lift)
        if days_per_week == 1
            # override. Weird bug, but we can have different lifts per day??
            lift_value = 1
        end
        Logging.info( "Lyft: ", lift_value)
        gc()
        return true, employees, lift_value
    else
        Logging.info( "Infeasible")
        gc()
        return false, {}, 0
    end
    # ok
    gc()
    return true, schedule
end

function validate_schedule(employees, env, weekly_schedule)
    days_per_week = size(env["coverage"], 1)

    day_length = Int[]
    days_per_week = size(env["coverage"], 1)

    for day in 1:days_per_week
        push!(day_length, length(env["coverage"][day]))
    end

    sum_hours = 0

    sum_coverage_by_hour = Any[]
    for day in 1:days_per_week
        push!(sum_coverage_by_hour, zeros(Int, day_length[day]))
    end

    for e in keys(weekly_schedule)
        num_shifts = length(weekly_schedule[e])

        # employees have shift_count within bounds
        if num_shifts < employees[e]["shift_count_min"]
            return false, "Employee $e has too few shifts"
        end

        if num_shifts > employees[e]["shift_count_max"]
            return false, "Employee $e has too many shifts"
        end

        employee_sum_hours = 0

        last_shift = null
        for i in 1:num_shifts
            shift_day = weekly_schedule[e][i]["day"]
            shift_start = weekly_schedule[e][i]["start"]
            shift_length = weekly_schedule[e][i]["length"]

            # Sum hours
            sum_hours += shift_length
            employee_sum_hours += shift_length

            # Coverage
            for t in shift_start:(shift_start + shift_length - 1)
                sum_coverage_by_hour[shift_day][t] += 1
            end

            # employees have shift lengths within bounds
            if shift_length < employees[e]["shift_time_min"]
                return false, "Employee $e shift $i length too short on Day $shift_day - $shift_length"
            end

            if shift_length > employees[e]["shift_time_max"]
                return false, "Employee $e shift $i length too long on Day $shift_day - $shift_length"
            end

            # INTERSHIFT
            if shift_day > 1 # check for intershift with last shift
                # find shift for preceding day - may not be in order
                for j in 1:num_shifts
                    # Intershift is not supposed to cover more than one day
                    last_shift_day = shift_day - 1
                    if weekly_schedule[e][j]["day"] == last_shift_day
                        last_shift_start = weekly_schedule[e][j]["start"]
                        last_shift_length = weekly_schedule[e][j]["length"]
                        last_shift_end = last_shift_start + last_shift_length - 1
                        hours_left_in_day = length(env["coverage"][last_shift_day]) - last_shift_end

                        if "intershift" in keys(employees[e])
                            intershift = employees[e]["intershift"]
                        else
                            intershift = env["intershift"]
                        end

                        overflow = intershift - hours_left_in_day - env["time_between_coverage"]
                        if overflow > shift_start
                            return false, "Not enough time before day $shift_day (overflow $overflow / shift start $shift_start)"
                        end
                    end
                end
            end
        end

        # employees have total hours within range
        if employee_sum_hours < employees[e]["hours_min"]
            return false, "Employee $e works too few hours in the week"
        end

        if employee_sum_hours > employees[e]["hours_max"]
            return false, "Employee $e works too many hours in the week"
        end

    end
    sum_coverage = sum(sum(env["coverage"]))

    # Hour by hour coverage
    for day in 1:days_per_week
        for hour in 1:day_length[day]
            if env["coverage"][day][hour] > sum_coverage_by_hour[day][hour]
                return false, "Insufficient employees working on day $day at hour $hour"
            end
        end
    end

    return true, "balling"
end

function valid_schedules_net(schedule)
    # For use with schedule_by_day in parallel mode
    return schedule["success"]
end

function return_schedules(schedule)
    # For use with schedule_by_day in parallel mode
    return schedule["schedule"]
end

max_stars_repo_name: shdann/autoshift

max_issues_repo_name: shdann/autoshift

max_forks_repo_name: shdann/autoshift


text:
module Fluxes

using Adapt
using DocStringExtensions

export AbstractFlux,
    FluxLW, FluxSW, init_flux_sw, set_flux_to_zero!, add_to_flux!

abstract type AbstractFlux{FT<:AbstractFloat,FTA2D<:AbstractArray{FT,2}} end

"""
    FluxLW{FT,FTA2D}

Upward, downward and net longwave fluxes at each level.

# Fields
$(DocStringExtensions.FIELDS)
"""
struct FluxLW{FT<:AbstractFloat,FTA2D<:AbstractArray{FT,2}} <:
       AbstractFlux{FT,FTA2D}
    "upward flux `[W/m²]` `(nlev,ncol)`"
    flux_up::FTA2D
    "downward flux `[W/m²]` `(nlev,ncol)`"
    flux_dn::FTA2D
    "net flux `[W/m²]` `(nlev,ncol)`"
    flux_net::FTA2D
end
FluxLW(flux_up, flux_dn, flux_net) =
    FluxLW{eltype(flux_up),typeof(flux_up)}(flux_up, flux_dn, flux_net)
Adapt.@adapt_structure FluxLW

function FluxLW(
    ncol::Int,
    nlay::Int,
    ::Type{FT},
    ::Type{DA},
) where {FT<:AbstractFloat,DA}
    flux_up = DA{FT}(undef, nlay + 1, ncol)
    flux_dn = DA{FT}(undef, nlay + 1, ncol)
    flux_net = DA{FT}(undef, nlay + 1, ncol)
    return FluxLW{FT,typeof(flux_net)}(flux_up, flux_dn, flux_net)
end

"""
    FluxSW{FT,FTA2D}

Upward, downward and net shortwave fluxes at each level.

# Fields
$(DocStringExtensions.FIELDS)
"""
struct FluxSW{FT<:AbstractFloat,FTA2D<:AbstractArray{FT,2}} <:
       AbstractFlux{FT,FTA2D}
    "upward flux `[W/m²]` `(nlev,ncol)`"
    flux_up::FTA2D
    "downward flux `[W/m²]` `(nlev,ncol)`"
    flux_dn::FTA2D
    "net flux `[W/m²]` `(nlev,ncol)`"
    flux_net::FTA2D
    "direct downward flux `[W/m²]` `(nlev,ncol)`"
    flux_dn_dir::FTA2D
end
FluxSW(flux_up, flux_dn, flux_net, flux_dn_dir) =
    FluxSW{eltype(flux_up),typeof(flux_up)}(
        flux_up,
        flux_dn,
        flux_net,
        flux_dn_dir,
    )
Adapt.@adapt_structure FluxSW

function FluxSW(
    ncol::Int,
    nlay::Int,
    ::Type{FT},
    ::Type{DA},
) where {FT<:AbstractFloat,DA}
    flux_up = DA{FT}(undef, nlay + 1, ncol)
    flux_dn = DA{FT}(undef, nlay + 1, ncol)
    flux_net = DA{FT}(undef, nlay + 1, ncol)
    flux_dn_dir = DA{FT}(undef, nlay + 1, ncol)
    return FluxSW{FT,typeof(flux_net)}(flux_up, flux_dn, flux_net, flux_dn_dir)
end

"""
    set_flux_to_zero!(flux::FluxLW{FT}) where {FT<:AbstractFloat}

Set longwave flux to zero

"""
function set_flux_to_zero!(flux::FluxLW{FT}) where {FT<:AbstractFloat}
    flux.flux_up .= FT(0)
    flux.flux_dn .= FT(0)
    flux.flux_net .= FT(0)
    return nothing
end

"""
    set_flux_to_zero!(flux::FluxSW{FT}) where {FT<:AbstractFloat}

Set shortwave flux to zero

"""
function set_flux_to_zero!(flux::FluxSW{FT}) where {FT<:AbstractFloat}
    flux.flux_up .= FT(0)
    flux.flux_dn .= FT(0)
    flux.flux_net .= FT(0)
    flux.flux_dn_dir .= FT(0)
    return nothing
end

"""
    add_to_flux!(flux1::FluxLW, flux2::FluxLW)

add longwave flux2 to longwave flux1
flux1 .+= flux2

"""
function add_to_flux!(flux1::FluxLW, flux2::FluxLW)
    flux1.flux_up .+= flux2.flux_up
    flux1.flux_dn .+= flux2.flux_dn
    flux1.flux_net .+= flux2.flux_net
    return nothing
end

"""
    add_to_flux!(flux1::FluxSW, flux2::FluxSW)

add shortwave flux2 to shortwave flux1
flux1 .+= flux2

"""
function add_to_flux!(flux1::FluxSW, flux2::FluxSW)
    flux1.flux_up .+= flux2.flux_up
    flux1.flux_dn .+= flux2.flux_dn
    flux1.flux_net .+= flux2.flux_net
    flux1.flux_dn_dir .+= flux2.flux_dn_dir
    return nothing
end

end

ext: jl


text:
module PAM

export pam

function build_phase(D,k)
    N = size(D,1)
    total_dists = [sum(D[:,j]) for j in 1:N]

    # initialize medoids with index of object with shortest distance to all others
    medoids = Int[findmin(total_dists)[2]]
    for j in 1:k-1
        TD = Vector{Float64}(undef,N)
        for a in 1:N
            td = 0.0
            for i in 1:N
                td += reduce(min, (D[i,m] for m in medoids), init=D[a,i])
            end
            TD[a] = td
        end
        push!(medoids, findmin(TD)[2])
    end

    return medoids
end

function swap_phase(D,k,M)
    M = copy(M)
    Mⱼ = similar(M)
    # Perform clustering
    assignments = Int[findmin(view(D, i,M))[2] for i in axes(D,1)]
    cumulative = similar(assignments)
    while true
        # Find minimum sum for each cluster (i.e. find the best medoid)
        for i in 1:k
            cluster = assignments .== i
            cumsum!(cumulative, cluster)
            D_slice = view(D, cluster, cluster)
            distances = sum(@view D_slice[:,i] for i in 1:last(cumulative))
            smallest_distance_idx = findmin(distances)[2]
            Mⱼ[i] = findfirst(==(smallest_distance_idx), cumulative)::Int
        end
        if sort(M) == sort(Mⱼ)
            return (medoids=M,assignments=assignments)
        else
            M,Mⱼ = Mⱼ,M
        end
    end
end

function pam(D,k)
    M = build_phase(D,k)
    return swap_phase(D,k,M)
end

end


text:
using Bedgraph
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

import Bedgraph.Record

module Bag
using Bedgraph
import Bedgraph.Record

const chroms = ["chr19", "chr19", "chr19", "chr19", "chr19", "chr19", "chr19", "chr19", "chr19"]
const firsts = [49302000, 49302300, 49302600, 49302900, 49303200, 49303500, 49303800, 49304100, 49304400]
const lasts = [49302300, 49302600, 49302900, 49303200, 49303500, 49303800, 49304100, 49304400, 49304700]
const values = [-1.0, -0.75, -0.50, -0.25, 0.0, 0.25, 0.50, 0.75, 1.00]

const browser1 = "browser position chr19:49302001-49304701"
const browser2 = "browser hide all"
const browser3 = "browser pack refGene encodeRegions"
const browser4 = "browser full altGraph"

const comment1 = "#	300 base wide bar graph, autoScale is on by default == graphing"
const comment2 = "#	limits will dynamically change to always show full range of data"
const comment3 = "#	in viewing window, priority = 20 positions this as the second graph"
const comment4 = "#	Note, zero-relative, half-open coordinate system in use for bedGraph format"

# space separated lines.
const line1 = "chr19 49302000 49302300 -1.0"
const line2 = "chr19 49302300 49302600 -0.75"
const line3 = "chr19 49302600 49302900 -0.50"
const line4 = "chr19 49302900 49303200 -0.25"
const line5 = "chr19 49303200 49303500 0.0"
const line6 = "chr19 49303500 49303800 0.25"
const line7 = "chr19 49303800 49304100 0.50"
const line8 = "chr19 49304100 49304400 0.75"
const line9 = "chr19 49304400 49304700 1.00"

const line_other_space = "2R 8225773 8226043 -0.426032509896305"
const line_other = "2R	8225773	8226043	-0.426032509896305"

const line_other_chrom = "1 3006665 3006673 2"

# Varaiations of line 1.
const line1_2 = "chr19   49302000    49302300    -1.0" # tab separated.
const line1_3 = "chr19  49302000     49302300        -1.0" # mix of tabs and spaces.
const line1_4 = " chr19 49302000 49302300 -1.0" # space at start.
const line1_5 = "    chr19 49302000 49302300 -1.0" # tab at start.
const line1_6 = "chr19 49302000 49302300 -1.0 " # space at end.
const line1_7 = "chr19 49302000 49302300 -1.0    " # tab at end.

const cells1 = ["chr19", "49302000", "49302300", "-1.0"]

const record1 = Record("chr19", 49302000, 49302300, -1.0)

const parameter_line_min = "track type=bedGraph"
const parameter_line = "track type=bedGraph name=\"BedGraph Format\" description=\"BedGraph format\" visibility=full color=200, 100, 0 altColor=0, 100, 200 priority=20"
const parameter_line_4 = "track type=bedGraph name=track_label description=center_label"
const parameter_line_long = "track type=bedGraph name=track_label description=center_label visibility=display_mode color=r, g, b altColor=r, g, b priority=priority autoScale=on|off alwaysZero=on|off gridDefault=on|off maxHeightPixels=max:default:min graphType=bar|points viewLimits=lower:upper yLineMark=real-value yLineOnOff=on|off windowingFunction=maximum|mean|minimum smoothingWindow=off|2-16"

const file = joinpath(@__DIR__, "data.bedgraph")
const file_headerless = joinpath(@__DIR__, "data-headerless.bedgraph")

const header = [browser1, browser2, browser3, browser4, comment1, comment2, comment3, comment4, parameter_line]
const records = [record1, Record(line2), Record(line3), Record(line4), Record(line5), Record(line6), Record(line7), Record(line8), Record(line9)]

end # module bag

@testset "Bedgraph" begin

@testset "Record Constructor" begin
	@test_nowarn Record("chr1", 1.0, 1.0, 1.1)
	@test Record("chr1", 1.0, 1.0, 1.0) == Record("chr1", 1, 1, 1.0)
	@test Record("chr1", "1.0", "1", "1.1") == Record("chr1", 1, 1, 1.1)
end

@testset "I/O" begin

	@test isfile(Bag.file)
	@test isfile(Bag.file_headerless)

	# Seek test.
	open(Bag.file, "r") do io
		@info "" readline(io, keep = true)
	    seek(io, Record)
	    @test position(io) == 536 #TODO: determine why the returned position on windows machines is 545.
	    @test readline(io) == Bag.line1
	end

	# Check things for headerless Bag.files.
	open(Bag.file_headerless, "r") do io

		# Check that the first record of a headerless bedGraph file can be sought.
	    seek(io, Record)
	    @test position(io) == 0
	    @test readline(io) == Bag.line1 # IO position is at the start of line 2.

		# Check behaviour of consecutive calls to seek(io, Record).
		seek(io, Record)
		seek(io, Record)
		@test readline(io) == Bag.line2

	end

	open(Bag.file) do io
		seek(io, Record)
		@test read!(io, Vector{Bedgraph.Record}(undef, length(Bag.records))) ==  Bag.records
	end

	@test_nowarn Bedgraph.BedgraphHeader()
	@test_nowarn Bedgraph.BedgraphHeader{Vector{String}}()

	@test read(Bag.file, Vector{Bedgraph.Record}) ==  Bag.records
	@test read(Bag.file, Bedgraph.BedgraphHeader{Vector{String}}).data == Bag.header

	open(Bag.file, "r") do io # Note: reading records first to check seek.
	    @test read(io, Vector{Record}) == Bag.records
		@test read(io, Bedgraph.BedgraphHeader).data == Bag.header
		@test read(io, Bedgraph.BedgraphHeader{Vector{String}}).data == Bag.header
	end

	open(Bag.file_headerless, "r") do io # Note: reading records first to check seek.
		@test read(io, Vector{Record}) == Bag.records
		@test read(io, Bedgraph.BedgraphHeader).data == []
	    @test read(io, Bedgraph.BedgraphHeader{Vector{String}}).data == []
	end

	@test Bag.records == open(Bag.file, "r") do io
		records = Vector{Record}()

		while !eof(seek(io, Bedgraph.Record))
	        record = read(io, Bedgraph.Record) #Note: no protection.
	        push!(records, record)
	    end

	    return records
	end

	@test Bag.records == open(Bag.file, "r") do io
	    seek(io, Bedgraph.Record)
	    return read(io, Vector{Bedgraph.Record})
	end

	@test_nowarn mktemp() do path, io
		header = Bedgraph.generate_basic_header(Bag.records)
	    write(io, header, Bag.records)

		seekstart(io)
		@test Bag.records == read(io, Vector{Record})
	end

	@test Record{Int} <: Record

	@test_nowarn mktemp() do path, io
	    write(io, Bag.records[[1,5,9]])

		seekstart(io)
		@test Bag.records[[1,5,9]] == read(io, Vector{Record{Int}})
	end

	@test_nowarn mktemp() do path, io
		header = Bedgraph.generate_basic_header(Bag.records)
	    write(io, header, Bag.records)

		seekstart(io)
		@test Bag.records == read(io, Vector{Record{Float64}})
	end

	# Not able to capture write(::IO, ::AbstractVector).
	# @test_nowarn mktemp() do path, io
	# 	# test AbstractVector
	# 	records = view(Bag.records, 2:4)
	# 	header = Bedgraph.generateBasicHeader(records)
	#     write(io, header, records)
	#
	# 	seekstart(io)
	# 	@test records == read(io, Vector{Record})
	# end

	@test_nowarn mktemp() do path, io
		header = Bedgraph.generate_basic_header("chr19", Bag.records[1].first, Bag.records[end].last, bump_forward=false)
		write(io, header, Bag.records)

		seekstart(io)
		@test Bag.records == read(io, Vector{Record})
	end

	@test_nowarn mktempdir() do path
		outputfile = joinpath(path, "test.bedgraph")
		header = Bedgraph.generate_basic_header(Bag.records)
		write(outputfile, header, Bag.records)

		@test Bag.records == read(outputfile, Vector{Record})
	end

end #testset I/O

@testset "Matching" begin

	@test Bedgraph.is_comment(Bag.comment1)
	@test Bedgraph.is_browser(Bag.browser3)

	@test Bedgraph.is_like_record("1 2 3 4") == true
	@test Bedgraph.is_like_record(Bag.parameter_line) == false
	@test Bedgraph.is_like_record(Bag.parameter_line_4) == false
	@test Bedgraph.is_like_record(Bag.parameter_line_min) == false
	@test Bedgraph.is_like_record(Bag.parameter_line_long) == false
	@test Bedgraph.is_like_record(Bag.line1) == true
	@test Bedgraph.is_like_record(Bag.line1_2) == true
	@test Bedgraph.is_like_record(Bag.line1_3) == true
	@test Bedgraph.is_like_record(Bag.line1_4) == true
	@test Bedgraph.is_like_record(Bag.line1_5) == true
	@test Bedgraph.is_like_record(Bag.line1_6) == true
	@test Bedgraph.is_like_record(Bag.line1_7) == true

	@test Bedgraph.is_like_record(Bag.line_other_space) == true
	@test Bedgraph.is_like_record(Bag.line_other) == true
	@test Bedgraph.is_like_record(Bag.line_other_chrom) == true

end #testset Matching

@testset "Line Splitting" begin

	@test Bedgraph._split_line(Bag.line1)[1:4] == Bag.cells1
	@test Bedgraph._split_line(Bag.line1_2)[1:4] == Bag.cells1
	@test Bedgraph._split_line(Bag.line1_3)[1:4] == Bag.cells1
	@test Bedgraph._split_line(Bag.line1_4)[1:4] == Bag.cells1
	@test Bedgraph._split_line(Bag.line1_5)[1:4] == Bag.cells1
	@test Bedgraph._split_line(Bag.line1_6)[1:4] == Bag.cells1
	@test Bedgraph._split_line(Bag.line1_7)[1:4] == Bag.cells1

end #testset Parsing

@testset "Sorting" begin

	# Testing isless function.
	@test isless(Record(Bag.line1), Record(Bag.line2)) == true
	@test isless(Record(Bag.line2), Record(Bag.line1)) == false

	unordered = Bag.records[[9,4,5,6,8,7,1,3,2]]

	@test Bag.records != unordered
	@test Bag.records == sort(unordered)

	# Testing sort by concept.
	@test Bag.records == sort(unordered, by=x->[x.chrom, x.first, x.last])

end #testset Sorting

@testset "Conversion" begin

	@test Bag.record1 == convert(Record, Bag.line1)
	@test Bag.record1 == Record(Bag.line1)
	@test Bag.record1 == convert(Record, string(Bag.line1, " ", "extra_cell"))
	@test Bag.record1 == Record(string(Bag.line1, " ", "extra_cell"))

end #testset Conversion

@testset "Internal Helpers" begin

	@test Bedgraph._range(Bag.record1) == Bag.record1.first : Bag.record1.last - 1
	@test Bedgraph._range(Bag.record1, right_open=false) == (Bag.record1.first + 1 ) : Bag.record1.last

	@test Bedgraph._range(Bag.records) == Bag.record1.first : Record(Bag.line9).last - 1
	@test Bedgraph._range(Bag.records, right_open=false) == Bag.record1.first + 1 : Record(Bag.line9).last

	bumped_records = Bedgraph._bump_forward(Bag.records)
	@test bumped_records[1].first == (Bag.records[1].first + 1)
	@test bumped_records[1].last == (Bag.records[1].last + 1)

	bumped_records = Bedgraph._bump_back(Bag.records)
	@test bumped_records[1].first == (Bag.records[1].first - 1)
	@test bumped_records[1].last == (Bag.records[1].last - 1)

end #testset Internal Helpers

@testset "NamedTuple" begin

	record = Record(Bag.line1)

	nt = (
		chrom = Bag.chroms[1],
		first = Bag.firsts[1],
		last = Bag.lasts[1],
		value = Bag.values[1],
	)

	@test record == Record(nt)

	@test convert(NamedTuple{(:chrom, :first, :last, :value)}, record) == convert(NamedTuple, record) == nt

	# Check renaming of fields.
	@test convert(NamedTuple{(:chrom, :left, :right, :value)}, record) == (
		chrom = Bag.chroms[1],
		left = Bag.firsts[1],
		right = Bag.lasts[1],
		value = Bag.values[1],
	)
end

end # total testset


text:
# Precompile statements generated using SnoopCompile, via:
#   julia> using SnoopCompile, Test, FileIO
#   julia> tinf = @snoopi_deep begin include("test/rw.jl"); include("test/loadsave.jl"); end
#   julia> ttot, pcs = ttot, pcs = SnoopCompile.parcel(tinf; tmin=0.05)
#   julia> SnoopCompile.write("tmp", pcs)
# With the cutoff at tmin = 0.05, this accounts for about 40% of overall inference time
# spent in executing /test/rw.jl and /test/loadsave.jl
# A single CodecZlib argument call was filtered out manually.

function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing

    Base.precompile(Tuple{typeof(write),JLDFile{MmapIO},String,Float64})   # time: 0.8394732
    Base.precompile(Tuple{typeof(write),JLDFile{MmapIO},String,Base.ReshapedArray{Int, 2, UnitRange{Int}, Tuple{}}})   # time: 0.5689673
    Base.precompile(Tuple{Core.kwftype(typeof(jldopen)),NamedTuple{(:compress,), Tuple{Bool}},typeof(jldopen),String,Bool,Bool,Bool,Type{MmapIO}})   # time: 0.3199544
    Base.precompile(Tuple{Core.kwftype(typeof(write)),NamedTuple{(:compress,), Tuple{Nothing}},typeof(write),Group{JLDFile{MmapIO}},String,Base.ReshapedArray{Int, 2, UnitRange{Int}, Tuple{}},JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.2051841
    Base.precompile(Tuple{typeof(read_data),JLDFile{MmapIO},ReadRepresentation{ComplexF64, ComplexF64},Tuple{ReadDataspace, RelOffset, Int, UInt16},Vector{ReadAttribute}})   # time: 0.187065
    Base.precompile(Tuple{typeof(read_data),JLDFile{MmapIO},ReadRepresentation{BitMatrix, OnDiskRepresentation{(0, 8, 16), Tuple{Vector{UInt}, Int, Tuple{Int, Int}}, Tuple{RelOffset, Int, Tuple{Int, Int}}}()},Tuple{ReadDataspace, RelOffset, Int, UInt16},Vector{ReadAttribute}})   # time: 0.1830797
    Base.precompile(Tuple{Core.kwftype(typeof(write)),NamedTuple{(:compress,), Tuple{Nothing}},typeof(write),Group{JLDFile{IOStream}},String,Base.ReshapedArray{Int, 2, UnitRange{Int}, Tuple{}},JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.172313
    Base.precompile(Tuple{typeof(write),JLDFile{MmapIO},String,Dict{Symbol, String}})   # time: 0.1482486
    Base.precompile(Tuple{typeof(write),JLDFile{MmapIO},String,IdDict{Symbol, String}})   # time: 0.14119
    Base.precompile(Tuple{typeof(write),JLDFile{MmapIO},String,Vector{Vector{Float64}}})   # time: 0.141007
    Base.precompile(Tuple{typeof(constructrr),JLDFile{MmapIO},DataType,CompoundDatatype,Vector{ReadAttribute}})   # time: 0.1277466
    Base.precompile(Tuple{typeof(h5fieldtype),JLDFile{MmapIO},Type{UnitRange{Int}},Type{UnitRange{Int}},Type{Val{true}}})   # time: 0.1257442
    Base.precompile(Tuple{typeof(write_dataset),JLDFile{MmapIO},WriteDataspace{0, Tuple{}},CommittedDatatype,OnDiskRepresentation{(0, 8, 16), Tuple{Vector{UInt}, Int, Tuple{Int, Int}}, Tuple{RelOffset, Int, Tuple{Int, Int}}},BitMatrix,JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.1243035
    Base.precompile(Tuple{typeof(write),JLDFile{MmapIO},String,Dict{String, String},JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.1156945
    Base.precompile(Tuple{typeof(read_data),JLDFile{MmapIO},ReadRepresentation{DataType, OnDiskRepresentation{(0, 16), Tuple{String, Vector{Any}}, Tuple{Vlen{String}, Vlen{RelOffset}}}()},Tuple{ReadDataspace, RelOffset, Int, UInt16},Nothing})   # time: 0.115149
    Base.precompile(Tuple{typeof(write),JLDFile{MmapIO},String,Vector{Vector{U} where U<:Integer}})   # time: 0.10986
    Base.precompile(Tuple{typeof(getindex),Group{JLDFile{IOStream}},String})   # time: 0.1069834
    Base.precompile(Tuple{typeof(write),JLDFile{MmapIO},String,Vector{Array{Int}}})   # time: 0.1035001
    Base.precompile(Tuple{typeof(setindex!),JLDFile{MmapIO},Vector{Union{Missing, Bool}},String})   # time: 0.0995374
    Base.precompile(Tuple{typeof(close),JLDFile{MmapIO}})   # time: 0.0988005
    Base.precompile(Tuple{typeof(read_data),JLDFile{MmapIO},ReadRepresentation{Base.ImmutableDict{Symbol, String}, CustomSerialization{Array, RelOffset}},Tuple{ReadDataspace, RelOffset, Int, UInt16},Vector{ReadAttribute}})   # time: 0.0972386
    Base.precompile(Tuple{Core.kwftype(typeof(write)),NamedTuple{(:compress,), Tuple{Nothing}},typeof(write),Group{JLDFile{MmapIO}},String,Dict{String, String},JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.0919869
    Base.precompile(Tuple{typeof(write),JLDFile{MmapIO},String,Tuple{Int, Int, Tuple{Int, Int, Vector{Int}}, Vector{Int}}})   # time: 0.0905843
    Base.precompile(Tuple{typeof(write),JLDFile{MmapIO},String,Vector{Tuple{}}})   # time: 0.0843623
    Base.precompile(Tuple{typeof(write),JLDFile{MmapIO},String,Tuple{Int, String}})   # time: 0.0840445
    Base.precompile(Tuple{typeof(read),JLDFile{MmapIO},String})   # time: 0.0816563
    Base.precompile(Tuple{typeof(write),JLDFile{MmapIO},String,BigFloat})   # time: 0.0796746
    Base.precompile(Tuple{typeof(h5convert!),IndirectPointer,Type{CustomSerialization{Vector{Pair{Any, Any}}, RelOffset}},JLDFile{MmapIO},IdDict,JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.0794924
    Base.precompile(Tuple{typeof(write_dataset),JLDFile{IOStream},WriteDataspace{0, Tuple{}},CommittedDatatype,OnDiskRepresentation{(0, 8, 16), Tuple{Vector{UInt}, Int, Tuple{Int, Int}}, Tuple{RelOffset, Int, Tuple{Int, Int}}},BitMatrix,JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.0782433
    Base.precompile(Tuple{Core.kwftype(typeof(write)),NamedTuple{(:compress,), Tuple{Nothing}},typeof(write),Group{JLDFile{MmapIO}},String,Float64,JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.0779112
    Base.precompile(Tuple{typeof(write),JLDFile{MmapIO},String,Vector{Vector{Int}}})   # time: 0.0767249
    Base.precompile(Tuple{typeof(write_dataset),JLDFile{MmapIO},WriteDataspace{1, Tuple{WrittenAttribute{WriteDataspace{0, Tuple{}}, BasicDatatype, RelOffset}}},BasicDatatype,Type{RelOffset},Vector{Union{Missing, Float32, Float64, Int32}},JLDWriteSession{Dict{UInt, RelOffset}},Bool})   # time: 0.0764436
    Base.precompile(Tuple{typeof(write),JLDFile{MmapIO},String,Base.ReshapedArray{ComplexF64, 1, Base.ReinterpretArray{ComplexF64, 2, Float64, Matrix{Float64}, false}, Tuple{}}})   # time: 0.0764319
    Base.precompile(Tuple{typeof(write),JLDFile{MmapIO},String,BitMatrix})   # time: 0.0733186
    Base.precompile(Tuple{Core.kwftype(typeof(write)),NamedTuple{(:compress,), Tuple{Nothing}},typeof(write),Group{JLDFile{MmapIO}},String,Vector{Array{Int}},JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.0726416
    Base.precompile(Tuple{typeof(h5type),JLDFile{MmapIO},Type,Base.ReshapedArray{Int, 2, UnitRange{Int}, Tuple{}}})   # time: 0.0719121
    Base.precompile(Tuple{typeof(write),JLDFile{MmapIO},String,ComplexF32})   # time: 0.07121
    Base.precompile(Tuple{Core.kwftype(typeof(write)),NamedTuple{(:compress,), Tuple{Nothing}},typeof(write),Group{JLDFile{IOStream}},String,Vector{Array{Int}},JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.0702584
    Base.precompile(Tuple{typeof(write_dataset),JLDFile{MmapIO},WriteDataspace{0, Tuple{}},CommittedDatatype,OnDiskRepresentation{(0, 8), Tuple{TypeVar, Any}, Tuple{RelOffset, RelOffset}},Type,JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.0697
    Base.precompile(Tuple{Core.kwftype(typeof(write)),NamedTuple{(:compress,), Tuple{Nothing}},typeof(write),Group{JLDFile{MmapIO}},String,Dict{Symbol, String},JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.0691246
    Base.precompile(Tuple{typeof(write),JLDFile{MmapIO},String,SubArray{Int, 1, Vector{Int}, Tuple{UnitRange{Int}}, true}})   # time: 0.0687224
    Base.precompile(Tuple{typeof(write_ref),JLDFile{IOStream},Vector{Int},JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.0681966
    Base.precompile(Tuple{typeof(write),JLDFile{MmapIO},String,Base.ReshapedArray{Int, 4, UnitRange{Int}, Tuple{}}})   # time: 0.068107
    Base.precompile(Tuple{typeof(setindex!),JLDFile{MmapIO},Vector{Union{Float64, Int}},String})   # time: 0.0679686
    Base.precompile(Tuple{typeof(write_ref),JLDFile{MmapIO},LineNumberNode,JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.0676823
    Base.precompile(Tuple{typeof(read_array),JLDFile{MmapIO},ReadDataspace,ReadRepresentation{Union{Int, Vector}, RelOffset},Int,UInt16,RelOffset,Vector{ReadAttribute}})   # time: 0.0675551
    Base.precompile(Tuple{typeof(setindex!),JLDFile{MmapIO},Vector{Union{Missing, Float32, Float64, Int32}},String})   # time: 0.0671966
    Base.precompile(Tuple{typeof(write),JLDFile{MmapIO},String,Vector{Union{Nothing, Int}}})   # time: 0.0671739
    Base.precompile(Tuple{typeof(setindex!),JLDFile{MmapIO},Vector{Union{Int, Vector}},String})   # time: 0.0666241
    Base.precompile(Tuple{typeof(write),JLDFile{MmapIO},String,Complex{Int}})   # time: 0.066388
    Base.precompile(Tuple{typeof(write_ref),JLDFile{MmapIO},Tuple{Symbol},JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.0653194
    Base.precompile(Tuple{typeof(write),JLDFile{MmapIO},String,Base.GitVersionInfo})   # time: 0.0652343
    Base.precompile(Tuple{typeof(write),JLDFile{IOStream},String,Vector{Char}})   # time: 0.0641642
    Base.precompile(Tuple{Core.kwftype(typeof(write)),NamedTuple{(:compress,), Tuple{Nothing}},typeof(write),Group{JLDFile{MmapIO}},String,NamedTuple{(:ptr,), Tuple{Ptr{Float64}}},JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.0624435
    Base.precompile(Tuple{Core.kwftype(typeof(write)),NamedTuple{(:compress,), Tuple{Nothing}},typeof(write),Group{JLDFile{IOStream}},String,Dict{Symbol, String},JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.0622999
    Base.precompile(Tuple{typeof(jldopen),AbstractString})   # time: 0.061449
    Base.precompile(Tuple{typeof(write),JLDFile{MmapIO},String,Expr})   # time: 0.0613437
    Base.precompile(Tuple{typeof(write),JLDFile{MmapIO},String,NamedTuple{(:ptr,), Tuple{Ptr{Float64}}},JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.0605991
    Base.precompile(Tuple{Core.kwftype(typeof(write)),NamedTuple{(:compress,), Tuple{Nothing}},typeof(write),Group{JLDFile{IOStream}},String,Tuple{Int, Int, Tuple{Int, Int, Vector{Int}}, Vector{Int}},JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.0597551
    Base.precompile(Tuple{typeof(write_ref_mutable),JLDFile{IOStream},Matrix{Int},JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.0595954
    Base.precompile(Tuple{Core.kwftype(typeof(write)),NamedTuple{(:compress,), Tuple{Nothing}},typeof(write),Group{JLDFile{IOStream}},String,Vector{Vector{U} where U<:Integer},JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.0593596
    Base.precompile(Tuple{Core.kwftype(typeof(write)),NamedTuple{(:compress,), Tuple{Nothing}},typeof(write),Group{JLDFile{IOStream}},String,Vector{Tuple{}},JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.0589484
    Base.precompile(Tuple{typeof(read_data),JLDFile{MmapIO},ReadRepresentation{ComplexF32, ComplexF32},Tuple{ReadDataspace, RelOffset, Int, UInt16},Vector{ReadAttribute}})   # time: 0.0588968
    Base.precompile(Tuple{Core.kwftype(typeof(write)),NamedTuple{(:compress,), Tuple{Nothing}},typeof(write),Group{JLDFile{IOStream}},String,Base.ReshapedArray{ComplexF64, 1, Base.ReinterpretArray{ComplexF64, 2, Float64, Matrix{Float64}, false}, Tuple{}},JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.0586527
    Base.precompile(Tuple{typeof(write),JLDFile{MmapIO},String,Matrix{Int}})   # time: 0.058283
    Base.precompile(Tuple{typeof(setindex!),JLDFile{IOStream},Vector{UInt8},String})   # time: 0.058052
    Base.precompile(Tuple{typeof(write_ref),JLDFile{IOStream},LineNumberNode,JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.0580383
    Base.precompile(Tuple{typeof(constructrr),JLDFile{IOStream},DataType,CompoundDatatype,Vector{ReadAttribute}})   # time: 0.0569429
    Base.precompile(Tuple{typeof(write_heap_object),JLDFile{IOStream},Type,Vector{RelOffset},JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.056889
    Base.precompile(Tuple{typeof(write_ref),JLDFile{MmapIO},Matrix{Float64},JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.0561273
    Base.precompile(Tuple{typeof(write),JLDFile{MmapIO},String,Array{Int, 0}})   # time: 0.0560988
    Base.precompile(Tuple{Core.kwftype(typeof(write)),NamedTuple{(:compress,), Tuple{Nothing}},typeof(write),Group{JLDFile{IOStream}},String,Float64,JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.0560133
    Base.precompile(Tuple{typeof(write),JLDFile{MmapIO},String,Vector{Vector{U} where U>:Int}})   # time: 0.0542099
    Base.precompile(Tuple{Core.kwftype(typeof(write)),NamedTuple{(:compress,), Tuple{Nothing}},typeof(write),Group{JLDFile{MmapIO}},String,Expr,JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.0540752
    Base.precompile(Tuple{Core.kwftype(typeof(write)),NamedTuple{(:compress,), Tuple{Nothing}},typeof(write),Group{JLDFile{IOStream}},String,Vector{Vector{Float64}},JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.0536389
    Base.precompile(Tuple{typeof(write),JLDFile{MmapIO},String,Vector{Vector{U} where Int<:U<:Real}})   # time: 0.0536368
    Base.precompile(Tuple{Core.kwftype(typeof(write)),NamedTuple{(:compress,), Tuple{Nothing}},typeof(write),Group{JLDFile{IOStream}},String,Vector{Any},JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.0534297
    Base.precompile(Tuple{typeof(write),JLDFile{MmapIO},String,Base.ReshapedArray{Int, 3, UnitRange{Int}, Tuple{}}})   # time: 0.0530808
    Base.precompile(Tuple{Core.kwftype(typeof(write)),NamedTuple{(:compress,), Tuple{Nothing}},typeof(write),Group{JLDFile{MmapIO}},String,SubArray{Int, 1, Vector{Int}, Tuple{UnitRange{Int}}, true},JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.0526628
    Base.precompile(Tuple{typeof(write),JLDFile{MmapIO},String,Matrix{Float64}})   # time: 0.0525738
    Base.precompile(Tuple{typeof(write_dataset),JLDFile{MmapIO},WriteDataspace{0, Tuple{}},CommittedDatatype,OnDiskRepresentation{(0, 8), Tuple{Int, Union{Nothing, Symbol}}, Tuple{Int, RelOffset}},LineNumberNode,JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.0525689
    Base.precompile(Tuple{typeof(h5fieldtype),JLDFile{MmapIO},Type{Tuple{Int, Int, Int}},Type{Tuple{Int, Int, Int}},Type{Val{true}}})   # time: 0.0525214
    Base.precompile(Tuple{typeof(h5convert!),IndirectPointer,Type{RelOffset},JLDFile{MmapIO},Array{Float64},JLDWriteSession{Dict{UInt, RelOffset}}})   # time: 0.0522247
    Base.precompile(Tuple{typeof(write),JLDFile{MmapIO},String,Vector{DataType}})   # time: 0.0515368
    Base.precompile(Tuple{typeof(fileio_save),File{FileIO.DataFormat{:JLD2}, String},String,String,String,Vararg{Any}})   # time: 0.0513424
    Base.precompile(Tuple{typeof(h5type),JLDFile{IOStream},Type,Base.ReshapedArray{Int, 2, UnitRange{Int}, Tuple{}}})   # time: 0.050349
    Base.precompile(Tuple{typeof(read_data),JLDFile{IOStream},FixedLengthString{String},Tuple{ReadDataspace, RelOffset, Int, UInt16},Vector{ReadAttribute}})   # time: 0.0502741
    Base.precompile(Tuple{typeof(read_data),JLDFile{IOStream},ReadRepresentation{LineNumberNode, OnDiskRepresentation{(0, 8), Tuple{Int, Union{Nothing, Symbol}}, Tuple{Int, RelOffset}}()},Tuple{ReadDataspace, RelOffset, Int, UInt16},Vector{ReadAttribute}})   # time: 0.0500727
end


text:
# Particle IDs used for testing. The definition is the same as in
# the Python package 'particle' from the SciKit-HEP group:
# https://github.com/scikit-hep/particle/blob/master/tests/conftest.py

@enum PDGIDS begin
    # Gauge and Higgs bosons
    Gluon = 21
    Photon = 22
    Z0 = 23
    WMinus = -24
    HiggsBoson = 25
    ZPrime = 32
    # Charged leptons
    Electron = 11
    Positron = -11
    Muon = 13
    AntiMuon = -13
    Tau = 15
    TauPrime = 17
    # Neutrinos
    Nu_e = 12
    NuBar_tau = -16
    # Quarks
    DQuark = 1
    UQuark = 2
    SQuark = 3
    CQuark = 4
    BQuark = 5
    TQuark = 6
    BPrimeQuark = 7  # 4th generation
    TPrimeQuark = 8
    # Quarkonia
    jpsi = 443
    psi_2S = 100443
    Upsilon_1S = 553
    Upsilon_4S = 300553
    # Light hadrons
    Pi0 = 111
    PiPlus = 211
    eta = 221
    eta_prime = 331
    a_0_1450_plus = 10211
    KL = 130
    KS = 310
    KMinus = -321
    rho_770_minus = -213
    phi = 333
    omega = 223
    K1_1270_0 = 10313
    K1_1400_0 = 20313
    rho_1700_0 = 30113
    a2_1320_minus = -215
    omega_3_1670 = 227
    f_4_2300 = 9010229  # example of a not-well-known meson
    Proton = 2212
    AntiNeutron = -2112
    Lambda = 3122
    Sigma0 = 3212
    SigmaPlus = 3222
    SigmaMinus = 3112
    Xi0 = 3322
    AntiXiMinus = -3312
    OmegaMinus = 3334
    # Charm hadrons
    D0 = 421
    DPlus = 411
    DsPlus = 431
    LcPlus = 4122
    # Beauty hadrons
    B0 = 511
    BPlus = 521
    Bs = 531
    BcPlus = 541
    Lb = 5122
    # Top hadrons
    T0 = 621
    LtPlus = 6122
    # Special particles
    Graviton = 39
    Reggeon = 110
    Pomeron = 990
    Odderon = 9990
    # Supersymmetric particles
    Gluino = 1000021
    Gravitino = 1000039
    STildeL = 1000003
    CTildeR = 2000004
    # R-hadrons
    RPlus_TTildeDbar = 1000612
    R0_GTildeG = 1000993
    RPlusPlus_GTildeUUU = 1092224
    # Q-balls
    QBall1 = 10000150
    QBall2 = -10000200
    # Dyons
    DyonSameMagElecChargeSign = 4110010
    DyonOppositeMagElecChargeSign = 4120010
    # Di-quarks
    DD1 = 1103
    SD0 = 3101
    # Nuclei
    HydrogenNucleus = 1000010010
    Carbon12 = 1000060120
    # Pentaquarks
    AntiUCbarCUDPentaquark = -9422144
    # example of spin 3/2 u-cbar-c-u-d pentaquark decaying to J/psi proton
    UCbarCUDPentaquark = 9422144
    # Technicolor
    Pi0TC = 3000111
    PiMinusTC = -3000211
    # Composite quarks and leptons
    UQuarkStar = 4000002
    AntiElectronStar = -4000011
    # Generator specific pseudoparticles or concepts
    AntiCHadron = -84
    # Invalid ID
    Invalid1 = 0  # illegal ID
    Invalid2 = 99999999  # general form is a 7-digit number
end


text:
module Service

using Dates, ExpiringCaches
using ..DataModel, ..Schema # , ..Auth

function getSystemAlerts()
    resp = Schema.getSystemAlerts()
    return resp
end

# function createAlbum(obj)
#     @assert haskey(obj, :name) && !isempty(obj.name)
#     @assert haskey(obj, :artist) && !isempty(obj.artist)
#     @assert haskey(obj, :songs) && !isempty(obj.songs)
#     @assert haskey(obj, :year) && 1900 < obj.year < Dates.year(Dates.now())
#     album = Album(obj.name, obj.artist, obj.year, obj.songs)
#     Mapper.create!(album)
#     return album
# end

# @cacheable Dates.Hour(1) function getAlbum(id::Int64)::Album
#     Mapper.get(id)
# end

# function updateAlbum(id, updated)
#     album = Mapper.get(id)
#     album.name = updated.name
#     album.artist = updated.artist
#     album.year = updated.year
#     album.songs = updated.songs
#     Mapper.update(album)
#     delete!(ExpiringCaches.getcache(getAlbum), (id,))
#     return album
# end

# function deleteAlbum(id)
#     Mapper.delete(id)
#     delete!(ExpiringCaches.getcache(getAlbum), (id,))
#     return
# end

# function pickAlbumToListen()
#     albums = Mapper.getAllAlbums()
#     leastTimesPicked = minimum(x -> x.timespicked, albums)
#     leastPickedAlbums = filter(x -> x.timespicked == leastTimesPicked, albums)
#     pickedAlbum = rand(leastPickedAlbums)
#     pickedAlbum.timespicked += 1
#     Mapper.update(pickedAlbum)
#     delete!(ExpiringCaches.getcache(getAlbum), (pickedAlbum.id,))
#     @info "picked album = $(pickedAlbum.name) on thread = $(Threads.threadid())"
#     return pickedAlbum
# end

# function createUser(user)
#     @assert haskey(user, :username) && !isempty(user.username)
#     @assert haskey(user, :password) && !isempty(user.password)
#     user = User(user.username, user.password)
#     Mapper.create!(user)
#     return user
# end

# function loginUser(user)
#     persistedUser = Mapper.get(user)
#     if persistedUser.password == user.password
#         persistedUser.password = ""
#         return persistedUser
#     else
#         throw(Auth.Unauthenticated())
#     end
# end

end # module


text:
@testset "DictObject" begin
    # """Test creation and property access on DictObjects"""
    d = DictObject(Dict("args_working?"=>"yes"), a=123, b="foobar", nerds=["chris","jen","janey"])
    @test d["a"]==123
    #@test d["a"]==123
    @test d["b"]=="foobar"
    #@test d["b"]=="foobar"
    @test d["nerds"]==["chris","jen","janey"]
    @test hasattr(d,"nerds") == true
    #@test d["nerds"]==["chris","jen","janey"]
    @test hasattr(d,"qwerqwer") == false

    # println(keys(d))
    @test all([key in keys(d) for key in ["args_working?", "a", "b", "nerds"]]) == true
    println(d)
    d["new_key"] = "new value!"
    @test d["new_key"] == "new value!"
end


text:
# ---------------------------------------------------------------------------- #
#
#   mesherFEFLOA.jl
#
#   feflo.a-specific read and write
#
#   λυτέος
#   Fall 2017
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

"""
    MesherFEFLOA

MesherFEFLOA type:
Type for FEFLOA meshes.
"""
struct MesherFEFLOA <: Mesher

end

"""
    writeMsh( mesh::Mesh2D, mesher::MesherFEFLOA, flname::String )

Writes 2D mesh in `mesh` to `flname` in format used by feflo.a.
"""
function writeMsh( mesh::Mesh2D, mesher::MesherFEFLOA, flname::String )

    np = size(mesh.p,1)

    fid = open( flname, "w" )

    @printf( fid, "MeshVersionFormatted 2\n\n" )

    @printf(fid, "Dimension 2\n\n")

    # Vertices
    @printf(fid, "Vertices\n%i\n", np)
    for jj = 1:size(mesh.p,1)
        @printf(fid, "%16.15e\t%16.15e\t%i\n", mesh.p[jj,1], mesh.p[jj,2], 1)
    end
    @printf(fid,"\n")

    # Triangles
    @printf(fid, "Triangles\n%i\n",size(mesh.t,1))
    for jj = 1:size(mesh.t,1)
        @printf(fid, "%i\t%i\t%i\t%i\n", mesh.t[jj,1], mesh.t[jj,2], mesh.t[jj,3], 1)
    end
    @printf(fid,"\n")

    # Edges
    nbound = size(mesh.fb,1)
    @printf(fid, "Edges\n%i\n",nbound)
    for ii = 1:nbound
        indf = mesh.fb[ii,1]
        @printf(fid, "%i\t%i\t%i\n", mesh.f[indf,1], mesh.f[indf,2], mesh.fb[ii,2])
    end
    @printf(fid,"\n")

    close(fid)

    # TODO: Corners

end

"""
    writeMsh( mesh::Mesh3D, mesher::MesherFEFLOA, flname::String )

Writes 3D mesh in `mesh` to `flname` in format used by feflo.a.
"""
function writeMsh( mesh::Mesh3D, mesher::MesherFEFLOA, flname::String )

    np = size(mesh.p,1)

    fid = open( flname, "w" )

    @printf( fid, "MeshVersionFormatted 2\n\n" )

    @printf(fid, "Dimension 3\n\n")

    # Vertices
    @printf(fid, "Vertices\n%i\n", np)
    for jj = 1:size(mesh.p,1)
        @printf(fid, "%16.15e\t%16.15e\t%16.15e\t%i\n", mesh.p[jj,1], mesh.p[jj,2],mesh.p[jj,3], 1)
    end
    @printf(fid,"\n")

    # Triangles
    @printf(fid, "Tetrahedra\n%i\n",size(mesh.t,1))
    for jj = 1:size(mesh.t,1)
        @printf(fid, "%i\t%i\t%i\t%i\t%i\n", mesh.t[jj,1], mesh.t[jj,2], mesh.t[jj,3], mesh.t[jj,4], 1)
    end
    @printf(fid,"\n")

    # Edges
    nbound = size(mesh.fb,1)
    @printf(fid, "Triangles\n%i\n",nbound)
    for ii = 1:nbound
        indf = mesh.fb[ii,1]
        @printf(fid, "%i\t%i\t%i\t%i\n", mesh.f[indf,1], mesh.f[indf,2], mesh.f[indf,3], mesh.fb[ii,2])
    end
    @printf(fid,"\n")

    close(fid)

    # TODO: Corners

end

"""
    readFEFLOA_2D( flname::String )

Reads 2D feflo.a mesh from `flname`.
"""
function readFEFLOA_2D( flname::String )

    dim = 2

    # Open file
    fid = open( flname, "r" )

    # Check dimension
    x = readline(fid)
    while isempty(x) || x[1] != 'D'
        x = readline(fid)
    end
    xx = split( x )
    if parse(Int64, xx[end]) != dim
        error(" FEFLOA read-in: dimension should be 2")
    end

    # Read vertices
    while isempty(x) || x[1] != 'V'
        x = readline(fid)
    end
    x  = readline(fid)
    np = parse( Int64, x )
    p  = fill( 0.0, np, 2 )
    for il in 1:np
        x  = readline(fid)
        xx = split(x)
        p[il,1] = parse( Float64, xx[1] )
        p[il,2] = parse( Float64, xx[2] )
    end

    # Read triangles
    while isempty(x) || x[1] != 'T'
        x = readline(fid)
    end
    x  = readline(fid)
    nt = parse( Int64, x )
    t  = fill( 0, nt, 3 )
    for il in 1:nt
        x  = readline(fid)
        xx = split(x)
        t[il,1] = parse( Int64, xx[1] )
        t[il,2] = parse( Int64, xx[2] )
        t[il,3] = parse( Int64, xx[3] )
    end

    # Read boundary
    while isempty(x) || x[1] != 'E'
        x = readline(fid)
    end
    x   = readline(fid)
    nb  = parse( Int64, x )
    bel = fill( 0, nb, 3 )
    for il in 1:nb
        x  = readline(fid)
        xx = split(x)
        bel[il,1] = parse( Int64, xx[1] )
        bel[il,2] = parse( Int64, xx[2] )
        bel[il,3] = parse( Int64, xx[3] )
    end

    close(fid)

    # Add boundary tags
    tags = Vector{String}( nb )
    for bb in 1:nb
        tags[bb] = string("bound", bb)
    end

    return( p, t, bel, tags )

end

"""
    readFEFLOA_3D( flname::String )

Reads 3D feflo.a mesh from `flname`.
"""
function readFEFLOA_3D( flname::String )

    dim = 3

    # Open file
    fid = open( flname, "r" )

    # Check dimension
    x = readline(fid)
    while isempty(x) || x[1] != 'D'
        x = readline(fid)
    end
    xx = split( x )
    if parse(Int64, xx[end]) != dim
        error(" FEFLOA read-in: dimension should be 3")
    end

    # Read vertices
    while isempty(x) || x[1] != 'V'
        x = readline(fid)
    end
    x  = readline(fid)
    np = parse( Int64, x )
    p  = fill( 0.0, np, 3 )
    for il in 1:np
        x  = readline(fid)
        xx = split(x)
        p[il,1] = parse( Float64, xx[1] )
        p[il,2] = parse( Float64, xx[2] )
        p[il,3] = parse( Float64, xx[3] )
    end

    # Read triangles
    while isempty(x) || x[1] != 'T'
        x = readline(fid)
    end
    x  = readline(fid)
    nt = parse( Int64, x )
    t  = fill( 0, nt, 4 )
    for il in 1:nt
        x  = readline(fid)
        xx = split(x)
        t[il,1] = parse( Int64, xx[1] )
        t[il,2] = parse( Int64, xx[2] )
        t[il,3] = parse( Int64, xx[3] )
        t[il,4] = parse( Int64, xx[4] )
    end

    # Read boundary
    while isempty(x) || x[1] != 'T'
        x = readline(fid)
    end
    x   = readline(fid)
    nb  = parse( Int64, x )
    bel = fill( 0, nb, 4 )
    for il in 1:nb
        x  = readline(fid)
        xx = split(x)
        bel[il,1] = parse( Int64, xx[1] )
        bel[il,2] = parse( Int64, xx[2] )
        bel[il,3] = parse( Int64, xx[3] )
        bel[il,4] = parse( Int64, xx[4] )
    end

    close(fid)

    # Add boundary tags
    nbounds = maximum( bel[:,4] )
    tags = Vector{String}( nbounds )
    for bb in 1:nbounds
        tags[bb] = string("bound", bb)
    end

    return( p, t, bel, tags )

end

"""
    runMesher( mesher::MesherFEFLOA, flMesh::String, flMetric::String, flMeshNew::String )

Runs FEFLOA mesher.
"""
function runMesher( mesher::MesherFEFLOA, flMesh::String, flMetric::String, flMeshNew::String )

    if flMetric[end-3:end] != ".sol"
        mv( flMetric, string(flMeshNew,".sol") )
    end
    run( `fefloa -in $flMesh -met $flMetric -out $flMeshNew` )
    # mv( string(flMeshNew,".mesh"), flMeshNew, remove_destination=true )

end


text:
#!/usr/bin/env bash
    #=
    exec julia --project="$(realpath $(dirname $(dirname $0)))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \
    "${BASH_SOURCE[0]}" "$@"
    =#

module CodingTheory

ENV["HECKE_PRINT_BANNER"] = "false"
# ENV["NEMO_PRINT_BANNER"] = "false"

# using Nemo: isprime, factor, fmpz # for prime power function
using Hecke: ispower, isprime_power, isprime # for perfect power function
# using Primes: isprime, primes
using LinearAlgebra: I
using FLoops: @floop, ThreadedEx
using Polynomials, StaticArrays

# from utils.jl
export sizeof_perfect_code, sizeof_all_words, has_identity, has_identity_on_left

include("utils.jl")

# Abstract types
export FinitePolynomial, AbstractCode, NonStaticAbstractWord, AbstractWord, Word,
        Codewords, Alphabet, Messages, CodeUniverse, CodeUniverseIterator,
        UniverseParameters
export no_round, getindex, setindex!, firstindex, lastindex, size, length, rand,
        gensym, genalphabet, eltype, isword, isabstractword

# RREF
export rref, rref!

# Bounds, Messages, Distance, and Primes
export rate, sphere_covering_bound, sphere_packing_bound, hamming_bound,
        singleton_bound, gilbert_varshamov_bound, elias_bassalygo_bound,
        plotkin_bound, construct_ham_matrix, isperfect, ishammingperfect,
        isgolayperfect
export get_codewords_greedy, get_codewords_random, get_all_words, get_codewords
export hamming_distance, hamming_ball, code_distance, t_error_detecting,
        t_error_correcting, find_error_detection_max, find_error_correction_max
export isprimepower

# Algebra
export FinitePolynomial, list_polys, multiplication_table, list_span, islinear,
        isirreducible, normal_form!, normal_form, equivalent_code!, equivalent_code,
        generator!, generator, parity_check, syndrome, isincode

# Powers
export isprimepower, isperfectpower

# Levenshtein
export levenshtein, levenshtein!

include("abstract_types.jl")
include("rref.jl")
include("messages.jl") # implicitly exports distance.jl and powers.jl
include("bounds.jl")
include("algebra.jl")
include("levenshtein.jl")

end # end module


text:
@doc doc"""
    FixedRankMatrices{m,n,k,T} <: Manifold

The manifold of $m\times n$ real-valued (complex-valued) matrices of fixed rank $k$, i.e.
````math
\mathcal M = \{ x \in \mathbb R^{m\times n} : \operatorname{rank}(x) = k \}.
````
# Representation with 3 matrix factors

A point $x\in\mathcal M$ can be stored using orthonormal matrices
$U\in\mathbb R^{m\times k}$, $V\in\mathbb R^{n\times k}$ as well as the $k$ singular
values of $x = USV^\mathrm{T}$. In other words, $U$ and $V$ are from the manifolds
[`Stiefel`](@ref)`(m,k)` and [`Stiefel`](@ref)`(n,k)`, respectively; see
[`SVDMPoint`](@ref) for details

The tangent space $T_x\mathcal M$ at a point $x\in\mathcal M$ with $x=USV^\mathrm{T}$
is given by
````math
T_x\mathcal M = \bigl\{ UMV^\mathrm{T} + U_xV^\mathrm{T} + UV_x^\mathrm{T} :
    M \in \mathbb R^{k\times k},
    U_x \in \mathbb R^{m\times k},
    V_x \in \mathbb R^{n\times k}
    \text{ s.t. }
    U_x^\mathrm{T}U = 0_k,
    V_x^\mathrm{T}V = 0_k
\bigr\},
````
where $0_k$ is the $k\times k$ zero matrix. See [`UMVTVector`](@ref) for details.

The (default) metric of this manifold is obtained by restricting the metric
on $\mathbb R^{m\times n}$ to the tangent bundle. This implementation follows[^Vandereycken2013].

# Constructor
    FixedRankMatrics(m,n,k,t=ℝ)

Generate the manifold of `m`-by-`n` real-valued matrices of rank `k`.

[^Vandereycken2013]:
    > Bart Vandereycken: "Low-rank matrix completion by Riemannian Optimization,
    > SIAM Journal on Optiomoization, 23(2), pp. 1214–1236, 2013.
    > doi: [10.1137/110845768](https://doi.org/10.1137/110845768),
    > arXiv: [1209.3834](https://arxiv.org/abs/1209.3834).
"""
struct FixedRankMatrices{M,N,K,T} <: Manifold end
function FixedRankMatrices(m::Int, n::Int, k::Int, t::AbstractNumbers = ℝ)
    return FixedRankMatrices{m,n,k,t}()
end

@doc doc"""
    SVDMPoint <: MPoint

A point on a certain manifold, where the data is stored in a svd like fashion,
i.e. in the form $USV^\mathrm{T}$, where this structure stores $U$, $S$ and
$V^\mathrm{T}$. The storage might also be shortened to just $k$ singular values
and accordingly shortened $U$ (columns) and $V^\mathrm{T}$ (rows)

# Constructors
* `SVDMPoint(A)` for a matrix `A`, stores its svd factors (i.e. implicitly $k=\min\{m,n\}$)
* `SVDMPoint(S)` for an `SVD` object, stores its svd factors (i.e. implicitly $k=\min\{m,n\}$)
* `SVDMPoint(U,S,Vt)` for the svd factors to initialize the `SVDMPoint`` (i.e. implicitly $k=\min\{m,n\}$)
* `SVDMPoint(A,k)` for a matrix `A`, stores its svd factors shortened to the
  best rank $k$ approximation
* `SVDMPoint(S,k)` for an `SVD` object, stores its svd factors shortened to the
  best rank $k$ approximation
* `SVDMPoint(U,S,Vt,k)` for the svd factors to initialize the `SVDMPoint`,
  stores its svd factors shortened to the best rank $k$ approximation
"""
struct SVDMPoint{TU<:AbstractMatrix,TS<:AbstractVector,TVt<:AbstractMatrix} <: MPoint
    U::TU
    S::TS
    Vt::TVt
end
SVDMPoint(A::AbstractMatrix) = SVDMPoint(svd(A))
SVDMPoint(S::SVD) = SVDMPoint(S.U, S.S, S.Vt)
SVDMPoint(A::Matrix, k::Int) = SVDMPoint(svd(A), k)
SVDMPoint(S::SVD, k::Int) = SVDMPoint(S.U, S.S, S.Vt, k)
SVDMPoint(U, S, Vt, k::Int) = SVDMPoint(U[:, 1:k], S[1:k], Vt[1:k, :])
==(x::SVDMPoint, y::SVDMPoint) = (x.U == y.U) && (x.S == y.S) && (x.Vt == y.Vt)

@doc doc"""
    UMVTVector <: TVector

A tangent vector that can be described as a product $UMV^\mathrm{T}$, at least
together with its base point, see for example [`FixedRankMatrices`](@ref)

# Constructors
* `UMVTVector(U,M,Vt)` store umv factors to initialize the `UMVTVector`
* `UMVTVector(U,M,Vt,k)` store the umv factors after shortening them down to
  inner dimensions $k$, i.e. in $UMV^\mathrm{T}$, $M\in\mathbb R^{k\times k}$
"""
struct UMVTVector{TU<:AbstractMatrix,TM<:AbstractMatrix,TVt<:AbstractMatrix} <: TVector
    U::TU
    M::TM
    Vt::TVt
end

UMVTVector(U, M, Vt, k::Int) = UMVTVector(U[:, 1:k], M[1:k, 1:k], Vt[1:k, :])

# here the division in M corrects for the first factor in UMV + x.U*Vt + U*x.Vt, where x is the base point to v.
*(v::UMVTVector, s::Number) = UMVTVector(v.U * s, v.M * s, v.Vt * s)
*(s::Number, v::UMVTVector) = UMVTVector(s * v.U, s * v.M, s * v.Vt)
/(v::UMVTVector, s::Number) = UMVTVector(v.U / s, v.M / s, v.Vt / s)
\(s::Number, v::UMVTVector) = UMVTVector(s \ v.U, s \ v.M, s \ v.Vt)
+(v::UMVTVector, w::UMVTVector) = UMVTVector(v.U + w.U, v.M + w.M, v.Vt + w.Vt)
-(v::UMVTVector, w::UMVTVector) = UMVTVector(v.U - w.U, v.M - w.M, v.Vt - w.Vt)
-(v::UMVTVector) = UMVTVector(-v.U, -v.M, -v.Vt)
+(v::UMVTVector) = UMVTVector(v.U, v.M, v.Vt)
==(v::UMVTVector, w::UMVTVector) = (v.U == w.U) && (v.M == w.M) && (v.Vt == w.Vt)

@doc doc"""
    check_manifold_point(M::FixedRankMatrices{m,n,k},x; kwargs...)

Check whether the matrix or [`SVDMPoint`](@ref) `x` ids a valid point on the
[`FixedRankMatrices`](@ref)`{m,n,k}` `M`, i.e. is (or represents) an `m`-by`n` matrix of
rank `k`. For the [`SVDMPoint`](@ref) the internal representation also has to have the right
shape, i.e. `x.U` and `x.Vt` have to be unitary.
"""
function check_manifold_point(M::FixedRankMatrices{m,n,k}, x; kwargs...) where {m,n,k}
    r = rank(x; kwargs...)
    s = "The point $(x) does not lie on the manifold of fixed rank matrices of size ($(m),$(n)) witk rank $(k), "
    if size(x) != (m, n)
        return DomainError(size(x), string(s, "since its size is wrong."))
    end
    if r > k
        return DomainError(r, string(s, "since its rank is too large ($(r))."))
    end
    return nothing
end
function check_manifold_point(
    F::FixedRankMatrices{m,n,k},
    x::SVDMPoint;
    kwargs...,
) where {m,n,k}
    s = "The point $(x) does not lie on the manifold of fixed rank matrices of size ($(m),$(n)) witk rank $(k), "
    if (size(x.U) != (m, k)) || (length(x.S) != k) || (size(x.Vt) != (k, n))
        return DomainError(
            [size(x.U)..., length(x.S), size(x.Vt)...],
            string(
                s,
                "since the dimensions do not fit (expected $(n)x$(m) rank $(k) got $(size(x.U,1))x$(size(x.Vt,2)) rank $(size(x.S)).",
            ),
        )
    end
    if !isapprox(x.U' * x.U, one(zeros(n, n)); kwargs...)
        return DomainError(
            norm(x.U' * x.U - one(zeros(n, n))),
            string(s, " since U is not orthonormal/unitary."),
        )
    end
    if !isapprox(x.Vt' * x.Vt, one(zeros(n, n)); kwargs...)
        return DomainError(
            norm(x.Vt' * x.Vt - one(zeros(n, n))),
            string(s, " since V is not orthonormal/unitary."),
        )
    end
    return nothing
end

@doc doc"""
    check_tangent_vector(M:FixedRankMatrices{m,n,k}, x, v)

Check whether the tangent [`UMVTVector`](@ref) `v` is from the tangent space of
the [`SVDMPoint`](@ref) `x` on the [`FixedRankMatrices`](@ref) `M`, i.e. that
`v.U` and `v.Vt` are (columnwise) orthogonal to `x.U` and `x.Vt`, respectively,
and its dimensions are consistent with `x` and `M`, i.e. correspond to `m`-by-`n`
matrices of rank `k`.
"""
function check_tangent_vector(
    M::FixedRankMatrices{m,n,k},
    x::SVDMPoint,
    v::UMVTVector;
    kwargs...,
) where {m,n,k}
    c = check_manifold_point(M, x)
    c === nothing || return c
    if (size(v.U) != (m, k)) || (size(v.Vt) != (k, n)) || (size(v.M) != (k, k))
        return DomainError(
            cat(size(v.U), size(v.M), size(v.Vt), dims = 1),
            "The tangent vector $(v) is not a tangent vector to $(x) on the fixed rank matrices since the matrix dimensions to not fit (expected $(m)x$(k), $(k)x$(k), $(k)x$(n)).",
        )
    end
    if !isapprox(v.U' * x.U, zeros(k, k); kwargs...)
        return DomainError(
            norm(v.U' * x.U - zeros(k, k)),
            "The tangent vector $(v) is not a tangent vector to $(x) on the fixed rank matrices since v.U'x.U is not zero. ",
        )
    end
    if !isapprox(v.Vt * x.Vt', zeros(k, k); kwargs...)
        return DomainError(
            norm(v.Vt * x.Vt - zeros(k, k)),
            "The tangent vector $(v) is not a tangent vector to $(x) on the fixed rank matrices since v.V'x.V is not zero.",
        )
    end
end

@doc doc"""
    inner(M::FixedRankMatrices, x::SVDMPoint, v::UMVTVector, w::UMVTVector)

Compute the inner product of `v` and `w` in the tangent space of `x` on the
[`FixedRankMatrices`](@ref) `M`, which is inherited from the embedding, i.e. can be computed
using `dot` on the elements (`U`, `Vt`, `M`) of `v` and `w`.
"""
function inner(::FixedRankMatrices, x::SVDMPoint, v::UMVTVector, w::UMVTVector)
    return dot(v.U, w.U) + dot(v.M, w.M) + dot(v.Vt, w.Vt)
end

function isapprox(::FixedRankMatrices, x::SVDMPoint, y::SVDMPoint; kwargs...)
    return isapprox(x.U * Diagonal(x.S) * x.Vt, y.U * Diagonal(y.S) * y.Vt; kwargs...)
end
function isapprox(
    ::FixedRankMatrices,
    x::SVDMPoint,
    v::UMVTVector,
    w::UMVTVector;
    kwargs...,
)
    return isapprox(
        x.U * v.M * x.Vt + v.U * x.Vt + x.U * v.Vt,
        x.U * w.M * x.Vt + w.U * x.Vt + x.U * w.Vt;
        kwargs...,
    )
end

@doc doc"""
    manifold_dimension(M::FixedRankMatrices{m,n,k,𝔽})

Return the manifold dimension for the `𝔽`-valued [`FixedRankMatrices`](@ref) `M`
of dimension `m`x`n` of rank `k`, namely

````math
k(m + n - k) \dim_ℝ 𝔽,
````

where $\dim_ℝ 𝔽$ is the [`real_dimension`](@ref) of `𝔽`.
"""
function manifold_dimension(::FixedRankMatrices{m,n,k,𝔽}) where {m,n,k,𝔽}
    return (m + n - k) * k * real_dimension(𝔽)
end

@doc doc"""
    project_tangent(M, x, A)
    project_tangent(M, x, v)

Project the matrix $A\in\mathbb R^{m,n}$ or a [`UMVTVector`](@ref) `v` from the embedding or
another tangent space onto the tangent space at $x$ on the [`FixedRankMatrices`](@ref) `M`,
further decomposing the result into $v=UMV$, i.e. a [`UMVTVector`](@ref) following
Section 3 in [^Vandereycken2013].

[^Vandereycken2013]:
    > Bart Vandereycken: "Low-rank matrix completion by Riemannian Optimization,
    > SIAM Journal on Optiomoization, 23(2), pp. 1214–1236, 2013.
    > doi: [10.1137/110845768](https://doi.org/10.1137/110845768),
    > arXiv: [1209.3834](https://arxiv.org/abs/1209.3834).
"""
project_tangent(::FixedRankMatrices, ::Any...)

function project_tangent!(
    ::FixedRankMatrices,
    vto::UMVTVector,
    x::SVDMPoint,
    A::AbstractMatrix,
)
    av = A * (x.Vt')
    uTav = x.U' * av
    aTu = A' * x.U
    vto.M .= uTav
    vto.U .= A * x.Vt' - x.U * uTav
    vto.Vt .= (aTu - x.Vt' * uTav')'
    return vto
end
function project_tangent!(
    M::FixedRankMatrices,
    vto::UMVTVector,
    x::SVDMPoint,
    v::UMVTVector,
)
    return project_tangent!(M, vto, x, v.U * v.M * v.Vt)
end

@doc doc"""
    representation_size(M::FixedRankMatrices{m,n,k})

Return the element size of a point on the [`FixedRankMatrices`](@ref) `M`, i.e.
the size of matrices on this manifold $(m,n)$.
"""
@generated representation_size(::FixedRankMatrices{m,n}) where {m,n} = (m, n)

@doc doc"""
    retract(M, x, v, ::PolarRetraction)

Compute an SVD-based retraction on the [`FixedRankMatrices`](@ref) `M` by computing
````math
    y = U_kS_kV_k^\mathrm{T},
````
where $U_k S_k V_k^\mathrm{T}$ is the shortened singular value decomposition $USV=x+v$,
in the sense that $S_k$ is the diagonal matrix of size $k\times k$ with the $k$ largest
singular values and $U$ and $V$ are shortened accordingly.
"""
retract(::FixedRankMatrices, ::Any, ::Any, ::PolarRetraction)

function retract!(
    ::FixedRankMatrices{M,N,k},
    y::SVDMPoint,
    x::SVDMPoint,
    v::UMVTVector,
    ::PolarRetraction,
) where {M,N,k}
    s = svd(x.U * Diagonal(x.S) * x.Vt + (x.U * v.M * x.Vt + v.U * x.Vt + v.U * v.Vt))
    y.U .= s.U[:, 1:k]
    y.S .= s.S[1:k]
    y.Vt .= s.Vt[1:k, :]
    return y
end

allocate(x::SVDMPoint) = SVDMPoint(allocate(x.U), allocate(x.S), allocate(x.Vt))
function allocate(x::SVDMPoint, ::Type{T}) where T
    return SVDMPoint(allocate(x.U,T), allocate(x.S,T), allocate(x.Vt,T))
end
allocate(v::UMVTVector) = UMVTVector(allocate(v.U), allocate(v.M), allocate(v.Vt))
function allocate(v::UMVTVector, ::Type{T}) where T
    return UMVTVector(allocate(v.U,T), allocate(v.M,T), allocate(v.Vt,T))
end

number_eltype(x::SVDMPoint) = typeof(one(eltype(x.U)) + one(eltype(x.S)) + one(eltype(x.Vt)))
number_eltype(v::UMVTVector) = typeof(one(eltype(v.U)) + one(eltype(v.M)) + one(eltype(v.Vt)))

one(x::SVDMPoint) = SVDMPoint(
    one(zeros(size(x.U, 1), size(x.U, 1))),
    ones(length(x.S)),
    one(zeros(size(x.Vt, 2), size(x.Vt, 2))),
    length(x.S),
)
one(v::UMVTVector) = UMVTVector(
    one(zeros(size(v.U, 1), size(v.U, 1))),
    one(zeros(size(v.M))),
    one(zeros(size(v.Vt, 2), size(v.Vt, 2))),
    size(v.M, 1),
)

function copyto!(x::SVDMPoint, y::SVDMPoint)
    copyto!(x.U, y.U)
    copyto!(x.S, y.S)
    copyto!(x.Vt, y.Vt)
    return x
end
function copyto!(v::UMVTVector, w::UMVTVector)
    copyto!(v.U, w.U)
    copyto!(v.M, w.M)
    copyto!(v.Vt, w.Vt)
    return v
end

@doc doc"""
    zero_tangent_vector(M::FixedRankMatrices, x::SVDMPoint)

Return a [`UMVTVector`](@ref) representing the zero tangent vector in the tangent space of
`x` on the [`FixedRankMatrices`](@ref) `M`, for example all three elements of the resulting
structure are zero matrices.
"""
function zero_tangent_vector(::FixedRankMatrices{m,n,k}, x::SVDMPoint) where {m,n,k}
    v = UMVTVector(
        zeros(eltype(x.U), m, k),
        zeros(eltype(x.S), k, k),
        zeros(eltype(x.Vt), k, n),
    )
    return v
end

function zero_tangent_vector!(
    ::FixedRankMatrices{m,n,k},
    v::UMVTVector,
    x::SVDMPoint,
) where {m,n,k}
    v.U .= zeros(eltype(v.U), m, k)
    v.M .= zeros(eltype(v.M), k, k)
    v.Vt .= zeros(eltype(v.Vt), k, n)
    return v
end


text:
function SeisMute(in;offset=[0.],tmute=0.,vmute=1500.,taper=0.1,dt=0.001)

	nt,nx = size(in[:,:])
	out = copy(in[:,:])
	for it = 1:nt
		for ix = 1:nx
			t = sqrt(tmute^2 + (offset[ix]/vmute).^2)
			if ((it-1)*dt <= t )
					if ((it-1)*dt > t - taper)
						out[it,ix] *= 1 - (t - (it-1)*dt)/taper
					else
						out[it,ix] = 0.
					end
			end
		end
	end

	return out
end

function SeisMute(in,h::Array{Header,1};offset=[0],tmute=0.,vmute=1500.,taper=0.1,dt=0.001)

	nt,nx = size(in)
	offset = Float32[]
	for itrace = 1 : nx
		#push!(offset,sqrt((h[itrace].gx - h[itrace].sx)^2 + (h[itrace].gy - h[itrace].sy)^2 ))
		push!(offset,abs(h[itrace].h))
	end
	out = SeisMute(in,offset=offset,tmute=tmute,vmute=vmute,taper=taper,dt=h[1].d1)

	return out,h
end

function SeisMute(m::String,d::String,adj;tmute=0.,vmute=1500.,taper=0.1)

	@compat parameters = Dict(:tmute=>tmute,:vmute=>vmute,:taper=>taper)
	if (adj==true)
		SeisProcess(d,m,[SeisMute],[parameters];key=["imx"])
	else
		SeisProcess(m,d,[SeisMute],[parameters];key=["imx"])
	end

end


text:
using SIMDPoly
using BenchmarkTools
using PyPlot

NSAMPLES = 10_000

degrees = 3:40
simd_times, ref_times = [], []

for n in degrees
    poly1 = ntuple(x -> Float64(x), n)
    poly2 = 2 .* poly1
    packed = packpoly2x2c(poly1, poly2)
    our = @benchmark evalpolydiv2x2c(z, $packed) samples=NSAMPLES setup=(z=rand()+rand()*im)
    ref = @benchmark (evalpoly(z, $poly1) / evalpoly(z, $poly2)) samples=NSAMPLES setup=(z=rand()+rand()*im)
    push!(simd_times, minimum(our).time)
    push!(ref_times, minimum(ref).time)
    #println(minimum(our).time, " / ", minimum(ref).time)
end

ioff()
plot(degrees, ref_times, linewidth=3, label="reference (Base.Math.evalpoly)")
plot(degrees, simd_times, linewidth=3, label="SIMD (4 x f64 lanes)")
legend()
ylim(ymin=0)
grid(axis="y", alpha=0.5)
xlabel("N")
ylabel("time (ns)")
title("degree N polynomial divison, complex input")
savefig("bm_plot_complex.svg")


text:
using RoadRunner
using Test

ant_str = """
    compartment c1, c2, c3, c4;
    c1 = 1.1; c2 = 2.2; c3 = 3.3; c4 = 4.4
"""

rr = RoadRunner.loada(ant_str)

@testset "compartment" begin
    @test RoadRunner.getNumberOfCompartments(rr) == 4
    @test RoadRunner.getCompartmentIds(rr) == ["c1", "c2", "c3", "c4"]
    @test RoadRunner.getCompartmentByIndex(rr, 0) == 1.1
    @test RoadRunner.getCompartmentByIndex(rr, 1) == 2.2
    @test RoadRunner.getCompartmentByIndex(rr, 2) == 3.3
    @test RoadRunner.getCompartmentByIndex(rr, 3) == 4.4
end


text:
# Autogenerated wrapper script for oneAPI_Level_Zero_Headers_jll for any
export ze_api

JLLWrappers.@generate_wrapper_header("oneAPI_Level_Zero_Headers")
JLLWrappers.@declare_file_product(ze_api)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_file_product(
        ze_api,
        "include/level_zero/ze_api.h",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#=
    Top level test suite, which runs unit tests, Pkg dependencies, regression & comparison tests.
=#

using Test

@testset "All tests" begin

    # Integration test
    include("test_c_sw.jl")

    #Interpolation test
    include("interpolate_test.jl")

    # Unit tests (in progress)
    # include("unit_tests.jl")

    # Test Pkg dependencies (not set up)
    # include("PkgTests.jl")

end


text:
PlotlyBase.savefig(p::SyncPlot, a...; k...) = savefig(p.plot, a...; k...)
PlotlyBase.savefig(io::IO, p::SyncPlot, a...; k...) = savefig(io, p.plot, a...; k...)

# TODO: overload mime methods once they are implemented via ORCA server


text:
#-------------------------------------------------------------------------------------------------
#--
#-- Create a range of default parameters
#--
using DataFrames, DataFramesMeta

function createDefaultParameters()
    # Create a vector containing triplets of (parameter name, initial value, range)
    params = [  ["r₀",            BaseR₀,         (2.0, 5.0)],
                ["tₗ",            4,              (3.0, 6.0)],
                ["tᵢ",            4,              (3.0, 6.0)],
                ["tₕ",            4,              (4.0, 4.0)],
                ["tᵤ",            14,             (12.0, 17.0)],
                ["γₑ",            0.8,            (0.2, 1.0)],
                ["γᵢ",            1.0,            (1.0, 1.0)],
                ["γⱼ",            1.5,            (0.2, 2.0)],
                ["γₖ",            2.0,            (0.2, 3.0)],
                ["δₖ",            1.0,            (0.9, 1.1)],
                ["δₗ",            2.0,            (0.8, 3.0)],
                ["δᵤ",            1.0,            (1.0, 1.0)],
                ["infectedM",     10.0,           (1.0, 50.0)],
                ["infectiousM",   20.0,           (1.0, 50.0)],

                # Mitigation profile with 4 points
                ["mv0",           1.0,            (1.0, 1.0)],
                ["mv1",           0.8,            (0.1, 1.5)],
                ["mv2",           0.5,            (0.1, 1.5)],
                ["mv3",           0.5,            (0.1, 1.5)],
                ["mv4",           0.5,            (0.1, 1.5)],
                ["mv5",           0.5,            (0.1, 1.5)],
                ["mv6",           0.5,            (0.1, 1.5)],
                ["mv7",           0.5,            (0.1, 1.5)],
                ["mv8",           0.5,            (0.1, 1.5)],
                ["mv9",           0.5,            (0.1, 1.5)]]

                # Fixed parameters
#                ["BED_max",       1000.0,         [1000.0, 1000.0]],
#                ["ICU_max",       1000.0,         [1000.0, 1000.0]]]

    nParams = size(params)[1]
    pNames = [params[i][1] for i in 1:nParams]
    initialParam = [params[i][2] for i in 1:nParams]
    rangeParam = Tuple{Float64, Float64}[params[i][3] for i in 1:nParams]

    return(initialParam, rangeParam, pNames)

end

#-------------------------------------------------------------------------------------------------
#--
#-- Define masks indicating variables relating to the disease and to the countries
#-- Note that 2 variables are never optimised to provide an anchor for the others
#-- They are the infectiousness of a symptomatic individual γᵢ and the fatality rate from ICU δᵤ
#--
diseaseMask = [1,
               1, 1, 1, 1,     # t parameters
               1, 0, 1, 1,     # γ parameters
               1, 1, 0,        # δ parameters
               0, 0,
               0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

countryMask = ones(length(diseaseMask)) .- diseaseMask

#-------------------------------------------------------------------------------------------------
#--
#-- Freeze the range of parameters specific to a given country to flex the disease-specific ones
#--
function freezeCountryRange()
    # List of parameters to freeze
    listParams = ["γᵢ", "δᵤ",
                  "infectedM", "infectiousM",
                  "mv0", "mv1", "mv2", "mv3", "mv4",
                  "mv5", "mv6", "mv7", "mv8", "mv9"]

    # Get the names of the parameters
    _, _, pNames = createDefaultParameters()

    for (c, _) in COUNTRY_LIST

        # for each parameter
        for position in 1:length(pNames)

            value = countryData[c][:params][position]

            # If the name is to be relaxed (not in the list)
            if findfirst(pNames[position] .== listParams) == nothing
                # Sets the range at the value +/- 10%
                countryData[c][:range][position] = (value * (1-0.10), value * (1+0.10))
            else
                # Otherwise, sets the optimisation range at that value
                countryData[c][:range][position] = (value, value)
            end
        end
    end
end

function freezeDiseaseRange()
#-- This function is the same code as freezeCountryRange, but different list

    # List of parameters to freeze
    listParams = ["r₀", "tₗ", "tᵢ", "tₕ", "tᵤ",
                  "γₑ", "γⱼ", "γₖ",
                  "δₖ", "δₗ", ]

    # Get the names of the parameters
    _, _, pNames = createDefaultParameters()

    for (c, _) in COUNTRY_LIST

        # for each parameter
        for position in 1:length(pNames)

            value = countryData[c][:params][position]

            # If the name is to be relaxed (not in the list)
            if findfirst(pNames[position] .== listParams) == nothing
                # Sets the range at the value +/- 10%
                countryData[c][:range][position] = (value * (1-0.10), value * (1+0.10))
            else
                # Otherwise, sets the optimisation range at that value
                countryData[c][:range][position] = (value, value)
            end
        end
    end
end

function deathsLoss(p)
    sol = calculateSolution(country, p)

    # Extract total deaths profile
    forecastDeaths = calculateTotalDeaths(sol)
    #println(sum(forecastDeaths)); println();

    return sqrt(sum( (countryData[country][:cases].deaths .- forecastDeaths) .^ 2 ))
end

#--
#-- Calculate the sum of all the losses of all the countries to optimise disease params.
#--
function allCountriesLoss(p)

    # The parameters passed to the individual loss is created with a mask defined
    totalLoss = 0.0

    for (country, _) in COUNTRY_LIST

        countryParams = countryData[country][:params]
        totalLoss += deathsLoss((diseaseMask .* p) .+ (countryMask .* countryParams))
    end

    return totalLoss
end

function calculateSolution(country, params)

    # First date should the date of the last death reported
    startDate = first(countryData[country][:cases].time)
    startDays = first(countryData[country][:cases].t)

    # Final date should the date of the last death reported
    endDate = last(countryData[country][:cases].time)
    endDays = last(countryData[country][:cases].t)

    tSpan = (startDays, endDays)

    ## Country-specific constants
    # Those are constants which cannot be changed to improve the model.

    BED_max = countryData[country][:hospital_capacity]
    ICU_max = countryData[country][:ICU_capacity]

    Age_Pyramid = transpose( Matrix(countryData[country][:age_distribution]))
    Age_Pyramid_frac = Age_Pyramid / sum(Age_Pyramid)

    TotalConfirmedAtStart = @where(countryData[country][:cases], :time .== startDate)[!, :cases][1]
    ConfirmedAtStart = TotalConfirmedAtStart .* Age_Pyramid_frac

    TotalDeathsAtStart = @where(countryData[country][:cases], :time .== startDate)[!, :deaths][1]
    DeathsAtStart = TotalDeathsAtStart .* Age_Pyramid_frac

    # Deconstruct the parameters
    r₀, tₗ, tᵢ, tₕ, tᵤ, γₑ, γᵢ, γⱼ, γₖ, δₖ, δₗ, δᵤ, infectedM, infectiousM,
        mv0, mv1, mv2, mv3, mv4, mv5, mv6, mv7, mv8, mv9 = params

    mitigation = [(0, mv0), (7, mv1), (14, mv2), (21, mv3), (28, mv4),
                  (35, mv5), (42, mv6), (49, mv7), (56, mv8), (63, mv9)]

    ## Parameter vector
    # Those are parameters which can be changed to improve the model
    TotalInfected = infectedM * TotalConfirmedAtStart
    InfectedAtStart = TotalInfected .* Age_Pyramid_frac

    TotalInfectious = infectiousM .* TotalConfirmedAtStart
    InfectiousAtStart = TotalInfectious .* Age_Pyramid_frac

    model_params = [r₀,
                  tₗ, tᵢ, tₕ, tᵤ,
                  γₑ, γᵢ, γⱼ, γₖ,
                  δₖ, δₗ, δᵤ,
                  mitigation,
                  BED_max,
                  ICU_max]

    ## Compartment vector
    # Note that values are initialised at 1 to avoid division by zero
    S0 = Age_Pyramid .- InfectedAtStart .- InfectiousAtStart .- DeathsAtStart
    E0 = InfectedAtStart
    I0 = InfectiousAtStart
    J0 = ones(nAgeGroup)
    H0 = ones(nAgeGroup)
    C0 = ones(nAgeGroup)
    R0 = ones(nAgeGroup)
    D0 = DeathsAtStart
    K0 = ones(nAgeGroup)
    L0 = ones(nAgeGroup)

    # Everybody confirmed is in hospital. Assume 1 ICU bed to stay away from zeros.
    BED = [TotalInfectious]
    ICU = [1.0]

    P0 = vcat(S0, E0, I0, J0, H0, C0, R0, D0, K0, L0, BED, ICU)

    # Differential equation solver
    model = ODEProblem(epiDynamics!, P0, tSpan, model_params)

    # Note: progress steps might be too quick to see!
    sol = solve(model, Tsit5(); progress = false)

    return sol
end


text:
# 1 ellipsoid on a plane

#push!(LOAD_PATH, "/Users/tomek/Dropbox/tkwork/codes/edyna")
#import edyna
include("../edyna.jl")

edyna.ellipsoid([0, 0, 0.1], 0.01, 0.02, 0.03, [1. 0. 0.; 0. 1. 0.; 0. 0. 1.], 1000)

edyna.obstacle([[[-1, -1, 0], [-1, 1, 0], [1, -1, 0]], [[1, -1, 0], [-1, 1, 0], [1, 1, 0]]])

edyna.gravity(0, 0, -10)

edyna.run_rigid(1.0, 0.001, 0.01, "out/1_on_plane")

max_stars_repo_name: tkoziara/edyna

max_issues_repo_name: tkoziara/edyna

max_forks_repo_name: tkoziara/edyna


text:
using SeaPearl
using Test
using LightGraphs
using Flux
using Zygote
using Random
using DataStructures
using ReinforcementLearning
const RL = ReinforcementLearning

@testset "SeaPearl.jl" begin
    include("CP/CP.jl")
     #RL tests are launched in CP/valueselection
    #include("MOI_wrapper/MOI_wrapper.jl")
    include("trailer.jl")
    include("datagen/datagen.jl")
    include("experiment/experiment.jl")
end

revise_user = "You use Revise, you're efficient in your work, well done ;)"


text:
using Base.Test
using Ev3dev

include("mock.jl")

let
    mktempdir() do dir
        @show dir
        mock_robot(dir, [(:large_motor, "outA"),
                         (:large_motor, "outB"),
                         (:medium_motor, "outC"),
                         (:ultrasound_sensor, "inA")])
        run(`ls $(dir)`)

        brick = Brick(dir)
        large_motors = find_devices(LargeMotor, brick)
        @test length(large_motors) == 2

        medium_motors = find_devices(MediumMotor, brick)
        @test length(medium_motors) == 1

        right_motor = find_device_at_address(LargeMotor, brick, "outA")
        @test right_motor.io.address() == "outA"

        ultrasound_sensor = find_devices(UltrasoundSensor, brick)[1]

        @test scaled_values(ultrasound_sensor) == [65.0]
    end
end


text:
#=
KTrussExample:
- Julia version: 1.6.2
- Author: samuel
- Date: 2021-09-12
=#

using SuiteSparseGraphBLAS
using SparseArrays
using LinearAlgebra
using Test

include("/Users/samuel/Documents/Jupyter/GraphBlasAlgorithmsFinal/src/KTruss.jl")
dim = 10000
matrix = GBMatrix(sparse(rand(0:1, dim, dim)))
@time ktruss(matrix, 4, dim)