text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module GMP

export BigInt

import Base: *, +, -, /, <, <<, >>, >>>, <=, ==, >, >=, ^, (~), (&), (|), xor,
             binomial, cmp, convert, div, divrem, factorial, fld, gcd, gcdx, lcm, mod,
             ndigits, promote_rule, rem, show, isqrt, string, powermod,
             sum, trailing_zeros, trailing_ones, count_ones, base, tryparse_internal,
             bin, oct, dec, hex, isequal, invmod, prevpow2, nextpow2, ndigits0z, widen, signed, unsafe_trunc, trunc,
             iszero, big

if Clong == Int32
    const ClongMax = Union{Int8, Int16, Int32}
    const CulongMax = Union{UInt8, UInt16, UInt32}
else
    const ClongMax = Union{Int8, Int16, Int32, Int64}
    const CulongMax = Union{UInt8, UInt16, UInt32, UInt64}
end
const CdoubleMax = Union{Float16, Float32, Float64}

gmp_version() = VersionNumber(unsafe_string(unsafe_load(cglobal((:__gmp_version, :libgmp), Ptr{Cchar}))))
gmp_bits_per_limb() = Int(unsafe_load(cglobal((:__gmp_bits_per_limb, :libgmp), Cint)))

const GMP_VERSION = gmp_version()
const GMP_BITS_PER_LIMB = gmp_bits_per_limb()

# GMP's mp_limb_t is by default a typedef of `unsigned long`, but can also be configured to be either
# `unsigned int` or `unsigned long long int`. The correct unsigned type is here named Limb, and must
# be used whenever mp_limb_t is in the signature of ccall'ed GMP functions.
if GMP_BITS_PER_LIMB == 32
    const Limb = UInt32
elseif GMP_BITS_PER_LIMB == 64
    const Limb = UInt64
else
    error("GMP: cannot determine the type mp_limb_t (__gmp_bits_per_limb == $GMP_BITS_PER_LIMB)")
end

mutable struct BigInt <: Integer
    alloc::Cint
    size::Cint
    d::Ptr{Limb}
    function BigInt()
        b = new(zero(Cint), zero(Cint), C_NULL)
        ccall((:__gmpz_init,:libgmp), Void, (Ptr{BigInt},), &b)
        finalizer(b, cglobal((:__gmpz_clear, :libgmp)))
        return b
    end
end

function __init__()
    try
        if gmp_version().major != GMP_VERSION.major || gmp_bits_per_limb() != GMP_BITS_PER_LIMB
            error(string("The dynamically loaded GMP library (version $(gmp_version()) with __gmp_bits_per_limb == $(gmp_bits_per_limb()))\n",
                         "does not correspond to the compile time version (version $GMP_VERSION with __gmp_bits_per_limb == $GMP_BITS_PER_LIMB).\n",
                         "Please rebuild Julia."))
        end

        ccall((:__gmp_set_memory_functions, :libgmp), Void,
              (Ptr{Void},Ptr{Void},Ptr{Void}),
              cglobal(:jl_gc_counted_malloc),
              cglobal(:jl_gc_counted_realloc_with_old_size),
              cglobal(:jl_gc_counted_free))
    catch ex
        Base.showerror_nostdio(ex,
            "WARNING: Error during initialization of module GMP")
    end
end

widen(::Type{Int128})  = BigInt
widen(::Type{UInt128}) = BigInt
widen(::Type{BigInt})  = BigInt

signed(x::BigInt) = x

convert(::Type{BigInt}, x::BigInt) = x

function tryparse_internal(::Type{BigInt}, s::AbstractString, startpos::Int, endpos::Int, base_::Integer, raise::Bool)
    _n = Nullable{BigInt}()

    # don't make a copy in the common case where we are parsing a whole String
    bstr = startpos == start(s) && endpos == endof(s) ? String(s) : String(SubString(s,startpos,endpos))

    sgn, base, i = Base.parseint_preamble(true,Int(base_),bstr,start(bstr),endof(bstr))
    if !(2 <= base <= 62)
        raise && throw(ArgumentError("invalid base: base must be 2 ≤ base ≤ 62, got $base"))
        return _n
    end
    if i == 0
        raise && throw(ArgumentError("premature end of integer: $(repr(bstr))"))
        return _n
    end
    z = BigInt()
    if Base.containsnul(bstr)
        err = -1 # embedded NUL char (not handled correctly by GMP)
    else
        err = ccall((:__gmpz_set_str, :libgmp),
                    Int32, (Ptr{BigInt}, Ptr{UInt8}, Int32),
                    &z, pointer(bstr)+(i-start(bstr)), base)
    end
    if err != 0
        raise && throw(ArgumentError("invalid BigInt: $(repr(bstr))"))
        return _n
    end
    Nullable(sgn < 0 ? -z : z)
end

function convert(::Type{BigInt}, x::Union{Clong,Int32})
    z = BigInt()
    ccall((:__gmpz_set_si, :libgmp), Void, (Ptr{BigInt}, Clong), &z, x)
    return z
end
function convert(::Type{BigInt}, x::Union{Culong,UInt32})
    z = BigInt()
    ccall((:__gmpz_set_ui, :libgmp), Void, (Ptr{BigInt}, Culong), &z, x)
    return z
end

convert(::Type{BigInt}, x::Bool) = BigInt(UInt(x))

function unsafe_trunc(::Type{BigInt}, x::Union{Float32,Float64})
    z = BigInt()
    ccall((:__gmpz_set_d, :libgmp), Void, (Ptr{BigInt}, Cdouble), &z, x)
    return z
end

function convert(::Type{BigInt}, x::Union{Float32,Float64})
    isinteger(x) || throw(InexactError())
    unsafe_trunc(BigInt,x)
end

function trunc(::Type{BigInt}, x::Union{Float32,Float64})
    isfinite(x) || throw(InexactError())
    unsafe_trunc(BigInt,x)
end

convert(::Type{BigInt}, x::Float16) = BigInt(Float64(x))
convert(::Type{BigInt}, x::Float32) = BigInt(Float64(x))

function convert(::Type{BigInt}, x::Integer)
    if x < 0
        if typemin(Clong) <= x
            return BigInt(convert(Clong,x))
        end
        b = BigInt(0)
        shift = 0
        while x < -1
            b += BigInt(~UInt32(x&0xffffffff))<<shift
            x >>= 32
            shift += 32
        end
        return -b-1
    else
        if x <= typemax(Culong)
            return BigInt(convert(Culong,x))
        end
        b = BigInt(0)
        shift = 0
        while x > 0
            b += BigInt(UInt32(x&0xffffffff))<<shift
            x >>>= 32
            shift += 32
        end
        return b
    end
end

rem(x::BigInt, ::Type{Bool}) = ((x&1)!=0)
function rem{T<:Union{Unsigned,Signed}}(x::BigInt, ::Type{T})
    u = zero(T)
    for l = 1:min(abs(x.size), cld(sizeof(T),sizeof(Limb)))
        u += (unsafe_load(x.d,l)%T) << ((sizeof(Limb)<<3)*(l-1))
    end
    x.size < 0 ? -u : u
end

rem(x::Integer, ::Type{BigInt}) = convert(BigInt, x)

function convert(::Type{T}, x::BigInt) where T<:Unsigned
    if sizeof(T) < sizeof(Limb)
        convert(T, convert(Limb,x))
    else
        0 <= x.size <= cld(sizeof(T),sizeof(Limb)) || throw(InexactError())
        x % T
    end
end

function convert(::Type{T}, x::BigInt) where T<:Signed
    n = abs(x.size)
    if sizeof(T) < sizeof(Limb)
        SLimb = typeof(Signed(one(Limb)))
        convert(T, convert(SLimb, x))
    else
        0 <= n <= cld(sizeof(T),sizeof(Limb)) || throw(InexactError())
        y = x % T
        (x.size > 0) ⊻ (y > 0) && throw(InexactError()) # catch overflow
        y
    end
end

function (::Type{Float64})(n::BigInt, ::RoundingMode{:ToZero})
    ccall((:__gmpz_get_d, :libgmp), Float64, (Ptr{BigInt},), &n)
end

function (::Type{T})(n::BigInt, ::RoundingMode{:ToZero}) where T<:Union{Float16,Float32}
    T(Float64(n,RoundToZero),RoundToZero)
end

function (::Type{T})(n::BigInt, ::RoundingMode{:Down}) where T<:CdoubleMax
    x = T(n,RoundToZero)
    x > n ? prevfloat(x) : x
end
function (::Type{T})(n::BigInt, ::RoundingMode{:Up}) where T<:CdoubleMax
    x = T(n,RoundToZero)
    x < n ? nextfloat(x) : x
end

function (::Type{T})(n::BigInt, ::RoundingMode{:Nearest}) where T<:CdoubleMax
    x = T(n,RoundToZero)
    if maxintfloat(T) <= abs(x) < T(Inf)
        r = n-BigInt(x)
        h = eps(x)/2
        if iseven(reinterpret(Unsigned,x)) # check if last bit is odd/even
            if r < -h
                return prevfloat(x)
            elseif r > h
                return nextfloat(x)
            end
        else
            if r <= -h
                return prevfloat(x)
            elseif r >= h
                return nextfloat(x)
            end
        end
    end
    x
end

convert(::Type{Float64}, n::BigInt) = Float64(n,RoundNearest)
convert(::Type{Float32}, n::BigInt) = Float32(n,RoundNearest)
convert(::Type{Float16}, n::BigInt) = Float16(n,RoundNearest)

promote_rule(::Type{BigInt}, ::Type{<:Integer}) = BigInt

big(::Type{<:Integer})  = BigInt
big(::Type{<:Rational}) = Rational{BigInt}

# Binary ops
for (fJ, fC) in ((:+, :add), (:-,:sub), (:*, :mul),
                 (:fld, :fdiv_q), (:div, :tdiv_q), (:mod, :fdiv_r), (:rem, :tdiv_r),
                 (:gcd, :gcd), (:lcm, :lcm),
                 (:&, :and), (:|, :ior), (:xor, :xor))
    @eval begin
        function ($fJ)(x::BigInt, y::BigInt)
            z = BigInt()
            ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &x, &y)
            return z
        end
    end
end

/(x::BigInt, y::BigInt) = float(x)/float(y)

function invmod(x::BigInt, y::BigInt)
    z = zero(BigInt)
    ya = abs(y)
    if ya == 1
        return z
    end
    if (y==0 || ccall((:__gmpz_invert, :libgmp), Cint, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &x, &ya) == 0)
        throw(DomainError())
    end
    # GMP always returns a positive inverse; we instead want to
    # normalize such that div(z, y) == 0, i.e. we want a negative z
    # when y is negative.
    if y < 0
        z = z + y
    end
    # The postcondition is: mod(z * x, y) == mod(big(1), m) && div(z, y) == 0
    return z
end

# More efficient commutative operations
for (fJ, fC) in ((:+, :add), (:*, :mul), (:&, :and), (:|, :ior), (:xor, :xor))
    @eval begin
        function ($fJ)(a::BigInt, b::BigInt, c::BigInt)
            z = BigInt()
            ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &a, &b)
            ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &z, &c)
            return z
        end
        function ($fJ)(a::BigInt, b::BigInt, c::BigInt, d::BigInt)
            z = BigInt()
            ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &a, &b)
            ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &z, &c)
            ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &z, &d)
            return z
        end
        function ($fJ)(a::BigInt, b::BigInt, c::BigInt, d::BigInt, e::BigInt)
            z = BigInt()
            ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &a, &b)
            ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &z, &c)
            ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &z, &d)
            ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &z, &e)
            return z
        end
    end
end

# Basic arithmetic without promotion
function +(x::BigInt, c::CulongMax)
    z = BigInt()
    ccall((:__gmpz_add_ui, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Culong), &z, &x, c)
    return z
end
+(c::CulongMax, x::BigInt) = x + c

function -(x::BigInt, c::CulongMax)
    z = BigInt()
    ccall((:__gmpz_sub_ui, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Culong), &z, &x, c)
    return z
end
function -(c::CulongMax, x::BigInt)
    z = BigInt()
    ccall((:__gmpz_ui_sub, :libgmp), Void, (Ptr{BigInt}, Culong, Ptr{BigInt}), &z, c, &x)
    return z
end
+(x::BigInt, c::ClongMax) = c < 0 ? -(x, -(c % Culong)) : x + convert(Culong, c)
+(c::ClongMax, x::BigInt) = c < 0 ? -(x, -(c % Culong)) : x + convert(Culong, c)
-(x::BigInt, c::ClongMax) = c < 0 ? +(x, -(c % Culong)) : -(x, convert(Culong, c))
-(c::ClongMax, x::BigInt) = c < 0 ? -(x + -(c % Culong)) : -(convert(Culong, c), x)

function *(x::BigInt, c::CulongMax)
    z = BigInt()
    ccall((:__gmpz_mul_ui, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Culong), &z, &x, c)
    return z
end
*(c::CulongMax, x::BigInt) = x * c
function *(x::BigInt, c::ClongMax)
    z = BigInt()
    ccall((:__gmpz_mul_si, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Clong), &z, &x, c)
    return z
end
*(c::ClongMax, x::BigInt) = x * c

/(x::BigInt, y::Union{ClongMax,CulongMax}) = float(x)/y
/(x::Union{ClongMax,CulongMax}, y::BigInt) = x/float(y)

# unary ops
for (fJ, fC) in ((:-, :neg), (:~, :com))
    @eval begin
        function ($fJ)(x::BigInt)
            z = BigInt()
            ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}), &z, &x)
            return z
        end
    end
end

function <<(x::BigInt, c::UInt)
    c == 0 && return x
    z = BigInt()
    ccall((:__gmpz_mul_2exp, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Culong), &z, &x, c)
    return z
end

function >>(x::BigInt, c::UInt)
    c == 0 && return x
    z = BigInt()
    ccall((:__gmpz_fdiv_q_2exp, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Culong), &z, &x, c)
    return z
end

>>>(x::BigInt, c::UInt) = x >> c

trailing_zeros(x::BigInt) = Int(ccall((:__gmpz_scan1, :libgmp), Culong, (Ptr{BigInt}, Culong), &x, 0))
trailing_ones(x::BigInt) = Int(ccall((:__gmpz_scan0, :libgmp), Culong, (Ptr{BigInt}, Culong), &x, 0))

count_ones(x::BigInt) = Int(ccall((:__gmpz_popcount, :libgmp), Culong, (Ptr{BigInt},), &x))

function divrem(x::BigInt, y::BigInt)
    z1 = BigInt()
    z2 = BigInt()
    ccall((:__gmpz_tdiv_qr, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z1, &z2, &x, &y)
    z1, z2
end

function cmp(x::BigInt, y::BigInt)
    ccall((:__gmpz_cmp, :libgmp), Int32, (Ptr{BigInt}, Ptr{BigInt}), &x, &y)
end
function cmp(x::BigInt, y::ClongMax)
    ccall((:__gmpz_cmp_si, :libgmp), Int32, (Ptr{BigInt}, Clong), &x, y)
end
function cmp(x::BigInt, y::CulongMax)
    ccall((:__gmpz_cmp_ui, :libgmp), Int32, (Ptr{BigInt}, Culong), &x, y)
end
cmp(x::BigInt, y::Integer) = cmp(x,big(y))
cmp(x::Integer, y::BigInt) = -cmp(y,x)

function cmp(x::BigInt, y::CdoubleMax)
    isnan(y) && throw(DomainError())
    ccall((:__gmpz_cmp_d, :libgmp), Int32, (Ptr{BigInt}, Cdouble), &x, y)
end
cmp(x::CdoubleMax, y::BigInt) = -cmp(y,x)

function isqrt(x::BigInt)
    z = BigInt()
    ccall((:__gmpz_sqrt, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}), &z, &x)
    return z
end

function ^(x::BigInt, y::Culong)
    z = BigInt()
    ccall((:__gmpz_pow_ui, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Culong), &z, &x, y)
    return z
end

function bigint_pow(x::BigInt, y::Integer)
    if y<0; throw(DomainError()); end
    if x== 1; return x; end
    if x==-1; return isodd(y) ? x : -x; end
    if y>typemax(Culong)
       x==0 && return x

       #At this point, x is not 1, 0 or -1 and it is not possible to use
       #gmpz_pow_ui to compute the answer. Note that the magnitude of the
       #answer is:
       #- at least 2^(2^32-1) ≈ 10^(1.3e9) (if Culong === UInt32).
       #- at least 2^(2^64-1) ≈ 10^(5.5e18) (if Culong === UInt64).
       #
       #Assume that the answer will definitely overflow.

       throw(OverflowError())
    end
    return x^convert(Culong, y)
end

^(x::BigInt , y::BigInt ) = bigint_pow(x, y)
^(x::BigInt , y::Bool   ) = y ? x : one(x)
^(x::BigInt , y::Integer) = bigint_pow(x, y)
^(x::Integer, y::BigInt ) = bigint_pow(BigInt(x), y)
^(x::Bool   , y::BigInt ) = Base.power_by_squaring(x, y)

function powermod(x::BigInt, p::BigInt, m::BigInt)
    r = BigInt()
    ccall((:__gmpz_powm, :libgmp), Void,
          (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}),
          &r, &x, &p, &m)
    return m < 0 && r > 0 ? r + m : r # choose sign conistent with mod(x^p, m)
end

powermod(x::Integer, p::Integer, m::BigInt) = powermod(big(x), big(p), m)

function gcdx(a::BigInt, b::BigInt)
    if b == 0 # shortcut this to ensure consistent results with gcdx(a,b)
        return a < 0 ? (-a,-one(BigInt),zero(BigInt)) : (a,one(BigInt),zero(BigInt))
    end
    g = BigInt()
    s = BigInt()
    t = BigInt()
    ccall((:__gmpz_gcdext, :libgmp), Void,
        (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}),
        &g, &s, &t, &a, &b)
    if t == 0
        # work around a difference in some versions of GMP
        if a == b
            return g, t, s
        elseif abs(a)==abs(b)
            return g, t, -s
        end
    end
    g, s, t
end

function sum(arr::AbstractArray{BigInt})
    n = BigInt(0)
    for i in arr
        ccall((:__gmpz_add, :libgmp), Void,
            (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}),
            &n, &n, &i)
    end
    return n
end

function factorial(x::BigInt)
    x.size < 0 && return BigInt(0)
    z = BigInt()
    ccall((:__gmpz_fac_ui, :libgmp), Void, (Ptr{BigInt}, Culong), &z, x)
    return z
end

function binomial(n::BigInt, k::UInt)
    z = BigInt()
    ccall((:__gmpz_bin_ui, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Culong), &z, &n, k)
    return z
end
binomial(n::BigInt, k::Integer) = k < 0 ? BigInt(0) : binomial(n, UInt(k))

==(x::BigInt, y::BigInt) = cmp(x,y) == 0
==(x::BigInt, i::Integer) = cmp(x,i) == 0
==(i::Integer, x::BigInt) = cmp(x,i) == 0
==(x::BigInt, f::CdoubleMax) = isnan(f) ? false : cmp(x,f) == 0
==(f::CdoubleMax, x::BigInt) = isnan(f) ? false : cmp(x,f) == 0
iszero(x::BigInt) = x == Clong(0)

<=(x::BigInt, y::BigInt) = cmp(x,y) <= 0
<=(x::BigInt, i::Integer) = cmp(x,i) <= 0
<=(i::Integer, x::BigInt) = cmp(x,i) >= 0
<=(x::BigInt, f::CdoubleMax) = isnan(f) ? false : cmp(x,f) <= 0
<=(f::CdoubleMax, x::BigInt) = isnan(f) ? false : cmp(x,f) >= 0

<(x::BigInt, y::BigInt) = cmp(x,y) < 0
<(x::BigInt, i::Integer) = cmp(x,i) < 0
<(i::Integer, x::BigInt) = cmp(x,i) > 0
<(x::BigInt, f::CdoubleMax) = isnan(f) ? false : cmp(x,f) < 0
<(f::CdoubleMax, x::BigInt) = isnan(f) ? false : cmp(x,f) > 0

string(x::BigInt) = dec(x)
show(io::IO, x::BigInt) = print(io, string(x))

bin(n::BigInt) = base( 2, n)
oct(n::BigInt) = base( 8, n)
dec(n::BigInt) = base(10, n)
hex(n::BigInt) = base(16, n)

bin(n::BigInt, pad::Int) = base( 2, n, pad)
oct(n::BigInt, pad::Int) = base( 8, n, pad)
dec(n::BigInt, pad::Int) = base(10, n, pad)
hex(n::BigInt, pad::Int) = base(16, n, pad)

function base(b::Integer, n::BigInt)
    2 <= b <= 62 || throw(ArgumentError("base must be 2 ≤ base ≤ 62, got $b"))
    nd = ndigits(n, b)
    str = Base._string_n(n < 0 ? nd+1 : nd)
    ccall((:__gmpz_get_str,:libgmp), Ptr{UInt8}, (Ptr{UInt8}, Cint, Ptr{BigInt}), str, b, &n)
    return str
end

function base(b::Integer, n::BigInt, pad::Integer)
    s = base(b, n)
    buf = IOBuffer()
    if n < 0
        s = s[2:end]
        write(buf, '-')
    end
    for i in 1:pad-sizeof(s) # `s` is known to be ASCII, and `length` is slower
        write(buf, '0')
    end
    write(buf, s)
    String(buf)
end

function ndigits0z(x::BigInt, b::Integer=10)
    b < 2 && throw(DomainError())
    if ispow2(b) && 2 <= b <= 62 # GMP assumes b is in this range
        Int(ccall((:__gmpz_sizeinbase,:libgmp), Csize_t, (Ptr{BigInt}, Cint), &x, b))
    else
        # non-base 2 mpz_sizeinbase might return an answer 1 too big
        # use property that log(b, x) < ndigits(x, b) <= log(b, x) + 1
        n = Int(ccall((:__gmpz_sizeinbase,:libgmp), Csize_t, (Ptr{BigInt}, Cint), &x, 2))
        lb = log2(b) # assumed accurate to <1ulp (true for openlibm)
        q,r = divrem(n,lb)
        iq = Int(q)
        maxerr = q*eps(lb) # maximum error in remainder
        if r-1.0 < maxerr
            abs(x) >= big(b)^iq ? iq+1 : iq
        elseif lb-r < maxerr
            abs(x) >= big(b)^(iq+1) ? iq+2 : iq+1
        else
            iq+1
        end
    end
end
ndigits(x::BigInt, b::Integer=10) = x.size == 0 ? 1 : ndigits0z(x,b)

prevpow2(x::BigInt) = x.size < 0 ? -prevpow2(-x) : (x <= 2 ? x : one(BigInt) << (ndigits(x, 2)-1))
nextpow2(x::BigInt) = x.size < 0 ? -nextpow2(-x) : (x <= 2 ? x : one(BigInt) << ndigits(x-1, 2))

Base.checked_abs(x::BigInt) = abs(x)
Base.checked_neg(x::BigInt) = -x
Base.checked_add(a::BigInt, b::BigInt) = a + b
Base.checked_sub(a::BigInt, b::BigInt) = a - b
Base.checked_mul(a::BigInt, b::BigInt) = a * b
Base.checked_div(a::BigInt, b::BigInt) = div(a, b)
Base.checked_rem(a::BigInt, b::BigInt) = rem(a, b)
Base.checked_fld(a::BigInt, b::BigInt) = fld(a, b)
Base.checked_mod(a::BigInt, b::BigInt) = mod(a, b)
Base.checked_cld(a::BigInt, b::BigInt) = cld(a, b)
Base.add_with_overflow(a::BigInt, b::BigInt) = a + b, false
Base.sub_with_overflow(a::BigInt, b::BigInt) = a - b, false
Base.mul_with_overflow(a::BigInt, b::BigInt) = a * b, false

function Base.deepcopy_internal(x::BigInt, stackdict::ObjectIdDict)
    if haskey(stackdict, x)
        return stackdict[x]
    end
    y = BigInt()
    ccall((:__gmpz_set,:libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}), &y, &x)
    stackdict[x] = y
    return y
end

end # module

max_stars_repo_name: odow/julia

max_issues_repo_name: odow/julia

max_forks_repo_name: odow/julia


text:
################################################################################
#
#  String I/O
#
################################################################################

function Base.show(io::IO, L::QuadLat)
  println(io, "Quadratic lattice of rank $(rank(L)) and degree $(degree(L))")
  println(io, "over")
  print(IOContext(io, :compact => true), base_ring(L))
end

################################################################################
#
#  Construction
#
################################################################################

@doc Markdown.doc"""
    quadratic_lattice(K::NumField, B::PMat; gram_ambient_space = F) -> QuadLat

Given a number field $K$ and a pseudo-matrix $B$, returns the quadratic lattice
spanned by the pseudo-matrix $B$ inside the quadratic space over $K$ with gram
matrix $F$.

If $F$ is not supplied, the gram matrix of the ambient space will be the
identity matrix.
"""
quadratic_lattice(::NumField, ::PMat; gram_ambient_space = nothing)

# TODO: At the moment I assume that B is a pseudo-hnf (probably)
function quadratic_lattice(K::NumField, B::PMat; gram_ambient_space = nothing, gram = nothing)
  if gram_ambient_space === nothing && gram === nothing
    return QuadLat(K, identity_matrix(K, nrows(B)), pseudo_matrix(B))
  end
  if gram_ambient_space !== nothing && gram === nothing
    return QuadLat(K, gram_ambient_space, B)
  end
  if gram_ambient_space === nothing && gram !== nothing
    z = QuadLat{typeof(K), typeof(gram), typeof(B)}()
    z.pmat = B
    z.gram = gram
    z.base_algebra = K
    return z
  end
end

@doc Markdown.doc"""
    quadratic_lattice(K::NumField, B::MatElem; gram_ambient_space = F) -> QuadLat

Given a number field $K$ and a matrix $B$, returns the quadratic lattice
spanned by the rows of $B$ inside the quadratic space over $K$ with gram matrix
$F$.

If $F$ is not supplied, the gram matrix of the ambient space will be the
identity matrix.
"""
quadratic_lattice(::NumField, ::MatElem; gram_ambient_space = nothing)

function quadratic_lattice(K::NumField, B::MatElem; gram_ambient_space = nothing, gram = nothing)
  if gram_ambient_space === nothing && gram === nothing
    return QuadLat(K, identity_matrix(K, nrows(B)), pseudo_matrix(B))
  end
  if gram_ambient_space !== nothing && gram === nothing
    return QuadLat(K, gram_ambient_space, pseudo_matrix(B))
  end
  if gram_ambient_space === nothing && gram !== nothing
    P = pseudo_matrix(B)
    z = QuadLat{typeof(K), typeof(B), typeof(P)}()
    z.pmat = P
    z.gram = gram
    z.base_algebra = K
    return z
  end
end

function quadratic_lattice(K::NumField; generators::Vector = Vector{elem_type(K)}[], gram_ambient_space)
  if length(generators) == 0
    pm = pseudo_matrix(identity_matrix(K, nrows(gram_ambient_space)))
  else
    z = zero_matrix(K, length(generators), ncols(gram_ambient_space))
    for i in 1:length(generators)
      for j in 1:ncols(gram_ambient_space)
        z[i, j] = generators[i][j]
      end
    end
    pm = pseudo_hnf(pseudo_matrix(z), :lowerleft)
    i = 1
    while iszero_row(pm.matrix, i)
      i += 1
    end
    pm = sub(pm, i:nrows(pm), 1:ncols(pm))
  end
  L = quadratic_lattice(K, pm, gram_ambient_space = gram_ambient_space)
  L.generators = Vector{elem_type(K)}[map(K, v) for v in generators]
  return L
end

@doc Markdown.doc"""
    lattice(V::QuadSpace, B::PMat) -> QuadLat

Given a quadratic space $V$ and a pseudo-matrix $B$, returns the quadratic lattice
spanned by the pseudo-matrix $B$ inside $V$.
"""
function lattice(V::QuadSpace, B::PMat)
  K = base_ring(V)
  z = QuadLat{typeof(K), typeof(gram_matrix(V)), typeof(B)}()
  z.pmat = B
  z.gram = gram_matrix(V, matrix(B))
  z.base_algebra = K
  z.space = V
  return z
end

@doc Markdown.doc"""
    lattice(V::QuadSpace, B::MatElem) -> QuadLat

Given a quadratic space $V$ and a matrix $B$, returns the quadratic lattice
spanned by the rows of $B$ inside $V$.
"""
function lattice(V::QuadSpace, B::MatElem)
  K = base_ring(V)
  pmat = pseudo_matrix(B)
  z = QuadLat{typeof(K), typeof(gram_matrix(V)), typeof(pmat)}()
  z.pmat = pmat
  z.gram = gram_matrix(V, B)
  z.base_algebra = K
  z.space = V
  return z
end

@doc Markdown.doc"""
    lattice(V::QuadSpace) -> QuadLat

Given a quadratic space $V$, returns the quadratic lattice with trivial basis
matrix.
"""
lattice(V::QuadSpace) = lattice(V, identity_matrix(base_ring(V), rank(V)))

################################################################################
#
#  Rational span
#
################################################################################

# docstring in ../Lattices.jl

function rational_span(L::QuadLat)
  if isdefined(L, :rational_span)
    return L.rational_span
  else
    G = gram_matrix_of_rational_span(L)
    V = quadratic_space(base_field(L), G)
    L.rational_span = V
    return V
  end
end

################################################################################
#
#  Involution
#
################################################################################

# docstring in ../Lattices.jl

involution(L::QuadLat) = identity

################################################################################
#
#  Hasse invariant
#
################################################################################

function hasse_invariant(L::QuadLat, p)
  return hasse_invariant(rational_span(L), p)
end

################################################################################
#
#  Witt invariant
#
################################################################################

function witt_invariant(L::QuadLat, p::Union{NfAbsOrdIdl, InfPlc})
  return witt_invariant(rational_span(L), p)
end

################################################################################
#
#  New lattice in same ambient space
#
################################################################################

function lattice_in_same_ambient_space(L::QuadLat, m::PMat)
  return quadratic_lattice(base_field(L), m,
                           gram_ambient_space = gram_matrix(ambient_space(L)))
end

################################################################################
#
#  Norm
#
################################################################################

function norm(L::QuadLat)
  if isdefined(L, :norm)
    return L.norm::fractional_ideal_type(base_ring(L))
  end
  G = gram_matrix_of_rational_span(L)
  C = coefficient_ideals(L)
  K = nf(order(C[1]))
  n = sum(G[i, i] * C[i]^2 for i in 1:length(C)) + K(2) * scale(L)
  L.norm = n
  return n
end

################################################################################
#
#  Scale
#
################################################################################

function scale(L::QuadLat)
  if isdefined(L, :scale)
    return L.scale::fractional_ideal_type(base_ring(L))
  end
  G = gram_matrix_of_rational_span(L)
  C = coefficient_ideals(L)
  to_sum = [ G[i, j] * C[i] * involution(L)(C[j]) for j in 1:length(C) for i in 1:j]
  s = sum(to_sum)
  L.scale = s
  return s
end

################################################################################
#
#  Rescale
#
################################################################################

function rescale(L::QuadLat, a)
  if isone(a)
    return L
  end
  K = fixed_field(L)
  b = K(a)
  gramamb = gram_matrix(ambient_space(L))
  return quadratic_lattice(base_field(L), pseudo_matrix(L),
                           gram_ambient_space = b * gramamb)
end

################################################################################
#
#  Bad primes
#
################################################################################

@doc Markdown.doc"""
    bad_primes(L::AbsLat; even = false) -> Vector{NfOrdIdl}

Returns the prime ideals dividing the scale and volume of $L$. If `even == true`
also the prime ideals dividing $2$ are included.
"""

function bad_primes(L::QuadLat; even::Bool = false)
  f = factor(scale(L))
  ff = factor(volume(L))
  for (p, e) in ff
    f[p] = 0
  end
  if even
    for p in prime_decomposition(base_ring(L), 2)
      f[p[1]] = 0
    end
  end
  return collect(keys(f))
end

################################################################################
#
#  Dual lattice
#
################################################################################

function dual(L::QuadLat)
  G, B = _gram_schmidt(gram_matrix_of_rational_span(L), involution(L))
  C = coefficient_ideals(L)
  Gi = inv(G)
  new_bmat = Gi * B
  new_coeff = eltype(C)[inv(c) for c in C]
  pm = pseudo_matrix(new_bmat, new_coeff)
  return quadratic_lattice(base_field(L), pm,
                           gram_ambient_space = gram_matrix(ambient_space(L)))
end

################################################################################
#
#  Jordan decomposition
#
################################################################################

function jordan_decomposition(L::QuadLat, p)
  F = gram_matrix(ambient_space(L))
  even = isdyadic(p)
  if even
    pi = elem_in_nf(uniformizer(p))
  else
    pi = zero(nf(order(p)))
  end

  oldval = PosInf()
  blocks = Int[]
  exponents = Int[]
  S = local_basis_matrix(L, p)
  n = ncols(S)
  k = 1
  while k <= n
    G = S * F * transpose(S)
    X = Union{Int, PosInf}[ iszero(G[i, i]) ? inf : valuation(G[i, i], p) for i in k:n]
    m, ii = findmin(X)
    ii = ii + (k - 1)
    pair = (ii, ii)

    for i in k:n
      for j in (i + 1):n
        tmp = iszero(G[i, j]) ? inf : valuation(G[i, j], p)
        if tmp < m
          m = tmp
          pair = (i, j)
        end
      end
    end

    if m != oldval
      push!(blocks, k)
      oldval = m
      push!(exponents, m)
    end

    if even && pair[1] != pair[2]
      swap_rows!(S, pair[1], k)
      swap_rows!(S, pair[2], k +1)

      T12 = (sub(S, k:k, 1:ncols(S)) * F * transpose(sub(S, k+1:k+1, 1:ncols(S))))[1, 1]
      for l in 1:ncols(S)
        S[k, l] = S[k, l] * pi^valuation(T12, p)//T12
      end

      T = (i, j) -> (sub(S, i:i, 1:ncols(S)) * F * transpose(sub(S, j:j, 1:ncols(S))))[1, 1]
      T11 = T(k, k)
      T22 = T(k + 1, k + 1)
      T12 = T(k, k + 1)
      d = T11*T22 - T12^2
      for l in (k + 2):n
        tl = T12 * T(k + 1, l) - T22 * T(k, l)
        ul = T12 * T(k, l) - T11 * T(k + 1, l)
        for u in 1:ncols(S)
          S[l, u] = S[l, u] + (tl//d) * S[k, u] + (ul//d) * S[k + 1, u]
        end
      end
      k = k + 2
    else
      if pair[1] == pair[2]
        swap_rows!(S, pair[1], k)
      else
        for u in 1:ncols(S)
          S[pair[1], u] = S[pair[1], u] + S[pair[2], u]
        end
        swap_rows!(S, pair[1], k)
      end
      nrm = (sub(S, k:k, 1:ncols(S)) * F * transpose(sub(S, k:k, 1:ncols(S))))[1, 1]
      XX = sub(S, k:k, 1:ncols(S)) * F * transpose(S)
      for l in k+1:n
        for u in 1:ncols(S)
          S[l, u] = S[l, u] - XX[1,l]//nrm * S[k, u]
        end
      end
      k = k + 1
    end
  end

  push!(blocks, n+1)
  matrices = typeof(F)[ sub(S, blocks[i]:(blocks[i+1] - 1), 1:ncols(S)) for i in 1:(length(blocks) - 1)]
  return matrices, typeof(F)[ m * F * transpose(m) for m in matrices], exponents
end

################################################################################
#
#  Maximal integral lattices
#
################################################################################

function guess_max_det(L::QuadLat, p)
  m = rank(L)
  R = base_ring(L)
  n = div(m, 2)
  d = det(gram_matrix_of_rational_span(L))
  e = 2 * valuation(base_ring(L)(2), p)
  if isodd(m)
    v = mod(valuation(d, p), 2)
    v = witt_invariant(L, p) == 1 ? v - e * n : 2 - v - e * n
  else
    if isodd(div(m * (m - 1), 2))
      d = -d
    end
    qd = quadratic_defect(d, p)
    if qd isa PosInf
      v = witt_invariant(L, p) == 1 ? -e * n : 2 - e * n
    else
      vd = valuation(d, p)
      # I cannot use div(qd, 2), because qd might be negative and I need to round
      # toward 0, e.g., I need div(-3, 2) to be -2 and not -1.
      v = vd - 2 * Int(fdiv(fmpz(qd), 2)) + e * (1 - n)
      if iseven(vd) && qd == vd + e && witt_invariant(L, p) == -1
        v = -e*n + 2
      end
    end
  end
  return v
end

@doc Markdown.doc"""
    ismaximal_integral(L::QuadLat, p) -> Bool, QuadLat

Checks whether L is maximal integral at $p$. If not, the second return value is
a minimal integral overlattice at p.
"""
function ismaximal_integral(L::QuadLat, p)
  @req order(p) == base_ring(L) "blabla do not match"
  #if iszero(L)
  #  return true, L
  #end

  if valuation(norm(L), p) < 0
    # this is a weird case? Magma does not return a second argument
    return false, L
  end

  # don't know what this does
  if guess_max_det(L, p) == valuation(volume(L), p)
    return true, L
  end

  R = base_ring(L)
  K = nf(R)

  k, h = ResidueField(R, p)
  hext = extend(h, K)

  BM = local_basis_matrix(L, p, type = :submodule)

  G = 2 * gram_matrix(ambient_space(L), BM)

  Gmodp = map(hext, G)

  r, V = left_kernel(Gmodp)
  @assert r > 0
  local v::dense_matrix_type(K)
  if !isdyadic(p)
    T = map(y -> hext\y, V)
    H = inv(elem_in_nf(uniformizer(p))) * T * G * transpose(T)
    Hmod = map_entries(hext, H)
    ok, __v = _isisotropic_with_vector_finite(Hmod)
    @assert ok
    _v = matrix(k, 1, length(__v), __v)
    e = map_entries(x -> hext\x, _v * V)
    sp = (e * G * e')[1, 1]
    valv = iszero(sp) ? inf : valuation(sp, p)
    @assert valv >= 2
    v = e
  else
    val2 = valuation(R(2), p)
    PP = enumerate_lines(k, nrows(V))
    for x in PP
      @assert !iszero(x)
      xV = matrix(k, 1, length(x), x) * V
      e = elem_type(K)[ hext\(xV[1, i]) for i in 1:ncols(xV) ]
      v = matrix(K, 1, length(e), e)
      _z = (v * G * transpose(v))[1, 1]
      # Test if valv >= val2 + 2
      if iszero(_z)
        break
      else
        valv = valuation(_z, p)
        @assert valv >= 1
        if valv >= val2 + 2
          break
        end
      end
    end
  end
  pia = anti_uniformizer(p)
  _LL = _sum_modules(pseudo_matrix(L), pseudo_matrix(pia * v * BM))
  LL = lattice(ambient_space(L), _LL)
  @assert volume(L) ==  volume(LL) * p^2 && valuation(norm(LL), p) >= 0
  return false, LL
end

@doc Markdown.doc"""
    ismaximal_integral(L::QuadLat) -> Bool, QuadLat

Checks whether $L$ is maximal integral. If not, the second return valiue is a
minimal integral overlattice.
"""
function ismaximal_integral(L::QuadLat)
  #if iszero(L)
  #  return true, L
  #end

  if !isintegral(norm(L))
    # is L a minimal integral over-lattice? I don't think so
    return false, L
  end

  for p in bad_primes(L, even = true)
    ok, LL = ismaximal_integral(L, p)
    if !ok
      return false, LL
    end
  end
  return true, L
end

function maximal_integral_lattice(L::QuadLat, p)
  @req base_ring(L) == order(p) "Second argument must be an ideal of the base ring of L"
  @req valuation(norm(L), p) >= 0 "The normal of the lattice must be locally integral"

  ok, LL = ismaximal_integral(L, p)
  while !ok
    L = LL
    ok, LL = ismaximal_integral(L, p)
  end
  return L
end

@doc Markdown.doc"""
    ismaximal_integral(L::QuadLat, p) -> Bool, QuadLat

Checks whether $L$ is maximal at $p$, that is, whether no overlattice has the
same norm. If not, the second return value is a proper overlattice with the
same norm.
"""
function ismaximal(L::QuadLat, p)
  @req order(p) == base_ring(L) "Asdsads"
  #if iszero(L)
  #  return true, L
  #end
  v = valuation(norm(L), p)
  x = uniformizer(p)^(-v)
  ok, LL = ismaximal_integral(rescale(L, x), p)
  if ok
    return true, L
  else
    return false, rescale(LL, inv(elem_in_nf(x)))
  end
end

@doc Markdown.doc"""
    maximal_integral_lattice(V::QuadSpace) -> QuadLat

Returns a lattice $L$ of $V$ such that the norm of $L$ is integral and $L$ is
maximal with respect to this property.
"""
function maximal_integral_lattice(V::QuadSpace)
  K = base_ring(V)
  L = lattice(V, identity_matrix(K, rank(V)))
  n = norm(L)
  R = order(n)
  if !isone(norm(n))
    fa = factor(n)
    d = prod(typeof(n)[fractional_ideal(R, p)^(fld(e, 2)) for (p, e) in fa])
    # fld = fdiv = floored division
    L = lattice(V, _module_scale_ideal(inv(d), pseudo_matrix(L)))
    n = norm(L)
    @assert isintegral(n)
  end

  return maximal_integral_lattice(L)
end

@doc Markdown.doc"""
    maximal_integral_lattice(L::QuadLat) -> QuadLat

Returns a maximal integral lattice containing $L$.
"""
function maximal_integral_lattice(L::QuadLat)
  @req isintegral(norm(L)) "Lattice must be integral"
  for p in bad_primes(L, even = true)
    L = maximal_integral_lattice(L, p)
  end
  return L
end

ext: jl


text:
module OptiMimi

using NLopt
using ForwardDiff, DiffResults
using MathProgBase

import Mimi: Model

export problem, solution, unaryobjective, objevals, setparameters, nameindexes, sensitivity, uncertainproblem, setsolution

include("registerdiff.jl")
include("matrixconstraints.jl")
include("linproghouse.jl")
include("makeroom.jl")
include("metamimi.jl")
include("uncertainty.jl")

allverbose = false
objevals = 0

mutable struct OptimizationProblem
    model::Model
    components::Vector{Symbol}
    names::Vector{Symbol}
    opt::Opt
    constraints::Vector{Function}
end

mutable struct BlackBoxOptimizationProblem{T}
    algorithm::Symbol
    model::Model
    components::Vector{Symbol}
    names::Vector{Symbol}
    objective::Function
    lowers::Vector{T}
    uppers::Vector{T}
end

mutable struct LinprogOptimizationProblem{T}
    model::Model
    components::Vector{Symbol}
    names::Vector{Symbol}
    objective::Function
    objectiveconstraints::Vector{Function}
    matrixconstraints::Vector{MatrixConstraintSet}
    exlowers::Vector{T}
    exuppers::Vector{T}
end

BlackBoxAlgorithms = [:separable_nes, :xnes, :dxnes]

"""Returns (ii, len, isscalar) with the index of each symbol and its length."""
function nameindexes(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, chnl)
    for ii in 1:length(components)
        dims = getdims(model, components[ii], names[ii])
        if length(dims) == 0
            put!(chnl, (ii, 1, true)) # It's a scalar
        else
            put!(chnl, (ii, prod(dims), false))
        end
    end
end

"""
Initialize a model with the results from an optimization.
"""
function setparameters(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, xx::Vector)
    startindex = 1
    for (ii, len, isscalar) in Channel((chnl) -> nameindexes(model, components, names, chnl))
        if isscalar
            set_param!(model, components[ii], names[ii], xx[startindex])
        else
            shape = getdims(model, components[ii], names[ii])
            reshaped = reshape(collect(model.md.number_type, xx[startindex:(startindex+len - 1)]), tuple(shape...))
            set_param!(model, components[ii], names[ii], reshaped)
        end
        startindex += len
    end
end

function sensitivity(model::Model, component::Symbol, parameter::Symbol, objective::Function, points::Vector{Float64})
    results = []
    for point in points
        setparameters(model, [component], [parameter], [point])
        run(model)
        push!(results, objective(model))
    end

    results
end

"""Generate the form of objective function used by the optimization, taking parameters rather than a model."""
function unaryobjective(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, objective::Function)
    function my_objective(xx::Vector)
        if allverbose
            println(xx)
        end

        global objevals
        objevals += 1

        setparameters(model, components, names, xx)
        run(model)
        objective(model)
    end

    my_objective
end

"""Create an NLopt-style objective function which does not use its grad argument."""
function gradfreeobjective(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, objective::Function)
    myunaryobjective = unaryobjective(model, components, names, objective)
    function myobjective(xx::Vector, grad::Vector)
        myunaryobjective(xx)
    end

    myobjective
end

"""Create an NLopt-style objective function which computes an autodiff gradient."""
function autodiffobjective(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, objective::Function)
    myunaryobjective = unaryobjective(model, components, names, objective)
    function myobjective(xx::Vector, grad::Vector)
        out = DiffResults.GradientResult(xx)
        ForwardDiff.gradient!(out, myunaryobjective, xx)
        if any(isnan.(DiffResults.gradient(out)))
            error("objective gradient is NaN")
        end
        grad[:] = DiffResults.gradient(out)
        DiffResults.value(out)
    end

    myobjective
end

"""Create a 0 point."""
function make0(model::Model, components::Vector{Symbol}, names::Vector{Symbol})
    initial = Float64[]
    for (ii, len, isscalar) in Channel((chnl) -> nameindexes(model, components, names, chnl))
        append!(initial, [0. for jj in 1:len])
    end

    initial
end

"""Expand parameter constraints to full vectors for every numerical parameter."""
function expandlimits(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, lowers::Vector{T}, uppers::Vector{T}) where T <: Real
    my_lowers = T[]
    my_uppers = T[]

    ## Replace with eachname
    totalvars = 0
    for (ii, len, isscalar) in Channel((chnl) -> nameindexes(model, components, names, chnl))
        append!(my_lowers, [lowers[ii] for jj in 1:len])
        append!(my_uppers, [uppers[ii] for jj in 1:len])
        totalvars += len
    end

    my_lowers, my_uppers, totalvars
end

"""Setup an optimization problem."""
function problem(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, lowers::Vector{T}, uppers::Vector{T}, objective::Function; constraints::Vector{Function}=Function[], algorithm::Symbol=:LN_COBYLA_OR_LD_MMA) where T <: Real
    my_lowers, my_uppers, totalvars = expandlimits(model, components, names, lowers, uppers)

    if algorithm == :GUROBI_LINPROG
        # Make no changes to objective!
    elseif model.md.number_type == Number
        if algorithm == :LN_COBYLA_OR_LD_MMA
            algorithm = :LD_MMA
        end
        if string(algorithm)[2] == 'N'
            warn("Model is autodifferentiable, but optimizing using a derivative-free algorithm.")
            myobjective = gradfreeobjective(model, components, names, objective)
        else
            println("Using AutoDiff objective.")
            myobjective = autodiffobjective(model, components, names, objective)
        end
    else
        if algorithm == :LN_COBYLA_OR_LD_MMA
            algorithm = :LN_COBYLA
        elseif string(algorithm)[2] == 'D'
            warn("Model is non-differentiable, but requested a gradient algorithm; instead using LN_COBYLA.")
            algorithm = :LN_COBYLA
        end

        myobjective = gradfreeobjective(model, components, names, objective)
    end

    if algorithm == :GUROBI_LINPROG
        LinprogOptimizationProblem(model, components, names, objective, constraints, MatrixConstraintSet[], my_lowers, my_uppers)
    else
        if algorithm in BlackBoxAlgorithms
            if length(constraints) > 0
                warn("Functional constraints not supported for BBO algorithms.")
            end

            BlackBoxOptimizationProblem(algorithm, model, components, names, myobjective, my_lowers, my_uppers)
        else
            opt = Opt(algorithm, totalvars)
            lower_bounds!(opt, my_lowers)
            upper_bounds!(opt, my_uppers)
            xtol_rel!(opt, 1e-6)

            max_objective!(opt, myobjective)

            for constraint in constraints
                let this_constraint = constraint
                    function my_constraint(xx::Vector, grad::Vector)
                        setparameters(model, components, names, xx)
                        run(model)
                        this_constraint(model)
                    end

                    inequality_constraint!(opt, my_constraint)
                end
            end

            OptimizationProblem(model, components, names, opt, constraints)
        end
    end
end

"""Solve an optimization problem."""
function solution(optprob::OptimizationProblem, generator::Function; maxiter=Inf, verbose=false)
    global allverbose
    allverbose = verbose

    if verbose
        println("Selecting an initial point.")
    end

    attempts = 0
    initial = []
    valid = false
    while attempts < maxiter
        initial = generator()

        setparameters(optprob.model, optprob.components, optprob.names, initial)
        run(optprob.model)

        valid = true
        for constraint in optprob.constraints
            if constraint(optprob.model) >= 0
                valid = false
                break
            end
        end

        if valid
            break
        end

        attempts += 1
        if attempts % 1000 == 0
            println("Could not find initial point after $attempts attempts.")
        end
    end

    if !valid
        throw(DomainError("Could not find a valid initial value."))
    end

    if verbose
        println("Optimizing...")
    end
    if maxiter < Inf
        maxeval!(optprob.opt, maxiter)
    end
    (minf,minx,ret) = optimize(optprob.opt, initial)

    (minf, minx)
end

"""Solve an optimization problem."""
function solution(optprob::BlackBoxOptimizationProblem; maxiter=Inf, verbose=false)
    global allverbose
    allverbose = verbose

    if verbose
        println("Optimizing...")
    end
    res = bboptimize(p -> -optprob.objective(p, 0 .* p); SearchRange=collect(zip(optprob.lowers, optprob.uppers)), Method=optprob.algorithm)

    (best_fitness(res), best_candidate(res))
end

"""Setup an optimization problem."""
function problem(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, lowers::Vector{T}, uppers::Vector{T}, objective::Function, objectiveconstraints::Vector{Function}, matrixconstraints::Vector{MatrixConstraintSet}) where T <: Real
    my_lowers, my_uppers, totalvars = expandlimits(model, components, names, lowers, uppers)
    LinprogOptimizationProblem(model, components, names, objective, objectiveconstraints, matrixconstraints, my_lowers, my_uppers)
end

"""Solve an optimization problem."""
function solution(optprob::LinprogOptimizationProblem, verbose=false)
    global allverbose
    allverbose = verbose

    initial = make0(optprob.model, optprob.components, optprob.names)

    if verbose
        println("Optimizing...")
    end

    if optprob.model.md.number_type == Number
        myobjective = unaryobjective(optprob.model, optprob.components, optprob.names, optprob.objective)
        f, b, A = lpconstraints(optprob.model, optprob.components, optprob.names, myobjective, objectiveconstraints)
    else
        f, b, A = lpconstraints(optprob.model, optprob.components, optprob.names, optprob.objective, optprob.objectiveconstraints)
    end

    f, b, A = combineconstraints(f, b, A, optprob.model, optprob.components, optprob.names, optprob.matrixconstraints)
    exlowers, exuppers = combinelimits(optprob.exlowers, optprob.exuppers, optprob.model, optprob.components, optprob.names, optprob.matrixconstraints)

    # Use -f, because linprog *minimizes* objective
    @time sol = linprog(-f, A, '<', b, optprob.exlowers, optprob.exuppers)

    sol.sol
end

end # module


text:
using Base.Test

t = 0:0.1:500
x = @. exp(2π*(-0.001+0.05im)*t) #+ exp(2π*(-0.002+0.03im)*t)
plot(t,real.(x))

import Harminv
reload("Harminv")

d = Harminv.HarminvData(x, 0.0, 0.05, 100)
Harminv.solve_once!(d)
Harminv.solve_again!(d)
Harminv.compute_amplitudes(d)


text:
abstract type AbstractStopRemoteControl end

"""
Turn a boolean to false to cancel this check in the functions stop! and start!.
"""
struct StopRemoteControl <: AbstractStopRemoteControl
  unbounded_and_domain_x_check::Bool
  domain_check::Bool
  optimality_check::Bool
  infeasibility_check::Bool
  unbounded_problem_check::Bool
  tired_check::Bool
  resources_check::Bool
  stalled_check::Bool
  iteration_check::Bool
  main_pb_check::Bool
  user_check::Bool
  user_start_check::Bool

  cheap_check::Bool #`stop!` and `start!` stop whenever one check worked
  #not used now

end

function StopRemoteControl(;
  unbounded_and_domain_x_check::Bool = true, #O(n)
  domain_check::Bool = true, #O(n)
  optimality_check::Bool = true,
  infeasibility_check::Bool = true,
  unbounded_problem_check::Bool = true, #O(n)
  tired_check::Bool = true,
  resources_check::Bool = true,
  stalled_check::Bool = true,
  iteration_check::Bool = true,
  main_pb_check::Bool = true, #O(n)
  user_check::Bool = true,
  user_start_check::Bool = true,
  cheap_check::Bool = false,
)
  return StopRemoteControl(
    unbounded_and_domain_x_check,
    domain_check,
    optimality_check,
    infeasibility_check,
    unbounded_problem_check,
    tired_check,
    resources_check,
    stalled_check,
    iteration_check,
    main_pb_check,
    user_check,
    user_start_check,
    cheap_check,
  )
end

"""
Return a StopRemoteControl with the most expansive checks at false (the O(n))
by default in Stopping when it has a main_stp.
"""
function cheap_stop_remote_control(;
  unbounded_and_domain_x_check::Bool = false,
  domain_check::Bool = false,
  optimality_check::Bool = true,
  infeasibility_check::Bool = true,
  unbounded_problem_check::Bool = false,
  tired_check::Bool = true,
  resources_check::Bool = true,
  stalled_check::Bool = true,
  iteration_check::Bool = true,
  main_pb_check::Bool = false,
  user_check::Bool = true,
  user_start_check::Bool = true,
  cheap_check::Bool = true,
)
  return StopRemoteControl(
    unbounded_and_domain_x_check,
    domain_check,
    optimality_check,
    infeasibility_check,
    unbounded_problem_check,
    tired_check,
    resources_check,
    stalled_check,
    iteration_check,
    main_pb_check,
    user_check,
    user_start_check,
    cheap_check,
  )
end

import Base.show
function show(io::IO, src::StopRemoteControl)
  varlines = "$(typeof(src)) controls the following checks \n"
  for f in fieldnames(typeof(src))
    varlines = string(varlines, @sprintf("%28s: %s \n", f, getfield(src, f)))
  end
  println(io, varlines)
end


text:
using CalibrationErrors
using Random

Random.seed!(1234)

@testset "Trivial tests" begin
    ece = ECE(UniformBinning(10))

    @test calibrationerror(ece, ([0 1; 1 0], [2, 1])) == 0
    @test calibrationerror(ece, ([0 0.5 0.5 1; 1 0.5 0.5 0], [2, 2, 1, 1])) == 0
end


text:
# This file was generated, do not modify it. # hide
figure(figsize=(8, 6))
hist(res, color="blue", edgecolor="white", bins=50,
     density=true, alpha=0.5)

savefig(joinpath(@OUTPUT, "hist_residuals.svg")) # hide


text:
# some convenience functions
function structname(equ)
    if occursin('.', equ)
        return equ[findlast(isequal('.'), equ)+1:end]
    else
        return equ
    end
end

teststring(equ) = structname(string(equ))
teststring(equ, pert) = structname(string(equ)) * " + " * structname(string(pert))

# testing parameters
const t = 1.
const ξ = [1.05, 0.5, 0.5]

macro test_equilibrium(equilibrium_module, equilibrium_periodicity)
    module_name = string(equilibrium_module)
    module_test = Symbol(module_name * "Test")

    quote
        @eval module $module_test
            using Test
            using LinearAlgebra

            # import analytic field module
            import ElectromagneticFields
            import ElectromagneticFields.$equilibrium_module

            # inject code
            $equilibrium_module.@code

            function test_equilibrium(t, ξ)
                @test ξ¹(t, ξ...) == ξ¹(t,ξ)
                @test ξ²(t, ξ...) == ξ²(t,ξ)
                @test ξ³(t, ξ...) == ξ³(t,ξ)

                @test x¹(t, ξ...) == x¹(t,ξ)
                @test x²(t, ξ...) == x²(t,ξ)
                @test x³(t, ξ...) == x³(t,ξ)

                @test g₁₁(t, ξ...) == g₁₁(t, ξ)
                @test g₁₂(t, ξ...) == g₁₂(t, ξ)
                @test g₁₃(t, ξ...) == g₁₃(t, ξ)
                @test g₂₁(t, ξ...) == g₂₁(t, ξ)
                @test g₂₂(t, ξ...) == g₂₂(t, ξ)
                @test g₂₃(t, ξ...) == g₂₃(t, ξ)
                @test g₃₁(t, ξ...) == g₃₁(t, ξ)
                @test g₃₂(t, ξ...) == g₃₂(t, ξ)
                @test g₃₃(t, ξ...) == g₃₃(t, ξ)

                @test g¹¹(t, ξ...) == g¹¹(t,ξ)
                @test g¹²(t, ξ...) == g¹²(t,ξ)
                @test g¹³(t, ξ...) == g¹³(t,ξ)
                @test g²¹(t, ξ...) == g²¹(t,ξ)
                @test g²²(t, ξ...) == g²²(t,ξ)
                @test g²³(t, ξ...) == g²³(t,ξ)
                @test g³¹(t, ξ...) == g³¹(t,ξ)
                @test g³²(t, ξ...) == g³²(t,ξ)
                @test g³³(t, ξ...) == g³³(t,ξ)

                @test DF₁₁(t, ξ...) == DF₁₁(t, ξ)
                @test DF₁₂(t, ξ...) == DF₁₂(t, ξ)
                @test DF₁₃(t, ξ...) == DF₁₃(t, ξ)
                @test DF₂₁(t, ξ...) == DF₂₁(t, ξ)
                @test DF₂₂(t, ξ...) == DF₂₂(t, ξ)
                @test DF₂₃(t, ξ...) == DF₂₃(t, ξ)
                @test DF₃₁(t, ξ...) == DF₃₁(t, ξ)
                @test DF₃₂(t, ξ...) == DF₃₂(t, ξ)
                @test DF₃₃(t, ξ...) == DF₃₃(t, ξ)

                @test DF̄₁₁(t, ξ...) == DF̄₁₁(t, ξ)
                @test DF̄₁₂(t, ξ...) == DF̄₁₂(t, ξ)
                @test DF̄₁₃(t, ξ...) == DF̄₁₃(t, ξ)
                @test DF̄₂₁(t, ξ...) == DF̄₂₁(t, ξ)
                @test DF̄₂₂(t, ξ...) == DF̄₂₂(t, ξ)
                @test DF̄₂₃(t, ξ...) == DF̄₂₃(t, ξ)
                @test DF̄₃₁(t, ξ...) == DF̄₃₁(t, ξ)
                @test DF̄₃₂(t, ξ...) == DF̄₃₂(t, ξ)
                @test DF̄₃₃(t, ξ...) == DF̄₃₃(t, ξ)

                @test dg₁₁dx₁(t, ξ...) == dg₁₁dx₁(t, ξ)
                @test dg₁₂dx₁(t, ξ...) == dg₁₂dx₁(t, ξ)
                @test dg₁₃dx₁(t, ξ...) == dg₁₃dx₁(t, ξ)
                @test dg₂₁dx₁(t, ξ...) == dg₂₁dx₁(t, ξ)
                @test dg₂₂dx₁(t, ξ...) == dg₂₂dx₁(t, ξ)
                @test dg₂₃dx₁(t, ξ...) == dg₂₃dx₁(t, ξ)
                @test dg₃₁dx₁(t, ξ...) == dg₃₁dx₁(t, ξ)
                @test dg₃₂dx₁(t, ξ...) == dg₃₂dx₁(t, ξ)
                @test dg₃₃dx₁(t, ξ...) == dg₃₃dx₁(t, ξ)

                @test dg₁₁dx₂(t, ξ...) == dg₁₁dx₂(t, ξ)
                @test dg₁₂dx₂(t, ξ...) == dg₁₂dx₂(t, ξ)
                @test dg₁₃dx₂(t, ξ...) == dg₁₃dx₂(t, ξ)
                @test dg₂₁dx₂(t, ξ...) == dg₂₁dx₂(t, ξ)
                @test dg₂₂dx₂(t, ξ...) == dg₂₂dx₂(t, ξ)
                @test dg₂₃dx₂(t, ξ...) == dg₂₃dx₂(t, ξ)
                @test dg₃₁dx₂(t, ξ...) == dg₃₁dx₂(t, ξ)
                @test dg₃₂dx₂(t, ξ...) == dg₃₂dx₂(t, ξ)
                @test dg₃₃dx₂(t, ξ...) == dg₃₃dx₂(t, ξ)

                @test dg₁₁dx₃(t, ξ...) == dg₁₁dx₃(t, ξ)
                @test dg₁₂dx₃(t, ξ...) == dg₁₂dx₃(t, ξ)
                @test dg₁₃dx₃(t, ξ...) == dg₁₃dx₃(t, ξ)
                @test dg₂₁dx₃(t, ξ...) == dg₂₁dx₃(t, ξ)
                @test dg₂₂dx₃(t, ξ...) == dg₂₂dx₃(t, ξ)
                @test dg₂₃dx₃(t, ξ...) == dg₂₃dx₃(t, ξ)
                @test dg₃₁dx₃(t, ξ...) == dg₃₁dx₃(t, ξ)
                @test dg₃₂dx₃(t, ξ...) == dg₃₂dx₃(t, ξ)
                @test dg₃₃dx₃(t, ξ...) == dg₃₃dx₃(t, ξ)

                @test dg¹¹dx₁(t, ξ...) == dg¹¹dx₁(t,ξ)
                @test dg¹²dx₁(t, ξ...) == dg¹²dx₁(t,ξ)
                @test dg¹³dx₁(t, ξ...) == dg¹³dx₁(t,ξ)
                @test dg²¹dx₁(t, ξ...) == dg²¹dx₁(t,ξ)
                @test dg²²dx₁(t, ξ...) == dg²²dx₁(t,ξ)
                @test dg²³dx₁(t, ξ...) == dg²³dx₁(t,ξ)
                @test dg³¹dx₁(t, ξ...) == dg³¹dx₁(t,ξ)
                @test dg³²dx₁(t, ξ...) == dg³²dx₁(t,ξ)
                @test dg³³dx₁(t, ξ...) == dg³³dx₁(t,ξ)

                @test dg¹¹dx₂(t, ξ...) == dg¹¹dx₂(t,ξ)
                @test dg¹²dx₂(t, ξ...) == dg¹²dx₂(t,ξ)
                @test dg¹³dx₂(t, ξ...) == dg¹³dx₂(t,ξ)
                @test dg²¹dx₂(t, ξ...) == dg²¹dx₂(t,ξ)
                @test dg²²dx₂(t, ξ...) == dg²²dx₂(t,ξ)
                @test dg²³dx₂(t, ξ...) == dg²³dx₂(t,ξ)
                @test dg³¹dx₂(t, ξ...) == dg³¹dx₂(t,ξ)
                @test dg³²dx₂(t, ξ...) == dg³²dx₂(t,ξ)
                @test dg³³dx₂(t, ξ...) == dg³³dx₂(t,ξ)

                @test dg¹¹dx₃(t, ξ...) == dg¹¹dx₃(t,ξ)
                @test dg¹²dx₃(t, ξ...) == dg¹²dx₃(t,ξ)
                @test dg¹³dx₃(t, ξ...) == dg¹³dx₃(t,ξ)
                @test dg²¹dx₃(t, ξ...) == dg²¹dx₃(t,ξ)
                @test dg²²dx₃(t, ξ...) == dg²²dx₃(t,ξ)
                @test dg²³dx₃(t, ξ...) == dg²³dx₃(t,ξ)
                @test dg³¹dx₃(t, ξ...) == dg³¹dx₃(t,ξ)
                @test dg³²dx₃(t, ξ...) == dg³²dx₃(t,ξ)
                @test dg³³dx₃(t, ξ...) == dg³³dx₃(t,ξ)

                @test X(t, ξ...) == X(t,ξ)
                @test Y(t, ξ...) == Y(t,ξ)
                @test Z(t, ξ...) == Z(t,ξ)

                @test J(t, ξ...) == J(t,ξ)
                @test B(t, ξ...) == B(t,ξ)
                @test φ(t, ξ...) == φ(t,ξ)

                @test A₁(t, ξ...) == A₁(t,ξ)
                @test A₂(t, ξ...) == A₂(t,ξ)
                @test A₃(t, ξ...) == A₃(t,ξ)

                @test B₁(t, ξ...) == B₁(t,ξ)
                @test B₂(t, ξ...) == B₂(t,ξ)
                @test B₃(t, ξ...) == B₃(t,ξ)

                @test a₁(t, ξ...) == a₁(t,ξ)
                @test a₂(t, ξ...) == a₂(t,ξ)
                @test a₃(t, ξ...) == a₃(t,ξ)

                @test b₁(t, ξ...) == b₁(t,ξ)
                @test b₂(t, ξ...) == b₂(t,ξ)
                @test b₃(t, ξ...) == b₃(t,ξ)

                @test c₁(t, ξ...) == c₁(t,ξ)
                @test c₂(t, ξ...) == c₂(t,ξ)
                @test c₃(t, ξ...) == c₃(t,ξ)

                @test E₁(t, ξ...) == E₁(t,ξ)
                @test E₂(t, ξ...) == E₂(t,ξ)
                @test E₃(t, ξ...) == E₃(t,ξ)

                @test A¹(t, ξ...) == A¹(t,ξ)
                @test A²(t, ξ...) == A²(t,ξ)
                @test A³(t, ξ...) == A³(t,ξ)

                @test B¹(t, ξ...) == B¹(t,ξ)
                @test B²(t, ξ...) == B²(t,ξ)
                @test B³(t, ξ...) == B³(t,ξ)

                @test a¹(t, ξ...) == a¹(t,ξ)
                @test a²(t, ξ...) == a²(t,ξ)
                @test a³(t, ξ...) == a³(t,ξ)

                @test b¹(t, ξ...) == b¹(t,ξ)
                @test b²(t, ξ...) == b²(t,ξ)
                @test b³(t, ξ...) == b³(t,ξ)

                @test c¹(t, ξ...) == c¹(t,ξ)
                @test c²(t, ξ...) == c²(t,ξ)
                @test c³(t, ξ...) == c³(t,ξ)

                @test E¹(t, ξ...) == E¹(t,ξ)
                @test E²(t, ξ...) == E²(t,ξ)
                @test E³(t, ξ...) == E³(t,ξ)

                @test B₍₁₎(t, ξ...) == B₍₁₎(t,ξ)
                @test B₍₂₎(t, ξ...) == B₍₂₎(t,ξ)
                @test B₍₃₎(t, ξ...) == B₍₃₎(t,ξ)

                @test a₍₁₎(t, ξ...) == a₍₁₎(t,ξ)
                @test a₍₂₎(t, ξ...) == a₍₂₎(t,ξ)
                @test a₍₃₎(t, ξ...) == a₍₃₎(t,ξ)

                @test b₍₁₎(t, ξ...) == b₍₁₎(t,ξ)
                @test b₍₂₎(t, ξ...) == b₍₂₎(t,ξ)
                @test b₍₃₎(t, ξ...) == b₍₃₎(t,ξ)

                @test c₍₁₎(t, ξ...) == c₍₁₎(t,ξ)
                @test c₍₂₎(t, ξ...) == c₍₂₎(t,ξ)
                @test c₍₃₎(t, ξ...) == c₍₃₎(t,ξ)

                @test dA₁dx₁(t, ξ...) == dA₁dx₁(t,ξ)
                @test dA₁dx₂(t, ξ...) == dA₁dx₂(t,ξ)
                @test dA₁dx₃(t, ξ...) == dA₁dx₃(t,ξ)

                @test dA₂dx₁(t, ξ...) == dA₂dx₁(t,ξ)
                @test dA₂dx₂(t, ξ...) == dA₂dx₂(t,ξ)
                @test dA₂dx₃(t, ξ...) == dA₂dx₃(t,ξ)

                @test dA₃dx₁(t, ξ...) == dA₃dx₁(t,ξ)
                @test dA₃dx₂(t, ξ...) == dA₃dx₂(t,ξ)
                @test dA₃dx₃(t, ξ...) == dA₃dx₃(t,ξ)

                @test dB₁dx₁(t, ξ...) == dB₁dx₁(t,ξ)
                @test dB₁dx₂(t, ξ...) == dB₁dx₂(t,ξ)
                @test dB₁dx₃(t, ξ...) == dB₁dx₃(t,ξ)

                @test dB₂dx₁(t, ξ...) == dB₂dx₁(t,ξ)
                @test dB₂dx₂(t, ξ...) == dB₂dx₂(t,ξ)
                @test dB₂dx₃(t, ξ...) == dB₂dx₃(t,ξ)

                @test dB₃dx₁(t, ξ...) == dB₃dx₁(t,ξ)
                @test dB₃dx₂(t, ξ...) == dB₃dx₂(t,ξ)
                @test dB₃dx₃(t, ξ...) == dB₃dx₃(t,ξ)

                @test db₁dx₁(t, ξ...) == db₁dx₁(t,ξ)
                @test db₁dx₂(t, ξ...) == db₁dx₂(t,ξ)
                @test db₁dx₃(t, ξ...) == db₁dx₃(t,ξ)

                @test db₂dx₁(t, ξ...) == db₂dx₁(t,ξ)
                @test db₂dx₂(t, ξ...) == db₂dx₂(t,ξ)
                @test db₂dx₃(t, ξ...) == db₂dx₃(t,ξ)

                @test db₃dx₁(t, ξ...) == db₃dx₁(t,ξ)
                @test db₃dx₂(t, ξ...) == db₃dx₂(t,ξ)
                @test db₃dx₃(t, ξ...) == db₃dx₃(t,ξ)

                @test dBdx₁(t, ξ...) == dBdx₁(t,ξ)
                @test dBdx₂(t, ξ...) == dBdx₂(t,ξ)
                @test dBdx₃(t, ξ...) == dBdx₃(t,ξ)

                @test d²A₁dx₁dx₁(t, ξ...) == d²A₁dx₁dx₁(t,ξ)
                @test d²A₁dx₁dx₂(t, ξ...) == d²A₁dx₁dx₂(t,ξ)
                @test d²A₁dx₁dx₃(t, ξ...) == d²A₁dx₁dx₃(t,ξ)

                @test d²A₁dx₂dx₁(t, ξ...) == d²A₁dx₂dx₁(t,ξ)
                @test d²A₁dx₂dx₂(t, ξ...) == d²A₁dx₂dx₂(t,ξ)
                @test d²A₁dx₂dx₃(t, ξ...) == d²A₁dx₂dx₃(t,ξ)

                @test d²A₁dx₃dx₁(t, ξ...) == d²A₁dx₃dx₁(t,ξ)
                @test d²A₁dx₃dx₂(t, ξ...) == d²A₁dx₃dx₂(t,ξ)
                @test d²A₁dx₃dx₃(t, ξ...) == d²A₁dx₃dx₃(t,ξ)

                @test d²A₂dx₁dx₁(t, ξ...) == d²A₂dx₁dx₁(t,ξ)
                @test d²A₂dx₁dx₂(t, ξ...) == d²A₂dx₁dx₂(t,ξ)
                @test d²A₂dx₁dx₃(t, ξ...) == d²A₂dx₁dx₃(t,ξ)

                @test d²A₂dx₂dx₁(t, ξ...) == d²A₂dx₂dx₁(t,ξ)
                @test d²A₂dx₂dx₂(t, ξ...) == d²A₂dx₂dx₂(t,ξ)
                @test d²A₂dx₂dx₃(t, ξ...) == d²A₂dx₂dx₃(t,ξ)

                @test d²A₂dx₃dx₁(t, ξ...) == d²A₂dx₃dx₁(t,ξ)
                @test d²A₂dx₃dx₂(t, ξ...) == d²A₂dx₃dx₂(t,ξ)
                @test d²A₂dx₃dx₃(t, ξ...) == d²A₂dx₃dx₃(t,ξ)

                @test d²A₃dx₁dx₁(t, ξ...) == d²A₃dx₁dx₁(t,ξ)
                @test d²A₃dx₁dx₂(t, ξ...) == d²A₃dx₁dx₂(t,ξ)
                @test d²A₃dx₁dx₃(t, ξ...) == d²A₃dx₁dx₃(t,ξ)

                @test d²A₃dx₂dx₁(t, ξ...) == d²A₃dx₂dx₁(t,ξ)
                @test d²A₃dx₂dx₂(t, ξ...) == d²A₃dx₂dx₂(t,ξ)
                @test d²A₃dx₂dx₃(t, ξ...) == d²A₃dx₂dx₃(t,ξ)

                @test d²A₃dx₃dx₁(t, ξ...) == d²A₃dx₃dx₁(t,ξ)
                @test d²A₃dx₃dx₂(t, ξ...) == d²A₃dx₃dx₂(t,ξ)
                @test d²A₃dx₃dx₃(t, ξ...) == d²A₃dx₃dx₃(t,ξ)

                @test d²b₁dx₁dx₁(t, ξ...) == d²b₁dx₁dx₁(t,ξ)
                @test d²b₁dx₁dx₂(t, ξ...) == d²b₁dx₁dx₂(t,ξ)
                @test d²b₁dx₁dx₃(t, ξ...) == d²b₁dx₁dx₃(t,ξ)

                @test d²b₁dx₂dx₁(t, ξ...) == d²b₁dx₂dx₁(t,ξ)
                @test d²b₁dx₂dx₂(t, ξ...) == d²b₁dx₂dx₂(t,ξ)
                @test d²b₁dx₂dx₃(t, ξ...) == d²b₁dx₂dx₃(t,ξ)

                @test d²b₁dx₃dx₁(t, ξ...) == d²b₁dx₃dx₁(t,ξ)
                @test d²b₁dx₃dx₂(t, ξ...) == d²b₁dx₃dx₂(t,ξ)
                @test d²b₁dx₃dx₃(t, ξ...) == d²b₁dx₃dx₃(t,ξ)

                @test d²b₂dx₁dx₁(t, ξ...) == d²b₂dx₁dx₁(t,ξ)
                @test d²b₂dx₁dx₂(t, ξ...) == d²b₂dx₁dx₂(t,ξ)
                @test d²b₂dx₁dx₃(t, ξ...) == d²b₂dx₁dx₃(t,ξ)

                @test d²b₂dx₂dx₁(t, ξ...) == d²b₂dx₂dx₁(t,ξ)
                @test d²b₂dx₂dx₂(t, ξ...) == d²b₂dx₂dx₂(t,ξ)
                @test d²b₂dx₂dx₃(t, ξ...) == d²b₂dx₂dx₃(t,ξ)

                @test d²b₂dx₃dx₁(t, ξ...) == d²b₂dx₃dx₁(t,ξ)
                @test d²b₂dx₃dx₂(t, ξ...) == d²b₂dx₃dx₂(t,ξ)
                @test d²b₂dx₃dx₃(t, ξ...) == d²b₂dx₃dx₃(t,ξ)

                @test d²b₃dx₁dx₁(t, ξ...) == d²b₃dx₁dx₁(t,ξ)
                @test d²b₃dx₁dx₂(t, ξ...) == d²b₃dx₁dx₂(t,ξ)
                @test d²b₃dx₁dx₃(t, ξ...) == d²b₃dx₁dx₃(t,ξ)

                @test d²b₃dx₂dx₁(t, ξ...) == d²b₃dx₂dx₁(t,ξ)
                @test d²b₃dx₂dx₂(t, ξ...) == d²b₃dx₂dx₂(t,ξ)
                @test d²b₃dx₂dx₃(t, ξ...) == d²b₃dx₂dx₃(t,ξ)

                @test d²b₃dx₃dx₁(t, ξ...) == d²b₃dx₃dx₁(t,ξ)
                @test d²b₃dx₃dx₂(t, ξ...) == d²b₃dx₃dx₂(t,ξ)
                @test d²b₃dx₃dx₃(t, ξ...) == d²b₃dx₃dx₃(t,ξ)

                @test d²Bdx₁dx₁(t, ξ...) == d²Bdx₁dx₁(t,ξ)
                @test d²Bdx₁dx₂(t, ξ...) == d²Bdx₁dx₂(t,ξ)
                @test d²Bdx₁dx₃(t, ξ...) == d²Bdx₁dx₃(t,ξ)

                @test d²Bdx₂dx₁(t, ξ...) == d²Bdx₂dx₁(t,ξ)
                @test d²Bdx₂dx₂(t, ξ...) == d²Bdx₂dx₂(t,ξ)
                @test d²Bdx₂dx₃(t, ξ...) == d²Bdx₂dx₃(t,ξ)

                @test d²Bdx₃dx₁(t, ξ...) == d²Bdx₃dx₁(t,ξ)
                @test d²Bdx₃dx₂(t, ξ...) == d²Bdx₃dx₂(t,ξ)
                @test d²Bdx₃dx₃(t, ξ...) == d²Bdx₃dx₃(t,ξ)

                @test periodicity(t,ξ...) == $equilibrium_periodicity
                @test periodicity(t,ξ)    == $equilibrium_periodicity
                @test periodicity(ξ...)   == $equilibrium_periodicity
                @test periodicity(ξ)      == $equilibrium_periodicity

                # check internal consistency
                @test from_cartesian(t, to_cartesian(t,ξ)) ≈ ξ  atol=1E-14

                let g = g(t,ξ), ḡ = ḡ(t,ξ), DF = DF(t,ξ), DF̄ = DF̄(t,ξ),
                    a = a(t,ξ), b = b(t,ξ), c = c(t,ξ),
                    a⃗ = a⃗(t,ξ), b⃗ = b⃗(t,ξ), c⃗ = c⃗(t,ξ),
                    â = aₚ(t,ξ), b̂ = bₚ(t,ξ), ĉ = cₚ(t,ξ)

                    @test J(t,ξ) ≈ sqrt(det(DF' * DF))  atol=1E-12
                    @test ḡ         ≈ inv(g)          atol=1E-12
                    @test DF̄        ≈ inv(DF)         atol=1E-12
                    @test DF' * DF  ≈ g               atol=1E-12
                    @test DF  * DF̄  ≈ Array(I, 3, 3)  atol=1E-12
                    @test DF̄  * DF̄' ≈ ḡ               atol=1E-12

                    if $equilibrium_module != ElectromagneticFields.Singular
                        @test g * a⃗ ≈ a  atol=1E-14
                        @test g * b⃗ ≈ b  atol=1E-14
                        @test g * c⃗ ≈ c  atol=1E-14

                        @test ḡ * a ≈ a⃗  atol=1E-14
                        @test ḡ * b ≈ b⃗  atol=1E-14
                        @test ḡ * c ≈ c⃗  atol=1E-14

                        @test â ≈ DF * a⃗  atol=1E-14
                        @test b̂ ≈ DF * b⃗  atol=1E-14
                        @test ĉ ≈ DF * c⃗  atol=1E-14

                        @test â ≈ DF̄' * a  atol=1E-14
                        @test b̂ ≈ DF̄' * b  atol=1E-14
                        @test ĉ ≈ DF̄' * c  atol=1E-14

                        @test a⃗' * a ≈ 1  atol=1E-14
                        @test b⃗' * b ≈ 1  atol=1E-14
                        @test c⃗' * c ≈ 1  atol=1E-14

                        @test â' * â ≈ 1  atol=1E-14
                        @test b̂' * b̂ ≈ 1  atol=1E-14
                        @test ĉ' * ĉ ≈ 1  atol=1E-14

                        @test â' * b̂ ≈ 0  atol=1E-14
                        @test b̂' * ĉ ≈ 0  atol=1E-14
                        @test ĉ' * â ≈ 0  atol=1E-14

                        @test a' * ḡ * a ≈ 1  atol=1E-14
                        @test b' * ḡ * b ≈ 1  atol=1E-14
                        @test c' * ḡ * c ≈ 1  atol=1E-14

                        @test a' * ḡ * b ≈ 0  atol=1E-14
                        @test b' * ḡ * c ≈ 0  atol=1E-14
                        @test c' * ḡ * a ≈ 0  atol=1E-14

                        @test a⃗' * g * a⃗ ≈ 1  atol=1E-14
                        @test b⃗' * g * b⃗ ≈ 1  atol=1E-14
                        @test c⃗' * g * c⃗ ≈ 1  atol=1E-14

                        @test a⃗' * g * b⃗ ≈ 0  atol=1E-14
                        @test b⃗' * g * c⃗ ≈ 0  atol=1E-14
                        @test c⃗' * g * a⃗ ≈ 0  atol=1E-14
                    end
                end
            end
        end

        # run tests
        @testset "$(rpad($module_name,60))" begin
            import .$module_test
            $module_test.test_equilibrium(t, ξ)
        end
    end
end

# perturbation list (equilibrium, parameters, perturbation, parameters, module)
perts = (
    (SymmetricQuadratic,    EzCosZ,    (2.)),
    (ThetaPinch,            EzCosZ,    (2.)),
)

# test equilibria

@test_equilibrium ABC                                           [0., 0., 0.]
@test_equilibrium AxisymmetricTokamakCartesian                  [0., 0., 0.]
@test_equilibrium AxisymmetricTokamakCylindrical                [0., 0., 2π]
@test_equilibrium AxisymmetricTokamakToroidal                   [0., 2π, 2π]
@test_equilibrium AxisymmetricTokamakToroidalRegularization     [0., 2π, 2π]
@test_equilibrium PenningTrapUniform                            [0., 0., 0.]
@test_equilibrium PenningTrapBottle                             [0., 0., 0.]
@test_equilibrium PenningTrapAsymmetric                         [0., 0., 0.]
@test_equilibrium Singular                                      [0., 0., 0.]
@test_equilibrium SymmetricQuadratic                            [0., 0., 0.]
@test_equilibrium ThetaPinch                                    [0., 0., 0.]
# @test_equilibrium SolovevFRC                                    [0., 0., 2π]
@test_equilibrium SolovevITER                                   [0., 0., 2π]
@test_equilibrium SolovevITERwXpoint                            [0., 0., 2π]
@test_equilibrium SolovevNSTX                                   [0., 0., 2π]
@test_equilibrium SolovevNSTXwXpoint                            [0., 0., 2π]
@test_equilibrium SolovevNSTXwDoubleXpoint                      [0., 0., 2π]
@test_equilibrium SolovevSymmetric                              [0., 0., 0.]
println()

# test perturbations
# @testset "$(rpad(teststring(equ[1], equ[3]),60))" for equ in perts begin
#         equ_obj = equ[1].init(equ[2]..., perturbation=equ[3].init(equ[4]...))
#         # test_equilibrium(equ[1], t, x)
#     end
# end
# println()

# test correctness of some of the magnetic fields

function test_axisymmetric_tokamak_cartesian_equilibrium(equ_mod, t=0., x=[1.5, 0.0, 0.5])
    @test equ_mod.B¹(t,x) ≈ equ_mod.B₁(t,x)   atol=1E-16
    @test equ_mod.B²(t,x) ≈ equ_mod.B₂(t,x)   atol=1E-16
    @test equ_mod.B³(t,x) ≈ equ_mod.B₃(t,x)   atol=1E-16

    @test equ_mod.B₁(t,x) ≈ - equ_mod.B₀ / equ_mod.q₀ * (equ_mod.q₀ * equ_mod.R₀ * equ_mod.Y(t,x) + equ_mod.X(t,x) * equ_mod.Z(t,x) ) / equ_mod.R(t,x)^2   atol=1E-16
    @test equ_mod.B₂(t,x) ≈ + equ_mod.B₀ / equ_mod.q₀ * (equ_mod.q₀ * equ_mod.R₀ * equ_mod.X(t,x) - equ_mod.Y(t,x) * equ_mod.Z(t,x) ) / equ_mod.R(t,x)^2   atol=1E-16
    @test equ_mod.B₃(t,x) ≈ + equ_mod.B₀ / equ_mod.q₀ * (equ_mod.R(t,x) - equ_mod.R₀) / equ_mod.R(t,x)   atol=1E-16
end

function test_axisymmetric_tokamak_cylindrical_equilibrium(equ_mod, t=0., x=[1.5, 0.5, π/5])
    @test equ_mod.B¹(t,x) == + equ_mod.B₀ / equ_mod.q₀ * equ_mod.Z(t,x) / equ_mod.R(t,x)
    @test equ_mod.B²(t,x) == - equ_mod.B₀ / equ_mod.q₀ * (equ_mod.R(t,x) - equ_mod.R₀) / equ_mod.R(t,x)
    @test equ_mod.B³(t,x) == - equ_mod.B₀ * equ_mod.R₀ / equ_mod.R(t,x)^2

    @test equ_mod.B₁(t,x) == + equ_mod.B₀ / equ_mod.q₀ * equ_mod.Z(t,x) / equ_mod.R(t,x)
    @test equ_mod.B₂(t,x) == - equ_mod.B₀ / equ_mod.q₀ * (equ_mod.R(t,x) - equ_mod.R₀) / equ_mod.R(t,x)
    @test equ_mod.B₃(t,x) == - equ_mod.B₀ * equ_mod.R₀
end

function test_axisymmetric_tokamak_toroidal_equilibrium(equ_mod, t=0., x=[0.5, π/10, π/5])
    @test equ_mod.B¹(t,x) == 0
    @test equ_mod.B²(t,x) == - equ_mod.B₀ / equ_mod.q₀ / equ_mod.R(t,x)
    @test equ_mod.B³(t,x) ≈  - equ_mod.B₀ * equ_mod.R₀ / equ_mod.R(t,x)^2 atol=1E-14

    @test equ_mod.B₁(t,x) == 0
    @test equ_mod.B₂(t,x) == - equ_mod.B₀ / equ_mod.q₀ * equ_mod.r(t,x)^2 / equ_mod.R(t,x)
    @test equ_mod.B₃(t,x) ≈  - equ_mod.B₀ * equ_mod.R₀ atol=1E-14
end

function test_consistency_axisymmetric_tokamak_cylindrical_equilibrium(equ_cyl, equ_car, t=0., ξ=[1.5, 0.5, π/5])
    x  = equ_cyl.to_cartesian(t,ξ)
    DF = equ_cyl.DF(t,ξ)
    DF̄ = equ_cyl.DF̄(t,ξ)

    B_cyl = [equ_cyl.B¹(t,ξ), equ_cyl.B²(t,ξ), equ_cyl.B³(t,ξ)]
    B_car = [equ_car.B¹(t,x), equ_car.B²(t,x), equ_car.B³(t,x)]

    B̂_cyl = [equ_cyl.B₁(t,ξ), equ_cyl.B₂(t,ξ), equ_cyl.B₃(t,ξ)]
    B̂_car = [equ_car.B₁(t,x), equ_car.B₂(t,x), equ_car.B₃(t,x)]

    @test B_cyl' * B̂_cyl ≈ B_car' * B̂_car  atol=1E-12
    @test DF  * B_cyl ≈ - B_car  atol=1E-12
    @test DF̄' * B̂_cyl ≈ - B̂_car  atol=1E-12
end

function test_consistency_axisymmetric_tokamak_toroidal_equilibrium(equ_tor, equ_car, t=0., ξ=[0.5, π/10, π/5])
    x  = equ_tor.to_cartesian(t,ξ)
    DF = equ_tor.DF(t,ξ)
    DF̄ = equ_tor.DF̄(t,ξ)

    B_tor = [equ_tor.B¹(t,ξ), equ_tor.B²(t,ξ), equ_tor.B³(t,ξ)]
    B_car = [equ_car.B¹(t,x), equ_car.B²(t,x), equ_car.B³(t,x)]

    B̂_tor = [equ_tor.B₁(t,ξ), equ_tor.B₂(t,ξ), equ_tor.B₃(t,ξ)]
    B̂_car = [equ_car.B₁(t,x), equ_car.B₂(t,x), equ_car.B₃(t,x)]

    @test B_tor' * B̂_tor ≈ B_car' * B̂_car  atol=1E-12
    @test DF  * B_tor ≈ - B_car  atol=1E-12
    @test DF̄' * B̂_tor ≈ - B̂_car  atol=1E-12
end

function test_symmetric_quadratic_equilibrium(equ_mod, t=0., x=[1.0, 0.5, 0.5])
    @test equ_mod.B¹(t,x) == equ_mod.B₁(t,x)
    @test equ_mod.B²(t,x) == equ_mod.B₂(t,x)
    @test equ_mod.B³(t,x) == equ_mod.B₃(t,x)

    @test equ_mod.B₁(t,x) == 0
    @test equ_mod.B₂(t,x) == 0
    @test equ_mod.B₃(t,x) == equ_mod.B₀ * (1 + equ_mod.X(t,x)^2 + equ_mod.Y(t,x)^2)

    @test equ_mod.B(t,x)  == equ_mod.B₀ * (1 + equ_mod.X(t,x)^2 + equ_mod.Y(t,x)^2)

    @test equ_mod.b¹(t,x) == 0
    @test equ_mod.b²(t,x) == 0
    @test equ_mod.b³(t,x) == 1

    @test equ_mod.b₁(t,x) == 0
    @test equ_mod.b₂(t,x) == 0
    @test equ_mod.b₃(t,x) == 1
end

function test_theta_pinch_equilibrium(equ_mod, t=0., x=[1.0, 0.5, 0.5])
    @test equ_mod.B¹(t,x) == equ_mod.B₁(t,x)
    @test equ_mod.B²(t,x) == equ_mod.B₂(t,x)
    @test equ_mod.B³(t,x) == equ_mod.B₃(t,x)

    @test equ_mod.B₁(t,x) == 0
    @test equ_mod.B₂(t,x) == 0
    @test equ_mod.B₃(t,x) == equ_mod.B₀

    @test equ_mod.B(t,x)  == equ_mod.B₀

    @test equ_mod.b¹(t,x) == 0
    @test equ_mod.b²(t,x) == 0
    @test equ_mod.b³(t,x) == 1

    @test equ_mod.b₁(t,x) == 0
    @test equ_mod.b₂(t,x) == 0
    @test equ_mod.b₃(t,x) == 1
end

function test_abc_equilibrium(equ_mod, t=0., x=[1.0, 0.5, 0.5])
    @test equ_mod.B¹(t,x) == equ_mod.B₁(t,x)
    @test equ_mod.B²(t,x) == equ_mod.B₂(t,x)
    @test equ_mod.B³(t,x) == equ_mod.B₃(t,x)

    @test equ_mod.B₁(t,x) == equ_mod.A₁(t,x)
    @test equ_mod.B₂(t,x) == equ_mod.A₂(t,x)
    @test equ_mod.B₃(t,x) == equ_mod.A₃(t,x)

    @test equ_mod.B₁(t,x) == equ_mod.a₀ * sin(equ_mod.Z(t,x)) + equ_mod.c₀ * cos(equ_mod.Y(t,x))
    @test equ_mod.B₂(t,x) == equ_mod.b₀ * sin(equ_mod.X(t,x)) + equ_mod.a₀ * cos(equ_mod.Z(t,x))
    @test equ_mod.B₃(t,x) == equ_mod.c₀ * sin(equ_mod.Y(t,x)) + equ_mod.b₀ * cos(equ_mod.X(t,x))
end

@testset "$(rpad("Magnetic Fields",60))" begin
    test_axisymmetric_tokamak_cartesian_equilibrium(AxisymmetricTokamakCartesianTest)
    test_axisymmetric_tokamak_cylindrical_equilibrium(AxisymmetricTokamakCylindricalTest)
    test_axisymmetric_tokamak_toroidal_equilibrium(AxisymmetricTokamakToroidalTest)
    test_symmetric_quadratic_equilibrium(SymmetricQuadraticTest)
    test_theta_pinch_equilibrium(ThetaPinchTest)
    test_abc_equilibrium(ABCTest)
end

@testset "$(rpad("Consistency",60))" begin
    test_consistency_axisymmetric_tokamak_cylindrical_equilibrium(AxisymmetricTokamakCylindricalTest, AxisymmetricTokamakCartesianTest)
    test_consistency_axisymmetric_tokamak_toroidal_equilibrium(AxisymmetricTokamakToroidalTest, AxisymmetricTokamakCartesianTest)
end

println()


text:
# Practical Summary: Ch. 16, GMM

## DGP: Poisson
using Distributions
n = 30
x = [ones(n) rand(n)]
β⁰ = [1., 2.]
λ = β -> exp.(x*β)
y = rand.(Poisson.(λ(β⁰)))

## ML, for reference
using Econometrics, SpecialFunctions
model = β -> -λ(β) + y.*x*β - loggamma.(y .+ 1.)    # note: loggamma(y+1) = log(factorial(y))
mleresults(model, β⁰)                               # but won't overflow

## GMM
using Econometrics
moments1 = β -> x.*(y - λ(β))
gmmresults(moments1, β⁰)
# note that the GMM estimator is identical to the ML estimator. You should
# be able to prove that result analytically. Because this GMM estimator is
# equivalent to ML, it is asymptotically efficient. Adding moment conditions
# won't improve asymptotic efficiency, even if they're valid.

## Trying a different GMM estimator, based on E(y)=λ, so E(y/λ)=1.
using Econometrics
moments2 = β -> x.* (y./λ(β) .- 1.0)
gmmresults(moments2, β⁰)

## Try out overidentified GMM estimator, using both sets of moments
using Econometrics, ForwardDiff
moments3 = β -> [moments1(β) moments2(β)]
βhat, objv, V, D, W, convergence = gmmresults(moments3, β⁰)

## how to get D and Omega (though gmmresults will also give them)
avgmoments = β -> (1/n)*[x'*(y - λ(β)) x'*(y./λ(β) .- 1.)]
D = ForwardDiff.jacobian(avgmoments,βhat)'
m = moments(βhat) # the moment contributions, evaluated at estimate
Ωhat = NeweyWest(m)

## let's see how the Hansen-Sargan test can detect
# incorrect moments
η = 0.01 # if this is different from zero, moments4, and, thus moments5, will not be valid
moments4 = β -> x.* (y./λ(β) .- 1.0 .+ η)
moments5 = β -> [moments1(β) moments4(β)]
gmmresults(moments5, β⁰)

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
"""
    create_surfaces(setup::PlateSetup, cells)

Create surface triangulation of a collection of stacked plates.
"""
function create_surfaces(setup::PlateSetup{T}, _) where {T}
    (; width, depth, heights) = setup

    n = length(heights)
    h = sum(heights) / 2
    z = cumsum([-h; heights])
    x = width / 2
    y = depth / 2

    points = zeros(T, 3, 4(n + 1))
    for i = 1:(n+1)
        points[:, (1+4(i-1)):(4i)] = [
            -x x x -x
            -y -y y y
            z[i] z[i] z[i] z[i]
        ]
    end

    facets = fill(0, 3, 2 + 10n)
    facets[:, 1:2] = [
        1 1
        2 3
        3 4
    ]
    for i = 1:n
        inds = 2 + 10(i - 1) .+ (1:10)
        facets[:, inds] .=
            4(i - 1) .+ [
                1 1 2 2 3 3 4 4 6 6
                2 6 3 7 4 8 1 5 7 8
                6 5 7 6 8 7 5 8 8 5
            ]
    end

    facetmarkers = fill(0, 2 + 10n)

    # Bottom and first side walls
    facetmarkers[1:10] .= n

    # Interfaces and side walls after interface
    for i = 1:(n-1)
        facetmarkers[10i+1] = i
        facetmarkers[10i+2] = i
        facetmarkers[(10i+3):(10(i+1))] .= n + i
    end

    # Top
    facetmarkers[(end-1):end] .= 2n - 1

    regions = zeros(T, 3, n)
    regions[3, :] .= (z[1:(end-1)] .+ z[2:end]) ./ 2

    (; points, facets, facetmarkers, regions)
end


text:
# include("model.jl")
include("path.jl")
N = 101
s = range(0,10,length=101)
e = sin.(s*10)*0.2

line = StraightPath(10, pi/4)
x,y = localToGlobal(line, s, e)
plot(Vector(x),Vector(y))

arc = ArcPath(-0*pi/4, 10, pi/3)
total_length(arc)
position_change(arc)
s = range(0,total_length(arc), length=N)
e = sin.(s*10)*0.2
x,y = localToGlobal(arc, s, e)
plot(x,y, aspect_ratio=:equal)
rad2deg(final_angle(arc))

line = StraightPath(10, pi/2)
arc = ArcPath(line, 10, pi)
line2 = StraightPath(arc, 12)
paths = [line,arc,line2]
paths = (line,arc,line2)
path = DubinsPath(paths)
plot(path, aspect_ratio=:equal)

s = range(0,total_length(path), length=N)
e = sin.(s*10)*0.2
e = fill(-0.2, N)
x,y = localToGlobal(path, s, e)
plot(x,y, aspect_ratio=:equal)

using BenchmarkTools
s = 1.2
@btime curvature($path, 1.2)
@code_warntype curvature(path, 1.2)
path


text:
# Graph States

using LightGraphs
import Base: convert
export GraphState, stabilizers, isLCeq

type GraphState
    graph::Graph
end

# operators

function stabilizers(gs::GraphState)
    res = [sigmax for i in vertices(gs.graph)]
    for node in vertices(gs.graph)
        for itr in vertices(gs.graph)
            if itr in neighbors(gs.graph,node)
                res[node] = kron(res[node],sigmaz)
            else
                res[node] = kron(res[node],sigmao)
            end
        end
    end
    return res
end

function isLCeq(gsa::GraphState,gsb::GraphState)
    @assert nv(gsa) == nv(gsb) "Number of nodes does not match"
    n = nv(gsa)
    generator_matrix_left = [speye(n) adjacency_matrix(gsa)]
    generator_matrix_right = [speye(n) adjacency_matrix(gsb)]
    P = [spzeros(n,n) speye(n);speye(n) spzeros(n,n)]
    flag = generator_matrix_left*P*generator_matrix_right.'
    return isempty(nonzeros(flag))
end

function convert(::Type{ComputState},gs::GraphState)
    return eigs(gs,nev=1,which=:SR)[2][:,1]|>ComputState
end

max_stars_repo_name: Roger-luo/openQC

max_issues_repo_name: Roger-luo/openQC

max_forks_repo_name: Roger-luo/openQC


text:
@testset "Config" begin
    include("config.jl")
end
@testset "Unit" begin
    include("unit.jl")
end
@testset "Continuous Linear" begin
    include("contlinear.jl")
end
@testset "Continuous Conic" begin
    include("contconic.jl")
end
@testset "Continuous Quadratic" begin
    include("contquadratic.jl")
end
@testset "Integer Linear" begin
    include("intlinear.jl")
end
@testset "Integer Conic" begin
    include("intconic.jl")
end


text:
using Documenter, IRKGaussLegendre.jl

makedocs()


text:
module TreeFlow
import GraphIdx: Vec, Weights, Graph, EdgeGraph, enumerate_edges

"""
    extract_non_tree!(graph, y, λt, π, α, λ)

Iterate over all edges of `grapy` (weighted by `λ`) and do the following:
If the edge `e` is a tree edge, store the corresponding `λ[e]` in tree order.
If the edge `e` is not part of the tree, compute the flow `α[e]` into the
node input `y`.
"""
function extract_non_tree!(
    graph::Graph,
    parent::Vector{Int},
    y::Array,
    alpha::Array,
    tlam::W1,
    lambda::W2,
) where {W1, W2<:Weights}
    enumerate_typed_edges(graph, parent) do istree::Bool, i::Int, u::Int, v::Int
        if istree
            tlam[u] = lambda[i]
        else
            y[u] += alpha[i]
            y[v] -= alpha[i]
        end
    end
end

"""

For convinience: Transform `λ` from Array into a `GraphIdx.Vec`.
"""
extract_non_tree!(
    graph::Graph,
    parent::Vector{Int},
    y::Array,
    alpha::Array,
    tlam::W1,
    lambda::Array{T},
) where {W1, T} =
    extract_non_tree!(
        graph::Graph,
        parent::Vector{Int},
        y::Array,
        alpha::Array,
        tlam::W1,
        GraphIdx.Vec(lambda),
    )

"""
    enumerate_typed_edges(func, edges, π)

Call for every edge `(u, v)` enumerated by `i` the function
```julia
func(istree, i, u, v)
```

whereby `istree` tells whether the edge is a tree edge.
If it is a tree edge, `u` is the child of `v`.
"""
@inline function enumerate_typed_edges(
    func::Func, graph::Graph, parent::Vector{Int}
) where {Func<:Function}
    enumerate_edges(graph) do i::Int, u::Int, v::Int
        if parent[v] == u
            func(true, i, v, u)
        elseif parent[u] == v
            func(true, i, u, v)
        else
            func(false, i, u, v)
        end
    end
end

"""
    update_tree!(α, αt, selected, graph, parent)

Update the global dual `α` by a tree dual `αt`.
"""
function update_tree!(α, αt, selected, graph::EdgeGraph, parent)
    @assert selected[1] < 0
    for i in @view selected[2:end]
        local u::Int, v::Int = graph.edges[i]
        α[i] = if parent[v] == u
            -αt[v]
        elseif parent[u] == v
            +αt[u]
        else
            error("Should not happen")
        end
    end
end

end # module

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
function get_equal_constraints(sdp::SparseSDP)
    pairs = Set()
    indices = collect(keys(cons(sdp)))
    for i = 1:length(indices), ip = i+1:length(indices)
        ci = indices[i]
        cip = indices[ip]
        bm = cons(sdp)[ci]
        bmp = cons(sdp)[cip]
        if ci != cip && isapprox(bm, bmp) && isapprox(get(rhs(sdp), ci, 0.0), get(rhs(sdp), cip, 0.0))
            push!(pairs, (ci, cip))
        end
    end
    pairs
end

function remove_equal_constraints!(sdp::SparseSDP)
    eqs = get_equal_constraints(sdp)
    for (ci, _) in eqs
         delete!(cons(sdp), ci)
         delete!(rhs(sdp), ci)
     end
end

function find_removable_row(A::Matrix)
    m = size(A, 1)
    n = size(A, 2)
    r = rank(A)
    for i = 1:m
        B = copy(A)
        B[i, :] = 0.0
        if rank(B) == r
            return i
        end
    end
end

function remove_linear_dependence!(sdp::SparseSDP; debug::Bool=false, dorref=false)
    constraintindices = union(keys(rhs(sdp)), keys(cons(sdp)))
    m = length(constraintindices)

    elindices = Dict()
    n = 1
    debug && println("building dictionaries")
    for bm in values(cons(sdp))
        for (bi, b) in blocks(bm)
            for ((ri, ci), _) in b
                if !haskey(elindices, (bi,ri,ci))
                    elindices[(bi, ri, ci)] = n
                    n += 1
                end
            end
        end
    end

    debug && println("building matrix")
    A = zeros(n, m)
    for j = 1:m
        A[n, j] = get(rhs(sdp), constraintindices[j], 0.0)
        for (bi, b) in cons(sdp)[constraintindices[j]]
            for ((ri, ci), v) in b
                A[elindices[(bi, ri, ci)], j] = v
            end
        end
    end

#        v = []
#        for j = 1:m
#            w = vcat(v, [j])
#            if rank(A[:, w]) == length(w)
#                v = w
#            else
#                ci = constraintindices[j]
#                delete!(rhs(sdp), ci)
#                delete!(cons(sdp), ci)
#            end
#        end
#    println(size(A))
#    U,S,V = svd(A)
#    println(S)
#    debug && println("removing constraints")
#    if dorref
#        debug && println("rref")
#        @time B = rref(A)
#        i = 1
#        counter=0
#        for j = 1:m
#            if isapprox(B[i, j], 0.0, atol=10e-8)
#                counter+=1
#                ci = constraintindices[j]
#                delete!(rhs(sdp), ci)
#                delete!(cons(sdp), ci)
#            else
#                i += 1
#            end
#        end
#        println(counter)
#    else
    println(size(A))
    _, r, p = qr(A, Val{true})

    for i = size(r, 1):-1:1
        println(r[i,i])
        if abs(r[i, i]) > 1e-14
            break
        end
        ci = constraintindices[p[i]]
        delete!(rhs(sdp), ci)
        delete!(cons(sdp), ci)
    end
#        i = 1
#        counter = 0
#        for j = 1:size(r, 2)
#            println(i, " " ,  j)
#            if abs(r[i,j]) >= 10e-10
#                i += 1
#            else
#                println(r[i,j])
#                counter += 1
#                ci = constraintindices[p[j]]
#                delete!(rhs(sdp), ci)
#                delete!(cons(sdp), ci)
#            end
#            if i > size(r, 1)
#                break
#            end
#        end
#        println(counter)
end


text:
using SuperResolution
using Test

@testset "segmentSR.jl" begin
    @test fieldnames(SuperResolution.SegmentSR)==(:select, :idx, :ground_truth, :loc, :frame, :mol_ID, :MU, :K, :N, :bounds, :all_ground_truth, :all_loc, :all_frame, :all_mol_ID)
    # #my_f(x,y)=2x+3y
    # @test my_f(2,1)==7
    # @test my_f(2,3)==13
    # # Write your tests here.
end

max_stars_repo_name: PGSch/SuperResolution

max_issues_repo_name: PGSch/SuperResolution

max_forks_repo_name: PGSch/SuperResolution


text:
export IntRangeSet

import Base: push!, in, show, foreach, collect, union, union!, intersect

type IntRange{T<:Integer}
    lv::T
    rv::T
    lp::Nullable{IntRange{T}}
    rp::Nullable{IntRange{T}}
    lc::Nullable{IntRange{T}}
    rc::Nullable{IntRange{T}}
    IntRange(lv, rv=lv, lp=nothing, rp=nothing, lc=nothing, rc=nothing) = new(lv, rv, lp, rp, lc, rc)
end

type IntRangeSet{T<:Integer} # <: Base.AbstractSet{T}
    balanced::Bool
    root::Nullable{IntRange{T}}
    cache::Nullable{IntRange{T}}
    IntRangeSet() = new(true, nothing, nothing)
end

# NOTE: tree.cache only used and updated in `push!` methods

function push!{T}(tree::IntRangeSet{T}, x::T)::IntRangeSet{T}
    if tree.root.isnull
        node = IntRange{T}(x)
        tree.root = node
        tree.cache = node
    elseif between_parents(x, tree.cache.value)
        push!(tree, tree.cache.value, x)
    else
        push!(tree, tree.root.value, x)
    end

    tree
end

function push!{T}(tree::IntRangeSet{T}, x::UnitRange{T})::IntRangeSet{T}
    x.start <= x.stop || return tree
    if tree.root.isnull
        node = IntRange{T}(x.start, x.stop)
        tree.root = node
        tree.cache = node
    elseif between_parents(x, tree.cache.value)
        push!(tree, tree.cache.value, x)
    else
        push!(tree, tree.root.value, x)
    end

    tree
end

function push!{T}(tree::IntRangeSet{T}, node::IntRange{T}, x::T)::IntRangeSet{T}
    if x == node.lv - 1
        extend_left!(tree, node, x)
        tree.cache = node
    elseif x == node.rv + 1
        extend_right!(tree, node, x)
        tree.cache = node
    elseif x < node.lv
        if node.lc.isnull
            x = IntRange{T}(x, x, node.lp, node)
            node.lc = x
            tree.cache = x
            tree.balanced = false
        else
            push!(tree, node.lc.value, x)
        end
    elseif x > node.rv
        if node.rc.isnull
            x = IntRange{T}(x, x, node, node.rp)
            node.rc = x
            tree.cache = x
            tree.balanced = false
        else
            push!(tree, node.rc.value, x)
        end
    end

    tree
end

"""
six position relations:
+-------------------+
     3333333333
    2222    5555
1111     44     6666
       ******
+-------------------+
1: insert to left child
2: extend left
3: extend both side
4: just ingnore
5: extend right
6: insert to right child
"""
function push!{T}(tree::IntRangeSet{T}, node::IntRange{T}, x::UnitRange{T})::IntRangeSet{T}
    if x.start < node.lv
        if x.stop < node.lv - 1
            if node.lc.isnull
                x = IntRange{T}(x.start, x.stop, node.lp, node)
                node.lc = x
                tree.cache = x
                tree.balanced = false
            else
                push!(tree, node.lc.value, x)
            end
        else
            extend_left!(tree, node, x.start)
            if x.stop > node.rv
                extend_right!(tree, node, x.stop)
            end
            tree.cache = node
        end
    elseif x.stop > node.rv
        if x.start <= node.rv + 1
            extend_right!(tree, node, x.stop)
            tree.cache = node
        else
            if node.rc.isnull
                x = IntRange{T}(x.start, x.stop, node, node.rp)
                node.rc = x
                tree.cache = x
                tree.balanced = false
            else
                push!(tree, node.rc.value, x)
            end
        end
    end

    tree
end

function extend_left!{T}(tree::IntRangeSet{T}, node::IntRange{T}, x::T)::IntRange{T}
    if !node.lc.isnull
        rmlc = right_most(tree, node.lc.value) # here can be optimized: stop at the first node that rv + 1 >= x
        if x <= rmlc.rv + 1
            fuse_rp!(tree, rmlc)
        else
            node.lv = x
        end
    else
        node.lv = x
    end

    x < node.lv ? extend_left!(tree, node, x) : node
end

function extend_right!{T}(tree::IntRangeSet{T}, node::IntRange{T}, x::T)::IntRange{T}
    if !node.rc.isnull
        lmrc = left_most(tree, node.rc.value)
        if x >= lmrc.lv - 1
            fuse_lp!(tree, lmrc)
        else
            node.rv = x
        end
    else
        node.rv = x
    end

    x > node.rv ? extend_right!(tree, node, x) : node
end

# left subtree will be lost, as they are covered by the fused node
function fuse_lp!{T}(tree::IntRangeSet{T}, node::IntRange{T})::IntRange{T}
    # 1. hoist right child to the position of node
    if node.lp.value.rc.value == node # linked to lp directly
        node.lp.value.rc = node.rc
    else
        node.rp.value.lc = node.rc
    end
    # 2. inherit lp
    p = node.rc
    while !p.isnull
        p.value.lp = node.lp
        p = p.value.lc
    end
    # 3. adjust lp range
    node.lp.value.rv = node.rv
    # 4. track balanced property
    tree.balanced = false
    node.lp.value
end

function fuse_rp!{T}(tree::IntRangeSet{T}, node::IntRange{T})::IntRange{T}
    if node.rp.value.lc.value == node
        node.rp.value.lc = node.lc
    else
        node.lp.value.rc = node.lc
    end
    p = node.lc
    while !p.isnull
        p.value.rp = node.rp
        p = p.value.rc
    end
    node.rp.value.lv = node.lv
    tree.balanced = false
    node.rp.value
end

function swap_lc!{T}(tree::IntRangeSet{T}, node::IntRange{T})::IntRange{T}
    # 1. preserve a pointer to the right subtree
    temp = node.lc.value.rc
    # 2. hoist left child to the position of node
    if !node.lp.isnull && node.lp.value.rc.value == node
        node.lp.value.rc = node.lc
    elseif !node.rp.isnull
        node.rp.value.lc = node.lc
    else # root
        tree.root = node.lc
    end
    # 3. inherit
    node.lc.value.rc = node
    node.lc.value.rp = node.rp
    # 4. setup node itself
    node.lp = node.lc
    node.lc = temp
    # 5. return the node in the origin position
    node.lp.value
end

function swap_rc!{T}(tree::IntRangeSet{T}, node::IntRange{T})::IntRange{T}
    temp = node.rc.value.lc
    if !node.rp.isnull && node.rp.value.lc.value == node
        node.rp.value.lc = node.rc
    elseif !node.lp.isnull
        node.lp.value.rc = node.rc
    else
        tree.root = node.rc
    end
    node.rc.value.lc = node
    node.rc.value.lp = node.lp
    node.rp = node.rc
    node.rc = temp
    node.rp.value
end

function rebalance!{T}(tree::IntRangeSet{T})::Int
    depth = tree.root.isnull ? 0 : rebalance!(tree, tree.root.value)
    tree.balanced = true
    depth
end

# TODO: this is O(n) operation. Maybe we should store balance factor in each node?
function rebalance!{T}(tree::IntRangeSet{T}, node::IntRange{T})::Int
    ld = node.lc.isnull ? 0 : rebalance!(tree, node.lc.value)
    rd = node.rc.isnull ? 0 : rebalance!(tree, node.rc.value)
    if ld - rd > 1
        swap_lc!(tree, node)
        ld
    elseif rd - ld > 1
        swap_rc!(tree, node)
        rd
    else # already balance
        max(ld, rd) + 1
    end
end

"apply f to each UnitRange in IntRangeSet, order is guaranteed"
function foreach{T}(f::Function, tree::IntRangeSet{T})::Void
    traverse(tree) do node
        f(node.lv:node.rv)
    end
end

"get a list of UnitRange that is in IntRangeSet, order is guaranteed"
function collect{T}(tree::IntRangeSet{T})::Vector{UnitRange{T}}
    list = UnitRange{T}[]
    foreach(x->push!(list, x), tree)
    list
end

function union{T}(t1::IntRangeSet{T}, t2::IntRangeSet{T})::IntRangeSet{T}
    tree = IntRangeSet{T}()
    union!(tree, t1)
    union!(tree, t2)
    tree
end

function union{T}(ts::IntRangeSet{T}...)::IntRangeSet{T}
    tree = IntRangeSet{T}()
    for t in ts
        union!(tree, t)
    end
    tree
end

function union!{T}(t1::IntRangeSet{T}, t2)::IntRangeSet{T}
    t1.balanced || rebalance!(t1)
    foreach(x->push!(t1, x), t2)
    t1
end

function intersect{T}(t1::IntRangeSet{T}, t2::IntRangeSet{T})::IntRangeSet{T}
    tree = IntRangeSet{T}()
    a, b = collect(t1), collect(t2)
    i, j = 1, 1

    while true
        if i > length(a) || j > length(b)
            break
        end

        start = max(a[i].start, b[j].start)

        if a[i].stop < b[j].stop
            push!(tree, start:a[i].stop)
            i += 1
        else
            push!(tree, start:b[j].stop)
            j += 1
        end
    end

    tree
end

# TODO: intersect all in one pass
function intersect{T}(ts::IntRangeSet{T}...)::IntRangeSet{T}
    reduce(intersect, ts)
end

function in{T}(x::T, tree::IntRangeSet{T})::Bool
    tree.balanced || rebalance!(tree)
    tree.root.isnull ? false : in(tree, tree.root.value, x)
end

function in{T}(tree::IntRangeSet{T}, node::IntRange{T}, x::T)::Bool
    if x < node.lv
        node.lc.isnull ? false : in(tree, node.lc.value, x)
    elseif x > node.rv
        node.rc.isnull ? false : in(tree, node.rc.value, x)
    else
        true
    end
end

function between_parents{T}(x::T, node::IntRange{T})::Bool
    !node.lp.isnull && x <= node.lp.value.rv + 1 && return false
    !node.rp.isnull && x >= node.rp.value.lv - 1 && return false
    true
end

function between_parents{T}(x::UnitRange{T}, node::IntRange{T})::Bool
    !node.lp.isnull && x.start <= node.lp.value.rv + 1 && return false
    !node.rp.isnull && x.stop  >= node.rp.value.lv - 1 && return false
    true
end

function traverse{T}(f::Function, tree::IntRangeSet{T})::Void
    traverse(f, tree, tree.root)
end

function traverse{T}(f::Function, tree::IntRangeSet{T}, node::Nullable{IntRange{T}})::Void
    if !node.isnull
        node = node.value
        traverse(f, tree, node.lc)
        f(node)
        traverse(f, tree, node.rc)
    end
end

function left_most{T}(tree::IntRangeSet{T}, node::IntRange{T})::IntRange{T}
    while !node.lc.isnull
        node = node.lc.value
    end
    node
end

function right_most{T}(tree::IntRangeSet{T}, node::IntRange{T})::IntRange{T}
    while !node.rc.isnull
        node = node.rc.value
    end
    node
end

function show{T}(io::IO, tree::IntRangeSet{T})::Void
    println(io, "IntRangeSets{$T}:")
    if tree.root.isnull
        println(io, "  (empty)")
    else
        traverse(tree) do node
            println(io, "  ", node.lv, ':', node.rv)
        end
    end
end


text:
function clear!(market::AbstractMarket)
    sell_book = market.ask_limit_orders
    buy_book = market.bid_limit_orders

    trades = Array{Trade, 1}()

    return clear!(buy_book, sell_book, market=market)
end

"Clear limit order books"
function clear!(bid_book::Array{BidLimitOrder, 1}, ask_book::Array{AskLimitOrder, 1}; market)

    trades = Array{Trade, 1}()

    # Clear empty orders before the clearing
    # ie. if an investor has set quantity
    # of an order to zero, it will be cleaned now
    clear_empty!(ask_book)
    clear_empty!(bid_book)

    # The clearing price is defined here
    # in case of call market. Nothing
    # should be returned if the price
    # is defined per order pair basis
    clearing_price = get_trade_price(market)
    while maxprice(bid_book) >= minprice(ask_book)

        order_ask = get_best(ask_book)
        order_bid = get_best(bid_book)

        # trade_price is clearing_price if defined,
        # else calculated per order pairs
        trade_price = (
            ~isnothing(clearing_price) ?
            clearing_price
            : get_trade_price(order_bid, order_ask, market=market)
        )

        trade = trade!(
            order_bid, order_ask,
            price=trade_price,
            from=market.currency,
            to=market.traded_asset
        )

        clear_empty!(ask_book)
        clear_empty!(bid_book)

        push!(trades, trade)
    end
    if ~isempty(trades)
        market.last_price = trades[end].price
    end
    return trades
end

function get_trade_price(market::AbstractMarket)
    # Must be redefined by market type. If not defined, nothing returned
    return nothing
end