text:
# syntax: proto2
using ProtoBuf
import ProtoBuf.meta

const Version = (;[
    Symbol("_START_VERSION") => Int32(0),
    Symbol("IR_VERSION_2017_10_10") => Int32(1),
    Symbol("IR_VERSION_2017_10_30") => Int32(2),
    Symbol("IR_VERSION_2017_11_3") => Int32(3),
    Symbol("IR_VERSION_2019_1_22") => Int32(4),
    Symbol("IR_VERSION_2019_3_18") => Int32(5),
    Symbol("IR_VERSION_2019_9_19") => Int32(6),
    Symbol("IR_VERSION") => Int32(7),
]...)

mutable struct StringStringEntryProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}

    function StringStringEntryProto(; kwargs...)
        obj = new(meta(StringStringEntryProto), Dict{Symbol,Any}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct StringStringEntryProto
const __meta_StringStringEntryProto = Ref{ProtoMeta}()
function meta(::Type{StringStringEntryProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_StringStringEntryProto)
            __meta_StringStringEntryProto[] = target = ProtoMeta(StringStringEntryProto)
            allflds = Pair{Symbol,Union{Type,String}}[:key => AbstractString, :value => AbstractString]
            meta(target, StringStringEntryProto, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_StringStringEntryProto[]
    end
end
function Base.getproperty(obj::StringStringEntryProto, name::Symbol)
    if name === :key
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :value
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct TensorAnnotation <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}

    function TensorAnnotation(; kwargs...)
        obj = new(meta(TensorAnnotation), Dict{Symbol,Any}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct TensorAnnotation
const __meta_TensorAnnotation = Ref{ProtoMeta}()
function meta(::Type{TensorAnnotation})
    ProtoBuf.metalock() do
        if !isassigned(__meta_TensorAnnotation)
            __meta_TensorAnnotation[] = target = ProtoMeta(TensorAnnotation)
            allflds = Pair{Symbol,Union{Type,String}}[:tensor_name => AbstractString, :quant_parameter_tensor_names => Base.Vector{StringStringEntryProto}]
            meta(target, TensorAnnotation, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_TensorAnnotation[]
    end
end
function Base.getproperty(obj::TensorAnnotation, name::Symbol)
    if name === :tensor_name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :quant_parameter_tensor_names
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{StringStringEntryProto}
    else
        getfield(obj, name)
    end
end

const TensorProto_DataType = (;[
    Symbol("UNDEFINED") => Int32(0),
    Symbol("FLOAT") => Int32(1),
    Symbol("UINT8") => Int32(2),
    Symbol("INT8") => Int32(3),
    Symbol("UINT16") => Int32(4),
    Symbol("INT16") => Int32(5),
    Symbol("INT32") => Int32(6),
    Symbol("INT64") => Int32(7),
    Symbol("STRING") => Int32(8),
    Symbol("BOOL") => Int32(9),
    Symbol("FLOAT16") => Int32(10),
    Symbol("DOUBLE") => Int32(11),
    Symbol("UINT32") => Int32(12),
    Symbol("UINT64") => Int32(13),
    Symbol("COMPLEX64") => Int32(14),
    Symbol("COMPLEX128") => Int32(15),
    Symbol("BFLOAT16") => Int32(16),
]...)

const TensorProto_DataLocation = (;[
    Symbol("DEFAULT") => Int32(0),
    Symbol("EXTERNAL") => Int32(1),
]...)

mutable struct TensorProto_Segment <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}

    function TensorProto_Segment(; kwargs...)
        obj = new(meta(TensorProto_Segment), Dict{Symbol,Any}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct TensorProto_Segment
const __meta_TensorProto_Segment = Ref{ProtoMeta}()
function meta(::Type{TensorProto_Segment})
    ProtoBuf.metalock() do
        if !isassigned(__meta_TensorProto_Segment)
            __meta_TensorProto_Segment[] = target = ProtoMeta(TensorProto_Segment)
            allflds = Pair{Symbol,Union{Type,String}}[:_begin => Int64, :_end => Int64]
            meta(target, TensorProto_Segment, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_TensorProto_Segment[]
    end
end
function Base.getproperty(obj::TensorProto_Segment, name::Symbol)
    if name === :_begin
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :_end
        return (obj.__protobuf_jl_internal_values[name])::Int64
    else
        getfield(obj, name)
    end
end

mutable struct TensorProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}

    function TensorProto(; kwargs...)
        obj = new(meta(TensorProto), Dict{Symbol,Any}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct TensorProto
const __meta_TensorProto = Ref{ProtoMeta}()
function meta(::Type{TensorProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_TensorProto)
            __meta_TensorProto[] = target = ProtoMeta(TensorProto)
            fnum = Int[1,2,3,4,5,6,7,8,12,9,13,14,10,11]
            pack = Symbol[:float_data,:int32_data,:int64_data,:double_data,:uint64_data]
            allflds = Pair{Symbol,Union{Type,String}}[:dims => Base.Vector{Int64}, :data_type => Int32, :segment => TensorProto_Segment, :float_data => Base.Vector{Float32}, :int32_data => Base.Vector{Int32}, :string_data => Base.Vector{Array{UInt8,1}}, :int64_data => Base.Vector{Int64}, :name => AbstractString, :doc_string => AbstractString, :raw_data => Array{UInt8,1}, :external_data => Base.Vector{StringStringEntryProto}, :data_location => Int32, :double_data => Base.Vector{Float64}, :uint64_data => Base.Vector{UInt64}]
            meta(target, TensorProto, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, pack, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_TensorProto[]
    end
end
function Base.getproperty(obj::TensorProto, name::Symbol)
    if name === :dims
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Int64}
    elseif name === :data_type
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :segment
        return (obj.__protobuf_jl_internal_values[name])::TensorProto_Segment
    elseif name === :float_data
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Float32}
    elseif name === :int32_data
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Int32}
    elseif name === :string_data
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Array{UInt8,1}}
    elseif name === :int64_data
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Int64}
    elseif name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :doc_string
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :raw_data
        return (obj.__protobuf_jl_internal_values[name])::Array{UInt8,1}
    elseif name === :external_data
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{StringStringEntryProto}
    elseif name === :data_location
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :double_data
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Float64}
    elseif name === :uint64_data
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{UInt64}
    else
        getfield(obj, name)
    end
end

mutable struct SparseTensorProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}

    function SparseTensorProto(; kwargs...)
        obj = new(meta(SparseTensorProto), Dict{Symbol,Any}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct SparseTensorProto
const __meta_SparseTensorProto = Ref{ProtoMeta}()
function meta(::Type{SparseTensorProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_SparseTensorProto)
            __meta_SparseTensorProto[] = target = ProtoMeta(SparseTensorProto)
            allflds = Pair{Symbol,Union{Type,String}}[:values => TensorProto, :indices => TensorProto, :dims => Base.Vector{Int64}]
            meta(target, SparseTensorProto, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_SparseTensorProto[]
    end
end
function Base.getproperty(obj::SparseTensorProto, name::Symbol)
    if name === :values
        return (obj.__protobuf_jl_internal_values[name])::TensorProto
    elseif name === :indices
        return (obj.__protobuf_jl_internal_values[name])::TensorProto
    elseif name === :dims
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Int64}
    else
        getfield(obj, name)
    end
end

mutable struct TensorShapeProto_Dimension <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}

    function TensorShapeProto_Dimension(; kwargs...)
        obj = new(meta(TensorShapeProto_Dimension), Dict{Symbol,Any}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct TensorShapeProto_Dimension
const __meta_TensorShapeProto_Dimension = Ref{ProtoMeta}()
function meta(::Type{TensorShapeProto_Dimension})
    ProtoBuf.metalock() do
        if !isassigned(__meta_TensorShapeProto_Dimension)
            __meta_TensorShapeProto_Dimension[] = target = ProtoMeta(TensorShapeProto_Dimension)
            allflds = Pair{Symbol,Union{Type,String}}[:dim_value => Int64, :dim_param => AbstractString, :denotation => AbstractString]
            oneofs = Int[1,1,1]
            oneof_names = Symbol[Symbol("value")]
            meta(target, TensorShapeProto_Dimension, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, oneofs, oneof_names)
        end
        __meta_TensorShapeProto_Dimension[]
    end
end
function Base.getproperty(obj::TensorShapeProto_Dimension, name::Symbol)
    if name === :dim_value
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :dim_param
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :denotation
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct TensorShapeProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}

    function TensorShapeProto(; kwargs...)
        obj = new(meta(TensorShapeProto), Dict{Symbol,Any}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct TensorShapeProto
const __meta_TensorShapeProto = Ref{ProtoMeta}()
function meta(::Type{TensorShapeProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_TensorShapeProto)
            __meta_TensorShapeProto[] = target = ProtoMeta(TensorShapeProto)
            allflds = Pair{Symbol,Union{Type,String}}[:dim => Base.Vector{TensorShapeProto_Dimension}]
            meta(target, TensorShapeProto, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_TensorShapeProto[]
    end
end
function Base.getproperty(obj::TensorShapeProto, name::Symbol)
    if name === :dim
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{TensorShapeProto_Dimension}
    else
        getfield(obj, name)
    end
end

mutable struct OperatorSetIdProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}

    function OperatorSetIdProto(; kwargs...)
        obj = new(meta(OperatorSetIdProto), Dict{Symbol,Any}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct OperatorSetIdProto
const __meta_OperatorSetIdProto = Ref{ProtoMeta}()
function meta(::Type{OperatorSetIdProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_OperatorSetIdProto)
            __meta_OperatorSetIdProto[] = target = ProtoMeta(OperatorSetIdProto)
            allflds = Pair{Symbol,Union{Type,String}}[:domain => AbstractString, :version => Int64]
            meta(target, OperatorSetIdProto, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_OperatorSetIdProto[]
    end
end
function Base.getproperty(obj::OperatorSetIdProto, name::Symbol)
    if name === :domain
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :version
        return (obj.__protobuf_jl_internal_values[name])::Int64
    else
        getfield(obj, name)
    end
end

const AttributeProto_AttributeType = (;[
    Symbol("UNDEFINED") => Int32(0),
    Symbol("FLOAT") => Int32(1),
    Symbol("INT") => Int32(2),
    Symbol("STRING") => Int32(3),
    Symbol("TENSOR") => Int32(4),
    Symbol("GRAPH") => Int32(5),
    Symbol("SPARSE_TENSOR") => Int32(11),
    Symbol("FLOATS") => Int32(6),
    Symbol("INTS") => Int32(7),
    Symbol("STRINGS") => Int32(8),
    Symbol("TENSORS") => Int32(9),
    Symbol("GRAPHS") => Int32(10),
    Symbol("SPARSE_TENSORS") => Int32(12),
]...)

mutable struct AttributeProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}

    function AttributeProto(; kwargs...)
        obj = new(meta(AttributeProto), Dict{Symbol,Any}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct AttributeProto (has cyclic type dependency)
const __meta_AttributeProto = Ref{ProtoMeta}()
function meta(::Type{AttributeProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_AttributeProto)
            __meta_AttributeProto[] = target = ProtoMeta(AttributeProto)
            fnum = Int[1,21,13,20,2,3,4,5,6,22,7,8,9,10,11,23]
            allflds = Pair{Symbol,Union{Type,String}}[:name => AbstractString, :ref_attr_name => AbstractString, :doc_string => AbstractString, :_type => Int32, :f => Float32, :i => Int64, :s => Array{UInt8,1}, :t => TensorProto, :g => "GraphProto", :sparse_tensor => SparseTensorProto, :floats => Base.Vector{Float32}, :ints => Base.Vector{Int64}, :strings => Base.Vector{Array{UInt8,1}}, :tensors => Base.Vector{TensorProto}, :graphs => "Base.Vector{GraphProto}", :sparse_tensors => Base.Vector{SparseTensorProto}]
            meta(target, AttributeProto, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_AttributeProto[]
    end
end
function Base.getproperty(obj::AttributeProto, name::Symbol)
    if name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :ref_attr_name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :doc_string
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :_type
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :f
        return (obj.__protobuf_jl_internal_values[name])::Float32
    elseif name === :i
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :s
        return (obj.__protobuf_jl_internal_values[name])::Array{UInt8,1}
    elseif name === :t
        return (obj.__protobuf_jl_internal_values[name])::TensorProto
    elseif name === :g
        return (obj.__protobuf_jl_internal_values[name])::Any
    elseif name === :sparse_tensor
        return (obj.__protobuf_jl_internal_values[name])::SparseTensorProto
    elseif name === :floats
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Float32}
    elseif name === :ints
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Int64}
    elseif name === :strings
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Array{UInt8,1}}
    elseif name === :tensors
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{TensorProto}
    elseif name === :graphs
        return (obj.__protobuf_jl_internal_values[name])::Any
    elseif name === :sparse_tensors
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{SparseTensorProto}
    else
        getfield(obj, name)
    end
end

mutable struct ValueInfoProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}

    function ValueInfoProto(; kwargs...)
        obj = new(meta(ValueInfoProto), Dict{Symbol,Any}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct ValueInfoProto (has cyclic type dependency)
const __meta_ValueInfoProto = Ref{ProtoMeta}()
function meta(::Type{ValueInfoProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_ValueInfoProto)
            __meta_ValueInfoProto[] = target = ProtoMeta(ValueInfoProto)
            allflds = Pair{Symbol,Union{Type,String}}[:name => AbstractString, :_type => "TypeProto", :doc_string => AbstractString]
            meta(target, ValueInfoProto, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_ValueInfoProto[]
    end
end
function Base.getproperty(obj::ValueInfoProto, name::Symbol)
    if name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :_type
        return (obj.__protobuf_jl_internal_values[name])::Any
    elseif name === :doc_string
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct NodeProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}

    function NodeProto(; kwargs...)
        obj = new(meta(NodeProto), Dict{Symbol,Any}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct NodeProto (has cyclic type dependency)
const __meta_NodeProto = Ref{ProtoMeta}()
function meta(::Type{NodeProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_NodeProto)
            __meta_NodeProto[] = target = ProtoMeta(NodeProto)
            fnum = Int[1,2,3,4,7,5,6]
            allflds = Pair{Symbol,Union{Type,String}}[:input => Base.Vector{AbstractString}, :output => Base.Vector{AbstractString}, :name => AbstractString, :op_type => AbstractString, :domain => AbstractString, :attribute => Base.Vector{AttributeProto}, :doc_string => AbstractString]
            meta(target, NodeProto, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_NodeProto[]
    end
end
function Base.getproperty(obj::NodeProto, name::Symbol)
    if name === :input
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{AbstractString}
    elseif name === :output
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{AbstractString}
    elseif name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :op_type
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :domain
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :attribute
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{AttributeProto}
    elseif name === :doc_string
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct TrainingInfoProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}

    function TrainingInfoProto(; kwargs...)
        obj = new(meta(TrainingInfoProto), Dict{Symbol,Any}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct TrainingInfoProto (has cyclic type dependency)
const __meta_TrainingInfoProto = Ref{ProtoMeta}()
function meta(::Type{TrainingInfoProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_TrainingInfoProto)
            __meta_TrainingInfoProto[] = target = ProtoMeta(TrainingInfoProto)
            allflds = Pair{Symbol,Union{Type,String}}[:initialization => "GraphProto", :algorithm => "GraphProto", :initialization_binding => Base.Vector{StringStringEntryProto}, :update_binding => Base.Vector{StringStringEntryProto}]
            meta(target, TrainingInfoProto, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_TrainingInfoProto[]
    end
end
function Base.getproperty(obj::TrainingInfoProto, name::Symbol)
    if name === :initialization
        return (obj.__protobuf_jl_internal_values[name])::Any
    elseif name === :algorithm
        return (obj.__protobuf_jl_internal_values[name])::Any
    elseif name === :initialization_binding
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{StringStringEntryProto}
    elseif name === :update_binding
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{StringStringEntryProto}
    else
        getfield(obj, name)
    end
end

mutable struct ModelProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}

    function ModelProto(; kwargs...)
        obj = new(meta(ModelProto), Dict{Symbol,Any}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct ModelProto (has cyclic type dependency)
const __meta_ModelProto = Ref{ProtoMeta}()
function meta(::Type{ModelProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_ModelProto)
            __meta_ModelProto[] = target = ProtoMeta(ModelProto)
            fnum = Int[1,8,2,3,4,5,6,7,14,20]
            allflds = Pair{Symbol,Union{Type,String}}[:ir_version => Int64, :opset_import => Base.Vector{OperatorSetIdProto}, :producer_name => AbstractString, :producer_version => AbstractString, :domain => AbstractString, :model_version => Int64, :doc_string => AbstractString, :graph => "GraphProto", :metadata_props => Base.Vector{StringStringEntryProto}, :training_info => Base.Vector{TrainingInfoProto}]
            meta(target, ModelProto, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_ModelProto[]
    end
end
function Base.getproperty(obj::ModelProto, name::Symbol)
    if name === :ir_version
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :opset_import
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{OperatorSetIdProto}
    elseif name === :producer_name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :producer_version
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :domain
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :model_version
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :doc_string
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :graph
        return (obj.__protobuf_jl_internal_values[name])::Any
    elseif name === :metadata_props
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{StringStringEntryProto}
    elseif name === :training_info
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{TrainingInfoProto}
    else
        getfield(obj, name)
    end
end

mutable struct GraphProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}

    function GraphProto(; kwargs...)
        obj = new(meta(GraphProto), Dict{Symbol,Any}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct GraphProto (has cyclic type dependency)
const __meta_GraphProto = Ref{ProtoMeta}()
function meta(::Type{GraphProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_GraphProto)
            __meta_GraphProto[] = target = ProtoMeta(GraphProto)
            fnum = Int[1,2,5,15,10,11,12,13,14]
            allflds = Pair{Symbol,Union{Type,String}}[:node => Base.Vector{NodeProto}, :name => AbstractString, :initializer => Base.Vector{TensorProto}, :sparse_initializer => Base.Vector{SparseTensorProto}, :doc_string => AbstractString, :input => Base.Vector{ValueInfoProto}, :output => Base.Vector{ValueInfoProto}, :value_info => Base.Vector{ValueInfoProto}, :quantization_annotation => Base.Vector{TensorAnnotation}]
            meta(target, GraphProto, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_GraphProto[]
    end
end
function Base.getproperty(obj::GraphProto, name::Symbol)
    if name === :node
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{NodeProto}
    elseif name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :initializer
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{TensorProto}
    elseif name === :sparse_initializer
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{SparseTensorProto}
    elseif name === :doc_string
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :input
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{ValueInfoProto}
    elseif name === :output
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{ValueInfoProto}
    elseif name === :value_info
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{ValueInfoProto}
    elseif name === :quantization_annotation
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{TensorAnnotation}
    else
        getfield(obj, name)
    end
end

mutable struct TypeProto_Tensor <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}

    function TypeProto_Tensor(; kwargs...)
        obj = new(meta(TypeProto_Tensor), Dict{Symbol,Any}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end
const __meta_TypeProto_Tensor = Ref{ProtoMeta}()
function meta(::Type{TypeProto_Tensor})
    ProtoBuf.metalock() do
        if !isassigned(__meta_TypeProto_Tensor)
            __meta_TypeProto_Tensor[] = target = ProtoMeta(TypeProto_Tensor)
            allflds = Pair{Symbol,Union{Type,String}}[:elem_type => Int32, :shape => TensorShapeProto]
            meta(target, TypeProto_Tensor, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_TypeProto_Tensor[]
    end
end
function Base.getproperty(obj::TypeProto_Tensor, name::Symbol)
    if name === :elem_type
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :shape
        return (obj.__protobuf_jl_internal_values[name])::TensorShapeProto
    else
        getfield(obj, name)
    end
end

mutable struct TypeProto_Sequence <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}

    function TypeProto_Sequence(; kwargs...)
        obj = new(meta(TypeProto_Sequence), Dict{Symbol,Any}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct TypeProto_Sequence (has cyclic type dependency)
const __meta_TypeProto_Sequence = Ref{ProtoMeta}()
function meta(::Type{TypeProto_Sequence})
    ProtoBuf.metalock() do
        if !isassigned(__meta_TypeProto_Sequence)
            __meta_TypeProto_Sequence[] = target = ProtoMeta(TypeProto_Sequence)
            allflds = Pair{Symbol,Union{Type,String}}[:elem_type => "TypeProto"]
            meta(target, TypeProto_Sequence, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_TypeProto_Sequence[]
    end
end
function Base.getproperty(obj::TypeProto_Sequence, name::Symbol)
    if name === :elem_type
        return (obj.__protobuf_jl_internal_values[name])::Any
    else
        getfield(obj, name)
    end
end

mutable struct TypeProto_Map <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}

    function TypeProto_Map(; kwargs...)
        obj = new(meta(TypeProto_Map), Dict{Symbol,Any}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct TypeProto_Map (has cyclic type dependency)
const __meta_TypeProto_Map = Ref{ProtoMeta}()
function meta(::Type{TypeProto_Map})
    ProtoBuf.metalock() do
        if !isassigned(__meta_TypeProto_Map)
            __meta_TypeProto_Map[] = target = ProtoMeta(TypeProto_Map)
            allflds = Pair{Symbol,Union{Type,String}}[:key_type => Int32, :value_type => "TypeProto"]
            meta(target, TypeProto_Map, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_TypeProto_Map[]
    end
end
function Base.getproperty(obj::TypeProto_Map, name::Symbol)
    if name === :key_type
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :value_type
        return (obj.__protobuf_jl_internal_values[name])::Any
    else
        getfield(obj, name)
    end
end

mutable struct TypeProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}

    function TypeProto(; kwargs...)
        obj = new(meta(TypeProto), Dict{Symbol,Any}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct TypeProto (has cyclic type dependency)
const __meta_TypeProto = Ref{ProtoMeta}()
function meta(::Type{TypeProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_TypeProto)
            __meta_TypeProto[] = target = ProtoMeta(TypeProto)
            fnum = Int[1,4,5,6]
            allflds = Pair{Symbol,Union{Type,String}}[:tensor_type => TypeProto_Tensor, :sequence_type => TypeProto_Sequence, :map_type => TypeProto_Map, :denotation => AbstractString]
            oneofs = Int[1,1,1,1]
            oneof_names = Symbol[Symbol("value")]
            meta(target, TypeProto, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, oneofs, oneof_names)
        end
        __meta_TypeProto[]
    end
end
function Base.getproperty(obj::TypeProto, name::Symbol)
    if name === :tensor_type
        return (obj.__protobuf_jl_internal_values[name])::TypeProto_Tensor
    elseif name === :sequence_type
        return (obj.__protobuf_jl_internal_values[name])::TypeProto_Sequence
    elseif name === :map_type
        return (obj.__protobuf_jl_internal_values[name])::TypeProto_Map
    elseif name === :denotation
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

export Version, AttributeProto_AttributeType, AttributeProto, ValueInfoProto, NodeProto, TrainingInfoProto, ModelProto, StringStringEntryProto, TensorAnnotation, GraphProto, TensorProto_DataType, TensorProto_DataLocation, TensorProto_Segment, TensorProto, SparseTensorProto, TensorShapeProto_Dimension, TensorShapeProto, TypeProto_Tensor, TypeProto_Sequence, TypeProto_Map, TypeProto, OperatorSetIdProto, AttributeProto_AttributeType, AttributeProto, ValueInfoProto, NodeProto, TrainingInfoProto, ModelProto, GraphProto, TypeProto_Sequence, TypeProto_Map, TypeProto


text:
using ADCME
using LinearAlgebra
using MAT
using AdFem
using PyPlot
using SparseArrays

function Q_exact(x, y)
    # (2*x^2*y*(x - 1)*(2*y - 1) + 1.0*(y^2 + 1)^2*(x*y*(x - 1) + x*y*(y - 1) + x*(x - 1)*(y - 1) + y*(x - 1)*(y - 1)) - (y^2 + 1)*(2*x*(x - 1)*(x + (x^2 + 1)*(y^2 + 1)) + 2*y*(x + (x^2 + 1)*(y^2 + 1))*(y - 1) + y*(2*x - 1)*(y - 1)*(2*x*(y^2 + 1) + 1)))/(y^2 + 1)^2
    # (2*x^2*y*(x - 1)*(2*y - 1) + x*(y^2 + 1)^2*(-y^2*(2*x - 1)*(y - 1)^2 + 1.0*y*(x - 1) + 1.0*(x - 1)*(y - 1)) - (y^2 + 1)*(2*x*(x - 1)*(x + (x^2 + 1)*(y^2 + 1)) + 2*y*(x + (x^2 + 1)*(y^2 + 1))*(y - 1) + y*(2*x - 1)*(y - 1)*(2*x*(y^2 + 1) + 1)))/(y^2 + 1)^2
    (2*x^2*y*(x - 1)*(2*y - 1)*(x + y + 1) - x*y^2*(2*x - 1)*(y - 1)^2*(y^2 + 1)^2*(x + y + 1) + x*(x - 1)*(2*y - 1)*(y^2 + 1)^2 - (y^2 + 1)*(x + y + 1)*(2*x*(x - 1)*(x + (x^2 + 1)*(y^2 + 1)) + 2*y*(x + (x^2 + 1)*(y^2 + 1))*(y - 1) + y*(2*x - 1)*(y - 1)*(2*x*(y^2 + 1) + 1)))/((y^2 + 1)^2*(x + y + 1))
end

function T_exact(x, y)
    x * (1-x) * y * (1-y)
end

function u_exact(x, y)
    x * y * (1-y)
end

function v_exact(x, y)
    1 / (1 + x + y)
end

function k_exact(x, y)
    1 + x^2 + x / (1 + y^2)
end

function k_nn(x, y)
    out = fc([x y], [20,20,20,1])^2 + 0.1 # N x 1
    squeeze(out)
end

#---------------------------------------------
# grid setup
m = 20
n = 20
h = 1/n

#---------------------------------------------
xy = fem_nodes(m, n, h)
x, y = xy[:,1], xy[:,2]
T0 = @. T_exact(x, y)
u = @. u_exact(x,y)
v = @. v_exact(x,y)
Q = @. Q_exact(x,y)
k = @. k_nn(x, y); k=stack(k)
# k0 = @. k_nn(xy)

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

bd = bcnode("all", m, n, h)

ugauss = fem_to_gauss_points(u, m, n, h)
vgauss = fem_to_gauss_points(v, m, n, h)
kgauss = fem_to_gauss_points(k, m, n, h)
Qgauss = fem_to_gauss_points(Q, m, n, h)

Advection = constant(compute_fem_advection_matrix1(constant(ugauss), constant(vgauss), m, n, h))
Laplace = compute_fem_laplace_matrix1(kgauss, m, n, h)
A = Advection + Laplace
A, _ = fem_impose_Dirichlet_boundary_condition1(A, bd, m, n, h)
b = constant(compute_fem_source_term1(constant(Qgauss), m, n, h))
b = scatter_update(b, bd, zeros(length(bd)))
T_computed = A\b

T_data = matread("chtenergy_data.mat")["T"]
# loss = mean((T_computed .- T_data)^2)

sample_size = 40
idx = rand(1:(m+1)*(n+1), sample_size)
observed_data = T_data[idx]
loss = mean((T_computed[idx] .- observed_data)^2)

noise = false
noise_level = 0.01
if noise
    noise_ratio = (1 - noise_level) .+ 2 * noise_level * rand(Float64, size(observed_data)) # uniform on (1-noise_level, 1+noise_level)
    observed_data = observed_data .* noise_ratio
end

loss = loss * 1e10
# ---------------------------------------------------
# create a session and run
max_iter = 100
sess = Session(); init(sess)
loss_ = BFGS!(sess, loss, max_iter)
figure(); semilogy(loss_); savefig("chtenergy_loss.png")

figure(figsize=(14,4));
subplot(131)
visualize_scalar_on_fem_points(k_exact.(x,y), m, n, h); title("conductivity exact")
subplot(132)
visualize_scalar_on_fem_points(run(sess, k), m, n, h); title("conductivity prediction")
subplot(133)
visualize_scalar_on_fem_points(k_exact.(x,y).-run(sess, k), m, n, h); title("conductivity difference")
savefig("chtenergy_k.png")

ext: jl


text:
### TODO:
### 1. Test speeds without caching a, b, d; maybe switch
### 2. Separate ALF computation to a different module
### 3. Compute H directly inside D array
### 4. Test skipping all the complicated indexing tricks; use fancy indexing
### 5. Allow specifying the extent of recursion / iterating over ℓ matrices

### NOTES:
### 1. Caching coefficients provides a ~11% speedup at low ℓ, but that falls off
###    as ℓ increases, reaching ~4% around ℓ=512.
### 2. This is probably a much more significant advantage for ALFs.

export WignerMatrixCalculator, H!, d!, D!, Y!, ϵ, WignerHindex, _WignerHindex

using ..SphericalFunctions: complex_powers!
using Quaternionic: AbstractQuaternion, to_euler_phases!

include("indexing.jl")
include("calculator.jl")
include("Hrecursions.jl")

@inline ϵ(m) = (m <= 0 ? 1 : (isodd(m) ? -1 : 1))

"""Return sign of input, with sign(0)=1"""
@inline sign(m) = (m < 0 ? -1 : 1)

"""
    d!(d, wigner, expiβ)

Compute Wigner's d matrix dˡₘₚ,ₘ(β)

# Parameters

* expiβ : array_like
    Values of expi(i*β) on which to evaluate the d matrix.
* out : array_like, optional
    Array into which the d values should be written.  It should be an array of
    floats, with size `self.dsize`.  If not present, the array will be created.
    In either case, the array will also be returned.
* workspace : array_like, optional
    A working array like the one returned by Wigner.new_workspace().  If not
    present, this object's default workspace will be used.  Note that it is not
    safe to use the same workspace on multiple threads.

# Returns

* d : array
    This is a 1-dimensional array of floats; see below.

# See Also

* H : Compute a portion of the H matrix
* D : Compute the full Wigner 𝔇 matrix
* rotate : Avoid computing the full 𝔇 matrix and rotate modes directly
* evaluate : Avoid computing the full 𝔇 matrix and evaluate modes directly

# Notes

This function is the preferred method of computing the d matrix for large ell
values.  In particular, above ell≈32 standard formulas become completely
unusable because of numerical instabilities and overflow.  This function uses
stable recursion methods instead, and should be usable beyond ell≈1000.

The result is returned in a 1-dimensional array ordered as

    [
        d(ell, mp, m, β)
        for ell in range(ell_max+1)
        for mp in range(-min(ℓ, mp_max), min(ℓ, mp_max)+1)
        for m in range(-ell, ell+1)
    ]

"""
function d!(d, w::WignerMatrixCalculator, expiβ::Complex)
    ell_min = ℓₘᵢₙ(w)
    ell_max = ℓₘₐₓ(w)
    mp_max = m′ₘₐₓ(w)
    if mp_max < ell_max
        throw(DomainError("ℓₘₐₓ = $ℓₘₐₓ",
                          "Cannot compute full d matrix up to ℓₘₐₓ with m′ₘₐₓ only $(mp_max)"
        ))
    end

    H!(w, expiβ)

    for ℓ in ell_min:ell_max
        for m′ in -ℓ:ℓ
            for m in -ℓ:ℓ
                i_d = WignerDindex(ℓ, m′, m, ell_min)
                i_H = WignerHindex(ℓ, m′, m, mp_max)
                d[i_d] = ϵ(m′) * ϵ(-m) * w.Hwedge[i_H]
            end
        end
    end

    d
end

function d!(w::WignerMatrixCalculator, β::Real)
    d = zeros(T(w), Wignerdsize(w))
    d!(w, expiβ, out)
end

"""
    D!(𝔇, w, R)
    D!(w, R)

Compute Wigner's 𝔇 matrix

# Parameters

* 𝔇 : array_like, optional
    Array into which the 𝔇 values should be written.  It should be an array of
    complex, with size `self.Dsize`.  If not present, the array will be
    created.  In either case, the array will also be returned.
* workspace : optional
    A working array like the one returned by Wigner.new_workspace().  If not
    present, this object's default workspace will be used.  Note that it is not
    safe to use the same workspace on multiple threads.
* R : Quaternion
    Array to be interpreted as a quaternionic array (thus its final dimension
    must have size 4), representing the rotations on which the 𝔇 matrix will be
    evaluated.

# Returns

* D : array
    This is a 1-dimensional array of complex; see below.

# See Also

* H : Compute a portion of the H matrix
* d : Compute the full Wigner d matrix
* rotate : Avoid computing the full 𝔇 matrix and rotate modes directly
* evaluate : Avoid computing the full 𝔇 matrix and evaluate modes directly

# Notes

This function is the preferred method of computing the 𝔇 matrix for large ell
values.  In particular, above ell≈32 standard formulas become completely
unusable because of numerical instabilities and overflow.  This function uses
stable recursion methods instead, and should be usable beyond ell≈1000.

This function computes 𝔇ˡₘₚ,ₘ(R).  The result is returned in a 1-dimensional
array ordered as

    [
        𝔇(ell, mp, m, R)
        for ell in range(ell_max+1)
        for mp in range(-min(ℓ, mp_max), min(ℓ, mp_max)+1)
        for m in range(-ell, ell+1)
    ]

"""
function D!(𝔇, w::WignerMatrixCalculator, R::AbstractQuaternion)
    ell_min = ℓₘᵢₙ(w)
    ell_max = ℓₘₐₓ(w)
    mp_max = m′ₘₐₓ(w)
    if mp_max < ell_max
        throw(DomainError("ℓₘₐₓ = $(ell_max)",
                          "Cannot compute full d matrix up to ℓₘₐₓ with m′ₘₐₓ only $(mp_max)"
        ))
    end

    to_euler_phases!(w.z, R)
    H!(w, w.z[2])
    complex_powers!(w.zₐpowers, w.z[1])
    complex_powers!(w.zᵧpowers, w.z[3])

    # 𝔇ˡₘₚ,ₘ(R) = dˡₘₚ,ₘ(R) exp[iϕₐ(m-mp)+iϕₛ(m+mp)] = dˡₘₚ,ₘ(R) exp[i(ϕₛ+ϕₐ)m+i(ϕₛ-ϕₐ)mp]
    # exp[iϕₛ] = R̂ₛ = hat(R[0] + 1j * R[3]) = zp
    # exp[iϕₐ] = R̂ₐ = hat(R[2] + 1j * R[1]) = zm.conjugate()
    # exp[i(ϕₛ+ϕₐ)] = zp * zm.conjugate() = z[2] = zᵧ
    # exp[i(ϕₛ-ϕₐ)] = zp * zm = z[0] = zₐ
    for ell in ell_min:ell_max
        for mp in -ell:-1
            i_D = WignerDindex(ell, mp, -ell, ell_min)
            for m in -ell:-1
                i_H = WignerHindex(ell, mp, m, mp_max)
                𝔇[i_D] = ϵ(mp) * ϵ(-m) * w.Hwedge[i_H] * conj(w.zᵧpowers[-m+1]) * conj(w.zₐpowers[-mp+1])
                i_D += 1
            end
            for m in 0:ell
                i_H = WignerHindex(ell, mp, m, mp_max)
                𝔇[i_D] = ϵ(mp) * ϵ(-m) * w.Hwedge[i_H] * w.zᵧpowers[m+1] * conj(w.zₐpowers[-mp+1])
                i_D += 1
            end
        end
        for mp in 0:ell
            i_D = WignerDindex(ell, mp, -ell, ell_min)
            for m in -ell:-1
                i_H = WignerHindex(ell, mp, m, mp_max)
                𝔇[i_D] = ϵ(mp) * ϵ(-m) * w.Hwedge[i_H] * conj(w.zᵧpowers[-m+1]) * w.zₐpowers[mp+1]
                i_D += 1
            end
            for m in 0:ell
                i_H = WignerHindex(ell, mp, m, mp_max)
                𝔇[i_D] = ϵ(mp) * ϵ(-m) * w.Hwedge[i_H] * w.zᵧpowers[m+1] * w.zₐpowers[mp+1]
                i_D += 1
            end
        end
    end
    𝔇
end

function D!(w::WignerMatrixCalculator, R::AbstractQuaternion)
    𝔇 = zeros(Complex{T(w)}, WignerDsize(w))
    D!(𝔇, w, R)
end

@doc raw"""
    Y!(Y, wigner, s, R)
    Y!(wigner, s, R)

Evaluate (and write into `Y`, if present) the values of ``{}_{s}Y_{\ell,
m}(R)`` for the input value of `s`, for all ``(\ell, m)`` throughout the range
specified by `wigner`.  `R` is assumed to be a unit quaternion (which may be
`Rotor`, or simply a `Quaternion`).  If `R` does not have unit magnitude, the
output elements will be too large by a factor ``|R|^{\ell}``.  If `Y` is not
present, a new array will be created.

The spherical harmonics of spin weight ``s`` are related to Wigner's
``\mathfrak{D}`` matrix as
```math
\begin{aligned}
{}_{s}Y_{\ell, m}(R)
  &= (-1)^s \sqrt{\frac{2\ell+1}{4\pi}} \mathfrak{D}^{(\ell)}_{m, -s}(R) \\
  &= (-1)^s \sqrt{\frac{2\ell+1}{4\pi}} \bar{\mathfrak{D}}^{(\ell)}_{-s, m}(\bar{R}).
\end{aligned}
```
"""
function Y!(Y, w::WignerMatrixCalculator, s::Int, R::AbstractQuaternion)
    if length(Y) < Ysize(w)
        error("Input `Y` has length $(length(Y)); it should be at least $(Ysize(w))")
    end
    ell_min = ℓₘᵢₙ(w)
    ell_max = ℓₘₐₓ(w)
    mp_max = m′ₘₐₓ(w)
    if mp_max < abs(s)
        throw(DomainError("ℓₘₐₓ = $(ell_max)",
                          "Cannot compute sYlm for spin weight $s with m′ₘₐₓ only $(mp_max)"
        ))
    end

    to_euler_phases!(w.z, R)
    H!(w, w.z[2])
    complex_powers!(w.zₐpowers, w.z[1])
    complex_powers!(w.zᵧpowers, w.z[3])

    # Yˡₘₚ,ₘ(R) = dˡₘₚ,ₘ(R) exp[iϕₐ(m-mp)+iϕₛ(m+mp)] = dˡₘₚ,ₘ(R) exp[i(ϕₛ+ϕₐ)m+i(ϕₛ-ϕₐ)mp]
    # exp[iϕₛ] = R̂ₛ = hat(R[0] + 1j * R[3]) = zp
    # exp[iϕₐ] = R̂ₐ = hat(R[2] + 1j * R[1]) = zm.conjugate()
    # exp[i(ϕₛ+ϕₐ)] = zp * zm.conjugate() = z[2] = zᵧ
    # exp[i(ϕₛ-ϕₐ)] = zp * zm = z[0] = zₐ
    i_D = 1
    @inbounds for ell in ell_min:ell_max
        if ell < abs(s)
            for mp in -ell:ell
                Y[i_D] = 0
                i_D += 1
            end
        else
            factor = (isodd(s) ? -1 : 1) * √((2ell+1)/(4T(w)(π)))
            for mp in -ell:-1
                i_H = WignerHindex(ell, mp, -s, mp_max)
                if -s < 0
                    Y[i_D] = factor * ϵ(mp) * ϵ(s) * w.Hwedge[i_H] * conj(w.zᵧpowers[s+1]) * conj(w.zₐpowers[-mp+1])
                    # println((ell, mp, i_D, factor, ϵ(mp), ϵ(s), w.Hwedge[i_H], conj(w.zᵧpowers[s+1]), conj(w.zₐpowers[-mp+1])))
                else
                    Y[i_D] = factor * ϵ(mp) * ϵ(s) * w.Hwedge[i_H] * w.zᵧpowers[-s+1] * conj(w.zₐpowers[-mp+1])
                    # println((ell, mp, i_D, factor, ϵ(mp), ϵ(s), w.Hwedge[i_H], w.zᵧpowers[-s+1], conj(w.zₐpowers[-mp+1])))
                end
                i_D += 1
            end
            for mp in 0:ell
                i_H = WignerHindex(ell, mp, -s, mp_max)
                if -s < 0
                    Y[i_D] = factor * ϵ(mp) * ϵ(s) * w.Hwedge[i_H] * conj(w.zᵧpowers[s+1]) * w.zₐpowers[mp+1]
                    # println((ell, mp, i_D, factor, ϵ(mp), ϵ(s), w.Hwedge[i_H], conj(w.zᵧpowers[s+1]), w.zₐpowers[mp+1]))
                else
                    Y[i_D] = factor * ϵ(mp) * ϵ(s) * w.Hwedge[i_H] * w.zᵧpowers[-s+1] * w.zₐpowers[mp+1]
                    # println((ell, mp, i_D, factor, ϵ(mp), ϵ(s), w.Hwedge[i_H], w.zᵧpowers[-s+1], w.zₐpowers[mp+1]))
                end
                i_D += 1
            end
        end
    end
    Y
end

function Y!(w::WignerMatrixCalculator, s::Int, R::AbstractQuaternion)
    Y = zeros(Complex{T(w)}, Ysize(w))
    Y!(Y, w, s, R)
end


text:
"""Transit function to call the specified diffusion scheme."""
function diffusion!(    u::AbstractMatrix,
                        v::AbstractMatrix,
                        Diag::DiagnosticVars,
                        S::ModelSetup)

    if S.parameters.diffusion == "constant"
        diffusion_constant!(u,v,Diag,S)
    elseif S.parameters.diffusion == "Smagorinsky"
        diffusion_smagorinsky!(u,v,Diag,S)
    else
        throw(error("Diffusion scheme $(S.parameters.diffusion) is unsupported."))
    end
end

"""Biharmonic diffusion operator with constant viscosity coefficient
Viscosity = ν∇⁴ ⃗u. Although constant, the coefficient is actually inside
Viscosity = ∇⋅ν∇∇² ⃗u."""
function diffusion_constant!(   u::AbstractMatrix,
                                v::AbstractMatrix,
                                Diag::DiagnosticVars,
                                S::ModelSetup)

    stress_tensor!(u,v,Diag)
    viscous_tensor_constant!(Diag,S)

    @unpack LLu1,LLu2,LLv1,LLv2 = Diag.Smagorinsky
    @unpack S11,S12,S21,S22 = Diag.Smagorinsky

    ∂x!(LLu1,S11)
    ∂y!(LLu2,S12)
    ∂x!(LLv1,S21)
    ∂y!(LLv2,S22)
end

""" Smagorinsky-like biharmonic viscosity
Viscosity = ∇ ⋅ (cSmag Δ⁴ |D| ∇∇² ⃗u)
The Δ⁴-scaling is omitted as gradient operators are dimensionless."""
function diffusion_smagorinsky!(u::AbstractMatrix,
                                v::AbstractMatrix,
                                Diag::DiagnosticVars,
                                S::ModelSetup)

    @unpack dudx,dvdy,dvdx,dudy = Diag.Vorticity

    ∂x!(dudx,u)
    ∂y!(dvdy,v)
    ∂x!(dvdx,v)
    ∂y!(dudy,u)

    # biharmonic diffusion
    stress_tensor!(u,v,Diag)
    smagorinsky_coeff!(Diag,S)
    viscous_tensor_smagorinsky!(Diag)

    @unpack LLu1,LLu2,LLv1,LLv2 = Diag.Smagorinsky
    @unpack S11,S12,S21,S22 = Diag.Smagorinsky

    ∂x!(LLu1,S11)
    ∂y!(LLu2,S12)
    ∂x!(LLv1,S21)
    ∂y!(LLv2,S22)
end

"""νSmag = cSmag * |D|, where deformation rate |D| = √((∂u/∂x - ∂v/∂y)^2 + (∂u/∂y + ∂v/∂x)^2).
The grid spacing Δ is omitted here as the operators are dimensionless."""
function smagorinsky_coeff!(Diag::DiagnosticVars,S::ModelSetup)

    @unpack dudx,dudy,dvdx,dvdy = Diag.Vorticity
    @unpack DS,DS_q,DT,νSmag,νSmag_q = Diag.Smagorinsky
    @unpack ep = S.grid
    @unpack cSmag = S.constants

    # horizontal shearing strain squared
    m,n = size(DS_q)
    @boundscheck (m+ep,n) == size(dudy) || throw(BoundsError())
    @boundscheck (m,n) == size(dvdx) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            DS_q[i,j] = (dudy[i+ep,j] + dvdx[i,j])^2
        end
    end

    Ixy!(DS,DS_q)

    # horizontal tension squared
    m,n = size(DT)
    @boundscheck (m+ep,n+2) == size(dudx) || throw(BoundsError())
    @boundscheck (m+2,n) == size(dvdy) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            DT[i,j] = (dudx[i+ep,j+1] + dvdy[i+1,j])^2
        end
    end

    # viscosity = Smagorinsky coefficient times deformation rate
    m,n = size(νSmag)
    @boundscheck (m,n) == size(DS) || throw(BoundsError())
    @boundscheck (m,n) == size(DT) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            νSmag[i,j] = cSmag*sqrt(DS[i,j] + DT[i,j])
        end
    end

    Ixy!(νSmag_q,νSmag)
end

"""Biharmonic stress tensor ∇∇²(u,v) = (∂/∂x(∇²u), ∂/∂y(∇²u); ∂/∂x(∇²v), ∂/∂y(∇²v))"""
function stress_tensor!(u::AbstractMatrix,v::AbstractMatrix,Diag::DiagnosticVars)
    @unpack Lu,Lv,dLudx,dLudy,dLvdx,dLvdy = Diag.Laplace
    ∇²!(Lu,u)
    ∇²!(Lv,v)
    ∂x!(dLudx,Lu)
    ∂y!(dLudy,Lu)
    ∂x!(dLvdx,Lv)
    ∂y!(dLvdy,Lv)
end

"""Biharmonic stress tensor times Smagorinsky coefficient
νSmag * ∇∇² ⃗u = (S11, S12; S21, S22)."""
function viscous_tensor_smagorinsky!(Diag::DiagnosticVars)

    @unpack dLudx,dLudy,dLvdx,dLvdy = Diag.Laplace
    @unpack νSmag,νSmag_q,S11,S12,S21,S22 = Diag.Smagorinsky
    @unpack ep = Diag.Smagorinsky

    m,n = size(S11)
    @boundscheck (m+2-ep,n) == size(νSmag) || throw(BoundsError())
    @boundscheck (m,n) == size(dLudx) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            S11[i,j] = νSmag[i+1-ep,j] * dLudx[i,j]
        end
    end

    m,n = size(S12)
    @boundscheck (m,n) == size(νSmag_q) || throw(BoundsError())
    @boundscheck (m+ep,n) == size(dLudy) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            S12[i,j] = νSmag_q[i,j] * dLudy[i+ep,j]
        end
    end

    m,n = size(S21)
    @boundscheck (m,n) == size(νSmag_q) || throw(BoundsError())
    @boundscheck (m,n) == size(dLvdx) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            S21[i,j] = νSmag_q[i,j] * dLvdx[i,j]
        end
    end

    m,n = size(S22)
    @boundscheck (m,n+2) == size(νSmag) || throw(BoundsError())
    @boundscheck (m,n) == size(dLvdy) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            S22[i,j] = νSmag[i,j+1] * dLvdy[i,j]
        end
    end
end

"""Biharmonic stress tensor times constant viscosity coefficient
νB * ∇∇² ⃗u = (S11, S12; S21, S22)"""
function viscous_tensor_constant!(  Diag::DiagnosticVars,
                                    S::ModelSetup)

    @unpack dLudx,dLudy,dLvdx,dLvdy = Diag.Laplace
    @unpack S11,S12,S21,S22 = Diag.Smagorinsky
    @unpack ep = S.grid
    @unpack νB = S.constants

    m,n = size(S11)
    @boundscheck (m,n) == size(dLudx) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            S11[i,j] = νB * dLudx[i,j]
        end
    end

    m,n = size(S12)
    @boundscheck (m+ep,n) == size(dLudy) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            S12[i,j] = νB * dLudy[i+ep,j]
        end
    end

    m,n = size(S21)
    @boundscheck (m,n) == size(dLvdx) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            S21[i,j] = νB * dLvdx[i,j]
        end
    end

    m,n = size(S22)
    @boundscheck (m,n) == size(dLvdy) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            S22[i,j] = νB * dLvdy[i,j]
        end
    end
end

"""Update u with bottom friction tendency (Bu,Bv) and biharmonic viscosity."""
function add_drag_diff_tendencies!( u::Array{Tprog,2},
                                    v::Array{Tprog,2},
                                    Diag::DiagnosticVars{T,Tprog},
                                    S::ModelSetup{T,Tprog}) where {T,Tprog}

    @unpack Bu,Bv = Diag.Bottomdrag
    @unpack LLu1,LLu2,LLv1,LLv2 = Diag.Smagorinsky
    @unpack halo,ep,Δt_diff = S.grid

    m,n = size(u) .- (2*halo,2*halo)
    @boundscheck (m+2-ep,n+2) == size(Bu) || throw(BoundsError())
    @boundscheck (m,n+2) == size(LLu1) || throw(BoundsError())
    @boundscheck (m+2-ep,n) == size(LLu2) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            u[i+2,j+2] += Δt_diff*(Tprog(Bu[i+1-ep,j+1]) + Tprog(LLu1[i,j+1]) + Tprog(LLu2[i+1-ep,j]))
        end
    end

    m,n = size(v) .- (2*halo,2*halo)
    @boundscheck (m+2,n+2) == size(Bv) || throw(BoundsError())
    @boundscheck (m,n+2) == size(LLv1) || throw(BoundsError())
    @boundscheck (m+2,n) == size(LLv2) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
             v[i+2,j+2] += Δt_diff*(Tprog(Bv[i+1,j+1]) + Tprog(LLv1[i,j+1]) + Tprog(LLv2[i+1,j]))
        end
    end
end


text:
# Continuation of ex 7.13, Least Squares WRM

using NumericalMethodsforEngineers
using Test

@sym begin
  LagrangePolynomial(xi_, yi_) := Module([N, sum, num, den],
    begin
      N = Length(xi)
      sum = 0
      For( i=1, i <= N, Increment(i),
        begin
          num = 1
          den = 1
          For( j=1, j <= N, Increment(j),
            begin
              If(j != i, num *= (x-xi[j]))
              If(j != i, den *= (xi[i]-xi[j]))
            end
          )
          sum += yi[i] * num/den
        end
      ),
      Return(Simplify(sum))
    end
  )
end

@sym begin
  ClearAll(xi, yi, N1, Y, Ydotdot, C11, ytilde1, ytilde2)
  xi = [0, 1//2, 1]
  yi = [0, a, 1]
  Y(x_) := LagrangePolynomial(xi, yi)
  Ydotdot(x_) = D(Y(x), x, 2)
  R(x_) := Simplify(Expand(Ydotdot(x) - 3*x - 4*Y(x)))
  R(x_) = Simplify(Expand(R(x) ./ (a => -C1/4)))
  Rdot(x_) = D(R(x), C1)
  C11 = Solve(Integrate(R(x)*Rdot(x), [x, 0, 1]), C1)[1]
  SetJ(r, ToString(Simplify(R(x))))
  SetJ(C1, ToString(C11[1][2]))
  ytilde1 = Simplify(Y(x) ./ (a => -C1/4))
  ytilde2 = Simplify(ytilde1 ./ C11)
  SetJ(y, ToString(Simplify(Expand(ytilde2))));
end

println("\n\nExample 7.17: y'' = 3x + 4y, y(0)=0, y(1)=1")
println("by 1-point least squares Weighted Residual Method")
@sym Print("\nY(x) = $(Y(x))\n")
@sym Print("Y(x) = $(Simplify(Expand(Y(x) ./ (a => -C1/4))))\n")
@sym Println("R(x) = $(R(x))\n")
@sym Print("C1 = $(Solve(Integrate(R(x)*Rdot(x), [x, 0, 1]), C1)[1][1][2])\n")
@eval ytilde_1pt_leastsquares(x) = $(Meta.parse(y))
@eval rf_1pt_leastsquares(x, C1) = $(Meta.parse(r))
println()

rf_1pt_leastsquares_1(x) = rf_1pt_leastsquares(x, Meta.parse(C1))
println()

@test r == "4 + 2C1 + x + 4C1*x - 8x^2 - 4C1*x^2"
@test y == "(1/27)*x*(-8 + 35x)"
#@test (quadgk(rf_1pt_leastsquares_1, 0, 1))[1] < 5*eps()


text:
abstract type MLJTask <: MLJType end # `Task` already has meaning in Base

mutable struct UnsupervisedTask <: MLJTask
    X
    input_scitypes::NamedTuple
    input_scitype_union::DataType
    input_is_multivariate::Bool
end

function UnsupervisedTask(; data=nothing, ignore=Symbol[], verbosity=1)

    data != nothing || error("You must specify data=... ")

    if ignore isa Symbol
        ignore = [ignore, ]
    end

    names = schema(data).names |> collect

    issubset(ignore, names) || error("ignore=$ignore contains a column name "*
                                     "not encountered in supplied data.")

    features = filter(names) do ftr
        !(ftr in ignore)
    end

    if length(features) == 1
        input_is_multivariate = false
        X = selectcols(data, features[1])
    else
        input_is_multivariate = true
        X = selectcols(data, features)
    end

    if input_is_multivariate
        input_scitypes = scitypes(X)
        input_scitype_union = Union{input_scitypes...}
    else
        input_scitype_union = scitype_union(X)
        input_scitypes = (input=input_scitype_union,)
    end

    if Unknown <: input_scitype_union
        @warn "An input feature with unknown scitype has been encountered. "*
        "Check the input_scitypes field of your task. "

    end

    if verbosity > 0
        @info "input_scitype_union = $input_scitype_union"
    end

    return UnsupervisedTask(X, input_scitypes, input_scitype_union,
                            input_is_multivariate)
end

# X and y can be different views of a common object.
mutable struct SupervisedTask <: MLJTask
    X
    y
    is_probabilistic
    input_scitypes
    target
    input_scitype_union
    target_scitype_union
    input_is_multivariate::Bool
end

function SupervisedTask(X, y::AbstractVector;
                        is_probabilistic=false,
                        input_is_multivariate=true,
                        target=:target,
                        verbosity=1)

    # is_probabilistic != nothing ||
    #     error("You must specify is_probabilistic=true or is_probabilistic=false. ")

    target_scitype_union = scitype_union(y)

    if target_scitype_union <: Tuple
        types = target_scitype_union.parameters
        unknown_found = reduce(|, [Unknown <: t for t in types])
    else
        unknown_found = Unknown <: target_scitype_union
    end

    if input_is_multivariate
        input_scitypes = scitypes(X)
        input_scitype_union = Union{input_scitypes...}
    else
        input_scitype_union = scitype_union(X)
        input_scitypes = (input=input_scitype_union,)
    end

    unknown_found = unknown_found || (Unknown <: input_scitype_union)

    if unknown_found
        @warn "An Unknown scitype has been encountered. "*
        "Check the input_scitypes and target_scitype_union fields of "*
        "your task. "
    end

    if verbosity > 0
        @info "\nis_probabilistic = $is_probabilistic\ninput_scitype_union = "*
        "$input_scitype_union \ntarget_scitype_union = $target_scitype_union"
    end

    return SupervisedTask(X, y,
                          is_probabilistic,
                          input_scitypes,
                          target,
                          input_scitype_union, target_scitype_union,
                          input_is_multivariate)
end

function SupervisedTask(; data=nothing, is_probabilistic=false, target=nothing, ignore=Symbol[], verbosity=1)

    target_as_input = target

    data != nothing ||
        error("You must specify data=... or use SupervisedTask(X, y, ...).")

    is_probabilistic != nothing ||
        error("You must specify is_probabilistic=true or is_proabilistic=false. ")

    if ignore isa Symbol
        ignore = [ignore, ]
    end

    names = schema(data).names |> collect

    issubset(ignore, names) ||
        error("ignore=$ignore contains a column name not encountered "*
              "in supplied data.")

    target != nothing ||
        error("You must specify target=... (use Symbol or Vector{Symbol}) ")

    if target isa Vector
        issubset(target, names) ||
            error("One or more specified target columns missing from "*
                  "supplied data. ")
        yrowtable = Tables.rowtable(selectcols(data, target))
        y = values.(yrowtable)
    else
        target in schema(data).names ||
            error("Specificed target not found in data. ")
        y = selectcols(data, target)
        target = [target, ]
    end

    features = filter(names) do ftr
        !(ftr in target) && !(ftr in ignore)
    end

    if length(features) == 1
        X = selectcols(data, features[1])
        input_is_multivariate = false
    else
        X = selectcols(data, features)
        input_is_multivariate = true
    end

    return SupervisedTask(X, y;
                          is_probabilistic=is_probabilistic,
                          input_is_multivariate=input_is_multivariate,
                          target=target_as_input,
                          verbosity=verbosity)
end

## RUDIMENTARY TASK OPERATIONS

nrows(task::MLJTask) = nrows(task.X)
Base.eachindex(task::MLJTask) = Base.OneTo(nrows(task))
nfeatures(task::MLJTask) = length(schema(task.X).names)

X_(task::MLJTask) = task.X
y_(task::SupervisedTask) = task.y

X_and_y(task::SupervisedTask) = (X_(task), y_(task))

# make tasks callable:
(task::UnsupervisedTask)() = X_(task)
(task::SupervisedTask)() = X_and_y(task)


text:
export lsqr

"""
x,flag,err,iter,resvec = lsqr(A,b,;x=[],maxIter=20,atol=1e-10,btol=1e-10,condlim=1e6,out=1,doBidiag=false,
                              storeInterm::Bool=false)

Least Squares QR method (LSQR) for solving

   min_x || A x - b ||,

where A is assumed is either over or underdetermined.

Implementation follows and comments refer to Algorithm LSQR in:

Paige, C. C., & Saunders, M. A. (1982).
LSQR: An Algorithm for Sparse Linear Equations and Sparse Least Squares.
ACM Transactions on Mathematical Software (TOMS), 8(1), 43–71. doi:10.1145/355984.355989

Required input:

  A       - function in (flag,v,alpha,x).
            if flag=='F'
                 x <- A*v + alpha*x (overwrite x)
            elseif flag=='T'
                 x <- A'*v + alpha*x (overwriting x)
            end
  b       - right hand side vector

Optional input:

  x           - starting guess (optional)
  maxIter     - maximum number of iterations
  atol        - error tolerance for (estimated) residual norm (for compatible systems)
  btol        - error tolerance for (estimated) gradient norm (for incompatible systems)
  condlim     - for stopping based on (estimated) condition number
  out         - flag for output (-1: no output, 0: only errors, 1: final status, 2: residual norm at each iteration)
  doBidiag    - returns the bidiagonalization of A at the final iteration
  storeInterm - store and return intermediate iterates

Output:

 x       - approximate solution
 flag    - exit flag ( 0 : stopping based on atol or btol,
                      -1 : maxIter reached without converging
                      -2 : condition number became too large
                      -9 : right hand side was zero)
 his     - status at each iteration
 U,S,V   - bidiagonalization (only build if doBidiag==true)
"""
function lsqr(A::SparseMatrixCSC,b::Vector;kwargs...)
    # build function that is compatible with lsqr
    ATv = zeros(size(A,2)) # output for A'*x

    # Af(flag,x,a=0.0,v=ATv) = (flag=='F') ? A_mul_B!(1.0,A,x,a,v) : At_mul_B!(1.0,A,x,a,v)
    Af(flag,x,a=0.0,v=ATv) = (flag=='F') ? mul!(v,A,x,1.0,a) : mul!(v,transpose(A),x,1.0,a)

    return lsqr(Af,b;kwargs...)

end

function lsqr(A::Any,b::Vector;kwargs...)
    # build function that is compatible with lsqr
    Af(flag,x,a=0.0,v=0.0) = (flag=='F') ? A*x.+a*v : A'*x.+a*v

    return lsqr(Af,b;kwargs...)

end

function lsqr(A::Function,b::Vector;x=[],maxIter=20,atol=1e-10,btol=1e-10,condlim=1e6,out=0,doBidiag=false,storeInterm::Bool=false)
    m  = length(b)

    nres0  = BLAS.nrm2(m,b,1) # norm(b)
    if nres0==0; return zero(eltype(b)),-9,[0.0]; end

    # (1) initialize
    if !isempty(x)
        b = copy(b)
        b = A('F',-x,1.0,b)
    end

    nrmb  = BLAS.nrm2(m,b,1) # norm(b)
    beta  = nrmb

    u     = b/beta
    v     = A('T',u,0.0)# A'*u
    n     = length(v)
	x     = (isempty(x)) ? zeros(n) : x
	if storeInterm;   X = zeros(n,maxIter) end
	if (beta <= btol*nres0); return x,1; end
	alpha = BLAS.nrm2(n,v,1)
	if alpha==0; return x,2; end
    v    /= alpha
    w     = copy(v)

    phib  = beta
    rhob  = alpha

    # (2) perform iteration
    # maxIter = min(maxIter,n,m)
    nBk  = 0.0
    nDk  = 0.0
    nres = beta
    STOP = fill(false,3)
    if doBidiag
        U    = zeros(m,maxIter+1)
        U[:,1] = u
        V    = zeros(n,maxIter)
        bs   = zeros(maxIter)
        as   = zeros(maxIter)
    end
    his  = zeros(maxIter,4)
    if out>1
        @printf "iter\t|r_k|\t\t|A'*r_k|\t|A|\t\tcond(A)\t\tSTOP\n"
    end
    flag = -1
    iter = 1
    while iter<=maxIter
        # (3) continue the bidiagonalization
        u     = A('F',v,-alpha,u) # A*v - alpha*u
        beta  = BLAS.nrm2(m,u,1)
        BLAS.scal!(m,1/beta,u,1) # faster than: u    /= beta
        if doBidiag
            as[iter] = alpha
            bs[iter] = beta
            U[:,iter+1] = u
            V[:,iter]   = v
        end
        nBk  += alpha^2+beta^2
        v     = A('T',u,-beta,v) # A'*u - beta*v
        alpha = BLAS.nrm2(n,v,1)
        BLAS.scal!(n,1/alpha,v,1) # faster than: v    /= alpha

        # (4) construct and apply next orthogonal transformation
        rho   = sqrt(rhob^2+beta^2)
        c     = rhob/rho
        s     = beta/rho
        theta = s*alpha
        rhob  = -c*alpha
        phi   = c*phib
        phib *= s
        nDk  += dot(w,w)/(rho*rho)

        # (5) update x,w
        BLAS.axpy!(n,phi/rho,w,1,x,1) # x += (phi/rho)*w
        if storeInterm; X[:,iter] = x; end

        # the following two lines are equivalent to w  = v - (theta/rho)*w
        w = BLAS.scal!(n,-theta/rho,w,1)
        w = BLAS.axpy!(n,1.0,v,1,w,1)

        # (6) check for stopping
        nrmA = sqrt(nBk)
        nrmG = abs(phib*alpha*c)
        conA = sqrt(nBk*nDk)
        nres *= s
        his[iter,:] = [nres nrmG nrmA conA]

        STOP[1] = nres <= btol*nres0+ atol*nrmA*norm(x)
        STOP[2] = (nrmG/(nrmA*nres)) <= atol
        STOP[3] = conA > condlim

        if out>1
            @printf "%03d\t%1.3e\t%1.3e\t%1.3e\t%1.3e\t[%d,%d,%d]\n" iter nres nrmG nrmA conA 1.0*STOP[1] 1.0*STOP[2] 1.0*STOP[3]
        end
        iter += 1

        if any(STOP)
            break;
        end
    end

    if STOP[1]; flag = 1; end
    if STOP[2] && !STOP[1]; flag = 2; end
    if STOP[3]; flag = -2; end

    if out>=0
        if flag==-1
            println(@sprintf("lsqr iterated maxIter (=%d) times without converging to desired tolerance. Estimated residual norm is %1.3e and estimated norm of gradient is %1.3e.",iter-1,his[iter-1,1],his[iter-1,2]))
        elseif flag==-2
            println(@sprintf("lsqr stopped at iteration %d because the matrix seems to be too ill-conditioned. Estimated condition number is %1.2e.",iter-1,his[iter-1,4]))
        elseif flag==1 && out>=1
            println(@sprintf("lsqr achieved desired tolerance at iteration %d. The system should be compatible. Residual norm is %1.2e.",iter-1,his[iter-1,1]))
        elseif flag==2 && out>=1
            println(@sprintf("lsqr achieved desired tolerance at iteration %d. The system should be incompatible. Gradient norm is %1.2e.",iter-1,his[iter-1,1]))
        end
    end

    x = storeInterm ? X[:,1:iter-1] : x;

    if doBidiag
        U = U[:,1:iter]
        V = V[:,1:iter-1]
		I,J,V = SparseArrays.spdiagm_internal(-1=>bs[1:iter-1],0=>as[1:iter-1])
		B = sparse(I,J,V,iter,iter-1)
        # B = spdiagm((bs[1:iter-1],as[1:iter-1]),-1:0,iter,iter-1)
        return x, flag,his[1:iter-1,:],U,B,V
    else
        return x, flag,his[1:iter-1,:]
    end
end


text:
mutable struct Range <: Reference
    file::Union{File,Nothing}
    from::Union{Int,Nothing}
    to::Union{Int,Nothing}
    function Range(f::File,from::Int,to::Int)
        @assert from <= to
        @assert length(f.data) >= to
        r = new(f,from,to)
        push!(f.references,r)
        return r
    end
end

@post function resolveRange(r::Range)
    if !isdestroyed(r) && !((r.from <= r.to))
        deletefromreferences!(r.file,r)
        destroy!(r)
    end
end

iterate(r::Range) = (r.from,r.from)

function iterate(r::Range,state::Int)
    if state < r.to
        return (state+1,state+1)
    else
        return nothing
    end
end

eltype(r::Range) = Int

lastindex(r::Range) = length(r)

function getindex(r::Range,idx::UnitRange{Int})
    first(idx) < 1 && throw(BoundsError(r,idx))
    last(idx) > length(r) && throw(BoundsError(r,idx))
    Range(r.file,r.from-1+first(idx),r.from-1+last(idx))
end

function getindex(r::Range,idx::Int)
    Line(r.file,r.from-1+first(idx))
end

@alive r value(r::Range) = r.file.data[r.from:r.to]
@alive r function value(r::Range,idx::Int)
    idx > length(r) && throw(BoundsError(r,idx))
    return r.file.data[r.from+idx-1]
end

@alive r length(r::Range) = r.to-r.from+1

==(a::Range,b::Range) = a.from == b.from && a.to == b.to && a.file === b.file

function show(io::IO,r::Range)
    if r.file == nothing
        print(io,"*Deleted*")
    else
        print(io,"Range($(r.from):$(r.to)) #$(truncate(value(r,1),12))▿$(truncate(value(r,length(r)),12))")
    end
end

@alive r next(r::Range,i::Int=1) = Line(r.file,r.to+i)
@alive r prev(r::Range,i::Int=1) = Line(r.file,r.from-i)

+(l::Range,i::Int) = next(l,i)
-(l::Range,i::Int) = prev(l,i)
+(i::Int,l::Range) = next(l,i)
-(i::Int,l::Range) = prev(l,i)

@alive r copy(r::Range) = Range(r.file,r.from,r.to)

#=
Delete functions
=#
@alive r function delete!(r::Range)
    notify_delete(r.file,r)
    deletefromreferences!(r.file,r)
    deleteat!(r.file.data,r.from:r.to)
    destroy!(r)
    return nothing
end

@alive r function destroy!(r::Range)
    r.from = nothing
    r.to = nothing
    r.file = nothing
end

isdestroyed(r::Range) = r.from == nothing && r.to == nothing && r.file == nothing

@alive sub sup function in(sub::Range,sup::Range)
    sub.from >= sup.from && sub.to <= sup.to
end
@alive sub sup function in(sub::Line,sup::Range)
    sub.ln >= sup.from && sub.ln <= sup.to
end

@resolveRange self function notify_delete(self::Range,del::Range)
    del_len = length(del)
    if del.to < self.from
        self.from -= del_len
        self.to -= del_len
    elseif del in self
        self.to -= del_len
    elseif !(del.from > self.to)
        deletefromreferences!(self.file,self)
        destroy!(self)
    end
end

@resolveRange self function notify_delete(self::Range,del::Line)
    del_len = 1
    if del.ln < self.from
        self.from -= del_len
        self.to -= del_len
    elseif del in self
        self.to -= del_len
    elseif !(del.ln > self.to)
        deletefromreferences!(self.file,self)
        destroy!(self)
    end
end

function notify_delete(self::Line,del::Range)
    del_len = length(del)
    if self.ln > del.to
        self.ln -= del_len
    elseif self in del
        deletefromreferences!(self.file,self)
        destroy!(self)
    end
end

#=
Insert functions
=#

@alive r function insert!(r::Range,str::Vector{String})
    notify_insert(r.file,r,length(str))
    splice!(r.file.data,r.from:r.from-1,str)
    r.from += length(str)
    r.to += length(str)
    if length(str) == 1
        return Line(r.file,r.from-1)
    else
        return Range(r.file,r.from-length(str),r.from-1)
    end
end

@alive r insert!(r::Range,str::String) = insert!(r::Range,[str,])
@alive r from insert!(r::Range,from::T) where T<:Reference = insert!(r,value(from))

@alive l function insert!(l::Line,str::Vector{String})
    notify_insert(l.file,l,length(str))
    splice!(l.file.data,l.ln:l.ln-1,str)
    l.ln += length(str)
    if length(str) == 1
        return Line(l.file,l.ln-1)
    else
        return Range(l.file,l.ln-length(str),l.ln-1)
    end
end

function notify_insert(self::Range,ins::Int,lines::Int)
    if self.from >= ins
        self.from += lines
    end
    if self.to >= ins
        self.to += lines
    end
end

notify_insert(self::Range,ins::Range,lines::Int) = notify_insert(self,ins.from,lines)
notify_insert(self::Range,ins::Line,lines::Int) = notify_insert(self,ins.ln,lines)

function notify_insert(self::Line,ins::Range,lines::Int)
    if self.ln >= ins.from
        self.ln += lines
    end
end

#=
Append functions
=#

@alive r function append!(r::Range,str::Vector{String})
    notify_append(r.file,r,length(str))
    splice!(r.file.data,r.to+1:r.to,str)
    if length(str) == 1
        return Line(r.file,r.to+1)
    else
        return Range(r.file,r.to+1,r.to+length(str))
    end
end

@alive l function append!(l::Line,str::Vector{String})
    notify_append(l.file,l,length(str))
    if l.ln < length(l.file.data)
        splice!(l.file.data,l.ln+1:l.ln,str)
    else
        append!(l.file.data,str)
    end
    if length(str) == 1
        return Line(l.file,l.ln+1)
    else
        return Range(l.file,l.ln+1,l.ln+length(str))
    end
end

append!(r::Range,str::String) = append!(r::Range,[str,])
@alive r from append!(r::Range,from::T) where T<:Reference = append!(r,value(from))

notify_append(self::Range,ins::Range,lines::Int) = notify_insert(self,ins.from,lines)
notify_append(self::Range,ins::Line,lines::Int) = notify_insert(self,ins.ln,lines)

function notify_append(self::Line,ins::Range,lines::Int)
    if self.ln > ins.to
        self.ln += lines
    end
end

#=
Replace functions
=#

@alive r function replace!(r::Range,str::Vector{String})
    lines = length(r)-length(str)
    if lines == 0
        r.file.data[r.from:r.to] = str
    else
        splice!(r.file.data,r.from:r.to,str)
        notify_replace(r.file,r,lines)
    end
    if length(str) == 1
        ln = r.from
        file = r.file
        deletefromreferences!(file,r)
        destroy!(r)
        return Line(file,ln)
    end
    r.to = r.from + length(str)-1
    return r
end

@alive r replace!(r::Range,str::String) = replace!(r,[str,])

@alive l function replace!(l::Line,str::Vector{String})
    lines = 1-length(str)
    if lines == 0
        l.file.data[l.ln] = str[1]
    else
        splice!(l.file.data,l.ln,str)
        notify_replace(l.file,l,lines)
    end
    if length(str) > 1
        from = l.ln
        file = l.file
        deletefromreferences!(file,l)
        destroy!(l)
        return Range(file,from,from+length(str)-1)
    end
    return l
end

@resolveRange self function notify_replace(self::Range,rep::Range,lines::Int)
    if self.from >= rep.to
        self.from -= lines
        self.to -= lines
    elseif rep in self
        self.to -= lines
    elseif !(rep.from > self.to)
        deletefromreferences!(self.file,self)
        destroy!(self)
    end
end

@resolveRange self function notify_replace(self::Range,rep::Line,lines::Int)
    if rep.ln < self.from
        self.from -= lines
        self.to -= lines
    elseif rep in self
        self.to -= lines
    elseif !(rep.ln > self.to)
        deletefromreferences!(self.file,self)
        destroy!(self)
    end
end

function notify_replace(self::Line,rep::Range,lines::Int)
    if self.ln > rep.to
        self.ln -= lines
    elseif self in rep
        deletefromreferences!(self.file,self)
        destroy!(self)
    end
end

replace!(r::Range,from::T) where T<:Reference = replace!(r,value(from))

###############################################################################
#                                     Move                                    #
###############################################################################

movebeforeline(ref::Range) = ref.from
moveafterline(ref::Range) = ref.to

@alive self to function move!(self::Range,to::T) where T<:Reference
    # Insert duplicate
    dup = insert!(to,self)
    # Delete original but keep reference
    notify_delete(self.file,self)
    deleteat!(self.file.data,self.from:self.to)
    # Move Range pointer
    self.from = dup.from
    self.to = dup.to
    # Delete reference for dup
    deletefromreferences!(dup.file,dup)
    return self
end

@alive self to function moveafter!(self::Range,to::T) where T<:Reference
    # Insert duplicate
    dup = append!(to,self)
    # Delete original but keep reference
    notify_delete(self.file,self)
    deleteat!(self.file.data,self.from:self.to)
    # Move Range pointer
    self.from = dup.from
    self.to = dup.to
    # Delete reference for dup
    deletefromreferences!(dup.file,dup)
    return self
end

###############################################################################
#                                Misc functions                               #
###############################################################################

@alive r function map!(fun::Function,r::Range)
    for i in r
        r.file.data[i] = fun(r.file.data[i])
    end
end

@alive r map(fun::Function,r::Range) = map(fun,r.file.data[r.from:r.to])


text:
function mlp(; layers=(), last=softmax, actf=relu, winit=Gaussian(0,.01), binit=Constant(0))
    prog = quote
        x0 = input()
    end
    N = length(layers)
    x0 = x1 = :x0
    for n=1:N
        x0 = x1; x1 = symbol("x$n")
        op = :($x1 = mlplayer($x0; n=$(layers[n]), f=$(n<N ? actf : last), winit=$winit, binit=$binit))
        push!(prog.args, op)
    end
    return prog
end

mlplayer(; n=1, f=relu, winit=Gaussian(0,.01), binit=Constant(0)) =
quote
    x1 = input()
    w1 = par($n,0; init=$winit)
    x2 = dot(w1,x1)
    b2 = par(0; init=$binit)
    x3 = add(b2,x2)
    y3 = $f(x3)
end


text:
module Print

pluralize(n::Real, name, plural="s", singular="") = "$n $(name)$(n != 1 ? plural : singular)"

function and_join(strs::AbstractVector)
    if length(strs) == 0
        "nothing"
    elseif length(strs) == 1
        strs[1]
    elseif length(strs) == 2
        join(strs, " and ")
    else
        join(strs[1:end-1], ", ") * ", and " * last(strs)
    end
end

function algs(df)
    algs = union(df.algs...)
    unstable = union(df.unstable...)
    and_join([alg ∈ unstable ? "$alg (unstable)" : "$alg" for alg in algs])
end

function domain(df)
    and_join(collect(pluralize(length(unique(a)), n) * z for (a,n,z) in [
        (df.ContainerType, "container type", ""),
        (df.Type, "element type", ""),
        (df.len, "length", " up to $(maximum(df.len))"),
        (df.order, "order", ""),
        (df.source, "distribution", ""),
    ]))
end

function silly_split(str::AbstractString, width=last(displaysize(stdout)))
    last_break=0
    v = collect(str)
    while true
        i = last_break+width
        if i > length(v)
            return typeof(str)(v)
        end
        while true
            if i <= last_break
                return typeof(str)(v)
            end
            if v[i] == ' '
                break
            end
            i -= 1
        end
        v[i] = '\n'
        last_break = i
    end
end

function _print(df, passed_tests, time, benchmark_time, color)
    printstyled(silly_split("$(algs(df)) passed $(pluralize(passed_tests, "test")) spanning"*
    " $(domain(df)) in $(round(Integer, time))s. $(round(Integer, benchmark_time))s"*
    " ($(round(Integer, benchmark_time/time*100))%) spent in benchmarks.\n"), color=color)
end

end

max_stars_repo_name: LilithHafner/SortMark

max_issues_repo_name: LilithHafner/SortMark


text:
"""
    get_header(table::Table)

Return the header string of the `Table` including `======` as the second line
"""
function get_header(table::Table)
    setup = table.setup
    ln = ""
    sum_width = 0
    for (name, width) in zip(setup.names, setup.widths)
        sum_width += width
        padding = width - 2 - length(name)
        if padding < 0
            error("The width must be at least 2+length(name) for each column")
        end
        ln *= repeat(" ", fld(padding, 2) + 1)
        ln *= name
        ln *= repeat(" ", cld(padding, 2) + 1)
    end
    equals = repeat("=", sum(sum_width))
    header = "$ln\n$equals"
    return header
end

"""
    print_header(table::Table)

Print the header of the given `table`. Calls [`get_header`](@ref).
"""
print_header(table::Table) = println(get_header(table))


text:
# Function and operator calls.

mutable struct FunctionNode <: AbstractSQLNode
    name::Symbol
    args::Vector{SQLNode}

    FunctionNode(;
             name::Union{Symbol, AbstractString},
             args = SQLNode[]) =
        new(Symbol(name), args)
end

FunctionNode(name; args = SQLNode[]) =
    FunctionNode(name = name, args = args)

FunctionNode(name, args...) =
    FunctionNode(name = name, args = SQLNode[args...])

"""
    Fun(; name, args = [])
    Fun(name; args = [])
    Fun(name, args...)

A function or an operator invocation.

# Example

```jldoctest
julia> person = SQLTable(:person, columns = [:person_id, :year_of_birth]);

julia> q = From(person) |>
           Where(Fun.not(Get.person_id .> 2000));

julia> print(render(q))
SELECT "person_1"."person_id", "person_1"."year_of_birth"
FROM "person" AS "person_1"
WHERE (NOT ("person_1"."person_id" > 2000))
```
"""
Fun(args...; kws...) =
    FunctionNode(args...; kws...) |> SQLNode

dissect(scr::Symbol, ::typeof(Fun), pats::Vector{Any}) =
    dissect(scr, FunctionNode, pats)

PrettyPrinting.quoteof(n::FunctionNode, ctx::QuoteContext) =
    Expr(:call,
         Expr(:., nameof(Fun),
                  QuoteNode(Base.isidentifier(n.name) ? n.name : string(n.name))),
         quoteof(n.args, ctx)...)

label(n::FunctionNode) =
    n.name

# Notation for making function nodes.

struct FunClosure
    name::Symbol
end

Base.show(io::IO, f::FunClosure) =
    print(io, Expr(:., nameof(Fun),
                       QuoteNode(Base.isidentifier(f.name) ? f.name : string(f.name))))

Base.getproperty(::typeof(Fun), name::Symbol) =
    FunClosure(name)

Base.getproperty(::typeof(Fun), name::AbstractString) =
    FunClosure(Symbol(name))

(f::FunClosure)(args...) =
    Fun(f.name, args = SQLNode[args...])

(f::FunClosure)(; args = SQLNode[]) =
    Fun(f.name, args = args)

# Broadcasting notation.

struct FunStyle <: Base.BroadcastStyle
end

Base.BroadcastStyle(::Type{<:AbstractSQLNode}) =
    FunStyle()

Base.BroadcastStyle(::FunStyle, ::Base.Broadcast.DefaultArrayStyle{0}) =
    FunStyle()

Base.broadcastable(n::AbstractSQLNode) =
    n

Base.Broadcast.instantiate(bc::Base.Broadcast.Broadcasted{FunStyle}) =
    bc

Base.copy(bc::Base.Broadcast.Broadcasted{FunStyle}) =
    Fun(nameof(bc.f), args = SQLNode[bc.args...])

Base.convert(::Type{AbstractSQLNode}, bc::Base.Broadcast.Broadcasted{FunStyle}) =
    FunctionNode(nameof(bc.f), args = SQLNode[bc.args...])


text:
# ---
# title: 785. Is Graph Bipartite?
# id: problem785
# author: Indigo
# date: 2021-01-28
# difficulty: Medium
# categories: Depth-first Search, Breadth-first Search, Graph
# link: <https://leetcode.com/problems/is-graph-bipartite/description/>
# hidden: true
# ---
#
# Given an undirected `graph`, return `true` if and only if it is bipartite.
#
# Recall that a graph is _bipartite_ if we can split its set of nodes into two
# independent subsets A and B, such that every edge in the graph has one node in
# A and another node in B.
#
# The graph is given in the following form: `graph[i]` is a list of indexes `j`
# for which the edge between nodes `i` and `j` exists.  Each node is an integer
# between `0` and `graph.length - 1`.  There are no self edges or parallel
# edges: `graph[i]` does not contain `i`, and it doesn't contain any element
# twice.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/21/bi1.jpg)
#
#
#
#     Input: graph = [[1,3],[0,2],[1,3],[0,2]]
#     Output: true
#     Explanation: We can divide the vertices into two groups: {0, 2} and {1, 3}.
#
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/21/bi2.jpg)
#
#
#
#     Input: graph = [[1,2,3],[0,2],[0,1,3],[0,2]]
#     Output: false
#     Explanation: We cannot find a way to divide the set of nodes into two independent subsets.
#
#
#
#
#
# **Constraints:**
#
#   * `1 <= graph.length <= 100`
#   * `0 <= graph[i].length < 100`
#   * `0 <= graph[i][j] <= graph.length - 1`
#   * `graph[i][j] != i`
#   * All the values of `graph[i]` are **unique**.
#   * The graph is **guaranteed** to be **undirected**.
#
#
## @lc code=start
using LeetCode

function is_bipartite(graph::Vector{Vector{Int}})
    for edge in graph
        edge .+= 1
    end
    color = fill(0, length(graph))
    q = Queue{Int}()
    for i in 1:length(graph)
        if color[i] != 0
            continue
        end
        color[i] = 1
        enqueue!(q, i)
        while !isempty(q)
            root = dequeue!(q)
            for neighbor in graph[root]
                if color[neighbor] == 0
                    color[neighbor] = -color[root]
                    enqueue!(q, neighbor)
                elseif color[neighbor] == color[root]
                    return false
                end
            end
        end
    end
    return true
end
## @lc code=end


text:
using NCDatasets
using Formatting
using PyPlot

if length(ARGS) >= 1
    fn = ARGS[1]
end

const Re = 6371000.0
const cp = 1004.0

ds = Dataset(fn, "r")

lat = ds["lat"][:]
lon = ds["lon"][:]
lev = ds["lev"][:]
ilev = ds["ilev"][:]

deg2rad = π/180.0

cross_section = (lon[2] - lon[1]) * deg2rad * Re * cos.(lat * deg2rad)

mass = ((ilev[2:end] - ilev[1:end-1]) / 9.8 * 100.0)'
mass = repeat(mass, outer=(length(lat), ))

mht = view(ds["VT"][:], :, :, :, 1)
mht = view(sum(mht, dims=(1,)), 1, :, :)       # take sum over lon
mht = view(sum(mht .* mass, dims=(2,)), :, 1)  # weight with mass
mht .*= cross_section * cp
close(ds)

mht ./= 1e15

fig, ax = plt[:subplots](1, 1, figsize=(6,4))

ax[:plot]([-90, 90], [0, 0], linewidth=2, color="#cccccc", dashes=(5,3))
ax[:plot](lat, mht, linewidth=2, color="black")
ax[:set_xlabel]("Latitude [deg]")
ax[:set_ylabel]("Heat transport [PW]")

ax[:set_title](fn)
ax[:set_xlim](-90, 90)
ax[:set_ylim](-10, 10)

plt[:show]()

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
#!/usr/bin/env julia --project=@.

using Bukdu # pipeline Conn routes resources get plug Router
using Sockets
using FileIO

using JWebImageDemo

function get_server()
    if haskey(ENV, "ON_HEROKU")
        (port=parse(Int, ENV["PORT"]), host=Sockets.IPAddr(0,0,0,0))
    else
        (port=8080, host=ip"127.0.0.1")
    end
end

server = get_server()

struct WelcomeController <: ApplicationController
    conn::Conn
end

import HTTP, HTTP.Parsers

"""
https://www.ietf.org/rfc/rfc2616.txt
"""
function parse_multipart(data::Vector{UInt8})
    chunk_end = HTTP.Parsers.find_end_of_chunk_size(data)
    if chunk_end == 0
        return nothing
    end

    boundary = data[1:chunk_end - 2] # without \r\n
    header_end = HTTP.Parsers.find_end_of_header(data)

    data_start = header_end + 1
    data_end = length(data) - chunk_end - 2
    # check that ending is with same boundary + two additional hyphens
    if (data[data_end + 1: end] == [boundary..., 0x2d, 0x2d, 0x0d, 0x0a])
       return data[data_start : data_end]
    end
    return nothing
end

function process_image(c::WelcomeController)
    data = parse_multipart(c.conn.request.body)
    dump(c.conn.request)
    if data != nothing
        buf = IOBuffer(data)
        img = FileIO.load(Stream(format"PNG", buf))
        out_img = JWebImageDemo.permute_image(img, 5)
        out_buf = IOBuffer()
        FileIO.save(Stream(format"PNG", out_buf), out_img)
        return Render("image/png", take!(out_buf))
    end
    render(Text, "Please upload some image file")
end

routes() do
    plug(Plug.Static, at="/", from=normpath(@__DIR__, "../public"))

    post("/api/process_image", WelcomeController, process_image)
end

Bukdu.start(server.port; host=server.host)

Router.call(get, "/") #
# CLI.routes()

Base.JLOptions().isinteractive==0 && wait()

# Bukdu.stop()


text:
using ITensors

#
# DMRG calculation of the extended Hubbard model
# ground state wavefunction, and spin densities
#

# Usage:
#
# Run:
#
# $ julia exthubbard.jl input.jl
#

# Include the specified input file
# Include the file from the first positional
# argument. If none is specified, use
# input.jl by default.
include(get(ARGS, 1, "input.jl"))

sweeps = Sweeps(nsweep, sweeps_args)

@show sweeps
@show N, Npart

sites = siteinds("Electron", N; conserve_qns=true)

ampo = AutoMPO()
for i=1:N
  ampo .+= U,"Nupdn",i
end
for b=1:N-1
  ampo .+= -t1,"Cdagup",b,"Cup",b+1
  ampo .+= -t1,"Cdagup",b+1,"Cup",b
  ampo .+= -t1,"Cdagdn",b,"Cdn",b+1
  ampo .+= -t1,"Cdagdn",b+1,"Cdn",b
  ampo .+= V1,"Ntot",b,"Ntot",b+1
end
for b=1:N-2
  ampo .+= -t2,"Cdagup",b,"Cup",b+2
  ampo .+= -t2,"Cdagup",b+2,"Cup",b
  ampo .+= -t2,"Cdagdn",b,"Cdn",b+2
  ampo .+= -t2,"Cdagdn",b+2,"Cdn",b
end
H = MPO(ampo,sites)

state = ["Emp" for n=1:N]
p = Ref(Npart)
for i=N:-1:1
  if p[] > i
    println("Doubly occupying site $i")
    state[i] = "UpDn"
    p[] -= 2
  elseif p[] > 0
    println("Singly occupying site $i")
    state[i] = (isodd(i) ? "Up" : "Dn")
    p[] -= 1
  end
end
# Initialize wavefunction to be bond
# dimension 10 random MPS with number
# of particles the same as `state`
psi0 = randomMPS(sites,state,10)

# Check total number of particles:
@show flux(psi0)

# Start DMRG calculation:
energy,psi = dmrg(H,psi0,sweeps)

upd = fill(0.0,N)
dnd = fill(0.0,N)
for j=1:N
  orthogonalize!(psi,j)
  psidag_j = dag(prime(psi[j], "Site"))
  upd[j] = scalar(psidag_j * op(sites, "Nup", j) * psi[j])
  dnd[j] = scalar(psidag_j * op(sites, "Ndn", j) * psi[j])
end

println("Up Density:")
for j=1:N
  println("$j $(upd[j])")
end
println()

println("Dn Density:")
for j=1:N
  println("$j $(dnd[j])")
end
println()

println("Total Density:")
for j=1:N
  println("$j $(upd[j]+dnd[j])")
end
println()

println("\nGround State Energy = $energy")


text:
const OPTIM_TYPES = Union{Int64, Rational{Int64}}

"""
    AbstractOptimizer{T} where {T <: Union{Int64, Rational{Int64}}}

An `AbstractOptimizer` is used to dispatch to various backend optimization routines like
static, synchronous, asynchronous, and others.

The type parameter encodes some of the high level behavior:

* `Int64`: The limit supplied by the type is an absolute memory limit.

* `Rational{Int64}`: The limit supplied by the type is a ratio of Large Memory (i.e. PMEM)
    to Limited Memory (i.e. DRAM)

    This will dispatch to another routine that will perform a grid search on the ratio to
    find an input ratio to supply that will result in an actual output ratio closest to
    the original input.
"""
abstract type AbstractOptimizer{T <: OPTIM_TYPES} end

const BANDWIDTHS = (
    # Split out CPU based on cases.
    cpu_pmem_dram_sync = 29000,
    cpu_dram_pmem_sync = 12000,
    cpu_pmem_dram_async = 2000,
    cpu_dram_pmem_async = 2500,

    # GPU Basically same in all diractions
    gpu = 12000,
)

_bw_remote_local_sync(::nGraph.Backend{nGraph.CPU}) = BANDWIDTHS[:cpu_pmem_dram_sync]
_bw_local_remote_sync(::nGraph.Backend{nGraph.CPU}) = BANDWIDTHS[:cpu_dram_pmem_sync]
_bw_remote_local_async(::nGraph.Backend{nGraph.CPU}) = BANDWIDTHS[:cpu_pmem_dram_async]
_bw_local_remote_async(::nGraph.Backend{nGraph.CPU}) = BANDWIDTHS[:cpu_dram_pmem_async]

_bw_remote_local_sync(::nGraph.Backend{nGraph.GPU}) = BANDWIDTHS[:gpu]
_bw_local_remote_sync(::nGraph.Backend{nGraph.GPU}) = BANDWIDTHS[:gpu]
_bw_remote_local_async(::nGraph.Backend{nGraph.GPU}) = BANDWIDTHS[:gpu]
_bw_local_remote_async(::nGraph.Backend{nGraph.GPU}) = BANDWIDTHS[:gpu]

# Bound on the actual between requested ratio for a workload.
#
# We iterate on the solution untill the ratio of PMEM to DRAM is within this bound of
# the requested ratio.
const RATIO_TOLERANCE = Ref(0.05)

function checkmargin(actual, wanted, tol = RATIO_TOLERANCE[])
    # Handle the cases where the denominator of "wanted" is zero
    rwanted, ractual = getratio.((wanted, actual))

    if iszero(rwanted.den) || iszero(rwanted.num)
        return true
    else
        return abs(ractual / rwanted - 1) <= tol
    end
end

geterr(actual, wanted) = abs(getratio(actual) / getratio(wanted) - 1)

_getratio(x::AbstractOptimizer) = x.ratio
Utils.getratio(x::AbstractOptimizer{Rational{Int64}}) = _getratio(x)
Utils.getratio(x::AbstractOptimizer{Int}) = getlimit(x)
Utils.getratio(x::Number) = x

getlimit(x::AbstractOptimizer{Int64}) = _getratio(x)

_optimizer(::T, r) where {T <: AbstractOptimizer} = T(r)

_numerator(x::AbstractOptimizer{Rational{Int64}}) = getratio(x).num
_numerator(x::AbstractOptimizer{Int64}) = getlimit(x)

_denominator(x::AbstractOptimizer{Rational{Int64}}) = getratio(x).den
_denominator(x::AbstractOptimizer{Int64}) = one(Int64)

name(::T) where {T <: AbstractOptimizer} = error("Name of $T not implemented")

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
using Test, LinearMaps, Aqua

@testset "code quality" begin
    Aqua.test_all(LinearMaps)
end

include("linearmaps.jl")

include("transpose.jl")

include("functionmap.jl")

include("scaledmap.jl")

include("composition.jl")

include("linearcombination.jl")

include("wrappedmap.jl")

include("uniformscalingmap.jl")

include("numbertypes.jl")

include("blockmap.jl")

include("kronecker.jl")

include("conversion.jl")

include("left.jl")

include("fillmap.jl")

include("nontradaxes.jl")

include("embeddedmap.jl")