text:
" Looks up character against input classes, returns integer code as Char, or
  else the input if the character was not found. "
function table_lookup(chr::Char, table::Array{ASCIIString, 1})
  for pos in 1:length(table)
    if chr in table[pos]
        return Char('0'+(pos-1))
    end
  end
  return chr
end

" Replace all the 'find' patterns with the corresponding 'replace' strings. "
function replace_all(str::ByteString, from::Array{Regex,1}, to::Array{Base.SubstitutionString{ASCIIString}, 1}, display=false::Bool)
  nstr = str
  for pos in 1:length(from)
    nstr = replace(nstr, from[pos], to[pos])
    if display
      println(nstr, " = ", from[pos], " :: ", to[pos])
    end
  end
  return nstr
end

" 'Squashes' a string by reducing any repeated characters to only one instance. "
function squash(str::ByteString)
  lc = 0
  return filter(str) do x
          if x != lc
           lc = x
           return true
          end
          return false
         end
end

" Lowercase and strip non-alpha chars from a word. Naturally asciifies it. "
function prep(str::ByteString)
  return ascii(replace(lowercase(str), r"[^a-z]", ""))
end

"""
  `soundex(str)`

  Transforms a string into its Soundex code.

  The Russell Soundex code is designed primarily for use with English names and has some
  known drawbacks, including a sensitivity to the first letter of a name (Christina is `c623`
  and similar Kristina is `k623`) and loss of some audible differences (Kant and Knuth, `k530`).
  The resulting code is always 1-letter and 3-digits.
 """
function soundex(str::ByteString)
  lstr = prep(str)

  #Squash repetitions.
  lstr = squash(lstr)

  soundex_table = ["aehiouwy", "bfpv", "cgjkqsxz", "dt", "l", "mn", "r"]

  #Look up soundex coding for 2:end
  body = map(x -> table_lookup(x, soundex_table), lstr[2:end])

  #Remove 0's
  body = join(split(body, '0'))

  #Pad with 0's
  body = body * "000"

  #Join first letter with trimmed body.
  return join([string(lstr[1]), body[1:3]])
end

"""
  fuzzy_soundex(str)

  Transforms a string into its Fuzzy Soundex code.

  The Fuzzy Soundex was an attempt to improve the reliability of Soundex. Like
  Phonex, it introduces some multi-character replacements as a prelude to a lookup
  table. Unlike Soundex, Phonex or Phonix, Fuzzy Soundex has a 1-letter 4-digit key.
"""
function fuzzy_soundex(str::ByteString)

  lstr = prep(str)

  #replace digrams
  fs_find = [r"ca", r"c[ck]", r"ce",r"ch$",r"ch?l",r"ch?r",r"ci",r"co",r"^[ct][sz]",r"cu",r"cy",r"dg",r"gh",r"^gn",r"^[hw]r", r"^hw", r"^kn|ng", r"ma?c",r"nst",r"^nt",r"p[fh]",r"rd?t$",r"sch",r"ti[ao]",r"tch"]
  fs_repl = [s"ka", s"kk", s"se",s"kk", s"kl",s"kr",s"si",s"ko",s"ss",s"ku",s"sy",s"gg",s"hh",s"nn",s"rr",s"ww",s"nn",s"mk",s"nss",s"tt",s"ff",s"rr",s"sss",s"sio",s"chh"]
  lstr = replace_all(lstr, fs_find, fs_repl)

  #crazy-ass fuzzy soundex table
  fuzzy_soundex_table = ["aehiouwy", "bfpv", "", "dt", "l", "mn", "r", "gkjqx", "", "csz"]
  body = map(x -> table_lookup(x, fuzzy_soundex_table), lstr[2:end])

  #remove duplicates
  body = squash(body)

  #Remove 0's
  body = join(split(body, '0'))

  #Pad with 0's
  body = body * "000"

  #Join first letter with trimmed body.
  return join([string(lstr[1]), body[1:min(4,end)]])
end

"""
  `phonex(str)`

  Transforms a string into its Phonex code.

  Lait & Randell's Phonex encoding scheme can be viewed as an improved version
  of Soundex. Like Soundex, it produces a 1-letter 3-digit code, but a range
  of modifications make it more resiliant to encoding errors involving the first
  character of a word, as well as other issues caused by interactions between
  characters within the rest of the word. """
function phonex(str::ByteString)

  lstr = prep(str)

  # remove trailing s
  if lstr[end] == 's'
    lstr = lstr[1:end-1]
  end

  #Remove duplicates
  lstr = squash(lstr)

  #Find/replace leading patterns
  phonex_pre_find = [r"^kn", r"^ph", r"^wr", r"^h", r"^[eiouy]", r"^p", r"^v", r"^[kq]", r"^j", r"^z"]
  phonex_pre_repl = [s"n", s"f", s"r", s"", s"a", s"b", s"f", s"c", s"g",s"s"]
  lstr = replace_all(lstr, phonex_pre_find, phonex_pre_repl)

  #Preparatory find/replace for main coding.
  phonex_main_find = [r"[dt]c",r"[lr]([^aeiou$])",r"([mn])g",r"a|e|h|i|o|u|w|y",r"a+"]
  phonex_main_repl = [s"c",s"\1", s"\1",s"a",s""]
  body = replace_all(lstr[2:end], phonex_main_find, phonex_main_repl)

  #look up phonex coding for 2:end
  phonex_table = ["","bfpv", "cgjkqsxz", "dt", "l", "mn", "r"]
  body = map(x -> table_lookup(x, phonex_table), body)

  #Pad with 0's
  body = body * "000"

  #Join first letter with trimmed body.
  return join([string(lstr[1]), body[1:3]])
end

"""
  `phonix(str)`

  Transforms a string into its Phonix code.

  T. N. Gadd's Phonix encoding scheme is based on the Soundex scheme, and is
  similar to the related Phonex scheme. The Phonix scheme uses a large set of
  hand-crafted rules specific to English text and pronunciation. A performance
  penalty might be expected as a result of this large ruleset.
"""
function phonix(str::ByteString)

  lstr = prep(str)

  #Apply this giant list of rules.
  phonix_find = [r"dg", r"c([oau])", r"c[yi]", r"ce", r"^cl([aeiou])", r"ck", r"[gj]c$", r"^ch?r([aeiou])", r"^wr", r"nc", r"ct", r"ph", r"aa", r"sch", r"btl", r"ght", r"augh", r"([aeiou])lj([aeiou])", r"lough", r"^q", r"^kn", r"^gn|gn$", r"(\w)gn([^aeiou])", r"ghn", r"gne$", r"ghne", r"gnes$", r"^ps", r"^pt", r"^cz", r"([aeiou]})wz(\w)", r"(\w)cz(\w)", r"lz", r"rz", r"(\w)z([aeiou])", r"zz", r"([aeiou])z(\w)", r"hrough", r"ough", r"([aeiou]})q([aeiou])", r"([aeiou])j([aeiou])", r"^yj([aeiou])", r"^gh", r"([aeiou])e$", r"^cy", r"nx", r"^pf", r"dt$", r"([td])l$", r"yth", r"^ts?j([aeiou])", r"^ts([aeiou])", r"tch", r"([aeiou])wsk", r"^[pm]n([aeiou])", r"([aeiou])stl", r"tnt$", r"eaux$", r"exci", r"x", r"ned$", r"jr", r"ee$", r"zs", r"([aeiou])h?r([^aeiou]|$)", r"re$", r"lle", r"([^aeiou])le(s?)$", r"e$", r"es$", r"([aeiou])ss", r"([aeiou])mb$", r"mpts", r"mps", r"mpt", r"^[aeiou]"]
  phonix_repl = [s"g", s"k", s"si", s"se", s"kl", s"k", s"k", s"kr", s"r", s"nk", s"kt", s"f", s"ar", s"sh", s"tl", s"t", s"arf", s"ld", s"low", s"kw", s"n", s"n", s"n", s"n", s"n", s"ne", s"ns", s"s", s"t", s"c", s"z", s"ch", s"lsh", s"rsh", s"s", s"ts", s"ts", s"rew", s"of", s"kw", s"y", s"y", s"g", s"gh", s"s", s"nks", s"f", s"t", s"il", s"ith", s"ch", s"t", s"ch", s"vsike", s"n", s"sl", s"ent", s"oh", s"ecs", s"ecs", s"nd", s"dr", s"ea", s"s", s"ah", s"ar", s"le", s"ile\1", s"", s"s", s"as", s"m", s"mps", s"ms", s"mt", s"v"]
  lstr = replace_all(lstr, phonix_find, phonix_repl)

  #remove duplicates
  lstr = squash(lstr)

  #look up phonix encoding for 2:end
  phonix_table = ["aehiouwy", "bp", "cgjkq", "dt", "l", "mn", "r", "fv", "sxz"]
  body = map(x -> table_lookup(x, phonix_table), lstr[2:end])

  #remove vowels
  body = replace(body, r"0+", "")

  #Pad with 0's
  body = body * "000"

  #Join first letter with trimmed body.
  return join([string(lstr[1]), body[1:3]])
end

"""
  `nysiis(str[, len=6])`

  Transform a string according to the New York State Identification and
  Intelligence System (NYSIIS) encoding scheme.

  Taft's NYSIIS scheme is reputed to be fairly popular. It is designed for
  English names, using a strict application of 1,2 and 3-letter substitutions.
  While simple enough, it does not appear particularly robust to common
  variations (Peter/Pete = `patar`/`pat`; Christina/Kristina = `chrast`/`crasta`).
"""
function nysiis(str::ByteString, len=6)

  lstr = prep(str)

  #Find/replace some initial leading/closing patterns.
  nysiis_pre_find = [r"^mac", r"^kn", r"^k", r"p[hf]", r"^sch", r"e[ie]$", r"[drn]t$|[rn]d$"]
  nysiss_pre_repl = [s"mcc", s"n", s"c", s"ff", s"sss", s"y", s"d"]
  lstr = replace_all(lstr, nysiis_pre_find, nysiss_pre_repl)

  #Find/replace rest of rules.
  #Oddly, these duplicate some work which is done in the prefix.
  nysiis_find = [r"ev", r"[eiou]", r"q", r"z", r"m|kn", r"k", r"sch", r"ph", r"([^aeiou])h(.)|(.)h([^aeiou])", r"([^aeiou])h$", r"([aeiou])w", r"s$", r"a$", r"as$",r"ay$"]
  nysiis_repl = [s"af", s"a", s"g", s"s", s"n", s"c", s"sss", s"ff", s"\1\2", s"\1", s"\1a", s"", s"", s"", s"y"]
  body = replace_all(lstr[2:end], nysiis_find, nysiis_repl)

  #remove duplicates
  body = squash(body)

  #Add first letter back on
  return join([lstr[1], body[1:min(len-1,end)]])
end

"""
  `caverphone(str)`

  Transforms a string into its Caverphone representation.

  David Hood's Caverphone algorithm was designed to assist in matching between
  electoral rolls in New Zealand, and is optimised for local accents. The output
  of the algorithm is a six-character code, with the final characters being padded
  with '1's in the case where the transformed string would otherwise have been
  shorter.
"""
function caverphone(str::ByteString)

  lstr = prep(str)

  #Apply the caverphone replacements
  cp_find = [r"^([crt]|en)ough", r"^gn", r"mb$", r"cq", r"c([iey])", r"tch", r"[cqx]", r"v", r"dg", r"ti([ao])", r"d", r"ph", r"b", r"sh", r"z", r"[eioua]", r"^3",  r"3gh3", r"gh", r"g", r"([stkfmn])+", r"w(h?[3y])", r"w",r"^h", r"h", r"r([3y])", r"r", r"l([3y])", r"l", r"j", r"y3", r"y"]
  cp_repl = [ s"\g<1>u2f", s"2n", s"m2", s"2q", s"s\1", s"2ch", s"k", s"f", s"2g", s"si\1", s"t", s"fh", s"p", s"s2", s"s", s"3", s"a", s"3kh3", s"22", s"k", s"\1", s"W\1", s"2", s"a", s"2", s"R\1", s"2", s"L\1", s"2", s"y",s"Y3",s"2"]
  lstr = replace_all(lstr, cp_find, cp_repl)

  #re-lower the string for consistency
  lstr = lowercase(lstr)

  #remove 2/3
  lstr = replace(lstr, r"[23]+", "")

  #pad
  lstr = lstr * "111111"

  #return first 6
  return lstr[1:6]
end

"""
  `metaphone(str[, len=4])`

  Transforms a word into its Metaphone reresentation.

  Lawrence Phillips' Metaphone technique is applicable to a range of sound-alike
  words. It produces a code which uses a sixteen-character alphabet, with 'x' for
  'sh' sounds and 'ø' for 'th' sounds, and no vowels apart from where they are the
  first letter. In some cases, such as the word 'persuade', its representation may
  be unintuitive regarding pronunciation.

  For comparison purposes, 4 characters are usually used, but you may vary the
  returned length `len` if you wish for a longer representation.
"""
function metaphone(str::ByteString,len::Int=4)

  lstr = prep(str)

  #substitution rules
  metaphone_find = [r"^x",  r"x",  r"mb",  r"sch",  r"t?ch|sh",  r"cia|[st]i[ao]",  r"s?c([eiy])(.+)",  r"(.+)dg([eiy])",  r"d", r"gh([^aeiou$])",  r"gn$",  r"gned$",  r"c|g|q",  r"ph|v",  r"th",  r"^wh",  r"[wy]([^aeiou])",  r"z",  r"^[gk]n",  r"^pm"]
  metaphone_repl = [s"s",  s"ks",  s"m",  s"sk",  s"x",  s"xa",  s"s\1\2", s"\1j\2",  s"t", s"h\1",  s"n",  s"ned", s"k",  s"f", s"0", s"w", s"\1", s"s",  s"n",  s"m"]
  lstr = replace_all(lstr, metaphone_find, metaphone_repl)

  #remove duplicates
  lstr = squash(lstr)

  #remove vowels apart from the first
  lstr = replace(lstr, r"([^^])[aeiou]+", s"\1")

  #return first len letters (if that many)
  return lstr[1:min(len,end)]
end

"""
  `double_metaphone(str)`

  Transforms a word into its Double-Metaphone reresentation.

  One of the issues often noted with Metaphone (and indeed other phonetic
  encoding schemes) is the dubious applicability to non-English names. The
  Double-Metaphone scheme tries to account for this, using a very large set
  of hand-crafted rules to reflex a varied set of European and Asian
  pronunciations.

  The technique is known as the 'Double' because it uses two alternative
  strings to represent a word. If the two representations turn out to be
  equivalent, both are returned, and should used as alternate keys.

  As such, this function returns either an `ASCIIString[rep1, rep2]`, or a single
  `ASCIIString`.

  NB: This implementation is based on another implementation in a different
  language, and some of the logic is speculative. It has been tested,
  but probably not enough.
"""
function double_metaphone(str::ByteString)

  lstr = prep(str)

  #Table of search patterns
  s_find = [r"^ps|^x", r"^sugar", r"sh(o[le][mzk]|eim)", r"sh", r"si[oa]", r"^sm", r"^s([nlw])", r"sz|sc[eiy]", r"sche([rn])", r"sch([aeiou])", r"^sch([^aeiouw])", r"sc"]
  j_find = [r"^jose", r"^ja", r"j$", r"j", r"^[gkp]n",r"wr",r"mb"]
  c_find = [r"([^aeiou]a)ch([^ie])",r"^cae",r"^chia",r"chae",r"^ch(arac|aris|or[^e]|ym|ia|em)",r"([ao]r)ch(.[std])", r"([aoue]|^)ch([nrlmbhfvw])",r"^mc", r"wi[ct]z", r"cz", r"ccia",  r"(^a|u)cc([eih][^u])", r"cc([eih][^u])",r"c[ckgq]", r"ci([aeo])", r"c([iey])", r"ch",r"([^^][^m])c", r"c"]
  d_find = [r"dg[iey]", r"dg",r"d[dt]?", r"([^m][^ae])ier",r"([yi])sl"]
  g_find = [r"([^aeiou])gh", r"^ghi", r"^gh", r"([bhd].{1,3})gh", r"([trlcg].u)gh", r"([^^i])gh", r"gh", r"^ogn", r"gn(e?y)", r"gn", r"gli",r"^g[eiy]([spblyn])",r"([drm][aou])ng(er|e?y)", r"([or])gy", r"g(e?[ry])",r"gier",r"g([iey])t", r"g([iey])", r"g"]
  h_find = [r"([^^aeiou])h([^aeiou])", r"([ia])ll([eao])$", r"ph", r"pb", r"q"]
  dbmeta_find = vcat(j_find, s_find, c_find, d_find , g_find, h_find, [r"([oa]i)s", r"tion", r"tia|tch", r"^th([oa])m", r"th", r"v", r"^w[aeiou]", r"^wh", r"([aeiou])w(sk[iy]|$)?",  r"([ie][oa]u)x", r"zh", r"z([oia])"])

  #Table of replacements
  s_repl = [s"s", s"sugar", s"s\1", s"x", s"s", s"sm", s"s\1", s"s", s"xe\1", s"sk\1", s"x\1", s"sk"]
  j_repl = [s"hose", s"ja", s"j", s"j", s"n", s"r", s"m"]
  c_repl = [s"\1k\2", s"sae", s"kia", s"kae",s"k\1",s"\1k\2", s"\1k\2", s"mk", s"ts", s"s", s"xa", s"\1ks\2", s"x\1", s"k", s"s\1", s"s", s"x", s"\1x", s"k"]
  d_repl = [s"j",s"tk",s"t", s"\1ie", s"\1l"]
  g_repl = [s"\1k", s"ji", s"k", s"\1", s"\1f", s"\1k", s"", s"okn", s"kn\1", s"n", s"kli", s"k\1", s"\1nj\2", s"\1jy", s"k\1",s"jier", s"k\1t", s"j\1", s"k"]
  h_repl = [s"\1\2", s"\1l\2", s"f",s"p",s"k"]
  dbmeta_repl = vcat(j_repl, s_repl, c_repl, d_repl, g_repl, h_repl, [s"\1", s"xn", s"x", s"t\1m", s"0", s"f", s"a", s"a", s"",  s"", s"\1ks", s"j", s"ts\1"])

  #Table of alternative replacements
  s_alt = [s"s", s"xugar", s"s\1", s"x", s"x", s"xm", s"x\1", s"s", s"ske\1", s"sk\1", s"s\1", s"sk"]
  j_alt = [s"hose", s"a", s"", s"h", s"n",s"r",s"m"]
  c_alt = [s"\1k\2", s"sae", s"kia", s"xae", s"k\1", s"\1k\2", s"\1k\2", s"mk", s"fx", s"x", s"xa", s"\1ks\2", s"x\1", s"k", s"x\1",s"s", s"x", s"\1k", s"k"]
  d_alt = [s"j",s"tk",s"t", s"\1ier", s"\1l"]
  g_alt = [s"\1k", s"ji", s"k", s"\1", s"\1f", s"\1k", s"", s"on", s"kn\1", s"kn", s"li", s"j\1", s"\1nj\2", s"\1jy", s"k\1", s"jier", s"k\1t", s"k\1", s"k"]
  h_alt = [s"\1\2", s"\1\2", s"f",s"p",s"k"]
  dbmeta_altr = vcat(j_alt, s_alt, c_alt, d_alt, g_alt, h_alt, [s"s", s"xn", s"x", s"t\1m", s"t",s"f", s"f", s"a", s"f", s"\1ks", s"j", s"s\1"])

  #Apply replacements
  astr = replace_all(lstr, dbmeta_find, dbmeta_altr)
  lstr = replace_all(lstr, dbmeta_find, dbmeta_repl)

  lstr = squash(lstr)
  astr = squash(astr)

  #Handle leading and other vowels
  lstr = replace(lstr, r"^[aeiouy]", "a")
  astr = replace(astr, r"^[aeiouy]", "a")
  lstr = replace(lstr, r"([^^])[aeiouyw]+", s"\1")
  astr = replace(astr, r"([^^])[aeiouyw]+", s"\1")

  #If the results are the same, return only one string.
  cmbi = [lstr, astr]
  if cmbi[1] == cmbi[2]
    return cmbi[1]
  end
  #Else, return both options.
  return cmbi
end

"""
  `match_rating_encode(str)`

  Transform a string into the representation used for the Match Rating Approach.

  This is the string encoding system used within the Western Airlines 'Match Rating'
  algorithm. It is a small set of simple transforms, most of which appear in other
  schemes.

  - For the similarity measure between two strings, use `match_rating(onestr, twostr)`

  - For an automatic binary response about the closeness of strings, call
  `meets_match_rating(onestr, twostr)`
"""
function match_rating_encode(str::ByteString)
  lstr = prep(str)

  #replace non-leading vowels
  lstr = string(lstr[1], replace(lstr[2:end], r"[aeiou]+", ""))

  #remove duplicate consonants
  lstr = squash(lstr)

  #return 6 chars, 1:3 and end-3:end (if there are 6 chars).
  if length(lstr) > 6
   lstr = join([lstr[1:3],lstr[end-3:end]])
  end
  return lstr
end

"""
  `lein(str)`

  Transform a string into its Lein Technique representation.

  I could not find much on the origins of this technique, but
  from inspection of the algorithm it appears to be an early
  competitor of Soundex. Like Soundex, it uses no bigram information,
  and produces a 1-letter 3-number code.
"""
function lein(str::ByteString)
  lstr = prep(str)

  #replace vowels and 'hwy'
  lstr = join([lstr[1], replace(lstr[2:end], r"[aeiouhwy]+", "")])

  #remove duplicates and truncate
  lstr = squash(lstr)[1:min(4,end)]

  #perform table lookup
  lein_table = ["", "dt", "mn", "lr", "bfpv", "cjkgqs", "xz"]
  lstr = string(lstr[1], map(x -> table_lookup(x, lein_table), lstr[2:end]))

  #pad with 0s
  lstr = lstr * ("0" ^ max(0, 4 - length(lstr)))

  return lstr
end

"""
  `roger_root(str)`

  Transforms a string into its Roger Root representation.

  A rather mysterious encoding scheme, somewhat more exotic than the
  similarly-sourced `lein` algorithm. The Roger Root maps a combination
  of bigrams and single characters to one or two digit codes, resulting
  in a 5-digit code for each word.
"""
function roger_root(str::ByteString)

  lstr = prep(str)

  rr_find = [r"^[aiou]", r"^[bp]", r"^(c[eiy]|t?s|z)", r"^t?[sc]{1,2}h", r"^([ckqx]|d?g)", r"^[dt]", r"^([gp]?f|ph|v)", r"^g?m", r"^[kgp]?n", r"^h", r"^j", r"^l", r"^w?r", r"^w", r"^y", r"[bp]+", r"(c[iey]|t?s|z)+", r"(t?[sc]{1,2}h|j)+", r"([ckqx]|d?g)+", r"([dt])+", r"([fv]|ph)+", r"l+", r"m+", r"n+", r"r+"]
  rr_repl = [s"1", s"09", s"00", s"06", s"07", s"01", s"08", s"03", s"02", s"2", s"3", s"05", s"04", s"4", s"5", s"9", s"0", s"6", s"7", s"1", s"8", s"5", s"3", s"2", s"4"]

  lstr = replace_all(lstr, rr_find, rr_repl)

  lstr = replace(lstr, r"[aeiouhwy]+", "")

  lstr = lstr * "00000"

  return lstr[1:5]
end

"""
  `canada(str)`

  Transforms a string into its Census Modified Statistics Canada representation.

  This is an extremely simple code, merely removing vowels and squashing letters,
  with no replacements or mapping to digits.
"""
function canada(str::ByteString)

  lstr = prep(str)

  lstr = squash(lstr)

  body = replace(lstr[2:end], r"[aeiouy]+", "")

  return string(lstr[1], body)
end

export soundex, metaphone, phonex, phonix, nysiis, double_metaphone, match_rating_encode, fuzzy_soundex, caverphone, lein, roger_root, canada


text:
export viz_policy, percCorrect, get_belief

"""
Helper function to draw aircraft
"""
function getACString(theta,x,y,fill,draw,width=2.0)
    return "\\node[aircraft top,fill="*fill*",draw="*draw*", minimum width="*string(width)*"cm,rotate="*string(theta)*",scale = 0.35] at (axis cs:"*string(x)*", "*string(y)*") {};"
end

function get_belief(pstate::Vector{Float64}, grid::RectangleGrid,interp::Bool=false)
    belief = spzeros(NSTATES, 1)
    indices, weights = interpolants(grid, pstate)
    if !interp
        indices = indices[findmax(weights)[2]]
        weights = 1.0
    end
    for i = 1:length(indices)
        belief[indices[i]] = weights[i]
    end # for i
    return belief
end # function get_belief

function percCorrect(nnetPath,tablePath)

    nnet = read_network(nnetPath);
    Q = h5open(tablePath, "r") do file
        read(file, "y")
    end
    X = h5open(tablePath, "r") do file
        read(file, "X")
    end

    rnges = h5open(tablePath, "r") do file
        read(file, "ranges")
    end
    means = h5open(tablePath, "r") do file
        read(file, "means")
    end

    Q = Q.*rnges[end].+means[end]
    X = X.*nnet.ranges[1:end-1].+nnet.means[1:end-1]
    nnet_policy = argmax(evaluate_network_multiple(nnet,X),dims=[1])
    table_policy = argmax(Q,dims=[1])
    return sum(nnet_policy.==table_policy)/length(nnet_policy)
end

"""
Visualize the neural network and/or table policies
"""
function viz_policy(;nnetPath::AbstractString="",tablePath::AbstractString="",batch_size::Int=-1,nnetRect=true)
    if nnetPath=="" && tablePath==""
        println("Please provide at least either a path to the neural network file or training data file")
        return
    end

    # Use aircraft shapes tikz style package
    resetPGFPlotsPreamble()
    pushPGFPlotsPreamble("\\usepackage{aircraftshapes}")

    # Load network if given
    nnet = nothing
    if nnetPath!=""
        nnet = read_network(nnetPath);
    end

    # Load table if given
    Q = nothing
    if tablePath!=""
        Q = h5open(tablePath, "r") do file
            read(file, "y")
        end

        rnges = h5open(tablePath, "r") do file
            read(file, "ranges")
        end
        means = h5open(tablePath, "r") do file
            read(file, "means")
        end

        Q = Q.*rnges[end].+means[end]
    end

    grid  = RectangleGrid(RANGES,THETAS,PSIS,OWNSPEEDS,INTRSPEEDS)

    COC = RGB(1.,1.,1.) # white
    SR = RGB(.0,.0,.5) # navy
    SL = RGB(.0,.600,.0) # green
    WR = RGB(.5,.5,.5) # grey
    WL = RGB(.7,.9,.0) # neon green
    ra_colors = [SL,WL,COC,WR,SR]
    bg_colors = [COC]

    # Create scatter plot classes for color key
    sc_string = "{"
    for i=0:4
        define_color("ra_$i",  ra_colors[i+1])
        if i==2
            sc_string *= "ra_$i={mark=square, style={black, mark options={fill=ra_$i}, mark size=6}},"
        else
            sc_string *= "ra_$i={style={ra_$i, mark size=6}},"
        end
    end

    # Color key as a scatter plot
    sc_string=sc_string[1:end-1]*"}"
    xx = [-1.5,-1.5,-1.5, -1.5, -1.5]
    yy = [1.65,1.15,0.65, 0.15, -0.35]
    zz = ["ra_0","ra_1","ra_2","ra_3","ra_4"]
    sc = string(sc_string)

    # Create manipulatable plotting tool
    @manipulate for psi  = convert(Array{Int32,1},round.(rad2deg.(PSIS))),
        ownspeed = OWNSPEEDS,
        intrspeed = INTRSPEEDS,
        zoom = 1.0,
        nbin = [100,150,200,250,1000],
        savePlot = [false,true],
        xshift = 0.0,
        yshift = 0.0,
        xscale = 1.0,
        yscale = 1.0

        # Ensure that zoom and scale factors don't result in dividing by zero or negative numbers
        zoom = zoom <0.1 ? 0.1 : zoom
        xscale = xscale <0.1 ? 0.1 : xscale
        yscale = yscale <0.1 ? 0.1 : yscale

        # Compute neural network Q-values if a neural network is given
        if nnet != nothing

            #Load table with the inputs needed to plot the heat map
            inputsNet= zeros(5,nbin*nbin)
            ind = 1
            for i=range((-1*RANGEMAX/zoom/xscale- xshift)*1000,stop=(RANGEMAX/zoom/xscale- xshift)*1000,length=nbin)
                for j=range((-1*RANGEMAX/zoom/yscale - yshift)*1000,stop=(RANGEMAX/zoom/yscale- yshift)*1000,length=nbin)
                    r = sqrt(i^2+j^2)
                    th = atan(j,i)
                    if nnetRect
                        inputsNet[:,ind] = [i,j,deg2rad.(psi),ownspeed,intrspeed]
                    else
                        inputsNet[:,ind] = [r,th,deg2rad.(psi),ownspeed,intrspeed]
                    end
                    ind = ind+1
                end
            end

            #Calculate all of the Q values from the input array
            q_nnet = zeros(5,nbin*nbin);
            ind = 1
            batch_size = batch_size<0 ? nbin*nbin : batch_size

            # Compute network outputs in batches
            while ind+batch_size<nbin*nbin
                input = inputsNet[:,ind:(ind+batch_size-1)]
                q_nnet[:,ind:(ind+batch_size-1)] = evaluate_network_multiple(nnet,input)
                ind=ind+batch_size
            end
            input = inputsNet[:,ind:end]
            q_nnet[:,ind:end] = evaluate_network_multiple(nnet,input)
        end

        # Q Table Heat Map
        function get_heat1(x::Float64, y::Float64)
            r = sqrt(x^2+y^2)*1000
            if r>RANGEMAX*1000
                return 0.0
            end
            th = atan(y,x)
            bel = get_belief([r,th,deg2rad.(psi),ownspeed,intrspeed],grid,false)
            qvals = Q[:,bel.rowval[1]]
            return rad2deg.(ACTIONS[findmax(qvals)[2]])
        end # function get_heat

       ind = 1
       #Neural Net Heat Map
       function get_heat2(x::Float64, y::Float64)
            qvals2 = q_nnet[:,ind]
            ind += 1
            r = sqrt(x^2+y^2)*1000
            th = atan(y,x)
            if r>RANGEMAX*1000
                return 0.0
            end
            return rad2deg.(ACTIONS[findmax(qvals2)[2]])
       end # function get_heat2

        numPlots = 3
        numPlots = nnet==nothing ? numPlots-1 : numPlots
        numPlots = Q==nothing ? numPlots-1 : numPlots
        numPlots = savePlot ? numPlots-1 : numPlots

        g = GroupPlot(numPlots, 1, style="height={8cm}, width={8cm}",groupStyle = "horizontal sep=2cm")

        # Plot Q-table, if possible
        if Q != nothing
            push!(g, Axis([
                        Plots.Image(get_heat1, (-1*RANGEMAX/zoom/xscale - xshift, RANGEMAX/zoom/xscale- xshift),
                            (-1*RANGEMAX/zoom/yscale- yshift, RANGEMAX/zoom/yscale- yshift),
                            zmin = -3, zmax = 3,
                            xbins = nbin, ybins = nbin,
                            colormap = ColorMaps.RGBArrayMap(ra_colors), colorbar=false),
                        Plots.Command(getACString(0.0,0.0,0.0,"black","white")),
                        Plots.Command(getACString(psi,RANGEMAX*0.82/zoom/xscale - xshift,RANGEMAX*0.82/zoom/yscale- yshift,"red","black"))
                        ], xlabel="Downrange (kft)", ylabel="Crossrange (kft)", title="Training Data"))
        end

        # Plot neural network, if possible
        if nnet != nothing
            push!(g, Axis([
                        Plots.Image(get_heat2, (-1*RANGEMAX/zoom/xscale - xshift, RANGEMAX/zoom/xscale- xshift),
                            (-1*RANGEMAX/zoom/yscale- yshift, RANGEMAX/zoom/yscale- yshift),
                            zmin = -3, zmax = 3,
                            xbins = nbin, ybins = nbin,
                            colormap = ColorMaps.RGBArrayMap(ra_colors), colorbar=false),
                        Plots.Command(getACString(0.0,0.0,0.0,"black","white")),
                        Plots.Command(getACString(psi,RANGEMAX*0.82/zoom/xscale - xshift,RANGEMAX*0.82/zoom/yscale - yshift,"red","black"))
                        ], xlabel="Downrange (kft)", ylabel="Crossrange (kft)", title="Neural Network"))
        end

        # Create Color Key
        f = (x,y)->x # Dummy function for background white image
        push!(g, Axis([
            Plots.Image(f, (-2,2), (-2,2),colormap = ColorMaps.RGBArrayMap(bg_colors),colorbar=false),
            Plots.Scatter(xx, yy, zz, scatterClasses=sc),
            Plots.Node("SL ",0.25,0.915,style="black,anchor=west", axis="axis description cs"),
            Plots.Node("WL ",0.25,0.790,style="black,anchor=west", axis="axis description cs"),
            Plots.Node("COC",0.25,0.665,style="black,anchor=west", axis="axis description cs"),
            Plots.Node("WR ",0.25,0.540,style="black,anchor=west", axis="axis description cs"),
            Plots.Node("SR ",0.25,0.415,style="black,anchor=west", axis="axis description cs"),
            ],style="xshift=-1.4cm",hideAxis =true))

        # Save plot if desired. Don't return g if want to save, since returning g will cleanup and delete the png images.
        if savePlot
            PGFPlots.save("PolicyPlot.tex", g, include_preamble=true)
        else
            return g
        end
        return nothing
    end # for p_int, v0, v1, pa, ta, etc
end # function viz_policy

ext: jl

max_stars_repo_name: cherrywoods/HorizontalCAS

max_issues_repo_name: ahmed-irfan/HorizontalCAS

max_forks_repo_name: ahmed-irfan/HorizontalCAS


text:
@enum ParsingState None Delimiter EOF Newline
# at start of field: check if eof, remove leading whitespace, check if empty field
# returns `true` if result of initial parsing is a null field
# also return `b` which is last byte read
@inline function checknullstart(io::IO, opt::CSV.Options, state)
    eof(io) && (state[] = EOF; return 0x00, true)
    b = unsafe_read(io, UInt8)
    while b == CSV.SPACE || b == CSV.TAB || b == opt.quotechar
        eof(io) && (state[] = EOF; return b, true)
        b = unsafe_read(io, UInt8)
    end
    if b == opt.delim
        state[] = Delimiter
        return b, true
    elseif b == NEWLINE # check for empty field
        state[] = Newline
        return b, true
    elseif b == RETURN
        state[] = Newline
        !eof(io) && unsafe_peek(io) == NEWLINE && unsafe_read(io, UInt8)
        return b, true
    end
    return b, false
end
# check if we've successfully finished parsing a field by whether
# we've encountered a delimiter or newline break or reached eof
@inline function checkdone(io::IO, b::UInt8, opt::CSV.Options, state)
    b == opt.quotechar && !eof(io) && (b = unsafe_read(io, UInt8))
    if b == opt.delim
        state[] = Delimiter
        return b, true
    elseif b == NEWLINE
        state[] = Newline
        return b, true
    elseif b == RETURN
        state[] = Newline
        !eof(io) && unsafe_peek(io) == NEWLINE && unsafe_read(io, UInt8)
        return b, true
    elseif eof(io) && b == opt.quotechar
        state[] = EOF
        return b, true
    end
    return b, false
end

CSVError{T}(::Type{T}, b, row, col) = CSV.CSVError("error parsing a `$T` value on column $col, row $row; encountered '$(Char(b))'")

# as a last ditch effort, after we've trying parsing the correct type,
# we check if the field is equal to a custom null type
# otherwise we give up and throw an error
@inline function checknullend{T}(io::IO, ::Type{T}, b::UInt8, opt::CSV.Options, row, col, state)
    !opt.nullcheck && throw(CSVError(T, b, row, col))
    i = 1
    while true
        b == opt.null.data[i] || throw(CSVError(T, b, row, col))
        (eof(io) || i == length(opt.null)) && break
        b = unsafe_read(io, UInt8)
        i += 1
    end
    if !eof(io)
        b = unsafe_read(io, UInt8)
        b, done = checkdone(io, b, opt, state)
        done ? (return true) : throw(CSVError(T, b, row, col))
    end
    state[] = EOF
    return true
end
"""
`CSV.parsefield{T}(io::IO, ::Type{T}, opt::CSV.Options=CSV.Options(), row=0, col=0)` => `Nullable{T}`

`io` is an `IO` type that is positioned at the first byte/character of an delimited-file field (i.e. a single cell)
leading whitespace is ignored for Integer and Float types.
returns a `Nullable{T}` saying whether the field contains a null value or not (empty field, missing value)
field is null if the next delimiter or newline is encountered before any other characters.
Specialized methods exist for Integer, Float, String, Date, and DateTime.
For other types `T`, a generic fallback requires `parse(T, str::String)` to be defined.
the field value may also be wrapped in `opt.quotechar`; two consecutive `opt.quotechar` results in a null field
`opt.null` is also checked if there is a custom value provided (i.e. "NA", "\\N", etc.)
For numeric fields, if field is non-null and non-digit characters are encountered at any point before a delimiter or newline, an error is thrown
"""
function parsefield end

function parsefield{T<:Integer}(io::IO, ::Type{T}, opt::CSV.Options=CSV.Options(), row=0, col=0, state=Ref{ParsingState}(None))
    v = zero(T)
    b, null = CSV.checknullstart(io, opt, state)
    null && return Nullable{T}()
    negative = false
    if b == CSV.MINUS # check for leading '-' or '+'
        negative = true
        b = unsafe_read(io, UInt8)
    elseif b == CSV.PLUS
        b = unsafe_read(io, UInt8)
    end
    while CSV.NEG_ONE < b < CSV.TEN
        # process digits
        v *= 10
        v += b - CSV.ZERO
        eof(io) && (state[] = EOF; return Nullable{T}(ifelse(negative,-v,v)))
        b = unsafe_read(io, UInt8)
    end
    b, done::Bool = CSV.checkdone(io, b, opt, state)
    if done
        return Nullable{T}(ifelse(negative,-v,v))
    elseif CSV.checknullend(io, T, b, opt, row, col, state)
        return Nullable{T}()
    else
        throw(CSVError("couldn't parse Int"))
    end
end

const REF = Array(Ptr{UInt8},1)

function parsefield{T<:AbstractFloat}(io::IOBuffer, ::Type{T}, opt::CSV.Options=CSV.Options(), row=0, col=0, state=Ref{ParsingState}(None))
    b, null = CSV.checknullstart(io, opt, state)
    v = zero(T)
    null && return Nullable{T}()
    # subtract 1 because we just read a valid byte, so position(io) is +1 from where a digit should be
    ptr = pointer(io.data) + position(io) - 1
    v = convert(T, ccall(:jl_strtod_c, Float64, (Ptr{UInt8},Ptr{Ptr{UInt8}}), ptr, REF))
    io.ptr += REF[1] - ptr - 1 # Hopefully io.ptr doesn't change for IOBuffer?
    eof(io) && (state[] = EOF; return Nullable{T}(v))
    b = unsafe_read(io, UInt8)
    b, done = checkdone(io, b, opt, state)
    if done
        return Nullable{T}(v)
    elseif checknullend(io, T, b, opt, row, col, state)
        return Nullable{T}()
    else
        throw(CSVError("couldn't parse $T"))
    end
end

function getfloat(io, T, b, opt, row, col, buf, state)
    ptr = pointer(buf.data)
    v = convert(T, ccall(:jl_strtod_c, Float64, (Ptr{UInt8},Ptr{Ptr{UInt8}}), ptr, CSV.REF))
    if (CSV.REF[1] - ptr) != position(buf)
        reset(io)
        b, null = CSV.checknullstart(io, opt, state)
        CSV.checknullend(io, T, b, opt, row, col, state) && return Nullable{T}()
        throw(CSVError(T, b, row, col))
    end
    return Nullable{T}(v)
end

function parsefield{T<:AbstractFloat}(io::IO, ::Type{T}, opt::CSV.Options=CSV.Options(), row=0, col=0, state=Ref{ParsingState}(None))
    mark(io)
    b, null = CSV.checknullstart(io, opt, state)
    v = zero(T)
    null && return Nullable{T}()
    buf = IOBuffer(zeros(UInt8, 18), true, true)
    Base.write(buf, b)
    eof(io) && return getfloat(io, T, b, opt, row, col, buf, state)
    while true
        b = unsafe_read(io, UInt8)
        b, isdone = CSV.checkdone(io, b, opt, state)
        isdone && return getfloat(io, T, b, opt, row, col, buf, state)
        Base.write(buf, b)
        eof(io) && return getfloat(io, T, b, opt, row, col, buf, state)
    end
    throw(CSVError("couldn't parse $T"))
end

function parsefield{T<:AbstractString}(io::IOBuffer, ::Type{T}, opt::CSV.Options=CSV.Options(), row=0, col=0, state=Ref{ParsingState}(None))
    eof(io) && (state[] = EOF; return Nullable{WeakRefString{UInt8}}(WeakRefStrings.NULLSTRING, true))
    ptr = pointer(io.data) + position(io)
    len = 0
    nullcheck = opt.nullcheck # if null is "", then we don't need to byte match it
    nulllen = length(opt.null.data)
    @inbounds while !eof(io)
        b = unsafe_read(io, UInt8)
        if b == opt.quotechar
            ptr += 1
            while !eof(io)
                b = unsafe_read(io, UInt8)
                if b == opt.escapechar
                    b = unsafe_read(io, UInt8)
                    len += 1
                elseif b == opt.quotechar
                    break
                end
                (nullcheck && len+1 <= nulllen && b == opt.null.data[len+1]) || (nullcheck = false)
                len += 1
            end
        elseif b == opt.delim
            state[] = Delimiter
            break
        elseif b == b == NEWLINE
            state[] = Newline
            break
        elseif b == RETURN
            state[] = Newline
            !eof(io) && unsafe_peek(io) == NEWLINE && unsafe_read(io, UInt8)
            break
        else
            (nullcheck && len+1 <= nulllen && b == opt.null.data[len+1]) || (nullcheck = false)
            len += 1
        end
    end
    eof(io) && (state[] = EOF)
    return (len == 0 || nullcheck) ? Nullable{WeakRefString{UInt8}}(WeakRefStrings.NULLSTRING, true) : Nullable{WeakRefString{UInt8}}(WeakRefString(ptr, len))
end

function parsefield{T<:AbstractString}(io::IO, ::Type{T}, opt::CSV.Options=CSV.Options(), row=0, col=0, state=Ref{ParsingState}(None))
    eof(io) && (state[] = EOF; return Nullable{T}())
    buf = IOBuffer()
    len = 0
    nullcheck = opt.nullcheck # if null is "", then we don't need to byte match it
    nulllen = length(opt.null)
    @inbounds while !eof(io)
        b = unsafe_read(io, UInt8)
        if b == opt.quotechar
            while !eof(io)
                b = unsafe_read(io, UInt8)
                if b == opt.escapechar
                    Base.write(buf, b)
                    b = unsafe_read(io, UInt8)
                    len += 1
                elseif b == opt.quotechar
                    break
                end
                (nullcheck && len+1 <= nulllen && b == opt.null.data[len+1]) || (nullcheck = false)
                Base.write(buf, b)
                len += 1
            end
        elseif b == opt.delim
            state[] = Delimiter
            break
        elseif b == NEWLINE
            state[] = Newline
            break
        elseif b == RETURN
            state[] = Newline
            !eof(io) && unsafe_peek(io) == NEWLINE && unsafe_read(io, UInt8)
            break
        else
            (nullcheck && len+1 <= nulllen && b == opt.null.data[len+1]) || (nullcheck = false)
            Base.write(buf, b)
            len += 1
        end
    end
    eof(io) && (state[] = EOF)
    return (len == 0 || nullcheck) ? Nullable{T}() : Nullable{T}(convert(T, takebuf_string(buf)))
end

@inline itr(io,n,val) = (for i = 1:n; val *= 10; val += unsafe_read(io, UInt8) - CSV.ZERO; end; return val)

function parsefield(io::IO, ::Type{Date}, opt::CSV.Options=CSV.Options(), row=0, col=0, state=Ref{ParsingState}(None))
    b, null = CSV.checknullstart(io, opt, state)
    null && return Nullable{Date}()
    if opt.datecheck # optimize for default yyyy-mm-dd
        if CSV.NEG_ONE < b < CSV.TEN
            year = CSV.itr(io, 3, Int(b - CSV.ZERO))
            unsafe_read(io, UInt8)
            month = CSV.itr(io, 2, 0)
            unsafe_read(io, UInt8)
            day = CSV.itr(io, 2, 0)
            eof(io) && (state[] = EOF; return Nullable{Date}(Date(year, month, day)))
            b = unsafe_read(io, UInt8)
        end
        b, done = checkdone(io, b, opt, state)
        if done
            return Nullable{Date}(Date(year, month, day))
        elseif checknullend(io, Date, b, opt, row, col, state)
            return Nullable{Date}()
        else
            throw(CSVError("couldn't parse Date"))
        end
    elseif opt.dateformat == EMPTY_DATEFORMAT
        throw(ArgumentError("Can't parse a `Date` type with $EMPTY_DATEFORMAT; please provide a valid Dates.DateFormat or date format string"))
    else
        buf = IOBuffer()
        Base.write(buf, b)
        while true
            b = unsafe_read(io, UInt8)
            if b == opt.delim
                state[] = Delimiter
                break
            elseif b == CSV.NEWLINE
                state[] = Newline
                break
            elseif b == CSV.RETURN
                state[] = Newline
                !eof(io) && unsafe_peek(io) == CSV.NEWLINE && unsafe_read(io, UInt8)
                break
            elseif eof(io)
                state[] = EOF
                Base.write(buf, b)
                break
            elseif b == opt.quotechar
                b, done = checkdone(io, b, opt, state)
                break
            end
            Base.write(buf, b)
        end
        val = takebuf_string(buf)
        val == opt.null && return Nullable{Date}()
        return Nullable{Date}(Date(val, opt.dateformat)::Date)
    end
end

function parsefield(io::IO, ::Type{DateTime}, opt::CSV.Options=CSV.Options(), row=0, col=0, state=Ref{ParsingState}(None))
    b, null = CSV.checknullstart(io, opt, state)
    null && return Nullable{DateTime}()
    if opt.datecheck # optimize for default yyyy-mm-ddTHH:MM:SS
        if CSV.NEG_ONE < b < CSV.TEN
            year = CSV.itr(io, 3, Int(b - CSV.ZERO))
            unsafe_read(io, UInt8)
            month = CSV.itr(io, 2, 0)
            unsafe_read(io, UInt8)
            day = CSV.itr(io, 2, 0)
            eof(io) && (state[] = EOF; return Nullable{DateTime}(DateTime(year, month, day)))
            b = unsafe_read(io, UInt8) # read the `T`
            hour = CSV.itr(io, 2, 0)
            unsafe_read(io, UInt8)
            minute = CSV.itr(io, 2, 0)
            unsafe_read(io, UInt8)
            second = CSV.itr(io, 2, 0)
            eof(io) && (state[] = EOF; return Nullable{DateTime}(DateTime(year, month, day, hour, minute, second)))
            b = unsafe_read(io, UInt8)
        end
        b, done = checkdone(io, b, opt, state)
        if done
            return Nullable{DateTime}(DateTime(year, month, day, hour, minute, second))
        elseif checknullend(io, DateTime, b, opt, row, col, state)
            return Nullable{DateTime}()
        else
            throw(CSVError("couldn't parse DateTime"))
        end
    elseif opt.dateformat == EMPTY_DATEFORMAT
        throw(ArgumentError("Can't parse a `DateTime` type with $EMPTY_DATEFORMAT; please provide a valid Dates.DateFormat or date format string"))
    else
        buf = IOBuffer()
        Base.write(buf, b)
        while true
            b = unsafe_read(io, UInt8)
            if b == opt.delim
                state[] = Delimiter
                break
            elseif b == CSV.NEWLINE
                state[] = Newline
                break
            elseif b == CSV.RETURN
                state[] = Newline
                !eof(io) && unsafe_peek(io) == CSV.NEWLINE && unsafe_read(io, UInt8)
                break
            elseif eof(io)
                state[] = EOF
                Base.write(buf, b)
                break
            elseif b == opt.quotechar
                b, done = checkdone(io, b, opt, state)
                break
            end
            Base.write(buf, b)
        end
        val = takebuf_string(buf)
        val == opt.null && return Nullable{DateTime}()
        return Nullable{DateTime}(DateTime(val, opt.dateformat)::DateTime)
    end
end

# Generic fallback
function getgeneric(io, T, b, opt, row, col, buf, state)
    str = takebuf_string(buf)
    try
        val = parse(T, str)::T
    catch e
        reset(io)
        b, null = CSV.checknullstart(io, opt, state)
        CSV.checknullend(io, T, b, opt, row, col, state) && return Nullable{T}()
        throw(CSVError(T, b, row, col))
    end
    return Nullable{T}(val)
end

function parsefield{T}(io::IO, ::Type{T}, opt::CSV.Options=CSV.Options(), row=0, col=0, state=Ref{ParsingState}(None))
    mark(io)
    b, null = CSV.checknullstart(io, opt, state)
    null && return Nullable{T}()
    buf = IOBuffer()
    Base.write(buf, b)
    eof(io) && (state[] = EOF; return getgeneric(io, T, b, opt, row, col, buf, state))
    while true
        b = unsafe_read(io, UInt8)
        b, isdone = CSV.checkdone(io, b, opt, state)
        isdone && return getgeneric(io, T, b, opt, row, col, buf, state)
        Base.write(buf, b)
        eof(io) && return getgeneric(io, T, b, opt, row, col, buf, state)
    end
    throw(CSVError("couldn't parse $T"))
end


text:
"""
    Cell

Represents a single cell.
"""
mutable struct Cell
    mutations::Array{Int64,1}
    clonetype::Int64
end

mutable struct CloneTracker
    parenttype::Int64
    parentmodule::Int64
    time::Float64
    mutations::Array{Int64, 1}
    N0::Int64
    Ndivisions::Int64
    avdivisions::Float64
    size::Int64
end

struct Clone
    parenttype::Int64
    parentmodule::Int64
    time::Float64
    mutations::Int64
    N0::Int64
    Ndivisions::Int64
    avdivisions::Float64
    freq::Float64
    freqp::Float64
end
struct ModuleTracker
    Nvec::Array{Int64, 1}
    tvec::Array{Float64, 1}
    cells::Array{Cell, 1}
    subclones::Array{CloneTracker, 1}
    id::Int64
    parentid::Int64
end

Base.length(moduletracker::ModuleTracker) = moduletracker.Nvec[end]

abstract type SimulationInput end

struct BranchingInput <: SimulationInput
    numclones::Int64
    Nmax::Int64
    clonalmutations::Int64
    selection::Array{Float64,1}
    μ::Float64
    b::Float64
    d::Float64
    tevent::Array{Float64,1}
    fixedmu::Bool
    maxclonesize::Union{Int64, Nothing}
    ploidy::Int64
end

struct MoranInput <: SimulationInput
    N::Int64
    numclones::Int64
    tmax::Float64
    clonalmutations::Int64
    selection::Array{Float64,1}
    μ::Float64
    bdrate::Float64
    tevent::Array{Float64,1}
    fixedmu::Bool
    ploidy::Int64
end

struct BranchingMoranInput <: SimulationInput
    numclones::Int64
    Nmax::Int64
    tmax::Float64
    clonalmutations::Int64
    selection::Array{Float64,1}
    μ::Float64
    bdrate::Float64
    b::Float64
    d::Float64
    tevent::Array{Float64,1}
    fixedmu::Bool
    ploidy::Int64
end
struct MultilevelInput <: SimulationInput
    modulesize::Int64
    maxtime::Float64
    maxmodules::Int64
    clonalmutations::Int64
    μ::Float64
    bdrate::Float64
    b::Float64
    d::Float64
    fixedmu::Bool
    branchrate::Float64
    branchinitsize::Int64
    ploidy::Int64
end

struct SimulationResult
    subclones::Array{Clone, 1}
    tend::Float64
    trueVAF::Array{Float64,1}
    cells::Array{Cell, 1}
    Nvec::Array{Int64, 1}
    tvec::Array{Float64, 1}
end

struct SampledData
    VAF::Array{Float64,1}
    depth::Array{Int64,1}
end

struct Simulation{T<:SimulationInput}
    input::T
    output::ModuleTracker
end

struct VAFResult{T<:SimulationInput}
    read_depth::Float64
    cellularity::Float64
    detectionlimit::Float64
    input::T
    trueVAF::Array{Float64,1}
    sampledVAF::Array{Float64,1}
    subclonefreq::Array{Float64, 1}
    subclonefreqp::Array{Float64, 1}

end

struct MultiSimulation{T<:SimulationInput}
    input::T
    output::Array{ModuleTracker, 1}
    # mutations::Array{Int64, 1}
end

function MultiSimulation{T}(input::T) where T<:SimulationInput
    return MultiSimulation(
        input,
        ModuleTracker[]
    )
end

Base.length(multisim::MultiSimulation) = length(multisim.output)
Base.iterate(multisim::MultiSimulation) = iterate(multisim.output)
Base.iterate(multisim::MultiSimulation, state) = iterate(multisim.output, state)
Base.getindex(multisim::MultiSimulation, i) = getindex(multisim.output, i)
Base.setindex(multisim::MultiSimulation, v, i) = getindex(multisim.output, v, i)
Base.firstindex(multisim::MultiSimulation) = firstindex(multisim.output)
Base.lastindex(multisim::MultiSimulation) = lastindex(multisim.output)

const Population = MultiSimulation{MultilevelInput}

function get_simulation(multsim, i)
    return Simulation(multsim.input, multsim.output[i])
end

function BranchingInput(;numclones = 1, Nmax = 10000, ploidy = 2, μ = 10.0,
    clonalmutations = μ, selection = fill(0.0,numclones), b = log(2.0), d = 0.0,
    tevent = collect(1.0:0.5:(1+numclones)/2), fixedmu = false,
    maxclonesize = nothing)

    return BranchingInput(
            numclones,
            Nmax,
            clonalmutations,
            selection,
            μ,
            b,
            d,
            tevent,
            fixedmu,
            maxclonesize,
            ploidy
    )
end

function MoranInput(;numclones = 1, N = 10000, ploidy = 2, μ = 10.0, clonalmutations = μ,
    selection = fill(0.0,numclones), bdrate = log(2.0), tmax = 15.0,
    tevent = collect(1.0:0.5:(1+numclones)/2), fixedmu = false)

    return MoranInput(
            N,
            numclones,
            tmax,
            clonalmutations,
            selection,
            μ,
            bdrate,
            tevent,
            fixedmu,
            ploidy
    )
end

function BranchingMoranInput(;numclones = 1, Nmax = 10000, ploidy = 2, μ = 10.0,
    clonalmutations = μ, selection = fill(0.0,numclones), bdrate = log(2.0), b = log(2),
    d = 0, tmax = 15.0, tevent = collect(1.0:0.5:(1+numclones)/2), fixedmu = false)

    return BranchingMoranInput(
            numclones,
            Nmax,
            tmax,
            clonalmutations,
            selection,
            μ,
            bdrate,
            b,
            d,
            tevent,
            fixedmu,
            ploidy
    )

end

function MultilevelInput(;modulesize=200, ploidy=2, μ=10.0, clonalmutations=0,
    bdrate=log(2.0), b=log(2), d=0, maxtime=15, maxmodules=10000, fixedmu=false, branchrate=5,
    branchfraction=0.1, branchinitsize=nothing)

    return MultilevelInput(
            modulesize,
            maxtime,
            maxmodules,
            clonalmutations,
            μ,
            bdrate,
            b,
            d,
            fixedmu,
            branchrate,
            branchinitsize !== nothing ? branchinitsize : ceil(modulesize * branchfraction),
            ploidy,
    )
end

age(moduletracker::ModuleTracker) = moduletracker.tvec[end]
age(populationtracker::Vector{ModuleTracker}) = maximum(map(age, populationtracker))
age(simulation::Simulation) = age(simulation.output)
age(multisim::MultiSimulation) = maximum(age(output) for output in multisim.output)


text:
using Base.Meta

function get_names(ex::Expr)
	@assert isexpr(ex, :block)
	return [arg for arg in ex.args if !isa(arg, LineNumberNode)]
end

function code_constants(name::Symbol, ex::Expr, start::Int, exponential::Bool)
	names = get_names(ex)
	values = start .+ (0:length(names)-1)
	if (exponential)
		values = [1 << x for x in values]
	end
	table = collect(zip(names, values))
	defs = [:(const $name = $value) for (name, value) in zip(names, values)]
	return quote
		const $name = $(QuoteNode(table))
		$(defs...)
	end
end

macro constants(name::Symbol, ex::Expr)
	return esc(code_constants(name, ex, 0, false))
end

macro flags(name::Symbol, start::Int, ex::Expr)
	return esc(code_constants(name, ex, start, true))
end

max_stars_repo_name: toivoh/fmicrosynth

max_issues_repo_name: toivoh/fmicrosynth

max_forks_repo_name: toivoh/fmicrosynth


text:
@show [1 2 3; 3 2 1] .+ [2 1 2; 0 2 1]
@show [1 2 3; 2 1 2] .+ 1
@show [1 2 3; 2 2 1] .- [1 1 1; 2 1 0]
@show [1 2 1; 1 2 3] .* [3 2 1; 1 0 1]
@show [1 2 3; 3 2 1] .* 2
@show [9 8 6; 3 2 3] ./ [3 1 2; 2 1 2]
@show [3 2 2; 1 2 3] .^ [1 2 3; 2 1 2]

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
struct DescriptorSet <: LavaAbstraction
    handle::Vk.DescriptorSet
    layout::Vk.DescriptorSetLayout
end

Base.@kwdef struct GlobalDescriptorsConfig
    textures::Int = 2048
    storage_images::Int = 512
    samplers::Int = 2048
end

struct DescriptorSetBindingState{T}
    type::Vk.DescriptorType
    elements::Vector{T}
end

struct GlobalDescriptorSet
    set::DescriptorSet
    state::Dictionary{Vk.DescriptorType, Vector{Vk.DescriptorImageInfo}}
end

struct ResourceDescriptors
    pool::Vk.DescriptorPool
    gset::GlobalDescriptorSet
end

function ResourceDescriptors(device, config::GlobalDescriptorsConfig = GlobalDescriptorsConfig())
    pool = Vk.DescriptorPool(device, 1, [
        Vk.DescriptorPoolSize(Vk.DESCRIPTOR_TYPE_SAMPLED_IMAGE, 1),
        Vk.DescriptorPoolSize(Vk.DESCRIPTOR_TYPE_STORAGE_IMAGE, 1),
        Vk.DescriptorPoolSize(Vk.DESCRIPTOR_TYPE_SAMPLER, 1),
    ])

    layout = Vk.DescriptorSetLayout(device, [
        Vk.DescriptorSetLayoutBinding(0, Vk.DESCRIPTOR_TYPE_SAMPLED_IMAGE, Vk.SHADER_STAGE_ALL; descriptor_count = config.textures),
        Vk.DescriptorSetLayoutBinding(1, Vk.DESCRIPTOR_TYPE_STORAGE_IMAGE, Vk.SHADER_STAGE_ALL; descriptor_count = config.storage_images),
        Vk.DescriptorSetLayoutBinding(2, Vk.DESCRIPTOR_TYPE_SAMPLER, Vk.SHADER_STAGE_ALL; descriptor_count = config.samplers),
        Vk.DescriptorSetLayoutBinding(3, Vk.DESCRIPTOR_TYPE_COMBINED_IMAGE_SAMPLER, Vk.SHADER_STAGE_ALL; descriptor_count = config.textures),
    ], next = Vk.DescriptorSetLayoutBindingFlagsCreateInfo(repeat([Vk.DESCRIPTOR_BINDING_PARTIALLY_BOUND_BIT], 4)))

    set = DescriptorSet(first(unwrap(Vk.allocate_descriptor_sets(device, Vk.DescriptorSetAllocateInfo(pool, [layout])))), layout)
    gset = GlobalDescriptorSet(set, dictionary([
        Vk.DESCRIPTOR_TYPE_SAMPLED_IMAGE => [],
        Vk.DESCRIPTOR_TYPE_STORAGE_IMAGE => [],
        Vk.DESCRIPTOR_TYPE_SAMPLER => [],
        Vk.DESCRIPTOR_TYPE_COMBINED_IMAGE_SAMPLER => [],
    ]))
    ResourceDescriptors(pool, gset)
end


text:
#=
A permutation is an ordered arrangement of objects. For example, 3124 is one possible permutation of the digits 1, 2, 3 and 4. If all of the permutations are listed numerically or alphabetically, we call it lexicographic order. The lexicographic permutations of 0, 1 and 2 are:

012   021   102   120   201   210

What is the millionth lexicographic permutation of the digits 0, 1, 2, 3, 4, 5, 6, 7, 8 and 9?
=#

using Combinatorics
function calc()
  x = nthperm(collect(0:9),1_000_000)
  reduce(*,["$(i)" for i in x])
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
module FastbackDummyTest

using Fastback
using Dates

# create instrument
inst = Instrument("AAPL");
# inst2 = Instrument("AAPL"; data=Dict(:a => 1))

# synthetic data
N = 100_000;
prices = 1000.0 .+ 100cumsum(randn(N) .+ 0.01);
bids = prices .- 0.01;
asks = prices .+ 0.01;
dts = map(x -> DateTime(2000, 1, 1) + Minute(x) + Millisecond(123), 1:N);

# create trading account
acc = Account(10_000.0);

# plot data collectors
collect_balance, balance_curve = periodic_collector(Float64, Second(1));
collect_equity, equity_curve = periodic_collector(Float64, Second(1));
collect_open_orders, open_orders_curve = max_value_collector(Int64);
collect_drawdown, drawdown_curve = drawdown_collector(Percentage::DrawdownMode, (v, dt, equity) -> dt - v.last_dt >= Second(1));

# backtest random trading strategy
for i in 1:N
    ba = BidAsk(dts[i], bids[i], asks[i])

    if i == N
        # # close all orders at end of backtest
        # for pos in acc.open_positions
        #     execute_order!(acc, CloseOrder(pos), ba)
        # end
        execute_order!(acc, OpenOrder(inst, 100.0, Long; data=i), ba)
    else
        # randomly trade
        if rand() > 0.99 && hour(ba.dt) < 15
            pos_dir = rand() > 0.5 ? Long::TradeDir : Short::TradeDir
            pos_size = match_target_exposure(acc.equity, pos_dir, ba)
            execute_order!(acc, OpenOrder(inst, pos_size, pos_dir; data=i), ba)
        end

        # close positions after 10 minutes
        for pos in acc.open_positions
            if ba.dt - pos.open_dt >= Minute(10)
                execute_order!(acc, CloseOrder(pos), ba)
            end
        end
    end

    update_account!(acc, inst, ba)

    # collect data for analysis
    collect_balance(ba.dt, acc.balance)
    collect_equity(ba.dt, acc.equity)
    collect_open_orders(ba.dt, length(acc.open_positions))
    collect_drawdown(ba.dt, acc.equity)
end

# print account
show(acc)

# print all closed positions as pretty table
print_positions(acc.closed_positions; max_print=NaN)

# print all closed position as pretty table including custom data column renderer
print_positions(acc.closed_positions; max_print=NaN, data_renderer=(pos, data) -> data)

# # plots
# include("../backtesting/plots.jl");
# using Measures
# gr();

# title_text = "Backtest $(inst.symbol)"
# title_plot = scatter(1:2,
#     marker=0,
#     markeralpha=0,
#     annotations=(1.5, 1.5, text(title_text)),
#     foreground_color_subplot=:white,
#     axis=false,
#     grid=false,
#     leg=false);
#
#
# # plot 1
# l = @layout[
#     a{0.04h} ;  # title
#     b ;         # equity
#     c{0.3h} ;   # drawdown
#     d{0.15h} ;  # open orders
# ];
# plot(
#     title_plot,
#     plot_equity(equity_curve),
#     plot_drawdown(drawdown_curve),
#     plot_open_orders(open_orders_curve),
#     titleloc = :left,
#     titlefont = font(10),
#     margin=0mm,
#     layout=l,
#     size=(750, 1000))
#
#
# # plot 2
# l = @layout[
#     a{0.04h} ;  # title
#     b ;         # returns dist. by day
#     c           # returns dist. by hour
# ];
# plot(
#     title_plot,
#     violin_nominal_returns_by_day(acc.closed_positions),
#     violin_nominal_returns_by_hour(acc.closed_positions);
#     titleloc = :left,
#     titlefont = font(10),
#     margin=0mm,
#     layout=l,
#     size=(750, 1000))
#
#
# # plot 3
# l = @layout[
#     a{0.04h} ;  # title
#     b ;         # cum. returns by day
#     c           # cum. returns by hour
# ];
# plot(
#     title_plot,
#     plot_nominal_cum_returns_by_hour(acc.closed_positions),
#     plot_nominal_cum_returns_by_weekday(acc.closed_positions);
#     titleloc = :left,
#     titlefont = font(10),
#     margin=0mm,
#     layout=l,
#     size=(750, 1000))

end # module


text:
using NetCDF, Printf;

export read_nctiles

## read_nctiles function
#
#examples:
#  mygrid=GridSpec("LLC90")
#  fileName="nctiles_grid/GRID"
#  Depth=read_nctiles(fileName,"Depth",mygrid)
#  hFacC=read_nctiles(fileName,"hFacC",mygrid)

"""
    read_nctiles(fileName,fldName,mygrid)

Read model output from Netcdf / NCTiles file and convert to MeshArray instance.
"""
function read_nctiles(fileName::String,fldName::String,mygrid::gcmgrid)

if (mygrid.class!="llc")||(mygrid.ioSize!=[90 1170])
  error("non-llc90 cases not implemented yet");
end;

fileIn=@sprintf("%s.%04d.nc",fileName,1);
x = ncread(fileIn,fldName);
ndims=length(size(x));

#initialize f
if ndims==2;
  f0=Array{Float64}(undef,90,0);
  f00=Array{Float64}(undef,0,90);
elseif ndims==3;
  f0=Array{Float64}(undef,90,0,size(x,3));
  f00=Array{Float64}(undef,0,90,size(x,3));
elseif ndims==4;
  f0=Array{Float64}(undef,90,0,size(x,3),size(x,4));
  f00=Array{Float64}(undef,0,90,size(x,3),size(x,4));
end;
f=[f0,f0,f0,f00,f00];

#fill in f
for ff=1:13;
  #read one tile
  fileIn=@sprintf("%s.%04d.nc",fileName,ff);
  x = ncread(fileIn,fldName);
  #combine tiles
  if ff<=3;
    f[1]=cat(f[1],x;dims=2);
  elseif ff<=6;
    f[2]=cat(f[2],x;dims=2);
  elseif ff==7;
    f[3]=x;
  elseif ff<=10;
    f[4]=cat(f[4],x;dims=1);
  elseif ff<=13;
    f[5]=cat(f[5],x;dims=1);
  end;

end;

fld=MeshArray(mygrid,f)
return fld

end


text:
# HINT network from Kruse et al. (2020)
# Author: Gabrio Rizzuti, rizzuti.gabrio@gatech.edu
# Date: August 2020

using DrWatson
@quickactivate "UncertaintyQuantificationAVP"
using UncertaintyQuantificationAVP
using InvertibleNetworks
using LinearAlgebra, Test, Random
Random.seed!(11)

# Define network
nx = 1
ny = 1
n_in = 2^3
n_hidden = 64
batchsize = 16
L = 20

# Multi-scale and single scale network
H = NetworkGlow1D(n_in, batchsize, n_hidden, L; logdet=true)

###################################################################################################
# Invertibility

# Test layers
test_size = 10
X = randn(Float32, nx, ny, n_in, test_size)

# Forward-backward
Z, logdet = H.forward(X)
_, X_ = H.backward(0f0.*Z, Z)
@test isapprox(norm(X - X_)/norm(X), 0f0; atol=1f-3)

# Forward-inverse
Z, logdet = H.forward(X)
X_ = H.inverse(Z)
@test isapprox(norm(X - X_)/norm(X), 0f0; atol=1f-3)

###################################################################################################
# Gradient test

# Loss
function loss(H, X)
    Z, logdet = H.forward(X)
    f = -log_likelihood(Z) - logdet
    ΔZ = -∇log_likelihood(Z)
    ΔX = H.backward(ΔZ, Z)[1]
    return f, ΔX
end

# Gradient test w.r.t. input
H = NetworkGlow1D(n_in, batchsize, n_hidden, L; logdet=true)
X = randn(Float32, nx, ny, n_in, batchsize)
X0 = randn(Float32, nx, ny, n_in, batchsize)
dX = X - X0

f0, ΔX = loss(H, X0)
h = 0.1f0
maxiter = 6
err1 = zeros(Float32, maxiter)
err2 = zeros(Float32, maxiter)

print("\nGradient test cond. HINT net: input\n")
for j=1:maxiter
    global h
    f = loss(H, X0 + h*dX)[1]
    err1[j] = abs(f - f0)
    err2[j] = abs(f - f0 - h*dot(dX, ΔX))
    print(err1[j], "; ", err2[j], "\n")
    h = h/2f0
end

@test isapprox(err1[end] / (err1[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err2[end] / (err2[1]/4^(maxiter-1)), 1f0; atol=1f1)

# Gradient test w.r.t. parameters
H = NetworkGlow1D(n_in, batchsize, n_hidden, L; logdet=true)
X = randn(Float32, nx, ny, n_in, batchsize); H.forward(X)
H0 = NetworkGlow1D(n_in, batchsize, n_hidden, L; logdet=true)
X0 = randn(Float32, nx, ny, n_in, batchsize); H0.forward(X0)
θ0 = get_pars(H0)
θ = get_pars(H)
dθ = θ .- θ0

f0, _ = loss(H0, X0); Δθ = deepcopy(get_grads(H0))
h = 0.1f0
maxiter = 6
err1 = zeros(Float32, maxiter)
err2 = zeros(Float32, maxiter)

print("\nGradient test cond. Glow net: pars\n")
for j=1:maxiter
    global h
    set_pars!(H0, θ0+h*dθ)
    f = loss(H0, X0)[1]
    err1[j] = abs(f - f0)
    err2[j] = abs(f - f0 - h*dot(dθ, Δθ))
    print(err1[j], "; ", err2[j], "\n")
    h = h/2f0
end

@test isapprox(err1[end] / (err1[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err2[end] / (err2[1]/4^(maxiter-1)), 1f0; atol=1f1)


text:
# export FDSNevq, FDSNevt, FDSNsta
export FDSNsta

# =============================================================================
# No export

function fdsn_chp(chans::Union{String,Array{String,1},Array{String,2}}; v::Int64 = KW.v)
  # Parse channel config
  if isa(chans, String)
    C = parse_chstr(chans, fdsn = true)
  elseif isa(chans, Array{String,1})
    C = parse_charr(chans, fdsn = true)
  else
    C = copy(chans)
  end
  minreq!(C)
  v > 1 && println(stdout, "Most compact request form = ", C)
  return C
end

"""
    S = FDSNsta(chans, KW)

Retrieve station/channel info for formatted parameter file (or string) `chans`
into an empty SeisData structure.

Standard keywords: rad, reg, src, to, v

Other keywords:
* msr: get MultiStage (full) responses?
* s: Start time
* t: Termination (end) time
* xf: Name of XML file to save station metadata

See also: chanspec, parsetimewin, get_data!, SeisIO.KW
"""
function FDSNsta(chans="*"::Union{String,Array{String,1},Array{String,2}};
                  msr ::Bool              = false,          # MultiStageResp
                  rad ::Array{Float64,1}  = KW.rad,         # Search radius
                  reg ::Array{Float64,1}  = KW.reg,         # Search region
                  s   ::TimeSpec          = 0,              # Start
                  src ::String            = KW.src,         # Source server
                  t   ::TimeSpec          = (-600),         # End or Length (s)
                  to  ::Int               = KW.to,          # Read timeout (s)
                  v   ::Int64             = KW.v,           # Verbosity
                  xf  ::String            = "FDSNsta.xml"   # XML filename
                 )

  d0, d1 = parsetimewin(s, t)
  v > 0 && @info(tnote("Querying FDSN stations"))
  URL = string(fdsn_uhead(src), "station/1/query")
  BODY = "level=response\nformat=xml\n"
  wc = "*"

  # Add geographic search to BODY
  if !isempty(rad)
    BODY *= string( "latitude=", rad[1], "\n",
                    "longitude=", rad[2], "\n",
                    "minradius=", rad[3], "\n",
                    "maxradius=", rad[4], "\n")
  end
  if !isempty(reg)
    BODY *= string( "minlatitude=", reg[1], "\n",
                    "maxlatitude=", reg[2], "\n",
                    "minlongitude=", reg[3], "\n",
                    "maxlongitude=", reg[4], "\n" )
  end

  # Add channel search to BODY
  if chans == wc
    (isempty(reg) && isempty(rad)) && error("No query! Please specify a search radius, a rectangular search region, or some channels.")
    BODY *= string("* * * * ", d0, " ", d1, "\n")
  else
    C = fdsn_chp(chans, v=v)
    Nc = size(C,1)
    for i = 1:Nc
      str = ""
      for j = 1:4
        str *= (" " * (isempty(C[i,j]) ? wc : C[i,j]))
      end
      BODY *= string(str, " ", d0, " ", d1, "\n")
    end
  end
  if v > 1
    printstyled("request url:", color=:light_green)
    println(URL)
    printstyled("request body: \n", color=:light_green)
    println(BODY)
  end
  open(xf, "w") do io
    request("POST", URL, webhdr, BODY, response_stream=io)
  end

  # Build channel list
  v > 0 && @info(tnote("Building list of channels"))
  io = open(xf, "r")
  xsta = read(io, String)
  close(io)
  S = FDSN_sta_xml(xsta, msr=msr, v=v)
  return S
end

function FDSNget!(U::SeisData, chans::Union{String,Array{String,1},Array{String,2}};
  autoname  ::Bool              = false,          # Auto-generate file names?
  fmt       ::String            = KW.fmt,         # Request format
  msr       ::Bool              = false,          # MultiStageResp
  nd        ::Real              = KW.nd,          # Number of days per request
  opts      ::String            = KW.opts,        # User-defined options
  rad       ::Array{Float64,1}  = KW.rad,         # Search radius
  reg       ::Array{Float64,1}  = KW.reg,         # Search region
  s         ::TimeSpec          = 0,              # Start
  si        ::Bool              = KW.si,          # Station info?
  src       ::String            = KW.src,         # Source server
  t         ::TimeSpec          = (-600),         # End or Length (s)
  to        ::Int64             = KW.to,          # Read timeout (s)
  v         ::Int64             = KW.v,           # Verbosity
  w         ::Bool              = KW.w,           # Write to disk?
  xf        ::String            = "FDSNsta.xml",  # XML filename
  y         ::Bool              = KW.y            # Sync?
  )

  parse_err = false
  n_badreq = 0
  wc = "*"
  d0, d1 = parsetimewin(s, t)

  # (1) Time-space query for station info
  if si
    S = FDSNsta(chans,
                msr   = msr,
                rad   = rad,
                reg   = reg,
                s     = d0,
                src   = src,
                t     = d1,
                to    = to,
                v     = v,
                xf    = xf
                )
  end

  # (1a) Can we autoname the file? True iff S.n == 1
  (S.n == 1) || (autoname = false)

  # (2) Build ID strings for data query
  ID_str = Array{String,1}(undef,S.n)
  for i = 1:S.n
    ID_mat = split(S.id[i], ".")
    ID_mat[isempty.(ID_mat)] .= wc
    ID_str[i] = join(ID_mat, " ")
  end
  if v > 1
    printstyled("data query strings:\n", color=:light_green)
    for i = 1:length(ID_str)
      println(stdout, ID_str[i])
    end
  end

  # (3) Data query
  v > 0 && @info(tnote("Data query begins"))
  URL = string(fdsn_uhead(src), "dataselect/1/query")
  if occursin("ncedc", URL) || occursin("scedc", URL)
    BODY = ""
    if fmt != "miniseed"
      @warn(string("format ", fmt, " ignored; server only allows miniseed."))
    end
  else
    BODY = "format=" * fmt * "\n"
    if !isempty(opts)
      OPTS = split(opts, "&")
      for opt in OPTS
        BODY *= string(opt, "\n")
      end
    end
  end

  # Set the data source
  for i = 1:S.n
    S.src[i] = URL
  end

  # Create variables for query
  ts = tstr2int(d0)
  ti = round(Int64, nd*86400000000)
  te = tstr2int(d1)
  t1 = deepcopy(ts)
  rn = 0
  while t1 < te
    rn += 1
    os = rn > 1 ? 1 : 0
    t1 = min(ts + ti, te)
    s_str = int2tstr(ts + os)
    t_str = int2tstr(t1)
    qtail = string(" ", s_str, " ", t_str, "\n")
    QUERY = identity(BODY)
    for i = 1:S.n
      QUERY *= ID_str[i]*qtail
    end
    if v > 1
      printstyled(string("request url: ", URL, "\n"), color=:light_green)
      printstyled(string("request body: \n", QUERY), color=:light_green)
    end

    # Request via "POST"
    if w
      if fmt == "miniseed"
        ext = "mseed"
      else
        ext = fmt
      end

      # Generate filename
      yy = s_str[1:4]
      mm = s_str[6:7]
      dd = s_str[9:10]
      HH = s_str[12:13]
      MM = s_str[15:16]
      SS = s_str[18:19]
      nn = lpad(div(parse(Int64, s_str[21:26]), 1000), 3, '0')
      jj = lpad(md2j(yy, mm, dd), 3, '0')
      if autoname
        fname = join([yy, jj, HH, MM, SS, nn, S.id[1], "R", ext], '.')
      else
        fname = join([yy, jj, HH, MM, SS, nn, "FDSNWS", src, ext], '.')
      end
      safe_isfile(fname) && @warn(string("File ", fname, " contains an identical request; overwriting."))
      open(fname, "w") do io
        request("POST", URL, webhdr, QUERY, readtimeout=to, response_stream=io)
      end
      io = open(fname, "r")
      parsable = true
    else
      (R, parsable) = get_http_post(URL, QUERY, to)
      io = IOBuffer(R)
    end

    # Parse data (if we can)
    if parsable
      if fmt == "mseed" || fmt == "miniseed"
        parsemseed!(S, io, v, KW.nx_add, KW.nx_add)
      elseif fmt == "geocsv" || fmt == "geocsv.tspair"
        read_geocsv_tspair!(S, io)
      elseif fmt == "geocsv.slist"
        read_geocsv_slist!(S, io)
      else
        parse_err = true
        n_badreq += 1
        S += SeisChannel(id = string("...YYY", n_badreq),
                         misc = Dict{String,Any}( "url" => URL,
                                                  "body" => QUERY,
                                                  "data" => read(io) ) )
      end
    else
      # Should only happen with an error message (parsable to String) in io
      parse_err = true
      n_badreq += 1
      S += SeisChannel(id = string("XX...", n_badreq),
                       misc = Dict{String,Any}( "url" => URL,
                                                "body" => QUERY,
                                                "data" => readlines(io) ) )
    end
    close(io)
    ts += ti
  end

  append!(U,S)
  # Done!
  v > 0 && @info(tnote("Done FDSNget query."))
  return parse_err
end


text:
module CImGui

using CSyntax
using CSyntax.CEnum
using CSyntax.CEnum: Cenum, name_value_pairs

include("LibCImGui.jl")
using .LibCImGui

# missing macros
const IM_COL32_R_SHIFT = 0
const IM_COL32_G_SHIFT = 8
const IM_COL32_B_SHIFT = 16
const IM_COL32_A_SHIFT = 24
const IM_COL32_A_MASK  = 0xff000000
IM_COL32(R,G,B,A) = (ImU32(A)<<IM_COL32_A_SHIFT) | (ImU32(B)<<IM_COL32_B_SHIFT) | (ImU32(G)<<IM_COL32_G_SHIFT) | (ImU32(R)<<IM_COL32_R_SHIFT)
const IM_COL32_WHITE       = IM_COL32(255,255,255,255)  # opaque white = 0xffffffff
const IM_COL32_BLACK       = IM_COL32(0,0,0,255)        # opaque black
const IM_COL32_BLACK_TRANS = IM_COL32(0,0,0,0)          # transparent black = 0x00000000
const IMGUI_PAYLOAD_TYPE_COLOR_3F = "_COL3F"
const IMGUI_PAYLOAD_TYPE_COLOR_4F = "_COL4F"
const FLT_MAX = igGET_FLT_MAX()

include("helper.jl")
include("wrapper.jl")

const IMGUI_VERSION = unsafe_string(GetVersion())

include("backend/GLFW/GLFWBackend.jl")
using .GLFWBackend

include("backend/OpenGL/OpenGLBackend.jl")
using .OpenGLBackend

end # module


text:
using DataFrames
using LinearAlgebra
using MixedModels
using NamedArrays
using Random
using SparseArrays
using Statistics
using StatsModels
using Tables
using Test

using MixedModels: dataset, likelihoodratiotest

const fms = Dict(
    :dyestuff => [@formula(yield ~ 1 + (1|batch))],
    :dyestuff2 => [@formula(yield ~ 1 + (1|batch))],
    :d3 => [@formula(y ~ 1 + u + (1+u|g) + (1+u|h) + (1+u|i))],
    :insteval => [
        @formula(y ~ 1 + service + (1|s) + (1|d) + (1|dept)),
        @formula(y ~ 1 + service*dept + (1|s) + (1|d)),
    ],
    :kb07 => [
        @formula(rt_trunc ~ 1+spkr+prec+load+(1|subj)+(1|item)),
        @formula(rt_trunc ~ 1+spkr*prec*load+(1|subj)+(1+prec|item)),
        @formula(rt_trunc ~ 1+spkr*prec*load+(1+spkr+prec+load|subj)+(1+spkr+prec+load|item)),
    ],
    :pastes => [
        @formula(strength ~ 1 + (1|sample)),
        @formula(strength ~ 1 + (1|sample) + (1|batch)),
    ],
    :penicillin => [@formula(diameter ~ 1 + (1|plate) + (1|sample))],
    :sleepstudy => [
        @formula(reaction ~ 1 + days + (1|subj)),
        @formula(reaction ~ 1 + days + zerocorr(1+days|subj)),
        @formula(reaction ~ 1 + days + (1|subj) + (0+days|subj)),
        @formula(reaction ~ 1 + days + (1+days|subj)),
    ],
)

const fittedmodels = Dict{Symbol,Vector{LinearMixedModel}}();
function models(nm::Symbol)
    get!(fittedmodels, nm) do
        fit.(MixedModel, fms[nm], Ref(dataset(nm)))
    end
end

const io = IOBuffer()

@testset "Dyestuff" begin
    fm1 = only(models(:dyestuff))

    @test length(fm1.allterms) == 3
    @test size(fm1.reterms) == (1, )
    @test lowerbd(fm1) == zeros(1)
    @test fm1.lowerbd == zeros(1)
    @test fm1.optsum.initial == ones(1)
    fm1.θ = ones(1)
    @test fm1.θ == ones(1)
    fm1.optsum.feval = -1
    @test_logs (:warn, "Model has not been fit") show(fm1)

    @test objective(updateL!(setθ!(fm1, [0.713]))) ≈ 327.34216280955366
    @test_deprecated MixedModels.describeblocks(IOBuffer(), fm1)

    show(io, BlockDescription(fm1))
    @test countlines(seekstart(io)) == 3
    output = String(take!(io))
    @test startswith(output, "rows:")

    fit!(fm1);
    @test :θ in propertynames(fm1)
    @test objective(fm1) ≈ 327.3270598811428 atol=0.001
    @test fm1.θ ≈ [0.752580] atol=1.e-5
    @test fm1.λ ≈ [LowerTriangular(reshape(fm1.θ, (1,1)))] atol=1.e-5
    @test deviance(fm1) ≈ 327.32705988 atol=0.001
    @test aic(fm1) ≈ 333.3270598811394 atol=0.001
    @test bic(fm1) ≈ 337.5306520261259 atol=0.001
    @test fixef(fm1) ≈ [1527.5]
    @test dispersion_parameter(fm1)
    @test first(first(fm1.σs)) ≈ 37.26034462135931 atol=0.0001
    @test fm1.β ≈ [1527.5]
    @test dof(fm1) == 3
    @test nobs(fm1) == 30
    @test MixedModels.fixef!(zeros(1),fm1) ≈ [1527.5]
    @test coef(fm1) ≈ [1527.5]
    fm1β = fm1.βs
    @test fm1β isa NamedTuple
    @test isone(length(fm1β))
    @test first(values(fm1β)) ≈ 1527.5
    fm1σρ = fm1.σρs
    @test fm1σρ isa NamedTuple
    @test isone(length(fm1σρ))
    @test isone(length(getproperty(first(fm1σρ), :σ)))
    @test isempty(getproperty(first(fm1σρ), :ρ))
    @test fm1.σ == sdest(fm1)
    @test fm1.b == ranef(fm1)
    @test fm1.u == ranef(fm1, uscale=true)
    @test fm1.stderror == stderror(fm1)
    @test isone(length(fm1.pvalues))
    @test fm1.objective == objective(fm1)
    @test fm1.σ ≈ 49.510099986291145 atol=1.e-5
    @test fm1.X == ones(30,1)
    ds = MixedModels.dataset(:dyestuff)
    @test fm1.y == ds[!, :yield]
    @test cond(fm1) == ones(1)
    @test first(leverage(fm1)) ≈ 0.15650534392640486 rtol=1.e-5
    @test sum(leverage(fm1)) ≈ 4.695160317792145 rtol=1.e-5
    cm = coeftable(fm1)
    @test length(cm.rownms) == 1
    @test length(cm.colnms) == 4
    @test fnames(fm1) == (:batch,)
    @test response(fm1) == ds[!, :yield]
    rfu = ranef(fm1, uscale = true)
    rfb = ranef(fm1)
    @test abs(sum(rfu[1])) < 1.e-5
    cv = condVar(fm1)
    @test length(cv) == 1
    @test size(first(cv)) == (1, 1, 6)
    show(IOBuffer(), fm1.optsum)

    @test logdet(fm1) ≈ 8.06014522999825 atol=0.001
    @test varest(fm1) ≈ 2451.2501089607676 atol=0.001
    @test pwrss(fm1) ≈ 73537.49947885796 atol=0.001
    @test stderror(fm1) ≈ [17.69455188898009] atol=0.0001

    vc = VarCorr(fm1)
    show(io, vc)
    str = String(take!(io))
    @test startswith(str, "Variance components:")
    @test vc.s == sdest(fm1)

    fit!(fm1, REML=true)
    @test objective(fm1) ≈ 319.65427684225216 atol=0.0001
    @test_throws ArgumentError loglikelihood(fm1)
    @test dof_residual(fm1) ≥ 0

    print(io, fm1)
    @test startswith(String(take!(io)), "Linear mixed model fit by REML")

    fm1.optsum.maxfeval = 5
    @test_logs (:warn, "NLopt optimization failure: MAXEVAL_REACHED") fit!(fm1)
    fm1.optsum.maxfeval = -1

    vc = fm1.vcov
    @test isa(vc, Matrix{Float64})
    @test only(vc) ≈ 409.79495436473167 rtol=1.e-6
end

@testset "Dyestuff2" begin
    ds2 = MixedModels.dataset(:dyestuff2)
    fm = only(models(:dyestuff2))
    @test lowerbd(fm) == zeros(1)
    show(IOBuffer(), fm)
    @test fm.θ ≈ zeros(1)
    @test objective(fm) ≈ 162.87303665382575
    @test abs(std(fm)[1][1]) < 1.0e-9
    @test std(fm)[2] ≈ [3.653231351374652]
    @test stderror(fm) ≈ [0.6669857396443261]
    @test coef(fm) ≈ [5.6656]
    @test logdet(fm) ≈ 0.0
    @test issingular(fm)
    refit!(fm, float(MixedModels.dataset(:dyestuff)[!, :yield]))
    @test objective(fm) ≈ 327.3270598811428 atol=0.001
end

@testset "penicillin" begin
    fm = only(models(:penicillin))
    @test size(fm) == (144, 1, 30, 2)
    @test fm.optsum.initial == ones(2)
    @test lowerbd(fm) == zeros(2)

    @test objective(fm) ≈ 332.18834867227616 atol=0.001
    @test coef(fm) ≈ [22.97222222222222] atol=0.001
    @test fixef(fm) ≈ [22.97222222222222] atol=0.001
    @test coef(fm)[1] ≈ mean(dataset(:penicillin).diameter)
    @test stderror(fm) ≈ [0.7445960346851368] atol=0.0001
    @test fm.θ ≈ [1.5375772376554968, 3.219751321180035] atol=0.001
    @test first(std(fm)) ≈ [0.8455645948223015] atol=0.0001
    @test std(fm)[2] ≈ [1.770647779277388] atol=0.0001
    @test varest(fm) ≈ 0.3024263987592062 atol=0.0001
    @test logdet(fm) ≈ 95.74614821367786 atol=0.001

    @test_throws ArgumentError condVar(fm)

    rfu = ranef(fm, uscale=true)
    @test length(rfu) == 2
    @test first(first(rfu)) ≈ 0.523162392717432 rtol=1.e-4

    rfb = ranef(fm)
    @test length(rfb) == 2
    @test last(last(rfb)) ≈ -3.001823834230942 rtol=1.e-4

    show(io, BlockDescription(fm))
    @test countlines(seekstart(io)) == 4
    tokens = Set(split(String(take!(io)), r"\s+"))
    @test "Diag/Dense" in tokens
    @test "Diagonal" in tokens
end

@testset "pastes" begin
    fm = last(models(:pastes))
    @test size(fm) == (60, 1, 40, 2)
    @test fm.optsum.initial == ones(2)
    @test lowerbd(fm) == zeros(2)

    @test objective(fm) ≈ 247.99446586289676 atol=0.001
    @test coef(fm) ≈ [60.05333333333329] atol=0.001
    @test fixef(fm) ≈ [60.05333333333329] atol=0.001
    @test stderror(fm) ≈ [0.6421359883527029] atol=0.0001
    @test fm.θ ≈ [3.5268858714382905, 1.3299230213750168] atol=0.001
    @test first(std(fm)) ≈ [2.904069002535747] atol=0.001
    @test std(fm)[2] ≈ [1.095070371687089] atol=0.0001
    @test std(fm)[3] ≈ [0.8234088395243269] atol=0.0001
    @test varest(fm) ≈ 0.6780020742644107 atol=0.0001
    @test logdet(fm) ≈ 101.0381339953986 atol=0.001

    show(io, BlockDescription(fm))
    @test countlines(seekstart(io)) == 4
    tokens = Set(split(String(take!(io)), r"\s+"))
    @test "Sparse" in tokens
    @test "Diagonal" in tokens

    lrt = likelihoodratiotest(models(:pastes)...)
    @test length(lrt.deviance) == length(lrt.formulas) == length(lrt.models )== 2
    @test first(lrt.tests.pvalues) ≈ 0.5233767966395597 atol=0.0001
end

@testset "InstEval" begin
    fm1 = first(models(:insteval))
    @test size(fm1) == (73421, 2, 4114, 3)
    @test fm1.optsum.initial == ones(3)
    @test lowerbd(fm1) == zeros(3)

    @test objective(fm1) ≈ 237721.7687745563 atol=0.001
    ftd1 = fitted(fm1);
    @test size(ftd1) == (73421, )
    @test ftd1 == predict(fm1)
    @test first(ftd1) ≈ 3.17876 atol=0.0001
    resid1 = residuals(fm1);
    @test size(resid1) == (73421, )
    @test first(resid1) ≈ 1.82124 atol=0.00001

    @testset "PCA" begin
        @test length(fm1.rePCA) == 3
        pca = MixedModels.PCA(fm1)
        @test length(pca) == 3
        @test :covcor in propertynames(first(pca))
        str = String(take!(io))
        show(io, first(pca), stddevs=true, variances=true)
        str = String(take!(io))
        @test !isempty(findall("Standard deviations:", str))
        @test !isempty(findall("Variances:", str))
    end

    show(io, BlockDescription(fm1))
    @test countlines(seekstart(io)) == 5
    tokens = Set(split(String(take!(io)), r"\s+"))
    @test "Sparse" in tokens
    @test "Sparse/Dense" in tokens
    @test "Diag/Dense" in tokens

    fm2 = last(models(:insteval))
    @test objective(fm2) ≈ 237585.5534151694 atol=0.001
    @test size(fm2) == (73421, 28, 4100, 2)
end

@testset "sleep" begin
    fm = last(models(:sleepstudy))
    @test lowerbd(fm) == [0.0, -Inf, 0.0]
    A11 = fm.A[Block(1,1)]
    @test isa(A11, UniformBlockDiagonal{Float64})
    @test isa(fm.L[Block(1, 1)], UniformBlockDiagonal{Float64})
    @test size(A11) == (36, 36)
    a11 = view(A11.data, :, :, 1)
    @test a11 == [10. 45.; 45. 285.]
    @test size(A11.data, 3) == 18
    λ = first(fm.λ)
    b11 = LowerTriangular(view(fm.L[Block(1, 1)].data, :, :, 1))
    @test b11 * b11' ≈ λ'a11*λ + I rtol=1e-5
    @test count(!iszero, Matrix(fm.L[Block(1, 1)])) == 18 * 4
    @test rank(fm) == 2

    @test objective(fm) ≈ 1751.9393444647046
    @test fm.θ ≈ [0.929221307, 0.01816838, 0.22264487096] atol=1.e-6
    @test pwrss(fm) ≈ 117889.46144025437
    @test logdet(fm) ≈ 73.90322021999222 atol=0.001
    @test stderror(fm) ≈ [6.632257721914501, 1.5022354739749826] atol=0.0001
    @test coef(fm) ≈ [251.40510484848477,10.4672859595959]
    @test fixef(fm) ≈  [251.40510484848477,10.4672859595959]
    @test std(fm)[1] ≈ [23.780468100188497, 5.716827903196682] atol=0.01
    @test logdet(fm) ≈ 73.90337187545992 atol=0.001
    @test cond(fm) ≈ [4.175251] atol=0.0001
    @test loglikelihood(fm) ≈ -875.9696722323523
    @test sum(leverage(fm)) ≈ 28.611525700136877 rtol=1.e-5
    σs = fm.σs
    @test length(σs) == 1
    @test keys(σs) == (:subj,)
    @test length(σs.subj) == 2
    @test first(values(σs.subj)) ≈ 23.780468626896745 atol=0.0001
    @test last(values(first(σs))) ≈ 5.716827808126002 atol=0.0001
    @test fm.corr ≈ [1.0 -0.1375451787621904; -0.1375451787621904 1.0] atol=0.0001

    u3 = ranef(fm, uscale=true)
    @test length(u3) == 1
    @test size(first(u3)) == (2, 18)
    @test first(u3)[1, 1] ≈ 3.030300122575336 atol=0.001
    u3n = first(ranef(fm, uscale=true, named=true))
    @test u3n isa NamedArray

    b3 = ranef(fm)
    @test length(b3) == 1
    @test size(first(b3)) == (2, 18)
    @test first(first(b3)) ≈ 2.815819441982976 atol=0.001

    simulate!(fm)  # to test one of the unscaledre methods

    fmnc = zerocorr!(deepcopy(fm))
    @test fmnc.optsum.feval < 0
    @test size(fmnc) == (180,2,36,1)
    @test fmnc.θ == [fm.θ[1], fm.θ[3]]
    @test lowerbd(fmnc) == zeros(2)

    fmnc = models(:sleepstudy)[2]
    @test size(fmnc) == (180,2,36,1)
    @test fmnc.optsum.initial == ones(2)
    @test lowerbd(fmnc) == zeros(2)

    @testset "zerocorr PCA" begin
        @test length(fmnc.rePCA) == 1
        @test fmnc.rePCA.subj ≈ [0.5, 1.0]
        @test any(Ref(fmnc.PCA.subj.loadings) .≈ (I(2), I(2)[:, [2,1]]))
        @test show(IOBuffer(), MixedModels.PCA(fmnc)) === nothing
    end

    @test deviance(fmnc) ≈ 1752.0032551398835 atol=0.001
    @test objective(fmnc) ≈ 1752.0032551398835 atol=0.001
    @test coef(fmnc) ≈ [251.40510484848585, 10.467285959595715]
    @test fixef(fmnc) ≈ [251.40510484848477, 10.467285959595715]
    @test stderror(fmnc) ≈ [6.707710260366577, 1.5193083237479683] atol=0.001
    @test fmnc.θ ≈ [0.9458106880922268, 0.22692826607677266] atol=0.0001
    @test first(std(fmnc)) ≈ [24.171449463289047, 5.799379721123582]
    @test last(std(fmnc)) ≈ [25.556130034081047]
    @test logdet(fmnc) ≈ 74.46952585564611 atol=0.001
    ρ = first(fmnc.σρs.subj.ρ)
    @test ρ === -0.0   # test that systematic zero correlations are returned as -0.0

    MixedModels.likelihoodratiotest(fm, fmnc)
    slp = MixedModels.dataset(:sleepstudy)
    fmrs = fit(MixedModel, @formula(reaction ~ 1+days + (0+days|subj)), slp);
    @test objective(fmrs) ≈ 1774.080315280528 rtol=0.00001
    @test fmrs.θ ≈ [0.24353985679033105] rtol=0.00001

    fm_ind = fit(MixedModel, @formula(reaction ~ 1+days + (1|subj) + (0+days|subj)), slp)
    @test objective(fm_ind) ≈ objective(fmnc)
    @test coef(fm_ind) ≈ coef(fmnc)
    @test fixef(fm_ind) ≈ fixef(fmnc)
    @test stderror(fm_ind) ≈ stderror(fmnc)
    @test fm_ind.θ ≈ fmnc.θ
    @test std(fm_ind) ≈ std(fmnc)
    @test logdet(fm_ind) ≈ logdet(fmnc)

    # combining [ReMat{T,S1}, ReMat{T,S2}] for S1 ≠ S2
    slpcat = categorical!(deepcopy(slp), [:days])
    fm_cat = fit(MixedModel, @formula(reaction ~ 1+days+(1|subj)+(0+days|subj)),slpcat)
    @test fm_cat isa LinearMixedModel
    σρ = fm_cat.σρs
    @test σρ isa NamedTuple
    @test isone(length(σρ))
    @test first(keys(σρ)) == :subj
    @test keys(σρ.subj) == (:σ, :ρ)
    @test length(σρ.subj) == 2
    @test length(first(σρ.subj)) == 11
    @test length(σρ.subj.ρ) == 55
    @test iszero(σρ.subj.ρ[46])
    @test σρ.subj.ρ[46] === -0.0

    show(io, BlockDescription(first(models(:sleepstudy))))
    @test countlines(seekstart(io)) == 3
    @test "Diagonal" in Set(split(String(take!(io)), r"\s+"))

    show(io, BlockDescription(last(models(:sleepstudy))))
    @test countlines(seekstart(io)) == 3
    @test "BlkDiag" in Set(split(String(take!(io)), r"\s+"))

end

@testset "d3" begin
    fm = only(models(:d3))
    @test pwrss(fm) ≈ 5.30480294295329e6 rtol=1.e-4
    @test objective(fm) ≈ 884957.5540213 rtol = 1e-4
    @test coef(fm) ≈ [0.4991229873, 0.31130780953] atol = 1.e-4
    @test length(ranef(fm)) == 3

    show(io, BlockDescription(fm))
    tokens = Set(split(String(take!(io)), r"\s+"))
    @test "BlkDiag/Dense" in tokens
end

@testset "kb07" begin
    pca = last(models(:kb07)).PCA
    @test keys(pca) == (:subj, :item)
    show(io, models(:kb07)[2])
    tokens = Set(split(String(take!(io)), r"\s+"))
    @test "Corr." in tokens
    @test "-0.89" in tokens
end

@testset "simulate!" begin
    ds = dataset(:dyestuff)
    fm = only(models(:dyestuff))
    resp₀ = copy(response(fm))
    # type conversion of ints to floats
    simulate!(Random.MersenneTwister(1234321), fm, β=[1], σ=1)
    refit!(fm,resp₀)
    refit!(simulate!(Random.MersenneTwister(1234321), fm))
    @test deviance(fm) ≈ 339.0218639362958 atol=0.001
    refit!(fm, float(ds.yield))
    Random.seed!(1234321)
    refit!(simulate!(fm))
    @test deviance(fm) ≈ 339.0218639362958 atol=0.001
    simulate!(fm, θ = fm.θ)
    @test_throws DimensionMismatch refit!(fm, zeros(29))

    # two implicit tests
    # 1. type conversion of ints to floats
    # 2. test method for default RNG
    parametricbootstrap(1, fm, β=[1], σ=1)

    bsamp = parametricbootstrap(MersenneTwister(1234321), 100, fm, use_threads=false)
    @test isa(propertynames(bsamp), Vector{Symbol})
    @test length(bsamp.objective) == 100
    @test keys(first(bsamp.bstr)) == (:objective, :σ, :β, :se, :θ)
    @test isa(bsamp.σs, Vector{<:NamedTuple})
    @test length(bsamp.σs) == 100
    allpars = DataFrame(bsamp.allpars)
    @test isa(allpars, DataFrame)
    cov = shortestcovint(shuffle(1.:100.))
    # there is no unique shortest coverage interval here, but the left-most one
    # is currently returned, so we take that. If this behavior changes, then
    # we'll have to change the test
    @test first(cov) == 1.
    @test last(cov) == 95.

    bsamp_threaded = parametricbootstrap(MersenneTwister(1234321), 100, fm, use_threads=true)
    # even though it's bad practice with floating point, exact equality should
    # be a valid test here -- if everything is working right, then it's the exact
    # same operations occuring within each bootstrap sample, which IEEE 754
    # guarantees will yield the same result
    @test sort(bsamp_threaded.σ) == sort(bsamp.σ)
    @test sort(bsamp_threaded.θ) == sort(bsamp.θ)
    @test sort(columntable(bsamp_threaded.β).β) == sort(columntable(bsamp.β).β)
    @test sum(issingular(bsamp)) == sum(issingular(bsamp_threaded))
end

@testset "Rank deficient" begin
    rng = MersenneTwister(0);
    x = rand(rng, 100);
    data = (x = x, x2 = 1.5 .* x, y = rand(rng, 100), z = repeat('A':'T', 5))
    model = fit(MixedModel, @formula(y ~ x + x2 + (1|z)), data)
    @test length(fixef(model)) == 2
    @test rank(model) == 2
    @test length(coef(model)) == 3
    ct = coeftable(model)
    @test ct.rownms ==  ["(Intercept)", "x", "x2"]
    @test fixefnames(model) == ["(Intercept)", "x"]
    @test coefnames(model) == ["(Intercept)", "x", "x2"]
    @test last(coef(model)) == -0.0
    stde = MixedModels.stderror!(zeros(3), model)
    @test isnan(last(stde))
    vc = vcov(model)
    @test isnan(last(vc))

    @test size(StatsModels.modelmatrix(model), 2) == 3
    # check preserving of name ordering in coeftable and placement of
    # pivoted-out element
    fill!(data.x, 0)
    model2 = fit(MixedModel, @formula(y ~ x + x2 + (1|z)), data)
    ct = coeftable(model2)
    @test ct.rownms ==  ["(Intercept)", "x", "x2"]
    @test fixefnames(model2) == ["(Intercept)", "x2"]
    @test coefnames(model2) == ["(Intercept)", "x", "x2"]
    @test coef(model2)[2] == -0.0
    @test last(fixef(model)) ≈ (last(fixef(model2)) * 1.5)
    stde2 = MixedModels.stderror!(zeros(3), model2)
    @test isnan(stde2[2])

end

@testset "coeftable" begin
    ct = coeftable(only(models(:dyestuff)));
    @test [3,4] == [ct.teststatcol, ct.pvalcol]
end


text:
import Pkg; Pkg.activate(joinpath(@__DIR__, "../../"))  # CytofRepFam
include(joinpath(@__DIR__, "ParTemp.jl"))
using Distributions
using LinearAlgebra
import Random


text:
using AutoHashEquals
abstract type SymbolicExpr end
import Base: show

struct One <: SymbolicExpr; end
show(io::IO, o::One) = print(io, 1)

struct NegOne <: SymbolicExpr; end
show(io::IO, o::NegOne) = print(io, -1)
Base.:-(a::SymbolicExpr) = *(NegOne(), a)
Base.:-(a::SymbolicExpr, b::SymbolicExpr) = +(a,-(b))

struct Var <: SymbolicExpr
    name::Symbol
end
show(io::IO, v::Var) = print(io, v.name)

const superscripts = OffsetArray(['⁰','¹','²','³','⁴','⁵','⁶','⁷','⁸','⁹'], 0:9)
struct Pow{T<:SymbolicExpr} <: SymbolicExpr
    s::T
    n::Int
    Pow(p::Pow, n::Integer) = new{typeof(p.s)}(p.s, p.n*n)
    Pow(s, n::Integer) = new{typeof(s)}(s, n)
end
Base.:^(s::SymbolicExpr, n::Integer) = Pow(s, n)
show(io::IO, p::Pow) = print(io, p.s, reverse(map(x->superscripts[x], digits(p.n)))...)

with_parens(x, T) = isa(x, T) ? sprint(print, "(", x, ")") : sprint(print, x)

@auto_hash_equals struct Plus <: SymbolicExpr
    x::NTuple{N, SymbolicExpr} where N
    function Plus(x::NTuple{N, SymbolicExpr}) where {N}
        es = SymbolicExpr[]
        for e in x
            isa(e, Plus) ? append!(es, e.x) : push!(es, e)
        end
        new((es...,))
    end
end
Plus(x::NTuple{1, SymbolicExpr}) = x[1]
function show(io::IO, p::Plus)
    length(p.x) == 0 && return print(io, "0")
    print(io, p.x[1])
    for e in p.x[2:end]
        if isa(e, Times) && e.x[1] == NegOne()
            print(io, " - ")
            print(io, with_parens(Times(e.x[2:end]), Plus))
        else
            print(io, " + ")
            print(io, with_parens(e, Plus))
        end
    end
end
Base.:+(x::SymbolicExpr...) = Plus(x)
Base.zero(::Type{SymbolicExpr}) = Plus(())

@auto_hash_equals struct Times <: SymbolicExpr
    x::NTuple{N, SymbolicExpr} where {N}
    function Times(x::NTuple{N, SymbolicExpr}) where {N}
        es = SymbolicExpr[]
        parity = One()
        function process!(x)
            for e in x
                isa(e, One) && continue
                if isa(e, NegOne)
                    parity = (parity == One() ? NegOne() : One())
                    continue
                end
                isa(e, Times) ? process!(e.x) : push!(es, e)
            end
        end
        process!(x)
        parity == One() || pushfirst!(es, NegOne())
        new((es...,))
    end
end
Times(x::NTuple{1, SymbolicExpr}) = x[1]
show(io::IO, p::Times) = join(io, map(x->with_parens(x, Plus), p.x), " * ")
Base.:*(x::SymbolicExpr...) = Times(x)

const subscripts = OffsetArray(['₀' + i for i = 0:9], 0:9)
struct RootOfUnity <: SymbolicExpr
    n
end

function simplify(p::Pow{RootOfUnity})
    (p.n == 0 || p.n == p.s.n) && return One()
    p.n == 1 && return p.s
    p.n > p.s.n && return simplify(Pow(p.s, mod(p.n, p.s.n)))
    iseven(p.s.n) && p.n >= div(p.s.n, 2) &&
        return simplify(Times((NegOne(), Pow(p.s, mod(p.n, div(p.s.n, 2))))))
    return p
end

function Base.:(==)(a::Union{RootOfUnity, Pow{RootOfUnity}},
         b::Union{RootOfUnity, Pow{RootOfUnity}})
    base = lcm((isa(x, RootOfUnity) ? x.n : x.s.n for x in (a,b))...)
    _d(r::RootOfUnity) = div(base, r.n)
    _d(p::Pow{RootOfUnity}) = _d(p.s) * p.n
    _d(a) == _d(b)
end

function accum(itr)
    base = lcm((isa(x, RootOfUnity) ? x.n : x.s.n for x in itr)...)
    _d(r::RootOfUnity) = div(base, r.n)
    _d(p::Pow{RootOfUnity}) = _d(p.s) * p.n
    n = sum(_d, itr)
    simplify(Pow(RootOfUnity(base), n))
end

function simplify(t::Times)
    t = Times(map(simplify, t.x))
    isa(t, Times) || return t
    for (i,e) in pairs(t.x)
        isa(e, One) && return simplify(Times(t.x[filter(j->j != i, 1:length(t.x))]))
    end
    roots = filter(x->isa(x, Union{RootOfUnity, Pow{RootOfUnity}}), collect(t.x))
    if length(roots) > 1
        others = filter(x->!isa(x, Union{RootOfUnity, Pow{RootOfUnity}}), collect(t.x))
        return Times((others..., accum(roots)))
    end
    t
end

find_var(x) = nothing
find_var(x::Var) = x.name
find_var(t::Times) = something(map(find_var, t.x)...)
find_var(t::Plus) = something(map(find_var, t.x)...)

cancels(a::Times, b) = isa(a.x[1], NegOne) && Times(a.x[2:end]) == b

function simplify(t::Plus)
    t = Plus(map(simplify, t.x))
    # Sort the terms by variable name
    t = Plus(tuple(sort(collect(t.x), by=find_var)...))
    # See if adjacent terms cancel each other out
    for i = 1:length(t.x)-1
        term = t.x[i]
        next_term = t.x[i+1]
        if (isa(term, Times) && cancels(term, next_term)) ||
            (isa(next_term, Times) && cancels(next_term, term))
            return simplify(Plus(t.x[filter(j->!in(j, (i, i+1)), 1:length(t.x))]))
        end
    end
    return t
end
simplify(x) = x

expand(x) = x
expand(p::Plus) = Plus(map(expand, p.x))
function expand(t::Times)
    for (i,e) in pairs(t.x)
        if isa(e, Plus)
            other_factors = Times(t.x[filter(j->j != i, 1:length(t.x))])
            return Plus(((expand(Times((pe, other_factors))) for pe in e.x)...,))
        end
    end
    return t
end

show(io::IO, ω::RootOfUnity) = print(io, string("ω", subscripts[ω.n]))


text:
#! format: off
########################### ElectricLoad ####################################

get_variable_multiplier(_, ::Type{<:PSY.ElectricLoad}, ::AbstractLoadFormulation) = -1.0

########################### ActivePowerVariable, ElectricLoad ####################################

get_variable_binary(::ActivePowerVariable, ::Type{<:PSY.ElectricLoad}, ::AbstractLoadFormulation) = false
get_variable_lower_bound(::ActivePowerVariable, d::PSY.ElectricLoad, ::AbstractLoadFormulation) = 0.0
get_variable_upper_bound(::ActivePowerVariable, d::PSY.ElectricLoad, ::AbstractLoadFormulation) = PSY.get_active_power(d)

########################### ReactivePowerVariable, ElectricLoad ####################################

get_variable_binary(::ReactivePowerVariable, ::Type{<:PSY.ElectricLoad}, ::AbstractLoadFormulation) = false

get_variable_lower_bound(::ReactivePowerVariable, d::PSY.ElectricLoad, ::AbstractLoadFormulation) = 0.0
get_variable_upper_bound(::ReactivePowerVariable, d::PSY.ElectricLoad, ::AbstractLoadFormulation) = PSY.get_reactive_power(d)

########################### ReactivePowerVariable, ElectricLoad ####################################

get_variable_binary(::OnVariable, ::Type{<:PSY.ElectricLoad}, ::AbstractLoadFormulation) = true

get_multiplier_value(::TimeSeriesParameter, d::PSY.ElectricLoad, ::StaticPowerLoad) = -1*PSY.get_max_active_power(d)
get_multiplier_value(::ReactivePowerTimeSeriesParameter, d::PSY.ElectricLoad, ::StaticPowerLoad) = -1*PSY.get_max_reactive_power(d)
get_multiplier_value(::TimeSeriesParameter, d::PSY.ElectricLoad, ::AbstractControllablePowerLoadFormulation) = PSY.get_max_active_power(d)

########################Objective Function##################################################
proportional_cost(cost::Nothing, ::OnVariable, ::PSY.ElectricLoad, ::AbstractControllablePowerLoadFormulation)=1.0
proportional_cost(cost::PSY.OperationalCost, ::OnVariable, ::PSY.ElectricLoad, ::AbstractControllablePowerLoadFormulation)=PSY.get_fixed(cost)

objective_function_multiplier(::VariableType, ::AbstractControllablePowerLoadFormulation)=OBJECTIVE_FUNCTION_NEGATIVE

variable_cost(::Nothing, ::PSY.ElectricLoad, ::ActivePowerVariable, ::AbstractControllablePowerLoadFormulation)=1.0
variable_cost(cost::PSY.OperationalCost, ::ActivePowerVariable, ::PSY.ElectricLoad, ::AbstractControllablePowerLoadFormulation)=PSY.get_variable(cost)

#! format: on

function get_default_time_series_names(
    ::Type{<:PSY.ElectricLoad},
    ::Type{<:Union{FixedOutput, AbstractLoadFormulation}},
)
    return Dict{Type{<:TimeSeriesParameter}, String}(
        ActivePowerTimeSeriesParameter => "max_active_power",
        ReactivePowerTimeSeriesParameter => "max_active_power",
    )
end

function get_default_attributes(
    ::Type{U},
    ::Type{V},
) where {U <: PSY.ElectricLoad, V <: Union{FixedOutput, AbstractLoadFormulation}}
    return Dict{String, Any}()
end

get_initial_conditions_device_model(
    ::OperationModel,
    ::DeviceModel{T, <:AbstractLoadFormulation},
) where {T <: PSY.ElectricLoad} = DeviceModel(T, StaticPowerLoad)

####################################### Reactive Power Constraints #########################
"""
Reactive Power Constraints on Controllable Loads Assume Constant power_factor
"""
function add_constraints!(
    container::OptimizationContainer,
    T::Type{<:ReactivePowerVariableLimitsConstraint},
    U::Type{<:ReactivePowerVariable},
    devices::IS.FlattenIteratorWrapper{V},
    ::DeviceModel{V, W},
    ::Type{<:PM.AbstractPowerModel},
) where {V <: PSY.ElectricLoad, W <: AbstractControllablePowerLoadFormulation}
    time_steps = get_time_steps(container)
    constraint = add_constraints_container!(
        container,
        T(),
        V,
        [PSY.get_name(d) for d in devices],
        time_steps,
    )
    jump_model = get_jump_model(container)
    for t in time_steps, d in devices
        name = PSY.get_name(d)
        pf = sin(atan((PSY.get_max_reactive_power(d) / PSY.get_max_active_power(d))))
        reactive = get_variable(container, U(), V)[name, t]
        real = get_variable(container, ActivePowerVariable(), V)[name, t] * pf
        constraint[name, t] = JuMP.@constraint(jump_model, reactive == real)
    end
end

function add_constraints!(
    container::OptimizationContainer,
    T::Type{ActivePowerVariableLimitsConstraint},
    U::Type{<:VariableType},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    X::Type{<:PM.AbstractPowerModel},
) where {V <: PSY.ControllableLoad, W <: DispatchablePowerLoad}
    add_parameterized_upper_bound_range_constraints(
        container,
        ActivePowerVariableTimeSeriesLimitsConstraint,
        U,
        ActivePowerTimeSeriesParameter,
        devices,
        model,
        X,
    )
    return
end

function add_constraints!(
    container::OptimizationContainer,
    T::Type{ActivePowerVariableLimitsConstraint},
    U::Type{<:VariableType},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    X::Type{<:PM.AbstractPowerModel},
) where {V <: PSY.ControllableLoad, W <: InterruptiblePowerLoad}
    add_parameterized_upper_bound_range_constraints(
        container,
        ActivePowerVariableTimeSeriesLimitsConstraint,
        U,
        ActivePowerTimeSeriesParameter,
        devices,
        model,
        X,
    )
    return
end

############################## FormulationControllable Load Cost ###########################
function objective_function!(
    container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    ::DeviceModel{T, U},
    ::Type{<:PM.AbstractPowerModel},
) where {T <: PSY.ControllableLoad, U <: DispatchablePowerLoad}
    add_variable_cost!(container, ActivePowerVariable(), devices, U())
    return
end

function objective_function!(
    container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    ::DeviceModel{T, U},
    ::Type{<:PM.AbstractPowerModel},
) where {T <: PSY.ControllableLoad, U <: InterruptiblePowerLoad}
    add_proportional_cost!(container, OnVariable(), devices, U())
    return
end