text:
@testset "BAM" begin
    bamdir = path_of_format("BAM")

    @testset "AuxData" begin
        auxdata = BAM.AuxData(UInt8[])
        @test isempty(auxdata)

        buf = IOBuffer()
        write(buf, "NM", UInt8('s'), Int16(1))
        auxdata = BAM.AuxData(take!(buf))
        @test length(auxdata) == 1
        @test auxdata["NM"] === Int16(1)
        @test collect(auxdata) == ["NM" => Int16(1)]

        buf = IOBuffer()
        write(buf, "AS", UInt8('c'), Int8(-18))
        write(buf, "NM", UInt8('s'), Int16(1))
        write(buf, "XA", UInt8('f'), Float32(3.14))
        write(buf, "XB", UInt8('Z'), "some text\0")
        write(buf, "XC", UInt8('B'), UInt8('i'), Int32(3), Int32[10, -5, 8])
        auxdata = BAM.AuxData(take!(buf))
        @test length(auxdata) == 5
        @test auxdata["AS"] === Int8(-18)
        @test auxdata["NM"] === Int16(1)
        @test auxdata["XA"] === Float32(3.14)
        @test auxdata["XB"] == "some text"
        @test auxdata["XC"] == Int32[10, -5, 8]
        @test convert(Dict{String,Any}, auxdata) == Dict(
            "AS" => Int8(-18),
            "NM" => Int16(1),
            "XA" => Float32(3.14),
            "XB" => "some text",
            "XC" => Int32[10, -5, 8])
    end

    @testset "Record" begin
        record = BAM.Record()
        @test !isfilled(record)
        @test repr(record) == "XAM.BAM.Record: <not filled>"
        @test_throws ArgumentError BAM.flag(record)
    end

    @testset "Reader" begin
        reader = open(BAM.Reader, joinpath(bamdir, "ce#1.bam"))
        @test isa(reader, BAM.Reader)
        @test eltype(reader) === BAM.Record
        @test startswith(repr(reader), "XAM.BAM.Reader{IOStream}:")

        # header
        h = header(reader)
        @test isa(h, SAM.Header)

        # first record
        record = BAM.Record()
        read!(reader, record)
        @test BAM.ismapped(record)
        @test BAM.isprimary(record)
        @test ! BAM.ispositivestrand(record)
        @test BAM.refname(record) == "CHROMOSOME_I"
        @test BAM.refid(record) === 1
        @test BAM.hasnextrefid(record)
        @test BAM.nextrefid(record) === 0
        @test BAM.hasposition(record) === hasleftposition(record) === true
        @test BAM.position(record) === leftposition(record) === 2
        @test BAM.hasnextposition(record)
        @test BAM.nextposition(record) === 0
        @test rightposition(record) == 102
        @test BAM.hastempname(record) === hasseqname(record) === true
        @test BAM.tempname(record) == seqname(record) == "SRR065390.14978392"
        @test BAM.hassequence(record) === hassequence(record) === true
        @test BAM.sequence(record) == sequence(record) == dna"""
        CCTAGCCCTAACCCTAACCCTAACCCTAGCCTAAGCCTAAGCCTAAGCCT
        AAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAA
        """
        @test BAM.seqlength(record) === 100
        @test BAM.hasquality(record)
        @test eltype(BAM.quality(record)) == UInt8
        @test BAM.quality(record) == [Int(x) - 33 for x in "#############################@B?8B?BA@@DDBCDDCBC@CDCDCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"]
        @test BAM.flag(record) === UInt16(16)
        @test BAM.cigar(record) == "27M1D73M"
        @test BAM.alignment(record) == Alignment([
            AlignmentAnchor(  0,   1, OP_START),
            AlignmentAnchor( 27,  28, OP_MATCH),
            AlignmentAnchor( 27,  29, OP_DELETE),
            AlignmentAnchor(100, 102, OP_MATCH)])
        @test record["XG"] == 1
        @test record["XM"] == 5
        @test record["XN"] == 0
        @test record["XO"] == 1
        @test record["AS"] == -18
        @test record["XS"] == -18
        @test record["YT"] == "UU"
        @test keys(record) == ["XG","XM","XN","XO","AS","XS","YT"]
        @test values(record) == [1, 5, 0, 1, -18, -18, "UU"]
        @test eof(reader)
        close(reader)

        # Test conversion from byte array to record
        dsize = BAM.data_size(record)
        array = Vector{UInt8}(undef, BAM.FIXED_FIELDS_BYTES + dsize)
        GC.@preserve array record begin
            ptr = Ptr{UInt8}(pointer_from_objref(record))
            unsafe_copyto!(pointer(array), ptr, BAM.FIXED_FIELDS_BYTES)
            unsafe_copyto!(array, BAM.FIXED_FIELDS_BYTES + 1, record.data, 1, dsize)
        end
        new_record = convert(BAM.Record, array)
        @test record.l_read_name == new_record.l_read_name
        @test record.mapq == new_record.mapq
        @test record.bin == new_record.bin
        @test record.block_size == new_record.block_size
        @test record.flag == new_record.flag
        @test record.n_cigar_op == new_record.n_cigar_op
        @test record.l_seq == new_record.l_seq
        @test record.next_refid == new_record.next_refid
        @test record.next_pos == new_record.next_pos
        @test record.refid == new_record.refid
        @test record.pos == new_record.pos
        @test record.tlen == new_record.tlen
        @test record.data == new_record.data

        # rightposition (also implicitly alignlength)
        records = collect(open(BAM.Reader, joinpath(bamdir, "ce#5b.bam")))
        @test BAM.rightposition(records[6]) == rightposition(records[6]) == 83

        # iterator
        @test length(collect(open(BAM.Reader, joinpath(bamdir, "ce#1.bam")))) == 1
        @test length(collect(open(BAM.Reader, joinpath(bamdir, "ce#2.bam")))) == 2

        # IOStream
        @test length(collect(BAM.Reader(open(joinpath(bamdir, "ce#1.bam"))))) == 1
        @test length(collect(BAM.Reader(open(joinpath(bamdir, "ce#2.bam"))))) == 2
    end

    @testset "Read long CIGARs" begin
        function get_cigar_lens(rec::BAM.Record)
            cigar_ops, cigar_n = BAM.cigar_rle(rec)
            field_ops, field_n = BAM.cigar_rle(rec, false)
            cigar_l = length(cigar_ops)
            field_l = length(field_ops)
            return cigar_l, field_l
        end

        function check_cigar_vs_field(rec::BAM.Record)
            cigar = BAM.cigar(rec)
            field = BAM.cigar(rec, false)
            cigar_l, field_l = get_cigar_lens(rec)
            return cigar != field && cigar_l != field_l
        end

        function check_cigar_lens(rec::BAM.Record, field_len, cigar_len)
            cigar_l, field_l = get_cigar_lens(rec)
            return cigar_l == cigar_len && field_l == field_len
        end

        reader = open(BAM.Reader, joinpath(bamdir, "cigar-64k.bam"))
        rec = BAM.Record()
        read!(reader, rec)
        @test !check_cigar_vs_field(rec)
        read!(reader, rec)
        @test check_cigar_vs_field(rec)
        @test check_cigar_lens(rec, 2, 72091)
    end

    function compare_records(xs, ys)
        if length(xs) != length(ys)
            return false
        end
        for (a, b) in zip(xs, ys)
            if !(
		        a.block_size    == b.block_size &&
		        a.refid         == b.refid &&
		        a.pos           == b.pos &&
		        a.l_read_name   == b.l_read_name &&
		        a.mapq          == b.mapq &&
		        a.bin           == b.bin &&
		        a.n_cigar_op    == b.n_cigar_op &&
		        a.flag          == b.flag &&
		        a.l_seq         == b.l_seq &&
		        a.next_refid    == b.next_refid &&
		        a.next_pos      == b.next_pos &&
		        a.tlen          == b.tlen &&
                a.data[1:BAM.data_size(a)] == b.data[1:BAM.data_size(b)])
                return false
            end
        end
        return true
    end

    @testset "Round trip" begin
        for specimen in list_valid_specimens("BAM")
            filepath = joinpath(bamdir, filename(specimen))
            mktemp() do path, _
                # copy
                if hastags(specimen) && in("bai", tags(specimen))
                    reader = open(BAM.Reader, filepath, index=filepath * ".bai")
                else
                    reader = open(BAM.Reader, filepath)
                end

                header_original = header(reader)

                writer = BAM.Writer(BGZFStream(path, "w"), BAM.header(reader, fillSQ=isempty(findall(header(reader), "SQ"))))

                records = BAM.Record[]
                for record in reader
                    push!(records, record)
                    write(writer, record)
                end
                close(reader)
                close(writer)

                reader = open(BAM.Reader, path)

                @test header(reader) == header_original
                @test compare_records(collect(reader), records)

                close(reader)

            end
        end
    end

    @testset "In-Place-Reading Pattern" begin

        file_bam = joinpath(bamdir, "ce#5b.bam")

        records = open(collect, BAM.Reader, file_bam)

        reader = open(BAM.Reader, file_bam)
        record = BAM.Record()
        i = 0
        while !eof(reader)
            empty!(record) # Reset the record.
            read!(reader, record)

            i = i + 1
            @test records[i] == record
        end

        close(reader)

    end

    @testset "Random access" begin
        filepath = joinpath(bamdir, "GSE25840_GSM424320_GM06985_gencode_spliced.head.bam")
        reader = open(BAM.Reader, filepath, index=filepath * ".bai")

        @test isa(eachoverlap(reader, "chr1", 1:100), BAM.OverlapIterator)
        @test isa(eachoverlap(reader, GenomicFeatures.Interval("chr1", 1, 100)), BAM.OverlapIterator)

        # expected values are counted using samtools
        for (refname, interval, expected) in [
                ("chr1", 1_000:10000,      21),
                ("chr1", 8_000:10000,      20),
                ("chr1", 766_000:800_000, 142),
                ("chr1", 786_000:800_000, 1),
                ("chr1", 796_000:800_000, 0)]
            intsect = eachoverlap(reader, refname, interval)
            @test eltype(intsect) == BAM.Record
            @test count(_ -> true, intsect) == expected
            # check that the intersection iterator is stateless
            @test count(_ -> true, intsect) == expected
        end

        # randomized tests
        for n in 1:50
            refindex = 1
            refname = "chr1"
            range = randrange(1:1_000_000)
            seekstart(reader)
            # linear scan
            expected = filter(collect(reader)) do record
                BAM.compare_intervals(record, (refindex, range)) == 0
            end
            # indexed scan
            actual = collect(eachoverlap(reader, refname, range))
            @test compare_records(actual, expected)
        end
        close(reader)

        filepath = joinpath(bamdir, "R_12h_D06.uniq.q40.bam")
        reader = open(BAM.Reader, filepath, index=filepath * ".bai")
        @test isempty(collect(eachoverlap(reader, "chr19", 5823708:5846478)))
        close(reader)
    end
end


text:
# METHOD TEMPLATES
# basic method template: "This method has priority over the more general T:<Union{....} template"
# set_wfm(obj::INSTR{:Agilent33500B},func="sin",ch="1") = write(obj,"SOURce$ch:FUNC $func")
# more general method template, superseded by basic method template
# set_wfm(obj::T;func="sin",ch="1") where {T<:Union{F332x1,F335x2}} = write(obj,"SOURce$ch:FUNC $func")

## TRIGGER
"""
# Trigger source
`set_trig_sour(obj::T;ch="1",sour="BUS") where {T <: Union{DG1000Z}}`
- sour = {"INTernal", "EXTernal", "BUS"}
"""
set_trig_sour(obj::T;ch="1",sour="BUS") where {T <: Union{DG1000Z}} = write(obj, "TRIG$ch:SOUR $sour")

"""
# Trigger delay (del in seconds)
`set_trig_del(obj::T;ch="1",del=0) where {T <: Union{DG1000Z}}`
- del = delay in s
"""
set_trig_del(obj::T;ch="1",del=0) where {T <: Union{DG1000Z}} = write(obj, "TRIG$ch:DEL $del")

"""
# Trigger slope
`set_trig_slope(obj::T;ch="1",slp="positive") where {T <: Union{DG1000Z}}`
- slp = {"POSitive", "NEGative"}
"""
set_trig_slope(obj::T;ch="1",slp="positive") where {T <: Union{DG1000Z}} = write(obj, "TRIG$ch:SLOP $slp")

"""
# Output trigger state
"""
#set_out_trig_stat(obj::T;st="normal") where {T <: Union{DG1000Z}} = write(obj, "OUTP:TRIG $st")

"""
# Output trigger slope
"""
#set_out_trig_slope(obj::T;slp="positive") where {T <: Union{DG1000Z}} = write(obj, "OUTP:TRIG:SLOP $slp") # RG1022 has SWEep:TRIGger:TRIGOut and BURSt:TRIGger:TRIGOut

"""
# Software trigger
`send_soft_trig(obj::T) where {T <: Union{DG1000Z}}`
"""
send_soft_trig(obj::T) where {T <: Union{DG1000Z}} = write(obj, "*TRG")

## OUTPUT FUNCTION
"""
# Waveform
`set_wfm(obj::T;ch="1",func="sin") where {T <: Union{DG1000Z}}``
- func = {
    "SINusoid", "SQUare", "RAMP", "PULSe", "NOISe", "USER", "HARMonic", "CUSTom", "DC", "KAISER", " ROUNDPM", "SINC", "NEGRAMP", "ATTALT", "AMPALT", "STAIRDN", "STAIRUP", "STAIRUD", "CPULSE", " PPULSE", "NPULSE", "TRAPEZIA", "ROUNDHALF", "ABSSINE", "ABSSINEHALF", "SINETRA", " SINEVER", "EXPRISE", "EXPFALL", "TAN", "COT", "SQRT", "X2DATA", "GAUSS", "HAVERSINE", "LORENTZ ", "DIRICHLET", "GAUSSPULSE", "AIRY", "CARDIAC", "QUAKE", "GAMMA", "VOICE", "TV", "COMBIN", " BANDLIMITED", "STEPRESP", "BUTTERWORTH", "CHEBYSHEV1", "CHEBYSHEV2", "BOXCAR", " BARLETT", "TRIANG", "BLACKMAN", "HAMMING", "HANNING", "DUALTONE", "ACOS", "ACOSH", " ACOTCON", "ACOTPRO", "ACOTHCON", "ACOTHPRO", "ACSCCON", "ACSCPRO", "ACSCHCON", " ACSCHPRO", "ASECCON", "ASECPRO", "ASECH", "ASIN", "ASINH", "ATAN", "ATANH", "BESSELJ", " BESSELY", "CAUCHY", "COSH", "COSINT", "COTHCON", "COTHPRO", "CSCCON", "CSCPRO", " CSCHCON", "CSCHPRO", "CUBIC", "ERF", "ERFC", "ERFCINV", "ERFINV", "LAGUERRE", "LAPLACE", " LEGEND", "LOG", "LOGNORMAL", "MAXWELL", "RAYLEIGH", "RECIPCON", "RECIPPRO", "SECCON", " SECPRO", "SECH", "SINH", "SININT", "TANH", "VERSIERA", "WEIBULL", "BARTHANN", "BLACKMANH", " BOHMANWIN", "CHEBWIN", "FLATTOPWIN", "NUTTALLWIN", "PARZENWIN", "TAYLORWIN", " TUKEYWIN", "CWPUSLE", "LFPULSE", "LFMPULSE", "EOG", "EEG", "EMG", "PULSILOGRAM", "TENS1", " TENS2", "TENS3", "SURGE", "DAMPEDOSC", "SWINGOSC", "RADAR", "THREEAM", "THREEFM", " THREEPM", "THREEPWM", "THREEPFM", "RESSPEED", "MCNOSIE", "PAHCUR", "RIPPLE", " ISO76372TP1", "ISO76372TP2A", "ISO76372TP2B", "ISO76372TP3A", "ISO76372TP3B", " ISO76372TP4", "ISO76372TP5A", "ISO76372TP5B", "ISO167502SP", "ISO167502VR", "SRC", " IGNITION", "NIMHDISCHARGE", "GATEVIBR"
}
"""
set_wfm(obj::T;ch="1",func="sin") where {T <: Union{DG1000Z}} = write(obj, "SOURce$ch:FUNC $func")

"""
# Enable Arb Waveform output
`set_wfm_arb(obj::T;ch="1") where {T <: Union{DG1000Z}}``
"""
set_wfm_arb(obj::T;ch="1") where {T <: Union{DG1000Z}} = write(obj, "SOURce$ch:APPLy:ARBitrary")

"""
# Set Amplitude
`set_amplit(obj::T;ch="1",unit="vpp",volt=0.01) where {T <: Union{DG1000Z}}``
- unit = {"VPP", "VRMS", "DBM"}
- volt = voltage in <unit>
"""
set_amplit(obj::T;ch="1",unit="vpp",volt=0.01) where {T <: Union{DG1000Z}} = write(obj, "SOURce$ch:VOLT:UNIT $unit ; :SOURce$ch:VOLT $volt")

"""
# Set Offset
`set_offs(obj::T;ch="1",volt=0.01) where {T <: Union{DG1000Z}}`
- volt = offset voltage in Volts
"""
set_offs(obj::T;ch="1",volt=0.01) where {T <: Union{DG1000Z}} = write(obj, "SOURce$ch:VOLT:OFFset $volt")

"""
# Set High Level
`set_hilev(obj::T;ch="1",volt=0.01) where {T <: Union{DG1000Z}}`
- volt = High Level in Volts (= Offset + Amplitude/2)
"""
set_hilev(obj::T;ch="1",volt=0.01) where {T <: Union{DG1000Z}} = write(obj, "SOURce$ch:VOLT:HIGH $volt")
"""
# Set Low Level
`set_lolev(obj::T;ch="1",volt=0.01) where {T <: Union{DG1000Z}}`
- volt = Low Level in Volts (= Offset - Amplitude/2)
"""
set_lolev(obj::T;ch="1",volt=0.01) where {T <: Union{DG1000Z}} = write(obj, "SOURce$ch:VOLT:LOW $volt")

"""
# Set Frequency
`set_freq(obj::T;ch="1",fr=1000) where {T <: Union{DG1000Z}}`
- fr = Frequency in Hertz
"""
set_freq(obj::T;ch="1",fr=1000) where {T <: Union{DG1000Z}} = write(obj, "SOURce$ch:FREQ $fr")

"""
# Set Period
`set_per(obj::T;ch="1",per=0.001) where {T <: Union{DG1000Z}}`
- per = Period in Seconds
"""
set_per(obj::T;ch="1",per=0.001) where {T <: Union{DG1000Z}} = write(obj, "SOURce$ch:PER $per")

"""
# Set Phase
`set_phase(obj::T;ch="1",phase=90) where {T <: Union{DG1000Z}}`
- phase = Phase in Deg
"""
set_phase(obj::T;ch="1",phase=90) where {T <: Union{DG1000Z}} = write(obj, "SOURce$ch:PHAS $phase")

"""
# Set Symmetry
`set_sym(obj::T;ch="1",func="ramp",sym=50) where {T <: Union{DG1000Z}}`
- func = always RAMP for DG1000Z
- sym = 0% to 100%
"""
set_sym(obj::T;ch="1",func="ramp",sym=50) where {T <: Union{DG1000Z}} = write(obj, "SOURce$ch:FUNC:RAMP:SYMM $sym")

# Both Edge times
set_both_edg(obj::T;ch="1",func="pulse",dt=10e-9) where {T <: Union{DG1000Z}} = write(obj, "SOURce$ch:FUNC:$func:TRAN $dt")
# Lead Edge time
set_lead_edg(obj::T;ch="1",func="pulse",dt=10e-9) where {T <: Union{DG1000Z}} = write(obj, "SOURce$ch:FUNC:$func:TRAN:LEAD $dt")
# Trail Edge time
set_trail_edg(obj::T;ch="1",func="pulse",dt=10e-9) where {T <: Union{DG1000Z}} = write(obj, "SOURce$ch:FUNC:$func:TRAN:TRA $dt")
# Width
set_width(obj::T;ch="1",func="pulse",wd=1e-4) where {T <: Union{DG1000Z}} = write(obj, "SOURce$ch:FUNC:$func:WIDTh $wd")
#
## OUTPUT
# Output Load
set_load(obj::T;ch="1",load="INF") where {T <: Union{DG1000Z}} = write(obj, "OUTP$ch:LOAD $load")
# Output Polarity
set_pol(obj::T;ch="1",pol="normal") where {T <: Union{DG1000Z}} = write(obj, "OUTP$ch:POL $pol")

"""
# Set Output on/off
`set_outp(obj::T;ch="1",st="off") where {T <: Union{DG1000Z}}`
- st = {"on", "off"}
"""
set_outp(obj::T;ch="1",st="off") where {T <: Union{DG1000Z}} = write(obj, "OUTP$ch $st")

## SYNC
# Sync state
set_sync_stat(obj::T;ch="1",st="on") where {T <: Union{DG1000Z}} = write(obj, "OUTP$ch:SYNC $st")

## ARBITRARY WAVEFORMS
"""
# Set Arb mode
`set_arb_mode(obj::T;ch="1",mode="FREQ") where {T <: Union{DG1000Z}}``
- mode = {"FREQ", "SRATe"}
"""
set_arb_mode(obj::T;ch="1",mode="FREQ") where {T <: Union{DG1000Z}} = write(obj, "SOURce$ch:FUNC:ARB:MODE $mode")

"""
# Set sample rate for Arb mode
`set_sample_rate(obj::T;ch="1",srate=1e6) where {T <: Union{DG1000Z}}`
- srate = Sample rate in samples per second
"""
set_sample_rate(obj::T;ch="1",srate=1e6) where {T <: Union{DG1000Z}} = write(obj, "SOURce$ch:FUNC:ARB:SRAT $srate")

# Clear all arbs from volatile memory
# clear_arbs(obj::T;ch="1") where {T<:Union{DG1000Z}} = write(obj,"SOURce$ch:DATA:VOLATILE:CLEar") # Doesnt seem supported on the Rigol_DG1000Z

"""
# Load Arb data points to volatile memory
`set_arb_data(obj::T;ch="1",values="-1,-1,-0.5,-0.5,0,0.5,0.5,1.0,1.0") where {T <: Union{DG1000Z}}`
- values = -1.0 to +1.0
"""
set_arb_data(obj::T;ch="1",values="-1,-1,-0.5,-0.5,0,0.5,0.5,1.0,1.0") where {T <: Union{DG1000Z}} = write(obj, "SOURce$ch:DATA VOLATILE,$values")

"""
# Set number of data points for Arb mode
`set_arb_points(obj::T;ch="1",points=10) where {T <: Union{DG1000Z}}`
- points = 8 to 16384
"""
set_arb_points(obj::T;ch="1",points=10) where {T <: Union{DG1000Z}} = write(obj, "SOURce$ch:DATA:POIN VOLATILE,$points")

## Sweep

set_sweep_mode(obj::T; ch="1", st="off") where {T <: Union{DG1000Z}} = write(obj, "SOURce$ch:SWE:STAT $st")

"""
- time = 1ms to 500s (in seconds)
"""
set_sweep_time(obj::T; ch="1", time=1) where {T <: Union{DG1000Z}} = write(obj, "SOURce$ch:SWE:TIME $time")

set_sweep_start(obj::T; ch="1", freq::Integer=1000) where {T <: Union{DG1000Z}} = write(obj, ":SOUR$ch:FREQ:STAR $freq")

set_sweep_stop(obj::T; ch="1", freq::Integer=100) where {T <: Union{DG1000Z}} = write(obj, ":SOUR$ch:FREQ:STOP $freq")

set_sweep_center(obj::T; ch="1", freq::Integer=100) where {T <: Union{DG1000Z}} = write(obj, ":SOUR$ch:FREQ:CENT $freq")

set_sweep_span(obj::T; ch="1", span::Integer=100) where {T <: Union{DG1000Z}} = write(obj, ":SOUR$ch:FREQ:SPAN $span")

ext: jl


text:
function op!(
  o::ITensor, ::OpName"Id", ::SiteType"Generic", s1::Index, sn::Index...; eltype=Float64
)
  s = (s1, sn...)
  n = prod(dim.(s))
  t = itensor(Matrix(one(eltype) * I, n, n), prime.(s)..., dag.(s)...)
  return settensor!(o, tensor(t))
end

function op!(o::ITensor, ::OpName"I", st::SiteType"Generic", s::Index...; kwargs...)
  return op!(o, OpName("Id"), st, s...; kwargs...)
end

function op!(o::ITensor, ::OpName"F", st::SiteType"Generic", s::Index; kwargs...)
  return op!(o, OpName("Id"), st, s; kwargs...)
end

function default_random_matrix(eltype::Type, s::Index...)
  n = prod(dim.(s))
  return randn(eltype, n, n)
end

# Haar-random unitary
#
# Reference:
# Section 4.6
# http://math.mit.edu/~edelman/publications/random_matrix_theory.pdf
function op!(
  o::ITensor,
  ::OpName"RandomUnitary",
  ::SiteType"Generic",
  s1::Index,
  sn::Index...;
  eltype=ComplexF64,
  random_matrix=default_random_matrix(eltype, s1, sn...),
)
  s = (s1, sn...)
  Q, _ = NDTensors.qr_positive(random_matrix)
  t = itensor(Q, prime.(s)..., dag.(s)...)
  return settensor!(o, tensor(t))
end

function op!(o::ITensor, ::OpName"randU", st::SiteType"Generic", s::Index...; kwargs...)
  return op!(o, OpName("RandomUnitary"), st, s...; kwargs...)
end


text:
"""
    AbstractStepRangeLen

Supertype for `StepSRangeLen` and `StepMRangeLen`. It's subtypes should behave
identically to `StepRangeLen`.
"""
abstract type AbstractStepRangeLen{T,R,S} <: AbstractRange{T} end

function Base.StepRangeLen{T,R,S}(r::AbstractStepRangeLen) where {T,R,S}
    return StepRangeLen{T,R,S}(convert(R, r.ref), convert(S, r.step), length(r), r.offset)
end

Base.nbitslen(r::AbstractStepRangeLen) = nbitslen(eltype(r), length(r), r.offset)

"""
    StepSRangeLen

A static range `r` where `r[i]` produces values of type `T` (in the second form,
`T` is deduced automatically), parameterized by a `ref`erence value, a `step`,
and the `len`gth. By default `ref` is the starting value `r[1]`, but
alternatively you can supply it as the value of `r[offset]` for some other
index `1 <= offset <= len`. In conjunction with `TwicePrecision` this can be
used to implement ranges that are free of roundoff error.
"""
struct StepSRangeLen{T,Tr,Ts,R,S,L,F} <: AbstractStepRangeLen{T,R,S} end

function StepSRangeLen{T,R,S}(ref::R, step::S, len::Integer, offset::Integer = 1) where {T,R,S}
    len >= 0 || throw(ArgumentError("length cannot be negative, got $len"))
    1 <= offset <= max(1,len) || throw(ArgumentError("StepSRangeLen: offset must be in [1,$len], got $offset"))
    return StepSRangeLen{T,R,S,tp2val(ref),tp2val(step),len,offset}()
end

function StepSRangeLen{T,R1,S1}(ref::R2, step::S2, len::Integer, offset::Integer) where {T,R1,S1,R2,S2}
    return StepSRangeLen{T,R1,S1}(R1(ref), S1(step), len, offset)
end

function (::Type{<:StepSRangeLen{Float64}})(r::AbstractRange)
    return _convertSSRL(StepSRangeLen{Float64,TwicePrecision{Float64},TwicePrecision{Float64}}, r)
end

function Base.getproperty(r::StepSRangeLen, s::Symbol)
    if s === :ref
        return _ref(r)
    elseif s === :step
        return step_hp(r)
    elseif s === :len
        return length(r)
    elseif s === :offset
        return _offset(r)
    else
        error("type $(typeof(r)) has no property $s")
    end
end

# convert TPVal to TwicePrecision
_ref(::StepSRangeLen{T,Tr,Ts,R,S,L,F}) where {T,Tr<:TwicePrecision,Ts,R,S,L,F} = convert(Tr, R)

_offset(::StepSRangeLen{T,Tr,Ts,R,S,L,F}) where {T,Tr,Ts,R,S,L,F} = F
_ref(::StepSRangeLen{T,Tr,Ts,R,S,L,F}) where {T,Tr,Ts,R,S,L,F} = R

"""
    StepMRangeLen

A mutable range `r` where `r[i]` produces values of type `T` (in the second form,
`T` is deduced automatically), parameterized by a `ref`erence value, a `step`,
and the `len`gth. By default `ref` is the starting value `r[1]`, but
alternatively you can supply it as the value of `r[offset]` for some other
index `1 <= offset <= len`. In conjunction with `TwicePrecision` this can be
used to implement ranges that are free of roundoff error.
"""
mutable struct StepMRangeLen{T,R,S} <: AbstractStepRangeLen{T,R,S}
    ref::R
    step::S
    len::Int
    offset::Int

    function StepMRangeLen{T,R,S}(ref::R, step::S, len::Integer, offset::Integer = 1) where {T,R,S}
        len >= 0 || throw(ArgumentError("length cannot be negative, got $len"))
        1 <= offset <= max(1,len) || throw(ArgumentError("StepMRangeLen: offset must be in [1,$len], got $offset"))
        return new(ref, step, len, offset)
    end
end

function (::Type{StepMRangeLen{Float64}})(r::AbstractRange)
    return _convertSMRL(StepMRangeLen{Float64,TwicePrecision{Float64},TwicePrecision{Float64}}, r)
end
function StepMRangeLen{T,R1,S1}(ref::R2, step::S2, len::Integer, offset::Integer) where {T,R1,S1,R2,S2}
    return StepMRangeLen{T,R1,S1}(R1(ref), S1(step), len, offset)
end

_offset(r::StepMRangeLen) = getfield(r, :offset)
_ref(r::StepMRangeLen) = getfield(r, :ref)

function Base.setproperty!(r::StepMRangeLen, s::Symbol, val)
    if s === :ref
        return set_ref!(r, val)
    elseif s === :step
        return set_step!(r, val)
    elseif s === :len
        return set_length!(r, val)
    elseif s === :offset
        return set_offset!(r, val)
    else
        error("type $(typeof(r)) has no property $s")
    end
end

for (F,f) in ((:M,:m), (:S,:s))
    SR = Symbol(:Step, F, :RangeLen)
    frange = Symbol(f, :range)
    CSRL = Symbol(:_convert, :S, F, :RL)
    _CSRL = Symbol(:__convert, :S, F, :RL)
    floatfrange = Symbol(:float, f, :range)
    @eval begin
        $(SR){T,R,S}(r::$(SR){T,R,S}) where {T,R,S} = r
        function $(SR){T,R,S}(r::$(SR)) where {T,R,S}
            return $(SR){T,R,S}(
                convert(R, _ref(r)),
                convert(S, step(r)),
                length(r),
                _offset(r)
               )
        end
        function $(SR){T}(r::Union{StepRangeLen,AbstractStepRangeLen}) where {T}
            return $(SR)(
                convert(T, _ref(r)),
                convert(T, step(r)),
                length(r),
                _offset(r)
               )
        end

        $(SR){T,R,S}(r::AbstractRange) where {T,R,S} = $(SR){T,R,S}(R(first(r)), S(step(r)), length(r))
        $(SR){T,R,S}(r::StepRangeLen) where {T,R,S} = $(SR){T,R,S}(R(r.ref), S(r.step), r.len, r.offset)
        $(SR){T}(r::AbstractRange) where {T} = $(SR)(T(first(r)), T(step(r)), length(r))
        $(SR)(r::AbstractRange) = $(SR){eltype(r)}(r)

        function $(SR)(
            ref::TwicePrecision{T},
            step::TwicePrecision{T},
            len::Integer,
            offset::Integer=1
           ) where {T}
            return $(SR){T,TwicePrecision{T},TwicePrecision{T}}(ref, step, len, offset)
        end

        $(SR){T,R,S}(r::$(SR){T,R,S}) where {T<:AbstractFloat,R<:TwicePrecision,S<:TwicePrecision} = r
        function $(SR)(ref::R, step::S, len::Integer, offset::Integer = 1) where {R,S}
            return $(SR){typeof(ref+0*step),R,S}(ref, step, len, offset)
        end
        function $(SR){T}(ref::R, step::S, len::Integer, offset::Integer = 1) where {T,R,S}
            return $(SR){T,R,S}(ref, step, len, offset)
        end

        function $(SR){T,R,S}(r::$(SR)) where {T<:AbstractFloat,R<:TwicePrecision,S<:TwicePrecision}
            return $(CSRL)($(SR){T,R,S}, r)
        end

        function (::Type{<:$(SR){Float64}})(r::$(SR))
            return $(CSRL)($(SR){Float64,TwicePrecision{Float64},TwicePrecision{Float64}}, r)
        end
        $(SR){T}(r::$(SR)) where {T<:IEEEFloat} = $(CSRL)($(SR){T,Float64,Float64}, r)

        function $(SR){T}(r::AbstractRange) where {T<:IEEEFloat}
            return $(CSRL)($(SR){T,Float64,Float64}, r)
        end

        function $(CSRL)(::Type{<:$(SR){T,R,S}}, r::$(SR){<:Integer}) where {T,R,S}
            return $(SR){T,R,S}(R(_ref(r)), S(step(r)), length(r), _offset(r))
        end

        function $(CSRL)(::Type{<:$(SR){T,R,S}}, r::AbstractRange{<:Integer}) where {T,R,S}
            return $(SR){T,R,S}(R(first(r)), S(step(r)), length(r))
        end

        function $(CSRL)(::Type{<:$(SR){T,R,S}}, r::AbstractRange{U}) where {T,R,S,U}
            # if start and step have a rational approximation in the old type,
            # then we transfer that rational approximation to the new type
            f, s = first(r), step(r)
            start_n, start_d = rat(f)
            step_n, step_d = rat(s)
            if start_d != 0 && step_d != 0 &&
                    U(start_n/start_d) == f && U(step_n/step_d) == s
                den = lcm(start_d, step_d)
                m = maxintfloat(T, Int)
                if den != 0 && abs(f*den) <= m && abs(s*den) <= m &&
                        rem(den, start_d) == 0 && rem(den, step_d) == 0
                    start_n = round(Int, f*den)
                    step_n = round(Int, s*den)
                    return $(floatfrange)(T, start_n, step_n, length(r), den)
                end
            end
            return $(_CSRL)($(SR){T,R,S}, r)
        end

        function $(_CSRL)(::Type{<:$(SR){T,R,S}}, r::$(SR){U}) where {T,R,S,U}
            return $(SR){T,R,S}(R(_ref(r)), S(step(r)), length(r), _offset(r))
        end
        function $(_CSRL)(::Type{<:$(SR){T,R,S}}, r::AbstractRange{U}) where {T,R,S,U}
            return $(SR){T,R,S}(R(first(r)), S(step(r)), length(r))
        end
   end
end


text:
# MLStyle patches

@active GlobalRef(x) begin
    if x isa GlobalRef
        (x.mod, x.name)
    else
        nothing
    end
end

@active Argument(x) begin
    if x isa Argument
        Some(x.n)
    else
        nothing
    end
end

@active SSAValue(x) begin
    if x isa SSAValue
        Some(x.id)
    else
        nothing
    end
end

@active SlotNumber(x) begin
    if x isa SlotNumber
        Some(x.id)
    else
        nothing
    end
end


text:
using PoincareInvariants

export guiding_center_4d_surface_ode,
       guiding_center_4d_surface_iode,
       guiding_center_4d_surface_iode_λ

export guiding_center_4d_ode_poincare_invariant_2nd,
       guiding_center_4d_ode_poincare_invariant_2nd_trapezoidal,
       guiding_center_4d_iode_poincare_invariant_2nd,
       guiding_center_4d_iode_poincare_invariant_2nd_trapezoidal,
       guiding_center_4d_iode_poincare_invariant_2nd_λ

function f_surface(i, j, nx, ny)
    f_surface(i/nx, j/ny)
end

function initial_conditions_surface(nx, ny)
    q₀ = zeros(4, nx*ny)

    for j in 1:ny
        for i in 1:nx
            q₀[:,nx*(j-1)+i] .= f_surface(i, j, nx, ny)
        end
    end

    return q₀
end

guiding_center_4d_surface_ode_init(q₀) = guiding_center_4d_ode(q₀, (μ = μ_surface(),); periodic=false)
guiding_center_4d_surface_iode_init(q₀) = guiding_center_4d_iode(q₀, (μ = μ_surface(),); periodic=false)
guiding_center_4d_surface_iode_λ_init(q₀) = guiding_center_4d_iode_λ(q₀, (μ = μ_surface(),); periodic=false)

function guiding_center_4d_surface_ode(nx, ny)
    guiding_center_4d_surface_ode_init(initial_conditions_surface(nx, ny))
end

function guiding_center_4d_surface_iode(nx, ny)
    guiding_center_4d_surface_iode_init(initial_conditions_surface(nx, ny))
end

function guiding_center_4d_surface_iode_λ(nx, ny)
    guiding_center_4d_surface_iode_λ_init(initial_conditions_surface(nx, ny))
end

function guiding_center_4d_ode_poincare_invariant_2nd(Δt, nx, ny, ntime, nsave, DT=Float64)
    PoincareInvariant2nd(guiding_center_4d_surface_ode_init, f_surface, ω, Δt, 4, nx, ny, ntime, nsave, DT)
end

function guiding_center_4d_ode_poincare_invariant_2nd_trapezoidal(Δt, nx, ny, ntime, nsave, DT=Float64)
    PoincareInvariant2ndTrapezoidal(guiding_center_4d_surface_ode_init, f_surface, ω, Δt, 4, nx, ny, ntime, nsave, DT)
end

function guiding_center_4d_iode_poincare_invariant_2nd(Δt, nx, ny, ntime, nsave, DT=Float64)
    PoincareInvariant2nd(guiding_center_4d_surface_iode_init, f_surface, ω, Δt, 4, nx, ny, ntime, nsave, DT)
end

function guiding_center_4d_iode_poincare_invariant_2nd_trapezoidal(Δt, nx, ny, ntime, nsave, DT=Float64)
    PoincareInvariant2ndTrapezoidal(guiding_center_4d_surface_iode_init, f_surface, ω, (D²ϑ₁, D²ϑ₂, D²ϑ₃, D²ϑ₄), Δt, 4, nx, ny, ntime, nsave)
end

function guiding_center_4d_iode_poincare_invariant_2nd_λ(Δt, nx, ny, ntime, nsave, DT=Float64)
    PoincareInvariant2nd(guiding_center_4d_surface_iode_init, f_surface, ω, (D²ϑ₁, D²ϑ₂, D²ϑ₃, D²ϑ₄), Δt, 4, nx, ny, ntime, nsave, DT)
end


text:
function linesearch_main(img, process_image, dims, iters)
    prefscale   = 1.0
    marginalize = false

    slider_out = nothing
    x1         = nothing
    x2         = nothing
    function redraw_op!(canvas, slider)
        slider_out = deepcopy(slider)
        param = slider*x1 + (1-slider)*x2
        if(isempty(Reactive._messages.data))
            res   = process_image(param)
            res   = clamp.(res, 0, 1)
            copy!(canvas, res)
        end
    end

    block = Gtk.Condition()
    dct = create_ui(img, redraw_op!, block)

    function propose_image(x1_in::Vector, x2_in::Vector)
        push!(dct["gui"]["slider"], 0.5)
        x1 = x1_in
        x2 = x2_in
        Gtk.wait(block)
        slider_out
    end
    prefbo_linesearch(propose_image, iters, dims, prefscale, marginalize;
                      search_budget=100000)
end


text:
using ViZDoom
using Test

# include("basic.jl")


text:
@doc raw"""
    AbstractDecoratorManifold{𝔽} <: AbstractManifold{𝔽}

Declare a manifold to be an abstract decorator.
A manifold which is a subtype of is a __decorated manifold__, i.e. has

* certain additional properties or
* delegates certain properties to other manifolds.

Most prominently, a manifold might be an embedded manifold, i.e. points on a manifold ``\mathcal M``
are represented by (some, maybe not all) points on another manifold ``\mathcal N``.
Depending on the type of embedding, several functions are dedicated to the embedding.
For example if the embedding is isometric, then the [`inner`](@ref) does not have to be
implemented for ``\mathcal M`` but can be automatically implemented by deligation to ``\mathcal N``.

This is modelled by the `AbstractDecoratorManifold` and traits. These are mapped to functions,
which determine the types of transparencies.
"""
abstract type AbstractDecoratorManifold{𝔽} <: AbstractManifold{𝔽} end

"""
    AbstractTrait

An abstract trait type to build a sequence of traits
"""
abstract type AbstractTrait end

"""
    EmptyTrait <: AbstractTrait

A Trait indicating that no feature is present.
"""
struct EmptyTrait <: AbstractTrait end

"""
    IsExplicitDecorator <: AbstractTrait

Specify that a certain type should dispatch per default to its [`decorated_manifold`](@ref).

!!! note
    Any decorator _behind_ this decorator might not have any effect, since the function
    dispatch is moved to its field at this point. Therefore this decorator should always be
    _last_ in the [`TraitList`](@ref).
"""
struct IsExplicitDecorator <: AbstractTrait end

"""
    TraitList <: AbstractTrait

Combine two traits into a combined trait.  Note that this introduces a preceedence.
the first of the traits takes preceedence if a trait is implemented for both functions.

# Constructor

    TraitList(head::AbstractTrait, tail::AbstractTrait)
"""
struct TraitList{T1<:AbstractTrait,T2<:AbstractTrait} <: AbstractTrait
    head::T1
    tail::T2
end

function Base.show(io::IO, t::TraitList)
    return print(io, "TraitList(", t.head, ", ", t.tail, ")")
end

"""
    active_traits(f, args...)

Return the list of traits applicable to the given call of function `f``. This function should be
overloaded for specific function calls.
"""
@inline active_traits(f, args...) = EmptyTrait()

"""
    merge_traits(t1, t2, trest...)

Merge two traits into a nested list of traits. Note that this takes trait preceedence into account,
i.e. `t1` takes preceedence over `t2` is any operations.
It always returns either ab [`EmptyTrait`](@ref) or a [`TraitList`](@ref).

This means that for
* one argument it just returns the trait itself if it is list-like, or wraps the trait in a
    single-element list otherwise,
* two arguments that are list-like, it merges them,
* two arguments of which only the first one is list-like and the second one is not,
    it appends the second argument to the list,
* two arguments of which only the second one is list-like, it prepends the first one to
    the list,
* two arguments of which none is list-like, it creates a two-element list.
* more than two arguments it recursively performs a left-assiciative recursive reduction
    on arguments, that is for example `merge_traits(t1, t2, t3)` is equivalent to
    `merge_traits(merge_traits(t1, t2), t3)`
"""
merge_traits()

@inline merge_traits() = EmptyTrait()
@inline merge_traits(t::EmptyTrait) = t
@inline merge_traits(t::TraitList) = t
@inline merge_traits(t::AbstractTrait) = TraitList(t, EmptyTrait())
@inline merge_traits(t1::EmptyTrait, ::EmptyTrait) = t1
@inline merge_traits(::EmptyTrait, t2::AbstractTrait) = merge_traits(t2)
@inline merge_traits(t1::AbstractTrait, t2::EmptyTrait) = TraitList(t1, t2)
@inline merge_traits(t1::AbstractTrait, t2::TraitList) = TraitList(t1, t2)
@inline merge_traits(::EmptyTrait, t2::TraitList) = t2
@inline function merge_traits(t1::AbstractTrait, t2::AbstractTrait)
    return TraitList(t1, TraitList(t2, EmptyTrait()))
end
@inline merge_traits(t1::TraitList, ::EmptyTrait) = t1
@inline function merge_traits(t1::TraitList, t2::AbstractTrait)
    return TraitList(t1.head, merge_traits(t1.tail, t2))
end
@inline function merge_traits(t1::TraitList, t2::TraitList)
    return TraitList(t1.head, merge_traits(t1.tail, t2))
end
@inline function merge_traits(
    t1::AbstractTrait,
    t2::AbstractTrait,
    t3::AbstractTrait,
    trest::AbstractTrait...,
)
    return merge_traits(merge_traits(t1, t2), t3, trest...)
end

"""
    parent_trait(t::AbstractTrait)

Return the parent trait for trait `t`, that is the more general trait whose behaviour it
inherits as a fallback.
"""
@inline parent_trait(::AbstractTrait) = EmptyTrait()

@inline function trait(f::TF, args...) where {TF}
    bt = active_traits(f, args...)
    return expand_trait(bt)
end

"""
    expand_trait(t::AbstractTrait)

Expand given trait into an ordered [`TraitList`](@ref) list of traits with their parent
traits obtained using [`parent_trait`](@ref).
"""
expand_trait(::AbstractTrait)

@inline expand_trait(e::EmptyTrait) = e
@inline expand_trait(t::AbstractTrait) = merge_traits(t, expand_trait(parent_trait(t)))
@inline function expand_trait(t::TraitList)
    et1 = expand_trait(t.head)
    et2 = expand_trait(t.tail)
    return merge_traits(et1, et2)
end

"""
    next_trait(t::AbstractTrait)

Return the next trait to consider, which by default is no following trait (i.e. [`EmptyTrait`](@ref)).

Expecially for a a [`TraitList`](@ref) this function returns the (remaining) tail of
the remaining traits.
"""
next_trait(::AbstractTrait) = EmptyTrait()

@inline next_trait(t::TraitList) = t.tail

#! format: off
# turn formatting for for the following functions
# due to the if with returns inside (formatter puts a return upfront the if)
function _split_signature(sig::Expr)
    if sig.head == :where
        where_exprs = sig.args[2:end]
        call_expr = sig.args[1]
    elseif sig.head == :call
        where_exprs = []
        call_expr = sig
    else
        error("Incorrect syntax in $sig. Expected a :where or :call expression.")
    end
    fname = call_expr.args[1]
    if isa(call_expr.args[2], Expr) && call_expr.args[2].head == :parameters
        # we have keyword arguments
        callargs = call_expr.args[3:end]
        kwargs_list = call_expr.args[2].args
    else
        callargs = call_expr.args[2:end]
        kwargs_list = []
    end
    argnames = map(callargs) do arg
        if isa(arg, Expr) && arg.head === :kw # default val present
            arg = arg.args[1]
        end
        if isa(arg, Expr) && arg.head === :(::) # typed
            return arg.args[1]
        end
        return arg
    end
    argtypes = map(callargs) do arg
        if isa(arg, Expr) && arg.head === :kw # default val present
            arg = arg.args[1]
        end
        if isa(arg, Expr)
            return arg.args[2]
        else
            return Any
        end
    end

    kwargs_call = map(kwargs_list) do kwarg
        if kwarg.head === :...
            return kwarg
        else
            if isa(kwarg.args[1], Symbol)
                kwargname = kwarg.args[1]
            else
                kwargname = kwarg.args[1].args[1]
            end
            return :($kwargname = $kwargname)
        end
    end

    return (;
        fname = fname,
        where_exprs = where_exprs,
        callargs = callargs,
        kwargs_list = kwargs_list,
        argnames = argnames,
        argtypes = argtypes,
        kwargs_call = kwargs_call,
    )
end
#! format: on

macro invoke_maker(argnum, type, sig)
    parts = ManifoldsBase._split_signature(sig)
    kwargs_list = parts[:kwargs_list]
    callargs = parts[:callargs]
    fname = parts[:fname]
    where_exprs = parts[:where_exprs]
    argnames = parts[:argnames]
    argtypes = parts[:argtypes]
    kwargs_call = parts[:kwargs_call]

    return esc(
        quote
            function ($fname)($(callargs...); $(kwargs_list...)) where {$(where_exprs...)}
                return invoke(
                    $fname,
                    Tuple{
                        $(argtypes[1:(argnum - 1)]...),
                        $type,
                        $(argtypes[(argnum + 1):end]...),
                    },
                    $(argnames...);
                    $(kwargs_call...),
                )
            end
        end,
    )
end

macro trait_function(sig, opts = :())
    parts = ManifoldsBase._split_signature(sig)
    kwargs_list = parts[:kwargs_list]
    callargs = parts[:callargs]
    fname = parts[:fname]
    where_exprs = parts[:where_exprs]
    argnames = parts[:argnames]
    argtypes = parts[:argtypes]
    kwargs_call = parts[:kwargs_call]

    argnametype_exprs = [:(typeof($(argname))) for argname in argnames]

    block = quote
        @inline function ($fname)($(callargs...); $(kwargs_list...)) where {$(where_exprs...)}
            return ($fname)(trait($fname, $(argnames...)), $(argnames...); $(kwargs_call...))
        end
        @inline function ($fname)(
            t::TraitList,
            $(callargs...);
            $(kwargs_list...),
        ) where {$(where_exprs...)}
            return ($fname)(next_trait(t), $(argnames...); $(kwargs_call...))
        end
        @inline function ($fname)(
            t::TraitList{IsExplicitDecorator},
            $(callargs...);
            $(kwargs_list...),
        ) where {$(where_exprs...)}
            arg1 = decorated_manifold($(argnames[1]))
            argt1 = typeof(arg1)
            return invoke(
                $fname,
                Tuple{argt1,$(argnametype_exprs[2:end]...)},
                arg1,
                $(argnames[2:end]...);
                $(kwargs_call...),
            )
        end
    end
    if !(:no_empty in opts.args)
        block = quote
            $block
            @inline function ($fname)(
                ::EmptyTrait,
                $(callargs...);
                $(kwargs_list...),
            ) where {$(where_exprs...)}
                return invoke(
                    $fname,
                    Tuple{supertype($(argtypes[1])),$(argnametype_exprs[2:end]...)},
                    $(argnames...);
                    $(kwargs_call...),
                )
            end
        end
    end
    return esc(block)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ c96551a0-881e-11eb-0d92-61edbb6c86e6
begin
	using Pkg; Pkg.activate("MLJ_env", shared=true)
  using PlutoUI
	using Test

end

# ╔═╡ accc6f84-881d-11eb-3b49-e9cd25a3b6fd
md"""
## Linear Algebra

ref. from book **"Data Science from Scratch"**, Chap 4

$(html"
© Pascal, March 2021
")
"""

# ╔═╡ 136cdd48-8835-11eb-0482-737cf7bd5a3e
PlutoUI.TableOfContents(indent=true, depth=4, aside=true)

# ╔═╡ b22cf860-8834-11eb-02de-d3ad803d1db6
html"""

"""

# ╔═╡ b1e26cfa-8834-11eb-16f5-85efb64c6a4c
md"""
### Vectors
Points in some finite dimensional space. A vector has a direction and a magnitude.$(html"
")
Vectors are useful to represent data (which lie in some multi-finite-dimensional space).
We will want to perform operations on those vectors, so let defines and implements those operations. $(html"
")
We will be using the `Vector` type from `Julia`.
"""

# ╔═╡ 9b2a7694-881e-11eb-325c-f3edd9178641
begin
	import Base: +, -

	function +(v₁::Vector{T}, v₂::Vector{T})::Vector{T} where T <: Number
		"""
		Vector Additon is component-wise
		"""
		@assert length(v₁) == length(v₂)
		v₁ .+ v₂
	end

	function +(vc::Vector{Vector{T}})::Vector{T} where T <: Number
		"""
		Sum all corresponding components of a collection of vector to return
		one vector
		"""
		@assert length(vc) > 0 "expect the collection to have at leat one vector element"
		n = length(vc[1])
		@assert all(v -> length(v) == n, vc) "Expect all vectors to be of he same length"
		v = zeros(T, n)
		for v₀ ∈ vc
			v .+= v₀
		end
		v
	end

	function -(v₁::Vector{T}, v₂::Vector{T})::Vector{T} where T <: Number
		"""
		Vector Additon is component-wise
		"""
		@assert length(v₁) == length(v₂)
		v₁ .- v₂
	end
end

# ╔═╡ 21a74512-881f-11eb-31e6-974bcaf59141
begin
	@test +([1, 2, 3], [4, 5, 6]) == [5, 7, 9]
	@test [1, 2, 3] + [4, 5, 6] == [5, 7, 9]

	vᵣ = [1.0, -1.0, π, ℯ]
	@test vᵣ + [-1.0, 1.0, -π, -ℯ] == zeros(eltype(vᵣ), length(vᵣ))
end

# ╔═╡ 11ff42a0-8828-11eb-3ba5-dbd12b5ffdda
md"""
Sometimes we want to sum a collection of vector component wise. that is create a new vector whose first component is the sum of evry first comp[onent of thew collection and  so for for the other components.

Here the collection is itsef a vector in `Julia`'s term.
cf. celle above for implementation
"""

# ╔═╡ 7dc3365c-8828-11eb-3c57-69dae41619cd
begin
	@test +([[1, 2], [3, 4], [5, 6], [7, 8]]) == [16, 20]
	@test +([[1, 3, 2] for _ ∈ 1:3]) == [3, 9, 6]
end

# ╔═╡ e9b7d526-8829-11eb-34d1-e317fdc7a1d0
md"""
`Multiplication` by a scalar, `Mean` and `dot` product.
"""

# ╔═╡ 3a7741c2-882a-11eb-2ffd-65b97e75a0b2
begin
	import Base.*

	*(x::T, v::Vector{T}) where T <: Number = v .* x
	*(v::Vector{T}, x::T) where T <: Number = v .* x

	function μ(vc::Vector{Vector{T}})::Vector{AbstractFloat} where T <: Number
		"""Element-wise mean"""
		@assert length(vc) > 0
		n = length(vc[1])
		@assert all(v -> length(v) == n, vc)
		reduce(+, vc) / n
 	end

	function dot(v₁::Vector{T}, v₂::Vector{T})::T where T <: Number
		@assert length(v₁) == length(v₂)
		sum(v₁ .* v₂)
	end
end

# ╔═╡ b3f61d4e-881f-11eb-2dcb-4dd0fd7fba06
begin
	@test -([1, 2, 3], [4, 5, 6]) == [-3, -3, -3]
	@test [1, 2, 3] - [4, 5, 6] == [-3, -3, -3]

	# vᵣ = [1.0, -1.0, π, ℯ]
	@test vᵣ - [-1.0, 1.0, -π, -ℯ] == [2.0, -2.0, 2π, 2ℯ]
end

# ╔═╡ 19f5b176-882b-11eb-3713-73c85f599be1
begin
	@test *(2, [1, 2, 3]) == [2, 4, 6]
	@test *(2, [1, 2, 3]) == *([1, 2, 3], 2)
	@test 2 * [1, 2, 3] == [1, 2, 3] * 2  ==  [2, 4, 6]
end

# ╔═╡ 50444c92-882b-11eb-2893-bf66b63ee595
begin
	μ([[1, 2, 1, 2], [3, 4, 3, 4], [5, 6, 5, 6]]) == Float64[3., 4., 3., 4.]
end

# ╔═╡ d358476c-882d-11eb-26bd-c167cabd7b28
begin
	@test dot([1, 2, 3], [4, 5, 6]) == 32
end

# ╔═╡ d33d417e-882d-11eb-0aae-a1d10512920f
function sum_of_square(v::Vector{T})::T where T <: Number
	"""
	return vᵢ * vᵢ ∀ i ∈ 1:length(v)
	"""
	dot(v, v)
end

# ╔═╡ d3202a62-882d-11eb-2e24-6d874dce773d
begin
	@test sum_of_square([1, 2, 3]) == 14
	@test sum_of_square([4, 5, 6]) == 77
end

# ╔═╡ 9db8e630-882e-11eb-07e5-0d5036f7e8e2
md"""
Let's now implement L₂ norm and euclidean distance.
"""

# ╔═╡ ba1d63da-882e-11eb-1736-37a246585172
function norm(v::Vector{T})::AbstractFloat where T <: Number
	√ sum(v .* v)
end

# ╔═╡ 9da1fcfc-882e-11eb-021f-1f89149817a9
begin
	@test norm([3, 4]) == 5
	@test norm([1, 0, 2]) ≈ √5
end

# ╔═╡ 9d8634ae-882e-11eb-15ad-75b06c847379
function distance(v₁::Vector{T}, v₂::Vector{T})::AbstractFloat where T <: Number
	"""
	distance(v₁, v₂) ≡ √(Σ (v₁ᵢ - v₂ᵢ)²)
	"""
	@assert length(v₁) == length(v₂)
	√ sum((v₁ .- v₂).^2)  # == norm(v₁ .- v₂)
end

# ╔═╡ 9d69ec22-882e-11eb-0637-d7e29f7b9f37
begin
	@test distance([2, 6, 7, 7, 5, 13, 14, 17, 11, 8],
		[3, 5, 5, 3, 7, 12, 13, 19, 22, 7]) ≈ 12.409673645990857
end

# ╔═╡ 0539916a-8835-11eb-3407-41f17efc2ee4
html"""

"""

# ╔═╡ 73e7de7a-8830-11eb-05ef-b93f6c774587
md"""
### Matrices

A two-dimensional collection of numbers.
"""

# ╔═╡ 7f051430-8830-11eb-1e43-797df6b195ed
function shape(m::Matrix{T})::Tuple{Integer, Integer} where T <: Number
	size(m)
end

# ╔═╡ 7ee4e642-8830-11eb-3d15-813c4dbb6420
@test shape([1 2 3; 4 5 6]) == (2, 3) # 2 rows, 3 columns

# ╔═╡ 7ec91688-8830-11eb-18a0-d5a67e267f9b
function get_row(m::Matrix{T}, r::Integer)::Vector{T}  where T <: Number
	@assert 1 ≤ r ≤ size(m)[1]
	view(m, r, :)  ## no copy!
end

# ╔═╡ 7e9afe92-8830-11eb-3772-150255e27923
begin
	@test get_row([1 2 3; 4 5 6], 1) == [1, 2, 3]
	@test get_row([1 2 3; 4 5 6], 2) == [4, 5, 6]
	@test_throws AssertionError get_row([1 2 3; 4 5 6], 0)
end

# ╔═╡ 648af952-8831-11eb-340b-5bc7ef7aaa73
function get_column(m::Matrix{T}, c::Integer)::Vector{T}  where T <: Number
	@assert 1 ≤ c ≤ size(m)[2]
	view(m, :, c)  ## no copy!
end

# ╔═╡ 6470ad90-8831-11eb-0413-5d016bf66536
begin
	@test get_column([1 2 3; 4 5 6], 1) == [1, 4]
	@test get_column([1 2 3; 4 5 6], 2) == [2, 5]
	@test get_column([1 2 3; 4 5 6], 3) == [3, 6]
	@test_throws AssertionError get_column([1 2 3; 4 5 6], 0)
end

# ╔═╡ 6458987c-8831-11eb-2b89-efac8fabf803
md"""
We will want to be able to create a new matrix given its shape and a function generator.
"""

# ╔═╡ c832331c-8831-11eb-0b7c-9f978c82f890
function make_matrix(nrows::Integer, ncols::Integer, fn::Function;
		DT::DataType=Float64)::Matrix
	@assert 1 ≤ nrows && 1 ≤ ncols

	m = zeros(DT, (nrows, ncols))
	for c ∈ 1:ncols, r ∈ 1:nrows
		m[r, c] = fn(r, c)
	end
	m
end

# ╔═╡ be10ce56-8832-11eb-0742-d90fa7f6f669
function identity_matrix(n::Integer; DT::DataType=Float64)::Matrix
	"""Returns the n × n identity matrix"""
	make_matrix(n, n, (i, j) -> i == j ? one(DT) : zero(DT); DT)
end

# ╔═╡ c81a12dc-8831-11eb-27e9-07132c6310a1
@test identity_matrix(5) == [1. 0. 0. 0. 0.; 0. 1. 0. 0. 0.; 0. 0. 1. 0. 0.; 0. 0. 0. 1. 0.; 0. 0. 0. 0. 1.]

# ╔═╡ 3c69e83a-8835-11eb-218e-0bb108ad7b9c
html"""

Mar 2021, Corto Inc
"""

# ╔═╡ Cell order:
# ╟─accc6f84-881d-11eb-3b49-e9cd25a3b6fd
# ╟─136cdd48-8835-11eb-0482-737cf7bd5a3e
# ╟─b22cf860-8834-11eb-02de-d3ad803d1db6
# ╟─b1e26cfa-8834-11eb-16f5-85efb64c6a4c
# ╠═c96551a0-881e-11eb-0d92-61edbb6c86e6
# ╠═9b2a7694-881e-11eb-325c-f3edd9178641
# ╠═21a74512-881f-11eb-31e6-974bcaf59141
# ╠═b3f61d4e-881f-11eb-2dcb-4dd0fd7fba06
# ╟─11ff42a0-8828-11eb-3ba5-dbd12b5ffdda
# ╠═7dc3365c-8828-11eb-3c57-69dae41619cd
# ╟─e9b7d526-8829-11eb-34d1-e317fdc7a1d0
# ╠═3a7741c2-882a-11eb-2ffd-65b97e75a0b2
# ╠═19f5b176-882b-11eb-3713-73c85f599be1
# ╠═50444c92-882b-11eb-2893-bf66b63ee595
# ╠═d358476c-882d-11eb-26bd-c167cabd7b28
# ╠═d33d417e-882d-11eb-0aae-a1d10512920f
# ╠═d3202a62-882d-11eb-2e24-6d874dce773d
# ╟─9db8e630-882e-11eb-07e5-0d5036f7e8e2
# ╠═ba1d63da-882e-11eb-1736-37a246585172
# ╠═9da1fcfc-882e-11eb-021f-1f89149817a9
# ╠═9d8634ae-882e-11eb-15ad-75b06c847379
# ╠═9d69ec22-882e-11eb-0637-d7e29f7b9f37
# ╟─0539916a-8835-11eb-3407-41f17efc2ee4
# ╟─73e7de7a-8830-11eb-05ef-b93f6c774587
# ╠═7f051430-8830-11eb-1e43-797df6b195ed
# ╠═7ee4e642-8830-11eb-3d15-813c4dbb6420
# ╠═7ec91688-8830-11eb-18a0-d5a67e267f9b
# ╠═7e9afe92-8830-11eb-3772-150255e27923
# ╠═648af952-8831-11eb-340b-5bc7ef7aaa73
# ╠═6470ad90-8831-11eb-0413-5d016bf66536
# ╟─6458987c-8831-11eb-2b89-efac8fabf803
# ╠═c832331c-8831-11eb-0b7c-9f978c82f890
# ╠═be10ce56-8832-11eb-0742-d90fa7f6f669
# ╠═c81a12dc-8831-11eb-27e9-07132c6310a1
# ╟─3c69e83a-8835-11eb-218e-0bb108ad7b9c

max_stars_repo_path: Data Science From Scratch/04_Linear_Algebra.jl

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks


text:
using GNSSTools
using Statistics
using ProgressMeter
using PyPlot
pygui(true)

h_parms = [1.9e-6, 2e-6, 9.7e-8, 9.1e-11, 9.9e-15]./16.3676e6^2  # TCXO
CN0 = 45.  # dB⋅Hz
f_d = 0
iterations = 100
T = 1e-3
f_s₀ = 2.048e6
f_ss = Array(1:100)
t_length = T

# L1 C/A
sigtype = define_l1ca_code_type(; sig_freq=1*L1_freq)
# generatesignal!(signal)
l1ca_SNRs = zeros(length(f_ss))
p = Progress(length(f_ss)*iterations, 1, "Processing L1 C/A...")
for i in 1:length(f_ss)
    SNR = 0
    f_s = f_ss[i] * f_s₀
    replica = definesignal(sigtype, f_s, T; f_d=f_d, allocate_noise_vectors=false,
                           skip_noise_generation=true)
    generatereplica!(replica)
    signal = definesignal(sigtype, f_s, t_length;
                          prn=1, code_start_idx=1000, f_d=f_d, fd_rate=0,
                          CN0=CN0, include_phase_noise=false,
                          include_thermal_noise=true, include_adc=false,
                          include_carrier=true, receiver_h_parms=h_parms,
                          skip_noise_generation=true)
    for j in 1:iterations
        definesignal!(signal; new_thermal_noise=true, new_phase_noise=false)
        generatesignal!(signal)
        corr_result = fft_correlate(signal.data[1:replica.sample_num], replica.data)
        PS = maximum(abs2.(corr_result))
        PN = ((sum(abs2.(corr_result))) - PS) / (replica.sample_num - 1)
        SNR += PS/PN
        next!(p)
    end
    l1ca_SNRs[i] = 10log10(SNR/iterations)
end

fig = figure()
plot(f_ss.*f_s₀./1e6, l1ca_SNRs, label="L1 C/A")
# plot(Ts, l5_SNRs, label="L5")
xlabel("Sampling Rate (MHz)")
ylabel("Acquired SNR (dB)")
legend()
xscale("log")
savefig("figures/ch2_snr_vs_sampling_rate.pdf", dpi=300)


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

using Markdown
using InteractiveUtils

# ╔═╡ b56a7126-46ba-11eb-2fdf-7f93842fc501
md"## _BroadCasting_
In this lesson we will talk about one of the most handy, if not useful, features of Julia: **array broadcasting**. Furthermore we will deal with some important differences between Julia and other programming languages (like Python or Matlab) which make extensive use of operations on arrays and with arrays."

# ╔═╡ 1fa5aed2-46bb-11eb-3d2c-37d60ffb1698
md" ## Working with Functions and Arrays
**Introduction**

Many languages (such as Python with numpy and Matlab) make extensive use of optimised C or Fortran routines under the hood to perform fast mathematical operations, as such the user is encouraged to write vectorised code, so that these routines can perform faster for loops, as a big part of the CPU time is spent on calling the underlying compiled routine, and not computing the actual result. More or less what happens is that the user writes vectorised code which communicates to C code that can run fast for loops and the result is then returned to the user in the form of an array (or matrix).

In Julia, since for loops are already as fast as they can be (close to the speed of C) there is no need to write vectorised code, as the interpreter will directly compile your code in optimised machine code which will run as fast as possible on your machine. In Julia nothing happens under the hood (beside the compilation of the functions) and almost all the functions in Julia are written in Julia, just like all the functions in C are written in C.

**Operations with Arrays**

Julia by default deals with operations on arrays and matrices as one would do in mathematics.

Let’s start with an example: from a mathematical point of view, we don’t know how to compute the sin of an array, as the sine function is defined only on single (dimensionless) values. At the same time the exp can work both on single values and matrices (as the exponential of a matrix has a well-defined geometrical meaning). For the same reason, you cannot multiply two arrays together, unless their size is matching correctly (i.e. one array is a row array and the other one is a column array) and in this case the multiplication of two arrays becomes the well-defined geometrical product of two arrays (which can be a scalar or a matrix, depending on the order of the multiplication):
"

# ╔═╡ bcbc3b46-46bb-11eb-1532-c561f25a9943
a = [1,2,3]

# ╔═╡ c6571362-46bb-11eb-2a5f-51ffab5cd8ab
b = [4,5,6]

# ╔═╡ ce409166-46bb-11eb-1ea6-d52e8922a64f
a*b

# ╔═╡ d39e1610-46bb-11eb-2e3b-71b169116ea1
c = [4 5 6]

# ╔═╡ e41dbb1c-46bb-11eb-07d3-ebd94d38fe07
a*c

# ╔═╡ 0db631e8-46bc-11eb-3684-85af439fcfb1
c*a

# ╔═╡ 1168dac0-46bc-11eb-1898-696d8b3565a5
d = reshape([1,2,3,4,5,6,7,8,9],3,3)

# ╔═╡ 1e6c7ef4-46bc-11eb-16a6-7bf5d86af426
d*a

# ╔═╡ 23dc16ae-46bc-11eb-3145-25a33518b1b2
md" This makes perfectly sense from a mathematical point of view and operators behave how we would mathematically expect. Nonetheless, in programming it is often useful to write operations which work on an element by element basis, and for this reason broadcasting comes to our help.
"

# ╔═╡ 3cf1e0e2-46bc-11eb-280d-01cb5a9b5043
md" ## Broadcasting
In Julia, with broadcasting we indicate the action of mapping a function or an operation (which are the same in Julia) over an array or a matrix element by element.

The broadcasting notation for operators consists of adding a dot . before the operator (for example .*)

Considering the example we get:
"

# ╔═╡ 6aef8c38-46bc-11eb-1329-8ba6c0556433
a .* c

# ╔═╡ 6d66571a-46bc-11eb-0821-2fcb49952c0e
c .* a

# ╔═╡ 7844aa58-46bc-11eb-0ebf-051cad17621a
a .* d

# ╔═╡ 7b144928-46bc-11eb-1c61-b76df7fe294d
md"
Notice that when we broadcast the multiplication with a matrix and an array, the array gets multiplied “in the same direction” as it is written, in the sense that if a vector is a column it gets applied column by column etc.

We can use the broadcasting notation also to map a function over an n-dimensional array. There is no speed gain in doing so, as it will be exactly equivalent to writing a for loop, but its conciseness may be useful sometimes. So the core idea in Julia is to write functions that take single values and use broadcasting when needed, unless the functions must explicitly work on arrays (for example to compute the mean of a series of values, perform matrix operations, vector multiplications, etc).

To broadcast a function over an array it is sufficient to put a dot before the brackets `.()`
"

# ╔═╡ 8ade4caa-46bc-11eb-1ed3-938b46b63dfa
sin.(a)

# ╔═╡ a590b9a2-46bc-11eb-04b6-4703bf84cbdd
md"
It is also possible to map a function over several values using the map function, but there are no real advantages in doing so and the broadcasting syntax is often more flexible.

As we will see in the future, in Julia it is really easy to write parallel code using multi-processing and multi-threading, and it is particularly simple to write such parallel structures when we encounter a for loop or a map function (which will become a parallel map through pmap), thus it is convenient to get used to writing your code in such paradigm.

If you are new to programming or you come from a compiled language (like C++), you should find it easy to think in terms of for loops and you should think of broadcasting as a fast and concise way to compute the value of a function over several input values.
"

# ╔═╡ add08b24-46bc-11eb-0637-43c62f9bfc00
md" ## Conclusion
We have learned what is broadcasting and how it can be used to perform element by element operations between vectors and how to map a function over an array using the concise broadcasting syntax.

We have also pointed out how Julia does not gain in performance by using “vectorised” notations (while other languages do) because the core Julia operations are implemented directly in the Julia language and thus there is no need to call compiled routines coded in other languages under the hood to offload all the heavy work.

"

# ╔═╡ Cell order:
# ╟─b56a7126-46ba-11eb-2fdf-7f93842fc501
# ╟─1fa5aed2-46bb-11eb-3d2c-37d60ffb1698
# ╠═bcbc3b46-46bb-11eb-1532-c561f25a9943
# ╠═c6571362-46bb-11eb-2a5f-51ffab5cd8ab
# ╠═ce409166-46bb-11eb-1ea6-d52e8922a64f
# ╠═d39e1610-46bb-11eb-2e3b-71b169116ea1
# ╠═e41dbb1c-46bb-11eb-07d3-ebd94d38fe07
# ╠═0db631e8-46bc-11eb-3684-85af439fcfb1
# ╠═1168dac0-46bc-11eb-1898-696d8b3565a5
# ╠═1e6c7ef4-46bc-11eb-16a6-7bf5d86af426
# ╟─23dc16ae-46bc-11eb-3145-25a33518b1b2
# ╟─3cf1e0e2-46bc-11eb-280d-01cb5a9b5043
# ╠═6aef8c38-46bc-11eb-1329-8ba6c0556433
# ╠═6d66571a-46bc-11eb-0821-2fcb49952c0e
# ╠═7844aa58-46bc-11eb-0ebf-051cad17621a
# ╟─7b144928-46bc-11eb-1c61-b76df7fe294d
# ╠═8ade4caa-46bc-11eb-1ed3-938b46b63dfa
# ╟─a590b9a2-46bc-11eb-04b6-4703bf84cbdd
# ╟─add08b24-46bc-11eb-0637-43c62f9bfc00


text:
import MathOptInterface
const MOI = MathOptInterface

mutable struct VariableInfo
    lower_bound::Float64  # May be -Inf even if has_lower_bound == true
    has_lower_bound::Bool # Implies lower_bound == Inf
    upper_bound::Float64  # May be Inf even if has_upper_bound == true
    has_upper_bound::Bool # Implies upper_bound == Inf
    is_fixed::Bool        # Implies lower_bound == upper_bound and !has_lower_bound and !has_upper_bound.
    start::Float64
end
# The default start value is zero.
VariableInfo() = VariableInfo(-Inf, false, Inf, false, false, 0.0)

mutable struct Optimizer <: MOI.AbstractOptimizer
    inner::Union{IpoptProblem,Nothing}
    variable_info::Vector{VariableInfo}
    nlp_data::MOI.NLPBlockData
    sense::MOI.OptimizationSense
    objective::Union{MOI.SingleVariable,MOI.ScalarAffineFunction{Float64},MOI.ScalarQuadraticFunction{Float64},Nothing}
    linear_le_constraints::Vector{Tuple{MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}}}
    linear_ge_constraints::Vector{Tuple{MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}}}
    linear_eq_constraints::Vector{Tuple{MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}}}
    quadratic_le_constraints::Vector{Tuple{MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}}}
    quadratic_ge_constraints::Vector{Tuple{MOI.ScalarQuadraticFunction{Float64}, MOI.GreaterThan{Float64}}}
    quadratic_eq_constraints::Vector{Tuple{MOI.ScalarQuadraticFunction{Float64}, MOI.EqualTo{Float64}}}
    options
end

struct EmptyNLPEvaluator <: MOI.AbstractNLPEvaluator end
MOI.features_available(::EmptyNLPEvaluator) = [:Grad, :Jac, :Hess]
MOI.initialize(::EmptyNLPEvaluator, features) = nothing
MOI.eval_objective(::EmptyNLPEvaluator, x) = NaN
function MOI.eval_constraint(::EmptyNLPEvaluator, g, x)
    @assert length(g) == 0
    return
end
MOI.eval_objective_gradient(::EmptyNLPEvaluator, g, x) = nothing
MOI.jacobian_structure(::EmptyNLPEvaluator) = Tuple{Int64,Int64}[]
MOI.hessian_lagrangian_structure(::EmptyNLPEvaluator) = Tuple{Int64,Int64}[]
function MOI.eval_constraint_jacobian(::EmptyNLPEvaluator, J, x)
    @assert length(J) == 0
    return
end
function MOI.eval_hessian_lagrangian(::EmptyNLPEvaluator, H, x, σ, μ)
    @assert length(H) == 0
    return
end

empty_nlp_data() = MOI.NLPBlockData([], EmptyNLPEvaluator(), false)

Optimizer(;options...) = Optimizer(nothing, [], empty_nlp_data(), MOI.FEASIBILITY_SENSE, nothing, [], [], [], [], [], [], options)

MOI.supports(::Optimizer, ::MOI.NLPBlock) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{MOI.SingleVariable}) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.SingleVariable}, ::Type{MOI.LessThan{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.SingleVariable}, ::Type{MOI.GreaterThan{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.SingleVariable}, ::Type{MOI.EqualTo{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.ScalarAffineFunction{Float64}}, ::Type{MOI.LessThan{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.ScalarAffineFunction{Float64}}, ::Type{MOI.GreaterThan{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.ScalarAffineFunction{Float64}}, ::Type{MOI.EqualTo{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.ScalarQuadraticFunction{Float64}}, ::Type{MOI.LessThan{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.ScalarQuadraticFunction{Float64}}, ::Type{MOI.GreaterThan{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.ScalarQuadraticFunction{Float64}}, ::Type{MOI.EqualTo{Float64}}) = true

function MOI.copy_to(model::Optimizer, src::MOI.ModelLike; copy_names = false)
    return MOI.Utilities.default_copy_to(model, src, copy_names)
end

MOI.get(model::Optimizer, ::MOI.NumberOfVariables) = length(model.variable_info)

function MOI.get(model::Optimizer, ::MOI.ListOfVariableIndices)
    return [MOI.VariableIndex(i) for i in 1:length(model.variable_info)]
end

function MOI.set(model::Optimizer, ::MOI.ObjectiveSense,
                 sense::MOI.OptimizationSense)
    model.sense = sense
    return
end

function MOI.empty!(model::Optimizer)
    model.inner = nothing
    empty!(model.variable_info)
    model.nlp_data = empty_nlp_data()
    model.sense = MOI.FEASIBILITY_SENSE
    model.objective = nothing
    empty!(model.linear_le_constraints)
    empty!(model.linear_ge_constraints)
    empty!(model.linear_eq_constraints)
    empty!(model.quadratic_le_constraints)
    empty!(model.quadratic_ge_constraints)
    empty!(model.quadratic_eq_constraints)
end

function MOI.is_empty(model::Optimizer)
    return isempty(model.variable_info) &&
           model.nlp_data.evaluator isa EmptyNLPEvaluator &&
           model.sense == MOI.FEASIBILITY_SENSE &&
           isempty(model.linear_le_constraints) &&
           isempty(model.linear_ge_constraints) &&
           isempty(model.linear_eq_constraints) &&
           isempty(model.quadratic_le_constraints) &&
           isempty(model.quadratic_ge_constraints) &&
           isempty(model.quadratic_eq_constraints)
end

function MOI.add_variable(model::Optimizer)
    push!(model.variable_info, VariableInfo())
    return MOI.VariableIndex(length(model.variable_info))
end
function MOI.add_variables(model::Optimizer, n::Int)
    return [MOI.add_variable(model) for i in 1:n]
end

function check_inbounds(model::Optimizer, vi::MOI.VariableIndex)
    num_variables = length(model.variable_info)
    if !(1 <= vi.value <= num_variables)
        error("Invalid variable index $vi. ($num_variables variables in the model.)")
    end
end

function check_inbounds(model::Optimizer, var::MOI.SingleVariable)
    return check_inbounds(model, var.variable)
end

function check_inbounds(model::Optimizer, aff::MOI.ScalarAffineFunction)
    for term in aff.terms
        check_inbounds(model, term.variable_index)
    end
end

function check_inbounds(model::Optimizer, quad::MOI.ScalarQuadraticFunction)
    for term in quad.affine_terms
        check_inbounds(model, term.variable_index)
    end
    for term in quad.quadratic_terms
        check_inbounds(model, term.variable_index_1)
        check_inbounds(model, term.variable_index_2)
    end
end

function has_upper_bound(model::Optimizer, vi::MOI.VariableIndex)
    return model.variable_info[vi.value].has_upper_bound
end

function has_lower_bound(model::Optimizer, vi::MOI.VariableIndex)
    return model.variable_info[vi.value].has_lower_bound
end

function is_fixed(model::Optimizer, vi::MOI.VariableIndex)
    return model.variable_info[vi.value].is_fixed
end

function MOI.add_constraint(model::Optimizer, v::MOI.SingleVariable, lt::MOI.LessThan{Float64})
    vi = v.variable
    check_inbounds(model, vi)
    if isnan(lt.upper)
        error("Invalid upper bound value $(lt.upper).")
    end
    if has_upper_bound(model, vi)
        error("Upper bound on variable $vi already exists.")
    end
    if is_fixed(model, vi)
        error("Variable $vi is fixed. Cannot also set upper bound.")
    end
    model.variable_info[vi.value].upper_bound = lt.upper
    model.variable_info[vi.value].has_upper_bound = true
    return MOI.ConstraintIndex{MOI.SingleVariable, MOI.LessThan{Float64}}(vi.value)
end

function MOI.add_constraint(model::Optimizer, v::MOI.SingleVariable, gt::MOI.GreaterThan{Float64})
    vi = v.variable
    check_inbounds(model, vi)
    if isnan(gt.lower)
        error("Invalid lower bound value $(gt.lower).")
    end
    if has_lower_bound(model, vi)
        error("Lower bound on variable $vi already exists.")
    end
    if is_fixed(model, vi)
        error("Variable $vi is fixed. Cannot also set lower bound.")
    end
    model.variable_info[vi.value].lower_bound = gt.lower
    model.variable_info[vi.value].has_lower_bound = true
    return MOI.ConstraintIndex{MOI.SingleVariable, MOI.GreaterThan{Float64}}(vi.value)
end

function MOI.add_constraint(model::Optimizer, v::MOI.SingleVariable, eq::MOI.EqualTo{Float64})
    vi = v.variable
    check_inbounds(model, vi)
    if isnan(eq.value)
        error("Invalid fixed value $(gt.lower).")
    end
    if has_lower_bound(model, vi)
        error("Variable $vi has a lower bound. Cannot be fixed.")
    end
    if has_upper_bound(model, vi)
        error("Variable $vi has an upper bound. Cannot be fixed.")
    end
    if is_fixed(model, vi)
        error("Variable $vi is already fixed.")
    end
    model.variable_info[vi.value].lower_bound = eq.value
    model.variable_info[vi.value].upper_bound = eq.value
    model.variable_info[vi.value].is_fixed = true
    return MOI.ConstraintIndex{MOI.SingleVariable, MOI.EqualTo{Float64}}(vi.value)
end

macro define_add_constraint(function_type, set_type, array_name)
    quote
        function MOI.add_constraint(model::Optimizer, func::$function_type, set::$set_type)
            check_inbounds(model, func)
            push!(model.$(array_name), (func, set))
            return MOI.ConstraintIndex{$function_type, $set_type}(length(model.$(array_name)))
        end
    end
end

@define_add_constraint(MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64},
                       linear_le_constraints)
@define_add_constraint(MOI.ScalarAffineFunction{Float64},
                       MOI.GreaterThan{Float64}, linear_ge_constraints)
@define_add_constraint(MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64},
                       linear_eq_constraints)
@define_add_constraint(MOI.ScalarQuadraticFunction{Float64},
                       MOI.LessThan{Float64}, quadratic_le_constraints)
@define_add_constraint(MOI.ScalarQuadraticFunction{Float64},
                       MOI.GreaterThan{Float64}, quadratic_ge_constraints)
@define_add_constraint(MOI.ScalarQuadraticFunction{Float64},
                       MOI.EqualTo{Float64}, quadratic_eq_constraints)

function MOI.supports(::Optimizer, ::MOI.VariablePrimalStart,
                      ::Type{MOI.VariableIndex})
    return true
end
function MOI.set(model::Optimizer, ::MOI.VariablePrimalStart,
                 vi::MOI.VariableIndex, value::Real)
    check_inbounds(model, vi)
    model.variable_info[vi.value].start = value
    return
end

function MOI.set(model::Optimizer, ::MOI.NLPBlock, nlp_data::MOI.NLPBlockData)
    model.nlp_data = nlp_data
    return
end

function MOI.set(model::Optimizer, ::MOI.ObjectiveFunction,
                 func::Union{MOI.SingleVariable, MOI.ScalarAffineFunction,
                             MOI.ScalarQuadraticFunction})
    check_inbounds(model, func)
    model.objective = func
    return
end

# In setting up the data for Ipopt, we order the constraints as follows:
# - linear_le_constraints
# - linear_ge_constraints
# - linear_eq_constraints
# - quadratic_le_constraints
# - quadratic_ge_constraints
# - quadratic_eq_constraints
# - nonlinear constraints from nlp_data

linear_le_offset(model::Optimizer) = 0
linear_ge_offset(model::Optimizer) = length(model.linear_le_constraints)
linear_eq_offset(model::Optimizer) = linear_ge_offset(model) + length(model.linear_ge_constraints)
quadratic_le_offset(model::Optimizer) = linear_eq_offset(model) + length(model.linear_eq_constraints)
quadratic_ge_offset(model::Optimizer) = quadratic_le_offset(model) + length(model.quadratic_le_constraints)
quadratic_eq_offset(model::Optimizer) = quadratic_ge_offset(model) + length(model.quadratic_ge_constraints)
nlp_constraint_offset(model::Optimizer) = quadratic_eq_offset(model) + length(model.quadratic_eq_constraints)

# Convenience functions used only in optimize!

function append_to_jacobian_sparsity!(jacobian_sparsity, aff::MOI.ScalarAffineFunction, row)
    for term in aff.terms
        push!(jacobian_sparsity, (row, term.variable_index.value))
    end
end

function append_to_jacobian_sparsity!(jacobian_sparsity, quad::MOI.ScalarQuadraticFunction, row)
    for term in quad.affine_terms
        push!(jacobian_sparsity, (row, term.variable_index.value))
    end
    for term in quad.quadratic_terms
        row_idx = term.variable_index_1
        col_idx = term.variable_index_2
        if row_idx == col_idx
            push!(jacobian_sparsity, (row, row_idx.value))
        else
            push!(jacobian_sparsity, (row, row_idx.value))
            push!(jacobian_sparsity, (row, col_idx.value))
        end
    end
end

# Refers to local variables in jacobian_structure() below.
macro append_to_jacobian_sparsity(array_name)
    escrow = esc(:row)
    quote
        for (func, set) in $(esc(array_name))
            append_to_jacobian_sparsity!($(esc(:jacobian_sparsity)), func, $escrow)
            $escrow += 1
        end
    end
end

function jacobian_structure(model::Optimizer)
    num_nlp_constraints = length(model.nlp_data.constraint_bounds)
    if num_nlp_constraints > 0
        nlp_jacobian_sparsity = MOI.jacobian_structure(model.nlp_data.evaluator)
    else
        nlp_jacobian_sparsity = []
    end

    jacobian_sparsity = Tuple{Int64,Int64}[]
    row = 1
    @append_to_jacobian_sparsity model.linear_le_constraints
    @append_to_jacobian_sparsity model.linear_ge_constraints
    @append_to_jacobian_sparsity model.linear_eq_constraints
    @append_to_jacobian_sparsity model.quadratic_le_constraints
    @append_to_jacobian_sparsity model.quadratic_ge_constraints
    @append_to_jacobian_sparsity model.quadratic_eq_constraints
    for (nlp_row, column) in nlp_jacobian_sparsity
        push!(jacobian_sparsity, (nlp_row + row - 1, column))
    end
    return jacobian_sparsity
end

append_to_hessian_sparsity!(hessian_sparsity, ::Union{MOI.SingleVariable,MOI.ScalarAffineFunction}) = nothing

function append_to_hessian_sparsity!(hessian_sparsity, quad::MOI.ScalarQuadraticFunction)
    for term in quad.quadratic_terms
        push!(hessian_sparsity, (term.variable_index_1.value,
                                 term.variable_index_2.value))
    end
end

function hessian_lagrangian_structure(model::Optimizer)
    hessian_sparsity = Tuple{Int64,Int64}[]
    if !model.nlp_data.has_objective && model.objective !== nothing
        append_to_hessian_sparsity!(hessian_sparsity, model.objective)
    end
    for (quad, set) in model.quadratic_le_constraints
        append_to_hessian_sparsity!(hessian_sparsity, quad)
    end
    for (quad, set) in model.quadratic_ge_constraints
        append_to_hessian_sparsity!(hessian_sparsity, quad)
    end
    for (quad, set) in model.quadratic_eq_constraints
        append_to_hessian_sparsity!(hessian_sparsity, quad)
    end
    nlp_hessian_sparsity = MOI.hessian_lagrangian_structure(model.nlp_data.evaluator)
    append!(hessian_sparsity, nlp_hessian_sparsity)
    return hessian_sparsity
end

function eval_function(var::MOI.SingleVariable, x)
    return x[var.variable.value]
end

function eval_function(aff::MOI.ScalarAffineFunction, x)
    function_value = aff.constant
    for term in aff.terms
        # Note the implicit assumtion that VariableIndex values match up with
        # x indices. This is valid because in this wrapper ListOfVariableIndices
        # is always [1, ..., NumberOfVariables].
        function_value += term.coefficient*x[term.variable_index.value]
    end
    return function_value
end

function eval_function(quad::MOI.ScalarQuadraticFunction, x)
    function_value = quad.constant
    for term in quad.affine_terms
        function_value += term.coefficient*x[term.variable_index.value]
    end
    for term in quad.quadratic_terms
        row_idx = term.variable_index_1
        col_idx = term.variable_index_2
        coefficient = term.coefficient
        if row_idx == col_idx
            function_value += 0.5*coefficient*x[row_idx.value]*x[col_idx.value]
        else
            function_value += coefficient*x[row_idx.value]*x[col_idx.value]
        end
    end
    return function_value
end

function eval_objective(model::Optimizer, x)
    # The order of the conditions is important. NLP objectives override regular
    # objectives.
    if model.nlp_data.has_objective
        return MOI.eval_objective(model.nlp_data.evaluator, x)
    elseif model.objective !== nothing
        return eval_function(model.objective, x)
    else
        # No objective function set. This could happen with FEASIBILITY_SENSE.
        return 0.0
    end
end

function fill_gradient!(grad, x, var::MOI.SingleVariable)
    fill!(grad, 0.0)
    grad[var.variable.value] = 1.0
end

function fill_gradient!(grad, x, aff::MOI.ScalarAffineFunction{Float64})
    fill!(grad, 0.0)
    for term in aff.terms
        grad[term.variable_index.value] += term.coefficient
    end
end

function fill_gradient!(grad, x, quad::MOI.ScalarQuadraticFunction{Float64})
    fill!(grad, 0.0)
    for term in quad.affine_terms
        grad[term.variable_index.value] += term.coefficient
    end
    for term in quad.quadratic_terms
        row_idx = term.variable_index_1
        col_idx = term.variable_index_2
        coefficient = term.coefficient
        if row_idx == col_idx
            grad[row_idx.value] += coefficient*x[row_idx.value]
        else
            grad[row_idx.value] += coefficient*x[col_idx.value]
            grad[col_idx.value] += coefficient*x[row_idx.value]
        end
    end
end

function eval_objective_gradient(model::Optimizer, grad, x)
    if model.nlp_data.has_objective
        MOI.eval_objective_gradient(model.nlp_data.evaluator, grad, x)
    elseif model.objective !== nothing
        fill_gradient!(grad, x, model.objective)
    end
    return
end

# Refers to local variables in eval_constraint() below.
macro eval_function(array_name)
    escrow = esc(:row)
    quote
        for (func, set) in $(esc(array_name))
            $(esc(:g))[$escrow] = eval_function(func, $(esc(:x)))
            $escrow += 1
        end
    end
end

function eval_constraint(model::Optimizer, g, x)
    row = 1
    @eval_function model.linear_le_constraints
    @eval_function model.linear_ge_constraints
    @eval_function model.linear_eq_constraints
    @eval_function model.quadratic_le_constraints
    @eval_function model.quadratic_ge_constraints
    @eval_function model.quadratic_eq_constraints
    nlp_g = view(g, row:length(g))
    MOI.eval_constraint(model.nlp_data.evaluator, nlp_g, x)
    return
end

function fill_constraint_jacobian!(values, start_offset, x, aff::MOI.ScalarAffineFunction)
    num_coefficients = length(aff.terms)
    for i in 1:num_coefficients
        values[start_offset+i] = aff.terms[i].coefficient
    end
    return num_coefficients
end

function fill_constraint_jacobian!(values, start_offset, x, quad::MOI.ScalarQuadraticFunction)
    num_affine_coefficients = length(quad.affine_terms)
    for i in 1:num_affine_coefficients
        values[start_offset+i] = quad.affine_terms[i].coefficient
    end
    num_quadratic_coefficients = 0
    for term in quad.quadratic_terms
        row_idx = term.variable_index_1
        col_idx = term.variable_index_2
        coefficient = term.coefficient
        if row_idx == col_idx
            values[start_offset+num_affine_coefficients+num_quadratic_coefficients+1] = coefficient*x[col_idx.value]
            num_quadratic_coefficients += 1
        else
            # Note that the order matches the Jacobian sparsity pattern.
            values[start_offset+num_affine_coefficients+num_quadratic_coefficients+1] = coefficient*x[col_idx.value]
            values[start_offset+num_affine_coefficients+num_quadratic_coefficients+2] = coefficient*x[row_idx.value]
            num_quadratic_coefficients += 2
        end
    end
    return num_affine_coefficients + num_quadratic_coefficients
end

# Refers to local variables in eval_constraint_jacobian() below.
macro fill_constraint_jacobian(array_name)
    esc_offset = esc(:offset)
    quote
        for (func, set) in $(esc(array_name))
            $esc_offset += fill_constraint_jacobian!($(esc(:values)),
                                                     $esc_offset, $(esc(:x)),
                                                     func)
        end
    end
end

function eval_constraint_jacobian(model::Optimizer, values, x)
    offset = 0
    @fill_constraint_jacobian model.linear_le_constraints
    @fill_constraint_jacobian model.linear_ge_constraints
    @fill_constraint_jacobian model.linear_eq_constraints
    @fill_constraint_jacobian model.quadratic_le_constraints
    @fill_constraint_jacobian model.quadratic_ge_constraints
    @fill_constraint_jacobian model.quadratic_eq_constraints

    nlp_values = view(values, 1+offset:length(values))
    MOI.eval_constraint_jacobian(model.nlp_data.evaluator, nlp_values, x)
    return
end

function fill_hessian_lagrangian!(values, start_offset, scale_factor,
                                  ::Union{MOI.SingleVariable,
                                          MOI.ScalarAffineFunction,Nothing})
    return 0
end

function fill_hessian_lagrangian!(values, start_offset, scale_factor,
                                  quad::MOI.ScalarQuadraticFunction)
    for i in 1:length(quad.quadratic_terms)
        values[start_offset + i] = scale_factor*quad.quadratic_terms[i].coefficient
    end
    return length(quad.quadratic_terms)
end

function eval_hessian_lagrangian(model::Optimizer, values, x, obj_factor, lambda)
    offset = 0
    if !model.nlp_data.has_objective
        offset += fill_hessian_lagrangian!(values, 0, obj_factor,
                                          model.objective)
    end
    for (i, (quad, set)) in enumerate(model.quadratic_le_constraints)
        offset += fill_hessian_lagrangian!(values, offset, lambda[i+quadratic_le_offset(model)], quad)
    end
    for (i, (quad, set)) in enumerate(model.quadratic_ge_constraints)
        offset += fill_hessian_lagrangian!(values, offset, lambda[i+quadratic_ge_offset(model)], quad)
    end
    for (i, (quad, set)) in enumerate(model.quadratic_eq_constraints)
        offset += fill_hessian_lagrangian!(values, offset, lambda[i+quadratic_eq_offset(model)], quad)
    end
    nlp_values = view(values, 1 + offset : length(values))
    nlp_lambda = view(lambda, 1 + nlp_constraint_offset(model) : length(lambda))
    MOI.eval_hessian_lagrangian(model.nlp_data.evaluator, nlp_values, x, obj_factor, nlp_lambda)
end

function constraint_bounds(model::Optimizer)
    constraint_lb = Float64[]
    constraint_ub = Float64[]
    for (func, set) in model.linear_le_constraints
        push!(constraint_lb, -Inf)
        push!(constraint_ub, set.upper)
    end
    for (func, set) in model.linear_ge_constraints
        push!(constraint_lb, set.lower)
        push!(constraint_ub, Inf)
    end
    for (func, set) in model.linear_eq_constraints
        push!(constraint_lb, set.value)
        push!(constraint_ub, set.value)
    end
    for (func, set) in model.quadratic_le_constraints
        push!(constraint_lb, -Inf)
        push!(constraint_ub, set.upper)
    end
    for (func, set) in model.quadratic_ge_constraints
        push!(constraint_lb, set.lower)
        push!(constraint_ub, Inf)
    end
    for (func, set) in model.quadratic_eq_constraints
        push!(constraint_lb, set.value)
        push!(constraint_ub, set.value)
    end
    for bound in model.nlp_data.constraint_bounds
        push!(constraint_lb, bound.lower)
        push!(constraint_ub, bound.upper)
    end
    return constraint_lb, constraint_ub
end

function MOI.optimize!(model::Optimizer)
    # TODO: Reuse model.inner for incremental solves if possible.
    num_variables = length(model.variable_info)
    num_linear_le_constraints = length(model.linear_le_constraints)
    num_linear_ge_constraints = length(model.linear_ge_constraints)
    num_linear_eq_constraints = length(model.linear_eq_constraints)
    nlp_row_offset = nlp_constraint_offset(model)
    num_quadratic_constraints = nlp_constraint_offset(model) -
                                quadratic_le_offset(model)
    num_nlp_constraints = length(model.nlp_data.constraint_bounds)
    num_constraints = num_nlp_constraints + nlp_row_offset

    evaluator = model.nlp_data.evaluator
    features = MOI.features_available(evaluator)
    has_hessian = (:Hess in features)
    init_feat = [:Grad]
    has_hessian && push!(init_feat, :Hess)
    num_nlp_constraints > 0 && push!(init_feat, :Jac)

    MOI.initialize(evaluator, init_feat)
    jacobian_sparsity = jacobian_structure(model)
    hessian_sparsity = has_hessian ? hessian_lagrangian_structure(model) : []

    # Objective callback
    if model.sense == MOI.MIN_SENSE
        objective_scale = 1.0
    elseif model.sense == MOI.MAX_SENSE
        objective_scale = -1.0
    else # FEASIBILITY_SENSE
        # TODO: This could produce confusing solver output if a nonzero
        # objective is set.
        objective_scale = 0.0
    end

    eval_f_cb(x) = objective_scale * eval_objective(model, x)

    # Objective gradient callback
    function eval_grad_f_cb(x, grad_f)
        eval_objective_gradient(model, grad_f, x)
        Compat.rmul!(grad_f,objective_scale)
    end

    # Constraint value callback
    eval_g_cb(x, g) = eval_constraint(model, g, x)

    # Jacobian callback
    function eval_jac_g_cb(x, mode, rows, cols, values)
        if mode == :Structure
            for i in 1:length(jacobian_sparsity)
                rows[i] = jacobian_sparsity[i][1]
                cols[i] = jacobian_sparsity[i][2]
            end
        else
            eval_constraint_jacobian(model, values, x)
        end
    end

    if has_hessian
        # Hessian callback
        function eval_h_cb(x, mode, rows, cols, obj_factor,
            lambda, values)
            if mode == :Structure
                for i in 1:length(hessian_sparsity)
                    rows[i] = hessian_sparsity[i][1]
                    cols[i] = hessian_sparsity[i][2]
                end
            else
                obj_factor *= objective_scale
                eval_hessian_lagrangian(model, values, x, objective_scale*obj_factor, lambda)
            end
        end
    else
        eval_h_cb = nothing
    end

    x_l = [v.lower_bound for v in model.variable_info]
    x_u = [v.upper_bound for v in model.variable_info]

    constraint_lb, constraint_ub = constraint_bounds(model)

    model.inner = createProblem(num_variables, x_l, x_u, num_constraints,
                            constraint_lb, constraint_ub,
                            length(jacobian_sparsity),
                            length(hessian_sparsity),
                            eval_f_cb, eval_g_cb, eval_grad_f_cb, eval_jac_g_cb,
                            eval_h_cb)

    # Ipopt crashes by default if NaN/Inf values are returned from the
    # evaluation callbacks. This option tells Ipopt to explicitly check for them
    # and return Invalid_Number_Detected instead. This setting may result in a
    # minor performance loss and can be overwritten by specifying
    # check_derivatives_for_naninf="no".
    addOption(model.inner, "check_derivatives_for_naninf", "yes")

    if !has_hessian
        addOption(model.inner, "hessian_approximation", "limited-memory")
    end
    if num_nlp_constraints == 0 && num_quadratic_constraints == 0
        addOption(model.inner, "jac_c_constant", "yes")
        addOption(model.inner, "jac_d_constant", "yes")
        if !model.nlp_data.has_objective
            # We turn on this option if all constraints are linear and the
            # objective is linear or quadratic. From the documentation, it's
            # unclear if it may also apply if the constraints are at most
            # quadratic.
            addOption(model.inner, "hessian_constant", "yes")
        end
    end

    model.inner.x = [v.start for v in model.variable_info]

    for (name,value) in model.options
        sname = string(name)
        if match(r"(^resto_)", sname) != nothing
            sname = replace(sname, r"(^resto_)", "resto.")
        end
        addOption(model.inner, sname, value)
    end
    solveProblem(model.inner)
end

function MOI.get(model::Optimizer, ::MOI.TerminationStatus)
    if model.inner === nothing
        return MOI.OPTIMIZE_NOT_CALLED
    end
    status = ApplicationReturnStatus[model.inner.status]
    if status == :Solve_Succeeded || status == :Feasible_Point_Found
        return MOI.LOCALLY_SOLVED
    elseif status == :Infeasible_Problem_Detected
        return MOI.LOCALLY_INFEASIBLE
    elseif status == :Solved_To_Acceptable_Level
        return MOI.ALMOST_LOCALLY_SOLVED
    elseif status == :Search_Direction_Becomes_Too_Small
        return MOI.NUMERICAL_ERROR
    elseif status == :Diverging_Iterates
        return MOI.NORM_LIMIT
    elseif status == :User_Requested_Stop
        return MOI.INTERRUPTED
    elseif status == :Maximum_Iterations_Exceeded
        return MOI.ITERATION_LIMIT
    elseif status == :Maximum_CpuTime_Exceeded
        return MOI.TIME_LIMIT
    elseif status == :Restoration_Failed
        return MOI.NUMERICAL_ERROR
    elseif status == :Error_In_Step_Computation
        return MOI.NUMERICAL_ERROR
    elseif status == :Invalid_Option
        return MOI.INVALID_OPTION
    elseif status == :Not_Enough_Degrees_Of_Freedom
        return MOI.INVALID_MODEL
    elseif status == :Invalid_Problem_Definition
        return MOI.INVALID_MODEL
    elseif status == :Invalid_Number_Detected
        return MOI.INVALID_MODEL
    elseif status == :Unrecoverable_Exception
        return MOI.OTHER_ERROR
    elseif status == :NonIpopt_Exception_Thrown
        return MOI.OHTER_ERROR
    elseif status == :Insufficient_Memory
        return MOI.MEMORY_LIMIT
    else
        error("Unrecognized Ipopt status $status")
    end
end

# Ipopt always has an iterate available.
function MOI.get(model::Optimizer, ::MOI.ResultCount)
    return (model.inner !== nothing) ? 1 : 0
end

function MOI.get(model::Optimizer, ::MOI.PrimalStatus)
    if model.inner === nothing
        return MOI.NO_SOLUTION
    end
    status = ApplicationReturnStatus[model.inner.status]
    if status == :Solve_Succeeded
        return MOI.FEASIBLE_POINT
    elseif status == :Feasible_Point_Found
        return MOI.FEASIBLE_POINT
    elseif status == :Solved_To_Acceptable_Level
        # Solutions are only guaranteed to satisfy the "acceptable" convergence
        # tolerances.
        return MOI.NEARLY_FEASIBLE_POINT
    elseif status == :Infeasible_Problem_Detected
        return MOI.INFEASIBLE_POINT
    else
        return MOI.UNKNOWN_RESULT_STATUS
    end
end

function MOI.get(model::Optimizer, ::MOI.DualStatus)
    if model.inner === nothing
        return MOI.NO_SOLUTION
    end
    status = ApplicationReturnStatus[model.inner.status]
    if status == :Solve_Succeeded
        return MOI.FEASIBLE_POINT
    elseif status == :Feasible_Point_Found
        return MOI.FEASIBLE_POINT
    elseif status == :Solved_To_Acceptable_Level
        # Solutions are only guaranteed to satisfy the "acceptable" convergence
        # tolerances.
        return MOI.NEARLY_FEASIBLE_POINT
    elseif status == :Infeasible_Problem_Detected
        # TODO: What is the interpretation of the dual in this case?
        return MOI.UNKNOWN_RESULT_STATUS
    else
        return MOI.UNKNOWN_RESULT_STATUS
    end
end

function MOI.get(model::Optimizer, ::MOI.ObjectiveValue)
    if model.inner === nothing
        error("ObjectiveValue not available.")
    end
    scale = (model.sense == MOI.MAX_SENSE) ? -1 : 1
    return scale * model.inner.obj_val
end

# TODO: This is a bit off, because the variable primal should be available
# only after a solve. If model.inner is initialized but we haven't solved, then
# the primal values we return do not have the intended meaning.
function MOI.get(model::Optimizer, ::MOI.VariablePrimal, vi::MOI.VariableIndex)
    if model.inner === nothing
        error("VariablePrimal not available.")
    end
    check_inbounds(model, vi)
    return model.inner.x[vi.value]
end

macro define_constraint_primal(function_type, set_type, prefix)
    constraint_array = Symbol(string(prefix) * "_constraints")
    offset_function = Symbol(string(prefix) * "_offset")
    quote
        function MOI.get(model::Optimizer, ::MOI.ConstraintPrimal,
                         ci::MOI.ConstraintIndex{$function_type, $set_type})
            if model.inner === nothing
                error("ConstraintPrimal not available.")
            end
            if !(1 <= ci.value <= length(model.$(constraint_array)))
                error("Invalid constraint index ", ci.value)
            end
            return model.inner.g[ci.value + $offset_function(model)]
        end
    end
end

@define_constraint_primal(MOI.ScalarAffineFunction{Float64},
                          MOI.LessThan{Float64}, linear_le)
@define_constraint_primal(MOI.ScalarAffineFunction{Float64},
                          MOI.GreaterThan{Float64}, linear_ge)
@define_constraint_primal(MOI.ScalarAffineFunction{Float64},
                          MOI.EqualTo{Float64}, linear_eq)
@define_constraint_primal(MOI.ScalarQuadraticFunction{Float64},
                          MOI.LessThan{Float64}, quadratic_le)
@define_constraint_primal(MOI.ScalarQuadraticFunction{Float64},
                          MOI.GreaterThan{Float64}, quadratic_ge)
@define_constraint_primal(MOI.ScalarQuadraticFunction{Float64},
                          MOI.EqualTo{Float64}, quadratic_eq)

function MOI.get(model::Optimizer, ::MOI.ConstraintPrimal,
                 ci::MOI.ConstraintIndex{MOI.SingleVariable,
                                         MOI.LessThan{Float64}})
    if model.inner === nothing
        error("ConstraintPrimal not available.")
    end
    vi = MOI.VariableIndex(ci.value)
    check_inbounds(model, vi)
    if !has_upper_bound(model, vi)
        error("Variable $vi has no upper bound -- ConstraintPrimal not defined.")
    end
    return model.inner.x[vi.value]
end

function MOI.get(model::Optimizer, ::MOI.ConstraintPrimal,
                 ci::MOI.ConstraintIndex{MOI.SingleVariable,
                                         MOI.GreaterThan{Float64}})
    if model.inner === nothing
        error("ConstraintPrimal not available.")
    end
    vi = MOI.VariableIndex(ci.value)
    check_inbounds(model, vi)
    if !has_lower_bound(model, vi)
        error("Variable $vi has no lower bound -- ConstraintPrimal not defined.")
    end
    return model.inner.x[vi.value]
end

function MOI.get(model::Optimizer, ::MOI.ConstraintPrimal,
                 ci::MOI.ConstraintIndex{MOI.SingleVariable,
                                         MOI.EqualTo{Float64}})
    if model.inner === nothing
        error("ConstraintPrimal not available.")
    end
    vi = MOI.VariableIndex(ci.value)
    check_inbounds(model, vi)
    if !is_fixed(model, vi)
        error("Variable $vi is not fixed -- ConstraintPrimal not defined.")
    end
    return model.inner.x[vi.value]
end

function MOI.get(model::Optimizer, ::MOI.ConstraintDual,
                 ci::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},
                                         MOI.LessThan{Float64}})
    if model.inner === nothing
        error("ConstraintDual not available.")
    end
    @assert 1 <= ci.value <= length(model.linear_le_constraints)
    # TODO: Unable to find documentation in Ipopt about the signs of duals.
    # Rescaling by -1 here seems to pass the MOI tests.
    return -1 * model.inner.mult_g[ci.value + linear_le_offset(model)]
end

function MOI.get(model::Optimizer, ::MOI.ConstraintDual,
                 ci::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},
                                         MOI.GreaterThan{Float64}})
    if model.inner === nothing
        error("ConstraintDual not available.")
    end
    @assert 1 <= ci.value <= length(model.linear_ge_constraints)
    # TODO: Unable to find documentation in Ipopt about the signs of duals.
    # Rescaling by -1 here seems to pass the MOI tests.
    return -1 * model.inner.mult_g[ci.value + linear_ge_offset(model)]
end

function MOI.get(model::Optimizer, ::MOI.ConstraintDual,
                 ci::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},
                                         MOI.EqualTo{Float64}})
    if model.inner === nothing
        error("ConstraintDual not available.")
    end
    @assert 1 <= ci.value <= length(model.linear_eq_constraints)
    # TODO: Rescaling by -1 for consistency, but I don't know if this is covered
    # by tests.
    return -1 * model.inner.mult_g[ci.value + linear_eq_offset(model)]
end

function MOI.get(model::Optimizer, ::MOI.ConstraintDual,
                 ci::MOI.ConstraintIndex{MOI.SingleVariable,
                                         MOI.LessThan{Float64}})
    if model.inner === nothing
        error("ConstraintDual not available.")
    end
    vi = MOI.VariableIndex(ci.value)
    check_inbounds(model, vi)
    if !has_upper_bound(model, vi)
        error("Variable $vi has no upper bound -- ConstraintDual not defined.")
    end
    # MOI convention is for feasible LessThan duals to be nonpositive.
    return -1 * model.inner.mult_x_U[vi.value]
end

function MOI.get(model::Optimizer, ::MOI.ConstraintDual,
                 ci::MOI.ConstraintIndex{MOI.SingleVariable,
                                         MOI.GreaterThan{Float64}})
    if model.inner === nothing
        error("ConstraintDual not available.")
    end
    vi = MOI.VariableIndex(ci.value)
    check_inbounds(model, vi)
    if !has_lower_bound(model, vi)
        error("Variable $vi has no lower bound -- ConstraintDual not defined.")
    end
    return model.inner.mult_x_L[vi.value]
end

function MOI.get(model::Optimizer, ::MOI.ConstraintDual,
                 ci::MOI.ConstraintIndex{MOI.SingleVariable,
                                         MOI.EqualTo{Float64}})
    if model.inner === nothing
        error("ConstraintDual not available.")
    end
    vi = MOI.VariableIndex(ci.value)
    check_inbounds(model, vi)
    if !is_fixed(model, vi)
        error("Variable $vi is not fixed -- ConstraintDual not defined.")
    end
    return model.inner.mult_x_L[vi.value] - model.inner.mult_x_U[vi.value]
end

function MOI.get(model::Optimizer, ::MOI.NLPBlockDual)
    if model.inner === nothing
        error("NLPBlockDual not available.")
    end
    return -1 * model.inner.mult_g[(1 + nlp_constraint_offset(model)):end]
end


text:
using Test
using Once

a = "foo"

@once f()::String = begin
    global a
    a *= "bar"
    v = a
    a *= "bar"
    return v
end

@test f() == "foobar"
@test a == "foobarbar"
a = ""
@test f() == "foobar"
@test a == ""
@test f() == "foobar"
@test a == ""

a = 7

@once f2()::Int = begin
    global a
    a *= 2
    v = a
    a *= 2
    return v
end

@test f2() == 14
@test a == 28
a = 1
@test f2() == 14
@test a == 1
@test f2() == 14
@test a == 1


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Pangomm_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Pangomm")
JLLWrappers.@generate_main_file("Pangomm", UUID("9c53b654-4175-57d2-a160-8980ed551c15"))
end  # module Pangomm_jll