text:
using Test

using BioGenerics
using FormatSpecimens
using GenomicFeatures
using XAM

import BioAlignments: Alignment, AlignmentAnchor, OP_START, OP_MATCH, OP_DELETE
import BGZFStreams: BGZFStream
import BioGenerics.Exceptions: MissingFieldException
import BioSequences: @dna_str, @aa_str

# Generate a random range within `range`.
function randrange(range)
    x = rand(range)
    y = rand(range)
    if x < y
        return x:y
    else
        return y:x
    end
end

@testset "SAM" begin
    samdir = path_of_format("SAM")

    @testset "MetaInfo" begin
        metainfo = SAM.MetaInfo()
        @test !isfilled(metainfo)
        @test occursin("not filled", repr(metainfo))

        metainfo = SAM.MetaInfo("CO", "some comment (parens)")
        @test isfilled(metainfo)
        @test string(metainfo) == "@CO\tsome comment (parens)"
        @test occursin("CO", repr(metainfo))
        @test SAM.tag(metainfo) == "CO"
        @test SAM.value(metainfo) == "some comment (parens)"
        @test_throws ArgumentError keys(metainfo)
        @test_throws ArgumentError values(metainfo)

        metainfo = SAM.MetaInfo("HD", ["VN" => "1.0", "SO" => "coordinate"])
        @test isfilled(metainfo)
        @test string(metainfo) == "@HD\tVN:1.0\tSO:coordinate"
        @test occursin("HD", repr(metainfo))
        @test SAM.tag(metainfo) == "HD"
        @test SAM.value(metainfo) == "VN:1.0\tSO:coordinate"
        @test keys(metainfo) == ["VN", "SO"]
        @test values(metainfo) == ["1.0", "coordinate"]
        @test SAM.keyvalues(metainfo) == ["VN" => "1.0", "SO" => "coordinate"]
        @test haskey(metainfo, "VN")
        @test haskey(metainfo, "SO")
        @test !haskey(metainfo, "GO")
        @test metainfo["VN"] == "1.0"
        @test metainfo["SO"] == "coordinate"
        @test_throws KeyError metainfo["GO"]
    end

    @testset "Header" begin
        header = SAM.Header()
        @test isempty(header)
        push!(header, SAM.MetaInfo("@HD\tVN:1.0\tSO:coordinate"))
        @test !isempty(header)
        @test length(header) == 1
        push!(header, SAM.MetaInfo("@CO\tsome comment"))
        @test length(header) == 2
        @test isa(collect(header), Vector{SAM.MetaInfo})
    end

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

        record = SAM.Record("r001\t99\tchr1\t7\t30\t8M2I4M1D3M\t=\t37\t39\tTTAGATAAAGGATACTG\t*")
        @test isfilled(record)
        @test occursin(r"^XAM.SAM.Record:\n", repr(record))
        @test SAM.ismapped(record)
        @test SAM.isprimary(record)
        @test SAM.hastempname(record)
        @test SAM.tempname(record) == "r001"
        @test SAM.hasflag(record)
        @test SAM.flag(record) === UInt16(99)
        @test SAM.hasrefname(record)
        @test SAM.refname(record) == "chr1"
        @test SAM.hasposition(record)
        @test SAM.position(record) === 7
        @test SAM.hasmappingquality(record)
        @test SAM.mappingquality(record) === UInt8(30)
        @test SAM.hascigar(record)
        @test SAM.cigar(record) == "8M2I4M1D3M"
        @test SAM.hasnextrefname(record)
        @test SAM.nextrefname(record) == "="
        @test SAM.hasnextposition(record)
        @test SAM.nextposition(record) === 37
        @test SAM.hastemplength(record)
        @test SAM.templength(record) === 39
        @test SAM.hassequence(record)
        @test SAM.sequence(record) == dna"TTAGATAAAGGATACTG"
        @test !SAM.hasquality(record)
        @test_throws MissingFieldException SAM.quality(record)
    end

    @testset "Reader" begin
        reader = open(SAM.Reader, joinpath(samdir, "ce#1.sam"))
        @test isa(reader, SAM.Reader)
        @test eltype(reader) === SAM.Record

        # header
        h = header(reader)
        @test string.(findall(header(reader), "SQ")) == ["@SQ\tSN:CHROMOSOME_I\tLN:1009800"]

        # first record
        record = SAM.Record()
        read!(reader, record)
        @test SAM.ismapped(record)
        @test SAM.refname(record) == "CHROMOSOME_I"
        @test SAM.position(record) == leftposition(record) == 2
        @test SAM.rightposition(record) == rightposition(record) == 102
        @test SAM.tempname(record) == seqname(record) == "SRR065390.14978392"
        @test SAM.sequence(record) == sequence(record) == dna"CCTAGCCCTAACCCTAACCCTAACCCTAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAA"
        @test SAM.sequence(String, record)          ==    "CCTAGCCCTAACCCTAACCCTAACCCTAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAA"
        @test SAM.seqlength(record) == 100
        @test SAM.quality(record)         == (b"#############################@B?8B?BA@@DDBCDDCBC@CDCDCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC" .- 33)
        @test SAM.quality(String, record) ==   "#############################@B?8B?BA@@DDBCDDCBC@CDCDCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"
        @test SAM.flag(record) == 16
        @test SAM.cigar(record) == "27M1D73M"
        @test SAM.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 eof(reader)
        close(reader)

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

        # iterator
        @test length(collect(open(SAM.Reader, joinpath(samdir, "ce#1.sam")))) == 1
        @test length(collect(open(SAM.Reader, joinpath(samdir, "ce#2.sam")))) == 2

        # IOStream
        @test length(collect(SAM.Reader(open(joinpath(samdir, "ce#1.sam"))))) == 1
        @test length(collect(SAM.Reader(open(joinpath(samdir, "ce#2.sam"))))) == 2
    end

    @testset "Round trip" begin
        function compare_records(xs, ys)
            if length(xs) != length(ys)
                return false
            end
            for (x, y) in zip(xs, ys)
                if x.data[x.filled] != y.data[y.filled]
                    return false
                end
            end
            return true
        end
        for specimen in list_valid_specimens("SAM")
            filepath = joinpath(samdir, filename(specimen))
            mktemp() do path, io
                # copy
                reader = open(SAM.Reader, filepath)

                header_original = header(reader)

                writer = SAM.Writer(io, header_original)

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

                close(reader)
                close(writer)

                reader = open(SAM.Reader, path)

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

                close(reader)

            end
        end
    end
end

@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.bin_mq_nl == new_record.bin_mq_nl
        @test record.block_size == new_record.block_size
        @test record.flag_nc == new_record.flag_nc
        @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 (x, y) in zip(xs, ys)
            if !(
                x.block_size == y.block_size &&
                x.refid      == y.refid &&
                x.pos        == y.pos &&
                x.bin_mq_nl  == y.bin_mq_nl &&
                x.flag_nc    == y.flag_nc &&
                x.l_seq      == y.l_seq &&
                x.next_refid == y.next_refid &&
                x.next_pos   == y.next_pos &&
                x.tlen       == y.tlen &&
                x.data[1:BAM.data_size(x)] == y.data[1:BAM.data_size(y)])
                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 "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:
using HDF5
using Test

@testset "properties" begin

fn = tempname()
h5open(fn, "w";
       userblock = 1024,
       alignment = (0, sizeof(Int)),
       libver_bounds = (HDF5.H5F_LIBVER_EARLIEST, HDF5.H5F_LIBVER_LATEST),
      ) do hfile
    # generic
    g = g_create(hfile, "group")
    d = d_create(g, "dataset", datatype(Int), dataspace((500,500)),
                 alloc_time = HDF5.H5D_ALLOC_TIME_EARLY,
                 chunk = (5, 10),
                 track_times = false)
    attrs(d)["metadata"] = "test"

    flush(hfile)

    fcpl = get_create_properties(hfile)
    fapl = get_access_properties(hfile)
    gcpl = get_create_properties(hfile["group"])
    dcpl = get_create_properties(hfile["group/dataset"])
    acpl = get_create_properties(attrs(hfile["group/dataset"])["metadata"])

    # Retrievability of properties
    @test isvalid(fcpl)
    @test isvalid(fapl)
    @test isvalid(gcpl)
    @test isvalid(dcpl)
    @test isvalid(acpl)

    # Retrieving property values:
    @test fcpl[:userblock] == 1024
    @test fcpl[:track_times] == true

    @test fapl[:alignment] == (0, sizeof(Int))
    # value is H5FD_SEC2, but "constant" is runtime value not loadable by _read_const()
    @test HDF5.h5i_get_type(fapl[:driver]) == HDF5.H5I_VFL
    # Docs say h5p_get_driver_info() doesn't error, but it does print an error message...
    #   https://portal.hdfgroup.org/display/HDF5/H5P_GET_DRIVER_INFO
    HDF5.hiding_errors() do
        @test fapl[:driver_info] == C_NULL
    end
    @test fapl[:fclose_degree] == HDF5.H5F_CLOSE_STRONG
    @test fapl[:libver_bounds] == (HDF5.H5F_LIBVER_EARLIEST, HDF5.H5F_LIBVER_LATEST)

    @test gcpl[:local_heap_size_hint] == 0
    @test gcpl[:track_times] == true

    @test HDF5.UTF8_LINK_PROPERTIES[][:char_encoding] == HDF5.H5T_CSET_UTF8
    @test HDF5.UTF8_LINK_PROPERTIES[][:create_intermediate_group] == 1

    @test dcpl[:alloc_time] == HDF5.H5D_ALLOC_TIME_EARLY
    @test dcpl[:chunk] == (5, 10)
    @test dcpl[:layout] == HDF5.H5D_CHUNKED
    @test dcpl[:track_times] == false

    @test acpl[:char_encoding] == HDF5.H5T_CSET_UTF8

    nothing
end

rm(fn, force=true)

end

ext: jl


text:
const VectorCone = Union{MathOptInterface.SecondOrderCone,
                         MathOptInterface.RotatedSecondOrderCone,
                         MathOptInterface.PowerCone,
                         MathOptInterface.DualPowerCone,
                         MathOptInterface.ExponentialCone,
                         MathOptInterface.DualExponentialCone}
const PositiveSemidefiniteCone = Union{MathOptInterface.PositiveSemidefiniteConeTriangle,
                                       MathOptInterface.PositiveSemidefiniteConeScaled}
const LinearFunction = Union{MathOptInterface.SingleVariable,
                             MathOptInterface.VectorOfVariables,
                             MathOptInterface.ScalarAffineFunction,
                             MathOptInterface.VectorAffineFunction}
const AffineFunction = Union{MathOptInterface.ScalarAffineFunction,
                             MathOptInterface.VectorAffineFunction}
const VariableFunction = Union{MathOptInterface.ScalarAffineFunction,
                               MathOptInterface.VectorAffineFunction}

const ScalarLinearDomain = Union{MathOptInterface.LessThan{Float64},
                                 MathOptInterface.GreaterThan{Float64},
                                 MathOptInterface.EqualTo{Float64},
                                 MathOptInterface.Interval{Float64} }
const VectorLinearDomain = Union{MathOptInterface.Nonpositives,
                                 MathOptInterface.Nonnegatives,
                                 MathOptInterface.Reals,
                                 MathOptInterface.Zeros}
const LinearDomain = Union{ScalarLinearDomain,VectorLinearDomain}
################################################################################
# ADD CONSTRAINT ###############################################################
################################################################################

function MathOptInterface.addconstraint!(
    m   :: MosekModel,
    axb :: MathOptInterface.ScalarAffineFunction{Float64},
    dom :: D) where {D <: MathOptInterface.AbstractScalarSet}

    N = 1
    conid = allocateconstraints(m,N)
    addlhsblock!(m, conid, fill(1,length(axb.variables)),axb.variables, axb.coefficients)

    if length(m.c_constant) < length(m.c_block)
        append!(m.c_constant,zeros(Float64,length(m.c_block) - length(m.c_constant)))
    end
    conidxs = getindexes(m.c_block,conid)
    m.c_constant[conidxs] = axb.constant

    addbound!(m,conid,conidxs,Float64[axb.constant],dom)

    conref = MathOptInterface.ConstraintReference{MathOptInterface.ScalarAffineFunction{Float64},D}(UInt64(conid) << 1)
    select(m.constrmap,MathOptInterface.ScalarAffineFunction{Float64},D)[conref.value] = conid

    conref
end

function MathOptInterface.addconstraint!(
    m   :: MosekModel,
    axb :: MathOptInterface.VectorAffineFunction{Float64},
    dom :: D) where { D <: MathOptInterface.AbstractVectorSet }

    N = MathOptInterface.dimension(dom)
    conid = allocateconstraints(m,N)
    addlhsblock!(m,
                 conid,
                 axb.outputindex,
                 axb.variables,
                 axb.coefficients)
    conidxs = getindexes(m.c_block,conid)
    m.c_constant[conidxs] = axb.constant

    m.c_block_slack[conid]

    addbound!(m,conid,conidxs,axb.constant,dom)
    conref = MathOptInterface.ConstraintReference{MathOptInterface.VectorAffineFunction{Float64},D}(UInt64(conid) << 1)
    select(m.constrmap,MathOptInterface.VectorAffineFunction{Float64},D)[conref.value] = conid

    conref
end

dim(dom :: MathOptInterface.PositiveSemidefiniteConeTriangle) = floor(Int,-.5 + sqrt(.25+2*MathOptInterface.dimension(dom)))
dim(dom :: MathOptInterface.PositiveSemidefiniteConeScaled) = floor(Int,-.5 + sqrt(.25+2*MathOptInterface.dimension(dom)))

function MathOptInterface.addconstraint!(m   :: MosekModel,
                                         axb :: MathOptInterface.VectorAffineFunction{Float64},
                                         dom :: PSDCone) where { PSDCone <: PositiveSemidefiniteCone }
    M = MathOptInterface.dimension(dom)
    N = dim(dom)

    conid = allocateconstraints(m,M)
    addlhsblock!(m,
                 conid,
                 axb.outputindex,
                 axb.variables,
                 axb.coefficients)
    conidxs = getindexes(m.c_block,conid)
    m.c_constant[conidxs] = axb.constant

    addbound!(m,conid,conidxs,axb.constant,dom)
    conref = MathOptInterface.ConstraintReference{MathOptInterface.VectorAffineFunction{Float64},PSDCone}(UInt64(conid) << 1)
    select(m.constrmap,MathOptInterface.VectorAffineFunction{Float64},PSDCone)[conref] = conid
    conref
end

################################################################################
# Variable constraints

# We allow following. Each variable can have
# - at most most upper and one lower bound
# - belong to at most one non-semidefinite cone
# - any number of semidefinite cones, which are implemented as ordinary constraints
# This is when things get a bit funky; By default a variable has no
# bounds, i.e. "free". Adding a `GreaterThan` or `Nonnegatives`
# constraint causes it to have a defined lower bound but no upper
# bound, allowing a `LessThan` or ``Nonpositives` constraint to be
# added later. Adding a `Interval` constraint defines both upper and
# lower bounds. Adding a `Reals` constraint will effectively be the
# same as an interval ]-infty;infty[, in that it will define both
# upper and lower bounds, and not allow those to be set afterwards.

domain_type_mask(dom :: MathOptInterface.Reals)        = (boundflag_lower | boundflag_upper)
domain_type_mask(dom :: MathOptInterface.Interval)     = (boundflag_lower | boundflag_upper)
domain_type_mask(dom :: MathOptInterface.EqualTo)      = (boundflag_lower | boundflag_upper)
domain_type_mask(dom :: MathOptInterface.GreaterThan)  = boundflag_lower
domain_type_mask(dom :: MathOptInterface.Nonnegatives) = boundflag_lower
domain_type_mask(dom :: MathOptInterface.LessThan)     = boundflag_upper
domain_type_mask(dom :: MathOptInterface.Nonpositives) = boundflag_upper

domain_type_mask(dom :: MathOptInterface.SecondOrderCone)        = boundflag_cone
domain_type_mask(dom :: MathOptInterface.RotatedSecondOrderCone) = boundflag_cone
domain_type_mask(dom :: MathOptInterface.ExponentialCone)        = boundflag_cone
domain_type_mask(dom :: MathOptInterface.PowerCone)              = boundflag_cone

domain_type_mask(dom :: MathOptInterface.PositiveSemidefiniteConeTriangle) = 0
domain_type_mask(dom :: MathOptInterface.PositiveSemidefiniteConeScaled)   = 0

domain_type_mask(dom :: MathOptInterface.Integer) = boundflag_int
domain_type_mask(dom :: MathOptInterface.ZeroOne) = (boundflag_int | boundflag_upper | boundflag_lower)

is_positivesemidefinite_set(dom :: Union{MathOptInterface.PositiveSemidefiniteConeTriangle,MathOptInterface.PositiveSemidefiniteConeScaled}) = true
is_positivesemidefinite_set(dom) = false
is_conic_set(dom :: VectorCone) = true
is_conic_set(dom) = false

addvarconstr(m :: MosekModel, subj :: Vector{Int}, dom :: MathOptInterface.Reals) = nothing
function addvarconstr(m :: MosekModel, subj :: Vector{Int}, dom :: MathOptInterface.Zeros)
    bnd = zeros(Float64,length(subj))
    putvarboundslice(m.task,subj,fill(MSK_BK_FX,length(subj)),bnd,bnd)
end
addvarconstr(m :: MosekModel, subj :: Vector{Int}, dom :: MathOptInterface.Interval) = putvarbound(m.task,subj[1],MSK_BK_RA,dom.lower, dom.upper)
addvarconstr(m :: MosekModel, subj :: Vector{Int}, dom :: MathOptInterface.EqualTo) = putvarbound(m.task,subj[1],MSK_BK_FX,dom.value, dom.value)
addvarconstr(m :: MosekModel, subj :: Vector{Int}, dom :: MathOptInterface.Integer) = putvartype(m.task,subj[1],MSK_VAR_TYPE_INT)
function addvarconstr(m :: MosekModel, subj :: Vector{Int}, dom :: MathOptInterface.ZeroOne)
    putvartype(m.task,subj[1],MSK_VAR_TYPE_INT)
    putvarbound(m.task,subj[1],MSK_BK_RA,0.0,1.0)
end

function addvarconstr(m :: MosekModel, subj :: Vector{Int}, dom :: MathOptInterface.LessThan)
    bk,lo,up = getvarbound(m.task,subj[1])
    bk = if (bk == MSK_BK_FR) MSK_BK_UP else MSK_BK_RA end
    putvarbound(m.task,Int32(subj[1]),bk,lo,dom.upper)
end

function addvarconstr(m :: MosekModel, subj :: Vector{Int}, dom :: MathOptInterface.GreaterThan)
    bk,lo,up = getvarbound(m.task,subj[1])
    bk = if (bk == MSK_BK_FR) MSK_BK_LO else MSK_BK_RA end
    putvarbound(m.task,Int32(subj[1]),bk,dom.lower,up)
end

function addvarconstr(m :: MosekModel, subj :: Vector{Int}, dom :: MathOptInterface.Nonnegatives)
    bkx = zeros(Int32,length(subj))
    blx = zeros(Float64,length(subj))
    bux = zeros(Float64,length(subj))
    for (i,j) in enumerate(subj)
        bk,lo,up = getvarbound(m.task,j)
        bkx[i] = if (bk == MSK_BK_FR) MSK_BK_RA else MSK_BK_LO end
        bux[i] = up
    end
    putvarboundlist(m.task,subj,bkx,blx,bux)
end

function addvarconstr(m :: MosekModel, subj :: Vector{Int}, dom :: MathOptInterface.Nonpositives)
    bkx = zeros(Int32,length(subj))
    blx = zeros(Float64,length(subj))
    bux = zeros(Float64,length(subj))
    for (i,j) in enumerate(subj)
        bk,lo,up = getvarbound(m.task,j)
        bkx[i] = if (bk == MSK_BK_FR) MSK_BK_RA else MSK_BK_UP end
        blx[i] = lo
    end
    putvarboundlist(m.task,subj,bkx,blx,bux)
end

abstractset2ct(dom::MathOptInterface.SecondOrderCone) = MSK_CT_QUAD
abstractset2ct(dom::MathOptInterface.RotatedSecondOrderCone) = MSK_CT_RQUAD
function MathOptInterface.addconstraint!(
    m   :: MosekModel,
    xs  :: MathOptInterface.SingleVariable,
    dom :: D) where {D <: MathOptInterface.AbstractScalarSet}

    subj = getindexes(m.x_block, ref2id(xs.variable))

    mask = domain_type_mask(dom)
    if mask & m.x_boundflags[subj[1]] != 0
        error("Cannot put multiple bound sets of the same type on a variable")
    end

    xcid = allocatevarconstraints(m,1)

    xc_sub = getindexes(m.xc_block,xcid)[1]

    m.xc_bounds[xcid]  = mask
    m.xc_idxs[xc_sub] = subj[1]

    addvarconstr(m,subj,dom)

    m.x_boundflags[subj[1]] .|= mask

    conref = MathOptInterface.ConstraintReference{MathOptInterface.SingleVariable,D}(UInt64(xcid) << 1)

    select(m.constrmap,MathOptInterface.SingleVariable,D)[conref.value] = xcid
    conref
end

function MathOptInterface.addconstraint!(m   :: MosekModel,
                                         xs  :: MathOptInterface.VectorOfVariables,
                                         dom :: D) where { D <: PositiveSemidefiniteCone }
    subj = Vector{Int}(length(xs.variables))
    for i in 1:length(subj)
        getindexes(m.x_block, ref2id(xs.variables[i]),subj,i)
    end

    mask = domain_type_mask(dom)
    if any(mask .& m.x_boundflags[subj[1]] .> 0)
        error("Cannot put multiple bound sets of the same type to a variable")
    end

    N = MathOptInterface.dimension(dom)

    if N < 3
        error("Invalid dimension for semidefinite constraint")
    end

    NN = (N*(N+1)) >> 1

    if length(subj) != NN
        error("Mismatching variable length for semidefinite constraint")
    end

    id = allocateconstraints(m,NN)

    appendbarvars(m.task,Int32(N))
    barvaridx = getnumbarvar(m.task)

    subi = getindexes(m.c_block,id)
    m.c_block_slack[subi] = -barvaridx

    subii32 = convert(Vector{Int32},subi)
    putaijlist(m.task,
               subii32,
               convert(Vector{Int32},subj),
               ones(Float64,NN))

    addbound!(m,id,subi,zeros(Float64,NN),dom)

    #putconboundlist(m.task,subii32,fill(MSK_BK_FX,NN),zeros(Float64,NN),zeros(Float64,NN))
    #idx = 1
    #for j in 1:N
    #    for i in 1:N
    #        symmatidx = appendsparsesymmat(m.task,N,Int32[i],Int32[j],Float64[-1.0])
    #        putbaraij(m.task,subii32[idx],barvaridx,[symmatidx],Float64[1.0])
    #        idx += 1
    #    end
    #end

    # HACK: We need to return a negative to indicate that this is
    # not, in fact, a real variable constraint, but rather a real
    # constraint, but to resolve return value at compile time we
    # need to disguise it as a variable constraint.
    #id2cref{MathOptInterface.VectorOfVariables,D}(-id)
    conref = MathOptInterface.ConstraintReference{MathOptInterface.VectorOfVariables,MathOptInterface.PositiveSemidefiniteConeTriangle}((UInt64(-id) << 1) | 1)
    select(m.constrmap,MathOptInterface.VectorOfVariables,MathOptInterface.PositiveSemidefiniteConeTriangle)[conref.value] = id

    conref
end

function aux_setvardom(m :: MosekModel,
                       xcid :: Int,
                       subj :: Vector{Int},
                       dom :: D) where { D <: VectorCone }
    appendcone(m.task,abstractset2ct(dom),  0.0, subj)
    coneidx = getnumcone(m.task)
    m.conecounter += 1
    putconename(m.task,coneidx,"$(m.conecounter)")
    m.xc_coneid[xcid] = m.conecounter
end
aux_setvardom(m :: MosekModel, xcid :: Int, subj :: Vector{Int},dom :: D) where {D <: MathOptInterface.AbstractSet} = addvarconstr(m,subj,dom)

function MathOptInterface.addconstraint!{D <: MathOptInterface.AbstractSet}(m :: MosekModel, xs :: MathOptInterface.VectorOfVariables, dom :: D)
    subj = Vector{Int}(length(xs.variables))
    for i in 1:length(subj)
        getindexes(m.x_block, ref2id(xs.variables[i]),subj,i)
    end

    mask = domain_type_mask(dom)
    if any(mask .& m.x_boundflags[subj[1]] .> 0)
        error("Cannot multiple bound sets of the same type to a variable")
    end

    N = MathOptInterface.dimension(dom)
    xcid = allocatevarconstraint(m,N)
    xc_sub = getindexes(m.xc_block,xcid)

    m.xc_bounds[xcid]  = mask
    m.xc_idxs[xc_sub] = subj

    aux_setvardom(m,xcid,subj,dom)

    m.x_boundflags[subj] .|= mask

    conref = MathOptInterface.ConstraintReference{MathOptInterface.VectorOfVariables,D}(UInt64(id) << 1)
    select(m.constrmap,MathOptInterface.VectorOfVariables,D)[conref.value] = id
    conref
end

################################################################################
################################################################################

# Put the linear left-hand side
function addlhsblock!(m        :: MosekModel,
                      conid    :: Int,
                      conidxs  :: Vector{Int},
                      varidxs  :: Vector{MathOptInterface.VariableReference},
                      cofs     :: Vector{Float64})
    consubi = getindexes(m.c_block,conid)
    subj = Array{Int}(length(varidxs))
    for i in 1:length(subj)
        getindexes(m.x_block,ref2id(varidxs[i]),subj,i)
    end

    N = length(consubi)

    At = sparse(subj, conidxs, cofs, getnumvar(m.task), N)
    putarowlist(m.task,convert(Vector{Int32},consubi),At)
end

addbound!(m :: MosekModel, conid :: Int, conidxs :: Vector{Int}, constant :: Vector{Float64}, dom :: MathOptInterface.Reals)                = putconboundlist(m.task,convert(Vector{Int32},conidxs),fill(MSK_BK_FR,MathOptInterface.dimension(dom)),-constant,-constant)
addbound!(m :: MosekModel, conid :: Int, conidxs :: Vector{Int}, constant :: Vector{Float64}, dom :: MathOptInterface.Zeros)                = putconboundlist(m.task,convert(Vector{Int32},conidxs),fill(MSK_BK_FX,MathOptInterface.dimension(dom)),-constant,-constant)
addbound!(m :: MosekModel, conid :: Int, conidxs :: Vector{Int}, constant :: Vector{Float64}, dom :: MathOptInterface.Nonnegatives)         = putconboundlist(m.task,convert(Vector{Int32},conidxs),fill(MSK_BK_LO,MathOptInterface.dimension(dom)),-constant,-constant)
addbound!(m :: MosekModel, conid :: Int, conidxs :: Vector{Int}, constant :: Vector{Float64}, dom :: MathOptInterface.Nonpositives)         = putconboundlist(m.task,convert(Vector{Int32},conidxs),fill(MSK_BK_UP,MathOptInterface.dimension(dom)),-constant,-constant)
addbound!(m :: MosekModel, conid :: Int, conidxs :: Vector{Int}, constant :: Vector{Float64}, dom :: MathOptInterface.GreaterThan{Float64}) = putconbound(m.task,Int32(conidxs[1]),MSK_BK_LO,dom.lower-constant[1],dom.lower-constant[1])
addbound!(m :: MosekModel, conid :: Int, conidxs :: Vector{Int}, constant :: Vector{Float64}, dom :: MathOptInterface.LessThan{Float64})    = putconbound(m.task,Int32(conidxs[1]),MSK_BK_UP,dom.upper-constant[1],dom.upper-constant[1])
addbound!(m :: MosekModel, conid :: Int, conidxs :: Vector{Int}, constant :: Vector{Float64}, dom :: MathOptInterface.EqualTo{Float64})     = putconbound(m.task,Int32(conidxs[1]),MSK_BK_FX,dom.value-constant[1],dom.value-constant[1])
addbound!(m :: MosekModel, conid :: Int, conidxs :: Vector{Int}, constant :: Vector{Float64}, dom :: MathOptInterface.Interval{Float64})    = putconbound(m.task,Int32(conidxs[1]),MSK_BK_RA,dom.lower-constant[1],dom.upper-constant[1])

function addbound!(m :: MosekModel, conid :: Int, conidxs :: Vector{Int}, constant :: Vector{Float64}, dom :: D) where { D <: VectorCone }
    N = MathOptInterface.dimension(dom)
    nalloc = ensurefree(m.x_block,N)

    varid = newblock(m.x_block,N)
    numvar = getnumvar(m.task)
    if length(m.x_block) > numvar
        appendvars(m.task, length(m.x_block) - numvar)
    end
    subj = getindexes(m.x_block,varid)

    putaijlist(m.task,conidxs,subj,-ones(Float64,N))
    putvarboundlist(m.task,subj,fill(MSK_BK_FR,N),zeros(Float64,N),zeros(Float64,N))
    putconboundlist(m.task,convert(Vector{Int32},conidxs),fill(MSK_BK_FX,N),-constant,-constant)

    m.c_block_slack[conid] = varid

    appendcone(m.task,abstractset2ct(dom),0.0,subj)
    coneidx = getnumcone(m.task)
    m.conecounter += 1
    putconename(m.task,coneidx,"$(m.conecounter)")
    m.c_coneid[conid] = m.conecounter
end

function addbound!(m :: MosekModel, conid :: Int, conidxs :: Vector{Int}, constant :: Vector{Float64}, dom :: MathOptInterface.PositiveSemidefiniteConeTriangle)
    dim = dim(dom)
    n = MathOptInterface.dimension(dom)
    appendbarvars(m.task,Int32[dim])
    barvaridx = getnumbarvar(m.task)

    idx = 1
    for j in 1:dim
        for i in j:dim
            matrixid = appendsparsesymmat(m.task,Int32(dim), Int32[i], Int32[j], Float64[-1.0])
            putbaraij(m.task, conidxs[idx], barvaridx, Int[matrixid], Float64[1.0])
            idx += 1
        end
    end

    putconboundlist(m.task,convert(Vector{Int32},conidxs),fill(MSK_BK_FX,length(constant)),-constant,-constant)

    m.c_block_slack[conid] = -barvaridx
end

function addbound!(m :: MosekModel, conid :: Int, conidxs :: Vector{Int}, constant :: Vector{Float64}, dom :: MathOptInterface.PositiveSemidefiniteConeScaled)
    dim = dim(dom)
    n = MathOptInterface.dimension(dom)

    appendbarvars(m.task,Int32[dim])
    barvaridx = getnumbarvar(m.task)

    idx = 1
    for j in 1:dim
        for i in j:dim
            if i != j
                matrixid = appendsparsesymmat(m.task,Int32(dim), Int32[i], Int32[j], Float64[- sqrt(2.0)])
            else
                matrixid = appendsparsesymmat(m.task,Int32(dim), Int32[i], Int32[j], Float64[-1.0])
            end
            putbaraij(m.task, conidxs[idx], barvaridx, Int[matrixid], Float64[1.0])
            idx += 1
        end
    end

    putconboundlist(m.task,convert(Vector{Int32},conidxs),fill(MSK_BK_FX,length(constant)),-constant,-constant)

    m.c_block_slack[conid] = -barvaridx
end

################################################################################
##  MODIFY #####################################################################
################################################################################

MathOptInterface.canmodifyconstraint(m::MosekModel, c::MathOptInterface.ConstraintReference{F,D}, dom) where { F <: LinearFunction, D <: VectorCone } = false
MathOptInterface.canmodifyconstraint(m::MosekModel, c::MathOptInterface.ConstraintReference{F,D}, dom) where { F <: LinearFunction, D <: PositiveSemidefiniteCone } = false
MathOptInterface.canmodifyconstraint(m::MosekModel, c::MathOptInterface.ConstraintReference{F,D}, dom::D) where { F <: LinearFunction, D <: LinearDomain } = haskey(select(m.constrmap,F,D),c.value)
MathOptInterface.canmodifyconstraint(m::MosekModel, c::MathOptInterface.ConstraintReference{F,D}, func::F) where { F <: Union{MathOptInterface.SingleVariable,MathOptInterface.ScalarAffineFunction}, D <: ScalarLinearDomain } = haskey(select(m.constrmap,F,D),c.value)
MathOptInterface.canmodifyconstraint(m::MosekModel, c::MathOptInterface.ConstraintReference{F,D}, change::MathOptInterface.AbstractFunctionModification) where { F <: AffineFunction, D <: MathOptInterface.AbstractSet } = haskey(select(m.constrmap,F,D),c.value)

chgbound(bl::Float64,bu::Float64,k::Float64,dom :: MathOptInterface.LessThan{Float64})    = bl,dom.upper-k
chgbound(bl::Float64,bu::Float64,k::Float64,dom :: MathOptInterface.GreaterThan{Float64}) = dom.lower-k,bu
chgbound(bl::Float64,bu::Float64,k::Float64,dom :: MathOptInterface.EqualTo{Float64})     = dom.value-k,dom.value-k
chgbound(bl::Float64,bu::Float64,k::Float64,dom :: MathOptInterface.Interval{Float64})    = dom.lower-k,dom.upper-k

function MathOptInterface.modifyconstraint!(m::MosekModel,
                                            xcref::MathOptInterface.ConstraintReference{F,D},
                                            dom::D) where { F    <: MathOptInterface.SingleVariable,
                                                            D    <: ScalarLinearDomain }
    xcid = ref2id(xcref)
    j = m.xc_idxs[getindexes(m.xc_block,xcid)[1]]
    bk,bl,bu = getvarbound(m.task,j)
    bl,bu = chgbound(bl,bu,0.0,dom)
    putvarbound(m.task,j,bk,bl,bu)
end

function MathOptInterface.modifyconstraint!(m::MosekModel,
                                            cref::MathOptInterface.ConstraintReference{F,D},
                                            dom::D) where { F    <: MathOptInterface.ScalarAffineFunction,
                                                            D    <: ScalarLinearDomain }
    cid = ref2id(cref)
    i = getindexes(m.c_block,cid)[1] # since we are in a scalar domain
    bk,bl,bu = getconbound(m.task,i)
    bl,bu = chgbound(bl,bu,0.0,dom)
    putconbound(m.task,i,bk,bl,bu)
end

function MathOptInterface.modifyconstraint!(m   ::MosekModel,
                                            c   ::MathOptInterface.ConstraintReference{MathOptInterface.ScalarAffineFunction{Float64},D},
                                            func::MathOptInterface.ScalarConstantChange{Float64}) where {D <: MathOptInterface.AbstractSet}

    cid = ref2id(c)

    i = getindexes(m.c_block, cid)[1]
    bk,bl,bu = getconbound(m.task,i)
    bl += m.c_constant[i] - func.new_constant
    bu += m.c_constant[i] - func.new_constant
    m.c_constant[i] = func.new_constant
    putconbound(m.task,i,bk,bl,bu)
end

function MathOptInterface.modifyconstraint!(m   ::MosekModel,
                                            c   ::MathOptInterface.ConstraintReference{MathOptInterface.ScalarAffineFunction{Float64},D},
                                            func::MathOptInterface.ScalarCoefficientChange{Float64}) where {D <: MathOptInterface.AbstractSet}
    cid = ref2id(c)
    xid = ref2id(func.variable)

    i = getindexes(m.c_block,cid)[1]
    j = getindexes(m.x_block,xid)[1]

    putaij(m.task,i,j,func.new_coefficient)
end

function MathOptInterface.modifyconstraint!(m::MosekModel,
                                            c::MathOptInterface.ConstraintReference{MathOptInterface.VectorAffineFunction{Float64},D},
                                            func::MathOptInterface.VectorConstantChange{Float64}) where {D <: MathOptInterface.AbstractSet}
    cid = ref2id(c)
    assert(cid > 0)

    subi = getindexes(m.c_block, cid)
    bk = Vector{Int32}(length(subi))
    bl = Vector{Float64}(length(subi))
    bu = Vector{Float64}(length(subi))

    bk,bl,bu = getconboundlist(m.task,convert(Vector{Int32},subi))
    bl += m.c_constant[subi] - func.new_constant
    bu += m.c_constant[subi] - func.new_constant
    m.c_constant[subi] = func.new_constant

    putconboundlist(m.task,convert(Vector{Int32},subi),bk,bl,bu)
end

function MathOptInterface.modifyconstraint!(m::MosekModel,
                                            c::MathOptInterface.ConstraintReference{MathOptInterface.VectorAffineFunction{Float64},D},
                                            func::MathOptInterface.MultirowChange{Float64}) where {D <: MathOptInterface.AbstractSet}
    cid = ref2id(c)
    assert(cid > 0)

    subi = getindexes(m.c_block, cid)[func.rows]
    xid = ref2id(func.variable)
    j = getindexes(m.x_block,xid)[1]

    putaijlist(m.task,convert(Vector{Int32},subi),fill(j,length(subi)),fund.new_coefficients)
end

MathOptInterface.cantransformconstraint(m::MosekModel, c::MathOptInterface.ConstraintReference{F,D1}, newdom::D2) where {F <: VariableFunction, D1, D2 } = false
MathOptInterface.cantransformconstraint(m::MosekModel, c::MathOptInterface.ConstraintReference{MathOptInterface.VectorAffineFunction,D1}, newdom::D2) where {D1 <: VectorLinearDomain, D2 <: VectorLinearDomain} = false
function MathOptInterface.cantransformconstraint(m::MosekModel,
                                                 cref::MathOptInterface.ConstraintReference{MathOptInterface.ScalarAffineFunction{Float64},D1},
                                                 newdom::D2) where {D1 <: ScalarLinearDomain,
                                                                    D2 <: ScalarLinearDomain}
    haskey(select(m.constrmap,MathOptInterface.ScalarAffineFunction{Float64},D1),cref.value)
end

function MathOptInterface.transformconstraint!(m::MosekModel,
                                               cref::MathOptInterface.ConstraintReference{F,D},
                                               newdom::D) where {F  <: MathOptInterface.AbstractFunction, D <: MathOptInterface.AbstractSet}
    MathOptInterface.modifyconstraint(m,cref,newdom)
    cref
end

function MathOptInterface.transformconstraint!(m::MosekModel,
                                               cref::MathOptInterface.ConstraintReference{MathOptInterface.ScalarAffineFunction{Float64},D1},
                                               newdom::D2) where {D1 <: ScalarLinearDomain,
                                                                  D2 <: ScalarLinearDomain}
    const F = MathOptInterface.ScalarAffineFunction{Float64}

    cid = ref2id(cref)

    subi = getindexes(m.c_block,cid)

    addbound!(m,cid,subi,m.c_constant[subi], newdom)

    newcref = MathOptInterface.ConstraintReference{F,D2}(UInt64(cid) << 1)
    delete!(select(m.constrmap,F,D1), cref.value)
    select(m.constrmap,F, D2)[newcref.value] = cid
    newcref
end

function MathOptInterface.transformconstraint!(m::MosekModel,
                                               cref::MathOptInterface.ConstraintReference{MathOptInterface.VectorAffineFunction{Float64},D1},
                                               newdom::D2) where {D1 <: VectorLinearDomain,
                                                                  D2 <: VectorLinearDomain}
    const F = MathOptInterface.VectorAffineFunction{Float64}

    cid = ref2id(cref)

    subi = getindexes(m.c_block,cid)

    addbound!(m,cid,subi,m.c_constant[subi], newdom)

    newcref = MathOptInterface.ConstraintReference{F,D2}(UInt64(cid) << 1)
    delete!(select(m.constrmap,F,D1), cref.value)
    select(m.constrmap,F,D2)[newcref.value] = cid
    newcref
end

#MathOptInterface.transformconstraint!(m::MosekModel, c::MathOptInterface.ConstraintReference{F,D1}, newdom::D2) where {F <: MathOptInterface.VectorAffineFunction , D1 <: VectorLinearDomain, D2 <: VectorLinearDomain} = false

################################################################################
##  DELETE #####################################################################
################################################################################

MathOptInterface.candelete(
    m   ::MosekModel,
    cref::MathOptInterface.ConstraintReference{F,D}) where {F <: Union{MathOptInterface.ScalarAffineFunction,
                                                                       MathOptInterface.VectorAffineFunction,
                                                                       MathOptInterface.SingleVariable,
                                                                       MathOptInterface.VectorOfVariables},
                                                            D <: Union{MathOptInterface.LessThan,
                                                                       MathOptInterface.GreaterThan,
                                                                       MathOptInterface.EqualTo,
                                                                       MathOptInterface.Interval,
                                                                       MathOptInterface.Zeros,
                                                                       MathOptInterface.Nonpositives,
                                                                       MathOptInterface.Nonnegatives,
                                                                       MathOptInterface.Reals}} = MathOptInterface.isvalid(m,cref)

MathOptInterface.candelete(
    m   ::MosekModel,
    cref::MathOptInterface.ConstraintReference{F,D}) where {F <: MathOptInterface.AbstractFunction,
                                                            D <: Union{MathOptInterface.SecondOrderCone,
                                                                       MathOptInterface.RotatedSecondOrderCone,
                                                                       MathOptInterface.ExponentialCone,
                                                                       MathOptInterface.DualExponentialCone,
                                                                       MathOptInterface.PowerCone,
                                                                       MathOptInterface.DualPowerCone}} = false

MathOptInterface.candelete(
    m   ::MosekModel,
    cref::MathOptInterface.ConstraintReference{F,D}) where {F <: Union{MathOptInterface.SingleVariable,
                                                                       MathOptInterface.VectorOfVariables},
                                                            D <: Union{MathOptInterface.SecondOrderCone,
                                                                       MathOptInterface.RotatedSecondOrderCone,
                                                                       MathOptInterface.ExponentialCone,
                                                                       MathOptInterface.DualExponentialCone,
                                                                       MathOptInterface.PowerCone,
                                                                       MathOptInterface.DualPowerCone}} = true

function Base.delete!(
    m::MosekModel,
    cref::MathOptInterface.ConstraintReference{F,D}) where {F <: Union{MathOptInterface.ScalarAffineFunction,
                                                                       MathOptInterface.VectorAffineFunction},
                                                            D <: Union{MathOptInterface.LessThan,
                                                                       MathOptInterface.GreaterThan,
                                                                       MathOptInterface.EqualTo,
                                                                       MathOptInterface.Interval,
                                                                       MathOptInterface.Zeros,
                                                                       MathOptInterface.Nonpositives,
                                                                       MathOptInterface.Nonnegatives,
                                                                       MathOptInterface.Reals}}

    delete!(select(m.constrmap, F, D), cref.value)

    cid = ref2id(cref)
    subi = getindexes(m.c_block,cid)

    n = length(subi)
    subi_i32 = convert(Vector{Int32},subi)
    ptr = fill(Int64(0),n)
    putarowlist(m.task,subi_i32,ptr,ptr,Int32[],Float64[])
    b = fill(0.0,n)
    putconboundlist(m.task,subi_i32,fill(MSK_BK_FX,n),b,b)

    m.c_constant[subi] = 0.0
    deleteblock(m.c_block,cid)
end

function Base.delete!(
    m::MosekModel,
    cref::MathOptInterface.ConstraintReference{F,D}) where {F <: Union{MathOptInterface.VectorAffineFunction},
                                                            D <: Union{MathOptInterface.SecondOrderCone,
                                                                       MathOptInterface.RotatedSecondOrderCone,
                                                                       MathOptInterface.ExponentialCone,
                                                                       MathOptInterface.DualExponentialCone,
                                                                       MathOptInterface.PowerCone,
                                                                       MathOptInterface.DualPowerCone}}

    delete!(select(m.constrmap, F, D), cref.value)
    xcid = ref2id(cref)
    sub = getindexes(m.xc_block,xcid)

    subj = [ getindexes(m.x_block,i)[1] for i in sub ]
    N = length(subj)
    m.x_boundflags[subj] .&= ~m.xc_bounds[xcid]
    asgn,coneidx = getconenameindex(m.task,"$(m.xc_coneid[xcid])")
    m.xc_coneid[xcid] = 0
    removecone(m.task,coneidx)

    m.x_numxc[subj]  -= 1
    m.xc_idxs[sub]    = 0
    m.xc_bounds[xcid] = 0

    deleteblock(m.xc_block,xcid)
end

function Base.delete!(
    m::MosekModel,
    cref::MathOptInterface.ConstraintReference{F,D}) where {F <: Union{MathOptInterface.SingleVariable,
                                                                       MathOptInterface.VectorOfVariables},
                                                            D <: Union{MathOptInterface.LessThan,
                                                                       MathOptInterface.GreaterThan,
                                                                       MathOptInterface.EqualTo,
                                                                       MathOptInterface.Interval,
                                                                       MathOptInterface.Zeros,
                                                                       MathOptInterface.Nonpositives,
                                                                       MathOptInterface.Nonnegatives,
                                                                       MathOptInterface.Reals,
                                                                       MathOptInterface.ZeroOne,
                                                                       MathOptInterface.Integer}}
    delete!(select(m.constrmap, F, D), cref.value)

    xcid = ref2id(cref)
    sub = getindexes(m.xc_block,xcid)

    subj = [ getindexes(m.x_block,i)[1] for i in sub ]
    N = length(subj)

    m.x_boundflags[subj] .&= ~m.xc_bounds[xcid]
    if m.xc_bounds[xcid] & boundflag_int != 0
        for i in 1:length(subj)
            putvartype(m.task,subj[i],MSK_VAR_TYPE_CONT)
        end
    end

    if m.xc_bounds[xcid] & boundflag_lower != 0 && m.xc_bounds[xcid] & boundflag_upper != 0
        bnd = fill(0.0, length(N))
        putvarboundlist(m.task,convert(Vector{Int32},subj),fill(MSK_BK_FR,N),bnd,bnd)
    elseif m.xc_bounds[xcid] & boundflag_lower != 0
        bnd = fill(0.0, length(N))
        bk,bl,bu = getvarboundlist(m.task, convert(Vector{Int32},subj))

        for i in 1:N
            if MSK_BK_RA == bk[i] || MSK_BK_FX == bk[i]
                bk[i] = MSK_BK_UP
            else
                bk[i] = MSK_BK_FR
            end
        end

        putvarboundlist(m.task,convert(Vector{Int32},subj),bk,bl,bu)
    elseif m.xc_bounds[xcid] & boundflag_upper != 0
        bnd = fill(0.0, length(N))
        bk,bl,bu = getvarboundlist(m.task, subj)

        for i in 1:N
            if MSK_BK_RA == bk[i] || MSK_BK_FX == bk[i]
                bk[i] = MSK_BK_LO
            else
                bk[i] = MSK_BK_FR
            end
        end

        putvarboundlist(m.task,convert(Vector{Int32},subj),bk,bl,bu)
    else
        assert(false)
        # should not happen
    end

    m.x_numxc[subj] -= 1
    m.xc_idxs[sub] = 0
    m.xc_bounds[xcid] = 0

    deleteblock(m.xc_block,xcid)
end

################################################################################
################################################################################
################################################################################

function allocateconstraints(m :: MosekModel,
                             N :: Int)
    numcon = getnumcon(m.task)
    alloced = ensurefree(m.c_block,N)
    id = newblock(m.c_block,N)

    M = numblocks(m.c_block) - length(m.c_block_slack)
    if alloced > 0
        appendcons(m.task, alloced)
        append!(m.c_constant, zeros(Float64,alloced))
    end
    if M > 0
        append!(m.c_block_slack, zeros(Float64,M))
        append!(m.c_coneid, zeros(Float64,M))
    end
    id
end

function allocatevarconstraints(m :: MosekModel,
                                N :: Int)
    nalloc = ensurefree(m.xc_block,N)
    id = newblock(m.xc_block,N)

    M = numblocks(m.xc_block) - length(m.xc_bounds)
    if M > 0
        append!(m.xc_bounds,zeros(Float64,M))
        append!(m.xc_coneid,zeros(Float64,M))
    end
    if nalloc > 0
        append!(m.xc_idxs, zeros(Float64,nalloc))
    end

    id
end

function allocatevariable(m :: MosekModel,N :: Int)
    numvar = getnumvar(m.task)
    alloced = ensurefree(m.x_block,N)
    if alloced > 0
        appendvars(m.task, length(m.x_block) - numvar)
        append!(m.x_boundflags, zeros(Int,length(m.x_block) - numvar))
        append!(m.x_numxc, zeros(Int,length(m.x_block) - numvar))
    end
    newblock(m.x_block,N)
end

MathOptInterface.isvalid(m::MosekModel, ref::MathOptInterface.ConstraintReference{F,D}) where { F,D } = haskey(select(m.constrmap,F,D),ref.value)
MathOptInterface.isvalid(m::MosekModel, ref::MathOptInterface.VariableReference) = allocated(m.x_block,ref2id(ref))

function getvarboundlist(t::Mosek.Task, subj :: Vector{Int32})
    n = length(subj)
    bk = Vector{Boundkey}(n)
    bl = Vector{Float64}(n)
    bu = Vector{Float64}(n)
    for i in 1:n
        bki,bli,bui = getvarbound(t,subj[i])
        bk[i] = bki
        bl[i] = bli
        bu[i] = bui
    end
    bk,bl,bu
end

function getconboundlist(t::Mosek.Task, subj :: Vector{Int32})
    n = length(subj)
    bk = Vector{Boundkey}(n)
    bl = Vector{Float64}(n)
    bu = Vector{Float64}(n)
    for i in 1:n
        bki,bli,bui = getconbound(t,subj[i])
        bk[i] = bki
        bl[i] = bli
        bu[i] = bui
    end
    bk,bl,bu
end


text:
export mu_lambda

# Executes mu-lambda evolution relative to goal starting with a random population of chromosomes
# mu is the number of parents and lambda is the number of children in each generation
# mu and lambda are components of the parameters
# gens is the maximum number of generations run
function mu_lambda(p::Parameters, goal::Goal, gens::Integer)
    mu = p.mu
    lambda = p.lambda
    funcs = p.funcs
    fitfunc = p.fitfunc
    perfect = p.targetfitness

    pop = [random_chromosome(p) for _ in 1:(mu+lambda) ]
    fit = [fitness(c,goal,fitfunc) for c in pop ]
    perm = sortperm(fit, rev=true)

    for t = 1:gens
        if fit[perm[1]] == perfect
            return (pop[perm[1]], t)
        end
        mutpop = [mutate(x,funcs) for x in pop[perm][mu+1:mu+lambda]]
        mutfit = [fitness(x, goal, fitfunc) for x in mutpop]
        newpop = vcat(pop[perm][1:mu], mutpop)
        newfit = vcat(fit[perm][1:mu], mutfit)
        pop = newpop
        fit = newfit

        perm = sortperm(fit, rev=true)
    end

    return (pop[perm][1], gens)
end


text:
module GoogleSheetsCSVExporter

import HTTP

struct Sheet
	documentID
	gid
end

function parseURI(uri::HTTP.URI)
	documentID = split(uri.path, "/")[4]
	fragmentParams = uri.fragment |> HTTP.queryparams
	Sheet(documentID, fragmentParams["gid"])
end

function exportURI(sheet::Sheet)
	"https://docs.google.com/spreadsheets/d/$(sheet.documentID)/export?format=csv&gid=$(sheet.gid)"
end

function openURI(uri::String)
	let str
		HTTP.open("GET", uri) do http
			str = read(http, String)
		end
		str
	end
end

function fromURI(uri::String)
	HTTP.URI(uri) |> parseURI |> exportURI |> openURI |> IOBuffer
end

end # GoogleSheetsCSVExporter


text:
function get_msgtype(::Message{m}) :: LDAPMessageType where {m}
    m
end

function destroy!(chain::MessageChain)
    if chain.handle != C_NULL
        type_of_last_msg_in_chain = LDAPMessageType(ldap_msgfree(chain.handle))
        #@info("Freed msg iter with last msg of type $type_of_last_msg_in_chain")
        chain.handle = C_NULL
    end
    nothing
end

function search(ldap::LDAPConnection, base::AbstractString, scope::LDAPScope;
        filter::Union{Nothing, AbstractString}=nothing, attr_desc_only::Bool=false, size_limit::Integer=-1) :: MessageChain

    result = Ref{Ptr{Cvoid}}()
    err = ldap_search_ext_s(ldap.handle, base, Cint(scope), filter, attr_desc_only, size_limit, result)
    msg_chain = MessageChain(result[], ldap) # assures result is freed regardless of the search err result
    error_check(err)
    return msg_chain
end

count_messages(chain::MessageChain) = ldap_count_messages(chain.ldap.handle, chain.handle)
count_entries(chain::MessageChain) = ldap_count_entries(chain.ldap.handle, chain.handle)
count_references(chain::MessageChain) = ldap_count_references(chain.ldap.handle, chain.handle)
count_results(chain::MessageChain) = count_messages(chain) - count_entries(chain) - count_references(chain)

each_message(chain::MessageChain) = MessageIterator{AbstractMessage}(chain)
each_entry(chain::MessageChain) = MessageIterator{Entry}(chain)
each_reference(chain::MessageChain) = MessageIterator{Reference}(chain)

for triplet in [ (:AbstractMessage, :ldap_first_message, :ldap_next_message),
                 (:Entry,           :ldap_first_entry,   :ldap_next_entry),
                 (:Reference,       :ldap_first_reference, :ldap_next_reference)
               ]
    @eval begin
        function Base.iterate(itr::MessageIterator{$(triplet[1])})
            first_msg_handle = ($(triplet[2]))(itr.chain.ldap.handle, itr.chain.handle)
            if first_msg_handle == C_NULL
                return nothing
            else
                return Message(first_msg_handle, itr.chain), first_msg_handle
            end
        end

        function Base.iterate(itr::MessageIterator{$(triplet[1])}, state::Ptr{Cvoid})
            next_msg_handle = ($(triplet[3]))(itr.chain.ldap.handle, state)
            if next_msg_handle == C_NULL
                return nothing
            else
                return Message(next_msg_handle, itr.chain), next_msg_handle
            end
        end
    end
end

Base.length(itr::MessageIterator{AbstractMessage}) = count_messages(itr.chain)
Base.length(itr::MessageIterator{Entry}) = count_entries(itr.chain)
Base.length(itr::MessageIterator{Reference}) = count_references(itr.chain)

function get_dn(entry::Entry) :: String
    dn_cstring = ldap_get_dn(entry.chain.ldap.handle, entry.handle)
    result = unsafe_string(dn_cstring)
    ldap_memfree(dn_cstring)
    return result
end

Base.show(io::IO, entry::Entry) = show(io, "Entry($(get_dn(entry)))")

function parse_controls_ptr_vec(controls_ptr_vec::Ptr{Ptr{LDAPControl}}, free_controls::Bool) :: Vector{Control}
    result = Vector{Control}()

    if controls_ptr_vec == C_NULL
        return result
    end

    local i = 1
    while true
        control_ptr = unsafe_load(controls_ptr_vec, i)

        if control_ptr == C_NULL
            break
        end

        ldap_control = unsafe_load(control_ptr)

        oid = unsafe_string(ldap_control.oid)
        val = unsafe_string(ldap_control.value.val, ldap_control.value.len)
        iscritical = Char(ldap_control.iscritical)
        push!(result, Control(oid, val, iscritical))

        i += 1
    end

    free_controls && ldap_controls_free(controls_ptr_vec)
    return result
end

function get_controls(entry::Entry) :: Vector{Control}
    controls_ptr_vec_ref = Ref{Ptr{Ptr{LDAPControl}}}()
    err = ldap_get_entry_controls(entry.chain.ldap.handle, entry.handle, controls_ptr_vec_ref)
    error_check(err)
    return parse_controls_ptr_vec(controls_ptr_vec_ref[], true)
end

function parse_result(result::Result) :: ParsedResult
    err_ref = Ref{Cint}()
    match_str_ref = Ref{Ptr{UInt8}}()
    err_msg_ref = Ref{Ptr{UInt8}}()
    referralsp_ref = Ref{Ptr{Ptr{UInt8}}}()
    controls_ref = Ref{Ptr{Ptr{LDAPControl}}}()

    err = ldap_parse_result(result.chain.ldap.handle, result.handle, err_ref, match_str_ref, err_msg_ref, referralsp_ref, controls_ref, false)
    error_check(err)

    local match_str::String = ""
    if match_str_ref[] != C_NULL
        match_str = unsafe_string(match_str_ref[])
        ldap_memfree(match_str_ref[])
    end

    local error_message::String = ""
    if err_msg_ref[] != C_NULL
        error_message = unsafe_string(err_msg_ref[])
        ldap_memfree(err_msg_ref[])
    end

    referrals = parse_null_terminated_string_vector(referralsp_ref[])
    if referralsp_ref[] != C_NULL
        ldap_memvfree(referralsp_ref[])
    end

    controls_vec = parse_controls_ptr_vec(controls_ref[], true)

    return ParsedResult(
            err_ref[],
            match_str,
            error_message,
            referrals,
            controls_vec
        )
end


text:
# http://www.gdal.org/classGDALMajorObject.html

abstract type GDALMajorObjectH end

abstract type GDALDatasetH <: GDALMajorObjectH end
abstract type GDALDriverH <: GDALMajorObjectH end
abstract type GDALRasterBandH <: GDALMajorObjectH end
abstract type OGRLayerH <: GDALMajorObjectH end

abstract type OGRDataSourceH <: GDALDatasetH end
abstract type VRTDatasetH <: GDALDatasetH end

abstract type OGRSFDriverH <: GDALDriverH end
abstract type VRTDriverH <: GDALDriverH end

abstract type GDALAsyncReaderH end
abstract type GDALColorTableH end
abstract type GDALContourGeneratorH end
abstract type GDALContourWriter end
abstract type GDALDerivedPixelFunc end
abstract type GDALMaskFunc end
abstract type GDALProgressFunc end
abstract type GDALRasterAttributeTableH end
abstract type GDALTransformerFunc end
abstract type GDALWarpOperationH end
abstract type GDALGridContext end

abstract type OGRCoordinateTransformationH end
abstract type OGRFeatureDefnH end
abstract type OGRFeatureH end
abstract type OGRFieldDefnH end
abstract type OGRGeometryH end
abstract type OGRGeomFieldDefnH end
abstract type OGRSpatialReferenceH end
abstract type OGRStyleMgrH end
abstract type OGRStyleTableH end
abstract type OGRStyleToolH end

abstract type CPLErrorHandler end
abstract type FilterFunc4ValuesType end
abstract type FilterFuncType end

abstract type VRTImageReadFunc end
abstract type VRTSourceH end
abstract type VRTSimpleSourceH end
abstract type VRTAveragedSourceH end
abstract type VRTComplexSourceH end
abstract type VRTFilteredSourceH end
abstract type VRTKernelFilteredSourceH end
abstract type VRTAverageFilteredSourceH end
abstract type VRTFuncSourceH end
abstract type VRTWarpedDatasetH end
abstract type VRTRasterBandH end
abstract type VRTSourcedRasterBandH end
abstract type VRTWarpedRasterBandH end
abstract type VRTDerivedRasterBandH end
abstract type VRTRawRasterBandH end

abstract type CPLFileFinder end
abstract type CPLVirtualMemCachePageCbk end
abstract type CPLVirtualMemFreeUserData end
abstract type CPLVirtualMemUnCachePageCbk end
abstract type VSIStatBufL end
abstract type _CPLXMLNode end
abstract type CPLVirtualMem end

# GDAL utils options
abstract type GDALBuildVRTOptions end
abstract type GDALBuildVRTOptionsForBinary end
abstract type GDALDEMProcessingOptions end
abstract type GDALDEMProcessingOptionsForBinary end
abstract type GDALGridOptions end
abstract type GDALGridOptionsForBinary end
abstract type GDALInfoOptions end
abstract type GDALInfoOptionsForBinary end
abstract type GDALNearblackOptions end
abstract type GDALNearblackOptionsForBinary end
abstract type GDALRasterizeOptions end
abstract type GDALRasterizeOptionsForBinary end
abstract type GDALTranslateOptions end
abstract type GDALTranslateOptionsForBinary end
abstract type GDALVectorTranslateOptions end
abstract type GDALVectorTranslateOptionsForBinary end
abstract type GDALWarpAppOptions end
abstract type GDALWarpAppOptionsForBinary end


text:
using Test

@testset "Exercise 3" begin

# Test for section on RV model
@testset "rv_model_const" begin
   @test isdefined(Main, :rv_model_const)
   @test rv_model_const(42)(1.0) ≈ 42.0
   @test rv_model_const(pi)(17.0) ≈ 3.14159  atol=1e-5
end; # rv_model_const

@testset "rv_model_circ" begin
   @test isdefined(Main, :rv_model_circ)
   # INSERT YOUR TESTS HERE
end; # rv_model_cric

# Test for section on probability
@testset "gauss_std_normal" begin
   @test isdefined(Main, :gauss_std_normal)
   @test gauss_std_normal(0) ≈ 1/sqrt(2pi)
   @test gauss_std_normal(3) ≈ exp(-4.5)/sqrt(2pi)
end; # gauss_std_normal

@testset "likelihood" begin
   @test isdefined(Main, :likelihood)
   @test likelihood(rv_model_const(-1.0),ones(10),ones(10),ones(10)) ≈ exp(-10*4/2)*(2pi)^(-5)
end; # likelihood

@testset "simulate_rvs_and_calc_likelihood" begin
   @test isdefined(Main, :simulate_rvs_and_calc_likelihood)
   @test exp(-120/2)*(2pi)^(-50) < simulate_rvs_and_calc_likelihood(100) < exp(-80/2)*(2pi)^(-50)
   @test simulate_rvs_and_calc_likelihood(600) == 0.0
end; # simulate_rvs_and_calc_likelihoodlikelihood

@testset "log_likelihood" begin
   @test isdefined(Main, :log_likelihood)
   @test log_likelihood(rv_model_const(-1.0),ones(10),ones(10),ones(10)) ≈ -10*4/2-10/2*log(2pi)
end; # log_likelihood

@testset "simulate_rvs_and_calc_log_likelihood" begin
   @test isdefined(Main, :simulate_rvs_and_calc_log_likelihood)
   @test -120/2-50*log(2pi) < simulate_rvs_and_calc_log_likelihood(100) < -80/2-50*log(2pi)
   @test -650/2-300*log(2pi) < simulate_rvs_and_calc_log_likelihood(600) < -550/2-300*log(2pi)
end; # simulate_rvs_and_calc_log_likelihood

end; # Exercise 3


text:
type FlatBreeder <: Breeder
  selection::Selection
  mutation::Mutation
  crossover::Crossover

  FlatBreeder(s::Selection, m::Mutation, c::Crossover) =
    new(s, m, c)
end

type FlatBreederDefinition <: BreederDefinition
  selection::SelectionDefinition
  mutation::MutationDefinition
  crossover::CrossoverDefinition

  FlatBreederDefinition() = new()
  FlatBreederDefinition(s::SelectionDefinition, m::MutationDefinition, c::CrossoverDefinition) =
    new(s, m, c)
end

"""
Composes a flat breeder from its definition.
"""
compose!(d::FlatBreederDefinition, sp::Species) =
  FlatBreeder(compose!(d.selection),
              compose!(d.mutation, sp),
              compose!(d.crossover, sp))

"""
TODO: Document flat breeders.
"""
function flat(f::Function)
  def = FlatBreederDefinition()
  f(def)
  def
end

"""
Performs the breeding process for each of the demes within the population
contained by a given state.
"""
function breed!{F}(b::FlatBreeder, sp::Species, buffer::IndividualCollection{F}, members::IndividualCollection{F}, n::Int)
  n_crossover = num_required(b.crossover, n)
  select(b.selection, buffer, sp, members, n_crossover)
  operate!(b.crossover, buffer, n_crossover)
  operate!(b.mutation, buffer, n)
  return
end


text:
function split!(quad::Quadtree, index::Int64)
    global level

    if index == 1
        output = [
                    ((quad.x1 + quad.x0) / 2),
                    (quad.y0),
                    (quad.x1 - quad.x0) / 2,
                    (quad.y1 - quad.y0) / 2
                ]*1.0
    end
    if index == 2
        output = [
                    (quad.x0),
                    (quad.y0),
                    (quad.x1 - quad.x0) / 2,
                    (quad.y1 - quad.y0) / 2
                ]*1.0
    end
    if index == 3
        output = [
                    (quad.x0),
                    (quad.y1 + quad.y0) / 2,
                    (quad.x1 - quad.x0) / 2,
                    (quad.y1 - quad.y0) / 2
                ]*1.0
    end
    if index == 4
        output = [
                    (quad.x1 + quad.x0) / 2,
                    (quad.y1 + quad.y0) / 2,
                    (quad.x1 - quad.x0) / 2,
                    (quad.y1 - quad.y0) / 2
                ]*1.0
    end

    if index == -1
        return
    end

    return output
end


text:
abstract type GANLoss end

# Loss that includes the generator calls out to loss that just takes in generated values
Lᴰ(t::GANLoss, G, D, z, x; wD = 1f0, yG = (), yD = ()) = Lᴰ(t, D, x, G(z, yG...), wD = wD, yG = yG, yD = yD)

## Binary Cross Entropy Loss
struct GAN_BCELoss <: GANLoss end
const LBCE = Flux.Losses.logitbinarycrossentropy
Lᴰ(t::GAN_BCELoss, D, x, xtilde; wD = 1f0, yG = (), yD = ()) = LBCE(D(vcat(x, yD...)), wD) + LBCE(D(vcat(xtilde, yG...)), 0f0)
Lᴳ(t::GAN_BCELoss, G, D, z; yG = ()) = LBCE(D(vcat(G(vcat(z, yG...)), yG...)), 1f0)

## Least Squares Loss
struct GAN_LSLoss <: GANLoss end
Lᴰ(t::GAN_LSLoss, D, x, xtilde; wD = 1f0, yG = (), yD = ()) = Flux.mse(sigmoid.(D(vcat(x, yD...))), wD) + Flux.mse(sigmoid.(D(vcat(xtilde, yG...))), 0f0)
Lᴳ(t::GAN_LSLoss, G, D, z; yG = ()) = Flux.mse(D(vcat(G(vcat(z, yG...)), yG...)), 1f0)

## Hinge Loss
struct GAN_HingeLoss <: GANLoss end
Lᴰ(t::GAN_HingeLoss, D, x, xtilde; wD = 1f0, yG = (), yD = ()) = mean(relu.(1f0 .+ (-2f0 .* wD .+ 1f0) .* D(vcat(x, yD...)))) + mean(relu.(1f0 .+ D(vcat(xtilde, yG...))))
Lᴳ(t::GAN_HingeLoss, G, D, z; yG = ()) = -mean(D(vcat(G(vcat(z, yG...)), yG...)))

## WLoss
struct GAN_WLoss <: GANLoss end
Lᴰ(t::GAN_WLoss, D, x, xtilde; wD = 1f0, yG = (), yD = ()) = mean(D(vcat(xtilde, yG...))) -  mean((wD .* 2f0 .- 1f0) .* D(vcat(x, yD...)))
Lᴳ(t::GAN_WLoss, G, D, z; yG = ()) = -mean(D(vcat(G(vcat(z, yG...)), yG...)))

## WGAN-GP
@with_kw struct GAN_WLossGP <: GANLoss
	λ::Float32 = 10f0 # gp parameter
	noise_range::Float32 = 1f0 # Range of the ϵ parameter that picks a point to evalutate the gradient (smaller is closer to "real" data)
end

function Lᴰ(t::GAN_WLossGP, D, x, xtilde; wD = 1f0, yG = (), yD = ())
	loss = mean(D(vcat(xtilde, yG...))) -  mean((wD .* 2f0 .- 1f0) .* D(vcat(x, yD...)))

	# Compute gradient penalty
	loss += (length(yG) == 1) ?  t.λ*gradient_penalty(D, vcat(x, yD...), vcat(xtilde, yG...)) : t.λ*gradient_penalty(D, x, xtilde)
	loss
end

Lᴳ(t::GAN_WLossGP, G, D, z; yG = ()) = -mean(D(vcat(G(vcat(z, yG...)), yG...)))

GANLosses = Dict("BCE" => GAN_BCELoss(), "LS" => GAN_LSLoss(), "Hinge" => GAN_HingeLoss(), "W" => GAN_WLoss(), "WGP" => GAN_WLossGP())


text:
using Test
using ImagineInterface
import Unitful: μm, s

sample_rate=10000s^-1
ocpi2 = rigtemplate("ocpi-2"; sample_rate = sample_rate)
pmin=100.0μm; pmax=180.0μm
z_spacing = 9.99μm #The space between slices in the z-stack.
stk_rate=5s^-1
reset_time = 0.06s
z_pad = 0.0μm
t_fullstk=1/stk_rate
t_imstk=t_fullstk-reset_time
n_stk_frames=round(Int,(pmax-pmin-2*z_pad)/z_spacing)
frac_realscan=z_spacing*n_stk_frames/(pmax-pmin) # fraction of "meaningful scan" in the forward scanning phase
exp_time=(t_imstk*frac_realscan)/n_stk_frames*0.90
flash_frac=2
stk = gen_unidirectional_stack(pmin, pmax, z_spacing, t_imstk, reset_time, exp_time, sample_rate, flash_frac; z_pad = z_pad)
@test size(stk["positioner"]) == size(stk["camera"]) == size(stk["laser"])

max_forks_repo_name: HolyLab/ImagineInterface


text:
using OndaEDF
using Documenter

makedocs(modules=[OndaEDF],
         sitename="OndaEDF",
         authors="Beacon Biosignals and other contributors",
         pages=["API Documentation" => "index.md"])

deploydocs(repo="github.com/beacon-biosignals/OndaEDF.jl.git",
           push_preview=true)


text:
using Celeste
using Celeste.Model
using Celeste.SensitiveFloats

import Celeste.PSF: get_psf_at_point, fit_psf,
       get_psf_transform, initialize_psf_params, transform_psf_params!,
       unwrap_psf_params, wrap_psf_params,
       unconstrain_psf_params, constrain_psf_params,
       transform_psf_sensitive_float!,
       PsfOptimizer, fit_raw_psf_for_celeste, trim_psf,
       BivariateNormalDerivatives
import Celeste.Model: eval_psf

using ForwardDiff
using StaticArrays

"""
Evaluate the sum of squared difference between the raw_psf and the psf
represented by psf_params, as well as the derivatives and hessians with
respect to unconstrained parameters.

Returns:
  - A sensitive float for the sum of squared differences.
"""
function evaluate_psf_fit{NumType <: Number}(
    psf_params::Vector{Vector{NumType}}, raw_psf::Matrix{Float64},
    calculate_gradient::Bool)

  K = length(psf_params)
  x_mat = PSF.get_x_matrix_from_psf(raw_psf)

  # TODO: allocate these outside?
  bvn_derivs = BivariateNormalDerivatives{NumType}(NumType)
  log_pdf = SensitiveFloat{NumType}(length(PsfParams), 1, true, true)
  pdf = SensitiveFloat{NumType}(length(PsfParams), 1, true, true)

  pixel_value = SensitiveFloat{NumType}(length(PsfParams), K, true, true)
  squared_error = SensitiveFloat{NumType}(length(PsfParams), K, true, true)

  PSF.evaluate_psf_fit!(
      psf_params, raw_psf, x_mat, bvn_derivs,
      log_pdf, pdf, pixel_value, squared_error, calculate_gradient)

  squared_error
end

function load_raw_psf(; x::Float64=500., y::Float64=500.)
  run_num = 3900
  camcol_num = 6
  field_num = 269
  b = 3

  psf_filename =
    @sprintf("%s/%s/%s/%s/psField-%06d-%d-%04d.fit",
        datadir, run_num, camcol_num, field_num,
                 run_num, camcol_num, field_num)
  psf_fits = FITSIO.FITS(psf_filename)
  raw_psf_comp = SDSSIO.read_psf(psf_fits, band_letters[b])
  close(psf_fits)

  eval_psf(raw_psf_comp, x, y)
end

function test_transform_psf_params()
  K = 2
  psf_params = initialize_psf_params(K, for_test=true)
  psf_params_original = deepcopy(psf_params)
  psf_params_free = deepcopy(psf_params)
  psf_transform = PSF.get_psf_transform(psf_params)

  transform_psf_params!(psf_params, psf_params_free, psf_transform, true)
  transform_psf_params!(psf_params, psf_params_free, psf_transform, false)

  psf_params_free_2 = unconstrain_psf_params(psf_params, psf_transform)
  psf_params_2 = constrain_psf_params(psf_params_free, psf_transform)

  for k=1:K
    @test_approx_eq psf_params[k] psf_params_original[k]
    @test_approx_eq psf_params_free[k] psf_params_free_2[k]
    @test_approx_eq psf_params[k] psf_params_2[k]
  end
end

function test_psf_fit()
  raw_psf = load_raw_psf()

  # Initialize params
  K = 2
  psf_params = initialize_psf_params(K, for_test=true)
  psf_param_vec = wrap_psf_params(psf_params)

  function pixel_value_wrapper_sf{NumType <: Number}(
      psf_param_vec::Vector{NumType}, calculate_gradient::Bool)

    local psf_params = unwrap_psf_params(psf_param_vec)
    bvn_derivs = BivariateNormalDerivatives{NumType}(NumType)
    log_pdf = SensitiveFloat{NumType}(length(PsfParams), 1, true, true)
    pdf = SensitiveFloat{NumType}(length(PsfParams), 1, true, true)

    local pixel_value = SensitiveFloat{NumType}(length(PsfParams), K, true, true)

    local sigma_vec, sig_sf_vec, bvn_vec
    sigma_vec, sig_sf_vec, bvn_vec = PSF.get_sigma_from_params(psf_params)

    # sigma_vec = Array(Matrix{NumType}, K)
    # sig_sf_vec = Array(GalaxySigmaDerivs{NumType}, K)
    #
    # for k = 1:K
    #   sigma_vec[k] = PSF.get_bvn_cov(psf_params[k][psf_ids.e_axis],
    #                                   psf_params[k][psf_ids.e_angle],
    #                                   psf_params[k][psf_ids.e_scale])
    #   sig_sf_vec[k] = GalaxySigmaDerivs(
    #     psf_params[k][psf_ids.e_angle],
    #     psf_params[k][psf_ids.e_axis],
    #     psf_params[k][psf_ids.e_scale], sigma_vec[k], calculate_gradient)
    #
    # end

    clear!(pixel_value)
    PSF.evaluate_psf_pixel_fit!(
        x, psf_params, sigma_vec, sig_sf_vec, bvn_vec,
        bvn_derivs, log_pdf, pdf, pixel_value, calculate_gradient)

    pixel_value
  end

  function pixel_value_wrapper_value{NumType <: Number}(psf_param_vec::Vector{NumType})
    pixel_value_wrapper_sf(psf_param_vec, false).v[]
  end

  x = @SVector [1.0, 2.0]

  sigma_vec = Array(Matrix{Float64}, K)
  for k = 1:K
    sigma_vec[k] = PSF.get_bvn_cov(psf_params[k][psf_ids.e_axis],
                                    psf_params[k][psf_ids.e_angle],
                                    psf_params[k][psf_ids.e_scale])
  end

  println("Testing single pixel value")
  psf_components = PsfComponent[
    PsfComponent(psf_params[k][psf_ids.weight], SVector{2,Float64}(psf_params[k][psf_ids.mu]), SMatrix{2,2,Float64,4}(sigma_vec[k]))
                  for k = 1:K ]

  psf_rendered = get_psf_at_point(psf_components, rows=[ x[1] ], cols=[ x[2] ])[1]
  @test_approx_eq psf_rendered pixel_value_wrapper_value(psf_param_vec)

  pixel_value = deepcopy(pixel_value_wrapper_sf(psf_param_vec, true))

  ad_grad = ForwardDiff.gradient(pixel_value_wrapper_value, psf_param_vec)
  ad_hess = ForwardDiff.hessian(pixel_value_wrapper_value, psf_param_vec)

  @test_approx_eq ad_grad pixel_value.d[:]
  @test_approx_eq ad_hess[:] pixel_value.h[:]

  # Test the whole least squares function.
  println("Testing psf least squares")

  # Fewer pixels for quick testing.  Also, ForwardDiff.hessian runs into strange
  # problems on the whole image.
  keep_pixels = 20:30

  function evaluate_psf_fit_wrapper_sf{NumType <: Number}(
        psf_param_vec::Vector{NumType}, calculate_gradient::Bool)
    local psf_params = unwrap_psf_params(psf_param_vec)
    local squared_error =
      evaluate_psf_fit(psf_params, raw_psf[keep_pixels, keep_pixels], calculate_gradient)
    squared_error
  end

  function evaluate_psf_fit_wrapper_value{NumType <: Number}(psf_param_vec::Vector{NumType})
    local squared_error = evaluate_psf_fit_wrapper_sf(psf_param_vec, false)
    squared_error.v[]
  end

  squared_error = deepcopy(evaluate_psf_fit_wrapper_sf(psf_param_vec, true))

  ad_grad = ForwardDiff.gradient(evaluate_psf_fit_wrapper_value, psf_param_vec)
  ad_hess = ForwardDiff.hessian(evaluate_psf_fit_wrapper_value, psf_param_vec)

  @test_approx_eq ad_grad squared_error.d[:]
  @test_approx_eq ad_hess[:] squared_error.h[:]
end

function test_transform_psf_sensitive_float()
  raw_psf = load_raw_psf()

  K = 2
  psf_params = initialize_psf_params(K, for_test=true)
  psf_transform = PSF.get_psf_transform(psf_params)
  psf_params_free = unconstrain_psf_params(psf_params, psf_transform)
  psf_params_free_vec = wrap_psf_params(psf_params_free)[:]

  # Fewer pixels for quick testing.  Also, ForwardDiff.hessian runs into strange
  # problems on the whole image.
  keep_pixels = 20:30

  function psf_fit_for_optim{NumType <: Number}(
      psf_params_free_vec::Vector{NumType}, calculate_gradient::Bool)

    local sf_free = SensitiveFloat{NumType}(length(PsfParams), K, true, true)
    local psf_params_free = unwrap_psf_params(psf_params_free_vec)
    local psf_params = constrain_psf_params(psf_params_free, psf_transform)
    local sf = evaluate_psf_fit(
      psf_params, raw_psf[keep_pixels, keep_pixels], calculate_gradient)
    transform_psf_sensitive_float!(
      psf_params, psf_transform, sf, sf_free, calculate_gradient)

    sf_free
  end

  function psf_fit_for_optim_val{NumType <: Number}(
      psf_params_free_vec::Vector{NumType})

    psf_fit_for_optim(psf_params_free_vec, false).v[]
  end

  sf_free = deepcopy(psf_fit_for_optim(psf_params_free_vec, true))

  expected_value =
    evaluate_psf_fit(psf_params, raw_psf[keep_pixels, keep_pixels], false).v[]
  ad_grad = ForwardDiff.gradient(psf_fit_for_optim_val, psf_params_free_vec)
  ad_hess = ForwardDiff.hessian(psf_fit_for_optim_val, psf_params_free_vec)

  @test_approx_eq expected_value sf_free.v[]
  @test_approx_eq sf_free.d[:] ad_grad
  @test_approx_eq sf_free.h ad_hess
end

function test_psf_optimizer()
  raw_psf = load_raw_psf()

  K = 2
  psf_params = initialize_psf_params(K, for_test=false)
  psf_transform = get_psf_transform(psf_params)
  psf_optimizer = PsfOptimizer(psf_transform, K)

  nm_result = fit_psf(psf_optimizer, raw_psf, psf_params)
  psf_params_fit =
    constrain_psf_params(unwrap_psf_params(Optim.minimizer(nm_result)), psf_transform)

  # Could this test be tighter?
  @test 0.0 < Optim.minimum(nm_result) < 1e-3

  celeste_psf = fit_raw_psf_for_celeste(raw_psf, K)[1]
  rendered_psf = get_psf_at_point(celeste_psf)

  @test_approx_eq Optim.minimum(nm_result) sum((raw_psf - rendered_psf) .^ 2)

  # Make sure that re-using the optimizer gets the same results.
  raw_psf_10_10 = load_raw_psf(x=10., y=10.)
  celeste_psf_10_10_v1, psf_params_10_10_v1 =
    fit_raw_psf_for_celeste(raw_psf_10_10, K)
  celeste_psf_10_10_v2, psf_params_10_10_v2 =
    fit_raw_psf_for_celeste(raw_psf_10_10, psf_optimizer, psf_params)
  for k=1:K
    @test_approx_eq psf_params_10_10_v1[k] psf_params_10_10_v2[k]
    for field in fieldnames(celeste_psf_10_10_v1[k])
      @test_approx_eq getfield(celeste_psf_10_10_v1[k], field) getfield(celeste_psf_10_10_v2[k], field)
    end
  end
end

function test_trim_psf()
    raw_psf = load_raw_psf()
    trim_percent = 0.95
    trimmed_psf = trim_psf(raw_psf; trim_percent=trim_percent)
    @test sum(abs(trimmed_psf)) >= trim_percent * sum(abs(raw_psf))
    @test size(trimmed_psf, 1) < size(raw_psf, 1)
    @test size(trimmed_psf, 2) < size(raw_psf, 2)
end

test_transform_psf_sensitive_float()
test_transform_psf_params()
test_psf_fit()
test_psf_optimizer()
test_trim_psf()


text:
function (kr::SingleRelPerms)(s)
    water_rel_perm(s, kr), oil_rel_perm(s, kr)
end

function (kr::SingleRelPerms)(s, θ)
    water_rel_perm(s, kr), oil_rel_perm(s, kr)
end

function (kr::MixedRelPerms)(s, θ)
    water_rel_perm(s, θ, kr), oil_rel_perm(s, θ, kr)
end

function water_rel_perm(sw::T, swr::P, sor::P, krw0::P, nw::P) where {T,P <: Float64}

    S  = (sw - swr) / (1 - sor - swr)

    if swr <= sw <= 1 - sor
        krw = krw0 * S^nw
    elseif 1 - sor < sw <= 1.0
        krw = krw0
    elseif sw < swr
        krw = 0.0
    end

end

function water_rel_perm(sw::AbstractVector{T}, swr::P, sor::P, krw0::P, nw::P) where {T,P <: Float64}
    water_rel_perm.(sw, swr, sor, krw0, nw)
end

function oil_rel_perm(sw::T, swr::P, sor::P, kro0::P, no::P) where {T,P <: Float64}

    S  = (sw - swr) / (1 - sor - swr)

    if swr <= sw <= 1 - sor
        kro = kro0 * (1 - S)^no
    elseif 1 - sor < sw <= 1.0
        kro = 0.0
    elseif sw < swr
        kro = kro0
    end
end

function oil_rel_perm(sw::AbstractVector{T}, swr::P, sor::P, kro0::P, no::P) where {T,P <: Float64}
    oil_rel_perm.(sw, swr, sor, kro0, no)
end

function oil_rel_perm(sw::Union{T, AbstractArray{T}}, kr::SingleRelPerms) where {T <: Float64}
    oil_rel_perm.(sw, kr.swr, kr.sor, kr.kro0, kr.no)
end

function oil_rel_perm(sw::Union{T, AbstractArray{T}}, θ::T, kr::MixedRelPerms) where {T <: Float64}
    swr = θ  *  kr.ww.swr + (1-θ)  *  kr.ow.swr
    sor = θ  *  kr.ww.sor + (1-θ)  *  kr.ow.sor
    kro0 = θ  *  kr.ww.kro0 + (1-θ)  *  kr.ow.kro0
    no = θ  *  kr.ww.no + (1-θ)  *  kr.ow.no

    oil_rel_perm.(sw, swr, sor, kro0, no)
end

function oil_rel_perm(sw::Union{T, AbstractArray{T}}, θ:: AbstractArray{T}, kr::MixedRelPerms) where {T <: Float64}
    kro(s,θ) = oil_rel_perm(s, θ, kr)
    kro.(sw, θ)
end

function water_rel_perm(sw::Union{T, AbstractArray{T}}, kr::SingleRelPerms) where {T <: Float64}
    water_rel_perm.(sw, kr.swr, kr.sor, kr.krw0, kr.nw)
end

function water_rel_perm(sw::Union{T, AbstractArray{T}}, θ::T, kr::MixedRelPerms) where {T <: Float64}
    swr = θ  *  kr.ww.swr + (1-θ)  *  kr.ow.swr
    sor = θ  *  kr.ww.sor + (1-θ)  *  kr.ow.sor
    krw0 = θ  *  kr.ww.krw0 + (1-θ)  *  kr.ow.krw0
    nw = θ  *  kr.ww.nw + (1-θ)  *  kr.ow.nw

    water_rel_perm.(sw, swr, sor, krw0, nw)
end

function water_rel_perm(sw::Union{T, AbstractArray{T}}, θ:: AbstractArray{T}, kr::MixedRelPerms) where {T <: Float64}
    krw(s,θ) = water_rel_perm(s, θ, kr)
    krw.(sw, θ)
end

function kro_derivative(sw::T, swr::P, sor::P, kro0::P, no::P) where {T,P <: Float64}
    derivative.(sw ->  oil_rel_perm(sw, swr, sor, kro0, no), sw)
end

function krw_derivative(sw::AbstractVector{T}, swr::P, sor::P, krw0::P, nw::P) where {T,P <: Float64}
    derivative.(sw ->  water_rel_perm(sw, swr, sor, krw0, nw), sw)
end

function fractional_flow(sw::T, swr::P, sor::P, krw0::P, kro0::P, nw::P, no::P, μw::P, μo::P) where {T,P <: Float64}
    krw = water_rel_perm(sw, swr, sor, krw0, nw)
    kro = oil_rel_perm(sw, swr, sor, kro0, no)
    λw = krw / μw
    λo = kro / μo

    λw ./ (λw .+ λo)
end

function fractional_flow(sw::AbstractVector{T}, swr::P, sor::P, krw0::P, kro0::P, nw::P, no::P, μw::P, μo::P) where {T,P <: Float64}
    fractional_flow.(sw, swr, sor, krw0, kro0, nw, no, μw, μo)
end

function fractional_flow(sw::Union{T,AbstractVector{T}}, kr::RelPerms, μw::P, μo::P) where {T,P <: Float64}
    swr = kr.swr
    sor = kr.sor
    krw0 = kr.krw0
    nw = kr.nw
    kro0 = kr.kro0
    no = kr.no

    krw = water_rel_perm(sw, swr, sor, krw0, nw)
    kro = oil_rel_perm(sw, swr, sor, kro0, no)
    λw = krw / μw
    λo = kro / μo

    λw ./ (λw .+ λo)
end

function fw_derivative(sw::Union{T,AbstractVector{T}}, swr::P, sor::P, krw0::P, kro0::P, nw::P, no::P, μw::P, μo::P) where {T,P <: Float64}
    derivative.(sw ->  fractional_flow(sw, swr, sor, krw0, kro0, nw, no, μw, μo), sw)
end

function fw_derivative(sw::Union{T,AbstractVector{T}}, kr::RelPerms, μw::P, μo::P) where {T,P <: Float64}
    swr = kr.swr
    sor = kr.sor
    krw0 = kr.krw0
    nw = kr.nw
    kro0 = kr.kro0
    no = kr.no

    derivative.(sw ->  fractional_flow(sw, swr, sor, krw0, kro0, nw, no, μw, μo), sw)
end