text:
using KUparser

# CONLL-U FORMAT:
# 1. ID: Word index, integer starting at 1 for each new sentence; may be a range for multiword tokens; may be a decimal number for empty nodes.
# 2. FORM: Word form or punctuation symbol.
# 3. LEMMA: Lemma or stem of word form.
# 4. UPOSTAG: Universal part-of-speech tag.
# 5. XPOSTAG: Language-specific part-of-speech tag; underscore if not available.
# 6. FEATS: List of morphological features from the universal feature inventory or from a defined language-specific extension; underscore if not available.
# 7. HEAD: Head of the current word, which is either a value of ID or zero (0).
# 8. DEPREL: Universal dependency relation to the HEAD (root iff HEAD = 0) or a defined language-specific subtype of one.
# 9. DEPS: Enhanced dependency graph in the form of a list of head-deprel pairs.
# 10. MISC: Any other annotation.

# Only changes the HEAD and DEPREL columns of the inputfile
function writeconllu(sentences, parses, inputfile)
    if length(sentences) != length(parses); error(); end
    v = sentences[1].vocab
    s = p = nothing
    ns = nw = nl = 0
    for line in eachline(inputfile)
        nl += 1
        if ismatch(r"^\d+\t", line)
            # info("$nl word")
            if s == nothing
                s = sentences[ns+1]
                p = parses[ns+1]
            end
            f = split(line, '\t')
            nw += 1
            if f[1] != "$nw"; error(); end
            if f[2] != v.words[s.word[nw]]; error(); end
            f[7] = string(p.head[nw])
            f[8] = v.deprels[p.deprel[nw]]
            print(join(f, "\t"))
        else
            if line == "\n"
                # info("$nl blank")
                if s == nothing; error(); end
                if nw != length(s.word); error(); end
                ns += 1; nw = 0
                s = p = nothing
            else
                # info("$nl non-word")
            end
            print(line)
        end
    end
    if ns != length(sentences); error(); end
end

function readconllu(filename::AbstractString, vocab=nothing)
    if vocab == nothing
        vocab = Vocab()
        vocab.postags = UPOSTAG
        vocab.deprels = UDEPREL
    end
    words = Dict{String,WordId}()
    for id in 1:length(vocab.words)
        words[vocab.words[id]] = id
    end
    # info(filename)
    sentences = Sentence[]
    s = Sentence(); s.vocab = vocab
    nw = 0
    for line in eachline(filename)
        line = chomp(line)
        fields = split(line, '\t') # TODO: this makes fields[2] a substring preventing garbage collection.
        if line == ""
            if nw > 0
                nw = 0
                push!(sentences,s)
                s = Sentence(); s.vocab = vocab
            end
        elseif line[1] == '#'
            continue
        elseif ismatch(r"^\d+$", fields[1])
            nw += 1
            if fields[1] != "$nw"; error("Token out of order in [$line]"); end

            if !haskey(words, fields[2])
                id = 1 + length(words)
                wform = String(fields[2])
                words[wform] = id
                push!(vocab.words, wform)
            end
            push!(s.word, words[fields[2]])

            if fields[4]=="_"
                pt = 0
            else
                pt = findfirst(vocab.postags, fields[4])
                if pt == 0; Base.warn_once("Bad postag: $(fields[4])"); end
            end
            push!(s.postag, pt)

            if fields[7]=="_"
                hd = 0
            else
                hd = parse(Position,fields[7])
                if !(0 <= hd <= typemax(eltype(s.head))); Base.warn_once("Bad head: $hd"); hd=0; end
            end
            push!(s.head, hd)

            if fields[8]=="_"
                dr = 0
            else
                dr = fields[8]
                # Get rid of language specific extensions
                drcolon = search(dr, ':')
                if drcolon > 0; dr = dr[1:drcolon-1]; end
                dr = findfirst(vocab.deprels, dr)
                if dr == 0; Base.warn_once("Bad deprel: $(fields[8])"); end
            end
            push!(s.deprel, dr)

        elseif ismatch(r"^\d+[-.]\d+$", fields[1])
            # skip
        else
            error("Cannot parse [$line]")
        end # if line == ""
    end # for line in eachline(filename)
    return sentences
end

function readvectors(filename::AbstractString, vocab::Vocab)
    lcwords = Dict{String,Vector{WordId}}()
    for id in 1:length(vocab.words)
        lcword = lowercase(vocab.words[id])
        if !haskey(lcwords, lcword); lcwords[lcword] = WordId[]; end
        push!(lcwords[lcword], id)
    end
    # info(filename)
    open(`xzcat $filename`) do io
        firstline = split(readline(io)) # first line gives count and dim
        wdims = parse(Int,firstline[2])
        # TODO: We are leaving unknown word vectors as zero!
        vocab.wvecs = zeros(WVtype, wdims, length(vocab.words))
        for line in eachline(io)
            n = search(line, ' ')
            w = SubString(line,1,n-1)
            if !haskey(lcwords,w); continue; end
            a = map(s->parse(Float32,s), split(SubString(line,n+1)))
            if wdims != length(a); throw(DimensionMismatch()); end
            for i in lcwords[w]
                vocab.wvecs[:,i] = a
            end
        end
    end
end

# Universal POS tags (17)
const UPOSTAG = String[ # Dict{String,UInt8}(
"ADJ"   ,# => 1, # adjective
"ADP"   ,# => 2, # adposition
"ADV"   ,# => 3, # adverb
"AUX"   ,# => 4, # auxiliary
"CCONJ" ,# => 5, # coordinating conjunction
"DET"   ,# => 6, # determiner
"INTJ"  ,# => 7, # interjection
"NOUN"  ,# => 8, # noun
"NUM"   ,# => 9, # numeral
"PART"  ,# => 10, # particle
"PRON"  ,# => 11, # pronoun
"PROPN" ,# => 12, # proper noun
"PUNCT" ,# => 13, # punctuation
"SCONJ" ,# => 14, # subordinating conjunction
"SYM"   ,# => 15, # symbol
"VERB"  ,# => 16, # verb
"X"     ,# => 17, # other
]

# Universal Dependency Relations (37)
const UDEPREL = String[ # Dict{String,UInt8}(
"root"       ,# => 1, # root
"acl"        ,# => 2, # clausal modifier of noun (adjectival clause)
"advcl"      ,# => 3, # adverbial clause modifier
"advmod"     ,# => 4, # adverbial modifier
"amod"       ,# => 5, # adjectival modifier
"appos"      ,# => 6, # appositional modifier
"aux"        ,# => 7, # auxiliary
"case"       ,# => 8, # case marking
"cc"         ,# => 9, # coordinating conjunction
"ccomp"      ,# => 10,# clausal complement
"clf"        ,# => 11, # classifier
"compound"   ,# => 12, # compound
"conj"       ,# => 13, # conjunct
"cop"        ,# => 14, # copula
"csubj"      ,# => 15, # clausal subject
"dep"        ,# => 16, # unspecified dependency
"det"        ,# => 17, # determiner
"discourse"  ,# => 18, # discourse element
"dislocated" ,# => 19, # dislocated elements
"expl"       ,# => 20, # expletive
"fixed"      ,# => 21, # fixed multiword expression
"flat"       ,# => 22, # flat multiword expression
"goeswith"   ,# => 23, # goes with
"iobj"       ,# => 24, # indirect object
"list"       ,# => 25, # list
"mark"       ,# => 26, # marker
"nmod"       ,# => 27, # nominal modifier
"nsubj"      ,# => 28, # nominal subject
"nummod"     ,# => 29, # numeric modifier
"obj"        ,# => 30, # object
"obl"        ,# => 31, # oblique nominal
"orphan"     ,# => 32, # orphan
"parataxis"  ,# => 33, # parataxis
"punct"      ,# => 34, # punctuation
"reparandum" ,# => 35, # overridden disfluency
"vocative"   ,# => 36, # vocative
"xcomp"      ,# => 37, # open clausal complement
]

# Universal Features (21) with number of possible values
# Abbr: abbreviation (1) Yes (all features also have None as an extra choice)
# Animacy: animacy (4) Anim,Hum,Inan,Nhum
# Aspect: aspect (6) Hab,Imp,Iter,Perf,Prog,Prosp
# Case: case (32)
# Definite: definiteness or state (5) Com,Cons,Def,Ind,Spec
# Degree: degree of comparison (5) Abs,Cmp,Equ,Pos,Sup
# Evident: evidentiality (2) Fh,Nfh
# Foreign: is this a foreign word? (1) Yes
# Gender: gender (4) Com,Fem,Masc,Neut
# Mood: mood (12) Adm,Cnd,Des,Imp,Ind,Jus,Nec,Opt,Pot,Prp,Qot,Sub
# NumType: numeral type (7) Card,Dist,Frac,Mult,Ord,Range,Sets
# Number: number (10) CollCountDualGrpaGrplInvPaucPlurPtanSingTri
# Person: person (5) 0,1,2,3,4
# Polarity: polarity (2) Neg,Pos
# Polite: politeness (4) ElevFormHumbInfm
# Poss: possessive (1) Yes
# PronType: pronominal type (11) ArtDemEmpExcIndIntNegPrsRcpRelTot
# Reflex: reflexive (1) Yes
# Tense: tense (5) FutImpPastPqpPres
# VerbForm: form of verb or deverbative (8) ConvFinGdvGerInfPartSupVnoun
# Voice: voice (8) ActAntipCauDirInvMidPassRcp

# conll17 word vectors
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Ancient_Greek/grc.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Arabic/ar.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Basque/eu.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Bulgarian/bg.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Catalan/ca.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/ChineseT/zh.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Croatian/hr.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Czech/cs.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Danish/da.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Dutch/nl.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/English/en.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Estonian/et.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Finnish/fi.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/French/fr.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Galician/gl.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/German/de.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Greek/el.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Hebrew/he.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Hindi/hi.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Hungarian/hu.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Indonesian/id.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Irish/ga.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Italian/it.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Japanese/ja.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Kazakh/kk.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Korean/ko.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Latin/la.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Latvian/lv.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Norwegian-Bokmaal/no_bokmaal.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Norwegian-Nynorsk/no_nynorsk.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Old_Church_Slavonic/cu.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Persian/fa.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Polish/pl.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Portuguese/pt.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Romanian/ro.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Russian/ru.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Slovak/sk.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Slovenian/sl.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Spanish/es.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Swedish/sv.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Turkish/tr.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Ukrainian/uk.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Urdu/ur.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Uyghur/ug.vectors.xz
# /mnt/ai/data/nlp/conll17/word-embeddings-conll17/Vietnamese/vi.vectors.xz

nothing


text:
# -*- coding: utf-8 -*-
using ModelingToolkit
using MacroTools
import MacroTools: postwalk
using Petri
using SemanticModels.ModelTools.PetriModels
using Test

# SIR  <- IR  -> SEIR
#  |       |      |
#  v       v      v
# SIRS <- IRS -> SEIRS

# +
using Catlab.WiringDiagrams
using Catlab.Doctrines
import Catlab.Doctrines.⊗
import Catlab.Graphics: to_graphviz
import Catlab.Graphics.Graphviz: run_graphviz
⊗(a::WiringDiagram, b::WiringDiagram) = otimes(a,b)
import Base: ∘
∘(a::WiringDiagram, b::WiringDiagram) = compose(b, a)
⊚(a,b) = b ∘ a
S, E, I, R, D= Ob(FreeSymmetricMonoidalCategory, :S, :E, :I, :R, :D)

inf  = WiringDiagram(Hom(:infection, S ⊗ I, I⊗I))
expo = WiringDiagram(Hom(:exposure, S⊗I, E⊗I))
rec  = WiringDiagram(Hom(:recovery, I,   R))
wan  = WiringDiagram(Hom(:waning,   R,   S))

sir_wire  = inf ⊚ (rec ⊗ rec)

sir = model(PetriModel, sir_wire)

dump(sir)

# +
@variables S, E, I, R

dump(S)
# -

sir = Petri.Model([S, I, R], [(I, R), (S+I, 2I)])

ir = Petri.Model([S, I, R], [(I, R)])

seir = Petri.Model([S, I, R], [(I, R), (S+I, I+E), (E, I)])

# +
irs = Petri.Model([S, I, R], [(I, R), (R, S)])

dump(sir)
dump(irs)
# -

sirs = Petri.Model([S, I, R], [(I, R), (S+I, 2*I), (R, S)])

rule = PetriModels.Span(sir, ir, seir)
sirs′ = PetriModels.pushout(irs, sir)
@test sirs′.Δ == sirs.Δ
seirs = PetriModels.solve(PetriModels.DPOProblem(rule, irs))
@test all(Set(seirs.Δ) .== Set([(S+I, I+E),
                             (E, I),
                             (I, R),
                             (R, S)]))
l = sir
c = ir
r = seir
c′ = irs

l′ = PetriModels.pushout(l, c′)
@test l′.Δ == sirs.Δ
@test PetriModels.dropdown(l,c,l′).Δ == c′.Δ
@test PetriModels.pushout(r, c′).Δ == seirs.Δ

function Δ(m::Petri.Model, ctx=:state)
    function updateblock(exp, sym)
        return postwalk(exp) do x
            if typeof(x) == Expr && x.head == :call
                if length(x.args) == 1
                    var = x.args[1]
                    # push!(args, var)
                    e = Expr(sym, :($ctx.$var), 1)
                    # @show "adding guard"
                    if sym == :-=
                        return quote
                            ($ctx.$var > 0 || return nothing ) && $e
                        end
                    end
                    return e
                end
                if length(x.args) >= 1 && x.args[1] == :(*)
                    op = x
                    try
                        # @info "trying"
                        # @show x
                        branch = x.args[3].args[2]
                        # @show branch
                        # @show branch.head
                        # @show branch.args[1]
                        if branch.head == :&&
                            # @info "&& found"
                            op = branch.args[2]
                            # @show op
                            op.args[end] = x.args[2]
                            # @show x
                            return x.args[3]
                        end
                    catch
                        # @info "catching: there was no branch"
                        changevalue = x.args[2]
                        statename = x.args[3].args[1]
                        # e = Expr(sym, statename, changevalue)
                        # return e
                        x.args[3].args[2] = changevalue
                        return x.args[3]
                    end
                end
                if length(x.args) >= 1 && x.args[1] == :(+)
                    # @show x
                    return quote
                        $(x.args[2:end]...)
                    end
                end
            end
        return x
        end
    end

    head(x) = try
        x.head
    catch
        nothing
    end

    function poolconditions(decrements)
        if decrements.head == :block
            steps = postwalk(MacroTools.striplines, decrements).args
            checks = Expr[]
            events = Expr[]
            map(steps) do s
                postwalk(s) do x
                    if head(x) == :&&
                        push!(checks, x.args[1])
                        push!(events, x.args[2])
                    end
                    return x
                end
            end
            decrements = quote $(checks...); $(events...)  end
        end
        return decrements
    end
    δf = map(m.Δ) do δ
        q = quote end
        # input states get decremented
        parents = δ[1]
        children = δ[2]

        exp1 = convert(Expr, parents)
        decrements = updateblock(exp1, :-=) |> poolconditions

        exp2 = convert(Expr, children)
        increments = updateblock(exp2, :+=)

        push!(q.args, decrements)
        push!(q.args, increments)

        sym = gensym("δ")
        @show MacroTools.striplines(q)
        :($sym(state) = $(q) )
    end
end

@show "SIR"
Δ(l, :state)
# @show "IR"
# funckit(c, :state)
@show "SEIR"
# funckit(r, :state)

@show "SIRS"
# funckit(l′, :state)
# @show "IRS"
# funckit(c′, :state)
@show "SEIRS"
# funckit(seirs, :state)

exprs = Δ(sirs, :state)
m = Petri.Model([S, I, R], exprs, [
    quote
    λ_2(state) = state.γ * state.I
    end,
    quote
    λ_1(state) = state.β * state.S * state.I / +(state.S, state.I, state.R)
    end,
    quote
    λ_3(state) = state.μ * state.R
    end],
                [
                    quote b_2(state) = state.I > 0 end,
                    quote b_1(state) = state.S > 0 && state.I > 0 end,
                    quote b_3(state) = state.R > 0 end]
                )

# +
mutable struct SIRState{T,F}
    S::T
    I::T
    R::T
    β::F
    γ::F
    μ::F
end

mutable struct ParamSIR{T, P}
    S::T
    I::T
    R::T
    params::P
end
# -

p = Petri.Problem(eval(m), SIRState(100, 1, 0, 0.5, 0.15, 0.05), 150)
#@show Petri.solve(p)

# @show Petri.funckit(Petri.Problem(l, missing, 10), :state)
function test_1()
    no_transitions = Tuple{Operation, Operation}[]
    @variables A, B, C, D
    states = [A, B, C, D]
    l = Petri.Model(states, [(A, B)])
    c = Petri.Model(states, no_transitions)
    r = Petri.Model(states, [(A, B + C)])
    rule = PetriModels.Span(l, c, r)
    c′ = Petri.Model(states, [(B, A)])
    r′ = PetriModels.solve(PetriModels.DPOProblem(rule, c′))
    @test r′.Δ == [(A, B+C), (B, A)]

    l′ = PetriModels.pushout(l, c′)
    @test l′.Δ == [(A, B), (B, A)]
    @test PetriModels.dropdown(l,c,l′).Δ == [(B, A)]
end

function test_2()
    no_transitions = Tuple{Operation, Operation}[]
    @variables A, B, C, D
    states = [A, B, C, D]
    l = Petri.Model(states, [(A, B), (B,C)])
    c = Petri.Model(states, no_transitions)
    r = Petri.Model(states, [(A, C)])
    rule = PetriModels.Span(l, c, r)
    c′ = Petri.Model(states, [(C, D)])
    r′ = PetriModels.solve(PetriModels.DPOProblem(rule, c′))
    @test r′.Δ == [(A, C), (C, D)]

    l′ = PetriModels.pushout(l, c′)
    @test l′.Δ == [(A, B), (B, C), (C, D)]
    @test PetriModels.dropdown(l,c,l′).Δ == [(C, D)]
    @test PetriModels.pushout(rule.r, c′).Δ == [(A, C), (C, D)]
end
test_1()
test_2()

exprs2 = Δ(Petri.Model([S,I,R], [(2S+I, 3I)]))

function Λ(m::Petri.Model{G, S, D, L, B}) where {G, S, D, L, B}
    head(x) = try
        x.head
    catch
        nothing
    end
    function ratecomp(exp, ctx)

        args = Dict{Symbol, Int}()
        postwalk(convert(Expr, exp)) do x
            if typeof(x) == Expr && x.head == :call
                if length(x.args) == 1
                    var = x.args[1]
                    args[var] = 1
                    return :($ctx.$var)
                end
                if x.args[1] == :(*)
                    args[var] = x.args[2]
                end
            end
        return x
        end
        return args
    end

    δf = map(enumerate(m.Δ)) do (i, δ)
        # input states are used to calc the rates
        parents = δ[1]

        # exp1 = convert(Expr, parents)
        ctx = :state
        rates = ratecomp(parents, ctx)
        q = :(*())
        map(collect(keys(rates))) do s
            r = rates[s]
            push!(q.args, :($ctx.$s / $r))
        end
        term = :($ctx.params[$i])
        push!(q.args, term)

        sym = gensym("λ")
        @show MacroTools.striplines(q)
        :($sym(state) = $(q) )
    end
end

function funckit(m::Petri.Model)
    return Petri.Model(m.g, m.S, Δ(m), Λ(m), missing)
end

function Petri.eval(m::Petri.Model{G, Z, D, L, Missing}) where {G, Z, D, L}
    Petri.Model(m.g, m.S, eval.(m.Δ), eval.(m.Λ), missing)
end

sir′ = funckit(sir)
m = sir′
m′ = eval(m)
@show m′
@show typeof(m′)
p = Petri.Problem(m′, ParamSIR(100, 1, 0, [0.15, 0.55/101]), 250)
soln = Petri.solve(p)
@test soln.S <= 10
@test soln.S + soln.I + soln.R == 101

m′ = eval(funckit(sirs))
@show m′
@show typeof(m′)
p = Petri.Problem(m′, ParamSIR(100, 1, 0, [0.15, 0.55/101, 0.15]), 250)
sirs_soln = Petri.solve(p)
@test sirs_soln.S >= 10
@test sirs_soln.S + sirs_soln.I + sirs_soln.R == 101
@show soln
@show sirs_soln

m′ = Petri.eval(funckit(seirs))
# @code_native m′.Δ[1](ParamSIR(100, 1, 0, [ 0.15, 0.55/101, 0.15, 0.1 ]))

ext: jl


text:
module CloudScat
using StaticArrays
using LinearAlgebra
using Interpolations
using Dates
using HDF5
using Parameters
using Formatting
using Logging
using Dates
using ProgressMeter
using Random
using Interpolations
using FastGaussQuadrature
using DelimitedFiles

export Params, Point, World, Observer, Homogeneous, VariableNR
export NoAbsorption, StratifiedAbsorption
export time, photo, image, stime

const Point = SVector{3, Float64}

include("constants.jl")
include("parameters.jl")
include("mieparams.jl")

include("composition.jl")
include("absorption.jl")
include("geometry.jl")
#include("phasefuncs.jl")
include("rayleigh.jl")

using .MieParams

const co = constants
const SAMPLES_PATH = normpath(joinpath(@__DIR__, "..", "samples"))
const UTIL_PATH = normpath(joinpath(@__DIR__, "..", "util"))

# For the dynamical dispatch we define types of collisions
abstract type ScatteringType end
struct Mie <: ScatteringType
    g::Float64
    ω0::Float64
end
struct Rayleigh <: ScatteringType end
struct Null <: ScatteringType end
struct Isotropic <: ScatteringType end

const NonNull = Union{Mie, Rayleigh, Isotropic}

# Phase functions
p(s::Mie, μ) = (1 / 4π) * (1 - s.g^2) / (1 + s.g^2 - 2 * s.g * μ)^(3/2)
p(s::Rayleigh, μ) = (3 / 16π) * (1 + μ^2)
p(s::Isotropic, μ) = (1 / 4π)

# Sampling from the phase functions.  Here ξ is a random number uniformly
# distributed in [0, 1).
sample(s::Mie, ξ) = (1 + s.g^2 - ((1 - s.g^2) / (1 + s.g * (2 * ξ - 1)))^2) / (2 * s.g)
function sample(s::Rayleigh, ξ)
    z = 2 * (2ξ - 1)
    w = z + sqrt(z^2 + 1)
    w^(1//3) - w^(-1//3)
end

# Test absorption of the particle
function reweight(s::Mie, w, wmin)
    w1 = s.ω0 * w
    if w1 < wmin
        w1 = trand() < w1 ? 1.0 : 0.0
    end
    w1
end
reweight(s::Rayleigh, w, wmin) = w
reweight(s::Isotropic, w, wmin) = w
reweight(s::Null, w, wmin) = w

# Whether a particle can scatter
scatters(s::Null) = false
scatters(s::Union{Rayleigh, Mie, Isotropic}) = true

struct World{Tc,Td,Tcomp,Tabs <: AbstractAbsorber,N}
    cloud::Tc
    domain::Td
    comp::Tcomp
    absorber::Tabs
    quadrule::Tuple{SVector{N,Float64}, SVector{N,Float64}}

    function World(cloud::Tc, domain::Td, comp::Tcomp, absorb::Tabs,
                   n::Int) where {Tc, Td, Tcomp, Tabs}
        x, w = gausslegendre(n)
        quadrule = (SVector{n}(x), SVector{n}(w))
        new{Tc, Td, Tcomp, Tabs, n}(cloud, domain, comp, absorb, quadrule)
    end

    World(cloud, domain, comp) = World(cloud, domain, comp, NoAbsorption(), 3)
    World(cloud, domain, comp, absorb) = World(cloud, domain, comp, absorb, 3)
end

"""
    Observer

Description of an observer.  An observer is located a some position `r` and
collects photons, building

 * a timeline of arriving photons (`obs`) with a time resolution `δt`
 * an image of the arriving directions of the photons (`img`) with cosine-angle
   resolution `δμ` and field-of-view `μmax`.

"""
struct Observer
    r::SVector{3, Float64}   # Location of the observer
    δt::Float64              # Time-resolution of the detector
    δu::Float64              # Spatial resoution of the camera (in cosine space)
    umax::Float64            # Field-of-view in cosine space (i.e. max. observable mu)
    obs::Array{Float64, 2}   # Observations
    img::Array{Float64, 3}   # Image array
end

"""
    Observer(obsdata)

Create a observer structure from a dictionary of observer properties.
"""
function Observer(;position, fov, pixels, tsample, nsamples)
    sposition = SVector{3, Float64}(position)
    umax = tan(deg2rad(fov)) / sqrt(2)
    δu = 2 * umax / pixels
    nt = Threads.nthreads()

    obs = Observer(sposition,
                   tsample,
                   δu,
                   umax,
                   zeros(nsamples, nt),
                   zeros(pixels, pixels, nt))
end

"""
    projectpx(obs, px, H)

  Project pixel coordinates into real space.

  Finds a pair `[x, y]` of positions that, at an altitude `H` would be projected
  into pixels `px = [xp, yp]` of the camera.  This is useful to simulate
  observations with the same geometry as an existing observation.
"""
function projectpx(obs::Observer, px::AbstractVector, H)
    @. -(obs.r[3] - H) * (obs.umax - px * obs.δu) + obs.r[[1, 2]]
end

@with_kw struct Results
    time::Vector{Float64}
    photo::Vector{Float64}
    image::Array{Float64, 2}
    delay::Float64
end

time(r::Results) = r.time
stime(r::Results) = time(r) .- r.delay
photo(r::Results) = r.photo
image(r::Results) = r.image

function main(params::Params, world::World, observers::Vector{Observer};
              saveto::Union{String,Nothing}=nothing, savecallbacks=(),
              verbose=true)

    function fmt(level, _module, group, id, file, line)
        return (:blue, format("{:<23}:", string(Dates.now())), "")
    end
    logger = ConsoleLogger(meta_formatter=fmt)

    with_logger(logger) do
        if verbose
            @info "CloudScat (c) Alejandro Luque IAA-CSIC, 2020"
            @info "Mie solver by Olli Wilkman (https://github.com/dronir/MieScatter.jl)"
            # Print a list of parameters
            @info "Input parameters:" params

            @info "Running with $(Threads.nthreads()) thread(s) (JULIA_NUM_THREADS)"

            # Run the MC simulation
            @info "Running the MC simulation"
        end
        run!(world, observers, params; verbose=verbose)

        # Build the results object to return; one for each observer
        results = map(observers) do obs
            rsource = 0.5 * (params.source_a + params.source_b)
            Results(time=[obs.δt * (i - 0.5) for i in 1:size(obs.obs, 1)],
                    # The dropdims(sum(...)) is to sum over each thread
                    photo=dropdims(sum(obs.obs, dims=2), dims=2),
                    image=dropdims(sum(obs.img, dims=3), dims=3),
                    delay=norm(obs.r - rsource) / co.c)
        end

        !isnothing(saveto) && save(saveto, observers, world, params,
                                   savecallbacks)

        results
    end
end

"""
    save(fname, p, observers, params)

Save the population and observations into a h5 file called `fname`.
"""
function save(fname, observers::Vector{Observer}, world::World, params::Params,
              callbacks=())

    h5open(fname, "w") do file
        g = create_group(file, "parameters")
        for field in fieldnames(Params)
            attributes(g)[String(field)] = getfield(params, field)
        end

        args = ()

        for (i, obs) in enumerate(observers)
            g = create_group(file, format("obs{:05d}", i))
            attributes(g)["altitude"] = obs.r[3]
            attributes(g)["shift"] = obs.r[1]

            rsource = 0.5 * (params.source_a + params.source_b)
            delay = norm(obs.r - rsource) / co.c
            @show delay
            attributes(g)["delay"] = delay

            g["t", shuffle = (), deflate = 3] =
                [obs.δt * (i - 0.5) for i in 1:size(obs.obs, 1)]

            # The dropdims(sum(...)) is to sum over each thread
            g["timeline", shuffle = (), deflate = 3] =
                dropdims(sum(obs.obs, dims=2), dims=2)

            g["image", shuffle = (), deflate = 3] =
                dropdims(sum(obs.img, dims=3), dims=3)

            g["depth", shuffle = (), deflate = 3] =
                depthimg(obs, world, params)

        end

        for cb in callbacks
            cb(file, observers, world, params)
        end

    end

    @info "Output written in $fname"
end

""" Compute and save droplet density in a grid. `cloud` is a geometrical
    shape defining the cloud extent, `comp` is a composition.
"""
function cloudgrid(cloud, comp, axes)
    ax, ay, az = axes

    n = zeros((length(ax), length(ay), length(az)))
    for (k, z) in enumerate(az), (j, y) in enumerate(ay), (i, x) in enumerate(ax)
        r = @SVector [x, y, z]
        inside(cloud, r) || continue

        loc = probe(comp, r)
        n[i, j, k] = density(comp, r, loc)
    end

    n
end

"""
Construct a saving callback to save a the save droplet density in a grid
"""
function savecloudgrid(axes...)
    (file, observers, world, params) -> begin
        g = create_group(file, "cloudgrid")

        n = cloudgrid(world.cloud, world.comp, axes)

        write(g, "density", n)
        x, y, z = axes
        write(g, "x", collect(x))
        write(g, "y", collect(y))
        write(g, "z", collect(z))
    end
end

"""
    newphoton(params)

Initialize a a new photon..
"""
function newphoton(params::Params)
    @unpack source_a, source_b = params

    ξ = trand()
    r = (1 - ξ) * source_a + ξ * source_b
    μ = randsphere()
    t = 0.0
    w = 1.0

    r, μ, t, w
end

"""
    run!(p, world, observers, interps, params)

Run the MC simulation on a photon population and a collection of observers.
"""
function run!(world::World, observers::Vector{Observer}, params::Params;
              verbose=true)
    @unpack N = params

    # This is to allow quick debug runs without changing an input file.
    N ÷= parse(Int, get(ENV, "CLOUDSCAT_DEBUG_REDUCTION", "1"))

    set_zero_subnormals(true)

    if verbose
        prog = Progress(N, 5)
        prog_lock = Threads.SpinLock()
    end

    count = Threads.Atomic{Int}(0)

    Threads.@threads for tid in 1:Threads.nthreads()
        @inbounds for i in getrange(tid, N)
            r, μ, t, w = newphoton(params)

            for o in observers
                observeone!(o, world, Isotropic(), r, μ, t, w, params)
            end

            iterate(r, μ, t, w, world, observers, params)

            if verbose
                lock(prog_lock)
                ProgressMeter.next!(prog)
                unlock(prog_lock)
            end
        end
    end
end

"""
    iterate(r, μ, t, w, world, observers, params)

Iterate over steps for a single particle with starting position `r`,
direction `μ`, time `t` and weight `w`.  The particles lives in `world` and
may be detected by a list of `observers`.  `params` are the simulation parameters.
"""
function iterate(r, μ, t, w, world::World, observers::Vector{Observer},
                  params::Params)
    @unpack max_iter = params

    for it in 1:max_iter
        (w > 0) || break

        # Propagate
        l = travel_null(r, μ, world, params)

        r1 = r + μ * l
        t += l / co.c

        # Check if inside the domain
        if !inside(world.domain, r1)
            w = 0.0
            break
        end

        # Absorption during the propagation
        w *= transmittance(world.absorber, r, r1, μ)
        r = r1

        # Choose scattering type
        scat = choosescat(r, world, params)

        # Accumulate observations
        if r[3] > params.zmin
            for o in observers
                observeone!(o, world, scat, r, μ, t, w, params)
            end
        end

        # Scatter and check if absorbed
        μ, w = scatterone(μ, w, scat, params)
    end
end

"""
    observeone!(o, w, scat, r, μ, t, weight, params)

Check one particle and one observer and add the particle's contribution to
the timeline and image.

NOTE: Absorption is not considered here: it would simply add a factor ω₀
"""
@inline @fastmath function observeone!(o::Observer, w::World,
                                       s::NonNull,
                                       r, μ, t, weight, params::Params)
    @unpack N, νray_ground, H = params
    # Distance to the observer
    sobs = norm(o.r - r)

    # Director to the observer
    μobs = (o.r - r) / sobs

    μscat = μ ⋅ μobs

    # Optical depth to the observer
    ν(r::Point) = miecollrate(r, w, params)

    τmie = optpath(w.cloud, ν, r, o.r, w.quadrule)

    # Optical depth to the observer from Rayleigh scattering
    τray = νray_ground * H / μobs[3] * (exp(-r[3] / H) - exp(-o.r[3] / H))

    τ = τmie + τray

    # Transmittance due to background absorption
    trans = transmittance(w.absorber, r, o.r, μobs)

    # The probability dens that the photon reaches the cloud top w/o scattering
    # NB: p(μ) gives the scattering probablity per unit solid angle,
    # imagine that the photon is a super-particle representing N real photons. Now
    # - the number of photons dispersed in a solid angle dΩ is then N p(μ) dΩ,
    # - the surface at the observer subtended by this solid angle is s^2 dΩ.
    # So the number of photons per unit surface at the observer is
    # N p(μ) / s^2.  To this we have to add an attenuation factor exp(-s / Λ).
    f = p(s, μscat) * exp(-τ) * trans / (sobs^2 * N)

    # Now find time of observation if the particle manages to escape
    tobs = t + sobs / co.c

    tid = Threads.threadid()

    # Update the observation curve
    ind = 1 + Int64(fld(tobs, o.δt))
    if ind <= size(o.obs, 1)
        @inbounds o.obs[ind, tid] += weight * f / o.δt
    end

    # Now update the image at the given pixels
    # The minus are here because to be more user-friendly we invert the
    # image formed in the sensor.
    px = 1 + Int64(fld(-μobs[1] / μobs[3] + o.umax, o.δu))
    py = 1 + Int64(fld(-μobs[2] / μobs[3] + o.umax, o.δu))

    if 0 < px <= size(o.img)[1] && 0 < py <= size(o.img)[2]
        # The solid angle subtended by the pixel is roughly
        # δu^2 * (cos θ)^3, where θ is the angle with the vertical.
        # Actually cos θ = μobs[3].  The units would be here
        # photons / sr / source photon.
        @inbounds o.img[px, py, tid] += weight * f / (o.δu^2 * μobs[3]^3)
    end
end

observeone!(o::Observer, w, ::Null, r, μ, t, weight, params::Params) = nothing

""" Compute a range appropriate for this thread.
    Obtained from https://stackoverflow.com/questions/52593588/julia-why-doesnt-shared-memory-multi-threading-give-me-a-speedup
"""
function getrange(tid, n)
    nt = Threads.nthreads()
    d, r = divrem(n, nt)
    from = (tid - 1) * d + min(r, tid - 1) + 1
    to = from + d - 1 + (tid <= r ? 1 : 0)
    from:to
end

"""
    randsphere()

Sample points from the unitary sphere.
"""
function randsphere()
    ϕ = 2π * trand()
    sinϕ, cosϕ = sincos(ϕ)

    u = 2 * trand() - 1
    v = sqrt(1 - u^2)

    @SVector [v * cosϕ, v * sinϕ, u]
end

"""
    travel(r, μ, params)

Sample the travel distance a single particle given its position `r` and
its direction `μ`.  This function considers the interfaces in the cloud geometry
and limits the particle movement up to the next interface.  It also considers
different collision rates in and out of the cloud.

DEPRECATED: DOES NOT WORK
"""
@fastmath function travel(r, μ, world::World, params::Params; ϵ=1e-3)
    @unpack νraymax = params
    νmax = νraymax + νmiemax(world.comp)

    # Select the optical depth to travel
    τ = -log(trand())

    # To check whether we use collision rates in or outside the cloud
    # we test a point slightly ahead of the particle in its path.  This is
    # to allow cases where the particle is exactly at an interface.
    rp = r + (ϵ / νmax) * μ

    # Still to do: use local collision rates.
    ν = inside(world.cloud, rp) ? νmax : νraymax

    l = τ / ν
    rend = r + l * μ

    minimizer = MinAccumulator(1.0)
    interfaces!(minimizer, world.cloud, r, rend)
    l * minimizer.value
end

"""
    travel_null(r, μ, params)

Sample the travel distance a single particle given its position `r` and
its direction `μ`.
"""
@fastmath function travel_null(r, μ, world::World, params::Params; ϵ=1e-3)
    @unpack νraymax = params
    νmax = νraymax + νmiemax(world.comp)

    # Select the optical depth to travel
    τ = -log(trand())
    τ / νmax
end

"""
    choosescat(z, params)

Choose the type of scattering event, depending on the altitude `z`.
"""
@inline function choosescat(r, world::World, params::Params)::ScatteringType
    @unpack νraymax, nair, H, σray, νray_ground = params
    νmax = νraymax + νmiemax(world.comp)

    inside(world.domain, r) || return Null()

    if inside(world.cloud, r)
        loc = probe(world.comp, r)
        Qext = mie_qext(world.comp, r, loc)
        νMie =  Qext * π * radius(world.comp, r, loc)^2 * density(world.comp, r, loc)
    else
        νMie = 0.0
    end

    νRay = νray_ground * exp(-r[3] / H)

    ξ = trand() * νmax

    if ξ <= νMie
        g = mie_g(world.comp, r, loc)
        ω0 = mie_ω0(world.comp, r, loc)
        return Mie(g, ω0)
    elseif ξ <= (νMie + νRay)
        return Rayleigh()
    else
        return Null()
    end
end

"""
Given a unitary vector μ, sample a scattering angle from the HG
distribution and finds a new vector that forms that angle with μ,
which is returned.

NB: This function breaks down if μ is directed exactly along z.
"""
@inline @fastmath function scatterone(μ, w, scat::ScatteringType, params)
    @unpack weight_min = params
    scatters(scat) || (return μ, w)

    # Return false if the particle is absorbed
    w = reweight(scat, w, weight_min)

    ϕ = 2π * trand()
    cosθ = sample(scat, trand())

    turn(μ, cosθ, ϕ), w
end

"""
    miecollrate(r, w, params)

Computes Mie collision rate at a given point.
"""
function miecollrate(r::Point, w::World, params::Params)
    inside(w.cloud, r) || return zero(Float64)
    loc = probe(w.comp, r)

    # Is the computation of ω0, g optimized out?
    Qext = mie_qext(w.comp, r, loc)
    Qext * π * radius(w.comp, r, loc)^2 * density(w.comp, r, loc)
end

"""
    turn(μ, cosθ, ϕ)

Find a unitary vector deviated with inclination θ and azimuth
ϕ with respect to μ.  For performance, the passed parameter is cosθ instead of θ.
"""
@inline @fastmath function turn(μ, cosθ, ϕ)
    sinϕ, cosϕ = sincos(ϕ)

    sinθ = sqrt(1 - cosθ^2)

    s = sqrt(1 - μ[3]^2)

    b(μ1, μ2, μ3) = μ1 * μ3 * cosϕ + μ2 * sinϕ

    μ1x = sinθ * b(μ[1], -μ[2], μ[3]) / s + μ[1] * cosθ
    μ1y = sinθ * b(μ[2],  μ[1], μ[3]) / s + μ[2] * cosθ
    μ1z = -s * sinθ * cosϕ + μ[3] * cosθ

    @SVector [μ1x, μ1y, μ1z]
end

const INTERFACE_LISTS = Vector{Float64}[]

@noinline function init_interface_list()
    # Allocate the list on the thread's own heap lazily
    # instead of the master thread heap to minimize memory conflict.
    tid = Threads.threadid()
    INTERFACE_LISTS[tid] = Float64[0.0]
    INTERFACE_LISTS[tid]
end

@inline function get_interface_list()
    @inbounds begin
        tid = Threads.threadid()
        isassigned(INTERFACE_LISTS, tid) || return init_interface_list()
        resize!(INTERFACE_LISTS[tid], 1)
        INTERFACE_LISTS[tid]
    end
end

"""
Compute the optical path between points `a` and `b` in a given geometry
`geom` and with a position-dependent collision rate `ν`.
"""
function optpath(geom, ν::Function, a::Point, b::Point, quadrule)
    v = ListAccumulator(get_interface_list())

    interfaces!(v, geom, a, b)
    sort!(v.list)

    push!(v.list, 1.0)

    f(s) = ν((1 - s) * a + s * b)

    L = norm(a - b)
    # path = L * quadgk(f, 0., v..., 1.0,
    #                   atol=1e-3, rtol=1e-2, order=3)[1]

    path = L * quadgauss(f, v.list, quadrule)
    path
end

@fastmath function quadgauss(f::Function, x, quadrule)
    xq, wq = quadrule
    x1 = x[1]

    I = zero(x1)
    for ix in 2:length(x)
        h = 0.5 * (x[ix] - x1)
        Ip = zero(x1)

        for i in 1:length(wq)
            @inbounds Ip += wq[i] * f(x1 + (1 + xq[i]) * h)
        end
        x1 = x[ix]
        I += h * Ip
    end
    I
end

function depthimg(o::Observer, w::World, params::Params)
    # This must be a large distance to be sure that the full cloud is closer
    # than this to the observer.
    INF_DIST = 1e6
    nx, ny, _ = size(o.img)
    dimg = zeros(nx, ny)

    for i in 1:nx, j in 1:ny
        μ0 = @SVector [(i + 0.5) * o.δu - o.umax, (j + 0.5) * o.δu - o.umax , -1.0]
        μ = μ0 ./ norm(μ0)
        rfar = o.r .+ INF_DIST .* μ
        dimg[i, j] = optpath(w.cloud, r -> miecollrate(r, w, params), rfar, o.r,
                             w.quadrule)
    end

    dimg
end

"""
Build a function that returns a fixed number `ν0` inside the geometry
`geom` and zero otherwise.
"""
function constinside(geom, ν0)
    function ν(r::Point)
        inside(geom, r) ? ν0 : zero(ν0)
    end
end

const RNG = MersenneTwister[]

@noinline function init_thread_rng()
    # Allocate the random number generator on the thread's own heap lazily
    # instead of the master thread heap to minimize memory conflict.
    RNG[Threads.threadid()] = MersenneTwister(Threads.threadid())
end

@inline function trand()
    @inbounds begin
        tid = Threads.threadid()
        isassigned(RNG, tid) || init_thread_rng()
        return rand(RNG[tid])
    end
end

function __init__()
    nth = Threads.nthreads()
    resize!(RNG, nth)
    resize!(INTERFACE_LISTS, nth)

    # for tid in 1:nth
    #     RNG[tid] = MersenneTwister(tid)
    # end
end

end # module


text:
@testset "Small special tests" begin
    @testset "Fix variable" begin
        m = Model(CSJuMPTestOptimizer())
        @variable(m, 1 <= x <= 9, Int)
        @variable(m, y == 2, Int)
        # should just return optimal with any 1-9 for x and y is fixed
        optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test 1 <= JuMP.value(x) <= 9 && length(CS.values(m, x)) == 1
        @test JuMP.value(y) == 2

        m = Model(CSJuMPTestOptimizer())
        @variable(m, 1 <= x <= 9, Int)
        @variable(m, y == 2, Int)
        @constraint(m, x + y == 10)
        optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.value(x) == 8
        @test JuMP.value(y) == 2
    end

    @testset "LessThan constraints JuMP" begin
        m = Model(CSJuMPTestOptimizer(; branch_strategy = :ABS))
        @variable(m, 1 <= x[1:5] <= 9, Int)
        @constraint(m, sum(x) <= 25)
        @constraint(m, sum(x) >= 20)
        weights = [1, 2, 3, 4, 5]
        @objective(m, Max, sum(weights .* x))
        optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.value.(x) == [1, 1, 5, 9, 9]
        @test JuMP.objective_value(m) == 99

        # minimize
        m = Model(CSJuMPTestOptimizer())
        @variable(m, 1 <= x[1:5] <= 9, Int)
        @constraint(m, sum(x) <= 25)
        @constraint(m, sum(x) >= 20)
        weights = [1, 2, 3, 4, 5]
        @objective(m, Min, sum(weights .* x))
        optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.value.(x) == [9, 8, 1, 1, 1]
        @test JuMP.objective_value(m) == 37

        # minimize with negative and positive real weights
        m = Model(CSJuMPTestOptimizer(; branch_strategy = :ABS))
        @variable(m, 1 <= x[1:5] <= 9, Int)
        @constraint(m, sum(x) <= 25)
        weights = [-0.1, 0.2, -0.3, 0.4, 0.5]
        @constraint(m, sum(x[i] for i = 1:5 if weights[i] > 0) >= 15)
        @objective(m, Min, sum(weights .* x))
        optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.value.(x) == [1, 9, 9, 5, 1]
        @test JuMP.objective_value(m) ≈ 1.5
    end

    @testset "Knapsack problems" begin
        m = Model(CSJuMPTestOptimizer())

        @variable(m, 1 <= x[1:5] <= 9, Int)
        @constraint(m, sum(x) <= 25)
        @constraint(m, x[2] + 1.2 * x[4] - x[5] <= 12)
        weights = [1.2, 3.0, -0.3, -5.2, 2.7]
        @objective(m, Max, dot(weights, x))

        optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL
        x_vals = JuMP.value.(x)
        @test x_vals[1] ≈ 5
        @test x_vals[2] ≈ 9
        @test x_vals[3] ≈ 1
        @test x_vals[4] ≈ 1
        @test x_vals[5] ≈ 9
        @test JuMP.objective_value(m) ≈ 51.8

        # less variables in the objective
        m = Model(CSJuMPTestOptimizer())

        @variable(m, 1 <= x[1:5] <= 9, Int)
        @constraint(m, sum(x) <= 25)
        @constraint(m, x[2] + 1.2 * x[4] - x[5] <= 12)
        @constraint(m, -x[3] - 1.2 * x[4] + x[5] <= 12)
        weights = [1.2, 3.0, -0.3, -5.2, 2.7]
        @objective(m, Max, x[3] + 2.7 * x[4] - x[1])

        optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL
        x_vals = JuMP.value.(x)
        @test sum(x_vals) <= 25
        @test x_vals[2] + 1.2 * x_vals[4] - x_vals[5] <= 12
        @test JuMP.objective_value(m) ≈ 32.3

        # minimize
        m = Model(CSJuMPTestOptimizer())

        @variable(m, 1 <= x[1:5] <= 9, Int)
        @constraint(m, sum(x) >= 25)
        @constraint(m, x[2] + 1.2 * x[4] >= 12)
        weights = [1.2, 3.0, 0.3, 5.2, 2.7]
        @objective(m, Min, dot(weights, x))

        optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL
        x_vals = JuMP.value.(x)
        @test x_vals[1] ≈ 3
        @test x_vals[2] ≈ 9
        @test x_vals[3] ≈ 9
        @test x_vals[4] ≈ 3
        @test x_vals[5] ≈ 1
        @test JuMP.objective_value(m) ≈ 51.6

        # minimize only part of the weights and some are negative
        m = Model(CSJuMPTestOptimizer())

        @variable(m, 1 <= x[1:5] <= 9, Int)
        @constraint(m, sum(x) >= 25)
        @constraint(m, x[2] + 1.2 * x[4] - x[1] >= 12)
        @constraint(m, x[5] <= 7)
        @objective(m, Min, 3 * x[2] + 5 * x[1] - 2 * x[3])

        optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL
        x_vals = JuMP.value.(x)
        @test x_vals[1] ≈ 1
        @test x_vals[2] ≈ 3
        @test x_vals[3] ≈ 9
        @test x_vals[4] ≈ 9
        @test sum(x_vals) >= 25
        @test JuMP.objective_value(m) ≈ -4

        m = Model(CSJuMPTestOptimizer())
        @variable(m, 1 <= x[1:5] <= 9, Int)
        @constraint(m, sum(x) <= 25)
        @constraint(m, -x[1] - x[2] - x[3] + x[4] + x[5] >= 5)
        weights = [-1, 2, 3, 4, 5]
        @objective(m, Min, sum(weights[1:3] .* x[1:3]))

        optimize!(m)

        x_vals = JuMP.value.(x)
        @test sum(x_vals) <= 25
        @test -x_vals[1] - x_vals[2] - x_vals[3] + x_vals[4] + x_vals[5] >= 5
        @test JuMP.objective_value(m) ≈ -3
    end

    @testset "Not supported constraints" begin
        m = Model(CSJuMPTestOptimizer())
        # must be an Integer upper bound
        @variable(m, 1 <= x[1:5] <= NaN, Int)
        @test_throws ErrorException optimize!(m)

        m = Model(CSJuMPTestOptimizer())
        # must be an Integer lower bound
        @variable(m, NaN <= x[1:5] <= 2, Int)
        @test_throws ErrorException optimize!(m)

        m = Model(CSJuMPTestOptimizer())
        @variable(m, 1 <= x[1:5] <= 2, Int)

        m = Model(CSJuMPTestOptimizer())
        @variable(m, 1 <= x[1:5] <= 2, Int)

        m = Model(CSJuMPTestOptimizer())
        @variable(m, 1 <= x[1:5] <= 2, Int)
    end

    @testset "Bipartite matching" begin
        match = CS.bipartite_cardinality_matching([2, 1, 3], [1, 2, 3], 3, 3)
        @test match.weight == 3
        @test match.match == [2, 1, 3]

        # no perfect matching
        match = CS.bipartite_cardinality_matching(
            [1, 2, 3, 4, 1, 2, 3, 3],
            [1, 1, 2, 2, 2, 2, 3, 4],
            4,
            4,
        )
        @test match.weight == 3
        # 4 is zero and the rest should be different
        @test match.match[4] == 0
        @test allunique(match.match)

        # more values than indices
        match = CS.bipartite_cardinality_matching(
            [1, 2, 3, 4, 1, 2, 3, 3, 2, 1, 2],
            [1, 1, 2, 2, 2, 2, 3, 4, 5, 5, 6],
            4,
            6,
        )
        @test match.weight == 4
        # all should be matched to different values
        @test allunique(match.match)
        # no unmatched vertex
        @test count(i -> i == 0, match.match) == 0

        # more values than indices with matching_init
        m = 4
        n = 6
        l = [1, 2, 3, 4, 1, 2, 3, 3, 2, 1, 2, 0, 0]
        r = [1, 1, 2, 2, 2, 2, 3, 4, 5, 5, 6, 0, 0]
        # don't use the zeros
        l_len = length(l) - 2
        matching_init = CS.MatchingInit(
            l_len,
            zeros(Int, m),
            zeros(Int, n),
            zeros(Int, m + 1),
            zeros(Int, m + n),
            zeros(Int, m + n),
            zeros(Int, m + n),
            zeros(Bool, m),
            zeros(Bool, n),
        )
        match = CS.bipartite_cardinality_matching(l, r, m, n; matching_init = matching_init)
        @test match.weight == 4
        # all should be matched to different values
        @test allunique(match.match)
        # no unmatched vertex
        @test count(i -> i == 0, match.match) == 0
    end

    @testset "Not equal" begin
        m = Model(CSJuMPTestOptimizer())

        @variable(m, 1 <= x <= 10, Int)
        @variable(m, 1 <= y <= 1, Int)
        @variable(m, 1 <= z <= 10, Int)
        @constraint(m, x != 2 - 1) # != 1
        @constraint(m, 2x != 4) # != 2
        @constraint(m, π / 3 * x != π) # != 3
        @constraint(m, 2.2x != 8.8) # != 4
        @constraint(m, 4x + 5y != 25) # != 5
        @constraint(m, 4x + π * y != 10) # just some random stuff
        @constraint(m, x + y + z - π != 10)
        @constraint(m, x + y + z + 2 != 10)
        @objective(m, Min, x)
        optimize!(m)

        @test JuMP.objective_value(m) == 6
        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.value(x) == 6
        @test JuMP.value(y) == 1
        # the values should be fixed
        @test length(CS.values(m, x)) == 1
        @test length(CS.values(m, y)) == 1
        @test length(CS.values(m, z)) == 1
        @test JuMP.value(x) + JuMP.value(y) + JuMP.value(z) + 2 != 10
    end

    @testset "Integers basic" begin
        m = Model(CSJuMPTestOptimizer())
        @variable(m, x, CS.Integers([1, 2, 4]))
        @variable(m, y, CS.Integers([2, 3, 5, 6]))
        @constraint(m, x == y)
        @objective(m, Max, x)
        optimize!(m)
        @test JuMP.value(x) ≈ 2
        @test JuMP.value(y) ≈ 2
        @test JuMP.objective_value(m) ≈ 2

        m = Model(optimizer_with_attributes(CS.Optimizer, "backtrack" => false))
        @variable(m, x, CS.Integers([1, 2, 4]))
        optimize!(m)
        com = CS.get_inner_model(m)
        @test !CS.has(com.search_space[1], 3)
        @test sort(CS.values(com.search_space[1])) == [1, 2, 4]

        m = Model(optimizer_with_attributes(CS.Optimizer, "backtrack" => false))
        @variable(m, y, CS.Integers([2, 5, 6, 3]))
        optimize!(m)
        com = CS.get_inner_model(m)
        @test !CS.has(com.search_space[1], 1)
        @test !CS.has(com.search_space[1], 4)
        @test sort(CS.values(com.search_space[1])) == [2, 3, 5, 6]
    end

    @testset "Biggest cube square number up to 100" begin
        m = Model(CSJuMPTestOptimizer())
        @variable(m, x, CS.Integers([i^2 for i = 1:20 if i^2 < 100]))
        @variable(m, y, CS.Integers([i^3 for i = 1:20 if i^3 < 100]))
        @constraint(m, x == y)
        @objective(m, Max, x)
        optimize!(m)
        @test JuMP.value(x) ≈ 64
        @test JuMP.value(y) ≈ 64
        @test JuMP.objective_value(m) ≈ 64
    end

    @testset "Pythagorean triples" begin
        m = Model(optimizer_with_attributes(
            CS.Optimizer,
            "all_solutions" => true,
            "branch_strategy" => :ABS,
            "activity.decay" => 0.999,
            "activity.max_probes" => 20,
            "logging" => [],
        ))
        @variable(m, x[1:3], CS.Integers([i^2 for i in 1:50]))
        @constraint(m, x[1] + x[2] == x[3])
        @constraint(m, x[1] <= x[2])
        optimize!(m)
        com = CS.get_inner_model(m)
        @test is_solved(com)
        @test MOI.get(m, MOI.ResultCount()) == 20
    end

    @testset "Infeasible by fixing variable to outside domain" begin
        m = Model(CSJuMPTestOptimizer())
        @variable(m, x, CS.Integers([1, 2, 4]))
        @constraint(m, x == 3)
        optimize!(m)
        @test JuMP.termination_status(m) == MOI.INFEASIBLE
    end

    @testset "Infeasible by fixing variable to two values" begin
        m = Model(CSJuMPTestOptimizer())
        @variable(m, x, CS.Integers([1, 2, 4]))
        @constraint(m, x == 1)
        @constraint(m, x == 2)
        optimize!(m)
        @test JuMP.termination_status(m) == MOI.INFEASIBLE
    end

    @testset "5 variables all equal" begin
        m = Model(optimizer_with_attributes(
            CS.Optimizer,
            "all_solutions" => true,
            "logging" => [],
        ))

        @variable(m, 5 <= x <= 10, Int)
        @variable(m, 2 <= y <= 15, Int)
        @variable(m, 1 <= z <= 7, Int)
        @variable(m, 2 <= a <= 9, Int)
        @variable(m, 6 <= b <= 10, Int)
        @constraint(m, x == y)
        # should not result in linking to x -> y -> x ...
        @constraint(m, y == x)
        @constraint(m, x == y)

        @constraint(m, y == z)
        @constraint(m, a == z)
        @constraint(m, b == y)
        optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.value(x) ==
              JuMP.value(y) ==
              JuMP.value(z) ==
              JuMP.value(a) ==
              JuMP.value(b)
        @test JuMP.value(x; result = 2) ==
              JuMP.value(y; result = 2) ==
              JuMP.value(z; result = 2) ==
              JuMP.value(a; result = 2) ==
              JuMP.value(b; result = 2)
        @test JuMP.value(x) == 6 || JuMP.value(x) == 7
        @test JuMP.value(x; result = 2) == 6 || JuMP.value(x; result = 2) == 7
        @test JuMP.value(x) != JuMP.value(x; result = 2)
    end

    @testset "x[1] <= x[1]" begin
        model = Model(CSJuMPTestOptimizer())
        n = 4
        @variable(model, 1 <= x[1:n] <= n, Int)
        for i in 1:n
            @constraint(model, x[1] <= x[i])
        end

        optimize!(model)
        status = JuMP.termination_status(model)
        @test status == MOI.OPTIMAL
        com = CS.get_inner_model(model)
        @test is_solved(com)
    end

    @testset "x[1] <= x[1] - 1 " begin
        model = Model(CSJuMPTestOptimizer())
        n = 4
        @variable(model, 1 <= x[1:n] <= n, Int)
        for i in 1:n
            @constraint(model, x[1] <= x[i] - 1)
        end

        optimize!(model)
        status = JuMP.termination_status(model)
        @test status == MOI.INFEASIBLE
    end

    @testset "x[1] == x[1]" begin
        model = Model(CSJuMPTestOptimizer())
        n = 4
        @variable(model, 1 <= x[1:n] <= n, Int)
        for i in 1:n
            @constraint(model, x[1] == x[i])
        end

        optimize!(model)
        status = JuMP.termination_status(model)
        @test status == MOI.OPTIMAL
        com = CS.get_inner_model(model)
        @test is_solved(com)
    end

    @testset "x[1] == x[1] - 1 " begin
        model = Model(CSJuMPTestOptimizer())
        n = 4
        @variable(model, 1 <= x[1:n] <= n, Int)
        for i in 1:n
            @constraint(model, x[1] == x[i] - 1)
        end

        optimize!(model)
        status = JuMP.termination_status(model)
        @test status == MOI.INFEASIBLE
    end

    @testset "x[1] == x[1] - 1 in reified" begin
        model = Model(CSCbcJuMPTestOptimizer())
        n = 4
        @variable(model, 1 <= x[1:n] <= n, Int)
        @variable(model, b[1:n], Bin)
        for i in 1:n
            @constraint(model, b[i] := {x[1] == x[i] - 1})
        end
        @objective(model, Max, sum(b))
        optimize!(model)
        status = JuMP.termination_status(model)
        @test status == MOI.OPTIMAL
        @test JuMP.objective_value(model) ≈ 3
        @test JuMP.value(b[1]) ≈ 0
    end

    @testset "x[1] == x[1] - 1 in indicator" begin
        model = Model(CSCbcJuMPTestOptimizer())
        n = 4
        @variable(model, 1 <= x[1:n] <= n, Int)
        @variable(model, b[1:n], Bin)
        for i in 1:n
            @constraint(model, b[i] => {x[1] == x[i] - 1})
        end
        @objective(model, Max, sum(b))
        optimize!(model)
        status = JuMP.termination_status(model)
        @test status == MOI.OPTIMAL
        @test JuMP.objective_value(model) ≈ 3
        @test JuMP.value(b[1]) ≈ 0
    end

    @testset "Infeasible all different in indicator" begin
        model = Model(CSCbcJuMPTestOptimizer(; branch_strategy = :ABS))
        n = 2
        @variable(model, 1 <= x[1:n] <= n - 1, Int)
        @variable(model, b, Bin)
        @constraint(model, b => {x in CS.AllDifferent()})
        @objective(model, Max, b)
        optimize!(model)
        status = JuMP.termination_status(model)
        @test status == MOI.OPTIMAL
        @test JuMP.objective_value(model) ≈ 0
        @test JuMP.value(b) ≈ 0
    end

    @testset "Infeasible all different in reified" begin
        model = Model(CSJuMPTestOptimizer())
        n = 4
        @variable(model, 1 <= x[1:n] <= n - 1, Int)
        @variable(model, b, Bin)
        @constraint(model, b := {x in CS.AllDifferent()})
        @objective(model, Max, b)
        optimize!(model)
        status = JuMP.termination_status(model)
        @test status == MOI.OPTIMAL
        @test JuMP.objective_value(model) ≈ 0
        @test JuMP.value(b) ≈ 0
    end

    @testset "Table dim not matching" begin
        model = Model(CSJuMPTestOptimizer())
        n = 4
        @variable(model, 1 <= x[1:n] <= n - 1, Int)
        @variable(model, b, Bin)
        @test_throws ArgumentError @constraint(model, b => {x in CS.TableSet([
            10 20
            11 5
        ])})
    end

    @testset "Infeasible table in indicator" begin
        model = Model(CSJuMPTestOptimizer())
        n = 2
        @variable(model, 1 <= x[1:n] <= n - 1, Int)
        @variable(model, b, Bin)
        @constraint(model, b => {x in CS.TableSet([
            10 20
            11 5
        ])})
        @objective(model, Max, b)
        optimize!(model)
        status = JuMP.termination_status(model)
        @test status == MOI.OPTIMAL
        @test JuMP.objective_value(model) ≈ 0
        @test JuMP.value(b) ≈ 0
    end

    @testset "Infeasible table in reified" begin
        model = Model(CSJuMPTestOptimizer())
        n = 2
        @variable(model, 1 <= x[1:n] <= n - 1, Int)
        @variable(model, b, Bin)
        @constraint(model, b := {x in CS.TableSet([
            10 20
            11 5
        ])})
        @objective(model, Max, b)
        optimize!(model)
        status = JuMP.termination_status(model)
        @test status == MOI.OPTIMAL
        @test JuMP.objective_value(model) ≈ 0
        @test JuMP.value(b) ≈ 0
    end

    @testset "Magic Square 5x5" begin
        n = 5
        model = Model(optimizer_with_attributes(
            CS.Optimizer,
            "traverse_strategy" => :BFS,
            "logging" => [],
            "branch_split" => :InHalf,
            "time_limit" => 3,
        ))

        # The total for each row, column, and the two main diaginals
        s = round(Int, n * (n^2 + 1) / 2)
        @variable(model, 1 <= x[1:n, 1:n] <= n^2, Int)
        @constraint(model, x[:] in CS.AllDifferent())

        for i in 1:n
            # Rows
            @constraint(model, sum(x[i, :]) == s)

            # Columns
            @constraint(model, sum(x[:, i]) == s)
        end

        # diagonals
        @constraint(model, sum([x[i, i] for i in 1:n]) == s)
        @constraint(model, s == sum([x[i, n - i + 1] for i in 1:n]))

        optimize!(model)
        status = JuMP.termination_status(model)
        @test status == MOI.OPTIMAL
        sol = convert.(Int, JuMP.value.(x))
        for i in 1:n
            @test sum(sol[i, :]) == 65
            @test sum(sol[:, i]) == 65
        end
        @test sum([sol[i, i] for i in 1:n]) == 65
        @test sum([sol[i, n - i + 1] for i in 1:n]) == 65
        @test allunique(sol)
    end

    @testset "Indicator with > and active" begin
        cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
        model = Model(optimizer_with_attributes(
            CS.Optimizer,
            "lp_optimizer" => cbc_optimizer,
            "logging" => []
        ))
        @variable(model, b >= 1, Bin)
        @variable(model, 0 <= x[1:4] <= 5, Int)
        @constraint(model, b => {sum([0.4,0.5,0.7,0.8] .* x) > 9})
        @objective(model, Min, sum([0.4,0.5,0.7,0.8] .* x))
        optimize!(model)
        CS.get_inner_model(model)
        @test JuMP.termination_status(model) == MOI.OPTIMAL
        @test JuMP.objective_value(model) ≈ 9.1
        @test JuMP.value(x[1]) ≈ 2.0
        @test JuMP.value(x[2]) ≈ 3.0
        @test JuMP.value(x[3]) ≈ 4.0
        @test JuMP.value(x[4]) ≈ 5.0
    end

    @testset "AllDifferent except 0" begin
        m = Model(optimizer_with_attributes(
            CS.Optimizer,
            "logging" => []
        ))
        @variable(m, 0 <= x[1:5] <= 3, Int)
        len = length(x)
        for i in 2:len, j in 1:i-1
            b = @variable(m, binary=true)
            @constraint(m, b := {x[i] != 0 && x[j] != 0})
            @constraint(m, b => {x[i] != x[j]} )
        end

        @objective(m, Min, sum(x))
        optimize!(m)
        CS.get_inner_model(m)
        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.objective_value(m) ≈ 0
        @test JuMP.value(x[1]) ≈ 0
        @test JuMP.value(x[2]) ≈ 0
        @test JuMP.value(x[3]) ≈ 0
        @test JuMP.value(x[4]) ≈ 0
        @test JuMP.value(x[5]) ≈ 0

        # maximize
        m = Model(optimizer_with_attributes(
            CS.Optimizer,
            "logging" => []
        ))
        @variable(m, 0 <= x[1:5] <= 3, Int)
        len = length(x)
        for i in 2:len, j in 1:i-1
            b = @variable(m, binary=true)
            @constraint(m, b := {x[i] != 0 && x[j] != 0})
            @constraint(m, b => {x[i] != x[j]} )
        end

        @objective(m, Max, sum(x))
        optimize!(m)
        CS.get_inner_model(m)
        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.objective_value(m) ≈ 3+2+1+0+0
    end

    @testset "Try to fulfill b with several conditions" begin
        # inner part is infeasible
        m = Model(optimizer_with_attributes(
            CS.Optimizer,
            "logging" => []
        ))
        @variable(m, 0 <= x[1:5] <= 3, Int)
        @variable(m, b, Bin)
        @constraint(m, b => {x[1] > x[2] && x[2] > x[3] && x[3] > x[4] && x[4] > x[5]})

        @objective(m, Max, b)
        optimize!(m)
        CS.get_inner_model(m)
        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.objective_value(m) ≈ 0
        @test JuMP.value(b) ≈ 0

        # inner part is feasible
        m = Model(optimizer_with_attributes(
            CS.Optimizer,
            "logging" => []
        ))
        @variable(m, 0 <= x[1:5] <= 4, Int)
        @variable(m, b, Bin)
        @constraint(m, b => {x[1] > x[2] && x[2] > x[3] && x[3] > x[4] && x[4] > x[5]})

        @objective(m, Max, b)
        optimize!(m)
        CS.get_inner_model(m)
        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.objective_value(m) ≈ 1
        @test JuMP.value(b) ≈ 1
        vx = JuMP.value.(x)
        @test vx[1] > vx[2] && vx[2] > vx[3] && vx[3] > vx[4] && vx[4] > vx[5]

        # inner part is feasible
        m = Model(optimizer_with_attributes(
            CS.Optimizer,
            "logging" => []
        ))
        @variable(m, 0 <= x[1:5] <= 4, Int)
        @variable(m, b, Bin)
        @constraint(m, b => {x[1] >= x[2] && (x[3] <= x[2] && (x[3] >= x[4] && x[4] > x[5])) && x in CS.AllDifferent()})

        @objective(m, Max, b)
        optimize!(m)
        CS.get_inner_model(m)
        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.objective_value(m) ≈ 1
        @test JuMP.value(b) ≈ 1
        vx = JuMP.value.(x)
        @test vx[1] > vx[2] && vx[2] > vx[3] && vx[3] > vx[4] && vx[4] > vx[5]

        # inner part is feasible have AndSet && AndSet
        m = Model(optimizer_with_attributes(
            CS.Optimizer,
            "logging" => []
        ))
        @variable(m, 0 <= x[1:5] <= 4, Int)
        @variable(m, b, Bin)
        @constraint(m, b => {(x[1] >= x[2] && x[3] <= x[2]) && (x[3] >= x[4] && x[4] > x[5] && x in CS.AllDifferent())})

        @objective(m, Max, b)
        optimize!(m)
        CS.get_inner_model(m)
        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.objective_value(m) ≈ 1
        @test JuMP.value(b) ≈ 1
        vx = JuMP.value.(x)
        @test vx[1] > vx[2] && vx[2] > vx[3] && vx[3] > vx[4] && vx[4] > vx[5]
    end

    @testset "Table in reified where setting activator to false" begin
        m = Model(optimizer_with_attributes(
            CS.Optimizer,
            "logging" => []
        ))
        @variable(m, 1 <= x[1:5] <= 4, Int)
        @variable(m, b, Bin)
        @constraint(m, b => {x[1] > x[2] && [x[1],x[2]] in CS.TableSet([1 2; 1 3])})

        @objective(m, Max, 10b+sum(x))
        optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.objective_value(m) ≈ 20
        @test JuMP.value(b) ≈ 0
    end

    @testset "Small test case for || in reified" begin
        m = Model(optimizer_with_attributes(
            CS.Optimizer,
            "logging" => []
        ))
        @variable(m, 1 <= x[1:5] <= 4, Int)
        @variable(m, b, Bin)
        @constraint(m, b := {(x[1] > x[2] || x in CS.AllDifferent()) })

        @objective(m, Max, 100b+sum(x))
        optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL

        @test JuMP.value(b) ≈ 1.0
        @test JuMP.value.(x) ≈ [4,3,4,4,4]
        @test JuMP.objective_value(m) ≈ sum([4,3,4,4,4]) + 100
    end

    @testset "Small test case for || and && in reified" begin
        m = Model(optimizer_with_attributes(
            CS.Optimizer,
            "logging" => []
        ))
        @variable(m, 1 <= x[1:5] <= 4, Int)
        @variable(m, b, Bin)
        table = [
            4 0;
            3 2;
        ]
        @constraint(m, b := {(x[1] > x[2] || x in CS.AllDifferent()) && x[4:5] in CS.TableSet(table) })

        @objective(m, Max, 100b+sum(x))
        optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL

        @test JuMP.value(b) ≈ 1.0
        @test JuMP.value.(x) ≈ [4,3,4,3,2]
        @test JuMP.objective_value(m) ≈ sum([4,3,4,3,2]) + 100
    end

    @testset "Small test case for || and &&" begin
        m = Model(optimizer_with_attributes(
            CS.Optimizer,
            "logging" => []
        ))
        @variable(m, 1 <= x[1:5] <= 4, Int)
        table = [
            4 0;
            3 2;
        ]
        @constraint(m, (x[1] > x[2] || x in CS.AllDifferent()) && x[4:5] in CS.TableSet(table))

        @objective(m, Max, sum(x))
        optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL

        @test JuMP.value.(x) ≈ [4,3,4,3,2]
        @test JuMP.objective_value(m) ≈ sum([4,3,4,3,2])
    end

    @testset "Small test case for && and TableSet" begin
        m = Model(optimizer_with_attributes(
            CS.Optimizer,
            "logging" => [],
            "all_solutions" => true
        ))
        @variable(m, x[1:2], Bin)
        @variable(m, y[1:2], Bin)
        xor = [
            0 1;
            1 0;
        ]
        and = [
            1 1;
        ]
        @constraint(m, (x in CS.TableSet(xor) || y in CS.TableSet(and)))

        optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL

        # x = 0 1 | y = 00 / 01 / 10 / 11
        # x = 1 0 | y = 00 / 01 / 10 / 11
        # y = 1 1 | x = 00 / 11 <- without double counting
        num_sols = MOI.get(m, MOI.ResultCount())
        @test num_sols == 10

        sets = [
            [0, 1, 0, 0],
            [0, 1, 0, 1],
            [0, 1, 1, 0],
            [0, 1, 1, 1],
            [1, 0, 0, 0],
            [1, 0, 0, 1],
            [1, 0, 1, 0],
            [1, 0, 1, 1],
            [0, 0, 1, 1],
            [1, 1, 1, 1]
        ]
        for i in 1:10
            found = false
            for j in 1:10
                if JuMP.value.([x...,y...]; result=j) == sets[i]
                    found = true
                    break
                end
            end
            @test found
        end
    end

    @testset "Monks and Doors" begin
        model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> true,
                "logging"=>[],
                "traverse_strategy"=>:BFS,
                "branch_split"=>:InHalf,
                "branch_strategy" => :ABS,
        ))

        num_doors = 4
        num_monks = 8
        @variable(model, doors[1:num_doors], Bin)
        da,db,dc,dd = doors
        door_names = ["A","B","C","D"]

        @variable(model, monks[1:num_monks], Bin)
        m1,m2,m3,m4,m5,m6,m7,m8 = monks

        # Monk 1: Door A is the exit.
        # M1 #= A (Picat constraint)
        @constraint(model, m1 == da)

        #  Monk 2: At least one of the doors B and C is the exit.
        # M2 #= 1 #<=> (B #\/ C)
        # @constraint(model, m2 := { db + dc >= 1})
        @constraint(model, m2 := { [db, dc] in CS.TableSet([
            1 0;
        ]) || dc == 1})

        #  Monk 3: Monk 1 and Monk 2 are telling the truth.
        # M3 #= 1 #<=> (M1 #/\ M2)
        # @constraint(model, m3 := { m1 + m2 == 2})
        # @constraint(model, m3 := { m1 == 1 && m2 == 1})
        @constraint(model, m3 := { m1 && m2})

        #  Monk 4: Doors A and B are both exits.
        # M4 #= 1 #<=> (A #/\ B)
        # @constraint(model, m4 := { da + db == 2})
        @constraint(model, m4 := { da && db })

        #  Monk 5: Doors A and C are both exits.
        # M5 #= 1 #<=> (A #/\ C)
        # @constraint(model, m5 := { da + dc == 2})
        @constraint(model, m5 := { da && dc })

        #  Monk 6: Either Monk 4 or Monk 5 is telling the truth.
        # M6 #= 1 #<=> (M4 #\/ M5)
        # @constraint(model, m6 := { m4 + m5 == 1})
        @constraint(model, m6 := { m4 || m5 })

        #  Monk 7: If Monk 3 is telling the truth, so is Monk 6.
        # M7 #= 1 #<=> (M3 #=> M6)
        @constraint(model, m7 := { [m3, m6] in CS.TableSet([1 1; 0 0; 0 1;]) && [m3, m6] in CS.TableSet([1 1; 0 0; 0 1;]) })
        # @constraint(model, m7 := { m3 => m6})  # This don't work!
        # @constraint(model, m7 := { m3 => {m6 == 1}}) # This don't work either.

        #  Monk 8: If Monk 7 and Monk 8 are telling the truth, so is Monk 1.
        # M8 #= 1 #<=> ((M7 #/\ M8) #=> (M1))
        b1 = @variable(model, binary=true)
        # @constraint(model, b1 := {m7 + m8 == 2})
        @constraint(model, b1 := {m7 && [m7,m8] in CS.TableSet([1 1; 2 2; 3 3;])})
        @constraint(model, m8 := {b1 <= m1})

        # Exactly one door is an exit.
        # (A + B + C + D) #= 1
        @constraint(model, da + db + dc + dd == 1)

        optimize!(model)

        status = JuMP.termination_status(model)
        @test JuMP.termination_status(model) == MOI.OPTIMAL

        @test MOI.get(model, MOI.ResultCount()) == 1
        doors_val = convert.(Integer,JuMP.value.(doors))
        monks_val = convert.(Integer,JuMP.value.(monks))
        @test doors_val == [1, 0, 0, 0]
        @test monks_val == [1, 0, 0, 0, 0, 0, 1, 1]

        com = CS.get_inner_model(model)
        for var in com.search_space
            @test issorted(var.init_vals)
        end
    end
end


text:
export MCVanilla, Vegas, Domain
export integral, ∫

using QuickTypes: @qstruct
using ArgCheck
using LinearAlgebra, Statistics
using StaticArrays
using Base.Threads: @threads
using Setfield: @settable

import Random

abstract type MCAlgorithm end

@settable @qstruct MCVanilla{R}(
        neval::Int64=10^6,
        rng::R=GLOBAL_PAR_RNG,
    ) do
        @argcheck neval > 2
end <: MCAlgorithm

"""
    Domain

Encodes a domain of integration, i.e. a product of finte intervals.
"""
struct Domain{N,T}
    lower::SVector{N,T}
    upper::SVector{N,T}
end

# function cartesian_product(d1::Domain, d2::Domain)
#     Domain(vcat(d1.lower, d2.lower), vcat(d1.upper, d2.upper))
# end

function pointtype(::Type{Domain{N,T}}) where {N,T}
    SVector{N,T}
end

function volume(dom::Domain)
    prod(dom.upper - dom.lower)
end

function Base.ndims(dom::Domain{N,T}) where {N,T}
    N
end

function Domain(interval::NTuple{2,Number})
    (a,b) = float.(promote(interval...))
    Domain(@SVector[a], @SVector[b])
end

function Domain(dom::Domain)
    dom
end

function Domain(lims)
    lower = SVector(map(float ∘ first, lims))
    upper = SVector(map(float ∘ last, lims))
    Domain(lower, upper)
end

function uniform(rng, dom::Domain)
    V = pointtype(typeof(dom))
    Δ = (dom.upper - dom.lower)
    dom.lower .+ rand(rng, V) .* Δ
end

"""

    integral(f, dom [,alg])

Compute the integral of the function `f` over a domain `dom`, via the algorithm `alg`.
"""
function integral(f, dom, alg=Vegas())
    f2, dom2, alg2 = canonicalize(f, dom, alg)
    integral_kernel(f2, dom2, alg2)
end

function canonicalize(f, dom, alg)
    f, Domain(dom), alg
end

function canonicalize(f, I::NTuple{2,Number}, alg)
    f_v = f ∘ first
    dom = Domain(I)
    f_v, dom, alg
end

const ∫ = integral

function integral_kernel(f, dom::Domain, alg::MCVanilla)
    mc_kernel(f, alg.rng, dom, neval=alg.neval)
end
function draw(rng, dom::Domain)
    vol = volume(dom)
    x = uniform(rng, dom)
    (position=x, weight = volume(dom))
end

"""
    mc_kernel(f, rng::AbstractRNG, dom; neval)

Monte Carlo integration of function `f` over `dom`.
`dom` must support the following methods:
* volume(dom): Return the volume of dom
* draw(rng, dom): Return an object with properties `position::SVector`, `weight::Real`.
"""
function mc_kernel(f, rng::AbstractRNG, dom; neval)
    N = neval
    x = uniform(rng, dom)
    y = float(f(x)) * volume(dom)
    sum = y
    sum2 = y.^2
    for _ in 1:(N-1)
        s = draw(rng, dom)
        y = s.weight * f(s.position)
        sum += y
        sum2 += y.^2
    end

    mean = sum/N
    var_f = (sum2/N) - mean .^2
    var_f = max(zero(var_f), var_f)
    var_fmean = var_f / N
    std = sqrt.(var_fmean)
    (value = mean, std=std, neval=N)
end

function mc_kernel(f, p::ParallelRNG, dom; neval)
    rngs = p.rngs
    res1 = mc_kernel(f, rngs[1], dom, neval=2)
    T = typeof(res1)
    nthreads = length(rngs)
    results = Vector{T}(undef, nthreads)
    neval_i = ceil(Int, neval / nthreads)
    @threads for i in 1:nthreads
        res = mc_kernel(f, rngs[i], dom, neval=neval_i)
        results[i] = res
    end
    fuseall(results)
end

function fuseall(results)
    N = sum(res->res.neval, results)
    value = sum(res->res.value * res.neval/N   , results)
    var   = sum(res->(res.std * res.neval/N).^2, results)
    (value=value, std=sqrt.(var), neval=N)
end

# function fuse(res1, res2)
#     var1 = res1.std .^ 2
#     var2 = res2.std .^ 2
#     w1, w2 = fusion_weights(var1, var2)
#     (
#         value = @.(w1 * res1.value + w2*res2.value),
#         std = @.(sqrt(w1^2*var1 + w2^2*var2))
#     )
# end
#
# function fusion_weights(var1::AbstractVector, var2::AbstractVector)
#     pairs = fusion_weights_scalar.(var1, var2)
#     first.(pairs), last.(pairs)
# end
#
# function fusion_weights(var1, var2)
#     fusion_weights_scalar(var1, var2)
# end
#
# function fusion_weights_scalar(var1, var2)
#     z = zero(var1)
#     o = one(var1)
#     if iszero(var1)
#         (o, z)
#     elseif iszero(var2)
#         (z,o)
#     else
#         p1 = o/var1
#         p2 = o/var2
#         p  = p1 + p2
#         (p1/p, p2/p)
#     end
# end


text:
using H5SectionsArrays
using HDF5
using Base.Test

a = rand(UInt8, 200,200,3)

# create a fake dataset
tempDir = tempname()
mkdir( tempDir )

open(joinpath(tempDir, "registry.txt"), "w") do f_registry
    for z in 1:size(a,3)
        write(f_registry, "1,$(z)_aligned 0 100 -100 200 200 true \n")
        f = h5open(joinpath(tempDir, "1,$(z)_aligned.h5"), "w")
        f["img"] = a[:,:,z]
        f["offset"] = [100,-100]
        f["size"] = [200,200]
        close(f)
    end
end

@testset "test cutout" begin
    # cutout the chunk
    ba = H5SectionsArray(tempDir)
    b = ba[101:300, -99:100, 1:3]

    @test all(a.==b)

    # clean it up
    rm(tempDir; recursive=true)

end


text:
using DSP
mutable struct Measure
    bpm::Integer
    fs::Float64
    beats::Integer
    samples::AbstractArray
end
function format(m::Measure)
    return string("bpm: $(m.bpm)
          \nfs: $(m.fs)
          \nbeats: $(m.beats)
          \nsamples: $(length(m.samples))
                        ")
end
function addToMeasure!(m::Measure, signal::AbstractVector, beat)
    beatLength = (60/m.bpm) * m.fs
    measureLength = Int(ceil(beatLength * m.beats - beatLength*beat))
    actualMeasureLength = Int(ceil(beatLength * m.beats))
    if(length(signal) > measureLength)
        throw(DomainError("
        error: signal is too large
        signal length: $(length(signal))
        measureLength: $measureLength
        "))
    end
    signalToPush = Float64[]
    append!(signalToPush, zeros(Int(ceil(beatLength * beat))))
    append!(signalToPush, signal/rms(signal))
    append!(signalToPush, zeros(actualMeasureLength-length(signalToPush)))
    push!(m.samples, signalToPush)
end
function renderMeasure(m::Measure)
    length(m.samples) > 0 || throw(DomainError(length(m.samples),"samples must not be empty"))
    if(length(m.samples)==1)
        return m.samples[1]
    else
        return (sum(m.samples))./maximum(sum(m.samples))
    end
end
#notes added in sequential order
function addNotes(m::Measure, notes::AbstractVector, beat, instrument)
    maxNotes = Int(ceil(m.beats/beat))
    if(length(notes) > maxNotes)
        throw(DomainError("
        error: too many Notes!
        number of notes entered: $(length(notes))
        Max notes allowed: $maxNotes
        "))
    end

    for i in 1:length(notes)
        notes[i]===nothing||(
        addToMeasure!(m, instrument(notes[i]), (i-1)*beat),
                println(notes[i]))
    end

end

max_stars_repo_name: justinbroce/julia-audio-stuffs

max_issues_repo_name: justinbroce/julia-audio-stuffs

max_forks_repo_name: justinbroce/julia-audio-stuffs


text:
const max_methods = 4  # maximum number of MethodInstances tracked for a particular :call statement

"""
    framecode, lenv = get_call_framecode(fargs, parentframe::FrameCode, idx::Int)

Return the framecode and environment for a call specified by `fargs = [f, args...]` (see [`prepare_args`](@ref)).
`parentframecode` is the caller, and `idx` is the program-counter index.
If possible, `framecode` will be looked up from the local method tables of `parentframe`.
"""
function get_call_framecode(fargs::Vector{Any}, parentframe::FrameCode, idx::Int; enter_generated::Bool=false)
    nargs = length(fargs)  # includes f as the first "argument"
    # Determine whether we can look up the appropriate framecode in the local method table
    if isassigned(parentframe.methodtables, idx)  # if this is the first call, this may not yet be set
        # The case where `methodtables[idx]` is a `Compiled` has already been handled in `bypass_builtins`
        d_meth = d_meth1 = parentframe.methodtables[idx]::DispatchableMethod
        local d_methprev
        depth = 1
        while true
            # TODO: consider using world age bounds to handle cache invalidation
            # Determine whether the argument types match the signature
            sig = d_meth.sig.parameters::SimpleVector
            if length(sig) == nargs
                # If this is generated, match only if `enter_generated` also matches
                fi = d_meth.frameinstance
                if fi isa FrameInstance
                    matches = !is_generated(scopeof(fi.framecode)) || enter_generated == fi.enter_generated
                else
                    matches = !enter_generated
                end
                if matches
                    for i = 1:nargs
                        if !isa(fargs[i], sig[i])
                            matches = false
                            break
                        end
                    end
                end
                if matches
                    # Rearrange the list to place this method first
                    # (if we're in a loop, we'll likely match this one again on the next iteration)
                    if depth > 1
                        parentframe.methodtables[idx] = d_meth
                        d_methprev.next = d_meth.next
                        d_meth.next = d_meth1
                    end
                    if fi isa Compiled
                        return Compiled(), nothing
                    else
                        fi = fi::FrameInstance
                        return fi.framecode, fi.sparam_vals
                    end
                end
            end
            depth += 1
            d_methprev = d_meth
            d_meth = d_meth.next
            d_meth === nothing && break
            d_meth = d_meth::DispatchableMethod
        end
    end
    # We haven't yet encountered this argtype combination and need to look it up by dispatch
    fargs[1] = f = to_function(fargs[1])
    ret = prepare_call(f, fargs; enter_generated=enter_generated)
    ret === nothing && return f(fargs[2:end]...), nothing
    is_compiled = isa(ret[1], Compiled)
    local framecode
    if is_compiled
        d_meth = DispatchableMethod(nothing, Compiled(), ret[2])
    else
        framecode, args, env, argtypes = ret
        # Store the results of the method lookup in the local method table
        fi = FrameInstance(framecode, env, is_generated(scopeof(framecode)) && enter_generated)
        d_meth = DispatchableMethod(nothing, fi, argtypes)
    end
    if isassigned(parentframe.methodtables, idx)
        d_meth.next = parentframe.methodtables[idx]
        # Drop the oldest d_meth, if necessary
        d_methtmp = d_meth.next
        depth = 2
        while d_methtmp.next !== nothing
            depth += 1
            d_methtmp = d_methtmp.next
            depth >= max_methods && break
        end
        if depth >= max_methods
            d_methtmp.next = nothing
        end
    else
        d_meth.next = nothing
    end
    parentframe.methodtables[idx] = d_meth
    if is_compiled
        return Compiled(), nothing
    else
        return framecode, env
    end
end


text:
function quick_sort!(arr::Vector{T}, l::Int = 1, r::Int = length(arr)) where {T}
    if l >= r
        return
    end
    mid = arr[(l+r)>>1]
    l_pos = l
    r_pos = r
    while l_pos <= r_pos
        while arr[l_pos] < mid
            l_pos += 1
        end
        while arr[r_pos] > mid
            r_pos -= 1
        end
        if l_pos <= r_pos
            t = arr[l_pos]
            arr[l_pos] = arr[r_pos]
            arr[r_pos] = t
            l_pos += 1
            r_pos -= 1
        end
    end
    if l < r_pos
        quick_sort!(arr, l, r_pos)
    end
    if l_pos < r
        quick_sort!(arr, l_pos, r)
    end
end

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia