text:
"""
Replacements pairs for all languages.
Source: <https://github.com/sindresorhus/transliterate/blob/master/replacements.js>
"""
const replacements = (
    # General (named Latin in the source)
    "la" => (
        "À" => "A",
        "Á" => "A",
        "Â" => "A",
        "Ã" => "A",
        "Ä" => "Ae",
        "Å" => "A",
        "Æ" => "AE",
        "Ç" => "C",
        "È" => "E",
        "É" => "E",
        "Ê" => "E",
        "Ë" => "E",
        "Ì" => "I",
        "Í" => "I",
        "Î" => "I",
        "Ï" => "I",
        "Ð" => "D",
        "Ñ" => "N",
        "Ò" => "O",
        "Ó" => "O",
        "Ô" => "O",
        "Õ" => "O",
        "Ö" => "Oe",
        "Ő" => "O",
        "Ø" => "O",
        "Ù" => "U",
        "Ú" => "U",
        "Û" => "U",
        "Ü" => "Ue",
        "Ű" => "U",
        "Ý" => "Y",
        "Þ" => "TH",
        "ß" => "ss",
        "à" => "a",
        "á" => "a",
        "â" => "a",
        "ã" => "a",
        "ä" => "ae",
        "å" => "a",
        "æ" => "ae",
        "ç" => "c",
        "è" => "e",
        "é" => "e",
        "ê" => "e",
        "ë" => "e",
        "ì" => "i",
        "í" => "i",
        "î" => "i",
        "ï" => "i",
        "ð" => "d",
        "ñ" => "n",
        "ò" => "o",
        "ó" => "o",
        "ô" => "o",
        "õ" => "o",
        "ö" => "oe",
        "ő" => "o",
        "ø" => "o",
        "ù" => "u",
        "ú" => "u",
        "û" => "u",
        "ü" => "ue",
        "ű" => "u",
        "ý" => "y",
        "þ" => "th",
        "ÿ" => "y",
        "ẞ" => "SS",
    ),

    # German
    "de" => (
        "ß" => "ss",
        "ä" => "ae",
        "Ä" => "Ae",
        "ö" => "oe",
        "Ö" => "Oe",
        "ü" => "ue",
        "Ü" => "Ue",
    ),

    # Vietnamese
    "vi" => (
        "à" => "a",
        "À" => "A",
        "á" => "a",
        "Á" => "A",
        "â" => "a",
        "Â" => "A",
        "ã" => "a",
        "Ã" => "A",
        "è" => "e",
        "È" => "E",
        "é" => "e",
        "É" => "E",
        "ê" => "e",
        "Ê" => "E",
        "ì" => "i",
        "Ì" => "I",
        "í" => "i",
        "Í" => "I",
        "ò" => "o",
        "Ò" => "O",
        "ó" => "o",
        "Ó" => "O",
        "ô" => "o",
        "Ô" => "O",
        "õ" => "o",
        "Õ" => "O",
        "ù" => "u",
        "Ù" => "U",
        "ú" => "u",
        "Ú" => "U",
        "ý" => "y",
        "Ý" => "Y",
        "ă" => "a",
        "Ă" => "A",
        "Đ" => "D",
        "đ" => "d",
        "ĩ" => "i",
        "Ĩ" => "I",
        "ũ" => "u",
        "Ũ" => "U",
        "ơ" => "o",
        "Ơ" => "O",
        "ư" => "u",
        "Ư" => "U",
        "ạ" => "a",
        "Ạ" => "A",
        "ả" => "a",
        "Ả" => "A",
        "ấ" => "a",
        "Ấ" => "A",
        "ầ" => "a",
        "Ầ" => "A",
        "ẩ" => "a",
        "Ẩ" => "A",
        "ẫ" => "a",
        "Ẫ" => "A",
        "ậ" => "a",
        "Ậ" => "A",
        "ắ" => "a",
        "Ắ" => "A",
        "ằ" => "a",
        "Ằ" => "A",
        "ẳ" => "a",
        "Ẳ" => "A",
        "ẵ" => "a",
        "Ẵ" => "A",
        "ặ" => "a",
        "Ặ" => "A",
        "ẹ" => "e",
        "Ẹ" => "E",
        "ẻ" => "e",
        "Ẻ" => "E",
        "ẽ" => "e",
        "Ẽ" => "E",
        "ế" => "e",
        "Ế" => "E",
        "ề" => "e",
        "Ề" => "E",
        "ể" => "e",
        "Ể" => "E",
        "ễ" => "e",
        "Ễ" => "E",
        "ệ" => "e",
        "Ệ" => "E",
        "ỉ" => "i",
        "Ỉ" => "I",
        "ị" => "i",
        "Ị" => "I",
        "ọ" => "o",
        "Ọ" => "O",
        "ỏ" => "o",
        "Ỏ" => "O",
        "ố" => "o",
        "Ố" => "O",
        "ồ" => "o",
        "Ồ" => "O",
        "ổ" => "o",
        "Ổ" => "O",
        "ỗ" => "o",
        "Ỗ" => "O",
        "ộ" => "o",
        "Ộ" => "O",
        "ớ" => "o",
        "Ớ" => "O",
        "ờ" => "o",
        "Ờ" => "O",
        "ở" => "o",
        "Ở" => "O",
        "ỡ" => "o",
        "Ỡ" => "O",
        "ợ" => "o",
        "Ợ" => "O",
        "ụ" => "u",
        "Ụ" => "U",
        "ủ" => "u",
        "Ủ" => "U",
        "ứ" => "u",
        "Ứ" => "U",
        "ừ" => "u",
        "Ừ" => "U",
        "ử" => "u",
        "Ử" => "U",
        "ữ" => "u",
        "Ữ" => "U",
        "ự" => "u",
        "Ự" => "U",
        "ỳ" => "y",
        "Ỳ" => "Y",
        "ỵ" => "y",
        "Ỵ" => "Y",
        "ỷ" => "y",
        "Ỷ" => "Y",
        "ỹ" => "y",
        "Ỹ" => "Y",
    ),

    # Arabic
    "ar" => (
        "ء" => "e",
        "آ" => "a",
        "أ" => "a",
        "ؤ" => "w",
        "إ" => "i",
        "ئ" => "y",
        "ا" => "a",
        "ب" => "b",
        "ة" => "t",
        "ت" => "t",
        "ث" => "th",
        "ج" => "j",
        "ح" => "h",
        "خ" => "kh",
        "د" => "d",
        "ذ" => "dh",
        "ر" => "r",
        "ز" => "z",
        "س" => "s",
        "ش" => "sh",
        "ص" => "s",
        "ض" => "d",
        "ط" => "t",
        "ظ" => "z",
        "ع" => "e",
        "غ" => "gh",
        "ـ" => "_",
        "ف" => "f",
        "ق" => "q",
        "ك" => "k",
        "ل" => "l",
        "م" => "m",
        "ن" => "n",
        "ه" => "h",
        "و" => "w",
        "ى" => "a",
        "ي" => "y",
        "َ" => "a",
        "ُ" => "u",
        "ِ" => "i",
        "٠" => "0",
        "١" => "1",
        "٢" => "2",
        "٣" => "3",
        "٤" => "4",
        "٥" => "5",
        "٦" => "6",
        "٧" => "7",
        "٨" => "8",
        "٩" => "9",
    ),

    # Persian / Farsi
    "fa" => (
        "چ" => "ch",
        "ک" => "k",
        "گ" => "g",
        "پ" => "p",
        "ژ" => "zh",
        "ی" => "y",
        "۰" => "0",
        "۱" => "1",
        "۲" => "2",
        "۳" => "3",
        "۴" => "4",
        "۵" => "5",
        "۶" => "6",
        "۷" => "7",
        "۸" => "8",
        "۹" => "9",
    ),

    # Pashto
    "ps" => (
        "ټ" => "p",
        "ځ" => "z",
        "څ" => "c",
        "ډ" => "d",
        "ﺫ" => "d",
        "ﺭ" => "r",
        "ړ" => "r",
        "ﺯ" => "z",
        "ږ" => "g",
        "ښ" => "x",
        "ګ" => "g",
        "ڼ" => "n",
        "ۀ" => "e",
        "ې" => "e",
        "ۍ" => "ai",
    ),

    # Urdu
    "ur" => (
        "ٹ" => "t",
        "ڈ" => "d",
        "ڑ" => "r",
        "ں" => "n",
        "ہ" => "h",
        "ھ" => "h",
        "ے" => "e",
    ),

    # Russian
    "ru" => (
        "А" => "A",
        "а" => "a",
        "Б" => "B",
        "б" => "b",
        "В" => "V",
        "в" => "v",
        "Г" => "G",
        "г" => "g",
        "Д" => "D",
        "д" => "d",
        "Е" => "E",
        "е" => "e",
        "Ж" => "Zh",
        "ж" => "zh",
        "З" => "Z",
        "з" => "z",
        "И" => "I",
        "и" => "i",
        "Й" => "J",
        "й" => "j",
        "К" => "K",
        "к" => "k",
        "Л" => "L",
        "л" => "l",
        "М" => "M",
        "м" => "m",
        "Н" => "N",
        "н" => "n",
        "О" => "O",
        "о" => "o",
        "П" => "P",
        "п" => "p",
        "Р" => "R",
        "р" => "r",
        "С" => "S",
        "с" => "s",
        "Т" => "T",
        "т" => "t",
        "У" => "U",
        "у" => "u",
        "Ф" => "F",
        "ф" => "f",
        "Х" => "H",
        "х" => "h",
        "Ц" => "Cz",
        "ц" => "cz",
        "Ч" => "Ch",
        "ч" => "ch",
        "Ш" => "Sh",
        "ш" => "sh",
        "Щ" => "Shch",
        "щ" => "shch",
        "Ъ" => "",
        "ъ" => "",
        "Ы" => "Y",
        "ы" => "y",
        "Ь" => "",
        "ь" => "",
        "Э" => "E",
        "э" => "e",
        "Ю" => "Yu",
        "ю" => "yu",
        "Я" => "Ya",
        "я" => "ya",
        "Ё" => "Yo",
        "ё" => "yo",
    ),

    # Romanian
    "ro" => (
        "ă" => "a",
        "Ă" => "A",
        "ș" => "s",
        "Ș" => "S",
        "ț" => "t",
        "Ț" => "T",
        "ţ" => "t",
        "Ţ" => "T",
    ),

    # Turkish
    "tr" => (
        "ş" => "s",
        "Ş" => "S",
        "ç" => "c",
        "Ç" => "C",
        "ğ" => "g",
        "Ğ" => "G",
        "ı" => "i",
        "İ" => "I",
    ),

    # Armenian
    "hy" => (
        "ա" => "a",
        "բ" => "b",
        "գ" => "g",
        "դ" => "d",
        "ե" => "ye",
        "զ" => "z",
        "է" => "e",
        "ը" => "u",
        "թ" => "t",
        "ժ" => "zh",
        "ի" => "i",
        "լ" => "l",
        "խ" => "kh",
        "ծ" => "ts",
        "կ" => "k",
        "հ" => "h",
        "ձ" => "dz",
        "ղ" => "r",
        "ճ" => "j",
        "մ" => "m",
        "յ" => "j",
        "ն" => "n",
        "շ" => "sh",
        "ո" => "vo",
        "չ" => "ch",
        "պ" => "p",
        "ջ" => "j",
        "ռ" => "r",
        "ս" => "s",
        "վ" => "v",
        "տ" => "t",
        "ր" => "re",
        "ց" => "ts",
        "ու" => "u",
        "ւ" => "v",
        "փ" => "p",
        "ք" => "q",
        "օ" => "o",
        "ֆ" => "f",
        "և" => "yev",
    ),

    # Georgian
    "ka" => (
        "ა" => "a",
        "ბ" => "b",
        "გ" => "g",
        "დ" => "d",
        "ე" => "e",
        "ვ" => "v",
        "ზ" => "z",
        "თ" => "t",
        "ი" => "i",
        "კ" => "k",
        "ლ" => "l",
        "მ" => "m",
        "ნ" => "n",
        "ო" => "o",
        "პ" => "p",
        "ჟ" => "zh",
        "რ" => "r",
        "ს" => "s",
        "ტ" => "t",
        "უ" => "u",
        "ფ" => "ph",
        "ქ" => "q",
        "ღ" => "gh",
        "ყ" => "k",
        "შ" => "sh",
        "ჩ" => "ch",
        "ც" => "ts",
        "ძ" => "dz",
        "წ" => "ts",
        "ჭ" => "tch",
        "ხ" => "kh",
        "ჯ" => "j",
        "ჰ" => "h",
    ),

    # Czech
    "cz" => (
        "č" => "c",
        "ď" => "d",
        "ě" => "e",
        "ň" => "n",
        "ř" => "r",
        "š" => "s",
        "ť" => "t",
        "ů" => "u",
        "ž" => "z",
        "Č" => "C",
        "Ď" => "D",
        "Ě" => "E",
        "Ň" => "N",
        "Ř" => "R",
        "Š" => "S",
        "Ť" => "T",
        "Ů" => "U",
        "Ž" => "Z",
    ),

    # Dhivehi
    "dv" => (
        "ހ" => "h",
        "ށ" => "sh",
        "ނ" => "n",
        "ރ" => "r",
        "ބ" => "b",
        "ޅ" => "lh",
        "ކ" => "k",
        "އ" => "a",
        "ވ" => "v",
        "މ" => "m",
        "ފ" => "f",
        "ދ" => "dh",
        "ތ" => "th",
        "ލ" => "l",
        "ގ" => "g",
        "ޏ" => "gn",
        "ސ" => "s",
        "ޑ" => "d",
        "ޒ" => "z",
        "ޓ" => "t",
        "ޔ" => "y",
        "ޕ" => "p",
        "ޖ" => "j",
        "ޗ" => "ch",
        "ޘ" => "tt",
        "ޙ" => "hh",
        "ޚ" => "kh",
        "ޛ" => "th",
        "ޜ" => "z",
        "ޝ" => "sh",
        "ޞ" => "s",
        "ޟ" => "d",
        "ޠ" => "t",
        "ޡ" => "z",
        "ޢ" => "a",
        "ޣ" => "gh",
        "ޤ" => "q",
        "ޥ" => "w",
        "ަ" => "a",
        "ާ" => "aa",
        "ި" => "i",
        "ީ" => "ee",
        "ު" => "u",
        "ޫ" => "oo",
        "ެ" => "e",
        "ޭ" => "ey",
        "ޮ" => "o",
        "ޯ" => "oa",
        "ް" => "",
    ),

    # Greek
    "el" => (
        "α" => "a",
        "β" => "v",
        "γ" => "g",
        "δ" => "d",
        "ε" => "e",
        "ζ" => "z",
        "η" => "i",
        "θ" => "th",
        "ι" => "i",
        "κ" => "k",
        "λ" => "l",
        "μ" => "m",
        "ν" => "n",
        "ξ" => "ks",
        "ο" => "o",
        "π" => "p",
        "ρ" => "r",
        "σ" => "s",
        "τ" => "t",
        "υ" => "y",
        "φ" => "f",
        "χ" => "x",
        "ψ" => "ps",
        "ω" => "o",
        "ά" => "a",
        "έ" => "e",
        "ί" => "i",
        "ό" => "o",
        "ύ" => "y",
        "ή" => "i",
        "ώ" => "o",
        "ς" => "s",
        "ϊ" => "i",
        "ΰ" => "y",
        "ϋ" => "y",
        "ΐ" => "i",
        "Α" => "A",
        "Β" => "B",
        "Γ" => "G",
        "Δ" => "D",
        "Ε" => "E",
        "Ζ" => "Z",
        "Η" => "I",
        "Θ" => "TH",
        "Ι" => "I",
        "Κ" => "K",
        "Λ" => "L",
        "Μ" => "M",
        "Ν" => "N",
        "Ξ" => "KS",
        "Ο" => "O",
        "Π" => "P",
        "Ρ" => "R",
        "Σ" => "S",
        "Τ" => "T",
        "Υ" => "Y",
        "Φ" => "F",
        "Χ" => "X",
        "Ψ" => "PS",
        "Ω" => "O",
        "Ά" => "A",
        "Έ" => "E",
        "Ί" => "I",
        "Ό" => "O",
        "Ύ" => "Y",
        "Ή" => "I",
        "Ώ" => "O",
        "Ϊ" => "I",
        "Ϋ" => "Y",
    ),

    # Latvian
    "lv" => (
        "ā" => "a",
        "ē" => "e",
        "ģ" => "g",
        "ī" => "i",
        "ķ" => "k",
        "ļ" => "l",
        "ņ" => "n",
        "ū" => "u",
        "Ā" => "A",
        "Ē" => "E",
        "Ģ" => "G",
        "Ī" => "I",
        "Ķ" => "K",
        "Ļ" => "L",
        "Ņ" => "N",
        "Ū" => "U",
        "č" => "c",
        "š" => "s",
        "ž" => "z",
        "Č" => "C",
        "Š" => "S",
        "Ž" => "Z",
    ),

    # Lithuanian
    "lt" => (
        "ą" => "a",
        "č" => "c",
        "ę" => "e",
        "ė" => "e",
        "į" => "i",
        "š" => "s",
        "ų" => "u",
        "ū" => "u",
        "ž" => "z",
        "Ą" => "A",
        "Č" => "C",
        "Ę" => "E",
        "Ė" => "E",
        "Į" => "I",
        "Š" => "S",
        "Ų" => "U",
        "Ū" => "U",
    ),

    # Macedonian
    "mk" => (
        "Ќ" => "Kj",
        "ќ" => "kj",
        "Љ" => "Lj",
        "љ" => "lj",
        "Њ" => "Nj",
        "њ" => "nj",
        "Тс" => "Ts",
        "тс" => "ts",
    ),

    # Ukrainian
    "uk" => (
        "Є" => "Ye",
        "І" => "I",
        "Ї" => "Yi",
        "Ґ" => "G",
        "є" => "ye",
        "і" => "i",
        "ї" => "yi",
        "ґ" => "g"
    ),

    # Polish
    "pl" => (
        "ą" => "a",
        "ć" => "c",
        "ę" => "e",
        "ł" => "l",
        "ń" => "n",
        "ś" => "s",
        "ź" => "z",
        "ż" => "z",
        "Ą" => "A",
        "Ć" => "C",
        "Ę" => "E",
        "Ł" => "L",
        "Ń" => "N",
        "Ś" => "S",
        "Ź" => "Z",
        "Ż" => "Z",
    ),

    # Hungarian
    # Conflicts with German and Latin.
    "hu" => (
        "ä" => "a",
        "Ä" => "A",
        "ö" => "o",
        "Ö" => "O",
        "ü" => "u",
        "Ü" => "U",
        "ű" => "u",
        "Ű" => "U",
    ),

    # Serbian
    # Conflicts with Vietnamese.
    "sr" => (
        "ђ" => "dj",
        "Ђ" => "Dj",
        "љ" => "lj",
        "Љ" => "Lj",
        "ј" => "j",
        "Ј" => "Ј",
        "Њ" => "Nj",
        "њ" => "nj",
        "ћ" => "c",
        "Ћ" => "C",
        "ц" => "c",
        "Ц" => "C",
        "џ" => "dz",
        "Џ" => "Dz",
        "đ" => "dj",
        "Đ" => "Dj"
    ),

    # Slovak
    # Conflicts with German and Latin.
    "sk" => (
        "ä" => "a",
        "Ä" => "A",
        "ľ" => "l",
        "ĺ" => "l",
        "ŕ" => "r",
        "Ľ" => "L",
        "Ĺ" => "L",
        "Ŕ" => "R",
    ),

    # Swedish
    # Conflicts with German and Latin.
    "sw" => (
        "å" => "o",
        "Å" => "o",
        "ä" => "a",
        "Ä" => "A",
        "ë" => "e",
        "Ë" => "E",
        "ö" => "o",
        "Ö" => "O",
    ),

    # Danish
    "da" => (
        "Æ" => "Ae",
        "Ø" => "Oe",
        "Å" => "Aa",
        "æ" => "ae",
        "ø" => "oe",
        "å" => "aa"
    )
)


text:
export
   AbstractSzegedy,
   Szegedy,
   sqrtstochastic

"""
    AbstractSzegedy

Type representing the abstract Szegedy model. Description of the default
parameter can be found in https://arxiv.org/abs/1611.02238, where two oracle
operator case is chosen.  Default representation of `AbstractSzegedy` is
`Szegedy`.
"""
abstract type AbstractSzegedy <: QWModelDiscr end

"""
    Szegedy{G, T}(graph, sqrtstochastic)

Default representation of `AbstractSzegedy`. Parameter `sqrtstochastic` needs to
be an element-wise square root of stochastic matrix.

    Szegedy(graph[, stochastic, checkstochastic])

Constructors of `AbstractSzegedy`. Parameter `stochastic` needs to be a
stochastic matrix.  Flag `checkstochastic` decides about checking the stochastic
properties.

Matrix `stochastic` defaults to the uniform walk operator, and `checkstochastic`
deafults to `false` in case of default `stochastic`. If matrix `stochastic` is
provided by the user, the default value of `stochastic` is `true`.
"""
struct Szegedy{G<:AbstractGraph, T<:Number} <: AbstractSzegedy
    graph::G
    sqrtstochastic::SparseMatrixCSC{T,Int}

    function Szegedy{G, T}(graph::G,
                           sqrtstochastic::SparseMatrixCSC{T}) where {G<:AbstractGraph, T<:Number}
        new{G, T}(graph, sqrtstochastic)
    end
end

function Szegedy(graph::G,
                 stochastic::SparseMatrixCSC{T},
                 checkstochastic::Bool=true) where {G<:AbstractGraph, T<:Number}
    if checkstochastic
        graphstochasticcheck(graph, stochastic)
    end
    Szegedy{G, T}(graph, sqrt.(stochastic))
end

function Szegedy(graph::AbstractGraph)
    Szegedy(graph, default_stochastic(graph), false)
end

"""
    sqrtstochastic(szegedy)

Returns the `sqrtstochastic` element of `szegedy`.
"""
sqrtstochastic(szegedy::AbstractSzegedy) = szegedy.sqrtstochastic

"""
    QWSearch(szegedy, marked[, penalty])

Creates `QWSearch` according to `AbstractSzegedy` model. By default parameter
`penalty` is set to 0. Evolution operators are constructed according to the
definition from https://arxiv.org/abs/1611.02238.

    QWSearch(qws[; marked, penalty]; check)

Update quantum walk search to new subset of marked elements and new penalty. By
default `marked` and `penalty` are the same as in `qws`.
"""
function QWSearch(szegedy::AbstractSzegedy,
                  marked::Vector{Int},
                  penalty::Real=0)
   walk_operators = szegedy_walk_operators(szegedy)
   oracle_operators = szegedy_oracle_operators(szegedy, marked)
   parameters = Dict{Symbol,Any}()
   parameters[:operators] = walk_operators .* oracle_operators

   QWSearch(szegedy, parameters, marked, penalty)
end

function QWSearch(qws::QWSearch{<:Szegedy};
                  marked::Vector{Int}=qws.marked,
                  penalty::Real=qws.penalty)
   oldmarked = qws.marked
   local corr_oracles
   if Set(marked) != Set(oldmarked)
    corr_oracles = szegedy_oracle_operators(model(qws), oldmarked) .*
                   szegedy_oracle_operators(model(qws), marked)
   else
      corr_oracles = fill(I, length(parameters(qws)[:operators]))#speye.(parameters(qws)[:operators])
   end
   QWSearch(model(qws),
            Dict(:operators => parameters(qws)[:operators].*corr_oracles),
            marked,
            penalty)
end

"""
    QWEvolution(szegedy)

Create `QWEvolution` according to `AbstractSzegedy` model. By default, the
constructed operator is of type `SparseMatrixCSC`.
"""
function QWEvolution(szegedy::AbstractSzegedy)
   parameters = Dict{Symbol,Any}(:operators => szegedy_walk_operators(szegedy))

   QWEvolution(szegedy, parameters)
end

"""
    check_qwmodel(szegedy, parameters)

Private function for checking the existance of `:operators`, its type, and the
dimensionality of its elements.
"""
function check_qwmodel(szegedy::AbstractSzegedy,
                       parameters::Dict{Symbol})
   @assert :operators in keys(parameters) "Parameters should contain key operators"
   @assert length(parameters[:operators]) == 2 "Only two operators allowed"
   @assert all(typeof(i) <: AbstractMatrix{<:Number} for i=parameters[:operators]) "Parameters should be a list of SparseMatrixCSC{<:Number}"
   order = nv(szegedy.graph)
   @assert all(size(i) == (order, order).^2 for i=parameters[:operators]) "Operators sizes mismatch"
end

"""
    check_qwdynamics(QWSearch, szegedy, marked, parameters)

Check whetver combination of `szegedy`, `marked` and `parameters` produces valid
`QWSearch` object. It checks where `parameters` consists of key `:operators` with
corresponding value being list of `AbstractMatrix`. Furthermore operators
needs to be square of size equals to square of `graph(szegedy).` order.
"""
function check_qwdynamics(::Type{QWSearch},
                          szegedy::AbstractSzegedy,
                          parameters::Dict{Symbol},
                          marked::Vector{Int})
   check_qwmodel(szegedy, parameters)
end

"""
    check_qwdynamics(QWEvolution, szegedy, parameters)

Check whetver combination of `szegedy`, and `parameters` produces a
valid `QWEvolution` object. It checks where `parameters` consists of key
`:operators` with corresponding value being a list of `AbstractMatrix` objects.
Furthermore operators need to be square of size equals to square of the order of
`graph(szegedy)`.
"""
function check_qwdynamics(::Type{QWEvolution},
                          szegedy::AbstractSzegedy,
                          parameters::Dict{Symbol})
   check_qwmodel(szegedy, parameters)
end

include("szegedy_stochastic.jl")
include("szegedy_operators.jl")
include("szegedy_evolution.jl")

ext: jl


text:
using ApproxFun, LowRankApprox, SingularIntegralEquations, Test

@testset "LowRankMatrix" begin
    L = LowRankApprox._LowRankMatrix(1 ./(1:10),1 ./(1:10.0))
    A = rand(10,10)

    @test isa(L⊕A,LowRankMatrix)

    @test rank(L+L) == 2rank(L⊕L)
    # @test rank(L-L) == 2rank(L⊖L)
    @test rank(2L+1+Matrix(L)) ≥ rank(L⊕Matrix(L))
end


text:
# HSA runtime and ROCm External Libraries
## copied from CUDAdrv/src/CUDAdrv.jl
if !parse(Bool, get(ENV, "JULIA_AMDGPU_DISABLE_ARTIFACTS", "false"))
    using hsa_rocr_jll
    using hsakmt_roct_jll
end

const ext = joinpath(@__DIR__, "ext.jl")
if !isfile(ext)
    @warn "Didn't find $ext, please build AMDGPU.jl"
    const hsa_configured = false
    const ext_libs_configured = false
else
    include(ext)
end

# Default (non-functional) values for critical variables,
# making it possible to load the package at all times.
if !hsa_configured
    const libhsaruntime_version = v"0.0"
    const libhsaruntime_vendor = "none"
    const libhsaruntime_path = nothing
end
if !ext_libs_configured
    const librocblas = nothing
    const librocsparse = nothing
    const librocalution = nothing
    const librocfft = nothing
    const librocrand = nothing
    const libmiopen = nothing
    const libhip = nothing
end

# ROCm-Device-Libs
device_libs_deps = joinpath(@__DIR__, "deps.jl")
isfile(device_libs_deps) && include(device_libs_deps)
const device_libs_path = joinpath(@__DIR__, "usr", "lib")

const configured = hsa_configured


text:
"""
    save(filename::AbstractString, f, args...; kwargs...)
    save(io::IO, f, args...; kwargs...)

Serialize the result of `modelproto(f, args...; kwargs...)` to a file with path `filename` or to `io`.

See [`modelproto`](@ref) for description of arguments.
"""
save(filename::AbstractString, f, args...; modelname=filename, kwargs...) = save(filename, modelproto(f, args...; modelname=modelname, kwargs...))
save(io::IO, f, args...; kwargs...) = save(io, modelproto(f, args...; kwargs...))

"""
    save(filename::AbstractString, mp::ONNX.ModelProto)
    save(io::IO, mp::ONNX.ModelProto)

Serialize the given [`ONNX.ModelProto`](@ref) to a file with path `filename` or to `io`.
"""
save(filename::AbstractString, mp::ONNX.ModelProto) = open(io -> save(io, mp), filename, "w")
save(io::IO, mp::ONNX.ModelProto) = ONNX.writeproto(io, mp)

"""
    modelproto(f; namestrat=name_runningnr(), posthook=validate, kwargs...)
    modelproto(f, inshapes::Tuple...; namestrat = name_runningnr(), posthook=validate, kwargs...)
    modelproto(f, indata::Pair{String, <:Any}...; modelname="model", namestrat = name_runningnr(), posthook=validate, kwargs...)

Return an [`ONNX.ModelProto`](@ref) from `f`.

Argument `inshapes` are size tuples representing the shape of each input. An attempt to infer sizes will be made if not provided.
Argument `indata` are pairs mapping names to size tuples. Names will be created automatically if not provided.

Argument `modelname` is a string which will be used as the name of the model. Must be non-empty to be valid ONNX.

Argument `namestrat` determines how nodes in the graph shall be named. Other keyword arguments are passed to the `ModelProto`.

Argument `posthook` will be called with the created `ONNX.ModelProto` as argument before returning it.

Other keyword arguments will be passed to `ONNX.ModelProto`.
"""
modelproto(f; kwargs...) = modelproto(f, infer_inshapes(f)...; kwargs...)
modelproto(f, inshapes::Union{Tuple, Missing}...;kwargs...) = modelproto(f, ("data_" .* string.(0:length(inshapes)-1) .=> inshapes)...; kwargs...)
function modelproto(f, indata::Pair{String, <:Any}...; modelname="model", namestrat = default_namestrat(f), posthook=validate, kwargs...)
    mp = modelproto(;kwargs...)
    mp.graph = graphproto(f, indata...; namestrat=namestrat)
    mp.graph.name=modelname
    posthook(mp)
    return mp
end

"""
    modelproto(g::CompGraph; outshape = shape, namestrat=default_namestrat(g); , posthook=validate, kwargs...)

Return an [`ONNX.ModelProto`](@ref) from `g`.

Argument `outshape` is a function which returns a size tuple representing the shape of the output of a given `AbstractVertex`.

Argument `namestrat` determines how nodes in the graph shall be named. Other keyword arguments are passed to the `ModelProto`.

Argument `posthook` will be called with the created `ONNX.ModelProto` as argument before returning it.

Other keyword arguments will be passed to `ONNX.ModelProto`.
"""
function modelproto(g::CompGraph; modelname="model", outshape = shape, namestrat=default_namestrat(g), posthook=validate, kwargs...)
    mp = modelproto(;kwargs...)
    mp.graph = graphproto(g, outshape, namestrat)
    mp.graph.name = modelname
    posthook(mp)
    return mp
end

function infer_inshapes(c::Chain)
    sh = infer_inshapes(first(c))
    sh isa FluxLayer && length(c) > 1 && return infer_inshapes(Chain(Base.tail(c.layers)...))
    return sh
end
infer_inshapes(sc::SkipConnection) = infer_inshapes(sc.layers)
infer_inshapes(l) = infer_inshapes(layertype(l), l)
infer_inshapes(lt::FluxTransparentLayer, ::Any) = lt
infer_inshapes(lt::FluxParLayer, l) = tuple(shape(lt, nin(l)...))

function infer_inshapes(::Any, f)
    ml = methods(f);
    for m in ml.ms
        m.sig isa DataType && return Tuple(infer_shape.(m.sig.types[2:end]))
    end
    return     @show(ntuple(i -> missing, ml.mt.max_args))
end
infer_shape(::Type{<:Any}) = missing
infer_shape(::Type{<:AbstractArray{T,N}}) where {T,N} = ntuple(i -> missing, N)

modelproto(;kwargs...) = ONNX.ModelProto(;
    ir_version=6,
    opset_import=[ONNX.OperatorSetIdProto(version=11)],
    producer_name="ONNXNaiveNASflux.jl",
    producer_version=string(Pkg.Types.Context().env.project.version), # TODO: Ugh....
    kwargs...)

"""
    AbstractProbe

Abstract base class for probes used for serialization.

Idea is that probe will "record" all seen operations based on how methods for those operations on the probe is defined.
"""
abstract type AbstractProbe end

# Called by several activation functions
Base.oftype(::AbstractProbe, x) = x

"""
    ProtoProbe <: AbstractProbe
    ProtoProbe(name, shape, nextname, graph)

Probe which builds an [`ONNX.GraphProto`](@ref) from seen operations.
"""
struct ProtoProbe{N,F,P,S} <: AbstractProbe
    name::N
    shape::S
    nextname::F
    graph::P
end
Base.Broadcast.broadcastable(p::ProtoProbe) = Ref(p)
NaiveNASlib.name(p::ProtoProbe) = p.name
nextname(p::ProtoProbe) = p.nextname
add!(p::ProtoProbe, n) = add!(p.graph, n)
shape(p::ProtoProbe) = p.shape
add_output!(p::ProtoProbe) = push!(p.graph.output, ONNX.ValueInfoProto(name(p), shape(p)))
Base.ndims(p::AbstractProbe) = length(shape(p))
function inputprotoprobe!(gp, name, shape, namestrat)
    push!(gp.input, ONNX.ValueInfoProto(name, shape))
    ProtoProbe(name, shape, namestrat, gp)
end

"""
    newnamestrat(p::ProtoProbe, f, pname=p.name)

Return a new `ProtoProbe` from `p` with `nextname = f` and name `pname`.
"""
newnamestrat(p::ProtoProbe, f, pname=p.name) = ProtoProbe(pname, p.shape, f, p.graph)

"""
    newfrom(p::ProtoProbe, outname::AbstractString, fshape=identity)

Return a new `ProtoProbe` with name `outname`. Argument `fshape` is used to determine a new shape (typically a function).
"""
newfrom(p::ProtoProbe, outname::AbstractString, fshape) = ProtoProbe(outname, nextshape(p, fshape), p.nextname, p.graph)
nextshape(p::AbstractProbe, f::Function) = f(shape(p))

add!(gp::ONNX.GraphProto, np::ONNX.NodeProto) = push!(gp.node, np)

add!(gp::ONNX.GraphProto, tp::ONNX.TensorProto) = push!(gp.initializer, tp)

"""
    Used to get activation functions as [`ONNX.AttributeProto`](@ref)s.
"""
struct ActivationAttributeProbe end

"""
    graphproto()

Return an [`ONNX.GraphProto`](@ref) with all fields initialized to empty arrays.
"""
graphproto(;kwargs...) = ONNX.GraphProto(;
node = ONNX.NodeProto[],
initializer = ONNX.TensorProto[],
input = ONNX.ValueInfoProto[],
output = ONNX.ValueInfoProto[],
value_info = ONNX.ValueInfoProto[],
kwargs...
)

"""
    graphproto(g::CompGraph, outshape = shape, namestrat=default_namestrat(g))

Return an [`ONNX.GraphProto`](@ref) from `g`.

Argument `outshape` is a function which returns the shape of an `AbstractVertex`.

Argument `namestrat` determines how nodes shall be named.
"""
graphproto(g::CompGraph, outshape = shape, namestrat=default_namestrat(g)) = graphproto(g, (recursename.(inputs(g), namestrat) .=> outshape.(inputs(g)))...;namestrat=namestrat)

"""
    graphproto(f, indata::Pair{String, <:Any}...; namestrat = name_runningnr())

Return an [`ONNX.GraphProto`](@ref) from `g`.

Argument indata are name => shape pairs for the input data.

Argument `namestrat` determines how nodes shall be named.
"""
function graphproto(f, indata::Pair{String, <:Any}...; namestrat = name_runningnr())
    gp = graphproto()
    pps = map(indata) do (name, shape)
        inputprotoprobe!(gp, name, shape, namestrat)
    end

    outpps = f(pps...)

    add_outputs!(gp, namestrat, outpps)

    return gp
end
add_outputs!(gp, ns, x) = add_outputs!(gp, ns, (x,))
add_outputs!(gp, ns, pps::NTuple{N, AbstractProbe}) where N  = add_output!.(pps)
function add_outputs!(gp, namestrat, pps::Tuple)
    # At least one of the outputs was not an AbstractProbe
    # This is probably because one of them is a constant
    # If there is at least one AbstractProbe we here assume that one contains the GraphProto for the non-constant ops
    anyprobeind = findfirst(x -> isa(x, AbstractProbe), pps)
    tempprobe = isnothing(anyprobeind) ? ProtoProbe("template", tuple(), namestrat, gp) : pps[anyprobeind]

    output_pps = constant.(pps, tempprobe, namestrat)
    add_output!.(output_pps)
end

# Only purpose is to snag the name in case this is the naming strategy
function (v::NaiveNASlib.MutationVertex)(pps::AbstractProbe...)
    ppsname = map(pps) do pp
        newnamestrat(pp, nextname(pp)(v))
    end
    ppout = base(v)(ppsname...)
    return newnamestrat(ppout, nextname(pps[1]))
end

actfun(::FluxLayer, l) = l.σ
function weightlayer(lt::FluxParLayer, l, pp, optype;attributes = ONNX.AttributeProto[])
    lname = recursename(l, nextname(pp))
    wname, bname = lname .* ("_weight", "_bias")

    add!(pp, ONNX.TensorProto(flipweights(lt, weights(l)), wname))
    inputnames = addbias!(lt, pp, bias(l), bname, [name(pp), wname])

    add!(pp, ONNX.NodeProto(
        input=inputnames,
        output=[lname],
        name=lname,
        attribute = attributes,
        op_type=optype))

    ppl = newfrom(pp, lname, s -> outshape(l, s))
    ppout = actfun(lt, l)(newnamestrat(ppl, f -> join([lname, genname(f)], "_"), lname))
    return newnamestrat(ppout, nextname(ppl))
end

function addbias!(lt, pp, b, name, inputnames)
    add!(pp, ONNX.TensorProto(b, name))
    return vcat(inputnames, name)
end
addbias!(lt, pp, ::Flux.Zeros, name, inputnames) = inputnames

function(l::Flux.Dense)(pp::AbstractProbe)
    ppl = pp
    if !ismissing(shape(pp)) && ndims(pp) == 3
        # Special case: Recurrent -> Dense. This is nothing special in flux as recurrent layers do 2D -> 2D
        # For it to be valid ONNX however we need to add a reshape so that time dimension becomes batch dimension
        outsize = shape(pp)[1]
        lname = recursename(l, nextname(pp))
        ppn = newnamestrat(pp, s -> lname * genname(s))
        ppn = reshape(ppn, nin(l)[], :)
        ppl = newnamestrat(ppn, s -> lname)
    end
    ppout = weightlayer(layertype(l), l, ppl, "Gemm")
    return newnamestrat(ppout, nextname(ppl))
end

(l::Flux.Conv)(pp::AbstractProbe) = weightlayer(layertype(l), l, pp, "Conv"; attributes = attribs(l))

attribs(l) = attribs(layertype(l), l)
attribs(::FluxConvolutional{N}, l) where N = ONNX.AttributeProto.([ "pads", "strides", "dilations"], [padexpand(Val(N), l.pad), reverse(l.stride), reverse(l.dilation)])
attribs(l::Union{MaxPool{N}, MeanPool{N}}) where N = ONNX.AttributeProto.(["kernel_shape", "pads", "strides"],  [reverse(l.k), padexpand(Val(N), l.pad), reverse(l.stride)])

# Interleave padding! (1,2) => [2,1,2,1], (1,1,2,2,3,3) => (3,2,1,3,2,1)
padexpand(::Val{N}, x::NTuple{N}) where N =  repeat(reverse(collect(x)), 2)
padexpand(::Val{N}, x::NTuple{M}) where {N,M} = vcat(collect(x[end-1:-2:1]), collect(x[end:-2:2]))

function(l::Flux.BatchNorm)(pp::AbstractProbe)
    lname = recursename(l, nextname(pp))
    γname, βname, μname, σ²name = lname .* ("_scale", "_bias", "_mean", "_var")

    add!(pp, ONNX.NodeProto(
        input=[name(pp), γname, βname, μname, σ²name],
        output=[lname],
        name=lname,
        attribute = ONNX.AttributeProto.(["epsilon", "momentum"], [l.ϵ, l.momentum]),
        op_type="BatchNormalization"))
    add!(pp, ONNX.TensorProto(l.γ, γname))
    add!(pp, ONNX.TensorProto(l.β, βname)) # Bias not optional for batchnorm
    add!(pp, ONNX.TensorProto(l.μ, μname))
    add!(pp, ONNX.TensorProto(l.σ², σ²name))

    ppout = actfun(layertype(l), l)(newnamestrat(pp, f -> join([lname, genname(f)], "_"), lname))
    return newnamestrat(ppout, nextname(pp))
end
actfun(::FluxBatchNorm, l) = l.λ

# Dropdims because ONNX expects recurrent layers to output tensors of shape [seq_length, num_directions, batch_size, hidden_size] where num_directions is 2 in case of bidirectional and 1 otherwise
# Flux.Recur is not bidirectional so we'll just assume the user wants to also drop num_directions so that recurrent layers can be stacked without hassle.
# Override Flux.Recur with some other method to circumvent this behaviour if not wanted
(m::Flux.Recur)(pps::AbstractProbe...) = dropdims(m.cell(m.state, pps...), dims=3)

(l::Flux.RNNCell)(h, pp::AbstractProbe) = recurrent_node(l, pp, "RNN")
(l::Flux.LSTMCell)(h, pp::AbstractProbe) = recurrent_node(l, pp, "LSTM")

function recurrent_node(l, pp, optype)
    lname = recursename(l, nextname(pp))
    wname, rname, bname = lname .* ("_W", "_R", "_B")

    hsize = size(l.Wh, 2)
    hsattrib = ONNX.AttributeProto("hidden_size", hsize)

    inputnames = [name(pp), wname, rname]

    # Flux weights are of shape [hidden_size, input_size]
    # ONNX wants them on the form [num_directions, hidden_size, input_size] (where num_directions is 2 for bidirectional else 1)
    # To spice things up a bit, all julia arrays are saved in reverse order, i.e we need to create a TensorProto from an array with the arrangement [input_size, hidden_size, num_directions].
    # First transpose the weights into [input_size, hidden_size], then reshape by adding 1 extra dimension
    Wi = permutedims(flipweights(layertype(l), l.Wi, hsize))
    add!(pp, ONNX.TensorProto(reshape(Wi, size(Wi)...,1), wname))
    Wh = permutedims(flipweights(layertype(l), l.Wh, hsize))
    add!(pp, ONNX.TensorProto(reshape(Wh, size(Wh)..., 1), rname))

    if !isa(l.b, Flux.Zeros)
        # ONNX has a separate bias for the recurrent part and wants the concatenation of input and recurrent biases.
        # We'll just hard code it to zeros. Doesn't matter which part is which as they are just added together in the ONNX expression for RNNs.
        b = flipweights(layertype(l), reshape(l.b, :, 1), hsize)
        add!(pp, ONNX.TensorProto(vcat(b, zeros(eltype(b), size(b))), bname))
        push!(inputnames, bname)
    end

    add!(pp, ONNX.NodeProto(
        input=inputnames,
        output=[lname],
        name=lname,
        attribute = push!(activation_attrib(l), hsattrib),
        op_type=optype))

    return newfrom(pp, lname, s -> outshape(l, s))
end

activation_attrib(l) = l.σ(ActivationAttributeProbe())
activation_attrib(l::Flux.LSTMCell) = ONNX.AttributeProto[] #Only default values supported by Flux

Base.tanh(::ActivationAttributeProbe) = rnnactattribs("Tanh")
Flux.elu(::ActivationAttributeProbe, α=1f0) = rnnactattribs("Elu", α)

rnnactattribs(op::AbstractString, α=0f0, β=0f0) = rnnactattribs([op], [α], [β])
rnnactattribs(ops::AbstractVector, αs, βs) = ONNX.AttributeProto.(["activations", "activation_alpha", "activation_beta"], [ops, αs, βs])

function attribfun(fhshape, optype, pps::AbstractProbe...; attributes = ONNX.AttributeProto[], lname = recursename(lowercase(optype), nextname(pps[1])))
    add!(pps[1], ONNX.NodeProto(
    input = collect(name.(pps)),
    output = [lname],
    name=lname,
    attribute = attributes,
    op_type= optype))
    return newfrom(pps[1], lname, fhshape)
end

Flux.relu(pp::AbstractProbe) = attribfun(identity, "Relu", pp)
Flux.elu(pp::AbstractProbe, α=1f0) = attribfun(identity, "Elu", pp; attributes = [ONNX.AttributeProto("alpha", α)])
Flux.selu(pp::AbstractProbe) = attribfun(identity, "Selu", pp)
Flux.selu(pp::AbstractProbe, γ, α) = attribfun(identity, "Selu", pp; attributes = ONNX.AttributeProto.(["gamma", "alpha"], [γ, α]))
Flux.softmax(pp::AbstractProbe; dims=1) =  onnxsoftmax(pp, np_axis = flux2numpydim(dims[end], ndims(pp)))
onnxsoftmax(pp::AbstractProbe; np_axis=1) =  attribfun(identity, "Softmax", pp; attributes=[ONNX.AttributeProto("axis", np_axis)])

(l::Flux.MaxPool)(pp::AbstractProbe) = attribfun(s -> outshape(l, s), "MaxPool", pp; attributes = attribs(l))
(l::Flux.MeanPool)(pp::AbstractProbe) = attribfun(s -> outshape(l, s), "AveragePool", pp; attributes = attribs(l))
(l::Flux.Dropout)(pp::AbstractProbe) = attribfun(identity, "Dropout", pp; attributes = [ONNX.AttributeProto("ratio", l.p)])

globalmeanpool(pp::AbstractProbe, wrap) = globalpool(pp, wrap, "GlobalAveragePool")
globalmaxpool(pp::AbstractProbe, wrap) = globalpool(pp, wrap, "GlobalMaxPool")

function globalpool(pp::AbstractProbe, wrap, type)
     gpp = attribfun(s -> ismissing(s) ? s : (1, 1, s[3:end]...), type, pp)
     ppnext = newnamestrat(gpp, f -> join([gpp.name, genname(f)], "_"), gpp.name)
     wpp = wrap(ppnext)
     return newnamestrat(wpp, nextname(gpp))
end

# Generate explicit combinations as I couldn't figure out how to avoid type piracy with varargs: https://discourse.julialang.org/t/extend-a-varargs-function-for-mixed-types/38233
function generate_elemwise(fm::Module, f, optype, argperms, m=@__MODULE__)
    for argtypes in argperms
        args = ntuple(i -> Symbol(:x, i), length(argtypes))
        sig = map(zip(args, argtypes)) do (a, at)
            isnothing(at) && return a
            :($a::$at)
        end

        @eval m $fm.$f($(sig...)) = elemwisefun($optype, $(args...))
    end
end

"""
    override_broadcast(f::F, argperms) where F

Prevent broadcasting of `f` when invoked with any combination of argument types in argperms.

Needed because broadcasting happens inside several ONNX operations.

For example, `[1,2,3] .+ 4` shall translate to `Add([1,2,3], 4)`, not as `Add(1, 4)`, `Add(2, 4)` and `Add(3, 4)`. One way to accomplish this is to override broadcasting when an `AbstractProbe` is one of the inputs.
"""
function override_broadcast(f::F, argperms, m=@__MODULE__) where F
    Broadcast.Broadcasted
    for Args in (Tuple{args...} for args in argperms)
        @eval m begin
            Base.Broadcast.Broadcasted(f::$F, args::$Args, axes=nothing) = $f(unwrap_broadcast.(args)...)
            Base.Broadcast.Broadcasted{Style}(f::$F, args::$Args, axes=nothing) where Style = $f(unwrap_broadcast.(args)...)
        end
    end
end
unwrap_broadcast(x) = x
unwrap_broadcast(x::Ref) = x.x

dummyfun(x,y) = "dummy $x"

argpermutations(n, args...) = Iterators.product(ntuple(_ -> args, n)...)
argpermswith(t, n::Integer, args...) = (a for a in argpermutations(n, t, args...) if t in a)

function gen_broadcastable_elemwise(f, optype, n=2)
    fs = Symbol(f)
    fm = which(ONNXNaiveNASflux, fs)
    generate_elemwise(fm, fs, optype, argpermswith(AbstractProbe, n, nothing))
    override_broadcast(f, argpermswith(Base.RefValue{<:AbstractProbe}, n, AbstractArray))
end

gen_broadcastable_elemwise(+, "Add")
gen_broadcastable_elemwise(*, "Mul")

function elemwisefun(optype, args...)
    # This mess is only to make sure we first draw the name of the op so that any constants base their name of it
    anyprobe = args[findfirst(x -> isa(x, AbstractProbe), args)]
    oname = recursename(lowercase(optype), nextname(anyprobe))
    nf = name_runningnr()
    refprobe = newnamestrat(anyprobe, f -> join([oname, nf(f)], "_"))
    return attribfun(identity, optype, constant.(args, refprobe, nextname(anyprobe))...; lname=oname)
end

constant(x::AbstractProbe, ::AbstractProbe, ns) = x
function constant(x, pp::AbstractProbe, ns)
    cname = recursename("constant", nextname(pp))
    add!(pp, ONNX.NodeProto(
    input = [],
    output = [cname],
    name=cname,
    attribute = ONNX.AttributeProto.(["value"], [ONNX.TensorProto(x, cname * "_value")]),
    op_type= "Constant"))
    ppo = newfrom(pp, cname, identity)
    return newnamestrat(ppo, ns)
end

function axisfun(fshape, optype, pps::AbstractProbe...; dims, axname="axes")
    fname = recursename(lowercase(optype), nextname(pps[1]))

    attrib = if isempty(dims)
        ONNX.AttributeProto[]
    else
        pok = filter(p -> !ismissing(shape(p)), pps)
        @assert !isempty(pok) "Must have at least one shape to determine axis!"
        np_axis = flux2numpydim.(dims, ndims(pok[1]))
        [ONNX.AttributeProto(axname, np_axis)]
    end

    add!(pps[1], ONNX.NodeProto(
        input = collect(name.(pps)),
        output = [fname],
        name = fname,
        attribute = attrib,
        op_type = optype
    ))
    return newfrom(pps[1], fname, fshape)
end

scal2tup(x) = (x,)
scal2tup(x::Tuple) = x

Base.cat(pps::AbstractProbe...; dims) = axisfun("Concat", pps...; dims=dims, axname="axis") do s
    sumshape = aggshape.(+, vcat(shape.(pps)...)...)
    return ntuple(i -> i in dims ? sumshape[i] : s[i], length(s))
end
Statistics.mean(pp::AbstractProbe; dims=()) = axisfun("ReduceMean", pp; dims=scal2tup(dims)) do s
    return ntuple(i -> i in dims ? 1 : s[i], length(s))
end
Base.dropdims(pp::AbstractProbe; dims) = axisfun(s -> rmdims(s, dims), "Squeeze", pp; dims=scal2tup(dims))

reshape_keepshape(pp::AbstractProbe, shape) = reshape(pp, shape)
Base.reshape(pp::AbstractProbe, shape...) = reshape(pp, shape)
function Base.reshape(pp::AbstractProbe, shape::Tuple)
    fname = recursename("Reshape", nextname(pp))
    sname = fname .* "_shape"
    fluxshape = collect(Int, map(s -> s == Colon() ? -1 : s, shape))

    add!(pp, ONNX.NodeProto(
        input=[name(pp), sname],
        output=[fname],
        name=fname,
        op_type="Reshape"))
    add!(pp, ONNX.TensorProto(reverse(fluxshape), sname))

    fshape = function(s)
        return map(enumerate(fluxshape)) do (ind, new)
            new == -1 && return missing # CBA to figure out how to do this...
            new == 0 && return s[ind]
            return new
        end |> Tuple
    end

    return newfrom(pp, fname, fshape)
end
expanddims(p::AbstractProbe, x, dims) = p

Flux.flatten(pp::AbstractProbe) = flatten(pp, ndims(pp)-1)

function flatten(pp::AbstractProbe, dim)
    fname = recursename("Flatten", nextname(pp))

    add!(pp, ONNX.NodeProto(
        input=[name(pp)],
        output=[fname],
        name=fname,
        attribute = [ONNX.AttributeProto("axis", -dim)],
        op_type="Flatten"))

    fshape = function (s)
        dim == 0 && return (aggshape(*, s), 1)
        absdim = dim < 0 ? length(s) + dim : dim
        return (aggshape(*, s[1:absdim]...), aggshape(*, s[absdim+1:end]))
    end
    return newfrom(pp, fname, fshape)
end


text:
function _indextomask(A, is)
    bitmask = falses(size(A))
    bitmask[is] .= 1

    return bitmask
end
function _slicetomask(A, slices::Tuple)
    slicemask = falses(size(A))
    if all(!isempty, slices)
        slicemask[slices...] .= 1
    end

    return slicemask
end

mask(A::AbstractArray, bitmask::AbstractArray{Bool}) = MaskedArray(A, bitmask)
# must be Vector{<:Integer} not AbstractVector{<:Integer} to avoid method ambiguity
mask(A::AbstractArray, bitmask::Vector{<:Integer}) = mask(A, _indextomask(A, bitmask))
mask(A::AbstractArray, slices::Tuple) = MaskedSliceArray(A, slices)
mask(A::AbstractArray, slice, slices...) = mask(A, (slice, slices...))
mask(A::MaskedArray, bitmask::AbstractArray{Bool}) = mask(A.data, A.mask .* bitmask)
mask(A::MaskedSliceArray, slices::Tuple) =
    mask(A.data, map((s, s̄) -> union(s, s̄), A.slices, to_indices(A.data, slices)))
mask(A::MaskedArray, slices::Tuple) = mask(A.data, A.mask .* _slicetomask(A, slices))
mask(A::MaskedSliceArray, bitmask::AbstractArray{Bool}) =
    mask(A.data, bitmask .* _slicetomask(A, A.slices))

unmask(A::MaskedArray) = A.data
unmask(A::MaskedSliceArray) = A.data

freeze(A::MaskedArray) = A.data .* A.mask
freeze(A::MaskedSliceArray) = A.data .* _slicetomask(A, A.slices)


text:
module Common

using Statistics
using NumericalDataManipulation.CommonMath

function relative_Δ(a¹::Float64, a²::Float64, b¹::Float64, b²::Float64)
    return (a² - a¹) / (b² - b¹)
end

function relative_Δ(a¹::Float64, a²::Float64, b::Float64)
    return (a² - a¹) / b
end

function relative(a::Float64, b::Float64)
    @debug "Relation: a=$a, b=$b, rel=$(a / b)"
    return a / b
end

function issmall(a::Float64, b::Float64; ϵʳᵉˡ=1e-2::Float64)
    return relative(a, b) < ϵʳᵉˡ
end

function Δknot_values(knots::Vector{Float64})
    @assert(!isempty(knots), "An empty knot vector is given.")
    @assert(issorted(knots), "The knot vector must be sorted in an ascending order.")
    @assert(length(knots) ≥ 2, "The knot vector must have 2 knots at least.")
    Δknots = knots[2:end] - knots[1:end-1]
    return Δknots
end

function nonzero(V::Vector{Float64})
    return filter(v -> v ≠ 0.0, V)
end

function finity(a::Float64, b::Float64, c::Float64)
    return (finity(a), finity(b), finity(c))
end
function finity(triple::Tuple{Float64, Float64, Float64})
    return map(number -> finity(number), triple)
end
function finity(v::Vector{Float64})
    return map(number -> finity(number), v)
end

function finity(number::Float64)
    if number == ∞⁻
        return ∞⁻
    elseif number == ∞
        return ∞
    else
        return 1
    end
end

function relfinity(a::Float64, b::Float64, c::Float64; D::Float64 = 1.0, ϵʳᵉˡ::Float64=1e-2)
    return (relfinity(a; D = D, ϵʳᵉˡ = ϵʳᵉˡ),
        relfinity(b; D = D, ϵʳᵉˡ = ϵʳᵉˡ),
        relfinity(c; D = D, ϵʳᵉˡ = ϵʳᵉˡ))
end
function relfinity(triple::Tuple{Float64, Float64, Float64}; D::Float64 = 1.0, ϵʳᵉˡ::Float64=1e-2)
    return map(number -> relfinity(number; D = D, ϵʳᵉˡ = ϵʳᵉˡ), triple)
end
function relfinity(vector::Vector{Float64}; D::Float64 = 1.0, ϵʳᵉˡ::Float64=1e-2)
    return map(number -> relfinity(number; D = D, ϵʳᵉˡ = ϵʳᵉˡ), vector)
end
function relfinity(number::Float64; D::Float64 = 1.0, ϵʳᵉˡ::Float64=1e-2)
    if number == ∞⁻
        return ∞⁻
    elseif number == ∞
        return ∞
    else
        if issmall(number, D; ϵʳᵉˡ = ϵʳᵉˡ)
            return 0
        else
            return 1
        end
    end
end

function sum_until(upper_bound_reached::Function, V::Vector{Float64})
    Σ = 0.0; iᵉ = 0
    for (i, v) ∈ enumerate(V)
        if upper_bound_reached(v, Σ + v)
            break
        end
        if v == ∞ || v == ∞⁻ || v == NaN
            @warn "v = $v with index $i; vector: $V"
        end
        Σ += v; iᵉ = i
    end
    return (Σ, iᵉ)
end

flatten(arr::Array) = mapreduce(x -> isa(x, Array) ? flatten(x) : x, append!, arr,init=[])
flatten(arr::Tuple) = mapreduce(x -> isa(x, Tuple) ? flatten(x) : x, append!, arr,init=[])
flatten(arr::Base.KeySet) = mapreduce(x -> isa(x, Base.KeySet) ? flatten(x) : x, append!, arr,init=[])

export relative_Δ, relative, issmall
export Δknot_values
export nonzero, finity, relfinity
export sum_until
export flatten

end  # module Common


text:
mutable struct RecurrentRVPTests <: TrainingTests
    tests::Vector{RecurrentRVPTest}
    images::AbstractArray{Float32, 4}
    action_images::AbstractArray{Float32, 4}
    actual_value::AbstractArray{Float32, 2}
    actual_reward::AbstractArray{Float32, 2}
    period::Int
    warmup::Int
    function RecurrentRVPTests(tests; period=1000, warmup=0)
        new(
            tests,
            tensorize(get_image.(tests)),
            tensorize(get_action_image.(tests)),
            Flux.unsqueeze(get_value.(tests), 1),
            Flux.unsqueeze(get_reward.(tests), 1),
            period,
            warmup,
        )
    end
end

function gettestfn(tests::RecurrentRVPTests)
    recurrent_inference_fn = recurrentinferencefn(tests)
    function fn(n::Network)
        recurrent_output = recurrent_inference_fn(n)
        reward_mse = getrewardmse(tests, recurrent_output)
        value_mse = getvaluemse(tests, recurrent_output)
        acc, fails = getpolicyacc(tests, recurrent_output)
        reward_mse, value_mse, acc, fails
    end
end

runtests(tests::RecurrentRVPTests, net) = gettestfn(tests)(net)
Base.lastindex(x::RecurrentRVPTests) = x.tests[end]

function df(tests::RecurrentRVPTests)
    test_num = "Test #"
    optimal_action = "Actual (π*)"
    d = Dict(
        test_num=>[],
        optimal_action=>[],
    )
    for (idx,test) in enumerate(tests.tests)
        push!(d[test_num], idx)
        push!(d[optimal_action], join(test.actions, ","))
    end
    DataFrame(
        test_num=>d[test_num],
        optimal_action=>d[optimal_action],
    )
end

function Base.show(io::IO, m::MIME"text/plain", x::RecurrentRVPTests)
    len = length(x.tests)
    title_paren = "all $len"
    pretty_table(
        io,
        df(x),
        show_row_number=false,
        nosubheader=true,
        title="Recurrent RVP Tests ($title_paren)"
    )
end


text:
import Pkg
Pkg.activate(".")

using Revise

include("max_rectangle.jl")
include("plot_grid.jl")
using .MaxRectangle
using .PlotGrid

using MLDataUtils
using ProgressBars
using PyCall

np = pyimport("numpy")
model_selection = pyimport("sklearn.model_selection")
neighbors = pyimport("sklearn.neighbors")
preprocessing = pyimport("sklearn.preprocessing")

plt = pyimport("matplotlib.pyplot")

function distributed_sampling(data::Array{Float64}, dim::Int64)
    println("Length of considered remaining data: $(size(data, 1))")

    nbrs = neighbors.NearestNeighbors(n_neighbors=2, algorithm="kd_tree", n_jobs=-1).fit(data)
    distances, __ = nbrs.kneighbors(data)
    # Use the maximum distance between a point to his nearest neighbors as cell radius
    cell_radius = maximum(distances[:, 2])
    cellsize = cell_radius / sqrt(2) * 2

    @time grid = MaxRectangle.max_rectangle(data, cellsize, dim);
    PlotGrid.plot_grid(grid, 1, 1, false)
    PlotGrid.plot_rasterized_samples(data, grid, cellsize)
end

function main()
    random_seed::Int64 = 1234
    test_size::Float64 = 0.4
    data_dir::String = "/home/finkeldey/workspace/distribution-sampling-jl/data/"

    filenames = filter(filename->occursin("_features.npy", filename), readdir(data_dir))
    train_files, val_files = model_selection.train_test_split(
        filenames,
        test_size=test_size,
        random_state=random_seed
    )
    val_files, __ = model_selection.train_test_split(
        val_files,
        test_size=0.5,
        random_state=random_seed
    )

    features::Array{Array{Float64, 2}} = []
    for (idx, __) in ProgressBar(enumerate(train_files))
        x__ = np.load("$data_dir/$(train_files[idx])")[:, 1:2]
        push!(features, unique(x__[x__[:, 1] .> 0, :], dims=1))
    end
    data = vcat(features...)

    scaler = preprocessing.MinMaxScaler()
    data = scaler.fit_transform(data)
    dim::Int64 = size(data, 2)

    distributed_sampling(data, dim)
end

main()

max_stars_repo_name: effaeff/distribution-sampling-jl

max_issues_repo_name: effaeff/distribution-sampling-jl

max_forks_repo_name: effaeff/distribution-sampling-jl


text:
n_man(sem_fit::SemFit) = n_man(sem_fit.model)

n_man(model::AbstractSemSingle) = n_man(model.observed)


text:
MRANGE(::Type{Float64}) = 10000000
MRANGE(::Type{Float32}) = 10000
IntF(::Type{Float64}) = Int64
IntF(::Type{Float32}) = Int32

@testset "Accuracy (max error in ulp) for $T" for T in (Float32, Float64)
    println("Accuracy tests for $T")

    xx = map(T, vcat(-10:0.0002:10, -1000:0.1:1000))
    fun_table = Dict(SLEEF.exp => Base.exp)
    tol = 1
    test_acc(T, fun_table, xx, tol)

    xx = map(T, vcat(-10:0.0002:10, -1000:0.02:1000))
    fun_table = Dict(SLEEF.asinh => Base.asinh, SLEEF.atanh => Base.atanh)
    tol = 1
    test_acc(T, fun_table, xx, tol)

    xx = map(T, vcat(1:0.0002:10, 1:0.02:1000))
    fun_table = Dict(SLEEF.acosh => Base.acosh)
    tol = 1
    test_acc(T, fun_table, xx, tol)

    xx = T[]
    for i = 1:10000
        s = reinterpret(T, reinterpret(IntF(T), T(pi)/4 * i) - IntF(T)(20))
        e = reinterpret(T, reinterpret(IntF(T), T(pi)/4 * i) + IntF(T)(20))
        d = s
        while d <= e
            append!(xx, d)
            d = reinterpret(T, reinterpret(IntF(T), d) + IntF(T)(1))
        end
    end
    xx = append!(xx, -10:0.0002:10)
    xx = append!(xx, -MRANGE(T):200.1:MRANGE(T))

    fun_table = Dict(SLEEF.sin => Base.sin, SLEEF.cos => Base.cos, SLEEF.tan => Base.tan)
    tol = 1
    test_acc(T, fun_table, xx, tol)

    fun_table = Dict(SLEEF.sin_fast => Base.sin, SLEEF.cos_fast => Base.cos, SLEEF.tan_fast => Base.tan)
    tol = 4
    test_acc(T, fun_table, xx, tol)

    global sin_sincos_fast(x) = (SLEEF.sincos_fast(x))[1]
    global cos_sincos_fast(x) = (SLEEF.sincos_fast(x))[2]
    fun_table = Dict(sin_sincos_fast => Base.sin, cos_sincos_fast => Base.cos)
    tol = 4
    test_acc(T, fun_table, xx, tol)

    global sin_sincos(x) = (SLEEF.sincos(x))[1]
    global cos_sincos(x) = (SLEEF.sincos(x))[2]
    fun_table = Dict(sin_sincos => Base.sin, cos_sincos => Base.cos)
    tol = 1
    test_acc(T, fun_table, xx, tol)

    xx = map(T, vcat(-1:0.00002:1))
    fun_table = Dict(SLEEF.asin_fast => Base.asin, SLEEF.acos_fast => Base.acos)
    tol = 3
    test_acc(T, fun_table, xx, tol)

    fun_table = Dict(SLEEF.asin => asin, SLEEF.acos => Base.acos)
    tol = 1
    test_acc(T, fun_table, xx, tol)

    xx = map(T, vcat(-10:0.0002:10, -10000:0.2:10000, -10000:0.201:10000))
    fun_table = Dict(SLEEF.atan_fast => Base.atan)
    tol = 3
    test_acc(T, fun_table, xx, tol)

    fun_table = Dict(SLEEF.atan => Base.atan)
    tol = 1
    test_acc(T, fun_table, xx, tol)

    xx1 = map(Tuple{T,T}, [zip(-10:0.050:10, -10:0.050:10)...])
    xx2 = map(Tuple{T,T}, [zip(-10:0.051:10, -10:0.052:10)...])
    xx3 = map(Tuple{T,T}, [zip(-100:0.51:100, -100:0.51:100)...])
    xx4 = map(Tuple{T,T}, [zip(-100:0.51:100, -100:0.52:100)...])
    xx = vcat(xx1, xx2, xx3, xx4)

    fun_table = Dict(SLEEF.atan_fast => Base.atan)
    tol = 2.5
    test_acc(T, fun_table, xx, tol)

    fun_table = Dict(SLEEF.atan => Base.atan)
    tol = 1
    test_acc(T, fun_table, xx, tol)

    xx = map(T, vcat(0.0001:0.0001:10, 0.001:0.1:10000, 1.1.^(-1000:1000), 2.1.^(-1000:1000)))
    fun_table = Dict(SLEEF.log_fast => Base.log)
    tol = 3
    test_acc(T, fun_table, xx, tol)

    fun_table = Dict(SLEEF.log => Base.log)
    tol = 1
    test_acc(T, fun_table, xx, tol)

    xx = map(T, vcat(0.0001:0.0001:10, 0.0001:0.1:10000))
    fun_table = Dict(SLEEF.log10 => Base.log10, SLEEF.log2 => Base.log2)
    tol = 1
    test_acc(T, fun_table, xx, tol)

    xx = map(T, vcat(0.0001:0.0001:10, 0.0001:0.1:10000, 10.0.^-(0:0.02:300), -10.0.^-(0:0.02:300)))
    fun_table = Dict(SLEEF.log1p => Base.log1p)
    tol = 1
    test_acc(T, fun_table, xx, tol)

    xx1 = map(Tuple{T,T}, [(x,y) for x = -100:0.20:100, y = 0.1:0.20:100])[:]
    xx2 = map(Tuple{T,T}, [(x,y) for x = -100:0.21:100, y = 0.1:0.22:100])[:]
    xx3 = map(Tuple{T,T}, [(x,y) for x = 2.1, y = -1000:0.1:1000])
    xx = vcat(xx1, xx2, xx2)
    fun_table = Dict(SLEEF.pow => Base.:^)
    tol = 1
    test_acc(T, fun_table, xx, tol)

    xx = map(T, vcat(-10000:0.2:10000, 1.1.^(-1000:1000), 2.1.^(-1000:1000)))
    fun_table = Dict(SLEEF.cbrt_fast => Base.cbrt)
    tol = 2
    test_acc(T, fun_table, xx, tol)

    fun_table = Dict(SLEEF.cbrt => Base.cbrt)
    tol = 1
    test_acc(T, fun_table, xx, tol)

    xx = map(T, vcat(-10:0.0002:10, -120:0.023:1000, -1000:0.02:2000))
    fun_table = Dict(SLEEF.exp2 => Base.exp2)
    tol = 1
    test_acc(T, fun_table, xx, tol)

    xx = map(T, vcat(-10:0.0002:10, -35:0.023:1000, -300:0.01:300))
    fun_table = Dict(SLEEF.exp10 => Base.exp10)
    tol = 1
    test_acc(T, fun_table, xx, tol)

    xx = map(T, vcat(-10:0.0002:10, -1000:0.021:1000, -1000:0.023:1000,
        10.0.^-(0:0.02:300), -10.0.^-(0:0.02:300), 10.0.^(0:0.021:300), -10.0.^-(0:0.021:300)))
    fun_table = Dict(SLEEF.expm1 => Base.expm1)
    tol = 2
    test_acc(T, fun_table, xx, tol)

    xx = map(T, vcat(-10:0.0002:10, -1000:0.02:1000))
    fun_table = Dict(SLEEF.sinh => Base.sinh, SLEEF.cosh => Base.cosh, SLEEF.tanh => Base.tanh)
    tol = 1
    test_acc(T, fun_table, xx, tol)

    @testset "xilogb at arbitrary values" begin
        xd = Dict{T,Int}(T(1e-30) => -100, T(2.31e-11) => -36, T(-1.0) => 0, T(1.0) => 0,
                T(2.31e11) => 37,  T(1e30) => 99)
        for (i,j) in xd
            @test SLEEF.ilogb(i)  === j
        end
    end

end


text:
@testset "1658.minimum-operations-to-reduce-x-to-zero.jl" begin
    @test min_operations([1, 1, 14, 2, 3], 5) == 2
    @test min_operations([5, 6, 7, 8, 9], 4) == -1
    @test min_operations([3, 2, 20, 1, 1, 3], 10) == 5
    @test min_operations(
        [
            8828,
            9581,
            49,
            9818,
            9974,
            9869,
            9991,
            10000,
            10000,
            10000,
            9999,
            9993,
            9904,
            8819,
            1231,
            6309,
        ],
        134365,
    ) == 16
end


text:
#=
BasicPrim:
- Julia version: 1.6.2
- Author: samuel
- Date: 2021-09-10
=#

using SuiteSparseGraphBLAS
using SparseArrays

function stoppingCondition(m, n)

infinito = true

for i = 1:n
        if(m[i] == 0.0)
            return false
    end
end
    return infinito
end

#A is the input matrix, n is the number of nodes, m is the source matrix Output: minimum spanning tree cost
function b_prim(A, n, m)

    d = GBVector{Float64}(n)

    weight = 0.0

    d = A[1,:]
    print("INITIAL WEIGHT: ")
    print(weight)
    print("\n\n")
    while(stoppingCondition(m, n) == false)

        u = argmin(m'+d)
        m[u[2]] = Inf
        weight = weight + d[u[2]]
        print("WEIGHT: ")
        print(weight)
        d = emul(d, A[u[1],:],  BinaryOps.MIN)
        print("  ITERATION FINISHED")
        print("\n\n\n")
    end

    print("COMPLETED")
    print(" ALL THE NODES HAVE BEEN PROCESSED (INF)")
    print("\n\n\n")
    print(m')
    return weight

end


text:
julia> q = Quaternion(1,0,0,0)
julia> q  = Quaternion (1, 2, 3, 4)
       q1 = Quaternion(2, 3, 4, 5)
       q2 = Quaternion(3, 4, 5, 6)
       r = 7.

julia> norm(q)
5.477225575051661

julia> -q
-1 - 2i - 3j - 4k

julia> conj(q)
1 - 2i - 3j - 4k

julia> r + q, q + r
(8.0 + 2.0i + 3.0j + 4.0k,8.0 + 2.0i + 3.0j + 4.0k)

julia> q1 + q2
5 + 7i + 9j + 11k

julia> r*q, q*r
(7.0 + 14.0i + 21.0j + 28.0k,7.0 + 14.0i + 21.0j + 28.0k)

julia> q1*q2, q2*q1, q1*q2 != q2*q1
(-56 + 16i + 24j + 26k,-56 + 18i + 20j + 28k,true)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@everywhere using AccelerationBenchmark
using DataFrames, CSV, Optim, LineSearches, CUTEst, JLD
run_tests_more = true # Run the Moré et al. tests from the O-ACCEL paper
run_tests_cutest = true # Run CUTEst tests

savejld = true
savecsv = true

function gettestset()
    wmax = 20
    m = 10
    ag = InitialStatic(alpha=1.0, scaled=false)
    agscaled = InitialStatic(alpha=1e-4, scaled=true)

    # delta and sigma are the Wolfe condition coefficients
    lshz = HagerZhang(delta = 0.1, sigma = 0.9) # Defaults
    # lslbfgs = HagerZhang(delta = 1e-4, sigma = 0.9)
    # lsgd = HagerZhang(delta = 1e-4, sigma = 0.1)
    # lsacc = HagerZhang(delta = 1e-4, sigma = 0.5)
    lslbfgs = lshz
    lsgd = lshz
    lsacc = lshz

    preconls = Static(alpha=1.0)
    gdprecon = GradientDescent(linesearch = preconls,
                               alphaguess = agscaled)

    solvers = [
        LBFGS(alphaguess = ag, linesearch = lslbfgs, scaleinvH0 = true, m = m),
        OACCEL(alphaguess = ag, linesearch = lsacc,
               nlprecon = gdprecon, wmax = wmax),
        NGMRES(alphaguess = ag, linesearch = lsacc,
               nlprecon = gdprecon, wmax = wmax),
        OACCEL(alphaguess = ag, linesearch = lsacc,
               nlprecon = gdprecon, wmax = m),
        NGMRES(alphaguess = ag, linesearch = lsacc,
               nlprecon = gdprecon, wmax = m),
        GradientDescent(alphaguess = ag, linesearch = lsgd)
    ]
    solvernames = [
        "L-BFGS",
        "O-ACCEL-$wmax",
        "N-GMRES-$wmax",
        "O-ACCEL-$m",
        "N-GMRES-$m",
        "GD"
    ]
    #stoptype = :FunctionTolRelative
    #stoptol  = 1e-10
    stoptype = :GradientTolRelative
    stoptol  = 1e-8
    timelog  = 2 # timelog = 1 can cause very inaccurate timings due to compilation and garbage collection
    maxiter  = 10000
    timelimit = 65.0
    TestSetup(solvers,solvernames,stoptype,stoptol,timelog,maxiter,timelimit)
end

testset = gettestset()
savebase = "data/"

if run_tests_more
    tests = [(:A,100), (:A,200),
             (:B,100), (:B,200),
             (:C,100), (:C,200),
             (:D,500), (:D,1000),
             (:E,100), (:E,200),
             (:F,200), (:F,500),
             (:G,100), (:G,200),
             (:D, 50000), (:D, 100000),
             (:E, 50000), (:E, 100000)]
    seeds = 0:99

    OACCEL2017.runmany(tests,seeds, testset;
                       savejld = savejld, savecsv = savecsv,
                       savebase = savebase)
end

if run_tests_cutest
    min_var = 20
    max_var = 100000
    cutestnames = sort(CUTEst.select(contype=:unc, max_var=max_var,
                                     min_var=min_var))
    oruns = createmeasures(cutestnames, testset)
    mdf = createmeasuredataframe(oruns)

    if savejld
        save(savebase*"/cutest_$(min_var)_$(max_var).jld", "mdf", mdf)
    end
    if savecsv
        CSV.write(savebase*"/cutest_$(min_var)_$(max_var).csv", mdf)
    end
end


text:
include("utils.jl")

build_so("011")

x = [
    "abc",
    "def",
    "gh",
]

ns = Cint[
    3,
    3,
    2
]

y = ccall(
    (:foo, "binary/lib011.so"),
    Void,
    (Ptr{Ptr{Uint8}}, Cint, Ptr{Cint}),
    x, length(x), ns,
)

@assert x[1] == "bcd"
@assert x[2] == "efg"
@assert x[3] == "hi"

@assert y === nothing


text:
module SummaryTest

using OnlineStats, StatsBase, Base.Test

x = randn(500)
x1 = randn(500)
x2 = randn(501)
xs = hcat(x1, x)

@testset "Summary" begin
@testset "Mean / Means" begin
    o = Mean()
    fit!(o, x1, 5)
    @test mean(o) ≈ mean(x1)
    @test nobs(o) == 500
    o = Mean(x2)
    @test mean(o) ≈ mean(x2)
    @test mean(o) == value(o)
    @test center(o, x[1]) == x[1] - mean(o)

    o = Means(2)
    fit!(o, xs, 5)
    @test mean(o) ≈ vec(mean(xs, 1))
    fit!(o, xs, 1)
    o = Means(xs)
    o2 = Means(2)
    fit!(o2, xs)
    @test mean(o) ≈ mean(o2)
    @test center(o, vec(xs[1, :])) == vec(xs[1, :]) - mean(o)

    o = Mean(x1)
    o2 = Mean(x2)
    merge!(o, o2)
    @test mean(o) ≈ mean(vcat(x1, x2))
    merge!(o, o2, :mean)
    merge!(o, o2, :singleton)
end
@testset "Variance / Variances" begin
    o = Variance(x1)
    @test var(o) ≈ var(x1)
    o = Variance()
    fit!(o, x1)
    @test var(o)            ≈ var(x1)
    @test std(o)            ≈ std(x1)
    @test mean(o)           ≈ mean(x1)
    @test var(o)            == value(o)
    @test center(o, x[1])   == x[1] - mean(o)
    @test zscore(o, x[1])   == (x[1] - mean(o)) / std(o)
    o = Variance(x1)
    o2 = Variance(x2)
    merge!(o, o2)
    @test var(o) ≈ var(vcat(x1, x2))

    o = Variances(xs)
    @test var(o) ≈ vec(var(xs, 1))
    @test mean(o) ≈ vec(mean(xs, 1))
    @test std(o) ≈ vec(std(xs, 1))
    @test var(o) == value(o)
    o2 = Variances(2)
    fit!(o2, xs)
    @test var(o) ≈ var(o2)
    @test center(o, vec(xs[1, :])) == vec(xs[1, :]) - mean(o)
    @test zscore(o, vec(xs[1, :])) == (vec(xs[1, :]) - mean(o)) ./ std(o)
    @test zscore(Variances(5), ones(5)) == ones(5)

    o = Variances(xs)
    o2 = Variances(xs)
    merge!(o, o2)
    @test var(o) ≈ vec(var(vcat(xs, xs), 1))
end
@testset "CovMatrix" begin
    o = CovMatrix(xs)
    @test cov(o) ≈ cov(xs)
    @test cor(o) ≈ cor(xs)

    o2 = CovMatrix(2)
    fit!(o2, xs)
    @test cov(o) ≈ cov(o2)
    @test cor(o) ≈ cor(o2)
    @test value(o) ≈ value(o2)
    @test mean(o) ≈ mean(o2)
    @test mean(o2) ≈ vec(mean(xs, 1))
    @test var(o) ≈ vec(var(xs, 1))
    @test std(o) ≈ vec(std(xs, 1))
    fit!(o2, xs, 1)

    o = CovMatrix(xs)
    o2 = CovMatrix(xs)
    merge!(o, o2)
    @test cov(o) ≈ cov(vcat(xs, xs))
end
@testset "Extrema" begin
    o = Extrema(x1)
    @test extrema(o) == extrema(x1)
    o2 = Extrema()
    fit!(o2, x1)
    @test extrema(o) == extrema(o2)

    o = Extrema(x1)
    o2 = Extrema(x2)
    merge!(o, o2)
    @test extrema(o) == extrema(vcat(x1, x2))
end
@testset "QuantileSGD / QuantileMM" begin
    o = QuantileSGD(x1)
    fit!(o, x2, 2)

    o = QuantileMM(x1)
    fit!(o, x2, 2)
end
@testset "Moments" begin
    o = Moments(x1)
    @test mean(o) ≈ mean(x1)
    @test var(o) ≈ var(x1)
    @test std(o) ≈ std(x1)
    @test_approx_eq_eps skewness(o) skewness(x1) .01
    @test_approx_eq_eps kurtosis(o) kurtosis(x1) .1
end
@testset "Diff / Diffs" begin
    Diff()
    Diff(Float64)
    Diff(Float32)
    Diff(Int64)
    Diff(Int32)
    y = randn(100)
    o = Diff(y)
    @test typeof(o) == Diff{Float64}
    @test last(o) == y[end]
    @test diff(o) == y[end] - y[end-1]
    @test value(o) == diff(o)
    y = rand(Int, 100)
    o = Diff(y)
    @test typeof(o) == Diff{Int64}
    @test last(o) == y[end]
    @test diff(o) == y[end] - y[end-1]

    Diffs(10)
    Diffs(Int32, 10)
    y = randn(100, 10)
    o = Diffs(y)
    @test last(o) == vec(y[end, :])
    @test diff(o) == vec(y[end, :] - y[end-1, :])
    @test value(o) == diff(o)
end
@testset "Sum / Sums" begin
    Sum()
    Sum(Float64)
    Sum(Float32)
    Sum(Int64)
    Sum(Int32)
    y = randn(100)
    o = Sum(y)
    @test typeof(o) == Sum{Float64}
    @test sum(o) ≈ sum(y)
    @test value(o) == sum(o)
    y = rand(Int, 100)
    o = Sum(y)
    @test typeof(o) == Sum{Int64}
    @test sum(o) ≈ sum(y)

    Sums(10)
    Sums(Int32, 10)
    y = randn(100, 10)
    o = Sums(y)
    @test sum(o) ≈ vec(sum(y,1))
    @test value(o) == sum(o)
end
end # summary
end #module


text:
macro custom(expr)
    (Meta.isexpr(expr, :call, 3) && expr.args[1] === :~) ||
        error("incorrect macro usage")
    quote
        $(esc(expr.args[2])) = 0.0
    end
end

macro mymodel1(ex)
    # check if expression was modified by the DynamicPPL "compiler"
    if ex == :(y ~ Uniform())
        return esc(:(x ~ Normal()))
    else
        return esc(:(z ~ Exponential()))
    end
end

struct MyModelStruct{T}
    x::T
end
Base.:~(x, y::MyModelStruct) = y.x
macro mymodel2(ex)
    # check if expression was modified by the DynamicPPL "compiler"
    if ex == :(y ~ Uniform())
        # Just returns 42
        return :(4 ~ MyModelStruct(42))
    else
        return :(return -1)
    end
end

@testset "compiler.jl" begin
    @testset "model macro" begin
        @model function testmodel_comp(x, y)
            s ~ InverseGamma(2,3)
            m ~ Normal(0,sqrt(s))

            x ~ Normal(m, sqrt(s))
            y ~ Normal(m, sqrt(s))

            return x, y
        end
        testmodel_comp(1.0, 1.2)

        # check if drawing from the prior works
        @model function testmodel01(x = missing)
            x ~ Normal()
            return x
        end
        f0_mm = testmodel01()
        @test mean(f0_mm() for _ in 1:1000) ≈ 0. atol=0.1

        # Test #544
        @model function testmodel02(x = missing)
            if x === missing
                x = Vector{Float64}(undef, 2)
            end
            x[1] ~ Normal()
            x[2] ~ Normal()
            return x
        end
        f0_mm = testmodel02()
        @test all(x -> isapprox(x, 0; atol = 0.1), mean(f0_mm() for _ in 1:1000))

        @model function testmodel03(x = missing)
            x ~ Bernoulli(0.5)
            return x
        end
        f01_mm = testmodel03()
        @test mean(f01_mm() for _ in 1:1000) ≈ 0.5 atol=0.1

        # test if we get the correct return values
        @model function testmodel1(x1, x2)
            s ~ InverseGamma(2,3)
            m ~ Normal(0,sqrt(s))

            x1 ~ Normal(m, sqrt(s))
            x2 ~ Normal(m, sqrt(s))

            return x1, x2
        end
        f1_mm = testmodel1(1., 10.)
        @test f1_mm() == (1, 10)

        # alternatives with keyword arguments
        testmodel1kw(; x1, x2) = testmodel1(x1, x2)
        f1_mm = testmodel1kw(x1 = 1., x2 = 10.)
        @test f1_mm() == (1, 10)

        @model function testmodel2(; x1, x2)
            s ~ InverseGamma(2,3)
            m ~ Normal(0,sqrt(s))

            x1 ~ Normal(m, sqrt(s))
            x2 ~ Normal(m, sqrt(s))

            return x1, x2
        end
        f1_mm = testmodel2(x1=1., x2=10.)
        @test f1_mm() == (1, 10)

        @info "Testing the compiler's ability to catch bad models..."

        # Test for assertions in observe statements.
        @model function brokentestmodel_observe1(x1, x2)
            s ~ InverseGamma(2,3)
            m ~ Normal(0,sqrt(s))

            x1 ~ Normal(m, sqrt(s))
            x2 ~ x1 + 2

            return x1, x2
        end

        btest = brokentestmodel_observe1(1., 2.)
        @test_throws ArgumentError btest()

        @model function brokentestmodel_observe2(x)
            s ~ InverseGamma(2,3)
            m ~ Normal(0,sqrt(s))

            x = Vector{Float64}(undef, 2)
            x ~ [Normal(m, sqrt(s)), 2.0]

            return x
        end

        btest = brokentestmodel_observe2([1., 2.])
        @test_throws ArgumentError btest()

        # Test for assertions in assume statements.
        @model function brokentestmodel_assume1()
            s ~ InverseGamma(2,3)
            m ~ Normal(0,sqrt(s))

            x1 ~ Normal(m, sqrt(s))
            x2 ~ x1 + 2

            return x1, x2
        end

        btest = brokentestmodel_assume1()
        @test_throws ArgumentError btest()

        @model function brokentestmodel_assume2()
            s ~ InverseGamma(2,3)
            m ~ Normal(0,sqrt(s))

            x = Vector{Float64}(undef, 2)
            x ~ [Normal(m, sqrt(s)), 2.0]

            return x
        end

        btest = brokentestmodel_assume2()
        @test_throws ArgumentError btest()

        # Test missing input arguments
        @model function testmodel_missing1(x)
            x ~ Bernoulli(0.5)
            return x
        end
        @test_throws MethodError testmodel_missing1()

        # Test missing initialization for vector observation turned parameter
        @model function testmodel_missing2(x)
            x[1] ~ Bernoulli(0.5)
            return x
        end
        @test_throws MethodError testmodel_missing2(missing)()

        # Test use of internal names
        @model function testmodel_missing3(x)
            x[1] ~ Bernoulli(0.5)
            global varinfo_ = __varinfo__
            global sampler_ = __sampler__
            global model_ = __model__
            global context_ = __context__
            global rng_ = __rng__
            global lp = getlogp(__varinfo__)
            return x
        end
        model = testmodel_missing3([1.0])
        varinfo = VarInfo(model)
        @test getlogp(varinfo) == lp
        @test varinfo_ isa AbstractVarInfo
        @test model_ === model
        @test sampler_ === SampleFromPrior()
        @test context_ === DefaultContext()
        @test rng_ isa Random.AbstractRNG

        # disable warnings
        @model function testmodel_missing4(x)
            x[1] ~ Bernoulli(0.5)
            global varinfo_ = __varinfo__
            global sampler_ = __sampler__
            global model_ = __model__
            global context_ = __context__
            global rng_ = __rng__
            global lp = getlogp(__varinfo__)
            return x
        end false
        lpold = lp
        model = testmodel_missing4([1.0])
        varinfo = VarInfo(model)
        @test getlogp(varinfo) == lp == lpold

        # test DPPL#61
        @model function testmodel_missing5(z)
            m ~ Normal()
            z[1:end] ~ MvNormal(fill(m, length(z)), 1.0)
            return m
        end
        model = testmodel_missing5(rand(10))
        @test all(z -> isapprox(z, 0; atol = 0.2), mean(model() for _ in 1:1000))

        # test Turing#1464
        @model function gdemo(x)
            s ~ InverseGamma(2, 3)
            m ~ Normal(0, sqrt(s))
            for i in eachindex(x)
                x[i] ~ Normal(m, sqrt(s))
            end
        end
        x = [1.0, missing]
        VarInfo(gdemo(x))
        @test ismissing(x[2])

        # https://github.com/TuringLang/Turing.jl/issues/1464#issuecomment-731153615
        vi = VarInfo(gdemo(x))
        @test haskey(vi.metadata, :x)
        vi = VarInfo(gdemo(x))
        @test haskey(vi.metadata, :x)
    end
    @testset "nested model" begin
        function makemodel(p)
            @model function testmodel(x)
                x[1] ~ Bernoulli(p)
                global lp = getlogp(__varinfo__)
                return x
            end
            return testmodel
        end
        model = makemodel(0.5)([1.0])
        varinfo = VarInfo(model)
        @test getlogp(varinfo) == lp
    end
    @testset "user-defined variable name" begin
        @model f1() = x ~ NamedDist(Normal(), :y)
        @model f2() = x ~ NamedDist(Normal(), @varname(y[2][:,1]))
        @model f3() = x ~ NamedDist(Normal(), @varname(y[1]))
        vi1 = VarInfo(f1())
        vi2 = VarInfo(f2())
        vi3 = VarInfo(f3())
        @test haskey(vi1.metadata, :y)
        @test vi1.metadata.y.vns[1] == VarName(:y)
        @test haskey(vi2.metadata, :y)
        @test vi2.metadata.y.vns[1] == VarName(:y, ((2,), (Colon(), 1)))
        @test haskey(vi3.metadata, :y)
        @test vi3.metadata.y.vns[1] == VarName(:y, ((1,),))
    end
    @testset "custom tilde" begin
        @model demo() = begin
            $(@custom m ~ Normal())
            return m
        end
        model = demo()
        @test all(iszero(model()) for _ in 1:1000)
    end
    @testset "docstring" begin
        "This is a test"
        @model function demo(x)
            m ~ Normal()
            x ~ Normal(m, 1)
        end

        s = @doc(demo)
        @test string(s) == "This is a test\n"

        # Verify that adding docstring didn't completely break execution of model
        m = demo(0.)
        @test m() isa Float64
    end
    @testset "type annotations" begin
        @model function demo_without(x)
            x ~ Normal()
        end
        @test isempty(VarInfo(demo_without(0.0)))

        @model function demo_with(x::Real)
            x ~ Normal()
        end
        @test isempty(VarInfo(demo_with(0.0)))
    end

    @testset "macros within model" begin
        # Macro expansion
        @model function demo1()
            @mymodel1(y ~ Uniform())
        end

        @test haskey(VarInfo(demo1()), @varname(x))

        # Interpolation
        # Will fail if:
        # 1. Compiler expands `y ~ Uniform()` before expanding the macros
        #    => returns -1.
        # 2. `@mymodel` is expanded before entire `@model` has been
        #    expanded => errors since `MyModelStruct` is not a distribution,
        #    and hence `tilde_observe` errors.
        @model function demo2()
            $(@mymodel2(y ~ Uniform()))
        end
        @test demo2()() == 42
    end
end


text:
using NCCL
using Test
using CUDAdrv
using CUDAnative

CUDAnative.initialize()

@show id = NCCL.UniqueID()
# initialise a communcator for nranks=1
# @show comm = NCCL.Communicator(1, id, 0)

@show comm = NCCL.Communicator(CUDAdrv.devices())


text:
export HermitianOperator, SymmetricOperator, AdjointOperator, TransposeOperator

struct HermitianOperator{T<:Number,B<:Operator} <: Operator{T}
    op::B
    uplo::Char
end

HermitianOperator(B::Operator{T}, uplo::Symbol=:U) where {T<:Number} = HermitianOperator{T,typeof(B)}(B, char_uplo(uplo))
convert(::Type{Operator{T}},A::HermitianOperator) where {T}=HermitianOperator(convert(Operator{T},A.op), A.uplo)

domainspace(P::HermitianOperator)=domainspace(P.op)
rangespace(P::HermitianOperator)=rangespace(P.op)
domain(P::HermitianOperator)=domain(P.op)
bandwidths(P::HermitianOperator) = (b = bandwidth(P.op, P.uplo == 'L' ? 1 : P.uplo == 'U' ? 2 : 0); (b, b))

function getindex(P::HermitianOperator,k::Integer,j::Integer)
    if P.uplo == 'L'
        if j > k
            conj(P.op[j,k])
        else
            P.op[k,j]
        end
    elseif P.uplo == 'U'
        if j < k
            conj(P.op[j,k])
        else
            P.op[k,j]
        end
    end
end

Hermitian(A::Operator, uplo::Symbol)=HermitianOperator(A, uplo)
Hermitian(A::Operator)=HermitianOperator(A)

struct SymmetricOperator{T<:Number,B<:Operator} <: Operator{T}
    op::B
    uplo::Char
end

SymmetricOperator(B::Operator{T}, uplo::Symbol=:U) where {T<:Number} = SymmetricOperator{T,typeof(B)}(B, char_uplo(uplo))
convert(::Type{Operator{T}},A::SymmetricOperator) where {T}=SymmetricOperator(convert(Operator{T},A.op), A.uplo)

domainspace(P::SymmetricOperator)=domainspace(P.op)
rangespace(P::SymmetricOperator)=rangespace(P.op)
domain(P::SymmetricOperator)=domain(P.op)
bandwidths(P::SymmetricOperator) = (b = bandwidth(P.op, P.uplo == 'L' ? 1 : P.uplo == 'U' ? 2 : 0); (b, b))

function getindex(P::SymmetricOperator,k::Integer,j::Integer)
    if P.uplo == 'L'
        if j > k
            P.op[j,k]
        else
            P.op[k,j]
        end
    elseif P.uplo == 'U'
        if j < k
            P.op[j,k]
        else
            P.op[k,j]
        end
    end
end

Symmetric(A::Operator, uplo::Symbol)=SymmetricOperator(A, uplo)
Symmetric(A::Operator)=SymmetricOperator(A)

struct AdjointOperator{T<:Number,B<:Operator} <: Operator{T}
    op::B
end

AdjointOperator(B::Operator{T}) where {T<:Number}=AdjointOperator{T,typeof(B)}(B)
convert(::Type{Operator{T}},A::AdjointOperator) where {T}=AdjointOperator(convert(Operator{T},A.op))

domainspace(P::AdjointOperator)=rangespace(P.op)
rangespace(P::AdjointOperator)=domainspace(P.op)
domain(P::AdjointOperator)=domain(P.op)
bandwidths(P::AdjointOperator) = reverse(bandwidths(P.op))

getindex(P::AdjointOperator,k::Integer,j::Integer) = conj(P.op[j,k])

function BandedMatrix(S::SubOperator{T,TO}) where {T,TO<:AdjointOperator}
    kr,jr=parentindices(S)
    adjoint(BandedMatrix(view(parent(S).op,jr,kr)))
end

adjoint(A::Operator)=AdjointOperator(A)

struct TransposeOperator{T<:Number,B<:Operator} <: Operator{T}
    op::B
end

TransposeOperator(B::Operator{T}) where {T<:Number}=TransposeOperator{T,typeof(B)}(B)
convert(::Type{Operator{T}},A::TransposeOperator) where {T}=TransposeOperator(convert(Operator{T},A.op))

domainspace(P::TransposeOperator)=rangespace(P.op)
rangespace(P::TransposeOperator)=domainspace(P.op)
domain(P::TransposeOperator)=domain(P.op)
bandwidths(P::TransposeOperator) = reverse(bandwidths(P.op))

getindex(P::TransposeOperator,k::Integer,j::Integer) = P.op[j,k]

function BandedMatrix(S::SubOperator{T,TO}) where {T,TO<:TransposeOperator}
    kr,jr=parentindices(S)
    transpose(BandedMatrix(view(parent(S).op,jr,kr)))
end

transpose(A::Operator)=TransposeOperator(A)


text:
module GridapGmshTests

using GridapGmsh
using Test

@testset "gmsh" begin
  include("gmshTests.jl")
end

@testset "GmshDiscreteModel" begin
 include("GmshDiscreteModelsTests.jl")
end

@testset "LagrangianFESpaces" begin
 include("LagrangianFESpacesTests.jl")
end

end # module


text:
# Copyright (c) 2020-2021, Bank of Canada
# All rights reserved.

# ----------------------------------------
# 1. Frequency types
# ----------------------------------------

"""
    abstract type Frequency end

The abstract supertype for all frequencies.
See also: [`Unit`](@ref) and [`YPFrequency`](@ref)
"""
abstract type Frequency end

"""
    struct Unit <: Frequency end

Represents a non-dimensional frequency (not associated with the calendar).
See also: [`YPFrequency`](@ref)
"""
struct Unit <: Frequency end

"""
    abstract type YPFrequency{N} <: Frequency end;

Represents a calendar frequency based on a number of periods in a year.
The type parameter `N` is the number of periods and must be a positive integer.
See also: [`Yearly`](@ref), [`Quarterly`](@ref), [`Monthly`](@ref)
"""
abstract type YPFrequency{N} <: Frequency end;

struct Yearly <: YPFrequency{1} end
struct Quarterly <: YPFrequency{4} end
struct Monthly <: YPFrequency{12} end

"""
    Yearly, Quarterly, Monthly

Frequencies corresponding to 1, 4, and 12 periods per year.
See also: [`YPFrequency`](@ref)
"""
Yearly, Quarterly, Monthly

# ----------------------------------------
# 2. MIT (moment in time) and Duration
# ----------------------------------------

primitive type MIT{F <: Frequency} <: Signed 64 end
primitive type Duration{F <: Frequency} <: Signed 64 end

"""
    MIT{F <: Frequency}, Duration{F <: Frequency}

Two types representing a
moment in time (like 2020Q1 or 2020-01-01) and
duration (the quantity of time between two moments).

Both of these have a Frequency as a type parameter and both
internally are represented by integer values.

If you imagine a time axis of the given `Frequency`, `MIT` values are
ordinal (correspond to points) while `Duration` values are cardinal
(correspond to distances).
"""
MIT, Duration

# --------------------------
# conversions with Int
MIT{F}(x::Int) where F <: Frequency = reinterpret(MIT{F}, x)
Int(x::MIT) = reinterpret(Int, x)

Duration{F}(x::Int) where F <: Frequency = reinterpret(Duration{F}, x)
Int(x::Duration) = reinterpret(Int, x)

# -------------------------
# frequencyof()

"""
frequencyof(x), frequencyof(T)

Return the Frequency type of the given value `x` or type `T`.
"""
function frequencyof end

# throw an error, except for values and types that have a frequency
# Q: should we return `nothing` instead?
# A: No. We assume that `frequencyof` returns a subtype of `Frequency`.
frequencyof(::T) where T = frequencyof(T)
frequencyof(T::Type) = throw(ArgumentError("$(T) does not have a frequency."))
frequencyof(::Type{MIT{F}}) where F <: Frequency = F
frequencyof(::Type{Duration{F}}) where F <: Frequency = F
# AbstractArray{<:MIT} cover MIT-ranges and vectors of MIT
frequencyof(::Type{<:AbstractArray{MIT{F}}}) where F <: Frequency = F
frequencyof(::Type{<:AbstractArray{Duration{F}}}) where F <: Frequency = F

# -------------------------
# YP-specific stuff

# construct from year and period
"""
    MIT{F}(year, period) where {F <: YPFrequency}

Construct an [`MIT`](@ref) instance from year and period. This is valid only for
frequencies subtyped from [`YPFrequency`](@ref).
"""
MIT{F}(y::Integer, p::Integer) where F <: YPFrequency{N} where N = MIT{F}(N * Int(y) + Int(p) - 1)

"""
    mit2yp(x::MIT)

Recover the year and period from a given [`MIT`](@ref) value. This is valid only
if the frequency is subtyped from [`YPFrequency`](@ref).
"""
function mit2yp end
mit2yp(x) = throw(ArgumentError("Value of type $(typeof(x)) cannot be represented as (year, period) pair. "))

# Careful with values near 0 and negative values
# the remainder `p` returned by divrem() has the same sign
# as the argument being divided. We need 0 <= p <= N-1, so
# in this case we add N to p and subtract 1 from y
# (since 1y = N*p we preserve the property that x = N*y+p)
@inline function mit2yp(x::MIT{<:YPFrequency{N}}) where {N}
    (y, p) = divrem(Int(x), N);
    p < 0 ? (y - 1, p + N + 1) : (y, p + 1)
end

@inline year(x) = mit2yp(x)[1]
@inline period(x) = mit2yp(x)[2]

"""
    year(x::MIT), period(x::MIT)

Return the year and period of an [`MIT`](@ref) value `x`. This only makes
sense if the frequency of `x` is subtyped from [`YPFrequency`](@ref).
"""
year, period

@inline mm(y::Integer, p::Integer) = MIT{Monthly}(y, p)
@inline qq(y::Integer, p::Integer) = MIT{Quarterly}(y, p)
@inline yy(y::Integer, p::Integer=1) = MIT{Yearly}(y, p)

"""
    mm(year, period), qq(year, period), yy(year, period=1)

IRIS type constructors for [`MIT`](@ref) values with frequencies that have
year and period.
"""
mm, qq, yy

"""
    pp(year, period; N)

Construct an [`MIT`](@ref) with frequency [`YPFrequency{N}`](@ref). For
[`Quarterly`](@ref), [`Monthly`](@ref), [`Yearly`](@ref), use [`qq`](@ref),
[`mm`](@ref), [`yy`](@ref) instead of this.

"""
@inline pp(y::Integer, p::Integer; N::Integer) = MIT{YPFrequency{N}}(y, p)

# -------------------------
# pretty printing

Base.show(io::IO, m::MIT{Unit}) = print(io, Int(m), 'U')
function Base.show(io::IO, m::MIT{F}) where F <: YPFrequency{N} where N
    if isconcretetype(F)
    periodletter = first("$(F)")
    else
        periodletter =  N == 1 ? 'Y' :
                        N == 4 ? 'Q' :
                        N == 12 ? 'M' : 'P';
    end
    print(io, year(m), periodletter)
    if N > 1
        # print(io, rpad(period(m), length(string(N))))
        print(io, period(m))
    end
end

Base.string(m::MIT{<:Frequency}) = repr(m)
Base.print(io::IO, m::MIT{<:Frequency}) = print(io, string(m))

Base.show(io::IO, d::Duration) = print(io, Int(d))
Base.string(d::Duration) = repr(d)
Base.print(io::IO, d::Duration) = print(io, string(d))

# -------------------------
# Arithmetic operations with MIT and Duration

# Arithmetic with MIT and Duration are limited to adding/subtracting and
# comparisons and even these are valid only with matching frequencies.
# We want to tightly control what works and what doesn't and also
# we want to provide meaningful error messages why operations
# that don't work are disabled.

Base.promote_rule(IT::Type{<:Integer}, MT::Type{<:MIT}) = throw(ArgumentError("Invalid arithmetic operation with $IT and $MT"))
Base.promote_rule(IT::Type{<:Integer}, DT::Type{<:Duration}) = throw(ArgumentError("Invalid arithmetic operation with $IT and $DT"))

mixed_freq_error(T1::Type, T2::Type) = throw(ArgumentError("Mixing frequencies not allowed: $(frequencyof(T1)) and $(frequencyof(T2))."))
mixed_freq_error(::T1, ::T2) where {T1,T2} = mixed_freq_error(T1, T2)
mixed_freq_error(T1::Type, T2::Type, T3::Type) = throw(ArgumentError("Mixing frequencies not allowed: $(frequencyof(T1)), $(frequencyof(T2)) and $(frequencyof(T3))."))
mixed_freq_error(::T1, ::T2, ::T3) where {T1,T2,T3} = mixed_freq_error(T1, T2, T3)

Base.promote_rule(T1::Type{<:MIT}, T2::Type{<:MIT}) = mixed_freq_error(T1, T2)
Base.promote_rule(::Type{MIT{F}}, T2::Type{MIT{F}}) where {F<:Frequency} = T1

# -------------------
# subtraction

# difference of two MIT is a Duration
Base.:(-)(l::MIT, r::MIT) = mixed_freq_error(l, r)
Base.:(-)(l::MIT{F}, r::MIT{F}) where F <: Frequency = Duration{F}(Int(l) - Int(r))
# difference of MIT and Duration is an MIT
Base.:(-)(l::MIT, r::Duration) = mixed_freq_error(l, r)
Base.:(-)(l::MIT{F}, r::Duration{F}) where F <: Frequency = MIT{F}(Int(l) - Int(r))
# difference of MIT and Integer is an MIT -- the Integer value is interpreted as a Duration of the appropriate frequency
Base.:(-)(l::MIT{F}, r::Integer) where F <: Frequency = MIT{F}(Int(l) - Int(r))
# Difference of two Duration is a Duration
Base.:(-)(l::Duration, r::Duration) = mixed_freq_error(l, r)
Base.:(-)(l::Duration{F}) where F <: Frequency = Duration{F}(-Int(l))
Base.:(-)(l::Duration{F}, r::Duration{F}) where F <: Frequency = Duration{F}(Int(l) - Int(r))
# difference of Duration and Integer is a Duration -- the Integer value is interpreted as a Duration of the same frequency
Base.:(-)(l::Duration{F}, r::Integer) where F <: Frequency = Duration{F}(Int(l) - Int(r))

Base.rem(x::Duration, y::Duration) = mixed_freq_error(x,y)
Base.rem(x::Duration{F}, y::Duration{F}) where {F<:Frequency} = Duration{F}(rem(Int(x),Int(y)))
Base.div(x::Duration, y::Duration, args...) = mixed_freq_error(x,y)
Base.div(x::Duration{F}, y::Duration{F}, args...) where {F<:Frequency} = Duration{F}(div(Int(x),Int(y),args...))

# -------------------
# Comparison for equality

# we can always compare two values of MIT and Duration for equality. But they are only equal if they are of the exact same type
Base.:(==)(l::Union{MIT,Duration}, r::Union{MIT,Duration}) = typeof(l) == typeof(r) && Int(l) == Int(r)
# Base.:(==)(l::Duration, r::Duration) = frequencyof(l) == frequencyof(r) && Int(l) == Int(r)
# Base.:(==)(::Duration, ::MIT) = false
# Base.:(==)(::MIT, ::Duration) = false
# when comparing to other Integer types, we compare the numerical values
Base.:(==)(l::Union{MIT,Duration}, r::Integer) = Int(l) == r
Base.:(==)(l::Integer, r::Union{MIT,Duration}) = l == Int(r)
# Base.:(==)(l::Duration, r::Integer) = Int(l) == r
# Base.:(==)(l::Integer, r::Duration) = l == Int(r)

# -------------------
# Comparison for order

# For MIT and Duration values, they can be ordered only if they are of the exact same type, otherwise we throw an ArgumentError.
Base.:(<)(l::MIT, r::MIT) = mixed_freq_error(l, r)
Base.:(<)(l::MIT{F}, r::MIT{F}) where F <: Frequency = Int(l) < Int(r)
Base.:(<)(l::Duration{F}, r::Duration{F}) where F <: Frequency = Int(l) < Int(r)
Base.:(<)(l::Duration, r::Duration) = mixed_freq_error(l, r)
Base.:(<)(l::MIT, r::Duration) = throw(ArgumentError("Illegal comparison of $(typeof(l)) and $(typeof(r))."))
Base.:(<)(l::Duration, r::MIT) = throw(ArgumentError("Illegal comparison of $(typeof(l)) and $(typeof(r))."))

# Comparison of Duration with Int is needed for indexing and iterating time series.
# Base.:(<)(l::Int, r::Duration) = l < Int(r)
# Base.:(<)(l::Duration, r::Int) = Int(l) < r

# <= is expressed as < or ==
Base.:(<=)(l::Union{MIT,Duration}, r::Union{MIT,Duration}) = (l < r) || (l == r)
# Base.:(<=)(l::Duration, r::Int) = (l < r) || (l == r)
# Base.:(<=)(l::Int, r::Duration) = (l < r) || (l == r)

# -------------------
# addition

# addition of two MIT is not allowed
Base.:(+)(::MIT, ::MIT) = throw(ArgumentError("Illegal addition of two `MIT` values."))
# addition of two Duration is valid only if they are of the same frequency
Base.:(+)(l::Duration, r::Duration) = mixed_freq_error(l, r)
Base.:(+)(l::Duration{F}, r::Duration{F}) where F <: Frequency = Duration{F}(Int(l) + Int(r))
# addition of MIT and Duration gives an MIT
Base.:(+)(l::MIT, r::Duration) = mixed_freq_error(l, r)
Base.:(+)(l::MIT{F}, r::Duration{F}) where F <: Frequency = MIT{F}(Int(l) + Int(r))
# addition of MIT or Duration with an Integer yields the same type as the MIT/Duration argument
Base.:(+)(l::Union{MIT,Duration}, r::Integer) = oftype(l, Int(l) + r)
Base.:(+)(l::Integer, r::Union{MIT,Duration}) = oftype(r, l + Int(r))

# NOTE: the rules above are meant to catch illegal arithmetic (where the units don't make sense).
# For indexing and iterating TSeries it's more convenient to return Int rather than Duration, however
# we choose to have the checks in place.

# -------------------
# one(x) is meant to be a dimensionless 1, so that's what we do
Base.one(::Union{MIT,Duration,Type{<:MIT},Type{<:Duration}}) = Int(1)

# -------------------
# Conversion to Float64 - that's needed for plotting
(T::Type{<:AbstractFloat})(x::MIT) = convert(T, Int(x))
# In the special case of YPFrequency we want the year to be the whole part and the period to be the fractional part.
(T::Type{<:AbstractFloat})(x::MIT{<:YPFrequency{N}}) where N = convert(T, ((y, p) = mit2yp(x); y + (p - 1) / N))
Base.promote_rule(::Type{<:MIT}, ::Type{T}) where T <: AbstractFloat = T

# ----------------------------------------
# 2.2 MIT{T} vector and dict support
# ----------------------------------------

# added so MIT can be used as dictionary keys
Base.hash(x::MIT{T}) where T <: Frequency = hash(("$T", Int(x)))

# # added for sorting Vector{MIT{T}} where T <: Frequency
# Base.sub_with_overflow(x::MIT{T}, y::MIT{T}) where T <: Frequency = begin
#     Base.checked_ssub_int(reinterpret(Int, x), reinterpret(Int, y))
# end

# ----------------------------------------
# 4 Convenience constants
# ----------------------------------------

struct _FConst{F <: Frequency} end
struct _FPConst{F <: Frequency,P} end
Base.:(*)(y::Integer, ::_FConst{F}) where F <: Frequency = MIT{F}(y)
Base.:(*)(y::Integer, ::_FPConst{F,P}) where {F <: Frequency,P} = MIT{F}(y, P)

Base.show(io::IO, C::_FConst) = print(io, 1C)
Base.show(io::IO, C::_FPConst) = print(io, 1C)

export Y, U, Q1, Q2, Q3, Q4, M1, M2, M3, M4, M5, M6, M7, M8, M9, M10, M11, M12
global const U = _FConst{Unit}()
global const Y = _FPConst{Yearly,1}()
global const Q1 = _FPConst{Quarterly,1}()
global const Q2 = _FPConst{Quarterly,2}()
global const Q3 = _FPConst{Quarterly,3}()
global const Q4 = _FPConst{Quarterly,4}()
global const M1 = _FPConst{Monthly,1}()
global const M2 = _FPConst{Monthly,2}()
global const M3 = _FPConst{Monthly,3}()
global const M4 = _FPConst{Monthly,4}()
global const M5 = _FPConst{Monthly,5}()
global const M6 = _FPConst{Monthly,6}()
global const M7 = _FPConst{Monthly,7}()
global const M8 = _FPConst{Monthly,8}()
global const M9 = _FPConst{Monthly,9}()
global const M10 = _FPConst{Monthly,10}()
global const M11 = _FPConst{Monthly,11}()
global const M12 = _FPConst{Monthly,12}()

# ----------------------------------------
# 5 Ranges of MIT
# ----------------------------------------

Base.:(:)(start::MIT, stop::MIT) = mixed_freq_error(start, stop)
Base.:(:)(start::MIT{F}, stop::MIT{F}) where F <: Frequency = UnitRange{MIT{F}}(start, stop)

Base.:(:)(::Int, ::MIT) = my_range_error()
Base.:(:)(::MIT, ::Int) = my_range_error()
my_range_error() = throw(ArgumentError("""Cannot mix Int and MIT in the same range.
    If you're using `begin` or `end` to index a TSeries make sure to use them at both ends.
    For example, instead of s[2:end] use s[begin+1:end].
"""))

Base.:(:)(start::MIT, step::Int, stop::MIT) = mixed_freq_error(start, stop)
Base.:(:)(start::MIT, step::Duration, stop::MIT) = mixed_freq_error(start, step, stop)
Base.:(:)(start::MIT{F}, step::Int, stop::MIT{F}) where {F<:Frequency} = start:Duration{F}(step):stop
Base.:(:)(start::MIT{F}, step::Duration{F}, stop::MIT{F}) where {F<:Frequency} = StepRange(start, step, stop)

Base.length(rng::UnitRange{<:MIT}) = convert(Int, last(rng) - first(rng) + 1)
Base.step(rng::UnitRange{<:MIT}) = convert(Int, 1)

Base.union(l::UnitRange{<:MIT}, r::UnitRange{<:MIT}) = mixed_freq_error(l, r)
Base.union(l::UnitRange{MIT{F}}, r::UnitRange{MIT{F}}) where F <: Frequency = min(first(l), first(r)):max(last(l), last(r))

# Base.issubset(l::UnitRange{<:MIT}, r::UnitRange{<:MIT}) = false
# Base.issubset(l::UnitRange{MIT{F}}, r::UnitRange{MIT{F}}) where F <: Frequency = first(r) <= first(l) && last(l) <= last(r)