text:
function get_ybus()
    #from matpower
    ybus = [
        (1379.7974868010331, -703.36748339344058), (-1379.7974868010326, 703.36748339344058), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.),
        (-1379.7974868010326, 703.36748339344058), (2149.4368630632848, -1322.9549474641478), (-258.13726441720388, 131.47721535006858), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (-511.5021118450482, 488.11024872063865), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.),
        (0., 0.), (-258.13726441720388, 131.47721535006858), (848.01838235020819, -474.03325448190787), (-347.72101776801628, 177.09070415433175), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (-242.16010016498808, 165.46533497750752), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.),
        (0., 0.), (0., 0.), (-347.72101776801628, 177.09070415433175), (681.65768847345453, -347.16970678365755), (-333.93667070543819, 170.0790026293258), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.),
        (0., 0.), (0., 0.), (0., 0.), (-333.93667070543819, 170.0790026293258), (446.07112747825659, -266.87883289493993), (-112.13445677281842, 96.799830265614105), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.),
        (0., 0.), (0., 0.), (0., 0.), (0., 0.), (-112.13445677281842, 96.799830265614105), (810.81081225833236, -653.40580767389326), (-71.786265476224997, 237.29348896482884), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (-626.89009000928888, 319.31248844345026), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.),
        (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (-71.786265476224997, 237.29348896482884), (274.89952940823412, -304.41736901807127), (-203.11326393200912, 67.123880053242416), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.),
        (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (-203.11326393200912, 67.123880053242416), (305.74544791639141, -140.85962398591522), (-102.63218398438232, 73.735743932672818), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.),
        (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (-102.63218398438232, 73.735743932672818), (204.81480857639477, -146.16404115885513), (-102.18262459201244, 72.428297226182309), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.),
        (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (-102.18262459201244, 72.428297226182309), (837.08726489889182, -315.40287284421652), (-734.9046403068794, 242.97457561803418), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.),
        (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (-734.9046403068794, 242.97457561803418), (1120.7984670263254, -370.57515184861893), (-385.89382671944588, 127.60057623058476), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.),
        (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (-385.89382671944588, 127.60057623058476), (453.3289876750257, -180.65753184311802), (-67.435160955579832, 53.056955612533272), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.),
        (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (-67.435160955579832, 53.056955612533272), (175.73097555203762, -195.60512077779214), (-108.29581459645779, 142.54816516525887), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.),
        (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (-108.29581459645779, 142.54816516525887), (259.62071319261491, -277.22988532565626), (-151.32489859615711, 134.68172016039742), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.),
        (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (-151.32489859615711, 134.68172016039742), (291.38961111396929, -236.96668363147336), (-140.06471251781215, 102.28496347107594), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.),
        (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (-140.06471251781215, 102.28496347107594), (184.74978127162765, -161.94594345980323), (-44.6850687538155, 59.660979988727298), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.),
        (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (-44.6850687538155, 59.660979988727298), (180.27011341680137, -165.98040041865715), (-135.58504466298587, 106.31942042992985), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.),
        (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (-135.58504466298587, 106.31942042992985), (135.58504466298587, -105.91942042992984), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.),
        (0., 0.), (-511.5021118450482, 488.11024872063865), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (570.30762971538297, -541.09854742765776), (-58.80551787033481, 52.988298707019162), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.),
        (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (-58.80551787033481, 52.988298707019162), (224.31234207269267, -246.34230257074509), (-165.50682420235785, 193.35400386372592), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.),
        (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (-165.50682420235785, 193.35400386372592), (247.7758849428638, -302.12927632370901), (-82.269060740505935, 108.77527245998309), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.),
        (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (-82.269060740505935, 108.77527245998309), (82.269060740505935, -108.77527245998309), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.),
        (0., 0.), (0., 0.), (-242.16010016498808, 165.46533497750752), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (352.09329719985436, -252.27338621838055), (-109.9331970348663, 86.808051240873013), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.),
        (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (-109.9331970348663, 86.808051240873013), (220.88165581050799, -173.62274367261969), (-110.9484587756417, 86.814692431746678), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.),
        (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (-110.9484587756417, 86.814692431746678), (110.9484587756417, -86.814692431746678), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.),
        (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (-626.89009000928888, 319.31248844345026), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (1074.7451893356217, -547.33723477801823), (-447.85509932633283, 228.02474633456791), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.),
        (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (-447.85509932633283, 228.02474633456791), (533.00728178818952, -303.10179328873727), (-85.152182461856654, 75.077046954169361), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.),
        (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (-85.152182461856654, 75.077046954169361), (198.45750094642182, -173.78595780940725), (-113.30531848456518, 98.708910855237903), (0., 0.), (0., 0.), (0., 0.), (0., 0.),
        (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (-113.30531848456518, 98.708910855237903), (364.06135511406887, -226.43390675676198), (-250.75603662950371, 127.72499590152407), (0., 0.), (0., 0.), (0., 0.),
        (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (-250.75603662950371, 127.72499590152407), (333.96709674337552, -209.96252764209243), (-83.21106011387181, 82.237531740568343), (0., 0.), (0., 0.),
        (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (-83.21106011387181, 82.237531740568343), (302.0744900671599, -337.3314942678345), (-218.86342995328809, 255.09396252726617), (0., 0.),
        (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (-218.86342995328809, 255.09396252726617), (356.39472533686063, -468.93294478625774), (-137.53129538357254, 213.83898225899159),
        (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (0., 0.), (-137.53129538357254, 213.83898225899159), (137.53129538357254, -213.2389822589916)
    ]

    Ybus = zeros(Complex, 33, 33)
    k = 1
    for i in 1:33
        for j in 1:33
            real = ybus[k][1]
            img = ybus[k][2]
            Ybus[i, j] = real + img * 1im
            k = k + 1
        end
    end
    return Ybus
end


text:
struct QNVal
  name::SmallString
  val::Int
  modulus::Int

  function QNVal(name, v::Int, m::Int=1)
    am = abs(m)
    if am > 1
      return new(SmallString(name), mod(v, am), m)
    end
    return new(SmallString(name), v, m)
  end
end

QNVal(v::Int, m::Int=1) = QNVal("", v, m)
QNVal() = QNVal("", 0, 0)

name(qv::QNVal) = qv.name
val(qv::QNVal) = qv.val
modulus(qv::QNVal) = qv.modulus
isactive(qv::QNVal) = modulus(qv) != 0
(qv1::QNVal < qv2::QNVal) = (name(qv1) < name(qv2))

function qn_mod(val::Int, modulus::Int)
  amod = abs(modulus)
  amod <= 1 && return val
  return mod(val, amod)
end

function -(qv::QNVal)
  return QNVal(name(qv), qn_mod(-val(qv), modulus(qv)), modulus(qv))
end

zero(qv::QNVal) = QNVal(name(qv), 0, modulus(qv))

(dir::Arrow * qv::QNVal) = QNVal(name(qv), Int(dir) * val(qv), modulus(qv))

(qv::QNVal * dir::Arrow) = (dir * qv)

function pm(qv1::QNVal, qv2::QNVal, fac::Int)
  if name(qv1) != name(qv2)
    error("Cannot add QNVals with different names \"$(name(qv1))\", \"$(name(qv2))\"")
  end
  if modulus(qv1) != modulus(qv2)
    error(
      "QNVals with matching name \"$(name(qv1))\" cannot have different modulus values "
    )
  end
  m1 = modulus(qv1)
  if m1 == 1 || m1 == -1
    return QNVal(name(qv1), val(qv1) + fac * val(qv2), m1)
  end
  return QNVal(name(qv1), Base.mod(val(qv1) + fac * val(qv2), abs(m1)), m1)
end

(qv1::QNVal + qv2::QNVal) = pm(qv1, qv2, +1)
(qv1::QNVal - qv2::QNVal) = pm(qv1, qv2, -1)

const ZeroVal = QNVal()

const maxQNs = 4
const QNStorage = SVector{maxQNs,QNVal}
const MQNStorage = MVector{maxQNs,QNVal}

"""
A QN object stores a collection of up to four
named values such as ("Sz",1) or ("N",0).
These values can include a third integer "m"
which makes them obey addition modulo m, for
example ("P",1,2) for a value obeying addition mod 2.
(The default is regular integer addition).

Adding or subtracting pairs of QN objects performs
addition and subtraction element-wise on each of
the named values. If a name is missing from the
collection, its value is treated as zero.
"""
struct QN
  data::QNStorage

  function QN()
    s = QNStorage(ntuple(_ -> ZeroVal, Val(maxQNs)))
    return new(s)
  end

  QN(s::QNStorage) = new(s)
end

QN(mqn::MQNStorage) = QN(QNStorage(mqn))
QN(mqn::NTuple{N,QNVal}) where {N} = QN(QNStorage(mqn))

function hash(obj::QN, h::UInt)
  # TODO: use an MVector or SVector
  # for performance here; put non-zero QNVals
  # to front and slice when passing to hash
  nzqv = QNVal[]
  for qv in obj.data
    if val(qv) != 0
      push!(nzqv, qv)
    end
  end
  return hash(nzqv, h)
end

"""
    QN(qvs...)

Construct a QN from a set of up to four
named value tuples.

Examples
```julia
q = QN(("Sz",1))
q = QN(("N",1),("Sz",-1))
q = QN(("P",0,2),("Sz",0)).
```
"""
function QN(qvs...)
  m = MQNStorage(ntuple(_ -> ZeroVal, Val(maxQNs)))
  for (n, qv) in enumerate(qvs)
    m[n] = QNVal(qv...)
  end
  Nvals = length(qvs)
  sort!(@view m[1:Nvals]; by=name, alg=InsertionSort)
  for n in 1:(length(qvs) - 1)
    if name(m[n]) == name(m[n + 1])
      error("Duplicate name \"$(name(m[n]))\" in QN")
    end
  end
  return QN(QNStorage(m))
end

"""
    QN(name,val::Int,modulus::Int=1)

Construct a QN with a single named value
by providing the name, value, and optional
modulus.
"""
QN(name, val::Int, modulus::Int=1) = QN((name, val, modulus))

"""
    QN(val::Int,modulus::Int=1)

Construct a QN with a single unnamed value
(equivalent to the name being the empty string)
with optional modulus.
"""
QN(val::Int, modulus::Int=1) = QN(("", val, modulus))

data(qn::QN) = qn.data

getindex(q::QN, n::Int) = getindex(data(q), n)

length(qn::QN) = length(data(qn))

lastindex(qn::QN) = length(qn)

isactive(qn::QN) = isactive(qn[1])

function nactive(q::QN)
  for n in 1:maxQNs
    !isactive(q[n]) && (return n - 1)
  end
  return maxQNs
end

function iterate(qn::QN, state::Int=1)
  (state > length(qn)) && return nothing
  return (qn[state], state + 1)
end

keys(qn::QN) = keys(data(qn))

"""
    val(q::QN,name)

Get the value within the QN q
corresponding to the string `name`
"""
function val(q::QN, name_)
  sname = SmallString(name_)
  for n in 1:maxQNs
    name(q[n]) == sname && return val(q[n])
  end
  return 0
end

"""
    modulus(q::QN,name)

Get the modulus within the QN q
corresponding to the string `name`
"""
function modulus(q::QN, name_)
  sname = SmallString(name_)
  for n in 1:maxQNs
    name(q[n]) == sname && return modulus(q[n])
  end
  return 0
end

"""
    zero(q::QN)

Returns a QN object containing
the same names as q, but with
all values set to zero.
"""
function zero(qn::QN)
  mqn = MQNStorage(undef)
  for i in 1:length(mqn)
    mqn[i] = zero(qn[i])
  end
  return QN(mqn)
end

function (dir::Arrow * qn::QN)
  mqn = MQNStorage(undef)
  for i in 1:length(mqn)
    mqn[i] = dir * qn[i]
  end
  return QN(mqn)
end

(qn::QN * dir::Arrow) = (dir * qn)

function -(qn::QN)
  mqn = MQNStorage(undef)
  for i in 1:length(mqn)
    mqn[i] = -qn[i]
  end
  return QN(mqn)
end

function (a::QN + b::QN)
  !isactive(b[1]) && return a

  ma = MQNStorage(data(a))
  for nb in 1:maxQNs
    !isactive(b[nb]) && break
    bname = name(b[nb])
    for na in 1:maxQNs
      aname = name(a[na])
      if !isactive(ma[na])
        ma[na] = b[nb]
        break
      elseif name(ma[na]) == bname
        ma[na] = ma[na] + b[nb]
        break
      elseif (bname < aname) && (na == 1 || bname > name(ma[na - 1]))
        for j in maxQNs:-1:(na + 1)
          ma[j] = ma[j - 1]
        end
        ma[na] = b[nb]
        break
      end
    end
  end
  return QN(QNStorage(ma))
end

(a::QN - b::QN) = (a + (-b))

function hasname(qn::QN, qv_find::QNVal)
  for qv in qn
    name(qv) == name(qv_find) && return true
  end
  return false
end

# Does not perform checks on if QN is already full, drops
# the last QNVal
# Rename insert?
function NDTensors.insertafter(qn::QN, qv::QNVal, pos::Int)
  return QN(NDTensors.insertafter(Tuple(qn), qv, pos)[1:length(qn)])
end

function addqnval(qn::QN, qv_add::QNVal)
  isactive(qn[end]) &&
    error("Cannot add QNVal, QN already contains maximum number of QNVals")
  for (pos, qv) in enumerate(qn)
    if qv_add < qv || !isactive(qv)
      return NDTensors.insertafter(qn, qv_add, pos - 1)
    end
  end
end

# Fills in the qns of qn1 that qn2 has but
# qn1 doesn't
function fillqns_from(qn1::QN, qn2::QN)
  # If qn1 has no non-trivial qns, fill
  # with qn2
  !isactive(qn1) && return zero(qn2)
  for qv2 in qn2
    if !hasname(qn1, qv2)
      qn1 = addqnval(qn1, zero(qv2))
    end
  end
  return qn1
end

# Make sure qn1 and qn2 have all of the same qns
function fillqns(qn1::QN, qn2::QN)
  qn1_filled = fillqns_from(qn1, qn2)
  qn2_filled = fillqns_from(qn2, qn1)
  return qn1_filled, qn2_filled
end

function isequal_assume_filled(qn1::QN, qn2::QN)
  for (qv1, qv2) in zip(qn1, qn2)
    modulus(qv1) != modulus(qv2) && error("QNVals must have same modulus to compare")
    qv1 != qv2 && return false
  end
  return true
end

function ==(qn1::QN, qn2::QN; assume_filled=false)
  if !assume_filled
    qn1, qn2 = fillqns(qn1, qn2)
  end
  return isequal_assume_filled(qn1, qn2)
end

function isless_assume_filled(qn1::QN, qn2::QN)
  for n in 1:length(qn1)
    val1 = val(qn1[n])
    val2 = val(qn2[n])
    val1 != val2 && return val1 < val2
  end
  return false
end

function isless(qn1::QN, qn2::QN; assume_filled=false)
  return <(qn1, qn2; assume_filled=assume_filled)
end

function <(qn1::QN, qn2::QN; assume_filled=false)
  if !assume_filled
    qn1, qn2 = fillqns(qn1, qn2)
  end
  return isless_assume_filled(qn1, qn2)
end

function have_same_qns(qn1::QN, qn2::QN)
  for n in 1:length(qn1)
    name(qn1[n]) != name(qn2[n]) && return false
  end
  return true
end

function have_same_mods(qn1::QN, qn2::QN)
  for n in 1:length(qn1)
    modulus(qn1[n]) != modulus(qn2[n]) && return false
  end
  return true
end

function show(io::IO, q::QN)
  print(io, "QN(")
  Na = nactive(q)
  for n in 1:Na
    v = q[n]
    n > 1 && print(io, ",")
    Na > 1 && print(io, "(")
    if name(v) != SmallString("")
      print(io, "\"$(name(v))\",")
    end
    print(io, "$(val(v))")
    if modulus(v) != 1
      print(io, ",$(modulus(v))")
    end
    Na > 1 && print(io, ")")
  end
  return print(io, ")")
end

function write(parent::Union{HDF5.File,HDF5.Group}, gname::AbstractString, q::QN)
  g = create_group(parent, gname)
  attributes(g)["type"] = "QN"
  attributes(g)["version"] = 1
  names = [String(name(q[n])) for n in 1:maxQNs]
  vals = [val(q[n]) for n in 1:maxQNs]
  mods = [modulus(q[n]) for n in 1:maxQNs]
  write(g, "names", names)
  write(g, "vals", vals)
  return write(g, "mods", mods)
end

function read(parent::Union{HDF5.File,HDF5.Group}, name::AbstractString, ::Type{QN})
  g = open_group(parent, name)
  if read(attributes(g)["type"]) != "QN"
    error("HDF5 group or file does not contain QN data")
  end
  names = read(g, "names")
  vals = read(g, "vals")
  mods = read(g, "mods")
  mqn = ntuple(n -> QNVal(names[n], vals[n], mods[n]), maxQNs)
  return QN(mqn)
end

ext: jl


text:
# ---
# title: 326. Power of Three
# id: problem326
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Math
# link: <https://leetcode.com/problems/power-of-three/description/>
# hidden: true
# ---
#
# Given an integer `n`, return _`true` if it is a power of three. Otherwise,
# return `false`_.
#
# An integer `n` is a power of three, if there exists an integer `x` such that
# `n == 3x`.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 27
#     Output: true
#
#
# **Example 2:**
#
#
#
#     Input: n = 0
#     Output: false
#
#
# **Example 3:**
#
#
#
#     Input: n = 9
#     Output: true
#
#
# **Example 4:**
#
#
#
#     Input: n = 45
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `-231 <= n <= 231 - 1`
#
#
#
# **Follow up:  **Could you do it without using any loop / recursion?
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
function run_maxent_ep_test(net, objider; checkall = false)

    # setup
    net = MetEP.fva_preprocess(net)
    net = MetNets.force_dims(net)
    M, N = size(net)
    objidx = MetNets.rxnindex(net, objider)

    # fbaout
    fbaout = MetLP.fba(net, objidx)
    fba_objval = MetNets.av(net, fbaout, objidx)

    # maxent ep
    epout = MetEP.maxent_ep(net; verbose = false)
    beta_vec = zeros(N)
    Δbeta = 1.1
    beta_vec[objidx] = 10
    ep_objval = MetNets.av(net, epout, objidx)

    obj_th = 1e-4
    prog = ProgressThresh(obj_th)
    for it in 1:5000
        epout_ = try
            beta_vec[objidx] *= Δbeta
            MetEP.maxent_ep(net;
                beta_vec, verbose = false,
                epsconv = 1e-4,
                solution = epout
            )
        catch err;
            beta = beta_vec[objidx]
            @warn("Error", it, beta, err)
            break
        end
        # !MetEP.is_converged(epout_) && break
        epout = epout_

        ep_objval = MetNets.av(net, epout, objidx)
        update!(prog, abs(ep_objval - fba_objval) / fba_objval)
        isapprox(ep_objval, fba_objval; rtol = obj_th) && break
    end
    finish!(prog)

    @show ep_objval
    @show fba_objval
    @test isapprox(ep_objval, fba_objval; rtol = obj_th)

    if checkall
        ep_av = MetNets.av(epout)
        fba_av = MetNets.av(fbaout)
        @show maximum(abs, ep_av - fba_av)

        rtol = 1e-2
        @test all(isapprox.(ep_av, fba_av; rtol = rtol))
    end

end

let
    @info("maxent_ep Toy Model")
    net = MetNets.toy_model()
    objider = MetNets.TOY_MODEL_BIOMASS_IDER
    checkall = true
    run_maxent_ep_test(net, objider; checkall)
    println()

    @info("maxent_ep EColi Core Model")
    net = MetNets.ecoli_core_model()
    objider = MetNets.ECOLI_MODEL_BIOMASS_IDER
    checkall = false # fba (max biomass) is degenerated
    run_maxent_ep_test(net, objider; checkall)
end


text:
#!/usr/bin/env julia
# Copyright 2017 Christian Henning
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#    http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using Distributions

include("./lib/Lists/src/Lists.jl")
using Lists

"""
    sigmoid(z::Float64)

Compute the sigmoid.
"""
function sigmoid(z::Float64)
    return 1.0 / (1.0 + exp(-z))
end

"""
The dataset considered for training and testing. It is always expected, that
the last dimension defines the batch size (cmp.
[MLDatasets.jl](https://github.com/JuliaML/MLDatasets.jl)). If no targets are
defined, then a self-supervised training is expected.
Inputs should be values between 0 and .
"""
type Dataset
    "Training input."
    train_in::Array{Float64, 2}
    "Training output/target."
    train_out::Any
    "Test input."
    test_in::Array{Float64, 2}
    "Test output/target."
    test_out::Any

    """
        Dataset(train_in::Array{Float64, 2}, test_in::Array{Float64, 2};
                train_out::Any = nothing, test_out::Any = nothing)

    Set fields of this type.
    """
    function Dataset(train_in::Array{Float64, 2}, test_in::Array{Float64, 2};
                     train_out::Any = nothing, test_out::Any = nothing)
        x = new()
        x.train_in = train_in
        x.train_out = train_out
        x.test_in = test_in
        x.test_out = test_out
        x
    end
end

"Type containing the hyperparameters of the learning process."
type Hyperparameters
    "Learning Rate."
    λ::Float64
    "Gamma parameter for exponential averaging."
    γ::Float64
    "Loss trade-off for sparsity regularizer."
    αₛ::Float64
    "Loss trade-off for participation regularizer."
    αₚ::Float64
    "Reward baseline."
    rᵦ::Float64
    "Momentum Learning: Influence of recent weight update."
    η::Float64
    """
    Number of epochs for training the network. Each epoch, a mini-batch will be
    presented to the network. Afterwards, the network weights will be updated.
    """
    num_epochs::Int
    "Batch size. Number of samples to process per epoch."
    batch_size::Int
    "Sample size. How often a single batch shall be presented during an epoch"
    sample_size::Int
    "Whether or not gradient clipping should be applied."
    clip_gradients::Bool
    """
    The value to which large gradient values should be clipped. This value
    applies to positive and negative gradients.
    """
    clipping_value::Float64
    "One can enable exponential rate decay."
    rate_decay::Bool
    "After how many epochs the learning rate should be halfed."
    rate_decay_interval::Int

    """
        Hyperparameters(;λ::Float64 = 1.0, γ::Float64 =0.99,
                        αₛ::Float64 = 0.0, αₚ::Float64 = 0.0
                        rᵦ::Float64 = 0.0, η::Float64 = 0.0,
                        num_epochs::Int = 1000, batch_size::Int = 1,
                        sample_size::Int = 1, clip_gradients::Bool = true,
                        clipping_value::Float64 = 1.0, rate_decay::Bool = true,
                        rate_decay_interval::Int = 1000)

    Set default values of type fields.
    """
    function Hyperparameters(;λ::Float64 = 1.0, γ::Float64 =0.99,
                             αₛ::Float64 = 0.0, αₚ::Float64 = 0.0,
                             rᵦ::Float64 = 0.0, η::Float64 = 0.0,
                             num_epochs::Int = 1000, batch_size::Int = 1,
                             sample_size::Int = 1, clip_gradients::Bool = true,
                             clipping_value::Float64 = 1.0,
                             rate_decay::Bool = true,
                             rate_decay_interval::Int = 1000)

        new(λ, γ, αₛ, αₚ, rᵦ, η, num_epochs, batch_size, sample_size,
            clip_gradients, clipping_value, rate_decay, rate_decay_interval)
    end
end

"""
Parameters of loss terms (including regularizers) that may differ within
layers.
"""
type LossParameters
    "Sparsity constraint."
    sparsity::Float64

    """
        LossParameters(;sparsity::Float64 = 0.5)

    Set default values of type fields.
    """
    function LossParameters(;sparsity::Float64 = 0.5)
        new(sparsity)
    end
end

"Struct that represents a single layer."
type Layer
    "The name of the layer."
    name::String
    "Size of layer (number of neurons)."
    size::Int
    "Function handle to forward function (computes forward pass)."
    forward_fcn::Function
    "Function handle to update function (computes weight update)."
    update_fcn::Function
    "Weights of layer."
    W::Array{Float64,2}
    "Recent update history. Used for momentum learning."
    ΔW::Array{Float64,2}
    "Struct of loss parameters defined for this layer."
    loss_params::LossParameters
    """
    Some loss terms might not only depend on the current forward computation
    but also on a state, that needs to be maintained. This is typically the
    case for a participation regularizer, that enforces a cell to be active
    every once in a while.
    """
    loss_memory::Dict{String,Any}
    "The regularizers added to this layer."
    regularizers::List{Function}

    """
        Layer(name::String, size::Int;
              loss_params::LossParameters = LossParameters())

    Set default values of type fields.
    """
    function Layer(name::String, size::Int;
                   loss_params::LossParameters = LossParameters())
        x = new()
        x.name = name
        x.size = size
        x.loss_params = loss_params
        x.loss_memory = Dict{String,Any}()
        x.regularizers = List{Function}()
        x
    end
end

"""
A `Network` struct shall contain the full architecture and state of a network.
"""
type Network
    "List of layers."
    layer::List{Layer}
    "Hyperparameters of the network."
    hparams::Hyperparameters
    """
    The global reward function used to compare output with target.
    Important, the actual reward is computed from this together with the
    specified regularizers.
    """
    global_reward_fcn::Function
    """
        Network(;hparams::Hyperparameters = Hyperparameters,
                global_reward_fcn::Function = mse_reward)

    Set default values of type fields.
    """
    function Network(;hparams::Hyperparameters = Hyperparameters(),
                     global_reward_fcn::Function = mse_reward)
        x = new()
        x.hparams = hparams
        x.layer = List{Layer}()
        x.global_reward_fcn = global_reward_fcn
        x
    end
end

"""
    mse_reward{T,S}(output::AbstractArray{T, 2},
                        target::AbstractArray{S, 2})

Compute the negative mean-squared error between `output` and `target`.

```math
- \\frac{1}{|output|} \\sum_i (output_i - target_i)^2
```
"""
function mse_reward{T,S}(output::AbstractArray{T, 2},
                         target::AbstractArray{S, 2})
    squeeze(-sum((output .- target).^2, 1) / size(output, 1), 1)
end

"""
    regularizer(x::Array{Float64, 2}, y::Array{Float64, 2},
                p::Array{Float64, 2}, layer::Layer, hparams::Hyperparameters)

This is the interface of a typical regularizer that may be implemented within
this framework. We assume, that all interesting regularizers can be expressed
on a layer-wise level.

# Arguments
- `x::Array{Float64, 2}`: the presynaptic output.
- `y::Array{Float64, 2}`: the postsynaptic output (current layer).
- `p::Array{Float64, 2}`: the postsynaptic p-values.
- `layer:Layer`: the postsynaptic layer.
- `hparams::Hyperparameters`: the training hyperparameters.
"""
function regularizer(x::AbstractArray{Bool, 2}, y::AbstractArray{Bool, 2},
                     p::Array{Float64, 2}, layer::Layer,
                     hparams::Hyperparameters)
    error("Not callable! Only abstract function interface.");
end

"""
    sparsity_regularizer(::Any, y::AbstractArray{Bool, 2}, ::Any, layer::Layer,
                         hparams::Hyperparameters)

A sparsity regularizer, that enforces the sparsity level defined in the
`LossParLossParameters` in this layer.

Across a batch the mean loss is returned.

```math
    - αₛ ((1-sparsity) - \\frac{1}{|y|} \\sum_i y_i)^2
```
"""
function sparsity_regularizer(::Any, y::AbstractArray{Bool, 2},
                              ::Any, layer::Layer, hparams::Hyperparameters)
    αₛ = hparams.αₛ
    part = 1.0 - layer.loss_params.sparsity
    # Note, that julia autmatically converts divisions into floats.
    squeeze(-αₛ * (part .- sum(y, 1) ./ size(y, 1)).^2, 1)
end

"""
    participation_regularizer(::Any, y::AbstractArray{Bool, 2}, ::Any,
                              layer::Layer, hparams::Hyperparameters)

A participation regularizer, that enforces neurons to participate on the
computation every once in a while (level of participation is defined by
``1-sparsity``). Therefore, this regularizer has to maintain a state of past
actactivation, which is kept as an exponentially smoothed trace.

```math
 - αₚ \\sum_i ((1-sparsity) - y_{trace,i})^2
```
"""
function participation_regularizer(::Any, y::AbstractArray{Bool, 2}, ::Any,
                                   layer::Layer, hparams::Hyperparameters)
    γ = hparams.γ
    αₚ = hparams.αₚ
    part = 1.0 - layer.loss_params.sparsity

    # Initialize trace, if not existing.
    if get(layer.loss_memory, "y_trace", -1) == -1
        layer.loss_memory["y_trace"] = part .* ones(size(y, 1))
    end
    y_trace = layer.loss_memory["y_trace"]

    y_trace_batch = γ .* y_trace .+ (1-γ) .* y
    # Update trace.
    y_trace .= squeeze(mean(y_trace_batch, 2), 2)

    squeeze(-αₚ * sum((part .- y_trace_batch).^2, 1) ./ length(y_trace), 1)
end

"""
    fully_connected_forward{T}(x::AbstractArray{T, 2},
                               W::Array{Float64, 2})

Computes the forward propagation of an input `x`  through a fully-connected
layer in a binomial network.

# Returns
Returns the p-values ``p = σ(W x)`` of the layer and the layer activation
``y = Bin(p)`` as an tuple:

    p::Array{Float64, 2}, y::BitArray{2}
"""
function fully_connected_forward{T}(x::AbstractArray{T, 2},
                                    W::Array{Float64, 2})
    p = sigmoid.(W * x)
    y = rand(size(p)) .< p

    p, y
end

"""
    fully_connected_update(x::Array{Float64, 2}, y::Array{Float64, 2},
                           p::Array{Float64, 2}, r::Array{Float64, 1},
                           hparams::Hyperparameters)

Compute the weight update `ΔW` of the current, fully-connected layer.

```math
ΔW = (r - rᵦ) * (y - p) * x^T
```

This method returns a 3-rank tensor ΔW, where the last dimension specifies the
sample from the batch that defines the weight update.
"""
function fully_connected_update(x::Array{Float64, 2}, y::Array{Float64, 2},
                                p::Array{Float64, 2}, r::Array{Float64, 1},
                                hparams::Hyperparameters)
    rᵦ = hparams.rᵦ
    ΔW = zeros(size(y, 1), size(x, 1), size(y, 2))
    for b = 1:size(y, 2)
        ΔW[:,:,b] = (r[b] - rᵦ) .* (y[:,b] .- p[:,b]) * x[:,b]'
    end

    ΔW
end

"""
    fully_connected_layer(name::String, in_size::Int, out_size::Int;
                          loss_params::LossParameters = LossParameters(),
                          initializer::Function = rand)

Create a fully-connected layer. The method returns a `layer` type.

# Arguments
- `name::String`: the name of the layer.
- `in_layer`: input layer to this one.
- `out_size`: output size.
- `loss_params`: loss parameters that shall apply to this layer.
- `initializer`: a function handle that is used to initialize the layer
    weights.

# Examples

```jldoctest
julia> init(dims...) = 0.1 * rand(dims)
init (generic function with 1 method)

julia> layer = fully_connected_layer("fc1", 784, 100, initializer=init)
Layer("fc1", ...
```
"""
function fully_connected_layer(name::String, in_layer::Layer, out_size::Int;
                               loss_params::LossParameters = LossParameters(),
                               initializer::Function = rand)
    layer = Layer(name, out_size, loss_params = loss_params)
    layer.W = initializer(out_size, in_layer.size)
    layer.ΔW = zeros(out_size, in_layer.size)
    forward_fcn(x) = fully_connected_forward(x, layer.W)
    layer.forward_fcn = forward_fcn
    layer.update_fcn = fully_connected_update

    layer
end

"""
    add_layer(network::Network, layer::Layer)

Add `layer` to the current end of the `network`.
"""
function add_layer(network::Network, layer::Layer)
    push!(network.layer, layer)
end

"""
    add_regularizer(layer::Layer, reg_fun::Function)

Add a regularizer function to a given layer. The `reg_fun` must have the same
interface as the function `regularizer(...)`.
"""
function add_regularizer(layer::Layer, reg_fun::Function)
    push!(layer.regularizers, reg_fun)
end

"""
    input_layer(name::String, size::Int)

Create an input layer. Every network has to start with an input layer. The
forward function of this layer is simply the identity function. An update
function does not exist for such a layer.
"""
function input_layer(name::String, size::Int)
    layer = Layer(name, size)
    layer.forward_fcn = identity
    update_fcn(args...) = error("Input layer has no weights to update")
    layer.update_fcn = update_fcn

    layer
end

"""
    train_network(network::Network, train_in::Array{Float64, 2};
                  train_out::Any = nothing, eval_interval::Int = 100,
                  test_in::Any = nothing, test_out::Any = nothing)

Train a `network` with the given training set. If no training output
`train_out` is provided, then self-supervision is assumed.

If you want the test accuracy to be printed, then please provide a test set.

# Returns

It returns a tuple containing two arrays. The first one maintains the training
loss development (one value each epoch). The second array contains the
evaluation of the test set (e.g., the test accuracy).
"""
function train_network(network::Network, train_in::Array{Float64, 2};
                       train_out::Any = nothing, eval_interval::Int = 100,
                       eval_accuracy::Bool = true, test_in::Any = nothing,
                       test_out::Any = nothing)
    num_samples = size(train_in, ndims(train_in))
    batch_indices = randperm(num_samples)
    batch_size = network.hparams.batch_size

    # Keep History of train loss and test evaluation values for later
    # assessment of the training progress.
    train_loss_dev = zeros(network.hparams.num_epochs)
    test_eval_dev = []

    hparams = network.hparams

    for epoch = 1:network.hparams.num_epochs
        # Apply exponential rate decay.
        if hparams.rate_decay && epoch % hparams.rate_decay_interval == 0
            hparams.λ /= 2
            @printf("Epoch %d, learning rate decayed to: %f\n",
                    epoch, hparams.λ)
        end

        # Determine current batch.
        bs_start = ((epoch-1) * batch_size) % num_samples + 1
        bs_end = bs_start + batch_size -1

        if bs_end > num_samples
            tmp_batch_inds = batch_indices[bs_start:num_samples]
            batch_indices = randperm(num_samples)
            batch_inds = vcat(tmp_batch_inds,
                              batch_indices[1:bs_end-num_samples])
        else
            batch_inds = batch_indices[bs_start:bs_end]
        end
        batch_in = train_in[:, batch_inds]
        if train_out == nothing
            batch_out = batch_in
        else
            batch_out = train_out[:, batch_inds]
        end

        # Process batch and update weights.
        train_loss = train_step(network, batch_in, batch_out)
        train_loss_dev[epoch] = train_loss

        # Report progress.
        if epoch % eval_interval == 0
            if test_in == nothing
                @printf("Epoch %d, loss of current train batch: %f\n",
                        epoch, train_loss)
            else
                if test_out == nothing
                    test_out = test_in
                end

                if eval_accuracy
                    test_acc = accuracy(network, test_in, test_out)
                    @printf("Epoch %d, current test accuracy: %f\n", epoch,
                            test_acc)
                    push!(test_eval_dev, test_acc)
                else
                    test_pred = predict(network, test_in)
                    test_loss = -mean(network.global_reward_fcn(test_pred,
                                                                test_out))
                    @printf("Epoch %d, global test loss: %f\n", epoch,
                            test_loss)
                    push!(test_eval_dev, test_loss)
                end
            end
        end
    end

    train_loss_dev, test_eval_dev
end

"""
    one_hot(labels::Array{Int, 1}; num_labels::Int = -1)

Encode a given vector of `labels` into an one hot array.

# Return
An array of type `Array{Bool, 2}`, where the last dimension has the same size
as `labels`.
"""
function one_hot(labels::Array{Int, 1}; num_labels::Int = -1)
    num_labels = (num_labels == -1 ? maximum(labels) : num_labels)
    encoded = zeros(Bool, num_labels, length(labels))
    for i = 1:length(labels)
        encoded[labels[i], i] = 1
    end
    encoded
end

"""
    train_step{T}(network::Network, batch_in::Array{Float64, 2},
                  batch_out::Array{T, 2})

Train the network according to the current batch. This includes the sampling of
the batch (i.e., a sufficient amount of forward sweeps (acc. to the
hyperparamers) until we have good estimates of the network activity.). The
weight update is the average update ``ΔW`` that would have been applied to the
each sample in the batch.
"""
function train_step{T}(network::Network, batch_in::Array{Float64, 2},
                       batch_out::Array{T, 2})
    layers = network.layer

    hparams = network.hparams
    sample_size = hparams.sample_size

    # During a forward sweep, we need to memorize the layer activities, as we
    # cannot compute ΔW without knowing the final reward.
    X = Array{Float64, 2}[]
    P = Array{Float64, 2}[]
    Y = Array{Float64, 2}[]
    # Reward is global (one number per sample in batch).
    R = zeros(size(batch_in, 2))

    # For each layer, keep a list of update matrices (one for each sample in
    # the batch).
    ΔW_batch = Array{Float64, 3}[]

    x = nothing
    p = nothing
    y = nothing

    # Sweep through the network `sample_size` times to get an average weight
    # update ΔW for each sample in the batch.
    for ss in 1:sample_size
        r = zeros(size(batch_in, 2))

        # Compute the activities an the final reward.
        for (lind, curr_layer) in enumerate(network.layer)
            # First layer is always input layer.
            if lind == 1
                y = curr_layer.forward_fcn(batch_in)
                continue
            end

            x = copy(y)

            p, y = curr_layer.forward_fcn(x)

            for reg_fcn in curr_layer.regularizers
                r .+= reg_fcn(x, y, p, curr_layer, hparams)
            end

            if ss == 1
                push!(X, x)
                push!(P, p)
                push!(Y, y)
            else
                X[lind-1] = x
                Y[lind-1] = y
                P[lind-1] = p
            end
        end
        r .+= network.global_reward_fcn(y, batch_out)
        R .+= r

        # Compute ΔW for each layer.
        for (lind, curr_layer) in enumerate(network.layer)
            # Input layers don't have weights.
            if lind == 1
                continue
            end

            x = X[lind-1]
            p = P[lind-1]
            y = Y[lind-1]

            ΔW = curr_layer.update_fcn(x, y, p, r, hparams)

            if ss == 1
                push!(ΔW_batch, ΔW)
            else
                ΔW_batch[lind-1] += ΔW
            end
        end
    end

    ΔW_batch /= sample_size

    ### FIXME dirty insertion
    for (lind, curr_layer) in enumerate(network.layer)
        # Input layers don't have weights.
        if lind == 1
            continue
        end

        y_trace = curr_layer.loss_memory["y_trace"]
        norm_unit = (y_trace .* (1 .- y_trace)) .< 1e-2
        norm_updates = norm_unit .* (0.5 .- y_trace)
        ΔW_batch[lind-1] .+= 1 .* norm_updates

        if lind == 4 && sum(norm_updates) > 0
            println("layer $(lind)")
            println(y_trace)
            println(norm_updates)
        end
    end
    ###

    ΔW_batch *= hparams.λ
    R /= sample_size

    # Go once through the network and update all weights, given the samples
    # weight update.
    for (lind, curr_layer) in enumerate(network.layer)
        # Input layers don't have weights.
        if lind == 1
            continue
        end

        # Use the mean update matrix ΔW across a batch.
        ΔW = squeeze(mean(ΔW_batch[lind-1], 3), 3)

        # If we are using momentum learning.
        if hparams.η > 0
            curr_layer.ΔW = hparams.η * curr_layer.ΔW + ΔW
            ΔW = curr_layer.ΔW
        end

        # Clip gradients.
        if hparams.clip_gradients
            ΔW[ΔW .> hparams.clipping_value] = hparams.clipping_value
            ΔW[ΔW .< -hparams.clipping_value] = -hparams.clipping_value

            if hparams.η > 0
                curr_layer.ΔW .= ΔW
            end
        end

        curr_layer.W .+= ΔW
    end

    train_loss = -mean(R)
end

"""
    accuracy{T}(network::Network, batch_in::Array{Float64, 2},
                batch_out::Array{T, 2})

This method calls the `predict` method to process the given input batch in the
network and subsequently compute an accuracy based on the output vector and the
one-hot encoded given output batch.
"""
function accuracy{T}(network::Network, batch_in::Array{Float64, 2},
                     batch_out::Array{T, 2})
    Y = predict(network, batch_in)

    correct_ind = mapslices(indmax, batch_out, 1) # argmax function
    predicted_ind = mapslices(indmax, Y, 1)

    sum(correct_ind .== predicted_ind) / length(correct_ind)
end

"""
    predict(network::Network, batch_in::Array{Float64, 2})

Process a batch with the given network and return the output activity for each
sample in the batch.
"""
function predict(network::Network, batch_in::Array{Float64, 2})
    layers = network.layer

    hparams = network.hparams
    sample_size = hparams.sample_size

    out_size = network.layer[end].size

    # We need the average output activity.
    Y = zeros(out_size, size(batch_in, 2))

    x = nothing
    p = nothing
    y = nothing

    for ss in 1:sample_size
        for (lind, curr_layer) in enumerate(network.layer)
            # First layer is always input layer.
            if lind == 1
                y = curr_layer.forward_fcn(batch_in)
                continue
            end

            x = copy(y)

            p, y = curr_layer.forward_fcn(x)
        end
        Y .+= y
    end

    Y /= sample_size
end

"""
    xavier_initializer(n_in::Int, n_out::Int)

This method returns a function that can be used as a weight initializer for a
layer. The weights sampled by this function follow the Xavier initialization
for a layer with `n_in` input neurons and `n_out` output neurons.

The distribution is a normal distribution with mean ``μ = 0`` and varianve
``σ = \\sqrt{\\frac{2}{n_{in} + n_{out}}}``.
"""
function xavier_initializer(n_in::Int, n_out::Int)
    #dist = Normal(0, 1/n_in) # As implemented in Caffe.
    dist = Normal(0, sqrt(2/(n_in+n_out))) # As implemented in Tensorflow.

    init(dims...) = rand(dist, dims...)
end

"""
    neural_bias_data_transform(data::Array{Float64, 2})

The neuron model implemented here does not consider neural biases. As a trick,
one can add an additional constant dimension to the input data, whose learned
corresponding weight is then considered as the bias.

# Return
A copy of the given data with the size of the first dimension being increased
by one, such that a 1 is added for all samples at the end of the first
dimension.
"""
function neural_bias_data_transform(data::Array{Float64, 2})
    batch_size = size(data)[2]
    tdata = vcat(data, ones(1, batch_size))

    tdata
end


text:
#This module defines the Spelliptic module and the corresponding operator
#overloading.
module Spelliptic
using IterativeSolvers
include("SortSparse.jl")
include("Factorization.jl")

export SpelMat
export SpelMatIter
export setCovStat!
export testAccuracy

#A "spelliptic" operator, that will contain the sparse fatorization of an
#elliptic integral operator/covariance matrix.
struct SpelMat
  #Size of the matrix
  N::Int64
  #Number of nonzero entries of the Factor
  nnz::Int64
  #An integer array containing the IDs of the dofs in order of the elimination
  P::Array{Int64,1}
  #An integer Array that contains as ith entry the number that the dof with Id i
  #has in the elimination ordering
  revP::Array{Int64,1}
  #An array of doubles that contains the "characteristic distance" of each of
  #the dofs, ordered in elimination order
  distances::Array{Float64,1}
  #The upper-triangular matrix containing the sparse factor
  U::SparseMatrixCSC{Float64, Int64}
end

#Constructors
#Constructor to generate a new SpelMat from a d \times N array of N points in
#d-dimensional space and an accuracy factor \rho.
function SpelMat( x::Array{Float64,2}, rho::Float64 )
  N::Int64 = size(x,2)
  P, revP, distances = sortPoints( x );
  ind, jnd = sparsityPattern( x, P, revP, distances, rho )
  nnz::Int64 = size( jnd, 1 )
  return SpelMat( N,
                  nnz,
                  P,
                  revP,
                  distances,
                  SparseMatrixCSC{Float64, Int64}( N,
                                                   N,
                                                   ind,
                                                   jnd,
                                                   ones(nnz))'
                 )
end

#Constructs a SpelMat, recycling the pattern, by reusing a given P, revP,
#and rho
function SpelMat( x::Array{Float64,2},
                  rho::Float64,
                  P::Array{Int64,1},
                  revP::Array{Int64, 1},
                  distances::Array{Float64, 1} )
  N::Int64 = size(x,2)
  ind, jnd = sparsityPattern( x, P, revP, distances, rho )
  nnz::Int64 = size( jnd, 1 )
  return SpelMat( N,
                  nnz,
                  P,
                  revP,
                  distances,
                  SparseMatrixCSC{Float64, Int64}( N,
                                                   N,
                                                   ind,
                                                   jnd,
                                                   ones(nnz))'
                 )
end

#Operator and function overloading:
#compute the log determinant
import Base.logdet
function logdet( K::SpelMat )
  return 2 * sum( log.( diag( K.U )))
end

#Operator overloading:
import Base.*
function *(K::SpelMat, v )
  return (((K.U * v[K.P])'*K.U)[K.revP])'
end

import Base.*
function *( v, K::SpelMat )
  return (((K.U * v[K.P]')'*K.U)[K.revP])
end

import Base.\
function \( K::SpelMat, v )
  return ((K.U \ ( v[ K.P, : ]' / K.U  )' )[ K.revP, : ] )
end

import Base./
function /( v, K::SpelMat )
  return (( K.U \ ( v[ : , K.P ] / K.U )' )[ : , K.revP ] )'
end

import Base.getindex
function getindex( K::SpelMat, i::Int64, j::Int64 )
  return dot( K.U[ :, K.revP[i]], K.U[ :, K.revP[j]])
end

function getindex( K::SpelMat, ivec, jvec )
  return K.U[ :, K.revP[ivec]]' * K.U[ :, K.revP[jvec]]
end

#Function to set the content of a SpelMat struct to the those given by a
#stationary covariance function.
#The variable sigma allows to add an optional nugget term to the covariance.
#Note that for larger nugget terms the approximation property will detoriate,
#Use SpelMatIter when working with large nugget terms.
function setCovStat!( K::SpelMat, x::Array{Float64,2}, covFunc, sigma::Float64 = 0. )
  for i = 1 : K.N
    for j = K.U.colptr[ i ] : ( K.U.colptr[ i + 1 ] - 1 )
      K.U.nzval[ j ] = norm( x[ :, K.P[ i ] ] - x[ :, K.P[ K.U.rowval[ j ] ] ] )
    end
  end
  K.U.nzval .= map( covFunc, K.U.nzval )
  if sigma != 0.
    for i = 1 : K.N
      for j = K.U.colptr[ i ] : ( K.U.colptr[ i + 1 ] - 1 )
        if K.U.rowval[ j ] == i
          K.U.nzval[ j ] += sigma
        end
      end
    end
  end
  @time icholU_high_level!( K.U )
end

function testAccuracy( K::SpelMat,
                      x::Array{Float64,2},
                      covFunc,
                      NTest::Int64,
                      sigma::Float64 = 0. )

  testI = randperm( K.N )[ 1 : NTest ]
  testJ = randperm( K.N )[ 1 : NTest ]
  errorMat = zeros( NTest, NTest )
  distMat = zeros( NTest, NTest )
  trueMat = zeros( NTest, NTest )
  appMat = zeros( NTest, NTest )

  for i = 1 : NTest
    for j = 1 : NTest
      trueMat[ i, j ] = covFunc( norm( x[ :, testI[ i ] ] -
             x[ :, testJ[ j ] ] ) )
      if ( testI[i] == testJ[j] )
        trueMat[ i , j ] += sigma
      end
      appMat[ i, j ] = K[ testI[ i ], testJ[ j ]]
      distMat[ i, j ] = norm( x[ :, testI[ i ]] - x[ :, testJ[ j ]] )
    end
  end
  errorMat = abs.( trueMat - appMat )
  return trueMat[:], appMat[:], errorMat[:], distMat[:]
end

###############################################################################
#An extension of SpelMat that treats the "nugget" by using iterative methods.
#The compression is only applied to the part that comes from a Green's function,
#then, conjugate gradient is used to invert the sum of this compressed matrix
#and its nugget.
mutable struct SpelMatIter
  K::SpelMat
  sigma::Float64
  tol::Float64
end

#Constructors
#Constructor to generate a new SpelMatIter from a d \times N array of N points in
#d-dimensional space and an accuracy factor \rho.
function SpelMatIter( x::Array{Float64,2}, rho::Float64 )
  return SpelMatIter( SpelMat( x, rho ), 0., 1e-3 )
end

#Constructs a SpelMatIter, recycling the pattern, by reusing a given P, revP,
#and rho
function SpelMatIter( x::Array{Float64,2},
                  rho::Float64,
                  P::Array{Int64,1},
                  revP::Array{Int64, 1},
                  distances::Array{Float64, 1} )

  return SpelMatIter( SpelMatIter( x::Array{Float64,2},
                      rho::Float64,
                      P::Array{Int64,1},
                      revP::Array{Int64, 1},
                      distances::Array{Float64, 1} ),
                      0., 1e-3 )
end

#Operator and function overloading:
function setCovStat!( K::SpelMatIter, x::Array{Float64,2}, covFunc, sigma::Float64, sigmaimp = 0.)
  setCovStat!( K.K, x, covFunc, sigmaimp )
  K.sigma = sigma
end

#Operator overloading:
import Base.*
function *(K::SpelMatIter, v )
  return (((K.K.U * v[K.K.P])'*K.K.U)[K.K.revP])' + K.sigma * v
end

import Base.*
function *( v, K::SpelMatIter )
  return (((K.K.U * v[K.K.P]')'*K.K.U)[K.K.revP]) + K.sigma * v
end

import Base.\
function \( K::SpelMatIter, v )
  return cg( K, v; verbose = false, tol = K.tol  )
end

import Base./
function /( v, K::SpelMatIter )
  return cg( K, v'; verbose = false, tol = K.tol )
end

import Base.getindex
function getindex( K::SpelMatIter, i::Int64, j::Int64 )
  return dot( K.K.U[ :, K.K.revP[i]], K.K.U[ :, K.K.revP[j]]) +
  K.sigma * convert( Float64, i==j )
end

function getindex( K::SpelMatIter, ivec, jvec )
  return K.K.U[ :, K.K.revP[ivec]]' * K.K.U[ :, K.K.revP[jvec]] +
  K.sigma * speye( K.K.N )[ivec, jvec]
end

import Base.eltype
function eltype( K::SpelMatIter )
  return eltype( K.K.U )
end

import Base.size
function size( K::SpelMatIter, d )
  return size( K.K.U, d )
end

import Base.A_mul_B!
function A_mul_B!( y, K::SpelMatIter, x )
  y .= K * x
end

#the end of Module SpellipticRL
end

max_stars_repo_name: f-t-s/Spelliptic

max_issues_repo_name: f-t-s/Spelliptic

max_forks_repo_name: f-t-s/Spelliptic


text:
using Base64

function testfn1(arg1, arg2; narg1=1, narg2=2)
    a1 = isa(arg1, Int) ? arg1 : parse(Int, arg1)
    a2 = isa(arg2, Int) ? arg2 : parse(Int, arg2)
    na1 = isa(narg1, Int) ? narg1 : parse(Int, narg1)
    na2 = isa(narg2, Int) ? narg2 : parse(Int, narg2)
    return (a1 * na1) + (a2 * na2)
end
testfn2(arg1, arg2; narg1=1, narg2=2) = testfn1(arg1, arg2; narg1=narg1, narg2=narg2)

testbinary(datalen::String) = testbinary(parse(Int,datalen))
testbinary(datalen::Int) = rand(UInt8, datalen)

testArray(x::Array{Float64, 2}) = sum(x) + x[1,2]

function testFile(;filename=nothing, filedata=nothing)
    filename = base64decode(filename)
    filedata = base64decode(filedata)
    #println("[", String(filename), "]")
    #println("[", String(filedata), "]")
    string(length(filename)) * "," * string(length(filedata))
end

testException() = error("testing exception handling")


text:
function GLONASSEphemeris(arg0::jint, arg1::jint, arg10::jdouble, arg11::jdouble, arg12::AbsoluteDate, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::jdouble, arg6::jdouble, arg7::jdouble, arg8::jdouble, arg9::jdouble)
    return GLONASSEphemeris((jint, jint, jdouble, jdouble, AbsoluteDate, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg10, arg11, arg12, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function GLONASSEphemeris(arg0::jint, arg1::jint, arg10::jdouble, arg11::jdouble, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::jdouble, arg6::jdouble, arg7::jdouble, arg8::jdouble, arg9::jdouble)
    return GLONASSEphemeris((jint, jint, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg10, arg11, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function get_date(obj::GLONASSEphemeris)
    return jcall(obj, "getDate", AbsoluteDate, ())
end

function get_delta_i(obj::GLONASSOrbitalElements)
    return jcall(obj, "getDeltaI", jdouble, ())
end

function get_delta_t(obj::GLONASSOrbitalElements)
    return jcall(obj, "getDeltaT", jdouble, ())
end

function get_delta_t_dot(obj::GLONASSOrbitalElements)
    return jcall(obj, "getDeltaTDot", jdouble, ())
end

function get_e(obj::GLONASSOrbitalElements)
    return jcall(obj, "getE", jdouble, ())
end

function get_lambda(obj::GLONASSOrbitalElements)
    return jcall(obj, "getLambda", jdouble, ())
end

function get_n4(obj::GLONASSEphemeris)
    return jcall(obj, "getN4", jint, ())
end

function get_na(obj::GLONASSEphemeris)
    return jcall(obj, "getNa", jint, ())
end

function get_pa(obj::GLONASSOrbitalElements)
    return jcall(obj, "getPa", jdouble, ())
end

function get_time(obj::GLONASSEphemeris)
    return jcall(obj, "getTime", jdouble, ())
end

function get_x(obj::GLONASSEphemeris)
    return jcall(obj, "getX", jdouble, ())
end

function get_x_dot(obj::GLONASSEphemeris)
    return jcall(obj, "getXDot", jdouble, ())
end

function get_x_dot_dot(obj::GLONASSEphemeris)
    return jcall(obj, "getXDotDot", jdouble, ())
end

function get_y(obj::GLONASSEphemeris)
    return jcall(obj, "getY", jdouble, ())
end

function get_y_dot(obj::GLONASSEphemeris)
    return jcall(obj, "getYDot", jdouble, ())
end

function get_y_dot_dot(obj::GLONASSEphemeris)
    return jcall(obj, "getYDotDot", jdouble, ())
end

function get_z(obj::GLONASSEphemeris)
    return jcall(obj, "getZ", jdouble, ())
end

function get_z_dot(obj::GLONASSEphemeris)
    return jcall(obj, "getZDot", jdouble, ())
end

function get_z_dot_dot(obj::GLONASSEphemeris)
    return jcall(obj, "getZDotDot", jdouble, ())
end


text:
function checkTagName(xmlel::XMLElement,tagname::String)
	return isequal(xmlel.tag.name,tagname)
end

"""
`getElements!(ret::Vector{XMLElement},xmlel::XMLElement,tagname::String)`

Pushes all `XMLElement`s that match the `tagname` to `ret`.
"""
function getElements!(ret::Vector{XMLElement},xmlel::XMLElement,tagname::String)
	for con in xmlel.content
		if typeof(con)==XMLElement
			if checkTagName(con,tagname)
				push!(ret,con)
			else
				getElements!(ret,con,tagname)
			end
		end
	end
	return nothing
end

"""
`getElements(xmlel::XMLElement,tagname::String)`

Gets all `XMLElement`s that match the `tagname` and returns a `Vector{XMLElement}`.
"""
function getElements(xmlel::XMLElement,tagname::String)
	ret = Vector{XMLElement}()
	getElements!(ret,xmlel,tagname)
	return ret
end

"""
`getChildrenbyTagName!(ret::Vector{XMLElement},xmlel::XMLElement,tagname::String)`

Pushes all `XMLElement`s that are direct children of `xmlel` and match the `tagname` to `ret`.
"""
function getChildrenbyTagName!(ret::Vector{XMLElement},xmlel::XMLElement,tagname::String)
	for con in xmlel.content
		if typeof(con)==XMLElement
			if checkTagName(con,tagname)
				push!(ret,con)
			end
		end
	end
	return nothing
end

"""
`getChildrenbyTagName(xmlel::XMLElement,tagname::String)`

Gets all `XMLElement`s that are direct children of `xmlel` and match the `tagname`.
Return a `Vector{XMLElement}`.
"""
function getChildrenbyTagName(xmlel::XMLElement,tagname::String)
	ret = Vector{XMLElement}()
	getChildrenbyTagName!(ret,xmlel,tagname)
	return ret
end

function copyElementWoContent(el::XMLElement)
	new_el = XMLElement(XMLTag(el.tag.name),Any[])
	foreach(x->push!(new_el.tag.attributes,XMLAttribute(x.key,x.val)),el.tag.attributes)
	return new_el
end

#function getAllPathesbyTag!(ret::Vector{String},xmlel::XMLElement,tagname::String,path="./")
#	for con in xmlel.content
#		if typeof(con)==XMLElement
#			_path = joinpath(path,con.tag.name)
#			for attribute in con.tag.attributes
#				_path = joinpath(_path,"@$(attribute.key)=$(attribute.val)")
#			end
#			if checkTagName(con,tagname)
#				push!(ret,_path)
#			else
#				getAllPathesbyTag!(ret,con,tagname,_path)
#			end
#		end
#	end
#	return nothing
#end
#
#function getAllPathesbyTag(xmlel::XMLElement,tagname::String,path="./")
#	ret = Vector{String}()
#	getAllPathesbyTag!(ret,xmlel,tagname,path)
#	return ret
#end

"""
`hasAttributekey(el::XMLElement,key::String))`

Return true if `el` has a attribute with keyname `key`. Returns false otherwise.
"""
function hasAttributekey(el::XMLElement,key::String)
	for attribute in el.tag.attributes
		if isequal(attribute.key,key)
			return true
		end
	end
	return false
end

"""
`getAttribute(el::XMLElement,key::String)`

Return attribute value if `el` has a attribute with keyname `key`. Returns emptystring otherwise.
"""
function getAttribute(el::XMLElement,key::String)
	for attribute in el.tag.attributes
		if isequal(attribute.key,key)
			return attribute.val
		end
	end
	return ""
end

"""
`setAttribute(el::XMLElement,key::String,val)`

Sets attribute value `val` to XML attribute `key` of `XMLElement` `el`.
Returns an `error("Attribute not found")` if `key` is not found. Returns `nothing` otherwise.
"""
function setAttribute(el::XMLElement,key::String,val)
	for attribute in el.tag.attributes
		if isequal(attribute.key,key)
			attribute.val = "\""*string(val)*"\""
			return nothing
		end
	end
	return error("Attribute not found")
end

Base.show(io::IO,tag::Union{XMLTag,XMLEmptyTag}) = print(io,"<"*tag.name*">")
function Base.show(io::IO,tag::Union{XMLTag,XMLEmptyTag},depth)
	#foreach(x->print("  "),1:depth)
	print(io,"<"*tag.name*">")
end
function Base.show(io::IO,el::XMLElement)
	depth = 0
	show(io,el.tag)
	for elc in el.content
		if typeof(elc) == XMLElement
			show(io,elc,depth+1)
		else
			print(io,elc)
		end
	end
	show(io,el.tag,depth)
end
function Base.show(io::IO,el::XMLElement,depth)
	show(io,el.tag,depth)
	for elc in el.content
		if typeof(elc) == XMLElement
			show(io,elc,depth+1)
		else
			#foreach(x->print("  "),1:depth)
			print(io,elc)
		end
	end
	show(io,el.tag,depth)
end

function Base.deepcopy(attr::XMLAttribute)
	return XMLAttribute(deepcopy(attr.key),deepcopy(attr.val))
end

function Base.deepcopy(tag::Tag) where {Tag<:AbstractXMLTag}
	return Tag(deepcopy(tag.name), map(deepcopy,tag.attributes))
end

function Base.deepcopy(el::XMLElement)
	return XMLElement(deepcopy(el.tag), map(deepcopy,el.content))
end

max_stars_repo_name: baxmittens/XMLParser

max_issues_repo_name: baxmittens/XMLParser

max_forks_repo_name: baxmittens/XMLParser


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiPolicyV1beta1RunAsUserStrategyOptions <: SwaggerModel
    ranges::Any # spec type: Union{ Nothing, Vector{IoK8sApiPolicyV1beta1IDRange} } # spec name: ranges
    rule::Any # spec type: Union{ Nothing, String } # spec name: rule

    function IoK8sApiPolicyV1beta1RunAsUserStrategyOptions(;ranges=nothing, rule=nothing)
        o = new()
        validate_property(IoK8sApiPolicyV1beta1RunAsUserStrategyOptions, Symbol("ranges"), ranges)
        setfield!(o, Symbol("ranges"), ranges)
        validate_property(IoK8sApiPolicyV1beta1RunAsUserStrategyOptions, Symbol("rule"), rule)
        setfield!(o, Symbol("rule"), rule)
        o
    end
end # type IoK8sApiPolicyV1beta1RunAsUserStrategyOptions

const _property_map_IoK8sApiPolicyV1beta1RunAsUserStrategyOptions = Dict{Symbol,Symbol}(Symbol("ranges")=>Symbol("ranges"), Symbol("rule")=>Symbol("rule"))
const _property_types_IoK8sApiPolicyV1beta1RunAsUserStrategyOptions = Dict{Symbol,String}(Symbol("ranges")=>"Vector{IoK8sApiPolicyV1beta1IDRange}", Symbol("rule")=>"String")
Base.propertynames(::Type{ IoK8sApiPolicyV1beta1RunAsUserStrategyOptions }) = collect(keys(_property_map_IoK8sApiPolicyV1beta1RunAsUserStrategyOptions))
Swagger.property_type(::Type{ IoK8sApiPolicyV1beta1RunAsUserStrategyOptions }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiPolicyV1beta1RunAsUserStrategyOptions[name]))}
Swagger.field_name(::Type{ IoK8sApiPolicyV1beta1RunAsUserStrategyOptions }, property_name::Symbol) =  _property_map_IoK8sApiPolicyV1beta1RunAsUserStrategyOptions[property_name]

function check_required(o::IoK8sApiPolicyV1beta1RunAsUserStrategyOptions)
    (getproperty(o, Symbol("rule")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiPolicyV1beta1RunAsUserStrategyOptions }, name::Symbol, val)
end


text:
function code_js(ex, A)
  cinfo = code_lowered(ex, A)[1]
  restructure(cinfo.code)
end

macro code_js(ex)
  isexpr(ex, :call) || error("@code_wasm f(xs...)")
  :(code_js($(esc(ex.args[1])), Tuple{$(map(_ -> Any, ex.args[2:end])...)}))
end


text:
using SHA, BinaryBuilder, Pkg, Pkg.BinaryPlatforms, Pkg.Artifacts
using BinaryBuilder: CompilerShard, BinaryBuilderBase

host_platform = Linux(:x86_64; libc=:musl)

# Test if something is older than a reference, or doesn't exist
function is_outdated(test, reference)
    if !isfile(test)
        return true
    end
    return stat(test).mtime < stat(reference).mtime
end

function unpacked_to_squashfs(unpacked_hash::Base.SHA1, name, version; platform=host_platform, target=nothing)
    path = artifact_path(unpacked_hash)
    squashfs_cs = CompilerShard(name, version, platform, :squashfs; target=target)
    art_name = BinaryBuilderBase.artifact_name(squashfs_cs)
    squash_hash = create_artifact() do target_dir
        target_squashfs = joinpath(target_dir, art_name)
        success(`mksquashfs $(path) $(target_squashfs) -force-uid 0 -force-gid 0 -comp xz -b 1048576 -Xdict-size 100% -noappend`)
    end
    @info("$(art_name) hash: $(bytes2hex(squash_hash.bytes))")
    return squash_hash
end

# Generate artifacts for a given path, creating both an unpacked version and a squashfs version
function generate_artifacts(path::AbstractString, name, version; platform=host_platform, target=nothing)
    if !isdir(path)
        error("$path is not a directory!")
    end

    # First, the unpacked version
    squashfs_cs = CompilerShard(name, version, platform, :unpacked; target=target)
    art_name = BinaryBuilderBase.artifact_name(squashfs_cs)
    unpacked_hash = create_artifact() do target
        rm(target; force=true, recursive=true)
        cp(path, target)
    end
    @info("$(art_name) hash: $(bytes2hex(unpacked_hash.bytes))")

    # Next, the squashfs version
    squashfs_hash = unpacked_to_squashfs(unpacked_hash, name, version; platform=platform, target=target)
    return (unpacked_hash, squashfs_hash)
end

function publish_artifact(repo::AbstractString, tag::AbstractString, hash::Base.SHA1, filename::AbstractString)
    mktempdir() do dir
        tarball_hash = archive_artifact(hash, joinpath(dir, "$(filename).tar.gz"))
        if repo == "local"
            @info "Skipping upload to GitHub because local build was requested!"
        else
            BinaryBuilder.upload_to_github_releases(repo, tag, dir)
        end
        return tarball_hash
    end
end

function get_next_shard_tag(cs)
    artifacts_toml = joinpath(dirname(dirname(pathof(BinaryBuilderBase))), "Artifacts.toml")
    meta = artifact_meta(BinaryBuilderBase.artifact_name(cs), artifacts_toml; platform=cs.host)
    if meta === nothing || !haskey(meta, "download")
        return "$(cs.name)-v$(cs.version)"
    end

    last_tag = basename(dirname(meta["download"][1]["url"]))
    if !startswith(last_tag, "$(cs.name)-v")
        error("Invalid last tag: $(last_tag)")
    end

    last_version = VersionNumber(last_tag[length(cs.name)+2:end])
    build_number = 0
    if isa(last_version.build, Tuple{<:UInt})
        build_number = last_version.build[1] + 1
    end
    next_version = VersionNumber(
        last_version.major,
        last_version.minor,
        last_version.patch,
        (),
        (build_number,),
    )
    return "$(cs.name)-v$(next_version)"
end

function upload_compiler_shard(repo, name, version, hash, archive_type; platform=host_platform, target=nothing)
    cs = CompilerShard(name, version, platform, archive_type; target=target)
    tag = get_next_shard_tag(cs)
    filename = BinaryBuilderBase.artifact_name(cs)
    tarball_hash = publish_artifact(repo, tag, hash, filename)

    return [
        ("https://github.com/$(repo)/releases/download/$(tag)/$(filename).tar.gz", tarball_hash),
    ]
end

# Given a name and the unpacked/squashfs hashes, slip this new compiler shard into
# the Artifacts.toml that belongs to BB
function insert_compiler_shard(name, version, hash, archive_type; platform=host_platform, download_info = nothing, target=nothing)
    cs = CompilerShard(name, version, platform, archive_type; target=target)
    artifacts_toml = joinpath(dirname(dirname(pathof(BinaryBuilderBase))), "Artifacts.toml")
    bind_artifact!(artifacts_toml, BinaryBuilderBase.artifact_name(cs), hash; platform=platform, download_info=download_info, lazy=true, force=true)
end

function upload_and_insert_shards(repo, name, version, unpacked_hash, squashfs_hash, platform; target=nothing)
    # Upload them both to GH releases on Yggdrasil
    unpacked_dl_info = upload_compiler_shard(repo, name, version, unpacked_hash, :unpacked; platform=platform, target=target)
    squashfs_dl_info = upload_compiler_shard(repo, name, version, squashfs_hash, :squashfs; platform=platform, target=target)

    # Insert these final versions into BB
    insert_compiler_shard(name, version, unpacked_hash, :unpacked; download_info=unpacked_dl_info, platform=platform, target=target)
    insert_compiler_shard(name, version, squashfs_hash, :squashfs; download_info=squashfs_dl_info, platform=platform, target=target)
end

function upload_and_insert_shards(repo, name, version, build_info; target=nothing)
    for platform in keys(build_info)
        unpacked_hash = build_info[platform][3]
        squashfs_hash = unpacked_to_squashfs(unpacked_hash, name, version; platform=platform, target=target)

        upload_and_insert_shards(repo, name, version, unpacked_hash, squashfs_hash, platform; target=target)
    end
end

macro flag(name)
    quote
        if isa($(esc(name)), AbstractString)
            $(esc(name))
        elseif $(esc(name)) === true
            $(string("--", name))
        else
            ""
        end
    end
end
function jrun(script::String, args::String...; verbose=verbose, debug=debug, deploy=false, register=false)
    run(`$(Base.julia_cmd()) --color=yes $(script) $(filter(x -> !isempty(x), [@flag(verbose), @flag(debug), @flag(deploy), @flag(register), args...]))`)
end

import BinaryBuilder: build_tarballs
function build_tarballs(project_path::String, args::String...; deploy=false, register=false)
    @info("Building $(project_path)...")
    cd(project_path) do
        jrun("build_tarballs.jl", args...; deploy=deploy, register=register)
    end
end

function find_deploy_arg(ARGS)
    dargs = ARGS[findall(arg->startswith(arg, "--deploy"), ARGS)]
    if length(dargs) > 1
        error("More than one deploy argument. Usage: --deploy|--deploy=local")
    end

    # No deployment
    if length(dargs) == 0
        return (ARGS, false, "")
    end

    ndARGS = filter(arg->!startswith(arg, "--deploy"), ARGS)
    if dargs[] == "--deploy"
        return (ndARGS, true, "JuliaPackaging/Yggdrasil")
    elseif dargs[] == "--deploy=local"
        return (ndARGS, true, "local")
    else
        error("--deploy argument must be --deploy or --deploy=local")
    end
end


text:
# Autogenerated wrapper script for SDL2_net_jll for aarch64-apple-darwin
export libSDL2_net

using SDL2_jll
JLLWrappers.@generate_wrapper_header("SDL2_net")
JLLWrappers.@declare_library_product(libSDL2_net, "@rpath/libSDL2_net-2.0.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header(SDL2_jll)
    JLLWrappers.@init_library_product(
        libSDL2_net,
        "lib/libSDL2_net.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using CSyntax.CRef
using Test

function cref_basic(x, y, z)
    z[] += 1
    return x, y, z
end

mutable struct FooCRef
    x
    y
end

@testset "CRef" begin
    x, y, z = 1, 2, 3
    @cref cref_basic(x, y, &z)
    @test z == 4

    f = FooCRef(10,20)
    @cref cref_basic(x, y, &f.x)
    @test f.x == 11
end


text:
module GPUCompiler

using LLVM
using LLVM.Interop

using DataStructures

using TimerOutputs

using Libdl

const to = TimerOutput()

timings() = (TimerOutputs.print_timer(to); println())

enable_timings() = (TimerOutputs.enable_debug_timings(GPUCompiler); return)

include("utils.jl")

# compiler interface and implementations
include("interface.jl")
include("error.jl")
include("native.jl")
include("ptx.jl")
include("gcn.jl")
include("spirv.jl")

include("runtime.jl")

# compiler implementation
include("irgen.jl")
include("optim.jl")
include("validation.jl")
include("rtlib.jl")
include("mcgen.jl")
include("debug.jl")
include("driver.jl")

# other reusable functionality
include("cache.jl")
include("execution.jl")
include("reflection.jl")

function __init__()
    TimerOutputs.reset_timer!(to)
    InitializeAllTargets()
    InitializeAllTargetInfos()
    InitializeAllAsmPrinters()
    InitializeAllAsmParsers()
    InitializeAllTargetMCs()

    return
end

end # module


text:
using Statistics
using Plots
using FFTW
using Optim
using Images, FileIO, ImageIO
using Printf
using Revise
using Profile
using LinearAlgebra
using JLD2
using Random
using Distributions
using FITSIO
using LineSearches
using Flux
using StatsBase
using Weave
using SparseArrays
using Distances

push!(LOAD_PATH, pwd()*"/../../../main")
using DHC_2DUtils
push!(LOAD_PATH, pwd()*"/../../../scratch_NM")
using Deriv_Utils_New
using Data_Utils
using Visualization
using ReconFuncs
include("../../../main/compute.jl")

#Functions
function J_hashindices(J_values, fhash)
    jindlist = []
    for jval in J_values
        push!(jindlist, findall(fhash["J_L"][:, 1].==jval))
    end
    return vcat(jindlist'...)
end

function J_S1indices(J_values, fhash)
    #Assumes this is applied to an object of length 2+Nf+Nf^2 or 2+Nf
    return J_hashindices(J_values, fhash) .+ 2
end

function compareS1(true_img, init_img, recon_img, filter_hash, dhc_args, tlist)
    JS1ind = J_S1indices([0, 1, 2, 3], filter_hash)
    strue = DHC_compute_wrapper(true_img, filter_hash; dhc_args...)
    s1true = strue[JS1ind]
    clim = (minimum(s1true), maximum(s1true))
    s1gaussian = DHC_compute_wrapper(gprecon, filter_hash; dhc_args...)[JS1ind]
    s1recon = DHC_compute_wrapper(recon_img, filter_hash; dhc_args...)[JS1ind]
    p1 = heatmap(s1true, title=tlist[1], clim=clim)
    p2 = heatmap(s1gaussian, title=tlist[2], clim=clim)
    p3 = heatmap(s1recon, title=tlist[3], clim=clim)
    p = plot(p1, p2, p3)
end

function plot_panel(true_img, gpinit, recon_img)
    clim = (minimum(true_img), maximum(true_img))
    kbins= convert(Array{Float64}, collect(1:32))
    JS1ind = J_S1indices([0, 1, 2, 3], filter_hash)
    true_ps = Data_Utils.calc_1dps(apodizer(true_img), kbins)
    gpps = Data_Utils.calc_1dps(apodizer(gpinit), kbins)
    recps = Data_Utils.calc_1dps(apodizer(recon_img), kbins)
    p1 = heatmap(recon_img, title="Recon", clim=clim)
    p3 = plot(log.(kbins), log.(true_ps), label="True")
    plot!(log.(kbins), log.(recps), label="Recon")
    plot!(log.(kbins), log.(gpps), label="Init")
    plot!(title="P(k)")
    xlabel!("lnk")
    ylabel!("lnP(k)")
    p2 = heatmap(true_img, title="True Img", clim=clim)
    p4= heatmap(gpinit, title="GPInit Img", clim=clim)

    residual = recon_img- true_img
    rlims = (minimum(residual), maximum(residual))
    p5 = heatmap(residual, title="Residual: Recon - True", clims=rlims, c=:bwr)
    p6 = heatmap(gpinit- true_img, title="Residual: GPInit - True", clims=rlims, c=:bwr)

    struesel = Data_Utils.fnlog(DHC_compute_wrapper(true_img, filter_hash, norm=false; dhc_args...))
    ssmoothsel = Data_Utils.fnlog(DHC_compute_wrapper(gpinit, filter_hash, norm=false; dhc_args...))
    sreconsel = Data_Utils.fnlog(DHC_compute_wrapper(recon_img, filter_hash, norm=false; dhc_args...))
    slims = (minimum(struesel[JS1ind]), maximum(struesel[JS1ind]))
    cg = cgrad([:blue, :white, :red])
    truephi, trueomg = round(struesel[2+filter_hash["phi_index"]], sigdigits=3), round(struesel[2+filter_hash["Omega_index"]], sigdigits=3)
    reconphi, reconomg = round(sreconsel[2+filter_hash["phi_index"]], sigdigits=3), round(sreconsel[2+filter_hash["Omega_index"]], sigdigits=3)
    smoothphi, smoothomg = round(ssmoothsel[2+filter_hash["phi_index"]], sigdigits=3), round(ssmoothsel[2+filter_hash["Omega_index"]], sigdigits=3)

    p7 = heatmap(struesel[JS1ind], title="True Coeffs ϕ=" * string(truephi) * "Ω=" * string(trueomg) , clims=slims, c=cg)
    p8 = heatmap(sreconsel[JS1ind], title="Recon Coeffs ϕ=" * string(reconphi) * "Ω=" * string(reconomg), clims=slims, c=cg)
    p9 = heatmap(ssmoothsel[JS1ind], title="GP Init ϕ=" * string(smoothphi) * "Ω=" * string(smoothomg), clims=slims, c=cg)
    p10 = heatmap(zeros(Nx, Nx))
    p = plot(p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, layout = (5, 2), size=(1800, 2400))
end

##Common
loaddf = load("../../../scratch_NM/StandardizedExp/Nx64/Data_1000.jld2")
trueimg = loaddf["true_img"]
Nx=64
pixall = sample(1:Nx^2, Integer(round(0.1*4096)), replace = false)
flatmask = falses(Nx*Nx)
flatmask[pixall] .= true
pixmask = reshape(flatmask, (Nx, Nx))

filter_hash = fink_filter_hash(1, 8, nx=Nx, t=1, wd=1, Omega=true)
dhc_args = Dict(:doS2=>false, :doS20=>true, :apodize=>false, :iso=>false)
strue = DHC_compute_wrapper(trueimg, filter_hash; dhc_args...)

#targets
Nf = length(filter_hash["filt_index"])
coeffmask = falses(2+Nf+Nf^2)
coeffmask[3+Nf:end] .= Diagonal(trues(Nf, Nf))[:]
starget = log.(strue[coeffmask])
scovinv = I
img_guess = fill(mean(trueimg), (Nx, Nx))
img_guess[pixmask] .= trueimg[pixmask]

#GP starting point
gr1d = collect(1:Nx)
xmat = reshape(repeat(gr1d, outer=[Nx]), (Nx, Nx))
ymat = xmat'
xfl = xmat[:]
yfl = ymat[:]
posmat = hcat(xfl, yfl)'
pwise_full = pairwise(Distances.Euclidean(), posmat, dims=2)
function SqExpBasisKernel(d; length=1.0, scale=1.0)
    return scale.^2 * exp(-d.^2/(2* (length.^2)))
end

covmat_full = map((x->SqExpBasisKernel(x; length=5.0, scale=0.5)), pwise_full) #L, scale, sort of optimized
cov_ss = covmat_full[flatmask, flatmask] + 0.01*I
testmask = (x->!x).(flatmask)
cov_sv = covmat_full[flatmask, testmask]
cov_vs = covmat_full[testmask, flatmask]

pred_gaussian = mean(trueimg) .+ cov_vs * inv(cov_ss) * (reshape(trueimg[flatmask], (length(trueimg[flatmask]), 1)) .- mean(trueimg))
gprecon = zeros(Nx^2)
gprecon[flatmask] .= trueimg[flatmask]
gprecon[testmask] .= pred_gaussian[:]
gprecon = reshape(gprecon, (Nx, Nx))
heatmap(gprecon)
using Feather
my_df = Dict([("field", trueimg), ("source_mask", pixmask), ("gprecon_julia", gprecon)])
using HDF5
h5open("fieldsrcs.h5", "w") do file
    write(file, "field", trueimg)
    #write(file, "source_mask", pixmask)
    write(file, "gprecon_julia", gprecon)
end
save(my_df, ".hdf5")

#
lval2 = 0.0
lval3 = 0.0
optim_settings = Dict([("iterations", 1000), ("norm", false), ("minmethod", ConjugateGradient())])
func_specific_params = Dict([(:reg_input=> img_guess), (:coeff_mask1=> coeffmask), (:target1=>starget), (:invcov1=>scovinv), (:coeff_mask2=> coeffmask), (:target2=>starget), (:invcov2=>scovinv), (:lambda2=>lval2), (:lambda3=>lval3), (:pixmask=>pixmask), (:func=>Data_Utils.fnlog), (:dfunc=>Data_Utils.fndlog)])

res, recon_img = image_recon_derivsum_custom(img_guess, filter_hash, dhc_args, ReconFuncs.Loss3Gaussian_transformed, ReconFuncs.dLoss3Gaussian_transformed_cp!; optim_settings=optim_settings, func_specific_params)
#coeff_mask1=nothing, target1=nothing, invcov1=nothing, reg_input=nothing, coeff_mask2=nothing, target2=nothing, invcov2=nothing, func=nothing, dfunc=nothing, lambda2=nothing, lambda3=nothing)
compareS1(loaddf["true_img"], img_guess, recon_img, filter_hash, dhc_args, ["True", "InitPointswMean", "Recon"])
heatmap(recon_img)

#Start with gprecon, no regularization
lval2 = 0.0
lval3 = 0.0
optim_settings = Dict([("iterations", 1000), ("norm", false), ("minmethod", ConjugateGradient())])
func_specific_params = Dict([(:reg_input=> img_guess), (:coeff_mask1=> coeffmask), (:target1=>starget), (:invcov1=>scovinv), (:coeff_mask2=> coeffmask), (:target2=>starget), (:invcov2=>scovinv), (:lambda2=>lval2), (:lambda3=>lval3), (:pixmask=>pixmask), (:func=>Data_Utils.fnlog), (:dfunc=>Data_Utils.fndlog)])

res, recon_img = image_recon_derivsum_custom(gprecon, filter_hash, dhc_args, ReconFuncs.Loss3Gaussian_transformed, ReconFuncs.dLoss3Gaussian_transformed_cp!; optim_settings=optim_settings, func_specific_params)
heatmap(recon_img)
compareS1(loaddf["true_img"], gprecon, recon_img, filter_hash, dhc_args, ["True", "GPRecon", "S2R-Recon"])
println("GP Recon Mean Abs Frac", mean(abs.((gprecon .- trueimg)./trueimg)))
println("S2r Recon Mean Abs Frac", mean(abs.((recon_img .- trueimg)./trueimg)))
println("GP Recon MSE", mean((gprecon .- trueimg).^2))
println("S2r Recon MSE", mean((recon_img .- trueimg).^2))
plot_panel(trueimg, gprecon, recon_img)
#decent-ish but not better than without

#Starting with the GP recon, regularizing with respect to the mean everywhere and the point values at those places: BAD
#Gets the coefficients right? Worse than without reg
lval2 = 0.0
lval3 = 1.0
optim_settings = Dict([("iterations", 1000), ("norm", false), ("minmethod", ConjugateGradient())])
func_specific_params = Dict([(:reg_input=> img_guess), (:coeff_mask1=> coeffmask), (:target1=>starget), (:invcov1=>scovinv), (:coeff_mask2=> coeffmask), (:target2=>starget), (:invcov2=>scovinv), (:lambda2=>lval2), (:lambda3=>lval3), (:pixmask=>pixmask), (:func=>Data_Utils.fnlog), (:dfunc=>Data_Utils.fndlog)])

res, recon_img = image_recon_derivsum_custom(gprecon, filter_hash, dhc_args, ReconFuncs.Loss3Gaussian_transformed, ReconFuncs.dLoss3Gaussian_transformed_cp!; optim_settings=optim_settings, func_specific_params)
heatmap(recon_img)
println("GP Recon Mean Abs Frac ", mean(abs.((gprecon .- trueimg)./trueimg)))
println("S2r Recon Mean Abs Frac ", mean(abs.((recon_img .- trueimg)./trueimg)))
println("GP Recon MSE ", mean((gprecon .- trueimg).^2))
println("S2r Recon MSE ", mean((recon_img .- trueimg).^2))

#Starting with the GP recon, regularizing with respect to GP recon at those places
lval2 = 0.0
lval3 = 1.0
optim_settings = Dict([("iterations", 1000), ("norm", false), ("minmethod", ConjugateGradient())])
func_specific_params = Dict([(:reg_input=> gprecon), (:coeff_mask1=> coeffmask), (:target1=>starget), (:invcov1=>scovinv), (:coeff_mask2=> coeffmask), (:target2=>starget), (:invcov2=>scovinv), (:lambda2=>lval2), (:lambda3=>lval3), (:pixmask=>pixmask), (:func=>Data_Utils.fnlog), (:dfunc=>Data_Utils.fndlog)])

res, recon_img = image_recon_derivsum_custom(gprecon, filter_hash, dhc_args, ReconFuncs.Loss3Gaussian_transformed, ReconFuncs.dLoss3Gaussian_transformed_cp!; optim_settings=optim_settings, func_specific_params)
println("GP Recon Mean Abs Frac ", mean(abs.((gprecon .- trueimg)./trueimg)))
println("S2r Recon Mean Abs Frac ", mean(abs.((recon_img .- trueimg)./trueimg)))
println("GP Recon MSE ", mean((gprecon .- trueimg).^2))
println("S2r Recon MSE ", mean((recon_img .- trueimg).^2))
plot_panel(trueimg, gprecon, recon_img)

#All S2R
#Starting with the GP recon, regularizing with respect to GP recon at those places
coeffmask = falses(2+Nf+Nf^2)
coeffmask[3+Nf:end] .= triu(trues(Nf, Nf))[:]
starget = log.(strue[coeffmask])
lval2 = 0.0
lval3 = 1.0
optim_settings = Dict([("iterations", 1000), ("norm", false), ("minmethod", ConjugateGradient())])
func_specific_params = Dict([(:reg_input=> gprecon), (:coeff_mask1=> coeffmask), (:target1=>starget), (:invcov1=>scovinv), (:coeff_mask2=> coeffmask), (:target2=>starget), (:invcov2=>scovinv), (:lambda2=>lval2), (:lambda3=>lval3), (:pixmask=>pixmask), (:func=>Data_Utils.fnlog), (:dfunc=>Data_Utils.fndlog)])

res, recon_img = image_recon_derivsum_custom(gprecon, filter_hash, dhc_args, ReconFuncs.Loss3Gaussian_transformed, ReconFuncs.dLoss3Gaussian_transformed_cp!; optim_settings=optim_settings, func_specific_params)
println("GP Recon Mean Abs Frac ", mean(abs.((gprecon .- trueimg)./trueimg)))
println("S2r Recon Mean Abs Frac ", mean(abs.((recon_img .- trueimg)./trueimg)))
println("GP Recon MSE ", mean((gprecon .- trueimg).^2))
println("S2r Recon MSE ", mean((recon_img .- trueimg).^2))
p = plot_panel(trueimg, gprecon, recon_img)
savefig(p, "../../../scratch_NM/Conpix/6-10/Case2_alls2r.png")

#Pick the wrong coefficient target: Can I take a random dust image's vector and make it look like that?
anotherdict = load("../../../scratch_NM/StandardizedExp/Nx64/Data_10.jld2")
heatmap(anotherdict["true_img"])
another_s = DHC_compute_wrapper(anotherdict["true_img"], filter_hash; dhc_args...)

coeffmask = falses(2+Nf+Nf^2)
coeffmask[3+Nf:end] .= triu(trues(Nf, Nf))[:]
starget = log.(another_s[coeffmask])
lval2 = 0.0
lval3 = 1.0
optim_settings = Dict([("iterations", 1000), ("norm", false), ("minmethod", ConjugateGradient())])
func_specific_params = Dict([(:reg_input=> gprecon), (:coeff_mask1=> coeffmask), (:target1=>starget), (:invcov1=>scovinv), (:coeff_mask2=> coeffmask), (:target2=>starget), (:invcov2=>scovinv), (:lambda2=>lval2), (:lambda3=>lval3), (:pixmask=>pixmask), (:func=>Data_Utils.fnlog), (:dfunc=>Data_Utils.fndlog)])

res, recon_img = image_recon_derivsum_custom(gprecon, filter_hash, dhc_args, ReconFuncs.Loss3Gaussian_transformed, ReconFuncs.dLoss3Gaussian_transformed_cp!; optim_settings=optim_settings, func_specific_params)
println("GP Recon Mean Abs Frac ", mean(abs.((gprecon .- trueimg)./trueimg)))
println("S2r Recon Mean Abs Frac ", mean(abs.((recon_img .- trueimg)./trueimg)))
println("GP Recon MSE ", mean((gprecon .- trueimg).^2))
println("S2r Recon MSE ", mean((recon_img .- trueimg).^2))
p = plot_panel(trueimg, gprecon, recon_img)

#All S2R, Norm
#Starting with the GP recon, regularizing with respect to GP recon at those places
dhc_args = Dict(:doS2=>false, :doS20=>true, :apodize=>false, :iso=>false, :norm=>true)

coeffmask = falses(2+Nf+Nf^2)
coeffmask[3+Nf:end] .= triu(trues(Nf, Nf))[:]
starget = log.(strue[coeffmask])
lval2 = 0.0
lval3 = 1.0
optim_settings = Dict([("iterations", 1000), ("norm", false), ("minmethod", ConjugateGradient())])

func_specific_params = Dict([(:reg_input=> gprecon), (:coeff_mask1=> coeffmask), (:target1=>starget), (:invcov1=>scovinv), (:coeff_mask2=> coeffmask), (:target2=>starget), (:invcov2=>scovinv), (:lambda2=>lval2), (:lambda3=>lval3),
    (:pixmask=>pixmask), (:func=>Data_Utils.fnlog), (:dfunc=>Data_Utils.fndlog)])

res, recon_img = image_recon_derivsum_custom(gprecon, filter_hash, dhc_args, ReconFuncs.Loss3Gaussian_transformed, ReconFuncs.dLoss3Gaussian_transformed_cp!; optim_settings=optim_settings, func_specific_params)
println("GP Recon Mean Abs Frac ", mean(abs.((gprecon .- trueimg)./trueimg)))
println("S2r Recon Mean Abs Frac ", mean(abs.((recon_img .- trueimg)./trueimg)))
println("GP Recon MSE ", mean((gprecon .- trueimg).^2))
println("S2r Recon MSE ", mean((recon_img .- trueimg).^2))
p = plot_panel(trueimg, gprecon, recon_img)
#savefig(p, "../../../scratch_NM/Conpix/6-10/Case2_alls2r.png")

#Resave coefficients
coeffdict = load("../../SavedCovMats/reg_apd_noiso_logcoeff.jld2")
using HDF5
sfdall = readsfd_fromsrc("../../data/dust10000.fits", 64, logbool=false)
h5open("../../SavedCovMats/reg_apd_noiso_logcoeffs_img.h5", "w") do file
    write(file, "coeffdbn", coeffdict["dbncoeffs"])
    write(file, "images", sfdall)
end

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
using LanguageServer, LanguageServer.SymbolServer

function main(args)
    if length(args) != 4
        error("Invalid number of arguments passed to julia language server.")
    end

    conn = stdout
    (outRead, outWrite) = redirect_stdout()

    if args[2] == "--debug=yes"
        ENV["JULIA_DEBUG"] = "all"
    elseif args[2] != "--debug=no"
        error("Invalid argument passed.")
    end

    symserver_store_path = joinpath(ARGS[4], "symbolstorev2")

    if !ispath(symserver_store_path)
        mkpath(symserver_store_path)
    end

    server = LanguageServerInstance(
        stdin,
        conn,
        args[1],
        args[3],
        nothing,
        symserver_store_path
    )

    run(server)
end

if abspath(PROGRAM_FILE) == @__FILE__
    main(ARGS)
end

max_stars_repo_name: wangl-cc/coc-julia

max_issues_repo_name: wangl-cc/coc-julia

max_forks_repo_name: wangl-cc/coc-julia


text:
using RainFARM
using Test
using Statistics

# write your own tests here

print("Testing main RainFARM function and aggregation:\n")
print("-----------------------------------------------\n")
nt = 2; nf = 8; ns = 64; nas = 8
prf = rand(ns, ns ,nt)
prf[2,2,1] = NaN
prl = agg(prf, nas, nt)
ww = ones(64, 64)
print("Testing fglob=false, fsmooth=false\n")
pr = rainfarm(prl, 1.7, nf, 1.; fglob=false, fsmooth=false, verbose=false)
pra = agg(pr, nas, nt)
eps1 = Statistics.mean((pra-prl).^2)
@test eps1 < 1e-20
print("Testing fglob=true, fsmooth=false, weights\n")
pr = rainfarm(prl, 1.7, nf, ww; fglob=true, fsmooth=false, verbose=false)
pra = agg(pr, nas, nt)
eps1 = (Statistics.mean(pra[:,:,1])-Statistics.mean(prl[:,:,1])).^2
@test eps1 < 1e-20
print("Testing fglob=false, fsmooth=true\n")
pr = rainfarm(prl, 1.7, nf, 1.; fglob=false, fsmooth=true, verbose=false)
pra = agg(pr, nas, nt)
eps1 = (Statistics.mean(pra[:,:,1])-Statistics.mean(prl[:,:,1])).^2
@test eps1 < 0.005
print("\nTesting auxiliary functions:\n")
print("----------------------------\n")
print("Testing fitslope\n")
fx = (1.:10.).^-2
sx = fitslope(fx)
eps1 = abs(sx-1)
@test eps1 < 1e-8
print("Testing lon_lat_fine\n")
lon = 0:7; lat = 0:7; nf = 2
(lon_f, lat_f) = lon_lat_fine(lon, lat, nf)
lon_c = collect(-0.25:0.5:7.25)
eps1 = sum((lon_c-lon_f).^2)
@test eps1 < 1e-20


text:
"""
    AlternatingUpdateSelector

Selects the composite update and manifold update every other iteration of the
`PartlySmoothOptimizer` method.
"""
struct AlternatingUpdateSelector <: AbstractUpdateSelector
end

function select_update!(::AlternatingUpdateSelector, state, optimizer, pb)
    state.previous_update = state.selected_update

    state.selected_update = optimizer.wholespace_update
    if state.previous_update == optimizer.wholespace_update && manifold_dimension(state.M)>0
        state.selected_update = optimizer.manifold_update
    end

    return
end

Base.summary(::AlternatingUpdateSelector) = "Alternate"


text:
using LazySets: STAR

for N in [Float64, Float32, Rational{Int}]

    # constructor with basis matrix
    S = Star(N[3, 3], N[1 0; 0 1], BallInf(N[0, 0], N(1)))

    # constructor vector-of-vectors basis
    W = Star(N[3, 3], [N[1, 0], N[0, 1]], BallInf(N[0, 0], N(1)))

    @test S ≈ W

    # getter functions
    @test center(S) == N[3, 3]
    @test basis(S) == N[1 0; 0 1]
    @test predicate(S) == BallInf(N[0, 0], N(1))
    @test dim(S) == 2

    # support function
    @test ρ(N[1, 0], S) == N(4)
    @test -ρ(N[-1, 0], S) == N(2)
    @test ρ(N[0, 1], S) == N(4)
    @test -ρ(N[0, -1], S) == N(2)

    # support vector
    @test σ(N[1, 1], S) == N[4, 4]
    @test σ(N[-1, -1], S) == N[2, 2]
    @test σ(N[1, -1], S) == N[4, 2]
    @test σ(N[-1, 1], S) == N[2, 4]

    # conversion from a polyhedral set
    P = HPolygon([HalfSpace(N[1, 1], N(-0.5)),
                  HalfSpace(N[-1.5, 0.5], N(4)),
                  HalfSpace(N[-0.5, -0.05], N(1.5)),
                  HalfSpace(N[1.5, -0.25], N(-1.5))])
    S = convert(Star, P)
    @test basis(S) == N[1 0; 0 1]
    @test center(S) == N[0, 0]
    @test predicate(S) == P

    # is_polyhedral
    @test is_polyhedral(S)

    # intersection with HalfSpace
    H = HalfSpace(N[0, 1], N(0))
    I = intersection(S, H)
    intersection!(S, H)
    addconstraint!(P, H)
    @test isequivalent(I, P)
    @test isequivalent(S, P)

    # check that we can intersect polyhedra that are axis-aligned
    B = BallInf(N[0, 0], N(1))
    S = Star(N[0, 0], N[1 0; 0 1], B)
    I = intersection(S, H)
    @test isequivalent(I, intersection(B, H))

    # concrete linear map (returns a star)
    M = N[2 0; 0 2]
    lm = linear_map(M, S)
    @test isa(lm, Star)
    @test isequivalent(lm, Star(N[0, 0], M, B))

    # concrete affine map (returns a star)
    v = N[1, 2]
    am = affine_map(M, S, v)
    @test isa(am, Star)
    @test isequivalent(am, Star(v, M, B))
end


text:
using Test
using NuFlux

const DATA_DIR = joinpath(@__DIR__, "../data")

@testset "load fluxtable" begin
    @test !isa(try NuFlux.readfluxfile(joinpath(DATA_DIR, "frj-ally-20-12-solmin.d")) catch e e end, Exception)
end