text:
using DataStreams, DataFrames

import Parsers: readbyte, peekbyte

substitute(::Type{Union{T, Missing}}, ::Type{T1}) where {T, T1} = Union{T1, Missing}
substitute(::Type{T}, ::Type{T1}) where {T, T1} = T1
substitute(::Type{Missing}, ::Type{T1}) where {T1} = Missing

"""
Represents the various configuration settings for delimited text file parsing.

Keyword Arguments:

 * `delim::Union{Char,UInt8}`: how fields in the file are delimited; default `','`
 * `quotechar::Union{Char,UInt8}`: the character that indicates a quoted field that may contain the `delim` or newlines; default `'"'`
 * `escapechar::Union{Char,UInt8}`: the character that escapes a `quotechar` in a quoted field; default `'\\'`
 * `missingstring::String`: indicates how missing values are represented in the dataset; default `""`
 * `dateformat::Union{AbstractString,Dates.DateFormat}`: how dates/datetimes are represented in the dataset; default `Base.Dates.ISODateTimeFormat`
 * `decimal::Union{Char,UInt8}`: character to recognize as the decimal point in a float number, e.g. `3.14` or `3,14`; default `'.'`
 * `truestring`: string to represent `true::Bool` values in a csv file; default `"true"`. Note that `truestring` and `falsestring` cannot start with the same character.
 * `falsestring`: string to represent `false::Bool` values in a csv file; default `"false"`
 * `internstrings`: whether strings should be interned, rather than creating a new object for each string field; default `true`
"""
struct Options{D}
    delim::UInt8
    quotechar::UInt8
    escapechar::UInt8
    missingstring::Vector{UInt8}
    null::Union{Vector{UInt8},Nothing} # deprecated
    missingcheck::Bool
    dateformat::D
    decimal::UInt8
    truestring::Vector{UInt8}
    falsestring::Vector{UInt8}
    internstrings::Bool
    # non-public for now
    datarow::Int
    rows::Int
    header::Union{Integer,UnitRange{Int},Vector}
    types
end

function Options(;delim=UInt8(','), quotechar=UInt8('"'), escapechar=UInt8('\\'), missingstring="", null=nothing, dateformat=nothing, decimal=UInt8('.'), truestring="true", falsestring="false", internstrings=true, datarow=-1, rows=0, header=1, types=Type[])
    Base.depwarn("CSV.Options is deprecated and will be removed in a future release; please use CSV.File and CSV.write directly", nothing)
    return Options(delim%UInt8, quotechar%UInt8, escapechar%UInt8,
            map(UInt8, collect(ascii(String(missingstring)))), null === nothing ? nothing : map(UInt8, collect(ascii(String(null)))), missingstring != "" || (null != "" && null != nothing),
            isa(dateformat, AbstractString) ? Dates.DateFormat(dateformat) : dateformat,
            decimal%UInt8, map(UInt8, collect(truestring)), map(UInt8, collect(falsestring)), internstrings, datarow, rows, header, types)
end
function Base.show(io::IO,op::Options)
    println(io, "    CSV.Options:")
    println(io, "        delim: '", Char(op.delim), "'")
    println(io, "        quotechar: '", Char(op.quotechar), "'")
    print(io, "        escapechar: '"); escape_string(io, string(Char(op.escapechar)), "\\"); println(io, "'")
    print(io, "        missingstring: \""); escape_string(io, isempty(op.missingstring) ? "" : String(collect(op.missingstring)), "\\"); println(io, "\"")
    println(io, "        dateformat: ", something(op.dateformat, ""))
    println(io, "        decimal: '", Char(op.decimal), "'")
    println(io, "        truestring: '$(String(op.truestring))'")
    println(io, "        falsestring: '$(String(op.falsestring))'")
    print(io, "        internstrings: $(op.internstrings)")
end

"""
A type that satisfies the `Data.Source` interface in the `DataStreams.jl` package.

A `CSV.Source` can be manually constructed in order to be re-used multiple times.

`CSV.Source(file_or_io; kwargs...) => CSV.Source`

Note that a filename string can be provided or any `IO` type. For the full list of supported
keyword arguments, see the docs for [`CSV.read`](@ref) or type `?CSV.read` at the repl

An example of re-using a `CSV.Source` is:
```julia
# manually construct a `CSV.Source` once, then stream its data to both a DataFrame
# and SQLite table `sqlite_table` in the SQLite database `db`
# note the use of `CSV.reset!` to ensure the `source` can be streamed from again
source = CSV.Source(file)
df1 = CSV.read(source, DataFrame)
CSV.reset!(source)
sq1 = CSV.read(source, SQLite.Sink, db, "sqlite_table")
```
"""
mutable struct Source{P, I, DF, D} <: Data.Source
    schema::Data.Schema
    parsinglayers::P
    io::I
    bools::Parsers.Trie
    dateformat::DF
    decimal::UInt8
    fullpath::String
    datapos::D # the position in the IOBuffer where the rows of data begins or columnpositions for transposed files
    pools::Vector{CategoricalPool{String, UInt32, CategoricalString{UInt32}}}
end

function Base.show(io::IO, f::Source)
    println(io, "CSV.Source: ", f.fullpath)
    println(io, f.io)
    show(io, f.schema)
end

"""
A type that satisfies the `Data.Sink` interface in the `DataStreams.jl` package.

A `CSV.Sink` can be manually constructed in order to be re-used multiple times.

`CSV.Sink(file_or_io; kwargs...) => CSV.Sink`

Note that a filename string can be provided or any `IO` type. For the full list of supported
keyword arguments, see the docs for [`CSV.write`](@ref) or type `?CSV.write` at the repl

An example of re-using a `CSV.Sink` is:
```julia
# manually construct a `CSV.Source` once, then stream its data to both a DataFrame
# and SQLite table `sqlite_table` in the SQLite database `db`
# note the use of `CSV.reset!` to ensure the `source` can be streamed from again
source = CSV.Source(file)
df1 = CSV.read(source, DataFrame)
CSV.reset!(source)
sq1 = CSV.read(source, SQLite.Sink, db, "sqlite_table")
```
"""
mutable struct Sink{D, B} <: Data.Sink
    options::Options{D}
    io::IOBuffer
    fullpath::Union{String, IO}
    datapos::Int # the position in the IOBuffer where the rows of data begins
    header::Bool
    colnames::Vector{String}
    cols::Int
    append::Bool
    quotefields::B
end

makedf(df::DateFormat) = df
makedf(df::String) = DateFormat(df)
makedf(df::Nothing) = df
makestr(d::UInt8) = string(Char(d))
makestr(d::Char) = string(d)
makestr(d::String) = d

readline!(io::IO; delim=",") = readline!(Parsers.Delimited(Parsers.Quoted(), delim; newline=true), io)
readline!(s::Source) = readline!(s.parsinglayers, s.io)
function countlines(io::IO, q::UInt8, e::UInt8)
    nl = 1
    b = 0x00
    while !eof(io)
        b = readbyte(io)
        if b === q
            while !eof(io)
                b = readbyte(io)
                if b === e
                    if eof(io)
                        break
                    elseif e === q && peekbyte(io) != q
                        break
                    end
                    b = readbyte(io)
                elseif b === q
                    break
                end
            end
        elseif b === UInt8('\n')
            nl += 1
        elseif b === UInt8('\r')
            nl += 1
            !eof(io) && peekbyte(io) === UInt8('\n') && readbyte(io)
        end
    end
    return nl - (b === UInt8('\n') || b === UInt8('\r'))
end
countlines(io::IO, q='"', e='\\') = countlines(io, UInt8(q), UInt8(e))

promote_type1(T::Type{<:Any}, ::Type{Any}) = T
promote_type1(::Type{Any}, T::Type{<:Any}) = T
# same types
promote_type1(::Type{T}, ::Type{T}) where {T} = T
# if we come across a Missing field, turn that column type into a Union{T, Missing}
promote_type1(T::Type{<:Any}, ::Type{Missing}) = Union{T, Missing}
promote_type1(::Type{Missing}, T::Type{<:Any}) = Union{T, Missing}
# these definitions allow Union{Int64, Missing} to promote to Union{Float64, Missing}
promote_type1(::Type{Union{T, Missing}}, ::Type{S}) where {T, S} = Union{promote_type1(T, S), Missing}
promote_type1(::Type{S}, ::Type{Union{T, Missing}}) where {T, S} = Union{promote_type1(T, S), Missing}
promote_type1(::Type{Union{T, Missing}}, ::Type{Union{S, Missing}}) where {T, S} = Union{promote_type1(T, S), Missing}
promote_type1(::Type{Union{String, Missing}}, ::Type{String}) = Union{String, Missing}
promote_type1(::Type{String}, ::Type{Union{String, Missing}}) = Union{String, Missing}
# basic promote type definitions from Base
promote_type1(::Type{Int64}, ::Type{Float64}) = Float64
promote_type1(::Type{Float64}, ::Type{Int64}) = Float64
promote_type1(::Type{Date}, ::Type{DateTime}) = DateTime
promote_type1(::Type{DateTime}, ::Type{Date}) = DateTime
# for cases when our current type can't widen, just promote to WeakRefString
promote_type1(::Type{<:Real}, ::Type{<:Dates.TimeType}) = String
promote_type1(::Type{<:Dates.TimeType}, ::Type{<:Real}) = String
promote_type1(::Type{T}, ::Type{String}) where T = String
promote_type1(::Type{Union{T, Missing}}, ::Type{String}) where T = Union{String, Missing}
promote_type1(::Type{String}, ::Type{T}) where T = String
promote_type1(::Type{String}, ::Type{Union{T, Missing}}) where T = Union{String, Missing}
# avoid ambiguity
promote_type1(::Type{Any}, ::Type{String}) = String
promote_type1(::Type{String}, ::Type{Any}) = String
promote_type1(::Type{String}, ::Type{String}) = String
promote_type1(::Type{String}, ::Type{Missing}) = Union{String, Missing}
promote_type1(::Type{Missing}, ::Type{String}) = Union{String, Missing}
promote_type1(::Type{Any}, ::Type{Missing}) = Missing
promote_type1(::Type{Missing}, ::Type{Missing}) = Missing

function detecttype(layers::Parsers.Delimited, io::IO, prevT, levels, row, col, bools, dateformat, dec, old)
    pos = position(io)
    result = Parsers.parse(layers, io, Tuple{Ptr{UInt8}, Int})
    Parsers.ok(result.code) || throw(Error(Parsers.Error(io, result), row, col))
    # @debug "res = $result"
    result.result === missing && return Missing
    # update levels
    incr!(levels, result.result::Tuple{Ptr{UInt8}, Int})

    if Int64 <: prevT || prevT == Missing
        seek(io, pos)
        res_int = Parsers.parse(layers, io, Int64)
        # @debug "res_int = $res_int"
        Parsers.ok(res_int.code) && return Int64
    end
    if Float64 <: prevT || Int64 <: prevT || prevT == Missing
        seek(io, pos)
        res_float = Parsers.parse(layers, io, Float64; decimal=dec)
        # @debug "res_float = $res_float"
        Parsers.ok(res_float.code) && return Float64
    end
    if Date <: prevT || DateTime <: prevT || prevT == Missing
        if dateformat === nothing
            # try to auto-detect TimeType
            seek(io, pos)
            res_date = Parsers.parse(layers, io, Date)
            # @debug "res_date = $res_date"
            Parsers.ok(res_date.code) && return Date
            seek(io, pos)
            res_datetime = Parsers.parse(layers, io, DateTime)
            # @debug "res_datetime = $res_datetime"
            Parsers.ok(res_datetime.code) && return DateTime
        else
            # use user-provided dateformat
            T = timetype(dateformat)
            # @debug "T = $T"
            seek(io, pos)
            res_dt = Parsers.parse(layers, io, T; dateformat=dateformat)
            # @debug "res_dt = $res_dt"
            Parsers.ok(res_dt.code) && return T
        end
    end
    if Bool <: prevT || prevT == Missing
        seek(io, pos)
        res_bool = Parsers.parse(layers, io, Bool; bools=bools)
        # @debug "res_bool = $res_bool"
        Parsers.ok(res_bool.code) && return Bool
    end
    return String
end

function datalayout(header, source, parsinglayers, delim, quotechar, escapechar, datarow, footerskip, rows, startpos, fs)
    rows = rows == 0 ? CSV.countlines(source, quotechar, escapechar) : rows
    seek(source, startpos)
    rows = fs == 0 ? -1 : max(-1, rows - datarow + 1 - footerskip) # rows now equals the actual number of rows in the dataset

    # figure out # of columns and header, either an Integer, AbstractRange, or Vector{String}
    # also ensure that `source` is positioned at the start of data
    if isa(header, Integer)
        # default header = 1
        if header <= 0
            # no header row in dataset; skip to data to figure out # of columns
            CSV.skipto!(parsinglayers, source, 1, datarow)
            datapos = position(source)
            row_vals = readsplitline(parsinglayers, source)
            seek(source, datapos)
            columnnames = ["Column$i" for i = eachindex(row_vals)]
        else
            CSV.skipto!(parsinglayers, source, 1, header)
            columnnames = [ismissing(x) ? "" : strip(x) for x in readsplitline(parsinglayers, source)]
            datarow != header+1 && CSV.skipto!(parsinglayers, source, header+1, datarow)
            datapos = position(source)
        end
    elseif isa(header, AbstractRange)
        CSV.skipto!(parsinglayers, source, 1, first(header))
        columnnames = [x for x in readsplitline(parsinglayers, source)]
        for row = first(header):(last(header)-1)
            for (i,c) in enumerate([x for x in readsplitline(parsinglayers, source)])
                columnnames[i] *= "_" * c
            end
        end
        datarow != last(header)+1 && CSV.skipto!(parsinglayers, source, last(header)+1, datarow)
        datapos = position(source)
    elseif fs == 0
        datapos = position(source)
        columnnames = header
        cols = length(columnnames)
    else
        CSV.skipto!(parsinglayers, source, 1, datarow)
        datapos = position(source)
        row_vals = readsplitline(parsinglayers, source)
        seek(source, datapos)
        if isempty(header)
            columnnames = ["Column$i" for i in eachindex(row_vals)]
        else
            length(header) == length(row_vals) || throw(ArgumentError("The length of provided header ($(length(header))) doesn't match the number of columns at row $datarow ($(length(row_vals)))"))
            columnnames = header
        end
    end
    return rows, columnnames, datapos
end

function datalayout_transpose(header, source, parsinglayers, delim, quotechar, escapechar, datarow, footerskip, fs)
    if isa(header, Integer) && header > 0
        # skip to header column to read column names
        row = skiptoheader!(parsinglayers, source, 1, header)
        # source now at start of 1st header cell
        columnnames = [strip(Parsers.parse(parsinglayers, source, String).result::String)]
        columnpositions = [position(source)]
        datapos = position(source)
        rows = countfields(source, parsinglayers)

        # we're now done w/ column 1, if EOF we're done, otherwise, parse column 2's column name
        cols = 1
        while !eof(source)
            # skip to header column to read column names
            row = skiptoheader!(parsinglayers, source, 1, header)
            cols += 1
            push!(columnnames, strip(Parsers.parse(parsinglayers, source, String).result::String))
            push!(columnpositions, position(source))
            readline!(parsinglayers, source)
        end
        seek(source, datapos)
    elseif isa(header, AbstractRange)
        # column names span several columns
        throw(ArgumentError("not implemented for transposed csv files"))
    elseif fs == 0
        # emtpy file, use column names if provided
        datapos = position(source)
        columnpositions = Int[]
        columnnames = header
    else
        # column names provided explicitly or should be generated, they don't exist in data
        # skip to datarow
        row = skiptoheader!(parsinglayers, source, 1, datarow)
        # source now at start of 1st data cell
        columnnames = [isa(header, Integer) || isempty(header) ? "Column1" : header[1]]
        columnpositions = [position(source)]
        datapos = position(source)
        rows = countfields(source, parsinglayers)
        # we're now done w/ column 1, if EOF we're done, otherwise, parse column 2's column name
        cols = 1
        while !eof(source)
            # skip to datarow column
            row = skiptoheader!(parsinglayers, source, 1, datarow)
            cols += 1
            push!(columnnames, isa(header, Integer) || isempty(header) ? "Column$cols" : header[cols])
            push!(columnpositions, position(source))
            readline!(parsinglayers, source)
        end
        seek(source, datapos)
    end
    rows = rows - footerskip # rows now equals the actual number of rows in the dataset
    return rows, columnnames, columnpositions
end

isascii(c) = Base.isascii(c)
isascii(c::UInt8) = c < 0x80

function Source(fullpath::Union{AbstractString,IO};
              delim=",",
              quotechar='"',
              escapechar='\\',
              missingstring=nothing,

              header::Union{Integer, UnitRange{Int}, Vector}=1, # header can be a row number, range of rows, or actual string vector
              datarow::Int=-1, # by default, data starts immediately after header or start of file
              types=Type[],
              typemap::Dict{Type, Type}=Dict{Type, Type}(),
              allowmissing::Symbol=:all,
              dateformat::Union{String, Dates.DateFormat, Nothing}=nothing,
              decimal='.',
              truestring=nothing,
              falsestring=nothing,
              categorical::Bool=true,
              strings::Union{Nothing, Symbol}=nothing,

              footerskip::Int=0,
              rows_for_type_detect::Int=100,
              rows::Int=0,
              limit::Union{Nothing, Int}=nothing,
              use_mmap::Bool=true,
              transpose::Bool=false,
              warn::Bool=true)
    # make sure character args are UInt8
    warn && Base.depwarn("CSV.Source is deprecated and will be removed in a future release; please use CSV.File instead (supports most of the same options)", nothing)
    isascii(delim) || throw(ArgumentError("non-ASCII characters not supported for delim argument: $delim"))
    isascii(quotechar) || throw(ArgumentError("non-ASCII characters not supported for quotechar argument: $quotechar"))
    isascii(escapechar) || throw(ArgumentError("non-ASCII characters not supported for escapechar argument: $escapechar"))

    # deprecations for CSV.File
    if strings === nothing
        strings = :intern # no need to depwarn
    elseif strings === :intern
        Base.depwarn("strings=:intern is deprecated, String values are interned by default now", nothing)
    elseif strings === :weakref
        Base.depwarn("strings=:weakref is deprecated, use typemap=Dict(String=>WeakRefString{UInt8}) instead", nothing)
        merge!(typemap, Dict(String=>WeakRefString{UInt8}))
    end
    if rows_for_type_detect != 100
        Base.depwarn("rows_for_type_detect is deprecated; CSV.File will randomly sample throughout an entire file to better infer column types", nothing)
    end
    if rows != 0
        Base.depwarn("the rows keyword argument is deprecated; to limit the number of rows read from a csv file, pass `limit=X` instead", nothing)
        limit = rows
    end
    if limit === nothing
        limit = 0
    end
    if truestring !== nothing
        Base.depwarn("truestring is deprecated in CSV.File in favor of allowing the passing of multiple strings via `truestrings=[...]`", nothing)
    else
        truestring = "true"
    end
    if falsestring !== nothing
        Base.depwarn("falsestring is deprecated in CSV.File in favor of allowing the passing of multiple strings via `falsestrings=[...]`", nothing)
    else
        falsestring = "false"
    end
    if missingstring !== nothing
        Base.depwarn("missingstring is deprecated in CSV.File in favor of allowing the passing of multiple strings via `missingstrings=[...]`", nothing)
    else
        missingstring = ""
    end

    # argument checks
    isa(fullpath, AbstractString) && (isfile(fullpath) || throw(ArgumentError("\"$fullpath\" is not a valid file")))
    header = (isa(header, Integer) && header == 1 && datarow == 1) ? -1 : header
    isa(header, Integer) && datarow != -1 && (datarow > header || throw(ArgumentError("data row ($datarow) must come after header row ($header)")))

    # open the file for property detection
    if isa(fullpath, IOBuffer)
        source = fullpath
        fs = bytesavailable(fullpath)
        fullpath = "<IOBuffer>"
    elseif isa(fullpath, IO)
        source = IOBuffer(Base.read(fullpath))
        fs = bytesavailable(source)
        fullpath = isdefined(fullpath, :name) ? fullpath.name : "__IO__"
    else
        source = open(fullpath, "r") do f
            IOBuffer(use_mmap ? Mmap.mmap(f) : Base.read(f))
        end
        fs = filesize(fullpath)
    end

    startpos = position(source)

    # BOM character detection
    if fs > 0 && Parsers.peekbyte(source) == 0xef
        Parsers.readbyte(source)
        Parsers.readbyte(source) == 0xbb || seek(source, startpos)
        Parsers.readbyte(source) == 0xbf || seek(source, startpos)
        startpos = position(source)
    end

    datarow = datarow == -1 ? (isa(header, Vector) ? 0 : last(header)) + 1 : datarow # by default, data starts on line after header

    # build Parsers parsinglayers
    bools = (truestring != "true" || falsestring != "false") ? Parsers.Trie([truestring=>true, falsestring=>false]) : Parsers.BOOLS
    df = makedf(dateformat)
    dec = decimal % UInt8
    d = makestr(delim)
    parsinglayers = Parsers.Delimited(Parsers.Quoted(Parsers.Strip(Parsers.Sentinel(missingstring), d == " " ? 0x00 : ' ', d == "\t" ? 0x00 : '\t'), quotechar, escapechar), d; newline=true)

    # data layout detection: figure out rows, columnnames, and datapos
    if transpose
        rows, columnnames, datapos = datalayout_transpose(header, source, parsinglayers, delim, quotechar, escapechar, datarow, footerskip, fs)
        originalcolumnpositions = [x for x in datapos]
    else
        rows, columnnames, datapos = datalayout(header, source, parsinglayers, delim, quotechar, escapechar, datarow, footerskip, limit, startpos, fs)
    end
    cols = length(columnnames)
    @debug "columnnames=$columnnames, cols=$cols, rows=$rows, datapos=$datapos"

    # Detect column types
    pools = CategoricalPool{String, UInt32, CatStr}[]
    if isa(types, Vector) && length(types) == cols
        # types might be a Vector{DataType}, which will be a problem if Unions are needed
        columntypes = convert(Vector{Type}, types)
    elseif isa(types, Dict) || isempty(types)
        @debug "detecting types"
        columntypes = Type[Any for x = 1:cols]
        levels = [Dict{String, Int}() for _ = 1:cols]
        lineschecked = 0
        while !eof(source) && lineschecked < min(rows < 0 ? rows_for_type_detect : rows, rows_for_type_detect)
            lineschecked += 1
            @debug "type detecting on row = $lineschecked..."
            for i = 1:cols
                transpose && seek(source, datapos[i])
                @debug "\tdetecting col = $i..."
                typ = CSV.detecttype(parsinglayers, source, columntypes[i], levels[i], lineschecked, i, bools, df, dec, true)::Type
                transpose && setindex!(datapos, position(source), i)
                @debug "$typ"
                columntypes[i] = CSV.promote_type1(columntypes[i], typ)
                @debug "...promoting to: $(columntypes[i])"
            end
        end
        if df === nothing && any(x->Base.nonmissingtype(x) <: Dates.TimeType, columntypes)
            df = any(x->Base.nonmissingtype(x) <: DateTime, columntypes) ? Dates.default_format(DateTime) : Dates.default_format(Date)
        end
        if categorical
            pools = Vector{CategoricalPool{String, UInt32, CatStr}}(undef, cols)
            for i = 1:cols
                T = columntypes[i]
                if length(levels[i]) / sum(values(levels[i])) < .67 && T !== Missing && Base.nonmissingtype(T) <: String
                    columntypes[i] = substitute(T, CategoricalArrays.catvaluetype(Base.nonmissingtype(T), UInt32))
                    pools[i] = CategoricalPool{String, UInt32}(collect(keys(levels[i])))
                end
            end
        end
    else
        throw(ArgumentError("$cols number of columns detected; `types` argument has $(length(types)) entries"))
    end
    if isa(types, Dict)
        if isa(types, Dict{String})
            colinds = indexin(keys(types), columnnames)
        else
            colinds = keys(types)
        end
        for (col, typ) in zip(colinds, values(types))
            columntypes[col] = typ
        end
        autocols = setdiff(1:cols, colinds)
    elseif isempty(types)
        autocols = collect(1:cols)
    else
        autocols = Int[]
    end
    if !isempty(typemap)
        for col = 1:cols
            columntypes[col] = get(typemap, columntypes[col]) do
                if haskey(typemap, Base.nonmissingtype(columntypes[col]))
                    Union{Missing, typemap[Base.nonmissingtype(columntypes[col])]}
                else
                    columntypes[col]
                end
            end
        end
    end
    if allowmissing != :auto
        if allowmissing == :all # allow missing values in all automatically detected columns
            for i = autocols
                T = columntypes[i]
                columntypes[i] = Union{Base.nonmissingtype(T), Missing}
            end
        elseif allowmissing == :none # disallow missing values in all automatically detected columns
            for i = autocols
                T = columntypes[i]
                columntypes[i] = Base.nonmissingtype(T)
            end
        else
            throw(ArgumentError("allowmissing must be either :all, :none or :auto"))
        end
    end
    sch = Data.Schema(columntypes, columnnames, ifelse(rows < 0, missing, rows))
    if transpose
        datapos = originalcolumnpositions
        seek(source, startpos)
    else
        seek(source, datapos)
    end
    return Source(sch, parsinglayers, source, bools, df, dec, String(fullpath), datapos, pools)
end

# Data.Source interface
"reset a `CSV.Source` to its beginning to be ready to parse data from again"
Data.reset!(s::CSV.Source) = (seek(s.io, s.datapos); return nothing)
Data.schema(source::CSV.Source) = source.schema
Data.accesspattern(::Type{<:CSV.Source}) = Data.Sequential
@inline Data.isdone(io::CSV.Source, row, col, rows, cols) = eof(io.io) || (!ismissing(rows) && row > rows)
@inline Data.isdone(io::Source, row, col) = Data.isdone(io, row, col, size(io.schema)...)
Data.streamtype(::Type{<:CSV.Source}, ::Type{Data.Field}) = true
Data.reference(source::CSV.Source) = source.io.data

@inline function parsefield(source::CSV.Source{P, I, DF, D}, ::Type{T}, row, col; kwargs...) where {P, I, DF, D, T}
    D === Vector{Int} && Parsers.fastseek!(source.io, source.datapos[col])
    r = Parsers.parse(source.parsinglayers, source.io, T; kwargs...)
    D === Vector{Int} && setindex!(source.datapos, position(source.io), col)
    if Parsers.ok(r.code)
        return r.result
    else
        throw(Error(Parsers.Error(source.io, r), row, col))
    end
end

@inline Data.streamfrom(source::Source, ::Type{Data.Field}, T, row, col::Int) = parsefield(source, Base.nonmissingtype(T), row, col)
@inline Data.streamfrom(source::Source, ::Type{Data.Field}, ::Union{Type{Bool}, Union{Bool, Missing}}, row, col::Int) = parsefield(source, Bool, row, col; bools=source.bools)
@inline Data.streamfrom(source::Source, ::Type{Data.Field}, ::Type{T}, row, col::Int) where {T <: Union{AbstractFloat, Missing}} = parsefield(source, Base.nonmissingtype(T), row, col; decimal=source.decimal)
@inline Data.streamfrom(source::Source, ::Type{Data.Field}, ::Type{T}, row, col::Int) where {T <: Union{Dates.TimeType, Missing}} = parsefield(source, Base.nonmissingtype(T), row, col; dateformat=source.dateformat)
@inline Data.streamfrom(source::Source, ::Type{Data.Field}, ::Type{Missing}, row, col::Int) = parsefield(source, Missing, row, col)

@inline function Data.streamfrom(source::Source, ::Type{Data.Field}, ::Union{Type{Union{CatStr, Missing}}, Type{CatStr}}, row, col::Int)
    str = parsefield(source, Tuple{Ptr{UInt8}, Int}, row, col)
    if str isa Missing
        return missing
    else
        @inbounds pool = source.pools[col]
        str::Tuple{Ptr{UInt8}, Int}
        i = get(pool, str, nothing)
        if i === nothing
            i = get!(pool, unsafe_string(str[1], str[2]))
            issorted(levels(pool)) || levels!(pool, sort(levels(pool)))
        end
        return CatStr(i, pool)
    end
end

function read(source::CSV.Source, sink::Union{Type, Nothing}=DataFrame, args...; append::Bool=false, transforms::Dict=Dict{Int,Function}())
    if append
        Base.depwarn("`CSV.read(source; append=true)` is deprecated in favor of sink-specific options; e.g. DataFrames supports `CSV.File(filename) |> x->append!(existing_df, x)` to append the rows of a csv file to an existing DataFrame", nothing)
    end
    if !isempty(transforms)
        Base.depwarn("`CSV.read(source; transforms=Dict(...)` is deprecated in favor of `CSV.File(filename) |> transform((a=x->x+1, b=x->string(\"custom_prefix\", x))) |> DataFrame`", nothing)
    end
    sink = Data.stream!(source, sink, args...; append=append, transforms=transforms)
    return Data.close!(sink)
end
function read(source::CSV.Source, sink::T; append::Bool=false, transforms::Dict=Dict{Int,Function}()) where {T}
    Base.depwarn("`CSV.read(source::CSV.Source, sink::$T)` is deprecated; use `CSV.read(source) |> sink` instead", nothing)
    if append
        Base.depwarn("`CSV.read(source; append=true)` is deprecated in favor of sink-specific options; e.g. DataFrames supports `CSV.File(filename) |> x->append!(existing_df, x)` to append the rows of a csv file to an existing DataFrame", nothing)
    end
    if !isempty(transforms)
        Base.depwarn("`CSV.read(source; transforms=Dict(...)` is deprecated in favor of `CSV.File(filename) |> transform((a=x->x+1, b=x->string(\"custom_prefix\", x))) |> DataFrame`", nothing)
    end
    sink = Data.stream!(source, sink; append=append, transforms=transforms)
    return Data.close!(sink)
end

function Sink(fullpath::Union{AbstractString, IO};
              delim::Char=',',
              quotechar::Char='"',
              escapechar::Char='\\',
              missingstring::AbstractString="",
              dateformat=nothing,
              header::Bool=true,
              colnames::Vector{String}=String[],
              append::Bool=false,
              quotefields::Bool=false)
    Base.depwarn("CSV.Sink is deprecated and will be removed in a future release; use `source |> CSV.write(file; kwargs...)` or `CSV.write(source, file; kwargs...)` directly.", nothing)
    delim = delim % UInt8; quotechar = quotechar % UInt8; escapechar = escapechar % UInt8
    dateformat = isa(dateformat, AbstractString) ? Dates.DateFormat(dateformat) : dateformat
    io = IOBuffer()
    options = CSV.Options(delim=delim, quotechar=quotechar, escapechar=escapechar, missingstring=missingstring, dateformat=dateformat)
    !append && header && !isempty(colnames) && writeheaders(io, colnames, options, Val{quotefields})
    return Sink(options, io, fullpath, position(io), !append && header && !isempty(colnames), colnames, length(colnames), append, Val{quotefields})
end

quoted(::Type{Val{true}},  val, q, e, d) =  string(q, replace(val, q=>string(e, q)), q)
quoted(::Type{Val{false}}, val, q, e, d) = (q in val || d in val) ? string(q, replace(val, q=>string(e, q)), q) : val

function writeheaders(io::IOBuffer, h::Vector{String}, options, quotefields)
    cols = length(h)
    q = Char(options.quotechar); e = Char(options.escapechar); d = Char(options.delim)
    for col = 1:cols
        Base.write(io, quoted(quotefields, h[col], q, e, d), ifelse(col == cols, UInt8('\n'), options.delim))
    end
    return nothing
end

# DataStreams interface
Data.streamtypes(::Type{CSV.Sink}) = [Data.Field]
Data.weakrefstrings(::Type{CSV.Sink}) = true

# Constructors
function Sink(sch::Data.Schema, T, append, file::Union{AbstractString, IO}; reference::Vector{UInt8}=UInt8[], kwargs...)
    sink = Sink(file; append=append, colnames=Data.header(sch), kwargs...)
    return sink
end

function Sink(sink, sch::Data.Schema, T, append; reference::Vector{UInt8}=UInt8[])
    sink.append = append
    sink.cols = size(sch, 2)
    !sink.header && !append && writeheaders(sink.io, Data.header(sch), sink.options, sink.quotefields)
    return sink
end

function Data.streamto!(sink::Sink, ::Type{Data.Field}, val, row, col::Int)
    q = Char(sink.options.quotechar); e = Char(sink.options.escapechar); d = Char(sink.options.delim)
    Base.write(sink.io, quoted(sink.quotefields, string(val), q, e, d), ifelse(col == sink.cols, UInt8('\n'), d))
    return nothing
end

function Data.streamto!(sink::Sink, ::Type{Data.Field}, val::Dates.TimeType, row, col::Int)
    v = Dates.format(val, sink.options.dateformat === nothing ? Dates.default_format(typeof(val)) : sink.options.dateformat)
    Base.write(sink.io, v, ifelse(col == sink.cols, UInt8('\n'), sink.options.delim))
    return nothing
end

const EMPTY_UINT8_ARRAY = UInt8[]
function Data.streamto!(sink::Sink, ::Type{Data.Field}, val::Missing, row, col::Int)
    Base.write(sink.io, sink.options.missingcheck ? sink.options.missingstring : EMPTY_UINT8_ARRAY, ifelse(col == sink.cols, UInt8('\n'), sink.options.delim))
    return nothing
end

function Data.close!(sink::CSV.Sink)
    io = isa(sink.fullpath, AbstractString) ? open(sink.fullpath, sink.append ? "a" : "w") : sink.fullpath
    Base.write(io, take!(sink.io))
    isa(sink.fullpath, AbstractString) && close(io)
    return sink
end

function write(sink::Sink, ::Type{T}, args...; append::Bool=false, transforms::Dict=Dict{Int,Function}()) where {T}
    Base.depwarn("CSV.write(file, source) is deprecated in favor of CSV.write(source, file); it will also now support piping like: `dataframe |> CSV.write(filename)`", nothing)
    sink = Data.stream!(T(args...), sink; append=append, transforms=transforms)
    return Data.close!(sink)
end

function write(sink::Sink, source; append::Bool=false, transforms::Dict=Dict{Int,Function}())
    Base.depwarn("CSV.write(file, source) is deprecated in favor of CSV.write(source, file); it will also now support piping like: `dataframe |> CSV.write(filename)`", nothing)
    sink = Data.stream!(source, sink; append=append, transforms=transforms)
    return Data.close!(sink)
end

function validate(fullpath::Union{AbstractString,IO}, sink::T; append::Bool=false, transforms::Dict=Dict{Int,Function}(), kwargs...) where {T}
    Base.depwarn("CSV.validate(filename, $T; kwargs...) is deprecated in favor of CSV.validate(filename; kwargs...)", nothing)
    validate(Source(fullpath; kwargs...))
    return
end


text:
#=
  "Cats and rats and elephants"
  https://www.allendowney.com/blog/2018/12/11/cats-and-rats-and-elephants/
  """
  A few weeks ago I posted 'Lions and Tigers and Bears',
  [https://www.allendowney.com/blog/2018/12/03/lions-and-tigers-and-bears/]
  which poses a Bayesian problem related to the Dirichlet distribution.  If you have
  not read it, you might want to start there.

  Now here’s the follow-up question:

    Suppose there are six species that might be in a zoo: lions and tigers and bears,
    and cats and rats and elephants. Every zoo has a subset of these species, and every
    subset is equally likely.

    One day we visit a zoo and see 3 lions, 2 tigers, and one bear. Assuming that every
    animal in the zoo has an equal chance to be seen, what is the probability that the
    next animal we see is an elephant?
  """

  Note: The approach in this model is the same as in the
  "Lions and Tigers and Bears" problem, just added
    - the three new animals
    - queries which calculates subsets and number of different animals

Distributions of variable animal7
lion       =>    3755  (0.375500)
tiger      =>    3677  (0.367700)
elephant   =>    2205  (0.220500)
bear       =>     323  (0.032300)
cat        =>      21  (0.002100)
rat        =>      19  (0.001900)

Distribution of unique combinations (unsorted)
Distributions of variable (num:0)
[2, 1, 0]	=>	4362 (0.4362)
[2, 0]	=>	3645 (0.3645)
[3, 2, 1, 0]	=>	1540 (0.154)
[1, 2, 0]	=>	324 (0.0324)
[3, 1, 0]	=>	29 (0.0029)
[2, 0, 1]	=>	25 (0.0025)
[1, 0, 2]	=>	18 (0.0018)
[3, 0, 1, 2]	=>	12 (0.0012)
[3, 1, 0, 2]	=>	11 (0.0011)
[4, 0, 2]	=>	5 (0.0005)
[3, 0, 2, 1]	=>	5 (0.0005)
[4, 0, 1]	=>	4 (0.0004)
[1, 0, 3]	=>	3 (0.0003)
[6, 0]	=>	3 (0.0003)
[1, 3, 0, 2]	=>	2 (0.0002)
[4, 1, 0]	=>	2 (0.0002)
[2, 0, 3, 1]	=>	2 (0.0002)
[3, 2, 0, 1]	=>	2 (0.0002)
[1, 3, 0]	=>	2 (0.0002)
[3, 0]	=>	1 (0.0001)
[1]	=>	1 (0.0001)
[2, 0, 1, 3]	=>	1 (0.0001)
[3, 0, 1]	=>	1 (0.0001)

  Cf lions__tiger_and_bears*.jl
     ~/blog/cats_rats_and_elephants.blog
     ~/webppl/cats_rats_and_elephants.wppl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function cats_rats_and_elephants()
    #  The animals.
    lion     = 1
    tiger    = 2
    bear     = 3
    cat      = 4
    rat      = 5
    elephant = 6

    #  Prior
    #  We draw 6 times with the multinomial distribution.
    #  What is the probability of different combinations of the number of each animal?
    alphas = [1/6, 1/6, 1/6, 1/6, 1/6, 1/6]; #  must sum to 1
    #  alphas = [2/12, 2/12, 2/12, 2/12, 2/12, 2/12];

    #  Draw 6 animals
    n = 6
    x ~ Multinomial(n,alphas)

    #  The probabilities to calculate ("aliased" for simplicity)
    probLion      = x[1]
    probTiger     = x[2]
    probBear      = x[3]
    probCat       = x[4]
    probRat       = x[5]
    probElephant  = x[6]

    #  Posterior: What is the probability of lion, tiger, and bear given the observations?
    n2 = n+1
    o = tzeros(n2)
    for i in 1:n2
        o[i] ~ Categorical(simplex([probLion,probTiger,probBear,probCat,probRat,probElephant]))
    end

    #  Observations
    true ~ Dirac(o[1] == lion)
    true ~ Dirac(o[2] == lion)
    true ~ Dirac(o[3] == lion)
    true ~ Dirac(o[4] == tiger)
    true ~ Dirac(o[5] == tiger)
    true ~ Dirac(o[6] == bear)

    animal7 ~ Dirac(o[7])

    # probbear7     ~ Dirac(o[7] == bear)
    # problion7     ~ Dirac(o[7] == lion)
    # probtiger7    ~ Dirac(o[7] == tiger)
    # probcat7      ~ Dirac(o[7] == cat)
    # probrat7      ~ Dirac(o[7] == rat)
    # probelephant7 ~ Dirac(o[7] == elephant)

    #  Which subset of the animals is the most probable?
    # xunique_ = unique(x)
    # xunique = sort(reshape(xunique_,length(xunique_),1),dims=1)
    xunique = unique(x) # unsorted

    #  Number of different animals in the zoo
    numAnimals ~ Dirac(length(xunique))

    return xunique

end

model = cats_rats_and_elephants()
num_chains = 4

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(15), 10_000)
# chns = sample(model, SMC(1000), 10_000)
# chns = sample(model, IS(), 10_000)

#
display(chns)
show_var_dist_pct(chns,:animal7,["lion","tiger","bear","cat","rat","elephant"])

println()
show_var_dist_pct(chns,:probbear7)

# println()
# show_var_dist_pct(chns,:problion7)

# println()
# show_var_dist_pct(chns,:probtiger7)

# println()
# show_var_dist_pct(chns,:probcat7)

# println()
# show_var_dist_pct(chns,:probrat7)

# println()
# show_var_dist_pct(chns,:probelephant7)

println()
show_var_dist_pct(chns,:numAnimals)

println("Distribution of unique combinations (unsorted)")
chains_params = Turing.MCMCChains.get_sections(chns, :parameters)
genq = generated_quantities(model, chains_params)
show_var_dist_pct(genq)

ext: jl

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
@testset "structures" begin
    primal_dual_map = Dualization.PrimalDualMap{Float64}()
    @test Dualization.is_empty(primal_dual_map)
    push!(primal_dual_map.primal_var_dual_con, VI(1) => CI{VI,MOI.EqualTo}(1))
    @test !Dualization.is_empty(primal_dual_map)

    # Constructors
    model = lp1_test()
    dp_f32 = Dualization.DualProblem{Float32}(model)
    @test typeof(dp_f32) == Dualization.DualProblem{Float32,TestModel{Float64}}

    dp_f64 = Dualization.DualProblem(model) # default is Float64
    @test typeof(dp_f64) == Dualization.DualProblem{Float64,TestModel{Float64}}
end


text:
module ProductizeBayes

import SHA
import StanSample

include("versionmodels/VersionedModel.jl")
include("versionmodels/VersionedFile.jl")

export VersionedModel,
        stan_sample,
        read_samples,
        read_summary,
        stan_summary,
        stan_generate_quantities,
        read_generated_quantities,
        diagnose
end # module


text:
####################################################################
function gpmodelVI(x, flux, σ2, idx; kernel=OU, logs=logs, maxiter = 10, seed = 1, show_trace = true)
####################################################################

    #--------------------------------------------------------
    # Auxiliary
    #--------------------------------------------------------

    rg = MersenneTwister(seed)

    JITTER = 1e-6

    Σdiag = Diagonal(σ2)

    numF  = length(unique(idx))

    @assert(all(sort(unique(idx)) .== collect(1:numF)))

    N = length(x)

    @assert(N == length(flux) == length(idx) == length(σ2))

    @printf("Running gpmodelVIexact_grad with %d and %d filters\n", N, numF)

    #--------------------------------------------------------
    # Calculate kernel matrix and keep fixed!
    #--------------------------------------------------------

    dist2 = pairwise(SqEuclidean(), reshape(x, 1, N), dims=2)

    K1 = zeros(N, N)

    calculatekernelmatrix!(K1, dist2, kernel, [0.0; logs])

    K = PDMat((K1 + K1') / 2 + JITTER * I) # This helps us speed calculations

    #--------------------------------------------------------------------
    # Define functions for calling lower bound.
    # Actual calculation takes place in function "complete_lower_bound"
    #--------------------------------------------------------------------

    v²ₐ = v²ᵦ = 1000.0

    # Convenient calls to lower bound

    function lb(μα, Σα, μβ, Σβ, μf, Σf)

        complete_lower_bound(; y=flux, idx=idx, S=Σdiag, μα=μα, Σα=Σα, v²ₐ=v²ₐ, μβ=μβ, Σβ=Σβ, v²ᵦ=v²ᵦ, μf=μf, Σf=Σf, K=K)

    end

    function lb_grad(μα, Σα, μβ, Σβ, μf, Σf)

        complete_lower_bound_grad(; y=flux, idx=idx, S=Σdiag, μα=μα, Σα=Σα, v²ₐ=v²ₐ, μβ=μβ, Σβ=Σβ, v²ᵦ=v²ᵦ, μf=μf, Σf=Σf, K=K)

    end

    # Unpacks vectorised parameters

    function unpack(param)

        @assert(length(param) == numF*4)

        local μα =               param[1+0*numF:1*numF]

        # the parametrisation used here has implications on the calculation
        # of the gradient, see complete_lower_bound.jl and exact_expectations_grad.jl
        local Σα = Diagonal(exp.(param[1+1*numF:2*numF]))

        local μβ =               param[1+2*numF:3*numF]

        # same comment as above applies here
        local Σβ = Diagonal(exp.(param[1+3*numF:4*numF]))

        return μα, Σα, μβ, Σβ

    end

    # Constructs necessary matrices

    function construct_matrices(μα, Σα, μβ, Σβ)

        local Ā  = Diagonal(μα[idx])

        local Σa = Diagonal(Σα.diag[idx])

        local b̄  = μβ[idx]

        return Ā, Σa, b̄

    end

    # Calculate optimal posterior q(f)

    function optimal_qf(Ā, Σa, b̄)

        # local Σf = ((Ā*inv(Σdiag)*Ā + inv(Σdiag)*Σa) + inv(K))\I

        local invDiagPart = (Ā*(Σdiag\Ā) + Σdiag\Σa)\I #
                                                       # This is faster/more stable
        local Σf = invDiagPart*((invDiagPart + K)\K)   #

        local μf = Σf * Ā * (Σdiag \ (flux - b̄))

        return μf, Σf

    end

    function fg!(F,G, param)

      # do common computations here

      local μα, Σα, μβ, Σβ = unpack(param)

      local Ā, Σa, b̄ = construct_matrices(μα, Σα, μβ, Σβ)

      local μf, Σf = optimal_qf(Ā, Σa, b̄)

      if G != nothing
        # code to compute gradient here
        # writing the result to the vector G
        local gradα, gradΣα, gradβ, gradΣβ = lb_grad(μα, Σα, μβ, Σβ, μf, Σf)

        copyto!(G, -1.0 * [gradα; gradΣα; gradβ; gradΣβ])

      end

      if F != nothing
        # value = ... code to compute objective function
        local value = -1.0 * lb(μα, Σα, μβ, Σβ, μf, Σf)
        return value
      end

    end

    #--------------------------------------------------------
    # Initialise GP hyperparameters with random search
    # and few iterations
    #--------------------------------------------------------

    a0 = [std(flux[findall(i.==idx)])  for i in unique(idx)]

    b0 = [mean(flux[findall(i.==idx)]) for i in unique(idx)]

    initsolutions   = [[a0.+5*randn(numF); log.(5*rand()*ones(numF)); b0.+5*randn(numF); log.(5*rand()*ones(numF))] for i=1:30]

    initialfitness  = @showprogress "\t Random restarts " map(initsolutions) do initsol₀

        Optim.optimize(Optim.only_fg!(fg!), initsol₀, LBFGS(), Optim.Options(iterations = 10)).minimum

    end

    solution = initsolutions[argmin(initialfitness)]

    # Use the LBFGS optimiser

    opt = Optim.Options(show_trace = show_trace, iterations = maxiter, show_every = 50, f_tol=0.0, g_tol=1e-8)

    result = Optim.optimize(Optim.only_fg!(fg!), solution, LBFGS(), opt)

    @printf("\t Minimum obj is %f\n\n", result.minimum)

    # retrieve result and return posterior of parameters describing
    # fitted line as object of type Distributions.MvNormal

    μα, Σα, μβ, Σβ = unpack(result.minimizer)

    #--------------------------------------
    # Predictive function
    #--------------------------------------

    μf, Σf = optimal_qf(construct_matrices(μα, Σα, μβ, Σβ)...)

    function predict(xtest)

        local Ntest = length(xtest)

        # dimensions: N × Ntest
        # local k = calculatekernelmatrix(x, xtest, kernel,  [0.0; logs])
        local k = calculatekernelmatrix(reshape(x, 1, N), reshape(xtest, 1, Ntest), kernel,  [0.0; logs])

        @assert(size(k, 1) == N); @assert(size(xtest, 1) == Ntest)

        # Ntest × 1
        local c = calculatekernelmatrix(xtest, xtest, kernel, [0.0; logs])

        # calculate mean
        local μpred = k'*(K\μf)

        # calculate covariance
        local Σpred = c - k'*(K\k) + k'*(K\(Σf/K))*k  + I*1e-4

        local fDensity = MvNormal(μpred, (Σpred.+Σpred')/2)

        local αDensity = MvNormal(μα, Σα)

        local βDensity = MvNormal(μβ, Σβ)

        local samples = map(1:numF) do f

            [rand(αDensity)[f]*rand(fDensity) .+ rand(βDensity)[f] for i=1:5_000]

        end

        samples#mean(samples), sqrt.(diag(cov(samples)))

        # return μpred, sqrt.(diag(Σpred))

    end

    return MvNormal([μα; μβ], Diagonal([Σα.diag; Σβ.diag])), predict, result.minimum

end


text:
using SimpleCaching
using Test

const testing_cache_dir = mktempdir(prefix = "SimpleCaching_test_")
const cached_type = "testing_object"

module A
    using SimpleCaching

    export heavy_computation

    function heavy_computation(cached_type, testing_cache_dir, x, s::Integer...)
        return @scachefast cached_type testing_cache_dir fill(x, s...)
    end
end

@testset "SimpleCaching.jl" begin
    # test @scache
    res1 = @scache cached_type testing_cache_dir fill(0.0, 10, 10, 10)
    res2 = @scache cached_type testing_cache_dir fill(0.0, 10, 10, 10)
    res3 = @scache cached_type testing_cache_dir fill(Float64(1 - 1), 10, 10, 10)
    res4 = @scache cached_type testing_cache_dir fill(res3[1,1,1], 10, 10, 10)

    res = fill(0.0, 10, 10, 10)

    @test res1 == res
    @test res2 == res
    @test res3 == res
    @test res4 == res

    # test @scachejld
    res1 = @scachejld cached_type testing_cache_dir fill(0.0, 20, 20, 20)
    res2 = @scachejld cached_type testing_cache_dir fill(0.0, 20, 20, 20)
    res3 = @scachejld cached_type testing_cache_dir fill(Float64(1 - 1), 20, 20, 20)
    res4 = @scachejld cached_type testing_cache_dir fill(res3[1,1,1], 20, 20, 20)

    res = fill(0.0, 20, 20, 20)

    @test res1 == res
    @test res2 == res
    @test res3 == res
    @test res4 == res

    # test using macro within another module
    using .A
    hc = heavy_computation(cached_type, testing_cache_dir, 0.0, 20, 20, 20)

    @test hc == res

    # test with local variables
    begin
        local n = 10

        @scache cached_type testing_cache_dir vcat(fill(1, n), fill(2, 2n))
    end
end


text:
# Configuration for nearly neutral simuation
function dfe_mixd( x::Int64 )
  dfe_mixed( x::Int64; adv_probability=0.01,
    alpha_disadv=dfe_disadv_alpha, beta_disadv=dfe_disadv_beta,
    alpha_adv=dfe_adv_alpha, beta_adv=dfe_adv_beta )
end

const nn_simtype = 2    # nearly neutral infinite sites
type_str = "M1"
const mu_list_flag = true
@everywhere const N_list     = [25,50,100,200,400,800,1600,3200,6400]     # popsize
const mu_list= [0.025,0.01,0.005,0.0025]   # Mutation rate as a multiple of 1.0/N
const ngens = 50000      # Generations after burn-in
const burn_in= 2.0    # generations of burn_in as a multiple of 1/mu
dfe = dfe_mixd
dfe_str = "dfe_mixd advprob=0.01 alpha_adv=0.2 alpha_disadv=0.5 beta_adv=0.005 beta_disadv=0.5"
dfe_adv_prob=0.01
dfe_adv_alpha=0.5
dfe_disadv_alpha=0.2
dfe_adv_beta=0.005
dfe_disadv_beta=0.5

max_stars_repo_name: um-tech-evolution/nearly-neutral

max_issues_repo_name: um-tech-evolution/nearly-neutral

max_forks_repo_name: um-tech-evolution/nearly-neutral


text:
function prometLinear(d::Number, q::Number, p::Number)::Float64
    if (d <= q)
        0
    elseif (q < d <= p)
        (d - q) / (p - q)
    else
        1
    end
end

function prometVShape(d::Number, q::Any, p::Number)::Float64
    if (d <= 0)
        0
    elseif (0 < d <= p)
        d / p
    else
        1
    end
end

function prometUsual(d::Number, q::Any, p::Any)::Float64
    if (d == 0)
        0
    else
        1
    end
end

function prometUsual(d::Number, q::Any, p::Any)::Float64
    if (d == 0)
        0
    else
        1
    end
end

function prometQuasi(d::Number, q::Number, p::Any)::Float64
    if (d <= q)
        0
    else
        1
    end
end

function prometLevel(d::Number, q::Number, p::Number)::Float64
    if (d <= q)
        0
    elseif (0 < d <= p)
        1 / 2
    else
       1
    end
end

function promethee(decisionMatrix::DataFrame, prefs::Array{Function,1}, weights::Array{Float64,1}, fns::Array{Function,1}, qs::Array, ps::Array)::PrometheeResult

    actionCount, criteriaCount = size(decisionMatrix)

    dValues = zeros(Float64, criteriaCount, actionCount, actionCount)

    for c in 1:criteriaCount, i in 1:actionCount, j in 1:actionCount
        dValues[c,i,j] = fns[c]([-1, 1]) * (decisionMatrix[i,c] - decisionMatrix[j,c])
    end

    pValues = zeros(Float64, criteriaCount, actionCount, actionCount)

    for c in 1:criteriaCount, i in 1:actionCount, j in 1:actionCount
        pValues[c,i,j] = prefs[c](dValues[c,i,j], qs[c], ps[c])
    end

    pValues .*= weights

    flows = sum(pValues, dims=1)[1,:,:]

    negSum = sum(flows, dims=1)[1,:] / (actionCount - 1)
    posSum = sum(flows, dims=2)[:,1] / (actionCount - 1)

    scores = negSum - posSum

    rankings = scores |> sortperm

    bestIndex = rankings |> first

    return PrometheeResult(decisionMatrix, weights, scores, rankings, bestIndex)

end

max_stars_repo_name: ersagunkuruca/JMcDM

max_issues_repo_name: ersagunkuruca/JMcDM

max_forks_repo_name: ersagunkuruca/JMcDM


text:
using DEMC
#include("../src/DEMC.jl")
using Random
using LinearAlgebra
using Distributions
Random.seed!(31953150)
# set up target distribution: Multivariate Normal
ndim = 5 # Number of dimensions
μ = rand(ndim) # mean of each dimension
A = rand(ndim, ndim)
Σ = A'*A .+ diagm(0 => 2*ones(ndim)) # variance covariance matrix
Σ = Σ./maximum(Σ)/100
distr = MvNormal(μ, Σ)

# log objective function
log_obj(mean) = logpdf(distr, mean)

# set up of DEMCz chain
Npar = length(μ)
blockindex = [1:Npar] # parameter blocks: here choose all parameters to be updated simultaenously
Nblocks = length(blockindex)
eps_scale = 1e-5*ones(Npar) # scale of random error around DE update
γ = 2.38 # scale of DE update, 2.38 is the "optimal" number for a normal distribution
N = 5 # number of chains
K = 10 # every K steps add current N draws to Z
Z = randn((10*ndim, ndim)) # initial distribution (completely off to make a difficult test case)

# Number of iterations in Chain
Ngen = 5000
mc, Z = DEMC.demcz_sample(log_obj, Z, N, K, 5000, Nblocks, blockindex, eps_scale, γ; verbose=false)
# run until almost converged, and start with draws from that chain
mc, Z = DEMC.demcz_sample(log_obj, Z[end-10*ndim:end,:], N, K, Ngen, Nblocks, blockindex, eps_scale, γ; prevrun=mc, verbose=false)

# drop first half of chain
Ntot = size(mc.chain,3)
keep = Int(Ntot-round(Ngen/2))+1:Ntot
Ngen_burned = length(keep)
chain_burned = mc.chain[:,:,keep]
logobj_burned = mc.log_obj[:, keep]
chainflat = DEMC.flatten_chain(chain_burned, N, Ngen_burned, Npar)'
bhat = mean(chainflat,dims=1)[:]
println("\n estimates: ", bhat, "\n dist to true: ", bhat - μ)
# covariance of estimates
b, Σb = DEMC.mean_cov_chain(chain_burned, N, Ngen_burned, Npar)

figure_path = "../img/normpdf_serial/"
accept_ratio, Rhat = DEMC.convergence_check(chain_burned, logobj_burned, figure_path; verbose = false)

@test all(Rhat.<1.1)
@test all(accept_ratio.> 0.1)
@test all(accept_ratio.< 0.45)


text:
# ---
# title: 154. Find Minimum in Rotated Sorted Array II
# id: problem154
# author: Indigo
# date: 2021-01-19
# difficulty: Hard
# categories: Array, Binary Search
# link: <https://leetcode.com/problems/find-minimum-in-rotated-sorted-array-ii/description/>
# hidden: true
# ---
#
# Suppose an array sorted in ascending order is rotated at some pivot unknown to
# you beforehand.
#
# (i.e.,  `[0,1,2,4,5,6,7]` might become  `[4,5,6,7,0,1,2]`).
#
# Find the minimum element.
#
# The array may contain duplicates.
#
# **Example 1:**
#
#
#
#     Input: [1,3,5]
#     Output: 1
#
# **Example 2:**
#
#
#
#     Input: [2,2,2,0,1]
#     Output: 0
#
# **Note:**
#
#   * This is a follow up problem to [Find Minimum in Rotated Sorted Array](https://leetcode.comhttps://leetcode.com/problems/find-minimum-in-rotated-sorted-array/description/).
#   * Would allow duplicates affect the run-time complexity? How and why?
#
#
## @lc code=start
using LeetCode
function find_min_154(nums::Vector{Int})
    left, right = 1, length(nums)
    while left < right
        mid = (left + right) ÷ 2
        if nums[mid] > nums[right]
            left = mid + 1
        elseif nums[mid] < nums[right]
            right = mid
        else
            right -= 1
        end
    end
    nums[left]
end
## @lc code=end


text:
const file=pwd()*"/file.txt"

function execute(cmd)
    @show cmd
    @show read(cmd)
    @show read(cmd,String)
    println("run(cmd);")
    run(cmd)
end

function main()
    cmd=pipeline(`cat $file`,`sort`,`uniq`)
    execute(cmd)
    cmd = pipeline(`cat $file`,`xargs echo`)
    execute(cmd)
    n=11
    cmd = pipeline(`seq $n`,`sort -nr`)
    execute(cmd)
    cmd = pipeline(`echo hello`, `sort`)
    execute(cmd)
    run(`echo hello '|' sort`)
    run(`echo hello \| sort`)
end

main()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
# Perform some tests on Step2

using Test
include("../src/Step2.jl")

@testset "Step2.jl" begin

    # Setting
    eH = 1.0
    eL = 0.1
    E = [eL; eH]
    e_size = size(E)[1]
    # Transition matrix: Pr(j|i) = Π(i,j) (columns sum to 1)
    πHH = 0.925 # e_H | e_H
    πHL = 0.5   # e_H | e_L
    Π = [1 - πHL  πHL
        1 - πHH πHH]
    q = 1.3
    a_lb = -2.0
    a_max = 5
    a_size = 100
    A = range(a_lb, a_max, length = a_size)
    A = [A; ];

    # Test 1: check first dimension of policy function is respected
    pol_func = ones(a_size - 10, e_size)
    @test_throws AssertionError step2(A, E, Π, pol_func)

    # Test 2: check second dimension of policy function is respected
    pol_func = ones(a_size, e_size - 1)
    @test_throws AssertionError step2(A, E, Π, pol_func)

end


text:
relu6(x) = min(max(zero(x),x), eltype(x)(6))
using Flux, CuArrays
m = Dense(10,5) |> gpu
x = rand(10) |> gpu

@show m(x)

m = Chain(
  Dense(28^2, 64),
  BatchNorm(64)) |>gpu

x = rand(28^2,3) |> gpu
@show typeof(x)
@show typeof(m.layers[1].W.data)
@show typeof(BatchNorm.σ)
@show m(x)


text:
module CoarseGrainingPrism3D
#####
using LinearAlgebra
using Profile
#####
#####

### GLOBAL CONSTANT
const K =3
const x = K+1
const y = K/2

### FUNDAMENTAL FUNCTION
#chop function
export numchop
function numchop(x::Real)
    if abs(x) <= 1000*eps()
        x = 0
    end
    return x
end
function numchop(x::ComplexF64)
    ansreal = real(x)
    ansim = imag(x)
    if abs(ansreal) <= 1000*eps()
        ansreal = 0
    end
    if abs(ansim) <= 1000*eps()
        ansim = 0
    end
    return x= ansreal + ansim*im
end
function numchop(x::Array)
    for i in 1:length(x)
        if x[i] != 0
            x[i] = numchop(x[i])
        end
    end
    return x
end

#Define delta_{ijkl} -> coupling rules
export delta
function delta(i::Real,j::Real,k::Real)
	sol = 0
	if i <= (j+k) && j<= (i+k) && i+j+k <= K &&2*(i+j+k)%2 == 0
		sol =1
	end
	return sol
end

#Define quantum numbers qn (this is real)
export qn
function qn(n::Real)
    sol = (exp(pi*n*im/(K+2)) - exp(-pi*n*im/(K+2))) / (exp(pi*im/(K+2)) - exp(-pi*im/(K+2)))
    return real(sol)
end

#Define qn factorial
export qnfact
function qnfact(n::Real)
    sol = 1
    for i in 1:n
        sol *= qn(i)
    end
    return sol
end

#Define square root of quantum dimension
export visqrt
function visqrt(i::Real)
    sol = ((-1+0im)^i )*sqrt(qn(2*i+1))
    return sol
end

#Define triangle equality
export trian
function trian(i::Real,j::Real,k::Real)
    sol = 0
    if delta(i,j,k) == 1
        sol = delta(i,j,k)*sqrt(qnfact(i+j-k)*qnfact(i-j+k)*qnfact(-i+j+k)/qnfact(i+j+k+1))
    end
    return sol
end

# Define Racah-Wigner six-j symbol
export RacWig6j
function RacWig6j(i::Real,j::Real,m::Real,k::Real,l::Real,n::Real)
    a = i+j+m; b = i+l+n; c = k+j+n; d = k+l+m;  e = i+j+k+l; f = i+k+m+n; g = j+l+m+n
    sol = 0
    if delta(i,j,m) != 0 && delta(i,l,n) != 0 && delta(k,j,n) != 0 && delta(k,l,m) != 0
        sumz = 0
        for z in max(a,b,c,d):min(e,f,g)
            sumz += (-1)^z *qnfact(z+1)/
                ((qnfact(e-z)*qnfact(f-z)*qnfact(g-z))* (qnfact(z-a)*qnfact(z-b)*qnfact(z-c)*qnfact(z-d)))
        end
        sol = trian(i,j,m)*trian(i,l,n)*trian(k,j,n)*trian(k,l,m)*sumz
    end
    return sol
end

#Define F-symbol
export Fsymb
function Fsymb(i::Real,j::Real,m::Real,k::Real,l::Real,n::Real)
    sol = 0
    if delta(i,j,m) != 0 && delta(i,l,n) != 0 && delta(k,j,n) != 0 && delta(k,l,m) != 0
        sol = (-1+0im)^(i+j+k+l)*sqrt(qn(2*m+1)*qn(2*n+1)) * RacWig6j(i,j,m,k,l,n)
    end
    return sol
end

#Define G-symbol
export Gsymb
function Gsymb(i::Real,j::Real,m::Real,k::Real,l::Real,n::Real)
    sol = 0
    if delta(i,j,m) != 0 && delta(i,l,n) != 0 && delta(k,j,n) != 0 && delta(k,l,m) != 0
        sol = Fsymb(i,j,m,k,l,n) /(visqrt(m)*visqrt(n))
    end
    return sol
end

### DATA
# Define amplitude for a prism function of 12 edges : usual 9 + 3 diag

export Prism
function Prism(β::Float64,α1::Float64,β1::Float64,α2::Float64,δ::Float64,a1::Float64,b::Float64,
                    a2::Float64,d::Float64,c1::Float64,c2::Float64,f::Float64)
    sol = 0
    if delta(c2,d,δ) !=0 && delta(c2,α2,a2) !=0 && delta(β,d,a2) !=0 && delta(β,α2,δ) !=0 &&
        delta(f,d,α1) !=0 && delta(f,c1,b) !=0 && delta(a1,d,b) !=0 && delta(a1,c1,α1) !=0 &&
        delta(α1,β1,δ) !=0 && delta(c2,β1,f) !=0

        dims = visqrt(β)*visqrt(α1)*visqrt(β1)*visqrt(α2)*visqrt(δ)*visqrt(a1)*visqrt(b)*visqrt(a2)*visqrt(d)*visqrt(c1)*visqrt(c2)*visqrt(f)
        sol = dims*Gsymb(c2,d,δ,β,α2,a2)*Gsymb(f,d,α1,a1,c1,b)*Gsymb(c2,d,δ,α1,β1,f)

    if imag(sol) == 0
        sol = real(sol)
    end
    end
    return sol
end

export SuperIndex
function SuperIndex(AllSpin)
	Size = length(AllSpin)
	ans = Int(2*sum(AllSpin[i]*x^(Size-i) for i in 1:Size)+1)
    return ans
end

export InverseSuperIndex
function InverseSuperIndex(y,Size)
    yt = y - 1
    j = zeros(0)
    for i in (Size-1):-1:0
        d,r = divrem(yt,x^i)
        push!(j,d/2)
        yt = r
    end
    return j
end

export PrAOpti
function PrAOpti()
    PrA = Array{Float64}(undef,0)
	s = Array{Int64}(undef,0)

	sol = 0
	dims = 0
    for β in 0.:0.5:y, d in 0.:0.5:y, a2 in 0.:0.5:y
        if delta(β,d,a2) !=0
            for α2 in 0.:0.5:y, δ in 0.:0.5:y
                if  delta(β,α2,δ) !=0
                    for c2 in 0.:0.5:y
                        if delta(c2,d,δ) !=0 && delta(c2,α2,a2) !=0
                            for f in 0.:0.5:y, β1 in 0.:0.5:y
                                if delta(c2,β1,f) !=0
                                    for α1 in 0.:0.5:y
                                        if delta(f,d,α1) !=0 && delta(α1,β1,δ) !=0
                                            for c1 in 0.:0.5:y, a1 in 0.:0.5:y
                                                if delta(a1,c1,α1) !=0
                                                    for b in 0.:0.5:y
                                                        if delta(f,c1,b) !=0 && delta(a1,d,b) !=0
                                                            dims = visqrt(β)*visqrt(α1)*visqrt(β1)*visqrt(α2)*visqrt(δ)*visqrt(a1)*visqrt(b)*visqrt(a2)*visqrt(d)*visqrt(c1)*visqrt(c2)*visqrt(f)
                                                            sol = numchop(dims*Gsymb(c2,d,δ,β,α2,a2)*Gsymb(f,d,α1,a1,c1,b)*Gsymb(c2,d,δ,α1,β1,f))
                                                            if sol !=0
                                                                push!(PrA,sol)
                                                                push!(s,SuperIndex([β,α1,β1,α2,δ,a1,b,a2,d,c1,c2,f]))
                                                            end
                                                        end
                                                    end
                                                end
                                            end
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
	return PrA, s
end

##### coarse-graining algorithm
# create function returning M_{AB}^C, where C fixed spins, A spins of tetra and B spins of pyramid
export BlocksPrism
function BlocksPrism(PrismData,c2::Float64,d::Float64,δ::Float64)
    PrA = PrismData[1]
    sIndex = PrismData[2]

    TempTetra = Array{Array}(undef,0)
    TempPyra = Array{Array}(undef,0)
    TempAmp = Array{Real}(undef,0)
    PosAmp = 1
    for i in sIndex
        (β,α1,β1,α2,δCut,a1,b,a2,dCut,c1,c2Cut,f) = InverseSuperIndex(i,12)
        if δCut == δ && dCut == d && c2Cut == c2
            push!(TempAmp,PrA[PosAmp])
            push!(TempTetra,[β α2 a2])
            push!(TempPyra,[c1 a1 α1 β1 b f])
        end
        PosAmp += 1
    end
    PosTetra = unique(TempTetra)
    PosPyra = unique(TempPyra)
    SizeA = length(PosTetra)
    SizeB = length(PosPyra)
    mat = zeros(SizeA,SizeB)
    for i in 1:length(TempTetra)
        qa = findall(x -> x == TempTetra[i], PosTetra)
        qb = findall(x -> x == TempPyra[i], PosPyra)
        mat[qa[1],qb[1]] = TempAmp[i]
    end

    if SizeA == 0
        mat = 0
    end
    return mat, PosTetra, PosPyra
end

export svdPrism
function svdPrism(PrismData,c2::Float64,d::Float64,δ::Float64)# parameters c2,d,δ
    mat, TetraIndex, PyraIndex = BlocksPrism(PrismData,c2,d,δ)

    U, s, V = svd(mat)
    V = adjoint(V)

    return U, numchop(s), V, TetraIndex, PyraIndex
end

#HYP : only take the firsts TRUNC singular values
export svdPrismTruncated
function svdPrismTruncated(PrismData,c2::Float64,d::Float64,δ::Float64,trunc)
	U, s, V, TetraIndex, PyraIndex = svdPrism(PrismData,c2,d,δ)

	TruncU = U[:,1:trunc]*sqrt.(s[1:trunc])
	TruncV = transpose(sqrt.(s[1:trunc]))*V[1:trunc,:]

	if real(sqrt(visqrt(c2)*visqrt(d)*visqrt(δ))) == 0
		TruncU = im*sqrt(visqrt(c2)*visqrt(d)*visqrt(δ))*TruncU
		TruncV = -im*sqrt(visqrt(c2)*visqrt(d)*visqrt(δ))*TruncV
	else
		TruncU = sqrt(visqrt(c2)*visqrt(d)*visqrt(δ))*TruncU
		TruncV = sqrt(visqrt(c2)*visqrt(d)*visqrt(δ))*TruncV
	end
    if length(TetraIndex) != 0
        (β,α2,a2) = TetraIndex[1]
        if sign(TruncU[1,1]) != sign(AmpTetraInPrism(c2,d,δ,β,α2,a2))
             TruncU = -(1)*TruncU
             TruncV = -(1)*TruncV
        end
    end

	return real(TruncU), real(TruncV), TetraIndex, PyraIndex, s[1:trunc]
end

export FullsvdPrismTruncated
function FullsvdPrismTruncated(PrismData,trunc)
	FullU = Array{Array}(undef,0)
	FullV = Array{Array}(undef,0)
	FullTetraindex = Array{Array}(undef,0)
	FullPyraindex = Array{Array}(undef,0)
	FullCutindex = Array{Array}(undef,0)
	Fulls = Array{Array}(undef,0)
	for c2 in 0:0.5:y, d in 0:0.5:y, δ in 0:0.5:y
		if delta(c2,d,δ) == 1
			temp = svdPrismTruncated(PrismData,c2,d,δ,trunc)
            if temp[1] !=[0]
    			push!(FullU,temp[1])
    			push!(FullV,temp[2])
    			push!(FullTetraindex,temp[3])
    			push!(FullPyraindex,temp[4])
    			push!(FullCutindex,[c2,d,δ])
    			push!(Fulls,temp[5])
            end
		end
	end

    if FullU[1][1] != 1
        FullU = (1/FullU[1][1])*FullU
        FullV = (1/FullV[1][1])*FullV
    end

	return FullU, FullV, FullTetraindex, FullPyraindex, FullCutindex, Fulls
end

export PrismUVTruncated
function PrismUVTruncated(PrismData,trunc)
	FullU, FullV, FullTetraindex, FullPyraindex, FullCutindex,s = FullsvdPrismTruncated(PrismData,trunc)

	AmpUV = Array{Float64}(undef,0)
	PosUVFaceExt = Array{Int64}(undef,0)
	PosUVFaceTrue = Array{Int64}(undef,0)
	PosUVFaceFalse = Array{Int64}(undef,0)
	SizeC = length(FullCutindex)
	for i in 1:SizeC, j in 1:SizeC
		(c2,d,δ) = FullCutindex[i]
		(c2t,dt,δt) = FullCutindex[j]
		#glue along dA/dB
		if d == dt
			TempTetra = FullTetraindex[i]
			TempPyra = FullPyraindex[j]
			for iA in 1:length(TempTetra), jB in 1:length(TempPyra)
				(β,α2,a2) = TempTetra[iA]
				(c1t,a1t,α1t,β1t,bt,ft) = TempPyra[jB]
				#and glue along uA/uB and vA/vB
				if β == bt && a2 == a1t
					push!(AmpUV,FullU[i][iA]*FullV[j][jB]/(visqrt(d)*visqrt(β)*visqrt(a2)))
					push!(PosUVFaceExt,SuperIndex([β1t,ft,β]))
					push!(PosUVFaceTrue,SuperIndex([c2,δt,c2t,δ,d]))
					push!(PosUVFaceFalse,SuperIndex([c2,α1t,c1t,α2,a2]))
				end
			end
		end
	end
	return AmpUV,PosUVFaceExt,PosUVFaceTrue,PosUVFaceFalse, s
end

export PrismEff
function PrismEff(PrismData,trunc)

	AmpUV,PosUVFaceExt,PosUVFaceTrue,PosUVFaceFalse, s = PrismUVTruncated(PrismData,trunc)
	FullAmp = Array{Real}(undef,0)
	BigSIndex = Array{Real}(undef,0)
	for i in 1:length(PosUVFaceTrue)
		# spins of first prism
		(Χ,c1,a2,c2,d) = InverseSuperIndex(PosUVFaceTrue[i],5) #shared spins
		(γ,e,f) = InverseSuperIndex(PosUVFaceExt[i],3)
		(χ,α1,a1,α2,β) = InverseSuperIndex(PosUVFaceFalse[i],5) # was face 2

		pos = findall(x->x == PosUVFaceTrue[i],PosUVFaceFalse)
		for j in pos
			# spins of second prism
			#(bA,j1A,j4B,j3B,uA)
			(γt,et,ft) = InverseSuperIndex(PosUVFaceExt[j],3)
			(Χ,c1t,a2t,c2t,dt) = InverseSuperIndex(PosUVFaceTrue[j],5)
			#super index for 18 spins prism
			TempIndex = SuperIndex([χ,α1,a1,α2,β,c1t,a2t,c2t,dt,γ,e,f,γt,et,ft,c1,c2,a2])
			# FullTot = vcat(FullTot,[TempIndex AmpUV[i]*AmpUV[j]])
			push!(FullAmp,AmpUV[i]*AmpUV[j])
			push!(BigSIndex,TempIndex)
		end
	end

    p = sortperm(BigSIndex)
    BigSIndex = BigSIndex[p]
    FullAmp = FullAmp[p]

    NewSindex = Array{Int64}(undef,0)
    NewAmp = Array{Float64}(undef,0)

    cnt = 0
    for i in 2:length(BigSIndex)
        cnt +=1
        if BigSIndex[i] == BigSIndex[i-1] && i != length(BigSIndex)
            FullAmp[i-cnt] += FullAmp[i]
        else
            push!(NewAmp,FullAmp[i-cnt])
            push!(NewSindex,BigSIndex[i-1])
            cnt = 0
        end
    end
    if BigSIndex[end] == BigSIndex[end-1]
        NewAmp[end] += FullAmp[end]
    else
        push!(NewAmp,FullAmp[end])
        push!(NewSindex,BigSIndex[end])
    end
    return NewAmp, NewSindex
end

####
"embedding map"
####

# 2-2 move for diag
# WARNING: Move done with Fsym, wich need to be either redefine at every step ? or always keep
# initial one ? but then we know something weird might happen -> do test !! IN THE FOLLOWING
# DONE WITH INITIAL Fsymb

export Move2_2
function Move2_2(SetAmpAndIndex,SizeSetSpins,PosMove,PosTrianglesOfMove)
	# SetIndexAndAmp: full set of amplitude and SuperIndex pos of non zero
	# SizeSetSpins: number of edges in polyhedra associated to SetIndexAndAmp
	# PosMove: Position of the edge/spin that we want to flip
	# PosTrianglesOfMove: Positions of the OTHER 4 edges, forming with PosMove the two ini triangles
	# 				that is PosTrianglesOfMove = [j1,j2 , h1,h2] and j1j2Posmove / h1h2Posmove ini triangles
	#				AND j1,h1 / j2,h2 share a vertex
	Amp,sIndex = SetAmpAndIndex

	TempAmp = Array{Real}(undef,0)
	TempsIndex = Array{Real}(undef,0)
	for i in 1:length(sIndex)
		SetSpins = InverseSuperIndex(sIndex[i],SizeSetSpins)
		for qMove in 0:0.5:y
            TempFsymb = Fsymb(SetSpins[PosTrianglesOfMove[1]],SetSpins[PosTrianglesOfMove[2]],SetSpins[PosMove],SetSpins[PosTrianglesOfMove[4]],SetSpins[PosTrianglesOfMove[3]],qMove)
            if numchop(TempFsymb) != 0
                TempSetSpins = copy(SetSpins)
				push!(TempsIndex,SuperIndex(insert!(deleteat!(TempSetSpins,PosMove),PosMove,qMove)))
				push!(TempAmp,TempFsymb*Amp[i])
			end
		end
	end

    p = sortperm(TempsIndex)
    TempsIndex = TempsIndex[p]
    TempAmp = TempAmp[p]

    NewsIndex = Array{Int64}(undef,0)
	NewData = Array{Float64}(undef,0)

	cnt = 0
	for i in 2:length(TempsIndex)
        cnt +=1
		if TempsIndex[i] == TempsIndex[i-1] && i != length(TempsIndex)
			TempAmp[i-cnt] += TempAmp[i]
		else
            if numchop(TempAmp[i-cnt]) != 0
                push!(NewData,TempAmp[i-cnt])
                push!(NewsIndex,TempsIndex[i-1])
            end
		    cnt = 0
		end
	end
    if TempsIndex[end] == TempsIndex[end-1]
        NewData[end] += TempAmp[end]
    else
        push!(NewData,TempAmp[end])
        push!(NewsIndex,TempsIndex[end])
    end

	return NewData,NewsIndex
end

export Move3_1
function Move3_1(SetAmpAndIndex,SizeSetSpins,PosTrianglesOfMove)
    # SetIndexAndAmp: full set of amplitude and SuperIndex pos of non zero
	# SizeSetSpins: number of edges in polyhedra associated to SetIndexAndAmp
	# PosTrianglesOfMove: Position of the 6 edges of the move, of the form
    #                   [j1 j2 j3 h1 h2 h3] where ji are disappearing spins, h1 h2 h3 keeping spins
    #                                triangles are: j1j3h1, j1j2h2, j2j3h3
	Amp,sIndex = SetAmpAndIndex

    PositionDelete = Array{Bool}(undef,0)
    for i in 1:SizeSetSpins
        if i in PosTrianglesOfMove[1:3]
            push!(PositionDelete,true)
        else
            push!(PositionDelete,false)
        end
    end

    TempAmp = Array{Real}(undef,0)
	TempsIndex = Array{Real}(undef,0)
    for i in 1:length(sIndex)
        SetSpins = InverseSuperIndex(sIndex[i],SizeSetSpins)
        TempDimFactor = (visqrt(SetSpins[PosTrianglesOfMove[1]])*visqrt(SetSpins[PosTrianglesOfMove[3]]))/visqrt(SetSpins[PosTrianglesOfMove[4]])
        TempFsymb =Fsymb(SetSpins[PosTrianglesOfMove[6]],SetSpins[PosTrianglesOfMove[5]],SetSpins[PosTrianglesOfMove[4]],SetSpins[PosTrianglesOfMove[1]],SetSpins[PosTrianglesOfMove[3]],SetSpins[PosTrianglesOfMove[2]])
        if numchop(TempFsymb) != 0
            TempSetSpins = copy(SetSpins)
            push!(TempsIndex,SuperIndex(deleteat!(TempSetSpins,PositionDelete)))
            push!(TempAmp,(1/(TempDimFactor*TempFsymb))*Amp[i])
        end
    end

   p = sortperm(TempsIndex)
   TempsIndex = TempsIndex[p]
   TempAmp = TempAmp[p]

   NewsIndex = Array{Int64}(undef,0)
   NewData = Array{Real}(undef,0)
   cntlast = 0
   cnt = 0
   for i in 2:length(TempsIndex)
       cnt+=1
       if TempsIndex[i] == TempsIndex[i-1] && i != length(TempsIndex)
           #TempAmp[i-cnt] += TempAmp[i]
       else
           if numchop(TempAmp[i-cnt]) != 0
               #push!(NewData,(1/cnt)*TempAmp[i-cnt])
               push!(NewData,TempAmp[i-cnt])
               push!(NewsIndex,TempsIndex[i-1])
               if i == length(TempsIndex)
                   cntlast=cnt
               end
           end
           cnt = 0
        end
    end
    if TempsIndex[end] == TempsIndex[end-1]
        #NewData[end] = cntlast*NewData[end]
        #NewData[end] += TempAmp[end]
        #NewData[end] = (1/(cntlast+1))*NewData[end]
    else
        push!(NewData,TempAmp[end])
        push!(NewsIndex,TempsIndex[end])
    end

    if NewData[1] != 1
        NewData = (1/NewData[1])*NewData
    end

	return NewData,NewsIndex
end

#### Embedding of PrismEff to space of Prism again
export EmbPrismEff
function EmbPrismEff(PrismData,trunc)
    SetAmpAndIndex = PrismEff(PrismData,trunc)

    #Move2_2 for big diagonal
    #Move2_2_BigDiag = Move2_2(SetAmpAndIndex,18,18,[10,11,14,15])
    Move2_2FirstTriangle = Move2_2(SetAmpAndIndex,18,16,[2,10,6,13])

    # Move 2_2 first triangle
    #Move2_2FirstTriangle = Move2_2(Move2_2_BigDiag,18,16,[2,10,6,13])
    Move2_2_BigDiag = Move2_2(Move2_2FirstTriangle,18,18,[10,11,14,15])

    # Move3_1 first triangle
    #Move3_1FirstTriangle = Move3_1(Move2_2FirstTriangle,18,[10,13,14,18,16,7])
    Move3_1FirstTriangle = Move3_1(Move2_2_BigDiag,18,[10,13,14,18,16,7])

    # Move 2_2 second triangle: SHOULD OR NOT ?????
    Move2_2SecondTriangle = Move2_2(Move3_1FirstTriangle,15,14,[4,11,8,12])

    # Move3_1 second triangle: SHOULD OR NOT ?????
    Move3_1SecondTriangle = Move3_1(Move2_2SecondTriangle,15,[10,11,12,15,3,14])

    EffPrism = Move3_1SecondTriangle

    return EffPrism
end

##### test function

#### data from initial tetra in prism
export AmpTetraInPrism
function AmpTetraInPrism(j4::Real,j5::Real,j6::Real,j1::Real,j2::Real,j3::Real)
	ans = 0
	if delta(j4,j5,j6)==1 && delta(j4,j2,j3)==1 && delta(j1,j5,j3)==1 && delta(j1,j2,j6)==1
		dims = visqrt(j1)*visqrt(j2)*visqrt(j3)*visqrt(j4)*visqrt(j5)*visqrt(j6)
		ans = dims*Gsymb(j4,j5,j6,j1,j2,j3)
	end
	return ans
end

export FullAmpTetraInPrism
function FullAmpTetraInPrism()
    Amp = Array{Real}(undef,0)
    Index = Array{Real}(undef,0)
    for j1 in 0:0.5:y, j2 in 0:0.5:y, j3 in 0:0.5:y, j4 in 0:0.5:y, j5 in 0:0.5:y, j6 in 0:0.5:y
        TempAmp = numchop(AmpTetraInPrism(j1,j2,j3,j4,j5,j6))
        if TempAmp!=0
            push!(Amp,TempAmp)
            push!(Index,SuperIndex([j1,j2,j3,j4,j5,j6]))
        end
    end
    return Amp,Index
end

export FullAmpPyramidInPrism
function FullAmpPyramidInPrism()
	Amp = Array{Real}(undef,0)
	Index = Array{Real}(undef,0)
	for j1 in 0:0.5:y, j2 in 0:0.5:y, j3 in 0:0.5:y, j4 in 0:0.5:y, j5 in 0:0.5:y, j6 in 0:0.5:y,
		a in 0:0.5:y, b in 0:0.5:y, c in 0:0.5:y
		TempAmp1 = numchop(AmpTetraInPrism(a,b,c,j1,j2,j3))
		TempAmp2 = numchop(AmpTetraInPrism(a,b,c,j4,j5,j6))
		if TempAmp1!=0 && TempAmp2!=0
			push!(Amp,TempAmp1*TempAmp2/(visqrt(a)*visqrt(b)*visqrt(c)))
			push!(Index,SuperIndex([a,b,c,j1,j2,j3,j4,j5,j6]))
		end
	end
	return Amp,Index
end

end # module


text:
include("../Inference.jl")
using Compat.Test

@testset "Inference" begin
    num_bit = 12
    psi = rand_state(num_bit)
    evidense = [1, 2, 3, 4, 5, 6, 7, 8, 9]
    fidelity = run_inference(psi, evidense)
    println("The final fidelity F = ", fidelity)
    @test isapprox(fidelity, 1, atol=1e-2)
end

max_stars_repo_name: Roger-luo/QCBM

max_issues_repo_name: Roger-luo/QCBM

max_forks_repo_name: Roger-luo/QCBM


text:
using Documenter, PairwiseListMatrices

makedocs(
    doctest=true,
    format=Documenter.HTML(
        prettyurls=get(ENV, "CI", nothing) == "true"
    ),
    sitename="PairwiseListMatrices",
    modules=[PairwiseListMatrices],
    pages=[
        "index.md",
        "api.md"
    ]
)

deploydocs(
    repo="github.com/diegozea/PairwiseListMatrices.jl.git",
    target="build",
    deps=nothing,
    make=nothing
)


text:
module SimpleDrawing
using Plots, LinearAlgebra

export newdraw, draw_circle, draw_arc, draw_segment, draw_point, draw_vector
export find_center, non_colinear_check
export finish, draw

"""
`newdraw()` clears the screen and sets up a blank drawing canvas with
no axis, grid, or legend, and sets the aspect ratio to 1.
"""
function newdraw()
    plot(aspectratio=1, legend=false, axis=false, grid=false)
end

function draw()
end

"""
`draw_circle(x,y,r)` (or `draw_circle(z,r)`) draws a circle centered
at `(x,y)` (at `z`) with radius `r`.
"""
function draw_circle(x::Real, y::Real, r::Real; opts...)
    f(t) = r*cos(t) + x
    g(t) = r*sin(t) + y
    plot!(f,g,0,2pi;opts...)
end

draw_circle(z::Complex, r::Real; opts...) =
    draw_circle(real(z), imag(z), r; opts...)

"""
`draw_circle(a,b,c)` draws a circle through the points given by
the three complex arguments.
"""
function draw_circle(a::Complex, b::Complex, c::Complex; opts...)
    z = find_center(a,b,c)
    r = abs(a-z)
    x,y = reim(z)
    draw_circle(x,y,r; opts...)
end

"""
`draw_arc(x,y,r,t1,t2)` draws the shorter arc centered at `(x,y)`
with radius `r` between angles `t1` and `t2`.
"""
function draw_arc(x::Real, y::Real, r::Real, t1::Real, t2::Real; opts...)
    f(t) = r*cos(t) + x
    g(t) = r*sin(t) + y

    plot!(f,g,t1,t2;opts...)
end

"""
`draw_arc(a,b,c)` where the arguments are `Complex`: Draw the
arc from `a` through `b` to `c`.
"""
function draw_arc(a::Complex, b::Complex, c::Complex; opts...)
    m2pi(x::Real) = mod(x,2π)

    z = find_center(a,b,c)
    x,y = reim(z)
    r = abs(b-z)
    ta = m2pi(angle(a-z))
    tb = m2pi(angle(b-z))
    tc = m2pi(angle(c-z))

    if ta > tc
        ta,tc = tc,ta
    end

    if ta<tb<tc
        return draw_arc(x,y,r,ta,tc; opts...)
    end

    return draw_arc(x,y,r,tc-2pi,ta; opts...)
end

"""
`draw_segment(x,y,xx,yy)` draws a line segment from `(x,y)` to
`(xx,yy)`. Also `draw_segment(z,zz)` for `Complex` arguments.
"""
function draw_segment(a::Real,b::Real,c::Real,d::Real; opts...)
    plot!([a,c],[b,d];opts...)
end

function draw_segment(a::Complex, b::Complex; opts...)
    x,y = reim(a)
    xx,yy = reim(b)
    draw_segment(x,y,xx,yy; opts...)
end

"""
`draw_point(x,y)` [or `draw_point(z)`] draws a point at coordinates `(x,y)`
[or complex location `z`].

`draw_point(list)` plots all the points in the one-dimensional `list` of
complex values.
"""
function draw_point(x::Real, y::Real; opts...)
    plot!([x],[y];marker=1,opts...)
end

draw_point(z::Complex; opts...) = draw_point(real(z),imag(z); opts...)

function draw_point(pts::Array{Complex{T},1};opts...) where T
    for p in pts
        draw_point(p;opts...)
    end
end

"""
`draw_vector` is used to draw vectors (line segments with an arrow at one end).
The variations are:
+ `draw_vector(x,y)` draws a vector from the origin to `(x,y)`.
+ `draw_vector(x,y,basex,basey)` draws a vector from `(basex,basey)` to
`(basex+x,basey+y)`.
+ `draw_vector(z)` draws a vector from the origin to the complex value `z`.
+ `draw_vector(z,basez)` draws a vector from the complex location `basez` to
`z+basez`.
"""
draw_vector(x::Real,y::Real; opts...) = draw_segment(0,0,x,y;arrow=:arrow, opts...)

draw_vector(z::Complex;opts...) = draw_segment(0+0im, z; arrow=:arrow, opts...)

draw_vector(x::Real,y::Real,basex::Real,basey::Real; opts...) =
    draw_segment(basex,basey,basex+x,basey+y; arrow=:arrow, opts...)

draw_vector(z::Complex, basez::Complex; opts...) =
    draw_segment(basez, basez+z; arrow=:arrow, opts...)

"""
`finish()` is used to clean up a drawing after various calls to `draw`.
It removes the axes and the grid lines, and sets the aspect ratio to one.
"""
function finish()
    plot!(aspectratio=1, legend=false, axis=false, grid=false)
end

export find_center, non_colinear_check

"""
`find_center(a,b,c)`: Given three points in the complex plane,
find the point `z` that is equidistant from all three. If the three
points are collinear then return `Inf + Inf*im`.
"""
function find_center(a::Complex, b::Complex, c::Complex)::Complex
    if !non_colinear_check(a,b,c)
        return Inf + im*Inf
    end
    A = collect(reim(a))
    B = collect(reim(b))
    C = collect(reim(c))
    AB = 0.5*(A+B)
    BC = 0.5*(B+C)
    M = [(A-B)'; (B-C)']
    rhs = [dot(AB,A-B), dot(BC,B-C)]

    Z = M\rhs
    return Z[1] + im*Z[2]
end

"""
`non_colinear_check(a,b,c)`: test if the complex numbers are distinct
and noncollinear.
"""
function non_colinear_check(a::Complex,b::Complex, c::Complex)::Bool
    if a==b || b==c || a==c
        return false
    end
    z = (b-a)/(c-a)
    return imag(z) != 0
end

include("cubic.jl")
include("curve.jl")

end


text:
baremodule UnsafeAtomics

abstract type Ordering end

function load end
function store! end
function cas! end
function modify! end

function add! end
function sub! end
function xchg! end
function and! end
function nand! end
function or! end
function xor! end
function max! end
function min! end

right(_, x) = x

module Internal

using Base.Sys: WORD_SIZE
using Base: llvmcall

using ..UnsafeAtomics: UnsafeAtomics, Ordering, right

if !@isdefined(⊼)
    ⊼(a, b) = ~(a & b)
end

include("orderings.jl")
include("core.jl")

end  # module Internal

const unordered = Internal.unordered
const monotonic = Internal.monotonic
const acquire = Internal.acquire
const release = Internal.release
const acq_rel = Internal.acq_rel
const seq_cst = Internal.seq_cst

# Julia names
const acquire_release = acq_rel
const sequentially_consistent = seq_cst

end  # baremodule UnsafeAtomics


text:
const STEPSPERKEY = 1
const SHIFTSTEPSPERKEY = 50

function default_mousemovecb(e::Engine, s::WindowState, ev::MouseMoveEvent)
    p = e.phys
    sim = getsim(p.model)
    ui = e.ui

    if ev.isdrag
        # Move free camera

        scaled_dx = ev.dx / s.width
        scaled_dy = ev.dy / s.height

        if s.right
            action = s.shift ? MJCore.mjMOUSE_MOVE_H : MJCore.mjMOUSE_MOVE_V
        elseif s.left
            action = s.shift ? MJCore.mjMOUSE_ROTATE_H : MJCore.mjMOUSE_ROTATE_V
        else
            action = MJCore.mjMOUSE_ZOOM
        end

        if iszero(p.pert[].active)
            mjv_moveCamera(sim.m, action, scaled_dx, scaled_dy, ui.scn, ui.cam)
        else
            mjv_movePerturb(sim.m, sim.d, action, scaled_dx, scaled_dy, ui.scn, p.pert)
        end
    end

    return
end

function default_buttoncb(e::Engine, s::WindowState, ev::ButtonEvent)
    p = e.phys
    pert = p.pert
    sim = getsim(p.model)
    ui = e.ui

    if isright(ev.button) || isleft(ev.button)
        # set/unset pertubation on previously selected body
        if ispress(ev.action) && s.control && pert[].select > 0
            newpert = s.right ? Int(MJCore.mjPERT_TRANSLATE) : Int(MJCore.mjPERT_ROTATE)
            # perturbation onset: reset reference
            iszero(pert[].active) && mjv_initPerturb(sim.m, sim.d, ui.scn, pert)
            pert[].active = newpert
        elseif isrelease(ev.action)
            # stop pertubation
            pert[].active = 0
        end
    end

    if ev.isdoubleclick
        # new select: find geom and 3D click point, get corresponding body

        # stop perturbation on new select
        pert[].active = 0

        # find the selected body
        selpnt = zeros(MVector{3,Float64})
        selgeom, selskin = Ref(Cint(0)), Ref(Cint(0))
        selbody = mjv_select(
            sim.m,
            sim.d,
            ui.vopt,
            s.width / s.height,
            s.x / s.width,
            (s.height - s.y) / s.height,
            ui.scn,
            selpnt,
            selgeom,
            selskin,
        )

        if isleft(ev.button) && nomod(s)
            if selbody >= 0
                # set body selection (including world body)

                # compute localpos
                tmp = selpnt - sim.d.xpos[:, selbody+1] # TODO
                pert[].localpos = reshape(sim.d.xmat[:, selbody+1], 3, 3) * tmp

                # record selection
                pert[].select = selbody
                pert[].skinselect = selskin[]
            else
                # no body selected, unset selection
                pert[].select = 0
                pert[].skinselect = -1
            end
        elseif isright(ev.button)
            # set camera to look at selected body (including world body)
            selbody >= 0 && (ui.cam[].lookat = selpnt)

            if s.control && selbody > 0
                # also set camera to track selected body (excluding world body)
                ui.cam[].type = Int(MJCore.mjCAMERA_TRACKING)
                ui.cam[].trackbodyid = selbody
                ui.cam[].fixedcamid = -1
            end
        end
    end

    return
end

function default_scrollcb(e::Engine, s::WindowState, ev::ScrollEvent)
    m = getsim(e.phys.model).m
    mjv_moveCamera(m, MJCore.mjMOUSE_ZOOM, 0.0, 0.05 * ev.dy, e.ui.scn, e.ui.cam)
    return
end

function handlers(e::Engine)
    return let e = e, ui = e.ui, p = e.phys
        [
            onevent(ButtonEvent) do s, ev
                default_buttoncb(e, s, ev)
            end,

            onevent(MouseMoveEvent) do s, ev
                default_mousemovecb(e, s, ev)
            end,

            onkey(GLFW.KEY_ESCAPE, what = "Quit") do s, ev
                ispress_or_repeat(ev.action) && (ui.shouldexit = true)
            end,

            onkey(GLFW.KEY_F1, what = "Show help message") do s, ev
                ispress_or_repeat(ev.action) && printhelp(e)
            end,

            onkey(GLFW.KEY_F2, what = "Toggle simulation info") do s, ev
                ispress_or_repeat(ev.action) && (ui.showinfo = !ui.showinfo)
            end,

            onkey(GLFW.KEY_F8, what = "Resize window to $RES_XGA (XGA)") do s, ev
                if ispress(ev.action) && e.ffmpeghandle === nothing
                    GLFW.SetWindowSize(s.window, RES_XGA...)
                end
            end,
            onkey(GLFW.KEY_F8, MOD_SHIFT, what = "Resize window to $RES_SXGA (SXGA)") do s, ev
                if ispress(ev.action) && e.ffmpeghandle === nothing
                    GLFW.SetWindowSize(s.window, RES_SXGA...)
                end
            end,

            onkey(GLFW.KEY_F9, what = "Resize window to $RES_HD (HD)") do s, ev
                if ispress(ev.action) && e.ffmpeghandle === nothing
                    GLFW.SetWindowSize(s.window, RES_HD...)
                end
            end,
            onkey(GLFW.KEY_F9, MOD_SHIFT, what = "Resize window to $RES_FHD (FHD)") do s, ev
                if ispress(ev.action) && e.ffmpeghandle === nothing
                    GLFW.SetWindowSize(s.window, RES_FHD...)
                end
            end,

            onkey(GLFW.KEY_F10, what = "Resize window to defaults") do s, ev
                if ispress(ev.action) && e.ffmpeghandle === nothing
                    GLFW.SetWindowSize(s.window, default_windowsize()...)
                end
            end,

            onkey(GLFW.KEY_F11, what = "Toggle fullscreen") do s, ev
                if ispress(ev.action) && e.ffmpeghandle === nothing
                    ismin = iszero(GLFW.GetWindowAttrib(s.window, GLFW.MAXIMIZED))
                    ismin ? GLFW.MaximizeWindow(s.window) : GLFW.RestoreWindow(s.window)
                end
            end,

            onscroll(what = "Zoom camera") do s, ev
                default_scrollcb(e, s, ev)
            end,

            onkey(GLFW.KEY_A, MOD_CONTROL, what = "Align camera scale") do s, ev
                ispress_or_repeat(ev.action) && alignscale!(ui, getsim(p.model))
            end,

            onkey(GLFW.KEY_V, MOD_CONTROL, what = "Toggle video recording") do s, ev
                if ispress(ev.action)
                    e.ffmpeghandle === nothing ? startrecord!(e) : stoprecord!(e)
                end
            end,

            onkey(GLFW.KEY_BACKSPACE, what = "Reset model") do s, ev
                ispress_or_repeat(ev.action) && reset!(p, mode(e))
            end,

            onkey(GLFW.KEY_R, MOD_CONTROL, what = "Toggle reverse") do s, ev
                if ispress_or_repeat(ev.action)
                    ui.reversed = !ui.reversed
                    p.timer.rate *= -1
                    resettime!(p)
                end
            end,

            onkey(GLFW.KEY_SPACE, what = "Pause") do s, ev
                ispress_or_repeat(ev.action) && setpause!(ui, p, !ui.paused)
            end,

            onevent(
                KeyEvent,
                when = describe(GLFW.KEY_RIGHT),
                what = "Step forward when paused (hold SHIFT for $(SHIFTSTEPSPERKEY) steps)"
            ) do s, ev
                if ui.paused && ev.key === GLFW.KEY_RIGHT && ispress_or_repeat(ev.action)
                    steps = s.shift ? 50 : 1
                    for _ = 1:steps
                        forwardstep!(p, mode(e))
                    end
                end
            end,

            onevent(
                KeyEvent,
                when = describe(GLFW.KEY_LEFT),
                what = "Step backwards when paused (hold SHIFT for $(SHIFTSTEPSPERKEY) steps)"
            ) do s, ev
                if ui.paused && ev.key === GLFW.KEY_LEFT && ispress_or_repeat(ev.action)
                    steps = s.shift ? 50 : 1
                    for _ = 1:steps
                        reversestep!(p, mode(e))
                    end
                end
            end,

            onkey(GLFW.KEY_ENTER, what = "Toggle speed mode") do s, ev
                if ispress_or_repeat(ev.action)
                    ui.speedmode = !ui.speedmode
                    setrate!(p.timer, ui.speedmode ? ui.speedfactor : 1)
                end
            end,

            onkey(GLFW.KEY_UP, MOD_SHIFT, what = "Increase sim rate in speedmode") do s, ev
                if ispress_or_repeat(ev.action)
                    ui.speedfactor *= 2
                    setrate!(p.timer, ui.speedfactor)
                end
            end,

            onkey(GLFW.KEY_DOWN, MOD_SHIFT, what = "Decrease sim rate in speedmode") do s, ev
                if ispress_or_repeat(ev.action)
                    ui.speedfactor /= 2
                    setrate!(p.timer, ui.speedfactor)
                end
            end,

            onkey(GLFW.KEY_RIGHT, MOD_CONTROL, what = "Cycle engine mode forward") do s, ev
                if ispress_or_repeat(ev.action)
                    switchmode!(e, inc(e.curmodeidx, 1, length(e.modes)))
                end
            end,

            onkey(GLFW.KEY_LEFT, MOD_CONTROL, what = "Cycle engine mode backwards") do s, ev
                if ispress_or_repeat(ev.action)
                    switchmode!(e, dec(e.curmodeidx, 1, length(e.modes)))
                end
            end,

            onkey(GLFW.KEY_MINUS, what = "Cycle label mode backwards") do s, ev
                if ispress_or_repeat(ev.action)
                    ui.vopt[].label = dec(ui.vopt[].label, 0, Int(MJCore.mjNLABEL) - 1)
                end
            end,

            onkey(GLFW.KEY_EQUAL, what = "Cycle label mode forward") do s, ev
                if ispress_or_repeat(ev.action)
                    ui.vopt[].label = inc(ui.vopt[].label, 0, Int(MJCore.mjNLABEL) - 1)
                end
            end,

            onkey(GLFW.KEY_LEFT_BRACKET, what = "Cycle frame mode backwards") do s, ev
                if ispress_or_repeat(ev.action)
                    ui.vopt[].frame = dec(ui.vopt[].frame, 0, Int(MJCore.mjNFRAME) - 1)
                end
            end,

            onkey(GLFW.KEY_RIGHT_BRACKET, what = "Cycle frame mode forward") do s, ev
                if ispress_or_repeat(ev.action)
                    ui.vopt[].frame = inc(ui.vopt[].frame, 0, Int(MJCore.mjNFRAME) - 1)
                end
            end,

            gen_mjflag_handlers(ui)...
        ]
    end
end

function gen_mjflag_handlers(ui::UIState)
    handlers = EventHandler[]
    let vopt = ui.vopt, scn = ui.scn
        for i=1:Int(MJCore.mjNVISFLAG)
            key = glfw_lookup_key(MJCore.mjVISSTRING[3, i])
            name = MJCore.mjVISSTRING[1, i]
            h = onkey(key, what = "Toggle $name Viz Flag") do s, ev
                ispress_or_repeat(ev.action) && (vopt[].flags = _toggle(vopt[].flags, i))
            end
            push!(handlers, h)
        end

        for i=1:Int(MJCore.mjNRNDFLAG)
            key = glfw_lookup_key(MJCore.mjRNDSTRING[3, i])
            name = MJCore.mjRNDSTRING[1, i]
            h = onkey(key, what = "Toggle $name Render Flag") do s, ev
                ispress_or_repeat(ev.action) && (scn[].flags = _toggle(scn[].flags, i))
            end
            push!(handlers, h)
        end

        n = MJCore.mjNGROUP

        h = onevent(KeyEvent, when = "[1-$n]", what = "Toggle Group Groups 1-$n") do s, ev
            i = Int(ev.key) - Int('0')
            if ispress_or_repeat(ev.action) && iszero(modbits(s)) && checkbounds(Bool, vopt[].geomgroup, i)
                vopt[].geomgroup = _toggle(vopt[].geomgroup, i)
            end
        end
        push!(handlers, h)

        h = onevent(KeyEvent, when = "SHIFT+[1-$n]", what = "Toggle Site Groups 1-$n") do s, ev
            i = Int(ev.key) - Int('0')
            if ispress_or_repeat(ev.action) && isshift(modbits(s)) && checkbounds(Bool, vopt[].sitegroup, i)
                vopt[].sitegroup = _toggle(vopt[].sitegroup, i)
            end
        end
        push!(handlers, h)
    end
    return handlers
end

@inline function _toggle(A::SVector{N,MJCore.mjtByte}, i::Integer) where {N}
    A = MVector(A)
    A[i] = ifelse(A[i] > 0, 0, 1)
    SVector(A)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct LocalNetworkGatewayPropertiesFormat <: SwaggerModel
    localNetworkAddressSpace::Any # spec type: Union{ Nothing, AddressSpace } # spec name: localNetworkAddressSpace
    gatewayIpAddress::Any # spec type: Union{ Nothing, String } # spec name: gatewayIpAddress
    fqdn::Any # spec type: Union{ Nothing, String } # spec name: fqdn
    bgpSettings::Any # spec type: Union{ Nothing, BgpSettings } # spec name: bgpSettings
    resourceGuid::Any # spec type: Union{ Nothing, String } # spec name: resourceGuid
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState

    function LocalNetworkGatewayPropertiesFormat(;localNetworkAddressSpace=nothing, gatewayIpAddress=nothing, fqdn=nothing, bgpSettings=nothing, resourceGuid=nothing, provisioningState=nothing)
        o = new()
        validate_property(LocalNetworkGatewayPropertiesFormat, Symbol("localNetworkAddressSpace"), localNetworkAddressSpace)
        setfield!(o, Symbol("localNetworkAddressSpace"), localNetworkAddressSpace)
        validate_property(LocalNetworkGatewayPropertiesFormat, Symbol("gatewayIpAddress"), gatewayIpAddress)
        setfield!(o, Symbol("gatewayIpAddress"), gatewayIpAddress)
        validate_property(LocalNetworkGatewayPropertiesFormat, Symbol("fqdn"), fqdn)
        setfield!(o, Symbol("fqdn"), fqdn)
        validate_property(LocalNetworkGatewayPropertiesFormat, Symbol("bgpSettings"), bgpSettings)
        setfield!(o, Symbol("bgpSettings"), bgpSettings)
        validate_property(LocalNetworkGatewayPropertiesFormat, Symbol("resourceGuid"), resourceGuid)
        setfield!(o, Symbol("resourceGuid"), resourceGuid)
        validate_property(LocalNetworkGatewayPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type LocalNetworkGatewayPropertiesFormat

const _property_map_LocalNetworkGatewayPropertiesFormat = Dict{Symbol,Symbol}(Symbol("localNetworkAddressSpace")=>Symbol("localNetworkAddressSpace"), Symbol("gatewayIpAddress")=>Symbol("gatewayIpAddress"), Symbol("fqdn")=>Symbol("fqdn"), Symbol("bgpSettings")=>Symbol("bgpSettings"), Symbol("resourceGuid")=>Symbol("resourceGuid"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_LocalNetworkGatewayPropertiesFormat = Dict{Symbol,String}(Symbol("localNetworkAddressSpace")=>"AddressSpace", Symbol("gatewayIpAddress")=>"String", Symbol("fqdn")=>"String", Symbol("bgpSettings")=>"BgpSettings", Symbol("resourceGuid")=>"String", Symbol("provisioningState")=>"ProvisioningState")
Base.propertynames(::Type{ LocalNetworkGatewayPropertiesFormat }) = collect(keys(_property_map_LocalNetworkGatewayPropertiesFormat))
Swagger.property_type(::Type{ LocalNetworkGatewayPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_LocalNetworkGatewayPropertiesFormat[name]))}
Swagger.field_name(::Type{ LocalNetworkGatewayPropertiesFormat }, property_name::Symbol) =  _property_map_LocalNetworkGatewayPropertiesFormat[property_name]

function check_required(o::LocalNetworkGatewayPropertiesFormat)
    true
end

function validate_property(::Type{ LocalNetworkGatewayPropertiesFormat }, name::Symbol, val)
end


text:
# Input data for generation of the cellular complex shown in figure of Data 2.1.2

V = [
    0.0 1.5 3.0 1.0 1.5 2.0 1.0 1.5 2.0 0.0 1.5 3.0;
    0.0 0.0 0.0 1.0 1.0 1.0 2.0 2.0 2.0 3.0 3.0 3.0
];
EV = [
    [1,  2], [2,  3], [ 4,  5], [ 5,  6],
    [7,  8], [8,  9], [10, 11], [11, 12],
    [1, 10], [4,  7], [ 6,  9], [ 3, 12],
    [2,  5], [8, 11]
];
FV = [
    [1, 2, 4, 5, 7, 8, 10, 11],
    [2, 3, 5, 6, 8, 9, 11, 12],
    [4, 5, 6, 7, 8, 9]
];

# The function K returns the Julia SparseArrays matrix providing
#  the characteristic matrix K𝑟 given as input an array CV
#  specifying each 𝑟-cell as array of vertex indices.

function K( CV )
    I = vcat( [ [k for h in CV[k]] for k=1:length(CV) ]...)
    J = vcat( CV...)
    X = Int8[1 for k=1:length(I)]
    return SparseArrays.sparse(I, J, X)
end

# Generation of the sparse binary matrices

copEV = K(EV);
copFV = K(FV);

Matrix(copEV)
Matrix(copFV)

if todisplay
    VV = [[k] for k in 1:size(V,2)];
    GL.VIEW( GL.numbering(1.0)( (V, [VV, EV]),GL.COLORS[1],0.1 ) );
end;

# Boundary matrix generation
# we construct the sparse matrix [𝜕2] from the product of characteristic matrices EV and the transposed FV

Matrix( copEV * copFV' )

copEF = (copEV * copFV') .÷ 2;

Matrix(copEF)

b1 = (copEF * sparse([1 1 1])') .% 2
b2 = (copEF * sparse([1 1 0])') .% 2

f1 = SparseArrays.findnz(b1)[1]
f2 = SparseArrays.findnz(b2)[1]

e1 = EV[f1]
e2 = EV[f2]

if todisplay
    GL.VIEW( GL.numbering(1.0)( (V, [VV, e1]),GL.COLORS[1],0.1 ) );
    GL.VIEW( GL.numbering(1.0)( (V, [VV, e2]),GL.COLORS[1],0.1 ) );
end;