text:
# #----------------------------------------------------------------------------------------------------------------------
# #----------------------------------------------------------------------------------------------------------------------
# This file contains functions and other snippets of code that are used in various calculations.
# #----------------------------------------------------------------------------------------------------------------------
# #----------------------------------------------------------------------------------------------------------------------

#######################################################################################################################
# CALCULATE REGIONAL CO₂ MITIGATION.
########################################################################################################################
# Description: This function calculates regional CO₂ mitigation levels as a function of a global carbon tax. It
#              uses the RICE2010 backstop price values and assumes a carbon tax of $0 in period 1.  If the number of
#              tax values is less than the total number of model time periods, the function assumes full decarbonization
#              (e.g. the tax = the backstop price) for all future periods without a specified tax.
#
# Function Arguments:
#
#       optimal_tax:    A vector of global carbon tax values to be optimized.
#       backstop_price: The regional backstop prices from RICE2010 (units must be in $1000s)
#       theta2:         The exponent on the abatement cost function (defaults to RICE2010 value).
#
# Function Output:
#
#       mitigation:     Regional mitigation rates resulting from the global carbon tax.
#----------------------------------------------------------------------------------------------------------------------

function mitigation_from_tax(optimal_tax::Array{Float64,1}, backstop_prices::Array{Float64,2}, theta2::Float64)

    # Initialize full tax vector with $0 tax in period 1 and maximum of the backstop price across all regions for remaining periods.
    full_tax = [0.0; maximum(backstop_prices, dims=2)[2:end]]

    # Set the periods being optimized to the optimized tax value (assuming full decarbonization for periods after optimization time frame).
    full_tax[2:(length(optimal_tax)+1)] = optimal_tax

    # Calculate regional mitigation rates from the full tax vector.
    mitigation = min.((max.(((full_tax ./ backstop_prices) .^ (1 / (theta2 - 1.0))), 0.0)), 1.0)

    return mitigation
end

#######################################################################################################################
# CREATE RICE OBJECTIVE FUNCTION.
########################################################################################################################
# Description: This function creates an objective function and an instance of RICE with user-specified parameter settings.
#              The objective function will take in a vector of global carbon tax values (cost-minimization) or regional
#              CO₂ mitigation rates (utilitarian) and returns the total economic welfare generated by that specifc
#              climate policy.
#
# Function Arguments:
#
#       run_utilitarian: A true/false indicator for whether or not to run the utilitarian optimization (true = run utilitarian).
#       ρ:                  Pure rate of time preference.
#       η:                  Elasticity of marginal utility of consumption.
#       backstop_prices:    The regional backstop prices from RICE2010 (units must be in dollars).
#       remove_negishi:     A true/false indicator for whether RICE should use a social welfare function with Negishi weights (true = remove Negishi weights).
#                           *Note: if using Negishi weights, RICE's discounting parameters default to η=1.5 and ρ=1.5%.
#
# Function Output:
#
#       rice_objective:     The objective function specific to user model settings.
#       m:                  An instance of RICE2010 consistent with user model settings.
#----------------------------------------------------------------------------------------------------------------------

function construct_rice_objective(run_utilitarian::Bool, ρ::Float64, η::Float64, backstop_prices::Array{Float64,2}, remove_negishi::Bool)

    # Get an instance of RICE given user settings.
    m = create_rice(ρ, η, remove_negishi)

    #--------------------------------------------------------------------------------------------------------
    # Create either a (i) cost-minimization or (ii) utilitarian objective function for this instance of RICE.
    #--------------------------------------------------------------------------------------------------------
    rice_objective = if run_utilitarian == false

        #---------------------------------------
        # Cost-minimzation price objective function.
        #---------------------------------------
        function(optimal_global_tax::Array{Float64,1})
            # Set the regional mitigation rates to the value implied by the global optimal carbon tax and return total welfare.
            update_param!(m, :MIU, mitigation_from_tax(optimal_global_tax, backstop_prices, 2.8))
            run(m)
            return m[:welfare, :UTILITY]
        end

    else

        #---------------------------------------
        # Utilitarian objective function.
        #---------------------------------------
        function(optimal_mitigation_vector::Array{Float64,1})
            # Number of periods with optimized rates (optimization program requires a vector).
            n_opt_periods = Int(length(optimal_mitigation_vector) / 12)
            # NOTE: for convenience, this objective directly optimizes the regional mitigation rates.
            # Initialze a regional mitigation array, assuming first period = 0% mitigation and periods after optimization achieve full decarbonization.
            optimal_regional_mitigation = vcat(zeros(1,12), ones(59,12))
            optimal_regional_mitigation[2:(n_opt_periods + 1), :] = reshape(optimal_mitigation_vector, (n_opt_periods, 12))
            # Set the optimal regional mitigation rates and return total welfare.
            update_param!(m, :MIU, optimal_regional_mitigation)
            run(m)
            return m[:welfare, :UTILITY]
        end

    end

    # Return the newly created objective function and the specific instance of RICE.
    return rice_objective, m
end

#######################################################################################################################
# OPTIMIZE RICE.
########################################################################################################################
# Description: This function takes an objective function (given user-supplied model settings), and optimizes it for the
#              cost-minimization (global carbon taxes) or utilitarian (regional carbon taxes) approach to find the
#              policy that maximizes global economic welfare. Note that the utilitarian objective function optimizes
#              on the decarbonization fraction (for more efficient code) and then calculates the corresponding regional
#              carbon tax values.
#
# Function Arguments:
#
#       optimization_algorithm:  The optimization algorithm to use from the NLopt package.
#       n_opt_periods:           The number of model time periods to optimize over.
#       stop_time:               The length of time (in seconds) for the optimization to run in case things do not converge.
#       tolerance:               Relative tolerance criteria for convergence (will stop if |Δf| / |f| < tolerance from one iteration to the next.)
#       backstop_price:          The regional backstop prices from RICE2010 (units must be in dollars).
#       run_utilitarian:      A true/false indicator for whether or not to run the utilitarian optimization (true = run utilitarian).
#       ρ:                       Pure rate of time preference.
#       η:                       Elasticity of marginal utility of consumption.
#       remove_negishi:          A true/false indicator for whether RICE should use a social welfare function with Negishi weights (true = remove Negishi weights).
#
# Function Output
#
#       optimized_policy_vector: The vector of optimized policy values returned by the optimization algorithm.
#       optimal_emissions:       Optimal emissions for all time periods.
#       optimal_mitigation:      Optimal mitigation rates for all time periods (2005-2595) resulting form the optimization.
#       optimal_tax:             The optimal tax for all periods and regions used to run the model.
#       opt_model:               An instance of RICE (with user-defined settings) set with the optimal CO₂ mitigation policy.
#       convergence_result:      Indicator for whether or not the optimization converged.
#----------------------------------------------------------------------------------------------------------------------

function optimize_rice(optimization_algorithm::Symbol, n_opt_periods::Int, stop_time::Int, tolerance::Float64, backstop_prices::Array{Float64,2}; run_utilitarian::Bool=true, ρ::Float64=0.008, η::Float64=1.5, remove_negishi::Bool=true)

    # -------------------------------------------------------------
    # Create objective function and values needed for optimization.
    #--------------------------------------------------------------

    # Create objective function and instance of RICE, given user settings.
    objective_function, optimal_model = construct_rice_objective(run_utilitarian, ρ, η, backstop_prices, remove_negishi)

    # Set objective function, upper bound, and number of optimzation objectives (will differ between cost-minimization and utilitarian approaches).
    if run_utilitarian == false
        # Number of objectives is equal to time periods being optimized.
        n_objectives = n_opt_periods
        # Upper bound is maximum of backstop price, across all regions.
        upper_bound = maximum(backstop_prices, dims=2)[2:(n_objectives+1)]
    else
        # Number of objectives is equal to time periods being optimizer × 12 regions.
        n_objectives = n_opt_periods * 12
        # Upper bound is 1.0 because we are optimizing the decarbonization rate (1 = full mitigation).
        upper_bound = ones(n_objectives)
    end

    # Create lower bound.
    lower_bound = zeros(n_objectives)

    # Create initial condition for algorithm (set at 50% of upper bound).
    starting_point = upper_bound/2

    # -------------------------------------------------------
    # Create an NLopt optimization object and optimize model.
    #--------------------------------------------------------
    opt = Opt(optimization_algorithm, n_objectives)

    # Set the bounds.
    lower_bounds!(opt, lower_bound)
    upper_bounds!(opt, upper_bound)

    # Assign the objective function to maximize.
    max_objective!(opt, (x, grad) -> objective_function(x))

    # Set termination time.
    maxtime!(opt, stop_time)

    # Set optimizatoin tolerance (will stop if |Δf| / |f| < tolerance from one iteration to the next).
    ftol_rel!(opt, tolerance)

    # Optimize model.
    maximum_objective_value, optimized_policy_vector, convergence_result = optimize(opt, starting_point)

    # Create optimal decarbonization rates for all time periods (approach to do so will differ between cost-minimization and utilitarian).
    if run_utilitarian == false
        optimal_mitigation = mitigation_from_tax(optimized_policy_vector, backstop_prices, 2.8)
    else
        optimal_mitigation = vcat(zeros(1,12), ones(59,12))
        optimal_mitigation[2:(n_opt_periods+1), :] = reshape(optimized_policy_vector, (n_opt_periods, 12))
    end

    # Run user-specified version of RICE with optimal mitigation policy.
    update_param!(optimal_model, :MIU, optimal_mitigation)
    run(optimal_model)

    # Create optimal tax rates for all time periods (approach to do so will differ between cost-minimization and utilitarian).
    if run_utilitarian == false
        optimal_tax = [0.0; maximum(backstop_prices, dims=2)[2:end]]
        optimal_tax[2:(length(optimized_policy_vector)+1)] = optimized_policy_vector
    else
        optimal_tax = optimal_model[:emissions, :CPRICE]
    end

    # Create optimal industiral emissions for all time periods.
    optimal_emissions = optimal_model[:emissions, :EIND]

    # Return results of optimization, optimal emissions, optimal mitigation rates, optimal taxes, and RICE run with optimal mitigation policies.
    return optimized_policy_vector, optimal_emissions, optimal_mitigation, optimal_tax, optimal_model, convergence_result
end

#######################################################################################################################
# LINEAR INTERPOLATION FUNCTION
########################################################################################################################
# Description: This function will carry out a linear interpolation between optimal tax rates for a given time spread to
#              produce a series of annual carbon tax values.
#
# Function Arguments:
#
#       years_with_data: Years along the model time horizon with optimized tax values (will interpolate between these).
#       data:            The data points to be interpolated between (in this case a vector of tax values).
#       spacing:         The number of years between each tax value.
#
# Function Output:
#
#       results:         A vector of interpolated values based on the original input data.
#----------------------------------------------------------------------------------------------------------------------

function interpolate(years_with_data::Array{Int,1}, data::Array{Float64,1}, spacing::Int)
    # Create a vector of the annual years to interpolate across.
    annual_years = collect(years_with_data[1]:years_with_data[end])
    # Create a vector for results and set last value (for convenience).
    results = zeros(length(annual_years))
    results[end] = data[end]
    # Loop through, linearly interpolating between each set of periods with data.
    for index in 1:(length(years_with_data)-1)
        for t = (spacing*(index-1)+1):(spacing*index)
            results[t] = (annual_years[t] - years_with_data[index])/spacing * data[index+1] + (years_with_data[index+1]-annual_years[t])/spacing*data[index]
        end
    end
    # Return interpolated results.
    return results
end

#######################################################################################################################
# CALCULATE INDEX FOR A GIVEN FUND YEAR
########################################################################################################################
# Description: Just for convenience, this function gives the index for a given year within the FUND model time horizon.
#
# Function Arguments:
#
#       index_year: The year to get an index for.
#
# Function Output:
#
#        output:    An index for the supplied year (within the time frame 1950 - final year).
#----------------------------------------------------------------------------------------------------------------------

function fund_year_index(index_year::Int)
    return index_year - 1950 + 1
end

#######################################################################################################################
# CREATE AN ANNUAL GLOBAL CARBON TAX ARRAY
########################################################################################################################
# Description: This function will take a vector of carbon tax values and create a linear interpolation across the extent
#              of the FUND model time horizon. It assumes an optimized tax value every 10 years from 2010-2200 (with linear
#              interpolation between). After 2200, the tax rate stays constant until 2300. The function will replicate
#              this vector for each FUND region so they all face the same carbon tax.
#
# Function Arguments:
#
#       tax:               A vector of tax values.
#       welfare_year:      The year to start calculating how a climate policy affects total welfare (FUND has a historic spin up period).
#       end_year:          The last year to run FUND.
#
# Function Output:
#
#       global_tax_vector: An interpolated tax vector for all time periods, replicated 16 times for each FUND region.
#----------------------------------------------------------------------------------------------------------------------

function create_global_tax_fund(tax::Array{Float64,1}, welfare_year::Int, end_year::Int)
    # Initialize global tax vector for entire length of FUND time horizon.
    global_tax_vector = zeros(length(1950:end_year))
    # Linearly interpolate between values from 2010-2200 (i.e. 20 tax values).
    global_tax_vector[fund_year_index(welfare_year):fund_year_index(2200)] = interpolate(collect(welfare_year:10:2200), tax[1:20], 10)
    # Assume tax rate remains constant after 2200.
    global_tax_vector[fund_year_index(2201):end] .= global_tax_vector[fund_year_index(2200)]
    # Repeat vector 16 times for each region.
    return repeat(global_tax_vector, outer=[1, 16])
end

#######################################################################################################################
# CREATE AN ANNUAL REGIONAL CARBON TAX ARRAY
########################################################################################################################
# Description: This function will take a vector of carbon tax values specific to each FUND region, and create a linear
#              interpolation across the extent of the FUND model time horizon. For each region, it assumes an optimized
#              tax value every 10 years from 2010-2200 (with linear interpolation between). After 2200, the tax rate stays
#              constant until 2300. The function will return an array, where each column is a unique regional tax series.
#
# Function Arguments:
#
#       tax:                 A vector of tax values for all regions.
#       welfare_year:        The year to start calculating how a climate policy affects total welfare (FUND has a historic spin up period).
#       end_year:            The last year to run FUND.
#
# Function Output:
#
#       regional_tax_matrix: A matrix of interpolated tax values for all time periods, with each column representing a FUND region.
#----------------------------------------------------------------------------------------------------------------------

function create_regional_tax_fund(tax::Array{Float64,1}, welfare_year::Int, end_year::Int)
    # Initialize regional tax array for entire length of FUND time horizon.
    regional_tax_matrix = zeros(length(1950:end_year), 16)
    # Reshape optimal tax vector to be a matrix (just makes indexing easier). Each row = a decade's tax value, each column = a region.
    opt_tax = reshape(tax, (20,16))

    # Loop through and create annual tax values for each region.
    for r = 1:16
        regional_tax_matrix[fund_year_index(welfare_year):fund_year_index(2200), r] = interpolate(collect(welfare_year:10:2200), opt_tax[1:20,r], 10)
        # Set all values after 2200 to a constant tax value.
        regional_tax_matrix[fund_year_index(2201):end, r] .= regional_tax_matrix[fund_year_index(2200), r]
    end

    return regional_tax_matrix
end

#######################################################################################################################
# CREATE REGIONAL CO₂ MITIGATION RATES GIVEN A CARBON TAX
########################################################################################################################
# Description: This function will take an instance of FUND run with an optimal carbon tax policy and return the resulting
#              CO₂ mitigation rates for each region. Note that a carbon tax in FUND has both a temporary and permanent
#              effect on emission reductions. The mitigation rates calculated here represent the percentage difference in
#              emissions for a given year relative to baseline case with no carbon tax or climate policy.
#
# Function Arguments:
#
#       optimal_model:    An instance of Mimi-FUND that has been run with an optimal carbon tax policy.
#       end_year:         The last year to run FUND.
#
# Function Output:
#
#       mitigation_rates: An array of CO₂ emission reduction rates caused by the carbon tax for each FUND region.
#----------------------------------------------------------------------------------------------------------------------

function regional_mitigation(optimal_model::Model, end_year::Int)

    # Get baseline version of FUND with no carbon tax but same parameter settings as optimal model.
    base_model = deepcopy(optimal_model)
    update_param!(base_model, :currtax, zeros(length(1950:end_year), 16))
    run(base_model)

    # Initialze an array to store regional mitigation rates.
    mitigation_rates = zeros(length(1950:end_year), 16)

    # Loop through each region and calculate mitigation relative to base case with no policy.
    # Note, first FUND period (1950) is a 'missing' value, so skip that.
    for r = 1:16
        mitigation_rates[2:end,r] = (base_model[:emissions, :emissionwithforestry][2:end,r] .- optimal_model[:emissions, :emissionwithforestry][2:end,r]) ./ base_model[:emissions, :emissionwithforestry][2:end,r]
    end

    return mitigation_rates
end

#######################################################################################################################
# CREATE FUND OBJECTIVE FUNCTION.
########################################################################################################################
# Description: This function creates an objective function and an instance of FUND with user-specified parameter settings.
#              The objective function will take in a vector of global or regional (utilitarian)
#              carbon tax values and returns the total economic welfare generated by that specifc climate policy.
#
# Function Arguments:
#
#       run_utilitarian: A true/false indicator for whether or not to run the utilitarian optimization (true = run utilitarian).
#       ρ:                  Pure rate of time preference.
#       η:                  Elasticity of marginal utility of consumption.
#       welfare_year:       The year to start calculating how a climate policy affects total welfare (FUND has a historic spin up period).
#       end_year:           The last year to run FUND.
#
# Function Output:
#
#       fund_objective:     The objective function specific to user model settings.
#       m:                  An instance of FUND consistent with user model settings.
#----------------------------------------------------------------------------------------------------------------------

function construct_fund_objective(run_utilitarian::Bool, ρ::Float64, η::Float64, welfare_year::Int, end_year::Int)

    # Get an instance of FUND given user settings.
    m = create_fund(ρ, η, welfare_year)

    #--------------------------------------------------------------------------------------------------------
    # Create either a (i) cost-minimization price or (ii) utilitarian objective function for this instance of FUND.
    #--------------------------------------------------------------------------------------------------------
    fund_objective = if run_utilitarian == false

        #---------------------------------------
        # Cost-minimization Price objective function.
        #---------------------------------------
        function(optimal_global_tax::Array{Float64,1})
            # Apply a global carbon tax to each reagion, run the model, and return total welfare.
            update_param!(m, :currtax, create_global_tax_fund(optimal_global_tax, welfare_year, end_year))
            run(m)
            return m[:fund_welfare, :UTILITY][end]
        end

    else

        #---------------------------------------
        # Utilitarian objective function.
        #---------------------------------------
        function regional_tax_fund_objective(regional_tax_vector::Array{Float64,1})
            # Create an array of unique optimal tax paths for each FUND region
            update_param!(m, :currtax, create_regional_tax_fund(regional_tax_vector, welfare_year, end_year))
            run(m)
            return m[:fund_welfare, :UTILITY][end]
        end
    end

    # Return the newly created objective function and the specific instance of RICE.
    return fund_objective, m
end

#######################################################################################################################
# OPTIMIZE FUND.
########################################################################################################################
# Description: This function takes an objective function (given user-supplied model settings), and optimizes it for the
#              cost-minimization (global carbon taxes) or utilitarian (regional carbon taxes) approach to find the
#              policy that maximizes global economic welfare.
#
# Function Arguments:
#
#       optimization_algorithm:  The optimization algorithm to use from the NLopt package.
#       n_opt_periods:           The number of model time periods to optimize over.
#       stop_time:               The length of time (in seconds) for the optimization to run in case things do not converge.
#       tolerance:               Relative tolerance criteria for convergence (will stop if |Δf| / |f| < tolerance from one iteration to the next.)
#       run_utilitarian:      A true/false indicator for whether or not to run the utilitarian optimization (true = run utilitarian).
#       ρ:                       Pure rate of time preference.
#       η:                       Elasticity of marginal utility of consumption.
#       welfare_year:            The year to start calculating how a climate policy affects total welfare (FUND has a historic spin up period).
#       end_year:                The last year to run FUND.
#
# Function Output
#
#       optimized_policy_vector: The vector of optimized policy values returned by the optimization algorithm.
#       optimal_mitigation:      Optimal mitigation rates for all time periods (2005-2595) resulting form the optimization.
#       optimal_tax:             The optimal tax for all periods and regions used to run the model.
#       opt_model:               An instance of FUND (with user-defined settings) set with the optimal CO₂ mitigation policy.
#       convergence_result:      Indicator for whether or not the optimization converged.
#----------------------------------------------------------------------------------------------------------------------

function optimize_fund(optimization_algorithm::Symbol, n_opt_periods::Int, stop_time::Int, tolerance::Float64; run_utilitarian::Bool=true, ρ::Float64=0.008, η::Float64=1.5, welfare_year::Int=2010, end_year::Int=2300)

    # -------------------------------------------------------------
    # Create objective function and values needed for optimization.
    #--------------------------------------------------------------

    # Create objective function and instance of FUND, given user settings.
    objective_function, optimal_model = construct_fund_objective(run_utilitarian, ρ, η, welfare_year, end_year)

    # Set number of optimzation objectives (will differ between cost-minimization and utilitarian approaches).
    if run_utilitarian == false
        # Number of objectives is equal to the number of tax values (one for every 10 years) being optimized.
        n_objectives = n_opt_periods
    else
        # Number of objectives is equal to number of tax values being optimized × 16 regions.
        n_objectives = n_opt_periods * 16
    end

    # Create lower and upper bounds (assume tax cannot go above $5000/ton CO₂).
    lower_bound = zeros(n_objectives)
    upper_bound = ones(n_objectives) .* 5000.0

    # Create initial condition for algorithm (start tax at $50 for all time periods and regions).
    starting_point = ones(n_objectives) .* 50.0

    # -------------------------------------------------------
    # Create an NLopt optimization object and optimize model.
    #--------------------------------------------------------
    opt = Opt(optimization_algorithm, n_objectives)

    # Set the bounds.
    lower_bounds!(opt, lower_bound)
    upper_bounds!(opt, upper_bound)

    # Assign the objective function to maximize.
    max_objective!(opt, (x, grad) -> objective_function(x))

    # Set termination time.
    maxtime!(opt, stop_time)

    # Set optimizatoin tolerance (will stop if |Δf| / |f| < tolerance from one iteration to the next).
    ftol_rel!(opt, tolerance)

    # Optimize model.
    maximum_objective_value, optimized_policy_vector, convergence_result = optimize(opt, starting_point)

    # Calculate the full tax series for all time periods and regions (approach to do so will differ between cost-minimization and utilitarian).
    if run_utilitarian == false
        optimal_tax = create_global_tax_fund(optimized_policy_vector, welfare_year, end_year)
    else
        optimal_tax = create_regional_tax_fund(optimized_policy_vector, welfare_year, end_year)
    end

    # Run user-specified version of FUND with optimal carbon tax.
    update_param!(optimal_model, :currtax, optimal_tax)
    run(optimal_model)

    # Create optimal regional decarbonization rates for all time periods resulting from the carbon tax.
    optimal_mitigation = regional_mitigation(optimal_model, end_year)

    # Return results of optimization, optimal mitigation rates, optimal taxes, and FUND run with optimal mitigation policies.
    return optimized_policy_vector, optimal_mitigation, optimal_tax, optimal_model, convergence_result
end

max_stars_repo_name: Environment-Research/Utilitarianism

max_issues_repo_name: Environment-Research/Utilitarianism

max_forks_repo_name: Environment-Research/Utilitarianism


text:
using .SearchProblem: AbstractSearchProblem, start_state, heuristic, cost, goal_test
using .GraphSearch
using .GraphSearch: AbstractTreeSearchAlgorithm, graph_search!, get_solution
using Random
using Statistics

abstract type AbstractWeightAdjustmentPolicy end
possible_weights(wap::AbstractWeightAdjustmentPolicy) = error("not implemented")
(wap::AbstractWeightAdjustmentPolicy)(;kwargs...)::Float64 = error("Not implemented")

struct ConstantWeightPolicy <: AbstractWeightAdjustmentPolicy
    weight::Float64
end
@inline possible_weights(cwp::ConstantWeightPolicy) = [cwp.weight]
@inline (cwp::ConstantWeightPolicy)(;kwargs...)::Float64 = cwp.weight

mutable struct AWAStar{S, A, WAP <: AbstractWeightAdjustmentPolicy} <: AbstractTreeSearchAlgorithm{S,A}
    weight::Float64
    walltime_limit::Float64
    nodes_budget::Int
    weight_adjustment_policy::WAP

    open_lists::MultipleOpenLists{S,A}
    closed_set::Set{S}
    path_costs::Dict{S, Float64}
    start_time::Float64
    stop_time::Union{Nothing,Float64}
    nodes_expended::Int
    interrupted::Bool

    nodes_expended_to_fist_solution::Union{Int,Nothing}
    num_solutions::Int
    all_solutions_costs::Vector{Float64}
    all_solutions_walltimes::Vector{Float64}
    all_solutions_nodes_expended::Vector{Int}
    all_solutions::Vector{Vector{A}}
    solution::Union{Vector{A},Nothing}
    solution_cost::Float64
    solution_walltime::Union{Float64,Nothing}
    solution_nodes_expended::Union{Int,Nothing}
    solution_w_trace::Union{Vector{Float64},Nothing}

    function AWAStar{S,A}(weight::Real, walltime_limit::Real, nodes_budget::Integer, wap::WAP) where {S,A,WAP<:AbstractWeightAdjustmentPolicy}
        awastar = new{S,A,WAP}(weight, walltime_limit, nodes_budget, wap)
        awastar.open_lists = MultipleOpenLists{S,A}(possible_weights(wap))
        awastar.closed_set = Set{S}()
        awastar.path_costs = Dict{S,Float64}()
        awastar.all_solutions_costs = Float64[]
        awastar.all_solutions_walltimes = Float64[]
        awastar.all_solutions_nodes_expended = Int[]
        awastar.all_solutions = Vector{A}[]
        return awastar
    end
end

@inline simulated_time(awastar::AWAStar, simulated_expansion_rate::Real)::Float64 = awastar.nodes_expended / simulated_expansion_rate

@inline function wall_time(awastar::AWAStar)::Float64
    if isnothing(awastar.stop_time)
        return time() - awastar.start_time
    else
        return awastar.stop_time - awastar.start_time
    end
end

@inline expansion_rate(awastar::AWAStar)::Float64 = awastar.nodes_expended / wall_time(awastar)

@inline quality(cost, optimal_cost)::Float64 = optimal_cost / cost

@inline function interrupt!(awastar::AWAStar)::Nothing
    awastar.interrupted = true
    nothing
end

function GraphSearch.init!(awastar::AWAStar{S,A}, sp::AbstractSearchProblem{S,A}) where {S,A}
    empty!(awastar.open_lists)
    empty!(awastar.closed_set)
    empty!(awastar.path_costs)
    awastar.nodes_expended = 0
    awastar.interrupted = false
    awastar.nodes_expended_to_fist_solution = nothing
    awastar.num_solutions = 0
    empty!(awastar.all_solutions_costs)
    empty!(awastar.all_solutions_walltimes)
    empty!(awastar.all_solutions_nodes_expended)
    empty!(awastar.all_solutions)
    awastar.solution = nothing
    awastar.solution_cost = Inf
    awastar.solution_walltime = nothing
    awastar.solution_nodes_expended = nothing
    awastar.solution_w_trace = nothing

    awastar.start_time = time()
    awastar.stop_time = nothing

    starting_state = start_state(sp)
    start_node = TreeSearchNode{S,A}(starting_state)
    start_node_g, start_node_h = start_node.path_cost, heuristic(sp, starting_state)
    @debug "AWA* Run" awastar.weight_adjustment_policy starting_state start_node_h
    push!(awastar.open_lists, starting_state, start_node, start_node_g, start_node_h)
    awastar.path_costs[starting_state] = start_node_g
end

function GraphSearch.stop_condition(awastar::AWAStar{S,A}, sp::AbstractSearchProblem{S,A})::Bool where {S,A}
    timeup::Bool = wall_time(awastar) >= awastar.walltime_limit
    budgetup::Bool = awastar.nodes_expended >= awastar.nodes_budget
    converged::Bool = false
    open_list_empty::Bool = false
    if length(awastar.open_lists) > 0
        _, g, h = peek(awastar.open_lists, 1.0)
        cost_lb::Float64 = g + h
        converged = cost_lb >= awastar.solution_cost
    else
        open_list_empty = true
    end
    stop = open_list_empty || converged || timeup || budgetup || awastar.interrupted
    if stop
        @debug "Stopping" open_list_empty converged timeup budgetup awastar.interrupted
    end
    return stop
end

function GraphSearch.node_expansion_policy(awastar::AWAStar{S,A}, sp::AbstractSearchProblem{S,A})::TreeSearchNode{S,A} where {S,A}
    node_f::Float64 = Inf
    local node::TreeSearchNode{S,A}
    while node_f >= awastar.solution_cost
        @assert length(awastar.open_lists) > 0 "Should not be possible"
        awastar.weight = awastar.weight_adjustment_policy(num_solutions=awastar.num_solutions)
        node, node_g, node_h = pop!(awastar.open_lists, awastar.weight)
        node.popped_using_w = awastar.weight
        node_f = node_g + node_h
    end
    return node
end

function GraphSearch.before_node_expansion!(awastar::AWAStar{S,A}, node::TreeSearchNode{S,A}, sp::AbstractSearchProblem{S,A}) where {S,A}
    push!(awastar.closed_set, node.state)
    return nothing
end

function GraphSearch.on_node_generation!(awastar::AWAStar{S,A}, child_node::TreeSearchNode{S,A}, sp::AbstractSearchProblem{S,A}) where {S,A}
    child_node_g, child_node_h = child_node.path_cost, heuristic(sp, child_node.state)
    child_node_f = child_node_g + child_node_h
    if child_node_f < awastar.solution_cost
        if goal_test(sp, child_node.state)
            awastar.path_costs[child_node.state] = child_node.path_cost
            if awastar.solution_cost == Inf
                awastar.nodes_expended_to_fist_solution = awastar.nodes_expended
            end
            awastar.solution_cost = child_node_f
            awastar.solution_walltime = wall_time(awastar)
            awastar.solution_nodes_expended = awastar.nodes_expended
            awastar.solution, awastar.solution_w_trace = get_solution(child_node)
            awastar.num_solutions += 1
            push!(awastar.all_solutions_walltimes, awastar.solution_walltime)
            push!(awastar.all_solutions_costs, awastar.solution_cost)
            push!(awastar.all_solutions_nodes_expended, awastar.solution_nodes_expended)
            push!(awastar.all_solutions, awastar.solution)
        elseif child_node.state in awastar.closed_set || haskey(awastar.open_lists, child_node.state)
            old_node_cost = awastar.path_costs[child_node.state]
            if old_node_cost > child_node.path_cost
                awastar.path_costs[child_node.state] = child_node.path_cost
                if child_node.state in awastar.closed_set
                    delete!(awastar.closed_set, child_node.state)
                else
                    delete!(awastar.open_lists, child_node.state)
                end
                push!(awastar.open_lists, child_node.state, child_node, child_node_g, child_node_h)
            end
        else
            awastar.path_costs[child_node.state] = child_node_g
            push!(awastar.open_lists, child_node.state, child_node, child_node_g, child_node_h)
        end
    end
    return nothing
end

function GraphSearch.on_node_expansion_finish!(awastar::AWAStar{S,A}, node::TreeSearchNode{S,A}, sp::AbstractSearchProblem{S,A}) where {S,A}
    awastar.nodes_expended += 1
    return nothing
end

function GraphSearch.on_stop!(awastar::AWAStar{S,A}, sp::AbstractSearchProblem{S,A}) where {S,A}
    awastar.stop_time = time()
    return nothing
end

function awastar_search_with_dynamic_weights(search_prob::AbstractSearchProblem{S,A}, weights_adjustment_policy::AbstractWeightAdjustmentPolicy, walltime_limit::Real, nodes_budget::Integer) where {S,A}
    awastar = AWAStar{S,A}(1, walltime_limit, nodes_budget, weights_adjustment_policy)
    return graph_search!(awastar, search_prob)
end

function awastar_search(search_prob::AbstractSearchProblem{S,A}, weight::Real, walltime_limit::Real, nodes_budget::Integer) where {S,A}
    return awastar_search_with_dynamic_weights(search_prob, ConstantWeightPolicy(weight), walltime_limit, nodes_budget)
end

ext: jl


text:
module TestShow
using Test
using Dates
using CategoricalArrays

@testset "show() for CategoricalPool{String} and its values" begin
    pool = CategoricalPool(["c", "b", "a"])
    opool = CategoricalPool(["c", "b", "a"], true)

    nv1 = CategoricalValue(pool, 1)
    nv2 = CategoricalValue(pool, 2)
    nv3 = CategoricalValue(pool, 3)

    ov1 = CategoricalValue(opool, 1)
    ov2 = CategoricalValue(opool, 2)
    ov3 = CategoricalValue(opool, 3)

    if VERSION >= v"1.6.0"
        @test sprint(show, pool) == "$CategoricalPool{String, UInt32}([\"c\", \"b\", \"a\"])"
        @test sprint(show, opool) == "$CategoricalPool{String, UInt32}([\"c\", \"b\", \"a\"]) with ordered levels"

        @test sprint(show, nv1) == repr(nv1) == "$CategoricalValue{String, UInt32} \"c\""
        @test sprint(show, nv2) == repr(nv2) == "$CategoricalValue{String, UInt32} \"b\""
        @test sprint(show, nv3) == repr(nv3) == "$CategoricalValue{String, UInt32} \"a\""

        @test sprint(show, ov1) == repr(ov1) == "$CategoricalValue{String, UInt32} \"c\" (1/3)"
        @test sprint(show, ov2) == repr(ov2) == "$CategoricalValue{String, UInt32} \"b\" (2/3)"
        @test sprint(show, ov3) == repr(ov3) == "$CategoricalValue{String, UInt32} \"a\" (3/3)"
    else
        @test sprint(show, pool) == "$CategoricalPool{String,UInt32}([\"c\", \"b\", \"a\"])"
        @test sprint(show, opool) == "$CategoricalPool{String,UInt32}([\"c\", \"b\", \"a\"]) with ordered levels"

        @test sprint(show, nv1) == repr(nv1) == "$CategoricalValue{String,UInt32} \"c\""
        @test sprint(show, nv2) == repr(nv2) == "$CategoricalValue{String,UInt32} \"b\""
        @test sprint(show, nv3) == repr(nv3) == "$CategoricalValue{String,UInt32} \"a\""

        @test sprint(show, ov1) == repr(ov1) == "$CategoricalValue{String,UInt32} \"c\" (1/3)"
        @test sprint(show, ov2) == repr(ov2) == "$CategoricalValue{String,UInt32} \"b\" (2/3)"
        @test sprint(show, ov3) == repr(ov3) == "$CategoricalValue{String,UInt32} \"a\" (3/3)"
    end

    @test sprint(show, nv1, context=:typeinfo=>typeof(nv1)) == "\"c\""
    @test sprint(show, nv2, context=:typeinfo=>typeof(nv2)) == "\"b\""
    @test sprint(show, nv3, context=:typeinfo=>typeof(nv3)) == "\"a\""

    @test sprint(show, ov1, context=:typeinfo=>typeof(ov1)) == "\"c\""
    @test sprint(show, ov2, context=:typeinfo=>typeof(ov2)) == "\"b\""
    @test sprint(show, ov3, context=:typeinfo=>typeof(ov3)) == "\"a\""

    @test sprint(show, nv1, context=:compact=>true) == "\"c\""
    @test sprint(show, nv2, context=:compact=>true) == "\"b\""
    @test sprint(show, nv3, context=:compact=>true) == "\"a\""

    @test sprint(show, ov1, context=:compact=>true) == "\"c\""
    @test sprint(show, ov2, context=:compact=>true) == "\"b\""
    @test sprint(show, ov3, context=:compact=>true) == "\"a\""

    @test sprint(print, nv1) == sprint(print, ov1) == "c"
    @test sprint(print, nv2) == sprint(print, ov2) == "b"
    @test sprint(print, nv3) == sprint(print, ov3) == "a"

    @test string(nv1) == string(ov1) == "c"
    @test string(nv2) == string(ov2) == "b"
    @test string(nv3) == string(ov3) == "a"

    @test String(nv1) == String(ov1) == "c"
    @test String(nv2) == String(ov2) == "b"
    @test String(nv3) == String(ov3) == "a"

    b = IOBuffer()
    @test write(b, nv1) == 1
    @test String(take!(b)) == "c"
    @test write(b, nv2) == 1
    @test String(take!(b)) == "b"
    @test write(b, nv3) == 1
    @test String(take!(b)) == "a"
    @test write(b, ov1) == 1
    @test String(take!(b)) == "c"
    @test write(b, ov2) == 1
    @test String(take!(b)) == "b"
    @test write(b, ov3) == 1
    @test String(take!(b)) == "a"
end

using JSON
@testset "JSON.lower" for pool in (CategoricalPool(["a"]),
                                   CategoricalPool([1]),
                                   CategoricalPool([1.0]))
    v = CategoricalValue(pool, 1)
    @test JSON.lower(v) == JSON.lower(unwrap(v))
    @test typeof(JSON.lower(v)) == typeof(JSON.lower(unwrap(v)))
end

using JSON3
using StructTypes
@testset "JSON3.write" begin
    v = CategoricalValue(CategoricalPool(["a"]), 1)
    @test JSON3.write(v) === "\"a\""

    v = CategoricalValue(CategoricalPool([1]), 1)
    @test JSON3.write(v) === "1"
    @test StructTypes.numbertype(typeof(v)) === Int

    v = CategoricalValue(CategoricalPool([2.0]), 1)
    @test JSON3.write(v) === "2.0"
    @test StructTypes.numbertype(typeof(v)) === Float64

    v = CategoricalValue(CategoricalPool([BigFloat(3.0,10)]), 1)
    @test JSON3.write(v) === "3.0"
    @test StructTypes.numbertype(typeof(v)) === BigFloat

    v = CategoricalValue(CategoricalPool([true,false]), 2)
    @test JSON3.write(v) == "false"
    @test StructTypes.numbertype(typeof(v)) === Bool
end

end


text:
struct CuDeviceArray{T,N} <: AbstractArray{T,N}
    shape::NTuple{N,Int}
    ptr::Ptr{T}
end

@inline function Base.setindex!(A::CuDeviceArray{T}, x, index::Integer) where {T}
    @boundscheck checkbounds(A, index)
    Base.unsafe_store!(A.ptr, x, index)
end

@inline function Base.getindex(A::CuDeviceArray{T}, index::Integer) where {T}
    @boundscheck checkbounds(A, index)
    Base.unsafe_load(A.ptr, index)::T
end

Base.size(g::CuDeviceArray) = g.shape

Base.checkbounds(::CuDeviceArray, I...) = nothing

for dim in (:x, :y, :z)
    # Thread index
    fn = Symbol("threadIdx_$dim")
    @eval @inline $fn() = Int((ccall($"llvm.nvvm.read.ptx.sreg.tid.$dim", llvmcall, UInt32, ()))+UInt32(1))

    # Block size (#threads per block)
    fn = Symbol("blockDim_$dim")
    @eval @inline $fn() =  Int(ccall($"llvm.nvvm.read.ptx.sreg.ntid.$dim", llvmcall, UInt32, ()))

    # Block index
    fn = Symbol("blockIdx_$dim")
    @eval @inline $fn() = Int((ccall($"llvm.nvvm.read.ptx.sreg.ctaid.$dim", llvmcall, UInt32, ()))+UInt32(1))

    # Grid size (#blocks per grid)
    fn = Symbol("gridDim_$dim")
    @eval @inline $fn() =  Int(ccall($"llvm.nvvm.read.ptx.sreg.nctaid.$dim", llvmcall, UInt32, ()))
end

@inline gridDim() =   (x=gridDim_x(),   y=gridDim_y(),   z=gridDim_z())

@inline blockIdx() =  (x=blockIdx_x(),  y=blockIdx_y(),  z=blockIdx_z())

@inline blockDim() =  (x=blockDim_x(),  y=blockDim_y(),  z=blockDim_z())

@inline threadIdx() = (x=threadIdx_x(), y=threadIdx_y(), z=threadIdx_z())

function kernel_vadd(a, b, c)
    i = (blockIdx().x-1) * blockDim().x + threadIdx().x
    c[i] = a[i] + b[i]

    return nothing
end

code_llvm(stdout, kernel_vadd, Tuple{CuDeviceArray{Float32,1},
                                     CuDeviceArray{Float32,1},
                                     CuDeviceArray{Float32,1}})

using TinyCassette
struct Ctx end

code_llvm(stdout, TinyCassette.execute,
          Tuple{Ctx,
                typeof(kernel_vadd),
                CuDeviceArray{Float32,1},
                CuDeviceArray{Float32,1},
                CuDeviceArray{Float32,1}})


text:
module TestStats

using Base.Test
using DataArrays
using DataFrames
using Dates

using TimeData

println("Running statistical function tests")

a = rand(5, 5)
tm = TimeData.Timematr(a)
tn = TimeData.Timenum(a)

#####################
## columnwise mean ##
#####################

expMean = mean(a, 1)
arrMean1 = array(mean(tm, 1))
arrMean2 = array(mean(tn, 1))

@test_approx_eq(expMean, arrMean1)
@test_approx_eq(expMean, arrMean2)

##################
## rowwise mean ##
##################

include("/home/chris/.julia/v0.3/TimeData/src/TimeData.jl")
expMean = mean(a, 2)
arrMean1 = TimeData.core(TimeData.rowmeans(tm))
arrMean2 = TimeData.core(TimeData.rowmeans(tn))

@test_approx_eq(expMean, arrMean1)
@test_approx_eq(expMean, arrMean2)

@test names(TimeData.rowmeans(tm)) == [:mean_values]

expProd = prod()

end


text:
"""
	number_of_nodes(s)

Return the number of nodes in the Stream Graph.
The number of nodes in a Stream Graph is defined as:
```math
n=\\sum_{v \\in V} n_v = \\frac{\\left| W \\right|}{\\left| T \\right|}
```
Where
```math
n_v=\\frac{\\left| T_v \\right|}{\\left| T \\right|}
```
is the contribution of node v.

### Reference
- Matthieu Latapy, Tiphaine Viard and Clémence Magnien Social
  Networks Analysis and Mining, 8: 61, 2018. "Stream Graphs and
  Link Streams for the Modeling of Interactions over Time".
  [(arXiv)](https://arxiv.org/pdf/1710.04073.pdf)
"""
function number_of_nodes(s::Union{StreamGraph,DirectedStreamGraph})
	length(s.W) != 0 ? sum([contribution(s,n) for (k,n) in s.W]) : 0.0
end

"""
	number_of_nodes(ls)

Return the number of nodes in the Link Stream.
The number of nodes in a Link Stream reduces to:
```math
n=\\left| V \\right|
```

### Reference
- Matthieu Latapy, Tiphaine Viard and Clémence Magnien Social
  Networks Analysis and Mining, 8: 61, 2018. "Stream Graphs and
  Link Streams for the Modeling of Interactions over Time".
  [(arXiv)](https://arxiv.org/pdf/1710.04073.pdf)
"""
number_of_nodes(ls::Union{LinkStream,DirectedLinkStream})=length(ls.V)

"""
	number_of_nodes(s)

Return the number of nodes in a state.
"""
number_of_nodes(s::State)=length(s.nodes)

"""
	number_of_nodes(tc)

Return the number of nodes at the cursor's current state.
"""
number_of_nodes(tc::TimeCursor)=length(nodes(tc))

"""
	number_of_nodes(tc,t)

Return the number of nodes in Gt using a time cursor.
"""
number_of_nodes(tc::TimeCursor,t::Float64)=length(nodes(tc,t))

"""
	number_of_nodes(tc,t0,t1)

Return the number of nodes which presence intersects with [t0;t1] using a time cursor.
"""
number_of_nodes(tc::TimeCursor,t0::Float64,t1::Float64)=length(nodes(tc,t0,t1))

"""
	number_of_links(s)

Return the number of links in the Stream.
The number of links in a Stream is defined as:
```math
m=\\sum_{uv \\in V \\otimes V} m_{uv} = \\frac{\\left| E \\right|}{\\left| T \\right|}
```
Where
```math
m_{uv}=\\frac{\\left| T_{uv} \\right|}{\\left| T \\right|}
```
is the contribution of the pair of nodes uv.

### Reference
- Matthieu Latapy, Tiphaine Viard and Clémence Magnien Social
  Networks Analysis and Mining, 8: 61, 2018. "Stream Graphs and
  Link Streams for the Modeling of Interactions over Time".
  [(arXiv)](https://arxiv.org/pdf/1710.04073.pdf)
"""
function number_of_links(s::AbstractStream)
	length(s.E) != 0 ? sum([contribution(s,l) for (k,v) in s.E for (kk,l) in v]) : 0.0
end

"""
	number_of_links(s)

Return the number of links in a state.
"""
number_of_links(s::State)=length(s.links)

"""
	number_of_links(tc)

Return the number of links in a cursor's current state.
"""
number_of_links(tc::TimeCursor)=number_of_links(tc.S)

"""
	number_of_links(tc,t)

Return the number of links present at a given time using a time cursor.
"""
number_of_links(tc::TimeCursor,t::Float64)=length(links(tc,t))

"""
	number_of_links(tc,t0,t1)

Return the number of links which presence intersects with [t0;t1] using a time cursor.
"""
number_of_links(tc::TimeCursor,t0::Float64,t1::Float64)=length(links(tc,t0,t1))

"""
	length(p)

Return the number of jumps in a path.
"""
length(p::AbstractPath)=length(p.jumps)


text:
# Boldify and unboldify words to create and read a dimension out of an ascii code.

bf_gap_uc = Int('𝐀') - Int('A')
bf_gap_lc = Int('𝐚') - Int('a')

boldify(s::AbstractString) =
    map(c -> Char(Int(c) + ('a' ≤ c ≤ 'z') * bf_gap_lc +
        ('A' ≤ c ≤ 'Z') * bf_gap_uc), s)

deboldify(s::AbstractString) =
    map(c -> Char(Int(c) - ('𝐚' ≤ c ≤ '𝐳') * bf_gap_lc -
    ('𝐀' ≤ c ≤ '𝐙') * bf_gap_uc), s)


text:
"""
    Semicircle(r)

The Wigner semicircle distribution with radius parameter `r` has probability
density function

```math
f(x; r) = \\frac{2}{\\pi r^2} \\sqrt{r^2 - x^2}, \\quad x \\in [-r, r].
```

```julia
Semicircle(r)   # Wigner semicircle distribution with radius r

params(d)       # Get the radius parameter, i.e. (r,)
```

External links

* [Wigner semicircle distribution on Wikipedia](https://en.wikipedia.org/wiki/Wigner_semicircle_distribution)
"""
struct Semicircle{T<:Real} <: ContinuousUnivariateDistribution
    r::T
    Semicircle{T}(r::T) where {T <: Real} = new{T}(r)
end

function Semicircle(r::T) where {T <: Real}
    @check_args(Semicircle, r > 0)
    return Semicircle{T}(r)
end

Semicircle(r::T, ::NoArgCheck) where {T <: Real} = Semicircle{T}(r)
Semicircle(r::Integer) = Semicircle(float(r))

@distr_support Semicircle -d.r +d.r

params(d::Semicircle) = (d.r,)

mean(d::Semicircle) = zero(d.r)
var(d::Semicircle) = d.r^2 / 4
skewness(d::Semicircle) = zero(d.r)
median(d::Semicircle) = zero(d.r)
mode(d::Semicircle) = zero(d.r)
entropy(d::Semicircle) = log(π * d.r) - oftype(d.r, 0.5)

function pdf(d::Semicircle, x::Real)
    xx, r = promote(x, float(d.r))
    if insupport(d, xx)
        return 2 / (π * r^2) * sqrt(r^2 - xx^2)
    else
        return oftype(r, 0)
    end
end

function logpdf(d::Semicircle, x::Real)
    xx, r = promote(x, float(d.r))
    if insupport(d, xx)
        return log(oftype(r, 2) / π) - 2 * log(r) + log(r^2 - xx^2) / 2
    else
        return oftype(r, -Inf)
    end
end

function cdf(d::Semicircle, x::Real)
    xx, r = promote(x, float(d.r))
    if insupport(d, xx)
        u = xx / r
        return (u * sqrt(1 - u^2) + asin(u)) / π + one(xx) / 2
    elseif x < minimum(d)
        return zero(r)
    else
        return one(r)
    end
end

@quantile_newton Semicircle


text:
#Path methods

Path(v::AbstractVector{Tuple{Int, Int}}) = Path([Point(p...) for p in v])
Path(v::AbstractVector{CartesianIndex{2}}) = Path([Point(p) for p in v])

function draw!(img::AbstractArray{T, 2}, path::Path, color::T) where T<:Colorant
    vertices = [CartesianIndex(p.y, p.x) for p in path.vertices]
    for i in 1:length(vertices)-1
		draw!(img, LineSegment(vertices[i], vertices[i+1]), color)
    end
	img
end

#Polygon methods

Polygon(v::AbstractVector{Tuple{Int, Int}}) = Polygon([Point(p...) for p in v])
Polygon(v::AbstractVector{CartesianIndex{2}}) = Polygon([Point(p) for p in v])

function draw!(img::AbstractArray{T, 2}, polygon::Polygon, color::T) where T<:Colorant
    draw!(img, Path(polygon.vertices), color)
    draw!(img, LineSegment(first(polygon.vertices), last(polygon.vertices)), color)
end

#RegularPolygon methods

RegularPolygon(point::CartesianIndex{2}, side_count::Int, side_length::T, θ::U) where {T<:Real, U<:Real} =
    RegularPolygon(Point(point), side_count, side_length, θ)

function draw!(img::AbstractArray{T, 2}, rp::RegularPolygon, color::T) where T<:Colorant
    n = rp.side_count
    ρ = rp.side_length/(2*sin(π/n))
    polygon = Polygon([ Point(round(Int, rp.center.x + ρ*cos(rp.θ + 2π*k/n)), round(Int, rp.center.y + ρ*sin(rp.θ + 2π*k/n))) for k in 1:n ])
    draw!(img, polygon, color)
end


text:
using Test
using Random
using ArgTools

import Tar
import Pkg.GitTools

const NON_STDLIB_TESTS = Main == @__MODULE__

if NON_STDLIB_TESTS
    using SimpleBufferStream

    using Tar_jll
    if isdefined(Tar_jll, :tar)
        const gtar = Tar_jll.tar
    end
end

tree_hash(path::AbstractString) = bytes2hex(GitTools.tree_hash(path))

function gen_file(file::String, size::Int)
    open(file, write=true) do io
        for i = 1:size; write(io, i % UInt8); end
    end
end

function make_test_tarball(tar_create::Function = Tar.create)
    root = mktempdir()
    let i = 0, paths = String[]
        a_lengths = [0, 10, 154, 155, 156, 255]
        @static if Sys.iswindows()
            filter!(a -> 2 + a + length(root) < 260, a_lengths)
        end
        for a in a_lengths
            dir = joinpath(root, "d"^a)
            push!(paths, normpath("..", dir))
            a > 0 && mkdir(dir)
            b_lengths = [10, 99, 100, 101, 255]
            # don't attempt to create file paths > 260 on windows
            @static if Sys.iswindows()
                filter!(b -> 4 + b + a + length(root) < 260, b_lengths)
            end
            for b in b_lengths
                for s in [0, 511, 512, 513, 1000]
                    f = rpad("$s-", b, "f")
                    x = rpad("$s-", b, "x")
                    file = joinpath(dir, f)
                    exec = joinpath(dir, x)
                    gen_file(file, s)
                    gen_file(exec, s)
                    chmod(file, 0o644)
                    chmod(exec, 0o755)
                    push!(paths, file)
                    push!(paths, exec)
                end
                dir′ = joinpath(dir, "s"^b)
                mkpath(dir′)
                push!(paths, dir′)
                link = joinpath(dir, "l"^b)
                target = relpath(paths[i += 1], link)
                symlink(target, link)
                push!(paths, link)
                broken = joinpath(dir, "b"^b)
                if target != "."
                    symlink(chop(target), broken)
                    push!(paths, broken)
                end
            end
        end
    end
    hash = tree_hash(root)
    tarball = tar_create(root)
    rm(root, force=true, recursive=true)
    return tarball, hash
end

const empty_tree_sha1 =
    "4b825dc642cb6eb9a060e54bf8d69288fbee4904"
const empty_tree_sha256 =
    "6ef19b41225c5369f1c104d45d8d85efa9b057b53b14b4b9b939dd74decc5321"

function test_empty_hashes(tarball::AbstractString)
    @test empty_tree_sha1 == Tar.tree_hash(tarball, skip_empty=true)
    @test empty_tree_sha1 == Tar.tree_hash(tarball, skip_empty=false)
    @test empty_tree_sha1 == Tar.tree_hash(tarball, algorithm="git-sha1")
    @test empty_tree_sha256 == Tar.tree_hash(tarball, algorithm="git-sha256")
    open(tarball) do io
        @test empty_tree_sha1 == Tar.tree_hash(io, skip_empty=true)
    end
    open(tarball) do io
        @test empty_tree_sha1 == Tar.tree_hash(io, skip_empty=false)
    end
    open(tarball) do io
        @test empty_tree_sha1 == Tar.tree_hash(io, algorithm="git-sha1")
    end
    open(tarball) do io
        @test empty_tree_sha256 == Tar.tree_hash(io, algorithm="git-sha256")
    end
end

function check_tree_hash(hash::AbstractString, root::AbstractString)
    @test tree_hash(root) == hash
    rm(root, recursive=true)
end

function make_test_dir(gen_skip::Bool=false)
    dir = mktempdir()
    touch(joinpath(dir, "file"))
    if gen_skip
        touch(joinpath(dir, "file.skip"))
    end
    mkdir(joinpath(dir, "dir"))

    # Create a file that actually has content
    fpath = joinpath(dir, "dir", "file")
    touch(fpath)
    open(fpath, write=true) do io
        write(io, rand(UInt8, 1000))
    end
    mkdir(joinpath(dir, "empty"))
    if gen_skip
        touch(joinpath(dir, "dir", "file.skip"))
    end
    if gen_skip
        mkdir(joinpath(dir, "dir.skip"))
        touch(joinpath(dir, "dir.skip", "file"))
    end
    if !Sys.iswindows()
        symlink("file", joinpath(dir, "link"))
        if gen_skip
            symlink("file", joinpath(dir, "link.skip"))
        end
    end
    return dir
end

const test_dir_paths = ["dir/file", "empty", "file", "link"]
Sys.iswindows() && pop!(test_dir_paths)

# uses Tar.list(callback, tarball) API
function tar_count(tarball; kwargs...)
    n = 0
    Tar.list(tarball; kwargs...) do hdr
        n += 1
    end
    return n
end

function tar_write_file(io::IO, path::String, data::String, mode::Integer=0o644)
    n = ncodeunits(data)
    Tar.write_header(io, Tar.Header(path, :file, mode, n, ""))
    Tar.write_data(io, IOBuffer(data), size=n)
end

function tar_write_link(io::IO, path::String, link::String, mode::Integer=0o755)
    Tar.write_header(io, Tar.Header(path, :symlink, mode, 0, link))
end

function tar_write_dir(io::IO, path::String, mode::Integer=0o755)
    Tar.write_header(io, Tar.Header(path, :directory, mode, 0, ""))
end

"""
    ChaosBufferStream(input::IO; chunksizes, sleepamnts)

Acts as a highly-inconvenient BufferStream, allowing tests to purposefully break
up their data flow reads/writes into chunks of random size (taken via `rand(chunksizes)`)
and with a pause of `rand(sleepamnts)` in between each chunk.

Usage example:

    t_hash = Tar.tree_hash(ChaosBufferStream(io))

Or slightly more complicated:

    http_io = BufferStream()
    @async HTTP.get(url; response_stream=http_io)
    t_hash = Tar.tree_hash(ChaosBufferStream(http_io; chunksizes=4000:8000, sleepamnts=[1e-3, 2e-3]))
"""
function ChaosBufferStream(input::IO; chunksizes = 1024:2048, sleepamnts = 1e-3:1e-4:5e-3)
    output = BufferStream()
    @async begin
        while !eof(input)
            sleep(rand(sleepamnts))
            chunk = read(input, rand(chunksizes))
            write(output, chunk)
        end
        close(output)
    end
    return output
end


text:
# # Entropy Maximization
#
# Here is a constrained entropy maximization problem:
#
# $$
# \begin{array}{ll}
#     \text{maximize}   & -\sum_{i=1}^n x_i \log x_i \\
#     \text{subject to} & \mathbf{1}' x = 1 \\
#                   & Ax \leq b
# \end{array}
# $$
#
# where $x \in \mathbf{R}^n$ is our optimization variable and $A \in \mathbf{R}^{m \times n}, b \in \mathbf{R}^{m}$.
#
# To solve this, we can simply use the `entropy` operation Convex.jl provides.

using Convex, SCS

n = 25;
m = 15;
A = randn(m, n);
b = rand(m, 1);

x = Variable(n);
problem = maximize(entropy(x), sum(x) == 1, A * x <= b)
solve!(problem, SCS.Optimizer(verbose=false))
problem.optval

#-

x.value


text:
"""
_compute_mutliplier: Additional function to estimate Lagrange multiplier of the problems
    (guarantee if LICQ holds)
"""
function _compute_mutliplier(pb    :: AbstractNLPModel,
                             x     :: Iterate,
                             gx    :: Iterate,
                             cx    :: Iterate,
                             Jx    :: MatrixType;
                             active_prec_c :: Float64 = 1e-6,
                			 active_prec_b :: Float64 = 1e-6)

 n  = length(x)
 nc = cx == nothing ? 0 : length(cx)

 #active res_bounds
 Ib = findall(x->(norm(x) <= active_prec_b),
			      min(abs.(x - pb.meta.lvar),
				      abs.(x - pb.meta.uvar)))
 if nc != 0
  #active constraints
  Ic = findall(x->(norm(x) <= active_prec_c),
                   min(abs.(cx-pb.meta.ucon),
                   abs.(cx-pb.meta.lcon)))

  Jc = hcat(Matrix(1.0I, n, n)[:,Ib], Jx'[:,Ic])
 else
  Ic = []
  Jc = hcat(Matrix(1.0I, n, n)[:,Ib])
 end

 mu, lambda = zeros(n), zeros(nc)
 if (Ib != []) || (Ic  != [])
  l = Jc \ (- gx)
  mu[Ib], lambda[Ic] = l[1:length(Ib)], l[length(Ib)+1:length(l)]
 end

 return mu, lambda
end


text:
@testset "Lasso small (h. split, $T)" for T in [Float32, Float64, ComplexF32, ComplexF64]

    using ProximalOperators
    using ProximalAlgorithms
    using LinearAlgebra
    using AbstractOperators: MatrixOp
    using RecursiveArrayTools: ArrayPartition

    A1 = T[  1.0  -2.0   3.0;
             2.0  -1.0   0.0;
            -1.0   0.0   4.0;
            -1.0  -1.0  -1.0]
    A2 = T[ -4.0  5.0;
            -1.0  3.0;
            -3.0  2.0;
             1.0  3.0]
    A = hcat(MatrixOp(A1), MatrixOp(A2))
    b = T[1.0, 2.0, 3.0, 4.0]

    m, n1 = size(A1)
    _, n2 = size(A2)
    n = n1 + n2

    R = real(T)

    lam = R(0.1)*norm(A'*b, Inf)
    @test typeof(lam) == R

    f = Translate(SqrNormL2(R(1)), -b)
    g = SeparableSum(NormL1(lam), NormL1(lam))

    x_star = ArrayPartition(
        T[-3.877278911564627e-01, 0, 0],
        T[2.174149659863943e-02, 6.168435374149660e-01]
    )

    TOL = R(1e-4)

    @testset "ForwardBackward" begin

        ## Nonfast/Nonadaptive

        x0 = ArrayPartition(zeros(T, n1), zeros(T, n2))
        solver = ProximalAlgorithms.ForwardBackward{R}(tol=TOL)
        x, it = solver(x0, f=f, A=A, g=g, L=opnorm([A1 A2])^2)
        @test eltype(x) == T
        @test norm(x - x_star, Inf) <= TOL
        @test it < 150

        # Nonfast/Adaptive

        x0 = ArrayPartition(zeros(T, n1), zeros(T, n2))
        solver = ProximalAlgorithms.ForwardBackward{R}(tol=TOL, adaptive=true)
        x, it = solver(x0, f=f, A=A, g=g)
        @test eltype(x) == T
        @test norm(x - x_star, Inf) <= TOL
        @test it < 300

        # Fast/Nonadaptive

        x0 = ArrayPartition(zeros(T, n1), zeros(T, n2))
        solver = ProximalAlgorithms.ForwardBackward{R}(tol=TOL, fast=true)
        x, it = solver(x0, f=f, A=A, g=g, L=opnorm([A1 A2])^2)
        @test eltype(x) == T
        @test norm(x - x_star, Inf) <= TOL
        @test it < 100

        # Fast/Adaptive

        x0 = ArrayPartition(zeros(T, n1), zeros(T, n2))
        solver = ProximalAlgorithms.ForwardBackward{R}(tol=TOL, adaptive=true, fast=true)
        x, it = solver(x0, f=f, A=A, g=g)
        @test eltype(x) == T
        @test norm(x - x_star, Inf) <= TOL
        @test it < 200
    end

    @testset "ZeroFPR" begin

        # ZeroFPR/Nonadaptive

        x0 = ArrayPartition(zeros(T, n1), zeros(T, n2))
        solver = ProximalAlgorithms.ZeroFPR{R}(tol=TOL)
        x, it = solver(x0, f=f, A=A, g=g, L=opnorm([A1 A2])^2)
        @test eltype(x) == T
        @test norm(x - x_star, Inf) <= TOL
        @test it < 20

        # ZeroFPR/Adaptive

        x0 = ArrayPartition(zeros(T, n1), zeros(T, n2))
        solver = ProximalAlgorithms.ZeroFPR{R}(adaptive=true, tol=TOL)
        x, it = solver(x0, f=f, A=A, g=g)
        @test eltype(x) == T
        @test norm(x - x_star, Inf) <= TOL
        @test it < 20

    end

    @testset "PANOC" begin

        # PANOC/Nonadaptive

        x0 = ArrayPartition(zeros(T, n1), zeros(T, n2))
        solver = ProximalAlgorithms.PANOC{R}(tol=TOL)
        x, it = solver(x0, f=f, A=A, g=g, L=opnorm([A1 A2])^2)
        @test eltype(x) == T
        @test norm(x - x_star, Inf) <= TOL
        @test it < 20

        ## PANOC/Adaptive

        x0 = ArrayPartition(zeros(T, n1), zeros(T, n2))
        solver = ProximalAlgorithms.PANOC{R}(adaptive=true, tol=TOL)
        x, it = solver(x0, f=f, A=A, g=g)
        @test eltype(x) == T
        @test norm(x - x_star, Inf) <= TOL
        @test it < 20

    end

end


text:
using MittagLefflerFunctions
using Polynomials
using PyPlot

α = 3/4
β = 1.0

E = MLQuad(α, β, 10, :hyperbola)

Nx = 300
Ny = 300

xmin, xmax = -6.0, 2.0
ymin, ymax = -4.0, 4.0
x = range(xmin, xmax, length=Nx)
y = range(ymin, ymax, length=Ny)

z = Complex{Float64}[ complex(x[j], y[i]) for i=1:Ny, j=1:Nx ]
modE = abs.(E.(z))

figure(1)
contourf(x, y, log10.(modE), levels=range(-2, 2, length=13))
colorbar(ticks=range(-2, 2, length=5))
xlabel(L"$\Re z$")
ylabel(L"$\Im z$")
grid(true)
savefig("figure02.pdf")


text:
hostname = readlines(`hostname`)[1]
cmd = readlines(pipeline(`cat /proc/$(getpid())/cmdline`,`xargs -0 echo`))[1]

if hostname=="ux305fa"
    ### Sync the code to hpc for example with rsync; also recreate empty folders; removes unneded files as well
    run(`rsync -a --delete-excluded '--include=*.jl' '--include=*/' '--exclude=*' . 'hpc05:~/data'`)
    ## ssh -c to run the same command we used
    run(`ssh hpc05 -t "export PATH=/home/jerdmanis/julia/bin:\$PATH; cd ~/data; $cmd"`)
    ### sync back created jld files
    run(`rsync -a '--include=*.jld' '--include=*/' '--exclude=*' 'hpc05:~/data/' .`)
    info("Sync finished")
    exit()
elseif hostname=="hpc05"
    ### remove all jld files
    using ClusterManagers

    run(`mkdir -p data`)
    run(`date`)
    println("How many cores you will need?")
    np = parse(Int,readline())
    println("How much time will you need hh:mm:ss?")
    walltime = readline()
    if walltime==""
        walltime="01:00:00"
    end
    println("How much memory will you need GB?")
    memory = readline()
    if memory==""
        memory="2"
    end

    #addprocs_pbs(np,queue="q1",res_list="walltime=$walltime;mem=$(memory)GB")
    addprocs_pbs(np,queue="q1",res_list="walltime=$walltime,mem=$(memory)GB")

    ### here I also could redefine println and pmap for a better laod balancing
    @everywhere import Base.println
    @everywhere println(x::AbstractString) = run(`echo "$x"`)

    @everywhere BLAS.set_num_threads(1)

    info("Cluster Initialized")
    run(`date`)
end

max_stars_repo_name: akels/WeylDisc

max_issues_repo_name: akels/WeylDisc

max_forks_repo_name: akels/WeylDisc


text:
using LinearAlgebra: norm, I, UniformScaling

"""
    scale(a::T, P::LazySet) where {T<:Real}

 Scale the set `P` by the scalar `a`. The scaling corresponds to applying a `linear_map` to `a` multiplied with an appropriately sized unit matrix and the set `P`.

 ### Input
 - `a` -- scalar
 - `P` -- lazy set

 ### Output
 A lazy set which correspoinds to the set `P` scaled by `a`.

 ### Note
 Enables to apply concrete multiplication operator `*ᶜ` to a scalar and a set (in addition to matrix and set).
"""
LazySets.scale(α::T, P::LazySet) where {T<:Real} =
    linear_map(α*Matrix{T}(I, dim(P), dim(P)), P)

"""
    linear_map(U::UniformScaling, P::LazySet)

 Concrete linear map of the set `P` and the uniform scaling `U`.

 ### Input
 - `U` -- uniform scaling
 - `P` -- lazy set

 ### Output
 A lazy set which corresponds to the linear map of `P` with `U`.
"""
LazySets.linear_map(U::UniformScaling, P::LazySet) =
    scale(U.λ, P)

"""
    ρ_exact(d::AbstractVector, P::LazySet)

 Compute support value for the lazy set `P` in the direction `d`. Avoid numerical errors for directions `d` with small norm by normalizing `d`. See  [`ρ`](@ref) for more detail.

 ### Input
 - `d` -- vector
 - `P` -- lazy set

 ### Output
 Support value for lazy set `P` in the direction `d`.
"""
ρ_exact(d::AbstractVector, P::LazySet) =  ρ(normalize(d), P) * norm(d)

"""
    σ_exact(d::AbstractVector, P::LazySet)

 Compute support vector for the lazy set `P` in the direction `d`. Avoid numerical errors for directions `d` with small norm by normalizing `d`. See  [`ρ`](@ref) for more detail.

 ### Input
 - `d` -- vector
 - `P` -- lazy set

 ### Output
 Support vector for lazy set `P` in the direction `d`.
"""
σ_exact(d::AbstractVector, P::LazySet) =  σ(normalize(d), P)

"""
    ρ_matrix(M::AbstractArray, P::LazySet)

Compute support value for the lazy set `P` in the direction of eachrow of matrix `M`.

### Input
- `M` -- matrix with rows as normal direction of constraints
- `P` -- lazy set

### Output
Vector of support values for lazy set `P` in the direction of each row of `M`.
"""
function ρ_matrix(M::AbstractArray, P::LazySet)
    num_constr = size(M, 1)
    res = zeros(num_constr)
    for i=1:num_constr
        res[i] = ρ_exact(view(M, i, :), P)
    end
    return res
end

"""
    canonical_length(X::LazySet)

Returns the support function of the given set along the positive and negative canonical directions.

### Inputs

- `X` - convex set

### Outputs

A matrix with `n = dims(X)` rows and two columns. Each row stands for
one dimension of `X` whereas the first column is the minimum and the second
column is the maximum value of the corresponding dimension.
"""
function canonical_length(X::LazySet{N}) where {N<:Real}
    dims = LazySets.dim(X)
    x = Matrix{N}(undef, dims, 2)
    for j = 1:dims
        ej = SingleEntryVector(j, dims, one(N))
        x[j,:] = [-ρ_exact(-ej, X), ρ_exact(ej, X)]
    end
    return x
end

"""
    marginal_enlargment(P::LazySet, ε=LazySets._TOL_F64.rtol)

 Increase the size of a lazy set `P` with available H-representation marginally. Hopefully by a Ball with a radius of `ε`.

 ### Input
 - `P` -- lazy set
 - `ε` -- (optional, default=LazySets.,_TOL_F64.rtol) marginal incrase

 ### Ouput
 Lazy set with marginally larger size than `P`.

 ### Note
 This function is useful, if we want to compute the minimal RPI set according to e.g. [`minimal_RPI_set_lazy_rakovic`](@ref) for an error set where the origin lies on the boundary and not in the interior. After marginally enlarging the error set, the algorithms for computing the minimal RPI set can be applied since the origin is now contained in the interior.
"""
function marginal_enlargment(P::LazySet, ε=LazySets._TOL_F64.rtol)
    F, g = tosimplehrep(P)
    g .+= ε/norm(g)
    return HPolytope(F, g)
end

max_stars_repo_name: mforets/InvariantSets

max_issues_repo_name: mforets/InvariantSets


text:
import IntervalArithmetic: interval
import Distributions
import StatsBase

const POTENTIAL_UVAL_TYPES = Union{T1, T2} where {T1<:Number, T2} where  T2 <: AbstractUncertainValue

"""
    UncertainScalarPopulation(values, probs)
    UncertainScalarPopulation(values, probs::Vector{Number})
    UncertainScalarPopulation(values, probs::Statsbase.AbstractWeights)

An `UncertainScalarPopulation`, which consists of some population members (`values`)
and some weights (`probs`) that indicate the relative importance of the
population members (for example during resampling).

## Fields

- **`values`**: The members of the population. Can be either numerical values, any
    type of uncertain value defined in this package (including populations), and
    `Measurement` instances from Measurements.jl.
- **`probs`**: The probabilities of sampling each member of the population.

## Constructors

- If `values` contains only scalar numeric values, then the `values` field
    will be of type `Vector{Number}`.
- If `values` contains one or more uncertain values, then the `values` field
    will be of type `Vector{AbstractUncertainValue}`

## Example

```julia

# Uncertain population consisting of CertainValues (scalars get promoted to
# CertainValue), theoretical distributions and KDE distributions
pop1 = UncertainScalarPopulation(
    [3.0, UncertainValue(Normal, 0, 1), UncertainValue(Gamma, 2, 3),
    UncertainValue(Uniform, rand(1000))], [0.5, 0.5, 0.5, 0.5])

# Uncertain population consisting of scalar values
pop2 = UncertainScalarPopulation([1, 2, 3], rand(3))
pop3 = UncertainScalarPopulation([1, 2, 3], Weights(rand(3)))

# Uncertain population consisting of uncertain populations
pop4 = UncertainScalarPopulation([pop1, pop2], [0.1, 0.5])

# Uncertain population consisting of uncertain populations, a scalar and
# a normal distribution. Assign random weights.
vals = [pop1, pop2, 2, UncertainValue(Normal, 0.3, 0.014)]
pop5 = UncertainScalarPopulation(vals, Weights(rand(4)))
```
"""
struct UncertainScalarPopulation{T, PW <: StatsBase.AbstractWeights} <: AbstractScalarPopulation{T, PW}
    values::Vector{T}
    probs::PW
end

"""
    UncertainScalarPopulation(values::Vector, probabilities::Vector{Float64})

Construct a population from a vector of values and a vector of probabilities associated
to those values."""
function UncertainScalarPopulation(values::Vector{T1}, probabilities::Vector{T2}) where {T1 <: Number, T2 <: Number}
    if length(values) != length(probabilities)
        throw(ArgumentError("Lengths of values and probability vectors do not match."))
    end
    UncertainScalarPopulation(values, StatsBase.weights(probabilities))
end
function UncertainScalarPopulation(values::VT, probabilities) where VT <: Vector{ELTYPE} where {ELTYPE<:POTENTIAL_UVAL_TYPES}
    if length(values) != length(probabilities)
        throw(ArgumentError("Lengths of values and probability vectors do not match."))
    end
    UncertainScalarPopulation(UncertainValue.(values), StatsBase.weights(probabilities))
end

"""
    ConstrainedUncertainScalarPopulation(values, probs)
    ConstrainedUncertainScalarPopulation(values, probs::Vector{Number})
    ConstrainedUncertainScalarPopulation(values, probs::Statsbase.AbstractWeights)

A `ConstrainedUncertainScalarPopulation`, which consists of some population
members (`values`)and some weights (`probs`) that indicate the relative importance of
the population members (for example during resampling). The uncertain values
for this type is meant to consist of constrained uncertain values
(generated by calling `constrain(uval, sampling_constraint`) on them.

This is just a convenience type to indicate that the population has been
constrained. It behaves identically to `UncertainScalarPopulation`.

There are different constructors for different types of `values`:

- If `values` contains only scalar numeric values, then the `values` field
    will be of type `Vector{Number}`.
- If `values` contains one or more uncertain values, then the `values` field
    will be of type `Vector{AbstractUncertainValue}`

"""
struct ConstrainedUncertainScalarPopulation{T, PW <: StatsBase.AbstractWeights} <: AbstractScalarPopulation{T, PW}
    values::Vector{T}
    probs::PW
end

"""
    ConstrainedUncertainScalarPopulation(values::Vector, probabilities::Vector{Float64})

Construct a constrained population from a vector of values and a vector of
probabilities associated to those values.
"""
function ConstrainedUncertainScalarPopulation(values::Vector{T1}, probabilities::Vector{T2}) where {T1 <: Number, T2 <: Number}
    if length(values) != length(probabilities)
        throw(ArgumentError("Lengths of values and probability vectors do not match."))
    end
    ConstrainedUncertainScalarPopulation(float.(values), StatsBase.weights(probabilities))
end
function ConstrainedUncertainScalarPopulation(values::VT, probabilities) where VT <: Vector{ELTYPE} where {ELTYPE<:POTENTIAL_UVAL_TYPES}
    if length(values) != length(probabilities)
        throw(ArgumentError("Lengths of values and probability vectors do not match."))
    end
    ConstrainedUncertainScalarPopulation(UncertainValue.(values), StatsBase.weights(probabilities))
end

export
UncertainScalarPopulation,
ConstrainedUncertainScalarPopulation


text:
# ---
# title: 912. Sort an Array
# id: problem912
# author: Indigo
# date: 2021-07-12
# difficulty: Medium
# categories:
# link: <https://leetcode.com/problems/sort-an-array/description/>
# hidden: true
# ---
#
# Given an array of integers `nums`, sort the array in ascending order.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [5,2,3,1]
#     Output: [1,2,3,5]
#
#
# **Example 2:**
#
#
#
#     Input: nums = [5,1,1,2,0,0]
#     Output: [0,0,1,1,2,5]
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 50000`
#   * `-50000 <= nums[i] <= 50000`
#
#
## @lc code=start
using LeetCode

sort_array(nums::Vector{Int}) = sort!(nums)
## @lc code=end


text:
# TODO bounds checking in constructors for all indices

abstract type AbstractPrimitive{J} <: ArrowVector{J} end
export AbstractPrimitive

"""
    Primitive{J} <: AbstractPrimitive{J}

A set of values which are stored contiguously in memory.  All other `ArrowVector` objects are built
from `Primitive`s.

## Constructors
    Primitive{J}(data::Vector{UInt8}, i::Integer, len::Integer)
    Primitive(data::Vector{UInt8}, i::Integer, x::AbstractVector)
    Primitive(v::AbstractVector)
    Primitive(Array, x::Abstractvector)

If `Array` is provided as the virst argument, the data will be allocated contiguously in a new buffer.

### Arguments
- `data`: a data buffer which the `Primitive` will refer to for accessing and storing values
- `i`: the index of the location (1-based) in `data` where the beginning of the value data is stored
- `x`: a vector which will be written into `data[i]` on construction
` `v`: existing reference data. constructors with will reference the original `v` as `data`
"""
struct Primitive{J} <: AbstractPrimitive{J}
    length::Int
    values_idx::Int64
    data::Vector{UInt8}
end
export Primitive

function Primitive{J}(data::AbstractVector{UInt8}, i::Integer, len::Integer) where J
    @boundscheck check_buffer_bounds(J, data, i, len)
    Primitive{J}(Int(len), Int(i), data)
end
function Primitive(data::AbstractVector{UInt8}, i::Integer, x::AbstractVector{J}) where J
    p = Primitive{J}(data, i, length(x))
    p[:] = x
    p
end
function Primitive{J}(data::AbstractVector{UInt8}, i::Integer, x::AbstractVector{K}) where {J,K}
    Primitive(data, i, convert(AbstractVector{J}, x))
end

# view of reinterpreted, will not include padding
function Primitive(v::AbstractVector{J}) where J
    # TODO in a future version this conversion will not be needed
    b = convert(Vector{UInt8}, reinterpret(UInt8, v))
    Primitive{J}(b, 1, length(v))
end
function Primitive{J}(v::AbstractVector{T}) where {J,T}
    Primitive(convert(AbstractVector{J}, v))
end

# create own buffer
function Primitive{J}(::Type{<:Array}, v::AbstractVector) where J
    b = Vector{UInt8}(undef, totalbytes(v))
    Primitive{J}(b, 1, v)
end
Primitive(::Type{<:Array}, v::AbstractVector) = Primitive(v)

# use undef buffer
function Primitive{J}(::UndefInitializer, i::Integer, len::Integer) where J
    b = Vector{UInt8}(undef, valuesbytes(J, len))
    Primitive{J}(b, i, len)
end

Primitive{J}(p::Primitive{J}) where J = Primitive{J}(p.length, p.values_idx, p.data)
Primitive{J}(p::Primitive{T}) where {J,T} = Primitive{J}(convert(AbstractVector{J}, p[:]))
Primitive(p::Primitive{J}) where J = Primitive{J}(p)

"""
    datapointer(A::Primitive)

Returns a pointer to the very start of the data buffer for `A` (i.e. does not depend on indices).
"""
datapointer(A::Primitive) = pointer(A.data)

valuespointer(A::Primitive) = datapointer(A) + A.values_idx - 1

#==================================================================================================
    NullablePrimitive
==================================================================================================#
"""
    NullablePrimitive{J} <: AbstractPrimitive{Union{J,Missing}}

A set of values stored contiguously in data with a bit mask specifying which values should be considered
null.  The bit mask and data needn't necessarily coexist within the same array.

## Constructors
    NullablePrimitive(bmask::Primitive{UInt8}, vals::Primitive{J})
    NullablePrimitive(data::Vector{UInt8}, bitmask_idx::Integer, values_idx::Integer, x::AbstractVector)
    NullablePrimitive(data::Vector{UInt8}, i::Integer, x::AbstractVector)
    NullablePrimitive(v::AbstractVector)
    NullablePrimitive(Array, x::AbstractVector)

If `Array` is passed to a constructor, the bit mask and values for the `NullablePrimitive` will be
contiguously allocated within a single array (bit mask first, then values).

### Arguments
- `bmask`: a `Primitive{UInt8}` containing the null bitmask for the `NullablePrimitive`
- `vals`: a `Primitive` containing the underlying data values for the `NullablePrimitive`
- `data`: a buffer in which the data for the `NullablePrimitive` will be stored
- `bitmask_idx`: the location within `data` where the null bit mask will be stored
- `values_idx`: the location within `data` where the values will be stored
- `x`, `v`: values to be stored in data
"""
struct NullablePrimitive{J} <: AbstractPrimitive{Union{J,Missing}}
    length::Int
    bitmask::Primitive{UInt8}
    values::Primitive{J}
end
export NullablePrimitive

# Primitive constructors
function NullablePrimitive{J}(bmask::Primitive{UInt8}, vals::Primitive{J}) where J
    NullablePrimitive{J}(length(vals), bmask, vals)
end
function NullablePrimitive(bmask::Primitive{UInt8}, vals::Primitive{J}) where J
    NullablePrimitive{J}(bmask, vals)
end

# buffer with location constructors
function NullablePrimitive{J}(data::Vector{UInt8}, bitmask_idx::Integer, values_idx::Integer,
                              len::Integer) where J
    bmask = Primitive{UInt8}(data, bitmask_idx, bitmaskbytes(len))
    vals = Primitive{J}(data, values_idx, len)
    NullablePrimitive{J}(bmask, vals)
end
function NullablePrimitive{J}(data::Vector{UInt8}, bitmask_idx::Integer, values_idx::Integer,
                              x::AbstractVector) where J
    bmask = Primitive{UInt8}(data, bitmask_idx, bitmaskpadded(x))
    vals = Primitive{J}(data, values_idx, length(x))
    setnonmissing!(vals, x)
    NullablePrimitive{J}(bmask, vals)
end
function NullablePrimitive(data::Vector{UInt8}, bitmask_idx::Integer, values_idx::Integer,
                           x::AbstractVector{J}) where J
    NullablePrimitive{J}(data, bitmask_idx, values_idx, x)
end
function NullablePrimitive(data::Vector{UInt8}, bitmask_idx::Integer, values_idx::Integer,
                           x::AbstractVector{Union{J,Missing}}) where J
    NullablePrimitive{J}(data, bitmask_idx, values_idx, x)
end
function NullablePrimitive(data::Vector{UInt8}, i::Integer, v::AbstractVector{T}
                          ) where {J,T<:Union{J,Union{J,Missing}}}
    NullablePrimitive(data, i, i+bitmaskbytes(v), v)
end
function NullablePrimitive{J}(data::Vector{UInt8}, i::Integer, v::AbstractVector{T}) where {J,T}
    NullablePrimitive(data, i, convert(AbstractVector{J}, v))
end

# contiguous new buffer constructors
function NullablePrimitive(::Type{<:Array}, v::AbstractVector{Union{J,Missing}}) where J
    b = Vector{UInt8}(undef, totalbytes(v))
    NullablePrimitive(b, 1, v)
end
function NullablePrimitive{J}(::Type{K}, v::AbstractVector{T}) where {J,K<:Array,T}
    NullablePrimitive(K, convert(AbstractVector{Union{J,Missing}}, v))
end
function NullablePrimitive(::Type{K}, v::AbstractVector{J}) where {K<:Array,J}
    NullablePrimitive(K, convert(AbstractVector{Union{J,Missing}}, v))
end

# new buffer constructors
function NullablePrimitive(v::AbstractVector{Union{J,Missing}}) where J
    bmask = Primitive(bitmaskpadded(v))
    vals = Primitive(replace_missing_vals(v))  # using first ensures exists
    NullablePrimitive(bmask, vals)
end
function NullablePrimitive(v::AbstractVector{J}) where J
    NullablePrimitive(convert(AbstractVector{Union{J,Missing}}, v))
end
function NullablePrimitive{J}(v::AbstractVector{K}) where {J,K}
    NullablePrimitive(convert(AbstractVector{Union{J,Missing}}, v))
end

function NullablePrimitive{J}(p::NullablePrimitive{J}) where J
    NullablePrimitive{J}(p.length, p.bitmask, p.values)
end
function NullablePrimitive{J}(p::NullablePrimitive{T}) where {J,T}
    NullablePrimitive{J}(convert(AbstractVector{J}, p[:]))
end
NullablePrimitive(p::NullablePrimitive{J}) where J = NullablePrimitive{J}(p)

#================================================================================================

    common interface

================================================================================================#
"""
    valuesbytes(A::AbstractVector)
    valuesbytes(::Type{C}, A::AbstractVector{<:AbstractString})

Computes the number of bytes needed to store the *values* of `A` (without converting the underlying
binary type). This does not include the number of bytes needed to store metadata such as a null
bitmask or offsets.

To obtain the number of values bytes needed to string data, one must input `C` the character encoding
type the string will be converted to (e.g. `UInt8`).
"""
valuesbytes(::Type{J}, len::Integer) where J = padding(sizeof(J)*len)
valuesbytes(A::AbstractVector{J}) where J = valuesbytes(J, len)
valuesbytes(A::AbstractVector{Union{J,Missing}}) where J = valuesbytes(J, len)
export valuesbytes

"""
    bitmaskbytes(A::AbstractVector)
    bitmaskbytes(::Type{Union{J,Missing}}, A::AbstractVector)

Compute the number of bytes needed to store a null bitmask for the data in `A`.  This is 0
unless `J <: Union{K,Missing}`. Note that this does not take into account scheme-dependent padding.
"""
bitmaskbytes(len::Integer) = padding(bytesforbits(len))
bitmaskbytes(A::AbstractVector) = 0
bitmaskbytes(::Type{Union{J,Missing}}, A::AbstractVector{J}) where J = bitmaskbytes(length(A))
function minbitmaskbytes(::Type{Union{J,Missing}}, A::AbstractVector{Union{J,Missing}}) where J
    bitmaskbytes(length(A))
end
bitmaskbytes(A::AbstractVector{Union{J,Missing}}) where J = bitmaskbytes(length(A))
export bitmaskbytes

minbitmaskbytes(len::Integer) = bytesforbits(len)
minbitmaskbytes(A::AbstractVector) = minbitmaskbytes(length(A))

"""
    totalbytes(A::AbstractVector)
    totalbytes(::Type{Union{J,Missing}}, A::AbstractVector)
    totalbytes(::Type{C}, A::AbstractVector)
    totalbytes(::Type{Union{J,Missing}}, ::Type{C}, A::AbstractVector)

Computes the minimum number of bytes needed to store `A` as an Arrow formatted primitive array or list.

To obtain the minimum bytes to store string data, one must input `C` the character encoding type the
string will be converted to (e.g. `UInt8`).
"""
totalbytes(A::AbstractVector) = bitmaskbytes(A) + valuesbytes(A)
function totalbytes(::Type{Union{J,Missing}}, A::AbstractVector{J}) where J
    bitmaskbytes(Union{J,Missing}, A) + valuesbytes(A)
end
function totalbytes(::Type{Union{J,Missing}}, A::AbstractVector{Union{J,Missing}}) where J
    bitmaskbytes(Union{J,Missing}, A) + valuesbytes(A)
end
export totalbytes

function unsafe_getvalue_contiguous(A::Union{Primitive{J},NullablePrimitive{J}},
                                    idx::AbstractVector{<:Integer}) where J
    ptr = convert(Ptr{J}, valuespointer(A)) + (first(idx)-1)*sizeof(J)
    unsafe_wrap(Array, ptr, length(idx))
end

"""
    unsafe_getvalue(A::ArrowVector, i)

Retrieve the value from memory location `i` using Julia 1-based indexing. `i` can be a single integer
index, an `AbstractVector` of integer indices, or an `AbstractVector{Bool}` mask.

This typically involves a call to `unsafe_load` or `unsafe_wrap`.
"""
function unsafe_getvalue(A::Union{Primitive{J},NullablePrimitive{J}},
                         i::Integer)::J where J
    unsafe_load(convert(Ptr{J}, valuespointer(A)), i)
end
function unsafe_getvalue(A::AbstractPrimitive, idx::UnitRange{<:Integer})
    unsafe_getvalue_contiguous(A, idx)
end
function unsafe_getvalue(A::AbstractPrimitive{T}, idx::AbstractVector{<:Integer}) where T
    T[unsafe_getvalue(A, i) for i ∈ idx]
end
function unsafe_getvalue(A::AbstractPrimitive{T}, idx::AbstractVector{Bool}) where T
    T[unsafe_getvalue(A, i) for i ∈ 1:length(A) if idx[i]]
end
unsafe_getvalue(A::Primitive, ::Colon) = unsafe_getvalue(A, 1:length(A))

_rawvalueindex_start(A::Primitive{J}, i::Integer) where J = A.values_idx + (i-1)*sizeof(J)
_rawvalueindex_stop(A::Primitive{J}, i::Integer) where J = A.values_idx + i*sizeof(J) - 1

function rawvalueindex_contiguous(A::Primitive{J}, idx::AbstractVector{<:Integer}) where J
    a = _rawvalueindex_start(A, first(idx))
    b = _rawvalueindex_stop(A, last(idx))
    a:b
end

function rawvalueindex(A::Primitive{J}, i::Integer) where J
    a = _rawvalueindex_start(A, i)
    b = _rawvalueindex_stop(A, i)
    a:b
end
function rawvalueindex(A::Primitive, idx::AbstractVector{<:Integer})
    vcat((rawvalueindex(A, i) for i ∈ idx)...)  # TODO inefficient use of vcat
end
rawvalueindex(A::Primitive, idx::UnitRange{<:Integer}) = rawvalueindex_contiguous(A, idx)
function rawvalueindex(A::Primitive, idx::AbstractVector{Bool})
    rawvalueindex(A, [i for i ∈ 1:length(A) if idx[i]])
end

"""
    rawvalues(A::Primitive, i)
    rawvalues(A::Primitive)

Gets the raw values for elements at locations `i` in the form of a `Vector{UInt8}`.
"""
rawvalues(A::Primitive, i::Union{<:Integer,AbstractVector{<:Integer}}) = A.data[rawvalueindex(A, i)]
rawvalues(A::Primitive, ::Colon) = rawvalues(A, 1:length(A))
rawvalues(A::Primitive) = rawvalues(A, :)

"""
    getvalue(A::ArrowVector, idx)

Get the values for indices `idx` from `A`.
"""
function getvalue(A::Union{Primitive{J},NullablePrimitive{J}}, i::Integer) where J
    reinterpret(J, rawvalues(A, i))[1]
end
function getvalue(A::Primitive{J}, i::AbstractVector{<:Integer}) where J
    reinterpret(J, rawvalues(A, i))
end
function getvalue(A::NullablePrimitive{J}, i::AbstractVector{<:Integer}) where J
    convert(Union{J,Missing}, reinterpret(J, rawvalues(A, i)))
end

"""
    unsafe_rawpaddedvalues(p::ArrowVector)

Retreive raw value data for `p` as a `Vector{UInt8}`.
"""
function unsafe_rawpaddedvalues(p::AbstractPrimitive)
    unsafe_rawpadded(valuespointer(p), valuesbytes(p))
end

function setvalue!(A::Primitive{J}, x::J, i::Integer) where J
    A.data[rawvalueindex(A, i)] = reinterpret(UInt8, [x])
end
function setvalue!(A::Primitive{J}, x::Vector{J}, idx::AbstractVector{<:Integer}) where J
    A.data[rawvalueindex(A, idx)] = reinterpret(UInt8, x)
end
function setvalue!(A::NullablePrimitive{J}, x::J, i::Integer) where J
    setvalue!(A.values, x, i)
end
function setvalue!(A::NullablePrimitive{J}, x::Vector{J}, idx::AbstractVector{<:Integer}) where J
    setvalue!(A.values, x, idx)
end

"""
    unsafe_setvalue!(A::ArrowVector{J}, x, i)

Set the value at location `i` to `x`.  If `i` is a single integer, `x` should be an element of type
`J`.  Otherwise `i` can be an `AbstractVector{<:Integer}` or `AbstractVector{Bool}` in which case
`x` should be an appropriately sized `AbstractVector{J}`.
"""
function unsafe_setvalue!(A::Union{Primitive{J},NullablePrimitive{J}}, x::J, i::Integer) where J
    unsafe_store!(convert(Ptr{J}, valuespointer(A)), x, i)
end
function unsafe_setvalue!(A::Union{Primitive{J},NullablePrimitive{J}}, v::AbstractVector{J},
                          idx::AbstractVector{<:Integer}) where J
    ptr = convert(Ptr{J}, valuespointer(A))
    for (x, i) ∈ zip(v, idx)
        unsafe_store!(ptr, x, i)
    end
end
function unsafe_setvalue!(A::Union{Primitive{J},NullablePrimitive{J}}, v::AbstractVector{J},
                          idx::AbstractVector{Bool}) where J
    ptr = convert(Ptr{J}, valuespointer(A))
    j = 1
    for i ∈ 1:length(A)
        if idx[i]
            unsafe_store!(ptr, v[j], i)
            j += 1
        end
    end
end
function unsafe_setvalue!(A::Union{Primitive{J},NullablePrimitive{J}}, v::Vector{J}, ::Colon) where J
    unsafe_copy!(convert(Ptr{J}, valuespointer(A)), pointer(v), length(v))
end

"""
    unsafe_construct(::Type{T}, A::Primitive, i::Integer, len::Integer)

Construct an object of type `T` using `len` elements from `A` starting at index `i` (1-based indexing).
This is mostly used by `AbstractList` objects to construct variable length objects such as strings
from primitive arrays.

Users must define new methods for new types `T`.
"""
function unsafe_construct(::Type{String}, A::Primitive{UInt8}, i::Integer, len::Integer)
    unsafe_string(convert(Ptr{UInt8}, valuespointer(A) + (i-1)), len)
end
# this is a fallback for other types and is safe
function unsafe_construct(::Type{T}, A::AbstractPrimitive, i::Integer, len::Integer) where T
    construct(T, A, i, len)
end
function unsafe_construct(::Type{T}, A::NullablePrimitive, i::Integer, len::Integer) where T
    throw(ErrorException("Unsafely constructing objects from NullablePrimitves is not supported."))
end

"""
    construct(::Type{T}, A::AbstractPrimitive{J}, i::Integer, len::Integer)

Construct an object of type `T` from `len` values in `A` starting at index `i`.
For this to work requires the existence of a constructor of the form `T(::Vector{J})`.
"""
function construct(::Type{T}, A::AbstractPrimitive, i::Integer, len::Integer) where T
    T(A[i:(i+len-1)])  # obviously this depends on the existence of this constructor
end

function setindex!(A::Primitive{J}, x, i::Integer) where J
    @boundscheck checkbounds(A, i)
    setvalue!(A, convert(J, x), i)  # should this conversion really be here?
end
# TODO inefficient in some cases because of conversion to Vector{J}
function setindex!(A::Primitive{J}, x::AbstractVector, idx::AbstractVector{<:Integer}) where J
    @boundscheck (checkbounds(A, idx); checkinputsize(x, idx))
    setvalue!(A, convert(Vector{J}, x), idx)
end
setindex!(A::Primitive, x::AbstractVector, ::Colon) = (A[1:end] = x)

function setindex!(A::NullablePrimitive{J}, x, i::Integer) where J
    @boundscheck checkbounds(A, i)
    o = setvalue!(A, convert(J, x), i)
    setnull!(A, false, i)  # important that this is last in case above fails
    o
end
function setindex!(A::NullablePrimitive{J}, x::Missing, i::Integer) where J
    @boundscheck checkbounds(A, i)
    setnull!(A, true, i)
    missing
end
# TODO this is horribly inefficient but really hard to do right for non-consecutive
function setindex!(A::NullablePrimitive, x::AbstractVector, idx::AbstractVector{<:Integer})
    @boundscheck (checkbounds(A, idx); checkinputsize(x, idx))
    for (ξ,i) ∈ zip(x, idx)
        @inbounds setindex!(A, ξ, i)
    end
end
function setindex!(A::NullablePrimitive, x::AbstractVector, idx::AbstractVector{Bool})
    @boundscheck (checkbounds(A, idx); checkinputsize(x, idx))
    j = 1
    for i ∈ 1:length(A)
        if idx[i]
            @inbounds setindex!(A, x[j], i)
            j += 1
        end
    end
    x
end
# TODO this probably isn't really much more efficient, should test
function setindex!(A::NullablePrimitive{J}, x::AbstractVector, ::Colon) where J
    @boundscheck checkinputsize(x, A)
    setnulls!(A, ismissing.(x))
    for i ∈ 1:length(A)
        !ismissing(x[i]) && setvalue!(A, convert(J, x[i]), i)
    end
    x
end


text:
#!/usr/bin/env julia
const MOD = 123
const YES = "yes"
const NO = "NO"

function solve(N::Int, M::Int, H::Matrix{String}, A::Vector{Int}, B::Vector{Float64}, Q::Int, X::Vector{Int})
    println("$N $M")
    @assert size(H) == (N-1, M-2)
    for i in 1:size(H,1)
        println(join(H[i, :], " "))
    end
    @assert length(A) == N - 1
    @assert length(B) == N - 1
    for (a, b) in zip(A, B)
        println("$a $b")
    end
    println(Q)
    @assert length(X) == M + Q
    foreach(println, X)

    println(YES)
    println(NO)
    println(MOD)
end

function main()
    tokens = Channel{String}(32)
    Task() do
        for line in eachline(@static VERSION < v"0.6" ? STDIN : stdin)
            for token in split(chomp(line))
                put!(tokens, token)
            end
        end
        close(tokens)
    end |> schedule
    N = parse(Int, take!(tokens))
    M = parse(Int, take!(tokens))
    H = similar(Matrix{String}, N-2+1, M-1-2+1)
    for i in 1:N-2+1
        for j in 1:M-1-2+1
            H[i,j] = take!(tokens)
        end
    end
    A = similar(Vector{Int}, N-2+1)
    B = similar(Vector{Float64}, N-2+1)
    for i in 1:N-2+1
        A[i] = parse(Int, take!(tokens))
        B[i] = parse(Float64, take!(tokens))
    end
    Q = parse(Int, take!(tokens))
    X = similar(Vector{Int}, M+Q)
    for i in 1:M+Q
        X[i] = parse(Int, take!(tokens))
    end
    solve(N, M, H, A, B, Q, X)
end

isempty(ARGS) && main()


text:
@testset "313.super-ugly-number.jl" begin
    primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
              73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151,
              157, 163, 167, 173, 179, 181, 191, 193, 197, 199]
    @test nth_super_ugly_number(1, primes) == 1
    @test nth_super_ugly_number(22, primes) == 22
    @test nth_super_ugly_number(222222, primes) == 1922232
end


text:
@testset "align_tensors (unprofiled)" begin

    @inferred align_tensors(A_TOST, B_TOST;method=ΛTAME_M(),matchingMethod = ΛTAME_GramMatching())
    @inferred align_tensors(A_TOST, B_TOST;method=ΛTAME_M(),matchingMethod = ΛTAME_rankOneMatching())
    @inferred align_tensors(A_TOST, B_TOST;method=LowRankTAME_M())
    @inferred align_tensors(A_TOST, B_TOST;method=TAME_M())
    @inferred align_tensors(A_TOST, B_TOST;method=ΛTAME_MultiMotif_M())

end

@testset "align_tensors_profiled" begin

    #for () in [(A_TOST, B_TOST), (A_UTOST,B_UTOST),(A_UST, B_UST)]
    @inferred align_tensors_profiled(A_TOST, B_TOST;method=ΛTAME_M(),matchingMethod = ΛTAME_GramMatching())
    @inferred align_tensors_profiled(A_TOST, B_TOST;method=ΛTAME_M(),matchingMethod = ΛTAME_rankOneMatching())
    @inferred align_tensors_profiled(A_TOST, B_TOST;method=LowRankTAME_M())
    @inferred align_tensors_profiled(A_TOST, B_TOST;method=TAME_M())
    #@inferred align_tensors_profiled(A_TOST, B_TOST;method=ΛTAME_MultiMotif_M())
        #not-implemnted right now
    #end

end

max_stars_repo_name: chuckcol/LambdaTAME

max_issues_repo_name: chuckcol/LambdaTAME

max_forks_repo_name: chuckcol/LambdaTAME


text:
export import_mesh, export_mesh
"""
    import_mesh(path::String)
    import_mesh(path::String, ::Type{T}=Float64) where {T<:AbstractFloat}

Import a mesh from file. The float type of the mesh may be specified with a second argument.
File type is inferred from the extension.
"""
function import_mesh(path::String, ::Type{T}) where {T<:AbstractFloat}
    @info "Reading "*path
    if endswith(path, ".inp")
        return read_abaqus(path, T)
    elseif endswith(path, ".xdmf")
        return read_xdmf(path, T)
    else
        error("Could not determine mesh file type from extension.")
    end
end

import_mesh(path::String) = import_mesh(path, Float64)

function export_mesh(path::String, mesh)
    @info "Writing "*path
    if endswith(path, ".xdmf")
        return write_xdmf(path, mesh)
    else
        error("Could not determine mesh file type from extension")
    end
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
using MetricsProject
using Test
using Random

@testset "MetricsProject.jl" begin
    # Tests for Berry1994
    Random.seed!(1)
    out = Berry1994(50, 10).comparison
    @test out.m_θhatⁱᵛ[1]  ≈ 3.906248115012579
    @test out.m_θhatⁱᵛ[2]  ≈ 1.695701084241484
    @test out.m_θhatⁱᵛ[3]  ≈ -0.8016480548335396

    # Tests for sdc (letting the default seed parameter = 0)
    out = sdc(500, 5).comparison
    @test out.θhat_t1ev[5] ≈ 1.566535306145989
    @test out.θhat_normal[5] ≈ 1.331458650920062
    # The estimatives must be invariant to re-scaling the mean utilities
    out2 = sdc(500, 5, μ = 10).comparison
    @test out.θhat_normal[5] ≈ out2.θhat_normal[5]

end