text:
using Dates

# Various functions to get the JSONtext of the data
function getdatapart(df, dffields, plottype::Symbol)

    # start the main string
    sb = StringBuilder()
    append!(sb, "[");

    # get the specific string for this type of data
    datasb = StringBuilder()
    numrows = length(df[!, 1]);

    # loop over rows and create a dictionary for each row
    if plottype == :multitrumpet #4D with 3 indices
        cols = (df[!, 1], df[!, 2], df[!, 3], df[!, 4])
        datastring = getdatapart_4d(cols, dffields, numrows, datasb)
    elseif plottype == :multiline || plottype == :trumpet #3D with 2 indices, one of which is time
        cols = (df[!, 1], df[!, 2], df[!, 3])
        datastring = getdatapart_3d_time(cols, dffields, numrows, datasb)
    elseif plottype == :multihistogram || plottype == :multibar || plottype == :multiscatter #3D with 2 indices, none of which is time
        cols = (df[!, 1], df[!, 2], df[!, 3])
        datastring = getdatapart_3d(cols, dffields, numrows, datasb)
    elseif plottype == :line  #2D with 1 index, one of which is time
        cols = (df[!, 1], df[!, 2])
        datastring = getdatapart_2d_time(cols, dffields, numrows, datasb)
    else # :bar and :histogram
        cols = (df[!, 1], df[!, 2])
        datastring = getdatapart_2d(cols, dffields, numrows, datasb)
    end

    append!(sb, datastring * "]");
    datapart = String(sb)

    return JSON.JSONText(datapart)
end

function getdatapart_4d(cols, dffields, numrows, datasb)
    for i = 1:numrows

        append!(datasb, "{
            \"" * dffields[1]  * "\":\"" * string(Date(cols[1][i])) * "\",
            \"" * dffields[2] * "\":\"" * string(cols[2][i]) * "\",
            \"" * dffields[3] * "\":\"" * string(cols[3][i]) * "\",
            \"" * dffields[4] * "\":\"" * string(cols[4][i]) * "\"}")

        if i != numrows
            append!(datasb, ",")
        end
    end
    return String(datasb)
end

function getdatapart_3d_time(cols, dffields, numrows, datasb)
    for i = 1:numrows

        append!(datasb, "{\"" * dffields[1]  * "\":\"" * string(Date(cols[1][i]))
            * "\",\"" * dffields[2] * "\":\"" * string(cols[2][i]) * "\",\""
            * dffields[3] * "\":\"" * string(cols[3][i]) * "\"}")

        if i != numrows
            append!(datasb, ",")
        end
    end
    return String(datasb)
end

function getdatapart_3d(cols, dffields, numrows, datasb)
    for i = 1:numrows

        append!(datasb, "{\"" * dffields[1]  * "\":\"" * string(cols[1][i])
            * "\",\"" * dffields[2] * "\":\"" * string(cols[2][i]) * "\",\""
            * dffields[3] * "\":\"" * string(cols[3][i]) * "\"}")

        if i != numrows
            append!(datasb, ",")
        end
    end
    return String(datasb)
end

function getdatapart_2d_time(cols, dffields, numrows, datasb)
    for i = 1:numrows
        append!(datasb, "{\"" * dffields[1]  * "\":\"" * string(Date(cols[1][i]))
            * "\",\"" * dffields[2] * "\":\"" * string(cols[2][i]) * "\"}")

        if i != numrows
            append!(datasb, ",")
        end
    end

    return String(datasb)
end

function getdatapart_2d(cols, dffields, numrows, datasb)
    for i = 1:numrows

        append!(datasb, "{\"" * dffields[1] * "\":\"" * string(cols[1][i]) *
            "\",\"" * dffields[2] * "\":\"" * string(cols[2][i]) * "\"}") #end of dictionary

        if i != numrows
            append!(datasb, ",")
        end
    end
    return String(datasb)
end


text:
## --- Parse a delimited string

    """
    ```julia
    delim_string_parse!(result, str, delim, [T];
        \toffset::Integer=0,
        \tmerge::Bool=false,
        \tundefval=NaN)
    ```
    Parse a delimited string `str` with delimiter `delim` into values of type `T`
    and return the answers in a pre-allocated `result` array provided as input.

    If `T` is not specified, it the `eltype` of the `result` array will be used by default.
    """
    function delim_string_parse!(result::Array, str::AbstractString, delim::Char, T::Type=eltype(result); offset::Integer=0, merge::Bool=false, undefval=NaN)

        # Make sure the output data type allows our chosen value for undefined data
        undefval = convert(T, undefval)

        # Ignore initial delimiter
        last_delim_pos = 0
        if ~isempty(str) && str[1] == delim
            last_delim_pos = 1
        end

        # Cycle through string parsing text betweeen delims
        delim_pos = 0
        n = offset
        if merge
            for i = 1:length(str)
                if str[i] == delim
                    delim_pos = i
                    if delim_pos > last_delim_pos+1
                        n += 1
                        parsed = nothing
                        if delim_pos > last_delim_pos+1
                            parsed = tryparse(T, str[(last_delim_pos+1):(delim_pos-1)])
                        end
                        result[n] = isnothing(parsed) ? undefval : parsed
                    end
                    last_delim_pos = delim_pos
                end
            end
        else
            for i = 1:length(str)
                if str[i] == delim
                    delim_pos = i
                    if delim_pos > last_delim_pos
                        n += 1
                        parsed = nothing
                        if delim_pos > last_delim_pos+1
                            parsed = tryparse(T, str[(last_delim_pos+1):(delim_pos-1)])
                        end
                        result[n] = isnothing(parsed) ? undefval : parsed
                        last_delim_pos = delim_pos
                    end
                end
            end
        end

        # Check for final value after last delim
        if length(str) > last_delim_pos
            n += 1
            parsed = tryparse(T, str[(last_delim_pos+1):length(str)])
            result[n] = isnothing(parsed) ? undefval : parsed
        end

        # Return the number of result values
        return n-offset
    end
    export delim_string_parse!

    """
    ```julia
    delim_string_parse(str, delim, T;
        \tmerge::Bool=false,
        \tundefval=NaN)
    ```
    Parse a delimited string `str` with delimiter `delim` into values of type `T`
    and return the answers as an array with eltype `T`
    """
    function delim_string_parse(str::AbstractString, delim::Char, T::Type=Float64; merge::Bool=false, undefval=NaN)

        # Allocate an array to hold our parsed results
        result = Array{T}(undef,ceil(Int,length(str)/2))

        # Parse the string
        n = delim_string_parse!(result, str, delim, T; merge=merge, undefval=undefval)

        # Return the result values
        return result[1:n]
    end
    export delim_string_parse

    """
    ```julia
    delim_string_function(f, str, delim, T;
        \tmerge::Bool=false,
    ```
    Parse a delimited string `str` with delimiter `delim` into substrings that will
    then be operated upon by function `f`. The results of `f` will be returned
    in an array with eltype `T`.
    """
    function delim_string_function(f::Function, str::AbstractString, delim::Char, T::Type; merge::Bool=false)

        # Max number of delimted values
        ndelims = 2
        for i = 1:length(str)
            if str[i] == delim
                ndelims += 1
            end
        end

        # Allocate output array
        result = Array{T}(undef,ceil(Int,ndelims))

        # Ignore initial delimiter
        last_delim_pos = 0
        if str[1] == delim
            last_delim_pos = 1
        end

        # Cycle through string parsing text betweeen delims
        delim_pos = 0
        n = 0
        if merge
            for i = 1:length(str)
                if str[i] == delim
                    delim_pos = i
                    if delim_pos > last_delim_pos+1
                        n += 1
                        if delim_pos > last_delim_pos+1
                            result[n] = f(str[(last_delim_pos+1):(delim_pos-1)])
                        end
                    end
                    last_delim_pos = delim_pos
                end
            end
        else
            for i = 1:length(str)
                if str[i] == delim
                    delim_pos = i
                    if delim_pos > last_delim_pos
                        n += 1
                        if delim_pos > last_delim_pos+1
                            result[n] = f(str[(last_delim_pos+1):(delim_pos-1)])
                        end
                        last_delim_pos = delim_pos
                    end
                end
            end
        end

        # Check for final value after last delim
        if length(str)>last_delim_pos
            n += 1
            result[n] = f(str[(last_delim_pos+1):length(str)])
        end

        # Return the result values
        return result[1:n]
    end
    export delim_string_function

    """
    ```julia
    parsedlm(str::AbstractString, delimiter::Char, T::Type=Float64; rowdelimiter::Char='\\n')
    ```
    Parse a string delimited by both row and column into a single (2-D) matrix. Default column delimiter is newline.
    Similar to `readdlm`, but operating on a string instead of a file.
    """
    function parsedlm(str::AbstractString, delimiter::Char, T::Type=Float64; rowdelimiter::Char='\n')
    	if T <: AbstractFloat
    		emptyval = T(NaN)
    	else
    		emptyval = zero(T)
    	end

    	# Count rows, and find maximum number of delimiters per row
    	numcolumns = maxcolumns = maxrows = 0
    	cₗ = delimiter
    	for c in str
    		(c == delimiter) && (numcolumns += 1)
    		if c == rowdelimiter
    			maxrows += 1
    			numcolumns += 1
    			# See if we have a new maximum, and reset the counters
    			(numcolumns > maxcolumns) && (maxcolumns = numcolumns)
    			numcolumns=0
    		end
    		cₗ = c
    	end
    	# If the last line isn't blank, add one more to the row counter
    	(cₗ != rowdelimiter) && (maxrows += 1)
    	maxchars = length(str)

    	# Allocate space for the imported array and fill with emptyval
    	parsedmatrix = fill(emptyval, maxrows, maxcolumns)

    	k = kₗ = 0 # Last delimiter position
    	for i = 1:maxrows
    		for j = 1:maxcolumns
    			while (k+1 <= maxchars) && (str[k+1] != delimiter) && (str[k+1] != rowdelimiter)
    				k += 1
    			end

    			# Otherwise, parse the string
    			parsed = tryparse(T, str[kₗ+1:k])
    			isnothing(parsed) || (parsedmatrix[i,j] = parsed)

                # If we're at the end of the string, move on
                (k+1 > maxchars) && break

    			# Step over the delimiter
    			kₗ = k += 1
    			# If we've hit a row delimiter, move to next row
    			(str[k] == rowdelimiter) && break
    		end
    	end
    	return parsedmatrix
    end
    export parsedlm

## --- Classifying imported datasets

    """
    ```julia
    plausiblynumeric(x)
    ```
    Return `true` if `x` can be parsed as a number, else `false`
    """
    function plausiblynumeric(x)
        if isa(x,Number)
            return true
        elseif isa(x,AbstractString) && tryparse(Float64,x) != nothing
            return true
        else
            return false
        end
    end
    export plausiblynumeric

    """
    ```julia
    nonnumeric(x)
    ```
    Return true for if `x` is not missing but cannot be parsed as a number
    """
    function nonnumeric(x)
        if isa(x,Number)
            false
        elseif isa(x,AbstractString) && (tryparse(Float64,x) != nothing || x == "")
            false
        else
            true
        end
    end
    export nonnumeric

## --- Transforming imported datasets

    """
    ```julia
    floatify(x, T::Type=Float64)
    ```
    Convert `x` to a floating-point number (default `Float64`) by any means necessary
    """
    function floatify(x, T::Type=Float64)
        if isa(x, Number)
            T(x)
        elseif isa(x,AbstractString) && tryparse(T,x) != nothing
            parse(T, x)
        else
            T(NaN)
        end
    end

    function _columnformat(x, standardize=true, floattype=Float64)
        if standardize
            if sum(plausiblynumeric.(x)) >= sum(nonnumeric.(x))
                return floatify.(x, floattype)
            else
                return string.(x)
            end
        else
            if all(xi -> isa(xi, AbstractString), x)
                return string.(x)
            elseif all(xi -> isa(xi, AbstractFloat), x)
                return float.(x)
            elseif all(xi -> isa(xi, Integer), x)
                return Integer.(x)
            else
                return x
            end
        end
    end

    function _sanitizevarname(s::AbstractString)
        s = replace(s, r"[\[\](){}]" => "") # Remove parentheses entirely
        s = replace(s, r"^([0-9])" => s"_\1") # Can't begin with a number
        s = replace(s, r"([\0-\x1F -/:-@\[-`{-~])" => s"_") # Everything else becomes an underscore
        return s
    end

    """
    ```julia
    elementify(data::Array, elements::Array=data[1,:];
        \timportas=:Dict,
        \tstandardize::Bool=true,
        \tfloattype=Float64,
        \tskipstart::Integer=1,
        \tskipnameless::Bool=true
    )
    ```
    Convert a flat array `data` into a dictionary (`importas=:Dict`) or named
    tuple (`importas=:Tuple`) with each column as a variable.
    Tuples are substantially more efficient, so should be favored where possible.
    """
    function elementify(data::Array, elements::Array=data[1,:];
            importas=:Dict,
            skipstart::Integer=1,
            standardize::Bool=true,
            floattype=Float64,
            skipnameless::Bool=true,
            sumduplicates::Bool=false
        )
        if importas == :Dict || importas==:dict
            # Output as dictionary
            if standardize
                # Constrain types somewhat for a modicum of type-stability
                if 1+skipstart == size(data,1)
                    result = Dict{String,Union{Vector{String}, String, Float64}}()
                else
                    result = Dict{String,Union{Vector{String}, Vector{Float64}}}()
                end
            else
                result = Dict{String, Any}()
            end

            # Process elements array
            elements = string.(elements)
            if skipnameless
                elements = elements[elements .!= ""]
            end
            result["elements"] = elements

            # Parse the input array, minus empty-named columns
            for i = 1:length(elements)
                if 1+skipstart == size(data,1)
                    column = data[end,i]
                else
                    column = data[(1+skipstart):end,i]
                end
                column_is_numeric = sum(plausiblynumeric.(column)) >= sum(nonnumeric.(column))

                if haskey(result,elements[i])
                    # If key already exists
                    if column_is_numeric && (standardize || (sum(plausiblynumeric.(result[elements[i]])) >= sum(nonnumeric.(result[elements[i]]))) )
                        # If either this column or the existing one is plausibly numeric, sum or average the two
                        if sumduplicates
                            result[elements[i]] = floatify.(result[elements[i]], floattype) + floatify.(column, floattype)
                        else
                            if 1+skipstart == size(data,1)
                                result[elements[i]] = nanmean( hcat(floatify.(result[elements[i]], floattype), floatify.(column, floattype)), dim=2 )[1]
                            else
                                result[elements[i]] = nanmean( hcat(floatify.(result[elements[i]], floattype), floatify.(column, floattype)), dim=2 )
                            end
                        end
                    elseif standardize
                        # If neither is numeric, but standardize is set, must return a string
                        result[elements[i]] = string.(result[elements[i]]) .* "|" .* string(lastcol)
                    else
                        # If neither is plausibly numeric, contatenate the columns and move on
                        result[elements[i]] = hcat(result[elements[i]], column)
                    end
                else
                    result[elements[i]] = _columnformat(column, standardize, floattype)
                end
            end

            # Return only unique elements, since dictionary keys must be unique
            result["elements"] = unique(elements)
            return result
        elseif importas==:Tuple || importas==:tuple || importas==:NamedTuple
            # Import as NamedTuple (more efficient future default)
            t = skipnameless ? elements .!= "" : isa.(elements, AbstractString)
            elements = _sanitizevarname.(elements[t])
            symbols = ((Symbol(e) for e ∈ elements)...,)
            values = [_columnformat(data[1+skipstart:end,i], standardize, floattype) for i in findall(t)]
            return NamedTuple{symbols}(values)
        end
    end
    export elementify

    """
    ```julia
    unelementify(dataset, elements;
        \tfloatout::Bool=false,
        \tfloattype=Float64,
        \tfindnumeric::Bool=false,
        \tskipnan::Bool=false,
        \trows=:
    )
    ```
    Convert a dict or named tuple of vectors into a 2-D array with variables as columns
    """
    function unelementify(dataset::Dict, elements::Array=sort(collect(keys(dataset)));
            floatout::Bool=false,
            floattype=Float64,
            findnumeric::Bool=false,
            skipnan::Bool=false,
            rows=:
        )

        # Find the elements in the input dict if they exist and aren't otherwise specified
        if any(elements .== "elements")
            elements = dataset["elements"]
        end

        # Figure out how many are numeric (if necessary), so we can export only
        # those if `findnumeric` is set
        if findnumeric
            is_numeric_element = Array{Bool}(undef,length(elements))
            for i = 1:length(elements)
                is_numeric_element[i] = sum(plausiblynumeric.(dataset[elements[i]])) > sum(nonnumeric.(dataset[elements[i]]))
            end
            elements = elements[is_numeric_element]
        end

        # Generate output array
        if floatout
            # Allocate output Array{Float64}
            result = Array{Float64}(undef, length(dataset[elements[1]][rows]), length(elements))

            # Parse the input dict. No column names if `floatout` is set
            for i = 1:length(elements)
                result[:,i] = floatify.(dataset[elements[i]][rows], floattype)
            end
        else
            # Allocate output Array{Any}
            result = Array{Any}(undef, length(dataset[elements[1]][rows])+1, length(elements))

            # Parse the input dict
            for i = 1:length(elements)
                # Column name goes in the first row, everything else after that
                result[1,i] = elements[i]
                result[2:end,i] .= dataset[elements[i]][rows]

                # if `skipnan` is set, replace each NaN in the output array with
                # an empty string ("") such that it is empty when printed to file
                # with dlmwrite or similar
                if skipnan
                    for n = 2:length(result[:,i])
                        if isa(result[n,i], AbstractFloat) && isnan(result[n,i])
                            result[n,i] = ""
                        end
                    end
                end
            end
        end
        return result
    end
    function unelementify(dataset::NamedTuple, elements=keys(dataset);
            floatout::Bool=false,
            floattype=Float64,
            findnumeric::Bool=false,
            skipnan::Bool=false,
            rows=:
        )
        # Figure out how many are numeric (if necessary), so we can export only
        # those if `findnumeric` is set
        if findnumeric
            elements = filter(x -> sum(plausiblynumeric.(dataset[x])) > sum(nonnumeric.(dataset[x])), elements)
        end

        # Generate output array
        if floatout
            # Allocate output Array{Float64}
            result = Array{floattype}(undef,length(dataset[elements[1]][rows]),length(elements))

            # Parse the input dict. No column names if `floatout` is set
            for i = 1:length(elements)
                result[:,i] = floatify.(dataset[elements[i]][rows], floattype)
            end
        else
            # Allocate output Array{Any}
            result = Array{Any}(undef,length(dataset[elements[1]][rows])+1,length(elements))

            # Parse the input dict
            for i = 1:length(elements)
                # Column name goes in the first row, everything else after that
                result[1,i] = string(elements[i])
                result[2:end,i] .= dataset[elements[i]][rows]

                # if `skipnan` is set, replace each NaN in the output array with
                # an empty string ("") such that it is empty when printed to file
                # with dlmwrite or similar
                if skipnan
                    for n = 2:length(result[:,i])
                        if isa(result[n,i], AbstractFloat) && isnan(result[n,i])
                            result[n,i] = ""
                        end
                    end
                end
            end
        end
        return result
    end
    export unelementify

## --- Concatenating / stacking datasets

    """
    ```julia
    concatenatedatasets(d1::AbstractDict, d2::AbstractDict)
    ```
    Vertically concatenate two Dict-based datasets, variable-by-variable
    """
    function concatenatedatasets(d1::AbstractDict, d2::AbstractDict)
        if isempty(keys(d1))
            result = d2
        elseif isempty(keys(d2))
            result = d1
        else
            # If 'elements' field doesn't exist, populate it
            if ~haskey(d1,"elements")
                d1["elements"] = sort(collect(keys(d1)))
            end
            if ~haskey(d2,"elements")
                d2["elements"] = sort(collect(keys(d2)))
            end

            # Find variable size
            s1 = size(d1[d1["elements"][1]])
            s2 = size(d2[d2["elements"][1]])

            # Combine datasets
            result = Dict()
            result["elements"] = d1["elements"] ∪ d2["elements"]
            for e in result["elements"]
                # Make any missing fields
                if haskey(d1, e) && ~haskey(d2, e)
                    if eltype(d1[e]) <: Number
                        d2[e] = fill(float(eltype(d1[e]))(NaN), s2)
                    else
                        d2[e] = fill("", s2)
                    end
                elseif haskey(d2, e) && ~haskey(d1, e)
                    if eltype(d2[e]) <: Number
                        d1[e] = fill(float(eltype(d2[e]))(NaN), s1)
                    else
                        d1[e] = fill("", s1)
                    end
                end

                # Combine fields
                result[e] = vcat(d1[e], d2[e])
            end
        end
        return result
    end
    export concatenatedatasets

## --- Renormalization of imported datasets

    """
    ```julia
    renormalize!(A::AbstractArray; dim, total=1.0)
    ```
    Normalize an array `A` in place such that it sums to `total`. Optionally may
    specify a dimension `dim` along which to normalize.
    """
    function renormalize!(A::AbstractArray; dim=:, total=1.0)
        current_sum = NaNStatistics._nansum(A, dim)
        A .*= total ./ current_sum
    end
    """
    ```julia
    renormalize!(dataset, [elements]; total=1.0)
    ```
    Normalize in-place a (i.e., compositional) `dataset` defined by a `Dict` or
    `NamedTuple` of one-dimensional numerical arrays, such that all the `elements`
    (i.e., variables -- by default all keys in the datset) sum to a given `total`
    (by default, `1.0`).

    Note that the arrays representing each element or variable are assumed to be
    of uniform length
    """
    function renormalize!(dataset::Union{Dict,NamedTuple}, elements=keys(dataset); total=1.0)
        # Note that this assumes all variables in the dataset are the same length!
        current_sum = zeros(size(dataset[first(keys(dataset))]))
        for e in elements
            current_sum .+= dataset[e] .|> x -> isnan(x) ? 0 : x
        end
        current_sum[current_sum .== 0] .= NaN

        for e in elements
            dataset[e] .*= total ./ current_sum
        end
        return dataset
    end
    export renormalize!

## --- High-level import/export functions

    """
    ```julia
    function importdataset(filepath, delim;
        \timportas=:Dict,
        \tstandardize::Bool=true,
        \tfloattype=Float64,
        \tskipstart::Integer=0,
        \tskipnameless::Bool=true,
        \tmindefinedcolumns::Integer=0
    )
    ```
    Import a delimited file specified by `filepath` with delimiter `delim` as a
    dataset in the form of either a `Dict` or a `NamedTuple`.

    Possible keyword arguments include:

        \timportas
    Specify the format of the imported dataset. Options include `:Dict` and `:Tuple`

        \tstandardize
    Convert columns to uniform type wherever possible. Boolean; `true` by default.

        \tfloattype
    Preferred floating-point type for numerical data. `Float64` by default.

        \tskipstart
    Ignore this many rows at the start of the input file (useful if input file has
    a header or other text before the column names). `0` by default.

        \tskipnameless
    Skip columns with no column name. Boolean; `true` by default

        \tmindefinedcolumns
    Skip rows with fewer than this number of delimiters. `0` by default.
    """
    function importdataset(filepath::AbstractString, delim::AbstractChar;
            importas=:Dict,
            standardize::Bool=true,
            floattype=Float64,
            skipstart::Integer=0,
            skipnameless::Bool=true,
            mindefinedcolumns::Integer=0
        )

        # Read file
        data = readdlm(filepath, delim, skipstart=skipstart)
        # Exclude rows with fewer than `mindefinedcolumns` columns
        if mindefinedcolumns > 0
            definedcolumns = vec(sum(.~ isempty.(data), dims=2))
            t = definedcolumns .>= mindefinedcolumns
            data = data[t,:]
        end

        return elementify(data,
            importas=importas,
            standardize=standardize,
            floattype=floattype,
            skipstart=1,
            skipnameless=skipnameless
        )
    end
    export importdataset

    """
    ```julia
    exportdataset(dataset, [elements], filepath, delim;
        \tfloatout::Bool=false,
        \tfindnumeric::Bool=false,
        \tskipnan::Bool=true,
        \tdigits::Integer,
        \tsigdigits::Integer
        \trows=:
    )
    ```
    Convert a dict or named tuple of vectors into a 2-D array with variables as columns
    Export a `dataset` (in the form of either a `Dict` or a `NamedTuple`),
    optionally specifying which `elements` to export, as a delimited ASCII text file
    with the name specified by `filepath` and delimiter `delim`.

    Possible keyword arguments include:

        \tdigits
        \tsigdigits
    Specify a number of absolute or significant digits to which to round the printed output.
    Default is no rounding.

        \tskipnan
    Leave `NaN`s as empty cells in the delimited output file. Boolean; `true` by default.

        \tfloatout
    Force all output to be represented as a floating-point number, or else `NaN`.
    Boolean; `false` by default.

        \tfindnumeric
    Export only numeric columns. Boolean; `false` by default.

        \trows
    specify which rows of the dataset to export. Default `:` exports all rows.
    """
    function exportdataset(dataset::Union{Dict,NamedTuple}, filepath::AbstractString, delim::AbstractChar;
            floatout::Bool=false,
            findnumeric::Bool=false,
            skipnan::Bool=true,
            digits::Integer=0,
            sigdigits::Integer=0,
            rows=:
        )

        # Convert dataset to flat 2d array
        data = unelementify(dataset,
            floatout=floatout,
            findnumeric=findnumeric,
            skipnan=skipnan,
            rows=rows
        )

        # Round output if applicable
        if sigdigits > 0
            data .= data .|> x -> isa(x, Number) ? round(x, sigdigits=sigdigits) : x
        elseif digits > 0
            data .= data .|> x -> isa(x, Number) ? round(x, digits=digits) : x
        end

        # Write to file
        return writedlm(filepath, data, delim)
    end
    # As above, but specifying which elements to export
    function exportdataset(dataset::Union{Dict,NamedTuple}, elements::Array, filepath::AbstractString, delim::AbstractChar;
            floatout::Bool=false,
            findnumeric::Bool=false,
            skipnan::Bool=true,
            digits::Integer=0,
            sigdigits::Integer=0,
            rows=:
        )

        # Convert dataset to flat 2d array
        data = unelementify(dataset, elements,
            floatout=floatout,
            findnumeric=findnumeric,
            skipnan=skipnan,
            rows=rows
        )

        # Round output if applicable
        if sigdigits > 0
            data .= data .|> x -> isa(x, Number) ? round(x, sigdigits=sigdigits) : x
        elseif digits > 0
            data .= data .|> x -> isa(x, Number) ? round(x, digits=digits) : x
        end

        # Write to file
        return writedlm(filepath, data, delim)
    end
    export exportdataset

## --- End of File

ext: jl


text:
using OrdinaryDiffEq, DiffEqDevTools, Test

function nonauto1(u, p, t)
  x, _ = u
  [
    t*x
    0
  ]
end

function nonauto2(u, p, t)
  _, y = u
  [
   y,
   t*y
  ]
end

function analytic(u0,p,t)
  x0, y0 = u0
  et = exp(t^2/2)
  [
   et*(x0 + t*y0),
   et*y0
  ]
end

u0 = [1.1, 2.2]
tspan = (0.0,1.0)
prob1 = ODEProblem(ODEFunction{true}((du,u,p,t)->du .= nonauto1(u,p,t) .+ nonauto2(u,p,t),
                              analytic=analytic),
                              u0,tspan)
prob2 = ODEProblem(ODEFunction{false}((u,p,t)-> nonauto1(u,p,t) .+ nonauto2(u,p,t),
                              analytic=analytic),
                              u0,tspan)
prob3 = SplitODEProblem(SplitFunction{true}((du,u,p,t)->du .= nonauto1(u,p,t),
                                       (du,u,p,t)->du .= nonauto2(u,p,t),
                              analytic=analytic),
                              u0,tspan)
prob4 = SplitODEProblem(SplitFunction{false}(nonauto1,
                                        nonauto2,
                              analytic=analytic),
                              u0,tspan)

testTol = 0.2

kwargs = (reltol=1e-16, abstol=1e-16)
for prob in [prob1, prob2, prob3, prob4]
  dts = 1 .//2 .^(12:-1:8)
  sim = test_convergence(dts,prob,KenCarp3(); kwargs...)
  @test sim.𝒪est[:l∞] ≈ 3 atol=testTol
  dts = 1 .//2 .^(9:-1:6)
  sim = test_convergence(dts,prob,CFNLIRK3(); kwargs...)
  @test sim.𝒪est[:l∞] ≈ 3 atol=testTol
  sim = test_convergence(dts,prob,KenCarp4(); kwargs...)
  @test sim.𝒪est[:l∞] ≈ 4 atol=testTol
  sim = test_convergence(dts,prob,KenCarp47(); kwargs...)
  @test sim.𝒪est[:l∞] ≈ 4 atol=testTol
  dts = 1 .//2 .^(7:-1:4)
  sim = test_convergence(dts,prob,KenCarp5(); kwargs...)
  @test sim.𝒪est[:l∞] ≈ 5 atol=testTol
  sim = test_convergence(dts,prob,KenCarp58(); kwargs...)
  @test sim.𝒪est[:l∞] ≈ 5 atol=testTol
end
for prob in [prob1, prob2]
  dts = 1 .//2 .^(7:-1:4)
  sim = test_convergence(dts,prob,ESDIRK54I8L2SA(); kwargs...)
  @test sim.𝒪est[:l∞] ≈ 5 atol=testTol
end

#=
prob = prob2
dts = 1 .//2 .^(9:-1:6)
sim = test_convergence(dts,prob,KenCarp4())
sim.𝒪est[:l∞]

using Plots
plot(sim)
=#


text:
export
    PreAllocatedStack

mutable struct PreAllocatedStack{T<:Integer, S}
    top::T
    stack::Vector{S}
end

push!(s::PreAllocatedStack) = (s.top += 1; s)
pop!(s::PreAllocatedStack) = (s.top -= 1; s)
getindex(s::PreAllocatedStack, idx) = s.stack[s.top + idx]


text:
module SubstituteM

export substituteGuess

# include("Utils.jl")
# include("FAfunc.jl")
# include("comparisons.jl")

using Main.UtilsM: removeExtras
using Main.FAM: fa
using Main.ComparisonM: faComparison

# Susbistution guess - ONLY WORKS WITH ENGLISH
function substituteGuess(text::String)

    engfreq = facomparison()
    textfreq = fa(uppercase(text), "ABCDEFGHIJKLMNOPQRSTUVWXYZ")
    sortedeng = sort(collect(engfreq), by = x -> x[2])
    sortedtext = sort(collect(textfreq), by = x -> x[2])
    for i in [1:sortedeng]

    end

end # function

end # module

max_stars_repo_name: Squalm/Cipher-Tools

max_issues_repo_name: Squalm/Cipher-Tools

max_forks_repo_name: Squalm/Cipher-Tools


text:
# """
#     responsetype(T)
# Return the numerical type of the time-response for composite LTISystem type `T`, e.g., `responsetype(ss(1.)) == Float64`
# If `T` is not a type, return `responsetype(typeof(T))`
# """
# responsetype{T<:Real}(::Type{T}) = T
# responsetype{T<:Complex}(::Type{T}) = promote_type(T.types...)
# responsetype(x) = typeof(x) # Fallback method
# responsetype{T}(::Type{T}) = T # Catch all

# Abstract type pyramid =============================================================

# NOTE: a lot that doesnt seem to be used...
#Base.promote_rule(::Type{StateSpace{T1}}, ::Type{StateSpace{T2}}) where {T1,T2} = StateSpace{promote_type(T1, T2)}

# NOTE: Is the below thing correct always?
Base.promote_rule(::Type{StateSpace{TE1,T1}}, ::Type{StateSpace{TE2,T2}}) where {TE1,TE2,T1,T2} =
    StateSpace{promote_type(TE1,TE2), promote_type(T1,T2)}

Base.promote_rule(::Type{StateSpace{TE1,T}}, ::Type{HeteroStateSpace{TE2,AT,BT,CT,DT}}) where {TE1,TE2,T,AT,BT,CT,DT} =
    HeteroStateSpace{promote_type(TE1, TE2),promote_type(Matrix{T},AT),promote_type(Matrix{T},BT),promote_type(Matrix{T},CT),promote_type(Matrix{T},DT)}

Base.promote_rule(::Type{TransferFunction{TE1,S1}}, ::Type{StateSpace{TE2,T2}}) where {TE1,TE2,T1,S1<:SisoTf{T1},T2} =
    StateSpace{promote_type(TE1, TE2), promote_type(T1,T2)}

Base.promote_rule(::Type{TransferFunction{TE1,S1}}, ::Type{DelayLtiSystem{T2}}) where {TE1,T1,S1<:SisoTf{T1},T2} =
    DelayLtiSystem{promote_type(T1,T2)}

Base.promote_rule(::Type{StateSpace{TE1,T1}}, ::Type{DelayLtiSystem{T2}}) where {TE1,T1,T2} =
    DelayLtiSystem{promote_type(T1,T2)}

Base.promote_rule(::Type{TransferFunction{TE1,S1}}, ::Type{TransferFunction{TE2,S2}}) where {TE1, TE2, S1, S2} =
    TransferFunction{promote_type(TE1,TE2),promote_type(S1, S2)}

#function Base.promote_rule{T<:StateSpace,P<:TransferFunction}(::Type{T}, ::Type{P})  #where T <: StateSpace where P <: TransferFunction
#    S = promote_type(primitivereal(P), primitivereal(T)) # TODO: this is not correct for P <: zpk, in which case the StateSpace system gets complex matrices
#    StateSpace{Matrix{S}}
#end

# Promotion of Number ==========================================================

Base.promote_rule(::Type{TransferFunction{TE,S}}, ::Type{T2}) where {TE, S<:SisoTf, T2<:Number} =
    TransferFunction{TE,promote_type(S, T2)}

# TODO Figure out a better way
function Base.promote_rule(::Type{StateSpace{TE, T1}}, ::Type{T2}) where {TE, T1, T2<:Number}
    StateSpace{TE, promote_type(T1, T2)}
end

function Base.promote_rule(::Type{SisoZpk{T1,C2}}, ::Type{T2}) where {T1, C2, T2<:Number}
    GainType = promote_type(T1, T2)
    return SisoZpk{GainType, complex(GainType)}
end
Base.promote_rule(::Type{SisoRational{T1}}, ::Type{T2}) where {T1, T2<:Number} = SisoRational{promote_type(T1,T2)}

# Promotion of Matrix ==========================================================

Base.promote_rule(::Type{TransferFunction{TE, SisoZpk{T1,TR1}}}, ::Type{MT}) where {TE, T1, TR1, T2, MT<:AbstractMatrix{T2}} =
    TransferFunction{TE, SisoZpk{T1, promote_type(TR1, T2)}}

Base.promote_rule(::Type{TransferFunction{TE, SisoRational{T1}}}, ::Type{MT}) where {TE, T1, T2, MT<:AbstractMatrix{T2}} =
    TransferFunction{TE, SisoRational{promote_type(T1, T2)}}

Base.promote_rule(::Type{StateSpace{TE, T1}}, ::Type{MT}) where {TE, T1, MT<:AbstractMatrix} =
    StateSpace{TE, promote_type(T,eltype(MT))}

Base.promote_rule(::Type{DelayLtiSystem{T1,S}}, ::Type{MT}) where {T1, S, MT<:AbstractMatrix} =
    DelayLtiSystem{promote_type(T1, eltype(MT)),S}


text:
# an idea: create one method named pivots with kwargs

function floorpivots{T,N}(ohlc::TimeArray{T,N})

    p  = lag(typical(ohlc))
    s1 = 2.*p .- lag(ohlc["High"])
    r1 = 2.*p .- lag(ohlc["Low"])
    s2 = p .- (r1 .- s1)
    r2 = (p .- s1) .+ r1
    s3 = p .- (r2 .- s1)
    r3 = (p .- s1) .+ r2

    TimeArray(s1.timestamp, [s3.values s2.values s1.values p.values r1.values r2.values r3.values],
              ["s3", "s2", "s1", "pivot", "r1", "r2", "r3"], ohlc.meta)

end

function woodiespivots{T,N}(ohlc::TimeArray{T,N})

    rng = lag(ohlc["High"]) .- lag(ohlc["Low"])

    p  = (lag(ohlc["High"]) .+ lag(ohlc["Low"]) .+ 2.*ohlc["Open"]) ./ 4
    s1 = 2.*p .- lag(ohlc["High"])
    r1 = 2.*p .- lag(ohlc["Low"])
    s2 = p .- rng
    r2 = p .+ rng
    s3 = s1 .- rng
    r3 = r1 .+ rng
    # s4 = s3 - range # can't get an answer that matches online calculators
    # r4 = r3 + range

    TimeArray(s1.timestamp, [s3.values s2.values s1.values p.values r1.values r2.values r3.values],
              ["s3", "s2", "s1", "pivot", "r1", "r2", "r3"], ohlc.meta)
end
#
# function camarillapivots{T,V}(hi::Array{SeriesPair{T,V},1},
#                               lo::Array{SeriesPair{T,V},1},
#                               cl::Array{SeriesPair{T,V},1})
#   #code here
#   ##  R4 = (H - L) * 1.1/2 + C
#   ##  R3 = (H - L) * 1.1/4 + C
#   ##  R2 = (H - L) * 1.1/6 + C
#   ##  R1 = (H - L) * 1.1/12 + C
#   ##  S1 = C - (H - L) * 1.1/12
#   ##  S2 = C - (H - L) * 1.1/6
#   ##  S3 = C - (H - L) * 1.1/4
#   ##  S4 = C - (H - L) * 1.1/2
# end
#
# function demarkpivots{T,V}(op::Array{SeriesPair{T,V},1},
#                            hi::Array{SeriesPair{T,V},1},
#                            lo::Array{SeriesPair{T,V},1},
#                            cl::Array{SeriesPair{T,V},1})
#   #code here
#   ## If Close < Open then X = (H + (L * 2) + C)
#   ## If Close > Open then X = ((H * 2) + L + C)
#   ## If Close = Open then X = (H + L + (C * 2))
#   ##
#   ## R1 = X / 2 - L
#   ## PP = X / 4
#   ## S1 = X / 2 - H
# end

function market_profile(x)
  #code here
end


text:
# this is basically log2_ext with a slightly modified range reduction.

function log21p_ext(x)
    u = 1.0 + x
    sig,n = frexp(u)
    au = 2sig; n -= 1

    # reduce to
    # k = 0, sqrt(3)/2 = 0.866 < au < sqrt(3/2) = 1.22
    # k = 1, sqrt(3/2) = 1.22 < au < sqrt(3) = 1.73

    if au < sqrt(3/2)
        bp = 1.0
        dp_h = 0.0
        dp_l = 0.0
    elseif au < sqrt(3)
        bp = 1.5
        # log2(1.5)
        dp_h = 5.84962487220764160156e-01
        dp_l = 1.35003920212974897128e-08
    else
        bp = 1.0
        dp_h = 0.0
        dp_l = 0.0

        n += 1
        au *= 0.5
    end

    # exact 1.0 + x = u + c
    if x >= 1.0
        # u <= 2x
        c = 1.0-(u-x)
    else
        # 1) x > -0.5: 0.5 <= u <= 2.0
        # 2) x <= -0.5: u = x+1.0 exact, hence so is u-1.0
        c = x-(u-1.0)
    end

    # x <= -0.5: c = 0
    # x < 0.5: |c| < eps(1.0)
    # x < 0x1p53: |c| in {-eps(x), 0, eps(x)}
    # otherwise: c = 1

    # scale to au:  1 + x = (au + ac)*2^n
    ac = ldexp(c,-n)
    # log(x) = 2s + 2/3 s^3 + 2/5 s^5...
    f = au - bp # usual f, reduced to -0.28 < u < 0.23
    f1 = f + ac
    f2 = ac + (f-f1)
    g = au + bp
    g_h = trunc32(g) # according to comment
    g_l = ac + (au + (bp - g_h))
    rg = 1.0/g

    ss = f1*rg # ss = (x-1)/(x+1) or (x-1.5)/(x+1.5) < 0.125 = 2^-3
    s_h = trunc32(ss)
    s_l = (((f1-s_h*g_h)-s_h*g_l) + f2)*rg # division correction?

    z = ss*ss # < 2^-6
    # log(au) = 2s + 2/3*(s^3 + r) = 2/3s (3 + s^2 + r)
    r = z*z*@horner(z,
                    5.99999999999994648725e-01,
                    4.28571428578550184252e-01,
                    3.33333329818377432918e-01,
                    2.72728123808534006489e-01,
                    2.30660745775561754067e-01,
                    2.06975017800338417784e-01) # < 2^-12
    # (s_h + s_l)^2 = s_h^2 + 2 s_l s_h + s_l s_l
    #   = s_h^2 + s_l (s_h + ss)
    r += s_l*(s_h+ss)
    s2 = s_h*s_h
    t_h = 3.0 + s2 + r
    t_h = trunc32(t_h)
    t_l = r-((t_h-3.0)-s2)
    # u+v = ss*(1+...)
    u = s_h*t_h
    v = s_l*t_h+t_l*ss

    # 2/(3log2)
    cp    =  9.61796693925975554329e-01 # /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */
    cp_h  =  9.61796700954437255859e-01 # /* 0x3FEEC709, 0xE0000000 =(float)cp */
    cp_l  = -7.02846165095275826516e-09 # /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h*/

    # 2/(3log2)*(ss+...)
    p_h = trunc32(u+v)
    p_l = v-(p_h-u)
    z_h = cp_h*p_h
    z_l = cp_l*p_h + p_l*cp + dp_l

    t = float64(n)
    t1 = trunc32(((z_h+z_l) + dp_h) + t)
    t2 = z_l - (((t1 - t) - dp_h) - z_h)
    Double64(t1,t2)
end

function exp2m1(p::Double64)
    z = p.lo + p.hi
    # check for under/overflow

    # 2^(p_l+p_h)

    zr = round(z)
    n = int(zr)
    p_h = p.hi - zr
    p_l = p.lo

    lg2  =  6.93147180559945286227e-01 # /* 0x3FE62E42, 0xFEFA39EF */
    lg2_h  =  6.93147182464599609375e-01 # /* 0x3FE62E43, 0x00000000 */
    lg2_l  = -1.90465429995776804525e-09 # /* 0xBE205C61, 0x0CA86C39 */

    pt = trunc32(p_l + p_h)
    u = pt*lg2_h
    v = (p_l-(pt-p_h))*lg2+pt*lg2_l
    z = u+v
    w = v-(z-u)

    expm1_kernel(n,Double64(z,w))
end

# 2^n - 1

# |z| < 0.5*log(2) = 0.35
# u ~ z
# r ~ -z^2/2

pow1p(x,y) = exp2(mul_ext(log21p_ext(x),y))
powm1(x,y) = exp2m1(mul_ext(log2_ext(x),y))
pow1pm1(x,y) = exp2m1(mul_ext(log21p_ext(x),y))


text:
using LogNormals
using Test, Distributions, LinearAlgebra, StatsBase, Missings, Random
using OffsetArrays, RecursiveArrayTools

function boot_sem_cor()
    #using FillArrays
    ##using BandedMatrices
    ##using Bootstrap
    nsum = 100
    probgap = 0.6
    acf0 = [1,0.4,0.1]
    Sigma = cormatrix_for_acf(nsum, acf0);
    dmn = MvNormal(ones(nsum), Symmetric(Sigma));
    #always use the same gaps
    igap = sample(1:nsum,trunc(Int,nsum*probgap), replace=false)
    b = allowmissing(rand(dmn));
    b[igap] .= missing;
    sem_cor(b, acf0, ExactMissing())
    nboot = 10_000
    resboot = map(1:nboot) do iboot
        b = allowmissing(rand(dmn));
        b[igap] .= missing
        sum(skipmissing(b)), mean(skipmissing(b)), sem_cor(b, acf0, ExactMissing())
    end;
    sums_boot, means_boot, semcor_boot = vectuptotupvec(resboot);
    mean(sums_boot), std(sums_boot)
    mean(means_boot), std(means_boot)
    dvn = SimpleDistributionVector(Fill(Normal(1,1), nsum)...);
    dvn[ismissing.(b)] .= missing;
    dsum = sum(dvn, AutoCorrelationFunction(acf0), SkipMissing())
    dmean = mean(dvn, AutoCorrelationFunction(acf0), SkipMissing())
    @test dmean.σ ≈ std(means_boot) atol=0.01
    # mean over Normals with missings works
    function sem_cor_benchmark(x, acfe, ms::MissingStrategy=ExactMissing())
        n = length(x)
        n <= 1 && return(std(x))
        nmiss = count(ismissing.(x))
        nfin = n - nmiss
        neff = effective_n_cor(x, acfe)
        σ2uncorr = var(skipmissing(x))
        # BLUE Var(x) for correlated: Zieba11 eq.(1)
        σ2 = σ2uncorr*(nfin-1)*neff/(nfin*(neff-1))
        if ms === ExactMissing() && (nmiss != 0)
            σ = Fill(√σ2, n)
            dv = ParamDistributionVector(Normal{nonmissingtype(eltype(x))}, x, σ)
            acfes = n < length(acfe) ? view(acfe,1:n) : acfe
            dm = mean(dv, AutoCorrelationFunction(acfes), SkipMissing())
            #@show acfes, dm, length(dv), count(ismissing,dv), count(ismissing,x)
            return(dm.σ)
        else
            return(√(σ2/neff))
        end
    end
    res_sem_cor = sem_cor(b, acf0)
    @test res_sem_cor ≈ std(means_boot) atol=0.02
    res_sem_cor2a = sem_cor_benchmark(b, acf0, ms=ExactMissing())
    @test res_sem_cor2a ≈ std(means_boot) atol=0.02
    res_sem_cor2b = sem_cor_benchmark(b, acf0, ms=SkipMissing())
    @test res_sem_cor2b ≈ std(means_boot) atol=0.02
    #using BenchmarkTools: @btime
    @btime sem_cor_benchmark($b, $acf0, ms=ExactMissing())
    @btime sem_cor_benchmark($b, $acf0, ms=SkipMissing())
    # neff based method much faster
    function tmpf()
        #using StatsPlots
        plot(density(means_boot))
        vline!([mean(means_boot), quantile(means_boot,[0.025,0.975])...])
        #
        # plot(density(means_boot), label = "density mean")
        # vline!([mean(means_boot), quantile(means_boot,[0.025,0.975])...], label = "stats mean")
        density(semcor_boot, label = "density sem_cor")
        vline!([mean(semcor_boot), quantile(semcor_boot,[0.025,0.975])...], label = "stats sem_cor")
        vline!([std(means_boot)], label = "std(means_bootstrap)") # true mu and sigma
        vline!([dmean.σ], label = "sum_normals - true distribs") # true mu and sigma
    end
end

function inspect_lags_autocorrelation()
    nsum = 1_000
    acf0 = [1,0.4,0.1]
    Sigma = cormatrix_for_acf(nsum, acf0);
    dmn = MvNormal(ones(nsum), Symmetric(Sigma));
    a = rand(dmn);
    probgap = 0.4
    clustersize = 2
    #clustersize = 50
    nboot = 1_000
    resboot = map(1:nboot) do iboot
        igap0 = sample(0:(nsum-clustersize),trunc(Int,nsum*probgap/clustersize), replace=true);
        b = allowmissing(a);
        # outer product see https://stackoverflow.com/a/44592419
        #igap0.+Base.OneTo(clustersize)'
        b[igap0.+Base.OneTo(clustersize)'] .= missing;
        ae = autocor(b, 1:2, ms=ExactMissing())
        af = autocor(b, 1:2, ms=SkipMissing())
        (ae, af)
    end;
    ae, af = VectorOfArray.(vectuptotupvec(resboot));
    mean(ae; dims=2)
    mean(af; dims=2)
    autocor(a, 1:2)

end

function benchmark_count_forlags()
    #using BenchmarkTools: @btime
    function count_forlags_slower(pred, x,lags)
        cnt = zeros(size(lags))
        ax = OffsetArrays.no_offset_view(axes(x,1))
        lx = length(x)
        for i in 1:lx
            for (ik,k) in enumerate(lags)
                i > lx-k && continue
                if (pred(x[ax[i]])::Bool || pred(x[ax[i+k]])::Bool); cnt[ik] += 1; end
            end
        end
        cnt
    end
    x = repeat([1,2,missing,missing,5], 500)
    lags = 0:4
    @btime count_forlags($ismissing, $x, $lags)
    @btime count_forlags_slower($ismissing, $x, $lags)
    # the version with calling subfunction is faster
    @code_warntype count_forlags(ismissing, x, lags)
    @code_warntype LogNormals.count_forlags2(ismissing, x, lags)
    @code_warntype LogNormals.count_forlag(ismissing, x, 2)
end


text:
baremodule AsyncFinalizers

function onfinalize end
function unsafe_unwrap end

module Internal

using Base.Threads: atomic_cas!, atomic_or!, atomic_fence

using ..AsyncFinalizers: AsyncFinalizers

if isfile(joinpath(@__DIR__, "config.jl"))
    include("config.jl")
else
    include("default-config.jl")
end

include("utils.jl")
include("bags.jl")
include("shims.jl")
include("runtime.jl")

function __init__()
    reinit()
    atexit(onexit)
end

# Use README as the docstring of the module:
@doc let path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    replace(read(path, String), r"^```julia"m => "```jldoctest README")
end AsyncFinalizers

end  # module Internal

end  # baremodule AsyncFinalizers


text:
"""
    function burnin(;
        players::Vector{Player},
        seasons::Int64,
        strength_method::String,
        verbose::Bool
    )

Generate an initial "burn in" ranking for the players.

Arguments
=========
- `players`:
    an array of players in the simulation.
- `seasons`:
    the number of seasons to simulate.
- `strength_method`:
    the method to determine probability team 1 beats team 2.
"""

function burnin(;
    players::Vector{Player},
    seasons::Int64,
    strength_method::String,
    verbose::Bool
)
    tournaments = generate_doubles_tournaments()
    sort!(tournaments, by = t -> t.date, alg=InsertionSort)

    # Simulate tournaments for each season
    for _ in 1:seasons
        # Simulate each tournament
        for tournament in tournaments
            play_doubles_tournament(
                tournament = tournament,
                players = players,
                ranking_scheme = "Burn In",
                strength_method = strength_method,
                verbose = verbose
            )
            updateRanking(
                players = players,
                ranking_scheme_to_apply = "ATP",
                ranking_scheme_to_update = "Burn In",
                last_tournament_id = tournament.id,
                verbose = verbose
            )
        end
    end
end

max_stars_repo_name: Danial-Hussain/Tennis-Ranking

max_issues_repo_name: Danial-Hussain/Tennis-Ranking

max_forks_repo_name: Danial-Hussain/Tennis-Ranking


text:
module Cubed

using Random
using ProgressMeter
using CUDA
using CSV, DataFrames

include("types.jl")
export System, Dynamic, Potential
include("msd.jl")
export msd!, savemsd
include("kernels.jl")
export gpu_energy!, gpu_ermak!
include("utils.jl")
export init!
include("move.jl")
export move!
include("simulation.jl")
export simulate

end # module


text:
abstract type AbstractChargeDriftModel{T <: SSDFloat} end
abstract type AbstractTemperatureModel{T <: SSDFloat} end

include("ElectricFieldChargeDriftModel/ElectricFieldChargeDriftModel.jl")
include("ADLChargeDriftModel/ADLChargeDriftModel.jl")


text:
using Makie

scene = Scene()

f(t, v, s) = (sin(v + t) * s, cos(v + t) * s, (cos(v + t) + sin(v)) * s)
t = Node(Base.time()) # create a life signal
limits = FRect3D(Vec3f0(-1.5, -1.5, -3), Vec3f0(3, 3, 6))
p1 = meshscatter!(scene, lift(t-> f.(t, range(0, stop = 2pi, length = 50), 1), t), markersize = 0.05)[end]
p2 = meshscatter!(scene, lift(t-> f.(t * 2.0, range(0, stop = 2pi, length = 50), 1.5), t), markersize = 0.05)[end]

lines = lift(p1[1], p2[1]) do pos1, pos2
    map((a, b)-> (a, b), pos1, pos2)
end
linesegments!(scene, lines, linestyle = :dot, limits = limits)
# record a video
N = 150
record(scene, "record_video.gif", 1:N) do i
    push!(t, Base.time())
end


text: Dict{String, Any}("a.b" => Dict{String, Any}("À" => Dict{String, Any}()), "À" => Dict{String, Any}("value" => "latin capital letter A with grave", "type" => "string"), "\"quoted\"" => Dict{String, Any}("quote" => Dict{String, Any}("value" => "true", "type" => "bool")), "\"" => Dict{String, Any}("value" => "just a quote", "type" => "string"), "\n" => Dict{String, Any}("value" => "newline", "type" => "string"), "backsp\b\b" => Dict{String, Any}())


text:
module Probabilistic

using LogicCircuits
using ..Utils

export

# ProbCircuits
ProbΔNode, ProbΔ, ProbΔ, ProbLeafNode, ProbInnerNode,
ProbLiteral, Prob⋀, Prob⋁, ProbCache, variable, num_parameters, compute_log_likelihood,
log_proba,
log_likelihood, estimate_parameters, log_likelihood_per_instance, marginal_log_likelihood_per_instance,
initial_mixture_model, estimate_parameters_from_aggregates, compute_ensemble_log_likelihood,
expectation_step, maximization_step, expectation_step_batch, train_mixture_with_structure, check_parameter_integrity,
ll_per_instance_per_component, ll_per_instance_for_ensemble,estimate_parameters_cached,
sample,
MPE, MAP,prob_origin,

# ProbFlowCircuits
marginal_pass_up, marginal_pass_down, marginal_pass_up_down,

# Mixtures
Mixture, AbstractFlatMixture, FlatMixture, FlatMixtureWithFlow,component_weights,FlatMixtureWithFlows,
log_likelihood, log_likelihood_per_instance, log_likelihood_per_instance_component,
init_mixture_with_flows, reset_mixture_aggregate_flows, aggregate_flows, estimate_parameters,
AbstractMetaMixture, MetaMixture,AbstractFlatMixture,AbstractMixture, components, num_components,

# EM Learner
train_mixture,

# Bagging
bootstrap_samples_ids, learn_mixture_bagging, learn_mixture_bagging2,
init_bagging_samples, train_bagging,

# VtreeLearner
MetisContext, metis_top_down, BlossomContext, blossom_bottom_up!,
test_top_down, test_bottom_up!,learn_vtree_bottom_up,

# MutualInformation
mutual_information, DisCache, conditional_entropy, sum_entropy_given_x,

# Clustering
clustering,

# Queries
pr_constraint, psdd_entropy, psdd_kl_divergence

include("Clustering.jl")
include("ProbCircuits.jl")
include("ProbFlowCircuits.jl")
include("MutualInformation.jl")
include("Mixtures.jl")
include("Bagging.jl")
include("EMLearner.jl")
include("VtreeLearner.jl")
include("Queries.jl")

end


text:
module Complex #hide
# # Complex numbers
# Trivial example to demonstrate the capability of handling complex numbers
using PetscWrap
using SlepcWrap

# Start by checking that you're using a PETSc/SLEPc build with complex number
(PetscScalar <: Real) && error("You must configure PETSc/SLEPc with complex numbers to run this example")

# Initialize SLEPc
SlepcInitialize()

# Create matrix
A = create_matrix(4, 4; auto_setup = true)

# Get rows handled by the local processor
A_rstart, A_rend = get_range(A)

# Create diag matrix with complex numbers
J = I = A_rstart:A_rend
V = im .* I
set_values!(A, I, J, V)
assemble!(A)

# Now we set up the eigenvalue solver
eps = create_eps(A; auto_setup = true)

# Then we solve
solve!(eps)

# Show converged eivenvalues
@show get_eigenvalues(eps)

# Finally, let's free the memory
destroy!(A)
destroy!(eps)

# And call finalize when you're done
SlepcFinalize()

end #hide


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule ZeroMQ_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("ZeroMQ")
JLLWrappers.@generate_main_file("ZeroMQ", UUID("8f1865be-045e-5c20-9c9f-bfbfb0764568"))
end  # module ZeroMQ_jll


text:
#=
Let d(n) be defined as the sum of proper divisors of n (numbers less than n which divide evenly into n).
If d(a) = b and d(b) = a, where a ≠ b, then a and b are an amicable pair and each of a and b are called amicable numbers.

For example, the proper divisors of 220 are 1, 2, 4, 5, 10, 11, 20, 22, 44, 55 and 110; therefore d(220) = 284. The proper divisors of 284 are 1, 2, 4, 71 and 142; so d(284) = 220.

Evaluate the sum of all the amicable numbers under 10000.
=#

function sum_proper_divisors(n)
  i = 1
  sum = 0
  while (i < n)
    if n%i==0
      sum += i
    end
    i += 1
  end
  sum
end

function calc()
  sum = 0
  for i in 1:9_999
    a = sum_proper_divisors(i)
    b = sum_proper_divisors(a)
    if i == b && a != b
      sum += i
    end
  end
  sum
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
"Do one iteration and get the current value of a Xorshift RNG object."
xorshift_next

let doc_xorshift64(r) = """
```julia
$r <: AbstractXorshift64
$r([seed])
```

$r RNG. The `seed` will be automatically convert to an `UInt64` number. A zero seed is not acceptable.
"""
    @doc doc_xorshift64("Xorshift64") Xorshift64
    @doc doc_xorshift64("Xorshift64Star") Xorshift64Star
end
@deprecate Xorshift64 Xoshiro256StarStar
@deprecate Xorshift64Star Xoshiro256StarStar

let doc_xorshift128(r) = """
```julia
$r <: AbstractXorshift128
$r([seed])
```

$r RNG. The `seed` can be a `Tuple` of two `UInt64`s, or an `Integer` which will be automatically convert to
an `UInt128` number. Zero seeds are not acceptable.
"""
    @doc doc_xorshift128("Xorshift128") Xorshift128
    @doc doc_xorshift128("Xorshift128Star") Xorshift128Star
    @doc doc_xorshift128("Xorshift128Plus") Xorshift128Plus
end
@deprecate Xorshift128 Xoshiro256StarStar
@deprecate Xorshift128Star Xoshiro256StarStar
@deprecate Xorshift128Plus Xoshiro256StarStar

let doc_xorshift1024(r) = """
```julia
$r <: AbstractXorshift1024
$r([seed...])
```

$r RNG. The `seed` can be a `Tuple` of 16 `UInt64`s, or several (no more than 16) `Integer`s which will all
be automatically converted to `UInt64` numbers. Zero seeds are not acceptable.
"""
    @doc doc_xorshift1024("Xorshift1024") Xorshift1024
    @doc doc_xorshift1024("Xorshift1024Star") Xorshift1024Star
    @doc doc_xorshift1024("Xorshift1024Plus") Xorshift1024Plus
end
@deprecate Xorshift1024 Xoshiro256StarStar
@deprecate Xorshift1024Star Xoshiro256StarStar
@deprecate Xorshift1024Plus Xoshiro256StarStar
# TODO: implement Xoroshiro1024Star / Xoroshiro1024StarStar and Xoroshiro512StarStar / Xoroshiro512Plus

let doc_xoroshiro64(r) = """
```julia
$r <: AbstractXoroshiro64
$r([seed])
```

$r RNG. The `seed` can be a `Tuple` of two `UInt32`s, or an `Integer` which will be automatically convert to
an `UInt64` number. Zero seeds are not acceptable.
"""
    @doc doc_xoroshiro64("Xoroshiro64Star") Xoroshiro64Star
    @doc doc_xoroshiro64("Xoroshiro64StarStar") Xoroshiro64StarStar
end

let doc_xoroshiro128(r) = """
```julia
$r <: AbstractXoroshiro128
$r([seed])
```

$r RNG. The `seed` can be a `Tuple` of two `UInt64`s, or an `Integer` which will be automatically convert to
an `UInt128` number. Zero seeds are not acceptable.
"""
    @doc doc_xoroshiro128("Xoroshiro128") Xoroshiro128
    @doc doc_xoroshiro128("Xoroshiro128Star") Xoroshiro128Star
    @doc doc_xoroshiro128("Xoroshiro128Plus") Xoroshiro128Plus
    @doc doc_xoroshiro128("Xoroshiro128StarStar") Xoroshiro128StarStar
end
@deprecate Xoroshiro128 Xoroshiro128Plus
@deprecate Xoroshiro128Star Xoroshiro128StarStar

let doc_xoshiro128(r) = """
```julia
$r <: AbstractXoshiro128
$r([seed])
```

$r RNG. The `seed` can be a `Tuple` of four `UInt32`s, or an `Integer` which will be automatically convert to
an `UInt128` number. Zero seeds are not acceptable.
"""
    @doc doc_xoshiro128("Xoshiro128Plus") Xoshiro128Plus
    @doc doc_xoshiro128("Xoshiro128StarStar") Xoshiro128StarStar
end

let doc_xoshiro256(r) = """
```julia
$r <: AbstractXoshiro256
$r([seed])
```

$r RNG. The `seed` can be a `Tuple` of four `UInt64`s, or an `Integer` which will be automatically convert to
an `UInt64` number (and then is initialized with SplitMix64). Zero seeds are not acceptable.
"""
    @doc doc_xoshiro256("Xoshiro256Plus") Xoshiro256Plus
    @doc doc_xoshiro256("Xoshiro256StarStar") Xoshiro256StarStar
end


text:
# Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors    #src
# This Source Code Form is subject to the terms of the Mozilla Public License   #src
# v.2.0. If a copy of the MPL was not distributed with this file, You can       #src
# obtain one at https://mozilla.org/MPL/2.0/.                                   #src

# # The urban planning problem

# An "urban planning" problem based on an [example from puzzlor](http://www.puzzlor.com/2013-08_UrbanPlanning.html).

using JuMP
import GLPK
import Test

function example_urban_plan()
    model = Model(GLPK.Optimizer)
    ## x is indexed by row and column
    @variable(model, 0 <= x[1:5, 1:5] <= 1, Int)
    ## y is indexed by R or C, the points, and an index in 1:5. Note how JuMP
    ## allows indexing on arbitrary sets.
    rowcol = ["R", "C"]
    points = [5, 4, 3, -3, -4, -5]
    @variable(model, 0 <= y[rowcol, points, 1:5] <= 1, Int)
    ## Objective - combine the positive and negative parts
    @objective(model, Max, sum(
          3 * (y["R", 3, i] + y["C", 3, i])
        + 1 * (y["R", 4, i] + y["C", 4, i])
        + 1 * (y["R", 5, i] + y["C", 5, i])
        - 3 * (y["R", -3, i] + y["C", -3, i])
        - 1 * (y["R", -4, i] + y["C", -4, i])
        - 1 * (y["R", -5, i] + y["C", -5, i])
        for i in 1:5)
    )
    ## Constrain the number of residential lots
    @constraint(model, sum(x) == 12)
    ## Add the constraints that link the auxiliary y variables to the x variables
    for i = 1:5
        @constraints(model, begin
            ## Rows
            y["R", 5, i] <= 1 / 5 * sum(x[i, :]) # sum = 5
            y["R", 4, i] <= 1 / 4 * sum(x[i, :]) # sum = 4
            y["R", 3, i] <= 1 / 3 * sum(x[i, :]) # sum = 3
            y["R", -3, i] >= 1 - 1 / 3 * sum(x[i, :]) # sum = 2
            y["R", -4, i] >= 1 - 1 / 2 * sum(x[i, :]) # sum = 1
            y["R", -5, i] >= 1 - 1 / 1 * sum(x[i, :]) # sum = 0
            ## Columns
            y["C", 5, i] <= 1 / 5 * sum(x[:, i]) # sum = 5
            y["C", 4, i] <= 1 / 4 * sum(x[:, i]) # sum = 4
            y["C", 3, i] <= 1 / 3 * sum(x[:, i]) # sum = 3
            y["C", -3, i] >= 1 - 1 / 3 * sum(x[:, i]) # sum = 2
            y["C", -4, i] >= 1 - 1 / 2 * sum(x[:, i]) # sum = 1
            y["C", -5, i] >= 1 - 1 / 1 * sum(x[:, i]) # sum = 0
        end)
    end
    ## Solve it
    optimize!(model)
    Test.@test termination_status(model) == MOI.OPTIMAL
    Test.@test primal_status(model) == MOI.FEASIBLE_POINT
    Test.@test objective_value(model) ≈ 14.0
    return
end

example_urban_plan()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
println("Real states:")
foreach(println, StateNamePuzzle.solve(StateNamePuzzle.realnames))

println("\nReal + fictitious state:")
foreach(println, StateNamePuzzle.solve(vcat(StateNamePuzzle.realnames,
    StateNamePuzzle.fictitious)))


text:
type Convergents
  q::PartialQuotients
  accuracy::Int
end

function Convergents(q::PartialQuotients; accuracy = 10)
  if q.accuracy < accuracy && q.periodlength == 0
    # To do : add a warn
    accuracy = q.accuracy
  end
  return Convergents(q, accuracy)
end

function pos2index(period, pos)
  (period == 0) && return pos
  (pos <= period + 1) && return pos
  return 2 + (pos - 2) % period
end

# state = (n, h_{n}, h_{n-1}, k_{n}, k_{n-1})
function Base.start(c::Convergents)
  return (0, big(1), big(0), big(0), big(1))
end

function Base.done(c::Convergents, state)
  return state[1] == c.accuracy
end

function Base.next(c::Convergents, state)
  n, h1, h2, k1, k2 = state
  n += 1
  i = pos2index(c.q.periodlength, n)
  h = c.q.partialquotients[i] * h1 + h2
  k = c.q.partialquotients[i] * k1 + k2
  return h // k, (n, h, h1, k, k1)
end


text:
# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Script with functions to import and transform raw data
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

using DataFrames, CSV, XLSX, LinearAlgebra, Statistics

# other scripts
dir_home = "x:/VIVES/1-Personal/Florian/git/2022_01_comext/src/"
include(dir_home * "functions.jl")

# location of data input/output (io)
dir_io = "C:/Users/u0148308/data/comext/"
dir_dropbox = "C:/Users/u0148308/Dropbox/BREXIT/"

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Import data into Julia and do initial cleaning
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# initialize DataFrame
# Notes:
#   - column names and types need to correspond to the final output of the "inital_cleaning" function
df = DataFrame(DECLARANT_ISO=String[], PARTNER_ISO=String[], TRADE_TYPE=String[], PRODUCT_NC=String[], FLOW=String[], PERIOD=Int64[],
               VALUE_IN_EUROS=Union{Missing, Float64}[], QUANTITY_IN_KG=Union{Missing, Float64}[],
               PRODUCT_SITC=String[], PRODUCT_CPA2008=String[], PRODUCT_CPA2002=String[], PRODUCT_CPA2_1=String[],
               PRODUCT_BEC=String[], PRODUCT_BEC5=String[], PRODUCT_SECTION=String[])

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

# for now just use 3 months otherwise Julia crashes since ~30GB of data
append!(df, initial_cleaning("2001", "01"))
append!(df, initial_cleaning("2020", "02"))
append!(df, initial_cleaning("2020", "03"))

# possibly need to sort before grouping, seems that grouping sorts incorrectly sometimes?
sort!(df, [:DECLARANT_ISO, :PARTNER_ISO, :PRODUCT_NC, :FLOW, :TRADE_TYPE, :PERIOD]) # do sorting just in case (unexplainable problem below)

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Unit prices
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# compute UNIT_PRICE
transform!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] => ByRow((v,q) -> v/q) => :UNIT_PRICE)

# compute UNIT_PRICE_CHANGE
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "PRODUCT_NC", "FLOW"]
gdf = groupby(df, cols_grouping)
df = transform(gdf, :UNIT_PRICE => pct_change => :UNIT_PRICE_CHANGE)

# compute MOM_PRICE_CHANGE
gdf = groupby(df, cols_grouping)
df = transform(gdf, [:PERIOD, :UNIT_PRICE] => mom_change => :MOM)

# export data
df_export = ifelse.(ismissing.(df), NaN, df)[1:30_000,:] # excel cannot deal with so many rows
XLSX.writetable(dir_io * "clean/" * "df_prices.xlsx", df_export, overwrite=true)

CSV.write(dir_io * "clean/" * "df_prices.csv", df)

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Share of Belgium
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# Notes:
#   - how to treat missing values? skipmissing for now
#   - what constitutes as total?
#       + for DECLARANT_ISO the total are all EU ctrys, termed EU
#           -  change in composition, most importantly GB exiting will increase shares
#       + for PARTNER_ISO the total is the entire world, termed WORLD
#       + the total will thus be the BE share of EU imports/exports to the entire world (intra + extra)

# compute EU_WORLD (EU imports/exports to the entire world)
cols_grouping = ["PRODUCT_NC", "FLOW", "PERIOD"]
prod_class_other = ["PRODUCT_SITC", "PRODUCT_CPA2002", "PRODUCT_CPA2008", "PRODUCT_CPA2_1", "PRODUCT_BEC", "PRODUCT_BEC5"]
cols_grouping = [cols_grouping; prod_class_other]

gdf = groupby(df, cols_grouping)
df_sum_EU_WORLD = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> (x -> sum(skipmissing(x))), renamecols=false)
df_sum_EU_WORLD.TRADE_TYPE .= "total"

# compute EU_EXTRA and EU_INTRA (EU imports/exports to extra/intra EU countries)
cols_grouping = ["TRADE_TYPE", "PRODUCT_NC", "FLOW", "PERIOD"]
cols_grouping = [cols_grouping; prod_class_other]
gdf = groupby(df, cols_grouping)
df_sum_EU_EXTRA_INTRA = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> (x -> sum(skipmissing(x))), renamecols=false)

# compute BE_WORLD (BE imports/exports to the entire world)
df_BE = subset(df, :DECLARANT_ISO => ByRow(x -> x == "BE"))
cols_grouping = ["PRODUCT_NC", "FLOW", "PERIOD"]
cols_grouping = [cols_grouping; prod_class_other]

gdf = groupby(df_BE, cols_grouping)
df_sum_BE_WORLD = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> (x -> sum(skipmissing(x))), renamecols=false)
df_sum_BE_WORLD.TRADE_TYPE .= "total"

# compute BE_EXTRA and BE_INTRA (BE imports/exports to extra/intra EU countries)
cols_grouping = ["TRADE_TYPE", "PRODUCT_NC", "FLOW", "PERIOD"]
cols_grouping = [cols_grouping; prod_class_other]
gdf = groupby(df_BE, cols_grouping)
df_sum_BE_EXTRA_INTRA = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> (x -> sum(skipmissing(x))), renamecols=false)

# merge tables
df_sum_EU = vcat(df_sum_EU_EXTRA_INTRA, df_sum_EU_WORLD)
df_sum_BE = vcat(df_sum_BE_EXTRA_INTRA, df_sum_BE_WORLD)

# join, rename and compute shares in percentages (?)
cols_join = ["TRADE_TYPE", "PRODUCT_NC", "FLOW", "PERIOD"]
cols_join = [cols_join; prod_class_other]
df_join = leftjoin(df_sum_EU, df_sum_BE, on=cols_join, makeunique=true)
rename!(df_join, ["VALUE_IN_EUROS", "QUANTITY_IN_KG", "VALUE_IN_EUROS_1", "QUANTITY_IN_KG_1"] .=> ["VALUE_EU", "QUANITY_EU", "VALUE_BE", "QUANITY_BE"])
df_join.VALUE_SHARE = df_join.VALUE_BE ./ df_join.VALUE_EU .* 100
df_join.QUANTITY_SHARE = df_join.QUANITY_BE ./ df_join.QUANITY_EU .* 100

# possibly need to sort before grouping, seems that grouping sorts incorrectly sometimes?
sort!(df_join, [:PRODUCT_NC, :FLOW, :TRADE_TYPE, :PERIOD]) # dont understand why I need the sorting here? otherwise some MOM are not taken?

# compute SHARE_CHANGE
cols_grouping = ["TRADE_TYPE", "PRODUCT_NC", "FLOW"]
gdf = groupby(df_join, cols_grouping)
df_share_BE = transform(gdf, [:VALUE_SHARE, :QUANTITY_SHARE] .=> pct_change .=> [:VALUE_SHARE_CHANGE, :QUANTITY_SHARE_CHANGE])

# compute MOM_CHANGE
gdf = groupby(df_share_BE, cols_grouping)
df_share_BE = transform(gdf, [:PERIOD, :VALUE_SHARE] => mom_change => :MOM_VALUE,
                             [:PERIOD, :QUANTITY_SHARE] => mom_change => :MOM_QUANTITY)

# export data
df_export = ifelse.(ismissing.(df_share_BE), NaN, df_share_BE)[1:15_000,:] # excel cannot deal with missing (but still an issue opening the file)
XLSX.writetable(dir_io * "clean/" * "df_share_BE.xlsx", df_export, overwrite=true)

CSV.write(dir_io * "clean/" * "df_share_BE.csv", df_share_BE)

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


text:
struct Suppressed{T}
    item::T
end
Base.show(io::IO, x::Suppressed) = print(io, "<suppressed ", x.item, '>')

function print_var(io::IO, name::Symbol, val)
    print("  | ")
    if val === nothing
        @assert false
    else
        val = something(val)
        T = typeof(val)
        try
            val = repr(val)
            if length(val) > 150
                val = Suppressed("$(length(val)) bytes of output")
            end
        catch
            val = Suppressed("printing error")
        end
        println(io, name, "::", T, " = ", val)
    end
end

print_locdesc(io::IO, frame::JuliaStackFrame) = println(io, locdesc(frame))

function print_locals(io::IO, frame::JuliaStackFrame)
    for i = 1:length(frame.locals)
        if !isa(frame.locals[i], Nothing)
            # #self# is only interesting if it has values inside of it. We already know
            # which function we're in otherwise.
            val = something(frame.locals[i])
            if frame.code.code.slotnames[i] == Symbol("#self#") && (isa(val, Type) || sizeof(val) == 0)
                continue
            end
            print_var(io, frame.code.code.slotnames[i], frame.locals[i])
        end
    end
    if frame.code.scope isa Method
        for (sym, value) in zip(sparam_syms(frame.code.scope), frame.sparams)
            print_var(io, sym, value)
        end
    end
end

function print_frame(io::IO, num::Integer, frame::JuliaStackFrame)
    print(io, "[$num] ")
    print_locdesc(io, frame)
    print_locals(io, frame)
end


text:
# Replication of Conesa et al "Taxing Capital"
using Parameters, QuantEcon, BasisMatrices,LinearAlgebra,Optim,DataFrames,Gadfly,SparseArrays,Arpack,Roots

#### Step 0: Create a struct for the model parameters
@with_kw mutable struct TaxMod
    # Demographics
    jr::Int64 = 46 # Retirement age
    J::Int64 = 81 # Maximum age
    nn::Float64 = 0.011 # Population growth
    ψ::Vector{Float64} = zeros(0) # survival probabilities

    # Preferences
    β::Float64 = 1.00093 # Discount Factor
    σ::Float64 = 4.0 # Risk aversion
    γ::Float64 = 0.377 # Consumption share

    # Labor Productivity Process
    var_α::Float64 = 0.14 # Variance types
    ρ::Float64 = 0.98 # persistence
    var_η::Float64 = 0.0289 # variance shock
    Nη::Int64 = 7 # Number of states
    η::Vector{Float64} = zeros(0) # productivity shock states
    Π::Matrix{Float64} = zeros(0,0) # Transition matrix
    ϵ::Matrix{Float64} = zeros(0,0) # ability x age efficiency matrix

    # Technology
    α::Float64 = 0.36 # capital share
    δ::Float64 = 0.0833 # depreciation
    TFP::Float64 = 1.0 # scale parameter Z
    w::Float64 = 1.0 # wages
    r::Float64 = 0.05 # interest rate

    # Government Policy
    τc::Float64 = 0.05 # consumption tax
    τk::Float64 = 0.36 # capital tax
    κ0::Float64 = 0.258 # marginal tax
    κ1::Float64 = 0.768 # tax progressivity
    κ2::Float64 = 0.0 # deduction
    τp::Float64 = 0.124 # payroll tax
    b::Float64 = 0.5 # Social security replacement rate
    ȳ::Float64 = 87000.0/37748.0

    # Grid sizes
    ns::Int64 = 7 #X Also number of states?
    na::Int64 = 10 # asset grid
    nl::Int64 = 4  # leisure grid
    nty::Int64 = 2
    maxit::Int64 = 10000

    # Other
    a̲::Float64 = 0.0 # borrowing constraint
    umin::Float64 = -1.0E+2 # minimum utility
    penscale::Float64 = 10000000 # penalty on utility for bad c or l
    l̅::Float64 = 0.99 # maximum value for labor supply

    # Spending stuff
    G::Float64 = 6.670
    Tr::Float64 = 2.0
    SS::Float64 = 2.0

    #Solution
    agrid::Vector{Float64} = zeros(0)
    lgrid::Vector{Float64} = zeros(0)
    k::Int = 2 #type of interpolation
    Vf::Array{Interpoland} = Interpoland[] # value function given a
    Vlf::Array{Interpoland} = Interpoland[] # Value function given a,l
    cf::Array{Interpoland} = Interpoland[]

end;

"""
    U(TM::TaxMod, c, l)

Calculates household utility given consumption, leisure, and preferences
"""
function U(TM::TaxMod, c, l)
    @unpack σ,γ,penscale,umin = TM

    if c <= 0.0 # Setting very low utility for negative consumption
		U = umin - penscale * c^2
	elseif l < 0.0 # Setting very low utility for negative labor supply
		U = umin - penscale * l^2
	elseif l >= 1.0 # Setting very low utility for labor supply over 1
		U = umin - penscale*(l-1.0)^2
	else # Real utility function
		U = (1/(1-σ)) * (((c^γ)*((1-l)^(1-γ)))^(1-σ))
	end

    return U
end

"""
    marginal_utility(TM::TaxMod, c, l)

Calculates marginal utility
"""
function marginal_utility(TM::TaxMod, c, l)
    @unpack σ,γ = TM

    if c > 0.0
        mu = (γ*(c^γ*(1.0-l)^(1.0-γ))^(1.0-σ))/c
	else # Very large MU if c is zero or negative
		mu = 1000000.0+ abs(c)^2.0
	end

end

"""
    tax_gs(y, k0, k1, k2)

Gouveia and Strauss tax function
"""
function tax_gs(y, k0, k1, k2)

    tax_gs = 0.0

    if k1 >= 0.001 && y > 0.0
        tax_gs = k0*(y - (y^(-k1)+k2)^(-1/k1))
    elseif k1 < 0.001 && y > 0.0 #
        tax_gs = k0*y + k2
    else # No taxes if income is negative/zero
        tax_gs = 0.0
    end

    return tax_gs
end

"""
    Tax3()

This tax function is separable in capital and labor income
"""

function Tax3(l0,l1,l2,ca0,lear,caear)

    Tax3=0.0

    # Labor Income Tax
    if l1 >= 0.000001 && lear >= 0.00000001
        Tax3 = l0*( lear - ( lear^(-l1)+l2 )^(-1.0/l1) )
    elseif l1 < 0.000001 && lear >= 0.00000001
        Tax3 = l0*lear+l2
    else
        Tax3 = 0.0
    end

    # Capital Income Tax
    Tax3 = Tax3 + ca0*caear

end

"""
    marginal_tax_gs(k0,k1,k2,cap_earn)
"""
function marginal_tax_gs(k0,k1,k2,cap_earn)

    mt = 0.0
    if k1 >= 0.001 && cap_earn > 0.0
        mt = k0*(1 - (1 + k2*cap_earn^k1)^(-1/k1 - 1))
    elseif k1 < 0.001 && cap_earn > 0.0
        mt = k0
    end

    return mt
end

"""
    marginal_tax3(ca0)
"""
function marginal_tax3(ca0)
    return ca0
end

### SEE PARAMS for basefun subroutine
tm = TaxMod()

##### Step 1: Create a grid of individual asset holdings. define survival rates, age distribution, and labor efficiency grids

"""
    l(c,j)

Static FOC relating consumption and leisure
"""
function l(TM::TaxMod, c)
    @unpack jr,γ,τc,ϵ,w,η = TM

    for s in 1:ns
        for ty in 1:nty
            for j in 1:J
                if j >= jr
                    l[s,ty,j] = 0
                else
  #                  l[s,ty,j] = 1 - ((1-γ)/γ)*(1 + τc)/(w*)

                cf[s]= Interpoland(abasis,c)
            end
        end
    end

end

"""
    setup_grid!()

This creates the grids for assets and labor. Also loads demographic
characteristics: age efficiency units from Hansen 1993 and population
from Bell and Miller 2002.
"""
a′grid = zeros(tm.na)
lgrid = zeros(tm.nl)

function setup_grid!(TM::TaxMod, a̅ = 75, curve = 2.5)
    @unpack na,nl,a̲,l̅,J,nty,nn,ns,ρ,var_η,var_α,Vf,cf,k = tm

    # Asset grid: with curvature
    tm.agrid = (a̅ - a̲).*LinRange(0,1,na).^curve .+ a̲
    # Labor grid: evenly spaced
    lgrid = LinRange(0,l̅,nl)

    # Loading demographic data from other script
    include("demographics.jl")
    #ephansen,pop,surv,Nu,mu,ep,measty,topop = age_eff_pop(J,nty,nn)
    TM.ψ = age_eff_pop(J,nty,nn).surv
    ephansen = age_eff_pop(J,nty,nn).ephansen

    # Labor Productivity Markov Chain
    mc = rouwenhorst(ns,ρ,var_η)
    TM.Π = mc.p
    TM.η = exp.(mc.state_values)

    # Fixed effects
    TM.ϵ = zeros(nty,J)
    TM.ϵ[1,1:J] = exp(-sqrt(var_α)) .* ephansen[1:J]
    TM.ϵ[2,1:J] = exp(sqrt(var_α)) .* ephansen[1:J]

    #First guess of interpolation functions
    abasis = Basis(SplineParams(a′grid,0,k))
    a = nodes(abasis)[1]

    Vf = TM.Vf = Array{Interpoland}(undef,ns,nty,J)
    cf = TM.cf = Array{Interpoland}(undef,ns,nty,J)
    lf = Array{Interpoland}(undef,ns,nty,J)

    for s in 1:ns
        for
            c = @. r̄*a + w̄*HH.ϵ[s]
            V = U(HH,c,l)./(1-β)
            l =

            Vf[s]= Interpoland(abasis,V)
            cf[s]= Interpoland(abasis,c)

        end
    end

    for s in 1:ns
        for ty in 1:nty
            for j in 1:J

                c[:,s,ty,j]= @. r*a + w*TM.ϵ[s]
                V = U.(c[:,s,ty,j],0.5)./(1-β)

                Vf[s]= Interpoland(abasis,V)
                cf[s]= Interpoland(abasis,c)
            end
        end
    end

end

setup_grid!(tm)
##### Step 2: Solving household problem

# Start with grid of next period's assets
# Guess consumption tommorrow
# Use to get consumption today w/ EE
# Use that to get labor using FOC
# Use c, labor, a′ to get assets today
#

"""
    iterate_endogenousgrid(HH,a′grid,cf′)

Iterates on Euler equation using endogenous grid method

c needs to be na x ns

"""
function iterate_endogenousgrid(tm::TaxMod,a′grid,cf′)
    @unpack ns,nty,J,jr,ψ,β,τk,τc,σ,κ0,κ1,κ2,w,ϵ,η,Π = tm

    # Finding initial guess for c′ at a′grid values for each ability type and period
    c′ = zeros(length(a′grid),ns,nty,J)
    for s in 1:ns
        for ty in 1:nty
            for j in 1:J
                c′[:,s,ty,j]= cf′[s,ty,j](a′grid)
            end
        end
    end

    c = zeros(length(a′grid),ns,nty,J)
    # Finding c from c′ using Euler Equation
    for i in 1:length(a′grid)
        for s in 1:ns
            for ty in 1:nty
                for j in 1:(J-1)
                    if j == J
                        c[i,s,ty,j] = 0
                    elseif j >= jr
                        c[i,s,ty,j] = 0
                    elseif j < jr
                        c[i,s,ty,j] = ((β * ψ[j] * ((1+r(1-τk))/(1+τc)) * c′[i,s,ty,j]^(-σ) *
                            ((1 - marginal_tax_gs(κ0,κ1,κ2,w*ϵ[ty,j]) - τss)*ϵ[ty,j]*η[s])/
                            ((1 - marginal_tax_gs(κ0,κ1,κ2,w*ϵ[ty,j+1]) - τss)*ϵ[ty,j+1]) ./ η) * Π[s,:])^(1/(σ))
                    end
                end
            end
        end

                    (1+r̄)*(c′).^(-γ)*Π' #RHS of Euler Equation
    c = EERHS.^(-1/γ)

    #compute implies assets
    a = ((c .+ a′grid) .- w̄ .*ϵ')./(1+r̄)

    cf = Vector{Interpoland}(undef,Nϵ)
    for s in 1:Nϵ
        if a[1,s]> a̲
            c̲ = r̄*a̲ + w̄*ϵ[s]
            cf[s]= Interpoland(Basis(SplineParams([a̲; a[:,s]],0,1)),[c̲;c[:,s]])
        else
            cf[s]= Interpoland(Basis(SplineParams(a[:,s],0,1)),c[:,s])
        end
    end
    return cf
end;


text:
module TestSimpleStateEvents

using TinyModia
using DifferentialEquations
@usingModiaPlot

SimpleStateEvents = Model(
    fmax = 1.5,
    m    = 1.0,
    k    = 1.0,
    d    = 0.1,
    s    = Var(init = 2.0),
    v    = Var(init = 0.0),
    equations = :[
        sPos = positive(s)
        f = if sPos; 0.0 else fmax end   # or: f = sPos ? 0.0 : fmax
        v = der(s)
        m*der(v) + d*v + k*s = f
    ]
)

model = @instantiateModel(SimpleStateEvents)

simulate!(model, Tsit5(), stopTime = 10, log=false, logEvents=false, requiredFinalStates = [1.1756992201144405, -0.5351307571761175])

plot(model, ["s", "v", "sPos", "f"])

end


text:
## ArbVector basic definitions
function clear!(v::arb_vec_struct)
    ccall(@libarb(_arb_vec_clear), Cvoid, (Ptr{arb_struct}, Clong), v.entries, v.n)
end
Base.cconvert(::Type{Ptr{arb_struct}}, v::Union{ArbVector,ArbRefVector}) = v.arb_vec
Base.unsafe_convert(::Type{Ptr{arb_struct}}, v::arb_vec_struct) = v.entries

# AbstractVector interface for ArbVector
Base.size(v::arb_vec_struct) = (v.n,)
Base.size(v::Union{ArbVector,ArbRefVector}) = size(v.arb_vec)

Base.getindex(v::arb_vec_struct, i::Integer) = v.entries + (i - 1) * sizeof(arb_struct)
Base.@propagate_inbounds function Base.getindex(v::ArbVector, i::Integer)
    @boundscheck checkbounds(v, i)
    x = Arb(prec = precision(v))
    x[] = v.arb_vec[i]
    x
end
Base.@propagate_inbounds function Base.getindex(v::ArbRefVector, i::Integer)
    @boundscheck checkbounds(v, i)
    return ArbRef(v.arb_vec[i], precision(v), cstruct(v))
end

"""
    ref(v::ArbVector, i)

Similar to `v[i]` but instead of an `Arb` returns an `ArbRef` which still shares the
memory with the `i`-th entry of `v`.
"""
Base.@propagate_inbounds function ref(v::Union{ArbVector,ArbRefVector}, i::Integer)
    @boundscheck checkbounds(v, i)
    return ArbRef(v.arb_vec[i], precision(v), cstruct(v))
end

Base.setindex!(v::arb_vec_struct, x, i::Integer) = (set!(v[i], x); x)
Base.@propagate_inbounds function Base.setindex!(
    v::Union{ArbVector,ArbRefVector},
    x,
    i::Integer,
)
    @boundscheck checkbounds(v, i)
    ref(v, i)[] = x
    return x
end

## AcbVector basic definitions
function clear!(v::acb_vec_struct)
    ccall(@libarb(_acb_vec_clear), Cvoid, (Ptr{acb_struct}, Clong), v.entries, v.n)
end
Base.cconvert(::Type{Ptr{acb_struct}}, v::Union{AcbVector,AcbRefVector}) = v.acb_vec
Base.unsafe_convert(::Type{Ptr{acb_struct}}, v::acb_vec_struct) = v.entries

# AbstractVector interface for ACbVector
Base.size(v::acb_vec_struct) = (v.n,)
Base.size(v::Union{AcbVector,AcbRefVector}) = size(v.acb_vec)

Base.getindex(v::acb_vec_struct, i::Integer) = v.entries + (i - 1) * sizeof(acb_struct)
Base.@propagate_inbounds function Base.getindex(v::AcbVector, i::Integer)
    @boundscheck checkbounds(v, i)
    x = Acb(prec = precision(v))
    x[] = v.acb_vec[i]
    x
end
Base.@propagate_inbounds function Base.getindex(v::AcbRefVector, i::Integer)
    @boundscheck checkbounds(v, i)
    return AcbRef(v.acb_vec[i], precision(v), cstruct(v))
end

"""
    ref(v::AcbVector, i)

Similar to `v[i]` but instead of an `Acb` returns an `AcbRef` which still shares the
memory with the `i`-th entry of `v`.
"""
Base.@propagate_inbounds function ref(v::Union{AcbVector,AcbRefVector}, i::Integer)
    @boundscheck checkbounds(v, i)
    return AcbRef(v.acb_vec[i], precision(v), cstruct(v))
end

Base.setindex!(v::acb_vec_struct, x, i::Integer) = (set!(v[i], x); x)
Base.@propagate_inbounds function Base.setindex!(
    v::Union{AcbVector,AcbRefVector},
    x,
    i::Integer,
)
    @boundscheck checkbounds(v, i)
    ref(v, i)[] = x
    return x
end

## Common methods
const Vectors = Union{ArbVector,ArbRefVector,AcbVector,AcbRefVector}

# General constructor
for T in [:ArbVector, :ArbRefVector, :AcbVector, :AcbRefVector]
    @eval function $T(v::AbstractVector; prec::Integer = _precision(first(v)))
        V = $T(length(v); prec = prec)
        @inbounds for (i, vᵢ) in enumerate(v)
            V[i] = vᵢ
        end
        return V
    end
end

# Arithmetic
for (jf, af) in [(:+, :add!), (:-, :sub!)]
    @eval function Base.$jf(
        v::T,
        w::T,
    ) where {T<:Union{ArbVector,ArbRefVector,AcbVector,AcbRefVector}}
        @boundscheck (length(v) == length(w) || throw(DimensionMismatch()))
        u = T(length(v); prec = max(precision(v), precision(w)))
        $af(u, v, w)
        u
    end
end
Base.:(-)(v::Vectors) = neg!(similar(v), v)

Base.copy(v::Vectors) = copy!(similar(v), v)
Base.copy!(w::T, v::T) where {T<:Vectors} = set!(w, v)
Base.copyto!(w::T, v::T) where {T<:Vectors} = set!(w, v)