text:
# result type inference

#################################################
#
#   Auxiliary functions
#
#################################################

## fptype

# @compat fptype{T<:Union{Bool,Int8,Int16,UInt8,UInt16}}(::Type{T}) = Float32
fptype{T<:Integer}(::Type{T}) = Float64
fptype{T<:AbstractFloat}(::Type{T}) = T

fptype{T<:Integer}(::Type{Complex{T}}) = Complex{fptype(T)}
fptype{T<:AbstractFloat}(::Type{Complex{T}}) = Complex{T}

## signed & unsigned type

signedtype{T<:Signed}(::Type{T}) = T
signedtype(::Type{UInt8}) = Int8
signedtype(::Type{UInt16}) = Int16
signedtype(::Type{UInt32}) = Int32
signedtype(::Type{UInt64}) = Int64
signedtype(::Type{UInt128}) = Int128

unsignedtype{T<:Unsigned}(::Type{T}) = T
unsignedtype(::Type{Int8}) = UInt8
unsignedtype(::Type{Int16}) = UInt16
unsignedtype(::Type{Int32}) = UInt32
unsignedtype(::Type{Int64}) = UInt64
unsignedtype(::Type{Int128}) = UInt128

## arithtype (unary)

if VERSION >= v"0.4-dev"
    arithtype{T<:Union{Signed,Unsigned}}(::Type{T}) = T
else
    arithtype{T<:Signed}(::Type{T}) = Int64
    arithtype{T<:Unsigned}(::Type{T}) = UInt64
end
arithtype(::Type{Bool}) = Int64
arithtype(::Type{Int128}) = Int128
arithtype(::Type{UInt128}) = UInt128

arithtype(::Type{Float32}) = Float32
arithtype(::Type{Float64}) = Float64
arithtype(::Type{Complex64}) = Complex64
arithtype(::Type{Complex128}) = Complex128
arithtype{T<:Integer}(::Type{Complex{T}}) = Complex{arithtype(T)}

## arithtype (binary)

arithtype{T<:Real}(::Type{T},::Type{T}) = arithtype(T)
arithtype{T1<:Real,T2<:Real}(::Type{T1}, ::Type{T2}) = arithtype(promote_type(T1, T2))
arithtype{T1<:Real,T2<:Real}(::Type{T1}, ::Type{Complex{T2}}) = Complex{arithtype(T1,T2)}
arithtype{T1<:Real,T2<:Real}(::Type{Complex{T1}}, ::Type{T2}) = Complex{arithtype(T1,T2)}
arithtype{T1<:Real,T2<:Real}(::Type{Complex{T1}}, ::Type{Complex{T2}}) = Complex{arithtype(T1,T2)}

#################################################
#
#   result type inference
#
#################################################

## arithmetics

# negate

result_type{T<:Number}(::Negate, ::Type{T}) = arithtype(T)

# abs

result_type{T<:Number}(::AbsFun, ::Type{T}) = arithtype(T)
result_type(::AbsFun, ::Type{Bool}) = Bool
result_type{T<:Unsigned}(::AbsFun, ::Type{T}) = T
result_type{T<:Real}(::AbsFun, ::Type{Complex{T}}) = fptype(T)

# abs2

result_type{T<:Number}(::Abs2Fun, ::Type{T}) = arithtype(T)
result_type(::Abs2Fun, ::Type{Bool}) = Bool
result_type{T<:Real}(::Abs2Fun, ::Type{Complex{T}}) = arithtype(T)

# sqr

result_type(::SqrFun, ::Type{Bool}) = Bool
result_type{T<:Number}(::SqrFun, ::Type{T}) = arithtype(T)

# real & imag

@compat result_type{T<:Real}(::Union{RealFun,ImagFun}, ::Type{T}) = T
@compat result_type{T<:Real}(::Union{RealFun,ImagFun}, ::Type{Complex{T}}) = T

# sign & signbit

result_type{T<:Real}(::SignFun, ::Type{T}) = T
result_type{T<:Real}(::SignbitFun, ::Type{T}) = Bool

# add & subtract

@compat result_type{T1<:Number,T2<:Number}(::Union{Add,Subtract}, ::Type{T1}, ::Type{T2}) = arithtype(T1, T2)

# multiply

result_type{T1<:Real,T2<:Real}(::Multiply, ::Type{T1}, ::Type{T2}) = arithtype(T1,T2)
result_type(::Multiply, ::Type{Bool}, ::Type{Bool}) = Bool
result_type{T<:Real}(::Multiply, ::Type{Bool}, ::Type{T}) = T
result_type{T<:Real}(::Multiply, ::Type{T}, ::Type{Bool}) = T

result_type{T1<:Real,T2<:Real}(::Multiply, ::Type{Complex{T1}}, ::Type{Complex{T2}}) =
    Complex{arithtype(T1, T2)}
result_type{T1<:Real,T2<:Real}(op::Multiply, ::Type{T1}, ::Type{Complex{T2}}) =
    Complex{result_type(op, T1, T2)}
result_type{T1<:Real,T2<:Real}(op::Multiply, ::Type{Complex{T1}}, ::Type{T2}) =
    Complex{result_type(op, T1, T2)}

# divide

result_type{T1<:Real,T2<:Real}(::Divide, ::Type{T1}, ::Type{T2}) = promote_type(T1, T2)
result_type{T1<:Integer,T2<:Integer}(::Divide, ::Type{T1}, ::Type{T2}) =
    promote_type(fptype(T1), fptype(T2))

result_type{T1<:Real,T2<:Complex}(op::Divide, ::Type{T1}, ::Type{T2}) =
    result_type(Multiply(), T1, fptype(T2))

result_type{T1<:Real,T2<:Real}(op::Divide, ::Type{Complex{T1}}, ::Type{T2}) =
    Complex{result_type(op, T1, T2)}

if VERSION >= v"0.4-dev"
    result_type{T1<:Real,T2<:Real}(op::Divide, ::Type{Complex{T1}},
                                   ::Type{Complex{T2}}) =
        Complex{result_type(op, T1, T2)}
else
    result_type{T1<:Complex,T2<:Complex}(op::Divide, ::Type{T1}, ::Type{T2}) =
        result_type(Multiply(), T1, fptype(T2))
end

# rdivide

result_type{T1<:Number,T2<:Number}(::RDivide, ::Type{T1}, ::Type{T2}) =
    result_type(Divide(), T2, T1)

# rcp

result_type{T<:Number}(::RcpFun, ::Type{T}) = fptype(T)

# pow

result_type(::Pow, ::Type{Bool}, ::Type{Bool}) = Bool
result_type{T<:Real}(::Pow, ::Type{T}, ::Type{Bool}) = T

result_type{T2<:Integer}(::Pow, ::Type{Bool}, ::Type{T2}) = Bool
result_type{T1<:Real,T2<:Integer}(::Pow, ::Type{T1}, ::Type{T2}) = arithtype(T1)
result_type{T1<:Real,T2<:Real}(::Pow, ::Type{T1}, ::Type{T2}) = fptype(promote_type(T1, T2))

result_type{T1<:Real}(::Pow, ::Type{Complex{T1}}, ::Type{Bool}) = Complex{T1}
result_type{T1<:Real,T2<:Integer}(::Pow, ::Type{Complex{T1}}, ::Type{T2}) = Complex{arithtype(T1)}

result_type{T1<:Real,T2<:Real}(::Pow, ::Type{Complex{T1}}, ::Type{T2}) =
    Complex{fptype(promote_type(T1, T2))}
result_type{T1<:Real,T2<:Real}(::Pow, ::Type{T1}, ::Type{Complex{T2}}) =
    Complex{fptype(promote_type(T1, T2))}
result_type{T1<:Real,T2<:Real}(::Pow, ::Type{Complex{T1}}, ::Type{Complex{T2}}) =
    Complex{fptype(promote_type(T1, T2))}

result_type{T<:Real, Tx<:Real}(::FixPow{T}, ::Type{Tx}) = result_type(Pow(), Tx, T)
result_type{T<:Real, Tx<:Real}(::FixAbsPow{T}, ::Type{Tx}) = result_type(Pow(), result_type(Abs(), Tx), T)

# max & min

@compat result_type{T<:Real}(::Union{MaxFun,MinFun}, ::Type{T}, ::Type{T}) = T
@compat result_type{T1<:Real,T2<:Real}(::Union{MaxFun,MinFun}, ::Type{T1}, ::Type{T2}) = promote_type(T1, T2)

# quotient & module

@compat result_type{T1<:Real,T2<:Real}(::Union{DivFun,RemFun}, ::Type{T1}, ::Type{T2}) = promote_type(T1, T2)

if VERSION >= v"0.4-dev"
    @compat result_type{T1<:Signed,T2<:Unsigned}(op::Union{DivFun,RemFun}, ::Type{T1}, ::Type{T2}) =
        signedtype(result_type(op, unsignedtype(T1), T2))
    @compat result_type{T1<:Unsigned,T2<:Signed}(op::Union{DivFun,RemFun}, ::Type{T1}, ::Type{T2}) =
        result_type(op, T1, unsignedtype(T2))
else
    @compat result_type{T1<:Signed,T2<:Unsigned}(::Union{DivFun,RemFun}, ::Type{T1}, ::Type{T2}) =
        signedtype(promote_type(T1, T2))
    @compat result_type{T1<:Unsigned,T2<:Signed}(::Union{DivFun,RemFun}, ::Type{T1}, ::Type{T2}) =
        unsignedtype(promote_type(T1, T2))
end

result_type{T1<:Real,T2<:Real}(::FldFun, ::Type{T1}, ::Type{T2}) =
    arithtype(T1, T2)
@compat result_type{T1<:Union{Unsigned,Bool}, T2<:Union{Unsigned,Bool}}(::FldFun, ::Type{T1}, ::Type{T2}) =
    promote_type(T1, T2)

result_type{T1<:Signed,T2<:Unsigned}(::FldFun, ::Type{T1}, ::Type{T2}) =
    result_type(DivFun(), T1, T2)
result_type{T1<:Unsigned,T2<:Signed}(::FldFun, ::Type{T1}, ::Type{T2}) =
    result_type(DivFun(), T1, T2)

result_type{T1<:Real,T2<:Real}(::ModFun, ::Type{T1}, ::Type{T2}) = promote_type(T1, T2)
result_type{T1<:Signed,T2<:Unsigned}(::ModFun, ::Type{T1}, ::Type{T2}) =
    unsignedtype(result_type(DivFun(), T1, T2))
result_type{T1<:Unsigned,T2<:Signed}(::ModFun, ::Type{T1}, ::Type{T2}) =
    signedtype(result_type(DivFun(), T1, T2))

# fma

result_type{T1<:Number,T2<:Number,T3<:Number}(::FMA, ::Type{T1}, ::Type{T2}, ::Type{T3}) =
    result_type(Add(), T1, result_type(Multiply(), T2, T3))

# ifelse

result_type{T<:Number}(::IfelseFun, ::Type{Bool}, ::Type{T}, ::Type{T}) = T

# logical

result_type(::Not, ::Type{Bool}) = Bool
result_type(::And, ::Type{Bool}, ::Type{Bool}) = Bool
result_type(::Or, ::Type{Bool}, ::Type{Bool}) = Bool

# bitwise

result_type{T<:Integer}(::BitwiseNot, ::Type{T}) = T

@compat result_type{T1<:Integer,T2<:Integer}(::Union{BitwiseAnd,BitwiseOr,BitwiseXor}, ::Type{T1}, ::Type{T2}) =
    promote_type(T1, T2)

# comparison

@compat result_type{T1<:Real,T2<:Real}(::Union{LT,GT,LE,GE}, ::Type{T1}, ::Type{T2}) = Bool
@compat result_type{T1<:Number,T2<:Number}(::Union{EQ,NE}, ::Type{T1}, ::Type{T2}) = Bool

# rounding

@compat result_type{T<:Real}(::Union{FloorFun, CeilFun, TruncFun, RoundFun}, ::Type{T}) = T
@compat result_type{T<:Integer}(::Union{IfloorFun, IceilFun, ItruncFun, IroundFun}, ::Type{T}) = T
@compat result_type{T<:AbstractFloat}(::Union{IfloorFun, IceilFun, ItruncFun, IroundFun}, ::Type{T}) = Int64

# number classification

@compat result_type{T<:Real}(::Union{IsnanFun, IsinfFun, IsfiniteFun}, ::Type{T}) = Bool

# algebraic functions

@compat result_type{T<:Number}(::Union{SqrtFun,CbrtFun,RsqrtFun,RcbrtFun}, ::Type{T}) = fptype(T)
result_type{T1<:Real,T2<:Real}(::HypotFun, ::Type{T1}, ::Type{T2}) = promote_type(fptype(T1), fptype(T2))

# exponential & logarithm

@compat result_type{T<:Number}(::Union{ExpFun,Exp2Fun,Exp10Fun}, ::Type{T}) = fptype(T)
@compat result_type{T<:Number}(::Union{LogFun,Log2Fun,Log10Fun}, ::Type{T}) = fptype(T)
@compat result_type{T<:Real}(::Union{Expm1Fun,Log1pFun}, ::Type{T}) = fptype(T)

@compat result_type{T<:Real}(::Union{XlogxFun,LogisticFun,LogitFun,SoftplusFun,InvsoftplusFun}, ::Type{T}) = fptype(T)
result_type{T1<:Real,T2<:Real}(::XlogyFun, ::Type{T1}, ::Type{T2}) = fptype(promote_type(T1, T2))
result_type{T1<:Real,T2<:Real}(::LogsumexpFun, ::Type{T1}, ::Type{T2}) = fptype(promote_type(T1, T2))

# trigonometric functions

@compat result_type{T<:Number}(::Union{SinFun,CosFun,TanFun,CotFun,SecFun,CscFun}, ::Type{T}) = fptype(T)
@compat result_type{T<:Number}(::Union{AsinFun,AcosFun,AtanFun,AcotFun,AsecFun,AcscFun}, ::Type{T}) = fptype(T)
@compat result_type{T<:Number}(::Union{SincFun,CoscFun,SinpiFun,CospiFun}, ::Type{T}) = fptype(T)
@compat result_type{T<:Number}(::Union{SindFun,CosdFun,TandFun,CotdFun,SecdFun,CscdFun}, ::Type{T}) = fptype(T)
@compat result_type{T<:Number}(::Union{AsindFun,AcosdFun,AtandFun,AcotdFun,AsecdFun,AcscdFun}, ::Type{T}) = fptype(T)

result_type{T<:Integer}(::SincFun, ::Type{T}) = T
if VERSION < v"0.4-dev"
    result_type{T<:Integer}(::SinpiFun, ::Type{T}) = T
    result_type{T<:Integer}(::CospiFun, ::Type{T}) = arithtype(T)
end

result_type{T1<:Real,T2<:Real}(::Atan2Fun, ::Type{T1}, ::Type{T2}) = promote_type(fptype(T1), fptype(T2))

# hyperbolic functions

@compat result_type{T<:Number}(::Union{SinhFun,CoshFun,TanhFun,CothFun,SechFun,CschFun}, ::Type{T}) = fptype(T)
@compat result_type{T<:Number}(::Union{AsinhFun,AcoshFun,AtanhFun,AcothFun,AsechFun,AcschFun}, ::Type{T}) = fptype(T)

# erf & friends

@compat result_type{T<:Number}(::Union{ErfFun,ErfcFun,ErfiFun,ErfcxFun}, ::Type{T}) = fptype(T)
@compat result_type{T<:Number}(::Union{ErfinvFun,ErfcinvFun}, ::Type{T}) = fptype(T)

# gamma, beta, & friends

@compat result_type{T<:Number}(::Union{GammaFun,LgammaFun,EtaFun,ZetaFun}, ::Type{T}) = fptype(T)
@compat result_type{T<:Integer}(::Union{EtaFun,ZetaFun}, ::Type{T}) = Float64

@compat result_type(::Union{GammaFun,LgammaFun}, ::Type{Complex64}) = Complex128

result_type{T<:AbstractFloat}(::DigammaFun, ::Type{T}) = T
result_type{T<:Integer}(::DigammaFun, ::Type{T}) = Float64

@compat result_type{T1<:Real,T2<:Real}(::Union{BetaFun,LbetaFun}, ::Type{T1}, ::Type{T2}) =
    promote_type(fptype(T1),fptype(T2))

@compat result_type{T1<:Integer,T2<:Integer}(::Union{BetaFun,LbetaFun}, ::Type{T1}, ::Type{T2}) = Float64

# airy & friends

@compat result_type{T<:Number}(::Union{AiryFun,AiryprimeFun,AiryaiFun,AiryaiprimeFun,AirybiFun,AirybiprimeFun},::Type{T}) =
    fptype(T)

# bessel & friends

@compat result_type{T<:Number}(::Union{Besselj0Fun, Besselj1Fun, Bessely0Fun, Bessely1Fun}, ::Type{T}) = fptype(T)
@compat result_type{T<:Number}(::Union{BesseliFun, BesseljFun, BesselkFun, BesselyFun}, ::Type{T}) = fptype(T)
@compat result_type{T<:Number}(::Union{Hankelh1Fun, Hankelh2Fun}, ::Type{T}) = fptype(T)


text:
import ComponentArrays

"""
    namedtuple(vi::DynamicPPL.TypedVarInfo)

Similar to `DynamicPPL.tonamedtuple` but values of return-value are only the values
of the different variables rather than a tuple of values and symbols.
"""
namedtuple(vi::DynamicPPL.TypedVarInfo) = _namedtuple(vi.metadata)
@generated function _namedtuple(metadata::NamedTuple{names}, start = 0) where {names}
    expr = Expr(:tuple)
    start = :(1)
    for f in names
        length = :(length(metadata.$f.vals))
        push!(expr.args, :(metadata.$f.vals[1:$length]))
    end
    return :(NamedTuple{$names}($expr))
end

# Constructor for `ComponentArray` from a `VarInfo`
function ComponentArrays.ComponentArray(varinfo::DynamicPPL.VarInfo)
    return _tocomponentarray(varinfo.metadata)
end

function ComponentArrays.ComponentArray(varinfo::DynamicPPL.Model)
    varinfo = DynamicPPL.VarInfo(model)
    return ComponentArrays.ComponentArray(varinfo)
end

@generated function _tocomponentarray(metadata::NamedTuple{names}) where {names}
    # Adapted from https://github.com/TuringLang/DynamicPPL.jl/blob/873d5e94fe778514d52e24338a3cd1bd0535eefa/src/varinfo.jl#L338-L348
    # and thus ought should allow usage such as
    #
    #    varinfo[sampler] .= ComponentArrays.getdata(ComponentArray(varinfo))
    #
    # i.e. have the same order as `varinfo[sampler]`.
    ranges = Expr(:tuple)
    expr = Expr(:vcat)
    start = :(1)
    for f in names
        length = :(length(metadata.$f.vals))
        finish = :($start + $length - 1)
        r = :($start:$finish)
        push!(ranges.args, r)
        push!(expr.args, :(metadata.$f.vals[1:$length]))
        start = :($start + $length)
    end

    return quote
        ax = $(ComponentArrays).Axis(NamedTuple{$names}($ranges))
        return $(ComponentArrays).ComponentArray($expr, ax)
    end
end

# Make it so that `NamedBijector` is compatible with `ComponentArray`
vec_mby(x) = x
vec_mby(x::Real) = [x]
vec_mby(x::AbstractArray) = vec(x)

@generated function (b::Bijectors.NamedBijector{names1})(
    x::ComponentArrays.ComponentVector
) where {names1}
    ax = first(ComponentArrays.getaxes(x))
    indexmap = ComponentArrays.indexmap(ax)

    expr = Expr(:call, :vcat)
    for (k, v) in pairs(indexmap)
        if k in names1
            push!(expr.args, :(vec_mby(b.bs.$k(x.$k))))
        else
            push!(expr.args, :(vec_mby(x.$k)))
        end
    end

    # Re-uses the implementation for `Stacked`, thus acting directly on the vector
    # rather than accessing the fields via the named fields.
    return :(ComponentArrays.ComponentArray($expr, $ax))
end

# Just works, but impl in Bijectors.jl unnecessarily specializes on `x::NamedTuple`
@generated function Bijectors.logabsdetjac(b::Bijectors.NamedBijector{names}, x) where {names}
    exprs = [:($(Bijectors).logabsdetjac(b.bs.$n, x.$n)) for n in names]
    return :(+($(exprs...)))
end

ext: jl


text:
module ESS

function all_help_topics()
    Base.Help.init_help()
    ## show all categories
    for cat = Base.Help.CATEGORY_LIST
        if !isempty(Base.Help.CATEGORY_DICT[cat])
            println()
            show(cat); println();
            for func = Base.Help.CATEGORY_DICT[cat]
                print("  ")
                show(func)
            end
        end
    end
end

## modified version of function show(io::IO, m::Method)
function fun_args(m::Method)
    tv = m.tvars
    io = OUTPUT_STREAM::IO
    if !isa(tv,Tuple)
        tv = (tv,)
    end
    if !isempty(tv)
        Base.show_delim_array(io, tv, '{', ',', '}', false)
    end
    li = m.func.code
    e = Base.uncompressed_ast(li)
    argnames = e.args[1]
    decls = map(Base.argtype_decl_string, argnames, {m.sig...})
    print(io, "(")
    print_joined(io, decls, ",", ",")
    print(io, ")")
end

## modified versionof show(io::IO, mt::MethodTable)
function fun_args(f::Function)
    mt = f.env
    print("(list nil nil '(")
    d = mt.defs
    while !is(d,())
        print("\"")
        fun_args(d)
        print("\" ")
        d = d.next
        if is(d,())
            print("))")
        end
    end
end

function fun_args(s::String)
    fun_args(eval(parse(s)))
end

end

max_stars_repo_name: jottos/emacs

max_issues_repo_name: jottos/emacs

max_forks_repo_name: jottos/emacs


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    #precompile(Tuple{typeof(Base.Printf.ini_dec),Base.GenericIOBuffer{Array{UInt8,1}},Float64,Int64,String,Int64,Int64,Char,Array{UInt8,1}})
    precompile(Tuple{typeof(Base.__cat),Array{Any,2},Tuple{Int64,Int64},Tuple{Bool,Bool},String,Vararg{Any,N} where N})
    precompile(Tuple{typeof(Base.__cat),Array{Float64,1},Tuple{Int64},Tuple{Bool},Array{Float64,1},Vararg{Any,N} where N})
    precompile(Tuple{typeof(Base.diff_names),NTuple{5,Symbol},Tuple{Symbol}})
    precompile(Tuple{typeof(Base.merge_names),NTuple{5,Symbol},Tuple{Symbol}})
    precompile(Tuple{typeof(Base.merge_types),NTuple{4,Symbol},Type{NamedTuple{(:clip, :first),Tuple{Nothing,Bool}}},Type{NamedTuple{(:W, :Vd),Tuple{Float64,Int64}}}})
    precompile(Tuple{typeof(Base.merge_types),NTuple{4,Symbol},Type{NamedTuple{(:first,),Tuple{Bool}}},Type{NamedTuple{(:F, :D, :par),Tuple{String,String,Tuple{Symbol,Int64}}}}})
    precompile(Tuple{typeof(Base.merge_types),NTuple{4,Symbol},Type{NamedTuple{(:first,),Tuple{Bool}}},Type{NamedTuple{(:pos, :B, :Vd),Tuple{NamedTuple{(:anchor,),Tuple{String}},String,Int64}}}})
    precompile(Tuple{typeof(Base.merge_types),Tuple{Symbol,Symbol},Type{NamedTuple{(:clip, :first),Tuple{Nothing,Bool}}},Type{NamedTuple{(),Tuple{}}}})
    precompile(Tuple{typeof(Base.merge_types),Tuple{Symbol,Symbol},Type{NamedTuple{(:first, :show),Tuple{Bool,Bool}}},Type{NamedTuple{(:show,),Tuple{Bool}}}})
    precompile(Tuple{typeof(Base.merge_types),Tuple{Symbol,Symbol},Type{NamedTuple{(:first,),Tuple{Bool}}},Type{NamedTuple{(:show,),Tuple{Bool}}}})
    precompile(Tuple{typeof(Base.merge_types),Tuple{Symbol},Type{NamedTuple{(:first,),Tuple{Bool}}},Type{NamedTuple{(),Tuple{}}}})
    precompile(Tuple{typeof(copyto!),Array{Float64,1},Tuple{Float64,Int64,Int64}})
    precompile(Tuple{typeof(haskey),Dict{Symbol,Any},String})
    precompile(Tuple{typeof(map),Type,NTuple{160,Char}})
    precompile(Tuple{typeof(map),Type,NTuple{320,Char}})
    precompile(Tuple{typeof(map),Type,NTuple{80,Char}})
    precompile(Tuple{typeof(merge),NamedTuple{(:first, :show),Tuple{Bool,Bool}},Base.Iterators.Pairs{Symbol,Bool,Tuple{Symbol},NamedTuple{(:show,),Tuple{Bool}}}})

    precompile(Tuple{Core.kwftype(typeof(GMT.grdimage)),NamedTuple{(:first, :show),Tuple{Bool,Bool}},typeof(grdimage),String,GMT.GMTgrid})
    precompile(Tuple{Core.kwftype(typeof(GMT.imshow)),NamedTuple{(:show,),Tuple{Bool}},typeof(imshow),String})
    precompile(Tuple{typeof(GMT.GMTJL_Set_Object),Ptr{Nothing},GMT.GMT_RESOURCE,GMT.GMTcpt})
    precompile(Tuple{typeof(GMT.GMTJL_Set_Object),Ptr{Nothing},GMT.GMT_RESOURCE,GMT.GMTgrid})
    precompile(Tuple{typeof(GMT.add_opt),Dict{Symbol,Any},String,String,Array{Any,2}})
    precompile(Tuple{typeof(GMT.add_opt_cpt),Dict{Symbol,Any},String,Array{Symbol,2},Char,Int64,Array{Float64,2}})
    precompile(Tuple{typeof(GMT.common_shade),Dict{Symbol,Any},String,GMT.GMTgrid,GMT.GMTcpt,Nothing,Nothing,String})
    precompile(Tuple{typeof(GMT.finish_PS_module),Dict{Symbol,Any},String,String,Bool,Bool,Bool,Array{Float64,2},Vararg{Any,N} where N})
    precompile(Tuple{typeof(GMT.finish_PS_module),Dict{Symbol,Any},String,String,Bool,Bool,Bool,GMT.GMTgrid,Vararg{Any,N} where N})
    precompile(Tuple{typeof(GMT.get_cpt_set_R),Dict{Symbol,Any},String,String,String,Int64,GMT.GMTgrid,Nothing,Nothing,String})
    precompile(Tuple{typeof(GMT.get_marker_name),Dict{Symbol,Any},Array{Symbol,2},Bool,Bool,Array{Float64,2}})
    precompile(Tuple{typeof(GMT.make_color_column),Dict{Symbol,Any},String,String,Int64,Int64,Int64,Bool,Bool,Array{Float64,2},Nothing})
    precompile(Tuple{typeof(GMT.put_in_legend_bag),Dict{Symbol,Any},String,Array{Float64,2}})
    precompile(Tuple{typeof(GMT.round_wesn),Array{Float64,2}})
    precompile(Tuple{typeof(convert),Type{Ptr{GMT.GMT_MATRIX}},Ptr{Nothing}})
    precompile(Tuple{typeof(plot),Array{Float64,1},Array{Float64,1}})
    precompile(Tuple{typeof(setproperty!),GMT.GMT_GRID_HEADER,Symbol,NTuple{320,UInt8}})
    precompile(Tuple{typeof(setproperty!),GMT.GMT_GRID_HEADER,Symbol,NTuple{80,UInt8}})
end


text:
@inline function cv_max(x::T,xL::T,xU::T,a::S) where {S,T<:AbstractFloat}
        ca::T = convert(T,a)
        if (xU<=ca)
          return ca, zero(T)
        elseif (ca<=xL)
          return x, one(T)
        else
          val::T = ca + (xU-ca)*(max(0.0,((x-ca)/(xU-ca))))^(MC_param.mu+1)
          dval::T = ((xU-ca)/(xU-ca))*(MC_param.mu+1)*(max(0.0,((x-ca)/(xU-ca))))^(MC_param.mu)
          return val, dval
        end
end
@inline function cc_max(x::T,xL::T,xU::T,a::S) where {S,T<:AbstractFloat}
  ca::T = convert(T,a)
  if (xU<=ca)
    d = zero(T)
  elseif (ca<=xL)
    d = one(T)
  else
    d = ((xU-ca)/(xU-ca))*(MC_param.mu+1)*(max(0.0,((x-ca)/(xU-ca))))^(MC_param.mu)
  end
  return line_seg(x,xL,max(xL,a),xU,max(xU,a)),dline_seg(x,xL,max(xL,a),xU,max(xU,a),d)
end
function cc_max_NS(x::T,lo::T,hi::T,c::S) where {S,T<:AbstractFloat}
  return line_seg(x,lo,max(lo,c),hi,max(hi,c)),dline_seg(x,lo,max(lo,c),hi,max(hi,c),((x>c) ? 1.0 : 0.0))
end
function cv_max_NS(x::T,xL::T,xU::T,c::S) where {S,T<:AbstractFloat}
  return max(x,c),((x>c) ? 1.0 : 0.0)
end

for i in union(int_list, float_list)
	eval( quote
	char = $i
@inline function max(x::SMCg{N,T},c::$i) where {N,T<:AbstractFloat}
  eps_min::T = x.Intv.lo
  eps_max::T = x.Intv.hi
	midcc::T,cc_id::Int64 = mid3(x.cc,x.cv,eps_max)
	midcv::T,cv_id::Int64 = mid3(x.cc,x.cv,eps_min)
  if (MC_param.mu >= 1)
     #println("ran me max 1")
	   cc::T,dcc::T = cc_max(midcc,x.Intv.lo,x.Intv.hi,c)
     cv::T,dcv::T = cv_max(midcv,x.Intv.lo,x.Intv.hi,c)
     gcc1::T,gdcc1::T = cc_max(x.cv,x.Intv.lo,x.Intv.hi,c)
 	   gcv1::T,gdcv1::T = cv_max(x.cv,x.Intv.lo,x.Intv.hi,c)
 	   gcc2::T,gdcc2::T = cc_max(x.cc,x.Intv.lo,x.Intv.hi,c)
 	   gcv2::T,gdcv2::T = cv_max(x.cc,x.Intv.lo,x.Intv.hi,c)
 	   cv_grad::SVector{N,T} = max(zero(T),gdcv1)*x.cv_grad + min(zero(T),gdcv2)*x.cc_grad
 	   cc_grad::SVector{N,T} = min(zero(T),gdcc1)*x.cv_grad + max(zero(T),gdcc2)*x.cc_grad
  else
    #println("ran me max 2")
    cc,dcc = cc_max_NS(midcc,x.Intv.lo,x.Intv.hi,c)
    cv,dcv = cv_max_NS(midcv,x.Intv.lo,x.Intv.hi,c)
    cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
  end

  return SMCg{N,T}(cc, cv, cc_grad, cv_grad, max(x.Intv),x.cnst,x.IntvBox,x.xref)
end
@inline max(c::$i,x::SMCg{N,T}) where {N,T<:AbstractFloat} = max(x,c)
@inline min(c::$i,x::SMCg{N,T}) where {N,T<:AbstractFloat} = -max(-x,-c)
@inline min(x::SMCg{N,T},c::$i) where {N,T<:AbstractFloat} = -max(-x,-c)
         end )
 end
# defines functions on which bivariant maximum mapping from Khan 2016
@inline function psil_max(x::T,y::T,lambda::Interval{T},nu::Interval{T},
                          f1::SMCg{N,T},f2::SMCg{N,T}) where {N,T<:AbstractFloat}
   if (nu.hi<=lambda.lo)
     val::T = x
   elseif (lambda.hi<=nu.lo)
     val = y
   elseif ((nu.lo<=lambda.lo)&&(lambda.lo<nu.hi))
     val = x+(nu.hi-lambda.lo)*max(zero(T),((y-x)/(nu.hi-lambda.lo)))^(MC_param.mu+1)
   else
     val =  y + (lambda.hi-nu.lo)*max(zero(T),(x-y)/(lambda.hi-nu.lo))^(MC_param.mu+1)
   end
   if (nu.hi <= lambda.lo)
     grad_val::SVector{N,T} = f1.cv_grad
   elseif (lambda.hi <= nu.lo)
     grad_val = f1.cc_grad
   else
     grad_val = max(zero(T),psil_max_dx(x,y,lambda,nu))*f1.cv_grad +
                min(zero(T),psil_max_dx(x,y,lambda,nu))*f1.cc_grad +
                max(zero(T),psil_max_dy(x,y,lambda,nu))*f2.cv_grad +
                min(zero(T),psil_max_dy(x,y,lambda,nu))*f2.cc_grad
   end
   return val,grad_val
end
@inline function thetar(x::T,y::T,lambda::Interval{T},nu::Interval{T}) where {T<:AbstractFloat}
    return (max(lambda.lo,nu.lo) + max(lambda.hi,nu.hi)-max(lambda.lo,nu.hi) +
    max(lambda.hi,nu.lo))*max(zero(T),((lambda.hi-x)/(lambda.hi-lambda.lo)-(y-nu.lo)/(nu.hi-nu.lo)))^(MC_param.mu+1)
end
function psil_max_dx(x::T,y::T,lambda::Interval{T},nu::Interval{T}) where {T<:AbstractFloat}
  if (nu.lo <= lambda.lo < nu.hi)
    return one(T)-(MC_param.mu+1)*max(zero(T),(y-x)/(nu.hi-lambda.lo))^MC_param.mu
  else
    return (MC_param.mu+1)*max(zero(T),(x-y)/(lambda.hi-nu.lo))^MC_param.mu
  end
end
function psil_max_dy(x::T,y::T,lambda::Interval{T},nu::Interval{T}) where {T<:AbstractFloat}
  if (nu.lo <= lambda.lo < nu.hi)
    return (MC_param.mu+1)*max(zero(T),(y-x)/(nu.hi-lambda.lo))^MC_param.mu
  else
    return one(T)-(MC_param.mu+1)*max(zero(T),(x-y)/(lambda.hi-nu.lo))^MC_param.mu
  end
end

@inline function psir_max(x::T,y::T,xgrad::SVector{N,T},ygrad::SVector{N,T},
                          lambda::Interval{T},nu::Interval{T}) where {N,T<:AbstractFloat}
    if (nu.hi<=lambda.lo)
      return x,xgrad
    elseif (lambda.hi<=nu.lo)
      return y,ygrad
    else
      val::T = max(lambda.hi,nu.hi)-(max(lambda.hi,nu.hi)-max(lambda.lo,nu.hi))*
          ((lambda.hi-x)/(lambda.hi-lambda.lo))-
          (max(lambda.hi,nu.hi)-max(lambda.hi,nu.lo))*((nu.hi-y)/(nu.hi-nu.lo))
          +thetar(x,y,nu,lambda)
      coeff = [(max(lambda.hi,nu.hi)-max(lambda.lo,nu.hi))/(lambda.hi-lambda.lo)
               (max(lambda.hi,nu.hi)-max(nu.lo,lambda.hi))/(nu.hi-nu.lo)
               (MC_param.mu+1)*(max(lambda.hi,nu.lo)+max(lambda.lo,nu.hi)-max(lambda.lo,nu.lo)-max(lambda.hi,nu.hi))
               max(0,((lambda.hi-x)/(lambda.hi-lambda.lo))-((y-nu.lo)/(nu.hi-nu.lo)))^MC_param.mu
               1/(lambda.hi-lambda.lo)
               1/(nu.hi-nu.lo)
               ]
      grad_val::SVector{N,T} = coeff[1]*xgrad + coeff[2]*ygrad +
                 coeff[3]*coeff[4]*(coeff[5]*xgrad+coeff[6]*ygrad)
      return val,grad_val
    end
end

@inline function max(x::SMCg{N,T},y::SMCg{N,T}) where {N,T<:AbstractFloat}
    if (MC_param.mu >= 1)
      cc::T = zero(T)
      cv::T = zero(T)
      temp_mid::T = zero(T)
      if ((y.Intv.hi<=x.Intv.lo)||(x.Intv.hi<=y.Intv.lo))
        cv,cv_grad::SVector{N,T} = psil_max(x.cv,y.cv,x.Intv,y.Intv,x,y)
      elseif ((y.Intv.lo<=x.Intv.lo) & (x.Intv.lo<y.Intv.hi))
        temp_mid,blank = mid3(x.cv,x.cc,y.cv-(y.Intv.hi-x.Intv.lo)*(MC_param.mu+1)^(-1/MC_param.mu))
        cv,cv_grad = psil_max(temp_mid,y.cv,x.Intv,y.Intv,x,y)
      elseif ((x.Intv.lo<y.Intv.lo) & (y.Intv.lo<x.Intv.hi))
        temp_mid,blank = mid3(y.cv,y.cc,x.cv-(x.Intv.hi-y.Intv.lo)*(MC_param.mu+1)^(-1/MC_param.mu))
        cv,cv_grad = psil_max(x.cv,temp_mid,x.Intv,y.Intv,x,y)
      end
      cc,cc_grad::SVector{N,T} = psir_max(x.cc,y.cc,x.cc_grad,y.cv_grad,x.Intv,y.Intv)
      return SMCg{N,T}(cc, cv, cc_grad, cv_grad, max(x.Intv,y.Intv),(x.cnst && y.cnst),x.IntvBox,x.xref)
    elseif (x.Intv.hi <= y.Intv.lo)
      cc = y.cc
      cc_grad = y.cnst ? zeros(y.cc_grad) : y.cc_grad
    elseif (x.Intv.lo >= y.Intv.hi)
      cc = x.cc
      cc_grad = x.cnst ? zeros(x.cc_grad) : x.cc_grad
    elseif (MC_param.multivar_refine)
      maxLL::T = max(x.Intv.lo,y.Intv.lo)
      maxLU::T = max(x.Intv.lo,y.Intv.hi)
      maxUL::T = max(x.Intv.hi,y.Intv.lo)
      maxUU::T = max(x.Intv.hi,y.Intv.hi)
      thin1::Bool = (diam(x.Intv) == zero(T))
      thin2::Bool = (diam(y.Intv) == zero(T))
      r11::T = thin1 ? zero(T) : (maxUL-maxLL)/diam(x.Intv)
      r21::T = thin1 ? zero(T) : (maxLU-maxUU)/diam(x.Intv)
      r12::T = thin2 ? zero(T) : (maxLU-maxLL)/diam(y.Intv)
      r22::T = thin2 ? zero(T) : (maxUL-maxUU)/diam(y.Intv)
      cc1::T = maxLL + r11*(x.cc-x.Intv.lo) + r12*(y.cc-y.Intv.lo)
      cc2::T = maxUU - r21*(x.cc-x.Intv.hi) - r22*(y.cc-y.Intv.hi)
      if (cc1 <= cc2)
        cc = cc1
        cc_grad =  (x.cnst ? zeros(y.cc_grad) : r11*x.cc_grad) + (y.cnst ? zeros(x.cc_grad) : r12*y.cc_grad)
      else
        cc = cc2
        cc_grad = -(x.cnst ? zeros(y.cc_grad) : r21*x.cc_grad) - (y.cnst ? zeros(x.cc_grad) : r22*y.cc_grad)
      end
    else
      ccMC::SMCg{N,T} = (x+y+abs(x-y))/2
      cc = ccMC.cc
      cc_grad = ccMC.cc_grad
    end
    cv = max(x.cv,y.cv)
    cv_grad = (x.cv > y.cv) ? (x.cnst ? zeros(x.cv_grad): x.cv_grad) :
                              (y.cnst ? zeros(y.cv_grad): y.cv_grad)
    cnst = y.cnst ? x.cnst : (x.cnst ? y.cnst : (x.cnst || y.cnst) )

    return SMCg{N,T}(cc, cv, cc_grad, cv_grad, max(x.Intv,y.Intv),cnst,x.IntvBox,x.xref)
end

@inline function maxcv(x::SMCg{N,T},y::SMCg{N,T}) where {N,T<:AbstractFloat}
        cv::T = zero(T)
        temp_mid::T = zero(T)
        if ((y.Intv.hi<=x.Intv.lo)||(x.Intv.hi<=y.Intv.lo))
            cv = psil_max(x.cv,y.cv,x.Intv,y.Intv)
        elseif ((y.Intv.lo<=x.Intv.lo) & (x.Intv.lo<y.Intv.hi))
          temp_mid = mid3(x.cv,x.cc,y.cv-(y.Intv.hi-x.Intv.lo)*(MC_param.mu+1)^(-1.0/MC_param.mu))
          cv = psil_max(temp_mid,y.cv,x.Intv,y.Intv)
        elseif ((x.Intv.lo<y.Intv.lo) & (y.Intv.lo<x.Intv.hi))
          temp_mid = mid3(y.cv,y.cc,x.cv-(x.Intv.hi-y.Intv.lo)*(MC_param.mu+1)^(-1.0/MC_param.mu))
          cv = psil_max(x.cv,temp_mid,x.Intv,y.Intv)
        end
        return cv
end
@inline function maxcc(x,y::SMCg{N,T}) where {N,T<:AbstractFloat}
        cc::T = psir_max(x.cc,y.cc,x.Intv,y.Intv)
end
@inline mincv(x,y::SMCg{N,T}) where {N,T<:AbstractFloat} = - maxcc(-x,-y)
@inline min(x::SMCg{N,T},y::SMCg{N,T}) where {N,T<:AbstractFloat} = -max(-x,-y)
@inline max(x::SMCg{N,T},y::Interval{T}) where {N,T<:AbstractFloat} = max(x,SMCg{N,T}(y))
@inline max(y::Interval{T},x::SMCg{N,T}) where {N,T<:AbstractFloat} = max(x,SMCg{N,T}(y))
@inline min(x::SMCg{N,T},y::Interval{T}) where {N,T<:AbstractFloat} = min(x,SMCg{N,T}(y))
@inline min(y::Interval{T},x::SMCg{N,T}) where {N,T<:AbstractFloat} = min(x,SMCg{N,T}(y))


text:
"""
    expand_poly(x, [degree = 5], [obsdim]) -> Matrix

Perform a polynomial basis expansion of the given `degree` for
the vector `x`.

The optional parameter `obsdim` can be used to specify which
dimension of the resulting `Matrix` should denote the
observations. By default the each column will denote an
observation, which means that the returning `Matrix` will have
the size `(degree, length(x))`.

```julia
julia> expand_poly(1:5, degree = 3)
3×5 Array{Float64,2}:
 1.0  2.0   3.0   4.0    5.0
 1.0  4.0   9.0  16.0   25.0
 1.0  8.0  27.0  64.0  125.0

julia> expand_poly(1:5, 3); # same but type-stable
```

If you want each row to denote an observation you can set
`obsdim = 1`.

```julia
julia> expand_poly(1:5, degree = 3, obsdim = 1)
5×3 Array{Float64,2}:
 1.0   1.0    1.0
 2.0   4.0    8.0
 3.0   9.0   27.0
 4.0  16.0   64.0
 5.0  25.0  125.0

julia> expand_poly(1:5, 3, ObsDim.First()); # same but type-stable
```
"""
function expand_poly(x::AbstractVector; degree::Integer = 5, obsdim = ObsDim.Last())
    expand_poly(x, degree, convert(LearnBase.ObsDimension, obsdim))
end

function expand_poly(x::AbstractVector, degree::Integer, ::ObsDim.Last)
    expand_poly(x, degree, ObsDim.Constant{2}())
end

function expand_poly{T<:Number}(x::AbstractVector{T}, degree::Integer, ::ObsDim.Constant{2} = ObsDim.Constant{2}())
    n = length(x)
    # FIXME: consider support for colorants (i.e. don't hardcode Float64)
    x_vec = convert(Vector{Float64}, x)
    X = zeros(Float64, (degree, n))
    @inbounds for i = 1:n
        for d = 1:degree
            X[d, i] += x_vec[i]^(d)
        end
    end
    X
end

function expand_poly{T<:Number}(x::AbstractVector{T}, degree::Integer, ::ObsDim.Constant{1})
    n = length(x)
    # FIXME: consider support for colorants (i.e. don't hardcode Float64)
    x_vec = convert(Vector{Float64}, x)
    X = zeros(Float64, (n, degree))
    @inbounds for d = 1:degree
        for i = 1:n
            X[i, d] += x_vec[i]^(d)
        end
    end
    X
end


text:
@testset "Test for SISO data" begin
    # Generate Input Data
    len = 8
    srand(1)
    u = rand(1, len)-0.5;

    # Define Dynamic Model
    mdl = Linear(5, 1)
    yterms = [[1, 2]]
    uterms = [[0, 1, 2]]
    diffeq = DifferenceEquation(mdl, yterms=yterms, uterms=uterms)

    # Simulate system for a given parameter
    Θ = [-0.1, 0.3, 1, 0.4, 0.3]
    y0 = reshape([0.0, 0.0], 1, 2)
    ϕ = [Θ; vec(y0)]
    yaux = simulate(diffeq, u, Θ, y0)
    y = [y0 yaux]

    @testset "Test one-step-ahead prediction" begin
        pred = one_step_ahead(mdl, yterms, uterms, IdData(y, u))

        @testset "Test Basic Construction" begin
            @test pred.dynmdl.yterms == Vector{Vector{Int}}(0)
            @test pred.dynmdl.uterms == [[1, 2], [0, 1, 2]]
            @test pred.u == [y; u]
            @test pred.y0 == Matrix{Float64}(0, 0)
            @test pred.length == size(yaux, 2)
            @test pred.start == 3
        end

        @testset "Test Prediction" begin
            @test predict(pred, Θ) ≈ yaux
        end

        @testset "Test Params Derivatives" begin

            function cal_Θ(theta)
                vec(predict(pred, theta))
            end

            dΘ_numeric = Calculus.finite_difference_jacobian(cal_Θ, Θ)

            ys = y_buffer(pred)
            dΘ = dΘ_buffer(pred)
            predict!(pred, ys, Θ; dΘ=dΘ)

            @test dΘ ≈ reshape(dΘ_numeric, size(dΘ))
        end

    end

    @testset "Test free-run simulation" begin
        # Test Basic Construction
        pred = free_run_simulation(mdl, yterms, uterms, IdData(y, u))

        @testset "Test Basic Construction" begin
            @test pred.dynmdl.yterms == [[1, 2]]
            @test pred.dynmdl.uterms == [[0, 1, 2]]
            @test pred.u == u
            @test pred.y0 == y0
            @test pred.length == size(yaux, 2)
            @test pred.start == 3
        end

        @testset "Test Prediction" begin
            @test predict(pred, Θ) ≈ yaux
        end

        @testset "Test Params Derivatives" begin

            function cal_Θ(theta)
                vec(predict(pred, theta))
            end

            dΘ_numeric = Calculus.finite_difference_jacobian(cal_Θ, Θ)

            ys = y_buffer(pred)
            dΘ = dΘ_buffer(pred)
            predict!(pred, ys, Θ; dΘ=dΘ)

            @test dΘ ≈ reshape(dΘ_numeric, size(dΘ))
        end

        @testset "Test Initial Conditions Derivatives" begin

            function cal_y0(x0)
                vec(predict(pred, Θ, reshape(x0, size(y0))))
            end

            dy0_numeric = Calculus.finite_difference_jacobian(cal_y0, vec(y0))

            ys = y_buffer(pred)
            dy0 = dy0_buffer(pred)
            predict!(pred, ys, Θ; dy0=dy0)

            @test dy0 ≈ reshape(dy0_numeric, size(dy0))
        end
    end
end

@testset "Test for MIMO data" begin
    # Generate Input Data
    len = 100
    srand(1)
    u = rand(2, len)-0.5;

    # Define Dynamic Model
    mdl = Linear(12, 2)
    yterms = [[1, 2, 5],[1, 2, 3]]
    uterms = [[2, 3, 4], [5, 6, 7]]
    diffeq = DifferenceEquation(mdl, yterms=yterms, uterms=uterms)

    # Simulate System for a given parameter
    Θ = rand(12*2)
    y0 = zeros(2, 5)
    ϕ = [Θ; vec(y0)]
    yaux = simulate(diffeq, u[:, 1:(end-2)], Θ, y0)
    y = [zeros(2, 2) y0 yaux]

    @testset "Test one-step-ahead prediction" begin
        pred = one_step_ahead(mdl, yterms, uterms, IdData(y, u))

        @testset "Test Basic Construction" begin
            @test pred.dynmdl.yterms == Vector{Vector{Int}}(0)
            @test pred.dynmdl.uterms == [[1, 2, 5],
                                         [1, 2, 3],
                                         [2, 3, 4],
                                         [5, 6, 7]]
            @test pred.u == [y; u][:, 1:(end-1)]
            @test pred.y0 == Matrix{Float64}(0, 0)
            @test pred.length == size(yaux, 2)
            @test pred.start == 8
        end

        @testset "Test Prediction" begin
            @test predict(pred, Θ) ≈ yaux
        end

        @testset "Test Params Derivatives" begin

            function cal_Θ(theta)
                vec(predict(pred, theta))
            end

            dΘ_numeric = Calculus.finite_difference_jacobian(cal_Θ, Θ)

            ys = y_buffer(pred)
            dΘ = dΘ_buffer(pred)
            predict!(pred, ys, Θ; dΘ=dΘ)

            @test dΘ ≈ reshape(dΘ_numeric, size(dΘ))
        end

    end

    @testset "Test free-run simulation" begin
        pred = free_run_simulation(mdl, yterms, uterms, IdData(y, u))

        @testset "Test Basic Construction" begin
            @test pred.dynmdl.yterms == [[1, 2, 5], [1, 2, 3]]
            @test pred.dynmdl.uterms == [[2, 3, 4], [5, 6, 7]]
            @test pred.u == u[:, 1:(end-2)]
            @test pred.y0 == y0
            @test pred.length == size(yaux, 2)
            @test pred.start == 8
        end

        @testset "Test Prediction" begin
            @test predict(pred, Θ) ≈ yaux
        end

        @testset "Test Params Derivatives" begin

            function cal_Θ(theta)
                vec(predict(pred, theta))
            end

            dΘ_numeric = Calculus.finite_difference_jacobian(cal_Θ, Θ)

            ys = y_buffer(pred)
            dΘ = dΘ_buffer(pred)
            predict!(pred, ys, Θ; dΘ=dΘ)

            @test dΘ ≈ reshape(dΘ_numeric, size(dΘ))
        end

        @testset "Test Initial Conditions Derivatives" begin

            function cal_y0(x0)
                vec(predict(pred, Θ, reshape(x0, size(y0))))
            end

            dy0_numeric = Calculus.finite_difference_jacobian(cal_y0, vec(y0))

            ys = y_buffer(pred)
            dy0 = dy0_buffer(pred)
            predict!(pred, ys, Θ; dy0=dy0)

            @test dy0 ≈ reshape(dy0_numeric, size(dy0))
        end
    end
end


text:
# 1D CA using Agents.jl
module CA1D
using Agents

mutable struct Cell <: AbstractAgent
  id::Int
  pos::Tuple{Int, Int}
  status::String
end

"""
    build_model(;rules::Dict, ncols=101)

Builds a 1D cellular automaton. `rules` is a dictionary with this format: `Dict("000" => "0")`. `ncols` is the number of cells in a 1D CA.
"""
function build_model(;rules::Dict, ncols::Integer=101)
  nv=(ncols,1)
  space = GridSpace(nv)
  properties = Dict(:rules => rules)
  model = ABM(Cell, space; properties = properties, scheduler=by_id)
  for n in 1:ncols
    add_agent!(Cell(n, (n,1), "0"), (n,1), model)
  end
  return model
end

function ca_step!(model)
  agentnum = nagents(model)
  new_status = Array{String}(undef, agentnum)
  for agid in 1:agentnum
    agent = model.agents[agid]
    center = agent.status
    if agent.id == agentnum
      # right = "0"  # this is for when there space is not toroidal
      right =  model.agents[1].status
      left = model.agents[agent.id-1].status
    elseif agent.id == 1
      # left = "0"  # this is for when there space is not toroidal
      left = model.agents[agentnum].status
      right = model.agents[agent.id+1].status
    else
      left = model.agents[agent.id-1].status
      right = model.agents[agent.id+1].status
    end
    rule = left*center*right
    new_status[agid] = model.properties[:rules][rule]
  end
  for ss in 1:agentnum
    model.agents[ss].status = new_status[ss]
  end
end

"""
    ca_run(model::ABM, runs::Integer)

Runs a 1D cellular automaton.
"""
function ca_run(model::ABM, runs::Integer)
  data = run!(model, dummystep, ca_step!, runs; agent_properties=[:pos, :status], when=1:runs)
end

end  # module


text:
module LintFilesTask

using Genie, Util, Lint, Logger

function description()
  """
  Lints the files in the indicated dir
  """
end

function run_task!()
  dir = joinpath(Pkg.dir("Genie"), "src")
  for filename in Util.walk_dir(dir)
    @show filename
    try
      for message in lintfile(filename)
        color = if iserror(message)
                  :red
                elseif iswarning(message)
                  :orange
                elseif isinfo(message)
                  :blue
                end
        print_with_color(color, string(message) * "\n\n")
      end
    catch ex
      Logger.log(string(ex), :err)
    end
  end
end

end

max_stars_repo_name: anirudhmurali/hello-julia-genie

max_issues_repo_name: anirudhmurali/hello-julia-genie

max_forks_repo_name: anirudhmurali/hello-julia-genie


text:
using AstroImages, FITSIO, Images, Random, Widgets, WCS, JLD

using Test, WCS
using SHA: sha256

import AstroImages: _float, render, _brightness_contrast, brightness_contrast

@testset "Conversion to float and fixed-point" begin
    @testset "Float" begin
        for T in (Float16, Float32, Float64)
            @test _float(T(-9.8)) === T(-9.8)
            @test _float(T(12.3)) === T(12.3)
        end
    end
    @testset "Integers" begin
        for (UIT, SIT) in ((UInt8,  Int8),
                           (UInt16, Int16),
                           (UInt32, Int32),
                           (UInt64, Int64))
            N = sizeof(UIT) * 8
            NT = Normed{UIT, N}
            maxint = UIT(big(2) ^ (N - 1))
            @test _float(typemin(UIT)) === _float(typemin(SIT)) === NT(0)
            @test _float(UIT(85)) === reinterpret(NT, UIT(85))
            @test _float(SIT(-85)) === _float(UIT(-85 + big(maxint)))
            @test _float(SIT(115)) === _float(UIT(115) + maxint)
            @test _float(typemax(UIT)) === _float(typemax(SIT)) === NT(1)
        end
    end
end

@testset "FITS and images" begin
    fname = tempname() * ".fits"
    for T in [UInt8, Int8, UInt16, Int16, UInt32, Int32, Int64,
              Float32, Float64]
        data = reshape(T[1:100;], 5, 20)
        FITS(fname, "w") do f
            write(f, data)
        end
        @test load(fname, 1)[1] == data
        @test load(fname, (1, 1))[1] == (data, data)
        img = AstroImage(fname)
        rendered_img = colorview(img)
        @test iszero(minimum(rendered_img))
    end
    rm(fname, force = true)
end

@testset "Control contrast" begin
    @test @inferred(_brightness_contrast(Gray, ones(Float32, 2, 2), 100, 100)) isa
          Array{Gray{Float32},2}
    Random.seed!(1)
    M = rand(2, 2)
    @test @inferred(_brightness_contrast(Gray, M, 100, 100)) ≈
        Gray.([0.48471895908315565  0.514787046406301;
               0.5280458879184159   0.39525854250610226]) rtol = 1e-12
    @test @inferred(_brightness_contrast(Gray, M, 0,   255)) == Gray.(M)
    @test @inferred(_brightness_contrast(Gray, M, 255,   0)) == Gray.(ones(size(M)))
    @test @inferred(_brightness_contrast(Gray, M,   0,   0)) == Gray.(zeros(size(M)))
    @test brightness_contrast(AstroImage(M)) isa Widgets.Widget{:manipulate,Any}
end

@testset "default handler" begin
    fname = tempname() * ".fits"
    @testset "less dimensions than 2" begin
        data = rand(2)
        FITS(fname, "w") do f
            write(f, data)
        end
        @test_throws ErrorException AstroImage(fname)
    end

    @testset "no ImageHDU" begin
        ## Binary table
        indata = Dict{String, Array}()
        i = length(indata) + 1
        indata["col$i"] = [randstring(10) for j=1:20]  # ASCIIString column
        i += 1
        indata["col$i"] = ones(Bool, 20)  # Bool column
        i += 1
        indata["col$i"] = reshape([1:40;], (2, 20))  # vector Int64 column
        i += 1
        indata["col$i"] = [randstring(5) for j=1:2, k=1:20]  # vector ASCIIString col
        indata["vcol"] = [randstring(j) for j=1:20]  # variable length column
        indata["VCOL"] = [collect(1.:j) for j=1.:20.] # variable length

        FITS(fname, "w") do f
            write(f, indata; varcols=["vcol", "VCOL"])
            @test_throws MethodError AstroImage(f)
        end
    end

    @testset "Opening AstroImage in different ways" begin
        data = rand(2,2)
        wcs = WCSTransform(2;)
        FITS(fname, "w") do f
            write(f, data)
        end
        f = FITS(fname)
        @test AstroImage(fname, 1) isa AstroImage
        @test AstroImage(Gray ,fname, 1) isa AstroImage
        @test AstroImage(Gray, f, 1) isa AstroImage
        @test AstroImage(data, wcs) isa AstroImage
        @test AstroImage((data,data), (wcs,wcs)) isa AstroImage
        @test AstroImage(Gray, data, wcs) isa AstroImage
        close(f)
    end

    @testset "Image HDU is not at 1st position" begin
        ## Binary table
        indata = Dict{String, Array}()
        i = length(indata) + 1
        indata["col$i"] = [randstring(10) for j=1:20]  # ASCIIString column
        i += 1
        indata["col$i"] = ones(Bool, 20)  # Bool column
        i += 1
        indata["col$i"] = reshape([1:40;], (2, 20))  # vector Int64 column
        i += 1
        indata["col$i"] = [randstring(5) for j=1:2, k=1:20]  # vector ASCIIString col
        indata["vcol"] = [randstring(j) for j=1:20]  # variable length column
        indata["VCOL"] = [collect(1.:j) for j=1.:20.] # variable length

        FITS(fname, "w") do f
            write(f, indata; varcols=["vcol", "VCOL"])
            write(f, rand(2, 2))
        end

        @test @test_logs (:info, "Image was loaded from HDU 3") AstroImage(fname) isa AstroImage
    end
    rm(fname, force = true)
end

@testset "Utility functions" begin
   @test size(AstroImage((rand(10,10), rand(10,10)))) == ((10,10), (10,10))
   @test length(AstroImage((rand(10,10), rand(10,10)))) == 2
end

@testset "multi image AstroImage" begin
    data1 = rand(10,10)
    data2 = rand(10,10)
    fname = tempname() * ".fits"
    FITS(fname, "w") do f
        write(f, data1)
        write(f, data2)
    end

    img = AstroImage(fname, (1,2))
    @test length(img.data) == 2
    @test img.data[1] == data1
    @test img.data[2] == data2

    f = FITS(fname)
    img = AstroImage(Gray, f, (1,2))
    @test length(img.data) == 2
    @test img.data[1] == data1
    @test img.data[2] == data2
    close(f)
end

@testset "multi wcs AstroImage" begin
    fname = tempname() * ".fits"
    f = FITS(fname, "w")
    inhdr = FITSHeader(["CTYPE1", "CTYPE2", "RADESYS", "FLTKEY", "INTKEY", "BOOLKEY", "STRKEY", "COMMENT",
                        "HISTORY"],
                       ["RA---TAN", "DEC--TAN", "UNK", 1.0, 1, true, "string value", nothing, nothing],
                       ["",
                        "",
                        "",
                        "floating point keyword",
                        "",
                        "boolean keyword",
                        "string value",
                        "this is a comment",
                        "this is a history"])

    indata = reshape(Float32[1:100;], 5, 20)
    write(f, indata; header=inhdr)
    write(f, indata; header=inhdr)
    close(f)

    img = AstroImage(fname, (1,2))
    f = FITS(fname)
    @test length(img.wcs) == 2
    @test WCS.to_header(img.wcs[1]) === WCS.to_header(WCS.from_header(read_header(f[1], String))[1])
    @test WCS.to_header(img.wcs[2]) === WCS.to_header(WCS.from_header(read_header(f[2], String))[1])

    img = AstroImage(Gray, f, (1,2))
    @test length(img.wcs) == 2
    @test WCS.to_header(img.wcs[1]) === WCS.to_header(WCS.from_header(read_header(f[1], String))[1])
    @test WCS.to_header(img.wcs[2]) === WCS.to_header(WCS.from_header(read_header(f[2], String))[1])
    close(f)
end

@testset "multi file AstroImage" begin
    fname1 = tempname() * ".fits"
    f = FITS(fname1, "w")
    inhdr = FITSHeader(["CTYPE1", "CTYPE2", "RADESYS", "FLTKEY", "INTKEY", "BOOLKEY", "STRKEY", "COMMENT",
                        "HISTORY"],
                    ["RA---TAN", "DEC--TAN", "UNK", 1.0, 1, true, "string value", nothing, nothing],
                    ["",
                        "",
                        "",
                        "floating point keyword",
                        "",
                        "boolean keyword",
                        "string value",
                        "this is a comment",
                        "this is a history"])

    indata1 = reshape(Int[1:100;], 5, 20)
    write(f, indata1; header=inhdr)
    close(f)

    fname2 = tempname() * ".fits"
    f = FITS(fname2, "w")
    indata2 = reshape(Int[1:100;], 5, 20)
    write(f, indata2; header=inhdr)
    close(f)

    fname3 = tempname() * ".fits"
    f = FITS(fname3, "w")
    indata3 = reshape(Int[1:100;], 5, 20)
    write(f, indata3; header=inhdr)
    close(f)

    img = AstroImage((fname1, fname2, fname3))
    f1 = FITS(fname1)
    f2 = FITS(fname2)
    f3 = FITS(fname3)

    @test length(img.data) == length(img.wcs) == 3
    @test img.data[1] == indata1
    @test img.data[2] == indata2
    @test img.data[3] == indata3
    @test WCS.to_header(img.wcs[1]) == WCS.to_header(img.wcs[2]) ==
        WCS.to_header(img.wcs[3]) == WCS.to_header(WCS.from_header(read_header(f1[1], String))[1])
    @test eltype(eltype(img.data)) == Int

    img = AstroImage(Gray, (f1, f2, f3), (1,1,1))
    @test length(img.data) == length(img.wcs) == 3
    @test img.data[1] == indata1
    @test img.data[2] == indata2
    @test img.data[3] == indata3
    @test WCS.to_header(img.wcs[1]) == WCS.to_header(img.wcs[2]) ==
        WCS.to_header(img.wcs[3]) == WCS.to_header(WCS.from_header(read_header(f1[1], String))[1])
    @test eltype(eltype(img.data)) == Int
    close(f1)
    close(f2)
    close(f3)
    rm(fname1, force = true)
    rm(fname2, force = true)
    rm(fname3, force = true)
end

include("plots.jl")
include("ccd2rgb.jl")


text:
using FileIO, Images

img = load("image.jpg")

sharpenkernel = reshape([-1.0, -1.0, -1.0, -1.0,  9.0, -1.0, -1.0, -1.0, -1.0], (3,3))

imfilt = imfilter(img, sharpenkernel)

save("imagesharper.png", imfilt)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   The Jacchia-Bowman 2008 Atmospheric Model, a product of the Space
#   Environment Technologies.
#
#   The source code is available at the following git:
#
#       http://sol.spacenvironment.net/jb2008/code.html
#
#   For more information about the model, see:
#
#       http://sol.spacenvironment.net/~JB2006/
#       http://sol.spacenvironment.net/~JB2008/
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Bowman, B. R., Tobiska, W. K., Marcos, F. A., Huang, C. Y., Lin, C. S.,
#       Burke, W. J (2008). A new empirical thermospheric density model JB2008
#       using new solar and geomagnetic indices. AIAA/AAS Astrodynamics
#       Specialist Conference, Honolulu, Hawaii.
#
#   [2] Bowman, B. R., Tobiska, W. K., Marcos, F. A., Valladares, C (2007). The
#       JB2006 empirical thermospheric density model. Journal of Atmospheric and
#       Solar-Terrestrial Physics, v. 70, p. 774-793.

#   [3] Jacchia, L. G (1970). New static models of the thermosphere and
#       exosphere with empirical temperature profiles. SAO Special Report #313.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

include("./jb2008_priv.jl")

export jb2008

"""
    function jb2008(JD::Number, glat::Number, glon::Number, h::Number)
    function jb2008(JD::Number, glat::Number, glon::Number, h::Number, F10::Number, F10ₐ::Number, S10::Number, S10ₐ::Number, M10::Number, M10ₐ::Number, Y10::Number, Y10ₐ::Number, DstΔTc::Number)

Compute the atmospheric density using the Jacchia-Bowman 2008 (JB2008) model.

If the space indices are not provided (first call), then they will be obtained
from the online database. In this case, the function `init_space_indices()` must
be called first and the function will throw an exception if the selected `JD` is
outside of the available data.

This model is a product of the **Space Environment Technologies**, more
information can be seen in the websites:

http://sol.spacenvironment.net/jb2006/

http://sol.spacenvironment.net/jb2008/

# Args

* `JD`: Julian day.
* `glat`: Geocentric latitude [rad].
* `glon`: Geocentric longitude [rad].
* `h`: Altitude [m].

* `F10`: 10.7-cm solar flux \\[10⁻²² W/(M² Hz)] (Tabular time 1 day earlier).
* `F10ₐ`: 10.7-cm averaged solar flux, 81-day centered on input time (Tabular
          time 1 day earlier).
* `S10`: EUV index (26-34 nm) scaled to F10.7 (Tabular time 1 day earlier).
* `S10ₐ`: EUV 81-day averaged centered index (Tabular time 1 day earlier).
* `M10`: MG2 index scaled to F10.7 (Tabular time 2 days earlier).
* `M10ₐ`: MG2 81-day averaged centered index (Tabular time 2 days earlier).
* `Y10`: Solar X-ray & Lya index scaled to F10.7 (Tabular time 5 days earlier).
* `Y10ₐ`: Solar X-ray & Lya 81-day averaged centered index (Tabular time 5 days
          earlier).
* `DstΔTc`: Temperature variation related to the Dst.

# Returns

An instance of the structure `JB2008_Output` with the computed values.

"""
function jb2008(JD::Number, glat::Number, glon::Number, h::Number)
    # Get the data in the desired Julian Day considering the tabular time of the
    # model.
    F10    = get_space_index(Val{:F10obs},  JD-1)
    F10ₐ   = get_space_index(Val{:F10Mobs}, JD-1)
    S10    = get_space_index(Val{:S10},     JD-1)
    S10ₐ   = get_space_index(Val{:S81a},    JD-1)
    M10    = get_space_index(Val{:M10},     JD-2)
    M10ₐ   = get_space_index(Val{:M81a},    JD-2)
    Y10    = get_space_index(Val{:Y10},     JD-5)
    Y10ₐ   = get_space_index(Val{:Y81a},    JD-5)
    DstΔTc = get_space_index(Val{:DstΔTc},  JD)

    jb2008(JD, glat, glon, h, F10, F10ₐ, S10, S10ₐ, M10, M10ₐ, Y10, Y10ₐ, DstΔTc)
end

function jb2008(JD::Number, glat::Number, glon::Number, h::Number, F10::Number,
                F10ₐ::Number, S10::Number, S10ₐ::Number, M10::Number,
                M10ₐ::Number, Y10::Number, Y10ₐ::Number, DstΔTc::Number)

    # Constants
    # ==========================================================================

    T₁    = 183.0       # Temperature at the lower bound [K].
    z₁    = 90.0        # Altitude of the lower bound [km].
    zx    = 125.0       # Altitude of the inflection point [km].
    Rstar = 8314.32     # Rstar is the universal gas-constant (mks)
                        # [joules/(K.kmol)].
    ρ₁    = 3.46e-6     # Density at `z₁` [kg/m³].
    A     = 6.02257e26  # Avogadro's constant (mks) [molecules/kmol].

    # Assumed sea-level composition.
    Mb₀  = 28.960
    q₀N₂ = 0.78110
    q₀O₂ = 0.20955
    q₀Ar = 9.3400e-3
    q₀He = 1.2890e-5

    # Molecular weights of each specie [kg/kmol].
    MN₂ = 28.0134
    MO₂ = 31.9988
    MO  = 15.9994
    MAr = 39.9480
    MHe =  4.0026
    MH  =  1.00797

    # Thermal diffusion coefficient for each specie.
    α_N₂ =  0.0
    α_O₂ =  0.0
    α_O  =  0.0
    α_Ar =  0.0
    α_He = -0.38
    α_H₂ =  0.0

    # Values used to stablish height step sizes in the integration process
    # between 90km to 105km, 105km to 500km, and above 500km.
    R1 = 0.010
    R2 = 0.025
    R3 = 0.075

    # Preliminaries
    # ==========================================================================

    # Convert the altitude from [m] to [km].
    h /= 1000

    # Compute the Sun position represented in the inertial reference frame
    # (MOD).
    s_i = sun_position_i(JD)

    # Compute the Sun declination [rad].
    δs = atan( s_i[3], sqrt(s_i[1]*s_i[1] + s_i[2]*s_i[2]) )

    # Compute the Sun right ascension [rad].
    Ωs = atan( s_i[2], s_i[1] )

    # Compute the right ascension of the selected location w.r.t. the inertial
    # reference frame.
    Ωp = glon + JDtoGMST(JD)

    # Compute the hour angle at the selected location, which is the angle
    # measured at the XY plane between the right ascension of the selected
    # position and the right ascension of the Sun.
    H = Ωp - Ωs

    # Compute the local solar time.
    lst = mod( (H + π)*12/π, 24 )

    # Algorithm
    # ==========================================================================

    # Eq. 2 [1], Eq. 14 [3]
    # --------------------------------------------------------------------------
    #
    # Nighttime minimum of the global exospheric temperature distribution when
    # the planetary geomagnetic index Kp is zero.
    ΔF10 =  F10 - F10ₐ
    ΔS10 =  S10 - S10ₐ
    ΔM10 = M10 - M10ₐ
    ΔY10 =  Y10 - Y10ₐ

    Wt = (F10ₐ/240)^(1/4)
    (Wt > 1) && (Wt = 1.0)
    Fsₐ = F10ₐ*Wt + S10ₐ*(1-Wt)
    Tc = 392.4 + 3.227Fsₐ + 0.298ΔF10 + 2.259ΔS10 + 0.312ΔM10 + 0.178ΔY10

    # Eq. 15 [3]
    # --------------------------------------------------------------------------
    η = abs(glat - δs)/2
    θ = abs(glat + δs)/2

    # Eq. 16 [3]
    # --------------------------------------------------------------------------
    τ = H - 0.64577182 + 0.10471976sin(H + 0.75049158)

    # Eq. 17 [3]
    # --------------------------------------------------------------------------
    m = 2.5
    n = 3.0
    R = 0.31
    C = cos(η)^m
    S = sin(θ)^m

    # NOTE: The original equation in [3] does not have the `abs` as in the
    # source-code of JB2008.
    Tl = Tc*( 1 + R*( S + (C-S)*abs( cos(τ/2) )^n ) )

    # Compute the correction to `Tc` considering the local solar time and
    # latitude.
    ΔTc = _jb2008_ΔTc(F10, lst, glat, h)

    # Compute the local exospheric temperature with the geomagnetic storm
    # effect.
    T_exo = Tl + DstΔTc
    T∞    = T_exo + ΔTc

    # Eq. 9 [3]
    # --------------------------------------------------------------------------
    #
    # Temperature at the inflection point `z = 125 km`.
    a  =  444.3807
    b  =    0.02385
    c  = -392.8292
    k  =   -0.0021357
    Tx = a + b*T∞ + c*exp(k*T∞)

    # Eq. 5 [3]
    # --------------------------------------------------------------------------
    #
    # From 90 to 105 km, for a given temperature profile `T[k]`, the density
    # ρ is computed by integrating the barometric equation:
    #
    #                  - _ -     _
    #                 |  M´ |    M´ g
    #   d lnρ´ = d ln | --- | - ------ dz ,
    #                 |  T  |    R* T
    #                  -   -
    #
    # which can be rewritten as:
    #
    #         -     -       _
    #        |  ρ´T  |      M´ g
    #   d ln | ----- | = - ------ dz .
    #        |   _   |      R* T
    #        |   M´  |
    #         -     -
    #
    #
    # Here, we will compute the following integral:
    #
    #
    #      / z₂
    #     |     _
    #     |     M´g
    #     |   ------- dz
    #     |      T
    #     |
    #    /   z₁
    #
    # in which z₁ is the minimum value between `h` and 105 km. The integration
    # will be computed by Newton-Cotes 4th degree method.

    z₂ = min(h, 105.0)

    int, z₂ = _jb2008_int(z₁, z₂, R1, Tx, T∞, _jb2008_δf1)

    Mb₁ = _jb2008_M(z₁)
    Tl₁ = _jb2008_T(z₁, Tx, T∞)
    Mb₂ = _jb2008_M(z₂)
    Tl₂ = _jb2008_T(z₂, Tx, T∞)

    # `Mbj` and `Tlj` contains, respectively, the mean molecular mass and local
    # temperature at the boundary of the integration interval.
    #
    # The `1000` factor is to convert `Rstar` to the appropriate units.
    ρ = ρ₁*(Mb₂/Mb₁)*(Tl₁/Tl₂)*exp(-1000int/Rstar)

    # Eq. 2 [3]
    # --------------------------------------------------------------------------
    NM = A*ρ
    N  = NM/Mb₂

    # Eq. 3 [3]
    # --------------------------------------------------------------------------
    NMoMb₀  = NM/Mb₀
    log_nN₂ = log(q₀N₂*NMoMb₀)
    log_nAr = log(q₀Ar*NMoMb₀)
    log_nHe = log(q₀He*NMoMb₀)

    # Eq. 4 [3]
    # --------------------------------------------------------------------------
    log_nO₂ = log(NMoMb₀ * (1 + q₀O₂) - N)
    log_nO  = log(2(N - NMoMb₀))
    log_nH  = 0.0

    Tz = 0.0
    if h <= 105
        Tz = Tl₂
        log_nH₂ = log_nHe - 25
    else
        # Eq.6 [3]
        # ----------------------------------------------------------------------
        #
        # From 100 km to 500 km, neglecting the thermal diffusion coefficient,
        # the eq. 6 in [3] can be written as:
        #
        #                          _
        #               1+αᵢ       M´ g
        #   d ln ( n . T    ) = - ------ dz ,
        #                          R* T
        #
        # where n number density for the i-th specie. This equations must be
        # integrated for each specie we are considering.
        #
        # Here, we will compute the following integral:
        #
        #
        #      / z₃
        #     |
        #     |     g
        #     |   ----- dz
        #     |     T
        #     |
        #    /   z₂
        #
        # in which z₃ is the minimum value between `h` and 500 km. The
        # integration will be computed by Newton-Cotes 4th degree method.

        z₃ = min(h, 500.0)

        int_1, z₃ = _jb2008_int(z₂, z₃, R1, Tx, T∞, _jb2008_δf2)

        Tl₃ = _jb2008_T(z₃, Tx, T∞)

        # If `h` is lower than 500 km, then keep integrating until 500 km due to
        # the hydrogen. Otherwise, continue the integration until `h`. Hence, we
        # will compute the following integral:
        #
        #      / z₄
        #     |
        #     |     g
        #     |   ----- dz
        #     |     T
        #     |
        #    /   z₃
        #
        # in which z₄ is the maximum value between `h` and 500 km. The
        # integration will be computed by Newton-Cotes 4th degree method.

        z₄ = max(h, 500.0)

        int_2, z₄ = _jb2008_int(z₃, z₄, (h <= 500) ? R2 : R3, Tx, T∞, _jb2008_δf2)

        Tl₄ = _jb2008_T(z₄, Tx, T∞)

        if h <= 500
            T500km    = Tl₄
            Tz        = Tl₃
            log_TfoTi = log(Tl₃/Tl₂)
            H_sign    = +1

            #          g
            # goRT = ------
            #         R*T
            goRT      = 1000int_1/Rstar
        else
            T500km    = Tl₃
            Tz        = Tl₄
            log_TfoTi = log(Tl₄/Tl₂)
            H_sign    = -1

            #          g
            # goRT = ------
            #         R*T
            goRT      = 1000(int_1 + int_2)/Rstar
        end

        log_nN₂ += -(1 + α_N₂)*log_TfoTi - goRT*MN₂
        log_nO₂ += -(1 + α_O₂)*log_TfoTi - goRT*MO₂
        log_nO  += -(1 + α_O )*log_TfoTi - goRT*MO
        log_nAr += -(1 + α_Ar)*log_TfoTi - goRT*MAr
        log_nHe += -(1 + α_He)*log_TfoTi - goRT*MHe

        # Eq. 7 [3]
        # ----------------------------------------------------------------------
        #
        # The equation of [3] will be converted from `log10` to `log`.
        # Furthermore, the units will be converted to [1/cm³] to [1/m³].

        log10_nH_500km = @evalpoly(log10(T∞), 73.13, -39.40, 5.5)
        log_nH_500km   = log(10) * ( log10_nH_500km + 6 )
        #                                             ^
        #                                             |
        #                         This factor converts from [1/cm³] to [1/m³]

        # Compute the hydrogen density based on the value at 500km.
        log_nH = log_nH_500km + H_sign *
            ( log(Tl₄/Tl₃) + ( 1000/Rstar ) * int_2 * MH )
    end

    # Eq. 24 [3] - Seasonal-latitudinal variation
    # ----------------------------------------------------------------------
    #
    # TODO: The term related to the year in the source-code of JB 2008 is
    # different from [3]. This must be verified.

    #         | Modified JD  |
    Φ = mod( ( JD - 2400000.5 - 36204 )/365.2422, 1 )

    Δlog₁₀ρ = 0.02 * ( h - 90 ) * sign(glat) * exp( -0.045( h - 90 ) ) *
        sin(glat)^2 * sin( 2π*Φ + 1.72 )

    # Convert from `log10` to `log`.
    Δlogρ = log(10)*Δlog₁₀ρ

    # Eq. 23 [3] - Semiannual variation
    # ----------------------------------------------------------------------

    if h < 2000
        # Compute the year given the selected Julian Day.
        year, month, day, = JDtoDate(JD)

        # Compute the day of the year.
        doy = JD - DatetoJD(year, 1, 1, 0, 0, 0) + 1

        # Use the new semiannual model from [1].
        Fz, Gz, Δsalog₁₀ρ = _jb2008_semiannual(doy, h, F10ₐ, S10ₐ, M10ₐ)

        (Fz < 0) && (Δsalog₁₀ρ = 0.0)

        # Convert from `log10` to `log`.
        Δsalogρ = log(10)*Δsalog₁₀ρ
    else
        Δsalogρ = 0.0
    end

    # Compute the total variation.
    Δρ = Δlogρ + Δsalogρ

    # Apply to the number densities.
    log_nN₂ += Δρ
    log_nO₂ += Δρ
    log_nO  += Δρ
    log_nAr += Δρ
    log_nHe += Δρ
    log_nH  += Δρ

    # Compute the high altitude exospheric density correction factor.
    #
    # The source-code computes this correction after computing the mass
    # density and mean molecular weight. Hence, the correction is applied
    # only to the total density. However, we will apply the correction to
    # each specie. TODO: Verify if this is reasonable.
    log_FρH  = log(_jb2008_highaltitude(h, F10ₐ))
    log_nN₂ += log_FρH
    log_nO₂ += log_FρH
    log_nO  += log_FρH
    log_nAr += log_FρH
    log_nHe += log_FρH
    log_nH  += log_FρH

    # Compute the mass density and mean molecular weight and convert number
    # density logs from natural to common.
    sum_n  = 0.0
    sum_mn = 0.0

    for (log_n, M) in [ (log_nN₂, MN₂), (log_nO₂, MO₂), (log_nO, MO),
                        (log_nAr, MAr), (log_nHe, MHe), (log_nH, MH) ]
        n = exp(log_n)
        sum_n  += n
        sum_mn += n*M
    end

    nN₂ = exp(log_nN₂)
    nO₂ = exp(log_nO₂)
    nO  = exp(log_nO)
    nAr = exp(log_nAr)
    nHe = exp(log_nHe)
    nH  = exp(log_nH)
    ρ   = sum_mn/A

    # Create and return the output structure.
    JB2008_Output{Float64}(nN₂, nO₂, nO, nAr, nHe, nH, ρ, T_exo, Tz)
end


text:
struct Node{R <: AbstractFloat, N}
    ind::CartesianIndex{N}
    val::R
end

isless(a::Node, b::Node) = isless(a.val, b.val)


text:
using LightGraphs

"""
    degreeisgreater(graph::AbstractGraph, vertex, deg)

Check if the degree of `vertex` is greater than `deg`.
"""
function degreeisgreater(graph::AbstractGraph, vertex, deg)
    return degree(graph)[vertex] >= deg
end

"""
    degreeislesser(graph::AbstractGraph, vertex, deg)

Check if the degree of `vertex` is lesser than `deg`.
"""
function degreeislesser(graph::AbstractGraph, vertex, deg)
    return degree(graph)[vertex] <= deg
end

"""
    degreeis(graph::AbstractGraph, vertex, degrees)

Check if the degree of `vertex` is contained in `degrees`.
"""
function degreeis(graph::AbstractGraph, vertex, degrees)
    return degree(graph)[vertex] in degrees
end

"""
    getmaxsize(arrays)

Return the maximum size of subarrays contained in `arrays`.
"""
function getmaxsize(arrays)
    try
        return maximum([length(array) for array in arrays])
    catch err
        if isa(err, ArgumentError)
            return -1
        end
    end
end

"""
    getminsize(arrays)

Return the minimum size of subarrays contained in `arrays`.
"""
function getminsize(arrays)
    try
        return minimum([length(array) for array in arrays])
    catch err
        if isa(err, ArgumentError)
            return -1
        end
    end
end

"""
    getsubset(arrays, values)

Return the subset of arrays contained in `arrays` wich contains `values`.
"""
function getsubset(arrays, values)
    return filter(array -> issubset(values, array), arrays)
end

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
abstract type NewtonParams <: SolverParams end

abstract type PreallocatedDataNewton{T <: Real, S} <: PreallocatedData{T, S} end

abstract type PreallocatedDataNewtonKrylov{T <: Real, S} <: PreallocatedDataNewton{T, S} end

include("K3Krylov.jl")
include("K3SKrylov.jl")
include("K3_5Krylov.jl")

abstract type PreallocatedDataNewtonStructured{T <: Real, S} <: PreallocatedDataNewton{T, S} end

# utils for K3_5 gpmr
include("K3_5gpmr_utils.jl")
include("K3SStructured.jl")
include("K3_5Structured.jl")


text:
# # An unstable Bickley jet in Shallow Water model
#
# This example uses Oceananigans.jl's `ShallowWaterModel` to simulate
# the evolution of an unstable, geostrophically balanced, Bickley jet
# The example is periodic in ``x`` with flat bathymetry and
# uses the conservative formulation of the shallow water equations.
# The initial conditions superpose the Bickley jet with small-amplitude perturbations.
# See ["The nonlinear evolution of barotropically unstable jets," J. Phys. Oceanogr. (2003)](https://doi.org/10.1175/1520-0485(2003)033<2173:TNEOBU>2.0.CO;2)
# for more details on this problem.
#
# The mass transport ``(uh, vh)`` is the prognostic momentum variable
# in the conservative formulation of the shallow water equations,
# where ``(u, v)`` are the horizontal velocity components and ``h``
# is the layer height.
#
# ## Install dependencies
#
# First we make sure that we have all of the packages that are required to
# run the simulation.
#
# ```julia
# using Pkg
# pkg"add Oceananigans, NCDatasets, Polynomials, CairoMakie"
# ```

using Oceananigans
using Oceananigans.Models: ShallowWaterModel

# ## Two-dimensional domain
#
# The shallow water model is a two-dimensional model and thus the number of vertical
# points `Nz` must be set to one.  Note that ``L_z`` is the mean depth of the fluid.

Lx, Ly, Lz = 2π, 20, 1
Nx, Ny = 128, 128

grid = RectilinearGrid(size = (Nx, Ny),
                       x = (0, Lx), y = (-Ly/2, Ly/2),
                       topology = (Periodic, Bounded, Flat))

# ## Physical parameters
#
# We choose non-dimensional parameters

const U = 1.0         # Maximum jet velocity

f = 1           # Coriolis parameter
g = 9.8         # Gravitational acceleration
Δη = f * U / g  # Maximum free-surface deformation as dictated by geostrophy

# ## Building a `ShallowWaterModel`
#
# We build a `ShallowWaterModel` with the `WENO5` advection scheme and
# 3rd-order Runge-Kutta time-stepping,

model = ShallowWaterModel(
                          timestepper = :RungeKutta3,
                          advection = WENO5(),
                          grid = grid,
                          gravitational_acceleration = g,
                          coriolis = FPlane(f=f))

# Use `architecture = GPU()` to run this problem on a GPU.

# ## Background state and perturbation
#
# The background velocity ``ū`` and free-surface ``η̄`` correspond to a
# geostrophically balanced Bickely jet with maximum speed of ``U`` and maximum
# free-surface deformation of ``Δη``,

h̄(x, y, z) = Lz - Δη * tanh(y)
ū(x, y, z) = U * sech(y)^2

# The total height of the fluid is ``h = L_z + \eta``. Linear stability theory predicts that
# for the parameters we consider here, the growth rate for the most unstable mode that fits
# our domain is approximately ``0.139``.

# The vorticity of the background state is

ω̄(x, y, z) = 2 * U * sech(y)^2 * tanh(y)

# The initial conditions include a small-amplitude perturbation that decays away from the
# center of the jet.

small_amplitude = 1e-4

 uⁱ(x, y, z) = ū(x, y, z) + small_amplitude * exp(-y^2) * randn()
uhⁱ(x, y, z) = uⁱ(x, y, z) * h̄(x, y, z)

# We first set a "clean" initial condition without noise for the purpose of discretely
# calculating the initial 'mean' vorticity,

ū̄h(x, y, z) = ū(x, y, z) * h̄(x, y, z)

set!(model, uh = ū̄h, h = h̄)

# We next compute the initial vorticity and perturbation vorticity,

uh, vh, h = model.solution

## Build velocities
u = uh / h
v = vh / h

## Build and compute mean vorticity discretely
ω = Field(∂x(v) - ∂y(u))
compute!(ω)

## Copy mean vorticity to a new field
ωⁱ = Field{Face, Face, Nothing}(model.grid)
ωⁱ .= ω

## Use this new field to compute the perturbation vorticity
ω′ = Field(ω - ωⁱ)

# and finally set the "true" initial condition with noise,

set!(model, uh = uhⁱ)

# ## Running a `Simulation`
#
# We pick the time-step so that we make sure we resolve the surface gravity waves, which
# propagate with speed of the order ``\sqrt{g L_z}``. That is, with `Δt = 1e-2` we ensure
# that `` \sqrt{g L_z} Δt / Δx,  \sqrt{g L_z} Δt / Δy < 0.7``.

simulation = Simulation(model, Δt = 1e-2, stop_time = 150)

# ## Prepare output files
#
# Define a function to compute the norm of the perturbation on the cross channel velocity.
# We obtain the `norm` function from `LinearAlgebra`.

using LinearAlgebra: norm

perturbation_norm(args...) = norm(v)

# Build the `output_writer` for the two-dimensional fields to be output.
# Output every `t = 1.0`.

simulation.output_writers[:fields] = NetCDFOutputWriter(model, (; ω, ω′),
                                                        filename = joinpath(@__DIR__, "shallow_water_Bickley_jet_fields.nc"),
                                                        schedule = TimeInterval(1),
                                                        overwrite_existing = true)

# Build the `output_writer` for the growth rate, which is a scalar field.
# Output every time step.

simulation.output_writers[:growth] = NetCDFOutputWriter(model, (; perturbation_norm),
                                                        filename = joinpath(@__DIR__, "shallow_water_Bickley_jet_perturbation_norm.nc"),
                                                        schedule = IterationInterval(1),
                                                        dimensions = (; perturbation_norm = ()),
                                                        overwrite_existing = true)

# And finally run the simulation.

run!(simulation)

# ## Visualize the results

# Load required packages to read output and plot.

using NCDatasets, Printf, CairoMakie
nothing # hide

# Define the coordinates for plotting.

x, y = xnodes(ω), ynodes(ω)
nothing # hide

# Read in the `output_writer` for the two-dimensional fields and then create an animation
# showing both the total and perturbation vorticities.

ds = NCDataset(simulation.output_writers[:fields].filepath, "r")

times = ds["time"][:]

n = Observable(1)

ω = @lift ds["ω"][:, :, 1, $n]
ω′ = @lift ds["ω′"][:, :, 1, $n]

ω′_lims = @lift (-maximum(abs, ds["ω′"][:, 1, :, $n]), maximum(abs, ds["ω′"][:, 1, :, $n]))

title = @lift @sprintf("t = %.1f", times[$n])

fig = Figure(resolution = (800, 440))

axis_kwargs = (xlabel = "x",
               ylabel = "y",
               aspect = AxisAspect(1),
               limits = ((0, Lx), (-Ly/2, Ly/2)))

ax_ω = Axis(fig[2, 1]; title = "total vorticity, ω", axis_kwargs...)
ax_ω′ = Axis(fig[2, 3]; title = "perturbation vorticity, ω - ω̄", axis_kwargs...)

hm_ω = heatmap!(ax_ω, x, y, ω, colorrange = (-1, 1), colormap = :balance)
Colorbar(fig[2, 2], hm_ω)

hm_ω′ = heatmap!(ax_ω′, x, y, ω′, colormap = :balance)
Colorbar(fig[2, 4], hm_ω′)

fig[1, :] = Label(fig, title, textsize=24, tellwidth=false)

# Finally, we record a movie.

frames = 1:length(times)

record(fig, "shallow_water_Bickley_jet.mp4", frames, framerate=12) do i
       @info "Plotting frame $i of $(frames[end])..."
       n[] = i
end
nothing #hide

# ![](shallow_water_Bickley_jet.mp4)

# It's always good practice to close the NetCDF files when we are done.

close(ds)
nothing # hide

# Read in the `output_writer` for the scalar field (the norm of ``v``-velocity).

ds2 = NCDataset(simulation.output_writers[:growth].filepath, "r")

     t = ds2["time"][:]
norm_v = ds2["perturbation_norm"][:]

close(ds2)
nothing # hide

# We import the `fit` function from `Polynomials.jl` to compute the best-fit slope of the
# perturbation norm on a logarithmic plot. This slope corresponds to the growth rate.

using Polynomials: fit

I = 6000:7000

degree = 1
linear_fit_polynomial = fit(t[I], log.(norm_v[I]), degree, var = :t)

# We can get the coefficient of the ``n``-th power from the fitted polynomial by using `n`
# as an index, e.g.,

constant, slope = linear_fit_polynomial[0], linear_fit_polynomial[1]

# We then use the computed linear fit coefficients to construct the best fit and plot it
# together with the time-series for the perturbation norm for comparison.

best_fit = @. exp(constant + slope * t)

lines(t, norm_v;
      linewidth = 4,
      label = "norm(v)",
      axis = (yscale = log10,
              limits = (nothing, (1e-3, 30)),
              xlabel = "time",
              ylabel = "norm(v)",
               title = "growth of perturbation norm"))

lines!(t[I], 2 * best_fit[I]; # factor 2 offsets fit from curve for better visualization
       linewidth = 4,
       label = "best fit")

axislegend()

fig # hide

# The slope of the best-fit curve on a logarithmic scale approximates the rate at which instability
# grows in the simulation. Let's see how this compares with the theoretical growth rate.

println("Numerical growth rate is approximated to be ", round(slope, digits=3), ",\n",
        "which is very close to the theoretical value of 0.139.")


text:
#InspectDR: Base functionnality and types for Gtk layer
#-------------------------------------------------------------------------------

import Gtk: getproperty, setproperty!, signal_connect, @guarded

#Create module aliasese to access constants:
#TODO: Figure out why constants are Int32 instead of Int???
import Gtk: GdkKeySyms
import Gtk: GtkPositionType, GdkEventType, GdkEventMask
import Gtk: GConstants.GtkShadowType

#==Extensions
===============================================================================#
#TODO: Why can't I just access these functions directly from Gtk/gen/gbox3???

function draw_value(scale::Gtk.GtkScale, draw_value_::Bool)
	ccall((:gtk_scale_set_draw_value,Gtk.libgtk),Void,(Ptr{Gtk.GObject},Cint),scale,draw_value_)
	return scale
end
function value_pos(scale::Gtk.GtkScale,pos::Int)
	ccall((:gtk_scale_set_value_pos,Gtk.libgtk),Void,(Ptr{Gtk.GObject},Cint),scale,pos)
	return scale
end

function can_focus(widget::Gtk.GtkWidget,can_focus_::Bool)
	ccall((:gtk_widget_set_can_focus,Gtk.libgtk),Void,(Ptr{Gtk.GObject},Cint),widget,can_focus_)
	return widget
end

function window_close(window::Gtk.GtkWindow)
	ccall((:gtk_window_close,Gtk.libgtk),Void,(Ptr{Gtk.GObject},),window)
	return
end

function focus(window::Gtk.GtkWindow,widget)
	ccall((:gtk_window_set_focus,Gtk.libgtk),Void,(Ptr{Gtk.GObject},Ptr{Gtk.GObject}),window,widget)
	return window
end
function focus(widget)
	ccall((:gtk_widget_grab_focus,Gtk.libgtk),Void,(Ptr{Gtk.GObject},),widget)
	return widget
end
#=
function activate_key(wnd::Gtk.GtkWindow,event::GdkEventKey)
	ccall((:gtk_window_activate_key,Gtk.libgtk),Cint,(Ptr{Gtk.GObject},Ptr{Gtk.GObject}),wnd,event)
	return wnd
end
#gboolean gtk_window_activate_key (GtkWindow *window, GdkEventKey *event);
=#

#==Constants
===============================================================================#
const XAXIS_SCALEMAX = 1000
const XAXIS_POS_STEPRES = 1/500

#==Display types
===============================================================================#
#Generic type used to spawn new InspectDR display windows:
immutable GtkDisplay <: Display
end

#==Main types
===============================================================================#

abstract InputState #Identifies current user input state.

type GtkSelection
	enabled::Bool
	bb::BoundingBox
	ext_start::PExtents2D #Exetents @ start of operation
	#Store ext_start to avoid accumulation of numerical errors.
end
GtkSelection() = GtkSelection(false, BoundingBox(0,0,0,0), PExtents2D())

type PlotWidget
	widget::_Gtk.Box #Base widget
	canvas::_Gtk.Canvas #Actual plot area
	src::Plot
	graphbb::BoundingBox
	state::InputState

	#Scrollbars to control x-scale & position:
	w_xscale::_Gtk.Scale
	xscale::_Gtk.Adjustment
	w_xpos::_Gtk.Scale
	xpos::_Gtk.Adjustment

	#Display image (Cached):
	bufsurf::Cairo.CairoSurface
#	bufbb::BoundingBox

	sel::GtkSelection

	#Restrict h/v motion:
	hallowed::Bool
	vallowed::Bool

	#Event handlers:
	eh_plothover::NullOr{HandlerInfo}
end

#Supports multiplot:
type GtkPlot
	wnd::_Gtk.Window
	grd::_Gtk.Grid #Holds subplot widgets
	subplots::Vector{PlotWidget}
	properties::Multiplot #Properties: ".subplots" ignored
	status::_Gtk.Label
end

#==Mutators
===============================================================================#
function settitle(::Type{GtkPlot}, wnd::_Gtk.Window, title::DisplayStringArg)
	if length(title)> 0
		title = "InspectDR - $(title)"
	else
		title = "InspectDR"
	end
	Gtk.setproperty!(wnd, :title, title)
end

function settitle(gplot::GtkPlot, title::DisplayStringArg)
	gplot.properties.title = DisplayString(title)
	settitle(GtkPlot, gplot.wnd, gplot.properties.title)
end

#==Helper functions
===============================================================================#
function copy_properties(mp::Multiplot)
	return Multiplot(mp.title, [], mp.ncolumns,
		mp.wplot, mp.hplot, mp.htitle, mp.fnttitle
	)
end

#==Main functions
===============================================================================#

function invalidbuffersize(pwidget::PlotWidget)
	return width(pwidget.canvas) != width(pwidget.bufsurf) ||
		height(pwidget.canvas) != height(pwidget.bufsurf)
end

#Render PlotWidget widget to buffer:
#-------------------------------------------------------------------------------
function render(pwidget::PlotWidget)
	#Create new buffer large enough to match canvas:
	#TODO: Is crating surfaces expensive?  This solution might be bad.
	if invalidbuffersize(pwidget)
		Cairo.destroy(pwidget.bufsurf)
		#TODO: use RGB surface? Gtk.cairo_surface_for() appears to generate ARGB surface (slower?)
		#pwidget.bufsurf = Cairo.CairoRGBSurface(width(pwidget.canvas),height(pwidget.canvas))
		pwidget.bufsurf = Gtk.cairo_surface_for(pwidget.canvas) #create similar
	end

	w = width(pwidget.canvas); h = height(pwidget.canvas)
	bb = BoundingBox(0, w, 0, h)
	ctx = CairoContext(pwidget.bufsurf)

	_reset(ctx)
	clear(ctx, bb)
	render(ctx, pwidget.src, bb)
	pwidget.graphbb = graphbounds(bb, pwidget.src.layout, pwidget.src.axes)
	Cairo.destroy(ctx)
end

#==IO functions
===============================================================================#
#_write() GtkPlot: Auto-coumpute w/h
function _write(path::AbstractString, mime::MIME, gplot::GtkPlot)
	mplot = copy_properties(gplot.properties)
	for s in gplot.subplots
		push!(mplot.subplots, s.src)
	end
	_write(path, mime, mplot)
end

write_png(path::AbstractString, gplot::GtkPlot) = _write(path, MIMEpng(), gplot)
write_svg(path::AbstractString, gplot::GtkPlot) = _write(path, MIMEsvg(), gplot)
write_eps(path::AbstractString, gplot::GtkPlot) = _write(path, MIMEeps(), gplot)
write_pdf(path::AbstractString, gplot::GtkPlot) = _write(path, MIMEpdf(), gplot)

#Last line


text:
using LinearAlgebra
using Distances, StaticArrays, StatsFuns

if true
x = rand(10000, 100000)
x = x'*x+UniformScaling(1.0)
end

#cholesky 10000x10000 takes 5-6 sec
# matmul two 100000 matrices takes 30 secr
# matvec takes 0.1 sec
#cholesky solve 0.15 sec

#forming kernel matrix takes 2.217 seconds when s is 1D, 3 seconds when s is 7D
#0.5 seconds to take dot product of 10,000^2 matrix with another matrix
#.4 seconds to square entries of 10,000^2 matrix

preallocated = zeros(10000,10000)

#15000x15000 case
#x*x 50 sec O(n^3)
#x*z 0.2 sec
# x hadamard product x: 1 second
#cholesky(x) solve: .3 seconds
#cholesky: 17 seconds O(n^3)

# ------------takeaway here is the math adds up-------------
#fastK(s, s):  7 seconds
#pairwise(s, s): 2 seconds
#exponentiate the result^: 2seconds
#divide by pi: 1.2 seconds

#0 pairwise: 2 seconds
#1.squaring operation ~1 second
#2. scale by - theta operation ~1 second
# 3. exponentiate operation ~2seconds
#divide by constant operation ~1 sec

#20000 x 20000 case
#fastK:  17 - 33 sec
#cholesky: 61 sec
#cholesky solve: 0.5 sec
#matvec with kernel matrix: 0.3 sec
#dot product k .* k: 1.3 sec

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
#
println("Eb) Three-step cascade computation and simulation for the photo-ionization of Si^- and its subsequent decay: AverageSCA model.")
using JLD2
#
println("aa")
setDefaults("method: continuum, asymptotic Coulomb")    ## setDefaults("method: continuum, Galerkin")
setDefaults("method: normalization, pure sine")         ## setDefaults("method: normalization, pure Coulomb")    setDefaults("method: normalization, pure sine")
#
if true
setDefaults("print summary: open", "zzz-Cascade-computation-photoionization.sum")

name = "Photoionization of Si- "
grid = Radial.Grid(false)
wa   = Cascade.Computation(Cascade.Computation(); name=name, nuclearModel=Nuclear.Model(14.), grid=grid, approach=Cascade.AverageSCA(),
                           scheme=Cascade.PhotonIonizationScheme([Photo], 1, [30.0, 80.0]),
                           initialConfigs=[Configuration("1s^2 2s^2 2p^6 3s^2 3p^3")] )
println(wa)
@show name
wb = perform(wa; output=true)
setDefaults("print summary: close", "")

elseif  false    ## Stepwise decay cascade
setDefaults("print summary: open", "zzz-Cascade-computation-following-decay.sum")

using JLD2
JLD2.@load "zzz-cascade-ionizing-computations-2020-02-02T20.jld"
iniMultiplets = results["generated multiplets:"]

name = "Si- (1s^-1) decay cascade"
grid = Radial.Grid(false)
wa   = Cascade.Computation(Cascade.Computation(); name=name, nuclearModel=Nuclear.Model(14.), grid=grid, approach=Cascade.AverageSCA(),
                           scheme=Cascade.StepwiseDecayScheme([Auger, Radiative], 1, 0, Shell[], Shell[]),
                           initialMultiplets=iniMultiplets )
println(wa)
@show name
wb = perform(wa; output=true)
setDefaults("print summary: close", "")

else
setDefaults("print summary: open", "zzz-Cascade-simulation.sum")

JLD2.@load "zzz-cascade-ionizing-computations-2020-02-02T20.jld"
resIon  = results
JLD2.@load "zzz-cascade-decay-computations-2020-02-02T20.jld"
resDecay = results

data = [resIon, resDecay]
name = "Simulation after Si- 1s and 2s ionization and subsequent decay"

wc   = Cascade.Simulation(Cascade.Simulation(), name=name, properties=Cascade.AbstractSimulationProperty[Cascade.IonDistribution()], ## , Cascade.FinalLevelDistribution()],
                          settings=Cascade.SimulationSettings(0., 0., 0., 0., 30., [(78, 2.0), (79, 1.0), (80, 0.5)]), computationData=data )
println(wc)
wd = perform(wc; output=true)
setDefaults("print summary: close", "")

end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: migrationhub_config
using AWS.Compat
using AWS.UUIDs

"""
    create_home_region_control(home_region, target)
    create_home_region_control(home_region, target, params::Dict{String,<:Any})

This API sets up the home region for the calling account only.

# Arguments
- `home_region`: The name of the home region of the calling account.
- `target`: The account for which this command sets up a home region control. The Target is
  always of type ACCOUNT.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"DryRun"`: Optional Boolean flag to indicate whether any effect should take place. It
  tests whether the caller has permission to make the call.
"""
function create_home_region_control(
    HomeRegion, Target; aws_config::AbstractAWSConfig=global_aws_config()
)
    return migrationhub_config(
        "CreateHomeRegionControl",
        Dict{String,Any}("HomeRegion" => HomeRegion, "Target" => Target);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_home_region_control(
    HomeRegion,
    Target,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return migrationhub_config(
        "CreateHomeRegionControl",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}("HomeRegion" => HomeRegion, "Target" => Target),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    describe_home_region_controls()
    describe_home_region_controls(params::Dict{String,<:Any})

This API permits filtering on the ControlId and HomeRegion fields.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"ControlId"`: The ControlID is a unique identifier string of your HomeRegionControl
  object.
- `"HomeRegion"`: The name of the home region you'd like to view.
- `"MaxResults"`: The maximum number of filtering results to display per page.
- `"NextToken"`: If a NextToken was returned by a previous call, more results are
  available. To retrieve the next page of results, make the call again using the returned
  token in NextToken.
- `"Target"`: The target parameter specifies the identifier to which the home region is
  applied, which is always of type ACCOUNT. It applies the home region to the current ACCOUNT.
"""
function describe_home_region_controls(; aws_config::AbstractAWSConfig=global_aws_config())
    return migrationhub_config(
        "DescribeHomeRegionControls"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function describe_home_region_controls(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return migrationhub_config(
        "DescribeHomeRegionControls",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    get_home_region()
    get_home_region(params::Dict{String,<:Any})

Returns the calling account’s home region, if configured. This API is used by other AWS
services to determine the regional endpoint for calling AWS Application Discovery Service
and Migration Hub. You must call GetHomeRegion at least once before you call any other AWS
Application Discovery Service and AWS Migration Hub APIs, to obtain the account's Migration
Hub home region.

"""
function get_home_region(; aws_config::AbstractAWSConfig=global_aws_config())
    return migrationhub_config(
        "GetHomeRegion"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function get_home_region(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return migrationhub_config(
        "GetHomeRegion", params; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end


text:
# Address map interfaces
function fill_array!(val::T, arr::Vector{T}, f_ind::Int) where T
    if length(arr) < f_ind
        resize!(arr, 2 * f_ind)
    end
    arr[f_ind] = val
    1
end

function fill_array!(val::Vector{T}, arr::Vector{T}, f_ind::Int) where T
    if length(arr) < f_ind + length(val)
        resize!(arr, 2 * (f_ind + length(val)))
    end
    arr[f_ind : f_ind + length(val) - 1] = val
    length(val)
end

function fill_array!(val::Value, arr::Vector{T}, f_ind::Int) where T
    fill_array!(get_value(val), arr, f_ind)
end

function fill_array!(val::Choice, arr::Vector{T}, f_ind::Int) where T
    fill_array!(get_value(val), arr, f_ind)
end

function fill_array!(val::T, arr::Vector{K}, f_ind::Int) where {K, T <: AddressMap}
    sorted_keys = sort(collect(keys(shallow_iterator(val))))
    idx = f_ind
    for k in sorted_keys
        n = fill_array!(get_sub(val, k), arr, idx)
        idx += n
    end
    idx - f_ind
end

function array(gs::K, ::Type{T}) where {T, K}
    arr = Vector{T}(undef, 32)
    n = fill_array!(gs, arr, 1)
    resize!(arr, n)
    arr
end

function from_array(::T, arr::Vector{T}, f_ind::Int) where T
    (1, arr[f_ind])
end

function from_array(val::Vector{T}, arr::Vector{T}, f_ind::Int) where {T}
    n_read = length(val)
    (n_read, arr[f_ind : f_ind + n_read - 1])
end

function from_array(::SelectAll, arr::Vector{K}, f_ind::Int) where K
    (1, Value(arr[f_ind]))
end

function from_array(v::Value, arr::Vector{K}, f_ind::Int) where K
    shape = get_value(v)
    n, v = from_array(shape, arr, f_ind)
    (n, Value(v))
end

function from_array(c::Choice, arr::Vector{K}, f_ind::Int) where K
    shape = get_value(c)
    n, v = from_array(shape, arr, f_ind)
    (n, Value(v))
end

function from_array(schema::T, arr::Vector{K}, f_ind::Int) where {K, T <: DynamicMap}
    sel = DynamicMap{Value}()
    sorted_keys = sort(collect(keys(shallow_iterator(schema))))
    idx = f_ind
    for k in sorted_keys
        (n, v) = from_array(get_sub(schema, k), arr, idx)
        idx += n
        set_sub!(sel, k, v)
    end
    (idx - f_ind, sel)
end

function target(schema::K, arr::Vector) where K
    (n, sel) = from_array(schema, arr, 1)
    n != length(arr) && error("Dimension error: length of arr $(length(arr)) must match $n.")
    sel
end


text:
module Pseudoseq

export Sequencing

include("PuzzleMaker.jl")
include("sequencing/Sequencing.jl")

end # module


text:
"""
    ControlAffineSystem

Models a nonlinear control affine system.

# Fields
- `n::Int`: state dimension
- `m::Int`: control dimension
- `f`: function f(x) that models the drift dynamics
- `g`: function g(x) that models the control directions
- `x0`: initial condition
- `x`: current state of system
- `xs`: system trajectory
"""
mutable struct ControlAffineSystem <: System
    n::Int
    m::Int
    f
    g
    x0
    x
    xs
end

"""
    ControlAffineSystem(n::Int, m::Int, f, g)
    ControlAffineSystem(n::Int, m::Int, f, g, x0)

Construct a nonlinear control affine system.
"""
ControlAffineSystem(n::Int, m::Int, f, g) = ControlAffineSystem(n, m, f, g, zeros(n), missing, missing)
ControlAffineSystem(n::Int, m::Int, f, g, x0) = ControlAffineSystem(n, m, f, g, x0, missing, missing)

"""
    initialize!(Σ::ControlAffineSystem)

Set current state to initial condition and initialize system trajectory.
"""
function initialize!(Σ::ControlAffineSystem)
    Σ.x = Σ.x0
    Σ.xs = [Σ.x]

    return Σ
end

"""
    step!(Σ::ControlAffineSystem, u, dt::Float64)

Integrate the dynamics of the system under control input u for dt seconds.
"""
function step!(Σ::ControlAffineSystem, u, dt::Float64)
    rhs(x, u, t) = Σ.f(x) + Σ.g(x)u
    sol = solve(ODEProblem(rhs, Σ.x, (0, dt), u))
    Σ.x = Σ.n == 1 ? sol[end] : sol[:, end]
    push!(Σ.xs, Σ.x)

    return Σ
end


text:
import Base: +, *, -, ==, getindex, setindex!
# import Base: +, *

abstract type AbstractSpline{N} end
struct Spline{N} <: AbstractSpline{N}
    tabpol::Vector{Polynomials.Polynomial{Rational{N}}}
    function Spline(tabpol::Vector{Polynomials.Polynomial{Rational{N}}}) where {N<:Signed}
        return new{N}(tabpol)
    end
end

function Base.getindex(sp::AbstractSpline{N}, index::Integer) where {N<:Signed}
    i = index + 1
    return if 1 <= i <= size(sp.tabpol, 1)
        sp.tabpol[i]
    else
        zero(Polynomials.Polynomial{Rational{N}})
    end
end
# function Base.getindex(sp::AbstractSpline{N}, index::AbstractRange) where{N<:Signed}
#     return Base.getindex.((sp::AbstractSpline{N},), index::AbstractRange)
# end
# function Base.setindex!(sp::AbstractSpline{N}, pol::Polynomials.Polynomial{Rational{N}}, index) where{N<:Signed}
#     sp.tabpol[index-1] = pol
# end
Base.size(sp::AbstractSpline, dim = 1) = size(sp.tabpol, 1)
function +(a::Spline{N}, b::Spline{N}) where {N<:Signed}
    sizenew = max(size(a.tabpol, 1), size(b.tabpol, 1))
    tabpolnew = zeros(Polynomials.Polynomial{Rational{N}}, sizenew)
    for i = 1:sizenew
        tabpolnew[i] += a[i-1] + b[i-1]
    end
    return (Spline(tabpolnew))
end
function -(a::Spline{N}, b::Spline{N}) where {N<:Signed}
    sizenew = max(size(a.tabpol, 1), size(b.tabpol, 1))
    tabpolnew = zeros(Polynomials.Polynomial{Rational{N}}, sizenew)
    for i = 1:sizenew
        tabpolnew[i] += a[i-1] - b[i-1]
    end
    return (Spline(tabpolnew))
end
function ==(a::Spline{N}, b::Spline{N}) where {N<:Signed}
    return a.tabpol == b.tabpol
end

function *(a::Spline{N}, pol::Polynomials.Polynomial{Rational{N}}) where {N<:Signed}
    tabpol2 = deepcopy(a.tabpol)
    for i = 1:size(tabpol2, 1)
        tabpol2[i] *= pol
    end
    return Spline(tabpol2)
end

function decal(a::Spline{N}, n) where {N<:Signed}
    return if n == 0
        a
    else
        tabpol = deepcopy(a.tabpol)
        poldec = Polynomials.Polynomial([-n // 1, one(Rational{N})])
        for i = 1:size(a.tabpol, 1)
            tabpol[i] = a.tabpol[i](poldec)
        end
        Spline(vcat(zeros(Polynomials.Polynomial{Rational{N}}, n), tabpol))
    end
end
w(p, j) = Polynomials.Polynomial([-j // p, 1 // p])
function _getbspline(n::N, j) where {N<:Signed}
    #   println("_getbspline($n=$n , j=$j ) N=$N")
    if n == zero(N)
        ret = decal(Spline(ones(Polynomials.Polynomial{Rational{N}}, 1)), j)
    else
        n1 = _getbspline(n - 1, j)
        n2 = decal(n1, 1)
        ret = n1 * w(n, j) + n2 * (1 - w(n, j + 1))
    end
    #    println("n=$n j=$j N=$N ret=$ret")
    return ret
end
function getbspline(n, j)
    return if (n > 8)
        _getbspline(big(n), j)
    else
        _getbspline(n, j)
    end
end
function (f::Spline{N})(x) where {N<:Signed}
    i = Int64((floor(x))) # it's different from ceil(x)
    if 0 <= i < size(f, 1)
        return f[i](x)
    else
        return zero(x)
    end
end
# struct SplineInt{N} <: AbstractSpline{N}
#     fact_order::N
#     tabpol::Vector{Polynomials.Polynomial{N}}
# end
# function SplineInt(order)
#     sp = getbspline(order, 0)
#     N = order <= 13 ? Int64 : BigInt
#     fact_order = factorial(big(order))
#     return SplineInt{N}(N(fact_order), map( x->Polynomial(N.(fact_order*coeffs(sp[x]))), 0:order))
# end
# function (f::SplineInt{N})(x::T) where{N<:Signed, T <: AbstractFloat}
#     i = Int64((floor(x))) # it's different from ceil(x)
#     if 0 <= i < size(f,1)
#         return f[i](x)/f.fact_order
#     else
#         return zero(x)/f.fact_order
#     end
# end
# function (f::SplineInt{N})(x::Union{Rational,Int}) where{N<:Signed}
#     i = Int64((floor(x))) # it's different from ceil(x)
#     if 0 <= i < size(f,1)
#         return f[i](x)//f.fact_order
#     else
#         return zero(x)//f.fact_order
#     end
# end


text:
module NVML

using ..APIUtils

using ..CUDA

using CEnum

using Memoize

using Libdl

@memoize function libnvml()
    if Sys.iswindows()
        # the NVSMI dir isn't added to PATH by the installer
        nvsmi = joinpath(ENV["ProgramFiles"], "NVIDIA Corporation", "NVSMI")
        if isdir(nvsmi)
            joinpath(nvsmi, "nvml.dll")
        else
            # let's just hope for the best
            "nvml"
        end
    else
        "libnvidia-ml.so.1"
    end
end
@memoize has_nvml() = Libdl.dlopen(libnvml(); throw_error=false) !== nothing

# core library
include("libnvml_common.jl")
include("error.jl")
include("libnvml.jl")
include("libnvml_deprecated.jl")

# wrappers
include("system.jl")
include("device.jl")

end