text:
module ArrowJSON

using Mmap
using StructTypes, JSON3, Tables, SentinelArrays, Arrow

# read json files as "table"
# write to arrow stream/file
# read arrow stream/file back

abstract type Type end
Type() = Null("null")
StructTypes.StructType(::Base.Type{Type}) = StructTypes.AbstractType()

children(::Base.Type{T}) where {T} = Field[]

mutable struct Int <: Type
    name::String
    bitWidth::Int64
    isSigned::Base.Bool
end

Int() = Int("", 0, true)
Type(::Base.Type{T}) where {T <: Integer} = Int("int", 8 * sizeof(T), T <: Signed)
StructTypes.StructType(::Base.Type{Int}) = StructTypes.Mutable()
function juliatype(f, x::Int)
    T = x.bitWidth == 8 ? Int8 : x.bitWidth == 16 ? Int16 :
        x.bitWidth == 32 ? Int32 : x.bitWidth == 64 ? Int64 : Int128
    return x.isSigned ? T : unsigned(T)
end

struct FloatingPoint <: Type
    name::String
    precision::String
end

Type(::Base.Type{T}) where {T <: AbstractFloat} = FloatingPoint("floatingpoint", T == Float16 ? "HALF" : T == Float32 ? "SINGLE" : "DOUBLE")
StructTypes.StructType(::Base.Type{FloatingPoint}) = StructTypes.Struct()
juliatype(f, x::FloatingPoint) = x.precision == "HALF" ? Float16 : x.precision == "SINGLE" ? Float32 : Float64

struct FixedSizeBinary <: Type
    name::String
    byteWidth::Int64
end

Type(::Base.Type{NTuple{N, UInt8}}) where {N} = FixedSizeBinary("fixedsizebinary", N)
children(::Base.Type{NTuple{N, UInt8}}) where {N} = Field[]
StructTypes.StructType(::Base.Type{FixedSizeBinary}) = StructTypes.Struct()
juliatype(f, x::FixedSizeBinary) = NTuple{x.byteWidth, UInt8}

struct Decimal <: Type
    name::String
    precision::Int32
    scale::Int32
end

Type(::Base.Type{Arrow.Decimal{P, S}}) where {P, S} = Decimal("decimal", P, S)
StructTypes.StructType(::Base.Type{Decimal}) = StructTypes.Struct()
juliatype(f, x::Decimal) = Arrow.Decimal{x.precision, x.scale}

mutable struct Timestamp <: Type
    name::String
    unit::String
    timezone::Union{Nothing ,String}
end

Timestamp() = Timestamp("", "", nothing)
unit(U) = U == Arrow.Meta.TimeUnit.SECOND ? "SECOND" :
          U == Arrow.Meta.TimeUnit.MILLISECOND ? "MILLISECOND" :
          U == Arrow.Meta.TimeUnit.MICROSECOND ? "MICROSECOND" : "NANOSECOND"
Type(::Base.Type{Arrow.Timestamp{U, TZ}}) where {U, TZ} = Timestamp("timestamp", unit(U), TZ === nothing ? nothing : String(TZ))
StructTypes.StructType(::Base.Type{Timestamp}) = StructTypes.Mutable()
unitT(u) = u == "SECOND" ? Arrow.Meta.TimeUnit.SECOND :
           u == "MILLISECOND" ? Arrow.Meta.TimeUnit.MILLISECOND :
           u == "MICROSECOND" ? Arrow.Meta.TimeUnit.MICROSECOND : Arrow.Meta.TimeUnit.NANOSECOND
juliatype(f, x::Timestamp) = Arrow.Timestamp{unitT(x.unit), x.timezone === nothing ? nothing : Symbol(x.timezone)}

struct Duration <: Type
    name::String
    unit::String
end

Type(::Base.Type{Arrow.Duration{U}}) where {U} = Duration("duration", unit(U))
StructTypes.StructType(::Base.Type{Duration}) = StructTypes.Struct()
juliatype(f, x::Duration) = Arrow.Duration{unit%(x.unit)}

struct Date <: Type
    name::String
    unit::String
end

Type(::Base.Type{Arrow.Date{U, T}}) where {U, T} = Date("date", U == Arrow.Meta.DateUnit.DAY ? "DAY" : "MILLISECOND")
StructTypes.StructType(::Base.Type{Date}) = StructTypes.Struct()
juliatype(f, x::Date) = Arrow.Date{x.unit == "DAY" ? Arrow.Meta.DateUnit.DAY : Arrow.Meta.DateUnit.MILLISECOND, x.unit == "DAY" ? Int32 : Int64}

struct Time <: Type
    name::String
    unit::String
    bitWidth::Int64
end

Type(::Base.Type{Arrow.Time{U, T}}) where {U, T} = Time("time", unit(U), 8 * sizeof(T))
StructTypes.StructType(::Base.Type{Time}) = StructTypes.Struct()
juliatype(f, x::Time) = Arrow.Time{unitT(x.unit), x.unit == "SECOND" || x.unit == "MILLISECOND" ? Int32 : Int64}

struct Interval <: Type
    name::String
    unit::String
end

Type(::Base.Type{Arrow.Interval{U, T}}) where {U, T} = Interval("interval", U == Arrow.Meta.IntervalUnit.YEAR_MONTH ? "YEAR_MONTH" : "DAY_TIME")
StructTypes.StructType(::Base.Type{Interval}) = StructTypes.Struct()
juliatype(f, x::Interval) = Arrow.Interval{x.unit == "YEAR_MONTH" ? Arrow.Meta.IntervalUnit.YEAR_MONTH : Arrow.Meta.IntervalUnit.DAY_TIME, x.unit == "YEAR_MONTH" ? Int32 : Int64}

struct UnionT <: Type
    name::String
    mode::String
    typIds::Vector{Int64}
end

Type(::Base.Type{Arrow.UnionT{T, typeIds, U}}) where {T, typeIds, U} = UnionT("union", T == Arrow.Meta.UnionMode.Dense ? "DENSE" : "SPARSE", collect(typeIds))
children(::Base.Type{Arrow.UnionT{T, typeIds, U}}) where {T, typeIds, U} = Field[Field("", fieldtype(U, i), nothing) for i = 1:fieldcount(U)]
StructTypes.StructType(::Base.Type{UnionT}) = StructTypes.Struct()
juliatype(f, x::UnionT) = Arrow.UnionT{x.mode == "DENSE" ? Arrow.Meta.UnionMode.DENSE : Arrow.Meta.UnionMode.SPARSE, Tuple(x.typeIds), Tuple{(juliatype(y) for y in f.children)...}}

struct List <: Type
    name::String
end

Type(::Base.Type{Vector{T}}) where {T} = List("list")
children(::Base.Type{Vector{T}}) where {T} = [Field("item", T, nothing)]
StructTypes.StructType(::Base.Type{List}) = StructTypes.Struct()
juliatype(f, x::List) = Vector{juliatype(f.children[1])}

struct LargeList <: Type
    name::String
end

StructTypes.StructType(::Base.Type{LargeList}) = StructTypes.Struct()
juliatype(f, x::LargeList) = Vector{juliatype(f.children[1])}

struct FixedSizeList <: Type
    name::String
    listSize::Int64
end

Type(::Base.Type{NTuple{N, T}}) where {N, T} = FixedSizeList("fixedsizelist", N)
children(::Base.Type{NTuple{N, T}}) where {N, T} = [Field("item", T, nothing)]
StructTypes.StructType(::Base.Type{FixedSizeList}) = StructTypes.Struct()
juliatype(f, x::FixedSizeList) = NTuple{x.listSize, juliatype(f.children[1])}

struct Struct <: Type
    name::String
end

Type(::Base.Type{NamedTuple{names, types}}) where {names, types} = Struct("struct")
children(::Base.Type{NamedTuple{names, types}}) where {names, types} = [Field(names[i], fieldtype(types, i), nothing) for i = 1:length(names)]
StructTypes.StructType(::Base.Type{Struct}) = StructTypes.Struct()
juliatype(f, x::Struct) = NamedTuple{Tuple(Symbol(x.name) for x in f.children), Tuple{(juliatype(y) for y in f.children)...}}

struct Map <: Type
    name::String
    keysSorted::Base.Bool
end

Type(::Base.Type{Dict{K, V}}) where {K, V} = Map("map", false)
children(::Base.Type{Dict{K, V}}) where {K, V} = [Field("entries", Arrow.KeyValue{K, V}, nothing)]
StructTypes.StructType(::Base.Type{Map}) = StructTypes.Struct()
juliatype(f, x::Map) = Dict{juliatype(f.children[1].children[1]), juliatype(f.children[1].children[2])}

Type(::Base.Type{Arrow.KeyValue{K, V}}) where {K, V} = Struct("struct")
children(::Base.Type{Arrow.KeyValue{K, V}}) where {K, V} = [Field("key", K, nothing), Field("value", V, nothing)]

struct Null <: Type
    name::String
end

Type(::Base.Type{Missing}) = Null("null")
StructTypes.StructType(::Base.Type{Null}) = StructTypes.Struct()
juliatype(f, x::Null) = Missing

struct Utf8 <: Type
    name::String
end

Type(::Base.Type{<:String}) = Utf8("utf8")
StructTypes.StructType(::Base.Type{Utf8}) = StructTypes.Struct()
juliatype(f, x::Utf8) = String

struct LargeUtf8 <: Type
    name::String
end

StructTypes.StructType(::Base.Type{LargeUtf8}) = StructTypes.Struct()
juliatype(f, x::LargeUtf8) = String

struct Binary <: Type
    name::String
end

Type(::Base.Type{Vector{UInt8}}) = Binary("binary")
children(::Base.Type{Vector{UInt8}}) = Field[]
StructTypes.StructType(::Base.Type{Binary}) = StructTypes.Struct()
juliatype(f, x::Binary) = Vector{UInt8}

struct LargeBinary <: Type
    name::String
end

StructTypes.StructType(::Base.Type{LargeBinary}) = StructTypes.Struct()
juliatype(f, x::LargeBinary) = Vector{UInt8}

struct Bool <: Type
    name::String
end

Type(::Base.Type{Base.Bool}) = Bool("bool")
StructTypes.StructType(::Base.Type{Bool}) = StructTypes.Struct()
juliatype(f, x::Bool) = Base.Bool

StructTypes.subtypekey(::Base.Type{Type}) = :name

const SUBTYPES = @eval (
    int=Int,
    floatingpoint=FloatingPoint,
    fixedsizebinary=FixedSizeBinary,
    decimal=Decimal,
    timestamp=Timestamp,
    duration=Duration,
    date=Date,
    time=Time,
    interval=Interval,
    union=UnionT,
    list=List,
    largelist=LargeList,
    fixedsizelist=FixedSizeList,
    $(Symbol("struct"))=Struct,
    map=Map,
    null=Null,
    utf8=Utf8,
    largeutf8=LargeUtf8,
    binary=Binary,
    largebinary=LargeBinary,
    bool=Bool
)

StructTypes.subtypes(::Base.Type{Type}) = SUBTYPES

const Metadata = Union{Nothing, Vector{NamedTuple{(:key, :value), Tuple{String, String}}}}
Metadata() = nothing

mutable struct DictEncoding
    id::Int64
    indexType::Type
    isOrdered::Base.Bool
end

DictEncoding() = DictEncoding(0, Type(), false)
StructTypes.StructType(::Base.Type{DictEncoding}) = StructTypes.Mutable()

mutable struct Field
    name::String
    nullable::Base.Bool
    type::Type
    children::Vector{Field}
    dictionary::Union{DictEncoding, Nothing}
    metadata::Metadata
end

Field() = Field("", true, Type(), Field[], nothing, Metadata())
StructTypes.StructType(::Base.Type{Field}) = StructTypes.Mutable()
Base.copy(f::Field) = Field(f.name, f.nullable, f.type, f.children, f.dictionary, f.metadata)

function juliatype(f::Field)
    T = juliatype(f, f.type)
    return f.nullable ? Union{T, Missing} : T
end

function Field(nm, ::Base.Type{T}, dictencodings) where {T}
    S = Arrow.maybemissing(T)
    type = Type(S)
    ch = children(S)
    if dictencodings !== nothing && haskey(dictencodings, nm)
        dict = dictencodings[nm]
    else
        dict = nothing
    end
    return Field(nm, T !== S, type, ch, dict, nothing)
end

mutable struct Schema
    fields::Vector{Field}
    metadata::Metadata
end

Schema() = Schema(Field[], Metadata())
StructTypes.StructType(::Base.Type{Schema}) = StructTypes.Mutable()

struct Offsets{T} <: AbstractVector{T}
    data::Vector{T}
end

Base.size(x::Offsets) = size(x.data)
Base.getindex(x::Offsets, i::Base.Int) = getindex(x.data, i)

mutable struct FieldData
    name::String
    count::Int64
    VALIDITY::Union{Nothing, Vector{Int8}}
    OFFSET::Union{Nothing, Offsets}
    TYPE_ID::Union{Nothing, Vector{Int8}}
    DATA::Union{Nothing, Vector{Any}}
    children::Vector{FieldData}
end

FieldData() = FieldData("", 0, nothing, nothing, nothing, nothing, FieldData[])
StructTypes.StructType(::Base.Type{FieldData}) = StructTypes.Mutable()

function FieldData(nm, ::Base.Type{T}, col, dictencodings) where {T}
    if dictencodings !== nothing && haskey(dictencodings, nm)
        refvals = DataAPI.refarray(col.data)
        if refvals !== col.data
            IT = eltype(refvals)
            col = (x - one(T) for x in refvals)
        else
            _, de = dictencodings[nm]
            IT = de.indexType
            vals = unique(col)
            col = Arrow.DictEncoder(col, vals, Arrow.encodingtype(length(vals)))
        end
        return FieldData(nm, IT, col, nothing)
    end
    S = Arrow.maybemissing(T)
    len = Arrow._length(col)
    VALIDITY = OFFSET = TYPE_ID = DATA = nothing
    children = FieldData[]
    if S <: Pair
        return FieldData(nm, Vector{Arrow.KeyValue{Arrow._keytype(S), Arrow._valtype(S)}}, (Arrow.KeyValue(k, v) for (k, v) in pairs(col)))
    elseif S !== Missing
        # VALIDITY
        VALIDITY = Int8[!ismissing(x) for x in col]
        # OFFSET
        if S <: Vector || S == String
            lenfun = S == String ? x->ismissing(x) ? 0 : sizeof(x) : x->ismissing(x) ? 0 : length(x)
            tot = sum(lenfun, col)
            if tot > 2147483647
                OFFSET = String[String(lenfun(x)) for x in col]
                pushfirst!(OFFSET, "0")
            else
                OFFSET = Int32[ismissing(x) ? 0 : lenfun(x) for x in col]
                pushfirst!(OFFSET, 0)
            end
            OFFSET = Offsets(OFFSET)
            push!(children, FieldData("item", eltype(S), Arrow.flatten(skipmissing(col)), dictencodings))
        elseif S <: NTuple
            if Arrow.getT(S) == UInt8
                DATA = [ismissing(x) ? Arrow.ArrowTypes.default(S) : String(collect(x)) for x in col]
            else
                push!(children, FieldData("item", Arrow.getT(S), Arrow.flatten(coalesce(x, Arrow.ArrowTypes.default(S)) for x in col), dictencodings))
            end
        elseif S <: NamedTuple
            for (nm, typ) in zip(Arrow.getnames(S), Arrow.getT(S))
                push!(children, FieldData(String(nm), typ, (getfield(x, nm) for x in col), dictencodings))
            end
        elseif S <: Arrow.UnionT
            U = eltype(S)
            tids = Arrow.typeids(S) === nothing ? (0:fieldcount(U)) : Arrow.typeids(S)
            TYPE_ID = [x === missing ? 0 : tids[Arrow.isatypeid(x, U)] for x in col]
            if Arrow.unionmode(S) == Arrow.Meta.UnionMode.Dense
                offs = zeros(Int32, fieldcount(U))
                OFFSET = Int32[]
                for x in col
                    idx = x === missing ? 1 : Arrow.isatypeid(x, U)
                    push!(OFFSET, offs[idx])
                    offs[idx] += 1
                end
                for i = 1:fieldcount(U)
                    SS = fieldtype(U, i)
                    push!(children, FieldData("$i", SS, Arrow.filtered(i == 1 ? Union{SS, Missing} : Arrow.maybemissing(SS), col), dictencodings))
                end
            else
                for i = 1:fieldcount(U)
                    SS = fieldtype(U, i)
                    push!(children, FieldData("$i", SS, Arrow.replaced(SS, col), dictencodings))
                end
            end
        elseif S <: KeyValue
            push!(children, FieldData("key", Arrow.keyvalueK(S), (x.key for x in col), dictencodings))
            push!(children, FieldData("value", Arrow.keyvalueV(S), (x.value for x in col), dictencodings))
        end
    end
    return FieldData(nm, len, VALIDITY, OFFSET, TYPE_ID, DATA, children)
end

mutable struct RecordBatch
    count::Int64
    columns::Vector{FieldData}
end

RecordBatch() = RecordBatch(0, FieldData[])
StructTypes.StructType(::Base.Type{RecordBatch}) = StructTypes.Mutable()

mutable struct DictionaryBatch
    id::Int64
    data::RecordBatch
end

DictionaryBatch() = DictionaryBatch(0, RecordBatch())
StructTypes.StructType(::Base.Type{DictionaryBatch}) = StructTypes.Mutable()

mutable struct DataFile <: Tables.AbstractColumns
    schema::Schema
    batches::Vector{RecordBatch}
    dictionaries::Vector{DictionaryBatch}
end

Base.propertynames(x::DataFile) = (:schema, :batches, :dictionaries)

function Base.getproperty(df::DataFile, nm::Symbol)
    if nm === :schema
        return getfield(df, :schema)
    elseif nm === :batches
        return getfield(df, :batches)
    elseif nm === :dictionaries
        return getfield(df, :dictionaries)
    end
    return Tables.getcolumn(df, nm)
end

DataFile() = DataFile(Schema(), RecordBatch[], DictionaryBatch[])
StructTypes.StructType(::Base.Type{DataFile}) = StructTypes.Mutable()

parsefile(file) = JSON3.read(Mmap.mmap(file), DataFile)

# make DataFile satisfy Tables.jl interface
function Tables.partitions(x::DataFile)
    if isempty(x.batches)
        # special case empty batches by producing a single DataFile w/ schema
        return (DataFile(x.schema, RecordBatch[], x.dictionaries),)
    else
        return (DataFile(x.schema, [x.batches[i]], x.dictionaries) for i = 1:length(x.batches))
    end
end

Tables.columns(x::DataFile) = x

function Tables.schema(x::DataFile)
    names = map(x -> x.name, x.schema.fields)
    types = map(x -> juliatype(x), x.schema.fields)
    return Tables.Schema(names, types)
end

Tables.columnnames(x::DataFile) =  map(x -> Symbol(x.name), x.schema.fields)

function Tables.getcolumn(x::DataFile, i::Base.Int)
    field = x.schema.fields[i]
    type = juliatype(field)
    return ChainedVector(ArrowArray{type}[ArrowArray{type}(field, length(x.batches) > 0 ? x.batches[j].columns[i] : FieldData(), x.dictionaries) for j = 1:length(x.batches)])
end

function Tables.getcolumn(x::DataFile, nm::Symbol)
    i = findfirst(x -> x.name == String(nm), x.schema.fields)
    return Tables.getcolumn(x, i)
end

struct ArrowArray{T} <: AbstractVector{T}
    field::Field
    fielddata::FieldData
    dictionaries::Vector{DictionaryBatch}
end
ArrowArray(f::Field, fd::FieldData, d) = ArrowArray{juliatype(f)}(f, fd, d)
Base.size(x::ArrowArray) = (x.fielddata.count,)

function Base.getindex(x::ArrowArray{T}, i::Base.Int) where {T}
    @boundscheck checkbounds(x, i)
    S = Base.nonmissingtype(T)
    if x.field.dictionary !== nothing
        fielddata = x.dictionaries[findfirst(y -> y.id == x.field.dictionary.id, x.dictionaries)].data.columns[1]
        field = copy(x.field)
        field.dictionary = nothing
        idx = x.fielddata.DATA[i] + 1
        return ArrowArray(field, fielddata, x.dictionaries)[idx]
    end
    if T === Missing
        return missing
    elseif S <: UnionT
        U = eltype(S)
        tids = Arrow.typeids(S) === nothing ? (0:fieldcount(U)) : Arrow.typeids(S)
        typeid = tids[x.fielddata.TYPE_ID[i]]
        if Arrow.unionmode(S) == Arrow.Meta.UnionMode.DENSE
            off = x.fielddata.OFFSET[i]
            return ArrowArray(x.field.children[typeid+1], x.fielddata.children[typeid+1], x.dictionaries)[off]
        else
            return ArrowArray(x.field.children[typeid+1], x.fielddata.children[typeid+1], x.dictionaries)[i]
        end
    end
    x.fielddata.VALIDITY[i] == 0 && return missing
    if S <: Vector{UInt8}
        return copy(codeunits(x.fielddata.DATA[i]))
    elseif S <: String
        return x.fielddata.DATA[i]
    elseif S <: Vector
        offs = x.fielddata.OFFSET
        A = ArrowArray{eltype(S)}(x.field.children[1], x.fielddata.children[1], x.dictionaries)
        return A[(offs[i] + 1):offs[i + 1]]
    elseif S <: Dict
        offs = x.fielddata.OFFSET
        A = ArrowArray(x.field.children[1], x.fielddata.children[1], x.dictionaries)
        return Dict(y.key => y.value for y in A[(offs[i] + 1):offs[i + 1]])
    elseif S <: Tuple
        if Arrow.getT(S) == UInt8
            A = x.fielddata.DATA
            return Tuple(map(UInt8, collect(A[i][1:x.field.type.byteWidth])))
        else
            sz = x.field.type.listSize
            A = ArrowArray{Arrow.getT(S)}(x.field.children[1], x.fielddata.children[1], x.dictionaries)
            off = (i - 1) * sz + 1
            return Tuple(A[off:(off + sz - 1)])
        end
    elseif S <: NamedTuple
        data = (ArrowArray(x.field.children[j], x.fielddata.children[j], x.dictionaries)[i] for j = 1:length(x.field.children))
        return NamedTuple{Arrow.getnames(S)}(Tuple(data))
    elseif S == Int64 || S == UInt64
        return parse(S, x.fielddata.DATA[i])
    elseif S <: Arrow.Decimal
        str = x.fielddata.DATA[i]
        return S(parse(Int128, str))
    elseif S <: Arrow.Date || S <: Arrow.Time
        val = x.fielddata.DATA[i]
        return Arrow.storagetype(S) == Int32 ? S(val) : S(parse(Int64, val))
    elseif S <: Arrow.Timestamp
        return S(parse(Int64, x.fielddata.DATA[i]))
    else
        return S(x.fielddata.DATA[i])
    end
end

# take any Tables.jl source and write out arrow json datafile
function DataFile(source)
    fields = Field[]
    metadata = nothing # TODO?
    batches = RecordBatch[]
    dictionaries = DictionaryBatch[]
    dictencodings = Dict{String, Tuple{Base.Type, DictEncoding}}()
    dictid = Ref(0)
    for (i, tbl1) in Tables.partitions(source)
        tbl = Arrow.toarrowtable(tbl1)
        if i == 1
            sch = Tables.schema(tbl)
            for (nm, T, col) in zip(sch.names, sch.types, Tables.Columns(tbl))
                if col isa Arrow.DictEncode
                    id = dictid[]
                    dictid[] += 1
                    codes = DataAPI.refarray(col.data)
                    if codes !== col.data
                        IT = Type(eltype(codes))
                    else
                        IT = Type(Arrow.encodingtype(length(unique(col))))
                    end
                    dictencodings[String(nm)] = (T, DictEncoding(id, IT, false))
                end
                push!(fields, Field(String(nm), T, dictencodings))
            end
        end
        # build record batch
        len = Tables.rowcount(tbl)
        columns = FieldData[]
        for (nm, T, col) in zip(sch.names, sch.types, Tables.Columns(tbl))
            push!(columns, FieldData(String(nm), T, col, dictencodings))
        end
        push!(batches, RecordBatch(len, columns))
        # build dictionaries
        for (nm, (T, dictencoding)) in dictencodings
            column = FieldData(nm, T, Tables.getcolumn(tbl, nm), nothing)
            recordbatch = RecordBatch(len, [column])
            push!(dictionaries, DictionaryBatch(dictencoding.id, recordbatch))
        end
    end
    schema = Schema(fields, metadata)
    return DataFile(schema, batches, dictionaries)
end

function Base.isequal(df::DataFile, tbl::Arrow.Table)
    Tables.schema(df) == Tables.schema(tbl) || return false
    i = 1
    for (col1, col2) in zip(Tables.Columns(df), Tables.Columns(tbl))
        if !isequal(col1, col2)
            @show i
            return false
        end
        i += 1
    end
    return true
end

end

max_stars_repo_name: StefanKarpinski/arrow

max_issues_repo_name: StefanKarpinski/arrow

max_forks_repo_name: StefanKarpinski/arrow


text:
module BayesianIntegral

using LinearAlgebra
using Random

include("integrate_given_hyperparameters.jl")
export gaussian_kernel, marginal_gaussian_kernel, K_matrix
export gaussian_kernel_hyperparameters, bayesian_integral_gaussian_exponential
export marginal_likelihood_gaussian_derivatives, log_likelihood
export correlation_vector_of_a_point
include("calibrate_weights.jl")
export calibrate_by_ML_with_SGD, RProp_params, calibrate_by_ML_with_Rprop

end

ext: jl


text:
################################################################################
###  Introduction to Scientific Programming and Machine Learning with Julia  ###
###                                                                          ###
### Run each script on a new clean Julia session                             ###
### GitHub: https://github.com/sylvaticus/IntroSPMLJuliaCourse               ###
### Licence (apply to all material of the course: scripts, videos, quizes,..)###
### Creative Commons By Attribution (CC BY 4.0), Antonello Lobianco          ###
################################################################################

# # 0102 Types and objects

# ## Some stuff to set-up the environment..

cd(@__DIR__)
using Pkg
Pkg.activate(".")
## If using a Julia version different than 1.7 please uncomment and run the following line (reproductibility guarantee will hower be lost)
## Pkg.resolve()
## Pkg.instantiate() # run this if you didn't in Segment 01.01
using Random
Random.seed!(123)

# ## Types

## 1 in not 1.0:
a = 1
b = 1.0
typeof(a) # type is inferred !
typeof(b)
## Convert type (cast)
a = 1
b = convert(Float64,a)
typeof(b)

# Type hierarchy in Julia:
#
# - `Any`
#   - `AbstractString`     (We'll see what all these "abstract" mean....)
#     - `String`
#     - ...
#   - `AbstractArray`
#     - `Array`
#     - ....
#   - `Number`
#     - `Complex`
#     - `Real`
#       - `Rational`
#       - `Integer`
#         - `Unsigned`
#           - `UInt64`
#           - ...
#         - `Signed`
#           - `Int32`
#           - `Int64`
#           - `BigInt`
#           - ...
#         - `Bool`
#       - `FixedPoints`
#         - ...
#       - `AbstractIrrational`
#         - `Irrational`
#       - `AbstractFloat`
#         - `Float32`
#         - `Float64`
#         - `BigFloat`
#         - ...
#   - ...
# Complete Number hierarchy: [https://upload.wikimedia.org/wikipedia/commons/d/d9/Julia-number-type-hierarchy.svg](https://upload.wikimedia.org/wikipedia/commons/d/d9/Julia-number-type-hierarchy.svg)

## Everythong is an object, i.e. of some "type"
c = typeof(a)
typeof(c)
d = sin
typeof(d) <: Function
typeof(+) <: Function

## Operators are just functions:
1 + 2
+(1,2) # this call the function "+"
import Base.+
## +(a,b,c) = a*b*c  # Defining my new crazy addition operation with 3 arguments
10+20+30            # This call it
10+20               # The addition with two parameters remains the same
10+20+30+40         # Also this one remains with the standard addition..

# !!! warning
#     After you tested this singular interpretation of the summation function, please restart Julia or few things will continue to work: with great power comes great responsability, and if we change the meaning of such a inner feature of the language, repercussions will be deep.

# ## Objects and variables

k = 10    # "create" an object Int64 in memory and binds (assign) it to the `k` identifier (the variable name)
typeof(k)
sizeof(k)  # bytes (1 byte is 8 bits)
bitstring(k)
0*2^0+1*2^1+0*2^2+1*2^3
m = k      # name binding: it binds (assign) the entity (object) referenced by a to the b identifier (the variable name)
m == k     # are the two objects equal ?
m === k    # are the two identifiers binding the same identical object in memory ?

# ## Mutability property of Julia objects

k = 10
v = [1,2]
p = 'z'
g = "hello"
ismutable(k)
ismutable(v)
ismutable(p)
ismutable(g)
## mutable objects are stored in memory "directly", while for mutable objects it is its memory address to be stored

# ## Three different ways to "copy" objects...

a = [[[1,2],3],4] # First element of a is said to be "mutable", second one is not:
ismutable(a[1])
ismutable(a[2])
b = a            # binding to a new identifier
c = copy(a)      # create a new copy of a and binds it to c
d = deepcopy(a)  # copy all the references recursively and assign this new object to d

c == a           # are the two objects equal ?
c === a          # are the two identifiers binding the same identical object in memory ?
a[2] = 40        # rebinds a[2] to an other objects and at the same time mutates object a:
b
c
d
a[1][2] = 30     # rebinds a[1][2] and at the same time mutates both a and a[1]
b
c
d
a[1][1][2] = 20
b
c
d
a = 5            # rebinds a:
b
c
d
# !!! note
#     Consider these memory isues when we'll discuss calling a function by reference/value !

max_stars_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_issues_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_forks_repo_name: sylvaticus/IntroSPMLJuliaCourse


text:
module LLVM

using Unicode
using Printf
using Libdl

## source code includes

include("base.jl")
include("version.jl")

# we don't embed the full path to LLVM, because the location might be different at run time.
const libllvm = basename(String(Base.libllvm_path()))
const libllvm_version = Base.libllvm_version

module API
using CEnum
using ..LLVM
using ..LLVM: libllvm

const llvm_version = if version() < v"12"
    "11"
elseif version().major == 12
    "12"
else
    "13"
end
const libdir = joinpath(@__DIR__, "..", "lib")

if !isdir(libdir)
    error("""
    The LLVM API bindings for v$llvm_version do not exist.
    You might need a newer version of LLVM.jl for this version of Julia.""")
end
import LLVMExtra_jll: libLLVMExtra

include(joinpath(libdir, llvm_version, "libLLVM_h.jl"))
include(joinpath(libdir, "libLLVM_extra.jl"))
include(joinpath(libdir, "libLLVM_julia.jl"))
end # module API

# LLVM API wrappers
include("support.jl")
include("types.jl")
include("passregistry.jl")
include("init.jl")
include("core.jl")
include("linker.jl")
include("irbuilder.jl")
include("analysis.jl")
include("moduleprovider.jl")
include("pass.jl")
include("passmanager.jl")
include("execution.jl")
include("buffer.jl")
include("target.jl")
include("targetmachine.jl")
include("datalayout.jl")
include("ir.jl")
include("bitcode.jl")
include("transform.jl")
include("debuginfo.jl")
include("dibuilder.jl")
include("jitevents.jl")
include("utils.jl")

has_orc_v1() = v"8" <= LLVM.version() < v"12"
if has_orc_v1()
    include("orc.jl")
end

has_orc_v2() = v"12" <= LLVM.version()
if has_orc_v2()
    include("orcv2.jl")
end

include("interop.jl")

include("deprecated.jl")

## initialization

function __init__()
    @debug "Using LLVM $libllvm_version at $(Base.libllvm_path())"

    # sanity checks
    if libllvm_version != Base.libllvm_version
        # this checks that the precompilation image isn't being used
        # after having upgraded Julia and the contained LLVM library.
        @error """LLVM.jl was precompiled for LLVM $libllvm_version, whereas you are now using LLVM $(Base.libllvm_version).
                  Please re-compile LLVM.jl."""
    end
    if version() !== runtime_version()
        # this is probably caused by a combination of USE_SYSTEM_LLVM
        # and an LLVM upgrade without recompiling Julia.
        @error """Julia was compiled for LLVM $(version()), whereas you are now using LLVM $(runtime_version()).
                  Please re-compile Julia and LLVM.jl (but note that USE_SYSTEM_LLVM is not a supported configuration)."""
    end

    _install_handlers()
    _install_handlers(GlobalContext())
end

end


text:
using HyperGraphsTools.Bridges.Graphs
import Graphs as _Graphs
using HyperGraphs, SimpleTraits

const AbstractSimpleGraph = _Graphs.SimpleGraphs.AbstractSimpleGraph
const AbstractSimpleEdge = _Graphs.SimpleGraphs.AbstractSimpleEdge

# helpers that return Tuples of edge vertices for easier testing
@traitfn make_edge_tuples(x::T::IsOriented) where {T<:AbstractHyperGraph} = [(src(e)[1], tgt(e)[1]) for e in hyperedges(x)]
@traitfn make_edge_tuples(x::T::(!IsOriented)) where {T<:AbstractHyperGraph} = [Tuple(vx) for vx in vertices(hyperedges(x))]
make_edge_tuples(g::T) where {T<:AbstractSimpleGraph} = [(e.src, e.dst) for e in collect(_Graphs.edges(g))]

# setting up some Graphs graphs
n = rand(1:50)
comp_u = _Graphs.complete_graph(n)
comp_d = _Graphs.complete_digraph(n)
path_u = _Graphs.path_graph(n)
path_d = _Graphs.path_digraph(n)

## Graphs to HyperGraphs ##

# bringing the Graphs types into HyperGraphs
comp_u_x = convert(HyperGraph, comp_u)
comp_d_x = convert(ChemicalHyperGraph, comp_d)
path_u_x = convert(HyperGraph, path_u)
path_d_x = convert(ChemicalHyperGraph, path_d)

@test make_edge_tuples(comp_u_x) == make_edge_tuples(comp_u)
@test make_edge_tuples(comp_d_x) == make_edge_tuples(comp_d)
@test make_edge_tuples(path_u_x) == make_edge_tuples(path_u)
@test make_edge_tuples(path_d_x) == make_edge_tuples(path_d)

# testing incorrect pairs of objects and types (because of non-matching orientation)

# convert a SimpleGraph to a ChemicalHyperGraph (oriented)
@test_throws ErrorException convert(ChemicalHyperGraph, comp_u)
@test_throws ErrorException convert(ChemicalHyperGraph, path_u)

# convert a SimpleDiGraph to a HyperGraph (unoriented)
@test_throws ErrorException convert(HyperGraph, comp_d)
@test_throws ErrorException convert(HyperGraph, path_d)

## HyperGraphs to Graphs ##

# example HyperGraph and its conversion to a Graphs.SimpleGraph
x = HyperGraph([HyperEdge([1, 2]), HyperEdge([2, 3]), HyperEdge([3, 4])])
g = convert(_Graphs.SimpleGraph, x)

@test make_edge_tuples(x) == make_edge_tuples(g)
@test_throws ErrorException convert(_Graphs.SimpleGraph, HyperGraph(HyperEdge([1, 1]))) # loop
@test_throws ErrorException convert(_Graphs.SimpleGraph, HyperGraph([HyperEdge([1, 2]), HyperEdge([1, 2])])) # parallel edge
@test_throws ErrorException convert(_Graphs.SimpleGraph, HyperGraph([HyperEdge([1, 2, 3])])) # cardinality >2
@test_throws MethodError convert(_Graphs.SimpleGraph, HyperGraph{String}()) # type is not Int
@test_throws ErrorException convert(_Graphs.SimpleGraph, HyperGraph([HyperEdge([1, 2]), HyperEdge([4, 5])])) # non-consecutive vertices


text:
using GLEPlot, Test, Colors
const G = GLEPlot

import Base: (//)

(//)(s1::String, s2::String) = strip(s1) == strip(s2)

isapproxstr(s1::AbstractString, s2::AbstractString) =
    isequal(map(s->replace(s, r"\s|\n"=>""), String.((s1, s2)))...)

function m2b(a,b)
    ismissing(a) && ismissing(b) && return true
    isinf(a) && isinf(b) && return true
    isnan(a) && isnan(b) && return true
    return (a == b)
end

function checkzip(z, v::VecOrMat)
    for (zi, vi) ∈ zip(z, eachrow(v))
      all(m2b(zi[j], vi[j]) for j ∈ eachindex(vi)) || begin
            @show zi
            @show vi
            return false
        end
    end
    return true
end


text:
import Base: ndims, size, getindex, reducedim

struct DomainBlocks{N} <: AbstractArray{ArrayDomain{N}, N}
    start::NTuple{N, Int}
    cumlength::Tuple
end
Base.@deprecate_binding BlockedDomains DomainBlocks

ndims(x::DomainBlocks{N}) where {N} = N
size(x::DomainBlocks) = map(length, x.cumlength)
function _getindex(x::DomainBlocks{N}, idx::Tuple) where N
    starts = map((vec, i) -> i == 0 ? 0 : getindex(vec,i), x.cumlength, map(x->x-1, idx))
    ends = map(getindex, x.cumlength, idx)
    ArrayDomain(map(UnitRange, map(+, starts, x.start), map((x,y)->x+y-1, ends, x.start)))
end

function getindex(x::DomainBlocks{N}, idx::Int) where N
    if N == 1
        _getindex(x, (idx,))
    else
        _getindex(x, ind2sub(x, idx))
    end
end

getindex(x::DomainBlocks, idx::Int...) = _getindex(x,idx)

Base.IndexStyle(::Type{<:DomainBlocks}) = IndexCartesian()

function Base.ctranspose(x::DomainBlocks{2})
    DomainBlocks(reverse(x.start), reverse(x.cumlength))
end
function Base.ctranspose(x::DomainBlocks{1})
    DomainBlocks((1, x.start[1]), ([1], x.cumlength[1]))
end

function (*)(x::DomainBlocks{2}, y::DomainBlocks{2})
    if x.cumlength[2] != y.cumlength[1]
        throw(DimensionMismatch("Block distributions being multiplied are not compatible"))
    end
    DomainBlocks((x.start[1],y.start[2]), (x.cumlength[1], y.cumlength[2]))
end

function (*)(x::DomainBlocks{2}, y::DomainBlocks{1})
    if x.cumlength[2] != y.cumlength[1]
        throw(DimensionMismatch("Block distributions being multiplied are not compatible"))
    end
    DomainBlocks((x.start[1],), (x.cumlength[1],))
end

merge_cumsums(x,y) = vcat(x, y+x[end])

function Base.cat(idx::Int, x::DomainBlocks, y::DomainBlocks)
    N = max(ndims(x), ndims(y))
    get_i(x,y, i) = length(x) <= i ? x[i] : length(y) <= i ? y[i] : Int[]
    for i=1:N
        i == idx && continue
        if get_i(x,y,i) != get_i(y,x,i)
            throw(DimensionMismatch("Blocked domains being concatenated have different distributions along dimension $i"))
        end
    end
    output = Any[x.cumlength...]
    output[idx] = merge_cumsums(x.cumlength[idx], y.cumlength[idx])
    DomainBlocks(x.start, (output...))
end

Base.hcat(xs::DomainBlocks...) = cat(2, xs...)
Base.vcat(xs::DomainBlocks...) = cat(1, xs...)

function reducedim(xs::DomainBlocks, dim::Int)
    DomainBlocks(xs.start,
        setindex(xs.cumlength,dim, [1]))
end
function reducedim(dom::DomainBlocks, dim::Tuple)
    reduce(reducedim, dom, dim)
end

cumulative_domains(x::DomainBlocks) = x


text:
using Test, MathOptInterface

const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities

import Cbc
const OPTIMIZER = Cbc.Optimizer()
MOI.set(OPTIMIZER, MOI.Silent(), true)

@testset "SolverName" begin
    @test MOI.get(OPTIMIZER, MOI.SolverName()) == "COIN Branch-and-Cut (Cbc)"
end

@testset "supports_default_copy_to" begin
    @test !MOIU.supports_allocate_load(OPTIMIZER, false)
    @test !MOIU.supports_allocate_load(OPTIMIZER, true)
    @test !MOIU.supports_default_copy_to(OPTIMIZER, false)
    @test !MOIU.supports_default_copy_to(OPTIMIZER, true)
end

const CACHE = MOIU.UniversalFallback(MOIU.Model{Float64}())
const CACHED = MOIU.CachingOptimizer(CACHE, OPTIMIZER)
const BRIDGED = MOI.Bridges.full_bridge_optimizer(CACHED, Float64)

const CONFIG = MOIT.TestConfig(duals = false, infeas_certificates = false)

@testset "basic_constraint_tests" begin
    MOIT.basic_constraint_tests(CACHED, CONFIG)
end

@testset "Unit" begin
    MOIT.unittest(BRIDGED, CONFIG, [
        "number_threads", # FIXME implement `MOI.NumberOfThreads`
        "solve_time", # FIXME implement `MOI.SolveTime`
        "solve_unbounded_model", # INFEASIBLE_OR_UNBOUNDED instead of DUAL_INFEASIBLE
        "delete_soc_variables", "solve_qcp_edge_cases", "solve_qp_edge_cases"  # No quadratics
    ])
end

@testset "ModelLike" begin
    @test MOI.get(CACHED, MOI.SolverName()) == "COIN Branch-and-Cut (Cbc)"
    @testset "default_objective_test" begin
         MOIT.default_objective_test(CACHED)
     end
     @testset "default_status_test" begin
         MOIT.default_status_test(CACHED)
     end
    @testset "nametest" begin
        MOIT.nametest(CACHED)
    end
    @testset "validtest" begin
        MOIT.validtest(CACHED)
    end
    @testset "emptytest" begin
        # Requires VectorOfVariables
        # MOIT.emptytest(CACHED)
    end
    @testset "orderedindicestest" begin
        MOIT.orderedindicestest(CACHED)
    end
    @testset "copytest" begin
        # Requires VectorOfVariables
        # MOIT.copytest(CACHED, MOIU.CachingOptimizer(
        #     ModelForCachingOptimizer{Float64}(),
        #     Cbc.Optimizer()
        # ))
    end
end

@testset "Continuous Linear" begin
    MOIT.contlineartest(BRIDGED, CONFIG,
                        # FIXME starting values
                        ["partial_start"])
end

@testset "Integer Linear" begin
    MOIT.intlineartest(BRIDGED, CONFIG, [
        "indicator1", "indicator2", "indicator3", "indicator4",
        # FIXME starting values
        "knapsack"
    ])
end

@testset "Test params" begin
    knapsack_model = MOIU.Model{Float64}()
    MOIU.loadfromstring!(knapsack_model, """
        variables: x, y
        maxobjective: x + y
        c1: x in ZeroOne()
        c2: y in ZeroOne()
        c3: x + y <= 1.0
    """)
    model = Cbc.Optimizer(maxNodes = 0, presolve = "off", cuts = "off",
                          heur = "off", logLevel = 0)
    MOI.copy_to(model, knapsack_model)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.NODE_LIMIT
    # We also check that options are not destroyed on `empty!`.
    MOI.empty!(model)
    MOI.copy_to(model, knapsack_model)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.NODE_LIMIT
end


text:
export abm_plot, abm_play, abm_video

"""
    abm_plot(model::ABM; kwargs...) → fig, abmstepper
    abm_plot!(ax::Axis/Axis3, model::ABM; kwargs...) → abmstepper
Plot an agent based model by plotting each individual agent as a marker and using
the agent's position field as its location on the plot. Requires `Agents`.

Return the overarching `fig` object, as well as a struct `abmstepper` that can be used
to interactively animate the evolution of the ABM and combine it with other subplots.
The figure is not displayed by default, you need to either return `fig` as a last statement
in your functions or simply call `display(fig)`.
Notice that models with `DiscreteSpace` are plotted starting from 0 to n, with
n the space size along each dimension.

To progress the ABM plot `n` steps simply do:
```julia
Agents.step!(abmstepper, model, agent_step!, model_step!, n)
```
You can still call this function with `n=0` to update the plot for a new `model`,
without doing any stepping. From `fig` you can obtain the plotted axis (to e.g. turn
off ticks, etc.) using `ax = content(fig[1, 1])`.
See [Sugarscape](@ref) for an example of using `abmstepper` to make an animation of
evolving the ABM and a heatmap in parallel with only a few lines of code.

## Agent related keywords
* `ac, as, am`: These three keywords decided the color, size, and marker, that
  each agent will be plotted as. They can each be either a constant or a *function*,
  which takes as an input a single argument and ouputs the corresponding value.
  For example:
  ```julia
  # ac = "#338c54"
  ac(a) = a.status == :S ? "#2b2b33" : a.status == :I ? "#bf2642" : "#338c54"
  # as = 10
  as(a) = 10*randn() + 1
  # am = :diamond
  am(a) = a.status == :S ? :circle : a.status == :I ? :diamond : :rect
  ```
  Notice that for 2D models, `am` can be/return a `Polygon` instance, which plots each agent
  as an arbitrary polygon. It is assumed that the origin (0, 0) is the agent's position when
  creating the polygon. In this case, the keyword `as` is meaningless, as each polygon has its
  own size. Use the functions `scale, rotate2D` to transform this polygon.

  3D models currently do not support having different markers. As a result, `am` cannot be a function.
  It should be a `Mesh` or 3D primitive (such as `Sphere` or `Rect3D`).
* `scheduler = model.scheduler`: decides the plotting order of agents
  (which matters only if there is overlap).
* `offset = nothing`: If not `nothing`, it must be a function taking as an input an
  agent and outputting an offset position tuple to be added to the agent's position
  (which matters only if there is overlap).
* `scatterkwargs = ()`: Additional keyword arguments propagated to the `scatter!` call.

## Preplot related keywords
* `heatarray = nothing` : A keyword that plots a heatmap over the space.
  Its values can be standard data accessors given to functions like `run!`, i.e.
  either a symbol (directly obtain model property) or a function of the model.
  The returned data must be a matrix of the same size as the underlying space.
  For example `heatarray = :temperature` is used in the [Daisyworld](@ref) example.
  But you could also define `f(model) = create_some_matrix_from_model...` and set
  `heatarray = f`. The heatmap will be updated automatically during model evolution
  in videos and interactive applications.
* `heatkwargs = NamedTuple()` : Keywords given to `Makie.heatmap` function
  if `heatarray` is not nothing.
* `static_preplot!` : A function `f(ax, model)` that plots something after the heatmap
  but before the agents. Notice that you can still make objects of this plot be visible
  above the agents using a translation in the third dimension like below:
  ```julia
  function static_preplot!(ax, model)
      obj = scatter!(ax, [50 50]; color = :red) # Show position of teacher
      hidedecorations!(ax) # hide tick labels etc.
      translate!(obj, 0, 0, 5) # be sure that the teacher will be above students
  end
  ```

## Figure related keywords
These only matter for `abm_plot` and not for `abm_plot!`.
* `resolution = (600, 600)`: Resolution of the figure.
* `backgroundcolor = DEFAULT_BG`: Background color of the figure.
* `axiskwargs = NamedTuple()`: Keyword arguments given to the main axis creation for e.g.
  setting `xticksvisible = false`.
* `aspect = DataAspect()`: The aspect ratio behavior of the axis.
"""
function abm_plot(model;
        resolution = (600,600),
        backgroundcolor = DEFAULT_BG, axiskwargs = NamedTuple(),
        aspect = model.space isa Agents.OpenStreetMapSpace ? AxisAspect(1) : DataAspect(),
        kwargs...
    )
    fig = Figure(; resolution, backgroundcolor)
    ax = fig[1,1][1,1] = agents_space_dimensionality(model) == 3 ?
        Axis3(fig; axiskwargs...) : Axis(fig; axiskwargs...)
    ax isa Axis && (ax.aspect = aspect)
    abmstepper = abm_plot!(ax, model; kwargs...)
    return fig, abmstepper
end

function abm_plot!(ax, model;
        ac = JULIADYNAMICS_COLORS[1], as = 10, am = :circle, offset = nothing,
        heatarray = nothing, heatkwargs = NamedTuple(), add_colorbar = true,
        static_preplot! = default_static_preplot, scatterkwargs = NamedTuple(),
        scheduler = model.scheduler,
        kwargs...
    )
    abmstepper = abm_init_stepper(model;
    ac, as, am, scheduler, offset, heatarray)
    abm_init_plot!(ax, model, abmstepper;
        heatkwargs, add_colorbar, static_preplot!, scatterkwargs
    )
    # temporarily disable inspector for poly plots
    if user_used_polygons(am, abmstepper.markers)
        inspector = DataInspector(ax.parent.scene)
        inspector.plot.enabled = false
    end
    return abmstepper
end

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

"""
    abm_play(model, agent_step! [, model_step!]; kwargs...) → fig, abmstepper
Launch an interactive application that plots an agent based model and can animate
its evolution in real time. Requires `Agents`.

The agents are plotted exactly like in [`abm_plot`](@ref), while the two functions
`agent_step!, model_step!` decide how the model will evolve, as in the standard
approach of Agents.jl and its `step!` function.

The application has three buttons:

* "step": advances the simulation once for `spu` steps.
* "run": starts/stops the continuous evolution of the model.
* "reset": resets the model to its original configuration.

Two sliders control the animation speed: "spu" decides how many model steps should be done
before the plot is updated, and "sleep" the `sleep()` time between updates.

## Keywords
* `ac, am, as, scheduler, offset, aspect, scatterkwargs`: propagated to [`abm_plot`](@ref).
* `spu = 1:100`: The values of the "spu" slider.
"""
function abm_play(model, agent_step!, model_step! = Agents.dummystep;
        spu = 1:100, kwargs...)
    fig, abmstepper = abm_plot(model; resolution = (600, 700), kwargs...)
    abm_play!(fig, abmstepper, model, agent_step!, model_step!; spu)
    display(fig)
    return fig, abmstepper
end

function abm_play!(fig, abmstepper, model, agent_step!, model_step!; spu)
    # preinitialize a bunch of stuff
    model0 = deepcopy(model)
    modelobs = Observable(model) # only useful for resetting
    speed, slep, step, run, reset, = abm_controls_play!(fig, model, spu, false)

    # Clicking the step button
    on(step) do clicks
        n = speed[]
        Agents.step!(abmstepper, model, agent_step!, model_step!, n)
    end

    # Clicking the run button
    isrunning = Observable(false)
    on(run) do clicks; isrunning[] = !isrunning[]; end
    on(run) do clicks
        @async while isrunning[]
            n = speed[]
            model = modelobs[] # this is useful only for the reset button
            Agents.step!(abmstepper, model, agent_step!, model_step!, n)
            slep[] == 0 ? yield() : sleep(slep[])
            isopen(fig.scene) || break # crucial, ensures computations stop if closed window.
        end
    end

    # Clicking the reset button
    on(reset) do clicks
        modelobs[] = deepcopy(model0)
        Agents.step!(abmstepper, modelobs[], agent_step!, model_step!, 0)
    end

    return nothing
end

function abm_controls_play!(fig, model, spu, add_update = false)
    controllayout = fig[2, 1] = GridLayout(tellheight = true)
    spusl = labelslider!(fig, "spu =", spu; tellwidth = true)
    if model.space isa Agents.ContinuousSpace
        _s, _v = 0:0.01:1, 0
    else
        _s, _v = 0:0.01:2, 1
    end
    slesl = labelslider!(fig, "sleep =", _s, sliderkw = Dict(:startvalue => _v))
    controllayout[1, :] = spusl.layout
    controllayout[2, :] = slesl.layout
    step = Button(fig, label = "step")
    run = Button(fig, label = "run")
    reset = Button(fig, label = "reset")
    if add_update
        update = Button(fig, label = "update")
        controllayout[3, :] = MakieLayout.hbox!(step, run, reset, update; tellwidth = false)
        upret = update.clicks
    else
        upret = nothing
        controllayout[3, :] = MakieLayout.hbox!(step, run, reset; tellwidth = false)
    end
    return spusl.slider.value, slesl.slider.value, step.clicks, run.clicks, reset.clicks, upret
end

##########################################################################################
"""
    abm_video(file, model, agent_step! [, model_step!]; kwargs...)
This function exports the animated time evolution of an agent based model into a video
saved at given path `file`, by recording the behavior of [`abm_play`](@ref) (without sliders).
The plotting is identical as in [`abm_plot`](@ref) and applicable keywords are propagated.

## Keywords
* `spf = 1`: Steps-per-frame, i.e. how many times to step the model before recording a new
  frame.
* `framerate = 30`: The frame rate of the exported video.
* `frames = 300`: How many frames to record in total, including the starting frame.
* `title = ""`: The title of the figure.
* `showstep = true`: If current step should be shown in title.
* `kwargs...`: All other keywords are propagated to [`abm_plot`](@ref).
"""
function abm_video(file, model, agent_step!, model_step! = Agents.dummystep;
        spf = 1, framerate = 30, frames = 300,  title = "", showstep = true,
        resolution = (600,600), colorscheme = JULIADYNAMICS_COLORS,
        backgroundcolor = DEFAULT_BG, axiskwargs = NamedTuple(), aspect = DataAspect(),
        as = 10, am = :circle, offset = nothing,
        heatarray = nothing, heatkwargs = NamedTuple(), add_colorbar = true,
        static_preplot! = default_static_preplot, scatterkwargs = NamedTuple(),
        kwargs...
    )
    ac = get(kwargs, :ac, colorscheme[1])
    scheduler = get(kwargs, :scheduler, model.scheduler)

    # add some title stuff
    s = Observable(0) # counter of current step
    if title ≠ "" && showstep
        t = lift(x -> title*", step = "*string(x), s)
    elseif showstep
        t = lift(x -> "step = "*string(x), s)
    else
        t = title
    end
    axiskwargs = (title = t, titlealign = :left, axiskwargs...)

    fig, abmstepper = abm_plot(model;
        resolution, colorscheme, backgroundcolor, axiskwargs, aspect,
        ac, as, am, scheduler, offset,
        heatarray, heatkwargs, add_colorbar,
        static_preplot!, scatterkwargs
    )

    record(fig, file; framerate) do io
        for j in 1:frames-1
            recordframe!(io)
            Agents.step!(abmstepper, model, agent_step!, model_step!, spf)
            s[] += spf; s[] = s[]
        end
        recordframe!(io)
    end
    return nothing
end


text:
using Documenter, Main.QuNet

makedocs(
    sitename = "QuNet",
    modules = [QuNet]
)

# deploydocs(
#     deps = Deps.pip("pygments", "mkdocs", "python-markdown-math")
# )

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet


text:
import ClimaCore.Geometry, LinearAlgebra, UnPack
import ClimaCore: Fields, Domains, Topologies, Meshes, DataLayouts, Operators, Geometry, Spaces

using CLIMAParameters
using CLIMAParameters.Planet: ρ_cloud_liq, ρ_cloud_ice, cp_l, cp_i, T_0, LH_f0
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

using RecursiveArrayTools
using OrdinaryDiffEq:
    ODEProblem, solve, SSPRK33, Rosenbrock23, Tsit5, SSPRK432, Feagin14, TsitPap8, CarpenterKennedy2N54
using DifferentialEquations
using UnPack
using LandHydrology
using LandHydrology.SoilHeatParameterizations
using LandHydrology.SoilWaterParameterizations

const FT = Float64

abstract type BC{FT <: AbstractFloat} end

mutable struct FluxBC{FT} <: BC{FT}
    top_heat_flux::FT
    top_water_flux::FT
    btm_heat_flux::FT
    btm_water_flux::FT
end

function compute_soil_rhs!(dY, Y, t, p)

    sp = p[1]
    param_set = p[2]
    zc = p[3]
    @unpack top_heat_flux, btm_heat_flux, top_water_flux, btm_water_flux = p[4]
    @unpack ν, vgn, vgα, vgm, ksat, θr, ρc_ds, κ_sat_unfrozen, κ_sat_frozen = sp
    ϑ_l = Y.x[1]
    θ_i = Y.x[2]
    ρe_int = Y.x[3]

    dϑ_l = dY.x[1]
    dθ_i = dY.x[2]
    dρe_int = dY.x[3]

    # Compute center values of everything
    θ_l = ϑ_l
    ρc_s = volumetric_heat_capacity.(θ_l, θ_i, ρc_ds, Ref(param_set))
    T = temperature_from_ρe_int.(ρe_int, θ_i, ρc_s, Ref(param_set))
    T_sfc = parent(T)[end]
    p[5] = T_sfc
    κ_dry = k_dry(param_set, sp)
    S_r = relative_saturation.(θ_l, θ_i, ν)
    kersten = kersten_number.(θ_i, S_r, Ref(sp))
    κ_sat = saturated_thermal_conductivity.(θ_l, θ_i, κ_sat_unfrozen, κ_sat_frozen)
    κ = thermal_conductivity.(κ_dry, kersten, κ_sat)
    ρe_int_l = volumetric_internal_energy_liq.(T, Ref(param_set))

    cs = axes(θ_i)

    S = effective_saturation.(θ_l; ν = ν, θr = θr)
    K = hydraulic_conductivity.(S; vgm = vgm, ksat = ksat)
    ψ = matric_potential.(S; vgn = vgn, vgα = vgα, vgm = vgm)
    h = ψ .+ zc

    interpc2f = Operators.InterpolateC2F()
    gradc2f_heat = Operators.GradientC2F()
    gradf2c_heat =
        Operators.GradientF2C(top = Operators.SetValue(top_heat_flux), bottom = Operators.SetValue(btm_heat_flux))

    gradc2f_water = Operators.GradientC2F()
    gradf2c_water =
        Operators.GradientF2C(top = Operators.SetValue(top_water_flux), bottom = Operators.SetValue(btm_water_flux))

    @. dϑ_l = -gradf2c_water(-interpc2f(K) * gradc2f_water(h)) #Richards equation
    @. dρe_int = -gradf2c_heat(-interpc2f(κ) * gradc2f_heat(T) - interpc2f(ρe_int_l * K) * gradc2f_water(h))
    dθ_i = Fields.zeros(eltype(θ_i), cs)

    return dY

end

# General composition
ν = FT(0.395);
ν_ss_quartz = FT(0.92)
ν_ss_minerals = FT(0.08)
ν_ss_om = FT(0.0)
ν_ss_gravel = FT(0.0);

#Water specific
Ksat = FT(4.42 / 3600 / 100) # m/s
S_s = FT(1e-3) #inverse meters
vg_n = FT(1.89)
vg_α = FT(7.5); # inverse meters
vg_m = FT(1) - FT(1) / vg_n
θ_r = FT(0)

# Heat specific
κ_quartz = FT(7.7) # W/m/K
κ_minerals = FT(2.5) # W/m/K
κ_om = FT(0.25) # W/m/K
κ_liq = FT(0.57) # W/m/K
κ_ice = FT(2.29); # W/m/K
ρp = FT(2700); # kg/m^3
κ_solid = k_solid(ν_ss_om, ν_ss_quartz, κ_quartz, κ_minerals, κ_om)
κ_sat_frozen = ksat_frozen(κ_solid, ν, κ_ice)
κ_sat_unfrozen = ksat_unfrozen(κ_solid, ν, κ_liq);
ρc_ds = FT((1 - ν) * 1.926e06); # J/m^3/K
a = FT(0.24)
b = FT(18.1)
κ_dry_parameter = FT(0.053)

#collect all params
msp = SoilParams{FT}(
    ν,
    vg_n,
    vg_α,
    vg_m,
    Ksat,
    θ_r,
    S_s,
    ν_ss_gravel,
    ν_ss_om,
    ν_ss_quartz,
    ρc_ds,
    κ_solid,
    ρp,
    κ_sat_unfrozen,
    κ_sat_frozen,
    a,
    b,
    κ_dry_parameter,
)

#Simulation and domain info
t0 = FT(0)
tf = FT(60 * 60 * 72)
dt = FT(0.02)

n = 50

zmax = FT(0)
zmin = FT(-1)
domain = Domains.IntervalDomain(zmin, zmax, x3boundary = (:bottom, :top))
mesh = Meshes.IntervalMesh(domain, nelems = n)

cs = Spaces.CenterFiniteDifferenceSpace(mesh)
fs = Spaces.FaceFiniteDifferenceSpace(cs)
zc = Fields.coordinate_field(cs)

# Boundary conditions
top_water_flux = FT(0)
top_heat_flux = FT(0)
bottom_water_flux = FT(0)
bottom_heat_flux = FT(0)
bc = FluxBC(top_heat_flux, top_water_flux, bottom_heat_flux, bottom_water_flux)

# Parameter structure

# initial conditions
T_max = FT(289.0)
T_min = FT(288.0)
c = FT(20.0)
T = @. T_min + (T_max - T_min) * exp(-(zc - zmax) / (zmin - zmax) * c)
Tsfc = parent(T)[end]
p = [msp, param_set, zc, bc, Tsfc]
θ_i = Fields.zeros(FT, cs)

theta_max = FT(ν * 0.5)
theta_min = FT(ν * 0.4)
θ_l = @. theta_min + (theta_max - theta_min) * exp(-(zc - zmax) / (zmin - zmax) * c)

ρc_s = volumetric_heat_capacity.(θ_l, θ_i, ρc_ds, Ref(param_set))
ρe_int = volumetric_internal_energy.(θ_i, ρc_s, T, Ref(param_set))

Y = ArrayPartition(θ_l, θ_i, ρe_int)

function ∑land_tendencies!(dY, Y, p, t)
    # Intermediate step to be added if needed
    compute_soil_rhs!(dY, Y, t, p)
end

land_prob = ODEProblem(∑land_tendencies!, Y, (t0, tf), p)
algorithm = CarpenterKennedy2N54()

land_simulation() = init(land_prob, algorithm, dt = dt, saveat = 1 * dt) # dt is the land model step

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
# EXCLUDE FROM TESTING
using VectorEngine

function do_sigsegv(addr::Int64)
   @veshow(unsafe_load(convert(Ptr{Int64}, addr)))
   return
end

vesig = VectorEngine.vefunction(do_sigsegv, Tuple{Int64})

vesig(0)

synchronize()


text:
# Printing, show and parse
# ------------------------

Base.summary(seq::BioSequence{<:DNAAlphabet}) = string(length(seq), "nt ", "DNA Sequence")
Base.summary(seq::BioSequence{<:RNAAlphabet}) = string(length(seq), "nt ", "RNA Sequence")
Base.summary(seq::BioSequence{<:AminoAcidAlphabet}) = string(length(seq), "aa ", "Amino Acid Sequence")
Base.summary(seq::BioSequence{<:CharAlphabet}) = string(length(seq), "char ", "Char Sequence")

function Base.print(io::IO, seq::BioSequence; width::Integer = 0)
    col = 1
    for x in seq
        if width > 0 && col > width
            write(io, '\n')
            col = 1
        end
        print(io, x)
        col += 1
    end
end

Base.show(io::IO, seq::BioSequence) = showcompact(io, seq)

function Base.show(io::IO, ::MIME"text/plain", seq::BioSequence)
    println(io, summary(seq), ':')
    showcompact(io, seq)
end

function showcompact(io::IO, seq::BioSequence)
    # don't show more than this many characters
    # to avoid filling the screen with junk
    if isempty(seq)
        print(io, "< EMPTY SEQUENCE >")
    else
        width = displaysize()[2]
        if length(seq) > width
            half = div(width, 2)
            for i in 1:half-1
                print(io, seq[i])
            end
            print(io, '…')
            for i in lastindex(seq)-half+2:lastindex(seq)
                print(io, seq[i])
            end
        else
            for x in seq
                print(io, convert(Char, x))
            end
        end
    end
end

function string_compact(seq::BioSequence)
    buf = IOBuffer()
    showcompact(buf, seq)
    return String(take!(buf))
end

Base.parse(::Type{S}, str::AbstractString) where {S<:BioSequence} = convert(S, str)


text:
# ========================
# iterator with no element
# ========================

struct EmptyIterator{T}
end

function Base.length(::EmptyIterator)
    return 0
end

function Base.eltype(::Type{EmptyIterator{T}}) where {T}
    return T
end

function Base.iterate(::EmptyIterator, state=nothing)
    return nothing
end

# ==============================
# iterator over a single element
# ==============================

struct SingletonIterator{T}
    element::T
end

function Base.length(::SingletonIterator)
    return 1
end

function Base.eltype(::Type{SingletonIterator{T}}) where {T}
    return T
end

function Base.iterate(it::SingletonIterator, state=0)
    if state == 0
        element = it.element
        state = nothing
        return (element, state)
    end
    return nothing
end

# ======================
# iterator over a vector
# ======================

struct VectorIterator{T}
    vector::T
end

function Base.length(it::VectorIterator)
    return length(it.vector)
end

function Base.eltype(::Type{VectorIterator{T}}) where {T}
    return eltype(T)
end

function Base.iterate(it::VectorIterator, state=1)
    if state > length(it)
        return nothing
    end
    element = it.vector[state]
    state += 1
    return (element, state)
end


text:
export read_traceheader

"""
# Info

Use: fileheader = read_fileheader(s::IO; bigendian::Bool = true)

Returns a binary trace header formed from the current position in the stream 's'.
"""
function read_traceheader(s::IO, th_byte2sample::Dict{String,Int32}; bigendian::Bool = true)

    # Initialize binary file header
    traceheader = BinaryTraceHeader()

    if bigendian

        # Read all header values and put in fileheader
        for k in keys(th_byte2sample)

            # Mark start of trace header
            m = mark(s)

            # Seek to this header value location
            skip(s, th_byte2sample[k])

            sym = Symbol(k)
            updateTH_b!(s, traceheader, getfield(traceheader, sym), sym)

            # Return to mark
            reset(s)
        end

    else

        # Read all header values and put in fileheader
        for k in keys(th_byte2sample)

            # Mark start of trace header
            m = mark(s)

            # Seek to this header value location
            skip(s, th_byte2sample[k])

            sym = Symbol(k)
            updateTH!(s, traceheader, getfield(traceheader, sym), sym)

            # Return to mark
            reset(s)
        end

    end

    # Seek to end of trace
    skip(s, 240)

    return traceheader
end

"""
Use: fileheader = read_traceheader(s::IO, keys = Array{String,1}; bigendian::Bool = true)

Returns a binary trace header formed from the current position in the stream 's', only reading
header values denoted in 'keys'.
"""
function read_traceheader(s::IO, keys::Array{String,1}, th_byte2sample::Dict{String, Int32};
                                bigendian::Bool = true)

    # Initialize binary file header
    traceheader = BinaryTraceHeader()

    if bigendian

        # Read all header values and put in fileheader
        for k in keys

            # Mark start of trace header
            m = mark(s)

            # Seek to this header value location
            skip(s, th_byte2sample[k])

            sym = Symbol(k)
            updateTH_b!(s, traceheader, getfield(traceheader, sym), sym)

            # Return to mark
            reset(s)
        end

    else

        # Read all header values and put in fileheader
        for k in keys

            # Mark start of trace header
            m = mark(s)

            # Seek to this header value location
            skip(s, th_byte2sample[k])

            sym = Symbol(k)
            updateTH!(s, traceheader, getfield(traceheader, sym), sym)

            # Return to mark
            reset(s)
        end

    end

    # Seek to end of trace
    skip(s, 240)

    return traceheader
end

function updateTH_b!(s::IO, th::BinaryTraceHeader, field::FT, sym::Symbol) where {FT<:Number}
    val = bswap(read(s, FT))::FT
    setfield!(th, sym, val)
end

function updateTH!(s::IO, th::BinaryTraceHeader, field::FT, sym::Symbol) where {FT<:Number}
    val = read(s, FT)::FT
    setfield!(th, sym, val)
end


text:
struct Image <: AbstractImage
    channels::AbstractArray{Float32, 3}
end

channels(x::Image) = x.channels
Base.size(image::Image) = size(image.channels)
tensorize(image::Image) = Flux.unsqueeze(channels(image), 4)
tensorize(images::Vector{Image}) = reduce((x,y)->cat(x, y; dims=4), map(tensorize, images))

function Base.show(io::IO, x::Image)
    show(io, "text/plain", x.channels)
    println(io, "\n")
end


text:
# TODO:
#   - give variables sensible names
#   - give functions sensible names
#       - follow Julia convention
#   - remove superfluous type annotaitons
#   - Multithreading...
#       - Each root node is independent

# Notes
# canonical labelling -> graph: (where nodes is number of nodes in subgraph)
# a = BitArray(64,[nodes])
# a.chunks = nauty.canonical_labelling(g)[1]
# Array{Int64,2}(a[end-nodes+1:end,:])
include("nauty.jl")

module kavosh
    import LightGraphs
    import nauty
    const lg = LightGraphs
    # IterTools is 10x-2x faster than Combinatorics
    # Make sure it uses revolving door algorithm
    import IterTools
    const it = IterTools

    # Find all subgraphs of size k in G
    # mild breaking change below: verbose now a keyword argument
    function getsubgraphs(G,k;norm=true, verbose=false)::Dict{Array{UInt64,1},Float64}
        # No speedup compared to []
        answers = Dict{Array{UInt64,1},Int64}()
        Visited = zeros(Bool,lg.nv(G))
        # For each node u
        for u in lg.vertices(G)
            if verbose #&& (mod(u,10) == 0)
                print("\r", round(u/lg.nv(G)*100), "% done")
            end
            # 3x speedup compared to Dict(). UInt8 uses ~5% less memory but is ~7% slower.
            S = Dict{Int64,Array{Int64,1}}()
            # "Global" variable Visited
            Visited .= false
            Visited[u] = true
            # S are parents?
            S[1] = [u]
            Enumerate_Vertex(G,u,S,k-1,1,Visited,answers)
        end
        if !norm return answers end
        normalisation = sum(values(answers))
        d = Dict()
        for (key,v) in answers
            d[key] = v/normalisation
        end
        d
    end

    # Take graph, root vertex, "Selection": the vertices that are part of the current motif, the number of nodes left to choose
    function Enumerate_Vertex(G::GraphType,u::Int64,S,Remainder::Int64,i::Int64,Visited::Array{Bool,1},answers)::Void where GraphType <: lg.SimpleGraphs.AbstractSimpleGraph
        # If there are no more nodes to choose, terminate
        s = copy(S) # Stops shorter trees from accidentally sharing data. Must be a neater way of doing this.
        if Remainder == 0
            # Next step: olieshomegrowncannonlabeller(lg.adjacency_matrix(lg.induced_subgraph(G,temp)))
            # This vcat line makes programme ~20% slower
            #temp = vcat(values(s)...)
            #push!(answers,temp)
            temp = vcat(values(s)...)
            # Most of memory usage is in the G[temp] call
            k = nauty.canonical_form(G[temp])[1]
            # Human readable alternative
            #k = nauty.label_to_adj(nauty.canonical_form(G[temp])[1],3)
            answers[k] = get(answers,k,0) + 1
            return
        else
            # Find vertices that could be part of unique motifs
            ValList = Validate(G,s[i],u,Visited)
            # Make sure we don't try to pick more than we can
            n = min(length(ValList),Remainder)
            # Pick k nodes from our current depth
            for k = 1:n
                for combination in it.subsets(ValList,k)
                    # set the current selection at the current depth to them nodes
                    # Might get a performance boost if s was just a list of numbers, and Parents was stored separately.
                    s[i+1] = combination
                    # and Remainder-k from other depths
                    Enumerate_Vertex(G,u,s,Remainder-k,i+1,Visited,answers)
                    # Repeat for all combinations of k nodes at current depth
                end
            end
            # Tidy up
            for v in ValList
                Visited[v] = false
            end
        end
    end

    # Take graph, selected vertices of previous layer, and the root vertex, return vertices that could form unique motifs
    # This is the bit where the labels are considered. Only labels bigger than root are considered. This is to stop double counting.
    function Validate(G::GraphType,Parents::Array{Int64,1},u::Int64,Visited::Array{Bool,1})::Array{Int64,1} where GraphType <: lg.SimpleGraphs.AbstractSimpleGraph
        ValList = Array{Int64,1}()
        # For all of the immediate neighbours of the parents
        for v in Parents
            # Original paper has this as the neighbours of U. 100% sure that it's a typo.
            for w in lg.all_neighbors(G,v)
                # If the label of the neighbour is greater than the parent, and the neighbour has not yet been visited,
                # I think perhaps it should be the root node, actually
                if (u < w) && !Visited[w]
                    # Mark it as visited and add it to our list of candidates
                    Visited[w] = true
                    push!(ValList, w)
                end
            end
        end
        return ValList
    end

    # Questions: what's the difference between ValList and Visited?
end


text:
##
using Printf
using JLD2
using DataDeps
using OceanTurb
using CairoMakie

include("modified_pacalowski_philander_model.jl")
# include("modified_diffusivity_model.jl")

ds = jldopen(joinpath(pwd(), "Data", "three_layer_constant_fluxes_hr192_Qu1.0e-03_Qb0.0e+00_f1.0e-04_Nh256_Nz128__statistics.jld2"))
OUTPUT_NAME = "modified_diffusivity_tstep60_wind_mixing_higher_diffusivity_test.mp4"
SIMULATION_NAME = "Wind Mixing test"

## Load LES grid information, boundary conditions, and initial conditions
Nz = ds["grid/Nz"]
Lz = ds["grid/Lz"]
u₀ = ds["timeseries/u/0"][:]
v₀ = ds["timeseries/v/0"][:]
T₀ = ds["timeseries/T/0"][:]
Fu = ds["boundary_conditions/u_top"]
Fθ = ds["boundary_conditions/θ_top"]
f₀ = ds["parameters/coriolis_parameter"]

## Construct OceanTurb models

# ΔRis = [0.5:0.1:1.0...]
ΔRis = [0.1, 1.]
constants = OceanTurb.Constants(Float64, f=f₀)
parameters = cat(PacanowskiPhilander.Parameters(), [ModifiedPacanowskiPhilanderParameters(ΔRi = ΔRi) for ΔRi in ΔRis]..., dims=1)
models = cat(PacanowskiPhilander.Model(N=Nz, L=Lz, stepper=:BackwardEuler, constants=constants, parameters=parameters[1]),
                [ModifiedPacanowskiPhilanderModel(N=Nz, L=Lz, stepper=:BackwardEuler, constants=constants, parameters=parameters[i+1]) for i in 1:length(ΔRis)]..., dims=1)
model_names = cat("Oceananigans.jl LES", "Pacanowski-Philander", ["Modified Pacanowski-Philander ΔRi = $ΔRi" for ΔRi in ΔRis]..., dims=1)
# pp_parameters = PacanowskiPhilander.Parameters()
# mpp_parameters = ModifiedPacanowskiPhilanderParameters()
# pp_model = PacanowskiPhilander.Model(N=Nz, L=Lz, stepper=:BackwardEuler, constants=constants, parameters=pp_parameters)
# mpp_model = ModifiedPacanowskiPhilanderModel(N=Nz, L=Lz, stepper=:BackwardEuler, constants=constants, parameters=mpp_parameters)

for model in models
    model.bcs[1].top = OceanTurb.FluxBoundaryCondition(Fu)
    model.bcs[3].top = OceanTurb.FluxBoundaryCondition(Fθ)

    model.solution[1].data[1:Nz] .= u₀
    model.solution[2].data[1:Nz] .= v₀
    model.solution[3].data[1:Nz] .= T₀
end

Δt = 60.0
times = [ds["timeseries/t/$i"] for i in keys(ds["timeseries/t"])]
Nt = length(times)

U_solutions = [zeros(Nz, Nt) for i in 1:length(models)]
V_solutions = [zeros(Nz, Nt) for i in 1:length(models)]
T_solutions = [zeros(Nz, Nt) for i in 1:length(models)]

U′W′_solutions = [zeros(Nz+1, Nt) for i in 1:length(models)]
V′W′_solutions = [zeros(Nz+1, Nt) for i in 1:length(models)]
W′T′_solutions = [zeros(Nz+1, Nt) for i in 1:length(models)]

Ri_solutions = [zeros(Nz+1, Nt) for i in 1:length(models)]

function get_diffusive_flux(field_index, model)
    flux = FaceField(model.grid)
    field = model.solution[field_index]
    K = model.timestepper.eqn.K[field_index]
    for i in interiorindices(flux)
        @inbounds flux[i] = - K(model, i) * ∂z(field, i)
    end
    return flux
end

function get_richardson_number_profile(model)
    Ri = FaceField(model.grid)
    for i in interiorindices(Ri)
        @inbounds Ri[i] = local_richardson(model, i)
    end
    return Ri
end

for i in 1:length(models)
    @info "Model $i/$(length(models))"
    model = models[i]
    U_solution = U_solutions[i]
    V_solution = V_solutions[i]
    T_solution = T_solutions[i]
    U′W′_solution = U′W′_solutions[i]
    V′W′_solution = V′W′_solutions[i]
    W′T′_solution = W′T′_solutions[i]
    Ri_solution = Ri_solutions[i]

    for n in 1:Nt
        OceanTurb.run_until!(model, Δt, times[n])
        @info "Time = $(times[n])"

        U_solution[:, n] .= model.solution[1][1:Nz]
        V_solution[:, n] .= model.solution[2][1:Nz]
        T_solution[:, n] .= model.solution[3][1:Nz]

        U′W′_solution[:, n] .= get_diffusive_flux(1, model)[1:Nz+1]
        V′W′_solution[:, n] .= get_diffusive_flux(2, model)[1:Nz+1]
        W′T′_solution[:, n] .= get_diffusive_flux(3, model)[1:Nz+1]

        U′W′_solution[Nz+1, n] = Fu

        Ri_solution[:, n] = get_richardson_number_profile(model)[1:Nz+1]
    end
end

# Get rid of ∞ and super large values.
Ri_solutions = [clamp.(Ri_solution, -1, 2) for Ri_solution in Ri_solutions]

## Plot!
##
U_LES_solution = zeros(Nz, Nt)
V_LES_solution = zeros(Nz, Nt)
T_LES_solution = zeros(Nz, Nt)

U′W′_LES_solution = zeros(Nz+1, Nt)
V′W′_LES_solution = zeros(Nz+1, Nt)
W′T′_LES_solution = zeros(Nz+1, Nt)
Ri_LES_solution = zeros(Nz+1, Nt)

for (n, iter) in zip(1:length(times), keys(ds["timeseries/t"]))
    U_LES_solution[:, n] = ds["timeseries/u/$iter"]
    V_LES_solution[:, n] = ds["timeseries/v/$iter"]
    T_LES_solution[:, n] = ds["timeseries/T/$iter"]

    U′W′_LES_solution[:, n] = ds["timeseries/wu/$iter"]
    V′W′_LES_solution[:, n] = ds["timeseries/wv/$iter"]
    W′T′_LES_solution[:, n] = ds["timeseries/wT/$iter"]
end

U′W′_LES_solution[end, :] .= Fu
W′T′_LES_solution[end, :] .= Fθ

zc = models[1].grid.zc
zf = models[1].grid.zf

frame = Node(1)
plot_title = @lift @sprintf("Modified Diffusivity %s: time = %s", SIMULATION_NAME, prettytime(times[$frame]))

U_models = [@lift U_solution[:, $frame] for U_solution in U_solutions]
V_models = [@lift V_solution[:, $frame] for V_solution in V_solutions]
T_models = [@lift T_solution[:, $frame] for T_solution in T_solutions]

U′W′_models = [@lift U′W′_solution[:, $frame] for U′W′_solution in U′W′_solutions]
V′W′_models = [@lift V′W′_solution[:, $frame] for V′W′_solution in V′W′_solutions]
W′T′_models = [@lift W′T′_solution[:, $frame] for W′T′_solution in W′T′_solutions]
Ri_models = [@lift Ri_solution[:, $frame] for Ri_solution in Ri_solutions]

U_LES = @lift U_LES_solution[:, $frame]
V_LES = @lift V_LES_solution[:, $frame]
T_LES = @lift T_LES_solution[:, $frame]

U′W′_LES = @lift U′W′_LES_solution[:, $frame]
V′W′_LES = @lift V′W′_LES_solution[:, $frame]
W′T′_LES = @lift W′T′_LES_solution[:, $frame]

fig = Figure(resolution=(1920, 1080))

colors = ["dodgerblue2", "forestgreen", "darkorange", "navyblue", "hotpink2", "saddlebrown", "gray27", "lawngreen"]

ax_U = fig[1, 1] = Axis(fig, xlabel="U (m/s)", ylabel="z (m)")
U_lines = cat(lines!(ax_U, U_LES, zc, linewidth=3, color="crimson"),
                [lines!(ax_U, U_models[i], zc, linewidth=3, color=colors[i]) for i in 1:length(U_models)]..., dims=1)
xlims!(ax_U, -1, 1)
ylims!(ax_U, -Lz, 0)

ax_V = fig[1, 2] = Axis(fig, xlabel="V (m/s)", ylabel="z (m)")
V_lines = cat(lines!(ax_V, V_LES, zc, linewidth=3, color="crimson"),
                [lines!(ax_V, V_models[i], zc, linewidth=3, color=colors[i]) for i in 1:length(V_models)]..., dims=1)
xlims!(ax_V, -1, 1)
ylims!(ax_V, -Lz, 0)

ax_T = fig[1, 3] = Axis(fig, xlabel="T (°C)", ylabel="z (m)")
T_lines = cat(lines!(ax_T, T_LES, zc, linewidth=3, color="crimson"),
                [lines!(ax_T, T_models[i], zc, linewidth=3, color=colors[i]) for i in 1:length(T_models)]..., dims=1)
ylims!(ax_T, -Lz, 0)

ax_U′W′ = fig[2, 1] = Axis(fig, xlabel="U′W′ (m²/s²)", ylabel="z (m)")
U′W′_lines = cat(lines!(ax_U′W′, U′W′_LES, zf, linewidth=3, color="crimson"),
                    [lines!(ax_U′W′, U′W′_models[i], zf, linewidth=3, color=colors[i]) for i in 1:length(U′W′_models)]..., dims = 1)
xlims!(ax_U′W′, extrema(U′W′_LES_solution))
ylims!(ax_U′W′, -Lz, 0)

ax_V′W′ = fig[2, 2] = Axis(fig, xlabel="V′W′ (m²/s²)", ylabel="z (m)")
V′W′_lines = cat(lines!(ax_V′W′, V′W′_LES, zf, linewidth=3, color="crimson"),
                    [lines!(ax_V′W′, V′W′_models[i], zf, linewidth=3, color=colors[i]) for i in 1:length(V′W′_models)]..., dims = 1)
xlims!(ax_V′W′, all(V′W′_LES_solution .== 0) ? (-1, 1) : extrema(V′W′_LES_solution))
ylims!(ax_V′W′, -Lz, 0)

ax_W′T′ = fig[2, 3] = Axis(fig, xlabel="W′T′ (m/s ⋅ °C)", ylabel="z (m)")
W′T′_lines = cat(lines!(ax_W′T′, W′T′_LES, zf, linewidth=3, color="crimson"),
                    [lines!(ax_W′T′, W′T′_models[i], zf, linewidth=3, color=colors[i]) for i in 1:length(W′T′_models)]..., dims = 1)
xlims!(ax_W′T′, extrema(W′T′_LES_solution))
ylims!(ax_W′T′, -Lz, 0)

ax_Ri = fig[2, 4] = Axis(fig, xlabel="Richardson number", ylabel="z (m)")
Ri_lines = [lines!(ax_Ri, Ri_models[i], zf, linewidth=3, color=colors[i]) for i in 1:length(Ri_models)]
xlims!(ax_Ri, extrema(Ri_solutions[end]))
ylims!(ax_Ri, -Lz, 0)

legend = fig[1, 4] = Legend(fig, U_lines, model_names)

supertitle = fig[0, :] = Label(fig, plot_title, textsize=30)
trim!(fig.layout)

record(fig, joinpath(pwd(), "Output/$OUTPUT_NAME"), 1:length(times), framerate=15) do n
    @info "Animating Pacanowski-Philander wind-mixing frame $n/$(length(times))..."
    frame[] = n
end
##


text:
################################################################################
#                                   Si₃N₄                                      #
################################################################################
export Si₃N₄

"""
This code first creates a symbolic representation of the
Sellmeier Equation model for the index of refraction.
Equation form is based on Luke, Okawachi, Lamont, Gaeta and Lipson,
"Broadband mid-infrared frequency comb generation in a Si3N4 microresonator"
Optics Letters Vol. 40, Issue 21, pp. 4823-4826 (2015)
https://doi.org/10.1364/OL.40.004823
valid from 0.31–5.504 um
Thermo-optic coefficients from
Xue, et al.
"Thermal tuning of Kerr frequency combs in silicon nitride microring resonators"
Opt. Express 24.1 (2016) http://doi.org/10.1364/OE.24.000687
The symbolic index model and its derivatives in turn are used to generate
numerical functions for the SiO₂ index, group index and GVD as a function
of temperature and wavelength.
Variable units are lm in [um] and T in [deg C]
"""

p_n²_Si₃N₄= (
    A₀ = 1,
    B₁ = 3.0249,
    C₁ = (0.1353406)^2,         #                           [μm²]
    B₂ = 40314,
    C₂ = (1239.842)^2,          #                           [μm²]
    dn_dT = 2.96e-5,            # thermo-optic coefficient  [K⁻¹]
    T₀ = 24.5,                  # reference temperature     [°C]
)

function make_Si₃N₄(;p_n²=p_n²_Si₃N₄)
	@variables λ, T
	n = sqrt(n²_sym_fmt1( λ ; p_n²...)) + p_n².dn_dT  *  ( T - p_n².T₀  )
	n² = n^2
	ε 	= diagm([n², n², n²])
	ng = ng_model(n,λ)
	gvd = gvd_model(n,λ)
	models = Dict([
		:n		=>	n,
		:ng		=>	ng,
		:gvd	=>	gvd,
		:ε 		=> 	diagm([n², n², n²]),
	])
	defaults =	Dict([
		:λ		=>		0.8,		# μm
		:T		=>		p_n².T₀,	# °C

	])
	Material(models, defaults, :Si₃N₄, colorant"firebrick1")
end

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

Si₃N₄ = make_Si₃N₄()


text:
"""
    adj = read_graph(fn)

Read a graph from a file in IJ or IJV format.
That is, each line of the file should represent one edge.
Each edge should be specified by the indices of its vertices,
separated by a whitespace or a comma.  If the graph is weighted,
the weight should follow the second index.  For example, the unweighted
complete graph on 3 vertices would appear as the file

```
1 2
1 3
2 3
```

A weighted path on 3 vertices with edge weights `1.5` and `2.5` would be

```
1, 2, 1.5
2, 3, 2.5
```

The function tries to detect the delimiter type (comma or whitespace)
from the first line of the file.  The format must be consistent.
Vertex indices start at `1`.
"""
function read_graph(fn::AbstractString)

    # first, detect the delimiter type
    fh = open(fn)
    ln = readline(fh)
    close(fh)

    if occursin(",", ln)
        data = readdlm(fn,',')
    else
        data = readdlm(fn)
    end

    n = maximum(data[:,1:2])

    edlist = convert(Array{Int64,2}, data[:,1:2])

    if size(data,2) == 3
        wts = convert(Array{Float64,1},data[:,3])
    else
        wts = ones(size(data,1))
    end

    a = sparse(edlist[:,1],edlist[:,2],wts,n,n)
    a = a + a'

    return a
end

"""Writes the upper portion of a matrix in ijv format, one row for each edge,
separated by commas.  Only writes the upper triangular portion.
The result can be read from Matlab like this:

```
>> dl = dlmread('graph.txt');
>> a = sparse(dl(:,1),dl(:,2),dl(:,3));
>> n = max(size(a))
>> a(n,n) = 0;
>> a = a + a';
```
"""
function writeIJV(filename::AbstractString, mat)

  (ai,aj,av) = findnz(triu(mat))
  fh = open(filename,"w")
  for i in 1:length(ai)
    write(fh, "$(ai[i]),$(aj[i]),$(av[i])\n")
  end
  close(fh)

end #write IJV


text:
innerproduct(f, g, measure; options...) = integral(x->conj(f(x))*g(x), measure)

applymeasure(μ::Measure, f::Function; options...) = integral(f, μ)

"A measure on a general domain with a general weight function `dσ = w(x) dx`."
struct GenericWeight{T} <: Weight{T}
    support          ::  Domain{T}
    weightfunction
end

name(m::GenericWeight) = "Measure with generic weight function"

unsafe_weightfun(m::GenericWeight, x) = m.weightfunction(x)

support(m::GenericWeight) = m.support

strings(m::GenericWeight) = (name(m), (string(m.support),), (string(m.weightfunction),))

name(m::LebesgueDomain) = "Lebesgue measure"
name(m::LegendreWeight) = "Legendre weight"
name(m::JacobiWeight) = "Jacobi weight (α = $(m.α), β = $(m.β))"
name(m::LaguerreWeight) = m.α == 0 ? "Laguerre weight" : "Generalized Laguerre measure (α = $(m.α))"
name(m::HermiteWeight) = "Hermite weight"
name(m::ChebyshevTWeight) = "Chebyshev weight of the first kind"
name(m::ChebyshevUWeight) = "Chebyshev weight of the second kind"

const FourierWeight = DomainIntegrals.LebesgueUnit
name(m::FourierWeight) = "Fourier (Lebesgue) measure"

######################################################
# Generating new measures from existing measures
######################################################

using DomainIntegrals: MappedWeight, ProductWeight

import DomainIntegrals: supermeasure, mappedmeasure

name(m::MappedWeight) = "Mapped measure"
strings(m::MappedWeight) = (name(m), strings(forward_map(m)), strings(supermeasure(m)))

apply_map(μ::Measure, m) = MappedWeight(m, μ)

function stencilarray(m::ProductWeight)
    A = Any[]
    push!(A, component(m,1))
    for i = 2:length(components(m))
        push!(A," ⊗ ")
        push!(A, component(m,i))
    end
    A
end

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

#############################
# Compatibility of measures
#############################

# By default, measures are compatible only when they are equal.
iscompatible(m1::M, m2::M) where {M <: Measure} = m1==m2
iscompatible(m1::Weight, m2::Weight) = false


text:
## using KronLinInv

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

function test2D( )

    ## 2D problem, so set nx = 1
    nx = 1
    ny = 20
    nz = 30
    nxobs = 1
    nyobs = 18
    nzobs = 24
    sigmaobs  = [1.0, 0.1, 0.1]
    corlenobs = [0.0, 1.4, 1.4]
    sigmam    = [1.0, 0.8, 0.8]
    corlenm   = [0.0, 2.5, 2.5]

    # run test
    result = test_mean_cov(nx,ny,nz,nxobs,nyobs,nzobs,
                           sigmaobs,corlenobs,sigmam,corlenm)

    return result
end

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

function test3D( )

    ## 2D problem, so set nx = 1
    nx = 8
    ny = 9
    nz = 12
    nxobs = 6
    nyobs = 7
    nzobs = 11
    sigmaobs  = [0.1, 0.1, 0.1]
    corlenobs = [1.4, 1.4, 1.4]
    sigmam    = [0.8, 0.8, 0.8]
    corlenm   = [2.5, 2.5, 2.5]

    # run test
    result = test_mean_cov(nx,ny,nz,nxobs,nyobs,nzobs,
                           sigmaobs,corlenobs,sigmam,corlenm)

    return result
end

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

function test_mean_cov(nx::Integer,ny::Integer,nz::Integer,
                       nxobs::Integer,nyobs::Integer,nzobs::Integer,
                       sigmaobs::Array{Float64,1},corlenobs::Array{Float64,1},
                       sigmam::Array{Float64,1},corlenm::Array{Float64,1})

    ##############################
    ## Setup the problem
    ##############################

    ## Covariance matrices
    # covariance on observed data
    Cd1,Cd2,Cd3 = mkCovSTAT(sigmaobs,nxobs,nyobs,nzobs,corlenobs,"gaussian")
    # covariance on model parameters
    Cm1,Cm2,Cm3 = mkCovSTAT(sigmam,nx,ny,nz,corlenm,"gaussian")

    ## Create the covariance matrix structure
    Covs = CovMats(Cd1,Cd2,Cd3,Cm1,Cm2,Cm3)

    ## Forward model operator
    G1 = rand(nxobs,nx)
    G2 = rand(nyobs,ny)
    G3 = rand(nzobs,nz)

    ## Create the forward operator structure
    Gfwd = FwdOps(G1,G2,G3)

    ##############################
    ## Setup the synthetic test
    ##############################

    ## Create a reference model
    refmod = rand(nx*ny*nz)

    ## Create a reference model
    mprior = copy(refmod) #0.5 .* ones(nx*ny*nz)

    ## Create some "observed" data
    ##   (without noise because it's just a test of the algo)
    dobs = kron(G1,kron(G2,G3)) * refmod

    ##############################
    ## Solve the inverse problem
    ##############################

    ## Calculate the required factors
    klifac = calcfactors(Gfwd,Covs)

    ## Calculate the posterior mean model
    postm = posteriormean(klifac,Gfwd,mprior,dobs)

    # serial version
    postm2 = KronLinInv.posteriormean_serial(klifac,Gfwd,mprior,dobs)

    ## Calculate the posterior covariance
    npts = nx*ny*nz
    astart, aend = 1,div(npts,3)
    bstart, bend = 1,div(npts,3)
    postC = blockpostcov(klifac,astart,aend,bstart,bend)

    ## serial version
    postC2 = KronLinInv.blockpostcov_serial(klifac,astart,aend,bstart,bend)

    ##############################
    ## Check results
    ##############################

    if (postm ≈ refmod) && (postm2 ≈ refmod)
        return true
    end
    return false
end

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

function mkCovSTAT(sigma::Array{Float64,1},nx::Integer,ny::Integer,nz::Integer,
                   corrlength::Array{Float64,1},kind::String)
    #
    #   Stationaly covariance model
    #
    #println("Creating covariance matrix...")

    function cgaussian(dist,corrlength)
        if maximum(dist)==0.0
            return 1.0
        else
            @assert(corrlength>0.0)
            return exp.(-(dist./corrlength).^2)
        end
    end
    function cexponential(dist,corrlength)
        if maximum(dist)==0.0
            return 1.0
        else
            @assert(corrlength>0.0)
            return exp.(-(dist./corrlength))
        end
    end

    npts = nx*ny*nz
    x = [float(i) for i=1:nx]
    y = [float(i) for i=1:ny]
    z = [float(i) for i=1:nz]
    covmat_x = zeros(nx,nx)
    covmat_y = zeros(ny,ny)
    covmat_z = zeros(nz,nz)

    if kind=="gaussian"
        calccovfun = cgaussian
    elseif kind=="exponential"
        calccovfun = cexponential
    else
        println("Error, no or wrong cov 'kind' specified")
        exit()
    end

    for i=1:nx
        covmat_x[i,:] .= sigma[1]^2 .*
            calccovfun(sqrt.((x.-x[i]).^2),corrlength[1])
    end
    for i=1:ny
        covmat_y[i,:] .= sigma[2]^2 .*
            calccovfun(sqrt.(((y.-y[i])).^2),corrlength[2])
    end
    for i=1:nz
        covmat_z[i,:] .= sigma[3]^2 .*
            calccovfun(sqrt.(((z.-z[i])).^2),corrlength[3])
    end

    return covmat_x,covmat_y,covmat_z
 end

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


text:
module AutoGainCtrl

# offline part
export agc!
export agc

# online part
export setagc
export onlineagc!
export onlineagc

include("offline-agc.jl")
include("online-agc.jl")

end #module

max_stars_repo_name: sonosole/AGC

max_issues_repo_name: sonosole/AutoGainCtrl

max_forks_repo_name: sonosole/AutoGainCtrl


text:
include("common.jl")
import Base: &, adjoint

mutable struct Equation
	lhs#::EX #relaxed this to allow arbitrary substitution, for example a tensor Ten(A) where A is replaced with an array. Do note that raw arrays should never be present in expressions, unless you are a practitioner of Chaos Magic (ie UNDEFINED results).
	rhs#::EX
	divisions
end
Equation(ex1,ex2)=Equation(ex1,ex2,Set())
Equation(ex1)=Equation(ex1,0,Set())
function tosym(expr)
	if isa(expr,Symbol)
		return QuoteNode(:($expr))
	elseif isa(expr,Expr)
		if expr.head==:$
			expr=esc(expr.args[1])
		elseif expr.head==:vcat
			for s in 1:length(expr.args)
				for p in 1:length(expr.args[s].args)
					expr.args[s].args[p]=tosym(expr.args[s].args[p])
				end
			end
		else
			if expr.head==:vect
				s1=1
			else
				s1=2
			end
			for s in s1:length(expr.args)
				expr.args[s]=tosym(expr.args[s])
			end
		end
	end
	return expr
end
macro equ(eq)
	neq=Expr(:comparison)
	push!(neq.args,tosym(eq.args[1]),≖,tosym(eq.args[2]))
	neq
end
macro equs(eqs...)
	l=length(eqs)
	eqa=Expr(:vect)
	for i in 1:l
		neq=Expr(:comparison)
		push!(neq.args,tosym(eqs[i].args[1]),≖,tosym(eqs[i].args[2]))
		push!(eqa.args,neq)
	end
	return eqa
end
function print(io::IO,eq::Equation)
	print(io,eq.lhs)
	print(io," = ")
	print(io,eq.rhs)
end
≖(a,b)=Equation(a,b)
complexity(eq::Equation)=complexity(eq.lhs)+complexity(eq.rhs)
isless(eq1::Equation,eq2::Equation)=complexity(eq1)<complexity(eq2)
mutable struct EquationChain
	expressions::Vector
end
EquationChain(a...)=EquationChain([a...])
start(eqc::EquationChain)=(1,eqc.expressions)
function next(eqc::EquationChain,state)
	return (state[2][state[1]],(state[1]+1,state[2]))
end
done(eqc::EquationChain,state)=state[1]>length(state[2])
getindex(eqc::EquationChain,i::Integer)=getindex(eqc.expressions,i)
setindex!(eqc::EquationChain,a)=setindex!(eqc.expressions,a)
length(eqc::EquationChain)=length(eqc.expressions)
push!(eqc::EquationChain,a)=push!(eqc.expressions,a)
function print(io::IO,eqc::EquationChain)
	for exi in 1:length(eqc)-1
		print(io,eqc[exi])
		print(io," ≖ ")
	end
	print(io,eqc[length(eqc)])
end
≖(a::EX,b::Equation)=EquationChain(a,b.lhs,b.rhs)
≖(b::Equation,a::EX)=EquationChain(b.lhs,b.rhs,a)
≖(eqc::EquationChain,a)=push!(eqc,a)
-(eq::Equation)=simplify(-eq.lhs≖-eq.rhs)
+(eq::Equation,ex::EX)=simplify(eq.lhs+ex≖eq.rhs+ex)
-(eq::Equation,ex::EX)=simplify(eq.lhs-ex≖eq.rhs-ex)
*(eq::Equation,ex::EX)=simplify(eq.lhs*ex≖eq.rhs*ex)
/(eq::Equation,ex::EX)=simplify(eq.lhs/ex≖eq.rhs/ex)
+(eq1::Equation,eq2::Equation)=simplify(eq1.lhs+eq2.lhs≖eq1.rhs+eq2.rhs)
-(eq1::Equation,eq2::Equation)=simplify(eq1.lhs-eq2.lhs≖eq1.rhs-eq2.rhs)
*(eq1::Equation,eq2::Equation)=simplify(eq1.lhs*eq2.lhs≖eq1.rhs*eq2.rhs)
/(eq1::Equation,eq2::Equation)=simplify(eq1.lhs/eq2.lhs≖eq1.rhs/eq2.rhs)
equation(ex::EX)=Equation(ex,0,Any[])
equation(ex1::EX,ex2::EX)=Equation(ex1,ex2,Any[])
adjoint(eq::Equation)=Equation(eq.rhs,eq.lhs)
==(eq1::Equation,eq2::Equation)=eq1.lhs==eq2.lhs&&eq1.rhs==eq2.rhs
function (&)(eq1::Equation,eq2::Equation)
	eq1=simplify(eq1);eq2=simplify(eq2)
	neq=simplify(Equation(replace(eq1.lhs,Dict(eq2.lhs=>eq2.rhs)),replace(eq1.rhs,Dict(eq2.lhs=>eq2.rhs))))
	if neq==eq1
		neq=Equation(eq1.lhs,simplify(eq1.rhs&eq2))
	end
	return neq
end
function (&)(eq::Equation,eqa::Array{Equation})
	for teq in eqa
		eq=eq&teq
	end
	return simplify(eq)
end
(&)(eq::Equation,fun::Function)=fun(eq)
(&)(ex::EX,fun::Function)=fun(ex)
function (&)(ex::Expression,eq::Equation)
	ex=simplify(ex);eq=simplify(eq)
	if isa(eq.lhs,Symbol)
		return simplify(replace(ex,Dict(eq.lhs=>eq.rhs)))
	end
	m=matches(ex,eq)
	if !isempty(m)
		return simplify(m[1])
	else
		tms=Tuple[]
		_terms=dcterms(ex)
		for t in 1:length(_terms)
			tm=matches(simplify(expression(ex[t])),eq)
			if !isempty(tm)
				_terms[t]=tm[1][1]
			else
				for f in 1:length(_terms[t])
					if isa(_terms[t][f],Component)
						ttm=matches(_terms[t][f],eq)
						if !isempty(ttm)
							_terms[t][f]=ttm[1]
						end
					end
				end
			end
		end
		return simplify(expression(_terms))
	end
end
function (&)(ex::Component,eq::Equation)
	ex=simplify(ex);eq=simplify(eq)
	if isa(eq.lhs,Symbol)
		return simplify(replace(ex,Dict(eq.lhs=>eq.rhs)))
	end
	m=matches(ex,eq)
	if !isempty(m)
		return simplify(m[1])
	else
		return ex
	end
end
function (&)(ex::Symbol,eq::Equation)
	if isa(eq.lhs,Symbol)&&eq.lhs==ex
		return eq.rhs
	end
#	m=matches(ex,eq)
#	if !isempty(m)
#		return simplify(m[1])
#	end
	return ex
end
(&)(x::Number,eq::Equation)=x
function (&)(ex::Union{Ex,Equation},eqa::Array)
	for teq in eqa
		ex=ex&teq
	end
	return simplify(ex)
end
function (&)(eq::Equation,sym::Symbol)
	if isa(eq.lhs,Component)
		s=solve(eq,sym,typeof(eq.lhs))
	else
		s=solve(eq,sym)
	end
	if isa(s,Nothing)
		error("Please implement code to solve $eq for $sym")
	end
#	if isa(s.lhs,Component) #how to avoid infinite recursion?
#		s=s&sym
#	end
	return s
end
function solve(c::Component,sym::Symbol,t::Type)
	print_with_color(:green,"Please implement a solver for $(typeof(c))")
	c
end
function equivalent(eq1::Equation,eq2::Equation)
	m=matches(eq2)
	for eq in m
		if eq1==eq
			return true
		end
	end
	return false
end
simplify!(eq::Equation)=begin;eq.lhs=simplify!(eq.lhs);eq.rhs=simplify!(eq.rhs);eq;end #the ! functions are not complete
function simplify(eq::Equation)
	lhs,rhs=simplify(eq.lhs),simplify(eq.rhs)
#	if (isa(rhs,Symbol)&&!isa(lhs,Symbol))||(isa(lhs,Number)&&!isa(rhs,Number)) #this is convenient sometimes but causes breakage
#		return Equation(rhs,lhs)
#	else
		return Equation(lhs,rhs)
#	end
end
function simplify!(eqa::Array{Equation})
	for eq in 1:length(eqa)
		eqa[eq]=simplify(eqa[eq])
	end
	return eqa
end
simplify(eqa::Array{Equation})=simplify!(deepcopy(eqa))
function pushallunique!(a1::Array,a2::Array)
	for d in a2
		if !(d∈a1)
			push!(a1,d)
		end
	end
	return a1
end
function pushallunique!(a1::Array,d)
	if !(d∈a1)
		push!(a1,d)
	end
	return a1
end
pushallunique!(a1::Array,b::Bool)=Void
function uniquefilter(a::Array{T}) where {T}
	pushallunique!(T[],a)
end
include("matchers.jl")
matchfuns=Function[]
push!(matchfuns,quadratic)
function matches(eq::Equation)
	m=Equation[]
	if eq.lhs==0
		return m #it only moves from left to right
	end
	eq=simplify(eq)
	terms=dcterms(eq.lhs)
	for term in 1:length(terms)
		teq=deepcopy(eq)
		tt=deepcopy(terms)
		if typeof(teq.rhs)==Expression
			nterm=Factor[]
			push!(nterm,-1)
			for fac in terms[term]
				push!(nterm,fac)
			end
			push!(teq.rhs,nterm)
		else
			teq.rhs=Expression(Term[[teq.rhs],[-1;terms[term]]])
		end
		if !isa(tt,X)
			deleteat!(tt,term)
		end
		teq.lhs=expression(tt)
		push!(m,teq)
	end
	for teq in m
		teq.rhs=sumnum(componify(teq.rhs))
		teq.lhs=sumnum(componify(teq.lhs))
	end
	return m
end
matches(ex::EX)=matches(equation(ex))
function matches(eq::Equation,fun::Function)
	m=Equation[]
	pushallunique!(m,fun(eq))
	return m
end
matches(ex::EX,fun::Function)=matches(equation(ex),fun)
function matches(eq::Equation,all::Bool)
	if eq.lhs==0
		return false #it only moves from left to right
	end
	eq=simplify(eq)
	m=Equation[]
	for fun in matchfuns
		pushallunique!(m,fun(eq))
	end
	pushallunique!(m,matches(eq,Div))
	pushallunique!(m,matches(eq,Sqrt))
	pushallunique!(m,matches(eq))
	return m
end
matches(ex::EX,all::Bool)=matches(equation(ex),all)
function matches(eqa::Array{Equation})
	neqa=deepcopy(eqa)
	for eq in eqa
		m=matches(eq)
		for teq in m
			push!(neqa,teq)
		end
	end
	return neqa
end
function matches!(eqa::Array{Equation})
	leqa=length(eqa)
	for eq in 1:leqa
		m=matches(eqa[eq])
		for teq in m
			if teq!=false&&indin(eqa,teq)==0
				push!(eqa,teq)
			end
		end
	end
	return eqa
end
function matches(eq::Equation,recursions::Integer)
	m=matches(eq)
	for r in 1:recursions
		matches!(m)
	end
	return m
end
matches(ex::EX,rec::Integer)=matches(equation(ex),rec)
include("div.jl")
include("sqrt.jl")
include("pow.jl")
include("log.jl")
include("der.jl")
include("vec.jl")
function evaluate(eq::Equation,symdic::Dict)
	for key in keys(symdic)
		if symdic[key]==0&&key∈eq.divisions
			error("Assigning zero to a value that the equation has been divided by.")
		end
	end
	return (evaluate(eq.lhs,symdic),evaluate(eq.rhs,symdic))
end
function solve(eq::Equation,sym::Symbol)
	eq=simplify(eq)
	indsl=expandindices(indsin(eq.lhs,sym))
	indsh=expandindices(indsin(eq.rhs,sym))
	if length(indsl)==0||length(indsh)==0
		if length(indsl)==0
			inds=indsh
			lhs=eq.rhs
			rhs=eq.lhs
		else
			inds=indsl
			lhs=eq.lhs
			rhs=eq.rhs
		end
		if length(inds)==1
			inds=inds[1]
			for termi in 1:length(lhs)
				if termi==inds[1]
					continue
				end
				rhs=rhs-expression(lhs[termi])
			end
			lhs=lhs[inds[1]]
			for fac in lhs
				if fac==lhs[inds[2]]
					continue
				end
				rhs=rhs/fac
			end
			lhs=lhs[inds[2]]
			return simplify(Equation(lhs,rhs))
		end
	end
end
function solve(eq::Equation,rec::Integer=1)
	seq=simplify(eq)
	mat=matches(seq,rec)
	sol=Equation[]
	for m in mat
		if isa(m.lhs,Symbol)
			push!(sol,m)
		end
	end
	return sol
end
function solve(eq::Equation,op)
	seq=simplify(eq)
	mat=matches(seq,op)
	sol=Equation[]
	for m in mat
		if length(m.lhs)==1
			push!(sol,m)
		end
	end
	return sol
end
solve(ex::Ex)=solve(equation(ex))
solve(ex::Ex,a)=solve(equation(ex),a)

include("units.jl")
include("printers.jl")


text:
module Conv
	function conv_func(t::DataType)
		return (x -> convert(t, x))
	end
	function conv_func(typestr::String)
		return eval(parse(" x -> convert($(typestr), x)"))
	end
	function conv(t::DataType, v::Any)
		return convert(t, v)
	end
	function conv(typestr::String, v::Any)
		return conv_func(typestr)(v)
	end
	function conv_type(typestr::String)
		return eval(parse("$typestr"))
	end
#
	type Converter
		f::Function
		t::DataType
	end
#
	function converter(t::DataType)
		return Converter(conv_func(t), t)
	end
	function converter(typestr::String)
		return converter(conv_type(typestr))
	end
	function vector(conv::Converter, V::Vector)
		n = length(V)
		ret = Array(conv.t, n)
		for i = 1:n
			ret[i] = conv.f(V[i])
		end
		return ret
	end
	function matrix(conv::Converter, M::Matrix)
		ret = Array(conv.t, (size(M,1), size(M,2)))
		for i = 1:length(M)
			ret[i] = conv.f(M[i])
		end
		return ret
	end
	function matrix(conv::Converter, M::SparseMatrixCSC)
		I,J,V = findnz(M)
		Vc = Array(conv.t, length(V))
		for i = 1:length(V)
			Vc[i] = conv.f(V[i])
		end
		return sparse(I,J,Vc)
	end
end

max_stars_repo_name: jadnohra/mlcp-lib

max_issues_repo_name: jadnohra/mlcp-lib

max_forks_repo_name: jadnohra/mlcp-lib


text:
# This file investigates mean square error (MSE) with respect to free variables in hidden 1D interpolation

using FractalTools
using Makie

# Construct interpolation data
# f(x) = [FractalTools.sinusoid(x), FractalTools.sinusoid(x)]  # For regular data
f(x) = [FractalTools.wen(x), FractalTools.weierstrass(x)]  # For irregular data
line = [[0.], [1.]]
pts = getdata(f, line, 11)

# Construct test data
tpts = getdata(line, 101)
npts = length(tpts)

# Compute errors
fvals = getindex.(map(pnt -> f(pnt[1]), tpts), 1)
freevars = 0.001 : 0.01 : 0.5
mse = map(freevars) do freevar
    interp = interpolate(pts, HInterp1D(fill(freevar, 2, 2)))
    ivals = getindex.(map(pnt -> interp(pnt[1]), tpts), 1)
    sum((fvals - ivals).^2) / npts
end

# Plot mse
fig = Figure()
ax = fig[1, 1] = Axis(fig, xlabel="Free Variable", ylabel="MSE", title="1D Hidden Interpolation MSE")
stem!(ax, freevars, mse, color=:black)
save(joinpath(@__DIR__, "hinterp1d_mse.png"), fig)
display(fig)


text:
function isvalid(subject::Dict{String, String}, display = false)
    if display
        println("isvalid: $subject")
    end
    typeof(iban(subject["value"])) == Dict{String,String}
end

@testset "iban_random         [one]" begin
    @test iban_random(
        CountryCode="BR"
    )["CountryCode"] == "BR"

    @test iban_random(
        CountryCode="BR",
        BankCode="18731592"
    )["BankCode"] == "18731592"

    @test iban_random(
        CountryCode="BR",
        BankCode="18731592",
        BranchCode = "89800"
    )["BranchCode"] == "89800"

    @test iban_random(
        CountryCode="BR",
        BankCode="18731592",
        BranchCode = "89800",
        AccountNumber = "6427460610"
    )["AccountNumber"] == "6427460610"

    @test iban_random(
        CountryCode="BR",
        BankCode="18731592",
        BranchCode = "89800",
        AccountNumber = "6427460610",
        AccountType = "X"
    )["AccountType"] == "X"

    @test iban_random(
        CountryCode="BR",
        BankCode="18731592",
        BranchCode = "89800",
        AccountNumber = "6427460610",
        AccountType = "X",
        OwnerAccountType = "m",
    )["OwnerAccountType"] == "m"
end

@testset "iban_random         [all]" begin
     @test isvalid(iban_random()) == true
     @test isvalid(iban_random(CountryCode="BR", BankCode="18731592")) == true

      # all supported country codes
      map(cc ->
        (@test isvalid(iban_random(CountryCode=cc)) == true),
        supported_countries()
      )
end


text:
using PkgBenchmark

#--track-allocation=all
config = BenchmarkConfig(id = nothing,
    juliacmd = `julia -O3`,
    env = Dict("JULIA_NUM_THREADS" => 4))

results = benchmarkpkg("Manifolds", config)
export_markdown("benchmark/results.md", results)


text:
# OpenCL argument utilities

# Argument accessors

# __agocl_global_offset_x - OpenCL Global Offset X
# __agocl_global_offset_y - OpenCL Global Offset Y
# __agocl_global_offset_z - OpenCL Global Offset Z
# __agocl_printf_addr - OpenCL address of printf buffer
# __agocl_queue_addr - OpenCL address of virtual queue used by enqueue_kernel
# __agocl_aqlwrap_addr - OpenCL address of AqlWrap struct used by enqueue_kernel
# __agocl_multigrid - Pointer argument used for Multi-grid synchronization