text:
using CUSPARSE
using CUDArt
using Base.Test

m = 25
n = 35
k = 10
blockdim = 5

#############
# test_sctr #
#############

function test_sctr!(elty)
    x = sparsevec(rand(1:m,k), rand(elty,k), m)
    y = zeros(elty,m)
    d_x = CudaSparseVector(x)
    d_y = CudaArray(y)
    d_y = CUSPARSE.sctr!(d_x,d_y,'O')
    h_y = to_host(d_y)
    if VERSION >= v"0.5.0-dev+742"
        y[x.nzind]  += x.nzval
    else
        y[x.rowval] += x.nzval
    end
    @test h_y ≈ y
end

function test_sctr(elty)
    x = sparsevec(rand(1:m,k), rand(elty,k), m)
    d_x = CudaSparseVector(x)
    d_y = CUSPARSE.sctr(d_x,'O')
    h_y = to_host(d_y)
    y = zeros(elty,m)
    if VERSION >= v"0.5.0-dev+742"
        y[x.nzind]  += x.nzval
    else
        y[x.rowval] += x.nzval
    end
    @test h_y ≈ y
end

types = [Float32,Float64,Complex64,Complex128]
for elty in types
    tic()
    test_sctr!(elty)
    println("sctr! took ", toq(), " for ", elty)
    tic()
    test_sctr(elty)
    println("sctr took ", toq(), " for ", elty)
end


text:
const DERIVED_NOUNS = Dict(
    Symbol("ACT PCPL") => ActiveParticle(),
    Symbol("PASS PCPL") => PassiveParticle(),
    Symbol("VN") => VerbalNoun()
)
const VERB_ASPECTS = Dict(
    Symbol("PERF") => Perfect(),
    Symbol("IMPF") => Imperfect(),
    Symbol("IMPV") => Imperative()
)
const VERB_MOODS = Dict(
    Symbol("IND") => Indicative(),
    Symbol("SUBJ") => Subjunctive(),
    Symbol("JUS") => Jussive()
)
const VERB_VOICES = Dict(
    Symbol("ACT") => Active(),
    Symbol("PASS") => Passive()
)
const VERB_FORMS = Dict(
    Symbol("(I)") => VerbFormI(),
    Symbol("(II)") => VerbFormII(),
    Symbol("(III)") => VerbFormIII(),
    Symbol("(IV)") => VerbFormIV(),
    Symbol("(V)") => VerbFormV(),
    Symbol("(VI)") => VerbFormVI(),
    Symbol("(VII)") => VerbFormVII(),
    Symbol("(VIII)") => VerbFormVIII(),
    Symbol("(IX)") => VerbFormIX(),
    Symbol("(X)") => VerbFormX(),
    Symbol("(XI)") => VerbFormXI(),
    Symbol("(XII)") => VerbFormXII()
)
const NOMINAL_CASES = Dict(
    Symbol("NOM") => Nominative(),
    Symbol("ACCC") => Accusative(Symbol("ACC"), "Accusative case", "منصوب"),
    Symbol("GEN") => Genetive()
)
const NOMINAL_STATES = Dict(
    Symbol("DEF") => Definite(),
    Symbol("INDEF") => Indefinite()
)
const NOUNS = Dict(
    Symbol("N") => Noun(),
    Symbol("PN") => ProperNoun()
)
const DERIVED_NOMINALS = Dict(
    Symbol("ADJ") => Adjective(),
    Symbol("IMPN") => ImperativeVerbalNoun()
)
const PRONOUNS = Dict(
    Symbol("PRON") => Personal(),
    Symbol("DEM") => Demonstrative(),
    Symbol("REL") => Relative()
)
const ADVERBS = Dict(
    Symbol("T") => Time(),
    Symbol("LOC") => Location()
)
const PREPOSITIONS = Dict(
    Symbol("P") => Preposition(),
    Symbol("Al+") => Preposition(Symbol("P"), "Determiner prefix ('the')", "حرف جر"),
    Symbol("bi+") => Preposition(Symbol("P"), "Preposition prefix ('by', 'with', 'in')", "حرف جر"),
    Symbol("ka+") => Preposition(Symbol("P"), "Preposition prefix ('like' or 'thus')", "حرف جر"),
    Symbol("ta+") => Preposition(Symbol("P"), "particle of oath prefix used as a preposition ('by Allah')", "حرف جر"),
    Symbol("w:P+") => Preposition(Symbol("P"), "Particle of oath prefix used as a preposition ('by the pen')", "حرف جر"),
    Symbol("l:P+") => Preposition(Symbol("P"), "Lam as a prefixed preposition", "حرف جر")
)
const LAM_PREFIXES = Dict(
    Symbol("l:EMPH+") => EmphaticLam(),
    Symbol("l:IMPV+") => ImperativeLam(),
    Symbol("l:PRP+") => PurposeLam()
)
const CONJUNCTIONS = Dict(
    Symbol("CONJ") => Coordinating(),
    Symbol("SUB") => Subordinating(),
    Symbol("w:CONJ+") => Coordinating(Symbol("CONJ"), "Waw as a conjunctin prefix ('and')", "حرف عطف"),
    Symbol("f:CONJ+") => Coordinating(Symbol("CONJ"), "Fa as a conjunction prefix ('and')", "حرف عطف")
)
const SUFFIXES = Dict(
    Symbol("+VOC") => Vocative(Symbol("VOC"), "Vocative suffix", "حرف نداء"),
    Symbol("+n:EMPH") => EmphaticNun(Symbol("+n:EMPH"), "Emphatic nun", "لام التوكيد")
)
const PARTICLES = Dict(
    Symbol("sa+") => Future(Symbol("FUT"), "Prefixed particle indicating the future ('they will')", "حرف استقبال"),
    Symbol("ya+") => Vocative(Symbol("VOC"), "A vocative prefix usually translated as 'O'", "حرف نداء"),
    Symbol("ha+") => Vocative(Symbol("VOC"), "A vocative prefix usually translated as 'Lo!'", "حرف نداء"),
    Symbol("A:INTG+") => Interogative(Symbol("INTG"), "Alif as a prefixed interrogative particle ('is?', 'did?', 'do?')", "حرف استفهام"),
    Symbol("A:EQ+") => Equalization(Symbol("EQ"), "Alif as a prefixed equalization particle ('whether')", "حرف تسوية"),
    Symbol("w:REM+") => Resumption(Symbol("REM"), "Waw as a resumption prefix ('then' or 'so')", "حرف استئنافية"),
    Symbol("w:CIRC+") => Circumstantial(Symbol("CIRC"), "Waw as a circumstantial prefix ('while')", "حرف حال"),
    Symbol("w:SUP+") => Supplemental(Symbol("SUP"), "Waw as a supplemental perfix ('then' pr 'so')", "حرف زائد"),
    Symbol("w:COM+") => Comitative(Symbol("COM"), "Waw as a comitative perfix ('with')", "واو المعية"),
    Symbol("f:REM+") => Resumption(Symbol("REM"), "Fa as a resumption prefix ('then' or 'so')", "حرف استئنافية"),
    Symbol("f:RSLT+") => Result(Symbol("RSLT"), "Fa as a result prefix ('then')", "حرف واقع في جواب الشرط"),
    Symbol("f:SUP+") => Supplemental(Symbol("SUP"), "Fa as a supplemental prefix ('then' or 'so')", "حرف زائد"),
    Symbol("f:CAUS+") => Cause(Symbol("CAUS"), "Fa as a cause prefix ('then' or 'so')", "حرف سببية"),
    Symbol("ACC") => Accusative(),
    Symbol("AMD") => Amendment(),
    Symbol("ANS") => Answer(),
    Symbol("AVR") => Aversion(),
    Symbol("CAUS") => Cause(),
    Symbol("CERT") => Certainty(),
    Symbol("CIRC") => Circumstantial(),
    Symbol("COM") => Comitative(),
    Symbol("COND") => Conditional(),
    Symbol("EQ") => Equalization(),
    Symbol("EXH") => Exhortation(),
    Symbol("EXL") => Explanation(),
    Symbol("EXP") => Exceptive(),
    Symbol("FUT") => Future(),
    Symbol("INC") => Inceptive(),
    Symbol("INT") => Interpretation(),
    Symbol("INTG") => Interogative(),
    Symbol("NEG") => Negative(),
    Symbol("PREV") => Preventive(),
    Symbol("PRO") => Prohibition(),
    Symbol("REM") => Resumption(),
    Symbol("RES") => Restriction(),
    Symbol("RET") => Retraction(),
    Symbol("RSLT") => Result(),
    Symbol("SUP") => Supplemental(),
    Symbol("SUR") => Surprise(),
    Symbol("VOC") => Vocative()
)
const DISCONNECTED = Dict(Symbol("INL") => DisconnectedLetters())
const PERSONS = Dict(
    Symbol("1") => FirstPerson(),
    Symbol("2") => SecondPerson(),
    Symbol("3") => ThirdPerson()
)
const GENDERS = Dict(
    Symbol("M") => Masculine(),
    Symbol("F") => Feminine()
)
const NUMBERS = Dict(
    Symbol("S") => Singular(),
    Symbol("D") => Dual(),
    Symbol("PL") => Plural()
)
const PARTOFSPEECH = Dict(
    VERB_VOICES..., VERB_ASPECTS..., VERB_MOODS..., VERB_FORMS..., NOMINAL_CASES...,
    NOMINAL_STATES..., NOUNS..., PRONOUNS..., DERIVED_NOMINALS..., CONJUNCTIONS..., ADVERBS...,
    DERIVED_NOUNS..., PREPOSITIONS..., DISCONNECTED..., LAM_PREFIXES..., PARTICLES...,
    PERSONS..., GENDERS..., NUMBERS..., SUFFIXES..., Symbol("V") => Verb()
)
const SIMPLE_ENCODING = Dict(
    Symbol(Char(0x0621)) => Symbol("Hamza"),
    Symbol(Char(0x0622)) => Symbol("Alif+Maddah"),
    Symbol(Char(0x0623)) => Symbol("HamzaAbove"),
    Symbol(Char(0x0624)) => Symbol("HamzaAbove"),
    Symbol(Char(0x0625)) => Symbol("HamzaBelow"),
    Symbol(Char(0x0626)) => Symbol("HamzaAbove"),
    Symbol(Char(0x0627)) => Symbol("Alif"),
    Symbol(Char(0x0628)) => Symbol("Ba"),
    Symbol(Char(0x0629)) => Symbol("TaMarbuta"),
    Symbol(Char(0x062A)) => Symbol("Ta"),
    Symbol(Char(0x062B)) => Symbol("Tha"),
    Symbol(Char(0x062C)) => Symbol("Jeem"),
    Symbol(Char(0x062D)) => Symbol("HHa"),
    Symbol(Char(0x062E)) => Symbol("Kha"),
    Symbol(Char(0x062F)) => Symbol("Dal"),
    Symbol(Char(0x0630)) => Symbol("Thal"),
    Symbol(Char(0x0631)) => Symbol("Ra"),
    Symbol(Char(0x0632)) => Symbol("Zain"),
    Symbol(Char(0x0633)) => Symbol("Seen"),
    Symbol(Char(0x0634)) => Symbol("Sheen"),
    Symbol(Char(0x0635)) => Symbol("Sad"),
    Symbol(Char(0x0636)) => Symbol("DDad"),
    Symbol(Char(0x0637)) => Symbol("TTa"),
    Symbol(Char(0x0638)) => Symbol("DTha"),
    Symbol(Char(0x0639)) => Symbol("Ain"),
    Symbol(Char(0x063A)) => Symbol("Ghain"),
    Symbol(Char(0x0640)) => Symbol("Tatweel"),
    Symbol(Char(0x0641)) => Symbol("Fa"),
    Symbol(Char(0x0642)) => Symbol("Qaf"),
    Symbol(Char(0x0643)) => Symbol("Kaf"),
    Symbol(Char(0x0644)) => Symbol("Lam"),
    Symbol(Char(0x0645)) => Symbol("Meem"),
    Symbol(Char(0x0646)) => Symbol("Noon"),
    Symbol(Char(0x0647)) => Symbol("Ha"),
    Symbol(Char(0x0648)) => Symbol("Waw"),
    Symbol(Char(0x0649)) => Symbol("AlifMaksura"),
    Symbol(Char(0x064A)) => Symbol("Ya"),
    Symbol(Char(0x064B)) => Symbol("Fathatan"),
    Symbol(Char(0x064C)) => Symbol("Dammatan"),
    Symbol(Char(0x064D)) => Symbol("Kasratan"),
    Symbol(Char(0x064E)) => Symbol("Fatha"),
    Symbol(Char(0x064F)) => Symbol("Damma"),
    Symbol(Char(0x0650)) => Symbol("Kasra"),
    Symbol(Char(0x0651)) => Symbol("Shadda"),
    Symbol(Char(0x0652)) => Symbol("Sukun"),
    Symbol(Char(0x0653)) => Symbol("Maddah"),
    Symbol(Char(0x0654)) => Symbol("HamzaAbove"),
    Symbol(Char(0x0670)) => Symbol("AlifKhanjareeya"),
    Symbol(Char(0x0671)) => Symbol("HamzatWasl"),
    Symbol(Char(0x06DC)) => Symbol("SmallHighSeen"),
    Symbol(Char(0x06DF)) => Symbol("SmallHighRoundedZero"),
    Symbol(Char(0x06E0)) => Symbol("SmallHighUprightRectangularZero"),
    Symbol(Char(0x06E2)) => Symbol("SmallHighMeemIsolatedForm"),
    Symbol(Char(0x06E3)) => Symbol("SmallLowSeen"),
    Symbol(Char(0x06E5)) => Symbol("SmallWaw"),
    Symbol(Char(0x06E6)) => Symbol("SmallYa"),
    Symbol(Char(0x06E8)) => Symbol("SmallHighNoon"),
    Symbol(Char(0x06EA)) => Symbol("EmptyCenterLowStop"),
    Symbol(Char(0x06EB)) => Symbol("EmptyCenterHighStop"),
    Symbol(Char(0x06EC)) => Symbol("RoundedHighStopWithFilledCenter"),
    Symbol(Char(0x06ED)) => Symbol("SmallLowMeem")
);
const BW_ENCODING = Dict(
    Symbol(Char(0x0621)) => Symbol('\''),
    Symbol(Char(0x0622)) => Symbol('('),
    Symbol(Char(0x0623)) => Symbol('>'),
    Symbol(Char(0x0624)) => Symbol('&'),
    Symbol(Char(0x0625)) => Symbol('<'),
    Symbol(Char(0x0626)) => Symbol('}'),
    Symbol(Char(0x0627)) => Symbol('A'),
    Symbol(Char(0x0628)) => Symbol('b'),
    Symbol(Char(0x0629)) => Symbol('p'),
    Symbol(Char(0x062A)) => Symbol('t'),
    Symbol(Char(0x062B)) => Symbol('v'),
    Symbol(Char(0x062C)) => Symbol('j'),
    Symbol(Char(0x062D)) => Symbol('H'),
    Symbol(Char(0x062E)) => Symbol('x'),
    Symbol(Char(0x062F)) => Symbol('d'),
    Symbol(Char(0x0630)) => Symbol('*'),
    Symbol(Char(0x0631)) => Symbol('r'),
    Symbol(Char(0x0632)) => Symbol('z'),
    Symbol(Char(0x0633)) => Symbol('s'),
    Symbol(Char(0x0634)) => Symbol('$'),
    Symbol(Char(0x0635)) => Symbol('S'),
    Symbol(Char(0x0636)) => Symbol('D'),
    Symbol(Char(0x0637)) => Symbol('T'),
    Symbol(Char(0x0638)) => Symbol('Z'),
    Symbol(Char(0x0639)) => Symbol('E'),
    Symbol(Char(0x063A)) => Symbol('g'),
    Symbol(Char(0x0640)) => Symbol('_'),
    Symbol(Char(0x0641)) => Symbol('f'),
    Symbol(Char(0x0642)) => Symbol('q'),
    Symbol(Char(0x0643)) => Symbol('k'),
    Symbol(Char(0x0644)) => Symbol('l'),
    Symbol(Char(0x0645)) => Symbol('m'),
    Symbol(Char(0x0646)) => Symbol('n'),
    Symbol(Char(0x0647)) => Symbol('h'),
    Symbol(Char(0x0648)) => Symbol('w'),
    Symbol(Char(0x0649)) => Symbol('Y'),
    Symbol(Char(0x064A)) => Symbol('y'),
    Symbol(Char(0x064B)) => Symbol('F'),
    Symbol(Char(0x064C)) => Symbol('N'),
    Symbol(Char(0x064D)) => Symbol('K'),
    Symbol(Char(0x064E)) => Symbol('a'),
    Symbol(Char(0x064F)) => Symbol('u'),
    Symbol(Char(0x0650)) => Symbol('i'),
    Symbol(Char(0x0651)) => Symbol('~'),
    Symbol(Char(0x0652)) => Symbol('o'),
    Symbol(Char(0x0653)) => Symbol('^'),
    Symbol(Char(0x0654)) => Symbol('#'),
    Symbol(Char(0x0670)) => Symbol('`'),
    Symbol(Char(0x0671)) => Symbol('{'),
    Symbol(Char(0x06DC)) => Symbol(':'),
    Symbol(Char(0x06DF)) => Symbol('@'),
    Symbol(Char(0x06E0)) => Symbol('\"'),
    Symbol(Char(0x06E2)) => Symbol('['),
    Symbol(Char(0x06E3)) => Symbol(';'),
    Symbol(Char(0x06E5)) => Symbol(','),
    Symbol(Char(0x06E6)) => Symbol('.'),
    Symbol(Char(0x06E8)) => Symbol('!'),
    Symbol(Char(0x06EA)) => Symbol('-'),
    Symbol(Char(0x06EB)) => Symbol('+'),
    Symbol(Char(0x06EC)) => Symbol('%'),
    Symbol(Char(0x06ED)) => Symbol(']')
);
const AR_DIACS_REGEX = Regex(
    string(
        Char(0x064B)[1], "|",
        Char(0x064C)[1], "|",
        Char(0x064D)[1], "|",
        Char(0x064E)[1], "|",
        Char(0x064F)[1], "|",
        Char(0x0640)[1], "|", # move this to special diac
        Char(0x0650)[1], "|",
        Char(0x0651)[1], "|",
        Char(0x0652)[1], "|",
        Char(0x0653)[1], "|",
        Char(0x0654)[1]
    )
);
const AR_DIACS = Symbol.(split(AR_DIACS_REGEX.pattern, "|"));
const SP_DEDIAC_MAPPING = Dict(
    Symbol(Char(0x0622)) => Symbol(Char(0x0627)),
    Symbol(Char(0x0623)) => Symbol(Char(0x0627)),
    Symbol(Char(0x0670)) => Symbol(Char(0x0627)),
    Symbol(Char(0x0671)) => Symbol(Char(0x0627)),
    Symbol(Char(0x0625)) => Symbol(Char(0x0627)),
    Symbol(Char(0x0624)) => Symbol(Char(0x0648)),
    Symbol(Char(0x0626)) => Symbol(Char(0x064A)),
    Symbol(Char(0x0649)) => Symbol(Char(0x064A)),
    Symbol(Char(0x0629)) => Symbol(Char(0x0647))
);
const SP_DEDIAC_KEYS = collect(keys(SP_DEDIAC_MAPPING));
const SP_REGEX_CHARS = ['\\', '^', '$', '.', '|', '?', '*', '+', ')', '(', ']', '[', '}', '{'];
const EN_CHAPTER_LABELS = ["The Opening", "The Cow", "The Family of 'Imran", "Women", "The Feast",
    "Livestock", "The Heights", "Battle Gains", "Repentance", "Jonah",
    "Hud", "Joseph", "Thunder", "Abraham", "Al-Hijr", "The Bee", "The Night Journey",
    "The Cave", "Mary", "Ta Ha", "The Prophets", "The Pilgrimage", "The Believers",
    "The Light", "The Differentiator", "The Poets", "The Ants", "The Story", "The Spider",
    "The Byzantines", "Luqman", "Bowing Down In Worship", "The Joint Forces",
    "Sheba", "The Creator", "Ya Sin", "Ranged In Rows", "Saad", "The Throngs",
    "The Foregiver", "[Verses] Made Clear", "Consultation", "Ornaments Of Gold",
    "Smoke", "Kneeling", "The Sand Dunes", "Muhammad", "Triumph", "The Private Rooms",
    "Qaaf", "Scattering [Winds]", "The Mountain", "The Star", "The Moon", "The Lord Of Mercy",
    "That Which Is Coming", "Iron", "The Dispute", "The Gathering [Of Forces]", "The Women Tested",
    "Solid Lines", "The Day of Congregation", "The Hypocrites", "Mutual Neglect",
    "Divorce", "Prohibition", "Control", "The Pen", "The Inevitable Hour", "The Ways Of Ascent",
    "Noah", "The Jinn", "Enfolded", "Wrapped In His Cloak", "The Resurrection", "Man",
    "[Winds] Sent Forth", "The Announcement", "The Forceful Chargers", "He Frowned",
    "The Rolling Up", "The Breaking", "Those Who Give Short Measure", "The Cracking",
    "The Towering Constellations", "The Night-Comer", "The Most High", "The Overwhelming Event",
    "Daybreak", "The City", "The Sun", "The Night", "The Morning Brightness", "Relief",
    "The Fig", "The Clinging Form", "The Night Of Glory", "Clear Evidence", "The Earthquake",
    "The Charging Steeds", "The Crashing Blow", "Competing For More", "The Fading Day",
    "The Backbiter", "The Elephant", "Quraysh", "Common Kindness", "Abundance", "The Disbelievers",
    "Help", "Palm Fibre", "Purity of Faith", "Daybreak", "People"];
const AR_CHAPTER_LABELS = ["ٱلْفَاتِحَة" ,"ٱلْبَقَرَة", "آلِ عِمْرَان", "ٱلنِّسَاء", "ٱلْمَائِدَة" , "ٱلْأَنْعَام", "ٱلْأَعْرَاف", "ٱلْأَنْفَال",
    "ٱلتَّوْبَة", "يُونُس", "هُود", "يُوسُف", "ٱلرَّعْد", "إِبْرَاهِيم", "ٱلْحِجْر", "ٱلنَّحْل",
    "ٱلْإِسْرَاء", "ٱلْكَهْف", "مَرْيَم", "طه", "ٱلْأَنْبِيَاء", "ٱلْحَجّ", "ٱلْمُؤْمِنُون", "ٱلنُّور",
    "ٱلْفُرْقَان", "ٱلشُّعَرَاء", "ٱلنَّمْل", "ٱلْقَصَص", "ٱلْعَنْكَبُوت", "ٱلرُّوم", "لُقْمَان",
    "ٱلسَّجْدَة", "ٱلْأَحْزَاب", "سَبَأ", "فَاطِر", "يس", "ٱلصَّافَّات", "ص", "ٱلزُّمَر", "غَافِر",
    "فُصِّلَت", "ٱلشُّورىٰ", "ٱلْزُّخْرُف", "ٱلدُّخَان", "ٱلْجَاثِيَة", "ٱلْأَحْقَاف", "مُحَمَّد", "ٱلْفَتْح",
    "ٱلْحُجُرَات", "ق", "ٱلذَّارِيَات", "ٱلطُّور", "ٱلنَّجْم", "ٱلْقَمَر", "ٱلرَّحْمَٰن", "ٱلْوَاقِعَة", "ٱلْحَدِيد",
    "ٱلْمُجَادِلَة", "ٱلْحَشْر", "ٱلْمُمْتَحَنَة", "ٱلصَّفّ", "ٱلْجُمُعَة", "ٱلْمُنَافِقُون", "ٱلتَّغَابُن", "ٱلطَّلَاق",
    "ٱلتَّحْرِيم", "ٱلْمُلْك", "ٱلْقَلَم", "ٱلْحَاقَّة", "ٱلْمَعَارِج", "نُوح", "ٱلْجِنّ", "ٱلْمُزَّمِّل",
    "ٱلْمُدَّثِّر", "ٱلْقِيَامَة", "ٱلْإِنْسَان", "ٱلْمُرْسَلَات", "ٱلنَّبَأ", "ٱلنَّازِعَات", "عَبَسَ", "ٱلتَّكْوِير",
    "ٱلْإِنْفِطَار", "ٱلْمُطَفِّفِين", "ٱلْإِنْشِقَاق", "ٱلْبُرُوج", "ٱلطَّارِق", "ٱلْأَعْلَىٰ", "ٱلْغَاشِيَة",
    "ٱلْفَجْ", "ٱلْبَلَد", "ٱلشَّمْس", "ٱللَّيْل", "ٱلضُّحَىٰ", "ٱلشَّرْح", "ٱلتِّين", "ٱلْعَلَق",
    "ٱلْقَدْر", "ٱلْبَيِّنَة", "ٱلزَّلْزَلَة", "ٱلْعَادِيَات", "ٱلْقَارِعَة", "ٱلتَّكَاثُر", "ٱلْعَصْر",
    "ٱلْهُمَزَة", "ٱلْفِيل", "قُرَيْش", "ٱلْمَاعُون", "ٱلْكَوْثَر", "ٱلْكَافِرُون", "ٱلنَّصْر", "ٱلْمَسَد",
    "ٱلْإِخْلَاص", "ٱلْفَلَق", "ٱلنَّاس"];

ext: jl


text:
#--------------------------------------------------------------------------------------------------
using Test

#--------------------------------------------------------------------------------------------------

@testset "Hensel" begin
    @testset "sundry tests" begin
        include("sunry_tests.jl")
    end
    @testset "mfill tests" begin
        include("mfill_tests.jl")
    end
    @testset "mexpansion tests" begin
        include("mexpansion_tests.jl")
    end
end

#--------------------------------------------------------------------------------------------------

max_stars_repo_name: ayarodionov/Mahler

max_issues_repo_name: ayarodionov/Mahler

max_forks_repo_name: ayarodionov/Mahler


text:
@testset "Results" begin
    rho = 0.8
    nrho = 50
    p = OrnsteinZernike.Parameters(1.0, 6.0, rho, 2^9, nrho, 1.0)
    pot = OrnsteinZernike.HardSphere()
    st = OrnsteinZernike.Structure(p, pot)
    cls = OrnsteinZernike.PercusYevick()

    inter = OrnsteinZernike.Interaction(p, st, cls)

    result = OrnsteinZernike.solve(inter)
    if 1.0 ∈ result.r
        first_contact = result.gr[result.r .== 1.0][1]
    else
        first_contact = maximum(result.gr)
    end
    ground_truth = 3.47 # For Percus-Yevick, this is the best it can get
    @show first_contact
    @test isapprox(first_contact, ground_truth; atol=0.01)
end


text:
function is_primal_infeasible!(δy, ws)

	settings = ws.settings
	# calculate unscaled norm of δy
	norm_δy = scaled_norm(ws.sm.E, δy.data, Inf)::eltype(δy)

	# make sure norm is unequal to zero before continuing
	if norm_δy > settings.eps_prim_inf
		# test condition A'δy = 0
		#A_δy = ws.p.A' * δy
		A_δy  = ws.utility_vars.vec_n
		mul!(A_δy, ws.p.A', δy.data)
		# unscale: A_δy = Dinv * A_δy
		lmul!(ws.sm.Dinv, A_δy)

		if norm(A_δy, Inf) <= settings.eps_prim_inf * norm_δy
			# scale δy to unit size
			@. δy *= (-1 / norm_δy)
			δyt_b = dot(δy, ws.p.b)
			# notice the in-place function. This is faster as we are using minus_unit_δy as workspace
			sF = support_function!(δy, ws.p.C, settings.eps_prim_inf) - δyt_b
			if sF <= settings.eps_prim_inf
				return true
			end
		end
	end
	return false
end

function is_dual_infeasible!(δx, ws)
	settings = ws.settings
	# calculate unscaled norm of δx
	norm_δx = scaled_norm(ws.sm.D, δx, Inf)::eltype(δx)

	if norm_δx > settings.eps_dual_inf
		# test condition <q,δx> < 0
		if  dot(ws.p.q, δx) / (norm_δx * ws.sm.c[]) < -settings.eps_dual_inf

			# test condition Pδx == 0
			P_δx = ws.utility_vars.vec_n
			# P_δx = P * δx
			mul!(P_δx, ws.p.P, δx)

			# unscale P_δx = Dinv * P_δx
			lmul!(ws.sm.Dinv, P_δx)

			if norm(P_δx, Inf) / (norm_δx * ws.sm.c[]) <= settings.eps_dual_inf

				# test condition Ax in Ktilde_b∞
				A_δx = ws.utility_vars.vec_m
				mul!(A_δx, ws.p.A, δx)

				# unscale A_δx = Einv * P_δx
				lmul!(ws.sm.Einv, A_δx)

				#normalize A_δx and create a SplitVector
				A_δx .*= 1/norm_δx
				A_δx_split = SplitVector(A_δx, ws.p.C)
				if in_pol_recc!(A_δx_split, ws.p.C, settings.eps_dual_inf)
					return true
				end
			end
		end
	end
	return false
end


text:
"""
Return the sum of ``val1`` and ``val2`` as two float64s, an integer part
and the fractional remainder.  If ``factor`` is not 1.0 then multiply the
sum by ``factor``.  If ``divisor`` is not 1.0 then divide the sum by
``divisor``.
The arithmetic is all done with exact floating point operations so no
precision is lost to rounding error.  This routine assumes the sum is less
than about 1e16, otherwise the ``frac`` part will be greater than 1.0.
Returns
-------
day, frac : float64
    Integer and fractional part of val1 + val2.
"""
function day_frac(val1, val2; factor=1.0, divisor=1.0)
    # Add val1 and val2 exactly, returning the result as two float64s.
    # The first is the approximate sum (with some floating point error)
    # and the second is the error of the float64 sum.
    sum12, err12 = two_sum(val1, val2)

    if factor != 1.0
        sum12, carry = two_product(sum12, factor)
        carry += err12 * factor
        sum12, err12 = two_sum(sum12, carry)
    end

    if divisor != 1.0
        q1 = sum12 / divisor
        p1, p2 = two_product(q1, divisor)
        d1, d2 = two_sum(sum12, -p1)
        d2 += err12
        d2 -= p2
        q2 = (d1 + d2) / divisor  # 3-part float fine here; nothing can be lost
        sum12, err12 = two_sum(q1, q2)
    end

    # get integer fraction
    day = Int(round(sum12))
    extra, frac = two_sum(sum12, -day)
    frac += extra + err12
    day, frac
end

"""
Add ``a`` and ``b`` exactly, returning the result as two float64s.
The first is the approximate sum (with some floating point error)
and the second is the error of the float64 sum.
Using the procedure of Shewchuk, 1997,
Discrete & Computational Geometry 18(3):305-363
http://www.cs.berkeley.edu/~jrs/papers/robustr.pdf
Returns
-------
sum, err : float64
    Approximate sum of a + b and the exact floating point error
"""
function two_sum(a, b)
    x = a + b
    eb = x - a
    eb = b - eb
    ea = x - b
    ea = a - ea
    x, ea + eb
end

"""
Multiple ``a`` and ``b`` exactly, returning the result as two float64s.
The first is the approximate product (with some floating point error)
and the second is the error of the float64 product.
Uses the procedure of Shewchuk, 1997,
Discrete & Computational Geometry 18(3):305-363
http://www.cs.berkeley.edu/~jrs/papers/robustr.pdf
Returns
-------
prod, err : float64
    Approximate product a * b and the exact floating point error
"""
function two_product(a, b)
    x = a * b
    ah, al = split(a)
    bh, bl = split(b)
    y1 = ah * bh
    y = x - y1
    y2 = al * bh
    y -= y2
    y3 = ah * bl
    y -= y3
    y4 = al * bl
    y = y4 - y
    x, y
end

"""
Split float64 in two aligned parts.
Uses the procedure of Shewchuk, 1997,
Discrete & Computational Geometry 18(3):305-363
http://www.cs.berkeley.edu/~jrs/papers/robustr.pdf
"""
function split(a)
    c = 134217729.0 * a  # 2**27+1.
    abig = c - a
    ah = c - abig
    al = a - ah
    ah, al
end


text:
# parsing utils
macro anyeq(var, vals...)
    ret = e = Expr(:||)
    for (i, v) in enumerate(vals)
        x = :($var == $v)
        push!(e.args, x)
        i >= length(vals) - 1 && continue
        ne = Expr(:||)
        push!(e.args, ne)
        e = ne
    end
    return esc(ret)
end

@inline lower(c) = Char(UInt32(c) | 0x20)
@inline ismark(c) = @anyeq(c, '-', '_', '.', '!', '~', '*', '\'', '(', ')')
@inline isalpha(c) = 'a' <= lower(c) <= 'z'
@inline isnum(c) = '0' <= c <= '9'
@inline isalphanum(c) = isalpha(c) || isnum(c)
@inline isuserinfochar(c) = isalphanum(c) || ismark(c) || @anyeq(c, '%', ';', ':', '&', '=', '+', '$', ',')
@inline ishex(c) =  isnum(c) || ('a' <= lower(c) <= 'f')
@inline ishostchar(c) = isalphanum(c) || @anyeq(c, '.', '-', '_', '~')
@inline isheaderchar(c) = c == CR || c == LF || c == Char(9) || (c > Char(31) && c != Char(127))

"""
    escapelines(string)

Escape `string` and insert '\n' after escaped newline characters.
"""

function escapelines(s)
    s = Base.escape_string(s)
    s = replace(s, "\\n", "\\n\n    ")
    return string("    ", strip(s))
end


text:
function Euler9()

	for a in (1:1000)
		for b in (a:1000)

			cc = a^2+b^2
			if mod(cc,sqrt(cc))==0
				if a + b + sqrt(cc) == 1000
					println(int(a*b*sqrt(cc)))
					break
				end
			end
		end
	end

end

@time Euler9()

max_stars_repo_name: gribeill/ProjectEuler

max_issues_repo_name: gribeill/ProjectEuler

max_forks_repo_name: gribeill/ProjectEuler


text:
try
    close(yingshi)
end
yingshi = UDPSocket();
bind(yingshi,ip"73.108.52.5",13000); # change this ip to the server (or julia code)

gains=[30.0,0.15,3.0];
function tune(gains)
  MsgOutString = ' ';
  for j in 1:length(gains)
      MsgOutString = string(MsgOutString,' ',gains[j]);
  end
  MsgOutString = string(MsgOutString," \n");
  send(yingshi,ip"141.212.141.245",36896,MsgOutString);
  # change this to the ip where you are running Simulink!
  print(MsgOutString)
end

#tune([90.0,20.,0.008]) # good gains
tune([90.0,0.0,0.0]) # good gains


text:
using Pkg

if !in("BayesNets", keys(Pkg.installed()))
	Pkg.add("BayesNets")
end

if !in("LightGraphs", keys(Pkg.installed()))
	Pkg.add("LightGraphs")
end

if !in("DataFrames", keys(Pkg.installed()))
	Pkg.add("DataFrames")
end

if !in("CSV", keys(Pkg.installed()))
	Pkg.add("CSV")
end

if !in("JLD2", keys(Pkg.installed()))
	Pkg.add("JLD2")
end

max_issues_repo_name: sisl/em-model-hobbydrone


text:
module Lu
#
  importall Shared_types
  using Shared_funcs
  using Conv
#
  :Dense
  :Sparse
#
  function MatrixTypesT(T) return Union(Matrix{T}, SparseMatrixCSC{T}) end
  MatrixTypes = Union(Matrix, SparseMatrixCSC)
#
  type Problem{T}
    conv::Conv.Converter
    params::Params
    r::Int
    c::Int
    A::MatrixTypesT(T)
  #
    Problem() = new()
  end
#
type Solution{T}
  solved::Bool
  L::MatrixTypesT(T)
  U::MatrixTypesT(T)
  p::Vector{Int}
  q::Vector{Int}
  rs::Vector{T}
  #
  function Solution(params::Params) x = new(); x.solved = false; return x; end
end
#
  function is_dense(prob::Problem) return issparse(prob.A) == false end
  function is_dense(sol::Solution) return issparse(sol.L) == false end
  function get_form(prob::Problem) return (is_dense(prob) ? :Dense : :Sparse) end
  function get_t(prob::Problem) return prob.conv.t end
  function get_r(prob::Problem) return prob.r end
  function get_c(prob::Problem) return prob.c end
  function comp_density(prob::Problem) return Shared_funcs.comp_density(prob.A) end
  function construct_Problem(params::Params, A::MatrixTypes)
    conv = Conv.converter(params["type"])
    prob = Problem{conv.t}()
    prob.conv = conv; prob.params = params;
    return prob
  end
  function fill_problem{T}(params::Params, prob::Problem{T}, A::MatrixTypes)
    prob.r, prob.c = size(A)
    prob.A = Conv.matrix(prob.conv, A)
  end
  function create_Problem(params::Params, A::MatrixTypes)
    prob = construct_Problem(params, A)
    fill_problem(params, prob, A)
    return prob
  end
  function construct_solution(type_t::DataType, params::Params)
    return Solution{type_t}(params)
  end
  function calc_solution_distance(prob::Problem, sol::Solution)
    if (sol.solved == false) return Inf end
    #println(dense(sol.L)); println(dense(sol.U));
    if (all(isfinite(sol.L)) == false) return Inf end
    if (all(isfinite(sol.U)) == false) return Inf end
    if (istril(sol.L) == false) return Inf end
    if (istriu(sol.U) == false) return Inf end
    nA = prob.A
    nA = isdefined(sol, :rs) ? scale(sol.rs, nA) : nA
    nA = nA[isdefined(sol, :p) ? sol.p : range(1,end),isdefined(sol, :q) ? sol.q : range(1,end)]
    normp = 1 #2-norm not yet implemented for sparse matrices
    #println(sol.p); println(sol.q); println(sol.rs);
    #println(dense(nA)); println(dense(sol.L*sol.U));
    return norm(nA - sol.L*sol.U, normp)
  end
end

max_stars_repo_name: jadnohra/mlcp-lib

max_issues_repo_name: jadnohra/mlcp-lib

max_forks_repo_name: jadnohra/mlcp-lib


text:
# Constraint Manager holds a list of constaints and can run the necessary
# calculus

include("constraints.jl")

using SparseArrays

abstract type constraintManager end

"""
    constraintManager_Base

Contains a list of constraints and the dual variables (lambda) that correspond
with each constraint.

See also [`constraintManager_Dynamics`](@ref)
"""
struct constraintManager_Base <: constraintManager
    cList::Array{constraint, 1}
    lambdaList
end

#############################################################
###      Update the dual variable of the constraints      ###
#############################################################
@doc raw"""
    updateDual(cM::constraintManager, yCurr, penalty::Float64)

Updates (in place) each λ in the lambdaList stores in the
[`constraintManager`](@ref).

```math
λ ⟵ λ + ρ c(y)
```

for each constraint

"""
function updateDual!(cM::constraintManager, yCurr, penalty::Float64)
    for (i, c) in enumerate(cM.cList)
        # Obtain the dual variable matching this constraint
        lambda = cM.lambdaList[i]

        # Evaluate the current constraint
        if (typeof(c) == AL_simpleCone) | (typeof(c) == AL_Multiple_simpleCone)
            cVal = getNormToProjVals(c, yCurr, lambda)
        else
            cVal = getNormToProjVals(c, yCurr)
        end
        # cVal = getNormToProjVals(c, yCurr)

        lambda = lambda + penalty * cVal
    end
end

#############################################################
###            Evaluate a List of Constraints             ###
#############################################################

function isCone(c::constraint)
    if (typeof(c) == AL_simpleCone) |
        (typeof(c) == AL_Multiple_simpleCone) |
        (typeof(c) == AL_simpleAngleCone) |
        (typeof(c) == AL_Multiple_AngleCone)

        return true
    end
    return false
end

@doc raw"""
    evalConstraints(yCurr, cM::constraintManager, penalty::Float64)

Evaluates each constraint and sums the result.

Recall that the constraint terms are given as
```math
ϕ(x) = f(x) + ρ Σ ||cᵢ(x)||² + Σ λᵢ^{\top}cᵢ(x)
```

This function returns the following part:
```math
ρ Σ ||cᵢ(x)||² + Σ λᵢ^{\top}cᵢ(x)
```

See also: [`constraintManager`](@ref)
"""
function evalConstraints(cM::constraintManager, yCurr, penalty::Float64)

    total_eval = 0

    # println("Evaluating the constraints")

    for (i, c) in enumerate(cM.cList)
        # Obtain the dual variable matching this constraint
        lambda = cM.lambdaList[i]
        # Evaluate the current constraint
        if isCone(c)
            cVal = getNormToProjVals(c, yCurr, lambda)
        else
            cVal = getNormToProjVals(c, yCurr)
        end
        # cVal = getNormToProjVals(c, yCurr)
        # println("Lambda = $(size(lambda)) and cVal = $(size(cVal))")
        # println("So, ||c(x)|| = $(norm(cVal)) -> $(size(norm(cVal)))")
        # println("And λ'c = $(lambda' * cVal) -> $(size(lambda' * cVal))")
        # println()

        # Evaluate and add the full constraint term
        total_eval += penalty * (norm(cVal)^2) + (lambda' * cVal)[1]
    end

    return total_eval
end

"""
    getConstraintViolationList(cM::constraintManager, yList,
                               penalty::Float64 = 1.0)

Takes a list of trajectories and calculates the TOTAL constraint violation
at each trajectory.
"""
function getConstraintViolationList(cM::constraintManager, yList,
                                    penalty::Float64 = 1.0)
    constVio = zeros(size(yList, 1))

    for (ind, y) in enumerate(yList)
        constVio[ind] = evalConstraints(cM, y, penalty)
    end

    return constVio
end

#############################################################
###     Evaluate the Gradient of a List of Constraints    ###
#############################################################
@doc raw"""
    evalConstraints(yCurr, cM::constraintManager, penalty::Float64)

Evaluates the gradient of each constraint and sums the result.

Recall that the constraint terms are given as
```math
∇ϕ(x) = ∇f(x) + Σ J(c)'(ρ * cᵢ(x) + λᵢ)
```

This function returns the following part:
```math
Σ J(c)'(ρ * cᵢ(x) + λᵢ)
```

See also: [`constraintManager`](@ref)
"""
function evalGradConstraints(cM::constraintManager, yCurr, penalty::Float64)

    if size(yCurr, 1) == 1
        # If this is a one dimensional optimization problem, the gradient
        # is actually just a derivative. Hence the vector notation will cause
        # an error.
        total_grad = 0
    else
        # otherwise the vector notation is necessary
        total_grad = zeros(size(yCurr))
    end

    for (i, c) in enumerate(cM.cList)
        # Select the correct dual variables
        lambda = cM.lambdaList[i]
        # Evaluate the constraint
        if (typeof(c) == AL_simpleCone) | (typeof(c) == AL_Multiple_simpleCone)
            cVal = getNormToProjVals(c, yCurr, lambda)
        else
            cVal = getNormToProjVals(c, yCurr)
        end
        # cVal = getNormToProjVals(c, yCurr)
        # Evaluate the jacobian of the constraint
        cJacob = getGradC(c, yCurr)
        # println("cJacob:")
        # display(cJacob)
        # println("cVal:")
        # println(cVal)

        # Evaluate and add the full constraint term
        total_grad += cJacob' * (penalty * cVal + lambda)
    end

    return total_grad

end

#############################################################
###     Evaluate the Hessian of a List of Constraints     ###
#############################################################
@doc raw"""
    evalHessConstraints(yCurr, cM::constraintManager, penalty::Float64)

Evaluates the Hessian of each constraint and sums the result.

Recall that the constraint terms are given as
```math
H(ϕ(x)) = H(f(x)) + H(ρ Σ ||cᵢ(x)||² + Σ λᵢ^{\top}cᵢ(x))
```

This function returns the following part:
```math
H(ρ Σ ||cᵢ(x)||² + Σ λᵢ^{\top}cᵢ(x)) = Σ H(ρ ||cᵢ(x)||² + λᵢ^{\top}cᵢ(x))
```

See also: [`constraintManager`](@ref)
"""
function evalHessConstraints(cM::constraintManager, yCurr, rho = 1)

    total_hess = spzeros(size(yCurr, 1), size(yCurr, 1))

    for c in cM.cList
        total_hess += Symmetric(getHessC_ALTerm(c, yCurr, rho))
        #getHessC(c)#getHessC_ALTerm(c, yCurr, rho)
    end

    return total_hess

end

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
import PackageCompiler
PackageCompiler.create_sysimage(["Plots", "Markdown", "DataFrames"]; sysimage_path = "DSSysimage.so", precompile_execution_file="precompilation.jl")

max_stars_repo_name: Jblew/towards-julia-notebook

max_issues_repo_name: Jblew/towards-julia-notebook

max_forks_repo_name: Jblew/towards-julia-notebook


text:
# truncation.jl
#
# Implements truncation schemes for truncating a tensor with svd, leftorth or rightorth

abstract TruncationScheme

immutable NoTruncation <: TruncationScheme
end
notrunc()=NoTruncation()

immutable TruncationError <: TruncationScheme
    epsilon::Real
end
truncerr(epsilon::Real)=TruncationError(epsilon)
Base.eps(t::TruncationError)=t.epsilon

immutable TruncationDimension <: TruncationScheme
    D::Integer
    epsilon::Real
end
truncdim(D::Int)=TruncationDimension(D,0)
dim(t::TruncationDimension)=t.D
Base.eps(t::TruncationDimension)=t.epsilon

Base.&(trunc1::TruncationDimension,trunc2::TruncationError)=TruncationDimension(dim(trunc1),max(eps(trunc1),eps(trunc2)))
Base.&(trunc2::TruncationError,trunc1::TruncationDimension)=TruncationDimension(dim(trunc1),max(eps(trunc1),eps(trunc2)))

immutable TruncationSpace{S<:ElementarySpace} <: TruncationScheme
    space::S
    epsilon::Real
end
truncspace(space::ElementarySpace) = TruncationSpace(space,0)
space(s::TruncationSpace) = s.space
dim(s::TruncationSpace) = dim(s.space)
Base.eps(t::TruncationSpace)=t.epsilon

Base.&(trunc1::TruncationSpace,trunc2::TruncationError)=TruncationSpace(space(trunc1),max(eps(trunc1),eps(trunc2)))
Base.&(trunc2::TruncationError,trunc1::TruncationSpace)=TruncationSpace(space(trunc1),max(eps(trunc1),eps(trunc2)))


text:
function MOI.copy!(dest::LinQuadOptimizer, src::MOI.ModelLike; copynames=false)
    return MOIU.defaultcopy!(dest, src, copynames)
end

function MOI.copy!(dest::MOI.ModelLike,    src::LinQuadOptimizer; copynames=false)
    return MOIU.defaultcopy!(dest, src, copynames)
end

function MOI.copy!(dest::LinQuadOptimizer, src::LinQuadOptimizer; copynames=false)
    return MOIU.defaultcopy!(dest, src, copynames)
end

MOI.canget(::LinQuadOptimizer, ::MOI.ListOfVariableAttributesSet) = true
function MOI.get(::LinQuadOptimizer, ::MOI.ListOfVariableAttributesSet)
    return MOI.AbstractVariableAttribute[]
end

MOI.canget(::LinQuadOptimizer, ::MOI.ListOfModelAttributesSet) = true
function MOI.get(model::LinQuadOptimizer, ::MOI.ListOfModelAttributesSet)
    return [MOI.ObjectiveSense(),
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        # TODO(odow): why is `nothing` returned?
        (MOI.get(model, MOI.Name()) != "") ? MOI.Name() : nothing
    ]
end

function MOI.canget(::LinQuadOptimizer, ::MOI.ListOfConstraintAttributesSet{F,S}) where {F, S}
    return true
end
function MOI.get(::LinQuadOptimizer, ::MOI.ListOfConstraintAttributesSet{F,S}) where {F, S}
    return MOI.AbstractConstraintAttribute[]
end


text:
module PPG2KP

# Constant used for the problems that need no change in the enumeration
# process (and often share the same method bodies).
const SIMILAR_4 = Union{Val{:G2KP}, Val{:G2MKP}, Val{:G2OPP}, Val{:G2CSP}}

# Include submodules. Even if not used here, they need to be here to be
# avaliable for users to access/import.
include("./Heuristic.jl")
include("./Args.jl")
include("./Enumeration.jl")

# Everything that is exported by PPG2KP.Enumeration ends up being used here.
using .Enumeration # Where all the plate enumeration logic is defined.
export ByproductPPG2KP # re-export ByproductPPG2KP from Enumeration

# Import third-party libraries.
using JuMP
using UnPack # for @unpack
import TimerOutputs.@timeit
import MathOptInterface
using RandomNumbers.Xorshifts # for Xoroshiro128Plus (for heuristic rng)
const MOI = MathOptInterface
using DocStringExtensions # For better struct documentation

# Import globals/modules/functions defined in GuillotineModels (main module).
import ..TIMER # Global module timer for use with TimerOutputs.@timeit.
import ..throw_if_timeout_now
import ..to_pretty_str # for debug
import ..CutPattern # type returned by get_cut_pattern
import ..BuildStopReason, ..BUILT_MODEL, ..FOUND_OPTIMUM
using ..Utilities
import ..Data.G2KP
import ..Data.SLOPP
import ..Data.MHLOPPW
import ..Data.SSSCSP
using ..Utilities: expand # for get_cut_pattern

# Included at the start because ModelByproduct has a RotationAwareData field.
include("rotation.jl")

struct ModelByproduct{D, S, P}
	preprocess_byproduct :: ByproductPPG2KP{D, S, P}
	found_optimum :: Bool
	# TODO: `optimum_if_found` will need to become a `Vector` of
	# `CutPattern`s in the future if this feature is to be supported by
	# G2CSP and G2MKP.
	optimum_if_found :: CutPattern{D, S}
	rad :: Union{Nothing, RotationAwareData{D, S}}
end

# Convenience in case of FOUND_OPTIMUM. `found_optimum` is set
# to `true` and `optimum_if_found` stores the solution.
function ModelByproduct(
	bp :: ByproductPPG2KP{D, S, P}, optimum :: CutPattern{D, S}
) where {D, S, P}
	return ModelByproduct{D, S, P}(bp, true, optimum)
end

# Convenience in case of BUILT_MODEL. `found_optimum` is set
# to `false` and `optimum_if_found` stores an empty (but valid) solution.
function ModelByproduct(bp :: ByproductPPG2KP{D, S, P}) where {D, S, P}
	return ModelByproduct{D, S, P}(bp, false, CutPattern(bp.L, bp.W), nothing)
end
# Convenience in case of BUILT_MODEL and rotation was allowed.
function ModelByproduct(
	bp :: ByproductPPG2KP{D, S, P}, rad :: RotationAwareData{D, S}
) where {D, S, P}
	return ModelByproduct{D, S, P}(bp, false, CutPattern(bp.L, bp.W), rad)
end

@enum Dimension begin
	LENGTH = 1
	WIDTH = 2
end
Base.getindex(tuple :: Tuple, s :: Dimension) = tuple[Int(s)]

@enum Relation begin
	PARENT = 1
	FIRST_CHILD = 2
	SECOND_CHILD = 3
end
Base.getindex(tuple :: Tuple, s :: Relation) = tuple[Int(s)]

@enum BaseModel BECKER FURINI

# Include code that is focused on implementing a specific feature, but does
# not merit a module of their own. Non-exported methods of these files are not
# expected to be used in this module, but exported methods are exports from
# this module (ou specialization of parent modules).
include("./get_cut_pattern.jl")
include("./MIP_start.jl")

function _partition_by_bits(bits, list)
	@assert length(bits) == length(list) # yes, it is length, not axes, what we want
	accepted = similar(list)
	rejected = similar(list)
	indexes = LinearIndices(list)
	qt_accepted = qt_rejected = 0 # yes, this does not need to be a specific type
	for (istrue, value) in zip(bits, list)
		if istrue
			qt_accepted += 1
			accepted[indexes[qt_accepted]] = value
		else
			qt_rejected += 1
			rejected[indexes[qt_rejected]] = value
		end
	end
	resize!(accepted, qt_accepted)
	resize!(rejected, qt_rejected)
	return accepted, rejected
end

# Delete the vars from the model and from the byproduct.
# Returns byproduct because if it is :cuts_made, then the old object is
# invalid and a new one is necessary (the vectors are changed inplace, but
# the first_vertical_cut_idx cannot be updated because the struct is
# immutable). If it is :picuts, then the passed struct is updated
# to a valid state and then returned.
function _delete_vars!(bp, model, to_keep_bits, var_set_name_in_model)
	@assert var_set_name_in_model in (:picuts, :cuts_made)
	var_set_name_in_bp = var_set_name_in_model == :picuts ? :np : :cuts
	vars = model[var_set_name_in_model]
	to_keep_idxs, to_delete_idxs = _partition_by_bits(to_keep_bits, keys(vars))
	deleteat!(getfield(bp, var_set_name_in_bp), to_delete_idxs)
	to_keep_vars, to_delete_vars = _partition_by_bits(to_keep_bits, vars)
	JuMP.delete(model, to_delete_vars)
	model[var_set_name_in_model] = to_keep_vars
	# Update the variable names to make the debug to make more sense.
	for (i, var) in enumerate(to_keep_vars)
		set_name(var, "$(var_set_name_in_model)[$(i)]")
	end
	var_set_name_in_model == :picuts && return bp
	new_fvci = searchsortedfirst(to_keep_idxs, bp.first_vertical_cut_idx)
	return ByproductPPG2KP( # new_fvci was computed and need to be updated
		bp.cuts, new_fvci, bp.np, bp.pli_lwb, bp.d, bp.l, bp.w, bp.L, bp.W,
		bp.mirror_plates
	)
end

# Include all code related to finding and removing unreachable plates/vars.
# It does make use of _partition_by_bits and _delete_vars! but both are
# used elsewhere so they are kept in this file. From the methods defined
# in the include below only _purge_unreachable! is used here.
include("./unreachable.jl")

"""
Auxiliar struct used to build the PPG2KP model.

Aggregates the inverted indexes that indicate the variables with non-zero
coefficients in model constraints.

$TYPEDFIELDS

"""
struct VarInvIndexes{P}
# Note the field order, as all of them have the same type but different
# meanings, switching the order in the default constructor can be calamitous.
"""
If indexed by a piece index, returns the indexes of all extraction variables
(picuts/np) that extract that piece from some plate.
"""
pii2pair :: Vector{Vector{P}}
"""
If indexed by a plate index, returns the indexes of all extraction variables
(picuts/np) that extract some piece from that plate.
"""
pli2pair :: Vector{Vector{P}}
"""
If indexed by a plate index, returns the indexes of all cut variables
(cuts_made/cuts) that have that plate as either first or second child (i.e.,
one of the plates generated by the cut).
"""
child2cut :: Vector{Vector{P}}
"""
If indexed by a plate index, returns the indexes of all cut variables
(cuts_made/cuts) that have that plate as the parent of the cut (i.e., the plate
that is being cut).
"""
parent2cut :: Vector{Vector{P}}
end

"""
    VarInvIndexes(num_piece_types, num_plate_types, np, cuts)

Utility constructor that creates the struct from a ByproductPPG2KP.
"""
@timeit TIMER function VarInvIndexes(
	num_piece_types :: P, num_plate_types :: P,
	@nospecialize(np), @nospecialize(cuts)
) :: VarInvIndexes{P} where {P <: Integer}
	pii2pair = [Vector{P}() for _ = 1:num_piece_types]
	pli2pair = [Vector{P}() for _ = 1:num_plate_types]
	child2cut = [Vector{P}() for _ = 1:num_plate_types]
	parent2cut = [Vector{P}() for _ = 1:num_plate_types]

	# Initialize all inverse indexes.
	for i in eachindex(cuts)
		parent, fchild, schild = cuts[i]
		push!(parent2cut[parent], i)
		!iszero(fchild) && push!(child2cut[fchild], i)
		!iszero(schild) && push!(child2cut[schild], i)
	end

	for i in eachindex(np)
		pli, pii = np[i]
		push!(pli2pair[pli], i)
		push!(pii2pair[pii], i)
	end

	# Note the parameter order, as all of them have the same type but different
	# meanings, switching the order can be calamitous.
	return VarInvIndexes(pii2pair, pli2pair, child2cut, parent2cut)
end

"""
    VarInvIndexes(bp :: ByproductPPG2KP{D, S, P}) :: VarInvIndexes{P}

Utility constructor that creates the struct from a ByproductPPG2KP.
"""
function VarInvIndexes(
	bp :: ByproductPPG2KP{D, S, P}
) :: VarInvIndexes{P} where {D, S, P}
	num_piece_types = convert(P, length(bp.l))
	num_plate_types = convert(P, length(bp.pli_lwb))
	return VarInvIndexes(num_piece_types, num_plate_types, bp.np, bp.cuts)
end

#=
# TODO: change enumeration to also use a Dict? and then use
# something like create_normalized_arg_subset in a subset of PPG2KP args
@timeit TIMER function _gen_cuts_wo(
	::Type{P}, d :: Vector{D}, l :: Vector{S}, w :: Vector{S},
	L :: S, W :: S, options :: Dict{String, Any} = Dict{String, Any}()
) where {D, S, P}
	return gen_cuts(P, d, SortedLinkedLW(D, l, w), L, W;
		ignore_2th_dim = options["ignore-2th-dim"],
		ignore_d = options["ignore-d"],
		round2disc = options["round2disc"],
		no_cut_position = options["no-cut-position"],
		no_redundant_cut = options["no-redundant-cut"],
		no_furini_symmbreak = options["no-furini-symmbreak"],
		faithful2furini2016 = options["faithful2furini2016"],
		quiet = options["quiet"], verbose = options["verbose"]
	)
end
=#

@timeit TIMER function _gen_cuts_wo(
	::Type{P}, d :: Vector{D}, sllw :: SortedLinkedLW{D, S},
	L :: S, W :: S, options :: Dict{String, Any} = Dict{String, Any}()
) where {D, S, P}
	return gen_cuts(P, d, sllw, L, W,
		Val(options["allow-rotation"] & options["mirror-plates"]);
		ignore_2th_dim = options["ignore-2th-dim"],
		ignore_d = options["ignore-d"],
		round2disc = options["round2disc"],
		hybridize_with_restricted = options["hybridize-with-restricted"],
		aggressive_hybridization = options["aggressive-hybridization"],
		no_cut_position = options["no-cut-position"],
		no_redundant_cut = options["no-redundant-cut"],
		no_furini_symmbreak = options["no-furini-symmbreak"],
		faithful2furini2016 = options["faithful2furini2016"],
		quiet = options["quiet"], verbose = options["verbose"],
	)
end

# INTERNAL USE.
#
# (Only used in _build_base_model! to make the shared demand constraint.)
#
# Helper functions that make it easier to index collections using the
# Vector{Union{D, Tuple{D, D}}} index lists.
function _flatten_index(a, i :: Integer, f)
	return a[i]
end
function _flatten_index(a, i :: Tuple{I, I}, f) where {I <: Integer}
	return f(a[first(i)], a[last(i)])
end

# INTERNAL USE.
#
# (Only used in _build_base_model! to compute a cut upper bound.)
#
# Helper function that computes a loose upper bound on the number of original
# plates available.
_ub_plates(inst :: MHLOPPW) = only(inst.available)
function _ub_plates(inst :: SSSCSP{D, S, P}) where {D, S, P}
	@unpack L, W, l, w, d = inst
	# The simplest (but terribly loose) upper bound is just `sum(d)`
	# (i.e., there is a solution in which each piece is in a plate of its
	# own). The following is the second-simplest thing: we have 4 buckets,
	# for pieces in which li <= L÷8 and wi <= W÷8 (we can fit 64 of these
	# in a plate), in which li <= L÷4 and wi <= W÷4 (we can fit 16 of these
	# in a plate), in which li <= L÷2 and wi <= W÷2 (we can fit 4 of these
	# in a plate), and in which li > L÷2 or wi > W÷2 (we give a whole plate
	# for these). This way, the bound is not terrible in the presence of
	# small- and medium-sized pieces.
	le8qt64, le4qt16, le2qt4, gt2 = zero(D), zero(D), zero(D), zero(D)

	L8, L4, L2 = L ÷ 8, L ÷ 4, L ÷ 2
	W8, W4, W2 = W ÷ 8, W ÷ 4, W ÷ 2
	for (li, wi, di) in zip(l, w, d)
		if     li <= L8 && wi <= W8
			le8qt64 += di
		elseif li <= L4 && wi <= W4
			le4qt16 += di
		elseif li <= L2 && wi <= W2
			le2qt4  += di
		else # i.e., li > L2 || wi > W2
			gt2     += di
		end
	end

	# The trick below is: if we sum (divisor - 1) to the number about to be
	# divided we transform ÷ into a round-up (not round-down) integer division.
	final_ub = (
		((le8qt64 + 63) ÷ 64) +
		((le4qt16 + 15) ÷ 16) +
		((le2qt4  +  3) ÷  4) +
		gt2
	)

	@assert final_ub <= sum(d)

	return final_ub
end

@timeit TIMER function _build_base_model!(
	problem :: SIMILAR_4, inst, model, bp :: ByproductPPG2KP{D, S, P},
	inv_idxs :: VarInvIndexes{P},
	rad :: Union{Nothing, RotationAwareData{D, S}},
	options :: Dict{String, Any} = Dict{String, Any}();
	build_LP_not_MIP = false
	# `build_LP_not_MIP` is used internally by the pricing (the relaxation needs
	# to be solved to solve the problem as a MIP).
	# `options["relax2lp"]` is a generic option to be handled in the generic
	# (i.e., formulation-agnostic) part of the code.
) where {D, S, P}
	# shorten the names
	@unpack d, l, w, L, W, np, cuts, cut_extraction, pli_lwb = bp
	@unpack pii2pair, pli2pair, child2cut, parent2cut = inv_idxs
	num_piece_types = convert(D, length(d))
	num_plate_types = length(pli_lwb)

	# This probably should not be here (but in an outer function)
	# but for now this will suffice.
	pii2dci = Vector{D}[Vector{D}() for i in 1:num_piece_types]
	for (dci, pii) in pairs(cut_extraction)
		!iszero(pii) && push!(pii2dci[pii], dci)
	end
	#@show length(cuts)
	#@show length(cut_extraction)

	if problem == Val(:G2KP) || problem == Val(:G2MKP)
		if options["allow-rotation"]
				# If a rotation-aware piece (i.e., dummy) maps to two original
				# pieces (only happens if both are perfect rotations of each other)
				# then we can map the dummy index to the any of the original pieces
				# (here we choose the first) because both share the same profit value
				# (this is a requirement to be a perfect rotation).
				rai2opi = first.(rad.sdi2opi[rad.rai2sdi]) :: Vector{D}
				p = inst.p[rai2opi] :: Vector{P}
		else
				p = inst.p :: Vector{P}
		end
	end

	# Note: using `d` as upper bound of `picuts` and `dc_sells` works even
	# with rotation enabled because `d` is the `rad.shared_d` expanded for the
	# dummies.
	@variable(model,
		0 <= picuts[i = 1:length(np)] <= d[np[i][2]], integer = !build_LP_not_MIP
	)
	if options["hybridize-with-restricted"]
		# dc_sells: Means "double-cut sells". For each piece type `pii` (or dummy
		# piece type in the case of rotation enabled), we have `dc_sells[pii]`
		# storing the number of pieces obtained by means of double cuts that were
		# actually sold. This variable is necessary because we cannot assume the
		# piece obtained by a double cut is always sold.  When double cuts are
		# enabled (i.e., hybridize-with-restricted), the normal cuts replaced by
		# double cuts cease to exist, and the double cuts may be needed in some
		# situations just to trim waste, without selling the piece they extract
		# because that piece demand is already exhausted.
		@variable(model,
			0 <= dc_sells[i = 1:num_piece_types] <= d[i],
			integer = !build_LP_not_MIP
		)
	end

	# We prefer to leave `cuts_made` without an upper bound to risk
	# making a bound that will not work for some combination of parameters.
	# EDIT: we are adding back a loose upper bound because it seems like
	# Gurobi have a bug that may give wrong answers if upper bounds are not set.
	# The `ub_cut_in_single_plate` is just an upper bound on the amount of
	# plates p that fit in a single original plate, where p is the parent plate
	# of the cut (i.e., the plate that is being divided by the cut).
	ub_cut_in_single_plate = getindex.(pli_lwb, 3)[getindex.(cuts, 1)]
	ub_cut = if problem === Val(:G2MKP) || problem === Val(:G2CSP)
		# NOTE: _ub_plates assumes no-rotation because this is a valid upper bound
		# for both rotation and no-rotation variants.
		ub_cut_in_single_plate .* _ub_plates(inst)
	else
		ub_cut_in_single_plate
	end
	@variable(model,
		0 <= cuts_made[i = 1:length(cuts)] <= ub_cut[i],
		integer = !build_LP_not_MIP
	)

	if options["hybridize-with-restricted"]
		# Constraint that guarantee dc_sells only sell pieces actually extracted
		# by double cuts.
		@constraint(model,
			dc_link_con[i = 1:length(d)], dc_sells[i] <= sum(cuts_made[pii2dci[i]])
		)
	end

	# If the problem is a knapsack problem, then the objective is to maximize
	# the profit from the packed/extracted pieces.
	# If the problem is CSP, then the objective is to minimize the number of
	# original plates used (all pieces must be packed).
	# If it is OPP, then there is no objective function.
	if problem === Val(:G2KP) || problem === Val(:G2MKP)
		# If the objective is profit maximization we need to account for the
		# sold pieces obtained by double cut extractionsa (`dc_sells`).
		if options["hybridize-with-restricted"]
			@objective(model, Max,
				sum(p[pii] * sum(picuts[pii2pair[pii]]) for pii = 1:num_piece_types) +
				sum(p[pii] * dc_sells[pii] for pii = 1:num_piece_types)
			)
		else
			@objective(model, Max,
				sum(p[pii] * sum(picuts[pii2pair[pii]]) for pii = 1:num_piece_types)
			)
		end
	elseif problem === Val(:G2CSP)
		# We create a variable to represent the number of original plates used:
		# `b`. This variable is not strictly necessary, we could instead minimize
		# the cuts and extractions made over plates of plate type #1. However,
		# without the variable we will have one less constraint compared to the
		# other problems (there will not be the `plate_cons` #1 constraint). This
		# would add some unfortunate assymetry; either the list of constraint would
		# start at the second index, or to map a plate to its constraint we would
		# need to always subtract one of the plate index (when it is G2CSP,
		# otherwise not, therefore a mess). This seems to be the better design,
		# no code (even G2CSP-specific) really need to be aware of this variable.
		@variable(model, b, integer = !build_LP_not_MIP)
		#@objective(model, Min, sum([b])) # the sum prevents errors with GLPK
		@objective(model, Min, b)
	end

	# c1: # Either there is a limited amount of original plates available
	# (one if it is KP or OPP, and 'n' if it is MKP), or there is an
	# unlimited amount of original plates available (CSP) but we are minimizing
	# the amount actually used (variable `b`).
	#
	# The original plates may be used either to extract pieces directly
	# or to be cut into two new subplates.
	if problem === Val(:G2CSP)
		plate_number_one = @constraint(model,
			sum(picuts[pli2pair[1]]) + sum(cuts_made[parent2cut[1]]) <= b
		)
	else
		qt_original :: Int = problem === Val(:G2MKP) ? only(inst.available) : 1
		plate_number_one = @constraint(model,
			sum(picuts[pli2pair[1]]) + sum(cuts_made[parent2cut[1]]) <= qt_original
		)
	end

	# c2: for each subplate type that is not the original plate, such subplate
	# type will be available the number of times it was the child of a cut,
	# subtracted the number of times it had a piece extracted or used for
	# further cutting. This constraint is the basis of the model and similar
	# for all problems it supports.
	# NOTE: the range is 1:(num_plate_types-1) and not 2:num_plate_types
	# because if the range does not start at 1, the returned type is not
	# a vector but a special JuMP type.
	plate_numbers_2_to_m = @constraint(model, [ppli=1:(num_plate_types-1)],
		sum(picuts[pli2pair[ppli+1]]) + sum(cuts_made[parent2cut[ppli+1]]) <=
		sum(cuts_made[child2cut[ppli+1]])
	)
	# give the constraints a name
	set_name.(plate_numbers_2_to_m, "plate_cons" .* string.(2:num_plate_types))
	# put them in the expected key inside the model
	model[:plate_cons] = pushfirst!(plate_numbers_2_to_m, plate_number_one)

	# The demand constraint, yet simple, is modified in complementary ways by
	# basically every flag that changes the formulation, leading an exponential
	# explosion of slightly different demand constraints. To avoid 3+ levels of
	# nested `if`s, we try to build each part of the demand constraint
	# separatedly (depending only in the options that affect that single part)
	# and then combine the parts at the end.

	# The right-hand-side of the demand constraint is determined by the
	# (no-)rotation. The (no-)rotation also defines how we translate
	# demand constraint indexes to piece indexes (`di2pii`).
	if options["allow-rotation"]
		demand_con_rhs = rad.shared_demand
		di2pii = rad.sdi2rai :: fieldtype(RotationAwareData{D, S}, :sdi2rai)
	else
		demand_con_rhs = d
		# We use typeof to avoid type-unstability on binding `di2pii`.
		di2pii = fieldtype(RotationAwareData{D, S}, :sdi2rai)()
		append!(di2pii, eachindex(d))
	end
	# With the rhs and the translator we may build lhs expression outside
	# of any `if`.
	demand_con_lhs = @expression(model, [di=1:length(di2pii)],
		sum(picuts[_flatten_index(pii2pair, di2pii[di], vcat)])
	)
	# If we are hybridizing with restricted we need to update lhs to
	# also include the extractions/sells made by double cuts.
	if options["hybridize-with-restricted"]
		demand_con_lhs = @expression(model, [di=1:length(di2pii)],
			demand_con_lhs[di] + sum(dc_sells[collect(di2pii[di])])
		)
	end
	# Now, with lhs and rhs in hand, we will build contraints that only differ in
	# the sense (i.e., either `<=` or `>=`).
	# For knapsack problems, demand means an upper bound on how many of each
	# piece type may be extracted/packed, for OPP and CSP it means a lower bound
	# on the number of each piece type that *must* be extracted/packed.
	if problem === Val(:G2KP) || problem === Val(:G2MKP)
		@constraint(model,
			demand_con[di=1:length(di2pii)],
			demand_con_lhs[di] <= demand_con_rhs[di]
		)
	elseif problem === Val(:G2OPP) || problem === Val(:G2CSP)
		@constraint(model,
			demand_con[di=1:length(di2pii)],
			demand_con_lhs[di] >= demand_con_rhs[di]
		)
	end

	#= Disabled while we rework it in a problem-agnostic fashion.
	lb = options["lower-bound"]
	if !iszero(lb)
		@constraint(model, obj_lb_con,
			sum(p[pii]*sum(picuts[pii2pair[pii]]) for pii = 1:num_piece_types) >= (lb + 1)
		)
	end

	ub = options["upper-bound"]
	if ub < sum(d .* p)
		@constraint(model, obj_ub_con,
			sum(p[pii]*sum(picuts[pii2pair[pii]]) for pii = 1:num_piece_types) <= ub
		)
	end
	=#

	return
end

# HIGH LEVEL EXPLANATION OF THE MODEL
#
# Variables:
#
# `picuts[n, pii]`: Integer. The number of pieces `pii` generated from
#   subplates of type `n`.
# `cuts_made[n1, n2, n3]`: Integer. The number of subplates of type
#   `n1` that are cut into subplates `n2` and `n3` (horizontal and
#   vertical cuts are together for now).
#
# Objective function:
#
# Maximize the profit of the pieces cut.
#   sum(p[pii] * picuts[_, pii])
#
# Constraints:
#
# There is exactly one of the original plate, which may be used for cutting
# or extracting a piece.
#   sum(picuts[1, _]) + sum(cuts_made[1, _,  _]) <= 1
# The number of subplates available depends on the number of plates that have
# it as children.
#   sum(picuts[n1>1, _]) + sum(cuts_made[n1>1, _, _]) <=
#     sum(cuts_made[_, n2, n3])
#     where n2 == n1 or n3 == n1, doubling cuts_made[_, n2, n3] if n2 == n3
# The number of pieces of some type is always less than or equal to the demand.
#   sum(picuts[_, pii]) <= d[pii]
#
# Unnecessary constraints:
#
# The number of times a pair pli-pii appear is at most the min between:
# d[pii] and the number of subplates pli that fit in the original plate.
#   sum(picuts[n, pii]) <= min(d[pii], max_fits[n])
"""
    build_complete_model(model, d, p, l, w, L, W[, start]; options)

TODO: describe method.
Note: `model[:plate_cons]` is a container of constraints of type and set
`(GenericAffExpr{Float64,VariableRef}, LessThan{Float64})` that represent all
the plate types (i.e., it will be the same length as ByproductPPG2KP.pli_lwb
and each position in `model[:plate_cons]` will correspond to the plate
described in the same position of `ByproductPPG2KP.pli_lwb`).
"""
#=
function build_complete_model(
	model, d :: Vector{D}, p :: Vector{P}, l :: Vector{S}, w :: Vector{S},
	L :: S, W :: S, start :: Float64 = time(),
	options :: Dict{String, Any} = Dict{String, Any}()
) :: ByproductPPG2KP{D, S, P} where {D, S, P}
	bp = _gen_cuts_wo(P, d, l, w, L, W, options)
	build_complete_model(model, p, bp, start, options)
	return bp
end
=#
function build_complete_model(
	problem, instance, model, bp :: ByproductPPG2KP{D, S, P},
	rad :: Union{Nothing, RotationAwareData{D, S}},
	start :: Float64 = time(), options :: Dict{String, Any} = Dict{String, Any}()
) :: ByproductPPG2KP{D, S, P} where {D, S, P}
	inv_idxs = VarInvIndexes(bp)
	limit :: Float64 = options["building-time-limit"]
	throw_if_timeout_now(start, limit)
	@assert issorted(
		@view bp.cuts[1:(bp.first_vertical_cut_idx-1)];
		by = c -> c[PARENT]
	)
	@assert	issorted(
		@view bp.cuts[bp.first_vertical_cut_idx:end];
		by = c -> c[PARENT]
	)
	_build_base_model!(problem, instance, model, bp, inv_idxs, rad, options)

	return bp
end

# Include all the code implementing the pricings referred below.
include("./pricing.jl")

# Auxiliar function to get the geometric info from the homogeneous
# problems and pass it to the enumeration.
function _homo_geometry(
	i :: Union{G2KP{D, S, P}, SSSCSP{D, S, P}}
) where {D, S, P}
	return D, S, P, i.d, i.l, i.w, i.L, i.W
end
function _homo_geometry(i :: MHLOPPW{D, S, P}) where {D, S, P}
	return D, S, P, i.dub, i.l, i.w, only(i.L), only(i.W)
end
function _get_profit(
	i :: Union{G2KP{D, S, P}, MHLOPPW{D, S, P}}
) :: Vector{P} where {D, S, P}
	return i.p
end
function _get_profit(
	i :: SSSCSP{D, S, P}
) :: Vector{P} where {D, S, P}
	return zeros(P, length(i.d))
end

function _no_arg_check_build_model(
	problem, instance, model, options :: Dict{String, Any}
) :: Tuple{BuildStopReason, <: Any}
	# Define 'start', this is the point where the building-time-limit
	# starts to be counted.
	start :: Float64 = time()
	limit :: Float64 = options["building-time-limit"]
	verbose = options["verbose"] && !options["quiet"]
	# Enumerate plates and cuts.
	enumeration_time = @elapsed begin
		D, _, P, d, l, w, L, W = _homo_geometry(instance)
		sllw = SortedLinkedLW(D, l, w)
		if options["allow-rotation"]
			rad = build_RAD(L, W, sllw, d, _get_profit(instance))
			# Here we expand the shared demand to all dummy pieces it applies,
			# so it is invisible to _gen_cuts_wo that we are dealing with
			# a rotation variant.
			unshared_demand = rad.shared_demand[rad.rai2sdi]
			bp = _gen_cuts_wo(P, unshared_demand, rad.ra_sllw, L, W, options)
		else
			bp = _gen_cuts_wo(P, d, sllw, L, W, options)
		end
	end
	if verbose
		@show enumeration_time
		#println("START_DEBUG_T1a")
		#for (pp, fc, sc) in bp.cuts
		#	print(bp.pli_lwb[pp][1], ' ', bp.pli_lwb[pp][2], ' ')
		#	print(bp.pli_lwb[fc][1], ' ', bp.pli_lwb[fc][2], ' ')
		#	!iszero(sc) && print(bp.pli_lwb[sc][1], ' ', bp.pli_lwb[sc][2])
		#	println()
		#end
		#println("END_DEBUG_T1a")
		println("qt_pevars_after_preprocess = $(length(bp.np))")
		println("qt_cmvars_after_preprocess = $(length(bp.cuts))")
		println("qt_plates_after_preprocess = $(length(bp.pli_lwb))")
	end
	throw_if_timeout_now(start, limit)

	# Deal with pricing or its absence. Builds the model based on the enumerated
	# plates and cuts (or a subset of them).
	switch_method = options["Gurobi-LP-method-inside-furini-pricing"]
	pricing = options["pricing"]
	# The warnings are done here to avoid putting them inside a method
	# that will not be called anyway.
	current_solver_name = solver_name(model)
	if !options["quiet"] && switch_method > -2
		if pricing != "furini"
			@warn "The flag Gurobi-LP-method-inside-furini-pricing has a" *
				" non-default value, but the pricing is '$pricing'. This flag" *
				" will have no effect besides warnings."
		end
		if current_solver_name != "Gurobi"
			@warn "The flag Gurobi-LP-method-inside-furini-pricing has a" *
				" non-default value, but the solver is '$current_solver_name'." *
				" This flag will have no effect besides warnings."
		end
	end
	bm = (options["faithful2furini2016"] ? FURINI : BECKER) :: BaseModel
	if pricing != "none"
		problem !== Val(:G2KP) && @error(
			"The pricing was only implemented for the G2KP problem. Other" *
			" problems have no implementation of the pricing."
		)
		options["allow-rotation"] && @error(
			"The pricing code was not yet examined to check if it supports" *
			" rotation. The combination of any pricing method plus rotation" *
			" is disabled for now."
		)
		# It is important to note that the pricing and the bm (BaseModel) have the
		# same two values ("becker"/BECKER and "furini"/FURINI) but they are
		# orthogonal/independent. You can have a FURINI model with "becker" pricing
		# and vice-versa. It is just because each author has defined both a model
		# and an optional pricing method (to be executed over it) but both pricings
		# are compatible with both models, and we allow them to be mixed. The
		# pricing methods just need to know the underlying base model to be able to
		# MIP-start it corretly (initially at least, this comment may be out of
		# date).
		pricing_time = @elapsed begin
			if pricing == "furini"
				bsr, mbp = _furini_pricing!(model, bp, p, start, options)
				bp = mbp.preprocess_byproduct
			else
				@assert pricing == "becker"
				bp = _becker_pricing!(model, bp, p, start, options)
			end
		end
		verbose && @show pricing_time
		# Could be cleaner...
		if @isdefined(bsr) && bsr == FOUND_OPTIMUM
			@assert @isdefined(mbp)
			return bsr, mbp
		end
	else # in the case there is no pricing phase
		# Needs to build the mode here, as when there is pricing the pricing
		# procedure is responsible for building the model.
		build_complete_model(
			problem, instance, model, bp,
			(options["allow-rotation"] ? rad : nothing), start, options
		)
		if options["MIP-start"] == "guaranteed"
			problem !== Val(:G2KP) && @error(
				"The MIP-start was only implemented for the G2KP problem. Other" *
				" problems have no implementation of the MIP-start."
			)
			options["allow-rotation"] && @error(
				"The MIP-start was not yet examined to check if it supports" *
				" rotation. Using MIP-start plus rotation is disabled for now."
			)
			heuristic_lb_time = @elapsed begin
				(heuristic_lb, _, _), _ = mip_start_by_heuristic!(
					model, bp, p, options["heuristic-seed"], bm
				)
			end
			if verbose
				@show heuristic_lb_time
				@show heuristic_lb
			end
		end
	end
	if verbose
		println("length_pe_after_pricing = $(length(model[:picuts]))")
		println("length_cm_after_pricing = $(length(model[:cuts_made]))")
		println("length_pc_after_pricing = $(length(model[:plate_cons]))")
	end
	throw_if_timeout_now(start, limit) # after the pricings or warm-starts

	# Check if variables and constraints made useless by the deletion
	# of other variables will be kept or not.
	# TODO: if we ever start dealing with multiple original plates, we may
	# need to update _handle_unreachable! (the exception is if we go with
	# some maximum.((L, W)) workaround were the original plates are all
	# cut from a dummy root plate).
	purge_unreachable_time = @elapsed begin
		purge = !options["do-not-purge-unreachable"]
		bp = _handle_unreachable!(bp, model, verbose, purge)
		if verbose && purge
			println("qt_pevars_after_purge = $(length(model[:picuts]))")
			println("qt_cmvars_after_purge = $(length(model[:cuts_made]))")
			println("qt_plates_after_purge = $(length(model[:plate_cons]))")
		end
	end
	verbose && @show purge_unreachable_time
	throw_if_timeout_now(start, limit) # after handle_unreachable!

	options["allow-rotation"] && return BUILT_MODEL, ModelByproduct(bp, rad)

	return BUILT_MODEL, ModelByproduct(bp)
end

# Only used to define `build_model(Val{:PPG2KP}, ...)`.
import ..Utilities.Args.create_normalized_arg_subset
import ..Utilities.Args.accepted_arg_list
import ..build_model
#=
"""
    build_model(::Val{:G2KP}, ::Val{:PPG2KP}, instance::G2KP, model[, options])
    build_model(::Val{:G2MKP}, ::Val{:PPG2KP}, instance::MHLOPPW, model[, options])
    build_model(::Val{:G2CSP}, ::Val{:PPG2KP}, instance::SSSCSP, model[, options])
    build_model(::Val{:G2OPP}, ::Val{:PPG2KP}, instance::SSSCSP, model[, options])

Build a PPG2KP-style model for a G2KP `instance` inside `model`.

Changes `model` by adding variables and constraints. `options` takes
arguments described in `Utilities.Args.accepted_arg_list(::Val{:PPG2KP})`.
"""
=#
#=
@timeit TIMER function build_model(
	HOMO_PROBLEMS, ::Val{:PPG2KP}, instance :: G2KP{D, S, P}, model,
	options :: Dict{String, Any} = Dict{String, Any}()
) :: Tuple{BuildStopReason, <: Any} where {D, S, P}
	norm_options = create_normalized_arg_subset(
		options, accepted_arg_list(Val(:PPG2KP))
	)
	return _no_arg_check_build_model(model, norm_options)
end
=#

const HOMO_PAIRS = [
	(:G2KP, :G2KP), (:G2CSP, :SSSCSP), (:G2OPP, :SSSCSP)
]

for (problem, instance) in HOMO_PAIRS
	eval(quote
		@timeit TIMER function build_model(
			::Val{$(QuoteNode(problem))}, ::Val{:PPG2KP}, instance :: $instance{D, S, P}, model,
			options :: Dict{String, Any} = Dict{String, Any}()
		) :: Tuple{BuildStopReason, <: Any} where {D, S, P}
			norm_options = create_normalized_arg_subset(
				options, accepted_arg_list(Val(:PPG2KP))
			)
			return _no_arg_check_build_model(
				Val($(QuoteNode(problem))), instance, model, norm_options
			)
		end
	end)
end

"""
    build_model(::Val{:G2KP}, ::Val{:PPG2KP}, instance::SLOPP, model[, options])

Build a PPG2KP-style model for a G2KP `instance` inside `model`.

Changes `model` by adding variables and constraints. `options` takes
arguments described in `Utilities.Args.accepted_arg_list(::Val{:PPG2KP})`.

This is a convenience method that takes a SLOPP instance. It errors if
the `SLOPP` object has a non-zero value in the `dlb` field. It is
the same as calling with a `G2KP` instance but with the `dub` field
being used as the `d` field (of `G2KP`).
"""
@timeit TIMER function build_model(
	::Val{:G2KP}, ::Val{:PPG2KP}, instance :: SLOPP{D, S, P}, m,
	options :: Dict{String, Any} = Dict{String, Any}()
) :: Tuple{BuildStopReason, <: Any} where {D, S, P}
	@assert all(iszero, instance.dlb)
	@unpack p, l, w, L, W = instance
	d = instance.dub
	return build_model(
		Val(:G2KP), Val(:PPG2KP), G2KP(L, W, l, w, d, p), m, options
	)
end

"""
    build_model(::Val{:G2CSP}, ::Val{:PPG2KP}, instance::G2KP, model[, options])

Build a PPG2KP-style model for a G2CSP `instance` inside `model`.

Changes `model` by adding variables and constraints. `options` takes
arguments described in `Utilities.Args.accepted_arg_list(::Val{:PPG2KP})`.

This is a convenience method that takes a G2KP instance. It allows
easily testing G2CSP over the same instances used commonly for
Knapsack Problem (G2KP). If `options["quiet"]` is not passed,
this method will always show a warnig, as the instance type is
not ideal.
"""
@timeit TIMER function build_model(
	::Val{:G2CSP}, ::Val{:PPG2KP}, instance :: G2KP{D, S, P}, m,
	options :: Dict{String, Any} = Dict{String, Any}()
) :: Tuple{BuildStopReason, <: Any} where {D, S, P}
	@unpack L, W, l, w, d = instance
	if !options["quiet"]
		@warn "Attention: converting a G2KP instance to a SSSCSP instance to" *
			" solve it for the G2CSP (problem). Profits dropped."
	end
	return build_model(
		Val(:G2CSP), Val(:PPG2KP), SSSCSP{D, S, P}(L, W, l, w, d), m, options
	)
end

"""
    build_model(::Val{:G2KP}, ::Val{:PPG2KP}, instance::SSSCSP, model[, options])

Build a PPG2KP-style model for a G2OPP `instance` inside `model`.

Changes `model` by adding variables and constraints. `options` takes
arguments described in `Utilities.Args.accepted_arg_list(::Val{:PPG2KP})`.

THIS IS A HACK THAT SHOULD BE DELETED AFTER. It gets a SSSCSP instance, and
does not interpret it as a G2CSP/G2BPP instance, but instead as a G2OPP
instance and then it transforms it in a G2KP problem with unitaty profits,
so the G2OPP (a decision problem) is true if obj matches the sum(d), and
is false if obj is smaller than sum(d) (a value larger than sum(d) is
impossible).
"""
@timeit TIMER function build_model(
	::Val{:G2KP}, ::Val{:PPG2KP}, instance :: SSSCSP{D, S, P}, m,
	options :: Dict{String, Any} = Dict{String, Any}()
) :: Tuple{BuildStopReason, <: Any} where {D, S, P}
	@unpack L, W, l, w, d = instance
	if !options["quiet"]
		@warn "Attention: converting a SSSCSP (understood as a G2OPP) instance to a G2KP instance to" *
			" solve it for the G2KP (problem). Unitary profits added."
	end
	return build_model(
		Val(:G2KP), Val(:PPG2KP), G2KP{D, S, P}(L, W, l, w, d, ones(P, length(d))), m, options
	)
end

"""
    build_model(::Val{:G2MKP}, ::Val{:PPG2KP}, instance::MHLOPPW, model[, options])

Build a PPG2KP-style model for a G2MKP `instance` inside `model`.

Changes `model` by adding variables and constraints. `options` takes
arguments described in `Utilities.Args.accepted_arg_list(::Val{:PPG2KP})`.

For now, this is the only format accepted by G2MKP, even if the G2MKP
refers to the homogeneous variant and, therefore, does not accept
large objects of different dimensions (just multiple copies of the
same large objects). A method accepting a more adequate format
(i.e., that only allow a single type of large object) could be
implemented.
"""
@timeit TIMER function build_model(
	::Val{:G2MKP}, ::Val{:PPG2KP}, instance :: MHLOPPW{D, S, P}, m,
	options :: Dict{String, Any} = Dict{String, Any}()
) :: Tuple{BuildStopReason, <: Any} where {D, S, P}
	if (length(instance.L), length(instance.W)) != (1, 1)
		error(
			"The G2MKP problem refers to the homogeneous variant." *
			" Multiple original plates (i.e., large objects) of distinct" *
			" dimensions are not allowed."
		)
	end
	norm_options = create_normalized_arg_subset(
		options, accepted_arg_list(Val(:PPG2KP))
	)
	return _no_arg_check_build_model(
		Val(:G2MKP), instance, m, norm_options
	)
end

end # module


text:
"""
    fastest_exact(::CombinatorialInstance)

Returns the overall fastest algorithm to solve the combinatorial instance
(i.e. this function returns a `CombinatorialAlgorithm` instance). If no such
algorithm is known, then the return value is `nothing` (for instance, there is
no algorithm that is always faster than the others). The returned algorithm
is always exact, and never approximate.
"""
function fastest_exact(::CombinatorialInstance)
    return nothing
end

# ============================================================
# = First, generic-named algorithms (dynamic, greedy, etc.). =
# ============================================================

"""
    struct DynamicProgramming <: CombinatorialAlgorithm

Use dynamic programming to solve the instance.

Some DP algorithms have a common name, such as `BellmanFordAlgorithm`.
"""
struct DynamicProgramming <: CombinatorialAlgorithm end

"""
    struct GreedyAlgorithm <: CombinatorialAlgorithm

Use a greedy algorithm to solve the instance.

Some greedy algorithms have a common name, such as `PrimAlgorithm`.
"""
struct GreedyAlgorithm <: CombinatorialAlgorithm end

"""
    struct LagrangianAlgorithm <: CombinatorialAlgorithm

Use Lagrangian relaxation to solve the instance. This is especially useful
to build an algorithm for a problem with a new constraint based on other
algorithms. These techniques do not always yield good approximation guarantees.
"""
struct LagrangianAlgorithm <: CombinatorialAlgorithm end

"""
    struct LagrangianAlgorithm <: CombinatorialAlgorithm

Use Lagrangian relaxation to solve the instance, followed by a refinement step.
This is especially useful to build an algorithm for a problem with a new
constraint based on other algorithms. These techniques often yield an constant
additive approximation, based on Lagrangian relaxation.
"""
struct LagrangianRefinementAlgorithm <: CombinatorialAlgorithm end

"""
    struct LagrangianAlgorithm <: CombinatorialAlgorithm

Use Lagrangian relaxation to solve the instance, followed by a refinement step
and more iterations as required. This is especially useful to build an algorithm
for a problem with a new constraint based on other algorithms. These techniques
often yield an constant multiplicative approximation, based on Lagrangian
relaxation and refinement.
"""
struct IteratedLagrangianRefinementAlgorithm <: CombinatorialAlgorithm end

"""
    struct DefaultLinearFormulation <: CombinatorialAlgorithm

Use the default LP formulation to solve the instance. For problems where
there are multiple common formulations, there should be no default formulation
defined.
"""
struct DefaultLinearFormulation <: CombinatorialLinearFormulation
    solver
end

# ===============================================================
# = Then, the ones with a recognised name (think Bellman-Ford). =
# ===============================================================

"""
    struct BellmanFordAlgorithm <: CombinatorialAlgorithm

Use the Bellman-Ford algorithm to solve the instance. Bellman-Ford is a special
case of dynamic programming.

This algorithm is sometimes called the Bellman-Ford-Moore algorithm.
"""
struct BellmanFordAlgorithm <: CombinatorialAlgorithm end

"""
    struct HungarianAlgorithm <: CombinatorialAlgorithm

Use the Hungarian algorithm to solve the instance.

This algorithm is sometimes called the Kuhn-Munkres algorithm.
"""
struct HungarianAlgorithm <: CombinatorialAlgorithm end

"""
    struct PrimAlgorithm <: CombinatorialAlgorithm

Use Prim's algorithm to solve the instance. Prim is a special case of greedy
algorithm.
"""
struct PrimAlgorithm <: CombinatorialAlgorithm end

# ==============================================================
# = Functions about the approximation guarantees of algorithms =
# ==============================================================

"""
    approximation_ratio(::CombinatorialInstance, ::CombinatorialAlgorithm)

Returns the approximation ratio for this algorithm when run on the input instance. When the problem is minimising
a function (e.g., finding the minimum-cost path), the ratio is defined as one constant ``r \\le 1.0``
(ideally, the lowest) such that

``f(x^\\star) \\leq f(x) \\leq r \\cdot f(x^\\star),``

where ``f(x^\\star)`` is the cost of the optimum solution and ``f(x)`` the one of the returned solution. On the
contrary, for a maximisation problem, the definition is reversed (with ``r \\le 1.0``):

``r \\cdot f(x^\\star) \\geq f(x) \\geq f(x^\\star).``

For an exact algorithm, the ratio is always ``1.0``, for both minimisation and maximisation problems.

The returned ratio might be constant, if the algorithm provides a constant ratio; if the ratio is not constant (i.e.
instance-dependent), it may either be a worst-case value or a truly instance-dependent ratio. Depending on the
algorithm, this behaviour might be tuneable.

If the algorithm has no guarantee, it should return `NaN`.
"""
function approximation_ratio(::CombinatorialInstance, ::CombinatorialAlgorithm)
    return 1.0
end

"""
    approximation_term(::CombinatorialInstance, ::CombinatorialAlgorithm)

Returns the approximation term for this algorithm when run on the input instance. When the problem is minimising
a function (e.g., finding the minimum-cost path), the term is defined as one constant ``t \\ge 0.0``
(ideally, the lowest) such that

``f(x^\\star) \\leq f(x) \\leq f(x^\\star) + t,``

where ``f(x^\\star)`` is the cost of the optimum solution and ``f(x)`` the one of the returned solution. On the
contrary, for a maximisation problem, the definition is reversed (with ``t \\ge 0.0``):

``f(x^\\star) \\leq f(x) \\leq f(x^\\star) - t.``

For an exact algorithm, the term is always ``0.0``, for both minimisation and maximisation problems.

The returned term might be constant, if the algorithm provides a constant term; if the term is not constant (i.e.
instance-dependent), it may either be a worst-case value or a truly instance-dependent term. Depending on the
algorithm, this behaviour might be tuneable.

If the algorithm has no guarantee, it should return `NaN`.
"""
function approximation_term(::CombinatorialInstance, ::CombinatorialAlgorithm)
    return 0.0
end


text:
# ------------------------------------------------------------------
# Copyright (c) 2017, Júlio Hoffimann Mendes <juliohm@stanford.edu>
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    BallNeighborhood(domain, radius)

A ball neighborhood of a given `radius` on a spatial `domain`.
"""
struct BallNeighborhood{D<:AbstractDomain} <: AbstractNeighborhood{D}
  # input fields
  domain::D
  radius # we cannot use coordtype(D) here yet in Julia v0.6

  # state fields
  cube::CubeNeighborhood{D}

  function BallNeighborhood{D}(domain, radius) where D<:AbstractDomain
    @assert radius > 0 "radius must be positive"
    @assert typeof(radius) == coordtype(domain) "radius and domain coordinate type must match"

    # cube of same radius
    cube = CubeNeighborhood(domain, radius)

    new(domain, radius, cube)
  end
end

BallNeighborhood(domain::D, radius) where {D<:AbstractDomain} = BallNeighborhood{D}(domain, radius)

function (neigh::BallNeighborhood{<:RegularGrid})(location::Int)
  # retrieve domain
  ndomain = neigh.domain

  # get neighbors in cube of same radius
  cneighbors = neigh.cube(location)

  # coordinates of the center
  xₒ = coordinates(ndomain, location)

  # pre-allocate memory for neighbors coordinates
  x = MVector{ndims(ndomain),coordtype(ndomain)}(undef)

  # discard neighbors outside of sphere
  neighbors = Vector{Int}()
  for neighbor in cneighbors
    coordinates!(x, ndomain, neighbor)

    # compute ||x-xₒ||^2
    sum² = zero(eltype(x))
    for i=1:ndims(ndomain)
      sum² += (x[i] - xₒ[i])^2
    end

    sum² ≤ neigh.radius^2 && push!(neighbors, neighbor)
  end

  neighbors
end

function (neigh::BallNeighborhood{<:PointSet})(location::Int)
  # retrieve domain
  ndomain = neigh.domain

  # center in real coordinates
  xₒ = coordinates(ndomain, location)

  # pre-allocate memory for neighbors coordinates
  x = MVector{ndims(ndomain),coordtype(ndomain)}(undef)

  neighbors = Vector{Int}()
  for loc in 1:npoints(ndomain)
    coordinates!(x, ndomain, loc)
    norm(x .- xₒ) ≤ neigh.radius && push!(neighbors, loc)
  end

  neighbors
end


text:
"""
    VonMisesFisher(M; params...)

The von Mises-Fisher (vMF) or Langevin distribution on the `Sphere` or `Stiefel` manifold `M`.

Given a matrix ``X ∈ 𝔽^{n × k}`` with IID entries ``X_{ij} ∼ \\mathrm{Normal}(F_{ij}, 1)``
for ``F ∈ 𝔽^{n × k}``, the restriction of the corresponding distribution in ``𝔽^{n × k}``
to the `Stiefel(n, k, 𝔽)` manifold, that is, the matrices for which ``X^\\mathrm{H} X = I_k``,
is the vMF distribution on the `Stiefel` manifold. The vMF distribution can also be specified
for any submanifold of the `Stiefel` manifold, including the `Sphere` and the `Circle`.

# Parameterizations

    VonMisesFisher(M::AbstractSphere{𝔽}; params...)
    VonMisesFisher(n::Int[, 𝔽]; params...)

Construct the vMF distribution on `Sphere(n-1,𝔽)==` ``𝔽𝕊^{n-1}```.

Implemented parameterizations are:

  - `(μ, κ)`: the modal direction ``μ ∈ 𝔽𝕊^{n-1}`` and concentration ``κ ∈ ℝ⁺``
  - `(c,)`: ``c = κ μ ∈ 𝔽^n``, the mean of the normal distribution in the embedded space.

The density of the vMF distribution on ``𝔽𝕊^{n-1}`` with respect to the normalized
Hausdorff measure is

```math
p(x | μ, κ) = \\frac{Γ(ν + 1)κ^ν}{2^ν I_ν(κ)} \\exp(κ \\Re⟨μ, x⟩)),
```

where ``ν = n/2-1``,  ``⟨⋅,⋅⟩`` is the Frobenius inner product, and ``I_ν(z)``
is the modified Bessel function of the first kind.

    VonMisesFisher(M::Stiefel{n,k,𝔽}; params...)
    VonMisesFisher(n::Int, k::Int[, 𝔽]; params...)

Construct the matrix vMF distribution on `Stiefel(n, k, 𝔽)=` ``\\mathrm{St}(n, k, 𝔽)``.

Implemented parameterizations are:

  - `(F,)`: a parameter matrix ``F ∈ 𝔽^{n × k}``, the mean of the normal distribution in the
    embedded space.
  - `(U, D, V)`: The SVD decomposition of ``F = U D V``, where ``U ∈ \\mathrm{St}(n, k, 𝔽)`` and
    ``V ∈ \\mathrm{U}(k, 𝔽)``.
  - `(H, P)`: The polar decomposition of ``F = H P``, where ``H ∈ \\mathrm{St}(n, k, 𝔽)`` is the
    mode, and ``P ∈ 𝔽^{k × k}`` is a Hermitian positive definite matrix.

The density of the vMF distribution on `Stiefel(n, k, 𝔽)` with respect to the normalized
Hausdorff measure is

```math
p(x | F) = \\frac{\\exp(\\Re⟨F, x⟩)}{_0 F_1(\\frac{n}{2}; \\frac{1}{4} F^\\mathrm{H}F)},
```

where ``_0 F_1(b; B)`` is a hypergeometric function with matrix argument ``B``.
"""
struct VonMisesFisher{M,N,T} <: ParameterizedMeasure{N}
    manifold::M
    par::NamedTuple{N,T}
end
VonMisesFisher(M; params...) = VonMisesFisher(M, NamedTuple(params))
VonMisesFisher(p::Int, 𝔽=ℝ; params...) = VonMisesFisher(Sphere(p - 1, 𝔽); params...)
function VonMisesFisher(n::Int, k::Int, 𝔽::AbstractNumbers=ℝ; params...)
    return VonMisesFisher(Stiefel(n, k, 𝔽); params...)
end

"""
    VonMises(𝔽=ℝ; params...)

The von Mises distribution on the real or complex `Circle`.

This is just a convenient alias for [`VonMisesFisher(Circle(𝔽); params...)`](@ref).

# Constructors

    VonMises(μ=π, κ=1)
    VonMises(ℂ; μ=im, κ=1)
    VonMises(ℂ; c=3+4im)
"""
const VonMises{𝔽,N,T} = VonMisesFisher{Circle{𝔽},N,T}
VonMises(𝔽=ℝ; params...) = VonMisesFisher(Circle(𝔽); params...)

function Base.show(io::IO, mime::MIME"text/plain", μ::VonMisesFisher)
    return show_manifold_measure(io, mime, μ)
end

Manifolds.base_manifold(d::VonMisesFisher) = getfield(d, :manifold)

function MeasureTheory.basemeasure(μ::VonMisesFisher)
    return normalize(Hausdorff(base_manifold(μ)))
end

function MeasureTheory.logdensity(d::VonMises{ℝ,(:μ, :κ)}, x)
    κ = d.κ
    return κ * cos(only(x) - only(d.μ)) - logbesseli(0, κ)
end

function MeasureTheory.logdensity(d::VonMisesFisher{M,(:μ, :κ)}, x) where {M}
    p = manifold_dimension(base_manifold(d)) + 1
    κ = d.κ
    return κ * real(dot(d.μ, x)) - lognorm_vmf(p, κ)
end
function MeasureTheory.logdensity(d::VonMisesFisher{M,(:c,)}, x) where {M}
    p = manifold_dimension(base_manifold(d)) + 1
    c = d.c
    κ = norm(c)
    return real(dot(c, x)) - lognorm_vmf(p, κ)
end
function MeasureTheory.logdensity(d::VonMisesFisher{M,(:F,)}, x) where {M}
    n, _ = representation_size(base_manifold(d))
    F = d.F
    return real(dot(F, x)) - logpFq((), (n//2,), (F'F) / 4)
end
function MeasureTheory.logdensity(d::VonMisesFisher{M,(:U, :D, :V)}, x) where {M}
    n, _ = representation_size(base_manifold(d))
    D = Diagonal(d.D)
    return real(dot(D * d.V', d.U' * x)) - logpFq((), (n//2,), D .^ 2 ./ 4)
end
function MeasureTheory.logdensity(d::VonMisesFisher{M,(:H, :P)}, x) where {M}
    n, _ = representation_size(base_manifold(d))
    P = d.P
    return real(dot(d.H, P, x)) - logpFq((), (n//2,), (P^2) / 4)
end

StatsBase.mode(d::VonMisesFisher{<:Any,(:μ, :κ)}) = d.μ
StatsBase.mode(d::VonMisesFisher{<:Any,(:c,)}) = (c = d.c; c / norm(c))
StatsBase.mode(d::VonMisesFisher{<:Any,(:F,)}) = (F = svd(d.F); F.U * F.Vt)
StatsBase.mode(d::VonMisesFisher{<:Any,(:U, :D, :V)}) = d.U * d.V'
StatsBase.mode(d::VonMisesFisher{<:Any,(:H, :P)}) = d.H

# ₀F₁(p//2; κ²/4) = 2ᵛ Iᵥ(κ) / κᵛ / Γ(v + 1) for v = p/2 - 1
# Note that the usual vMF constant Cₚ(κ) is defined wrt the un-normalized Hausdorff
# measure, whereas we use the normalized Hausdorff measure here.
function lognorm_vmf(p, κ)
    iszero(κ) && return log(one(κ))
    ν = p//2 - 1
    r = logbesseli(ν, κ) + ν * (logtwo - log(κ))
    return r + loggamma(oftype(r, p//2))
end

#
# Random sampling
#

function Base.rand(rng::AbstractRNG, T::Type, d::VonMisesFisher)
    p = default_point(d, T)
    return Random.rand!(rng, p, d)
end

##
## von Mises
##

# Best, D. J., and Nicholas I. Fisher. "Efficient simulation of the von Mises distribution."
# Journal of the Royal Statistical Society: Series C (Applied Statistics) 28.2 (1979): 152-157.
# doi: 10.2307/2346732
function Base.rand(rng::AbstractRNG, T::Type, d::VonMises{ℝ,(:μ, :κ)})
    κ = T(d.κ)
    tκ = 2κ
    τ = 1 + sqrt(1 + tκ^2)
    ρ = (τ - sqrt(2τ)) / tκ
    r = (1 + ρ^2) / 2ρ
    f = zero(T)
    while true
        z = cospi(rand(rng, T))
        f = T((1 + r * z) / (r + z))
        c = κ * (r - f)
        u = rand(rng, T)
        (c * (2 - c) > u || log(c / u) + 1 ≥ c) && break
    end
    θ₀ = acos(f)
    θ = (rand(rng, (θ₀, -θ₀)))
    return mod2pi(θ + T(d.μ) + π) - π
end
function Base.rand(rng::AbstractRNG, T::Type, d::VonMises{ℂ,(:μ, :κ)})
    θ = rand(rng, T, VonMises(ℝ; μ=angle(d.μ), κ=d.κ))
    return cis(θ)
end
function Base.rand(rng::AbstractRNG, T::Type, d::VonMises{ℂ,(:c,)})
    c = d.c
    θ = rand(rng, T, VonMises(ℝ; μ=angle(c), κ=abs(c)))
    return cis(θ)
end

##
## von Mises-Fisher on the sphere
##

function Random.rand!(
    rng::AbstractRNG, p::AbstractArray, d::VonMisesFisher{<:AbstractSphere,(:μ, :κ)}
)
    n = manifold_dimension(base_manifold(d)) + 1
    return _rand_vmf_sphere!(rng, p, n, d.μ, d.κ)
end
function Random.rand!(
    rng::AbstractRNG, p::AbstractArray, d::VonMisesFisher{<:AbstractSphere,(:c,)}
)
    n = manifold_dimension(base_manifold(d)) + 1
    return _rand_vmf_sphere!(rng, p, n, d.c)
end

# Andrew T.A Wood (1994) Simulation of the von mises fisher distribution,
# Communications in Statistics - Simulation and Computation, 23:1, 157-164
# doi: 10.1080/0361091940881316
function _rand_vmf_sphere!(rng, p, n, μ, κ)
    eltype(p) <: Real && isone(n) && return _rand_vmf_0sphere!(rng, p, κ * μ)
    T = real(eltype(p))
    t = _rand_normal_scale_vmf_sphere(rng, n, T(κ))
    _rand_tangent_vmf_sphere!(rng, p)
    _combine_tangent_normal_sphere!(p, t)
    _reflect_from_xaxis_to_c!(p, μ, 1)
    return p
end
function _rand_vmf_sphere!(rng, p, n, c)
    eltype(p) <: Real && isone(n) && return _rand_vmf_0sphere!(rng, p, c)
    T = real(eltype(p))
    κ = T(norm(c))
    t = _rand_normal_scale_vmf_sphere(rng, n, κ)
    _rand_tangent_vmf_sphere!(rng, p)
    _combine_tangent_normal_sphere!(p, t)
    _reflect_from_xaxis_to_c!(p, c, κ)
    return p
end
function _rand_vmf_0sphere!(rng, p, c)
    p[1] = rand(rng, Bernoulli(; logitp=2 * c[1])) ? 1 : -1
    return p
end

# in the tangent-normal parameterization
# p = t x + √(1 - t²) [0; ξ], for x the x-axis, ξ ∈ 𝕊ⁿ⁻², and t ∈ [-1, 1],
# then ξ follows the normalized Hausdorff measure, and p(t) ∝ (1 - t^2)^((n-3)/2).
# so we draw t and ξ, compose p, then then transform it using the reflection that
# takes the x to μ.
# Method due to Wood, 1994. Adapted also for complex and quaternionic spheres.
function _rand_normal_scale_vmf_sphere(rng, n, κ)
    T = eltype(κ)
    twoκ = 2κ
    if n == 3
        # 2t+1 follows an exponential distribution truncated to [0, 1]
        # so we use inverse transform sampling
        u = rand(rng, T)
        return 1 + log(u + exp(-twoκ) * (1 - u)) / κ
    end
    m = T((n - 1)//2)
    a = twoκ / (n - 1)
    b = sqrt(a^2 + 1) - a
    x = (1 - b) / (1 + b)
    c = κ * x + (n - 1) * log1p(-x^2)
    βdist = Beta(m, m)

    z = rand(rng, T, βdist)
    t = (1 - (1 + b) * z) / (1 - (1 - b) * z)
    while κ * t + (n - 1) * log1p(-x * t) - c < log(rand(rng))
        z = rand(rng, T, βdist)
        t = (1 - (1 + b) * z) / (1 - (1 - b) * z)
    end
    return t
end
function _rand_tangent_vmf_sphere!(rng, p)
    randn!(rng, p)
    p[1] -= real(p[1])
    rdiv!(p, norm(p))
    return p
end
function _combine_tangent_normal_sphere!(p, t)
    rmul!(p, sqrt(1 - t^2))
    p[1] += t
    return p
end
# in-place apply Householder reflection p ↦ p - q 2𝕽⟨q,p⟩/‖q‖², for q=e₁-c/‖c‖
function _reflect_from_xaxis_to_c!(p, c, cnorm=norm(c))
    num = real(p[1]) - real(dot(c, p)) / cnorm
    den = cnorm - real(c[1])
    α = num / den
    p .+= c .* α
    p[1] -= α * cnorm
    return p
end

##
## von Mises-Fisher on the Stiefel manifold
##

function Random.rand!(
    rng::AbstractRNG, p::AbstractArray, d::VonMisesFisher{<:Stiefel{n,k},(:U, :D, :V)}
) where {n,k}
    dim𝔽 = real_dimension(number_system(base_manifold(d)))
    return _rand_vmf_stiefel!(rng, p, dim𝔽, n, k, d.U, d.D, d.V)
end
function Random.rand!(
    rng::AbstractRNG, p::AbstractArray, d::VonMisesFisher{<:Stiefel{n,k},(:F,)}
) where {n,k}
    U, D, V = svd(d.F)
    dim𝔽 = real_dimension(number_system(base_manifold(d)))
    return _rand_vmf_stiefel!(rng, p, dim𝔽, n, k, U, D, V)
end
function Random.rand!(
    rng::AbstractRNG, p::AbstractArray, d::VonMisesFisher{<:Stiefel{n,k},(:H, :P)}
) where {n,k}
    D, V = eigen(Hermitian(d.P))
    U = d.H * V
    dim𝔽 = real_dimension(number_system(base_manifold(d)))
    return _rand_vmf_stiefel!(rng, p, dim𝔽, n, k, U, D, V)
end

# Peter Hoff. Simulation of the Matrix Bingham—von Mises—Fisher Distribution,
# With Applications to Multivariate and Relational Data.
# Journal of Computational and Graphical Statistics. 18(2). 2009.
function _rand_vmf_stiefel!(rng, p, dim𝔽, n, k, U, D, V)
    if isone(k)
        _rand_vmf_sphere!(rng, vec(p), dim𝔽 * n, vec(U), D[1])
        rmul!(p, V[1]')
        return p
    end
    T = real(eltype(p))
    Z = fill!(similar(p), zero(T))
    Z₁ = @view Z[:, 1]
    U₁ = @view U[:, 1]
    y = similar(Z₁)
    z = similar(U₁)
    while true
        _rand_vmf_sphere!(rng, Z₁, dim𝔽 * n, U₁, D[1])
        lcrit = zero(T)
        for j in 2:k
            s = n - j + 1
            r = dim𝔽 * s
            @views begin
                N = _nullbasis(Z[:, 1:(j - 1)])
                Uⱼ = U[:, j]
                Zⱼ = Z[:, j]
                zⱼ = z[1:s]
                yⱼ = y[1:s]
            end
            Dⱼ = D[j]
            if Dⱼ > 0
                mul!(zⱼ, N', Uⱼ, Dⱼ, false)
                _rand_vmf_sphere!(rng, yⱼ, r, zⱼ)
                mul!(Zⱼ, N, yⱼ)
                nzⱼ = norm(zⱼ)
                ν = r//2 - 1
                lcrit += T(
                    logbesseli(ν, nzⱼ) - logbesseli(ν, Dⱼ) + ν * (log(Dⱼ) - log(nzⱼ))
                )
            else  # sample from uniform distribution, lcrit contribution is zero
                randn!(rng, yⱼ)
                mul!(Zⱼ, N, yⱼ, inv(norm(yⱼ)), false)
            end
        end
        log(rand(rng)) < lcrit && break
    end
    mul!(p, Z, V')
    return p
end

# basis N of null space of A, such that N'A=0
# compute basis of null space of matrix A
function _nullbasis(A)
    F = qr(A)
    rank = size(F.R, 1)
    return F.Q[:, (rank + 1):end]
end


text:
module Utils

using HTTP: Request, payload, queryparams
using HTTP.URIs: URI, splitpath
using JSON3

export json_payload, query_parameters, uri_parameter

# json_payload

"""
    json_payload(req::Request)::JSON3.Object

Retrieve the json payload from a request.
"""
json_payload(req::Request) = req |> payload |> IOBuffer |> JSON3.read

# query_parameters

"""
    query_parameters(req::Request)::Dict{String, String}

Retrieve the parameters encoded in the querystring from the request.
"""
query_parameters(req::Request) = req.target |> URI |> queryparams

# uri_parameter

"""
    uri_parameter(req::Request, segment_idx::Int)::String

Retrieve a uri parameter from a request.
"""
uri_parameter(req::Request, segment_idx::Int) = splitpath(req.target) |> x -> x[segment_idx]

end


text:
immutable ScreenSample
    "name of original file"
    filename::String
    "total number of reads"
    reads::Int
    "Number of called hits"
    hits::Int
    "Number of high quality reads that didn't match"
    no_match::Int
    "Number of reads with too many errors in their constant region"
    constant_error::Int
    "Number of reads with low quality sgRNA regions"
    low_quality::Int
    "Number of reads with a priming site starting at each position of the read"
    alt_staggers::Vector{Int}
    counts::Dict{String, Int}
    hit_qualities::Matrix{Int}
    all_qualities::Matrix{Int}

    function ScreenSample(filename, read_length, max_phred_quality)
        new(filename, 0, 0, 0, 0, 0, zeros(Int, read_length),
            Dict{String, Int}(),
            zeros(Int, max_phred_quality, read_length),
            zeros(Int, max_phred_quality, read_length))
    end
end

function Base.show(io::IO, s::ScreenSample)
    println(io, typeof(s))
    println(io, "\tFilename: $(s.filename)")
    println(io, "\tHits: ", format(s.hits, commas=true))
    println(io, "\tReads: ", format(s.reads, commas=true))
    println(io, "\t% reads matched: $(round(s.hits/s.reads*100, 4))%")
    println(io, "\t% unmatchable reads: $(round(s.no_match/(s.hits+s.no_match)*100, 4))%")
    println(io, "\t% reads had errors in the constant region: $(round(s.constant_error/s.reads*100, 4))%")
    println(io, "\t% reads had low quality sgrna regions: $(round(s.low_quality/s.reads*100, 4))%")
end


text:
# Code snippet handling
struct CodeSnippet
    code_lines::Vector{String}
    lang::Any
end
function init_code(lang)
    return CodeSnippet(Vector{String}(undef, 0), lang)
end
function push_code!(code, str; ind_lvl = 1, ind_str = "    ")
    # Indent only non-empty lines.
    indentation = repeat(ind_str, ind_lvl)
    indented_string = isempty(str) ? "" : indentation * str
    return push!(code.code_lines, indented_string)
end
function push_code_verbatim_string!(code, str)
    for line in split(str, "\n")
        push!(code.code_lines, line)
    end
end
function push_comment!(code, str; ind_lvl = 1, ind_str = "    ")
    # Convert empty comments to empty lines.
    return push_code!(
        code,
        isempty(str) ? "" : comment(code.lang, str),
        ind_lvl = ind_lvl,
        ind_str = ind_str,
    )
end
function to_string(code)
    return join(code.code_lines, "\n")
end


text:
@testset "largest_component" begin
    A = load_matrix_network("dfs_example")
    (Acc,p) = largest_component(A)
    @test size(Acc,1) == 5

    (Acc,p) = largest_component(A,true)
    @test size(Acc,1) == 6

    A = load_matrix_network("cores_example")
    (Acc1,p1) = largest_component(A)
    (Acc2,p2) = largest_component(A,true)
    @test Acc1 == Acc2
end