text:
@enum CSVType singledata pairdata assocdata groupdata

"""
    getfileextension(filepath)

A quick helper to get the file extension of any given path (without the dot).

# Examples
```julia-repl
julia> getfileextension("~/Desktop/text.txt")
"txt"
```
"""

function getfileextension(filepath::AbstractString)
    dotpos = findlast(isequal('.'), filepath)
    isnothing(dotpos) && return ""
    return filepath[dotpos+1:end]
end

"""
    getpaths(location; relativetodatabase=false)

Returns database paths that is optionally relative to Clapeyron.jl directory.
If path is a file, then return an Array containing a single path to that file.
If path is a directory, then return an Array containing paths to all csv files in that directory.

# Examples
```julia-repl
julia> getpaths("SAFT/PCSAFT"; relativetodatabase=true)
3-element Array{String,1}:
 "/home/user/.julia/packages/Clapeyron.jl/xxxxx/database/SAFT/PCSAFT/data_PCSAFT_assoc.csv"
 "/home/user/.julia/packages/Clapeyron.jl/xxxxx/database/SAFT/PCSAFT/data_PCSAFT_like.csv"
 "/home/user/.julia/packages/Clapeyron.jl/xxxxx/database/SAFT/PCSAFT/data_PCSAFT_unlike.csv"

```
"""
function getpaths(location::AbstractString; relativetodatabase::Bool=false)
    # We do not use realpath here directly because we want to make the .csv suffix optional.
    filepath = relativetodatabase ? normpath(dirname(pathof(Clapeyron)), "..", "database", location) : location
    isfile(filepath) && return [realpath(filepath)]
    isfile(filepath * ".csv") && return [realpath(filepath * ".csv")]
    if !isdir(filepath)
        relativetodatabase ? error("The path ", location, " does not exist in the Clapeyron database.") :
            error("The path ", location, " does not exist.")
    end
    files = joinpath.(filepath, readdir(filepath))
    typeof(files)
    return realpath.(files[isfile.(files) .& (getfileextension.(files) .== "csv")])
end

export getparams
function getparams(components::Array{String,1}, locations::Array{String,1}=String[]; userlocations::Array{String,1}=String[], asymmetricparams::Array{String,1}=String[], ignore_missingsingleparams::Bool=false, verbose::Bool=false)
    # Gets all parameters from database.
    # locations is a list of paths relative to the Clapeyron database directory.
    # userlocations is a list of paths input by the user.
    # If parameters exist in multiple files, Clapeyron gives priority to files in later paths.
    # asymmetricparams is a list of parameters for which matrix reflection is disabled.
    # ignore_missingsingleparams gives users the option to disable component existence check in single params.
    filepaths = string.(vcat([(getpaths.(locations; relativetodatabase=true)...)...], [(getpaths.(userlocations)...)...]))
    allcomponentsites = findsitesincsvs(components, filepaths; verbose=verbose)
    allparams, paramsourcecsvs, paramsources = createparamarrays(components, filepaths, allcomponentsites; verbose=verbose)
    return packageparams(allparams, components, allcomponentsites, paramsourcecsvs, paramsources; asymmetricparams=asymmetricparams, ignore_missingsingleparams=ignore_missingsingleparams)
end

function getparams(groups::GroupParam, locations::Array{String,1}=String[]; userlocations::Array{String,1}=String[], asymmetricparams::Array{String,1}=String[], ignore_missingsingleparams::Bool=false, verbose::Bool=false)
    # For GC.
    return getparams(groups.flattenedgroups, locations; userlocations=userlocations, asymmetricparams=asymmetricparams, ignore_missingsingleparams=ignore_missingsingleparams, verbose=verbose)
end

function getparams(components::String, locations::Array{String,1}=String[]; userlocations::Array{String,1}=String[], asymmetricparams::Array{String,1}=String[], ignore_missingsingleparams::Bool=false, verbose::Bool=false)
    return getparams([components], locations; userlocations=userlocations, asymmetricparams=asymmetricparams, ignore_missingsingleparams=ignore_missingsingleparams, verbose=verbose)
end

function packageparams(allparams::Dict, components::Array{String,1}, allcomponentsites::Array{Array{String,1},1}, paramsourcecsvs::Dict{String,Set{String}}, paramsources::Dict{String,Set{String}}; asymmetricparams::Array{String,1}=String[], ignore_missingsingleparams::Bool=false)
    # Package params into their respective Structs.
    output = Dict{String,ClapeyronParam}()
    for (param, value) ∈ allparams
        if value isa Array{<:Any,1}
            newvalue, ismissingvalues = defaultmissing(value)
            if !ignore_missingsingleparams && any(ismissingvalues)
                error("Missing values exist in single parameter ", param, ": ", value, ".")
            end
            output[param] = SingleParam(param, newvalue, ismissingvalues, components, allcomponentsites, collect(paramsourcecsvs[param]), collect(paramsources[param]))
        elseif value isa Array{<:Array,2}
            param ∉ asymmetricparams && mirrormatrix!(value)
            newvalue_ismissingvalues = defaultmissing.(value)
            newvalue = getindex.(newvalue_ismissingvalues, 1)
            ismissingvalues = getindex.(newvalue_ismissingvalues, 2)
            output[param] = AssocParam(param, newvalue, ismissingvalues, components, allcomponentsites, collect(paramsourcecsvs[param]), collect(paramsources[param]))
        elseif value isa Array{<:Any, 2}
            param ∉ asymmetricparams && mirrormatrix!(value)
            newvalue, ismissingvalues = defaultmissing(value)
            if (!ignore_missingsingleparams
                && !all([ismissingvalues[x,x] for x ∈ 1:size(ismissingvalues,1)])
                && any([ismissingvalues[x,x] for x ∈ 1:size(ismissingvalues,1)]))
                error("Partial missing values exist in diagonal of pair parameter ", param, ": ", [value[x,x] for x ∈ 1:size(ismissingvalues,1)], ".")
            end
            output[param] = PairParam(param, newvalue, ismissingvalues, components, allcomponentsites, collect(paramsourcecsvs[param]), collect(paramsources[param]))
        else
            error("Format for ", param, " is incorrect.")
        end
    end
    return output
end

function createparamarrays(components::Array{String,1}, filepaths::Array{String,1}, allcomponentsites::Array{Array{String,1},1}; verbose::Bool=false)
    # Returns Dict with all parameters in their respective arrays.
    checkfor_clashingheaders(filepaths)
    allparams = Dict{String,Any}()
    paramsourcecsvs = Dict{String,Set{String}}()
    paramsources = Dict{String,Set{String}}()
    # Read the filepaths in reverse in order to ensure that unused sources do not get added.
    for filepath ∈ reverse(filepaths)
        csvtype = readcsvtype(filepath)
        if csvtype == groupdata
            verbose && println("Skipping groupdata csv ", filepath)
            continue
        end
        headerparams = readheaderparams(filepath)
        verbose && println("Searching for ", string(csvtype), " headers ", headerparams, " for components ", components, " at ", filepath, "...")
        foundparams, paramtypes, sources = findparamsincsv(components, filepath, headerparams; verbose=verbose)
        foundcomponents = collect(keys(foundparams))
        foundparams = swapdictorder(foundparams)
        for headerparam ∈ headerparams
            if !haskey(allparams, headerparam)
                if ismissing(paramtypes[headerparam])
                    allparams[headerparam] = createemptyparamsarray(Any, csvtype, components, allcomponentsites)
                else
                    allparams[headerparam] = createemptyparamsarray(paramtypes[headerparam], csvtype, components, allcomponentsites)
                end
            end
            if !haskey(paramsourcecsvs, headerparam)
                paramsourcecsvs[headerparam] = Set{String}()
                paramsources[headerparam] = Set{String}()
            end
            if csvtype == singledata
                isempty(foundparams) && continue
                for (component, value) ∈ foundparams[headerparam]
                    currenttype = nonmissingtype(eltype(allparams[headerparam]))
                    if paramtypes[headerparam] != currenttype
                        if currenttype <: Number && paramtypes[headerparam] <: String || currenttype <: String && paramtypes[headerparam] <: Number
                            allparams[headerparam] = convert(Array{Any}, allparams[headerparam])
                        else
                            try
                                allparams[headerparam] = convert(Array{Union{Missing,paramtypes[headerparam]}}, allparams[headerparam])
                            catch e
                            end
                        end
                    end
                    idx = findfirst(isequal(component), components)
                    if allparams[headerparam] isa Array{<:Any,2}
                        !ismissing(allparams[headerparam][idx,idx]) && continue
                        allparams[headerparam][idx,idx] = value
                    else
                        !ismissing(allparams[headerparam][idx]) && continue
                        allparams[headerparam][idx] = value
                    end
                    push!(paramsourcecsvs[headerparam], filepath)
                    !ismissing(sources[component]) && push!(paramsources[headerparam], sources[component])
                end
            end
            if csvtype == pairdata
                if allparams[headerparam] isa Array{<:Any,1}
                    allparams[headerparam] = convertsingletopair(allparams[headerparam], outputmissing=true)
                end
                isempty(foundparams) && continue
                for (componentpair, value) ∈ foundparams[headerparam]
                    currenttype = nonmissingtype(eltype(allparams[headerparam]))
                    if paramtypes[headerparam] != currenttype
                        if currenttype <: Number && paramtypes[headerparam] <: String || currenttype <: String && paramtypes[headerparam] <: Number
                            allparams[headerparam] = convert(Array{Any}, allparams[headerparam])
                        else
                            try
                                allparams[headerparam] = convert(Array{Union{Missing,paramtypes[headerparam]}}, allparams[headerparam])
                            catch e
                            end
                        end
                    end
                    idx1 = findfirst(isequal(componentpair[1]), components)
                    idx2 = findfirst(isequal(componentpair[2]), components)
                    !ismissing(allparams[headerparam][idx1,idx2]) && continue
                    allparams[headerparam][idx1,idx2] = value
                    push!(paramsourcecsvs[headerparam], filepath)
                    !ismissing(sources[componentpair]) && push!(paramsources[headerparam], sources[componentpair])
                end
            end
            if csvtype == assocdata
                isempty(foundparams) && continue
                for (assocpair, value) ∈ foundparams[headerparam]
                    currenttype = nonmissingtype(eltype(first(allparams[headerparam])))
                    if currenttype != paramtypes[headerparam]
                        if currenttype <: Number && paramtypes[headerparam] <: String || currenttype <: String && paramtypes[headerparam] <: Number
                            allparams[headerparam] = convert(Array{Array{Any,2},2}, allparams[headerparam])
                        else
                            try
                                allparams[headerparam] = convert(Array{Array{Union{Missing,paramtypes[headerparam]},2},2}, allparams[headerparam])
                            catch e
                            end
                        end
                    end
                    idx1 = findfirst(isequal(assocpair[1][1]), components)
                    idx2 = findfirst(isequal(assocpair[1][2]), components)
                    idx21 = findfirst(isequal(assocpair[2][1]), allcomponentsites[idx1])
                    idx22 = findfirst(isequal(assocpair[2][2]), allcomponentsites[idx2])
                    !ismissing(allparams[headerparam][idx1,idx2][idx21,idx22]) && continue
                    allparams[headerparam][idx1,idx2][idx21,idx22] = value
                    push!(paramsourcecsvs[headerparam], filepath)
                    !ismissing(sources[assocpair]) && push!(paramsources[headerparam], sources[assocpair])
                end
            end
        end
    end
    return allparams, paramsourcecsvs, paramsources
end

function defaultmissing(array::Array; defaultvalue=nothing)
    # Returns a non-missing Array with missing values replaced by default values.
    # Also returns an Array of Bools to retain overwritten information.
    arraycopy = deepcopy(array)
    type = nonmissingtype(eltype(array))
    ismissingvalues = [ismissing(ai) for ai in arraycopy]
    if type != Any && type === eltype(array)
        return arraycopy, ismissingvalues
    end
    if isnothing(defaultvalue)
        if type === Union{} # If it only contains missing
            # Default to 0. unless specified with defaultvalue parameter
            # but make type Any.
            arraycopy = convert(Array{Any}, arraycopy)
            arraycopy[ismissingvalues] .= 0.
            type = Any
        elseif type === Any
            # This means that it has a mix of Number and String types.
            # If that was the case, convert it all to String.
            arraycopy[ismissingvalues] .= ""
            arraycopy = string.(arraycopy)
            type = String
        elseif type <: Number
            arraycopy[ismissingvalues] .= 0
        elseif type <: AbstractString
            arraycopy[ismissingvalues] .= ""
        else
            error("Unsupported type.")
        end
    else
        arraycopy = convert(Array{Any}, arraycopy)
        arraycopy[ismissingvalues] .= defaultvalue
        type = typeof(defaultvalue)
    end
    return convert(Array{type}, arraycopy), ismissingvalues
end

function swapdictorder(dict::Dict)
    # Swap the first two levels in a nested dictionary.
    isempty(dict) && return dict
    output = Dict()
    outerkeys = keys(dict)
    try
        innerkeys = keys(first(values(dict)))
        for innerkey ∈ innerkeys, outerkey ∈ outerkeys
            if !haskey(output, innerkey)
                output[innerkey] = Dict{Any,Any}(outerkey => dict[outerkey][innerkey])
            end
            push!(output[innerkey], outerkey => dict[outerkey][innerkey])
        end
    catch e
        error("The format of the nested dictionary is not correct.")
    end
    return output
end

function findparamsincsv(components::Array{String,1}, filepath::AbstractString, headerparams::Array{String,1}; columnreference::AbstractString="species", sitecolumnreference::AbstractString="site", sourcecolumnreference::AbstractString="source", verbose::Bool=false, normalisecomponents::Bool=true)
    # Returns a Dict with all matches in a particular file for one parameter.
    normalised_columnreference = normalisestring(columnreference)
    csvtype = readcsvtype(filepath)
    df = CSV.File(filepath; header=3, pool=0,lazystrings=true)
    csvheaders = String.(Tables.columnnames(df))
    normalised_components = normalisestring.(components; isactivated=normalisecomponents)
    normalised_csvheaders = normalisestring.(csvheaders)
    normalised_headerparams = normalisestring.(headerparams)
    normalised_headerparams ⊈ normalised_csvheaders && error("Headers ", setdiff(normalised_headerparams, normalised_csvheaders), " not present in csv header.")

    foundvalues = Dict()
    paramtypes = Dict(headerparams .=> [Tables.columntype(df, Symbol(x)) for x ∈ headerparams])
    sources = Dict()

    normalised_sourcecolumnreference = normalisestring(sourcecolumnreference)
    getsources = false
    if normalised_sourcecolumnreference ∈ normalised_csvheaders
        getsources = true
        sourcecolumn = Symbol(normalised_csvheaders[findfirst(isequal(normalised_sourcecolumnreference), normalised_csvheaders)])
    end

    if csvtype == singledata
        lookupcolumnindex = findfirst(isequal(normalised_columnreference), normalised_csvheaders)
        isnothing(lookupcolumnindex) && error("Header ", normalised_columnreference, " not found.")
        lookupcolumn = Symbol(csvheaders[lookupcolumnindex])
        for row ∈ Tables.rows(df)
            component = row[lookupcolumn]
            foundcomponentidx = findfirst(isequal(normalisestring(component; isactivated=normalisecomponents)), normalised_components)
            isnothing(foundcomponentidx) && continue
            verbose && print("Found component: ", component)
            component = components[foundcomponentidx]
            foundvalues[component] = Dict()
            for headerparam ∈ headerparams
                foundvalues[component][headerparam] = row[Symbol(headerparam)]
            end
            verbose && println(" with values ", foundvalues[component])
            if getsources
                source = row[sourcecolumn]
                sources[component] = source
            else
                sources[component] = missing
            end
        end
    elseif csvtype == pairdata
        lookupcolumnindex1 = findfirst(x -> x[1:end-1] == normalised_columnreference && x[end] == '1', normalised_csvheaders)
        lookupcolumnindex2 = findfirst(x -> x[1:end-1] == normalised_columnreference && x[end] == '2', normalised_csvheaders)
        isnothing(lookupcolumnindex1) && error("Header ", normalised_columnreference * "1", " not found.")
        isnothing(lookupcolumnindex2) && error("Header ", normalised_columnreference * "2", " not found.")
        lookupcolumn1 = Symbol(csvheaders[lookupcolumnindex1])
        lookupcolumn2 = Symbol(csvheaders[lookupcolumnindex2])
        for row ∈ Tables.rows(df)
            component1 = row[lookupcolumn1]
            component2 = row[lookupcolumn2]
            foundcomponentidx1 = findfirst(isequal(normalisestring(component1; isactivated=normalisecomponents)), normalised_components)
            foundcomponentidx2 = findfirst(isequal(normalisestring(component2; isactivated=normalisecomponents)), normalised_components)
            (isnothing(foundcomponentidx1) || isnothing(foundcomponentidx2)) && continue
            verbose && print("Found component pair: ", (component1, component2))
            componentpair = (components[foundcomponentidx1], components[foundcomponentidx2])
            foundvalues[componentpair] = Dict()
            for headerparam ∈ headerparams
                foundvalues[componentpair][headerparam] = row[Symbol(headerparam)]
            end
            verbose && println(" with values ", foundvalues[componentpair])
            if getsources
                source = row[sourcecolumn]
                sources[componentpair] = source
            else
                sources[componentpair] = missing
            end
        end
    elseif csvtype == assocdata
        lookupcolumnindex1 = findfirst(x -> x[1:end-1] == normalised_columnreference && x[end] == '1', normalised_csvheaders)
        lookupcolumnindex2 = findfirst(x -> x[1:end-1] == normalised_columnreference && x[end] == '2', normalised_csvheaders)
        isnothing(lookupcolumnindex1) && error("Header ", normalised_columnreference * "1", " not found.")
        isnothing(lookupcolumnindex2) && error("Header ", normalised_columnreference * "2", " not found.")
        lookupcolumn1 = Symbol(csvheaders[lookupcolumnindex1])
        lookupcolumn2 = Symbol(csvheaders[lookupcolumnindex2])

        normalised_sitecolumnreference = normalisestring(sitecolumnreference)
        lookupsitecolumnindex1 = findfirst(x -> x[1:end-1] == normalised_sitecolumnreference && x[end] == '1', normalised_csvheaders)
        lookupsitecolumnindex2 = findfirst(x -> x[1:end-1] == normalised_sitecolumnreference && x[end] == '2', normalised_csvheaders)
        isnothing(lookupcolumnindex1) && error("Header ", normalised_sitecolumnreference * "1", " not found.")
        isnothing(lookupcolumnindex2) && error("Header ", normalised_sitecolumnreference * "2", " not found.")
        lookupsitecolumn1 = Symbol(csvheaders[lookupsitecolumnindex1])
        lookupsitecolumn2 = Symbol(csvheaders[lookupsitecolumnindex2])

        for row ∈ Tables.rows(df)
            component1 = row[lookupcolumn1]
            component2 = row[lookupcolumn2]
            foundcomponentidx1 = findfirst(isequal(normalisestring(component1; isactivated=normalisecomponents)), normalised_components)
            foundcomponentidx2 = findfirst(isequal(normalisestring(component2; isactivated=normalisecomponents)), normalised_components)
            (isnothing(foundcomponentidx1) || isnothing(foundcomponentidx2)) && continue
            site1 = row[lookupsitecolumn1]
            site2 = row[lookupsitecolumn2]
            verbose && print("Found assoc pair: ", ((component1, component2),(site1, site2)))
            assocpair = ((components[foundcomponentidx1], components[foundcomponentidx2]), (site1, site2))
            foundvalues[assocpair] = Dict()
            for headerparam ∈ headerparams
                foundvalues[assocpair][headerparam] = row[Symbol(headerparam)]
            end
            verbose && println(" with values ", foundvalues[assocpair])
            if getsources
                source = row[sourcecolumn]
                sources[assocpair] = source
            else
                sources[assocpair] = missing
            end
        end
    else
        error("File is of type ", String(csvtype), " and cannot be read with this function.")
    end
    return foundvalues, paramtypes, sources
end

function normalisestring(str::AbstractString; isactivated::Bool=true, tofilter::Regex=r"[ \-\_]", changecase::Bool=true)
    !isactivated && return str
    normalisedstring = replace(str, tofilter => "")
    changecase && (normalisedstring = lowercase(normalisedstring))
    return normalisedstring
end

function findgroupsincsv(components::Array{String,1}, filepath::AbstractString; columnreference::AbstractString="species", groupcolumnreference::AbstractString="groups", verbose::Bool=false, normalisecomponents=true)
    # Returns a Dict with the group string that will be parsed in buildspecies.
    csvtype = readcsvtype(filepath)
    csvtype != groupdata && return Dict{String,String}()
    normalised_components = normalisestring.(components; isactivated=normalisecomponents)
    normalised_columnreference = normalisestring(columnreference)
    normalised_groupcolumnreference = normalisestring(groupcolumnreference)
    csvtype = readcsvtype(filepath)
    df = CSV.File(filepath; header=3,lazystrings=true)
    csvheaders = String.(Tables.columnnames(df))
    normalised_csvheaders = normalisestring.(csvheaders)

    normalised_columnreference ∉ normalised_csvheaders && error("Header ", normalised_columnreference, " not found.")
    normalised_groupcolumnreference ∉ normalised_csvheaders && error("Header ", normalised_groupcolumnreference, " not found.")

    foundgroups = Dict{String,String}()

    lookupcolumnindex = findfirst(isequal(normalised_columnreference), normalised_csvheaders)
    isnothing(lookupcolumnindex) && error("Header ", normalised_columnreference, " not found.")
    lookupcolumn = Symbol(csvheaders[lookupcolumnindex])

    lookupgroupcolumnindex = findfirst(isequal(normalised_groupcolumnreference), normalised_csvheaders)
    isnothing(lookupgroupcolumnindex) && error("Header ", normalised_groupcolumnreference, " not found.")
    lookupgroupcolumn = Symbol(csvheaders[lookupgroupcolumnindex])
    for row ∈ Tables.rows(df)
        component = row[lookupcolumn]
        foundcomponentidx = findfirst(isequal(normalisestring(component; isactivated=normalisecomponents)), normalised_components)
        isnothing(foundcomponentidx) && continue
        verbose && print("Found component: ", component)
        component = components[foundcomponentidx]
        foundgroups[component] = row[lookupgroupcolumn]
        verbose && println(" with groups ", foundgroups[component])
    end
    return foundgroups
end

function readcsvtype(filepath::AbstractString)
    # Searches for type from second line of CSV.
    keywords = ["like", "single", "unlike", "pair", "assoc", "group", "groups"]
    words = split(lowercase(rstrip(getline(filepath, 2), ',')), ' ')
    foundkeywords = intersect(words, keywords)
    isempty(foundkeywords) && error("Unable to determine type of database", filepath, ". Check that keyword is present on Line 2.")
    length(foundkeywords) > 1 && error("Multiple keywords found in database ", filepath, ": ", foundkeywords)
    first(foundkeywords) == "single" && return singledata
    first(foundkeywords) == "like" && return singledata
    first(foundkeywords) == "pair" && return pairdata
    first(foundkeywords) == "unlike" && return pairdata
    first(foundkeywords) == "assoc" && return assocdata
    first(foundkeywords) == "group" && return groupdata
    first(foundkeywords) == "groups" && return groupdata
end

function getline(filepath::AbstractString, selectedline::Int)
    # Simple function to return text from filepath at selectedline.
    open(filepath) do file
        linecount = 1
        for line ∈ eachline(file)
            linecount == selectedline && return line
            linecount += 1
        end
        error("Selected line number exceeds number of lines in file")
    end
end

function readheaderparams(filepath::AbstractString; headerline::Int = 3)
    # Returns array of filtered header strings at line 3.
    headers = split(getline(filepath, headerline), ',')
    ignorelist = ["source", "species", "dipprnumber", "smiles", "site"]
    return String.(filter(x -> normalisestring(x; tofilter=r"[ \-\_\d]") ∉ ignorelist, headers))
end

function checkfor_clashingheaders(filepaths::Array{String,1})
    # Raises an error if the header of any assoc parameter clashes with a non-assoc parameter.
    headerparams = []
    headerparams_assoc = []
    for filepath in filepaths
        csvtype = readcsvtype(filepath)
        if csvtype == singledata || csvtype == pairdata
            append!(headerparams, readheaderparams(filepath))
        elseif csvtype == assocdata
            append!(headerparams_assoc, readheaderparams(filepath))
        end
    end
    clashingheaders = intersect(headerparams, headerparams_assoc)
    !isempty(clashingheaders) && error("Headers ", clashingheaders, " appear in both loaded assoc and non-assoc files.")
end

function findsitesincsvs(components::Array{String,1}, filepaths::Array{String,1}; columnreference::AbstractString="species", sitecolumnreference::AbstractString="site", verbose::Bool=false, normalisecomponents::Bool=true)
    # Look for all relevant sites in the database.
    # Note that this might not necessarily include all sites associated with a component.
    normalised_components = normalisestring.(components; isactivated=normalisecomponents)
    normalised_columnreference = normalisestring(columnreference)
    normalised_sitecolumnreference = normalisestring(sitecolumnreference)
    sites = Dict(components .=> [Set{String}() for _ ∈ 1:length(components)])
    for filepath ∈ filepaths
        csvtype = readcsvtype(filepath)
        csvtype != assocdata && continue

        df = CSV.File(filepath; header=3,lazystrings=true)
        csvheaders = String.(Tables.columnnames(df))
        normalised_csvheaders = normalisestring.(String.(Tables.columnnames(df)))

        lookupcolumnindex1 = findfirst(x -> x[1:end-1] == normalised_columnreference && x[end] == '1', normalised_csvheaders)
        lookupcolumnindex2 = findfirst(x -> x[1:end-1] == normalised_columnreference && x[end] == '2', normalised_csvheaders)
        isnothing(lookupcolumnindex1) && error("Header ", normalised_columnreference * "1", " not found.")
        isnothing(lookupcolumnindex2) && error("Header ", normalised_columnreference * "2", " not found.")
        lookupcolumn1 = Symbol(csvheaders[lookupcolumnindex1])
        lookupcolumn2 = Symbol(csvheaders[lookupcolumnindex2])

        normalised_sitecolumnreference = normalisestring(sitecolumnreference)
        lookupsitecolumnindex1 = findfirst(x -> x[1:end-1] == normalised_sitecolumnreference && x[end] == '1', normalised_csvheaders)
        lookupsitecolumnindex2 = findfirst(x -> x[1:end-1] == normalised_sitecolumnreference && x[end] == '2', normalised_csvheaders)
        isnothing(lookupcolumnindex1) && error("Header ", normalised_sitecolumnreference * "1", " not found.")
        isnothing(lookupcolumnindex2) && error("Header ", normalised_sitecolumnreference * "2", " not found.")
        lookupsitecolumn1 = Symbol(csvheaders[lookupsitecolumnindex1])
        lookupsitecolumn2 = Symbol(csvheaders[lookupsitecolumnindex2])

        for row ∈ Tables.rows(df)
            component1 = row[Symbol(lookupcolumn1)]
            component2 = row[Symbol(lookupcolumn2)]
            foundcomponentidx1 = findfirst(isequal(normalisestring(component1; isactivated=normalisecomponents)), normalised_components)
            foundcomponentidx2 = findfirst(isequal(normalisestring(component2; isactivated=normalisecomponents)), normalised_components)
            (isnothing(foundcomponentidx1) || isnothing(foundcomponentidx2)) && continue
            push!(sites[components[foundcomponentidx1]], row[Symbol(lookupsitecolumn1)])
            push!(sites[components[foundcomponentidx2]], row[Symbol(lookupsitecolumn2)])
        end
    end
    output = Array{Array{String,1}}(undef, 0)
    for component ∈ components
        push!(output, collect(sites[component]))
    end
    verbose && println("Found sites for ", components, " are ", output, ".")
    return output
end

function createemptyparamsarray(datatype::Type, csvtype::CSVType, components::Array{String,1})
    allcomponentsites = [String[] for _ in 1:length(components)]
    return createemptyparamsarray(datatype, csvtype, components, allcomponentsites)
end

function createemptyparamsarray(datatype::Type, csvtype::CSVType, components::Array{String,1}, allcomponentsites::Array{Array{String,1},1})
    # Creates a missing array of the appropriate size.
    componentslength = length(components)
    csvtype == singledata && return (Array{Union{Missing,datatype}}(undef, componentslength) .= missing)
    csvtype == pairdata && return (Array{Union{Missing,datatype}}(undef, componentslength, componentslength) .= missing)
    if csvtype == assocdata
        output = Array{Array{Union{Missing,datatype},2},2}(undef, componentslength, componentslength)
        for i ∈ 1:componentslength, j ∈ 1:componentslength
            output[i,j] = (Array{Union{Missing,datatype}}(undef, length(allcomponentsites[i]), length(allcomponentsites[j])) .= missing)
        end
        return output
    end
end

function createemptyparamsarray(csvtype::CSVType, components::Array{String,1})
    allcomponentsites = [String[] for _ in 1:length(components)]
    return createemptyparamsarray(Any, csvtype, components, allcomponentsites)
end

function createemptyparamsarray(csvtype::CSVType, components::Array{String,1}, allcomponentsites::Array{Array{String,1},1})
    return createemptyparamsarray(Any, csvtype, components, allcomponentsites)
end

function convertsingletopair(params::Array{T,1}; outputmissing::Bool=false) where T
    # Returns a missing square matrix with its diagonal matrix replaced by the given parameters.
    paramslength = length(params)
    if outputmissing
        output = Array{Union{Missing,T}}(undef, paramslength, paramslength) .= missing
    else
        if T == Any
            output = Array{T}(undef, paramslength, paramslength) .= 0
        elseif T <: Number
            output = Array{T}(undef, paramslength, paramslength) .= 0
        elseif T <: AbstractString
            output = Array{T}(undef, paramslength, paramslength) .= ""
        else
            error("Data type ", T, " not supported.")
        end
    end
    for i = 1:paramslength
        output[i,i] = params[i]
    end
    return output
end

function mirrormatrix!(matrix::Array{T,2}) where T
    # Mirrors a square matrix.
    matrixsize = size(matrix)
    matrixsize[1] != matrixsize[2] && error("Matrix is not square.")
    for i ∈ 2:matrixsize[1], j ∈ 1:i-1
        lowervalue = matrix[i,j]
        uppervalue = matrix[j,i]
        if !ismissing(lowervalue) && !ismissing(uppervalue) && lowervalue != uppervalue
            error("Dissimilar non-zero entries exist across diagonal:", lowervalue, ", ", uppervalue)
        end
        !ismissing(lowervalue) && (matrix[j,i] = lowervalue)
        !ismissing(uppervalue) && (matrix[i,j] = uppervalue)
    end
    return matrix
end

function mirrormatrix!(matrix::Array{Array{T,2},2}) where T
    # Mirrors a square matrix.
    matrixsize = size(matrix)
    matrixsize[1] != matrixsize[2] && error("Matrix is not square.")
    [mirrormatrix!(matrix[i,i]) for i ∈ 1:matrixsize[1]]
    for i ∈ 2:matrixsize[1], j ∈ 1:i-1
        matrixsize2 = size(matrix[i,j])
        for a ∈ 1:matrixsize2[1], b ∈ 1:matrixsize2[2]
            lowervalue = matrix[i,j][a,b]
            uppervalue = matrix[j,i][b,a]
            if !ismissing(lowervalue) && !ismissing(uppervalue) && lowervalue != uppervalue
                error("Dissimilar non-zero entries exist across diagonal:", lowervalue, ", ", uppervalue)
            end
            !ismissing(lowervalue) && (matrix[j,i][b,a] = lowervalue)
            !ismissing(uppervalue) && (matrix[i,j][a,b] = uppervalue)
        end
    end
    return matrix
end

export GroupParam
function GroupParam(gccomponents::Array{<:Any,1}, grouplocations::Array{String,1}=String[]; usergrouplocations::Array{String,1}=String[], verbose::Bool=false)
    # The format for gccomponents is an arary of either the species name (if it
    # available in the Clapeyron database, or a tuple consisting of the species
    # name, followed by a list of group => multiplicity pairs.  For example:
    # gccomponents = ["ethane",
    #                ("hexane", ["CH3" => 2, "CH2" => 4]),
    #                ("octane", ["CH3" => 2, "CH2" => 6])]
    BuildSpeciesType = Union{Tuple{String, Array{Pair{String, Int64},1}}, String, Tuple{String}}
    any(.!(isa.(gccomponents, BuildSpeciesType))) && error("The format of the components is incorrect.")
    filepaths = string.(vcat([(getpaths.(grouplocations; relativetodatabase=true)...)...], [(getpaths.(usergrouplocations)...)...]))
    components = String[]
    allcomponentgroups = Array{Array{String,1},1}(undef, 0)
    allcomponentngroups = Array{Array{Int,1},1}(undef, 0)
    componentstolookup = String[]
    append!(componentstolookup, [x for x ∈ gccomponents[isa.(gccomponents, String)]])
    append!(componentstolookup, [first(x) for x ∈ gccomponents[isa.(gccomponents, Tuple{String})]])
    allfoundcomponentgroups = Dict{String,String}()
    groupsourcecsvs = String[]
    for filepath in filepaths
        csvtype = readcsvtype(filepath)
        if csvtype != groupdata
            verbose && println("Skipping ", csvtype, " csv at ", filepath)
            continue
        end
        verbose && println("Searching for groups for components ", componentstolookup, " at ", filepath, "...")
        merge!(allfoundcomponentgroups, findgroupsincsv(componentstolookup, filepath; verbose=verbose))
        append!(groupsourcecsvs, [filepath])
    end
    for gccomponent ∈ gccomponents
        if gccomponent isa Tuple{String, Array{Pair{String, Int64},1}}
            append!(components, [gccomponent[1]])
            groupsandngroups = gccomponent[2]
        elseif gccomponent isa String
            !haskey(allfoundcomponentgroups, gccomponent) && error("Predefined component ", gccomponent, " not found in any group input csvs.")
            append!(components, [gccomponent])
            groupsandngroups = eval(Meta.parse(allfoundcomponentgroups[gccomponent]))
        elseif gccomponent isa Union{String, Tuple{String}}
            !haskey(allfoundcomponentgroups, first(gccomponent)) && error("Predefined component ", gccomponent, " not found in any group input csvs.")
            append!(components, [first(gccomponent)])
            groupsandngroups = eval(Meta.parse(allfoundcomponentgroups[first(gccomponent)]))
        end
        componentgroups = String[]
        componentngroups = Int[]
        for (componentgroup, ncomponentgroup) ∈ groupsandngroups
            append!(componentgroups, [componentgroup])
            append!(componentngroups, [ncomponentgroup])
        end
        append!(allcomponentgroups, [componentgroups])
        append!(allcomponentngroups, [componentngroups])
    end
    flattenedgroups = unique([(allcomponentgroups...)...])
    allcomponentnflattenedgroups = [zeros(Int, length(flattenedgroups)) for _ in 1:length(components)]
    for i in 1:length(components)
        for (k, group) in enumerate(flattenedgroups)
            if group in allcomponentgroups[i]
                allcomponentnflattenedgroups[i][k] = allcomponentngroups[i][findfirst(isequal(group), allcomponentgroups[i])]
            else
                allcomponentnflattenedgroups[i][k] = 0
            end
        end
    end
    return GroupParam(components, allcomponentgroups, allcomponentngroups, flattenedgroups, allcomponentnflattenedgroups, groupsourcecsvs)
end


text:
#=
    The file contains the implementation of the FGLM algorithm,
    an algorithm for ordering conversion
=#

mutable struct NextMonomials
    # monomials to check
    monoms::Vector{Int}
    load::Int
end

Base.isempty(m::NextMonomials) = m.load == 0

function initialize_nextmonomials(ht)
    etmp = ht.exponents[1]
    etmp .= UInt16(0)
    vidx = insert_in_hash_table!(ht, etmp)
    monoms = Vector{Int}(undef, 2^64)
    monoms[1] = vidx
    load = 1
    NextMonomials(monoms, load)
end

function insertnexts!(
            m::NextMonomials, basis::Basis,
            ht::MonomialHashtable, monom::Int)
    etmp = ht.exponents[1]
    etmp .= UInt16(0)
    emonom = ht.exponents[monom]
    eprod = Vector{UInt16}(undef, ht.explen)
    for i in ht.nvars
        etmp[i] = UInt16(1)
        eprod = divisi
    end
end

function extractnext!(m::NextMonomials)
    # assuming m.monoms is sorted (reversed)
    monom = m.monoms[m.load]
    m.load -= 1
    monom
end

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

function fglm_f4!(
            ring::PolyRing,
            basis::Basis,
            ht::MonomialHashtable)

    newbasis = initialize_basis(ring, basis.ntotal)
    nextmonoms = initialize_nextmonomials(ht)
    matrix = initialize_double_matrix()

    while !empty(nextmonoms)
        monom = extractnext!(nextmonoms)

        monom_nf, ht = normal_form_f4!(ring, basis, ht, monom)

        exists = linear_relation!(matrix, monom_nf)

        if exists
            add_generator!(newbasis, matrix)
        else
            insertnexts!(nextmonoms, newbasis, ht, monom)
        end
    end

    newbasis, ht
end

function fglm_f4(
            ring::PolyRing,
            basisexps::Vector{Vector{Vector{UInt16}}},
            basiscoeffs::Vector{Vector{UInt64}},
            rng::Rng,
            tablesize::Int=2^16)
    basis, ht = initialize_structures(ring, basisexps,
                                        basiscoeffs, rng, tablesize)

    basis, ht = fglm_f4!(ring, basis, ht, tobereduced)

    export_basis_data(basis, ht)
end

ext: jl


text:
__precompile__()

module InverseLaplace

import SpecialFunctions
# I have to export everything in order for Documenter.jl to find
# the strings. A few hours of work would probably be enough to solve the problem.
export ILt, setNterms
# export optimize, opteval # Broken at the moment
export Weeks, WeeksErr, setparameters
export Talbot, GWR, ILT
export TransformPair, ILtPair, abserr, iltpair_power
export ilt, talbot, gwr

abstract type AbstractILt end

"""
   ft::Talbot = Talbot(func::Function, Nterms::Integer=32)

Return `ft`, which estimates the inverse Laplace transform of `func` with
the fixed Talbot algorithm. `ft(t)` evaluates the transform at `t`.  You may
want to tune `Nterms` together with `setprecision(BigFloat, x)`.

# Example

Compute the inverse transform of the transform of `cos` at argument `pi/2`.
```julia-repl
julia> ft = Talbot(s -> s / (s^2 + 1), 80);

julia> ft(pi / 2)
-3.5366510684573195e-5
```

Note that given `Float64` input, the precision of the returned value may not be satisfactory.
```julia-repl
julia> Float64(ft(big(pi) / 2))
2.114425886215604e-49
```

!!! note
    This function uses the fixed Talbot method. It evaluates the Laplace transform
    function for complex arguments. The GWR method is, in general, less accurate and
    less stable, but does not evaluate the Laplace transform function for complex arguments.
"""
struct Talbot{T<:Base.Callable} <: AbstractILt
    Laplace_space_func::T  # Laplace space function
    Nterms::Int
end

const talbot_default_num_terms = 32
Talbot(Laplace_space_func::Base.Callable) = Talbot(Laplace_space_func, talbot_default_num_terms)

"""
   ft::GWR = GWR(func::Function, Nterms::Integer=16)

Return `ft`, which estimates the inverse Laplace transform of `func` with
the GWR algorithm. `ft(t)` evaluates the transform at `t`.

# Example

Compute the inverse transform of the transform of `cos` at argument `pi / 2`.
```
julia> ft = GWR(s -> s / (s^2 + 1), 16);

julia> ft(pi / 2)
-0.001
```
"""
struct GWR{T<:Base.Callable} <: AbstractILt
    Laplace_space_func::T  # Laplace space function
    Nterms::Int
end

const gwr_default_num_terms = 16
GWR(Laplace_space_func::Base.Callable) = GWR(Laplace_space_func, gwr_default_num_terms)

"""
    ILT(function, Nterms=32)

This is an alias for the default `Talbot()` method.
"""
ILT(args...) = Talbot(args...)

function Base.show(io::IO, ailt::AbstractILt)
    print(io, string(typeof(ailt)), "(Nterms=", ailt.Nterms, ')')
end

# TODO: get rid of this in favor of above.
# Rely on broadcasting, as well.
struct ILt{T<:Base.Callable, V<:Base.Callable} <: AbstractILt
    Laplace_space_func::T
    iltfunc::V
    Nterms::Int
end

"""
    itrans = ILt(func, iltfunc=talbot, Nterms=32)

 * deprecated. Use, ILT, Talbot, GWR, or Weeks instead *

Return an object that estimates the inverse Laplace transform of
the function `func` using the algorithm implemented by function `iltfunc`.
`itrans(t)` estimates the inverse transform for argument `t`.  The
accuracy of the estimates depends strongly on the choice of `iltfunc`, `t`, `Nterms`,
and the precision of the data type of the argument to `func`. The default value
of `32` may give extremely inaccurate estimates.

`iltfunc` may be either `talbot` or `gwr`.

## Example

```julia
julia> itr = ILt( s -> 1/(1 + s^2), talbot);

julia> itr([ pi/4, pi/2, 3*pi/4, -pi])
4-element Array{Float64,1}:
  0.707107
  1.0
  0.707107
 -3.66676e-12
```
"""
ILt(func, iltfunc) = ILt(func, iltfunc, talbot_default_num_terms)

# Make this the default
"""
    ilt(func::Function, t::AbstractFloat, M::Integer=32)

`ilt` is an alias for the default inverse Laplace transform method `talbot`.
"""
ilt(args...) = talbot(args...)
ILt(func) = ILt(func, talbot, talbot_default_num_terms)

"""
    setNterms(ailt::AbstractILt, Nterms::Integer)

set the number of terms used in the inverse Laplace tranform `ailt`. If
`ailt` stores internal data, it will be recomputed, so that subsequent
calls `ailt(t)` reflect the new value of `Nterms`.
"""
setNterms(ailt::AbstractILt, N::Integer) = (ailt.Nterms = N)

## Make the ILT data types callable.
(ailt::ILt)(t) = ailt.iltfunc(ailt.Laplace_space_func, t, ailt.Nterms)
(ailt::ILt)(t,N) = ailt.iltfunc(ailt.Laplace_space_func, t, N)

(ailt::Talbot)(t) = talbot(ailt.Laplace_space_func, t, ailt.Nterms)
(ailt::Talbot)(t, n::Integer) = talbot(ailt.Laplace_space_func, t, n)

(ailt::GWR)(t) = gwr(ailt.Laplace_space_func, t, ailt.Nterms)
(ailt::GWR)(t, n::Integer) = gwr(ailt.Laplace_space_func, t, n)

include("fixed_talbot.jl")
include("gwr.jl")
include("weeks.jl")
include("pairtest.jl")

end # module


text:
# Taken and modified from https://github.com/vtjnash/Polynomial.jl
# Polynomial.jl was deprecated in favor of Polynomials.jl, which uses reversed
# indexing. The former was copied to this file, and modified/trimmed to just
# the bare functions required for TransferFunction support.

Base.eps{T}(::Type{T}) = zero(T)
Base.eps{F<:AbstractFloat}(x::Type{F}) = Base.eps(F)
Base.eps{T}(x::Type{Complex{T}}) = eps(T)

immutable Poly{T<:Number}
    a::Vector{T}
    nzfirst::Int #for effiencicy, track the first non-zero index
    function Poly(a::Vector{T})
        la = length(a)
        i = 0
        for i = 1:la
            if abs(a[i]) > 2*eps(T)  break  end
        end
        new(a, i)
    end
end

Poly{T<:Number}(a::Vector{T}) = Poly{T}(a)

Base.convert{T}(::Type{Poly{T}}, p::Poly) = Poly(convert(Vector{T}, p.a))
Base.promote_rule{T, S}(::Type{Poly{T}}, ::Type{Poly{S}}) = Poly{promote_type(T, S)}
Base.eltype{T}(::Poly{T}) = T

Base.length(p::Poly) = length(p.a) - p.nzfirst + 1
Base.endof(p::Poly) = length(p)
deg(p::Poly) = length(p) - 1

Base.getindex(p::Poly, i) = p.a[i - 1 + p.nzfirst]
Base.getindex(p::Poly, c::Colon) = p.a[p.nzfirst:end]
Base.setindex!(p::Poly, v, i) = (p.a[i - 1 + p.nzfirst] = v)

Base.copy(p::Poly) = Poly(copy(p.a[p.nzfirst:end]))

Base.zero{T}(p::Poly{T}) = Poly([zero(T)])
Base.zero{T}(::Type{Poly{T}}) = Poly([zero(T)])
Base.one{T}(p::Poly{T}) = Poly([one(T)])
Base.one{T}(::Type{Poly{T}}) = Poly([one(T)])

function Base.show(io::IO, p::Poly)
    print(io,"Poly(")
    print(io,p)
    print(io,")")
end

Base.print(io::IO, p::Poly) = print_poly(io, p)

function print_poly{T}(io::IO, p::Poly{T}, var=:x)
    n = length(p)
    if n == 1
        print(io, p[1])
    else
        for j = 1:n
            pj = p[j]
            magpj = abs(pj)
            if magpj > 2*eps(T)
                if j == 1
                    real(pj) < 0 && print(io, "-")    #Prepend - if first and negative
                else
                    real(pj) < 0 ? print(io," - ") : print(io," + ")
                end
                #Print pj if pj is the last coefficient, or pj is not identically 1
                if j == n || abs(magpj - 1) > 2*eps(T)
                    print(io, magpj)
                end
                exp = n-j
                if exp > 0
                    print(io, var)
                    if exp > 1
                        print(io, '^', exp)
                    end
                end
            end
        end
    end
end

*(c::Number, p::Poly) = Poly(c * p.a[p.nzfirst:end])
*(p::Poly, c::Number) = Poly(c * p.a[p.nzfirst:end])
/(p::Poly, c::Number) = Poly(p.a[p.nzfirst:end] / c)
./(p::Poly, c::Number) = /(p, c)
-(p::Poly) = Poly(-p.a[p.nzfirst:end])

-(p::Poly, c::Number) = +(p, -c)
+(c::Number, p::Poly) = +(p, c)
function +(p::Poly, c::Number)
    if length(p) < 1
        return Poly([c,])
    else
        p2 = copy(p)
        p2.a[end] += c;
        return p2;
    end
end
function -(c::Number, p::Poly)
    if length(p) < 1
        return Poly([c,])
    else
        p2 = -p;
        p2.a[end] += c;
        return p2;
    end
end

function +{T,S}(p1::Poly{T}, p2::Poly{S})
    R = promote_type(T,S)
    n = length(p1)
    m = length(p2)
    if n > m
        a = Array(R, n)
        for i = 1:m
            a[n-m+i] = p1[n-m+i] + p2[i]
        end
        for i = 1:n-m
            a[i] = p1[i]
        end
    else
        a = Array(R, m)
        for i = 1:n
            a[m-n+i] = p1[i] + p2[m-n+i]
        end
        for i = 1:m-n
            a[i] = p2[i]
        end
    end
    Poly(a)
end

function -{T,S}(p1::Poly{T}, p2::Poly{S})
    R = promote_type(T,S)
    n = length(p1)
    m = length(p2)
    if n > m
        a = Array(R, n)
        for i = 1:m
            a[n-m+i] = p1[n-m+i] - p2[i]
        end
        for i = 1:n-m
            a[i] = p1[i]
        end
    else
        a = Array(R, m)
        for i = 1:n
            a[m-n+i] = p1[i] - p2[m-n+i]
        end
        for i = 1:m-n
            a[i] = -p2[i]
        end
    end
    Poly(a)
end

function *{T,S}(p1::Poly{T}, p2::Poly{S})
    R = promote_type(T,S)
    n = length(p1)
    m = length(p2)
    if n == 0 || m == 0
        return Poly(R[])
    end
    a = zeros(R, n+m-1)
    for i = 1:length(p1)
        for j = 1:length(p2)
            a[i+j-1] += p1[i] * p2[j]
        end
    end
    Poly(a)
end

function ==(p1::Poly, p2::Poly)
    if length(p1) != length(p2)
        return false
    else
        return p1.a[p1.nzfirst:end] == p2.a[p2.nzfirst:end]
    end
end

function isapprox(p1::Poly, p2::Poly; kwargs...)
    if length(p1) != length(p2)
        return false
    else
        return isapprox(p1.a[p1.nzfirst:end], p2.a[p2.nzfirst:end]; kwargs...)
    end
end

function polyval{T}(p::Poly{T}, x::Number)
    R = promote_type(T, typeof(x))
    lenp = length(p)
    if lenp == 0
        return zero(R)
    else
        y = convert(R, p[1])
        for i = 2:lenp
            y = p[i] + x.*y
        end
        return y
    end
end

# compute the roots of a polynomial
function roots{T}(p::Poly{T})
    R = promote_type(T, Float64)
    num_zeros = 0
    if length(p) == 0
        return zeros(R, 0)
    end
    while abs(p[end-num_zeros]) <= 2*eps(T)
        if num_zeros == length(p)-1
            return zeros(R, 0)
        end
        num_zeros += 1
    end
    n = length(p)-num_zeros-1
    if n < 1
        return zeros(R, length(p)-1)
    end
    companion = zeros(R, n, n)
    a0 = p[end-num_zeros]
    for i = 1:n-1
        companion[1,i] = -p[end-num_zeros-i] / a0
        companion[i+1,i] = 1;
    end
    companion[1,end] = -p[1] / a0
    D,V = eig(companion)
    r = zeros(eltype(D),length(p)-1)
    r[1:n] = 1./D
    return r
end


text:
# Utils for MathOptInterface
#
##################################################
# Import legacy from LinQuadOptInterface to ease the integration
# of KNITRO quadratic and linear facilities.
"""
    canonical_quadratic_reduction(func::ScalarQuadraticFunction)

Reduce a ScalarQuadraticFunction into three arrays, returned in the following
order:
 1. a vector of quadratic row indices
 2. a vector of quadratic column indices
 3. a vector of quadratic coefficients

Warning: we assume in this function that all variables are correctly
ordered, that is no deletion or swap has occured.
"""
function canonical_quadratic_reduction(func::MOI.ScalarQuadraticFunction)
    quad_columns_1, quad_columns_2, quad_coefficients = (
        Int32[term.variable_index_1.value for term in func.quadratic_terms],
        Int32[term.variable_index_2.value for term in func.quadratic_terms],
        [term.coefficient for term in func.quadratic_terms]
    )
    # Take care of difference between MOI standards and KNITRO ones.
    for i in 1:length(quad_coefficients)
        if quad_columns_1[i] == quad_columns_2[i]
            quad_coefficients[i] *= .5
        end
    end
    # Take care that Julia is 1-indexed.
    quad_columns_1 .-= 1
    quad_columns_2 .-= 1
    return quad_columns_1, quad_columns_2, quad_coefficients
end

"""
    canonical_linear_reduction(func::Quad)

Reduce a ScalarQuadraticFunction into two arrays, returned in the following
order:
 1. a vector of linear column indices
 2. a vector of linear coefficients

Warning: we assume in this function that all variables are correctly
ordered, that is no deletion or swap has occured.
"""
function canonical_linear_reduction(func::MOI.ScalarQuadraticFunction)
    affine_columns = Int32[term.variable_index.value for term in func.affine_terms]
    affine_coefficients = [term.coefficient for term in func.affine_terms]
    affine_columns .-= 1
    return affine_columns, affine_coefficients
end
function canonical_linear_reduction(func::MOI.ScalarAffineFunction)
    affine_columns = Int32[term.variable_index.value for term in func.terms]
    affine_coefficients = [term.coefficient for term in func.terms]
    affine_columns .-= 1
    return affine_columns, affine_coefficients
end

function canonical_vector_affine_reduction(func::MOI.VectorAffineFunction)
    index_cols = Int32[]
    index_vars = Int32[]
    coefs = Float64[]

    for t in func.terms
        push!(index_cols, t.output_index)
        push!(index_vars, t.scalar_term.variable_index.value)
        push!(coefs, t.scalar_term.coefficient)
    end
    index_cols .-= 1
    index_vars .-= 1
    return index_cols, index_vars, coefs
end

# Convert Julia'Inf to KNITRO's Inf.
function check_value(val::Float64)
    if val > KN_INFINITY
        return KN_INFINITY
    elseif val < -KN_INFINITY
        return -KN_INFINITY
    end
    return val
end


text:
# Core functions for hidden Markov models

# Copyright (c) 2016 Jeffrey W. Miller
# This software is released under the MIT License.
# If you use this software in your research, please cite:
#
#   Jeffrey W. Miller (2016). Lecture Notes on Advanced Stochastic Modeling. Duke University, Durham, NC.

# This code assumes the following have been defined:
#   log_q = function such that log_q(x_t,phi[z_t]) = log p(x_t | Z_t = z_t), where phi is a vector of parameters.
#   qrnd = function such that qrnd(phi[z_t]) generates a sample from p(x_t|z_t), where phi is a vector of parameters.
#   phi_init = function that creates an array of phi parameters.
#   phi_max! = function that finds the maximal phi parameters for the given data x and weights gamma.

# The functions below involve the following variables:
#   m = number of hidden states
#   n = length of observed sequence
#   log_pi = length m vector with log_pi[s] = log p(Z_1=s).
#   log_T = m-by-m matrix with log_T[r,s] = log p(Z_t=s | Z_{t-1}=r) = log of transition probability r->s.
#   phi = length m vector with phi[s] = parameters for the emission distribution p(x_t | z_t=s).

using Distributions

# Generate a sequence of n hidden states z and observations x.
# OUTPUT:
#    x,z = vectors x_{1:n} and z_{1:n} sampled from p(x_{1:n},z_{1:n}).
function generate(n,log_pi,log_T,phi)
    m = length(log_pi)
    z = zeros(Int,n)
    D = [Categorical(exp(vec(log_T[r,:]))) for r=1:m]
    z[1] = rand(Categorical(exp(log_pi)))
    for t = 2:n
        z[t] = rand(D[z[t-1]])
    end
    x = [qrnd(phi[z[t]]) for t=1:n]
    return x,z
end

# Find most probable sequence of hidden states z given observations x.
# OUTPUT: z = a vector z_{1:n} maximizing p(z_{1:n} | x_{1:n}).
function viterbi(x,log_pi,log_T,phi)
    m = length(log_pi)  # number of hidden states
    n = length(x)  # length of observed sequence
    F = zeros(n,m)  # maxima
    A = zeros(n,m)  # maximizers

    # Compute F and A
    for s = 1:m
        F[1,s] = log_pi[s] + log_q(x[1],phi[s])
    end
    for j = 2:n
        for s = 1:m
            maxval = -Inf
            maxind = 0
            for r = 1:m
                newval = F[j-1,r] + log_T[r,s]
                if newval > maxval
                    maxval = newval
                    maxind = r
                end
            end
            F[j,s] = maxval +  log_q(x[j],phi[s])
            A[j,s] = maxind
        end
    end

    # Recover optimal z
    z = zeros(Int,n)
    z[n] = indmax(F[n,:])
    for t = n-1:-1:1
        z[t] = A[t+1,z[t+1]]
    end
    return z
end

# Compute log(sum(exp(x))) in a way that avoids numerical underflow/overflow.
function logsumexp(x)
    mx = maximum(x)
    if mx==-Inf
        return -Inf
    elseif mx==Inf
        return Inf
    else
        s = 0.0
        for xi in x
            s += exp(xi-mx)
        end
        return log(s)+mx
    end
end

# Run the forward algorithm to compute log(p(x_{1:j},z_j) for each j, z_j.
# OUTPUT: G = n-by-m matrix with G[j,s] = log(p(x_{1:j},z_j=s)).
function forward(x,log_pi,log_T,phi)
    m = length(log_pi)  # number of hidden states
    n = length(x)  # length of observed sequence
    G = zeros(n,m)  # G[j,s] = log(p(x_{1:j},z_j=s))
    l = zeros(m)  # temporary vector

    for s = 1:m
        G[1,s] = log_pi[s] + log_q(x[1],phi[s])
    end
    for j = 2:n
        for s = 1:m
            for r = 1:m
                l[r] = G[j-1,r] + log_T[r,s] + log_q(x[j],phi[s])
            end
            G[j,s] = logsumexp(l)
        end
    end
    log_m = logsumexp(G[n,:])  # log(p(x_{1:n}))
    return G,log_m
end

# Run the backward algorithm to compute log(p(x_{j+1:n}|z_j)) for each j, z_j.
# OUTPUT: H = n-by-m matrix with H[j,s] = log(p(x_{j+1:n}|z_j=s)).
function backward(x,log_pi,log_T,phi)
    m = length(log_pi)  # number of hidden states
    n = length(x)  # length of observed sequence
    H = zeros(n,m)  # H[j,r] = log(p(x_{j+1:n}|z_j=r))
    l = zeros(m)  # temporary vector

    for r = 1:m
        H[n,r] = 0.0
    end
    for j = n-1:-1:1
        for r = 1:m
            for s = 1:m
                l[s] = H[j+1,s] + log_T[r,s] + log_q(x[j+1],phi[s])
            end
            H[j,r] = logsumexp(l)
        end
    end
    return H
end

# Estimate the HMM parameters using the Baum-Welch algorithm.
function estimate(x,m,tolerance)
    n = length(x)  # length of observed sequence

    # randomly initialize parameters
    log_pi = log(rand(Dirichlet(m,1)))
    log_T = zeros(m,m)
    for r = 1:m; log_T[r,:] = log(rand(Dirichlet(m,1))); end
    phi = phi_init(m,x)

    # initialize EM variables
    gamma = zeros(n,m)  # gamma[j,r] = p(z_j=r | x)
    delta = zeros(n-1,m,m)  # delta[j,r,s] = p(z_j=r, z_{j+1}=s | x)
    log_m = -Inf
    log_m_old = 0.0
    n_iterations = 0

    println("\nEstimating parameters using Baum-Welch...")
    while abs(log_m - log_m_old) > tolerance
        log_m_old = log_m

        # E-step: compute gamma and delta
        G,log_m = forward(x,log_pi,log_T,phi)
        H = backward(x,log_pi,log_T,phi)
        for j = 1:n
            for r = 1:m
                gamma[j,r] = exp(G[j,r] + H[j,r] - log_m)
            end
        end
        for j = 1:n-1
            for r = 1:m
                for s = 1:m
                    delta[j,r,s] = exp(G[j,r] + log_T[r,s] +
                        log_q(x[j+1],phi[s]) + H[j+1,s] - log_m)
                end
            end
        end

        # M-step: update log_pi, log_T, and phi
        log_pi = vec(log(gamma[1,:] ./ sum(gamma[1,:])))
        C = squeeze(sum(delta,1),1)
        log_T = log(C./sum(C,2))
        phi_max!(phi,x,gamma)

        # check change in log likelihood
        @assert(log_m > log_m_old)

        n_iterations += 1
        println(n_iterations, " ", log_m)
    end
    return log_pi,log_T,phi,log_m
end

max_stars_repo_name: jwmi/HMM

max_issues_repo_name: jwmi/HMM

max_forks_repo_name: jwmi/HMM


text:
# Script to obtain sectoral WTO MFN tariffs

# In the paper they use EU MFN tariffs as a proxy for hard Brexit (take EU as reporter country)
# Tariffs are levied at product level (HS6), thus they need to be aggregated to the 56 sectors of the WIOD table.
# Use the correspondence of HS2007 and CPA2008 (obtained via Eurostat: RAMON)
# WTO tariff data was obtained via: http://tariffdata.wto.org/Default.aspx?culture=en-US

using DataFrames, CSV, Statistics, Plots, XLSX

# ---------- Correspondence tables -------------------------------------------------------------------------------------------------------

# HS2007 vis-a-vis CPA2008
raw_HS_CPA = CSV.read("data/HS2007_CPA2008.csv", DataFrame)[2:end, 1:2]
rename!(raw_HS_CPA, [:HS2007, :CPA2008])
transform!(raw_HS_CPA, [:HS2007, :CPA2008] .=> ByRow(string) .=> [:HS2007, :CPA2008], renamecols=false) # change type to string
raw_HS_CPA = replace.(raw_HS_CPA, r"[.]" => "") # take out dot of string
f(x) = x[1:2] # select first two elements of cell
raw_HS_CPA[:, :CPA2digits] = f.(raw_HS_CPA[:, :CPA2008]) # aggregates CPA2008 to two digit level

# CPA2008 (2 digits) vis-a-vis NACE Rev.2 (and the WIOD classification of sectors which aggregates some NACE)
raw_CPA_NACE = CSV.read("data/CPA2008_NACE2.csv", DataFrame)
transform!(raw_CPA_NACE, names(raw_CPA_NACE) .=> ByRow(string) .=> names(raw_CPA_NACE), renamecols=false) # change type to string
raw_CPA_NACE[:, :CPA2digits] = lpad.(raw_CPA_NACE[:, :CPA2digits], 2, '0') # pad with leading zeros so all have 2 digits

# combine HS2007, CPA2008, CPA2digits, NACE Rev. 2
c_table_HS_NACE = leftjoin(raw_HS_CPA, raw_CPA_NACE, on=:CPA2digits)

# HS2007 vis-a-vis HS2017
raw_HS_HS = CSV.read("data/HS2007_HS2017.csv", DataFrame)
raw_HS_HS = lpad.(raw_HS_HS[:, [:HS2017, :HS2007]], 6, '0') # pad with leading zeros so all have 6 digits

# add correspondence vis-a-vis HS2017 as tariff data from WTO uses HS2017
c_table_HS_NACE = leftjoin(c_table_HS_NACE, raw_HS_HS, on=:HS2007)
describe(c_table_HS_NACE)[:, :nmissing] # 62 of 5642 HS2007 codes do not have a corresponding HS2017 code (would say thats pretty good!)

XLSX.writetable("clean/table_correspondence_HS_NACE.xlsx", c_table_HS_NACE, overwrite=true)

# ---------- WTO tariffs -----------------------------------------------------------------------------------------------------------------

begin # vector for columns to take out
delete = ["Base nomenclature", "Number of subheadings", "Number of lines", "Number of AVDuties", "Number of non AVDuties", "List of non AVDuties",
 "List of unit values", "Exchange rate", "Import value"]
end

raw_tariffs = CSV.read("data/MFN_tariffs_2020.txt", DataFrame)[:, Not(delete)]
col_names = ["reporter", "year", "HS_level", "HS2017", "avtariff_avg", "avtariff_min", "avtariff_max", "tariff_free", "HS2017_long"] # new names
rename!(raw_tariffs, col_names) # rename columns
transform!(raw_tariffs, [:reporter, :HS2017, :HS2017_long] .=> ByRow(string) .=> [:reporter, :HS2017, :HS2017_long], renamecols=false) # change type to string

# depending on digit code pad with leading zeros (not absolutely necessary since we lose all codes with digits < 6 when matched with correspondence tables)
for i in 1:size(raw_tariffs, 1)
    raw_tariffs.HS2017[i] = lpad(raw_tariffs.HS2017[i], raw_tariffs.HS_level[i], "0")
end

# ---------- Correspondence: Tariffs to NACE Rev.2 -----------------------------------------------------------------------------------------------------------------

dropmissing!(c_table_HS_NACE) # cannot join on column with missing values
raw_tariffs = raw_tariffs[map(x->length(x), raw_tariffs[:, :HS2017]).==6, :] # only keep 6 digit HS codes

tariff_NACE = leftjoin(raw_tariffs, c_table_HS_NACE, on=:HS2017)
describe(tariff_NACE) # 172 of 11.330 missing due to matching, 623 of 11.330 missing due to data reported missing
# what to do with missing tariff data? either drop or assume zero tariffs, i.e. impacts the calculation of the average

# ---------- Replicate Figure 1: Average tariffs imposed by EU -----------------------------------------------------------------------------------------------------

dropmissing!(tariff_NACE) # still need to decide what to do with missing tariff data - for now just drop all of them

# data for averages with zero tariffs, maximum and minimum (always 0)
gdf = groupby(tariff_NACE, [:reporter, :WIOD, :Figure1])
d_tariff_NACE = combine(gdf, :avtariff_avg => mean => :avg, :avtariff_avg => maximum => :max) # we lose some sectors (non-tradables, for example public administration)
filter!(row -> row.reporter == "European Union", d_tariff_NACE)
sort!(d_tariff_NACE, :avg, rev=true)
#filter!(row -> row.reporter == "United States of America", d_tariff_NACE)

# data for sample average without zero tariffs
avg_zeros = filter(x -> x.avtariff_avg > 0 && x.reporter == "European Union", tariff_NACE)
gdf = groupby(avg_zeros, [:reporter, :WIOD, :Figure1])
avg_zeros = combine(gdf, :avtariff_avg => mean => :avg)

# Figure 1
begin
    figure1 = bar(d_tariff_NACE.max, xticks=(1:size(d_tariff_NACE,1),
     d_tariff_NACE.Figure1), ylims=(0,25), xrotation=90, label="min/max", color=:orange)
    xlabel!("NACE Rev. 2 activties")
    ylabel!("AV Applied Tariff (%)")
    title!("Sectoral tariff rates applied by the EU")
    scatter!(d_tariff_NACE.avg, label="average", color=:blue)
    plot!(fill(mean(d_tariff_NACE.avg), size(d_tariff_NACE,1)), width=2,
     color=:red, label="sample average of tariffs = $(round(mean(d_tariff_NACE.avg), digits=1)) ")
    plot!(fill(mean(avg_zeros.avg), size(avg_zeros, 1)), width=2,
     color=:green, label="sample average of non-zero tariffs = $(round(mean(avg_zeros.avg), digits=1))")
    figure1
    savefig(figure1, "images/figure1.png")
end

# ---------- Final table containinng all tariffs (MFN, NTB), countries, sectors -------------------------------------------------------------------------------------------

# also include non tariff barriers (NTBs) for "soft" brexit (2.77%) and "hard" brexit (8.31%)
# NTBs are based on Dhingra et al. (2017) and Berden et al. (2009)

WIOD = raw_CPA_NACE.WIOD
unique!(WIOD)

begin
    EU27 = ["AUT", "BEL", "BGR", "HRV", "CYP", "CZE", "DNK", "EST", "FIN", "FRA", "DEU", "GRC", "HUN", "IRL", "ITA", "LVA",
    "LTU", "LUX", "MLT", "NLD", "POL", "PRT", "ROU", "SVK", "SVN", "ESP", "SWE"]
end

N = size(EU27, 1) + 1 # EU27 + GBR
S = size(WIOD, 1)

country = repeat([EU27; "GBR"], inner=S)
WIOD = repeat(WIOD, outer=N)
d_ctry_WIOD = DataFrame(iso3 = country, WIOD = WIOD, NTB_soft=2.77, NTB_hard=8.31)

gdf = groupby(tariff_NACE, [:reporter, :WIOD])
avg_tariff_NACE = combine(gdf, :avtariff_avg => mean => :avg_tariff) # we lose some sectors (non-tradables, for example public administration)
filter!(row -> row.reporter == "European Union", avg_tariff_NACE)

all_tariffs = leftjoin(d_ctry_WIOD, avg_tariff_NACE[:, Not(:reporter)], on=:WIOD) # we do not have tariffs in 27 sectors (get lost when grouping data according to sectors)
# what to do? assume zero tariffs since not reported
all_tariffs = coalesce.(all_tariffs, 0.0) # assume zero tariffs for missing sectors
sort!(all_tariffs, [:iso3, :WIOD])

XLSX.writetable("clean/all_tariffs.xlsx", all_tariffs, overwrite=true)


text:
struct FailureBasedVariableSelection{TakeObjective} <: AbstractVariableSelection{TakeObjective} end

FailureBasedVariableSelection(;take_objective=true) = FailureBasedVariableSelection{take_objective}()

function (::FailureBasedVariableSelection{false})(cpmodel::CPModel)
    selectedVar = nothing
    minSize = Inf
    for (k, x) in branchable_variables(cpmodel)
        if x !== cpmodel.objective && !isbound(x) && length(x.domain) / cpmodel.statistics.infeasibleStatusPerVariable[parentId(x)] <= minSize
            selectedVar = x
            minSize = length(x.domain) / cpmodel.statistics.infeasibleStatusPerVariable[parentId(x)]
        end
    end
    if isnothing(selectedVar) && !isbound(cpmodel.objective)
        return cpmodel.objective
    end
    return selectedVar
end

function (::FailureBasedVariableSelection{true})(cpmodel::CPModel)
    selectedVar = nothing
    minSize = Inf
    for (k, x) in branchable_variables(cpmodel)
        if !isbound(x) && length(x.domain)/ cpmodel.statistics.infeasibleStatusPerVariable[parentId(x)] <= minSize
            selectedVar = x
            minSize = length(x.domain) / cpmodel.statistics.infeasibleStatusPerVariable[parentId(x)]
        end
    end
    return selectedVar
end


text:
include("1.2.1 Quadratic/.jl")
include("1.2.2 Quartic/.jl")
include("1.2.3 General/.jl")
include("1.2.4 Improper/.jl")

max_stars_repo_path: src/rules/1 Algebraic functions/1.2 Trinomial products/.jl


text:
# Henyey-Greenstein phase function.
phg(g, μ) = (1 / 4π) * (1 - g^2) / (1 + g^2 - 2 * g * μ)^(3/2)

# Sampling from the HG phase function.
# Here ξ in a random variable uniformly distributed between 0 and 1
μhg(g, ξ) = (1 + g^2 - ((1 - g^2) / (1 + g * (2 * ξ - 1)))^2) / (2 * g)

# Rayleigh phase function
pr(μ) = (3 / 16π) * (1 + μ^2)

# Hack to allow interchanging both phase functions
pr(g, μ) = pr(μ)

# Isotropic scattering
piso(μ) = (1 / 4π)
piso(g, μ) = piso(μ)

# Sample from the Rayleigh phase function.  The PDF notes are wrong.
function μr(ξ)
    z = 2 * (2ξ - 1)
    w = z + sqrt(z^2 + 1)
    w^(1//3) - w^(-1//3)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "AcousticsToolbox"
version = v"2020.11.4"

# Collection of sources required to complete build
sources = [
    ArchiveSource("http://oalib.hlsresearch.com/AcousticsToolbox/at_2020_11_4.zip", "6f9fe5d3d7fd99fc34fb2a2c3308f2e660a1cc0ea93f0651111af698acf442f0")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/at_*
rm -rf ../__MACOSX
perl -p -i -e 's/\-march=native//; s/\-ffast\-math//; s/\-mtune=native//;' Makefile
find . -name *.exe -exec rm {} \;
make
mkdir -p $bindir
find . -name *.exe -exec cp {} $bindir \;
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
platforms = expand_gfortran_versions(platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("field3d.exe", :field3d),
    ExecutableProduct("field.exe", :field),
    ExecutableProduct("kraken.exe", :kraken),
    ExecutableProduct("krakenc.exe", :krakenc),
    ExecutableProduct("bellhop3d.exe", :bellhop3d),
    ExecutableProduct("sparc.exe", :sparc),
    ExecutableProduct("scooter.exe", :scooter),
    ExecutableProduct("bounce.exe", :bounce),
    ExecutableProduct("bellhop.exe", :bellhop)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
using Test
using Distributed
using DistributedArrays

# Disable scalar indexing to avoid falling back on generic methods
# for AbstractArray
DistributedArrays.allowscalar(false)

# add at least 3 worker processes
if nworkers() < 3
    n = max(3, min(8, Sys.CPU_THREADS))
    addprocs(n; exeflags=`--check-bounds=yes`)
end
@assert nprocs() > 3
@assert nworkers() >= 3

@everywhere using Distributed
@everywhere using DistributedArrays
@everywhere using DistributedArrays.SPMD
@everywhere using Random
@everywhere using LinearAlgebra

@everywhere Random.seed!(1234 + myid())

const MYID = myid()
const OTHERIDS = filter(id-> id != MYID, procs())[rand(1:(nprocs()-1))]

function check_leaks()
    if length(DistributedArrays.refs) > 0
        sleep(0.1)  # allow time for any cleanup to complete and test again
        length(DistributedArrays.refs) > 0 && @warn("Probable leak of ", length(DistributedArrays.refs), " darrays")
    end
end

include("darray.jl")
include("spmd.jl")


text:
atol = 1e-5
rtol = 0.0
maxit = 100

println("=====================")
println("  Performance Tests ")
println("=====================")
for (f, df, init, _, name) in tests
   println("-------------------")
   println("  $name Test")
   println("-------------------")
   x, it_hist, ierr, x_hist = nsoli(init(), f,
                                 atol = atol, rtol = rtol, maxit = maxit)
   # improve the solution a bit
   y = copy(x)
   for i = 1:3
      y -= df(y) \ f(y)
   end
   println("     Error Message: ", ierr)
   println("      # iterations: ", length(x_hist) )
   println("          Residual: ", norm(f(x), Inf))
   println("     Residual xref: ", norm(f(y), Inf))
   println("  Error |x - xref|: ", norm(x - y, Inf))
end


text:
"""
    Base.rand(markovchain::MarkovChain, ntimesteps::Integer) -> Vector{Int}

Simulate `markovchain` for `ntimesteps` time-steps, returning one sample path (i.e.
trajectory) of length `ntimesteps`.

The time steps represented in the sample path are `0, 1, ..., ntimesteps-1`, for we start
with an initial state.

##### Complexity
- Time complexity: `O(n)`.
"""
function Base.rand(markovchain::MarkovChain, ntimesteps::Integer)
    rowmarginals = Vector{Categorical}(undef, length(markovchain.initialdist))
    for (index, row) in enumerate(eachrow(markovchain.transitionmatrix))
        rowmarginals[index] = Categorical(row[:,1])
    end

    # allocate memory
    realizationvec = Vector{Int}(undef, ntimesteps)
    realizationvec[1] = rand(Categorical(markovchain.initialdist))
    for step in 2:ntimesteps
        # sample from the distribution associated with the previous state's row
        realizationvec[step] = rand(rowmarginals[realizationvec[step-1]])
    end
    return realizationvec
end

max_stars_repo_name: LoganDGraham/StochasticProcesses

max_issues_repo_name: LoganDGraham/StochasticProcesses

max_forks_repo_name: LoganDGraham/StochasticProcesses


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _mag2flux{T<:AbstractFloat}(mag::T, zero_point::T, ABwave::T)
    if isnan(ABwave)
        return exp10(-0.4*(mag + zero_point))
    else
        return exp10(-0.4*(mag + 2.406 + 5*log10(float(ABwave))))
    end
end

"""
    mag2flux(mag[, zero_point, ABwave=number]) -> flux

### Purpose ###

Convert from magnitudes to flux expressed in erg/(s cm² Å).

### Explanation ###

This is the reverse of `flux2mag`.

### Arguments ###

* `mag`: the magnitude to be converted in flux.  It can be either a scalar or an
  array.
* `zero_point`: scalar giving the zero point level of the magnitude.  If not
 supplied then defaults to 21.1 (Code et al 1976).  Ignored if the `ABwave`
 keyword is supplied
* `ABwave` (optional numeric keyword): wavelength, scalar or array, in
 Angstroms.  If supplied, then the input `mag` is assumed to contain Oke AB
 magnitudes (Oke & Gunn 1983, ApJ, 266, 713;
 http://adsabs.harvard.edu/abs/1983ApJ...266..713O).

### Output ###

The flux.  It is of the same type, scalar or array, as `mag`.

If the `ABwave` keyword is set, then the flux is given by the expression

\$\$\\text{flux} = 10^{-0.4(\\text{mag} +2.406 + 4\\log_{10}(\\text{ABwave}))}\$\$

Otherwise the flux is given by

\$\$\\text{flux} =  10^{-0.4(\\text{mag} + \\text{zero point})}\$\$

### Example ###

``` julia
mag2flux(8.3)
# => 1.7378008287493692e-12
mag2flux(8.3, 12)
# => 7.58577575029182e-9
mag2flux(8.3, ABwave=12)
# => 3.6244115683017193e-7
```

### Notes ###

Code of this function is based on IDL Astronomy User's Library.
"""
mag2flux(mag::Real, zero_point::Real=21.1; ABwave::Real=NaN) =
    _mag2flux(promote(float(mag), float(zero_point), float(ABwave))...,)

function mag2flux{N<:Real}(mag::AbstractArray{N}, zero_point::Real=21.1;
                           ABwave::Real=NaN)
    flux = similar(mag, typeof(float(one(N))))
    for i in eachindex(mag)
        flux[i] = mag2flux(mag[i], zero_point, ABwave=ABwave)
    end
    return flux
end


text:
raw"""

@article{marching_cubes_jgt,
    author = {Thomas Lewiner and Hélio Lopes and Antônio Wilson Vieira and Geovan Tavares},
    title = {Efficient implementation of marching cubes cases with topological guarantees},
    year = {2003},
    month = {december},
    journal = {Journal of Graphics Tools},
    volume = {8},
    number = {2},
    pages = {1--15},
    publisher = {A.K.Peters},
    doi = {10.1080/10867651.2003.10487582},
    url = {\url{http://thomas.lewiner.org/pdfs/marching_cubes_jgt.pdf}}
}

Adapted to `Julia` by T Bltg (github.com/t-bltg).
"""

module MarchingCubes

import Base: RefValue
using StaticArrays

export MC, march, march_legacy

include("lut.jl")

const Vertex{F} = SVector{3,F}
const Normal{F} = SVector{3,F}
const Triangle = SVector{3,Int}

"""
    MC(vol::Array{F,3}, I::Type = Int; x = [], y = [], z = [])

# Description
Structure to hold temporaries and output of the marching cubes algorithm.
Vertices, normals and triangles are accessible using `m.vertices`, `m.normals` and `m.triangles`.
1-based indexing is assumed.
Optionally, pass `x`, `y`, `z` to normalize the vertices coordinates instead of `0:nx-1`, `0:ny-1` and `0:nz-1`.

# Arguments
    - `vol`: rank 3 array of floats (volume) on which the Marching Cubes algorithm is applied.
    - `I`: Int32, Int64, ...: vertices / normals / triangles index `Integer` type (defaults to `Int`).
    - `x` normalize vertex.x to the coordinates of `x`.
    - `y` normalize vertex.y to the coordinates of `y`.
    - `z` normalize vertex.z to the coordinates of `z`.
"""
struct MC{F,I}
    nx::Int  # height of the grid
    ny::Int  # width of the grid
    nz::Int  # depth of the grid
    vol::RefValue{Array{F,3}}  # implicit function values sampled on the grid
    cube::MVector{8,F}  # values of the implicit function on the active cube
    tv::MVector{3,I}  # triangle buffer
    vtx::MVector{3,F}  # vertex buffer
    nrm::MVector{3,F}  # normal buffer
    vert_indices::Array{I,4}  # pre-computed vertex indices
    triangles::Vector{Triangle}  # output triangles
    vertices::Vector{Vertex{F}}  # output vertex positions
    normals::Vector{Normal{F}}  # output vertex normals
    x::RefValue{Vector{F}}
    y::RefValue{Vector{F}}
    z::RefValue{Vector{F}}

    MC(
        vol::Array{F,3},
        I::Type{G} = Int;
        x::AbstractVector{F} = F[],
        y::AbstractVector{F} = F[],
        z::AbstractVector{F} = F[],
    ) where {F<:AbstractFloat,G<:Integer} = begin
        m = new{F,I}(
            size(vol)...,
            Ref(vol),
            zeros(F, 8),
            zeros(Int, 3),
            zeros(F, 3),
            zeros(F, 3),
            zeros(I, (3, size(vol)...)),
            Triangle[],
            Vertex[],
            Normal[],
            Ref(x),
            Ref(y),
            Ref(z),
        )
        sz = size(vol) |> prod
        sizehint!(m.triangles, nextpow(2, sz ÷ 6))
        sizehint!(m.vertices, nextpow(2, sz ÷ 2))
        sizehint!(m.normals, nextpow(2, sz ÷ 2))
        m
    end
end

"""
    lut_entry(vol, cb, i, j, k, iso)

# Description
Cube sign representation.
"""
lut_entry(vol::Array{F,3}, cb, i, j, k, iso) where {F} = begin
    lut = UInt8(0)
    @inbounds for p ∈ 0:7
        c = vol[i+((p⊻(p>>1))&1), j+((p>>1)&1), k+((p>>2)&1)] - iso
        abs(c) < eps(F) && (c = eps(F))
        c > 0 && (lut += (UInt8(1) << p))
        cb[p+1] = c
    end
    Int(lut) + 1
end

"""
    denormalize(m::MC)

# Description
Optional denormalization step of vertices to user coordinates.
"""
denormalize(m::MC) = begin
    if length(m.x[]) > 0 && length(m.y[]) > 0 && length(m.z[]) > 0
        mx, Mx = extrema(m.x[])
        my, My = extrema(m.y[])
        mz, Mz = extrema(m.z[])
        scl =
            @SVector([Mx - mx, My - my, Mz - mz]) ./
            @SVector([m.nx - 1, m.ny - 1, m.nz - 1])
        off = @SVector([mx, my, mz])
        @inbounds for (n, v) in enumerate(m.vertices)
            m.vertices[n] = Vertex(off .+ v .* scl)
        end
    end
    return
end

"""
    march_legacy(m::MC, isovalue::Number)

# Description
Original Marching Cubes algorithm.

# Arguments
    - `m`: Marching Cubes data structure.
    - `isovalue`: isosurface value.
"""
march_legacy(m::MC{F}, isovalue::Number = 0) where {F} = begin
    empty!(m.triangles)
    empty!(m.vertices)
    empty!(m.normals)

    vol = m.vol[]
    iso = F(isovalue)

    compute_intersection_points(m, vol, m.cube, iso)

    @inbounds for k ∈ 1:m.nz-1, j ∈ 1:m.ny-1, i ∈ 1:m.nx-1
        lut = lut_entry(vol, m.cube, i, j, k, iso)
        nt = 0
        while casesClassic[lut][3nt+1] > 0
            nt += 1
        end
        add_triangle(m, i, j, k, casesClassic[lut], nt)
    end

    denormalize(m)
    return
end

"""
    march(m::MC, isovalue::Number)

# Description
Topologically controlled Marching Cubes algorithm.

Arguments
    - `m`: Marching Cubes data structure.
    - `isovalue`: isosurface value.
"""
march(m::MC{F}, isovalue::Number = 0) where {F} = begin
    empty!(m.triangles)
    empty!(m.vertices)
    empty!(m.normals)

    vol = m.vol[]
    iso = F(isovalue)

    compute_intersection_points(m, vol, m.cube, iso)

    @inbounds for k ∈ 1:m.nz-1, j ∈ 1:m.ny-1, i ∈ 1:m.nx-1
        lut = lut_entry(vol, m.cube, i, j, k, iso)

        case = cases[lut][1]  # case of the active cube in [1; 16]
        cfg = cases[lut][2]  # configuration of the active cube
        subcfg = 1  # subconfiguration of the active cube

        if case == 1
            nothing
        elseif case == 2
            add_triangle(m, i, j, k, tiling1[cfg], 1)
        elseif case == 3
            add_triangle(m, i, j, k, tiling2[cfg], 2)
        elseif case == 4
            if test_face(m.cube, test3[cfg])
                add_triangle(m, i, j, k, tiling3_2[cfg], 4)  # 3.2
            else
                add_triangle(m, i, j, k, tiling3_1[cfg], 2)  # 3.1
            end
        elseif case == 5
            if test_interior(case, m.cube, cfg, subcfg, test4[cfg])
                add_triangle(m, i, j, k, tiling4_1[cfg], 2)  # 4.1.1
            else
                add_triangle(m, i, j, k, tiling4_2[cfg], 6)  # 4.1.2
            end
        elseif case == 6
            add_triangle(m, i, j, k, tiling5[cfg], 3)
        elseif case == 7
            if test_face(m.cube, test6[cfg][1])
                add_triangle(m, i, j, k, tiling6_2[cfg], 5)  # 6.2
            else
                if test_interior(case, m.cube, cfg, subcfg, test6[cfg][2])
                    add_triangle(m, i, j, k, tiling6_1_1[cfg], 3)  # 6.1.1
                else
                    add_triangle(m, i, j, k, tiling6_1_2[cfg], 9, add_c_vertex(m, i, j, k))  # 6.1.2
                end
            end
        elseif case == 8
            test_face(m.cube, test7[cfg][1]) && (subcfg += 1)
            test_face(m.cube, test7[cfg][2]) && (subcfg += 2)
            test_face(m.cube, test7[cfg][3]) && (subcfg += 4)
            if subcfg == 1
                add_triangle(m, i, j, k, tiling7_1[cfg], 3)
            elseif subcfg == 2 || subcfg == 3
                add_triangle(m, i, j, k, tiling7_2[cfg][subcfg-1], 5)
            elseif subcfg == 4
                add_triangle(m, i, j, k, tiling7_3[cfg][1], 9, add_c_vertex(m, i, j, k))
            elseif subcfg == 5
                add_triangle(m, i, j, k, tiling7_2[cfg][3], 5)
            elseif subcfg == 6 || subcfg == 7
                add_triangle(m, i, j, k, tiling7_3[cfg][subcfg-4], 9, add_c_vertex(m, i, j, k))
            elseif subcfg == 8
                if test_interior(case, m.cube, cfg, subcfg, test7[cfg][4])
                    add_triangle(m, i, j, k, tiling7_4_2[cfg], 9)
                else
                    add_triangle(m, i, j, k, tiling7_4_1[cfg], 5)
                end
            end
        elseif case == 9
            add_triangle(m, i, j, k, tiling8[cfg], 2)
        elseif case == 10
            add_triangle(m, i, j, k, tiling9[cfg], 4)
        elseif case == 11
            if test_face(m.cube, test10[cfg][1])
                if test_face(m.cube, test10[cfg][2])
                    add_triangle(m, i, j, k, tiling10_1_1_[cfg], 4)  # 10.1.1
                else
                    add_triangle(m, i, j, k, tiling10_2[cfg], 8, add_c_vertex(m, i, j, k))  # 10.2
                end
            else
                if test_face(m.cube, test10[cfg][2])
                    add_triangle(m, i, j, k, tiling10_2_[cfg], 8, add_c_vertex(m, i, j, k))  # 10.2
                else
                    if test_interior(case, m.cube, cfg, subcfg, test10[cfg][3])
                        add_triangle(m, i, j, k, tiling10_1_1[cfg], 4)  # 10.1.1
                    else
                        add_triangle(m, i, j, k, tiling10_1_2[cfg], 8)  # 10.1.2
                    end
                end
            end
        elseif case == 12
            add_triangle(m, i, j, k, tiling11[cfg], 4)
        elseif case == 13
            if test_face(m.cube, test12[cfg][1])
                if test_face(m.cube, test12[cfg][2])
                    add_triangle(m, i, j, k, tiling12_1_1_[cfg], 4)  # 12.1.1
                else
                    add_triangle(m, i, j, k, tiling12_2[cfg], 8, add_c_vertex(m, i, j, k))  # 12.2
                end
            else
                if test_face(m.cube, test12[cfg][2])
                    add_triangle(m, i, j, k, tiling12_2_[cfg], 8, add_c_vertex(m, i, j, k))  # 12.2
                else
                    if test_interior(case, m.cube, cfg, subcfg, test12[cfg][3])
                        add_triangle(m, i, j, k, tiling12_1_1[cfg], 4)  # 12.1.1
                    else
                        add_triangle(m, i, j, k, tiling12_1_2[cfg], 8)  # 12.1.2
                    end
                end
            end
        elseif case == 14
            test_face(m.cube, test13[cfg][1]) && (subcfg += 1)
            test_face(m.cube, test13[cfg][2]) && (subcfg += 2)
            test_face(m.cube, test13[cfg][3]) && (subcfg += 4)
            test_face(m.cube, test13[cfg][4]) && (subcfg += 8)
            test_face(m.cube, test13[cfg][5]) && (subcfg += 16)
            test_face(m.cube, test13[cfg][6]) && (subcfg += 32)
            if (sc = subcfg13[subcfg]) == 0  # 13.1
                add_triangle(m, i, j, k, tiling13_1[cfg], 4)
            elseif sc ≥ 1 && sc ≤ 6  # 13.2
                add_triangle(m, i, j, k, tiling13_2[cfg][sc], 6)
            elseif sc ≥ 7 && sc ≤ 18  # 13.3
                add_triangle(m, i, j, k, tiling13_3[cfg][sc-6], 10, add_c_vertex(m, i, j, k))
            elseif sc ≥ 19 && sc ≤ 22  # 13.4
                add_triangle(m, i, j, k, tiling13_4[cfg][sc-18], 12, add_c_vertex(m, i, j, k))
            elseif sc ≥ 23 && sc ≤ 26  # 13.5
                subcfg = sc - 22
                if test_interior(case, m.cube, cfg, subcfg, test13[cfg][6])
                    add_triangle(m, i, j, k, tiling13_5_1[cfg][subcfg], 6)
                else
                    add_triangle(m, i, j, k, tiling13_5_2[cfg][subcfg], 10)
                end
            elseif sc ≥ 27 && sc ≤ 38  # 13.3
                add_triangle(m, i, j, k, tiling13_3_[cfg][sc-26], 10, add_c_vertex(m, i, j, k))
            elseif sc ≥ 39 && sc ≤ 44  # 13.2
                add_triangle(m, i, j, k, tiling13_2_[cfg][sc-38], 6)
            elseif sc == 45  # 13.1
                add_triangle(m, i, j, k, tiling13_1_[cfg], 4)
            elseif sc > 45
                @error "Impossible case 13 ?"
            end
        elseif case == 15
            add_triangle(m, i, j, k, tiling14[cfg], 4)
        end
    end

    denormalize(m)
    return
end

"""
    compute_intersection_points(m::MC, vol, cb, iso)

# Description
Computes almost all the vertices of the mesh by interpolation along the cubes edges.
"""
compute_intersection_points(m::MC{F}, vol, cb, iso) where {F} = begin
    @inbounds for k ∈ axes(vol, 3), j ∈ axes(vol, 2), i ∈ axes(vol, 1)
        c0 = vol[i, j, k] - iso
        c1 = i < m.nx ? vol[i+1, j, k] - iso : c0
        c2 = j < m.ny ? vol[i, j+1, k] - iso : c0
        c3 = k < m.nz ? vol[i, j, k+1] - iso : c0
        m.cube[1] = abs(c0) < eps(F) ? eps(F) : c0
        m.cube[2] = abs(c1) < eps(F) ? eps(F) : c1
        m.cube[4] = abs(c2) < eps(F) ? eps(F) : c2
        m.cube[5] = abs(c3) < eps(F) ? eps(F) : c3
        if m.cube[1] < 0
            m.cube[2] > 0 &&
                (m.vert_indices[1, i, j, k] = add_x_vertex(m, vol, cb, i, j, k))
            m.cube[4] > 0 &&
                (m.vert_indices[2, i, j, k] = add_y_vertex(m, vol, cb, i, j, k))
            m.cube[5] > 0 &&
                (m.vert_indices[3, i, j, k] = add_z_vertex(m, vol, cb, i, j, k))
        else
            m.cube[2] < 0 &&
                (m.vert_indices[1, i, j, k] = add_x_vertex(m, vol, cb, i, j, k))
            m.cube[4] < 0 &&
                (m.vert_indices[2, i, j, k] = add_y_vertex(m, vol, cb, i, j, k))
            m.cube[5] < 0 &&
                (m.vert_indices[3, i, j, k] = add_z_vertex(m, vol, cb, i, j, k))
        end
    end
    return
end

"""
    test_interior(case, cb, cfg, subcfg, s)

# Description
Tests if the components of the tesselation of the cube should be connected through the interior of the cube.
"""
test_interior(case, cb::MVector{N,T}, cfg, subcfg, s) where {N,T} = begin
    test = 0
    At = Bt = Ct = Dt = T(0)
    @inbounds begin
        if case == 5 || case == 11
            a = (cb[5] - cb[1]) * (cb[7] - cb[3]) - (cb[8] - cb[4]) * (cb[6] - cb[2])
            b = (
                cb[3] * (cb[5] - cb[1]) + cb[1] * (cb[7] - cb[3]) -
                cb[2] * (cb[8] - cb[4]) - cb[4] * (cb[6] - cb[2])
            )
            t = -b / 2a
            (t < 0 || t > 1) && return s > 0

            At = cb[1] + (cb[5] - cb[1]) * t
            Bt = cb[4] + (cb[8] - cb[4]) * t
            Ct = cb[3] + (cb[7] - cb[3]) * t
            Dt = cb[2] + (cb[6] - cb[2]) * t
        elseif case == 7 || case == 8 || case == 13 || case == 14
            # reference edge of the triangulation
            edge = if case == 7
                test6[cfg][3]
            elseif case == 8
                test7[cfg][5]
            elseif case == 13
                test12[cfg][4]
            elseif case == 14
                tiling13_5_1[cfg][subcfg][1]
            else
                0
            end
            if edge == 1
                t = cb[1] / (cb[1] - cb[2])
                Bt = cb[4] + (cb[3] - cb[4]) * t
                Ct = cb[8] + (cb[7] - cb[8]) * t
                Dt = cb[5] + (cb[6] - cb[5]) * t
            elseif edge == 2
                t = cb[2] / (cb[2] - cb[3])
                Bt = cb[1] + (cb[4] - cb[1]) * t
                Ct = cb[5] + (cb[8] - cb[5]) * t
                Dt = cb[6] + (cb[7] - cb[6]) * t
            elseif edge == 3
                t = cb[3] / (cb[3] - cb[4])
                Bt = cb[2] + (cb[1] - cb[2]) * t
                Ct = cb[6] + (cb[5] - cb[6]) * t
                Dt = cb[7] + (cb[8] - cb[7]) * t
            elseif edge == 4
                t = cb[4] / (cb[4] - cb[1])
                Bt = cb[3] + (cb[2] - cb[3]) * t
                Ct = cb[7] + (cb[6] - cb[7]) * t
                Dt = cb[8] + (cb[5] - cb[8]) * t
            elseif edge == 5
                t = cb[5] / (cb[5] - cb[6])
                Bt = cb[8] + (cb[7] - cb[8]) * t
                Ct = cb[4] + (cb[3] - cb[4]) * t
                Dt = cb[1] + (cb[2] - cb[1]) * t
            elseif edge == 6
                t = cb[6] / (cb[6] - cb[7])
                Bt = cb[5] + (cb[8] - cb[5]) * t
                Ct = cb[1] + (cb[4] - cb[1]) * t
                Dt = cb[2] + (cb[3] - cb[2]) * t
            elseif edge == 7
                t = cb[7] / (cb[7] - cb[6])
                Bt = cb[6] + (cb[5] - cb[6]) * t
                Ct = cb[2] + (cb[1] - cb[2]) * t
                Dt = cb[3] + (cb[4] - cb[3]) * t
            elseif edge == 8
                t = cb[8] / (cb[8] - cb[5])
                Bt = cb[7] + (cb[6] - cb[7]) * t
                Ct = cb[3] + (cb[2] - cb[3]) * t
                Dt = cb[4] + (cb[1] - cb[4]) * t
            elseif edge == 9
                t = cb[1] / (cb[1] - cb[5])
                Bt = cb[4] + (cb[8] - cb[4]) * t
                Ct = cb[3] + (cb[7] - cb[3]) * t
                Dt = cb[2] + (cb[6] - cb[2]) * t
            elseif edge == 10
                t = cb[2] / (cb[2] - cb[6])
                Bt = cb[1] + (cb[5] - cb[1]) * t
                Ct = cb[4] + (cb[8] - cb[4]) * t
                Dt = cb[3] + (cb[7] - cb[3]) * t
            elseif edge == 11
                t = cb[3] / (cb[3] - cb[7])
                Bt = cb[2] + (cb[6] - cb[2]) * t
                Ct = cb[1] + (cb[5] - cb[1]) * t
                Dt = cb[4] + (cb[8] - cb[4]) * t
            elseif edge == 12
                t = cb[4] / (cb[4] - cb[8])
                Bt = cb[3] + (cb[7] - cb[3]) * t
                Ct = cb[2] + (cb[6] - cb[2]) * t
                Dt = cb[1] + (cb[5] - cb[1]) * t
            else
                @error "Invalid edge $edge"
            end
        else
            @error "Invalid ambiguous case $case"
        end
    end

    At ≥ 0 && (test += 1)
    Bt ≥ 0 && (test += 2)
    Ct ≥ 0 && (test += 4)
    Dt ≥ 0 && (test += 8)

    return if test == 6 || test == 8 || test == 9 || test == 12 || (test ≥ 0 && test ≤ 4)
        s > 0
    elseif test == 5
        (At * Ct - Bt * Dt < eps(T)) ? s > 0 : s < 0
    elseif test == 10
        (At * Ct - Bt * Dt ≥ eps(T)) ? s > 0 : s < 0
    else  # test == 7 || test == 11 || test ≥ 13
        s < 0
    end
end

"""
    test_face(cb, face)

# Description
Tests if the components of the tesselation of the cube should be connected by the interior of an ambiguous face.
"""
test_face(cb::MVector{N,T}, face) where {N,T} = begin
    @inbounds begin
        if (af = abs(face)) == 1
            A = cb[1]
            B = cb[5]
            C = cb[6]
            D = cb[2]
        elseif af == 2
            A = cb[2]
            B = cb[6]
            C = cb[7]
            D = cb[3]
        elseif af == 3
            A = cb[3]
            B = cb[7]
            C = cb[8]
            D = cb[4]
        elseif af == 4
            A = cb[4]
            B = cb[8]
            C = cb[5]
            D = cb[1]
        elseif af == 5
            A = cb[1]
            B = cb[4]
            C = cb[3]
            D = cb[2]
        elseif af == 6
            A = cb[5]
            B = cb[8]
            C = cb[7]
            D = cb[6]
        else
            @error "Invalid face code $face"
        end
    end
    abs(A * C - B * D) < eps(T) ? face ≥ 0 : face * A * (A * C - B * D) ≥ 0  # face and A invert signs
end

"""
    add_triangle(m::MC, i, j, k, tri, n, v12 = 0)

# Description
Routine to add a triangle to the mesh.

# Arguments
  - `trig` the code for the triangle as a sequence of edges index.
  - `n` the number of triangles to produce.
  - `v12` the index of the interior vertex to use, if necessary.
"""
add_triangle(m::MC, i, j, k, tri, n, v12 = 0) = begin
    @inbounds for t ∈ 1:3n
        tr = tri[t]
        id = (t - 1) % 3 + 1
        m.tv[id] = tv = if tr == 1
            m.vert_indices[1, i, j, k]
        elseif tr == 2
            m.vert_indices[2, i+1, j, k]
        elseif tr == 3
            m.vert_indices[1, i, j+1, k]
        elseif tr == 4
            m.vert_indices[2, i, j, k]
        elseif tr == 5
            m.vert_indices[1, i, j, k+1]
        elseif tr == 6
            m.vert_indices[2, i+1, j, k+1]
        elseif tr == 7
            m.vert_indices[1, i, j+1, k+1]
        elseif tr == 8
            m.vert_indices[2, i, j, k+1]
        elseif tr == 9
            m.vert_indices[3, i, j, k]
        elseif tr == 10
            m.vert_indices[3, i+1, j, k]
        elseif tr == 11
            m.vert_indices[3, i+1, j+1, k]
        elseif tr == 12
            m.vert_indices[3, i, j+1, k]
        else
            v12
        end
        tv == 0 && @warn "Invalid triangle $(length(m.triangles) + 1)"
        id == 3 && push!(m.triangles, Triangle(m.tv))
    end
    return
end

"""
    ∇x(args...; kwargs...)

# Description
Interpolates the horizontal gradient of the implicit function at the lower vertex of the specified cube.
"""
@inline ∇x(vol, i, j, k, nx) = @inbounds(
    i > 1 ?
    (i < nx ? (vol[i+1, j, k] - vol[i-1, j, k]) / 2 : (vol[i, j, k] - vol[i-1, j, k])) :
    (vol[i+1, j, k] - vol[i, j, k])
)

@inline ∇y(vol, i, j, k, ny) = @inbounds(
    j > 1 ?
    (j < ny ? (vol[i, j+1, k] - vol[i, j-1, k]) / 2 : (vol[i, j, k] - vol[i, j-1, k])) :
    (vol[i, j+1, k] - vol[i, j, k])
)

@inline ∇z(vol, i, j, k, nz) = @inbounds(
    k > 1 ?
    (k < nz ? (vol[i, j, k+1] - vol[i, j, k-1]) / 2 : (vol[i, j, k] - vol[i, j, k-1])) :
    (vol[i, j, k+1] - vol[i, j, k])
)

@inline norm(x) = @inbounds √(x[1]^2 + x[2]^2 + x[3]^2)

"""
    add_x_vertex(m::MC, vol, cb, i, j, k)

# Description
Adds a vertex on the current horizontal edge.
"""
add_x_vertex(m::MC{F}, vol, cb, i, j, k) where {F} = begin
    @inbounds begin
        u = cb[1] / (cb[1] - cb[2])
        m.nrm[1] = (1 - u) * ∇x(vol, i, j, k, m.nx) + u * ∇x(vol, i + 1, j, k, m.nx)
        m.nrm[2] = (1 - u) * ∇y(vol, i, j, k, m.ny) + u * ∇y(vol, i + 1, j, k, m.ny)
        m.nrm[3] = (1 - u) * ∇z(vol, i, j, k, m.nz) + u * ∇z(vol, i + 1, j, k, m.nz)
        (mag = norm(m.nrm)) > eps(F) && (m.nrm ./= mag)
    end
    push!(m.vertices, Vertex(i - 1 + u, j - 1, k - 1))
    push!(m.normals, Normal(m.nrm))
    length(m.vertices)
end

add_y_vertex(m::MC{F}, vol, cb, i, j, k) where {F} = begin
    @inbounds begin
        u = cb[1] / (cb[1] - cb[4])
        m.nrm[1] = (1 - u) * ∇x(vol, i, j, k, m.nx) + u * ∇x(vol, i, j + 1, k, m.nx)
        m.nrm[2] = (1 - u) * ∇y(vol, i, j, k, m.ny) + u * ∇y(vol, i, j + 1, k, m.ny)
        m.nrm[3] = (1 - u) * ∇z(vol, i, j, k, m.nz) + u * ∇z(vol, i, j + 1, k, m.nz)
        (mag = norm(m.nrm)) > eps(F) && (m.nrm ./= mag)
    end
    push!(m.vertices, Vertex(i - 1, j - 1 + u, k - 1))
    push!(m.normals, Normal(m.nrm))
    length(m.vertices)
end

add_z_vertex(m::MC{F}, vol, cb, i, j, k) where {F} = begin
    @inbounds begin
        u = cb[1] / (cb[1] - cb[5])
        m.nrm[1] = (1 - u) * ∇x(vol, i, j, k, m.nx) + u * ∇x(vol, i, j, k + 1, m.nx)
        m.nrm[2] = (1 - u) * ∇y(vol, i, j, k, m.ny) + u * ∇y(vol, i, j, k + 1, m.ny)
        m.nrm[3] = (1 - u) * ∇z(vol, i, j, k, m.nz) + u * ∇z(vol, i, j, k + 1, m.nz)
        (mag = norm(m.nrm)) > eps(F) && (m.nrm ./= mag)
    end
    push!(m.vertices, Vertex(i - 1, j - 1, k - 1 + u))
    push!(m.normals, Normal(m.nrm))
    length(m.vertices)
end

"""
    add_c_vertex(m::MC, i, j, k)

# Description
Add a vertex inside the current cube.
"""
add_c_vertex(m::MC{F}, i, j, k) where {F} = begin
    u = 0
    @inbounds begin
        m.vtx .= 0
        m.nrm .= 0
        # compute the average of the intersection points of the cube
        (n = m.vert_indices[1, i, j, k]) > 0 &&
            (u += 1; m.vtx .+= m.vertices[n]; m.nrm .+= m.normals[n])
        (n = m.vert_indices[2, i+1, j, k]) > 0 &&
            (u += 1; m.vtx .+= m.vertices[n]; m.nrm .+= m.normals[n])
        (n = m.vert_indices[1, i, j+1, k]) > 0 &&
            (u += 1; m.vtx .+= m.vertices[n]; m.nrm .+= m.normals[n])
        (n = m.vert_indices[2, i, j, k]) > 0 &&
            (u += 1; m.vtx .+= m.vertices[n]; m.nrm .+= m.normals[n])
        (n = m.vert_indices[1, i, j, k+1]) > 0 &&
            (u += 1; m.vtx .+= m.vertices[n]; m.nrm .+= m.normals[n])
        (n = m.vert_indices[2, i+1, j, k+1]) > 0 &&
            (u += 1; m.vtx .+= m.vertices[n]; m.nrm .+= m.normals[n])
        (n = m.vert_indices[1, i, j+1, k+1]) > 0 &&
            (u += 1; m.vtx .+= m.vertices[n]; m.nrm .+= m.normals[n])
        (n = m.vert_indices[2, i, j, k+1]) > 0 &&
            (u += 1; m.vtx .+= m.vertices[n]; m.nrm .+= m.normals[n])
        (n = m.vert_indices[3, i, j, k]) > 0 &&
            (u += 1; m.vtx .+= m.vertices[n]; m.nrm .+= m.normals[n])
        (n = m.vert_indices[3, i+1, j, k]) > 0 &&
            (u += 1; m.vtx .+= m.vertices[n]; m.nrm .+= m.normals[n])
        (n = m.vert_indices[3, i+1, j+1, k]) > 0 &&
            (u += 1; m.vtx .+= m.vertices[n]; m.nrm .+= m.normals[n])
        (n = m.vert_indices[3, i, j+1, k]) > 0 &&
            (u += 1; m.vtx .+= m.vertices[n]; m.nrm .+= m.normals[n])
        m.vtx ./= u
        (mag = norm(m.nrm)) > eps(F) && (m.nrm ./= mag)
    end
    push!(m.vertices, Vertex(m.vtx))
    push!(m.normals, Normal(m.nrm))
    length(m.vertices)
end

include("example.jl")

end


text:
using MultivariateStochasticVolatility
using Test

@testset "MultivariateStochasticVolatility.jl" begin
    # Write your tests here.
end