text:
const MSB = 0x80
const MASK7 = 0x7f
const MASK8 = 0xff
const MASK3 = 0x07

const WIRETYP_VARINT   = 0
const WIRETYP_64BIT    = 1
const WIRETYP_LENDELIM = 2
const WIRETYP_GRPSTART = 3   # deprecated
const WIRETYP_GRPEND   = 4   # deprecated
const WIRETYP_32BIT    = 5

"""
The abstract type from which all generated protobuf structs extend.
"""
abstract type ProtoType end

# TODO: wiretypes should become julia types, so that methods can be parameterized on them
const WIRETYPES = Dict{Symbol,Tuple}(
    :int32          => (WIRETYP_VARINT,     :write_varint,  :read_varint,   Int32),
    :int64          => (WIRETYP_VARINT,     :write_varint,  :read_varint,   Int64),
    :uint32         => (WIRETYP_VARINT,     :write_varint,  :read_varint,   UInt32),
    :uint64         => (WIRETYP_VARINT,     :write_varint,  :read_varint,   UInt64),
    :sint32         => (WIRETYP_VARINT,     :write_svarint, :read_svarint,  Int32),
    :sint64         => (WIRETYP_VARINT,     :write_svarint, :read_svarint,  Int64),
    :bool           => (WIRETYP_VARINT,     :write_bool,    :read_bool,     Bool),
    :enum           => (WIRETYP_VARINT,     :write_varint,  :read_varint,   Int32),

    :fixed64        => (WIRETYP_64BIT,      :write_fixed,   :read_fixed,    UInt64),
    :sfixed64       => (WIRETYP_64BIT,      :write_fixed,   :read_fixed,    Int64),
    :double         => (WIRETYP_64BIT,      :write_fixed,   :read_fixed,    Float64),

    :string         => (WIRETYP_LENDELIM,   :write_string,  :read_string,   AbstractString),
    :bytes          => (WIRETYP_LENDELIM,   :write_bytes,   :read_bytes,    Vector{UInt8}),
    :obj            => (WIRETYP_LENDELIM,   :writeproto,    :readproto,     Any),
    :map            => (WIRETYP_LENDELIM,   :write_map,     :read_map,      Dict),

    :fixed32        => (WIRETYP_32BIT,      :write_fixed,   :read_fixed,    UInt32),
    :sfixed32       => (WIRETYP_32BIT,      :write_fixed,   :read_fixed,    Int32),
    :float          => (WIRETYP_32BIT,      :write_fixed,   :read_fixed,    Float32)
)

aliaswiretypes(wtype::Symbol) = wiretypes(WIRETYPES[wtype][4])

wiretypes(::Type{Int32})                            = [:int32, :sint32, :enum, :sfixed32]
wiretypes(::Type{Int64})                            = [:int64, :sint64, :sfixed64]
wiretypes(::Type{UInt32})                           = [:uint32, :fixed32]
wiretypes(::Type{UInt64})                           = [:uint64, :fixed64]
wiretypes(::Type{Bool})                             = [:bool]
wiretypes(::Type{Float64})                          = [:double]
wiretypes(::Type{Float32})                          = [:float]
wiretypes(::Type{T}) where {T<:AbstractString}      = [:string]
wiretypes(::Type{Vector{UInt8}})                    = [:bytes]
wiretypes(::Type{Dict{K,V}}) where {K,V}            = [:map]
wiretypes(::Type)                                   = [:obj]
wiretypes(::Type{Vector{T}}) where {T}              = wiretypes(T)

wiretype(::Type{T}) where {T}                       = wiretypes(T)[1]

defaultval(::Type{T}) where {T<:Number}             = [zero(T)]
defaultval(::Type{T}) where {T<:AbstractString}     = [convert(T,"")]
defaultval(::Type{Bool})                            = [false]
defaultval(::Type{Vector{T}}) where {T}             = Any[T[]]
defaultval(::Type)                                  = []

function _write_uleb(io::IO, x::T) where T <: Integer
    nw = 0
    cont = true
    while cont
        byte = x & MASK7
        if (x >>>= 7) != 0
            byte |= MSB
        else
            cont = false
        end
        nw += write(io, UInt8(byte))
    end
    nw
end

# max number of 7bit blocks for reading n bytes
# d,r = divrem(sizeof(T)*8, 7)
# (r > 0) && (d += 1)
const _max_n = [2, 3, 4, 5, 6, 7, 8, 10]

function _read_uleb_base(io::IO, ::Type{T}) where T <: Integer
    res = zero(T)
    n = 0
    byte = UInt8(MSB)
    while (byte & MSB) != 0
        byte = read(io, UInt8)
        res |= (convert(T, byte & MASK7) << (7*n))
        n += 1
    end
    n, res
end

function _read_uleb(io::IO, ::Type{Int32})
    n, res = _read_uleb_base(io, Int32)

    # negative int32 are encoded in 10 bytes (ref: https://developers.google.com/protocol-buffers/docs/encoding)
    # > if you use int32 or int64 as the type for a negative number, the resulting varint is always ten bytes long
    #
    # but Julia can be tolerant like the C protobuf implementation (unlike python)

    if n > _max_n[sizeof(res < 0 ? Int64 : Int32)]
        @debug("overflow reading $T. returning 0")
        return Int32(0)
    end

    res
end

function _read_uleb(io::IO, ::Type{T}) where T <: Integer
    n, res = _read_uleb_base(io, T)
    # in case of overflow, consider it as missing field and return default value
    if n > _max_n[sizeof(T)]
        @debug("overflow reading $T. returning 0")
        return zero(T)
    end
    res
end

function _write_zigzag(io::IO, x::T) where T <: Integer
    nbits = 8*sizeof(x)
    zx = (x << 1) ⊻ (x >> (nbits-1))
    _write_uleb(io, zx)
end

function _read_zigzag(io::IO, ::Type{T}) where T <: Integer
    zx = _read_uleb(io, UInt64)
    # result is positive if zx is even
    convert(T, iseven(zx) ? (zx >>> 1) : -signed((zx+1) >>> 1))
end

##
# read and write field keys
_write_key(io::IO, fld::Int, wiretyp::Int) = _write_uleb(io, (fld << 3) | wiretyp)
function _read_key(io::IO)
    key = _read_uleb(io, UInt64)
    wiretyp = key & MASK3
    fld = key >>> 3
    (fld, wiretyp)
end

##
# read and write field values
write_varint(io::IO, x::T) where {T <: Integer} = _write_uleb(io, x)
write_varint(io::IO, x::Int32) = _write_uleb(io, x < 0 ? Int64(x) : x)
write_bool(io::IO, x::Bool) = _write_uleb(io, x ? 1 : 0)
write_svarint(io::IO, x::T) where {T <: Integer} = _write_zigzag(io, x)

read_varint(io::IO, ::Type{T}) where {T <: Integer} = _read_uleb(io, T)
read_bool(io::IO) = Bool(_read_uleb(io, UInt64))
read_bool(io::IO, ::Type{Bool}) = read_bool(io)
read_svarint(io::IO, ::Type{T}) where {T <: Integer} = _read_zigzag(io, T)

write_fixed(io::IO, x::UInt32) = _write_fixed(io, x)
write_fixed(io::IO, x::UInt64) = _write_fixed(io, x)
write_fixed(io::IO, x::Int32) = _write_fixed(io, reinterpret(UInt32, x))
write_fixed(io::IO, x::Int64) = _write_fixed(io, reinterpret(UInt64, x))
write_fixed(io::IO, x::Float32) = _write_fixed(io, reinterpret(UInt32, x))
write_fixed(io::IO, x::Float64) = _write_fixed(io, reinterpret(UInt64, x))
function _write_fixed(io::IO, ux::T) where T <: Unsigned
    N = sizeof(ux)
    for n in 1:N
        write(io, UInt8(ux & MASK8))
        ux >>>= 8
    end
    N
end

read_fixed(io::IO, typ::Type{UInt32}) = _read_fixed(io, convert(UInt32,0), 4)
read_fixed(io::IO, typ::Type{UInt64}) = _read_fixed(io, convert(UInt64,0), 8)
read_fixed(io::IO, typ::Type{Int32}) = reinterpret(Int32, _read_fixed(io, convert(UInt32,0), 4))
read_fixed(io::IO, typ::Type{Int64}) = reinterpret(Int64, _read_fixed(io, convert(UInt64,0), 8))
read_fixed(io::IO, typ::Type{Float32}) = reinterpret(Float32, _read_fixed(io, convert(UInt32,0), 4))
read_fixed(io::IO, typ::Type{Float64}) = reinterpret(Float64, _read_fixed(io, convert(UInt64,0), 8))
function _read_fixed(io::IO, ret::T, N::Int) where T <: Unsigned
    for n in 0:(N-1)
        byte = convert(T, read(io, UInt8))
        ret |= (byte << (8*n))
    end
    ret
end

function write_bytes(io::IO, data::Vector{UInt8})
    n = _write_uleb(io, sizeof(data))
    n += write(io, data)
    n
end

function read_bytes(io::IO)
    n = _read_uleb(io, UInt64)
    data = Vector{UInt8}(undef, n)
    read!(io, data)
    data
end
read_bytes(io::IO, ::Type{Vector{UInt8}}) = read_bytes(io)

write_string(io::IO, x::AbstractString) = write_string(io, String(x))
write_string(io::IO, x::String) = write_bytes(io, @static isdefined(Base, :codeunits) ? unsafe_wrap(Vector{UInt8}, x) : Vector{UInt8}(x))

read_string(io::IO) = String(read_bytes(io))
read_string(io::IO, ::Type{T}) where {T <: AbstractString} = convert(T, read_string(io))

##
# read and write protobuf structures

mutable struct ProtoMetaAttribs
    fldnum::Int             # the field number in the structure
    fld::Symbol
    ptyp::Symbol            # protobuf type
    occurrence::Int         # 0: optional, 1: required, 2: repeated
    packed::Bool            # if repeated, whether packed
    default::Array          # the default value, empty array if none is specified, first element is used if something is specified
    meta::Any               # the ProtoMeta if this is a nested type
end

mutable struct ProtoMeta
    jtype::Type
    symdict::Dict{Symbol,ProtoMetaAttribs}
    numdict::Dict{Int,ProtoMetaAttribs}
    ordered::Vector{ProtoMetaAttribs}
    oneofs::Vector{Int}
    oneof_names::Vector{Symbol}

    ProtoMeta(jtype::Type, ordered::Vector{ProtoMetaAttribs}, oneofs::Vector{Int}=Int[], oneof_names::Vector{Symbol}=Symbol[]) = _setmeta(new(), jtype, ordered, oneofs, oneof_names)
end

function _setmeta(meta::ProtoMeta, jtype::Type, ordered::Vector{ProtoMetaAttribs}, oneofs::Vector{Int}, oneof_names::Vector{Symbol})
    symdict = Dict{Symbol,ProtoMetaAttribs}()
    numdict = Dict{Int,ProtoMetaAttribs}()
    for attrib in ordered
        symdict[attrib.fld] = numdict[attrib.fldnum] = attrib
    end
    meta.jtype = jtype
    meta.symdict = symdict
    meta.numdict = numdict
    meta.ordered = ordered
    meta.oneofs = oneofs
    meta.oneof_names = oneof_names
    meta
end

function write_map(io::IO, fldnum::Int, dict::Dict)
    dmeta = mapentry_meta(typeof(dict))
    iob = IOBuffer()

    n = 0
    for key in keys(dict)
        @debug("write_map writing key: $key")
        val = dict[key]
        writeproto(iob, key, dmeta.ordered[1])
        @debug("write_map writing val: $val")
        writeproto(iob, val, dmeta.ordered[2])
        n += _write_key(io, fldnum, WIRETYP_LENDELIM)
        n += write_bytes(io, take!(iob))
    end
    n
end

function writeproto(io::IO, val, attrib::ProtoMetaAttribs)
    fld = attrib.fldnum
    meta = attrib.meta
    ptyp = attrib.ptyp
    wiretyp, write_fn, read_fn, jtyp = WIRETYPES[ptyp]
    iob = IOBuffer()
    wfn = eval(write_fn)

    n = 0
    wfn(iob, convert(jtyp, val), meta)
    n += _write_key(io, fld, WIRETYP_LENDELIM)
    n += write_bytes(io, take!(iob))
    n
end

function writeproto(io::IO, val::T, attrib::ProtoMetaAttribs) where T<:Number
    fld = attrib.fldnum
    ptyp = attrib.ptyp
    wiretyp, write_fn, read_fn, jtyp = WIRETYPES[ptyp]
    wfn = eval(write_fn)

    n = 0
    n += _write_key(io, fld, wiretyp)
    n += wfn(io, convert(jtyp, val))
    n
end

function writeproto(io::IO, val::T, attrib::ProtoMetaAttribs) where T<:AbstractString
    fld = attrib.fldnum
    ptyp = attrib.ptyp
    wiretyp, write_fn, read_fn, jtyp = WIRETYPES[ptyp]

    n = 0
    n += _write_key(io, fld, wiretyp)
    n += write_string(io, val)
    n
end

writeproto(io::IO, val::Dict, attrib::ProtoMetaAttribs) = write_map(io, attrib.fldnum, convert(attrib.meta.jtype, val))

function writeproto(io::IO, val::Vector{UInt8}, attrib::ProtoMetaAttribs)
    fld = attrib.fldnum
    ptyp = attrib.ptyp
    wiretyp, write_fn, read_fn, jtyp = WIRETYPES[ptyp]

    n = 0
    n += _write_key(io, fld, wiretyp)
    n += write_bytes(io, val)
    n
end

function writeproto(io::IO, val::Array{T}, attrib::ProtoMetaAttribs) where T
    ptyp = attrib.ptyp
    wiretyp, write_fn, read_fn, jtyp = WIRETYPES[ptyp]
    wfn = eval(write_fn)

    writeproto(io, val, attrib, wfn)
end

function writeproto(io::IO, val::Array{T}, attrib::ProtoMetaAttribs, wfn::F) where {T,F}
    fld = attrib.fldnum
    meta = attrib.meta
    ptyp = attrib.ptyp
    wiretyp, write_fn, read_fn, jtyp = WIRETYPES[ptyp]
    iob = IOBuffer()

    n = 0
    (attrib.occurrence == 2) || error("expected meta attributes of $(attrib.fld) to specify an array")
    if attrib.packed
        # write elements as a length delimited field
        if ptyp == :obj
            error("can not write object field $fld as packed")
        else
            for eachval in val
                wfn(iob, convert(jtyp, eachval))
            end
        end
        n += _write_key(io, fld, WIRETYP_LENDELIM)
        n += write_bytes(io, take!(iob))
    else
        # write each element separately
        # maps can not be repeated
        if ptyp == :obj
            for eachval in val
                wfn(iob, convert(jtyp, eachval), meta)
                n += _write_key(io, fld, WIRETYP_LENDELIM)
                n += write_bytes(io, take!(iob))
            end
        else
            for eachval in val
                n += _write_key(io, fld, wiretyp)
                n += wfn(io, convert(jtyp, eachval))
            end
        end
    end
    n
end

function writeproto(io::IO, obj, meta::ProtoMeta=meta(typeof(obj)))
    n = 0
    @debug("writeproto writing an obj with meta: $meta")
    for attrib in meta.ordered
        fld = attrib.fld
        if isfilled(obj, fld)
            @debug("writeproto writing field: $fld")
            n += writeproto(io, getfield(obj, fld), attrib)
        else
            @debug("field not set: $fld")
            (attrib.occurrence == 1) && error("missing required field $fld (#$(attrib.fldnum))")
        end
    end
    n
end

function read_lendelim_packed(io, fld, reader, jtyp::Type)
    iob = IOBuffer(read_bytes(io))
    while !eof(iob)
        val = reader(iob, jtyp)
        push!(fld, val)
    end
    nothing
end

function read_lendelim_obj(io, val, meta::ProtoMeta, reader)
    fld_buf = read_bytes(io)
    reader(IOBuffer(fld_buf), val, meta)
    val
end

instantiate(t::Type) = ccall(:jl_new_struct_uninit, Any, (Any,), t)
instantiate(t::T) where {T <: ProtoType} = T()

function skip_field(io::IO, wiretype::Integer)
    if wiretype == WIRETYP_LENDELIM
        read_bytes(io)
    elseif wiretype == WIRETYP_64BIT
        read_fixed(io, UInt64)
    elseif wiretype == WIRETYP_32BIT
        read_fixed(io, UInt32)
    elseif wiretype == WIRETYP_VARINT
        read_varint(io, UInt64)
    end
    nothing
end

function read_map(io, dict::Dict{K,V}) where {K,V}
    iob = IOBuffer(read_bytes(io))

    dmeta = mapentry_meta(Dict{K,V})
    key_wtyp, key_wfn, key_rfn, key_jtyp = WIRETYPES[dmeta.numdict[1].ptyp]
    val_wtyp, val_wfn, val_rfn, val_jtyp = WIRETYPES[dmeta.numdict[2].ptyp]
    key_val = Vector{Union{K,V}}(undef, 2)

    while !eof(iob)
        fldnum, wiretyp = _read_key(iob)
        @debug("reading map fldnum: $fldnum")

        fldnum = Int(fldnum)
        attrib = dmeta.numdict[fldnum]

        if fldnum == 1
            key_val[1] = read_field(iob, nothing, attrib, key_wtyp, K)
        elseif fldnum == 2
            key_val[2] = read_field(iob, nothing, attrib, val_wtyp, V)
        else
            skip_field(iob, wiretyp)
        end
    end
    @debug("read map key: $(key_val[1])=$(key_val[2])")
    dict[key_val[1]] = key_val[2]
    dict
end

function read_field(io, container, attrib::ProtoMetaAttribs, wiretyp, jtyp_specific)
    ptyp = attrib.ptyp
    fld = attrib.fld

    _wiretyp, write_fn, read_fn, jtyp = WIRETYPES[ptyp]
    rfn = eval(read_fn)
    isrepeat = (attrib.occurrence == 2)

    if jtyp_specific != nothing
        jtyp = jtyp_specific
    elseif ptyp == :obj
        jtyp = attrib.meta.jtype
    elseif ptyp == :map
        jtyp = attrib.meta.jtype
    end

    if isrepeat
        arr_val = ((container != nothing) && isdefined(container, fld)) ? convert(Vector{jtyp}, getfield(container, fld)) : jtyp[]
        # Readers should accept repeated fields in both packed and expanded form.
        # Allows compatibility with old writers when [packed = true] is added later.
        # Only repeated fields of primitive numeric types (isbitstype == true) can be declared "packed".
        # Maps can not be repeated
        if isbitstype(jtyp) && (wiretyp == WIRETYP_LENDELIM)
            read_lendelim_packed(io, arr_val, rfn, jtyp)
        elseif ptyp == :obj
            push!(arr_val, read_lendelim_obj(io, instantiate(jtyp), attrib.meta, rfn))
        else
            push!(arr_val, rfn(io, jtyp))
        end
        return arr_val
    else
        (wiretyp != _wiretyp) && !isrepeat && error("cannot read wire type $wiretyp as $ptyp")

        if ptyp == :obj
            val_obj = ((container != nothing) && isdefined(container, fld)) ? getfield(container, fld) : instantiate(jtyp)
            return read_lendelim_obj(io, val_obj, attrib.meta, rfn)
        elseif ptyp == :map
            val_map = ((container != nothing) && isdefined(container, fld)) ? convert(jtyp, getfield(container, fld)) : jtyp()
            return read_map(io, val_map)
        else
            @debug("reading type $jtyp")
            return rfn(io, jtyp)
        end
    end
end

function readproto(io::IO, obj, meta::ProtoMeta=meta(typeof(obj)))
    @debug("readproto begin: $(typeof(obj))")
    fillunset(obj)
    fldnums = collect(keys(meta.numdict))
    while !eof(io)
        fldnum, wiretyp = _read_key(io)
        @debug("reading fldnum: $(typeof(obj)).$fldnum")

        fldnum = Int(fldnum)
        # ignore unknown fields
        if !(fldnum in fldnums)
            @debug("skipping unknown field: $(typeof(obj)).$fldnum")
            skip_field(io, wiretyp)
            continue
        end
        attrib = meta.numdict[fldnum]

        val = read_field(io, obj, attrib, wiretyp, nothing)
        fld = attrib.fld
        fillset(obj, fld)
        if (attrib.occurrence == 2) || (attrib.ptyp == :obj) || (attrib.ptyp == :map)
            setfield!(obj, fld, convert(fld_type(obj, fld), val))
        else
            setfield!(obj, fld, val)
        end
    end

    # populate defaults
    fnames = fld_names(typeof(obj))
    fill = filled(obj)
    for attrib in meta.ordered
        fld = attrib.fld
        idx = something(findfirst(isequal(fld), fnames))
        # TODO: do not fill if oneof the fields in the oneof
        if !isfilled(obj, fld) && (length(attrib.default) > 0) && !_isset_oneof(fill, meta.oneofs, idx)
            default = attrib.default[1]
            setfield!(obj, fld, convert(fld_type(obj, fld), deepcopy(default)))
            @debug("readproto set default: $(typeof(obj)).$fld = $default")
            fillset_default(obj, fld)
        end
    end
    @debug("readproto end: $(typeof(obj))")
    obj
end

##
# helpers
oiddict() = @static isdefined(Base, :IdDict) ? IdDict() : ObjectIdDict()
const _metacache = oiddict() # dict of Type => ProtoMeta
const _mapentry_metacache = oiddict()
const _fillcache = Dict{UInt,BitArray{2}}()

const DEF_REQ = Symbol[]
const DEF_FNUM = Int[]
const DEF_VAL = Dict{Symbol,Any}()
const DEF_PACK = Symbol[]
const DEF_WTYPES = Dict{Symbol,Symbol}()
const DEF_ONEOFS = Int[]
const DEF_ONEOF_NAMES = Symbol[]
const DEF_FIELD_TYPES = Dict{Symbol,String}()

meta(typ::Type) = haskey(_metacache, typ) ? _metacache[typ] : meta(typ, DEF_REQ, DEF_FNUM, DEF_VAL, true, DEF_PACK, DEF_WTYPES, DEF_ONEOFS, DEF_ONEOF_NAMES, DEF_FIELD_TYPES)
function meta(typ::Type, required::Array, numbers::Array, defaults::Dict, cache::Bool=true, pack::Array=DEF_PACK, wtypes::Dict=DEF_WTYPES,
                oneofs::Vector{Int}=DEF_ONEOFS, oneof_names::Vector{Symbol}=DEF_ONEOF_NAMES, field_types::Dict{Symbol,String}=DEF_FIELD_TYPES)
    haskey(_metacache, typ) && return _metacache[typ]
    d = Dict{Symbol,Any}()
    for (k,v) in defaults
        d[k] = v
    end
    meta(typ, convert(Vector{Symbol}, required), convert(Vector{Int}, numbers), d, cache, convert(Vector{Symbol}, pack), wtypes, oneofs, oneof_names, field_types)
end
function meta(typ::Type, required::Vector{Symbol}, numbers::Vector{Int}, defaults::Dict{Symbol,Any}, cache::Bool=true, pack::Vector{Symbol}=DEF_PACK,
                wtypes::Dict=DEF_WTYPES, oneofs::Vector{Int}=DEF_ONEOFS, oneof_names::Vector{Symbol}=DEF_ONEOF_NAMES, field_types::Dict{Symbol,String}=DEF_FIELD_TYPES)
    haskey(_metacache, typ) && return _metacache[typ]

    m = ProtoMeta(typ, ProtoMetaAttribs[])
    cache && (_metacache[typ] = m)

    attribs = ProtoMetaAttribs[]
    names = fld_names(typ)
    types = typ.types
    for fldidx in 1:length(names)
        fldname = names[fldidx]
        fldtyp = (fldname in keys(field_types)) ? Core.eval(typ.name.module, (@static (VERSION < v"0.7.0-alpha") ? parse : Meta.parse)(field_types[fldname])) : types[fldidx]
        fldnum = isempty(numbers) ? fldidx : numbers[fldidx]
        isarr = (fldtyp <: Array) && !(fldtyp === Vector{UInt8})
        repeat = isarr ? 2 : (fldname in required) ? 1 : 0

        elemtyp = isarr ? eltype(fldtyp) : fldtyp
        wtyp = get(wtypes, fldname, wiretype(elemtyp))
        packed = (isarr && (fldname in pack))
        default = haskey(defaults, fldname) ? Any[defaults[fldname]] : defaultval(fldtyp)

        fldmeta = (wtyp == :obj) ? meta(elemtyp) :
                  (wtyp == :map) ? mapentry_meta(elemtyp) : nothing
        push!(attribs, ProtoMetaAttribs(fldnum, fldname, wtyp, repeat, packed, default, fldmeta))
    end
    _setmeta(m, typ, attribs, oneofs, oneof_names)
    m
end

function mapentry_meta(typ::Type{Dict{K,V}}) where {K,V}
    m = ProtoMeta(typ, ProtoMetaAttribs[])
    _mapentry_metacache[typ] = m

    attribs = ProtoMetaAttribs[]
    push!(attribs, ProtoMetaAttribs(1, :key, wiretype(K), 0, false, defaultval(K), nothing))

    isarr = (V <: Array) && !(V === Vector{UInt8})
    repeat = isarr ? 2 : 0
    packed = isarr
    wtyp = wiretype(V)
    vmeta = (wtyp == :obj) ? meta(V) :
            (wtyp == :map) ? mapentry_meta(V) : nothing
    push!(attribs, ProtoMetaAttribs(2, :value, wtyp, repeat, packed, defaultval(V), vmeta))

    _setmeta(m, typ, attribs, DEF_ONEOFS, DEF_ONEOF_NAMES)
    m
end

function _unset_oneofs(fill::BitArray{2}, oneofs::Vector{Int}, idx::Int)
    oneofidx = isempty(oneofs) ? 0 : oneofs[idx]
    if oneofidx > 0
        # unset other fields in the oneof group
        for uidx = 1:length(oneofs)
            if (oneofs[uidx] == oneofidx) && (uidx !== idx)
                fill[1:2,uidx] .= false
            end
        end
    end
end

function _isset_oneof(fill::BitArray{2}, oneofs::Vector{Int}, idx::Int)
    oneofidx = isempty(oneofs) ? 0 : oneofs[idx]
    if oneofidx > 0
        # find if any field in the oneof group is set
        for uidx = 1:length(oneofs)
            if oneofs[uidx] == oneofidx
                fill[1,uidx] && (return true)
            end
        end
    end
    false
end

fillunset(obj) = (fill!(filled(obj), false); nothing)
fillunset(obj, fld::Symbol) = _fillset(obj, fld, false, false)
fillset(obj, fld::Symbol) = _fillset(obj, fld, true, false)
fillset_default(obj, fld::Symbol) = _fillset(obj, fld, true, true)
function _fillset(obj, fld::Symbol, val::Bool, isdefault::Bool)
    fill = filled(obj)
    fnames = fld_names(typeof(obj))
    idx = something(findfirst(isequal(fld), fnames))
    fill[1,idx] = val
    (!val || isdefault) && (fill[2,idx] = val)
    val && _unset_oneofs(fill, meta(typeof(obj)).oneofs, idx)
    nothing
end

function filled(obj)
    oid = objectid(obj)
    haskey(_fillcache, oid) && return _fillcache[oid]

    fnames = fld_names(typeof(obj))
    fill = fill!(BitArray(undef, 2, length(fnames)), false)
    oneofs = meta(typeof(obj)).oneofs
    for idx in 1:length(fnames)
        if isdefined(obj, fnames[idx])
            fill[1,idx] = true
            _unset_oneofs(fill, oneofs, idx)
        end
    end
    if !isimmutable(obj)
        _fillcache[oid] = fill
        finalizer(obj->delete!(_fillcache, objectid(obj)), obj)
    end
    fill
end

isfilled(obj, fld::Symbol) = _isfilled(obj, fld, false)
isfilled_default(obj, fld::Symbol) = _isfilled(obj, fld, true)
function _isfilled(obj, fld::Symbol, isdefault::Bool)
    fnames = fld_names(typeof(obj))
    idx = something(findfirst(isequal(fld), fnames))
    filled(obj)[isdefault ? 2 : 1, idx]
end

function isfilled(obj)
    fill = filled(obj)
    flds = meta(typeof(obj)).ordered
    for idx in 1:length(flds)
        fld = flds[idx]
        if fld.occurrence == 1
            fill[1,idx] || (return false)
            (fld.meta != nothing) && !isfilled(getfield(obj, fld.fld)) && (return false)
        end
    end
    true
end

function which_oneof(obj, oneof::Symbol)
    m = meta(typeof(obj))
    fill = filled(obj)
    fnames = fld_names(typeof(obj))
    oneofs = m.oneofs
    oneof_idx = something(findfirst(isequal(oneof), m.oneof_names))

    for idx in 1:length(oneofs)
        (oneofs[idx] == oneof_idx) && fill[1,idx] && (return fnames[idx])
    end
    Symbol()
end

function show(io::IO, m::ProtoMeta)
    println(io, "ProtoMeta for $(m.jtype)")
    println(io, m.ordered)
end

##
# Enum Lookup

abstract type ProtoEnum end

function lookup(en::ProtoEnum, val)
    for name in fld_names(typeof(en))
        (val == getfield(en, name)) && return name
    end
    error("Enum $(typeof(en)) has no value: $val")
end


text:
using FastGaussQuadrature
include("chebyshev.jl")
include("pgrid.jl")

idx(ki, xi, order) = (ki - 1) * order + xi

struct CompositeGrid
    grid::Vector{Float64} # panel x finegrid
    wgrid::Vector{Float64}
    panel::Vector{Float64}
    x::Vector{Float64}
    w::Vector{Float64}
    # xp::Matrix{Float64} # the fine grid
    # wp::Matrix{Float64} # the weight of the fine grid
    Np::Int  # number of panel - 1
    order::Int  # number of fine grid
    type::Symbol
    function CompositeGrid(panel, order, type::Symbol)
        if type == :cheb
            x, w = barychebinit(order)
        elseif type == :gaussian
            x, w = gausslegendre(order)
        else
            throw("not implemented!")
        end

        # println("Quadrature type : ", type)
        # println("Quadrature points : ", x)
        # println("Quadrature weights: ", w)

        Np = length(panel) - 1

        grid = zeros(Float64, Np * order)
        wgrid = zeros(Float64, Np * order)
        # xp = zeros(Float64, (Np, order))
        # wp = zeros(Float64, (Np, order))

        for p in 1:Np
            a, b = panel[p], panel[p + 1]
            for xi in 1:order
                index = idx(p, xi, order)
                grid[index] = (a + b) / 2 + (b - a) / 2 * x[xi]
                wgrid[index] = (b - a) / 2 * w[xi]
            end
        end

        return new(grid, wgrid, panel, x, w, Np, order, type)
    end
end

function KPanel(Nk, kF, maxK, minK)
    panel = Grid.boseKUL(0.5 * kF, maxK, minK, Nk, 1).grid
    panel[1] = 0.0  # the kgrid start with 0.0
    return panel
end

function QPanel(Nk, kF, maxK, minK, k)
    panel = pGrid(k, kF, maxK, minK, [Nk,Nk,Nk,Nk]).grid

    # panel = Grid.boseKUL(0.5 * kF, maxK, minK, Nk, 1).grid
    # panel[1] = 0.0  # the kgrid start with 0.0

    return panel
end

"""
interpolate!(f, k::CompositeGrid, grid)

map f array in the grid k to a new grid

# Arguments
- f::Vector{Float64}: vector of data.
- k::CompositeGrid: the grid object that f is defined on
- grid::Vector{Float64}: new grid points
"""
function interpolate(f, k::CompositeGrid, grid)
    @assert k.type == :cheb
    order = k.order
    ff = zeros(eltype(f), length(grid))
    kpidx = 1 # panel index of the kgrid
    head, tail = idx(kpidx, 1, order), idx(kpidx, order, order)
    # extract all x in the kpidx-th k panel
    fx = @view f[head:tail] # all F in the same kpidx-th K panel
    x = @view k.grid[head:tail]
    w = @view k.wgrid[head:tail]

    for (qi, q) in enumerate(grid)
        # for a given q, one needs to find the k panel to do interpolation
        if q > k.panel[kpidx + 1]
            # if q is too large, move k panel to the next
            # println("before $q, $kpidx, $(k.panel[kpidx]) -> $(k.panel[kpidx + 1])")
            while q > kgrid.panel[kpidx + 1]
                kpidx += 1
            end
            # println("after $q, $kpidx, $(k.panel[kpidx]) -> $(k.panel[kpidx + 1])")
            head, tail = idx(kpidx, 1, order), idx(kpidx, order, order)
            fx = @view f[head:tail] # all F in the same kpidx-th K panel
            x = @view k.grid[head:tail]
            w = @view k.wgrid[head:tail]
            @assert kpidx <= k.Np
        end
        ff[qi] = barycheb(order, q, fx, w, x) # the interpolation is independent with the panel length
        # @assert k.panel[kpidx - 1][order] <= q <= k.panel[pidx + 1][1] "$q for kpidx=$kpidx with $x"
    end
    return ff
end

if abspath(PROGRAM_FILE) == @__FILE__
    grid=CompositeGrid([0.0, 1.0, 2.0], 4, :cheb)
    println(grid.grid)
end

ext: jl


text:
push!(LOAD_PATH, ".")

module RealIndividual

include("utility.jl")

using Individual
using Parallel

abstract type AbstractRealIndividual{GeneType <: Real} <: Individual.AbstractIndividual{GeneType} end

@define RealIndividual begin
  RealIndividual.Individual.@Individual
end

struct RealIndividualType{GeneType} <: AbstractRealIndividual{GeneType}
  @RealIndividual

  function RealIndividualType{GeneType}(x...) where {GeneType}
    args = build(x...)
    new(args[1], args[2], args[3])
  end
end

function build(x...)
  return Individual.build(x...)
end

function Individual.generateRandom!(this::AbstractRealIndividual)
  for i in eachindex(this)
    gene_lb, gene_ub = this.bounds[i]
    this[i] = gene_lb + Parallel.threadRand() * (gene_ub - gene_lb)
  end
end

end

max_stars_repo_name: Markhyz/Julia-Genetic-Algorithm

max_issues_repo_name: Markhyz/Julia-Genetic-Algorithm

max_forks_repo_name: Markhyz/Julia-Genetic-Algorithm


text:
@testset "ToyModel" begin

println("Testing TB Toy Model")

TB = TightBinding
at = (1,2,3) * bulk(:Al, pbc=false, cubic=true)
tbm = TightBinding.ToyTBModel(r0=2.5, rcut=4.9)

X = positions(at)
X[2] += [0.123, 0.234, 0.01]
set_positions!(at, X)
@show length(at)

println("-------------------------------------------")
println("Testing ToyTBModel: ")
@show length(at)
print("check that hamiltonian evaluates ... ")
H, M = hamiltonian(tbm, at)
println("ok.")
print("check that `densitymatrix` evaluates ... ")
Γ = densitymatrix(tbm, at)
println("ok.")
print("check that `energy` evaluates ... ")
E = energy(tbm, at)
println("ok : E = ", E, ".")
print("check that `forces` evaluates ... ")
frc = forces(tbm, at)
println("ok : |f|∞ = ", maximum(norm.(frc)), ".")
@test true

println("-------------------------------------------")
println("  Finite-difference test with ToyTBModel:  ")
println("-------------------------------------------")
@test fdtest(tbm, at, verbose=true)
println("-------------------------------------------")

end


text:
# source code representation cannot be recovered, but code_lowered yes
# see for example Symbolics.jl or Cassette.jl

# unary minus, plus

# some unary functions:
# sqrt, cbrt, exp, log, log2, log10

# some trig functions:
# sin    cos    tan    cot    sec    csc
# sinh   cosh   tanh   coth   sech   csch
# asin   acos   atan   acot   asec   acsc
# asinh  acosh  atanh  acoth  asech  acsch
# sinc   cosc

# Implement:

# unary minus, plus
# sqrt, cbrt, exp, log, log2, log10
# sin    cos    tan

# what not easy?
# simplify!
# rewriting rules
# trivial cases:
#   if x - x -> 0
#   if x / x -> 1
#   if x*x -> x^2
#   if x^i / x^j -> x^(i-j)
#   if x^0 -> 1
#   math operations on constants should be evaluated
#   if c1 + c2 -> c3(c1.value + c2.value)

# implement constants

module symbolics

function default_variable_value()
    return 0.0
end

abstract type AbstractVariable end
abstract type AbstractSymbolicExpression end
abstract type AbstractConstant end

# ============================================================================ #
# Constant
# ============================================================================ #

struct Constant <: AbstractConstant
    value::Float64

    function Constant(value)
        return new( value )
    end
end

function evaluate(x::Constant)
    return x.value
end

function to_string(x::Constant)
    return string(x.value)
end

# ============================================================================ #
# Variable
# ============================================================================ #

mutable struct Variable <: AbstractVariable
    name::String
    value::Float64

    function Variable(name, value = 0.0)
        return new(string(name), value)
    end
end

macro variable(varname, value = default_variable_value())
    lhs = esc(varname)
    rhs = :(     Variable( $(QuoteNode(varname)), $value)   )
    ex  = :($(lhs) = $(rhs))

    return ex
end

function evaluate(x::Variable)
    return x.value
end

function to_string(x::Variable)
    return x.name
end

# ============================================================================ #
# SymbolicExpression
# ============================================================================ #

mutable struct SymbolicExpression <: AbstractSymbolicExpression
    value::Float64
    head::Symbol
    args::Vector{Any}

    function SymbolicExpression(head, args...)
        return new(default_variable_value(), head, collect(args)) # [args...])
    end
end

for op in (:+, :-, :/, :*, :^)
    @eval function Base.$(op)(
        x::Tx, y::Ty) where {
            Tx <: Union{Constant, Variable, SymbolicExpression},
            Ty <: Union{Constant, Variable, SymbolicExpression}}

        ex = SymbolicExpression(Symbol(Base.$op), x, y)
        return ex
    end

    @eval function Base.$(op)(
        x::Tx, y::Ty) where {
            Tx <: Number,
            Ty <: Union{Constant, Variable, SymbolicExpression}}

        c = Constant(x)
        ex = SymbolicExpression(Symbol(Base.$op), c, y)
        return ex
    end

    @eval function Base.$(op)(
        x::Tx, y::Ty) where {
            Tx <: Union{Constant, Variable, SymbolicExpression},
            Ty <: Number}

        c = Constant(y)
        ex = SymbolicExpression(Symbol(Base.$op), x, c)
        return ex
    end
end

function Base.:(+)(
    x::Tx) where {Tx <: Union{Constant, Variable, SymbolicExpression}}
    return x
end

function Base.:(-)(
    x::Tx) where {Tx <: Union{Constant, Variable, SymbolicExpression}}
    ex = SymbolicExpression(:(-), x)
    return ex
end

for op in (:sqrt, :cbrt, :exp, :log, :log2, :log10, :sin, :cos, :tan)
    @eval function Base.$(op)(
        x::Tx) where {Tx <: Union{Constant, Variable, SymbolicExpression}}

        ex = SymbolicExpression(Symbol(Base.$op), x)
        return ex
    end
end

function evaluate(ex::SymbolicExpression)

    values = Vector{Any}(undef, length(ex.args))
    for (ii, arg) in enumerate(ex.args)
        values[ii] = evaluate(arg)
    end

    ex2call  = Expr(:call, ex.head, values...)
    val      = eval(ex2call)
    ex.value = val
    return val
end

function substitute!(
    ex::SymbolicExpression,
    old_var::T,
    new_ex::U) where {
        T  <: Union{Constant, Variable, SymbolicExpression},
        U <: Union{Constant, Variable, SymbolicExpression}}

    for (ii, ee) in enumerate(ex.args)
        if ee isa Variable || ee isa Constant
            if ee == old_var
                ex.args[ii] = new_ex
            end
        elseif ee isa SymbolicExpression
            substitute!(ee, old_var, new_ex)
        else
            error("something went wrong")
        end
    end

    return ex
end

# ============================================================================ #
# Operations
# ============================================================================ #

function is_constant(x)
    return false
end
function is_constant(x::Constant)
    return true
end

function is_variable(x)
    return false
end
function is_variable(x::Variable)
    return true
end

function is_expression(x)
    return false
end
function is_expression(x::SymbolicExpression)
    return true
end

function is_composite_expression(x::T) where {T <: Union{Constant, Variable}}
    return false
end
function is_composite_expression(ex::SymbolicExpression)
    return true
end

function is_commutative(ex::SymbolicExpression)
    return  is_plus_operator(ex.head) ||
            is_product_operator(ex.head)
end

function to_string(ex::SymbolicExpression)
# s expression style is easy, but not very nice
    if is_binary_operator(ex.head)
        __to_string_binary_op(ex)
    elseif is_unary_operator(ex.head)
        __to_string_unary_op(ex)
    else
        __show_k_ary_expression(ex)
    end
end

function __to_string_minus_op(ex::SymbolicExpression)
    total_str = "-" * ex.args[1].name
    return total_str
end

function __to_string_unary_op(ex::SymbolicExpression)
    if is_minus_operator(ex.head)
        return __to_string_minus_op(ex)
    end
    arg_str   = to_string(ex.args[1])
    total_str = string(ex.head) * "(" * arg_str * ")"
    return total_str
end

function __show_k_ary_expression(ex)
    op = string(ex.head)

end

function __to_string_binary_op(ex::SymbolicExpression)

    L = length(ex.args)
    if L < 2
       return __to_string_unary_op(ex)
    end

    op  = string(ex.head)
    lhs = ex.args[1]
    rhs = ex.args[2]

    # c op c
    # var op var
    # comp op comp -> hard case
    #
    # c op var | var op c
    #
    # c op comp

    # rule const op const
    if is_constant(ex.args[1]) && is_constant(ex.args[2])
        return __to_string_const_op_const(ex)
    end

    # rule: var op var
    if is_variable(ex.args[1]) && is_variable(ex.args[2])
        return __to_string_var_op_var(ex)
    end

    # rule: c op var | var op c
    if (is_variable(ex.args[1]) || is_constant(ex.args[1])) &&
        (is_variable(ex.args[2]) || is_constant(ex.args[2]))
        return __to_string_var_op_var(ex)
    end

    # c op comp
    if is_constant(ex.args[1]) && is_composite_expression(ex.args[2])
        return __to_string_const_op_comp(ex, true)
    end
    if is_constant(ex.args[2]) && is_composite_expression(ex.args[1])
        return __to_string_const_op_comp(ex, false)
    end

    # rule: constant */^ var or comp both way
    if (is_product_operator(ex.head) || is_divide_operator(ex.head) || is_power_operator(ex.head)) &&
        ( is_constant(ex.args[1]) || is_constant(ex.args[2]) ||
          is_variable(ex.args[1]) || is_variable(ex.args[2]) )
        return __to_string_constant_op_composite_expression(ex)
    end

    # rule: composite + composite
    if is_plus_operator(ex.head) & is_composite_expression(ex.args[1]) & is_composite_expression(ex.args[2])
        return __to_string_comp_plus_comp(ex)
    end

    # rule: composite - composite
    if is_minus_operator(ex.head) & is_composite_expression(ex.args[1]) & is_composite_expression(ex.args[2])
        return __to_string_comp_minus_comp(ex)
    end

    # rule: composite * composite
    if is_product_operator(ex.head) & is_composite_expression(ex.args[1]) & is_composite_expression(ex.args[2])
        return __to_string_comp_prod_comp(ex)
    end

    return __to_string_fallback(ex)
end

function __to_string_fallback(ex)

    op        = string(ex.head)
    lhs       = ex.args[1]
    rhs       = ex.args[2]
    lhs_str   = to_string(lhs)
    rhs_str   = to_string(rhs)
    ws        = " "

    left_paren   = "("
    right_paren  = ")"

    return left_paren * lhs_str * right_paren * ws * op * ws * left_paren * rhs_str * right_paren
end

function __to_string_const_op_const(ex)
    op        = string(ex.head)
    lhs       = ex.args[1]
    rhs       = ex.args[2]
    lhs_str   = to_string(lhs)
    rhs_str   = to_string(rhs)
    ws        = " "

    total_str = lhs_str * ws * op * ws * rhs_str

    return total_str
end

function __to_string_var_op_var(ex::SymbolicExpression)
    op        = string(ex.head)
    lhs       = ex.args[1]
    rhs       = ex.args[2]
    lhs_str   = to_string(lhs)
    rhs_str   = to_string(rhs)
    ws        = is_power_operator(ex.head) ? "" : " "
    total_str = lhs_str * ws * op * ws * rhs_str

    return total_str
end

function __to_string_const_op_comp(ex, lhs_is_const::Bool)

    op        = ex.head
    op_str    = string(op)
    lhs       = ex.args[1]
    rhs       = ex.args[2]
    lhs_str   = to_string(lhs)
    rhs_str   = to_string(rhs)
    ws        = " "
    left_paren   = "("
    right_paren  = ")"

    if is_plus_operator(op)
        # c + ex
        return lhs_str * ws * op_str * ws * rhs_str
    end
    if is_minus_operator(op)
        # c - ex
        # c - (a +- b)
        # paren if the next op is +-
        if is_plus_operator(rhs.head) || is_minus_operator(rhs.head)
            if lhs_is_const
                return lhs_str * ws * op_str * ws * left_paren * rhs_str * right_paren
            else
                return left_paren * lhs_str * right_paren * ws * op_str * ws * rhs_str
            end
        else
            return lhs_str * ws * op_str * ws * rhs_str
        end
    end
    if is_product_operator(op) || is_divide_operator(rhs.head)
        # dont paren if the next op is also * or /
        if is_product_operator(rhs.head) || is_divide_operator(rhs.head)
            return lhs_str * ws * op_str * ws * rhs_str
        else
            if lhs_is_const
                return lhs_str * ws * op_str * ws * left_paren * rhs_str * right_paren
            else
                return left_paren * lhs_str * right_paren * ws * op_str * ws * rhs_str
            end
        end
    end

end

function __to_string_constant_op_composite_expression(ex::SymbolicExpression)
    op  = string(ex.head)
    lhs = ex.args[1]
    rhs = ex.args[2]

    if is_constant(lhs) || is_variable(lhs)
        return __to_string_constant_op_composite_expression_constvarleft(ex)
    end

    if is_constant(rhs) || is_variable(rhs)
        return __to_string_constant_op_composite_expression_constvarright(ex)
    end

#=
    op  = string(ex.head)
    lhs = ex.args[1]
    rhs = ex.args[2]

    always_paren = is_power_operator(ex.head)
    rather_paren = is_composite_expression(lhs) || is_composite_expression(rhs)
    do_paren     = always_paren | rather_paren
    left_paren   = do_paren ? "(" : ""
    right_paren  = do_paren ? ")" : ""
    ws           = is_power_operator(ex.head) ? "" : " "

    lhs_str   = to_string(lhs)
    rhs_str   = to_string(rhs)

    if is_constant(lhs) || is_variable(lhs)
        total_str = lhs_str * ws * op * ws * left_paren * rhs_str * right_paren
    end

    if is_constant(rhs) || is_variable(rhs)
        total_str = left_paren * lhs_str * right_paren * ws * op * ws * rhs_str
    end

    return total_str
=#
end

function __to_string_constant_op_composite_expression_constvarleft(ex)

    # a * (b +- c)
    # a / (b +- c)
    # a ^ (b +-*/ c)
    #
    # a * b */ c
    # a / b */ c

    op        = ex.head
    op_str    = string(op)
    lhs       = ex.args[1]
    rhs       = ex.args[2]
    lhs_str   = to_string(lhs)
    rhs_str   = to_string(rhs)
    ws        = " "
    left_paren   = "("
    right_paren  = ")"

    b1 = is_plus_operator(rhs.head) || is_minus_operator(rhs.head)

    do_paren = false
    if (is_product_operator(op) ||  is_divide_operator(op)) && b1
        do_paren = true
    end

    omit_ws = false
    if is_power_operator(op)
        do_paren = true
        omit_ws  = true
    end

    ws          = omit_ws ? "" : ws
    left_paren  = do_paren ? left_paren : ""
    right_paren = do_paren ? right_paren : ""

    return lhs_str * ws * op_str * ws * left_paren * rhs_str * right_paren
end

function __to_string_constant_op_composite_expression_constvarright(ex)

    # (a +- b) * c
    # (a +- b) / c
    # (a +- b) ^ c
    #
    # a * b */ c

    op        = ex.head
    op_str    = string(op)
    lhs       = ex.args[1]
    rhs       = ex.args[2]
    lhs_str   = to_string(lhs)
    rhs_str   = to_string(rhs)
    ws        = " "
    left_paren   = "("
    right_paren  = ")"

    b1 = is_plus_operator(lhs.head) || is_minus_operator(lhs.head)

    do_paren = false
    if (is_product_operator(op) ||  is_divide_operator(op)) && b1
        do_paren = true
    end

    omit_ws = false
    if is_power_operator(op)
        do_paren = true
        omit_ws  = true
    end

    ws          = omit_ws ? "" : ws
    left_paren  = do_paren ? left_paren : ""
    right_paren = do_paren ? right_paren : ""

    return left_paren * lhs_str * right_paren * ws * op_str * ws * rhs_str
end

function __to_string_comp_plus_comp(ex)
    op        = string(ex.head)
    lhs       = ex.args[1]
    rhs       = ex.args[2]
    lhs_str   = to_string(lhs)
    rhs_str   = to_string(rhs)
    ws        = " "

    total_str = lhs_str * ws * op * ws * rhs_str

    return total_str
end

function __to_string_comp_minus_comp(ex)

    # a + b - (c +- d)
    # a + b - c */ d

    # if the rhs op is
    # +, - -> paren
    # *, / -> dont paren if the rhs.args are either const or var
    op        = string(ex.head)
    lhs       = ex.args[1]
    rhs       = ex.args[2]
    lhs_str   = to_string(lhs)
    rhs_str   = to_string(rhs)
    ws        = " "
    always_paren   = is_plus_operator(rhs.head) || is_minus_operator(rhs.head)
    possibly_paren = ( is_product_operator(rhs.head) ||
                       is_divide_operator(rhs.head) ) &&
                       !(is_constant(rhs.args[1]) || is_variable(rhs.args[1]))
    do_paren    = always_paren | possibly_paren
    left_paren  = do_paren ? "(" : ""
    right_paren = do_paren ? ")" : ""

    total_str = lhs_str * ws * op * ws * left_paren * rhs_str * right_paren

    return total_str
end

function __to_string_comp_prod_comp(ex)

    # special rule (not implemented yet)
    # if in one branch only * -> no paren is needed

    # (a +- b) * (c +- d)
    # a*/b * c*/d

    # if the rhs op is
    # +, - -> paren
    # *, / -> dont paren if the rhs.args are either const or var
    op        = string(ex.head)
    lhs       = ex.args[1]
    rhs       = ex.args[2]
    lhs_str   = to_string(lhs)
    rhs_str   = to_string(rhs)
    ws        = " "

    always_paren   = is_plus_operator(rhs.head) || is_minus_operator(rhs.head)
    possibly_paren = ( is_product_operator(rhs.head) ||
                       is_divide_operator(rhs.head) ) &&
                       !(is_constant(rhs.args[1]) || is_variable(rhs.args[1]))
    do_paren_right  = always_paren | possibly_paren

    always_paren   = is_plus_operator(lhs.head) || is_minus_operator(lhs.head)
    possibly_paren = ( is_product_operator(lhs.head) ||
                       is_divide_operator(lhs.head) ) &&
                       !(is_constant(lhs.args[1]) || is_variable(lhs.args[1]))
    do_paren_left  = always_paren | possibly_paren

    left_paren  = "("
    right_paren = ")"

    left_str  = do_paren_left ? (left_paren * lhs_str * right_paren) : lhs_str
    right_str = do_paren_right ? (left_paren * rhs_str * right_paren) : rhs_str
    total_str = left_str *  ws * op * ws * right_str

    return total_str
end

include("is_operators.jl")

end # module

module stest
using ..symbolics

symbolics.@variable x 10
symbolics.@variable y 0.2
symbolics.@variable z 3

e1 = x + y
e2 = x / y

e3 = z ^ (x + y)

e4 = -x
e5 = -x - y
e6 = sin(x + y - z * x)

# e3_val = symbolics.evaluate(e3)

# str1 = symbolics.to_string( stest.e1 )
# str2 = symbolics.to_string( stest.e2 )
# str3 = symbolics.to_string( stest.e3 )
# str4 = symbolics.to_string( stest.e4 )
# str5 = symbolics.to_string( stest.e5 )
# str6 = symbolics.to_string( stest.e6 )

# e7 = y * 10
# str7 = symbolics.to_string( stest.e7 )
# e8 = 10 * y
# str8 = symbolics.to_string( stest.e8 )

# e9 = 10 + (x - y)
# str9 = symbolics.to_string( stest.e9 )

# !!!
# this doesnt work correctly now
# e7 = y * 10 * (x * 2 - 3*z)
e7 = x * y * x * (x * 2 - 3*z)
str7 = symbolics.to_string( stest.e7 )

# str = symbolics.to_string( stest.e6 )

# x_val = symbolics.evaluate(x)
# y_val = symbolics.evaluate(x)

# e1_val = symbolics.evaluate(e1)

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
using Groups

@testset "sum of squares in FreeGroup *-algebra" begin
    StarAlgebras.star(g::Groups.GroupElement) = inv(g)
    F = Groups.FreeGroup(4)
    S = [Groups.gens(F); inv.(Groups.gens(F))]

    ID = one(F)
    RADIUS = 3
    @time E_R, sizes = Groups.wlmetric_ball(S, ID, radius = 2 * RADIUS)
    @test sizes == [9, 65, 457, 3201, 22409, 156865]

    b = StarAlgebras.Basis{UInt32}(E_R)

    tmstr = StarAlgebras.MTable{true}(b, table_size = (sizes[RADIUS], sizes[RADIUS]))

    RG = StarAlgebra(F, b, tmstr)

    g, h, k, l = S[1:4]

    length(b)

    G = (one(RG) - RG(g))
    @test G^2 == StarAlgebras.mul!(zero(G), G, G) == 2one(RG) - RG(g) - star(RG(g))
    @test star(G * G) == G * G

    @testset "Sums of hermitian squares" begin
        ∗ = StarAlgebras.star

        𝕀 = one(RG)

        G = (𝕀 - RG(g))
        H = (𝕀 - RG(h))
        K = (𝕀 - RG(k))
        L = (𝕀 - RG(l))
        GH = (𝕀 - RG(g * h))
        KL = (𝕀 - RG(k * l))

        X = (2𝕀 - ∗(RG(g)) - RG(h))
        Y = (2𝕀 - ∗(RG(g * h)) - RG(k))

        @test -(2𝕀 - RG(g * h) - ∗(RG(g * h))) + 2G^2 + 2H^2 == X^2
        @test (2𝕀 - RG(g * h) - ∗(RG(g * h))) == GH^2
        @test -(2𝕀 - RG(g * h * k) - ∗(RG(g * h * k))) + 2GH^2 + 2K^2 == Y^2
        @test -(2𝕀 - RG(g * h * k) - ∗(RG(g * h * k))) +
              2(GH^2 - 2G^2 - 2H^2) +
              4G^2 +
              4H^2 +
              2K^2 == Y^2

        @test GH^2 - 2G^2 - 2H^2 == -X^2
        @test -(2𝕀 - RG(g * h * k) - ∗(RG(g * h * k))) + 4G^2 + 4H^2 + 2K^2 == 2X^2 + Y^2

        @test GH^2 == 2G^2 + 2H^2 - (2𝕀 - ∗(RG(g)) - RG(h))^2
        @test KL^2 == 2K^2 + 2L^2 - (2𝕀 - ∗(RG(k)) - RG(l))^2

        @test -(2𝕀 - ∗(RG(g * h * k * l)) - RG(g * h * k * l)) + 2 * GH^2 + 2 * KL^2 ==
              (2𝕀 - ∗(RG(g * h)) - RG(k * l))^2

        @test -(2𝕀 - ∗(RG(g * h * k * l)) - RG(g * h * k * l)) +
              2(2G^2 + 2H^2 - (2𝕀 - ∗(RG(g)) - RG(h))^2) +
              2(2K^2 + 2L^2 - (2𝕀 - ∗(RG(k)) - RG(l))^2) ==
              (2𝕀 - ∗(RG(g * h)) - RG(k * l))^2

        @test -(2𝕀 - ∗(RG(g * h * k * l)) - RG(g * h * k * l)) +
              2(2G^2 + 2H^2) +
              2(2K^2 + 2L^2) ==
              (2𝕀 - ∗(RG(g * h)) - RG(k * l))^2 +
              2(2𝕀 - ∗(RG(g)) - RG(h))^2 +
              2(2𝕀 - ∗(RG(k)) - RG(l))^2

        @test -(2𝕀 - ∗(RG(g * h * k * l)) - RG(g * h * k * l)) +
              2(2𝕀 - ∗(RG(g * h * k)) - RG(g * h * k)) +
              2L^2 == (2𝕀 - ∗(RG(g * h * k)) - RG(l))^2

        @test 2𝕀 - ∗(RG(g * h * k)) - RG(g * h * k) ==
              2GH^2 + 2K^2 - (2𝕀 - ∗(RG(g * h)) - RG(k))^2

        @test -(2𝕀 - ∗(RG(g * h * k * l)) - RG(g * h * k * l)) +
              2(2GH^2 + 2K^2 - (2𝕀 - ∗(RG(g * h)) - RG(k))^2) +
              2L^2 == (2𝕀 - ∗(RG(g * h * k)) - RG(l))^2

        @test -(2𝕀 - ∗(RG(g * h * k * l)) - RG(g * h * k * l)) + 2(2GH^2 + 2K^2) + 2L^2 ==
              (2𝕀 - ∗(RG(g * h * k)) - RG(l))^2 + 2(2𝕀 - ∗(RG(g * h)) - RG(k))^2

        @test -(2𝕀 - ∗(RG(g * h * k * l)) - RG(g * h * k * l)) +
              8G^2 +
              8H^2 +
              4K^2 +
              2L^2 ==
              (2𝕀 - ∗(RG(g * h * k)) - RG(l))^2 +
              2(2𝕀 - ∗(RG(g * h)) - RG(k))^2 +
              4(2𝕀 - ∗(RG(g)) - RG(h))^2

        @test -(2𝕀 - ∗(RG(g * h * k * l)) - RG(g * h * k * l)) + 2GH^2 + 2KL^2 ==
              (2𝕀 - ∗(RG(g * h)) - RG(k * l))^2

        @test -(2𝕀 - ∗(RG(g * h * k * l)) - RG(g * h * k * l)) +
              2(2G^2 + 2H^2) +
              2(2K^2 + 2L^2) ==
              (2𝕀 - ∗(RG(g * h)) - RG(k * l))^2 +
              2(2𝕀 - ∗(RG(k)) - RG(l))^2 +
              2(2𝕀 - ∗(RG(g)) - RG(h))^2
    end
end


text:
using Cthulhu, OceanTurb

Δt = 1
N = 128
m_be = KPP.Model(N=N, stepper=:BackwardEuler)
m_fe = KPP.Model(N=N, stepper=:ForwardEuler)

m_fe.bcs.T.top = BoundaryCondition(Flux, 1.0)
m_be.bcs.T.top = BoundaryCondition(Flux, 1.0)

#@descend OceanTurb.calc_diffusive_lhs!(Δt, m_be.timestepper.lhs, m_be.timestepper.eqn.K, m_be.solution, m_be)
@descend OceanTurb.calc_explicit_rhs!(m_fe.timestepper.rhs, m_fe.timestepper.eqn, m_fe.solution, m_fe)

#OceanTurb.update!(m_be.bcs, m_be.timestepper.eqn, m_be.solution, m_be)

#@descend OceanTurb.calc_implicit_rhs!(m.timestepper.rhs, m.timestepper.eqn, m.solution, m)


text:
using Poly
using ISL
import ISL.API
using Printf

export run_polyhedral_model

"""
Run polyhedral model on a kernel. Does the following:
creates a context and space for ISL
create an ISL union set representing the domains and instructions in kernel
create an ISL union maps representing access relations in kernel
create an ISL schedule from intial implementation
create dependencies maps by analyzing the accesses and schedule
add loop optimal ordering to dependencies
create a new schedule from all dependencies
create an AST from the schedule
parse AST to julia code
return expression

if debug is set, ISL errors are printed to stderror
if verbose is
    0: no printing (unless debug is on, where ISL errors are printed)
    1: final Julia code is printed
    2: initial C code, loop orderings, new schedule, final C code, and final Julia code printed
    3: all of 2 plus domain, access relations, dependence relations, original schedule, and new schedule constraints printed
"""
function run_polyhedral_model(kernel::LoopKernel; debug=false, verbose=0, tile=-1, thread=false)::Expr
    """
    initialize context and set error printing options to warn/silence
    """
    context = ISL.API.isl_ctx_alloc()
    if debug
        ISL.API.isl_options_set_on_error(context, 0) # warn errors
    else
        ISL.API.isl_options_set_on_error(context, 1) # silence errors
    end

    """
    get the domain (instructions) of the kernel
    """
    instructions = ISL.API.isl_union_set_read_from_str(context, instructions_isl_rep(kernel))
    if verbose >= 3
        println("===DOMAIN===")
        ISL.API.isl_union_set_dump(instructions)
    end

    """
    construct the access relations
    """
    may_read, must_write = accesses_isl_rep(kernel)
    if verbose >= 3
        println("===ACCESS RELATIONS===")
        @show may_read
        @show must_write
    end
    may_read = ISL.API.isl_union_map_read_from_str(context, may_read)
    must_write = ISL.API.isl_union_map_read_from_str(context, must_write)

    """
    construct the original schedule
    """
    schedule = schedule_tree_isl_rep(context, kernel.domains[1], kernel)
    if verbose >= 3
        println("===ORIGINAL SCHEDULE===")
        ISL.API.isl_schedule_dump(schedule)
    end
    if verbose >= 2
        # show original code in C representation
        space = ISL.API.isl_set_read_from_str(context, "{:}")
        build = ISL.API.isl_ast_build_from_context(space)
        ast = ISL.API.isl_ast_build_node_from_schedule(build, ISL.API.isl_schedule_copy(schedule))
        c = ISL.API.isl_ast_node_get_ctx(ast)
        p = ISL.API.isl_printer_to_str(c)
        p = ISL.API.isl_printer_set_output_format(p, 4) # 4 = C code
        q = ISL.API.isl_printer_print_ast_node(p, ast)
        str = ISL.API.isl_printer_get_str(q)
        str = Base.unsafe_convert(Ptr{Cchar}, str)
        println("===ORIGINAL C CODE===")
        println(Base.unsafe_string(str))
    end

    """
    dependency analysis
    """
    # read after write deps
    access = ISL.API.isl_union_access_info_from_sink(ISL.API.isl_union_map_copy(may_read))
    access = ISL.API.isl_union_access_info_set_must_source(access, ISL.API.isl_union_map_copy(must_write))
    access = ISL.API.isl_union_access_info_set_schedule(access, ISL.API.isl_schedule_copy(schedule))
    flow = ISL.API.isl_union_access_info_compute_flow(access)
    raw_deps = ISL.API.isl_union_flow_get_must_dependence(flow)
    # write after read and write after write deps
    access = ISL.API.isl_union_access_info_from_sink(ISL.API.isl_union_map_copy(must_write))
    access = ISL.API.isl_union_access_info_set_must_source(access, must_write)
    access = ISL.API.isl_union_access_info_set_may_source(access, may_read)
    access = ISL.API.isl_union_access_info_set_schedule(access, schedule)
    flow = ISL.API.isl_union_access_info_compute_flow(access)
    waw_deps = ISL.API.isl_union_flow_get_must_dependence(flow)
    war_deps = ISL.API.isl_union_flow_get_may_dependence(flow)
    if verbose >= 3
        println("===DEPENDENCIES (RAW, WAW, WAR)===")
        ISL.API.isl_union_map_dump(raw_deps)
        ISL.API.isl_union_map_dump(waw_deps)
        ISL.API.isl_union_map_dump(war_deps)
    end

    """
    scheduling options
    """
    # maximize band size- better for tiling
    ISL.API.isl_options_set_schedule_maximize_band_depth(context, 1)
    # schedule weakly-connected components together
    ISL.API.isl_options_set_schedule_whole_component(context, 1)
    # avoid coalescing if possible
    ISL.API.isl_options_set_schedule_treat_coalescing(context, 1)
    # step over tile dim
    ISL.API.isl_options_set_tile_scale_tile_loops(context, 1)
    # don't start at 0 for tiles (cleans up code for 1-indexed Julia)
    ISL.API.isl_options_set_tile_shift_point_loops(context, 0)

    """
    construct schedule constraints from dependency analysis
    """
    # dependence constraints
    all_deps = ISL.API.isl_union_map_union(waw_deps, war_deps)
    all_deps = ISL.API.isl_union_map_union(all_deps, raw_deps)
    schedule_constraints = ISL.API.isl_schedule_constraints_on_domain(ISL.API.isl_union_set_copy(instructions))
    schedule_constraints = ISL.API.isl_schedule_constraints_set_validity(schedule_constraints, ISL.API.isl_union_map_copy(all_deps))

    # proximity constraints (keeps loops nested based on dependencies)
    schedule_constraints = ISL.API.isl_schedule_constraints_set_proximity(schedule_constraints, ISL.API.isl_union_map_copy(all_deps))

    # coincidence constraints (things to be scheduled together- doesn't usually do much)
    # schedule_constraints = ISL.API.isl_schedule_constraints_set_coincidence(schedule_constraints, ISL.API.isl_union_map_copy(all_deps))

    if verbose >= 3
        println("===SCHEDULE CONSTRAINTS===")
        ISL.API.isl_schedule_constraints_dump(schedule_constraints)
    end

    """
    compute the new schedule
    """
    schedule = ISL.API.isl_schedule_constraints_compute_schedule(schedule_constraints)
    if verbose >= 3
        println("===NEW SCHEDULE PRE TILING===")
        ISL.API.isl_schedule_dump(schedule)
    end

    """
    reordering of schedule loops
    """
    loop_orderings = get_best_nesting_orders(kernel)
    if verbose >= 2
        println("===LOOP ORDERINGS===")
        @show loop_orderings
    end
    reorder_schedule_loops(kernel, schedule, context, loop_orderings)
    if verbose >= 3
        println("===SCHEDULE AFTER REORDERING===")
        ISL.API.isl_schedule_dump(schedule)
    end

    """
    tiling of band nodes
    """
    if tile != 0
        schedule = tile_schedule(kernel, schedule, context, tile)
    end

    if verbose >= 2
        println("===NEW SCHEDULE===")
        ISL.API.isl_schedule_dump(schedule)
    end
    if verbose >= 3
        println("===DOMAIN AFTER TILING===")
        ISL.API.isl_union_set_dump(ISL.API.isl_schedule_get_domain(schedule))
    end

    """
    construct AST from new schedule
    """
    space = ISL.API.isl_set_read_from_str(context, "{:}")
    build = ISL.API.isl_ast_build_from_context(space)
    ast = ISL.API.isl_ast_build_node_from_schedule(build, schedule)

    if verbose >= 2
        # print C code
        c = ISL.API.isl_ast_node_get_ctx(ast)
        p = ISL.API.isl_printer_to_str(c)
        p = ISL.API.isl_printer_set_output_format(p, 4) # 4 = C code
        q = ISL.API.isl_printer_print_ast_node(p, ast)
        str = ISL.API.isl_printer_get_str(q)
        str = Base.unsafe_convert(Ptr{Cchar}, str)
        println("===FINAL C CODE===")
        println(Base.unsafe_string(str))
    end

    """
    parse back to Julia code
    """
    expr = parse_ast(ast, kernel, thread=thread, outermost=true)
    if verbose >= 1
        println("===FINAL JULIA CODE===")
        @show expr
    end
    return expr
end

"""
POLYHEDRAL MODEL GENERATION FUNCTIONS
Julia -> ISL
"""

"""
helper to turn gensyms into strings without special characters
"""
function sym_to_str(s::Symbol)::String
   str = string(s)
   str = replace.(str, r"#" => "")
   str = replace.(str, r":" => "")
   return str
end

"""
helper to return a unique identifier for a kernel existential quantifier
"""
function unique_identifier()::String
    id = gensym("id")
    return sym_to_str(id)
end

"""
helper to get string rep of params in kernel
"""
function get_params_str(kernel::LoopKernel)
    params = string(kernel.consts)
    params = replace.(params, r":" => "")
    params = replace.(params, r"Symbol" => "")
    return params
end

"""
format a condition on an instruction properly
"""
function format_instruction_cond(cond::Expr)::String
    stringcond = string(cond)
    stringcond = replace.(stringcond, r"%" => " mod ")
    stringcond = replace.(stringcond, r"==" => " = ")
    stringcond = replace.(stringcond, r"&&" => " and ")
    stringcond = replace.(stringcond, r"&" => " and ")
    stringcond = replace.(stringcond, "||" => " or ")
    stringcond = replace.(stringcond, "|" => " or ")
    return stringcond
end

"""
helper to get the domains related to instruction(s) in the ISL format
ex: 0 <= i <= n and 0 <= j <= n
"""
function get_instructions_domains(instructions::Vector{Instruction}, kernel::LoopKernel)::String

    conditions = ""
    count = 1
    for domain in kernel.domains
        for instruction in instructions
            if domain.iname in instruction.dependencies
                if count != 1
                    conditions = string(conditions, " and ")
                end
                count += 1
                instructioncond = ""
                if instruction.cond != :()
                    # instruction has separate domain condition
                    instructioncond = string(" and ", format_instruction_cond(instruction.cond))
                end
                if domain.step != 1
                    identifier = unique_identifier()
                    conditions = string(conditions, @sprintf("exists %s: %s = %s*%s + %s and %s <= %s <= %s", identifier, string(domain.iname), string(domain.step), identifier, string(domain.lowerbound), string(domain.lowerbound), string(domain.iname), string(domain.upperbound)), instructioncond)
                else
                    conditions = string(conditions, @sprintf("%s <= %s <= %s", string(domain.lowerbound), string(domain.iname), string(domain.upperbound)), instructioncond)
                end
            end
        end
    end
    return conditions
end

"""
get the related iteration spaces for an instruction
ex: [i, j]
if primes is true, returns [i', j']
"""
function get_related_domains(instruction::Instruction, kernel::LoopKernel, primes=false)
    count = 1
    domains = "["
    for domain in kernel.domains
        for iname in instruction.dependencies
            if domain.iname == iname
                if count != 1
                    domains = string(domains, ", ")
                end
                count += 1
                domains = string(domains, domain.iname)
                if primes
                    domains = string(domains, "\'")
                end
            end
        end
    end
    domains = string(domains, "]")
    return domains
end

"""
construct the instructions' domains representation in ISL syntax form
ex: (id = :mult) C[i, j] += A[i, k] * B[k, j] in matrix multiplication becomes
[n, m, r, A, B, C] -> {mult[i, j, k] : 0 <= i <= n and 0 <= j <= m and 0 <= k <= r}
"""
function instructions_isl_rep(kernel::LoopKernel)::String
    params = get_params_str(kernel)

    insn_map = "{"
    icount = 1
    for instruction in kernel.instructions
        count = 1
        name = string(sym_to_str(instruction.iname), get_related_domains(instruction, kernel))
        conditions = get_instructions_domains([instruction], kernel)

        if icount != 1
            insn_map = string(insn_map, "; ")
        end
        icount += 1

        if conditions != ""
            insn_map = string(insn_map, name, " : ", conditions)
        else
            insn_map = string(insn_map, name)
        end
    end
    insn_map = string(insn_map, "}")
    if insn_map == "{}"
        insn_map = "{:}"
    end
    return @sprintf("%s -> %s", params, insn_map)
end

"""
helpers for extracting accesses from an expression expr
"""
function expr_accesses(expr::Expr, kernel::LoopKernel)::Vector{Union{Expr, Symbol}}
    deps = []
    if expr.head == :call
        # RHS is operation on multiple terms, add all parts
        for arg in expr.args[2:end]
            append!(deps, expr_accesses(arg, kernel))
        end
    elseif expr.head == :macrocall
        # RHS is macro call on multiple terms, add all parts except macro and line number node
        for arg in expr.args[3:end]
            append!(deps, expr_accesses(arg, kernel))
        end
    else
        push!(deps, expr)
    end
    return deps
end

function expr_accesses(sym::Symbol, kernel::LoopKernel)::Vector{Union{Expr, Symbol}}
    if sym in kernel.consts
        return []
    end
    if sym in [d.iname for d in kernel.domains]
        return []
    end
    return [:($sym[])]
end

function expr_accesses(num::Number, kernel::LoopKernel)::Vector{Union{Expr, Symbol}}
    return []
end

function expr_accesses(l::LineNumberNode, kernel::LoopKernel)::Vector{Union{Expr, Symbol}}
    return []
end

"""
construct the access relations map in ISL syntax form

ex: (id = :mult) C[i, j] += A[i, k] * B[k, j] becomes
may_read:
[n, A, B, C] -> {
mult[i, j, k] -> A[k, i] : 0 <= i <= n and 0 <= j <= n and 0 <= k <= n;
mult[i, j, k] -> B[j, k] : 0 <= i <= n and 0 <= j <= n and 0 <= k <= n}
must_write:
[n, A, B, C] -> {mult[i, j, k] -> C[i, j] : 0 <= i <= n and 0 <= j <= n and 0 <= k <= n}
"""
function accesses_isl_rep(kernel::LoopKernel)::Tuple{String, String}
    params = get_params_str(kernel)

    may_read = string(params, " -> {")
    must_write = string(params, " -> {")
    rcount = 1
    wcount = 1
    mcount = 1
    for instruction in kernel.instructions
        ds = get_related_domains(instruction, kernel)
        reads = expr_accesses(instruction.body.args[2], kernel)
        writes = expr_accesses(instruction.body.args[1], kernel)
        head = instruction.body.head
        if head != :(=)
            # read from LHS also
            append!(reads, writes)
        end
        # must writes in LHS
        for part in writes
            if mcount != 1
                must_write = string(must_write, "; ")
            end
            mcount += 1
            must_write = string(must_write, sym_to_str(instruction.iname), ds, " -> ", part)
            # get domains
            conditions = get_instructions_domains([instruction], kernel)
            if conditions != ""
                must_write = string(must_write, " : ", conditions)
            end
        end
        # may reads in RHS and LHS
        for part in reads
            if rcount != 1
                may_read = string(may_read, "; ")
            end
            rcount += 1
            may_read = string(may_read, sym_to_str(instruction.iname), ds, " -> ", part)
            # get domains
            conditions = get_instructions_domains([instruction], kernel)
            if conditions != ""
                may_read = string(may_read, " : ", conditions)
            end
        end
    end

    may_read = string(may_read, "}")
    if may_read == "{}"
        may_read = "{:}"
    end
    must_write = string(must_write, "}")
    if must_write == "{}"
        must_write = "{:}"
    end
    return may_read, must_write
end

"""
Get the loop ordering dependencies for a kernel
For example, to specify the order j, k, i in a loop nest:
    "[n, r, m] -> {out253[i, j, k] -> out253[i', j', k']: j, k, i << j', k', i'}"
"""
function get_loop_ordering_deps_isl(kernel::LoopKernel, loop_orderings::Set{Vector{Symbol}})::String
    params = get_params_str(kernel)

    ordering_deps = "{"
    count = 0

    for order in loop_orderings
        for domain in kernel.domains
            if domain.iname in order
                for instruction in domain.instructions
                    if typeof(instruction) == Instruction
                        # found instruction to use
                        if count != 0
                            ordering_deps = string(ordering_deps, "; ")
                        end
                        count += 1
                        # sym[i, j] -> sym[i', j']
                        ordering_deps = string(ordering_deps, sym_to_str(instruction.iname), get_related_domains(instruction, kernel), " -> ", sym_to_str(instruction.iname), get_related_domains(instruction, kernel, true), ": ")
                        order_string = ""
                        order_string_prime = ""
                        ocount = 0
                        for iname in order
                            if iname in instruction.dependencies
                                if ocount != 0
                                    order_string = string(order_string, ", ")
                                    order_string_prime = string(order_string_prime, ", ")
                                end
                                ocount += 1
                                order_string = string(order_string, iname)
                                order_string_prime = string(order_string_prime, iname, "\'")
                            end
                        end
                        # i, j << i', j'
                        ordering_deps = string(ordering_deps, order_string, " << ", order_string_prime)
                        break
                    end
                end
            end
        end
    end

    ordering_deps = string(ordering_deps, "}")
    return @sprintf("%s -> %s", params, ordering_deps)
end

"""
construct a partial schedule for a grouping of instructions in domains
"""
function partial_schedule_isl(domains::Vector{Domain}, instructions::Vector{Instruction}, kernel::LoopKernel)::String

    mupa = "["

    dcount = 0
    for domain in domains
        set = ""

        icount = 0
        for instruction in instructions
            if domain.iname in instruction.dependencies
                if icount != 0
                    set = string(set, "; ")
                end
                icount += 1
                set = string(set, sym_to_str(instruction.iname), get_related_domains(instruction, kernel), "-> [(", domain.iname, ")]")
            end
        end

        if set != ""
            if dcount != 0
                mupa = string(mupa, ", ")
            end
            dcount += 1
            mupa = string(mupa, "{", set, "}")
        end
    end

    mupa = string(mupa, "]")

    return mupa
end

"""
construct an original schedule for dependency analysis
"""
function schedule_tree_isl_rep(context::Ptr{ISL.API.isl_ctx}, domain::Domain, kernel::LoopKernel, parent_doms=[])::Ptr{ISL.API.isl_schedule}

    new_parent_doms = [domain]
    append!(new_parent_doms, parent_doms)
    schedule = :nothing
    # navigate through original scheduling
    for instruction in domain.instructions
        if typeof(instruction) == Instruction
            # construct instruction schedule
            temp_kern = LoopKernel([instruction], kernel.domains, kernel.args, kernel.argtypes, kernel.consts)
            inst_dom = instructions_isl_rep(temp_kern)
            inst_dom = ISL.API.isl_union_set_read_from_str(context, inst_dom)
            inst_schedule = ISL.API.isl_schedule_from_domain(inst_dom)

            # put in sequence
            if schedule == :nothing
                schedule = inst_schedule
            else
                schedule = ISL.API.isl_schedule_sequence(schedule, inst_schedule)
            end
        else # Domain
            dom_schedule = schedule_tree_isl_rep(context, instruction, kernel, new_parent_doms)
            if schedule == :nothing
                schedule = dom_schedule
            else
                schedule = ISL.API.isl_schedule_sequence(schedule, dom_schedule)
            end
        end
    end

    # insert partial schedule
    insts = Instruction[]
    for instruction in kernel.instructions
        if domain.iname in instruction.dependencies
            push!(insts, instruction)
        end
    end
    partial = partial_schedule_isl([domain], insts, kernel)
    if partial != "[]"
        partial = ISL.API.isl_multi_union_pw_aff_read_from_str(context, partial)
        schedule = ISL.API.isl_schedule_insert_partial_schedule(schedule, partial)
    end

    if schedule == :nothing
        error("Mis-defined initial schedule! Check that instructions are properly nested")
    end

    return schedule
end

"""
PARSING FUNCTIONS
ISL -> Julia
"""

"""
parse an isl AST back into julia code
"""
function parse_ast(ast::Ptr{ISL.API.isl_ast_node}, kernel::LoopKernel; thread=false, outermost=false)::Expr
    node_type = ISL.API.isl_ast_node_get_type(ast)

    if node_type == ISL.API.isl_ast_node_for
        # for loop
        expr = parse_ast_for(ast, kernel, thread=thread, outermost=outermost)

    elseif node_type == ISL.API.isl_ast_node_if
        # if statement
        expr = parse_ast_if(ast, kernel, thread=thread)

    elseif node_type == ISL.API.isl_ast_node_block
        # compound node- multiple elements in block
        expr = parse_ast_block(ast, kernel, thread=thread)

    elseif node_type == ISL.API.isl_ast_node_mark
        # mark in the ast- comment

    elseif node_type == ISL.API.isl_ast_node_user
        # expression statement
        expr = parse_ast_user(ast, kernel, thread=thread)

    elseif node_type == ISL.API.isl_ast_node_error
        # error- quit

    end
end

"""
parse an ast node representing a for loop into an expression
"""
function parse_ast_for(ast::Ptr{ISL.API.isl_ast_node}, kernel::LoopKernel; thread=false, outermost=false)::Expr
    executes_once = ISL.API.isl_ast_node_for_is_degenerate(ast) # if loop executes only once

    # loop qualities
    iter = ISL.API.isl_ast_node_for_get_iterator(ast)
    id = ISL.API.isl_ast_expr_get_id(iter) # id of iterator
    ISL.API.isl_id_free(id)
    name = Base.unsafe_convert(Ptr{Cchar}, ISL.API.isl_id_get_name(id))
    name = Symbol(Base.unsafe_string(name)) # iterator symbol
    lb = parse_ast_expr(ISL.API.isl_ast_node_for_get_init(ast)) # initial condition
    body = parse_ast(ISL.API.isl_ast_node_for_get_body(ast), kernel) # body of loop

    if executes_once == ISL.API.isl_bool_true
        # loop only runs once, just set name = init and run body
        expr = quote
            $name = $init
            $body
        end
        return expr
    else
        step = parse_ast_expr(ISL.API.isl_ast_node_for_get_inc(ast)) # incremental step
        cond = parse_ast_expr(ISL.API.isl_ast_node_for_get_cond(ast)) # final condition
        ub = cond.args[3]
        if cond.args[1] == :(<)
            ub = :($ub - 1)
        elseif cond.args[1] != :(<=)
            error("unexpected conditional operator ", cond.head)
        end

        if outermost && thread
            if step == 1
                # for name = lb:up
                    # body
                expr = quote
                    Base.Threads.@threads for $name = $lb:$ub
                        $body
                    end
                end
                return expr

            else
                # for name = lb:step:up
                    # body
                expr = quote
                    Base.Threads.@threads for $name = $lb:$step:$ub
                        $body
                    end
                end
                return expr
            end
        else
            if step == 1
                # for name = lb:up
                    # body
                expr = quote
                    for $name = $lb:$ub
                        $body
                    end
                end
                return expr

            else
                # for name = lb:step:up
                    # body
                expr = quote
                    for $name = $lb:$step:$ub
                        $body
                    end
                end
                return expr
            end
        end
    end
end

"""
parse an ast node representing an if block into an expression
"""
function parse_ast_if(ast::Ptr{ISL.API.isl_ast_node}, kernel::LoopKernel; thread=false)::Expr
    cond = ISL.API.isl_ast_node_if_get_cond(ast)
    cond = parse_ast_expr(cond)
    if_body = ISL.API.isl_ast_node_if_get_then(ast)
    if_body = parse_ast(if_body, kernel, thread=thread)
    has_else = ISL.API.isl_ast_node_if_has_else(ast)
    if has_else == ISL.API.isl_bool_true
        # if cond
            # if_body
        # else
            # else_body
        # end
        else_body = ISL.API.isl_ast_node_if_get_else(ast)
        else_body = parse_ast(else_body, kernel, thread=thread)
        expr = quote
            if $cond
                $if_body
            else
                $else_body
            end
        end
        return expr
    else
        # if cond
            # if_body
        # end
        expr = quote
            if $cond
                $if_body
            end
        end
        return expr
    end
end

"""
parse an ast node representing a block into an expression
"""
function parse_ast_block(ast::Ptr{ISL.API.isl_ast_node}, kernel::LoopKernel; thread=false)::Expr
    children = ISL.API.isl_ast_node_block_get_children(ast)
    n = ISL.API.isl_ast_node_list_n_ast_node(children)
    exs = []
    for i=0:n-1
        child = ISL.API.isl_ast_node_list_get_at(children, i)
        push!(exs, parse_ast(child, kernel, thread=thread))
    end
    return quote $(exs...) end
end

"""
helpers for replacing symbols in expression
"""
function replace_expr_syms(ex::Expr, sym_dict::Dict{Symbol, Union{Symbol, Number, Expr}})::Expr
    new_args = []
    for arg in ex.args
        push!(new_args, replace_expr_syms(arg, sym_dict))
    end
    new_ex = Expr(ex.head, new_args...)
    return new_ex
end

function replace_expr_syms(sym::Symbol, sym_dict::Dict{Symbol, Union{Symbol, Number, Expr}})
    if haskey(sym_dict, sym)
        return sym_dict[sym]
    end
    return sym
end

function replace_expr_syms(num::Number, sym_dict::Dict{Symbol, Union{Symbol, Number, Expr}})
    return num
end

function replace_expr_syms(l::LineNumberNode, sym_dict::Dict{Symbol, Union{Symbol, Number, Expr}})
    return l
end

"""
parse an ast node representing an expression
"""
function parse_ast_user(ast::Ptr{ISL.API.isl_ast_node}, kernel::LoopKernel; thread=false)::Expr
    expr = ISL.API.isl_ast_node_user_get_expr(ast)
    ISL.API.isl_ast_expr_free(expr)
    first_expr = ISL.API.isl_ast_expr_op_get_arg(expr, 0) # first arg is "function" name (symbol cooresponding to instruction)
    n = ISL.API.isl_ast_expr_get_op_n_arg(expr)
    id = ISL.API.isl_ast_expr_id_get_id(first_expr)
    ISL.API.isl_ast_expr_free(first_expr)
    name = Base.unsafe_convert(Ptr{Cchar}, ISL.API.isl_id_get_name(id)) # name of identifier
    name = Base.unsafe_string(name)
    ISL.API.isl_id_free(id)

    # search for matching instruction
    for instruction in kernel.instructions
        if sym_to_str(instruction.iname) == name
            # found instruction
            # replace original iterators with new iterators
            old_ids = []
            for domain in kernel.domains
                for iname in instruction.dependencies
                    if domain.iname == iname
                        push!(old_ids, iname)
                    end
                end
            end
            new_ids = Dict{Symbol, Union{Symbol, Number, Expr}}()
            for i=1:Int(n)-1
                arg = ISL.API.isl_ast_expr_op_get_arg(expr, i)
                arg = parse_ast_expr(arg)
                new_ids[old_ids[i]] = arg
            end
            ex = replace_expr_syms(instruction.body, new_ids)
            return :(@inbounds $ex)
        end
    end

    return :()
end

"""
parse an ast expression condition
"""
function parse_ast_expr(expr::Ptr{ISL.API.isl_ast_expr})::Union{Symbol, Expr, Number}
    type = ISL.API.isl_ast_expr_get_type(expr)
    if type == ISL.API.isl_ast_expr_id
        id = ISL.API.isl_ast_expr_id_get_id(expr)
        name = Base.unsafe_convert(Ptr{Cchar}, ISL.API.isl_id_get_name(id)) # name of identifier
        name = Base.unsafe_string(name)
        ISL.API.isl_id_free(id)
        return :($(Symbol(name)))
    elseif type == ISL.API.isl_ast_expr_int
        val = ISL.API.isl_ast_expr_int_get_val(expr)
        num = ISL.API.isl_val_get_d(val)
        if Int(num) == num
            num = Int(num)
        end
        ISL.API.isl_val_free(val)
        return :($num)
    else
        if ISL.API.isl_ast_expr_get_op_n_arg(expr) == 1
            arg1 = parse_ast_expr(ISL.API.isl_ast_expr_op_get_arg(expr, 0))
            arg2 = :()
        else
            arg1 = parse_ast_expr(ISL.API.isl_ast_expr_op_get_arg(expr, 0))
            arg2 = parse_ast_expr(ISL.API.isl_ast_expr_op_get_arg(expr, 1))
        end

        op_type = ISL.API.isl_ast_expr_op_get_type(expr)
        if op_type == ISL.API.isl_ast_expr_op_and || op_type == ISL.API.isl_ast_expr_op_and_then
            return :($arg1 && $arg2)
        elseif op_type == ISL.API.isl_ast_expr_op_or || op_type == ISL.API.isl_ast_expr_op_or_else
            return :($arg1 || $arg2)
        elseif op_type == ISL.API.isl_ast_expr_op_max
            return :(max($arg1, $arg2))
        elseif op_type == ISL.API.isl_ast_expr_op_min
            return :(min($arg1, $arg2))
        elseif op_type == ISL.API.isl_ast_expr_op_minus
            return :(-$arg1)
        elseif op_type == ISL.API.isl_ast_expr_op_add
            return :($arg1 + $arg2)
        elseif op_type == ISL.API.isl_ast_expr_op_sub
            return :($arg1 - $arg2)
        elseif op_type == ISL.API.isl_ast_expr_op_mul
            return :($arg1 * $arg2)
        elseif op_type == ISL.API.isl_ast_expr_op_fdiv_q
            return :(Int(floor($arg1 / $arg2)))
        elseif op_type == ISL.API.isl_ast_expr_op_pdiv_q || op_type == ISL.API.isl_ast_expr_op_div
            return :(Int($arg1 / $arg2))
        elseif op_type == ISL.API.isl_ast_expr_op_pdiv_r || op_type == ISL.API.isl_ast_expr_op_zdiv_r
            return :($arg1 % $arg2)
        elseif op_type == ISL.API.isl_ast_expr_op_eq
            return :($arg1 == $arg2)
        elseif op_type == ISL.API.isl_ast_expr_op_le
            return :($arg1 <= $arg2)
        elseif op_type == ISL.API.isl_ast_expr_op_ge
            return :($arg1 >= $arg2)
        elseif op_type == ISL.API.isl_ast_expr_op_lt
            return :($arg1 < $arg2)
        elseif op_type == ISL.API.isl_ast_expr_op_gt
            return :($arg1 > $arg2)
        else
            # invalid
            return :()
        end
    end
end


text:
using BenchmarkTools

using ArnoldiMethod: partialschur, SM, LM
using LinearMaps
using LinearAlgebra
using SparseArrays
using DelimitedFiles

mymatrix(n) = spdiagm(-1 => fill(-1.0, n-1), 0 => fill(2.0, n), 1 => fill(-1.001, n-1))

struct ShiftAndInvert{TA,TB,TT}
    A_lu::TA
    B::TB
    temp::TT
end

function (M::ShiftAndInvert)(y,x)
    mul!(M.temp, M.B, x)
    # @show typeof(y), typeof(M.A_lu), typeof(M.temp)
    ldiv!(y, M.A_lu, M.temp)
end

function construct_linear_map(A,B)
    a = ShiftAndInvert(factorize(A),B,Vector{eltype(A)}(undef, size(A,1)))
    LinearMap(a, size(A,1), ismutating=true)
end

# For inverting in non-generealized case
function construct_linear_map(A)
    a = factorize(A)
    function f(y,x)
        ldiv!(y, a, x)
    end
    LinearMap{eltype(A)}(f, size(A,1))
end

function bencharnoldi()
    A = mymatrix(6000)
    A = construct_linear_map(mymatrix(6000))
    target = LM()

    @benchmark partialschur(B, mindim=11, maxdim=22, nev=10, tol=1e-10, restarts=100000, which=tar) setup = (B = $A; tar = $target)

end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ cb205c40-83a9-11eb-329c-e529ab101736
begin
	cd("..")
	using Pkg
	Pkg.activate(".")
end

# ╔═╡ 80030c2c-83ab-11eb-2a2e-39f1fcbcffc3
using ColabINAOE2021

# ╔═╡ 1c6fa178-83ad-11eb-154f-a11ce04426dd
using Plots

# ╔═╡ 5d4a918e-83b3-11eb-2445-f79669bef553
using LaTeXStrings

# ╔═╡ e90745a4-83ac-11eb-3c40-6f6d11de047e
using LinearAlgebra

# ╔═╡ 31a312be-83a8-11eb-365c-cd02c33f3b94
md" ## Energy exchange studies in the masses array evolution "

# ╔═╡ cdcb9f1a-83a8-11eb-14fe-67c783a8f2da
md"

In [1] we stipulate how to obtain the time evolution of initial amplitude arrays $\mathbf{Q}(0) = \sum_{l}q_{l}(0)\ket{l}$ under a sort of harmonic nearest-neighbour couplings given by coupling matrices $\mathbb{M}$. We can obtain closed-analytics solutions when the matrices $\mathbb{M}$ could be factorized as $\mathbb{A}\mathbb{B}\mathbb{A}^{\dagger}$ over a discrete and finite orthonormal space, like the Fock space of bras and kets.

-- Energy exchange: We define the energy exchange of one element _i_ to other element _j_ in an array as the function

$E_{i\rightarrow j} = k\left(q_{i}\dot{q}_{j} - \dot{q}_{i}q_{j}\right),$

for some restitution constant $k$. Because we have the explicit functions for three cases of arrays: one circular, and two lineals; we can obtain the energy exchange in such simple form, just obtaining the time-derivatives and composing the function.

To notice, is that the energy exchange depends on time, so we have specifically the energy exchange rate, that is, how much an element _j_ receive (or deliver) to an element _i_ in a timespan.

_____

[1] Urzúa, A. R., Ramos-Prieto, I., Soto-Eguibar, F., & Moya-Cessa, H. (2020). Dynamical analysis of mass–spring models using Lie algebraic methods. Physica A: Statistical Mechanics and Its Applications, 540, 123193. [doi](https://doi.org/10.1016/j.physa.2019.123193)
"

# ╔═╡ c00b788e-83ad-11eb-3456-1937a68cc057
pyplot()

# ╔═╡ 46a432d4-83ac-11eb-2410-93ca49a73b74
# Global Constants

begin
	N = 32
	dim = N + 1
	k = 1
end;

# ╔═╡ 7978a3de-83ac-11eb-3ec2-ef162c5ca306
t_list = range(0, stop = 200, length = 200)

# ╔═╡ b22933a6-83ac-11eb-36e0-53bda1f79066
begin
	v0 = zeros(dim)
	v0[1] = 1
	#v0[] = 1

	q0 = normalize(v0)
end;

# ╔═╡ 84833ef8-83b4-11eb-3b7d-0388d5b1a38e
md"### 1. Circular array "

# ╔═╡ f9f33a28-83b2-11eb-3b1c-f7db876a1e0a
plot([qc(1, 1, t, N, q0) for t in t_list], linewidth = 2, label = L"q_{1}(t)", xlabel = L"Time\; [t]", ylabel = L"Amplitude")

# ╔═╡ 0a5165a6-83af-11eb-3d16-91b1821ac83c
ec01 = wijct(1, 0, 1, t_list, N, q0);

# ╔═╡ 2e85e55c-83af-11eb-2b62-f7bcf8192c37
plot(ec01, linewidth = 2, label = L"E_{0\rightarrow 1}(t)", xlabel = L"Time\; [t]", ylabel = L"Energy")

# ╔═╡ 9616dc6a-83b4-11eb-3f55-3b3b13ae6a94
md" ### 2. Lineal array "

# ╔═╡ a103ac0a-83b4-11eb-0766-cdea31fc445e
plot([ql(1, 1, t, N, q0) for t in t_list], linewidth = 2, label = L"q_{1}(t)", xlabel = L"Time\; [t]", ylabel = L"Amplitude")

# ╔═╡ c646c242-83b4-11eb-009d-038d9bc9df7c
el01 = wijlt(1, 0, 1, t_list, N, q0);

# ╔═╡ ce72d0a0-83b4-11eb-11be-e38b4881c08d
plot(el01, linewidth = 2, label = L"E_{0\rightarrow 1}(t)", xlabel = L"Time\; [t]", ylabel = L"Energy")

# ╔═╡ 190f4c7e-83b5-11eb-2f75-7559d8264ede
md" ### 3. Kravchuk array "

# ╔═╡ 284ee190-83b5-11eb-0940-917c1c56f682
plot([qk(1, t, N, q0) for t in t_list], linewidth = 2, label = L"q_{1}(t)", xlabel = L"Time\; [t]", ylabel = L"Amplitude")

# ╔═╡ 3407b2f0-83b5-11eb-3003-6787a4694722
ek01 = wijkt(0, 1, t_list, N, q0);

# ╔═╡ a37ee8c8-83b6-11eb-32e4-51d7e3407b42
plot(ek01, linewidth = 2, label = L"E_{0\rightarrow 1}(t)", xlabel = L"Time\; [t]", ylabel = L"Energy")

# ╔═╡ e6099ea2-83b6-11eb-102c-45960096ab4e
[sum(ec01), sum(el01), sum(ek01)]

# ╔═╡ Cell order:
# ╟─31a312be-83a8-11eb-365c-cd02c33f3b94
# ╟─cdcb9f1a-83a8-11eb-14fe-67c783a8f2da
# ╠═cb205c40-83a9-11eb-329c-e529ab101736
# ╠═80030c2c-83ab-11eb-2a2e-39f1fcbcffc3
# ╠═1c6fa178-83ad-11eb-154f-a11ce04426dd
# ╠═5d4a918e-83b3-11eb-2445-f79669bef553
# ╠═c00b788e-83ad-11eb-3456-1937a68cc057
# ╠═e90745a4-83ac-11eb-3c40-6f6d11de047e
# ╠═46a432d4-83ac-11eb-2410-93ca49a73b74
# ╠═7978a3de-83ac-11eb-3ec2-ef162c5ca306
# ╠═b22933a6-83ac-11eb-36e0-53bda1f79066
# ╟─84833ef8-83b4-11eb-3b7d-0388d5b1a38e
# ╠═f9f33a28-83b2-11eb-3b1c-f7db876a1e0a
# ╠═0a5165a6-83af-11eb-3d16-91b1821ac83c
# ╠═2e85e55c-83af-11eb-2b62-f7bcf8192c37
# ╟─9616dc6a-83b4-11eb-3f55-3b3b13ae6a94
# ╠═a103ac0a-83b4-11eb-0766-cdea31fc445e
# ╠═c646c242-83b4-11eb-009d-038d9bc9df7c
# ╠═ce72d0a0-83b4-11eb-11be-e38b4881c08d
# ╟─190f4c7e-83b5-11eb-2f75-7559d8264ede
# ╠═284ee190-83b5-11eb-0940-917c1c56f682
# ╠═3407b2f0-83b5-11eb-3003-6787a4694722
# ╠═a37ee8c8-83b6-11eb-32e4-51d7e3407b42
# ╠═e6099ea2-83b6-11eb-102c-45960096ab4e


text:
# # Chapter 4 - 回帰分析
# ## 4.1 線形回帰分析
## #src
include("../src/data.jl")
include("../src/linear_regression.jl")

using Measures
using Plots

## #src
#-
## 1. データの作成
data = regression_data1()

## #src
#-
function linear_regression(data; normalize=false)
    ## 2. データを学習用と評価用に分割
    n_train = floor(Int64, size(data.X, 1) * 0.9)

    ## 学習データ(全体の90%)
    X = normalize ? (data.X .- mean(data.X; dims = 1)) ./ std(data.X; corrected = false, dims = 1) : data.X
    Y = normalize ? (data.Y .- mean(data.Y; dims = 1)) ./ std(data.Y; corrected = false, dims = 1) : data.Y
    x_train = X[1:n_train, :]
    y_train = Y[1:n_train]

    ## 評価データ(全体の10%)
    x_test = X[n_train + 1:end, :]
    y_test = Y[n_train + 1:end]

    ## 3. 線形モデルの学習
    model = LinearRegression(x_train, y_train)
    train!(model)
    println("モデルパラメーター:\nw=$(model.w), \nb=$(model.b)")
    println("平方平均二乗誤差=$(rmse(model, x_test, y_test))ドル")
    println("決定係数=$(r2(model, x_test, y_test))")
    return model
end

## #src
#-
## 4. 線形モデルの評価
model = linear_regression(data)

## #src
#-
## 5. 線形モデルのプロット
scatter(data.X, data.Y,
        label = "Data",
        legend = :topleft,
        xlabel = "Garage Area x (square feet)",
        ylabel = "Sale Price y (USD)",
        margin = 12mm)
xs = 0:100:6000
## 予測を計算する。xs は reshapeする必要あり
ys = predict(model, reshape(xs, :, 1))
plot!(xs, ys, linewidth = 2, label = "Linear Regression Model")

## #src
# ### 4.1.9 重回帰分析の実行
data2 = regression_data2()
model2 = linear_regression(data2)

## #src
# ### 4.1.11 標準化による説明変数の重要度分析
model3 = linear_regression(data2, normalize = true)

## #src
# ### 4.1.12 外れ値と正則化
scatter(data.X, data.Y,
        label = "Train data without outlier",
        legend = :topleft,
        xlabel = "Garage Area x (square feet)",
        ylabel = "Sale Price y (USD)",
        margin = 12mm)
data_outlier = regression_data3()
scatter!(data_outlier.X, data_outlier.Y,
        label = "Trains data with outlier")
scatter!(data.X[data.Y .< 700000, :], data.Y[data.Y .< 700000],
        label = "Train data")
plot!(xs, ys, linewidth = 2, label = "Linear Regression Model (No outlier)")
model_outlier = linear_regression(data_outlier)
ys_outlier = predict(model_outlier, reshape(xs, :, 1))
plot!(xs, ys_outlier, linewidth = 2, label = "Linear Regression Model (With outlier)")

## #src
#-
## L2正則平均二乗誤差の正則化
## 正則化項なしの場合
scatter(data.X, data.Y,
        label = "Data",
        legend = :topleft,
        xlabel = "Garage Area x (square feet)",
        ylabel = "Sale Price y (USD)",
        margin = 12mm)
plot!(xs, ys_outlier, linewidth = 2, label = "Linear Regression Model")

## #src
#-
## 正則化項あり(lamb=1)
model_reg = linear_regression(data_outlier)
train_regularized!(model_reg, 1.0)
scatter(data.X, data.Y,
        label = "Data",
        legend = :topleft,
        xlabel = "Garage Area x (square feet)",
        ylabel = "Sale Price y (USD)",
        margin = 12mm)
xs = 0:100:6000
ys_reg = predict(model_reg, reshape(xs, :, 1))
plot!(xs, ys_reg, linewidth = 2, label = "Linear Regression Model")


text:
function StenoGraph(ex)
    if ex.head == :block
        exs = filter(x -> !isa(x, LineNumberNode), ex.args)
        exs = StenoGraphs.variable_as_node!.(exs)
        exs = StenoGraphs.addition_to_vector!.(exs)
        vec = Expr(:call, :vcat, exs...)
        return esc(:(meld(vcat(StenoGraphs.unarrow.($vec)...))))
    else
        StenoGraph(Expr(:block, ex))
    end
end

macro StenoGraph(ex)
    StenoGraph(ex)
end


text:
module AbstractIterativeSolvers

using LinearAlgebra

export GMRES, GMRES_verbose, GMRES_hi

function GMRES(A,b,inner,tol,n,cond)
    nom = a -> sqrt(abs(inner(a,a)))
    H = zeros(Complex{Float64},n+1,n)
    bnorm = nom(b)
    x = 0.
    conv_history = []
    Q = [(1.0/bnorm)*b]
    for i = 1:n
       #tic()
       #println("Operator application: ")
       v = A(Q[i])
       #toc()
       #tic()
       #println("Inner products: ")
       for j = 1:i
           #tic()
           H[j,i] = inner(Q[j],v)
           #toc()
           v = cond(v - H[j,i]*Q[j])
       end
       v = cond(v)
       #println("Assembling Q:")
       H[i+1,i] = nom(v)
       Q = vcat(Q,[copy((1.0/H[i+1,i])*v)])
       #print("Arnoldi: ")
       #toc()
       #return v
       if i > 1
           # Solve H[1:i+1,1:i]*x = bnorm*e_1, using least squares
           # TODO: Implement Givens rotations
           rhs = zeros(Float64,i+1)
           rhs[1] = bnorm
           x = H[1:i+1,1:i]\rhs
           res = norm(H[1:i+1,1:i]*x-rhs)
           conv_history = vcat(conv_history,[i,res])
           print("iteration = ")
           print(i)
           print(", residual = ")
           println(res)
           if res < tol
               return  [Q,x,conv_history]
           end
       end
    end
    println("GMRES did not terminate")
    return [Q,x,conv_history]
end

function GMRES_hi(A,b,inner,tol,n,cond,prec,flag = false)
    nom = a -> sqrt(abs(inner(a,a)))
    H = zeros(Complex{BigFloat},n+1,n)
    bnorm = nom(b)
    x = BigFloat(0,prec)
    conv_history = []
    Q = [(1/bnorm)*b]
    for i = 1:n
       #tic()
       #println("Operator application: ")
       v = A(Q[i])
       #toc()
       #tic()
       #println("Inner products: ")
       for j = 1:i
           #tic()
           H[j,i] = inner(Q[j],v)
           #toc()
           v = cond(v - H[j,i]*Q[j])
       end
       v = cond(v)
       #println("Assembling Q:")
       H[i+1,i] = nom(v)
       Q = vcat(Q,[copy((1/H[i+1,i])*v)])
       #print("Arnoldi: ")
       #toc()
       #return v
       if i > 1
           # Solve H[1:i+1,1:i]*x = bnorm*e_1, using least squares
           # TODO: Implement Givens rotations
           rhs = zeros(BigFloat,i+1)
           rhs[1] = bnorm
           x = H[1:i+1,1:i]\rhs
           res = norm(H[1:i+1,1:i]*x-rhs)
           conv_history = vcat(conv_history,[i,res])
           if flag
               print("iteration = ")
               print(i)
               print(", residual = ")
               println(res)
           end
           if res < tol
               return  [Q,x,conv_history]
           end
       end
    end
    println("GMRES did not terminate")
    return [Q,x,conv_history]
end

function GMRES_verbose(A,b,inner,tol,n,cond)

    nom = a -> sqrt(abs(inner(a,a)))
    H = zeros(Complex{Float64},n+1,n)
    bnorm = nom(b)
    x = 0.
    Q = [(1.0/bnorm)*b]
    for i = 1:n
       #tic()
       println("Operator application: ")
       @time v = A(Q[i])
       #toc()
       #tic()
       println("Inner products: ")
       for j = 1:i
           #tic()
           @time H[j,i] = inner(Q[j],v)
           #toc()
           v = cond(v - H[j,i]*Q[j])
       end
       v = cond(v)
       println("Assembling Q:")
       @time H[i+1,i] = nom(v)
       @time Q = vcat(Q,[copy((1.0/H[i+1,i])*v)])
       #print("Arnoldi: ")
       #toc()
       #return v
       if i > 1
           # Solve H[1:i+1,1:i]*x = bnorm*e_1, using least squares
           # TODO: Implement Givens rotations
           rhs = zeros(Float64,i+1)
           rhs[1] = bnorm
           println("Solve least squares")
           @time x = H[1:i+1,1:i]\rhs
           res = norm(H[1:i+1,1:i]*x-rhs)
           print("iteration = ")
           print(i)
           print(", residual = ")
           println(res)
           if res < tol
               return [Q,x]
           end
       end
    end
    println("GMRES did not terminate")
    return[Q,x]
end

end # module


text:
import Voronoi

using Random
using Printf
include("Draw.jl")

function randf(start, finish, n)
   v = rand(n)
   return map(x -> start + x*(finish-start), v)
end
Random.seed!(100)

WIDTH = 100.0
HEIGHT = 100.0

function benchmarkLloyd(Z::Array{Tuple{Real, Real}, 1}, ϵ::Number)
   global WIDTH, HEIGHT

   points = deepcopy(Z)

   results = Array{Tuple{Real, Real, Real, Real}}([])

   times = @elapsed V, f, ∇f, ξ = Voronoi.Optimization.init(points, WIDTH, HEIGHT)
   evaluations = 1
   push!(results, (f, ξ, evaluations, times))

   while ξ > ϵ
      time = @elapsed V, f, ξ = Voronoi.Optimization.lloydIteration(points, V)

      evaluations += 1
      times += time

      push!(results, (f, ξ, evaluations, times))
   end

   return results
end

function benchmarkGradient(Z::Array{Tuple{Real, Real}, 1}, ϵ::Number)
   global WIDTH, HEIGHT

   points = deepcopy(Z)

   results = Array{Tuple{Real, Real, Real, Real}}([])

   times = @elapsed V, f, ∇f, ξ = Voronoi.Optimization.init(points, WIDTH, HEIGHT)
   evaluations = 1
   push!(results, (f, ξ, evaluations, times))
   μ = 0.5
   α = 0.01
   λ₀ = 0.01

   while ξ > ϵ
      d = ∇f

      time = @elapsed V, f, ∇f, ξ, evals = Voronoi.Optimization.lineSearch(points, d, α, λ₀, μ, ϵ, V, f, ∇f, ξ)

      evaluations += evals
      times += time

      push!(results, (f, ξ, evaluations, times))
   end

   return results
end

function benchmark(n, ϵ)
   global WIDTH, HEIGHT

   points = convert(Array{Tuple{Real, Real}}, collect(zip(randf(1, WIDTH-1, n), randf(1, HEIGHT-1, n))))

   lloyd_results = benchmarkLloyd(points, ϵ)
   gradient_results = benchmarkGradient(points, ϵ)

   lloyd_f = Array{Number}([])
   lloyd_ξ = Array{Number}([])
   lloyd_evals = Array{Number}([])
   lloyd_time = Array{Number}([])
   gradient_f = Array{Number}([])
   gradient_ξ = Array{Number}([])
   gradient_evals = Array{Number}([])
   gradient_time = Array{Number}([])

   i = 1
   t₀ = lloyd_results[1][4]
   for (fᵢ, ξᵢ, evalsᵢ, timeᵢ) in lloyd_results
      push!(lloyd_f, fᵢ)
      push!(lloyd_ξ, ξᵢ)
      push!(lloyd_evals, evalsᵢ)
      push!(lloyd_time, timeᵢ - t₀)
      i += 1
   end

   i = 1
   t₀ = gradient_results[1][4]
   for (fᵢ, ξᵢ, evalsᵢ, timeᵢ) in gradient_results
      push!(gradient_f, fᵢ)
      push!(gradient_ξ, ξᵢ)
      push!(gradient_evals, evalsᵢ)
      push!(gradient_time, timeᵢ - t₀)
      i += 1
   end

   #@printf("f = %.10g\t\t|∇f| = %.10g\t\ttime: %s seconds\n", f, ξ, time)
   #@printf("f = %.10g\t\t|∇f| = %.10g\t\ttime: %s seconds\n", f, ξ, time)

   Draw.init(0, 0)
   Draw.clear(string("n = ", n))
   Draw.plot(lloyd_time, lloyd_f, "xkcd:azure", "Algoritmo de Lloyd", "-")
   Draw.plot(gradient_time, gradient_f, "xkcd:tomato", "Método do Gradiente", ":")
   Draw.legend("Tempo (s)", "Valor da Função Objetivo")
   Draw.commit()
   println("Press Return to finish.")
   command = readline(stdin)

   Draw.init(0, 0)
   Draw.clear(string("n = ", n))
   Draw.plot(lloyd_time, lloyd_ξ, "xkcd:azure", "Algoritmo de Lloyd", "-")
   Draw.plot(gradient_time, gradient_ξ, "xkcd:tomato", "Método do Gradiente", ":")
   Draw.legend("Tempo (s)", "Norma do Gradiente")
   Draw.commit()
   println("Press Return to finish.")
   command = readline(stdin)
end

ϵ = 10^-2
benchmark(3, ϵ)
benchmark(30, ϵ)
benchmark(100, ϵ)
benchmark(300, ϵ)
benchmark(1000, ϵ)

max_stars_repo_name: arthursantana/coveringWithCircles

max_issues_repo_name: arthursantana/coveringWithCircles

max_forks_repo_name: arthursantana/coveringWithCircles