text:
export USSettlementCalendar, USNYSECalendar, USGovernmentBondCalendar, USNERCCalendar

abstract UnitedStatesCalendar <: WesternCalendar

type USSettlementCalendar <: UnitedStatesCalendar; end
type USNYSECalendar <: UnitedStatesCalendar; end
type USNERCCalendar <: UnitedStatesCalendar; end
type USGovernmentBondCalendar <: UnitedStatesCalendar; end

string(::USSettlementCalendar) = "US Settlement"
string(::USNYSECalendar) = "New York Stock Exchange"
string(::USNERCCalendar) = "North American Energy Reliability Council"
string(::USGovernmentBondCalendar) = "US Government Bond Market"

function isWorkingDay(c::USSettlementCalendar, dt::CalendarTime)
	d=day(dt)
	dd=dayofyear(dt)
	m=month(dt)
	w=dayofweek(dt)
	y=year(dt)
	em=easterMonday(c, y)
	if (isWeekend(c, w)
            # New Year's Day (possibly moved to Monday if on Sunday)
            || ((d == 1 || (d == 2 && w == Monday)) && m == January)
            # (or to Friday if on Saturday)
            || (d == 31 && w == Friday && m == December)
            # Martin Luther King's birthday (third Monday in January)
            || ((d >= 15 && d <= 21) && w == Monday && m == January)
            # Washington's birthday (third Monday in February)
            || ((d >= 15 && d <= 21) && w == Monday && m == February)
            # Memorial Day (last Monday in May)
            || (d >= 25 && w == Monday && m == May)
            # Independence Day (Monday if Sunday or Friday if Saturday)
            || ((d == 4 || (d == 5 && w == Monday) ||
                 (d == 3 && w == Friday)) && m == July)
            # Labor Day (first Monday in September)
            || (d <= 7 && w == Monday && m == September)
            # Columbus Day (second Monday in October)
            || ((d >= 8 && d <= 14) && w == Monday && m == October)
            # Veteran's Day (Monday if Sunday or Friday if Saturday)
            || ((d == 11 || (d == 12 && w == Monday) || (d == 10 && w == Friday)) && m == November)
            # Thanksgiving Day (fourth Thursday in November)
            || ((d >= 22 && d <= 28) && w == Thursday && m == November)
            # Christmas (Monday if Sunday or Friday if Saturday)
            || ((d == 25 || (d == 26 && w == Monday) ||
                 (d == 24 && w == Friday)) && m == December))
            return false
    end
    return true
end

function isWorkingDay(c::USNYSECalendar, dt::CalendarTime)
	d=day(dt)
	dd=dayofyear(dt)
	m=month(dt)
	w=dayofweek(dt)
	y=year(dt)
	em=easterMonday(c, y)

	if (isWeekend(c, w)
            # New Year's Day (possibly moved to Monday if on Sunday)
            || ((d == 1 || (d == 2 && w == Monday)) && m == January)
            # Washington's birthday (third Monday in February)
            || ((d >= 15 && d <= 21) && w == Monday && m == February)
            # Good Friday
            || (dd == em-3)
            # Memorial Day (last Monday in May)
            || (d >= 25 && w == Monday && m == May)
            # Independence Day (Monday if Sunday or Friday if Saturday)
            || ((d == 4 || (d == 5 && w == Monday) ||
                 (d == 3 && w == Friday)) && m == July)
            # Labor Day (first Monday in September)
            || (d <= 7 && w == Monday && m == September)
            # Thanksgiving Day (fourth Thursday in November)
            || ((d >= 22 && d <= 28) && w == Thursday && m == November)
            # Christmas (Monday if Sunday or Friday if Saturday)
            || ((d == 25 || (d == 26 && w == Monday) || (d == 24 && w == Friday)) && m == December)
            )
		return false
	end

    if (y >= 1998)
        if (# Martin Luther King's birthday (third Monday in January)
            ((d >= 15 && d <= 21) && w == Monday && m == January)
            # President Reagan's funeral
            || (y == 2004 && m == June && d == 11)
            # September 11, 2001
            || (y == 2001 && m == September && (11 <= d && d <= 14))
            # President Ford's funeral
            || (y == 2007 && m == January && d == 2)
            )
        	return false
    	end
    elseif (y <= 1980)
        if (# Presidential election days
            ((y % 4 == 0) && m == November && d <= 7 && w == Tuesday)
            # 1977 Blackout
            || (y == 1977 && m == July && d == 14)
            # Funeral of former President Lyndon B. Johnson.
            || (y == 1973 && m == January && d == 25)
            # Funeral of former President Harry S. Truman
            || (y == 1972 && m == December && d == 28)
            # National Day of Participation for the lunar exploration.
            || (y == 1969 && m == July && d == 21)
            # Funeral of former President Eisenhower.
            || (y == 1969 && m == March && d == 31)
            # Closed all day - heavy snow.
            || (y == 1969 && m == February && d == 10)
            # Day after Independence Day.
            || (y == 1968 && m == July && d == 5)
            # June 12-Dec. 31, 1968
            # Four day week (closed on Wednesdays) - Paperwork Crisis
            || (y == 1968 && dd >= 163 && w == Wednesday)
            ) return false
		end
    else
    	if (# Nixon's funeral
            (y == 1994 && m == April && d == 27))
    		return false
    	end
    end

    return true
end

function isWorkingDay(c::USGovernmentBondCalendar, dt::CalendarTime)
	d=day(dt)
	dd=dayofyear(dt)
	m=month(dt)
	w=dayofweek(dt)
	y=year(dt)
	em=easterMonday(c,y)

	if (isWeekend(c, w)
            # New Year's Day (possibly moved to Monday if on Sunday)
            || ((d == 1 || (d == 2 && w == Monday)) && m == January)
            # Martin Luther King's birthday (third Monday in January)
            || ((d >= 15 && d <= 21) && w == Monday && m == January)
            # Washington's birthday (third Monday in February)
            || ((d >= 15 && d <= 21) && w == Monday && m == February)
            # Good Friday
            || (dd == em-3)
            # Memorial Day (last Monday in May)
            || (d >= 25 && w == Monday && m == May)
            # Independence Day (Monday if Sunday or Friday if Saturday)
            || ((d == 4 || (d == 5 && w == Monday) || (d == 3 && w == Friday)) && m == July)
            # Labor Day (first Monday in September)
            || (d <= 7 && w == Monday && m == September)
            # Columbus Day (second Monday in October)
            || ((d >= 8 && d <= 14) && w == Monday && m == October)
            # Veteran's Day (Monday if Sunday or Friday if Saturday)
            || ((d == 11 || (d == 12 && w == Monday) ||
                 (d == 10 && w == Friday)) && m == November)
            # Thanksgiving Day (fourth Thursday in November)
            || ((d >= 22 && d <= 28) && w == Thursday && m == November)
            # Christmas (Monday if Sunday or Friday if Saturday)
            || ((d == 25 || (d == 26 && w == Monday) ||
                 (d == 24 && w == Friday)) && m == December))
            return false
    end
    return true
end

function isWorkingDay(c::USNERCCalendar, dt::CalendarTime)
	d=day(dt)
	dd=dayofyear(dt)
	m=month(dt)
	w=dayofweek(dt)
	y=year(dt)
	em=easterMonday(c,y)

	if (isWeekend(c, w)
            # New Year's Day (possibly moved to Monday if on Sunday)
            || ((d == 1 || (d == 2 && w == Monday)) && m == January)
            # Memorial Day (last Monday in May)
            || (d >= 25 && w == Monday && m == May)
            # Independence Day (Monday if Sunday)
            || ((d == 4 || (d == 5 && w == Monday)) && m == July)
            # Labor Day (first Monday in September)
            || (d <= 7 && w == Monday && m == September)
            # Thanksgiving Day (fourth Thursday in November)
            || ((d >= 22 && d <= 28) && w == Thursday && m == November)
            # Christmas (Monday if Sunday)
            || ((d == 25 || (d == 26 && w == Monday)) && m == December))
            return false
     end
     return true
end


text:
AdaBoostRegressor_ = SKEN.AdaBoostRegressor
@sk_reg mutable struct AdaBoostRegressor <: MMI.Deterministic
    base_estimator::Any    = nothing
    n_estimators::Int      = 50::(_ > 0)
    learning_rate::Float64 = 1.0::(_ > 0)
    loss::String           = "linear"::(_ in ("linear","square","exponential"))
    random_state::Any      = nothing
end
MMI.fitted_params(model::AdaBoostRegressor, (f, _, _)) = (
    estimators           = f.estimators_,
    estimator_weights    = f.estimator_weights_,
    estimator_errors     = f.estimator_errors_,
    feature_importances_ = f.feature_importances_
    )

# ----------------------------------------------------------------------------
AdaBoostClassifier_ = SKEN.AdaBoostClassifier
@sk_clf mutable struct AdaBoostClassifier <: MMI.Probabilistic
    base_estimator::Any    = nothing
    n_estimators::Int      = 50::(_ > 0)
    learning_rate::Float64 = 1.0::(_ > 0)
    algorithm::String      = "SAMME.R"::(_ in ("SAMME", "SAMME.R"))
    random_state::Any      = nothing
end
MMI.fitted_params(m::AdaBoostClassifier, (f, _, _)) = (
    estimators        = f.estimators_,
    estimator_weights = f.estimator_weights_,
    estimator_errors  = f.estimator_errors_,
    classes           = f.classes_,
    n_classes         = f.n_classes_
    )
metadata_model(AdaBoostClassifier,
    input   = Table(Continuous),
    target  = AbstractVector{<:Finite},
    weights = false,
    descr   = "Adaboost ensemble classifier."
    )

# ============================================================================
BaggingRegressor_ = SKEN.BaggingRegressor
@sk_reg mutable struct BaggingRegressor <: MMI.Deterministic
    base_estimator::Any      = nothing
    n_estimators::Int        = 10::(_>0)
    max_samples::Union{Int,Float64}  = 1.0::(_>0)
    max_features::Union{Int,Float64} = 1.0::(_>0)
    bootstrap::Bool          = true
    bootstrap_features::Bool = false
    oob_score::Bool          = false
    warm_start::Bool         = false
    n_jobs::Option{Int}      = nothing
    random_state::Any        = nothing
    verbose::Int             = 0
end
MMI.fitted_params(model::BaggingRegressor, (f, _, _)) = (
    estimators          = f.estimators_,
    estimators_samples  = f.estimators_samples_,
    estimators_features = f.estimators_features_,
    oob_score           = model.oob_score ? f.oob_score_ : nothing,
    oob_prediction      = model.oob_score ? f.oob_prediction_ : nothing
    )

# ----------------------------------------------------------------------------
BaggingClassifier_ = SKEN.BaggingClassifier
@sk_clf mutable struct BaggingClassifier <: MMI.Probabilistic
    base_estimator::Any      = nothing
    n_estimators::Int        = 10::(_>0)
    max_samples::Union{Int,Float64}  = 1.0::(_>0)
    max_features::Union{Int,Float64} = 1.0::(_>0)
    bootstrap::Bool          = true
    bootstrap_features::Bool = false
    oob_score::Bool          = false
    warm_start::Bool         = false
    n_jobs::Option{Int}      = nothing
    random_state::Any        = nothing
    verbose::Int             = 0
end
MMI.fitted_params(m::BaggingClassifier, (f, _, _)) = (
    base_estimator        = f.base_estimator_,
    estimators            = f.estimators_,
    estimators_samples    = f.estimators_samples_,
    estimators_features   = f.estimators_features_,
    classes               = f.classes_,
    n_classes             = f.n_classes_,
    oob_score             = m.oob_score ? f.oob_score_ : nothing,
    oob_decision_function = m.oob_score ? f.oob_decision_function_ : nothing
    )
metadata_model(BaggingClassifier,
    input   = Table(Continuous),
    target  = AbstractVector{<:Finite},
    weights = false,
    descr   = "Bagging ensemble classifier."
    )

# ============================================================================
GradientBoostingRegressor_ = SKEN.GradientBoostingRegressor
@sk_reg mutable struct GradientBoostingRegressor <: MMI.Deterministic
    loss::String                    = "ls"::(_ in ("ls","lad","huber","quantile"))
    learning_rate::Float64          = 0.1::(_>0)
    n_estimators::Int               = 100::(_>0)
    subsample::Float64              = 1.0::(_>0)
    criterion::String               = "friedman_mse"::(_ in ("mse","mae","friedman_mse"))
    min_samples_split::Union{Int,Float64} = 2::(_>0)
    min_samples_leaf::Union{Int,Float64}  = 1::(_>0)
    min_weight_fraction_leaf::Float64     = 0.0::(_≥0)
    max_depth::Int                 = 3::(_>0)
    min_impurity_decrease::Float64 = 0.0::(_≥0)
    init::Any                      = nothing
    random_state::Any              = nothing
    max_features::Union{Int,Float64,String,Nothing} = nothing::(_===nothing || (isa(_,String) && (_ in ("auto","sqrt","log2"))) || _ > 0)
    alpha::Float64                 = 0.9::(_>0)
    verbose::Int                   = 0
    max_leaf_nodes::Option{Int}    = nothing::(_===nothing || _>0)
    warm_start::Bool               = false
    presort::Union{Bool,String}    = "auto"::(isa(_, Bool) || _ == "auto")
    validation_fraction::Float64   = 0.1::(_>0)
    n_iter_no_change::Option{Int}  = nothing
    tol::Float64                   = 1e-4::(_>0)
end
MMI.fitted_params(m::GradientBoostingRegressor, (f, _, _)) = (
    feature_importances = f.feature_importances_,
    train_score         = f.train_score_,
    loss                = f.loss_,
    init                = f.init_,
    estimators          = f.estimators_,
    oob_improvement     = m.subsample < 1 ? f.oob_improvement_ : nothing
    )

# ----------------------------------------------------------------------------
GradientBoostingClassifier_ = SKEN.GradientBoostingClassifier
@sk_clf mutable struct GradientBoostingClassifier <: MMI.Probabilistic
    loss::String                    = "deviance"::(_ in ("deviance","exponential"))
    learning_rate::Float64          = 0.1::(_>0)
    n_estimators::Int               = 100::(_>0)
    subsample::Float64              = 1.0::(_>0)
    criterion::String               = "friedman_mse"::(_ in ("mse","mae","friedman_mse"))
    min_samples_split::Union{Int,Float64} = 2::(_>0)
    min_samples_leaf::Union{Int,Float64}  = 1::(_>0)
    min_weight_fraction_leaf::Float64     = 0.0::(_≥0)
    max_depth::Int                 = 3::(_>0)
    min_impurity_decrease::Float64 = 0.0::(_≥0)
    init::Any                      = nothing
    random_state::Any              = nothing
    max_features::Union{Int,Float64,String,Nothing} = nothing::(_===nothing || (isa(_,String) && (_ in ("auto","sqrt","log2"))) || _ > 0)
    verbose::Int                   = 0
    max_leaf_nodes::Option{Int}    = nothing::(_===nothing || _>0)
    warm_start::Bool               = false
    presort::Union{Bool,String}    = "auto"::(isa(_, Bool) || _ == "auto")
    validation_fraction::Float64   = 0.1::(_>0)
    n_iter_no_change::Option{Int}  = nothing
    tol::Float64                   = 1e-4::(_>0)
end
MMI.fitted_params(m::GradientBoostingClassifier, (f, _, _)) = (
    n_estimators        = f.n_estimators_,
    feature_importances = f.feature_importances_,
    train_score         = f.train_score_,
    loss                = f.loss_,
    init                = f.init_,
    estimators          = f.estimators_,
    oob_improvement     = m.subsample < 1 ? f.oob_improvement_ : nothing
    )
metadata_model(GradientBoostingClassifier,
    input   = Table(Continuous),
    target  = AbstractVector{<:Finite},
    weights = false,
    descr   = "Gradient boosting ensemble classifier."
    )

# ============================================================================
RandomForestRegressor_ = SKEN.RandomForestRegressor
@sk_reg mutable struct RandomForestRegressor <: MMI.Deterministic
    n_estimators::Int              = 100::(_ > 0)
    criterion::String              = "mse"::(_ in ("mae", "mse"))
    max_depth::Option{Int}         = nothing::(_ === nothing || _ > 0)
    min_samples_split::Union{Int,Float64} = 2::(_ > 0)
    min_samples_leaf::Union{Int,Float64}  = 1::(_ > 0)
    min_weight_fraction_leaf::Float64     = 0.0::(_ ≥ 0)
    max_features::Union{Int,Float64,String,Nothing} = "auto"::(_ === nothing || (isa(_, String) && (_ in ("auto","sqrt","log2"))) || _ > 0)
    max_leaf_nodes::Option{Int}    = nothing::(_ === nothing || _ > 0)
    min_impurity_decrease::Float64 = 0.0::(_ ≥ 0)
    bootstrap::Bool                = true
    oob_score::Bool                = false
    n_jobs::Option{Int}            = nothing
    random_state::Any              = nothing
    verbose::Int                   = 0
    warm_start::Bool               = false
end
MMI.fitted_params(model::RandomForestRegressor, (f, _, _)) = (
    estimators          = f.estimators_,
    feature_importances = f.feature_importances_,
    n_features          = f.n_features_,
    n_outputs           = f.n_outputs_,
    oob_score           = model.oob_score ? f.oob_score_ : nothing,
    oob_prediction      = model.oob_score ? f.oob_prediction_ : nothing
    )
metadata_model(RandomForestRegressor,
    input   = Table(Count,Continuous),
    target  = AbstractVector{Continuous},
    weights = false,
    descr   = "Random forest regressor."
    )

# ----------------------------------------------------------------------------
RandomForestClassifier_ = SKEN.RandomForestClassifier
@sk_clf mutable struct RandomForestClassifier <: MMI.Probabilistic
    n_estimators::Int              = 100::(_ > 0)
    criterion::String              = "gini"::(_ in ("gini","entropy"))
    max_depth::Option{Int}         = nothing::(_ === nothing || _ > 0)
    min_samples_split::Union{Int,Float64} = 2::(_ > 0)
    min_samples_leaf::Union{Int,Float64}  = 1::(_ > 0)
    min_weight_fraction_leaf::Float64     = 0.0::(_ ≥ 0)
    max_features::Union{Int,Float64,String,Nothing} = "auto"::(_ === nothing || (isa(_, String) && (_ in ("auto","sqrt","log2"))) || _ > 0)
    max_leaf_nodes::Option{Int}    = nothing::(_ === nothing || _ > 0)
    min_impurity_decrease::Float64 = 0.0::(_ ≥ 0)
    bootstrap::Bool                = true
    oob_score::Bool                = false
    n_jobs::Option{Int}            = nothing
    random_state::Any              = nothing
    verbose::Int                   = 0
    warm_start::Bool               = false
    class_weight::Any              = nothing
end
MMI.fitted_params(m::RandomForestClassifier, (f, _, _)) = (
    estimators            = f.estimators_,
    classes               = f.classes_,
    n_classes             = f.n_classes_,
    n_features            = f.n_features_,
    n_outputs             = f.n_outputs_,
    feature_importances   = f.feature_importances_,
    oob_score             = m.oob_score ? f.oob_score_ : nothing,
    oob_decision_function = m.oob_score ? f.oob_decision_function_ : nothing
    )
metadata_model(RandomForestClassifier,
    input   = Table(Count,Continuous),
    target  = AbstractVector{<:Finite},
    weights = false,
    descr   = "Random forest classifier."
    )

const ENSEMBLE_REG = Union{Type{<:AdaBoostRegressor}, Type{<:BaggingRegressor}, Type{<:GradientBoostingRegressor}}

MMI.input_scitype(::ENSEMBLE_REG)  = Table(Continuous)
MMI.target_scitype(::ENSEMBLE_REG) = AbstractVector{Continuous}

# ============================================================================
ExtraTreesRegressor_ = SKEN.ExtraTreesRegressor
@sk_reg mutable struct ExtraTreesRegressor <: MMI.Deterministic
    n_estimators::Int              = 100::(_>0)
    criterion::String              = "mse"::(_ in ("mae", "mse"))
    max_depth::Option{Int}         = nothing::(_ === nothing || _ > 0)
    min_samples_split::Union{Int,Float64}  = 2::(_ > 0)
    min_samples_leaf::Union{Int,Float64}   = 1::(_ > 0)
    min_weight_fraction_leaf::Float64      = 0.0::(_ ≥ 0)
    max_features::Union{Int,Float64,String,Nothing} = "auto"::(_ === nothing || (isa(_, String) && (_ in ("auto","sqrt","log2"))) || _ > 0)
    max_leaf_nodes::Option{Int}    = nothing::(_ === nothing || _ > 0)
    min_impurity_decrease::Float64 = 0.0::(_ ≥ 0)
    bootstrap::Bool                = true
    oob_score::Bool                = false
    n_jobs::Option{Int}            = nothing
    random_state::Any              = nothing
    verbose::Int                   = 0
    warm_start::Bool               = false
end
MMI.fitted_params(m::ExtraTreesRegressor, (f, _, _)) = (
    estimators          = f.estimators_,
    feature_importances = f.feature_importances_,
    n_features          = f.n_features_,
    n_outputs           = f.n_outputs_,
    oob_score           = m.oob_score ? f.oob_score_ : nothing,
    oob_prediction      = m.oob_score ? f.oob_prediction_ : nothing,
    )
metadata_model(ExtraTreesRegressor,
    input   = Table(Continuous),
    target  = AbstractVector{Continuous},
    weights = false,
    descr   = "Extra trees regressor, fits a number of randomized decision trees on various sub-samples of the dataset and uses averaging to improve the predictive accuracy and control over-fitting."
    )

# ----------------------------------------------------------------------------
ExtraTreesClassifier_ = SKEN.ExtraTreesClassifier
@sk_clf mutable struct ExtraTreesClassifier <: MMI.Probabilistic
    n_estimators::Int              = 100::(_>0)
    criterion::String              = "gini"::(_ in ("gini", "entropy"))
    max_depth::Option{Int}         = nothing::(_ === nothing || _ > 0)
    min_samples_split::Union{Int,Float64}  = 2::(_ > 0)
    min_samples_leaf::Union{Int,Float64}   = 1::(_ > 0)
    min_weight_fraction_leaf::Float64      = 0.0::(_ ≥ 0)
    max_features::Union{Int,Float64,String,Nothing} = "auto"::(_ === nothing || (isa(_, String) && (_ in ("auto","sqrt","log2"))) || _ > 0)
    max_leaf_nodes::Option{Int}    = nothing::(_ === nothing || _ > 0)
    min_impurity_decrease::Float64 = 0.0::(_ ≥ 0)
    bootstrap::Bool                = true
    oob_score::Bool                = false
    n_jobs::Option{Int}            = nothing
    random_state::Any              = nothing
    verbose::Int                   = 0
    warm_start::Bool               = false
    class_weight::Any              = nothing
end
MMI.fitted_params(m::ExtraTreesClassifier, (f, _, _)) = (
    estimators            = f.estimators_,
    classes               = f.classes_,
    n_classes             = f.n_classes_,
    feature_importances   = f.feature_importances_,
    n_features            = f.n_features_,
    n_outputs             = f.n_outputs_,
    oob_score             = m.oob_score ? f.oob_score_ : nothing,
    oob_decision_function = m.oob_score ? f.oob_decision_function_ : nothing,
    )
metadata_model(ExtraTreesClassifier,
    input   = Table(Continuous),
    target  = AbstractVector{<:Finite},
    weights = false,
    descr   = "Extra trees classifier, fits a number of randomized decision trees on various sub-samples of the dataset and uses averaging to improve the predictive accuracy and control over-fitting."
    )

ext: jl


text:
#Written By Johnathan Bizzano
export MultiVarPolynomial

struct MultiVarPolynomial{T, N} <: AbstractMathFun{T, N}
    summation::RegressableSummation{T, N}

    function MultiVarPolynomial{T, N}(powers::AbstractArray) where {T, N}
        new{T, N}(Functions.RegressableSummation{T, N}(GeneratePowers(powers, T, N)))
    end

    function GeneratePowers(iter::AbstractArray, T::Type, Var_Len)
        iter = adjoint(iter)
        arr = TermsArray{T, Var_Len}(undef, size(iter, 1))
        i = 1
        vars = [Functions.VarFunction{T}(v) for v in 1:Var_Len]
        for row in eachrow(iter)
            powCollection = MFCollection(undef, Var_Len)
            for col in 1:Var_Len
                powCollection[col] = Functions.PowFunction{T, Var_Len}(vars[col], row[col])
            end
            arr[i] = Functions.ProductFunction{T, Var_Len}(powCollection)
            i += 1
        end
        arr
    end

    (f::MultiVarPolynomial{T, N})(args::AbstractArray{T}; Data_Type::Type=T) where {T, N} = f.summation(args, Data_Type = Data_Type)
    Functions.train!(f::MultiVarPolynomial{T, N}, data::AbstractDataSet; Data_Type::Type=T) where {T, N} = Functions.train!(f.summation, data, Data_Type = Data_Type)
    Base.sizeof(x::MultiVarPolynomial) = sizeof(x.summation)
    Base.string(x::MultiVarPolynomial) = string(x.summation)
    Functions.simplify(x::MultiVarPolynomial) = Functions.simplify(x.summation)
    Functions.trainer(f::MultiVarPolynomial{T, N}; precision=10) where {T, N} = RegressableSummationTrainer{T, N}(f.summation, precision = precision)
end


text:
using MySQL
using DataFrames

"""
    get_semantic_occurrences_df(db, mesh, umls_concepts...)

Given a mesh descriptor and filtering UMLS concepts, returns a DataFrame with
PMIDs and related mesh descriptor filtered by the concepts.
"""
function get_semantic_occurrences_df(db::MySQL.Connection, mesh::String, umls_concepts::String...)

    concept_string = "'" * join(umls_concepts,"','") * "'"

    # select pmids from medline which use the mesh descriptor
    query_string = """ SELECT mh.pmid, mrc.str as descriptor
                    FROM pubmed_comorbidities.MH_SMALL mh

                    JOIN (select uid from pubmed_comorbidities.ALL_MESH where STR = '$mesh') md
                    ON md.uid <> mh.desc_uid

                    JOIN pubmed_comorbidities.ALL_MESH mrc
                    ON mrc.uid = mh.desc_uid

                    JOIN pubmed_comorbidities.MRSTY_FOR_MESH mrs
                    ON mrc.cui = mrs.cui

                    JOIN pubmed_comorbidities.MH_SMALL mh2
                    ON mh2.pmid = mh.pmid

                    join (select uid from pubmed_comorbidities.ALL_MESH where STR = '$mesh') md2
                    on md2.uid = mh2.desc_uid

                    WHERE mrs.sty in ($concept_string); """;

    println(query_string)

    articles_df = MySQL.query(db, query_string, DataFrame)

    unique!(articles_df)

    return articles_df

end


text:
__precompile__()

module RICoverAlg

# package code goes here

export preprocessinput,convertInputsToFloats, cancelZeroStates, checkForTriviality, transform_utilities
export do_cover_iterations,cancel_zero_alternatives,run_cover_algorithm,compute_factor
export risolve
# export checkOptimalityOfp

include("preprocessinput.jl")
include("cover_algorithm.jl")

# function checkOptimalityOfp(U,g,lmd,p)
#       Z = exp.(U./lmd)
#       ϵ_max = lmd * log(findmax(compute_factor(Z,g,p))[1])
#       payoff = lmd * ((log.(p'*Z))*g)[1]
#       return payoff, ϵ_max, payoff/(ϵ_max+payoff)
# end

function risolve(U, g, lmd;
                  initialIterations::Int  = 10,
                  aTolSetChoiceProbToZero::Float64 = sqrt(eps()),
                  aTolerance4Stopping::Float64 = 0.0,
                  maxIterations::Int = 9999)
      # solution = (optimalProbabilityVector, Maximal Payoff, Number of Iterations, Message)
      ϵ_max = 0
      Z, g, shiftterm, solution = preprocessinput(U, g, lmd)

      if solution[4] == "Run Cover Algorithm"
            ## running cover algorithm
            solution = run_cover_algorithm(Z, g,
                                          initialIterations  = initialIterations,
                                          aTolSetChoiceProbToZero = aTolSetChoiceProbToZero,
                                          aTolerance4Stopping = aTolerance4Stopping,
                                          maxIterations = maxIterations)
            ## computing optimality bounds
            ϵ_max = lmd * log(findmax(compute_factor(Z,g,solution[1]))[1])
      end

      payoff = lmd*(solution[2])+shiftterm
      mention_if_beyond_maxIterations = ""
      if solution[3] >= maxIterations
            mention_if_beyond_maxIterations = ", i.e., we reached maxIterations so that we aborted"
      end
      output_message = solution[4] * "; Number of iterations: $(solution[3])"*
                        mention_if_beyond_maxIterations*
                        "; Resulting payoff $payoff is $(payoff/(ϵ_max+payoff))"*
                        " optimal (maximal payoff <= $(ϵ_max+payoff))"
      return (solution[1],payoff,output_message)
end

end # module


text:
module TestCLBlast

using Compat
using Compat.Test
using Compat.Random

using CLBlast, OpenCL

const n_L1 = 64
const m_L2 = 60
const n_L2 = 50
const kl = 2
const ku = 3
const m_L3 = 6
const n_L3 = 5
const k_L3 = 4
@compat const eltypes = (Float32, Float64, ComplexF32, ComplexF64)

device, ctx, queue = cl.create_compute_context()

@time @testset "BLAS Level 1" begin include("L1_test.jl") end
@time @testset "BLAS Level 2" begin include("L2_test.jl") end
@time @testset "BLAS Level 3" begin include("L3_test.jl") end
@time @testset "Auxiliary Tests" begin include("auxiliary_test.jl") end

end # module


text:
using Compat: uninitialized, Nothing, Cvoid, AbstractDict
using Images
using ImageMagick
using DataStructures
using Distances

displayInYellow(s) = string("\x1b[93m",s,"\x1b[0m")
displayInRed(s) = string("\x1b[91m",s,"\x1b[0m")

include("lib_window.jl")
include("lib_opengl.jl")
include("lib_math.jl")
include("lib_time.jl")

"""
TODO
"""
function waitForFileReady(path::String, func::Function, tryCount=100, tryWait=0.1)
	result=false
	for i = 1:tryCount

		#try reading file
		if stat(path).size > 0
			open(path) do file
				 result=func(file)
			end
			if result break end
		end

		Libc.systemsleep(tryWait) #wait
	end
	result
end

"""
TODO
"""
function fileGetContents(path::String, tryCount=100, tryWait=0.1)
	content=nothing
	waitForFileReady(path,(x)->(content=readstring(x); content != nothing),tryCount,tryWait)
	content
end

TITLE = "Julia OpenGL"
STARTTIME = Dates.time()
PREVTIME = STARTTIME
FRAMES = 0
MAX_FRAMES = 0
FPS = 0
MAX_FPS = 0
ITERATION = 0
BLOCK_COUNT = 0
SIZE = 0

"""
TODO
"""
function UpdateCounters()
  UpdateTimers()
  showFrames()
end

prevTime = Ref(0.0)

"""
TODO
"""
function showFrames()
  global TITLE, TIMERS, FRAMES, MAX_FRAMES, FPS, MAX_FPS, ITERATION, BLOCK_COUNT, PREVTIME, RENDER_METHOD

  time = Dates.time() #GetTimer("FRAME_TIMER")

  ITERATION +=1
  if !OnTime(1.0, prevTime, time) FRAMES += 1; return end

  #FPS = FRAMES/(time - PREVTIME)
  #PREVTIME = time
  #if MAX_FPS < FPS MAX_FPS = FPS end
  #if FPS > 15 COUNT += 1 end
  #fpms = FPS > 0 ? (1000.0 / FPS) : 0
  #max_fmps = MAX_FPS > 0 ? (1000.0 / MAX_FPS) : 0
  #norm_fps = FPS/MAX_FPS

  if MAX_FRAMES < FRAMES MAX_FRAMES = FRAMES end
  const fps = FRAMES
  const max_fps = MAX_FRAMES
  const fpms = FRAMES > 0 ? (1000.0 / FRAMES) : 0
  const max_fmps = MAX_FRAMES > 0 ? (1000.0 / MAX_FRAMES) : 0
  const norm_fps = FRAMES / MAX_FRAMES

  GLFW.SetWindowTitle(window, "$(TITLE) - FPS $(round(fps, 2))[$(round(max_fps, 2))] | FMPS $(round(fpms, 2))[$(round(max_fmps, 2))] - Blocks $CHUNK_SIZE^3 ($BLOCK_COUNT) - IT $ITERATION")
  FRAMES = 0
end

include("cubeData.jl")
include("camera.jl")
include("frustum.jl")
include("chunk.jl")
include("mesh.jl")
include("texture.jl")

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
__precompile__()

module SCIP

if isfile(joinpath(dirname(@__FILE__),"..","deps","deps.jl"))
    include("../deps/deps.jl")
else
    error("SCIP.jl not properly installed. Please run Pkg.build(\"SCIP\")")
end

importall MathProgBase.SolverInterface

include("types.jl")
include("wrapper.jl")
include("mpb_interface.jl")
include("params.jl")

function CSIPversion()
    "$(_majorVersion()).$(_minorVersion()).$(_patchVersion())"
end

function __init__()
    version = VersionNumber(CSIPversion())
    if !(version.major == 0 && version.minor == 3 && version.patch == 4)
        depsdir = realpath(joinpath(dirname(@__FILE__),"..","deps"))
        error("Current CSIP version installed is $(version), but we require version 0.3.4. On Linux, delete the contents of the `$depsdir` directory except for `build.jl`, then rerun Pkg.build(\"SCIP\").")
    end
end

end


text:
module Catch22
using catch22_jll
using DimensionalData
using Libdl
import Statistics.mean, Statistics.std

include("Feature.jl")
include("FeatureSet.jl")
include("FeatureArray.jl")
include("metadata.jl")
include("testdata.jl")

catch22_jll.__init__() # Initialise the c-library

zscore(𝐱::AbstractVector) = (𝐱 .- mean(𝐱))./(std(𝐱))

"""
    _catch22(𝐱::AbstractArray{Float64}, fName::Symbol)
    _catch22(fName::Symbol, 𝐱::AbstractArray{Float64})
Evaluate the feature `fName` on the time series `𝐱`. If an array is supplied, features are calculated for its columns and returned as a Vector. See `Catch22.featuredescriptions` for a summary of the 22 available time series features.

# Examples
```julia-repl
𝐱 = Catch22.testdata[:test]
Catch22._catch22(𝐱, :DN_HistogramMode_5)
```
"""
function _catch22(𝐱::AbstractVector, fName::Symbol)::Float64
    if any(isinf.(𝐱)) || any(isnan.(𝐱))
        return NaN
    end
    𝐱 = zscore(𝐱)
    𝐱 = convert(Vector{Float64}, 𝐱)
    fType = featuretypes[fName]
    if fType <: AbstractFloat
        ccall(dlsym(dlopen(ccatch22), fName), Cdouble, (Ptr{Array{Cdouble}},Cint), 𝐱, Int(size(𝐱, 1)))
    elseif fType <: Integer
        ccall(dlsym(dlopen(ccatch22), fName), Cint, (Ptr{Array{Cdouble}},Cint), 𝐱, Int(size(𝐱, 1)))
    end
end
function _catch22(X::AbstractArray{Float64, 2}, fName::Symbol)::AbstractArray{Float64, 2}
    mapslices(𝐱 -> _catch22(𝐱, fName), X, dims=[1])
end

"""
    catch22(𝐱::Vector)
    catch22(X::Array)
    catch22[featurename::Symbol](X::Array)
Evaluate all features for a time series vector `𝐱` or the columns of an array `X`.
`catch22` is a FeatureSet, which means it can be indexed by feature names (as symbols) to return a subset of the available features.
`getnames(catch22)`, `getkeywords(catch22)` and `getdescriptions(catch22)`` will also return feature names, keywords and descriptions respectively.
Features are returned in a `FeatureArray`, in which array rows are annotated by feature names. A `FeatureArray` can be converted to a regular array with `Array(F)`.

# Examples
```julia-repl
𝐱 = Catch22.testdata[:test]
𝐟 = catch22(𝐱)

X = randn(100, 10)
F = catch22(X)
F = catch22[:DN_HistogramMode_5](X)
```
"""
catch22 = FeatureSet([(x -> _catch22(x, f)) for f ∈ featurenames], featurenames, featurekeywords, featuredescriptions)
export catch22

for f = featurenames
    eval(quote
        $f = catch22[$(Meta.quot(f))]; export $f
    end)
end

"""
    DN_HistogramMode_5(x::AbstractVector{Union{Float64, Int}}) # For example
An alternative to `catch22(...)``. All features, such as `DN_HistogramMode_5`, are exported as Features and can be evaluated by calling their names.

# Examples
```julia-repl
𝐱 = Catch22.testdata[:test]
f = DN_HistogramMode_5(𝐱)
```
"""
DN_HistogramMode_5;

end


text:
module Shazam
include("AudioDisplay.jl")
#import Pkg;
#Pkg.add("FFTW")
#Pkg. add("WAV")

using Plots
using FFTW
using WAV

export AB1Aufgabe1a
export AB1A1aPlay
export AB1A1bPlay
export AB1Aufgabe1b
export AB1A1c1Play
export AB1Aufgabe1c1
export AB1A1c2Play
export AB1Aufgabe1c2
export AB1A1dPlay
export AB1Aufgabe1d
export AB1A2aPlay
export AB1Aufgabe2a
export AB1A2bPlay
export AB1Aufgabe2b
export AB1ZAPlay
export AB1ZusatzaufgabeA
export AB1ZBPlay
export AB1ZusatzaufgabeB
export AB1ZCPlay
export AB1ZusatzaufgabeC
export AB1ZDPlay
export AB1ZusatzaufgabeD
export AB2A1Play
export AB2Aufgabe1
export AB2Aufgabe2b1
export AB2Aufgabe2b2
export AB2A31Play
export AB2Aufgabe31
export AB2A32Play
export AB2Aufgabe32
export AB2Zusatzaufgabe
export AB3Aufgabe1
export AB4Aufgabe1a
export AB4Aufgabe1b
export AB4Aufgabe1c

function AB1A1aPlay()
    eps = 1/44100
    x = range(0.0,stop=2,step=eps);
    y = zeros(Float64,length(x));
    for i = 1:length(y)
        y[i] = sin(2*pi*200*x[i]);
    end
    #display(y)
    wavwrite(y, "SinustonAB1A1a.wav", Fs=44100);
    #y, fs = wavread("SinustonAB1A1a.wav");
	#wavplay(y, fs)
	inline_audioplayer("SinustonAB1A1a.wav")
    plot(y)
end

function AB1Aufgabe1a(amplitude,frequenz)
    if amplitude!=1 || frequenz!=200
		@warn "Amplitude und/oder Frequenz falsch!"
    end
    x = 0:0.0001:0.1;
    y = zeros(1001);
    for i = 1:1001
        y[i] = amplitude*sin(2*pi*frequenz*x[i]);
    end
    labels = ["dein Ton"];
    plot(x,y, label = labels)
end

function AB1A1bPlay(amplitude_leiser, amplitude_lauter)
    eps = 1/44100
    x = range(0.0,stop=2,step=eps);
    w = zeros(Float64,length(x));
    for i = 1:length(w)
        w[i] = sin(2*pi*200*x[i]);
    end
    y = zeros(Float64,length(x));
    for i = 1:length(y)
        y[i] = amplitude_leiser*sin(2*pi*200*x[i]);
    end
    z = zeros(Float64,length(x));
    for i = 1:length(z)
        z[i] = amplitude_lauter*sin(2*pi*200*x[i]);
    end
    ton  = [w; y; z];
    #display(ton)
    wavwrite(ton, "SinustonAB1A1b.wav", Fs=44100);
    #ton, fs = wavread("SinustonAB1A1b.wav");
    #wavplay(ton, fs)
	inline_audioplayer("SinustonAB1A1b.wav")

end

function AB1Aufgabe1b(amplitude,amplitude_leiser,amplitude_lauter)
    if amplitude_leiser >= amplitude
		@warn "Der erste Ton ist nicht leiser als der Ton aus Teil a)!"
    end
    if amplitude_lauter <= amplitude
		@warn "Der zweite Ton ist nicht lauter als der Ton aus Teil a)!"
    end
    x = 0:0.0001:0.1;
    y = zeros(1001,3);
    for i = 1:1001
        y[i,1] = sin(2*pi*200*x[i]);
        y[i,2] = amplitude_leiser*sin(2*pi*200*x[i]);
        y[i,3] = amplitude_lauter*sin(2*pi*200*x[i]);
    end
    labels = ["Teil a)" "leiser" "lauter"];
    plot(x,y, label = labels)
end

function AB1A1c1Play(A)
    eps = 1/44100
    x = range(0.0,stop=2,step=eps);
    y = zeros(Float64,length(x));
    for i = 1:length(y)
        y[i] = A(x[i])*sin(2*pi*200*x[i]);
    end
    #display(ton)
    wavwrite(y, "SinustonAB1A1c1.wav", Fs=44100);
    #y, fs = wavread("SinustonAB1A1c1.wav");
    #wavplay(y, fs)
	inline_audioplayer("SinustonAB1A1c1.wav")
end

function AB1Aufgabe1c1(A)
    if A(0) != 0
		@warn "Die Amplitude hat zur Zeit t = 0 nicht den Wert 0!"
    end
    m1 = A(1)-A(0);
    m2 = A(2)-A(1);
    if m1 != m2
		@warn "Die Funktion A besitzt kein lineares Wachstum!"
    end
    labels = ["dein Ton"];
    x = 0:0.0001:0.1;
    y = zeros(1001);
    for i = 1:1001
        y[i] = A(x[i])*sin(2*pi*200*x[i]);
    end
    plot(x,y, label = labels)
end

function AB1A1c2Play(A, A2)
    eps = 1/44100
    x = range(0.0,stop=2,step=eps);
    y = zeros(Float64,length(x));
    for i = 1:length(y)
        y[i] = A(x[i])*sin(2*pi*200*x[i]);
    end
    z = zeros(Float64,length(x));
    for i = 1:length(z)
        z[i] = A2(x[i])*sin(2*pi*200*x[i]);
    end
    w = zeros(44100,1);
    ton = [y;w;z];
    #display(ton)
    wavwrite(ton, "SinustonAB1A1c2.wav", Fs=44100);
    #y, fs = wavread("SinustonAB1A1c2.wav");
    #wavplay(y, fs)
	inline_audioplayer("SinustonAB1A1c2.wav")
end

function AB1Aufgabe1c2(A,A2)
    if A2(0) != 0
		@warn "Die Amplitude hat bei t = 0 nicht den Wert 0!"
    end
    m1 = A2(1)-A2(0);
    m2 = A2(2)-A2(1);
    if m1 != m2
		@warn "Die Funktion A2 besitzt kein lineares Wachstum!"
    end
    if A(0.1) > A2(0.1)
		@warn "Die Funktion A2 wächst nicht schneller als die Funktion A!"
    end
    labels = ["Teil 1" "Teil 2"];
    x = 0:0.0001:0.1;
    y = zeros(1001,2);
    for i = 1:1001
        y[i,1] = A(x[i])*sin(2*pi*200*x[i]);
        y[i,2] = A2(x[i])*sin(2*pi*200*x[i]);
    end
    plot(x,y,label=labels)
end

function AB1A1dPlay(A)
    eps = 1/44100
    x = range(0.0,stop=1,step=eps);
    y = zeros(Float64,length(x));
    for i = 1:length(y)
        y[i] = A(x[i])*sin(2*pi*200*x[i]);
    end
    #display(y)
    wavwrite(y, "SinustonAB1A1d.wav", Fs=44100);
    #y, fs = wavread("SinustonAB1A1d.wav");
    #wavplay(y, fs)
	inline_audioplayer("SinustonAB1A1d.wav")
end

function AB1Aufgabe1d(A3)
    if A3(1) != 0
		@warn "Bei t = 1 hat die Amplitude nicht den Wert 0!"
    end
    m1 = A3(1)-A3(0);
    m2 = A3(2)-A3(1);
    if abs(m1 - m2) > 0.0000000001
		@warn "Die Funktion A3 besitzt kein lineares Abklingverhalten!"
    end
    labels = ["dein Ton"];
    x = 0:0.0001:1;
    y = zeros(length(x));
    for i = 1:length(y)
        y[i] = A3(x[i])*sin(2*pi*200*x[i]);
    end
    plot(x,y, label = labels)
end

function AB1A2aPlay(F_hoeher,F_tiefer)
    eps = 1/44100
    x = range(0.0,stop=2,step=eps);
    y = zeros(Float64,length(x));
    for i = 1:length(y)
        y[i] = sin(2*pi*200*x[i]);
    end
    z = zeros(Float64,length(x));
    for i = 1:length(z)
        z[i] = sin(F_hoeher*2*pi*200*x[i]);
    end
    w = zeros(Float64,length(x));
    for i = 1:length(w)
        w[i] = sin(F_tiefer*2*pi*200*x[i]);
    end
    ton = [y; z; w];
    #display(ton)
    wavwrite(ton, "SinustonAB1A2a.wav", Fs=44100);
    #ton, fs = wavread("SinustonAB1A2a.wav");
    #wavplay(ton, fs)
	inline_audioplayer("SinustonAB1A2a.wav")
end

function AB1Aufgabe2a(F_hoeher,F_tiefer)
    if F_hoeher >= 6 || F_hoeher <= 0
		@warn "F_hoeher liegt nicht zwischen 0 und 6!"
    end
    if F_tiefer >= 6 || F_tiefer <= 0
		@warn "F_tiefer liegt nicht zwischen 0 und 6!"
    end
    if F_hoeher <= 1
		@warn "Der erste Ton ist nicht höher als der Referenzton!"
    end
    if F_tiefer >= 1
		@warn "Der zweite Ton ist nicht tiefer als der Referenzton!"
    end
    x = 0:0.00001:0.01;
    y = zeros(length(x),3);
    for i = 1:length(x)
        y[i,1] = sin(2*pi*200*x[i]);
        y[i,2] = sin(F_hoeher*2*pi*200*x[i]);
        y[i,3] = sin(F_tiefer*2*pi*200*x[i]);
    end
    labels = ["Referenz" "hoeher" "tiefer"];
    plot(x,y, label = labels)
end

function AB1A2bPlay(F)
    eps = 1/44100
    x = range(0.0,stop=1,step=eps);
    y = zeros(Float64,length(x));
    for i = 1:length(y)
        y[i] = sin(F(x[i])*2*pi*200*x[i]);
    end
    #display(y)
    wavwrite(y, "SinustonAB1A2b.wav", Fs=44100);
    #y, fs = wavread("SinustonAB1A2b.wav");
    #wavplay(y, fs)
	inline_audioplayer("SinustonAB1A2b.wav")
end

function AB1Aufgabe2b(F)
    if F(0) != 1
		@warn "Zum Zeitpunkt t = 0 hat die Frequenz von g nicht den Wert 200 Hz!"
    end
    if F(0.1) != 2
		@warn "Zum Zeitpunkt t = 0.1 hat die Frequenz von g nicht den Wert 400Hz!"
    end
    m1 = F(1)-F(0);
    m2 = F(2)-F(1);
    if abs(m1-m2) > 0.00000001
		@warn "Die Funktion F besitzt kein lineares Wachstum!"
    end
    labels = ["dein Signal"];
    x = 0:0.0001:0.1;
    y = zeros(1001);
    for i = 1:1001
        y[i] = sin(F(x[i])*2*pi*200*x[i]);
    end
    plot(x,y, label = labels)
end

function AB1ZAPlay(A,a,b)
    eps = 1/44100
    a_ton = a - 1;
    b_ton = b + 1;
    x = range(a_ton,stop=b_ton,step=eps);
    y = zeros(Float64,length(x));
    for i = 1:length(y)
        y[i] = A(x[i])*sin(2*pi*200*x[i]);
    end
    #display(y)
    wavwrite(y, "SinustonAB1A2b.wav", Fs=44100);
    #y, fs = wavread("SinustonAB1A2b.wav");
    #wavplay(y, fs)
	inline_audioplayer("SinustonAB1A2b.wav")
end

function AB1ZusatzaufgabeA(A4,t_min,a,b)
    if t_min <= a || t_min >= b
		@warn " t_min liegt nicht zwischen a und b!"
    end
    min = 1;
    for i = -10:10
        if A4(t_min + i/100) < A4(t_min)
            min = 0;
        end
    end
    if min == 0
		@warn "Bei t = t_min ist die Amplitude nicht minimal!"
    end
    labels = ["dein Ton"];
    if t_min > a && t_min < b
        x = a:(b-a)/10000:b;
        y = zeros(length(x));
        for i = 1:length(y)
            y[i] = A4(x[i])*sin(2*pi*200*x[i]);
        end
        plot(x,y, label = labels)
    end
end

function AB1ZBPlay(A,a,b)
    eps = 1/44100
    x = range(a,stop=b,step=eps);
    y = zeros(Float64,length(x));
    for i = 1:length(y)
        y[i] = A(x[i])*sin(2*pi*50*x[i]);
    end
    #display(y)
    wavwrite(y, "SinustonAB1ZB.wav", Fs=44100);
    #y, fs = wavread("SinustonAB1ZB.wav");
    #wavplay(y, fs)
	inline_audioplayer("SinustonAB1ZB.wav")
end

function AB1ZusatzaufgabeB(A5,t_max,a,b)
    if t_max <= a || t_max >= b
		@warn "t_max liegt nicht zwischen a und b!"
    end
    if A5(a) != 0
		@warn "Bei t = a hat die Amplitude nicht den Wert 0!"
    end
    if A5(b) != 0
		@warn "Bei t = b hat die Amplitude nicht den Wert 0!"
    end
    max = 1;
    for i = -10:10
        if A5(t_max + i/100) > A5(t_max)
            max = 0;
        end
    end
    if max == 0
		@warn "Bei t = t_max ist die Amplitude nicht maximal!"
    end
    labels = ["dein Ton"];
    if t_max > a && t_max < b
        x = a:(b-a)/10000:b;
        y = zeros(length(x));
        for i = 1:length(y)
            y[i] = A5(x[i])*sin(2*pi*50x[i]);
        end
        plot(x,y, label = labels)
    end
end

function AB1ZCPlay(a,F)
    eps = 1/44100
    b = a + 1;
    x = range(a,stop=b,step=eps);
    y = zeros(Float64,length(x));
    for i = 1:length(y)
        y[i] = sin(F(x[i])*2*pi*200*x[i]);
    end
    #display(y)
    wavwrite(y, "SinustonAB1ZC.wav", Fs=44100);
    #y, fs = wavread("SinustonAB1ZC.wav");
    #wavplay(y, fs)
	inline_audioplayer("SinustonAB1ZC.wav")
end

function AB1ZusatzaufgabeC(a,b,F)
    if b <= a
		@warn "b muss größer sein als a!"
    end
    increase = 1;
    for i = 1:100
        if F(a+(b-a)/100*i) < F(a+(b-a)/100*(i+1))
            increase = 0;
        end
    end
    if increase == 0
		@warn "Die Frequenz fällt nicht auf dem Intervall [a,b]!"
    end
    labels = ["dein Signal"];
    if a < b
        x = a:(b-a)/10000:b;
        y = zeros(length(x));
        for i = 1:length(y)
            y[i] = sin(F(x[i])*2*pi*200*x[i]);
        end
        plot(x,y, label = labels)
    end
end

function AB1ZDPlay(a,b,F)
    eps = 1/44100
    x = range(a,stop=b,step=eps);
    y = zeros(Float64,length(x));
    for i = 1:length(y)
        y[i] = sin(F(x[i])*2*pi*200*x[i]);
    end
    #display(y)
    wavwrite(y, "SinustonAB1ZD.wav", Fs=44100);
    #y, fs = wavread("SinustonAB1ZD.wav");
    #wavplay(y, fs)
	inline_audioplayer("SinustonAB1ZD.wav")
end

function AB1ZusatzaufgabeD(a,b,F)
    if b <= a
		@warn "b muss größer sein als a!"
    end
    labels = ["dein Signal"];
    if a < b
        x = a:(b-a)/10000:b;
        y = zeros(length(x));
        for i = 1:length(y)
            y[i] = sin(F(x[i])*2*pi*200*x[i]);
        end
        plot(x,y, label = labels)
    end
end

function AB2A1Play(g)
    eps = 1/44100
    x = range(0.0,stop=2,step=eps);
    y = zeros(Float64,length(x));
    for i = 1:length(y)
        y[i] = g(x[i]);
    end
    #display(y)
    wavwrite(y, "SinustonAB2A1.wav", Fs=44100);
    #y, fs = wavread("SinustonAB2A1.wav");
    #wavplay(y, fs)
	inline_audioplayer("SinustonAB2A1.wav")
end

function AB2Aufgabe1(g)
    x = 0:0.0001:0.1;
    y = zeros(length(x));
    for i = 1:length(y)
        y[i] = g(x[i]);
    end
    labels = ["dein Dreiklang"];
    plot(x,y, label = labels)
end

function AB2Aufgabe2b1(g)
    samplingrate = 44100;
    x = 0:1/44100:6;
    y = zeros(length(x));
    for i = 1:length(y)
        y[i] = g(x[i]);
    end
    yfft = fft(y);
    f = zeros(length(yfft));
    for i = 1:length(yfft)
        f[i] = (i-1)*samplingrate/length(yfft);
    end
    fplot = f[1:5000];
    yplot = yfft[1:5000];
    yfind = abs.(yplot);
    m1, ind1 = findmax(yfind);
    yfind[ind1] = 0;
    m2, ind2 = findmax(yfind);
    yfind[ind2] = 0;
    m3, ind3 = findmax(yfind);
    labels = ["Frequenzspektrum"];
    plot(fplot,abs.(yplot)/m3, label = labels)
end

function AB2Aufgabe2b2(g)
    samplingrate = 44100;
    x = 0:1/44100:6;
    y = zeros(length(x));
    for i = 1:length(y)
        y[i] = g(x[i]);
    end
    yfft = fft(y);
    f = zeros(length(yfft));
    for i = 1:length(yfft)
        f[i] = (i-1)*samplingrate/length(yfft);
    end
    fplot = f[1:5000];
    yplot = yfft[1:5000];
    yfind = abs.(yplot);
    m1, ind1 = findmax(yfind);
    yfind[ind1] = 0;
    m2, ind2 = findmax(yfind);
    yfind[ind2] = 0;
    m3, ind3 = findmax(yfind);
    print("Frequenz 1 = ", f[ind1], "     Frequenz 2 = ", f[ind2], "     Frequenz 3 = ", f[ind3])
end

function AB2A31Play(y1,y2,y3)
    eps = 1/44100
    x = range(0.0,stop=2,step=eps);
    y = zeros(Float64,length(x));
    for i = 1:length(y)
        y[i] = y1(x[i]);
    end
    z = zeros(Float64,length(x));
    for i = 1:length(z)
        z[i] = y2(x[i]);
    end
    w = zeros(Float64,length(x));
    for i = 1:length(w)
        w[i] = y3(x[i]);
    end
    ton = [y; z; w];
    #display(ton)
    wavwrite(ton, "SinustonAB2A31.wav", Fs=44100);
    #ton, fs = wavread("SinustonAB2A31.wav");
    #wavplay(ton, fs)
	inline_audioplayer("SinustonAB2A31.wav")
end

function AB2Aufgabe31(y1,y2,y3)
    x = 0:0.00001:0.01;
    y = zeros(length(x),3);
    for i = 1:length(x)
        y[i,1] = y1(x[i]);
        y[i,2] = y2(x[i]);
        y[i,3] = y3(x[i]);
    end
    labels = ["y1" "y2" "y3"];
    plot(x,y, label = labels)
end

function AB2A32Play(S)
    eps = 1/44100
    x = range(0.0,stop=2,step=eps);
    y = zeros(Float64,length(x));
    for i = 1:length(y)
        y[i] = S(x[i]);
    end
    #display(y)
    wavwrite(y, "SinustonAB2A32.wav", Fs=44100);
    #y, fs = wavread("SinustonAB2A32.wav");
    #wavplay(y, fs)
	inline_audioplayer("SinustonAB2A32.wav")
end

function AB2Aufgabe32(S)
    x = 0:0.00001:0.01;
    y = zeros(length(x));
    for i = 1:length(x)
        y[i] = S(x[i]);
    end
    labels = ["S"];
    plot(x,y, label = labels)
end

function AB2Zusatzaufgabe(F)
    saeg(t) = mod2pi(t+pi)-pi;
    x = -2*pi:0.01:2*pi;
    y = zeros(length(x),2);
    for i = 1:length(x)
        y[i,1] = F(x[i]);
        y[i,2] = saeg(x[i]);
    end
    labels = ["F" "Saegezahnfunktion"];
    plot(x,y, label = labels)
end

function AB3Aufgabe1(pointFreq_8,T_delta_8,anchorFreq_3,pointFreq_10,anchorFreq_4,pointFreq_11_2,T_delta_11_2,anchorFreq_6,pointFreq_16_2,T_delta_16_2,pointFreq_17_2,T_delta_17_2,pointFreq_18,T_delta_18,pointFreq_19,T_delta_19,pointFreq_20,T_delta_20)
    if pointFreq_8 != 40
		@warn "Frequenz von Datenpunkt 8 ist falsch!"
    end
    if T_delta_8 != 4
		@warn "Abstand von Datenpunkt 8 zum Anchor Point 2 stimmt nicht!"
    end
    if anchorFreq_3 != 10
		@warn "Frequenz von Anchor Point 3 stimmt nicht!"
    end
    if pointFreq_10 != 50
		@warn "Frequenz von Datenpunkt 10 stimmt nicht!"
    end
    if anchorFreq_4 != 50
		@warn "Frequenz von Anchor Point 4 stimmt nicht!"
    end
    if pointFreq_11_2 != 10
		@warn "Frequenz von Datenpunkt 11 stimmt nicht!"
    end
    if T_delta_11_2 != 1
		@warn "Abstand von Datenpunkt 11 zu Anchor Point 4 stimmt nicht!"
    end
    if anchorFreq_6 != 30
		@warn "Frequenz von Anchor Point 6 stimmt nicht!"
    end
    if pointFreq_16_2 != 30
		@warn "Frequenz von Datenpunkt 16 stimmt nicht!"
    end
    if T_delta_16_2 != 0
		@warn "Abstand von Datenpunkt 16 zum Anchor Point 6 stimmt nicht!"
    end
    if pointFreq_17_2 != 20
		@warn "Frequenz von Datenpunkt 17 stimmt nicht!"
    end
    if T_delta_17_2 != 1
		@warn "Abstand von Datenpunkt 17 zu Anchor Point 6 stimmt nicht!"
    end
    if pointFreq_18 != 40
		@warn "Frequenz von Datenpunkt 18 stimmt nicht!"
    end
    if T_delta_18 != 2
		@warn "Abstand zwischen Datenpunkt 18 und Anchor Point 6 stimmt nicht!"
    end
    if pointFreq_19 != 10
		@warn "Frequenz von Datenpunkt 19 stimmt nicht!"
    end
    if T_delta_19 != 3
		@warn "Abstand zwischen Datenpunkt 19 und Anchor Point 6 stimmt nicht!"
    end
    if pointFreq_20 != 30
		@warn "Frequenz von Datenpunkt 20 stimmt nicht!"
    end
    if T_delta_20 != 4
		@warn "Abstand von Datenpunkt 20 zu Anchor Point 6 stimmt nicht!"
    end
end

function AB4Aufgabe1a(Spalte3,Spalte4)
    Spalte2 = [22; 22; 22; 49; 49; 34; 34; 34; 85; 49; 58; 58; 58; 58; 58; 19; 19; 19; 67; 67; 82; 82; 52; 70; 70; 25; 25; 25; 40; 40];
    if Spalte3 != [1; 1; 1; 1; 1; 4; 4; 4; 4; 4; 7; 7; 7; 7; 7; 10; 10; 10; 10; 10; 13; 13; 13; 13; 13; 16; 16; 16; 16; 16]
		@warn "Spalte3 ist falsch!"
    end
    Spalte4_temp = Spalte2 - Spalte3;
    if Spalte4 !=  Spalte4_temp
		@warn "Spalte4 ist falsch!"
    end
    print("Spalte4 = ", Spalte4)
end

function AB4Aufgabe1b(Spalte2,Spalte3,Spalte4)
    if Spalte2 != [52; 52; 52; 52; 52; 55; 55; 55; 55; 55; 58; 58; 58; 58; 58; 61; 61; 61; 61; 61; 64; 64; 64; 64; 64; 67; 67; 67; 67; 67];
		@warn "Spalte2 ist falsch!"
    end
    if Spalte3 != [1; 1; 1; 1; 1; 4; 4; 4; 4; 4; 7; 7; 7; 7; 7; 10; 10; 10; 10; 10; 13; 13; 13; 13; 13; 16; 16; 16; 16; 16]
		@warn "Spalte3 ist falsch!"
    end
    Spalte4_temp = Spalte2 - Spalte3;
    if Spalte4 !=  Spalte4_temp
		@warn "Spalte4 ist falsch!"
    end
    print("Spalte4 = ", Spalte4)
end

function AB4Aufgabe1c(N)
    if N != 2
		@warn "Falsch!"
    end
    if N == 2
        print("Richtig, Song02 ist der gesuchte Song!")
    end
end

end


text:
using CloseOpenIntervals
using Test

function mysum2(X = CartesianIndices((SafeCloseOpen(10),SafeCloseOpen(10))))
  s = 0
  for I in X
    s += sum(Tuple(I))
  end
  s
end
mysum2_allocated() = @allocated(mysum2())
const ArrayInterface = CloseOpenIntervals.ArrayInterface
@testset "CloseOpenIntervals.jl" begin
  function mysum(x, N)
    s = zero(eltype(x))
    @inbounds @fastmath for i ∈ CloseOpen(N)
      s += x[i+1]
    end
    s
  end
  function mysafesum(x, N)
    s = zero(eltype(x))
    @inbounds @fastmath for i ∈ SafeCloseOpen(N)
      s += x[i+1]
    end
    s
  end
  x = rand(128);
  for n ∈ 1:128
    @test mysum(x,n) ≈ mysafesum(x,n) ≈ sum(view(x,1:n))
    @test length(CloseOpen(n)) == n
  end
  @test @inferred(mysum(x,0)) == first(x)
  @test @inferred(mysafesum(x,0)) == 0.0
  @test @inferred(mysum(x, CloseOpenIntervals.StaticInt(128))) ≈ sum(x)
  @test @inferred(ArrayInterface.static_length(CloseOpen(CloseOpenIntervals.StaticInt(128)))) === CloseOpenIntervals.StaticInt(128)
  @test @inferred(eltype(CloseOpen(7))) === Int
  @test ArrayInterface.known_length(CloseOpen(CloseOpenIntervals.StaticInt(128))) == 128
  @test @inferred(mysum2()) == sum(0:9)*2*length(0:9)
  @test mysum2_allocated() == 0
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.10.3
#   kernelspec:
#     display_name: Julia 1.7.0
#     language: julia
#     name: julia-1.7
# ---

# %%
using Distributions

tstat(μ, X) = (mean(X) - μ)/√(var(X) / length(X))
pvalue(μ, X) = 2ccdf(TDist(length(X) - 1), abs(tstat(μ, X)))

function ci(X, α)
    X̄ = mean(X)
    c = quantile(TDist(length(X) - 1), 1 - α/2)
    d = √(var(X) / length(X))
    X̄ - c*d, X̄ + c*d
end

isininterval(x, int) = first(int) ≤ x ≤ last(int)
isinci(μ, X, α) = isininterval(μ, ci(X, α))

prob_trueisinci(dist, n, α; L=10^6) = mean(isinci(mean(dist), rand(dist, n), α) for _ in 1:L)

# %%
prob_trueisinci(Normal(), 20, 0.05)

# %%
prob_trueisinci(Uniform(), 20, 0.05)

# %%
prob_trueisinci(Exponential(), 20, 0.05)

# %%
prob_trueisinci(MixtureModel([Normal(), Normal(10, 1)], [0.95, 0.05]), 20, 0.05)

# %%
using Plots

# %%
normal = Normal(mean(Uniform()), std(Uniform()))
P = plot(x -> pdf(Uniform(), x), -1, 2; label="Uniform()")
plot!(x -> pdf(normal, x), -1, 2; label="normal approx.", ls=:dash)
title!("pdfs")
Q = plot(x -> cdf(Uniform(), x), -1, 2; label="Uniform()")
plot!(x -> cdf(normal, x), -1, 2; label="normal approx.", ls=:dash)
title!("cdfs"; legend=:bottomright)
plot(P, Q; size=(720, 300))

# %%
normal = Normal(mean(Exponential()), std(Exponential()))
P = plot(x -> pdf(Exponential(), x), -3, 6; label="Exponential()")
plot!(x -> pdf(normal, x), -3, 6; label="normal approx.", ls=:dash)
title!("pdfs")
Q = plot(x -> cdf(Exponential(), x), -3, 6; label="Exponential()")
plot!(x -> cdf(normal, x), -3, 6; label="normal approx.", ls=:dash)
title!("cdfs"; legend=:bottomright)
plot(P, Q; size=(720, 300))

# %%
mixnormal = MixtureModel([Normal(), Normal(10, 1)], [0.95, 0.05])
P = plot(x -> pdf(mixnormal, x), -5, 15; label="mixnormal")
plot!(x -> pdf(Normal(), x), -5, 15; label="Normal()", ls=:dash)
title!("pdfs")
Q = plot(x -> cdf(mixnormal, x), -5, 15; label="mixnormal")
plot!(x -> cdf(Normal(), x), -5, 15; label="Normal()", ls=:dash)
title!("cdfs"; legend=:bottomright)
plot(P, Q; size=(720, 300))

# %%
mixnormal = MixtureModel([Normal(), Normal(10, 1)], [0.95, 0.05])
normal = Normal(mean(mixnormal), std(mixnormal))
P = plot(x -> pdf(mixnormal, x), -5, 15; label="mixnormal")
plot!(x -> pdf(normal, x), -5, 15; label="normal approx.", ls=:dash)
title!("pdfs")
Q = plot(x -> cdf(mixnormal, x), -5, 15; label="mixnormal")
plot!(x -> cdf(normal, x), -5, 15; label="normal approx.", ls=:dash)
title!("cdfs"; legend=:bottomright)
plot(P, Q; size=(720, 300))

# %%
model = Normal()
@show model
μ, σ² = mean(model), var(model)
PP = []
for n in (5, 10, 20, 100)
    Z = [(mean(rand(model, n)) - μ)/√(σ²/n) for _ in 1:10^5]
    P = histogram(Z; norm=true, alpha=0.3, bin = -3.5:0.1:6, label="")
    plot!(x -> pdf(Normal(), x); label="", ls=:dash, lw=1.5)
    title!("sample size n = $n")
    push!(PP, P)
end
plot(PP...; size=(800, 600), layout=(2, 2), titlefontsize=11)

# %%
model = Uniform()
@show model
μ, σ² = mean(model), var(model)
PP = []
for n in (5, 10, 20, 100)
    Z = [(mean(rand(model, n)) - μ)/√(σ²/n) for _ in 1:10^5]
    P = histogram(Z; norm=true, alpha=0.3, bin = -3.5:0.1:6, label="")
    plot!(x -> pdf(Normal(), x); label="", ls=:dash, lw=1.5)
    title!("sample size n = $n")
    push!(PP, P)
end
plot(PP...; size=(800, 600), layout=(2, 2), titlefontsize=11)

# %%
model = Exponential()
@show model
μ, σ² = mean(model), var(model)
PP = []
for n in (5, 10, 20, 100)
    Z = [(mean(rand(model, n)) - μ)/√(σ²/n) for _ in 1:10^5]
    P = histogram(Z; norm=true, alpha=0.3, bin = -3.5:0.1:6, label="")
    plot!(x -> pdf(Normal(), x); label="", ls=:dash, lw=1.5)
    title!("sample size n = $n")
    push!(PP, P)
end
plot(PP...; size=(800, 600), layout=(2, 2), titlefontsize=11)

# %%
model = MixtureModel([Normal(), Normal(10, 1)], [0.95, 0.05])
@show model
μ, σ² = mean(model), var(model)
PP = []
for n in (5, 10, 20, 100)
    Z = [(mean(rand(model, n)) - μ)/√(σ²/n) for _ in 1:10^5]
    P = histogram(Z; norm=true, alpha=0.3, bin = -3.5:0.1:6, label="")
    plot!(x -> pdf(Normal(), x); label="", ls=:dash, lw=1.5)
    title!("sample size n = $n")
    push!(PP, P)
end
plot(PP...; size=(800, 600), layout=(2, 2), titlefontsize=11)

# %%

max_stars_repo_path: 0012/probability of true value being in CI.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
module KMarkovGenSeq

using FASTX, LinearAlgebra
import CUDA: copyto!, CuMatrix, CUBLAS
import CUDA

export count_freq!, db_freq, query_freq

function count_freq!(cnt::AbstractVector, seq::AbstractVector{UInt8}, len)
    mask = typemax(UInt) >> (64 - 2 * len)
    len += 1
    idx = UInt(0)
    flg = 0
    for i in 1:len-1
        @inbounds idx += seq[i]
        (seq[i] == 4) && (flg = len)
        idx <<= 2
        flg -= 1
    end

    for num in @view(seq[len:end])
        (num != 0x4) ? (idx += num) : (flg = len)
        (flg <= 0) && (@inbounds cnt[idx + 1] += 1)
        flg -= 1
        idx &= mask
        idx <<= 2
    end
end

"""
Process the counts of each k-strr in several database
sequences. And the -log frequency will be given for the origin
sequences and their reverse complementary sequences.

The first value to return is the descriptions of every database
sequences. The second is the -log frequencies.

Args:

- `dbseqs`: the forlder of the database sequences
- `len`: do K-Markov algo with K = `len`
"""
function db_freq(dbseqs, len)
    map = fill(0x4, 127)
    map[['A', 'G', 'C', 'T'] .|> UInt8] .= [0x0, 0x1, 0x2, 0x3]
    map[['a', 'g', 'c', 't'] .|> UInt8] .= [0x0, 0x1, 0x2, 0x3]
    res = zeros(Float32, 4 ^ (len + 1), length(dbseqs), 2)

    descriptions = String[] #empty array
    @sync @inbounds for (idx, dbseq) in dbseqs |> enumerate
        reader = open(FASTA.Reader, dbseq) #read one db seq
        record = FASTA.Record()
        read!(reader, record)
        (description(record) !== nothing) && push!(descriptions, description(record)) #store the descriptions
        @simd for i in record.sequence
            record.data[i] = map[record.data[i]]
        end

        Threads.@spawn begin
            count_freq!(@view(res[:, idx, 1]), @view(record.data[record.sequence]), len)
            @simd for i in record.sequence
                record.data[i] == 0x4 || (record.data[i] = 0x3 - record.data[i])
            end
            count_freq!(@view(res[:, idx, 2]), @view(record.data[record.sequence |> reverse]), len)
        end

        close(reader)
    end
    @inbounds Threads.@threads for idx in 1:length(dbseqs)
        for i in 1:4:4^(len + 1) # calculate the P(Oi | O_{i+1}) = F(O)
            res[i:i+3, idx, 1] ./= sum(@view(res[i:i+3, idx, 1])) # if all 0, gives NaN
            res[i:i+3, idx, 2] ./= sum(@view(res[i:i+3, idx, 2]))
        end
    end

    res .= .-log.(res) # if 0, gives -inf; if NaN, gives NaN
    res[isnan.(res)] .= 10 # treat inf and NaN as 10
    res[isinf.(res)] .= 10
    descriptions, res
end

"""
Perform K-Markov algorithm to predict where sequences come from.

Args:

- `query_seqs`: the file containing the query sequences
- `dbseqs`: the forlder of the database sequences
- `len`: do K-Markov algo with K = `len`
- `out_file`: the output file to give the predictions, default is `nothing` and no file will be given
- `use_gpu`: whether to use gpu
"""
function query_freq(query_seqs, dbseq_dir, len, out_file=nothing, use_gpu=false)
    batch_size = 100
    map = fill(0x4, 127)
    map[['A', 'G', 'C', 'T'] .|> UInt8] .= [0x0, 0x1, 0x2, 0x3]
    map[['a', 'g', 'c', 't'] .|> UInt8] .= [0x0, 0x1, 0x2, 0x3]

    query_freqs_h = zeros(Float32, 4 ^ (len + 1), batch_size)
    query_freqs_d = use_gpu ? CUDA.cu(query_freqs_h) : nothing

    dbseqs = dbseq_dir * "/" .* readdir(dbseq_dir)
    descriptions, db_freq_h = db_freq(dbseqs, len)
    db_freq_d = use_gpu ? CUDA.cu(db_freq_h) : nothing

    product_h = zeros(Float32, batch_size, length(dbseqs), 2)
    product_d = use_gpu ? CUDA.cu(product_h) : nothing

    reader = open(FASTA.Reader, query_seqs)

    total_query = 0
    cls_res = Int[]

    io = !isnothing(out_file) ? open(joinpath("output/", out_file), "w") : nothing
    while !eof(reader)
        cnt = 0
        # VERY IMPORTANT INITIALIZE!!!
        fill!(query_freqs_h, 0.0)
        @sync for idx in 1:batch_size
            eof(reader) && break

            record = FASTA.Record()
            read!(reader, record)
            @inbounds @simd for i in record.sequence
                record.data[i] = map[record.data[i]]
            end
            Threads.@spawn count_freq!(@view(query_freqs_h[:, idx]), @view(record.data[record.sequence]), len)
            cnt += 1
        end
        if db_freq_d !== nothing
            copyto!(query_freqs_d, query_freqs_h)
            cu_mat_mul!(@view(product_d[:, :, 1]), query_freqs_d, @view(db_freq_d[:, :, 1]))
            cu_mat_mul!(@view(product_d[:, :, 2]), query_freqs_d, @view(db_freq_d[:, :, 2]))
            copyto!(product_h, product_d)
        else
            mul!(@view(product_h[:, :, 1]), query_freqs_h', @view(db_freq_h[:, :, 1]))
            mul!(@view(product_h[:, :, 2]), query_freqs_h', @view(db_freq_h[:, :, 2]))
        end

        argmin_score = argmin(product_h, dims=(2, 3))

        append!(cls_res, argmin_score[i][2] for i in 1:cnt)
        if !isnothing(io)
            for i in 1:cnt
                println(io, i + total_query, "\t", descriptions[argmin_score[i][2]])
            end
        end
        total_query += cnt
    end
    !isnothing(io) && close(io)

    if use_gpu
        CUDA.unsafe_free!(product_d)
        CUDA.unsafe_free!(query_freqs_d)
        CUDA.unsafe_free!(db_freq_d)
    end

    total_query, cls_res, descriptions
end

"""
Wrapper of cuda lib functon cublasSgemm_v2 to do matrix multiplication.

Calculate C = A' ⋅ B

dims: m×n = (p×m)'p×n
"""

function cu_mat_mul!(C::CuMatrix, A::CuMatrix, B::CuMatrix)
    handle = CUBLAS.handle()
    row_A, col_A = size(A)
    row_B, col_B = size(B)
    CUBLAS.cublasSgemm_v2(
        handle,
        CUBLAS.CUBLAS_OP_T,
        CUBLAS.CUBLAS_OP_N,
        col_A, col_B, row_B,
        1.0,
        A, row_A,
        B, row_B,
        0.0,
        C, col_A
    )
end

end


text:
using Atum
using Atum.Euler

using Test
using Printf
using StaticArrays: SVector
using LinearAlgebra: norm

if !@isdefined integration_testing
  const integration_testing = parse(
    Bool,
    lowercase(get(ENV, "ATUM_INTEGRATION_TESTING", "false")),
  )
end

function vortex(law, x⃗, t)
  FT = eltype(law)
  γ = constants(law).γ
  c⃗ = SVector(FT(5), FT(0))
  u⃗₀ = SVector(FT(1), FT(0))
  β = FT(5)
  r⃗ = x⃗ - c⃗ - u⃗₀ * t
  f = β * exp(1 - r⃗' * r⃗)

  ρ = (1 - (γ - 1) * f ^ 2 / (16 * γ * π ^ 2)) ^ (1 / (γ - 1))
  u⃗ = u⃗₀ + f / 2π * SVector(-r⃗[2], r⃗[1])
  ρu⃗ = ρ * u⃗
  p = ρ ^ γ
  ρe = Euler.energy(law, ρ, ρu⃗, p)

  SVector(ρ, ρu⃗..., ρe)
end

function run(A, FT, N, K; volume_form=WeakForm())
  Nq = N + 1

  law = EulerLaw{FT, 2}()

  cell = LobattoCell{FT, A}(Nq, Nq)
  vx = range(FT(0), stop=FT(20), length=K+1)
  vy = range(FT(-5), stop=FT(5), length=K+1)
  grid = brickgrid(cell, (vx, vy); periodic=(true, true))

  dg = DGSEM(; law, grid, volume_form,
               surface_numericalflux = RusanovFlux())

  cfl = FT(1 // 4)
  dt = cfl * step(vy) / N
  timeend = FT(5)

  q = vortex.(Ref(law), points(grid), FT(0))

  @info @sprintf """Starting
  N           = %d
  K           = %d
  volume_form = %s
  norm(q)     = %.16e
  """ N K volume_form weightednorm(dg, q)

  odesolver = LSRK54(dg, q, dt)
  solve!(q, timeend, odesolver)

  qexact = vortex.(Ref(law), points(grid), timeend)
  errf = weightednorm(dg, q .- qexact)

  @info @sprintf """Finished
  norm(q)      = %.16e
  norm(q - qe) = %.16e
  """ weightednorm(dg, q) errf
  errf
end

let
  A = Array
  FT = Float64
  N = 4

  expected_error = Dict()

  #form, lev
  expected_error[WeakForm(), 1] = 2.1296858506446990e+00
  expected_error[WeakForm(), 2] = 5.9526168482660669e-01
  expected_error[WeakForm(), 3] = 6.4358678661161267e-02
  expected_error[WeakForm(), 4] = 2.3640146930597296e-03

  expected_error[FluxDifferencingForm(EntropyConservativeFlux()), 1] = 2.0789522602391144e+00
  expected_error[FluxDifferencingForm(EntropyConservativeFlux()), 2] = 3.5478518460567871e-01
  expected_error[FluxDifferencingForm(EntropyConservativeFlux()), 3] = 3.0653814065252447e-02
  expected_error[FluxDifferencingForm(EntropyConservativeFlux()), 4] = 1.3056958215461041e-03

  nlevels = integration_testing ? 4 : 1

  @testset for volume_form in (WeakForm(),
                               FluxDifferencingForm(EntropyConservativeFlux()))
    errors = zeros(FT, nlevels)
    for l in 1:nlevels
      K = 5 * 2 ^ (l - 1)
      errf = run(A, FT, N, K; volume_form)
      errors[l] = errf
      @test errors[l] ≈ expected_error[volume_form, l]
    end

    if nlevels > 1
      rates = log2.(errors[1:(nlevels-1)] ./ errors[2:nlevels])
      @info "Convergence rates\n" *
        join(["rate for levels $l → $(l + 1) = $(rates[l])" for l in 1:(nlevels - 1)], "\n")
      @test rates[end] ≈ N + 1 atol = 0.5
    end
  end
end


text:
export LinearRegressor, RidgeRegressor, LassoRegressor, ElasticNetRegressor,
       RobustRegressor, HuberRegressor, QuantileRegressor, LADRegressor,
       LogisticClassifier, MultinomialClassifier

const SymStr = Union{Symbol,String}

include("regressors.jl")
include("classifiers.jl")

const REG_MODELS = (LinearRegressor, RidgeRegressor, LassoRegressor,
                    ElasticNetRegressor, RobustRegressor, HuberRegressor,
                    QuantileRegressor, LADRegressor)
const CLF_MODELS = (LogisticClassifier, MultinomialClassifier)
const ALL_MODELS = (REG_MODELS..., CLF_MODELS...)

#= ==========
   REGRESSORS
   ========== =#

function MMI.fit(m::Union{REG_MODELS...}, verb::Int, X, y)
    Xmatrix = MMI.matrix(X)
    sch = MMI.schema(X)
    features = (sch === nothing) ? nothing : sch.names
    reg     = glr(m)
    solver  = m.solver === nothing ? _solver(reg, size(Xmatrix)) : m.solver
    # get the parameters
    θ = fit(reg, Xmatrix, y; solver=solver)
    # return
    return (θ, features), nothing, NamedTuple{}()
end

MMI.predict(m::Union{REG_MODELS...}, (θ, features), Xnew) = apply_X(MMI.matrix(Xnew), θ)

function MMI.fitted_params(m::Union{REG_MODELS...}, (θ, features))
    m.fit_intercept && return (coefs = coef_vec(θ[1:end-1], features), intercept = θ[end])
    return (coefs = coef_vec(θ, features), intercept = nothing)
end

#= ===========
   CLASSIFIERS
   =========== =#

function MMI.fit(m::Union{CLF_MODELS...}, verb::Int, X, y)
    Xmatrix  = MMI.matrix(X)
    sch = MMI.schema(X)
    features = (sch === nothing) ? nothing : sch.names
    yplain   = convert.(Int, MMI.int(y))
    classes  = MMI.classes(y[1])
    nclasses = length(classes)
    if nclasses < 2
        throw(DomainError("The target `y` needs to have two or more levels."))
    elseif nclasses == 2
        # recode to ± 1
        yplain[yplain .== 1] .= -1
        yplain[yplain .== 2] .= 1
        # force the binary case
        nclasses    = 0
    end
    # NOTE: here the number of classes is either 0 or > 2
    clf = glr(m, nclasses)
    solver = m.solver === nothing ? _solver(clf, size(Xmatrix)) : m.solver
    # get the parameters
    θ = fit(clf, Xmatrix, yplain, solver=solver)
    # return
    return (θ, features, classes, nclasses), nothing, NamedTuple{}()
end

function MMI.predict(m::Union{CLF_MODELS...}, (θ, features, classes, c), Xnew)
    Xmatrix = MMI.matrix(Xnew)
    preds = apply_X(Xmatrix, θ, c)
    if c > 2 # multiclass
        preds .= softmax(preds)
    else # binary (necessarily c==0)
        preds  .= sigmoid.(preds)
        preds   = hcat(1.0 .- preds, preds) # scores for -1 and 1
        return MMI.UnivariateFinite(classes, preds)
    end
    return MMI.UnivariateFinite(classes, preds)
end

function MMI.fitted_params(m::Union{CLF_MODELS...}, (θ, features, classes, c))
    # helper function to assemble the results
    _fitted_params(coefs, features, intercept) =
        (classes=classes, coefs=coef_vec(coefs, features), intercept=intercept)
    # multiclass or binary?
    if c > 2
        W = reshape(θ, :, c)
        if m.fit_intercept
            return _fitted_params(W, features, W[end, :])
        end
        return _fitted_params(W, features, nothing)
    end
    # single class (necessarily c==0)
    m.fit_intercept && return _fitted_params(θ[1:end-1], features, θ[end])
    return _fitted_params(θ, features, nothing)
end

@static VERSION < v"1.1" && (eachrow(A::AbstractVecOrMat) = (view(A, i, :) for i in axes(A, 1)))

coef_vec(W::AbstractMatrix, features) =
    [feature => coef for (feature, coef) in zip(features, eachrow(W))]
coef_vec(θ::AbstractVector, features) =
    [feature => coef for (feature, coef) in zip(features, θ)]
coef_vec(W::AbstractMatrix, ::Nothing) = W
coef_vec(θ::AbstractVector, ::Nothing) = θ

#= =======================
   METADATA FOR ALL MODELS
   ======================= =#

MMI.metadata_pkg.(ALL_MODELS,
    name="MLJLinearModels",
    uuid="6ee0df7b-362f-4a72-a706-9e79364fb692",
    url="https://github.com/alan-turing-institute/MLJLinearModels.jl",
    julia=true,
    license="MIT",
    is_wrapper=false)

descr_(M) = descr(M) *
    "\n→ based on [MLJLinearModels](https://github.com/alan-turing-institute/MLJLinearModels.jl)" *
    "\n→ do `@load $(MMI.name(M)) pkg=\"MLJLinearModels\" to use the model.`" *
    "\n→ do `?$(MMI.name(M))` for documentation."
lp_(M) = "MLJLinearModels.$(MMI.name(M))"

for M in REG_MODELS
    MMI.metadata_model(M,
        input=MMI.Table(MMI.Continuous),
        target=AbstractVector{MMI.Continuous},
        weights=false,
        descr=descr_(M), path=lp_(M))
end
for M in CLF_MODELS
    MMI.metadata_model(M,
        input=MMI.Table(MMI.Continuous),
        target=AbstractVector{<:MMI.Finite},
        weights=false,
        descr=descr_(M), path=lp_(M))
end


text:
##

using ACE
using Printf, Test, LinearAlgebra, StaticArrays
using ACE: evaluate, evaluate_d, Rn1pBasis, Ylm1pBasis,
           Product1pBasis, Scal1pBasis
using Random: shuffle
using ACEbase.Testing: fdtest, print_tf

##

maxdeg = 5
r0 = 1.0
rcut = 3.0
maxorder = 3
Bsel = SimpleSparseBasis(maxorder, maxdeg)

trans = trans = PolyTransform(1, r0)
Pk = ACE.scal1pbasis(:x, :k, maxdeg, trans, rcut)
RnYlm = ACE.Utils.RnYlm_1pbasis()

B1p = RnYlm * Pk
ACE.init1pspec!(B1p, Bsel)
length(B1p)

##

PosScalState{T} = ACE.State{(:rr, :x), Tuple{SVector{3, T}, T}}

Base.promote_rule(::Union{Type{S}, Type{PosScalState{S}}},
             ::Type{PosScalState{T}}) where {S, T} =
      PosScalState{promote_type(S, T)}

X = rand(PosScalState{Float64})
cfg = ACEConfig([ rand(PosScalState{Float64}) for _=1:10 ])

Rn = B1p.bases[1]
Ylm = B1p.bases[2]
Pk = B1p.bases[3]

ACE.gradtype(B1p, cfg)
ACE.valtype(B1p, cfg)

ACE.acquire_B!(Pk, cfg)

A = evaluate(B1p, cfg)
dA = evaluate_d(B1p, cfg)
A1, dA1 = ACE.evaluate_ed(B1p, cfg)

println(@test( A ≈ A1 ))
println(@test( dA ≈ dA1 ))

##

# gradient test

_vec2X(x) = PosScalState{eltype(x)}((rr = SVector{3}(x[1:3]), x = x[4]))
_vec2cfg(x) = ACEConfig( [_vec2X(x)] )
_X2vec(X) = [X.rr; [X.x]]

for ntest = 1:30
   x0 = randn(4)
   c = rand(length(B1p))
   F = x -> sum(ACE.evaluate(B1p, _vec2X(x)) .* c)
   dF = x -> _X2vec( sum(ACE.evaluate_d(B1p, _vec2cfg(x)) .* c))
   print_tf(@test fdtest(F, dF, x0; verbose=false))
end
println()

##


text:
# Autogenerated wrapper script for pprof_jll for armv7l-linux-gnueabihf
export pprof

JLLWrappers.@generate_wrapper_header("pprof")
JLLWrappers.@declare_executable_product(pprof)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        pprof,
        "bin/pprof",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
rand(::Type{RuleNode}, grammar::Grammar, typ::Symbol, max_depth::Int=10)
Generates a random RuleNode of return type typ and maximum depth max_depth.
"""
function Base.rand(::Type{RuleNode}, grammar::ContextGrammar, typ::Symbol, max_depth::Int=10,
	context::Union{GrammarContext,Nothing}=nothing, bin::Union{NodeRecycler,Nothing}=nothing)

	if context === nothing
		# first call, context not set yet
		init_node = RuleNode(0)
		context = GrammarContext(init_node)
	end

	rules = deepcopy(grammar[typ])
	# propagate constraints
	rules = propagate_contraints(grammar, context, rules)

	if max_depth <= 1
	    terminals = filter(r->isterminal(grammar,r), rules)
	    rule_index = !isempty(terminals) ? StatsBase.sample(terminals) : StatsBase.sample(rules)
	else
	    rule_index = StatsBase.sample(rules)
	end

	rulenode = iseval(grammar, rule_index) ?
	    RuleNode(bin, rule_index, Core.eval(grammar, rule_index)) :
	    RuleNode(bin, rule_index)

	if !grammar.isterminal[rule_index]
	    	for ch in child_types(grammar, rule_index)
			#create context for the child expansion
			new_context = GrammarContext(rulenode, deepcopy(context.nodeLocation))
			push!(new_context.nodeLocation, length(rulenode.children) + 1)
			push!(rulenode.children, rand(RuleNode, grammar, ch, max_depth-1, new_context, bin))
	    	end
	end

	return rulenode
end

"""
    rand(::Type{RuleNode}, grammar::Grammar, typ::Symbol, dmap::AbstractVector{Int}, max_depth::Int=10)
Generates a random RuleNode of return type typ and maximum depth max_depth guided by a minimum depth map dmap.
"""
function Base.rand(::Type{RuleNode}, grammar::ContextGrammar, typ::Symbol, dmap::AbstractVector{Int},
	max_depth::Int=10, context::Union{GrammarContext,Nothing}=nothing, bin::Union{NodeRecycler,Nothing}=nothing)

	if context === nothing
		# first call, context not set yet
		init_node = RuleNode(0)
		context = GrammarContext(init_node)
	end

    	rules = grammar[typ]
	    # propagate constraints
	rules = propagate_constraints(grammar, context, rules)
    	rule_index = StatsBase.sample(filter(r->dmap[r] ≤ max_depth, rules))

	rulenode = iseval(grammar, rule_index) ?
 		RuleNode(bin, rule_index, Core.eval(grammar, rule_index)) :
 		RuleNode(bin, rule_index)

    	if !grammar.isterminal[rule_index]
        	for ch in child_types(grammar, rule_index)
			# create context for the child expansion
			new_context = GrammarContext(rulenode, deepcopy(context.nodeLocation))
			push!(new_context.nodeLocation, length(rulenode.children) + 1)
            		push!(rulenode.children, rand(RuleNode, grammar, ch, dmap, max_depth-1, new_context, bin))
       	 	end
    	end
    return rulenode
end


text:
module LatinSquares
using JuMP
using MathProgBase
using Cbc
using ChooseOptimizer

SOLVER = CbcSolver

export set_latin_solver

include("ortho_latin.jl")
include("latin.jl")
include("latin_print.jl")

# package code goes here

end # module


text:
abstract type Importance <: InferenceType end
struct LikelihoodWeighting <: Importance end
struct ImportanceSampling <: Importance end
const LW = LikelihoodWeighting()
const IS = ImportanceSampling()

@doc raw"""
    likelihood_weighting_results()
Outer constructor for `SamplingResults`.
"""
likelihood_weighting_results() = NonparametricSamplingResults{LikelihoodWeighting}(LW, Array{Float64, 1}(), Array{Any, 1}(), Array{Trace, 1}())
likelihood_weighting_results(A::DataType, D::DataType) = NonparametricSamplingResults{LikelihoodWeighting}(LW, Array{Float64, 1}(), Array{Any, 1}(), Array{TypedTrace{A, D}, 1}())

@doc raw"""
    function lw_step(f, params...)

Perform one step of likelihood weighting -- draw a single proposal from the prior and compute
the log weight as equal to the likelihood. Returns a tuple (log weight, rval, trace).
"""
function lw_step(f, params...)
    t = trace()
    r_n = f(t, params...)
    logprob!(t)
    log_w = loglikelihood(t)
    (log_w, r_n, t)
end

@doc raw"""
    function lw_step(f, types::Tuple{DataType, DataType}, params...)

Perform one step of likelihood weighting -- draw a single proposal from the prior and compute
the log weight as equal to the likelihood. Returns a tuple (log weight, rval, trace).
"""
function lw_step(f, types::Tuple{DataType, DataType}, params...)
    t = trace(types[1], types[2])
    r_n = f(t, params...)
    logprob!(t)
    log_w = loglikelihood(t)
    (log_w, r_n, t)
end

@doc raw"""
    function likelihood_weighting(f, params...; nsamples :: Int = 1)

Given a stochastic function f and arguments to the function `params...`, executes `nsamples`
iterations of importance sampling by using the prior as a proposal distribution. The importance
weights are given by ``\log W_n = \ell(t_n)``. Returns an `SamplingResults` instance.
"""
function likelihood_weighting(f, params...; nsamples :: Int = 1)
    results = likelihood_weighting_results()
    for n in 1:nsamples
        log_w, r_n, t = lw_step(f, params...)
        push!(results.log_weights, log_w)
        push!(results.return_values, r_n)
        push!(results.traces, deepcopy(t))
    end
    results
end

@doc raw"""
    function likelihood_weighting(f, types::Tuple{DataType, DataType}, params...; nsamples :: Int = 1)

Given a stochastic function f and arguments to the function `params...`, executes `nsamples`
iterations of importance sampling by using the prior as a proposal distribution. The importance
weights are given by ``\log W_n = \ell(t_n)``. Returns an `SamplingResults` instance.
"""
function likelihood_weighting(f, types::Tuple{DataType, DataType}, params...; nsamples :: Int = 1)
    results = likelihood_weighting_results(types...)
    for n in 1:nsamples
        log_w, r_n, t = lw_step(f, types, params...)
        push!(results.log_weights, log_w)
        push!(results.return_values, r_n)
        push!(results.traces, deepcopy(t))
    end
    results
end

importance_sampling_results() = NonparametricSamplingResults{ImportanceSampling}(IS, Array{Float64, 1}(), Array{Any, 1}(), Array{Trace, 1}())
importance_sampling_results(A::DataType, D::DataType) = NonparametricSamplingResults{ImportanceSampling}(IS, Array{Float64, 1}(), Array{Any, 1}(), Array{TypedTrace{A, D}, 1}())

@doc raw"""
    function is_step(f, q; params = ())

Perform one step of importance sampling -- draw a single sample from the proposal `q`, replay
it through `f`, and record the log weight as ``\log W_n = \log p(x, z_n) - \log q(z_n)``. Returns
a tuple (log weight, rval, trace).
"""
function is_step(f, q; params = ())
    proposed_t = trace()
    q(proposed_t, params...)
    logprob!(proposed_t)
    replayed_t, g = replay(f, proposed_t)
    r_n = g(params...)
    logprob!(replayed_t)
    log_w = replayed_t.logprob_sum - proposed_t.logprob_sum
    (log_w, r_n, deepcopy(replayed_t))
end

@doc raw"""
    function is_step(f, q, types::Tuple{DataType,DataType}; params = ())

Perform one step of importance sampling -- draw a single sample from the proposal `q`, replay
it through `f`, and record the log weight as ``\log W_n = \log p(x, z_n) - \log q(z_n)``. Returns
a tuple (log weight, rval, trace).
"""
function is_step(f, q, types::Tuple{DataType,DataType}; params = ())
    proposed_t = trace(types[1], types[2])
    q(proposed_t, params...)
    logprob!(proposed_t)
    replayed_t, g = replay(f, proposed_t)
    r_n = g(params...)
    logprob!(replayed_t)
    log_w = replayed_t.logprob_sum - proposed_t.logprob_sum
    (log_w, r_n, deepcopy(replayed_t))
end

@doc raw"""
    function importance_sampling(f, q; params = (), nsamples :: Int = 1)

Given a stochastic function `f`, a proposal function `q`, and a tuple of `params` to pass to `f` and `q`,
compute `nsamples` iterations of importance sampling. `q` must have the same input signature
as `f`. Returns a `SamplingResults` instance.
"""
function importance_sampling(f, q; params = (), nsamples :: Int = 1)
    results = importance_sampling_results()
    for n in 1:nsamples
        log_w, r_n, t = is_step(f, q; params = params)
        push!(results.log_weights, log_w)
        push!(results.return_values, r_n)
        push!(results.traces, t)
    end
    results
end

@doc raw"""
    function importance_sampling(f, q, types::Tuple{DataType,DataType}; params = (), nsamples :: Int = 1)

Given a stochastic function `f`, a proposal function `q`, and a tuple of `params` to pass to `f` and `q`,
compute `nsamples` iterations of importance sampling. `q` must have the same input signature
as `f`. Returns a `SamplingResults` instance.
"""
function importance_sampling(f, q, types::Tuple{DataType,DataType}; params = (), nsamples :: Int = 1)
    results = importance_sampling_results(types...)
    for n in 1:nsamples
        log_w, r_n, t = is_step(f, q, types; params = params)
        push!(results.log_weights, log_w)
        push!(results.return_values, r_n)
        push!(results.traces, t)
    end
    results
end

@doc raw"""
    function log_evidence(r :: SamplingResults{LikelihoodWeighting})

Computes the log evidence (log partition function),

``
\log Z \equiv \log p(x) \approx -\log N_{\text{samples}} + \log \sum_{n=1}^{N_{\text{samples}}} W_n.
``
"""
function log_evidence(r :: SamplingResults{I}) where I <: Importance
    # Z \approx 1/L \sum W_l
    # \log Z \approx -\log L + \logsumexp(W)
    -log(length(r.log_weights)) + logsumexp(r.log_weights)
end

@doc raw"""
    function normalized_weights(r :: SamplingResults{LikelihoodWeighting})

Computes the normalized weights ``w_n`` from unnormalized weights ``W_n``:

``
w_n = W_n / p(x) = \exp\{\ell(t_n) - \log Z\}.
``
"""
function normalized_weights(r :: SamplingResults{I}) where I <: Importance
    log_normalizer = log_evidence(r)
    exp.(r.log_weights .- log_normalizer)
end

@doc raw"""
    function sample(r :: SamplingResults{LikelihoodWeighting}, k, n :: Int)

Draws `n` samples from the empirical marginal posterior at address `k`.
"""
function sample(r :: SamplingResults{I}, k, n :: Int) where I <: Importance
    v = r[k]
    weights = StatsBase.Weights(normalized_weights(r))
    if n == 1
        StatsBase.sample(v, weights)
    else
        StatsBase.sample(v, weights, n)
    end
end

function sample(r :: ParametricSamplingResults{I}, k, n :: Int) where I <: Importance
    v = getsampled(r, k)
    weights = StatsBase.Weights(normalized_weights(r))
    if n == 1
        StatsBase.sample(v, weights)
    else
        StatsBase.sample(v, weights, n)
    end
end

export LikelihoodWeighting, ImportanceSampling
export likelihood_weighting, importance_sampling, lw_step, is_step
export log_evidence, normalized_weights, sample


text:
using Unitful
using UnitfulBrewing
using Test

@testset "Dimensions and quantities" begin
    # new dimensions
    @test UnitfulBrewing.𝐂𝐋*UnitfulBrewing.𝐂𝐋 === UnitfulBrewing.𝐂𝐋^2 # Color
    @test UnitfulBrewing.𝐃𝐏*UnitfulBrewing.𝐃𝐏 === UnitfulBrewing.𝐃𝐏^2 # Diastatic Power
    @test UnitfulBrewing.𝐁𝐓*UnitfulBrewing.𝐁𝐓 === UnitfulBrewing.𝐁𝐓^2 # Bitterness
    @test UnitfulBrewing.𝐒𝐂*UnitfulBrewing.𝐒𝐂 === UnitfulBrewing.𝐒𝐂^2 # SugarContents

    # US Volumes not in Unitful
    @test @macroexpand(u"tsp") == u"tsp"
    @test @macroexpand(u"bbl") == u"bbl"
    @test uconvert(u"bbl", 42u"gal") == 1u"bbl"
    @test uconvert(u"gal", 16u"cup") == 1u"gal"
    @test uconvert(u"gal", 768u"tsp") == 1u"gal"

    # Imperial Volumes
    @test @macroexpand(u"ibbl") == u"ibbl"
    @test uconvert(u"ibbl", 36u"igal") == 1u"ibbl"
    @test uconvert(u"igal", 160u"ifloz") == 1u"igal"
    @test uconvert(u"ipt", 4u"gi") == 1u"ipt"

    # sugar content and gravity
    @test @macroexpand(u"°P") == u"°P"
    @test @macroexpand(u"sg") == u"sg"
    @test @macroexpand(u"gu") == u"gu"

    # diastatic power
    @test @macroexpand(u"°Lintner") == u"°Lintner"
    @test @macroexpand(u"°WK") == u"°WK"
    @test uconvert(u"°Lintner", 19u"°WK") == 10u"°Lintner"

    # color
    @test @macroexpand(u"SRM") == u"SRM"
    @test @macroexpand(u"°L") == u"°L"
    @test @macroexpand(u"EBC") == u"EBC"
    @test @macroexpand(u"Lovi") == u"°L"
    @test uconvert(u"srm", 20u"ebc") == (2000//197)u"SRM"
    @test uconvert(u"EBC", 10u"SRM") == (197//10)u"EBC"
    @test uconvert(u"°L", 0u"SRM") == (7600//13546)u"°L"
    @test uconvert(u"°L", 20u"SRM") == (207600//13546)u"°L"

    # pH
    @test [1,2,3]u"pwrH" == u"pH⁺" * [1,2,3]
    @test 3u"pH⁺" < 5u"pwrH"

    # time
    @test uconvert(u"week", 7u"day") == 1u"wk"
    @test uconvert(u"min", 60u"sec") == 1u"min"

    # Throw errors
    @test_throws LoadError @macroexpand(u"ton Lovi")
    @test_throws LoadError @macroexpand(u"lovibond")
end

@testset "Equivalences" begin
    # density and concentration
    @test uconvert(u"mg/l", 1u"ppm", DensityConcentration()) === 1u"mg/l"
    @test uconvert(u"kg/l", 10u"percent", DensityConcentration()) === (1//10)u"kg/l"
    @test uconvert(u"ppm", 1u"mg/l", DensityConcentration()) === 1.0u"ppm"

    # sugar contents and gravity with SugarGravity
    @test uconvert(u"sg", 10u"°P", SugarGravity()) ≈ 1.0401606425u"sg"
    @test uconvert(u"gu", 10u"°P", SugarGravity()) ≈ 40.160642u"gu"
    @test uconvert(u"°P", 1.040u"sg", SugarGravity()) ≈ 9.9615384615u"°P"
    @test uconvert(u"°P", 40u"gu", SugarGravity()) ≈ 9.9615384615u"°P"

    # sugar contents and gravity with SugarGravityQuad
    @test uconvert(u"sg", 10u"°P", SugarGravityQuad()) ≈ 1.040032121u"sg"
    @test uconvert(u"gu", 10u"°P", SugarGravityQuad()) ≈ 40.03u"gu" (atol = 0.01u"gu")
    @test uconvert(u"°P", 1.040u"sg", SugarGravityQuad()) ≈ 9.99224u"°P"
    @test uconvert(u"°P", 40u"gu", SugarGravityQuad()) ≈ 9.99u"°P" (atol = 0.01u"°P")
end


text:
module ML

function ptr(T)
  @assert isbitstype(T)
  Vector{T}([1])
end

const Cstr = Ptr{Cchar}

include("consts.jl")

const Env = Ptr{Nothing}
const Link = Ptr{Nothing}

function find_lib_ker()
    @static if Sys.isapple()
        # TODO: query OS X metadata for non-default installations
        # https://github.com/JuliaLang/julia/issues/8733#issuecomment-167981954
        mpath = "/Applications/Mathematica.app"
        if isdir(mpath)
            lib = joinpath(mpath,"Contents/Frameworks/mathlink.framework/mathlink")
            ker = joinpath(mpath,"Contents/MacOS/MathKernel")
            return lib, ker
        end
    end

    @static if Sys.islinux()
        archdir = Sys.ARCH == :arm ?    "Linux-ARM" :
                  Sys.ARCH == :x86_64 ? "Linux-x86-64" :
                                        "Linux"

        # alternatively, "math" or "wolfram" is often in PATH, so could use
        # echo \$InstallationDirectory | math | sed -n -e 's/Out\[1\]= //p'

        for mpath in ["/usr/local/Wolfram/Mathematica","/opt/Wolfram/WolframEngine"]
            if isdir(mpath)
                vers = readdir(mpath)
                ver = vers[indmax(map(VersionNumber,vers))]

                lib = Libdl.find_library(
                          ["libML$(Sys.WORD_SIZE)i4","libML$(Sys.WORD_SIZE)i3"],
                          [joinpath(mpath,ver,"SystemFiles/Links/MathLink/DeveloperKit",archdir,"CompilerAdditions")])
                ker = joinpath(mpath,ver,"Executables/MathKernel")
                return lib, ker
            end
        end
    end

    @static if Sys.iswindows()
        archdir = Sys.ARCH == :x86_64 ? "Windows-x86-64" :
                                        "Windows"

        #TODO: query Windows Registry, see RCall.jl
        mpath = "C:\\Program Files\\Wolfram Research\\Mathematica"
        if isdir(mpath)
            vers = readdir(mpath)
            ver = vers[indmax(map(VersionNumber,vers))]
            lib = Libdl.find_library(
                          ["libML$(Sys.WORD_SIZE)i4","libML$(Sys.WORD_SIZE)i3"],
                          [joinpath(mpath,ver,"SystemFiles\\Links\\MathLink\\DeveloperKit",archdir,"SystemAdditions")])
            ker = joinpath(mpath,ver,"math.exe")
            return lib, ker
        end
    end

    error("Could not find Mathematica installation")
end

const mlib,mker = find_lib_ker()

function Open(path = mker)
  # MLInitialize
  mlenv = ccall((:MLInitialize, mlib), Env, (Cstr,), C_NULL)
  mlenv == C_NULL && error("Could not MLInitialize")

  # MLOpenString
  # local link
  err = ptr(Cint)
  args = "-linkname '\"$path\" -mathlink' -linkmode launch"
  link = ccall((:MLOpenString, mlib), Link,
                (Env, Cstr, Ptr{Cint}),
                mlenv, args, err)
  err[1]==0 || mlerror(link, "MLOpenString")

  # Ignore first input packet
  @assert NextPacket(link) == Pkt.INPUTNAME
  NewPacket(link)

  return link
end

Close(link::Link) = ccall((:MLClose, mlib), Nothing, (Link,), link)

ErrorMessage(link::Link) =
  ccall((:MLErrorMessage, mlib), Cstr, (Link,), link) |> unsafe_string

for f in [:Error :ClearError :EndPacket :NextPacket :NewPacket]
  fstr = string("ML", f)
  @eval $f(link::Link) = ccall(($fstr, mlib), Cint, (Link,), link)
end

mlerror(link, name) = error("MathLink Error $(Error(link)) in $name: " * ErrorMessage(link))

# Put fns

PutFunction(link::Link, name::AbstractString, nargs::Int) =
  ccall((:MLPutFunction, mlib), Cint, (Link, Cstr, Cint),
    link, name, nargs) != 0 || mlerror(link, "MLPutFunction")

for (f, Tj, Tc) in [(:PutInteger64, Int64, Int64)
                    (:PutInteger32, Int32, Int32)
                    (:PutString, AbstractString, Cstr)
                    (:PutSymbol, Symbol, Cstr)
                    (:PutReal32, Float32, Float32)
                    (:PutReal64, Float64, Float64)]
  fstr = string("ML", f)
  @eval $f(link::Link, x::$Tj) =
          ccall(($fstr, mlib), Cint, (Link, $Tc), link, x) != 0 ||
            mlerror(link, $fstr)
end

# Get fns

GetType(link::Link) =
  ccall((:MLGetType, mlib), Cint, (Link,), link) |> Char

for (f, T) in [(:GetInteger64, Int64)
               (:GetInteger32, Int32)
               (:GetReal32, Float32)
               (:GetReal64, Float64)]
  fstr = string("ML", f)
  @eval function $f(link::Link)
    i = ptr($T)
    ccall(($fstr, mlib), Cint, (Link, Ptr{$T}), link, i) != 0 ||
      mlerror(link, $fstr)
    i[1]
  end
end

function GetString(link::Link)
  s = ptr(Cstr)
  ccall((:MLGetString, mlib), Cint, (Link, Ptr{Cstr}), link, s) != 0 ||
    mlerror(link, "GetString")
  r = s[1] |> unsafe_string
  ReleaseString(link, s)
  return r
end

function GetSymbol(link::Link)
  s = ptr(Cstr)
  ccall((:MLGetSymbol, mlib), Cint, (Link, Ptr{Cstr}), link, s) != 0 ||
    mlerror(link, "GetSymbol")
  r = s[1] |> unsafe_string |> unescape_string |> Symbol
  ReleaseSymbol(link, s)
  return r
end

function GetFunction(link::Link)
  name = ptr(Cstr)
  nargs = ptr(Cint)
  ccall((:MLGetFunction, mlib), Cint, (Link, Ptr{Cstr}, Ptr{Cint}),
    link, name, nargs) != 0 || mlerror(link, "MLGetFunction")
  r = name[1] |> unsafe_string |> Symbol, nargs[1]
  ReleaseString(link, name)
  return r
end

ReleaseString(link::Link, s) = ccall((:MLReleaseString, mlib), Nothing, (Link, Cstr), link, s[1])
ReleaseSymbol(link::Link, s) = ccall((:MLReleaseSymbol, mlib), Nothing, (Link, Cstr), link, s[1])

end


text:
# Chains of InvOrAdjOrVecc
# =================================================

# `*` concatinates chains
# ------------------------------------

function *(O1::A, O2::InvOrAdjOrVecc) where A<:InvOrAdjOrVecc
    tuple(O1, O2)
end
# the second argument type below is unionall
# ... this allows you to bipass for a specific InvOrAdjOrVecc

function *(O1::NTuple{N,InvOrAdjOrVecc}, O2::InvOrAdjOrVecc) where N
    Base.front(O1) * (Base.last(O1) * O2)
end

function *(O1::InvOrAdjOrVecc, O2::NTuple{N,InvOrAdjOrVecc}) where N
    (O1 * Base.first(O2)) * Base.tail(O2)
end

function *(O1::NTuple{N,InvOrAdjOrVecc}, O2::NTuple{M,InvOrAdjOrVecc}) where {N,M}
    tuple(O1..., O2...)
end

# `inv` broadcasts and reverses order
# ------------------------------------

function inv(O1::NTuple{N,InvOrAdjOrVecc}) where N
    tuple((inv(op) for op in reverse(O1))...)
end

pinv(O1::NTuple{N,InvOrAdjOrVecc}) where N = inv(O1)

# `adjoint` broadcasts and reverses order
function adjoint(O1::NTuple{N,InvOrAdjOrVecc}) where N
    tuple((adjoint(op) for op in reverse(O1))...)
end

# `\`
# ------------------------------------

function \(O1::A, O2::InvOrAdjOrVecc) where A<:InvOrAdjOrVecc
    inv(O1) * O2
end
# the second argument type below is unionall
# ... this allows you to bipass for a specific InvOrAdjOrVecc

function \(O1::NTuple{N,InvOrAdjOrVecc}, O2::InvOrAdjOrVecc) where N
    inv(O1) * O2
end

function \(O1::InvOrAdjOrVecc, O2::NTuple{N,InvOrAdjOrVecc}) where N
    inv(O1) * O2
end

function \(O1::NTuple{N,InvOrAdjOrVecc}, O2::NTuple{M,InvOrAdjOrVecc}) where {N,M}
    inv(O1) * O2
end

# `/`
# ------------------------------------

function /(O1::A, O2::InvOrAdjOrVecc) where A<:InvOrAdjOrVecc
    O1 * inv(O2)
end
# the second argument type below is unionall
# ... this allows you to bipass for a specific InvOrAdjOrVecc

function /(O1::NTuple{N,InvOrAdjOrVecc}, O2::InvOrAdjOrVecc) where N
    O1 * inv(O2)
end

function /(O1::InvOrAdjOrVecc, O2::NTuple{N,InvOrAdjOrVecc}) where N
    O1 * inv(O2)
end

function /(O1::NTuple{N,InvOrAdjOrVecc}, O2::NTuple{M,InvOrAdjOrVecc}) where {N,M}
    O1 * inv(O2)
end

# activate the lazy tuple when operating
# ------------------------------------

# function mul!(rf::AbstractVector, O1::NTuple{N,InvOrAdjOrVecc}, f::AbstractVector) where N
#     copyto!(rf, O1 * f) # is there a better way to do this?
# end

function *(O1::NTuple{N,InvOrAdjOrVecc}, f::AbstractVector) where N
    foldr(*, (O1..., f))
end

function \(O1::NTuple{N,InvOrAdjOrVecc}, f::AbstractVector) where N
    inv(O1) * f
end

# function ldiv!(O1::NTuple{N,InvOrAdjOrVecc}, f::AbstractVector) where N
#     mul!(f, inv(O1), copy(f))
# end

# ----------------------------------------------

function show(io::IO, ::MIME"text/plain", A::NTuple{N,InvOrAdjOrVecc}) where N
    println(io, "$N - tuple of VecchiaFactors or Adjoints")
    println(io, typeof(A))
end


text:
using DSP

lowpass = digitalfilter(Lowpass(0.27),Butterworth(4))
function find_peaks(x)#,stage=:final)
  x = filtfilt(lowpass,x)
  dx = diff(x)
  ddx = diff(dx)
  cross0 = find(dx[1:end-1].*dx[2:end] .< 0) .+ 1

  # stage == :cross && return cross0

  peaks = cross0[ddx[cross0 .- 1] .< 0]

  # stage == :max && return peaks

  if !isempty(peaks)
    peaks[x[peaks] .> 0.75maximum(x[peaks])]
  else
    Int[]
  end
end

function guess_source_count(window)
  allpeaks = map(1:size(window,1)) do i
    window[i,find_peaks(window[i,:])]
  end

  maxi = map(allpeaks) do peaks
    isempty(peaks) ? typemin(eltype(peaks)) : maximum(peaks)
  end |> indmax
  length(allpeaks[maxi])
end

function slicerange(x,range,axis)
  dim = axisdim(x,axis)
  args = Array{Any}(ndims(x))
  args .= Colon()
  args[dim] = range
  getindex(x,args...)
end

function map_window(fn,x,window,delta;showprogress=false)
  indices = (1:floor(Int,delta / Δt(x)):ntimes(x))[1:end-1]
  width = floor(Int,window / Δt(x))
  vals = map(indices) do i
    fn(slicerange(x,clamp.(i:i+width,1,ntimes(x)),Axis{:time}))
  end
  AxisArray(vals,Axis{:time}(linspace(0s,duration(x),length(vals))))
end

function source_bumps(cs;window=1s,delta=0.25s)
  indices = 1:floor(Int,delta / Δt(cs)):ntimes(cs)
  y = similar(cs,Axis{:time}(times(cs)[indices]),AxisArrays.axes(cs)[2:end]...)

  for (k,i) in enumerate(indices)
    win = cs[atindex(0s .. window,i)]
    meanwin = meanabs(win,1)
    y[Axis{:time}(k)] = filtfilt(lowpass,copy(meanwin'))'
  end

  y
end

function source_peaks(cs;window=1s,delta=0.25s,stage=:final)
  indices = 1:floor(Int,delta / Δt(cs)):ntimes(cs)
  y = similar(cs,Axis{:time}(times(cs)[indices]),AxisArrays.axes(cs)[2:end]...)

  for (k,i) in enumerate(indices)
    win = cs[atindex(0s .. window,i)]
    meanwin = meanabs(win,1)
    peaks = mapslices(meanwin,2) do col
      ixs = find_peaks(col,stage)
      col .= 0
      col[ixs] = 1.0

      col
    end
    y[Axis{:time}(k)] = peaks
  end

  y
end

function source_count_by_peaks(cs;window=1s,delta=0.25s,buildup=1s)
  map_window(cs[buildup .. last(times(cs))],window,delta) do win
    guess_source_count(meanabs(win,1))
  end
end

max_stars_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_issues_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_forks_repo_name: JuliaTagBot/AuditoryBistabilityLE


text:
#=================================
Methods using Affine Arithmetic
=================================#
using .AffineArithmetic: AFF

"""
    enclose_AffineArithmetic(f::Function, dom::Interval_or_IntervalBox)

Compute a range enclosure using interval arithmetic substitution.

### Input

- `f`   -- function
- `dom` -- hyperrectangular domain, either a unidimensional  `Interval` or
           a multidimensional `IntervalBox`

### Output

An interval representing the range enclosure (minimum and maximum) of `f` over
its domain `dom`.

### Algorithm

The result is obtained by substitution of `dom` on the function `f` and
application of the rules of interval arithmetic.

We refer to the documentation and source code of `IntervalArithmetic` for details.
"""
function enclose_AffineArithmetic(f::Function, dom::Interval_or_IntervalBox)
    _enclose_AffineArithmetic(f, dom)
end

# univariate
function _enclose_AffineArithmetic(f::Function, dom::Interval)
    x = AFF(dom, 1, 1)
    return interval(f(x))
end

# multivariate
function _enclose_AffineArithmetic(f::Function, dom::IntervalBox{N}) where {N}
    x = [AFF(dom[i], N, i) for i in 1:N]
    return interval(f(x...))
end