text:
function run_opf(system::PSY.System, model_type::Type, optimizer; kwargs...)
    pmi_data = get_pm_data(system; kwargs...)
    return run_opf(
        pmi_data,
        model_type,
        optimizer;
        (k for k in kwargs if first(k) ∉ ACCEPTED_PM_DATA_KWARGS)...,
    )
end

function run_mn_opf(system::PSY.System, model_type::Type, optimizer; kwargs...)
    pmi_data = get_pm_data(system; kwargs...)
    return run_mn_opf(
        pmi_data,
        model_type,
        optimizer;
        (k for k in kwargs if first(k) ∉ ACCEPTED_PM_DATA_KWARGS)...,
    )
end

function run_pf(system::PSY.System, model_type::Type, optimizer; kwargs...)
    pmi_data = get_pm_data(system; kwargs...)
    return run_pf(
        pmi_data,
        model_type,
        optimizer;
        (k for k in kwargs if first(k) ∉ ACCEPTED_PM_DATA_KWARGS)...,
    )
end

function run_ac_pf(system::PSY.System, optimizer; kwargs...)
    pmi_data = get_pm_data(system; kwargs...)
    return run_ac_pf(
        pmi_data,
        optimizer;
        (k for k in kwargs if first(k) ∉ ACCEPTED_PM_DATA_KWARGS)...,
    )
end

function run_dc_pf(system::PSY.System, optimizer; kwargs...)
    pmi_data = get_pm_data(system; kwargs...)
    return run_dc_pf(
        pmi_data,
        optimizer;
        (k for k in kwargs if first(k) ∉ ACCEPTED_PM_DATA_KWARGS)...,
    )
end

function compute_ac_pf(system::PSY.System; kwargs...)
    pmi_data = get_pm_data(system; kwargs...)
    return compute_ac_pf(
        pmi_data,
        (k for k in kwargs if first(k) ∉ ACCEPTED_PM_DATA_KWARGS)...,
    )
end

function compute_dc_pf(system::PSY.System; kwargs...)
    pmi_data = get_pm_data(system; kwargs...)
    return compute_dc_pf(
        pmi_data,
        (k for k in kwargs if first(k) ∉ ACCEPTED_PM_DATA_KWARGS)...,
    )
end

function run_model(system::PSY.System, model_type::Type, optimizer, build_method; kwargs...)
    pmi_data = get_pm_data(system; kwargs...)
    return run_model(
        pmi_data,
        model_type,
        optimizer,
        build_method;
        (k for k in kwargs if first(k) ∉ ACCEPTED_PM_DATA_KWARGS)...,
    )
end

function instantiate_model(system::PSY.System, model_type::Type, build_method; kwargs...)
    pmi_data = get_pm_data(system; kwargs...)
    return instantiate_model(
        pmi_data,
        model_type,
        build_method;
        (k for k in kwargs if first(k) ∉ ACCEPTED_PM_DATA_KWARGS)...,
    )
end

max_stars_repo_name: NREL-SIIP/PowerModelsInterface

max_issues_repo_name: NREL-SIIP/PowerModelsInterface

max_forks_repo_name: NREL-SIIP/PowerModelsInterface


text:
module FluxDash

#use packages
using CSV: Options
using Dash, DashBootstrapComponents, DashHtmlComponents, DashCoreComponents
using DashTable, DashDaq
using Core: Typeof
using Base64, Dates, CSV, DataFrames
using PlotlyJS
using Random, StatsBase
using Flux
using Flux: Data.DataLoader

#Initialize channel to get data during training
chn = Channel{Tuple{Float64, Int, Chain}}(Inf)
m = nothing #Initialized Chain
in_labels = nothing #Intialized label name vector corresponds to input data
out_labels = nothing #Intialized label name vector corresponds to output data
X_test = nothing #Intialized input test data
y_test = nothing #Intialized output test data
ep = 0
st = false
x_graph = []
y_graph = []

include("utils.jl")
include("render_dhc.jl")
include("generate_nn.jl")

function make_app()
	#CDN for fonts
	external_stylesheets = "https://cdnjs.cloudflare.com/ajax/libs/font-awesome/4.7.0/css/font-awesome.min.css"

	global chn, m, in_labels, out_labels, X_test, y_test,ep, st, x_graph, y_graph
	#Initialize channel to get data during training
	chn = Channel{Tuple{Float64, Int, Chain}}(Inf)
	m = nothing #Initialized Chain
	in_labels = nothing #Intialized label name vector corresponds to input data
	out_labels = nothing #Intialized label name vector corresponds to output data
	X_test = nothing #Intialized input test data
	y_test = nothing #Intialized output test data
	ep = 0
	st = false
	x_graph = []
	y_graph = []
	app = dash(external_stylesheets=[dbc_themes.BOOTSTRAP,external_stylesheets],
	    suppress_callback_exceptions=true)

	#setup dashboard
	app.layout = html_div([
	    dcc_store(id="dflx-data-memory"), #memory to store data
	    dcc_store(id="dflx-testdataX-memory"), #memory to store input labels
	    dcc_store(id="dflx-testdatay-memory"), #memory to store output labels
	    dcc_store(id="dflx-inlabels-memory"), #memory to store input test data
	    dcc_store(id="dflx-outlabels-memory"), #memory to store output test data
	    render_all() #To render all other components
	])

	#=
	This callback aims to upload data as CSV file and store data in memory.
	Moreover, this callback render interactive graph and table with loaded data
	=#
	callback!(
	    app,
	    Output("dflx-data-memory", "data"),
	    Output("dflx-output-data-upload", "children"),
	    Output("dflx-output-plot", "children"),
	    Output("dflx-nn-tab", "disabled"),
	    Output("dflx-input-name-list", "children"),
	    Output("dflx-output-name-list", "children"),
	    Input("dflx-upload-data", "contents"),
		Input("dflx-def-upload-data", "n_clicks"),
	    State("dflx-upload-data", "filename"),
		State("dflx-def-upload-data-drop", "value")
	) do contents, sel_clicks, filename, drop_val

		ctx1 =callback_context()
	    if !(ctx1 isa Nothing)
			if !isempty(ctx1.triggered)
				if ctx1.triggered[1].prop_id == "dflx-upload-data.contents"

					df = parse_contents(contents, filename)
				elseif ctx1.triggered[1].prop_id == "dflx-def-upload-data.n_clicks"
					if drop_val == "1"
						df = CSV.read(download("https://raw.githubusercontent.com/efmanu/fluxdash/master/FluxDash/datasets/real_estate.csv"), DataFrame)
						df[!,"id"] = 1:nrow(df)
						filename = "Real_Estate.csv"
					elseif drop_val == "2"
						df = CSV.read(download("https://raw.githubusercontent.com/plotly/datasets/master/solar.csv"), DataFrame)
						df[!,"id"] = 1:nrow(df)
						filename = "Solar.csv"
					elseif drop_val == "3"
						df = CSV.read(download("https://raw.githubusercontent.com/efmanu/fluxdash/master/FluxDash/datasets/kc_house_data.csv"), DataFrame)
						df[!,"id"] = 1:nrow(df)
						filename = "House_Price.csv"
					else
						throw(PreventUpdate())
					end
				end
				return_data = render_table(df, filename; n_rows = 5)
				return_graph = render_graph_title(df)
				row_names = names(df)
				row_names = row_names[row_names .!= "id"]
				return (
					df, return_data,
					return_graph, false,
					render_mult_graphtitle(row_names, id_val="dflx-nninput-dropdown"),
					render_mult_graphtitle(row_names, id_val="dflx-nnoutput-dropdown")
				)
			else
				throw(PreventUpdate())
			end
		else
			throw(PreventUpdate())
	    end
	end

	#=
	This callback aims to render graph interactively based on selected labels
	=#
	callback!(app,
	    Output("dflx-plot-graph", "figure"),
	    Input("dflx-graph-dropdown-x", "value"),
	    Input("dflx-graph-dropdown-y", "value"),
	    State("dflx-data-memory", "data")
	) do valuex, valuey, df
	    new_tp = DataFrame(;zip(Tuple(Symbol.(df.colindex.names)), Tuple(df.columns))...)
	    return  render_graph(new_tp, valuex, valuey)
	end

	#=
	This callback aims to render Neural network hideen layers
	=#
	callback!(app,
	    Output("dflx-nn-layer-count", "data"),
	    Output("dflx-hidden-layers", "children"),
	    Input("dflx-add-hidden-layer", "n_clicks"),
	    Input("dflx-remove-hidden-layer", "n_clicks"),
	    State("dflx-nn-layer-count", "data")
	) do nclicks_plus, nclicks_minus, layer_count
	    ctx = callback_context()
	    if isempty(ctx.triggered)
	        throw(PreventUpdate())
	    end
	    if layer_count isa Nothing
	        layer_count = 0
	    end
	    if split(ctx.triggered[1].prop_id,".")[1] == "dflx-add-hidden-layer"
	        layer_count += 1
	    else
	        layer_count -= 1
	    end
	    return layer_count, render_nn_layer(layer_count)
	end

	#=
	This callback aims to start neural network training
	Also this callback graph between number of  epoch vs loss
	=#
	callback!(app,
	    Output("dflx-training-updates", "children"),
	    Output("dflx-testdataX-memory", "data"),
	    Output("dflx-testdatay-memory", "data"),
	    Output("dflx-inlabels-memory", "data"),
	    Output("dflx-outlabels-memory", "data"),
	    Input("dflx-nnlayer-submit", "n_clicks"),
	    State("dflx-nninput-dropdown", "value"),
	    State("dflx-nnoutput-dropdown", "value"),
	    State("dflx-hidden-layers", "children"),
	    State("dflx-data-memory", "data"),
	    State("dflx-nntrain-percent", "value"),
	    State("dflx-nntrain-epoch", "value"),
	    State("dflx-nntrain-optimizer","value"),
	    State("dflx-nntrain-lrate","value"),
	    prevent_initial_call=true) do nclicks, labels_in, labels_out, child, df, tr_len, epval, opt, eta
	    if child isa Nothing
	        throw(PreventUpdate())
	    end
	    global in_labels
	    global out_labels
	    in_labels = labels_in
	    out_labels = labels_out
	    global X_test
	    global y_test
	    new_tp = DataFrame(;zip(Tuple(Symbol.(df.colindex.names)), Tuple(df.columns))...)
	    X_train, X_test, y_train, y_test = format_nn_data(new_tp, in_labels, out_labels; training_percent=tr_len)
	    hidden_outs = [ch.props.children[1].props.children[1].props.children[2].props.value for ch in child]
	    hidden_activations = [ch.props.children[1].props.children[1].props.children[3].props.value for ch in child]
		reset_vars()
	    train_nn(X_train, y_train, in_labels, out_labels, hidden_outs, chn, hidden_activations; training_percent = tr_len, epoc = epval, opt = opt, eta = eta)
	    return render_training(), X_test, y_test, in_labels, out_labels
	end

	#=
	This callback updates the loss vs epoch graph every 100 ms and stop updation if
	training finished. Also renders test results, if training is finished
	=#
	callback!(app,
	  Output("live-update-graph", "figure"),
	  Output("interval-component", "disabled"),
	  Output("dflx-pred-tab", "disabled"),
	  Output("dflx-testing-section", "children"),
	  Output("dflx-pred-tab", "children"),
	  Input("interval-component", "n_intervals"),
	  State("live-update-graph", "figure"),
	  State("dflx-nntrain-epoch", "value"),
	  prevent_initial_call=true) do n, stg, epval
	    global m, ep, st
		global x_graph, y_graph
	    global X_test
	    global y_test
	    if isready(chn)
	    	l,ep,m = take!(chn)
			if ep >= epval
				st = true
			end
			append!(x_graph, ep)
			append!(y_graph, l)
		end

	    if st
	        test_render = render_testing(m, X_test,y_test)
	        pred_render = render_prediction_componnets(in_labels)
	    else
	        test_render = dbc_progress(value=(ep/epval)*100)
	        pred_render = html_div()
	    end
	    return Dict(
	    "data" => [
	        Dict(
	            "x" => x_graph,
	            "y" => y_graph,
	            "mode" => "line",
	        ),
	    ]
	    ), st, !st, test_render, pred_render

	end

	#=
	This callback render predcted output by inputing user data to trained model
	=#
	callback!(app,
	    Output("pred_out", "children"),
	    Input((index = ALL, type = "preddone"), "value")
	) do in_values
	    global m
	    global out_labels
	    if (in_values isa Vector) && !(nothing in in_values)
	        prediction =  m(in_values)
	        return [
	            html_div([
	                dcc_input(
	                    id = (index = idx, type = "preddout"),
	                    value = val,
	                    disabled = true
	                )
	            ]) for (idx, val) in enumerate(prediction)]
	    else
	        return html_div("nothing")
	    end
	end
	return app
end

end # module

ext: jl

max_stars_repo_name: DEICIDIA/fluxdash

max_issues_repo_name: DEICIDIA/fluxdash

max_forks_repo_name: DEICIDIA/fluxdash


text:
using Formatting
using ArgParse
using JSON
using NCDatasets

function parse_commandline()

    s = ArgParseSettings()
    @add_arg_table s begin

        "--input-file"
            help = "CESM formated domain file."
            arg_type = String
            required = true

        "--output-file"
            help = "Output SCRIP file."
            arg_type = String
            default = ""

        "--center-lon"
            help = "Longitude variable name."
            arg_type = String
            required = true

        "--center-lat"
            help = "Latitude variable name."
            arg_type = String
            required = true

        "--corner-lon"
            help = "Longitude variable name."
            arg_type = String
            required = true

        "--corner-lat"
            help = "Latitude variable name."
            arg_type = String
            required = true

        "--mask"
            help = "Mask variable name."
            arg_type = String
            default = "mask"

        "--area"
            help = "Area variable name."
            arg_type = String
            default = "area"

        "--angle-unit"
            help = "Unit of lat, lon. `deg` or `rad`. Default is `deg`."
            arg_type = String
            default = "deg"

        "--mask-value"
            help = "Mask value. Active grids will be set to 1. Inactive grids will be set to 0."
            arg_type = Float64
            required = true

        "--all-active"
            action = :store_true

    end

    return parse_args(ARGS, s)
end

println("Running ", @__FILE__)

parsed = parse_commandline()
print(json(parsed, 4))

if parsed["output-file"] == ""
    parsed["output-file"] = format("{:s}.SCRIP.nc", basename(splitext(parsed["input-file"])[1]))
end

Dataset(parsed["input-file"]) do ds

    c = ( parsed["angle-unit"] == "deg" ) ? π/180.0 : 1.0

    global grid_center_lat = reshape(convert(Array{Float64}, replace(ds[parsed["center-lat"]][:], missing=>NaN)) * c, :)
    global grid_center_lon = reshape(convert(Array{Float64}, replace(ds[parsed["center-lon"]][:], missing=>NaN)) * c, :)
    global grid_corner_lat = reshape(convert(Array{Float64}, replace(ds[parsed["corner-lat"]][:], missing=>NaN)) * c, 4, :)
    global grid_corner_lon = reshape(convert(Array{Float64}, replace(ds[parsed["corner-lon"]][:], missing=>NaN)) * c, 4, :)
    global grid_imask      = reshape(convert(Array{Int64},   replace(ds[parsed["mask"]][:],       missing=>NaN)), :)
    global grid_area       = reshape(convert(Array{Float64}, replace(ds[parsed["area"]][:],       missing=>NaN)), :)

    global grid_dims = collect(Int64, size(ds[parsed["center-lat"]]))
    global grid_size = reduce(*, grid_dims)

    if parsed["all-active"]
        println("--all-active is on.")
        grid_imask .= 1.0
    else
        _grid_imask = copy(grid_imask)
        grid_imask[_grid_imask .== parsed["mask-value"]] .= 1.0
        grid_imask[_grid_imask .!= parsed["mask-value"]] .= 0.0
    end
end

println(grid_dims)

Dataset(parsed["output-file"], "c") do ds

    defDim(ds, "grid_size", grid_size)
    defDim(ds, "grid_corners", 4)
    defDim(ds, "grid_rank", 2)

    for (varname, vardata, varnctype, vardim, attrib) in [
        ("grid_dims",        grid_dims,       Int32,   ("grid_rank",),                 Dict()),
        ("grid_center_lat",  grid_center_lat, Float64, ("grid_size",),                 Dict("units" => "radians")),
        ("grid_center_lon",  grid_center_lon, Float64, ("grid_size",),                 Dict("units" => "radians")),
        ("grid_imask",       grid_imask,      Int32,   ("grid_size",),                 Dict("units" => "unitless")),
        ("grid_corner_lat",  grid_corner_lat, Float64, ("grid_corners", "grid_size",), Dict("units" => "radians")),
        ("grid_corner_lon",  grid_corner_lon, Float64, ("grid_corners", "grid_size",), Dict("units" => "radians")),
        ("grid_area",        grid_area,       Float64, ("grid_size",),                 Dict("units" => "radians^2", "long_name" => "area weights")),
    ]
        println("Doing varname:", varname)
        var = defVar(ds, varname, varnctype, vardim)

        for (k, v) in attrib
            var.attrib[k] = v
        end

        rng = []
        for i in 1:length(vardim)-1
            push!(rng, Colon())
        end

        push!(rng, 1:size(vardata)[end])
        var[rng...] = vardata

    end

end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
using OpenVR

eError = Ref{ETrackedPropertyError}()
ivrsystem = VR_Init(eError, OpenVR.VRApplication_Scene, C_NULL)

if ivrsystem != C_NULL

eError = Ref{OpenVR.ETrackedPropertyError}()
bufsize = UInt32(1024)
buf = zeros(UInt8,bufsize)
len = OpenVR.GetStringTrackedDeviceProperty(ivrsystem,OpenVR.k_unTrackedDeviceIndex_Hmd,OpenVR.Prop_TrackingSystemName_String,pointer(buf),bufsize,pointer_from_objref(eError))
trackingsystemname = unsafe_string(pointer(buf),len-1)

end


text:
export time_loop

function time_loop(n)
    global x, dr, F
    global v, storage_e
    global trajzmin, dhdea
    global trajtheta, dhdv, blk
    global eigvec_H, dhp
    global eigval_H, surfpinit
    global phi, psi, storage_aop, surfp, storage_op, storage_x, storage_v
    dr = x[1, :] - x[2, :]
    v_cond = (v[1, 3]*m_N + v[2, 3]*m_O)/(m_N + m_O)
    z_cond = (x[1, 3] + x[2, 3])/2.0
    if z_cond > z_end && v_cond > 0.0
        println("Boundary reached: End simulation")
        return
    end

    # Record minimial distance of z approach and orientational angle
    min_z_no = minimum(x[1:2, 3])
    if min_z_no <= trajzmin
        trajzmin = min_z_no
        trajtheta = acos((dr[3])/sqrt(dr' * dr))
    end

    # get state corresponding to current surface
    for i in 1:Ne
        phi[:, i] = eigvec_H[:, surfp[i]]
    end

    # Propagate electronic Hamiltonian
    psi = eigvec_H' * psi
    psi = @. exp(-1im * eigval_H * dt/hbar)

    psi = repeat(psi, Ne)
    psi = reshape(psi, (41, Ne))
    psi = @.psi * psi
    psi = eigvec_H * psi

    #recordings
    storage_aop[n] = sum(abs.(psi[1, :]).^2)
    storage_op_temp = sum(abs.((eigvec_H' * psi).^2), dims=2)
    storage_op[:, n] = storage_op_temp

    rtemp = sqrt(dr' * dr)
    storage_e[1, n] = 0.5 * (m_N + m_O)*sum(((m_N * v[1, :] + m_O * v[2, :])/(m_N + m_O)).^2)
    T_vib = 0.50 * μ * sum(((v[1, :] - v[2, :]).*(x[1, :] - x[2, :])/rtemp).^2)
    U_vib = F_n *(1.0 - exp(-δ_n * (rtemp -  r_0_NO)))^2
    E_vib = U_vib + T_vib
    storage_e[2, n] = E_vib
    T_tot = 0.5 * mass_arr[1] * sum(v[1, :].^2) + 0.5*mass_arr[2]*sum(v[2, :].^2)
    storage_e[3, 1] = T_tot - T_vib - storage_e[1, n]
    storage_e[4, 1] = sum(eigval_H[surfp]) - sum(eigval_H[surfpinit])

    #Propagate nuclear motion

    vdot = @. 1.0/m_spread * F
    x = @. x + v*dt + 0.5*vdot * dt^2
    vtemp = @. v + 0.5 *vdot * dt

    # Get adiabatic eigenvalues, eigenvectors
    hp = get_E_all(x) #Hamiltonian parameters
    dhp = get_F_all(x)

    # dhp = get_F_all(x, dhp)

    H = h0 + vm * hp[3] / sqrt(delta_E)
    H[1, 1] = H[1, 1] + hp[2] - hp[1]
    eig_vec_hold = eigvec_H

    eigval_H = eigvals(H)
    eigvec_H = eigvecs(H)

    for j in 1:Ms+1
        vHmaxloc = argmax(abs.(eig_vec_hold' * eigvec_H[:, j]))
        sign_b = sign(dot(eig_vec_hold[:, vHmaxloc], eigvec_H[:, j]))
        eigvec_H[:, j] = eigvec_H[:, j]*sign_b
    end
    # Get dH / d Ep matrices where Ep is parameter in Newns-Anderson Hamiltonian
    # (i.e. Ep = V or Ep = E_a = (EI-EN) )
    blk .= 0
    for j in 1:Ms+1
        dhdea[:, j] = @. eigvec_H[1, :] * eigvec_H[1, j]
    end
    dhdv = eigvec_H' * (vm * eigvec_H)/sqrt(delta_E)

    #calculate forces

    F .= 0
    for j in 1:Ne
        temp_arr = @. dhdea[surfp[j], surfp[j]] * (dhp[2, :] - dhp[1, :])
        temp_arr2 = reshape(temp_arr, (3, N+2))
        temp_arr3 = permutedims(temp_arr2, (2, 1))
        F = F + temp_arr3

        temp_arr = @. dhdv[surfp[j], surfp[j]] * dhp[3, :]
        temp_arr2 = reshape(temp_arr, (3, N+2))
        temp_arr3 = permutedims(temp_arr2, (2, 1))
        F = F + temp_arr3
    end
    temp_arr = reshape(dhp[1,:], (3, N+2))
    temp_arr2 = permutedims(temp_arr, (2, 1))
    F = F + temp_arr2

    #fix last layer of grid
    F[399:end, :] .= 0.0
    v[399:end, :] .= 0.0
    F = -F

    v = @. vtemp + 0.5 / m_spread * F * dt

    #storage
    storage_x[1:3, n] = x[1, :]
    storage_x[4:6, n] = x[2, :]

    storage_v[1:3, n] = v[1, :]
    storage_v[4:6, n] = v[2, :]
    return
end


text:
function lintref(ex::Expr, ctx::LintContext)
    sub1 = ex.args[1]
    guesstype(ex, ctx) # tickle the type checks on the expression
    lintexpr(sub1, ctx)
    sub1type = guesstype(sub1, ctx)
    for i=2:length(ex.args)
        if ex.args[i] == sub1
            if length(ex.args)==2 && sub1type <: Array{Int}
                msg(ctx, :I473, sub1, "value at position #$(i-1) is the referenced " *
                    "$(sub1). OK if it represents permutations")
            else # almost certain it's an error
                msg(ctx, :E434, sub1, "value at position #$(i-1) is the " *
                    "referenced $(sub1). Possible typo?")
            end
        end
        lintexpr(ex.args[i], ctx)
    end
end

function linttyped_vcat(ex::Expr, ctx::LintContext)
    lintref(ex, ctx)
end


text:
using JudiLing
using Test

@testset "test_combo tests" begin
    data_path = joinpath(@__DIR__, "data", "latin_mini.csv")
    data_output_dir = joinpath(@__DIR__, "out_data")
    output_dir = joinpath(@__DIR__, "out")
    JudiLing.test_combo(
        :train_only,
        val_sample_size = 100,
        data_path = data_path,
        data_prefix = "latin",
        data_output_dir = data_output_dir,
        n_features_columns = ["Lexeme","Person","Number","Tense","Voice","Mood"],
        n_grams_target_col = :Word,
        n_grams_tokenized = false,
        grams = 3,
        n_features_base = ["Lexeme"],
        n_features_inflections = ["Person","Number","Tense","Voice","Mood"],
        output_dir = output_dir,
        verbose = false
        )

    JudiLing.test_combo(
        :random_split,
        val_sample_size = 100,
        data_path = data_path,
        data_prefix = "latin",
        data_output_dir = data_output_dir,
        n_features_columns = ["Lexeme","Person","Number","Tense","Voice","Mood"],
        n_grams_target_col = :Word,
        n_grams_tokenized = false,
        if_combined = true,
        grams = 3,
        n_features_base = ["Lexeme"],
        n_features_inflections = ["Person","Number","Tense","Voice","Mood"],
        output_dir = output_dir,
        verbose = false
        )

    JudiLing.test_combo(
        :careful_split,
        val_sample_size = 100,
        data_path = data_path,
        data_prefix = "latin",
        data_output_dir = data_output_dir,
        n_features_columns = ["Lexeme","Person","Number","Tense","Voice","Mood"],
        n_grams_target_col = :Word,
        n_grams_tokenized = false,
        grams = 3,
        n_features_base = ["Lexeme"],
        n_features_inflections = ["Person","Number","Tense","Voice","Mood"],
        output_dir = output_dir,
        verbose = false
        )

    rm(data_output_dir, force = true, recursive = true)
    rm(output_dir, force = true, recursive = true)
end


text:
include("../include.jl")

results = Dict{String, Dict{String,problem_summary}}()
#results["IPOPT"] = load_results("other_solver_results/ipopt-results.txt");
#results["big_run"] = load("../results/OLD-august-21/big_run/summary.jld", "summary")
#results["IPOPT"] = convert_JuMP(load("../results/OLD-august-21/ipopt_test3/summary.jld", "summary"))

#results["IPOPT no perturb"] = convert_JuMP(load("../results/ipopt/large_no_perturb/summary.jld", "summary"))
#results["IPOPT"] = convert_JuMP(load("../results/ipopt/large/summary.jld", "summary"))
#results["static"] = load("../results/one_phase/large/summary.jld", "summary")
#results["large_dynamic"] = load("../results/one_phase/large_dynamic/summary.jld", "summary")
#results["dynamic"] = load("../results/one_phase/large_dynamic_1hr/summary.jld", "summary")
#results["one phase 2"] = load("../results/one_phase/plain/summary.jld", "summary")

#results["1"] = load("../results/one_phase/sept_1_corrections/summary.jld", "summary")
#results["2"] = load("../results/one_phase/sept_2_corrections/summary.jld", "summary")
#results["3"] = load("../results/one_phase/sept_3_corrections/summary.jld", "summary")
#results["4"] = load("../results/one_phase/sept_4_corrections/summary.jld", "summary")

# #=
#results["0"] = load("../results/one_phase/dual_ls_0/summary.jld", "summary")
#results["1"] = load("../results/one_phase/dual_ls_1/summary.jld", "summary")
#results["2"] = load("../results/one_phase/dual_ls_2/summary.jld", "summary")
# =#

# infeasible problems
#=
results["one phase"] = load("../results/one_phase/infeas-3/summary.jld", "summary")
results["IPOPT"] = convert_JuMP(load("../results/ipopt/infeas-3/summary.jld", "summary"))
=#

# #=
#results["max step"] = load("../results/one_phase/sept_max_step_stable/summary.jld", "summary")
#results["log barrier"] = load("../results/one_phase/sept_log_barrier_stable/summary.jld", "summary")
#results["filter"] = load("../results/one_phase/sept_3_corrections/summary.jld", "summary")
# =#

#=
results["regularizer"] = load("../results/one_phase/sept_no_regularizer/summary.jld", "summary")
results["no_regularizer"] = load("../results/one_phase/sept_3_corrections/summary.jld", "summary")
=#

#results["eta_affine"] = load("../results/one_phase/eta/affine/summary.jld", "summary")
#results["eta_mehrotra_stb"] = load("../results/one_phase/eta/mehrotra_stb/summary.jld", "summary")

#=
results["norm tol"] = load("../results/one_phase/norm_tol/summary.jld", "summary")
results["high tol"] = load("../results/one_phase/high_tol/summary.jld", "summary")
=#

#=
results["medium"] = convert_JuMP(load("../results/ipopt/medium/summary.jld", "summary"))
results["medium high tol"] = convert_JuMP(load("../results/ipopt/medium_high_tol/summary.jld", "summary"))
results["medium perturb"] = convert_JuMP(load("../results/ipopt/medium_perturb/summary.jld", "summary"))
results["medium perturb high tol"] = convert_JuMP(load("../results/ipopt/medium_perturb_high_tol/summary.jld", "summary"))
=#

# #=
#results["dynamic"] = load("../results/one_phase/sept_dynamic/summary.jld", "summary")
#results["static"] = load("../results/one_phase/sept_3_corrections/summary.jld", "summary")
# =#

# compare one phase and ipopt

#results["sept 2"] = load("../results/one_phase/sept_2_corrections/summary.jld", "summary")
#results["Oct21"] = load("../results/one_phase/Oct21/summary.jld", "summary")

#=
results["Oct21_penalty"] = load("../results/one_phase/Oct21_penalty/summary.jld", "summary")
results["Oct21_no_penalty"] = load("../results/one_phase/Oct21_no_penalty/summary.jld", "summary")
results["Oct21_penalty_full"] = load("../results/one_phase/Oct21_penalty_full/summary.jld", "summary")
results["Oct21_penalty_prox"] = load("../results/one_phase/Oct21_penalty_prox/summary.jld", "summary")
=#

#results["Oct22_a_penalty_only"] = load("../results/one_phase/Oct22/a_penalty_only/summary.jld", "summary")
#results["Oct22_a_penalty_only_no_prox"] = load("../results/one_phase/Oct22/a_penalty_only_no_prox/summary.jld", "summary")
#results["Oct22_penalty"] = load("../results/one_phase/Oct22/penalty/summary.jld", "summary")
#results["Oct22_penalty_prox"] = load("../results/one_phase/Oct22/penalty_no_prox/summary.jld", "summary")
#results["Oct23"] = load("../results/one_phase/Oct23/summary.jld", "summary")
#results["Oct23_MORE_TIME"] = load("../results/one_phase/Oct23_MORE_TIME/summary.jld", "summary")
#results["Oct26"] = load("../results/one_phase/Oct26/summary.jld", "summary")
#results["Oct28"] = load("../results/one_phase/Oct28/summary.jld", "summary")

#=
results["init0.1"] = load("../results/one_phase/init/0.1/summary.jld", "summary")
results["init1"] = load("../results/one_phase/init/1/summary.jld", "summary")
results["init10"] = load("../results/one_phase/init/10/summary.jld", "summary")
results["init100"] = load("../results/one_phase/init/100/summary.jld", "summary")
results["init1000"] = load("../results/one_phase/init/1000/summary.jld", "summary")
=#

min_num = 100
max_num = 10000
lrg_problems = lrg_problem_func(min_num, max_num)
problem_list = CUTEst.select(custom_filter=lrg_problems);
problem_list = convert(Array{String,1},problem_list);

results = restrict_problems(results, problem_list)

#results["stable_first/false"] = load("results/stable_first/false/summary.jld", "summary")
#results["stable_first/true"] = load("results/stable_first/true/summary.jld", "summary")

if true
#error_free_results = remove_errors(results, [:NaN_ERR,:Error])
#error_free_results = remove_errors(results, [:MAX_TIME])
#error_free_results = remove_errors(results, [:NaN_ERR, :MAX_TIME])
error_free_results = remove_errors(results, [])
overlapping_results = overlap(error_free_results)
#overlapping_results = restrict_to_set(overlapping_results,[:primal_infeasible])

overlapping_results = restrict_to_set(overlapping_results,[:Optimal,:primal_infeasible,:dual_infeasible])

elseif false
overlapping_results = overlap(results)
infeas_results = select(overlapping_results, [:primal_infeasible])
overlapping_results = union_results(infeas_results, results)
show_results(overlapping_results)
else
overlapping_results = overlap(results)
end

##
## FAILURES
##

println("kkt")
for (method_name, sum_data) in overlapping_results
    println(method_name, "=" ,tot(sum_data,[:Optimal]), " out of ", length(sum_data))
end

println("unbounded")
for (method_name, sum_data) in overlapping_results
    println(method_name, "=" ,tot(sum_data,[:dual_infeasible]), " out of ", length(sum_data))
end

println("primal infeasible")
for (method_name, sum_data) in overlapping_results
    println(method_name, "=" ,tot(sum_data,[:primal_infeasible]), " out of ", length(sum_data))
end

println("failures")
for (method_name, sum_data) in overlapping_results
    println(method_name, "=" ,tot_failures(sum_data), " out of ", length(sum_data))
end

@show shared_failures(overlapping_results)
@show list_combined_success(overlapping_results)

#outcomes_table(overlapping_results)

println("failures")
for (method_name, sum_data) in overlapping_results
    println(method_name)#
    @show failure_causes(sum_data)
    #, "=" ,tot_failures(sum_data), " out of ", length(sum_data))
end

print_failure_problems(overlapping_results);
its, best, ratios, times = compute_its_etc(overlapping_results,MAX_IT=3000);

##
## PLOT ITERATION RATIO CURVES
##

using PyPlot

plot_iteration_ratios(its, best, ratios)
plot_iterations(its, best, ratios, 10000)

##
## FUNCTION VALUES
##
problem_list = collect(keys(first(overlapping_results)[2]))
method_list = collect(keys(overlapping_results))

winners = Dict()
winners["tie"] = 0
winners["fail"] = 0
begin
  for method_name in method_list
      print(pd(""), pd(method_name))
      winners[method_name] = 0
  end
  print("\n")

  for problem_name in problem_list
      print(pd(problem_name))
      best = "fail"
      both_kkt = true
      fbest = 1e14
      for method_name in method_list
        info = overlapping_results[method_name][problem_name]
        fval = info.fval
        con_vio = info.con_vio
        print(rd(con_vio))
        print(rd(fval))

        if info.status == :Optimal
          if con_vio > 1e-6
            both_kkt = false
            println("CON VIOLATION!!")
          else
            @assert(con_vio < 1e-6)
            ftol = 1e-1
            ftol_scaled = ftol * (1.0 + abs(fbest))
            if fval <= fbest - ftol_scaled
              fbest = fval
              best = method_name
            elseif fbest - ftol_scaled < fval && fval <= fbest + ftol_scaled
              fbest = fval
              best = "tie"
            end
          end
        else
          both_kkt = false
        end
      end
      if both_kkt
        winners[best] += 1
      end
      print(pd(best))
      print("\n")
  end
end

winners

function best_fvals(method_name1, method_name2)
    better_fval = Dict()
    succeed = Dict()

    feas_tol = 1e-4
    fval_tol = 1e-1

    better_fval[method_name1] = 0
    better_fval[method_name2] = 0

    succeed[method_name1] = 0
    succeed[method_name2] = 0

    for problem_name in problem_list
      fval1 = overlapping_results[method_name1][problem_name].fval
      con_vio1 = overlapping_results[method_name1][problem_name].con_vio

      fval2 = overlapping_results[method_name2][problem_name].fval
      con_vio2 = overlapping_results[method_name2][problem_name].con_vio

      if con_vio1 < feas_tol && con_vio2 < feas_tol
        fbest = min(fval1,fval2)
        ftol_scaled = fval_tol * (1.0 + abs(fbest))

        if fval1 < fval2 - ftol_scaled
          better_fval[method_name1] += 1
        elseif fval2 < fval1 - ftol_scaled
          better_fval[method_name2] += 1
        end
      end

      if con_vio1 < feas_tol
        succeed[method_name1] += 1
      end

      if con_vio2 < feas_tol
        succeed[method_name2] += 1
      end
    end

    return better_fval, succeed
end

method_name1 = "IPOPT"
method_name2 = "IPOPT no perturb"
best_fvals(method_name1, method_name2)

##
## PLOT TIME CURVES
##
best = best_its(times)

min_y = 1.0

ratios = Dict()
y_vals = collect(1:length(best)) / length(best)
for (method_name, val) in its
  ratios[method_name] = times[method_name] ./ best
  semilogx(sort(ratios[method_name]), y_vals, label=method_name, basex=2)
  min_y = min(min_y,sum(ratios[method_name] .== 1.0) / length(best))
end
ax = gca()
ax[:set_xlim]([1.0,2^8.0])
ax[:set_ylim]([min_y,1.0])

#ax[:xaxis][:ticker] = 0.5
#title("Comparsion on 45 CUTEst problems")
xlabel("times ratio to best solver")
ylabel("proportion of problems")

legend()


text:
# Julia code to compute the correlation coefficient
using Distributions
using Plots

x = rand(MvNormal([0, 0], [3 1; 1 1]), 1000)
# x = rand(MvNormal([0, 0], [3 0; 0 3]), 1000)
# x = rand(MvNormal([0, 0], [3 2.9; 2.9 3]), 1000)
scatter(x[1, :], x[2, :])

σ₁ = std(x[1, :])
σ₂ = std(x[2, :])
μ₁ = mean(x[1, :])
μ₂ = mean(x[2, :])
Exy = mean(x[1, :] .* x[2, :])
ρ = (Exy - μ₁ * μ₂) / (σ₁ * σ₂)

####################################################################

# Julia code to compute a mean vector.
using Statistics

X = randn(100, 2)
mean(X, dims=1)

####################################################################

# Julia code to compute covariance matrix.
using Statistics

X = randn(100, 2)
cov(X)

####################################################################

# Julia code: Overlay random numbers with the Gaussian contour.
using Distributions
using Plots

p = MvNormal([0, 0], [0.25 0.3; 0.3 1])

X = rand(p, 1000)
x₁ = -2.5:0.01:2.5
x₂ = -3.5:0.01:3.5

f(x₁, x₂) = pdf(p, [x₁, x₂])

scatter(X[1, :], X[2, :], legend=false)
contour!(x₁, x₂, f, linewidth=2)

####################################################################

# Julia code: Gaussian(0,1) --> Gaussian(mu,sigma)
using Distributions

x = rand(MvNormal([0, 0], [1 0; 0 1]), 1000)
Σ = [3 -0.5; -0.5 1]
μ = [1, -2]
y = Σ^(1/2) * x .+ μ

####################################################################

# Julia code: Gaussian(mu,sigma) --> Gaussian(0,1)
using Distributions

y = rand(MvNormal([1, -2], [3 -0.5; -0.5 1]), 100)
μ = mean(y, dims=2)
Σ = cov(y')
x = Σ^(-1/2) * (y .- μ)

####################################################################

# Julia code to perform the principal component analysis
using LinearAlgebra
using Distributions

x = rand(MvNormal([0, 0], [2 -1.9; -1.9 2]), 1000)
Σ = cov(x')
S, U = eigen(Σ)
U[:, 1]
U[:, 2]


text:
using Documenter, ImageDraw

makedocs(sitename="Documentation",
            format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"))
deploydocs(repo = "github.com/JuliaImages/ImageDraw.jl.git")


text:
include("gkMCSetup.jl")

function makeRecycleLTD(ys::Vector{Float64}, N::Int64)
  rpio::RPIO{Float64Particle} = RPIO{Float64Particle}(N, length(ys), 1)
  function ltd(in::SVector{d, Float64})
    θ::GKθ = toGKθ(in)
    lp::Float64 = lglogprior(θ)
    if lp == -Inf return -Inf end
    model::Model = makeGKModel(θ, ys)
    return rpEstimate(model, rpio)
  end
end

logtargetRecycle = makeRecycleLTD(ys, N)
PRecycle = makeAMKernel(logtargetRecycle, d)

Random.seed!(12345)
setRNGs(54321)

recycleElapsed = @elapsed (chainRecycle = simulateChainProgress(PRecycle, t0, chainLength))
rar = PRecycle(:acceptanceRate)
recycleCovEstimate = PRecycle(:covEstimate)

@save "data/gkMCRecycle.jld" chainRecycle rar recycleCovEstimate recycleElapsed


text:
FEM.linquadinterpmod()
using FEM.LinQuadInterpMod

module LinQuadMod

using FEM.LinQuadInterpMod

import FEM: createinterp, creategps, createstorage, get_field,
            Bmatrix, doftypes, get_ndofs, get_geotype, get_ref_area, dNdxmatrix
import FEM: AbstractMaterialStatus, AbstractElemStorage, AbstractFElement, FENode2
import FEM: Vertex4, Point2, GaussPoint2, Du, Dv, GeoTrig

export LinQuad

type LinQuadStorage <: AbstractElemStorage
    B::Matrix{Float64}
    DeBe::Matrix{Float64}
    Ke::Matrix{Float64}
    ɛ::Vector{Float64}
    f_int::Vector{Float64}
    u_field::Vector{Float64}
end

LinQuadStorage() = LinQuadStorage(zeros(4, 8), zeros(4,8),
                                  zeros(8,8), zeros(4), zeros(8), zeros(8))

type LinQuad{T <: AbstractMaterialStatus} <: AbstractFElement{T}
    vertices::Vertex4
    gps::Vector{GaussPoint2}
    n::Int
    interp::LinQuadInterp
    storage::LinQuadStorage
    matstats::Vector{T}
    temp_matstats::Vector{T}
end

# Constructors
function LinQuad{T <: AbstractMaterialStatus}(vertices::Vertex4, n, interp::LinQuadInterp,
                 storage::LinQuadStorage, gps::Vector{GaussPoint2}, matstat::T)
    matstats = T[]
    temp_matstats = T[]
    for i in 1:length(gps)
        push!(matstats, copy(matstat))
        push!(temp_matstats, copy(matstat))
    end
    LinQuad(vertices, gps, n, interp, storage, matstats, temp_matstats)
end

function LinQuadv(v::Vector{Int}, n, interp::LinQuadInterp,
                 storage::LinQuadStorage, gps::Vector{GaussPoint2})
    LinQuad(Vertex4(v[1], v[2], v[3], v[4]), gps, n, interp, storage)
end

get_ref_area(::LinQuad) = 4.0
get_ndofs(::LinQuad) = 8
get_geoelem(ele::LinQuad) = GeoQuad(ele.n, ele.vertices)
get_geotype(::LinQuad) = GeoQuad

createstorage(::Type{LinQuad}) = LinQuadStorage()
createinterp(::Type{LinQuad}) = LinQuadInterp()

function creategps(::Type{LinQuad})
    p = 1 / sqrt(3)
    [GaussPoint2(Point2(-p, -p), 1.0);
     GaussPoint2(Point2( p, -p), 1.0);
     GaussPoint2(Point2( p,  p), 1.0);
     GaussPoint2(Point2(-p,  p), 1.0)]
end

# Avoids allocation (remove?)
const LINTRIG_DOFTYPES = [Du, Dv]
function doftypes(::LinQuad, ::Int)
    return LINTRIG_DOFTYPES
end

function Bmatrix(elem::LinQuad, gp::GaussPoint2, nodes::Vector{FENode2})
    dNdx = dNdxmatrix(elem.interp, gp.local_coords, elem.vertices, nodes)
    B = elem.storage.B
    for i in 1:4
        B[1, 2*i - 1] = dNdx[i, 1]
        B[2, 2*i]     = dNdx[i, 2]

        B[4, 2*i - 1] = dNdx[i, 2]
        B[4, 2*i]     = dNdx[i, 1]
    end
    return B
end

end


text:
using HorizonSideRoots
#include("")
r = Robot()


text:
# Deterministic Finite Automaton
# ==============================

struct DFANode
    edges::Vector{Tuple{Edge,DFANode}}
    final::Bool
    eof_actions::ActionList
    nfanodes::Set{NFANode}
end

function DFANode(final::Bool, eof_actions::ActionList, nodes::Set{NFANode})
    return DFANode(Tuple{Edge,DFANode}[], final, eof_actions, nodes)
end

function DFANode(final::Bool, nodes::Set{NFANode})
    return DFANode(final, ActionList(), nodes)
end

function Base.show(io::IO, node::DFANode)
    print(io, summary(node), "(<", length(node.edges), " edges", '@', objectid(node), ">)")
end

struct DFA
    start::DFANode
end

# Check if the DFA is really deterministic or not.
function validate(dfa::DFA)
    for s in traverse(dfa.start)
        for i in 1:lastindex(s.edges), j in 1:i-1
            ei = s.edges[i][1]
            ej = s.edges[j][1]
            if overlaps(ei, ej)
                error("found non-deterministic edges")
            end
        end
    end
end

function nfa2dfa(nfa::NFA, unambiguous::Bool=true)
    newnodes = Dict{Set{NFANode},DFANode}()
    new(S) = get!(newnodes, S, DFANode(nfa.final ∈ S, S))
    isvisited(S) = haskey(newnodes, S)
    S = epsilon_closure(nfa.start)
    start = new(S)
    unvisited = [S]
    while !isempty(unvisited)
        # TODO: support fail
        S = pop!(unvisited)
        S_actions = accumulate_actions(S)
        s′ = new(S)

        if s′.final
            union!(s′.eof_actions, S_actions[nfa.final])
        end

        # accumulate edges and preconditions
        labels = Vector{ByteSet}()
        preconds = Set{Symbol}()
        for s in S, (e, t) in s.edges
            if !iseps(e)
                push!(labels, e.labels)
                union!(preconds, precondition_names(e.precond))
            end
        end

        # append DFA edges and nodes
        pn = collect(preconds)
        for label in disjoint_split(labels)
            # This enumeration will not finish in reasonable time when there
            # are too many preconditions.
            edges = Dict{Tuple{DFANode,ActionList},Vector{UInt64}}()
            for pv in UInt64(0):UInt64((1 << length(pn)) - 1)
                T = Set{NFANode}()
                actions = ActionList()
                for s in S, (e, t) in s.edges
                    if !isdisjoint(e.labels, label) && satisfies(e, pn, pv)
                        push!(T, t)
                        union!(actions, e.actions)
                        union!(actions, S_actions[s])
                    end
                end
                if !isempty(T)
                    T = epsilon_closure(T)
                    if !isvisited(T)
                        push!(unvisited, T)
                    end
                    push!(get!(edges, (new(T), actions), UInt64[]), pv)
                end
            end
            for ((t′, actions), pvs) in edges
                pn′, pvs′ = remove_redundant_preconds(pn, pvs)
                for pv′ in pvs′
                    push!(s′.edges, (Edge(label, make_precond(pn′, pv′), actions), t′))
                end
            end
        end
    end
    # Each key represents a set of NFANodes that collapses to one DFANode.
    # If any set contain conflicting possible actions, raise an error.
    if false
        for nfaset in keys(newnodes)
            validate_nfanodes(nfaset)
        end
    end
    return DFA(start)
end

function epsilon_closure(node::NFANode)
    return epsilon_closure(Set([node]))
end

function epsilon_closure(nodes::Set{NFANode})
    closure = Set{NFANode}()
    unvisited = copy(nodes)
    while !isempty(unvisited)
        s = pop!(unvisited)
        push!(closure, s)
        for (e, t) in s.edges
            if iseps(e) && t ∉ closure
                push!(unvisited, t)
            end
        end
    end
    return closure
end

"Find the nodes in this set where no other node in the set has an edge to it"
function gettop(S::Set{NFANode})
    top = copy(S)
    for s in S
        for (e, t) in s.edges
            if iseps(e)
                delete!(top, t)
            end
        end
    end
    return top
end

"Find paths from top nodes through epsilon edges, keeping track of actions taken."
function get_epsilon_paths(tops::Set{NFANode})
    paths = Tuple{Union{Nothing, Edge}, Vector{Symbol}}[]
    heads = [(node, Symbol[]) for node in tops]
    visited = Set{NFANode}()
    while !isempty(heads)
        node, actions = pop!(heads)
        if iszero(length(node.edges))
            push!(paths, (nothing, actions))
        end
        for (edge, child) in node.edges
            if iseps(edge)
                if !in(node, visited)
                    push!(heads, (child, append!(copy(actions), [a.name for a in edge.actions])))
                end
            else
                push!(paths, (edge, actions))
            end
        end
        push!(visited, node)
    end
    return paths
end

function validate_paths(paths::Vector{Tuple{Union{Nothing, Edge}, Vector{Symbol}}})
    all(actions == paths[1][2] for (n, actions) in paths) && return nothing
    for i in 1:length(paths) - 1
        edge1, actions1 = paths[i]
        for j in i+1:length(paths)
            edge2, actions2 = paths[j]
            # If either ends with EOF, they don't have same conditions and we can continue
            (edge1 === nothing) ⊻ (edge2 === nothing) && continue
            actions1 == actions2 && continue
            eof = (edge1 === nothing) & (edge2 === nothing)
            !(eof || overlaps(edge1, edge2)) && continue
            byte = eof ? "EOF" : repr(first(intersect(edge1.labels, edge2.labels)))
            a1 = isempty(actions1) ? "nothing" : actions1
            a2 = isempty(actions2) ? "nothing" : actions2
            error("Ambiguous DFA: Input $byte can lead to actions $a1 or $a2")
        end
    end
end

function validate_nfanodes(nodes::StableSet{NFANode})
    # First get "tops", that's the starting epsilon nodes that cannot be
    # reached by another epsilon node. All paths lead from those
    tops = gettop(nodes)

    # Quick path: If only one path, everything is OK
    length(tops) == 1 && all(length(node.edges) == 1 for node in nodes) && return nothing

    # Now we transverse all possible epsilon-paths and keep track of the actions
    # taken along the way
    paths = get_epsilon_paths(tops)

    # If any two paths have different actions, and can be reached with the same
    # byte, the DFA's actions cannot be resolved, and we raise an error.
    validate_paths(paths)
end

function disjoint_split(sets::Vector{ByteSet})
    # TODO: maybe too slow when length(sets) is large
    cut(s1, s2) = (intersect(s1, s2), setdiff(s1, s2))
    m = typemax(UInt64)
    disjsets = [ByteSet(m, m, m, m)]
    disjsets′ = ByteSet[]
    for x in sets
        for y in disjsets
            y1, y2 = cut(y, x)
            if !all(isdisjoint(z, y1) for z in sets)
                push!(disjsets′, y1)
            end
            if !all(isdisjoint(z, y2) for z in sets)
                push!(disjsets′, y2)
            end
        end
        disjsets, disjsets′ = disjsets′, disjsets
        empty!(disjsets′)
    end
    return disjsets
end

function accumulate_actions(S::Set{NFANode})
    top = gettop(S)
    @assert !isempty(top)
    actions = Dict(s => ActionList() for s in S)
    visited = Set{NFANode}()
    unvisited = top
    while !isempty(unvisited)
        s = pop!(unvisited)
        push!(visited, s)
        for (e, t) in s.edges
            if iseps(e)
                @assert !isconditioned(e.precond)
                union!(actions[t], e.actions)
                union!(actions[t], actions[s])
                if t ∉ visited
                    push!(unvisited, t)
                end
            end
        end
    end
    return actions
end

function satisfies(edge::Edge, names::Vector{Symbol}, pv::UInt64)
    for (n, v) in edge.precond
        i = findfirst(isequal(n), names)
        @assert i !== nothing
        @assert 0 < i ≤ 64
        vi = bitat(pv, i)
        if !(v == BOTH || (v == TRUE && vi) || (v == FALSE && !vi))
            return false
        end
    end
    return true
end

function remove_redundant_preconds(names::Vector{Symbol}, pvs::Vector{UInt64})
    mask(n) = ((1 << n) - 1) % UInt64
    newnames = Symbol[]
    pvs = copy(pvs)
    left = length(names)
    for name in names
        sort!(pvs)
        fnd = findfirst(pv -> bitat(pv, left), pvs)
        k = ifelse(fnd === nothing, 0, fnd) # TODO: See if there is a more elegant way of doing this.
        if (k - 1) * 2 == length(pvs)
            redundant = true
            for i in 1:k-1
                m = mask(left - 1)
                if pvs[i] & m != pvs[i+k-1] & m
                    redundant = false
                    break
                end
            end
        else
            redundant = false
        end
        if redundant
            left -= 1
            for i in 1:lastindex(pvs)
                # remove the redundant bit
                pvs[i] = pvs[i] & mask(left)
            end
        else
            push!(newnames, name)
            for i in 1:lastindex(pvs)
                # circular left shift
                pvs[i] = ((pvs[i] << 1) & mask(left)) | bitat(pvs[i], left)
            end
        end
    end
    return newnames, unique(pvs)
end

function make_precond(names::Vector{Symbol}, pv::UInt64)
    precond = Precondition()
    for (i, n) in enumerate(names)
        push!(precond, n => bitat(pv, i) ? TRUE : FALSE)
    end
    return precond
end

function bitat(x::UInt64, i::Integer)
    return ((x >> (i - 1)) & 1) == 1
end

function reduce_nodes(dfa::DFA)
    equivalents = get_equivalent(collect(traverse(dfa.start)))
    newnodes = Dict{Set{DFANode},DFANode}()
    new(S) = get!(newnodes, S) do
        s = first(S)
        return DFANode(s.final, s.eof_actions, foldl((x, s) -> union(x, s.nfanodes), S, init=Set{NFANode}()))
    end
    isvisited(T) = haskey(newnodes, T)
    S = equivalents[dfa.start]
    start = new(S)
    unvisited = [S]
    while !isempty(unvisited)
        S = pop!(unvisited)
        s′ = new(S)
        for (e, t) in first(S).edges
            T = equivalents[t]
            if !isvisited(T)
                push!(unvisited, T)
            end

            # Here, add the old edge to the new DFA. If an equivalent edge already
            # exists, instead add the labels of the old edge to the existing one.
            existing_edge = false
            for (i, (e´, t´)) in enumerate(s′.edges)
                if t´ == new(T) && e´.precond == e.precond && e´.actions == e.actions
                    newlabels = union(e.labels, e´.labels)
                    newedge = Edge(newlabels, e.precond, e.actions)
                    s′.edges[i] = (newedge, t´)
                    existing_edge = true
                    break
                end
            end
            existing_edge || push!(s′.edges, (e, new(T)))
        end
    end
    return DFA(start)
end

function get_groupof(v::Vector{DFANode})
    # These are sets of the outgoing bytes from every node
    labels = Dict(s => foldl((x, y) -> union(x, y[1].labels), s.edges, init=ByteSet()) for s in v)

    # First we create groups of nodes that MAY be identical based on quick-to-compute
    # characteristics. This narrows down the number of comparisons needed later.
    # Nodes may be equal if the have the same outgoing bytes, same final node and same
    # EOF actions.
    groupof = Dict{DFANode,Vector{DFANode}}()
    for s1 in v
        unique = true
        for group in values(groupof)
            s2 = first(group)
            if s1.final == s2.final && labels[s1] == labels[s2] && s1.eof_actions == s2.eof_actions
                push!(group, s1)
                groupof[s1] = group
                unique = false
                break
            end
        end
        unique && (groupof[s1] = [s1])
    end
    return groupof
end

function equivalent_pairs(groupof, v)
    # Here, we check each pair in each group. If they share an edge with an overlapping byte,
    # compatible preconditions, but the edge have different actions or leads no a non-equivalent
    # node, they're different. Because a parent's nonequivalence relies on the status of
    # its children, we need to update the parents every time we update a node.
    equivalent_pairs = Set{Tuple{DFANode,DFANode}}()
    for group in values(groupof), s1 in group, s2 in group
        push!(equivalent_pairs, (s1, s2))
    end

    # Get a node => vector of parents of node dict
    parentsof = Dict{DFANode,Vector{DFANode}}()
    for s in v, (e,c) in s.edges
        if haskey(parentsof, c)
            push!(parentsof[c], s)
        else
            parentsof[c] = [s]
        end
    end
    unupdated = Set(v)
    while !isempty(unupdated)
        s1 = pop!(unupdated)
        group = groupof[s1]
        for s2 in group
            (s1 == s2 || (s1, s2) ∉ equivalent_pairs) && continue
            for (e1, t1) in s1.edges, (e2, t2) in s2.edges
                if overlaps(e1, e2) && ((t1, t2) ∉ equivalent_pairs || e1.actions != e2.actions)
                    delete!(equivalent_pairs, (s1, s2))
                    haskey(parentsof, s1) && union!(unupdated, parentsof[s1])
                    break
                end
            end
        end
    end
    return equivalent_pairs
end

function split_group(group::Vector{DFANode}, pairs)
    remaining = copy(group)
    result = Vector{DFANode}[]
    nonequivalents = DFANode[]
    while !isempty(remaining)
        s1 = first(remaining)
        equivalents = DFANode[]
        for s2 in remaining
            (s1, s2) ∈ pairs ? push!(equivalents, s2) : push!(nonequivalents, s2)
        end
        push!(result, equivalents)
        (nonequivalents, remaining) = (remaining, nonequivalents)
        empty!(nonequivalents)
    end
    return result
end

"Creates a DFANode => Set{DFANode} dict with equivalent nodes for every node in v"
function get_equivalent(v::Vector{DFANode})
    groupof = get_groupof(v)
    pairs = equivalent_pairs(groupof, v)
    groups = collect(keys(IdDict{Vector{DFANode},Nothing}(v => nothing for v in values(groupof))))
    result = Dict{DFANode,Set{DFANode}}()
    for group in groups
        for subgroup in split_group(group, pairs)
            S = Set(subgroup)
            for node in subgroup
                result[node] = S
            end
        end
    end
    return result
end

function overlaps(e1::Edge, e2::Edge)
    return !(isdisjoint(e1.labels, e2.labels) || conflicts(e1.precond, e2.precond))
end

function revoke_finals(p::Function, dfa::DFA)
    newnodes = Dict{DFANode,DFANode}()
    new(s) = get!(newnodes, s) do
        return DFANode(s.final && !p(s), s.eof_actions, s.nfanodes)
    end
    for s in traverse(dfa.start)
        s′ = new(s)
        for (e, t) in s.edges
            push!(s′.edges, (e, new(t)))
        end
    end
    return DFA(new(dfa.start))
end

function dfa2nfa(dfa::DFA)
    newnodes = Dict{DFANode,NFANode}()
    new(s) = get!(newnodes, s, NFANode())
    final = NFANode()
    for s in traverse(dfa.start)
        s′ = new(s)
        for (e, t) in s.edges
            push!(s′.edges, (e, new(t)))
        end
        if s.final
            push!(s′.edges, (Edge(eps, s.eof_actions), final))
        end
    end
    start = NFANode()
    push!(start.edges, (Edge(eps), new(dfa.start)))
    return NFA(start, final)
end

function remove_dead_nodes(dfa::DFA)
    backrefs = Dict(dfa.start => Set{DFANode}())
    for s in traverse(dfa.start), (_, t) in s.edges
        push!(get!(backrefs, t, Set{DFANode}()), s)
    end

    alive = Set{DFANode}()
    unvisited = [s for s in keys(backrefs) if s.final]
    while !isempty(unvisited)
        s = pop!(unvisited)
        push!(alive, s)
        for t in backrefs[s]
            if t ∉ alive
                push!(unvisited, t)
            end
        end
    end
    @assert dfa.start ∈ alive

    newnodes = Dict{DFANode,DFANode}()
    new(s) = get!(newnodes, s, DFANode(s.final, s.eof_actions, s.nfanodes))
    isvisited(s) = haskey(newnodes, s)
    unvisited = [dfa.start]
    while !isempty(unvisited)
        s = pop!(unvisited)
        s′ = new(s)
        for (e, t) in s.edges
            if t ∈ alive
                if !isvisited(t)
                    push!(unvisited, t)
                end
                push!(s′.edges, (e, new(t)))
            end
        end
    end

    return DFA(new(dfa.start))
end


text:
# helper module with misc DArray utilities
module joDAutils
    using DistributedArrays
    using JOLI: jo_convert, joAbstractLinearOperator

    function DArray5(init, dims, procs, idxs, cuts)
        dist = map(i->length(i)-1,cuts)
        np = prod(dist)
        procs = reshape(procs[1:np], ntuple(i->dist[i], length(dist)))
        id = DistributedArrays.next_did()
        return DArray(id, init, dims, procs, idxs, cuts)
    end
    function jo_x_mv!(A::joAbstractLinearOperator,in::DArray{ADDT,2},out::DArray{ARDT,2}) where {ADDT,ARDT}
        out[:L]=jo_convert(ARDT,A*in[:L])
        return nothing
    end
    function jo_x_mv!(F::Function,in::DArray{ADDT,2},out::DArray{ARDT,2}) where {ADDT,ARDT}
        out[:L]=jo_convert(ARDT,F(in[:L]))
        return nothing
    end
end
using .joDAutils

export dparts
"""
    julia> dparts(da)

return partitioning vector of DArray if partioned in single dimension

# Signature

    dparts(da::DArray{T,N})

# Arguments

- `da`: DArray

# Notes

- if DArray is quasi-distributed (over single worker), dparts returns size(da,N)

"""
function dparts(da::DArray{T,N}) where {T,N}
    chunks=map(i->length(i)-1,da.cuts)
    dim=findfirst(i->i>1,chunks); dim = dim==nothing ? N : dim
    ldim=findlast(i->i>1,chunks); ldim = ldim==nothing ? N : ldim
    dim==ldim || throw(joPAsetupException("joPAsetup: cannot return parts of a DArray partitioned in multiple dimensions"))
    parts=map(i->length(i[dim]),da.indices)
    return vec(parts)
end

export dcopy
"""
    julia> dcopy(dtr,[dst])

copy transpose(DArray) into a new DArray using predefined joPAsetup

# Signature

    dcopy(Dtr::Transpose{T,<:DArray{T,2}},dst::joPAsetup)
    dcopy(Dtr::Transpose{T,<:DArray{T,2}})

# Arguments

- `dtr`: transpose(DArray)
- `dst`: target joPAsetup

"""
function dcopy(Dtr::Transpose{T,<:DArray{T,2}},dst::joPAsetup) where T
    Dst=joPAsetup(parent(Dtr))
    Dst.dims==reverse(dst.dims) || throw(joPAsetupException("the sizes of original array and provided target distributor do not match"))
    Dst.procs==dst.procs || throw(joPAsetupException("the workers of original array and provided target distributor do not match"))

    D = parent(Dtr)
    joDAutils.DArray5(dst.dims, dst.procs, dst.idxs, dst.cuts) do I
        # @debug I rI=reverse(I)
        lp = Array{T}(undef, map(length, I))
        rp = convert(Array, D[reverse(I)...])
        transpose!(lp, rp)
    end
end
function dcopy(Dtr::Transpose{T,<:DArray{T,2}}) where T
    Dst=transpose(joPAsetup(parent(Dtr)))
    return dcopy(Dtr,Dst)
end

export dalloc
"""
    julia> dalloc(dims, [...])

Allocates a DistributedArrays.DArray without value assigment.

Use it to allocate quicker the array that will have all elements overwritten.

# Signature

    dalloc(dims::Dims, [...])

# Arguments

- optional trailing arguments are the same as those accepted by `DArray`.

"""
dalloc(dims::Dims, args...) = DArray(I->Array{joFloat}(map(length,I)), dims, args...)
dalloc(::Type{T}, dims::Dims, args...) where {T} = DArray(I->Array{T}(undef,map(length,I)), dims, args...)
dalloc(::Type{T}, d1::Integer, drest::Integer...) where {T} = dalloc(T, convert(Dims, tuple(d1, drest...)))
dalloc(d1::Integer, drest::Integer...) = dalloc(joFloat, convert(Dims, tuple(d1, drest...)))
dalloc(d::Dims) = dalloc(joFloat, d)

export dfill
"""
    julia> dfill(F, d; [DT])

Constructs a DistributedArrays.DArray, according to given distributor, filled with elements provided by anonymous function F.

# Signature

    dfill(F::Function,d::joPAsetup;DT::DataType=d.DT)

# Arguments

- `F`: anonymous function of the form `I->f(...,map(length,I)))`
- `d`: see help for joPAsetup
- `DT`: keyword argument to overwrite the type in joPAsetup

# Notes

- function F will be passed via `map(length,I)` the tuple with dimensions of local part
- one has to pass array type manualy to F

# Examples

- `dfill(I->ones(d.DT,map(length,I)),d)`: fill a distributed array with ones of type d.DT

"""
function dfill(F::Function,d::joPAsetup;DT::DataType=d.DT)
    id=DistributedArrays.next_did()
    procs = reshape(d.procs, ntuple(i->d.chunks[i], length(d.chunks)))
    return DArray(id, F, d.dims, procs, d.idxs, d.cuts)
end

"""
    julia> dfill(x, d; [DT])

Constructs a DistributedArrays.DArray, according to given distributor, filled with x.

# Signature

    dfill(x::Number,d::joPAsetup;DT::DataType=d.DT)

# Arguments

- `d`: see help for joPAsetup
- `DT`: keyword argument to overwrite the type in joPAsetup

# Examples

- `dfill(3.,d)`: fill a distributed array with d.DT(3.)
- `dfill(3.,d;DT=Float32)`: fill a distributed array with Float32(3.)

"""
function dfill(x::Number,d::joPAsetup;DT::DataType=d.DT)
    id=DistributedArrays.next_did()
    X=DT(x)
    init=I->fill(X,map(length,I))
    procs = reshape(d.procs, ntuple(i->d.chunks[i], length(d.chunks)))
    return DArray(id, init, d.dims, procs, d.idxs, d.cuts)
end

"""
    julia> dalloc(d; [DT])

Allocates a DistributedArrays.DArray, according to given distributor, without value assigment.

Use it to allocate quicker the array that will have all elements overwritten.

# Signature

    dalloc(d::joPAsetup;DT::DataType=d.DT)

# Arguments

- `d`: see help for joPAsetup
- `DT`: keyword argument to overwrite the type in joPAsetup

# Examples

- `dalloc(d)`: allocate an array
- `dalloc(d,DT=Float32)`: allocate array and overwite d.DT with Float32

"""
function dalloc(d::joPAsetup;DT::DataType=d.DT)
    id=DistributedArrays.next_did()
    init=I->Array{DT}(undef,map(length,I))
    procs = reshape(d.procs, ntuple(i->d.chunks[i], length(d.chunks)))
    return DArray(id, init, d.dims, procs, d.idxs, d.cuts)
end
#DArray(d::joPAsetup;DT::DataType=d.DT) = dalloc(d;DT=DT)

export dzeros
"""
    julia> dzeros(d; [DT])

Constructs a DistributedArrays.DArray, according to given distributor, filled with zeros.

# Signature

    dzeros(d::joPAsetup;DT::DataType=d.DT)

# Arguments

- `d`: see help for joPAsetup
- `DT`: keyword argument to overwrite the type in joPAsetup

# Examples

- `dzeros(d)`: allocate an array of zeros
- `dzeros(d,DT=Float32)`: allocate array of Float32 zeros

"""
function dzeros(d::joPAsetup;DT::DataType=d.DT)
    id=DistributedArrays.next_did()
    init=I->zeros(DT,map(length,I))
    procs = reshape(d.procs, ntuple(i->d.chunks[i], length(d.chunks)))
    return DArray(id, init, d.dims, procs, d.idxs, d.cuts)
end

export dones
"""
    julia> dones(d; [DT])

Constructs a DistributedArrays.DArray, according to given distributor, filled with ones.

# Signature

    dones(d::joPAsetup;DT::DataType=d.DT)

# Arguments

- `d`: see help for joPAsetup
- `DT`: keyword argument to overwrite the type in joPAsetup

# Examples

- `dones(d)`: allocate an array of ones
- `dones(d,DT=Float32)`: allocate array of Float32 ones

"""
function dones(d::joPAsetup;DT::DataType=d.DT)
    id=DistributedArrays.next_did()
    init=I->ones(DT,map(length,I))
    procs = reshape(d.procs, ntuple(i->d.chunks[i], length(d.chunks)))
    return DArray(id, init, d.dims, procs, d.idxs, d.cuts)
end

export drand
"""
    julia> drand(d; [DT], [RNG])

Constructs a DistributedArrays.DArray, according to given distributor, filled using built-in rand.

# Signature

    drand(d::joPAsetup;DT::DataType=d.DT,RNG::AbstractRNG=RandomDevice())

# Arguments

- `d`: see help for joPAsetup
- `DT`: keyword argument to overwrite the type in joPAsetup
- `RNG`: random-number generator function (see help for rand/randn)

# Examples

- `drand(d)`: allocate an array with rand
- `drand(d,DT=Float32)`: allocate array with rand of Float32
- `drand(d,DT=Float32,RNG=MersenneTwister(1234))`: allocate array with rand of Float32 using MersenneTwister() random device

"""
function drand(d::joPAsetup;DT::DataType=d.DT,RNG::AbstractRNG=RandomDevice())
    id=DistributedArrays.next_did()
    init=I->rand(RNG,DT,map(length,I))
    procs = reshape(d.procs, ntuple(i->d.chunks[i], length(d.chunks)))
    return DArray(id, init, d.dims, procs, d.idxs, d.cuts)
end

export drandn
"""
    julia> drandn(d; [DT], [RNG])

Constructs a DistributedArrays.DArray, according to given distributor, filled using built-in randn.

# Signature

    drandn(d::joPAsetup;DT::DataType=d.DT,RNG::AbstractRNG=RandomDevice())

# Arguments

- `d`: see help for joPAsetup
- `DT`: keyword argument to overwrite the type in joPAsetup
- `RNG`: random-number generator function (see help for rand/randn)

# Notes

- only float type are supported by randn (see help for randn)

# Examples

- `drandn(d)`: allocate an array with randn
- `drandn(d,DT=Float32)`: allocate array with randn of Float32
- `drandn(d,DT=Float32,RNG=MersenneTwister(1234))`: allocate array with randn of Float32 using MersenneTwister() random device

"""
function drandn(d::joPAsetup;DT::DataType=d.DT,RNG::AbstractRNG=RandomDevice())
    DT<:Integer && @warn "Cannot use Integer type in randn.\n\t Overwite joPAsetup's type using DT keyword\n\t or create Float joPAsetup.\n\t Falling back to joFloat!" key="JOLI:drandn:Integer" maxlog=1
    DT= (DT<:Integer) ? joFloat : DT
    id=DistributedArrays.next_did()
    init=I->randn(RNG,DT,map(length,I))
    procs = reshape(d.procs, ntuple(i->d.chunks[i], length(d.chunks)))
    return DArray(id, init, d.dims, procs, d.idxs, d.cuts)
end

export distribute
"""
    julia> distribute(A,d)

Distributes DArray according to given joPAsetup.

# Signature

    distribute(a::AbstractArray,d::joPAsetup)

# Arguments

- `A`: array to ditribute
- `d`: see help for joPAsetup

# Notes

- the type in joPAsetup is ignored here
- distributes over last non-singleton (worker-wise) dimension
- one of the dimensions must be large enough to hold at least one element on each worker

# Examples

- `distribute(A,d)`: distribute A using given distributor
- `distribute(A,joPAsetup(size(A)...))`: distribute A using default distributor settings

"""
function distribute(A::AbstractArray,d::joPAsetup)
    size(A)==d.dims || throw(joPAsetupException("joPAsetup: array size does not match dims of joPAsetup"))
    id=DistributedArrays.next_did()
    s = DistributedArrays.verified_destination_serializer(reshape(d.procs, size(d.idxs)), size(d.idxs)) do pididx
        A[d.idxs[pididx]...]
    end
    init = I->DistributedArrays.localpart(s)
    procs = reshape(d.procs, ntuple(i->d.chunks[i], length(d.chunks)))
    return DArray(id, init, d.dims, procs, d.idxs, d.cuts)
end


text:
g = Extrusion(5)
i = ConstantInflow(0)
o = Orifice(C = 0.7, A = 0.005, z = 0)
model = ReservoirModel(geometry = g, inflow = i, outlet = o, z₀ = 10.0)
sim = Simulation(model = model, Δt = 5.0, n = 500)
t, z, Qout = run(sim)

struct ContinuousCurve <: ReservoirGeometry end

surface_area(g::ContinuousCurve, z) = 100.0 * z^2

struct ContinuousHydrograph <: Inflow end

function discharge(i::ContinuousHydrograph, t)
    if t < 40_000.0
        2000.0 * (1.0 - t / 40_000.0)
    else
        0.0
    end
end

g = ContinuousCurve()
i = ContinuousHydrograph()
o = Weir(C = 1.0, B = 100.0, z = 95.0)
model = ReservoirModel(geometry = g, inflow = i, outlet = o, z₀ = 50.0)
sim = Simulation(model = model, Δt = 100.0, n = 1000)
t, z, Qout = run(sim)


text:
module EXT_TEST
using Compat
export EXT_TEST_hmmalign
export EXT_TEST_hmmsearch
export EXT_TEST_cmsearch

######################################################################
# FUNCTION:              EXT_TEST_hmmalign
######################################################################
function EXT_TEST_hmmalign()
	try
@compat		run(pipeline(`hmmalign -h`, DevNull))
	catch
		return false
	end
	return true
end

######################################################################
# FUNCTION:              EXT_TEST_hmmsearch
######################################################################
function EXT_TEST_hmmsearch()
	try
@compat		run(pipeline(`hmmsearch -h`, DevNull))
	catch
		return false
	end
	return true
end

######################################################################
# FUNCTION:              EXT_TEST_cmsearch
######################################################################
function EXT_TEST_cmsearch()
	try
@compat		run(pipeline(`cmsearch -h`, DevNull))
	catch
		return false
	end
	return true
end

#/module
end


text:
module Pade

using ..MittagLefflerFunctions: rΓ
using OffsetArrays
using LinearAlgebra: svd, lu

export MLPade

struct MLPade{T<:AbstractFloat}
    α :: T
    β :: T
    m :: Int64
    n :: Int64
    p :: Vector{T}
    q :: Vector{T}
end

function MLPade(α::T, β::T, m::Integer, n::Integer,
                method::Symbol) where T <: AbstractFloat
    A = pade_matrix(α, β, m, n)
    if method == :fix
        p, q = pade_fix(A)
    elseif method == :svd
        p, q, y = pade_svd(A)
    elseif method == :lu
        p, q = pade_lu(A)
    else
        error("Unknown error $method")
    end
    return MLPade(α, β, m, n, p, q)
end

function (E::MLPade{T})(x::T) where T <: AbstractFloat
    p, q = E.p, E.q
    if x < 0
        throw(DomainError(x, "argument must be greater than or equal to zero"))
    elseif x < 1
        R = evalpoly(x, p) / evalpoly(x, q)
    else
        rx = 1 / x
        R = rx * evalpoly(rx, p[end:-1:1]) / evalpoly(rx, q[end:-1:1])
    end
    return R
end

function pade_matrix(α::T, β::T,
                     m::Integer, n::Integer) where T <: AbstractFloat
    if mod(m+n,2) == 0
        throw(ArgumentError("m+n must be odd"))
    end
    r = div(m+n-1,2)
    a = OffsetArray{T}(undef, 0:m-1)
    b = OffsetArray{T}(undef, 0:n-1)
    pow = one(T)
    for k = 0:m-1
        a[k] = pow * rΓ(β+k*α)
        pow = -pow
    end
    b[0] = zero(T)
    pow = -one(T)
    for k = 1:n-1
        pow = -pow
        b[k] = pow * rΓ(β-k*α)
    end
    A = zeros(T, 2r, 2r+1)
    if m ≥ r+1
        row = 0
        for k = 0:r-1
            row += 1
            A[row,k+1] = one(T)
            for j = 0:k
                A[row,r+1+j] = -a[k-j]
            end
        end
        for k = r:m-1
            row += 1
            for j = 0:r
                A[row,r+1+j] = -a[k-j]
            end
        end
        for k = r-n+1:r-1
            row += 1
            A[row,k+1] = one(T)
            for j = k+1:r
                A[row,r+1+j] = -b[j-k]
            end
        end
    else
        row = 0
        for k = 0:m-1
            row += 1
            A[row,k+1] = one(T)
            for j = 0:k
                A[row,r+1+j] = -a[k-j]
            end
        end
        for k = -(r-m):-1
            row += 1
            for j = 0:r
                A[row,r+1+j] = -b[j-k]
            end
        end
        for k = 0:r-1
            row += 1
            A[row,k+1] = one(T)
            for j = k+1:r
                A[row,r+1+j] = -b[j-k]
            end
        end
    end
    return A
end

function pade_svd(α::T, β::T, m::Integer, n::Integer) where T <: AbstractFloat
    A = pade_matrix(α, β, m, n)
    p, q = pade_svd(A)
    return p, q
end

function pade_svd(A::Matrix{T}) where T <: AbstractFloat
    r = div(size(A,1), 2)
    F = svd(A, full=true)
    p = F.V[1:r,2r+1]
    q = F.V[r+1:2r+1,2r+1]
    q1 = q[1]
    p ./= q1
    q ./= q1
    return p, q, F.V[:,2r+1]
end

function pade_lu(A::Matrix{T}) where T <: AbstractFloat
    r = div(size(A,1), 2)
    F = lu(A)
    y = zeros(T, 2r+1)
    y[2r+1] = one(T)
    for i = 2r:-1:1
        s = zero(T)
        for j = i+1:2r+1
            s += F.U[i,j] * y[j]
        end
        y[i] = -s / F.U[i,i]
    end
    p = y[1:r]
    q = y[r+1:2r+1]
    μ = q[1]
    for k = 1:r
        p[k] /= μ
    end
    for k = 1:r+1
        q[k] /= μ
    end
    return p, q
end

function pade_fix(A::Matrix{T}) where T <: AbstractFloat
    r = div(size(A,1), 2)
    Atilde = [ A[:,1:r] A[:,r+2:2r+1] ]
    b = -A[:,r+1]
    y = Atilde \ b
    p = y[1:r]
    q = Vector{T}(undef, r+1)
    q[1] = one(T)
    q[2:r+1] = y[r+1:2r]
    return p, q
end

end # module


text:
module FluxFederated

# exceptions
include("exceptions.jl")
export ConnectionError

# crypt_utils
include("crypt_utils.jl")

export gen_salt, gen_key32, to_str, to_bytes, Bytes
export RSAKey, PublicRSAKey, PrivateRSAKey
export gen_RSAKeys, gen_prime, pack_int, unpack_int
export secure_rand, secure_randstring

# sock_utils
include("sock_utils.jl")

export handshake_client, handshake_server, secure_send, secure_receive
export secure_send_bytes, secure_send_int, secure_receive_bytes, secure_receive_int
export var_to_bytes, var_from_bytes, model_to_bytes, model_from_buffer, model_from_bytes

# serve
include("async_serve.jl")
export serve_model

end # module


text:
using Dates
start_time = Dates.now()
println("Starting imports: ", start_time)
# using StatsBase
# using Distributions
# using Gadfly
# using TimerOutputs
using Random
#import Cairo, Fontconfig
using JSON
end_time = Dates.now()
println("Finished imports: ", end_time - start_time)
include("./main.jl")

# Function to plot graph
# function plotgraph(xaxis::Array{Int64,1},yaxis::Array{Float64,1},x_label::String,y_label::String,heading::String)
#     e = plot(x=xaxis, y=yaxis, Guide.xlabel(x_label), Guide.ylabel(y_label), Guide.title(heading), Geom.line,
#         Theme(minor_label_font_size=14pt, major_label_font_size=14pt));
#     return(e)
# end

function write_results(filename, x_vals, scores, mem, time)
    f = open(filename, write=true)
    d = Dict(
        "x_vals" => x_vals,
        "scores" => scores,
        "mems" => mem,
        "times" => time
    )
    JSON.write(f, json(d))
end

#function for testing with increasing number of input sequences.
function msatest1()
    println("TEST msatest1")
    score_mat = Float64[]
    num_of_inp_seq = Int64[]
    bytes = Float64[]
    time = Float64[]

    #Vary the number of input sequences and record reading.
    for t in 10:20:500
        println("t=", t)
        sequences = generate_sequences(t, 50)
        MSA = @timed PSO_MSA(sequences,20)
        push!(score_mat,MSA[1])
        push!(num_of_inp_seq,t)
        push!(bytes, MSA[3])
        push!(time, MSA[2] - MSA[4])
    end

    #First reading considered as a dry run.
    popfirst!(score_mat)
    popfirst!(num_of_inp_seq)
    popfirst!(bytes)
    popfirst!(time)

    println(num_of_inp_seq)
    println(score_mat)
    println(bytes)

    write_results("msatest1.json", num_of_inp_seq, score_mat, bytes, time)

    #z1 = plotgraph(num_of_inp_seq, score_mat, "Number of input sequences", "Score", "Plot of Score v/s Number of input sequences")
    #z2 = plotgraph(num_of_inp_seq, bytes, "Number of input sequences", "Memory in bytes", "Plot of Memory v/s Number of input sequences")
    #z3 = plotgraph(num_of_inp_seq, time, "Number of input sequences", "Time in seconds", "Plot of Time v/s Number of input sequences")
    #display(z1)
    #display(z2)
    #display(z3)

    #draw(PNG("msatest1-z1.png", 8inch, 8inch), z1)
end

#function for testing with increasing length of input sequences.
function msatest2()
    println("TEST msatest2")
    score_mat = Float64[]
    inp_seq_len = Int64[]
    bytes = Float64[]
    time = Float64[]

    #Vary the length of input sequences and record reading.
    for N in 20:20:800
        println("N=", N)
        sequences = generate_sequences(10, N)
        MSA = @timed PSO_MSA(sequences, 20)
        push!(score_mat, MSA[1])
        push!(inp_seq_len, N)
        push!(bytes, MSA[3])
        push!(time, MSA[2]-MSA[4])
    end

    #First reading considered as a dry run.
    popfirst!(score_mat)
    popfirst!(inp_seq_len)
    popfirst!(bytes)
    popfirst!(time)

    write_results("msatest2.json", inp_seq_len, score_mat, bytes, time)

    # z1 = plotgraph(inp_seq_len, score_mat, "Input sequence length", "Score", "Plot of Score v/s Length of input sequence")
    # z2 = plotgraph(inp_seq_len, bytes, "Input sequence length", "Memory in bytes", "Plot of Memory v/s Length of input sequence")
    # z3 = plotgraph(inp_seq_len, time, "Input sequence length", "Time in seconds", "Plot of Time v/s Length of input sequence")
    # display(z1)
    # display(z2)
    # display(z3)
end

#function for testing with increasing number of iterations.
function msatest3()
    println("TEST msatest3")
    score_mat = Float64[]
    itr = Int64[]
    bytes = Float64[]
    time = Float64[]
    sequences = generate_sequences(10, 50)

    #Vary the number of iterations and record reading.
    for i in 0:5:120
        println("i=", i)
        Random.seed!(1)
        MSA = @timed PSO_MSA(sequences,i)
        push!(score_mat,MSA[1])
        push!(itr,i)
        push!(bytes, MSA[3])
        push!(time, MSA[2]-MSA[4])
    end

    #First reading considered as a dry run.
    popfirst!(score_mat)
    popfirst!(itr)
    popfirst!(bytes)
    popfirst!(time)

    write_results("msatest3.json", itr, score_mat, bytes, time)

    # z1 = plotgraph(itr, score_mat, "Number of iterations", "Score", "Plot of Score v/s Number of iterations")
    # z2 = plotgraph(itr, bytes, "Number of iterations", "Memory in bytes", "Plot of Memory v/s Number of iterations")
    # z3 = plotgraph(itr, time, "Number of iterations", "Time in seconds", "Plot of Time v/s Number of iterations")
    # display(z1)
    # display(z2)
    # display(z3)
end

#msatest1()
#msatest2()
msatest3()


text:
# Is this usable data?
function isconcrete{T<:Number}(x::T)
    !isna(x) && isfinite(x)
end

function isconcrete(x)
    !isna(x)
end

# How many concrete elements in an iterable
function concrete_length(xs)
    n = 0
    for x in xs
        if isconcrete(x)
            n += 1
        end
    end
    n
end

function nonzero_length(xs)
    n = 0
    for x in xs
        if x != 0
            n += 1
        end
    end
    n
end

# Create a new object of type T from a with missing values (i.e., those set to
# nothing) inherited from b.
function inherit{T}(a::T, b::T)
    c = copy(a)
    inherit!(c, b)
    c
end

function inherit!{T}(a::T, b::T)
    for field in T.names
        val = getfield(a, field)
        # TODO: this is a hack to let non-default labelers overide the defaults
        if val === nothing || val === string || val === fmt_float
            setfield(a, field, getfield(b, field))
        end
    end
    nothing
end

isnothing(u) = is(u, nothing)
issomething(u) = !isnothing(u)

negate(f) = x -> !f(x)

function has{T,N}(xs::AbstractArray{T,N}, y::T)
    for x in xs
        if x == y
            return true
        end
    end
    return false
end

function push!{T}(xs::Vector{T}, ys::T...)
    for y in ys
        push!(xs, y)
    end
end

Maybe(T) = Union(T, Nothing)

# Float64 -> String, trimming trailing zeros when appropriate.
# This is largely taken from cairo's function _cairo_dtostr.
function fmt_float(x::Float64)
    if x < 0.1
        a = @sprintf("%0.4f", x)
    else
        a = @sprintf("%f", x)
    end

    n = length(a)
    while a[n] == '0'
        n -= 1
    end

    if a[n] == '.'
        n -= 1
    end

    a[1:n]
end

function lerp(x::Float64, a::Float64, b::Float64)
    a + (b - a) * max(min(x, 1.0), 0.0)
end

# Generate a unique id, primarily for assigning IDs to SVG elements.
let next_unique_svg_id_num = 0
    global unique_svg_id
    function unique_svg_id()
        uid = @sprintf("id%d", next_unique_svg_id_num)
        next_unique_svg_id_num += 1
        uid
    end
end

# Remove any markup or whitespace from a string.
function escape_id(s::String)
    s = replace(s, r"<[^>]*>", "")
    s = replace(s, r" ", "_")
    s
end


text:
# Copyright 2017,2018 Eric S. Tellez
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using Test
using DelimitedFiles

function loadiris()
    url = "https://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data"
    filename = basename(url)
    if !isfile(filename)
        download(url, filename)
    end

    data = readdlm(filename, ',')
    X = data[:, 1:4]
    X = [Float64.(X[i, :]) for i in 1:size(X, 1)]
    y = String.(data[:, 5])
    X, y
end

function loadlinearreg()
    X = [rand(2) .+ i for i in 1:100]
    y = range(1, stop=100, length=100) .+ rand(100)
    X, y
end


text:
module Pinguu

using CSV
using DataFrames
using Dates
using Format
using Glob
using HTTP
using Indicators
using JSON
using Telegram, Telegram.API
using PyCall
using Query
using StatsBase
using Suppressor
using TypedTables
binance = pyimport("binance.client");
export binance

include("api/check_orders.jl")
include("api/check_why_there_is_no_sell.jl")
include("api/determine_closed.jl")
include("api/examine_order.jl")
include("api/place_order_command.jl")
include("api/trade.jl")
include("api/trade_helper.jl")
include("api/buy_helper.jl")
include("api/BNB_control.jl")
include("api/get_current_ticker.jl")

include("analysis/check_balances.jl")
include("analysis/convert.jl")
include("analysis/ema.jl")
include("analysis/klines.jl")
include("analysis/trigger.jl")
include("analysis/calc_order_fee.jl")
include("analysis/create_tax_report.jl")

include("interface/history.jl")
include("interface/notification.jl")

include("invest/cancel_invest.jl")
include("invest/check_liquid.jl")
include("invest/move_order_to_invest.jl")
include("invest/move_to_filled.jl")
include("invest/place_invest.jl")

include("tasks/scanner.jl")

include("stats/stats_fct.jl")
include("grid_trading/grid_trading.jl")

# Analysis export
export convert_klines, get_ema
export get_trigger_limits, place_triggers

# API export
export buy_coin_and_place_sell, buy_coin_and_place_sell_test, check_open_orders
export place_first_order!, check_open_sell!, check_open_buy!
export place_first_poly_order!, check_open_poly_sell!, check_open_poly_buy!

const base_dir = ENV["pinguu"]
const telegram = ENV["telegram"]
export base_dir
export telegram

end


text:
# This file was generated, do not modify it. # hide
linreg = @load LinearRegressor pkg="GLM"

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials