text:
# https://github.com/jheinen/GR.jl

const gr_linetype = Dict(
  :auto => 1, :solid => 1, :dash => 2, :dot => 3, :dashdot => 4,
  :dashdotdot => -1 )

const gr_markertype = Dict(
  :auto => 1, :none => -1, :ellipse => -1, :rect => -7, :diamond => -13,
  :utriangle => -3, :dtriangle => -5, :pentagon => -14, :hexagon => 3,
  :cross => 2, :xcross => 5, :star5 => 3 )

const gr_halign = Dict(:left => 1, :hcenter => 2, :right => 3)
const gr_valign = Dict(:top => 1, :vcenter => 3, :bottom => 5)

const gr_font_family = Dict(
  "times" => 1, "helvetica" => 5, "courier" => 9, "bookman" => 14,
  "newcenturyschlbk" => 18, "avantgarde" => 22, "palatino" => 26)

function gr_getcolorind(v)
  c = getColor(v)
  return convert(Int, GR.inqcolorfromrgb(c.r, c.g, c.b))
end

function gr_getaxisind(p)
  axis = get(p, :axis, :none)
  if axis in [:none, :left]
    return 1
  else
    return 2
  end
end

function gr_display(plt::Plot{GRPackage}, clear=true, update=true,
                    subplot=[0, 1, 0, 1])
  d = plt.plotargs

  clear && GR.clearws()

  mwidth, mheight, width, height = GR.inqdspsize()
  w, h = d[:size]
  viewport = zeros(4)
  if w > h
    ratio = float(h) / w
    msize = mwidth * w / width
    GR.setwsviewport(0, msize, 0, msize * ratio)
    GR.setwswindow(0, 1, 0, ratio)
    viewport[1] = subplot[1] + 0.1  * (subplot[2] - subplot[1])
    viewport[2] = subplot[1] + 0.95 * (subplot[2] - subplot[1])
    viewport[3] = ratio * (subplot[3] + 0.1  * (subplot[4] - subplot[3]))
    viewport[4] = ratio * (subplot[3] + 0.95 * (subplot[4] - subplot[3]))
  else
    ratio = float(w) / h
    msize = mheight * h / height
    GR.setwsviewport(0, msize * ratio, 0, msize)
    GR.setwswindow(0, ratio, 0, 1)
    viewport[1] = ratio * (subplot[1] + 0.1  * (subplot[2] - subplot[1]))
    viewport[2] = ratio * (subplot[1] + 0.95 * (subplot[2] - subplot[1]))
    viewport[3] = subplot[3] + 0.1  * (subplot[4] - subplot[3])
    viewport[4] = subplot[3] + 0.95 * (subplot[4] - subplot[3])
  end

  extrema = zeros(2, 4)
  num_axes = 1
  cmap = false
  axes_2d = true

  for axis = 1:2
    xmin = ymin = typemax(Float64)
    xmax = ymax = typemin(Float64)
    for p in plt.seriesargs
      if axis == gr_getaxisind(p)
        if axis == 2
          num_axes = 2
        end
        if p[:linetype] == :bar
          x, y = 1:length(p[:y]), p[:y]
        elseif p[:linetype] == :ohlc
          x, y = 1:size(p[:y], 1), p[:y]
        elseif p[:linetype] in [:hist, :density]
          x, y = Base.hist(p[:y])
        elseif p[:linetype] in [:heatmap, :hexbin]
          E = zeros(length(p[:x]),2)
          E[:,1] = p[:x]
          E[:,2] = p[:y]
          if isa(p[:nbins], Tuple)
            xbins, ybins = p[:nbins]
          else
            xbins = ybins = p[:nbins]
          end
          cmap = true
          x, y, H = Base.hist2d(E, xbins, ybins)
        else
          if p[:linetype] in [:contour, :surface]
            cmap = true
          end
          if p[:linetype] in [:surface, :wireframe, :path3d, :scatter3d]
            axes_2d = false
          end
          x, y = p[:x], p[:y]
        end
        xmin = min(minimum(x), xmin)
        xmax = max(maximum(x), xmax)
        if p[:linetype] == :ohlc
          for val in y
            ymin = min(val.open, val.high, val.low, val.close, ymin)
            ymax = max(val.open, val.high, val.low, val.close, ymax)
          end
        else
          ymin = min(minimum(y), ymin)
          ymax = max(maximum(y), ymax)
        end
      end
    end
    extrema[axis,:] = [xmin, xmax, ymin, ymax]
  end

  if num_axes == 2 || !axes_2d
    viewport[2] -= 0.05
  end
  if cmap
    viewport[2] -= 0.1
  end
  GR.setviewport(viewport[1], viewport[2], viewport[3], viewport[4])

  scale = d[:scale]
  for axis = 1:num_axes
    xmin, xmax, ymin, ymax = extrema[axis,:]
    if scale & GR.OPTION_X_LOG == 0
      xmin, xmax = GR.adjustlimits(xmin, xmax)
      majorx = 5
      xtick = GR.tick(xmin, xmax) / majorx
    else
      xtick = majorx = 1
    end
    if scale & GR.OPTION_Y_LOG == 0
      ymin, ymax = GR.adjustlimits(ymin, ymax)
      majory = 5
      ytick = GR.tick(ymin, ymax) / majory
    else
      ytick = majory = 1
    end
    if scale & GR.OPTION_FLIP_X == 0
      xorg = (xmin, xmax)
    else
      xorg = (xmax, xmin)
    end
    if scale & GR.OPTION_FLIP_Y == 0
      yorg = (ymin, ymax)
    else
      yorg = (ymax, ymin)
    end

    GR.setwindow(xmin, xmax, ymin, ymax)
    if axis == 1 && haskey(d, :background_color)
      GR.savestate()
      GR.setfillintstyle(GR.INTSTYLE_SOLID)
      GR.setfillcolorind(gr_getcolorind(d[:background_color]))
      GR.fillrect(xmin, xmax, ymin, ymax)
      GR.restorestate()
    end
    GR.setscale(scale)

    if axes_2d
      diag = sqrt((viewport[2] - viewport[1])^2 + (viewport[4] - viewport[3])^2)
      GR.setlinewidth(1)
      charheight = max(0.018 * diag, 0.01)
      GR.setcharheight(charheight)
      ticksize = 0.0075 * diag
      GR.grid(xtick, ytick, 0, 0, majorx, majory)
      if num_axes == 1
        GR.axes(xtick, ytick, xorg[1], yorg[1], majorx, majory, ticksize)
        GR.axes(xtick, ytick, xorg[2], yorg[2], -majorx, -majory, -ticksize)
      elseif axis == 1
        GR.axes(xtick, ytick, xorg[1], yorg[1], majorx, majory, ticksize)
      else
        GR.axes(xtick, ytick, xorg[2], yorg[2], -majorx, majory, -ticksize)
      end
    end
  end

  if get(d, :title, "") != ""
    GR.savestate()
    GR.settextalign(GR.TEXT_HALIGN_CENTER, GR.TEXT_VALIGN_TOP)
    GR.text(0.5, min(ratio, 1), d[:title])
    GR.restorestate()
  end
  if get(d, :xlabel, "") != ""
    GR.savestate()
    GR.settextalign(GR.TEXT_HALIGN_CENTER, GR.TEXT_VALIGN_BOTTOM)
    GR.text(0.5, 0, d[:xlabel])
    GR.restorestate()
  end
  if get(d, :ylabel, "") != ""
    GR.savestate()
    GR.settextalign(GR.TEXT_HALIGN_CENTER, GR.TEXT_VALIGN_TOP)
    GR.setcharup(-1, 0)
    GR.text(0, 0.5 * (viewport[3] + viewport[4]), d[:ylabel])
    GR.restorestate()
  end
  if get(d, :yrightlabel, "") != ""
    GR.savestate()
    GR.settextalign(GR.TEXT_HALIGN_CENTER, GR.TEXT_VALIGN_TOP)
    GR.setcharup(1, 0)
    GR.text(1, 0.5 * (viewport[3] + viewport[4]), d[:yrightlabel])
    GR.restorestate()
  end

  legend = false

  for p in plt.seriesargs
    GR.savestate()
    xmin, xmax, ymin, ymax = extrema[gr_getaxisind(p),:]
    GR.setwindow(xmin, xmax, ymin, ymax)
    if p[:linetype] in [:path, :line, :steppre, :steppost, :sticks, :hline, :vline, :ohlc]
      haskey(p, :linestyle) && GR.setlinetype(gr_linetype[p[:linestyle]])
      haskey(p, :linewidth) && GR.setlinewidth(p[:linewidth])
      haskey(p, :linecolor) && GR.setlinecolorind(gr_getcolorind(p[:linecolor]))
    end
    if p[:linetype] == :path
      if haskey(p, :fillcolor)
        GR.setfillcolorind(gr_getcolorind(p[:fillcolor]))
        GR.setfillintstyle(GR.INTSTYLE_SOLID)
      end
      if p[:fillrange] != nothing
        GR.fillarea([p[:x][1]; p[:x]; p[:x][length(p[:x])]], [p[:fillrange]; p[:y]; p[:fillrange]])
      end
      GR.polyline(p[:x], p[:y])
      legend = true
    end
    if p[:linetype] == :line
      GR.polyline(p[:x], p[:y])
      legend = true
    elseif p[:linetype] in [:steppre, :steppost]
      n = length(p[:x])
      x = zeros(2*n + 1)
      y = zeros(2*n + 1)
      x[1], y[1] = p[:x][1], p[:y][1]
      j = 2
      for i = 2:n
        if p[:linetype] == :steppre
          x[j], x[j+1] = p[:x][i-1], p[:x][i]
          y[j], y[j+1] = p[:y][i],   p[:y][i]
        else
          x[j], x[j+1] = p[:x][i],   p[:x][i]
          y[j], y[j+1] = p[:y][i-1], p[:y][i]
        end
        j += 2
      end
      GR.polyline(x, y)
      legend = true
    elseif p[:linetype] == :sticks
      x, y = p[:x], p[:y]
      for i = 1:length(y)
        GR.polyline([x[i], x[i]], [ymin, y[i]])
      end
      legend = true
    elseif p[:linetype] == :scatter || (p[:markershape] != :none && axes_2d)
      haskey(p, :markercolor) && GR.setmarkercolorind(gr_getcolorind(p[:markercolor]))
      haskey(p, :markershape) && GR.setmarkertype(gr_markertype[p[:markershape]])
      if haskey(d, :markersize)
        if typeof(d[:markersize]) <: Number
          GR.setmarkersize(d[:markersize] / 4.0)
          GR.polymarker(p[:x], p[:y])
        else
          c = p[:markercolor]
          GR.setcolormap(-GR.COLORMAP_GLOWING)
          for i = 1:length(p[:x])
            if isa(c, ColorGradient) && p[:zcolor] != nothing
              ci = round(Int, 1000 + p[:zcolor][i] * 255)
              GR.setmarkercolorind(ci)
            end
            GR.setmarkersize(d[:markersize][i] / 4.0)
            GR.polymarker([p[:x][i]], [p[:y][i]])
          end
        end
      else
        GR.polymarker(p[:x], p[:y])
      end
      legend = true
    elseif p[:linetype] == :bar
      y = p[:y]
      for i = 1:length(y)
        GR.setfillcolorind(gr_getcolorind(p[:fillcolor]))
        GR.setfillintstyle(GR.INTSTYLE_SOLID)
        GR.fillrect(i-0.4, i+0.4, max(0, ymin), y[i])
        GR.setfillcolorind(1)
        GR.setfillintstyle(GR.INTSTYLE_HOLLOW)
        GR.fillrect(i-0.4, i+0.4, max(0, ymin), y[i])
      end
    elseif p[:linetype] in [:hist, :density]
      h = Base.hist(p[:y])
      x, y = float(collect(h[1])), float(h[2])
      for i = 2:length(y)
        GR.setfillcolorind(gr_getcolorind(p[:fillcolor]))
        GR.setfillintstyle(GR.INTSTYLE_SOLID)
        GR.fillrect(x[i-1], x[i], ymin, y[i])
        GR.setfillcolorind(1)
        GR.setfillintstyle(GR.INTSTYLE_HOLLOW)
        GR.fillrect(x[i-1], x[i], ymin, y[i])
      end
    elseif p[:linetype] in [:hline, :vline]
      for xy in p[:y]
        if p[:linetype] == :hline
          GR.polyline([xmin, xmax], [xy, xy])
        else
          GR.polyline([xy, xy], [ymin, ymax])
        end
      end
    elseif p[:linetype] in [:heatmap, :hexbin]
      E = zeros(length(p[:x]),2)
      E[:,1] = p[:x]
      E[:,2] = p[:y]
      if isa(p[:nbins], Tuple)
        xbins, ybins = p[:nbins]
      else
        xbins = ybins = p[:nbins]
      end
      x, y, H = Base.hist2d(E, xbins, ybins)
      counts = round(Int32, 1000 + 255 * H / maximum(H))
      n, m = size(counts)
      GR.setcolormap(GR.COLORMAP_COOLWARM)
      GR.cellarray(xmin, xmax, ymin, ymax, n, m, counts)
      GR.setviewport(viewport[2] + 0.02, viewport[2] + 0.05,
                     viewport[3], viewport[4])
      GR.setspace(0, maximum(counts), 0, 90)
      diag = sqrt((viewport[2] - viewport[1])^2 + (viewport[4] - viewport[3])^2)
      charheight = max(0.016 * diag, 0.01)
      GR.setcharheight(charheight)
      GR.colormap()
    elseif p[:linetype] == :contour
      x, y, z = p[:x], p[:y], p[:z].surf
      zmin, zmax = minimum(z), maximum(z)
      if typeof(p[:levels]) <: Array
        h = p[:levels]
      else
        h = linspace(zmin, zmax, p[:levels])
      end
      GR.setspace(zmin, zmax, 0, 90)
      GR.setcolormap(GR.COLORMAP_COOLWARM)
      GR.contour(x, y, h, reshape(z, length(x) * length(y)), 1000)
      GR.setviewport(viewport[2] + 0.02, viewport[2] + 0.05,
                     viewport[3], viewport[4])
      l = round(Int32, 1000 + (h - minimum(h)) / (maximum(h) - minimum(h)) * 255)
      GR.setwindow(xmin, xmax, zmin, zmax)
      GR.cellarray(xmin, xmax, zmax, zmin, 1, length(l), l)
      ztick = 0.5 * GR.tick(zmin, zmax)
      diag = sqrt((viewport[2] - viewport[1])^2 + (viewport[4] - viewport[3])^2)
      charheight = max(0.016 * diag, 0.01)
      GR.setcharheight(charheight)
      GR.axes(0, ztick, xmax, zmin, 0, 1, 0.005)
    elseif p[:linetype] in [:surface, :wireframe]
      x, y, z = p[:x], p[:y], p[:z].surf
      zmin, zmax = GR.adjustrange(minimum(z), maximum(z))
      GR.setspace(zmin, zmax, 40, 70)
      xtick = GR.tick(xmin, xmax) / 2
      ytick = GR.tick(ymin, ymax) / 2
      ztick = GR.tick(zmin, zmax) / 2
      diag = sqrt((viewport[2] - viewport[1])^2 + (viewport[4] - viewport[3])^2)
      charheight = max(0.018 * diag, 0.01)
      ticksize = 0.01 * (viewport[2] - viewport[1])
      GR.savestate()
      GR.setlinewidth(1)
      GR.grid3d(xtick, 0, ztick, xmin, ymin, zmin, 2, 0, 2)
      GR.grid3d(0, ytick, 0, xmax, ymin, zmin, 0, 2, 0)
      GR.restorestate()
      z = reshape(z, length(x) * length(y))
      if p[:linetype] == :surface
        GR.setcolormap(GR.COLORMAP_COOLWARM)
        GR.gr3.surface(x, y, z, GR.OPTION_COLORED_MESH)
      else
        GR.setfillcolorind(0)
        GR.surface(x, y, z, GR.OPTION_FILLED_MESH)
      end
      GR.setlinewidth(1)
      GR.setcharheight(charheight)
      GR.axes3d(xtick, 0, ztick, xmin, ymin, zmin, 2, 0, 2, -ticksize)
      GR.axes3d(0, ytick, 0, xmax, ymin, zmin, 0, 2, 0, ticksize)
      if cmap
        GR.setviewport(viewport[2] + 0.07, viewport[2] + 0.1,
                       viewport[3], viewport[4])
        GR.colormap()
      end
    elseif p[:linetype] in [:path3d, :scatter3d]
      x, y, z = p[:x], p[:y], p[:z]
      zmin, zmax = GR.adjustrange(minimum(z), maximum(z))
      GR.setspace(zmin, zmax, 40, 70)
      xtick = GR.tick(xmin, xmax) / 2
      ytick = GR.tick(ymin, ymax) / 2
      ztick = GR.tick(zmin, zmax) / 2
      diag = sqrt((viewport[2] - viewport[1])^2 + (viewport[4] - viewport[3])^2)
      charheight = max(0.018 * diag, 0.01)
      ticksize = 0.01 * (viewport[2] - viewport[1])
      GR.savestate()
      GR.setlinewidth(1)
      GR.grid3d(xtick, 0, ztick, xmin, ymin, zmin, 2, 0, 2)
      GR.grid3d(0, ytick, 0, xmax, ymin, zmin, 0, 2, 0)
      GR.restorestate()
      if p[:linetype] == :scatter3d
        haskey(p, :markercolor) && GR.setmarkercolorind(gr_getcolorind(p[:markercolor]))
        haskey(p, :markershape) && GR.setmarkertype(gr_markertype[p[:markershape]])
        for i = 1:length(z)
          px, py = GR.wc3towc(x[i], y[i], z[i])
          GR.polymarker([px], [py])
        end
      else
        haskey(p, :linewidth) && GR.setlinewidth(p[:linewidth])
        GR.polyline3d(x, y, z)
      end
      GR.setlinewidth(1)
      GR.setcharheight(charheight)
      GR.axes3d(xtick, 0, ztick, xmin, ymin, zmin, 2, 0, 2, -ticksize)
      GR.axes3d(0, ytick, 0, xmax, ymin, zmin, 0, 2, 0, ticksize)
    elseif p[:linetype] == :ohlc
      y = p[:y]
      n = size(y, 1)
      ticksize = 0.5 * (xmax - xmin) / n
      for i in 1:n
        GR.polyline([i-ticksize, i], [y[i].open, y[i].open])
        GR.polyline([i, i], [y[i].low, y[i].high])
        GR.polyline([i, i+ticksize], [y[i].close, y[i].close])
      end
    elseif p[:linetype] in [:pie]
      println("TODO: add support for linetype $(p[:linetype])")
    end
    GR.restorestate()
  end

  if d[:legend] && legend
    GR.savestate()
    GR.selntran(0)
    GR.setscale(0)
    w = 0
    i = 0
    for p in plt.seriesargs
      if typeof(p[:label]) <: Array
        i += 1
        lab = p[:label][i]
      else
        lab = p[:label]
      end
      tbx, tby = GR.inqtext(0, 0, lab)
      w = max(w, tbx[3])
    end
    px = viewport[2] - 0.05 - w
    py = viewport[4] - 0.06
    GR.setfillintstyle(GR.INTSTYLE_SOLID)
    GR.setfillcolorind(0)
    GR.fillrect(px - 0.08, px + w + 0.02, py + 0.03, py - 0.03 * length(plt.seriesargs))
    GR.setlinetype(1)
    GR.setlinecolorind(1)
    GR.setlinewidth(1)
    GR.drawrect(px - 0.08, px + w + 0.02, py + 0.03, py - 0.03 * length(plt.seriesargs))
    haskey(d, :linewidth) && GR.setlinewidth(d[:linewidth])
    i = 0
    for p in plt.seriesargs
      if p[:linetype] in [:path, :line, :steppre, :steppost, :sticks]
        haskey(p, :linecolor) && GR.setlinecolorind(gr_getcolorind(p[:linecolor]))
        haskey(p, :linestyle) && GR.setlinetype(gr_linetype[p[:linestyle]])
        GR.polyline([px - 0.07, px - 0.01], [py, py])
      end
      if p[:linetype] == :scatter || p[:markershape] != :none
        haskey(p, :markercolor) && GR.setmarkercolorind(gr_getcolorind(p[:markercolor]))
        haskey(p, :markershape) && GR.setmarkertype(gr_markertype[p[:markershape]])
        if p[:linetype] in [:path, :line, :steppre, :steppost, :sticks]
          GR.polymarker([px - 0.06, px - 0.02], [py, py])
        else
          GR.polymarker([px - 0.06, px - 0.04, px - 0.02], [py, py, py])
        end
      end
      if typeof(p[:label]) <: Array
        i += 1
        lab = p[:label][i]
      else
        lab = p[:label]
      end
      GR.settextalign(GR.TEXT_HALIGN_LEFT, GR.TEXT_VALIGN_HALF)
      GR.text(px, py, lab)
      py -= 0.03
    end
    GR.selntran(1)
    GR.restorestate()
  end

  if haskey(d, :anns)
    GR.savestate()
    for ann in d[:anns]
      x, y, val = ann
      x, y = GR.wctondc(x, y)
      alpha = val.font.rotation
      family = lowercase(val.font.family)
      GR.setcharheight(0.7 * val.font.pointsize / d[:size][2])
      GR.setcharup(sin(val.font.rotation), cos(val.font.rotation))
      if haskey(gr_font_family, family)
        GR.settextfontprec(100 + gr_font_family[family], GR.TEXT_PRECISION_STRING)
      end
      GR.settextcolorind(gr_getcolorind(val.font.color))
      GR.settextalign(gr_halign[val.font.halign], gr_valign[val.font.valign])
      GR.text(x, y, val.str)
    end
    GR.restorestate()
  end

  update && GR.updatews()
end

function gr_display(subplt::Subplot{GRPackage})
  clear = true
  update = false
  l = enumerate(subplt.layout)
  nr = nrows(subplt.layout)
  for (i, (r, c)) in l
    nc = ncols(subplt.layout, r)
    if i == length(l)
      update = true
    end
    subplot = [(c-1)/nc, c/nc, 1-r/nr, 1-(r-1)/nr]
    gr_display(subplt.plts[i], clear, update, subplot)
    clear = false
  end
end

function _create_plot(pkg::GRPackage; kw...)
  isijulia() && GR.inline("svg")
  d = Dict(kw)
  Plot(nothing, pkg, 0, d, Dict[])
end

function _add_series(::GRPackage, plt::Plot; kw...)
  d = Dict(kw)
  push!(plt.seriesargs, d)
  plt
end

function _add_annotations{X,Y,V}(plt::Plot{GRPackage}, anns::AVec{@compat(Tuple{X,Y,V})})
  if haskey(plt.plotargs, :anns)
    append!(plt.plotargs[:anns], anns)
  else
    plt.plotargs[:anns] = anns
  end
end

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

function _before_update_plot(plt::Plot{GRPackage})
end

function _update_plot(plt::Plot{GRPackage}, d::Dict)
  scale = 0
  d[:xscale] == :log10 && (scale |= GR.OPTION_X_LOG)
  d[:yscale] == :log10 && (scale |= GR.OPTION_Y_LOG)
  get(d, :xflip, false) && (scale |= GR.OPTION_FLIP_X)
  get(d, :yflip, false) && (scale |= GR.OPTION_FLIP_Y)
  plt.plotargs[:scale] = scale

  for k in (:title, :xlabel, :ylabel)
    haskey(d, k) && (plt.plotargs[k] = d[k])
  end
end

function _update_plot_pos_size(plt::PlottingObject{GRPackage}, d::Dict)
end

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

function Base.getindex(plt::Plot{GRPackage}, i::Int)
  d = plt.seriesargs[i]
  d[:x], d[:y]
end

function Base.setindex!(plt::Plot{GRPackage}, xy::Tuple, i::Integer)
  d = plt.seriesargs[i]
  d[:x], d[:y] = xy
  plt
end

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

function _create_subplot(subplt::Subplot{GRPackage}, isbefore::Bool)
  true
end

function _expand_limits(lims, plt::Plot{GRPackage}, isx::Bool)
  # TODO: call expand limits for each plot data
end

function _remove_axis(plt::Plot{GRPackage}, isx::Bool)
  # TODO: if plot is inner subplot, might need to remove ticks or axis labels
end

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

function Base.writemime(io::IO, m::MIME"image/png", plt::PlottingObject{GRPackage})
  ENV["GKS_WSTYPE"] = "png"
  gr_display(plt)
  GR.emergencyclosegks()
  write(io, readall("gks.png"))
end

function Base.writemime(io::IO, m::MIME"image/svg+xml", plt::PlottingObject{GRPackage})
  isijulia() || return
  ENV["GKS_WSTYPE"] = "svg"
  gr_display(plt)
  GR.emergencyclosegks()
  write(io, readall("gks.svg"))
end

function Base.display(::PlotsDisplay, plt::Plot{GRPackage})
  gr_display(plt)
end

function Base.display(::PlotsDisplay, plt::Subplot{GRPackage})
  true
end


text:
const DO_NOT_PUSH_TO_THESE_DESTINATIONS = String[
    "MirrorUpdater",
    "MirrorUpdater.jl",
    "UnofficialJuliaMirror/MirrorUpdater",
    "UnofficialJuliaMirror/MirrorUpdater.jl",
    "UnofficialJuliaMirrorSnapshots/MirrorUpdater",
    "UnofficialJuliaMirrorSnapshots/MirrorUpdater.jl",
    "UnofficialJuliaMirrorBot/MirrorUpdater",
    "UnofficialJuliaMirrorBot/MirrorUpdater.jl",
    #
    "Snapshots",
    "Snapshots.jl",
    "UnofficialJuliaMirror/Snapshots",
    "UnofficialJuliaMirror/Snapshots.jl",
    "UnofficialJuliaMirrorSnapshots/Snapshots",
    "UnofficialJuliaMirrorSnapshots/Snapshots.jl",
    "UnofficialJuliaMirrorBot/Snapshots",
    "UnofficialJuliaMirrorBot/Snapshots.jl",
    #
    "RemoveLFS",
    "RemoveLFS.jl",
    "UnofficialJuliaMirror/RemoveLFS",
    "UnofficialJuliaMirror/RemoveLFS.jl",
    "UnofficialJuliaMirrorSnapshots/RemoveLFS",
    "UnofficialJuliaMirrorSnapshots/RemoveLFS.jl",
    "UnofficialJuliaMirrorBot/RemoveLFS",
    "UnofficialJuliaMirrorBot/RemoveLFS.jl",
    ]

ext: jl


text:
using Arpack

export eigenmode, eigenfrequency

"    eigenmode(d::AbstractDevice{1}, pol::Polarization, neff::Number, neigenvalues::Int)"
function eigenmode(d::AbstractDevice{1}, pol::Polarization, neff::Number, neigenvalues::Int)
    (ϵ₀, μ₀, c₀) = normalize_parameters(d);
    ω = d.ω[1];

    Tϵ = sparse(Diagonal((ϵ₀*d.ϵᵣ[:])));
    Tϵxinv = sparse(Diagonal((ϵ₀*grid_average(d.ϵᵣ, x̂)[:]).^-1));

    δxb = δ(x̂, Backward, d.grid);
    δxf = δ(x̂, Forward,  d.grid);

    if pol == TM
        A = ω^2*μ₀*Tϵ + δxf*δxb;
    end
    if pol == TE
        A = ω^2*μ₀*Tϵ + Tϵ*δxf*Tϵxinv*δxb;
    end

    estimatedβ = ω/c₀*neff;

    (β², vectors) = eigs(A, nev=neigenvalues, sigma=estimatedβ^2);

    β = sqrt.(β²);
    return (β, vectors)
end

"    eigenmode(d::AbstractDevice{2}, neff::Number, neigenvalues::Int)"
function eigenmode(d::AbstractDevice{2}, neff::Number, neigenvalues::Int)
    (ϵ₀, μ₀, c₀) = normalize_parameters(d);
    ω = d.ω[1];

    Tϵx    = sparse(Diagonal(ϵ₀*grid_average(d.ϵᵣ, x̂)[:]));
    Tϵxinv = sparse(Diagonal((ϵ₀*grid_average(d.ϵᵣ, x̂)[:]).^-1));
    Tϵy    = sparse(Diagonal(ϵ₀*grid_average(d.ϵᵣ, ŷ)[:]));
    Tϵyinv = sparse(Diagonal((ϵ₀*grid_average(d.ϵᵣ, ŷ)[:]).^-1));
    Tϵzinv = sparse(Diagonal((ϵ₀*grid_average(grid_average(d.ϵᵣ, x̂), ŷ)[:]).^-1));

    δxb = δ(x̂, Backward, d.grid);
    δxf = δ(x̂, Forward,  d.grid);
    δyb = δ(ŷ, Backward, d.grid);
    δyf = δ(ŷ, Forward,  d.grid);

    A = ω^2*μ₀*blockdiag(Tϵy, Tϵx) + blockdiag(Tϵy, Tϵx)*vcat(-δyf, δxf)*Tϵzinv*hcat(-δyb, δxb) + vcat(δxb, δyb)*hcat(δxf, δyf);

    estimatedβ = ω/c₀*neff;

    (β², vectors) = eigs(A, nev=neigenvalues, sigma=estimatedβ^2);

    β = sqrt.(β²);
    fields = Array{FieldAll}(undef, neigenvalues);
    for i in 1:neigenvalues
        hx = vectors[1:length(d.grid), i];
        hy = vectors[length(d.grid)+1:2*length(d.grid), i];
        hz = 1/(1im*β[i])*(δxf*hx+δyf*hy);
        ex = 1/(1im*ω)*Tϵxinv*(δyb*hz+1im*β[i]*hy);
        ey = 1/(1im*ω)*Tϵyinv*(-1im*β[i]*hx-δxb*hz);
        ez = 1/(1im*ω)*Tϵzinv*(δxb*hy-δyb*hx);
        fields[i] = FieldAll(d.grid, ω, ex, ey, ez, hx, hy, hz);
    end

    return (β, fields)
end

"    eigenfrequency(d::AbstractDevice, pol::Polarization, neigenvalues::Int; which::Symbol=:LM)"
function eigenfrequency(d::AbstractDevice, pol::Polarization, neigenvalues::Int; which::Symbol=:LM)
    (ϵ₀, μ₀, c₀) = normalize_parameters(d);
    ω₀ = d.ω[1];

    (Sxf, Sxb, Syf, Syb) = S_create(d.grid, ω₀);

    δxb = Sxb*δ(x̂, Backward, d.grid);
    δxf = Sxf*δ(x̂, Forward,  d.grid);
    δyb = Syb*δ(ŷ, Backward, d.grid);
    δyf = Syf*δ(ŷ, Forward,  d.grid);

    if pol == TM
        fields = Array{FieldTM}(undef, neigenvalues);

        Tϵᵣ⁻¹ = sparse(Diagonal(1 ./ d.ϵᵣ[:]));
        A = Tϵᵣ⁻¹*δxf*δxb + Tϵᵣ⁻¹*δyf*δyb;

        (ω²μϵ, ez) = eigs(A, nev=neigenvalues, sigma=-ω₀^2*μ₀*ϵ₀, which=which);
        ω = sqrt.(-ω²μϵ/μ₀/ϵ₀);

        for i = 1:neigenvalues
            hxi = -1/1im/ω[i]/μ₀*δyf*ez[:,i];
            hyi = 1/1im/ω[i]/μ₀*δxf*ez[:,i];
            fields[i] = FieldTM(d.grid, ω[i], ez[:,i], hxi, hyi)
        end

        return (ω, fields)
    end
    if pol == TE
        fields = Array{FieldTE}(undef, neigenvalues);

        Tϵx⁻¹ = sparse(Diagonal(1 ./ grid_average(ϵ₀*d.ϵᵣ, x̂)[:]));
        Tϵy⁻¹ = sparse(Diagonal(1 ./ grid_average(ϵ₀*d.ϵᵣ, ŷ)[:]));
        A = δxf*Tϵx⁻¹*δxb + δyf*Tϵy⁻¹*δyb;

        (ω²μ, hz) = eigs(A, nev=neigenvalues, sigma=-ω₀^2*μ₀, which=which);
        ω = sqrt.(-ω²μ/μ₀);

        for i = 1:neigenvalues
            exi = 1/1im/ω[i]*Tϵx⁻¹*δyb*hz[:,i];
            eyi = 1/1im/ω[i]*Tϵy⁻¹*(-δxb*hz[:,i]);
            fields[i] = FieldTE(d.grid, ω[i], hz[:,i], exi, eyi)
        end

        return (ω, fields)
    end
end


text:
immutable Raw <: Format end

# Raw Parser
# ==========

function _get_seqs(io::Union{IO, AbstractString})
    SEQS = ASCIIString[]
    for line in eachline(io)
        push!(SEQS, chomp(line))
    end
    SEQS
end

function parse(io::Union{IO, AbstractString}, format::Type{Raw}, output::Type{Matrix{Residue}};
               deletefullgaps::Bool=true, checkalphabet::Bool=false)
    SEQS = _get_seqs(io)
    _strings_to_msa(SEQS, deletefullgaps, checkalphabet)
end

function parse(io::Union{IO,AbstractString}, format::Type{Raw}, output::Type{AnnotatedMultipleSequenceAlignment};
               generatemapping::Bool=false, deletefullgaps::Bool=true, checkalphabet::Bool=false, keepinserts::Bool=false)
    SEQS = _get_seqs(io)
    annot = Annotations()
    if keepinserts
        _keepinserts!(SEQS, annot)
    end
    if generatemapping
        MSA, MAP = _to_msa_mapping(SEQS)
        setannotfile!(annot, "NCol", string(size(MSA,2)))
        setannotfile!(annot, "ColMap", join(vcat(1:size(MSA,2)), ','))
        for i in 1:length(SEQS)
            setannotsequence!(annot, string(i), "SeqMap", MAP[i])
        end
    else
        MSA = convert(Matrix{Residue}, SEQS)
    end
    msa = AnnotatedMultipleSequenceAlignment(IndexedArray(ASCIIString[ string(i) for i in 1:length(SEQS)]), MSA, annot)
    if checkalphabet
        deletenotalphabetsequences!(msa, SEQS)
    end
    if deletefullgaps
        deletefullgapcolumns!(msa)
    end
    msa
end

parse(io::Union{IO, AbstractString}, format::Type{Raw};
      deletefullgaps::Bool=true, checkalphabet::Bool=false) = parse(io, Raw, Matrix{Residue};
                                       deletefullgaps=deletefullgaps, checkalphabet=checkalphabet)

# Print Raw
# =========

function print(io::IO, msa::Matrix{Residue}, format::Type{Raw})
    for i in 1:nsequences(msa)
        seq = asciisequence(msa, i)
        println(io, seq)
    end
end

print(io::IO, msa::AbstractMultipleSequenceAlignment, format::Type{Raw}) = print(io, msa.msa, Raw)

print(msa::Matrix{Residue}, format::Type{Raw}) = print(STDOUT, msa, Raw)
print(msa::AbstractMultipleSequenceAlignment, format::Type{Raw}) = print(STDOUT, msa.msa, Raw)


text:
"""
```
impulse_responses(m, system)

impulse_responses(system, horizon)
```

Compute impulse responses for a single draw.

### Inputs

- `m::AbstractModel`: model object
- `system::System{S}`: state-space system matrices
- `horizon::Int`: number of periods ahead to forecast

where `S<:AbstractFloat`

### Outputs

- `states::Array{S, 3}`: matrix of size `nstates` x `horizon` x `nshocks` of
  state impulse response functions
- `obs::Array{S, 3}`: matrix of size `nobs` x `horizon` x `nshocks` of
  observable impulse response functions
- `pseudo::Array{S, 3}`: matrix of size `npseudo` x `horizon` x `nshocks` of
  pseudo-observable impulse response functions
"""
function impulse_responses{S<:AbstractFloat}(m::AbstractModel, system::System{S})
    horizon = impulse_response_horizons(m)
    impulse_responses(system, horizon)
end

function impulse_responses{S<:AbstractFloat}(system::System{S}, horizon::Int)
    # Setup
    nshocks      = size(system[:RRR], 2)
    nstates      = size(system[:TTT], 1)
    nobs         = size(system[:ZZ], 1)
    npseudo      = size(system[:ZZ_pseudo], 1)

    states = zeros(S, nstates, horizon, nshocks)
    obs    = zeros(S, nobs,    horizon, nshocks)
    pseudo = zeros(S, npseudo, horizon, nshocks)

    # Set constant system matrices to 0
    system = zero_system_constants(system)

    z0 = zeros(S, nstates)

    for i = 1:nshocks
        # Isolate single shock
        shocks = zeros(S, nshocks, horizon)
        shocks[i, 1] = -sqrt(system[:QQ][i, i]) # a negative 1 s.d. shock

        # Iterate state space forward
        states[:, :, i], obs[:, :, i], pseudo[:, :, i], _ = forecast(system, z0, shocks)
    end

    return states, obs, pseudo
end


text:
module TestPrimitivePrognosticConversion

using CLIMAParameters
using CLIMAParameters.Planet
using Test
using StaticArrays
using UnPack

using ClimateMachine
ClimateMachine.init()
const ArrayType = ClimateMachine.array_type()

const clima_dir = dirname(dirname(pathof(ClimateMachine)))

include(joinpath(clima_dir, "test", "Common", "Thermodynamics", "profiles.jl"))
using ClimateMachine.BalanceLaws
using ClimateMachine.ConfigTypes
using ClimateMachine.VariableTemplates
using ClimateMachine.Thermodynamics
using ClimateMachine.TemperatureProfiles
using ClimateMachine.Atmos
using ClimateMachine.BalanceLaws:
    prognostic_to_primitive!, primitive_to_prognostic!
const BL = BalanceLaws

struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

atol_temperature = 5e-1
atol_energy = cv_d(param_set) * atol_temperature

import ClimateMachine.BalanceLaws: vars_state

# Assign aux.ref_state different than state for general testing
function assign!(aux, bl, nt, st::Auxiliary)
    @unpack p, ρ, e_pot, = nt
    aux.ref_state.ρ = ρ / 2
    aux.ref_state.p = p / 2
    aux.orientation.Φ = e_pot
end

function assign!(state, bl, nt, st::Prognostic, aux)
    @unpack u, v, w, ρ, e_kin, e_pot, T, q_pt = nt
    state.ρ = ρ
    assign!(state, bl, compressibility_model(bl), nt, st, aux)
    state.ρu = SVector(state.ρ * u, state.ρ * v, state.ρ * w)
    param_set = parameter_set(bl)
    state.energy.ρe = state.ρ * total_energy(param_set, e_kin, e_pot, T, q_pt)
    assign!(state, bl, moisture_model(bl), nt, st)
end
assign!(state, bl, moisture::DryModel, nt, ::Prognostic) = nothing
assign!(state, bl, moisture::EquilMoist, nt, ::Prognostic) =
    (state.moisture.ρq_tot = state.ρ * nt.q_pt.tot)
function assign!(state, bl, moisture::NonEquilMoist, nt, ::Prognostic)
    state.moisture.ρq_tot = state.ρ * nt.q_pt.tot
    state.moisture.ρq_liq = state.ρ * nt.q_pt.liq
    state.moisture.ρq_ice = state.ρ * nt.q_pt.ice
end
# Assign prog.ρ = aux.ref_state.ρ in anelastic1D
assign!(state, bl, ::Compressible, nt, ::Prognostic, aux) = nothing
function assign!(state, bl, ::Anelastic1D, nt, ::Prognostic, aux)
    state.ρ = aux.ref_state.ρ
end

function assign!(state, bl, nt, st::Primitive, aux)
    @unpack u, v, w, ρ, p = nt
    state.ρ = ρ
    state.u = SVector(u, v, w)
    state.p = p
    assign!(state, bl, compressibility_model(bl), nt, st, aux)
    assign!(state, bl, moisture_model(bl), nt, st)
end
assign!(state, bl, moisture::DryModel, nt, ::Primitive) = nothing
assign!(state, bl, moisture::EquilMoist, nt, ::Primitive) =
    (state.moisture.q_tot = nt.q_pt.tot)
function assign!(state, bl, moisture::NonEquilMoist, nt, ::Primitive)
    state.moisture.q_tot = nt.q_pt.tot
    state.moisture.q_liq = nt.q_pt.liq
    state.moisture.q_ice = nt.q_pt.ice
end
# Assign prim.p = aux.ref_state.p in anelastic1D
assign!(state, bl, ::Compressible, nt, ::Primitive, aux) = nothing
function assign!(state, bl, ::Anelastic1D, nt, ::Primitive, aux)
    state.p = aux.ref_state.p
end

@testset "Prognostic-Primitive conversion (dry)" begin
    FT = Float64
    compressibility = (Anelastic1D(), Compressible())
    for comp in compressibility
        physics = AtmosPhysics{FT}(
            param_set;
            moisture = DryModel(),
            compressibility = comp,
        )
        bl = AtmosModel{FT}(
            AtmosLESConfigType,
            physics;
            init_state_prognostic = x -> x,
        )
        vs_prog = vars_state(bl, Prognostic(), FT)
        vs_prim = vars_state(bl, Primitive(), FT)
        vs_aux = vars_state(bl, Auxiliary(), FT)
        prog_arr = zeros(varsize(vs_prog))
        prim_arr = zeros(varsize(vs_prim))
        aux_arr = zeros(varsize(vs_aux))
        prog = Vars{vs_prog}(prog_arr)
        prim = Vars{vs_prim}(prim_arr)
        aux = Vars{vs_aux}(aux_arr)
        for nt in PhaseDryProfiles(param_set, ArrayType)
            assign!(aux, bl, nt, Auxiliary())

            # Test prognostic_to_primitive! identity
            assign!(prog, bl, nt, Prognostic(), aux)
            prog_0 = deepcopy(parent(prog))
            prim_arr .= 0
            prognostic_to_primitive!(bl, prim, prog, aux)
            @test !all(parent(prim) .≈ parent(prog)) # ensure not calling fallback
            primitive_to_prognostic!(bl, prog, prim, aux)
            @test all(parent(prog) .≈ prog_0)

            # Test primitive_to_prognostic! identity
            assign!(prim, bl, nt, Primitive(), aux)
            prim_0 = deepcopy(parent(prim))
            prog_arr .= 0
            primitive_to_prognostic!(bl, prog, prim, aux)
            @test !all(parent(prim) .≈ parent(prog)) # ensure not calling fallback
            prognostic_to_primitive!(bl, prim, prog, aux)
            @test all(parent(prim) .≈ prim_0)
        end
    end
end

@testset "Prognostic-Primitive conversion (EquilMoist)" begin
    FT = Float64
    compressibility = (Compressible(),) # Anelastic1D() does not converge
    for comp in compressibility
        physics = AtmosPhysics{FT}(
            param_set;
            moisture = EquilMoist(; maxiter = 5), # maxiter=3 does not converge
            compressibility = comp,
        )
        bl = AtmosModel{FT}(
            AtmosLESConfigType,
            physics;
            init_state_prognostic = x -> x,
        )
        vs_prog = vars_state(bl, Prognostic(), FT)
        vs_prim = vars_state(bl, Primitive(), FT)
        vs_aux = vars_state(bl, Auxiliary(), FT)
        prog_arr = zeros(varsize(vs_prog))
        prim_arr = zeros(varsize(vs_prim))
        aux_arr = zeros(varsize(vs_aux))
        prog = Vars{vs_prog}(prog_arr)
        prim = Vars{vs_prim}(prim_arr)
        aux = Vars{vs_aux}(aux_arr)
        err_max_fwd = 0
        err_max_bwd = 0
        for nt in PhaseEquilProfiles(param_set, ArrayType)
            assign!(aux, bl, nt, Auxiliary())

            # Test prognostic_to_primitive! identity
            assign!(prog, bl, nt, Prognostic(), aux)
            prog_0 = deepcopy(parent(prog))
            prim_arr .= 0
            prognostic_to_primitive!(bl, prim, prog, aux)
            @test !all(parent(prim) .≈ parent(prog)) # ensure not calling fallback
            primitive_to_prognostic!(bl, prog, prim, aux)
            @test all(parent(prog)[1:4] .≈ prog_0[1:4])
            @test isapprox(parent(prog)[5], prog_0[5]; atol = atol_energy)
            # @test all(parent(prog)[5] .≈ prog_0[5]) # fails
            @test all(parent(prog)[6] .≈ prog_0[6])
            err_max_fwd = max(abs(parent(prog)[5] .- prog_0[5]), err_max_fwd)

            # Test primitive_to_prognostic! identity
            assign!(prim, bl, nt, Primitive(), aux)
            prim_0 = deepcopy(parent(prim))
            prog_arr .= 0
            primitive_to_prognostic!(bl, prog, prim, aux)
            @test !all(parent(prim) .≈ parent(prog)) # ensure not calling fallback
            prognostic_to_primitive!(bl, prim, prog, aux)
            @test all(parent(prim)[1:4] .≈ prim_0[1:4])
            # @test all(parent(prim)[5] .≈ prim_0[5]) # fails
            @test isapprox(parent(prim)[5], prim_0[5]; atol = atol_energy)
            @test all(parent(prim)[6] .≈ prim_0[6])
            err_max_bwd = max(abs(parent(prim)[5] .- prim_0[5]), err_max_bwd)
        end
    end
    # We may want/need to improve this later, so leaving debug info:
    # @show err_max_fwd
    # @show err_max_bwd
end

@testset "Prognostic-Primitive conversion (NonEquilMoist)" begin
    FT = Float64
    compressibility = (Anelastic1D(), Compressible())
    for comp in compressibility
        physics = AtmosPhysics{FT}(
            param_set;
            moisture = NonEquilMoist(),
            compressibility = comp,
        )
        bl = AtmosModel{FT}(
            AtmosLESConfigType,
            physics;
            init_state_prognostic = x -> x,
        )
        vs_prog = vars_state(bl, Prognostic(), FT)
        vs_prim = vars_state(bl, Primitive(), FT)
        vs_aux = vars_state(bl, Auxiliary(), FT)
        prog_arr = zeros(varsize(vs_prog))
        prim_arr = zeros(varsize(vs_prim))
        aux_arr = zeros(varsize(vs_aux))
        prog = Vars{vs_prog}(prog_arr)
        prim = Vars{vs_prim}(prim_arr)
        aux = Vars{vs_aux}(aux_arr)
        for nt in PhaseEquilProfiles(param_set, ArrayType)
            assign!(aux, bl, nt, Auxiliary())

            # Test prognostic_to_primitive! identity
            assign!(prog, bl, nt, Prognostic(), aux)
            prog_0 = deepcopy(parent(prog))
            prim_arr .= 0
            prognostic_to_primitive!(bl, prim, prog, aux)
            @test !all(parent(prim) .≈ parent(prog)) # ensure not calling fallback
            primitive_to_prognostic!(bl, prog, prim, aux)
            @test all(parent(prog) .≈ prog_0)

            # Test primitive_to_prognostic! identity
            assign!(prim, bl, nt, Primitive(), aux)
            prim_0 = deepcopy(parent(prim))
            prog_arr .= 0
            primitive_to_prognostic!(bl, prog, prim, aux)
            @test !all(parent(prim) .≈ parent(prog)) # ensure not calling fallback
            prognostic_to_primitive!(bl, prim, prog, aux)
            @test all(parent(prim) .≈ prim_0)
        end
    end
end

@testset "Prognostic-Primitive conversion (array interface)" begin
    FT = Float64
    physics = AtmosPhysics{FT}(param_set; moisture = DryModel())
    bl = AtmosModel{FT}(
        AtmosLESConfigType,
        physics;
        init_state_prognostic = x -> x,
    )
    vs_prog = vars_state(bl, Prognostic(), FT)
    vs_prim = vars_state(bl, Primitive(), FT)
    vs_aux = vars_state(bl, Auxiliary(), FT)
    prog_arr = zeros(varsize(vs_prog))
    prim_arr = zeros(varsize(vs_prim))
    aux_arr = zeros(varsize(vs_aux))
    prog = Vars{vs_prog}(prog_arr)
    prim = Vars{vs_prim}(prim_arr)
    aux = Vars{vs_aux}(aux_arr)
    for nt in PhaseDryProfiles(param_set, ArrayType)
        assign!(aux, bl, nt, Auxiliary())

        # Test prognostic_to_primitive! identity
        assign!(prog, bl, nt, Prognostic(), aux)
        prog_0 = deepcopy(parent(prog))
        prim_arr .= 0
        BL.prognostic_to_primitive!(bl, prim_arr, prog_arr, aux_arr)
        @test !all(prog_arr .≈ prim_arr) # ensure not calling fallback
        BL.primitive_to_prognostic!(bl, prog_arr, prim_arr, aux_arr)
        @test all(parent(prog) .≈ prog_0)

        # Test primitive_to_prognostic! identity
        assign!(prim, bl, nt, Primitive(), aux)
        prim_0 = deepcopy(parent(prim))
        prog_arr .= 0
        BL.primitive_to_prognostic!(bl, prog_arr, prim_arr, aux_arr)
        @test !all(prog_arr .≈ prim_arr) # ensure not calling fallback
        BL.prognostic_to_primitive!(bl, prim_arr, prog_arr, aux_arr)
        @test all(parent(prim) .≈ prim_0)
    end
end

end


text:
module MyExamplePackage

greet() = print("Hello World!")

include("ExtraFile.jl")
export my_f
end # module


text:
push!(LOAD_PATH, "/home/debian/juliapackages")
using LabConnections.BeagleBone
BeagleBone.precompile_bb()
server = run_server()
println("Opening server")
while isopen(server)
    #To keep julia from returning, thus killing the process
    println("server kept open")
    sleep(10)
end


text:
# PREAMBLE

# PKG_SETUP

# ## Setup

using RoboDojo

# ## Initial conditions
q1 = nominal_configuration(hopper)
v1 = zeros(hopper.nq)

# ## Time
T = 100
h = 0.01

# ## Simulator
s = Simulator(hopper, T, h=h)

# ## Simulate
simulate!(s, q1, v1)

# ## Visualizer
vis = Visualizer()
render(vis)

# ## Visualize
visualize!(vis, s)


text:
using JLD2, MCMCChains
import KenyaCoVaccines
using StatsPlots, Distributions, Optim

## Group fitted models by urban (Nairobi and Mombasa) and non-urban

non_urban_filenames = ["modelfits/Nakuru_model.jld2",
    "modelfits/Uasin Gishu_model.jld2",
    "modelfits/Embu_model.jld2",
    "modelfits/Kisumu_model.jld2",
    "modelfits/Nyeri_model.jld2",
    "modelfits/Siaya_model.jld2",
    "modelfits/Kisii_model.jld2",
    "modelfits/Kilifi_model.jld2",
    "modelfits/Kwale_model.jld2"]

## Combine posterior distribution draws of detection parameters across groups

non_urban_parameter_chains = []

for filename in non_urban_filenames
    fit_dict = load(filename)
    model = fit_dict[first(keys(fit_dict))]
    push!(non_urban_parameter_chains, model.MCMC_results.chain)
end

## Define a posterior density for the hyperparameters of an underlying Gamma distribution
#The assumption is that the non-cities have specific parameters drawn from a hierarchy (See Gelman et al)

# log-posterior density
function log_post_density_for_detection_parameters(k, θ, post_array::AbstractArray, param_sym::Symbol, prior_k, prior_θ)
    # LPD=0.0
    LPD = logpdf(prior_k, k) + logpdf(prior_θ, θ)
    for chn in post_array
        param_draws = get(chn, param_sym)
        LPD += log(mean([pdf(Gamma(k, θ), p) for p in param_draws[param_sym]]))
    end
    LPD
end

###
### Fit the population distribution of detection parameters
###
obj_p_test(x) = -log_post_density_for_detection_parameters(exp(x[1]), exp(x[2]), non_urban_parameter_chains, :p_test, Exponential(10), Exponential(0.5 / 10))
obj_χ(x) = -log_post_density_for_detection_parameters(exp(x[1]), exp(x[2]), non_urban_parameter_chains, :χ, Exponential(10), Exponential(1.5 / 10))
##

fit1 = optimize(obj_p_test, log.([10.0, 0.1]), BFGS())
fit2 = optimize(obj_χ, log.([10.0, 0.1]), BFGS())

fitted_pop_paramsptest = exp.(fit1.minimizer)
fitted_pop_paramschi = exp.(fit2.minimizer)

## Save fitted priors --
## NB: Prior variance is doubled compared to fit
fitted_priors = (d_ptest = Gamma(fitted_pop_paramsptest[1] / 2, fitted_pop_paramsptest[2] * 2),
    d_χ = Gamma(fitted_pop_paramschi[1] / 2, fitted_pop_paramschi[2] * 2))

@save("data/non_urban_detection_fits.jld2", fitted_priors)

##Plots to examine the inferred priors

plot(fitted_priors.d_ptest, title = "Fitted prior for p_test (doubled variance)", lab = "")
plot(fitted_priors.d_χ, title = "Fitted prior for chi (doubled variance)", lab = "")


text:
function mdl_source!(voltage_r, voltage_i, current_r, current_i, device::PSY.Source)
    #Load device parameters
    V_R =
        device.ext[CONTROL_REFS][V_source_index] *
        cos(device.ext[CONTROL_REFS][θ_source_index])
    V_I =
        device.ext[CONTROL_REFS][V_source_index] *
        sin(device.ext[CONTROL_REFS][θ_source_index])
    R_th = PSY.get_R_th(device)
    X_th = PSY.get_X_th(device)
    Zmag = R_th^2 + X_th^2

    #update current
    current_r[1] += R_th * (V_R - voltage_r[1]) / Zmag + X_th * (V_I - voltage_i[1]) / Zmag #in system pu flowing out
    current_i[1] += R_th * (V_I - voltage_i[1]) / Zmag - X_th * (V_R - voltage_r[1]) / Zmag #in system pu flowing out

    return
end


text:
# Donut chart w/ VegaLite

using VegaLite, DataFrames

df = DataFrame(category=1:6, value=[4,6,10,3,7,8])

df |> @vlplot(mark={:arc, innerRadius=50}, theta=:value, color="category:n", view={stroke=nothing})

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
export argmin_x

"""
solve the x-minimization step in PARSDMM, i.e., x-minimization w.r.t. the augmented Lagrangian
"""
function argmin_x(
                  Q               ::Union{Array{TF,2},SparseMatrixCSC{TF,Int32},SparseMatrixCSC{TF,Int64}},
                  rhs             ::Vector{TF},
                  x               ::Vector{TF},
                  x_min_solver    ::String,
                  x_solve_tol_ref ::TF,
                  i               ::Integer,
                  log_PARSDMM,
                  Q_offsets=[],
                  Ax_out=zeros(TF,length(x)) ::Vector{TF}
                  ) where {TF<:Real}

    #Initialize
    flag    = 0
    reslres = 0.0
    iter    = 0

  if x_min_solver == "CG_normal" # for package: KrylovMethods
    if typeof(Q)==Array{TF,2}
      #set up multi-threaded matrix-vector product in compressed diagonal storage format
      Af(in) =  (fill!(Ax_out,0); CDS_MVp_MT(size(Q,1),size(Q,2),Q,Q_offsets,in,Ax_out); return Ax_out)

      #determine what relative residual CG needs to reach
      if i<3 #i is the PARSDMM iteration counter
        x_solve_tol_ref = TF(max(0.1*norm(Af(x)-rhs)/norm(rhs),10*eps(TF))) #10% of current relative residual
      else
        x_solve_tol_ref = TF(min(max(0.1*norm(Af(x)-rhs)/norm(rhs),10*eps(TF)),x_solve_tol_ref)) #10% of current relative residual, but not larger than previous relative residual
      end

      (x,flag,relres,iter) = cg(Af,rhs,tol=x_solve_tol_ref,maxIter=1000,x=x,out=0);
    else
      #CG with native Julia sparse CSC format MPVs
      (x,flag,relres,iter) = cg(Q,rhs,tol=x_solve_tol_ref,maxIter=1000,x=x,out=0);
    end

    end

    x_solve_tol_ref=TF(x_solve_tol_ref)
  return x,iter,relres,x_solve_tol_ref
end


text:
# ps3000setadvtriggerchannelproperties.jl

export ps3000setadvtriggerchannelproperties, TriggerChannelProperties, PS3000ThresholdMode

immutable TriggerChannelProperties
  thresholdmajor  :: Int16  # scaled in 16-bit ADC counts
  thresholdminor  :: Int16  # scaled in 16-bit ADC
  hysteresis      :: UInt16 # scaled in 16-bit ADC
  channel         :: Int16 # A = 0, B = 1, C = 2, D = 3 I assume.
  thresholdmode   :: Int16
end

baremodule PS3000ThresholdMode
  const level = 0x0000
  const window = 0x0001
end

function ps3000setadvtriggerchannelproperties(handle :: Int16,
                                              channelproperties :: TriggerChannelProperties,
                                              autotriggermilliseconds)
  nchannelproperties = 1
  ps_status = ccall((:ps3000SetAdvTriggerChannelProperties, ps3000driver), Int16,
    (Int16, Ref{TriggerChannelProperties}, Int16, Int32),
    handle, Ref{TriggerChannelProperties}(channelproperties),
     nchannelproperties, autotriggermilliseconds)
  if ps_status == 0
    error("ps3000 Error: one of the parameters is out of range")
  end
  return nothing
end

function ps3000setadvtriggerchannelproperties(handle :: Int16)
  # call with no conditions to turn trigger off
  nchannelproperties = 0
  autotriggermilliseconds = 0
  ps_status = ccall((:ps3000SetAdvTriggerChannelProperties, ps3000driver), Int16,
    (Int16, Ref{Void}, Int16, Int32),
    handle, Ref{Void}(C_NULL),
     nchannelproperties, autotriggermilliseconds)
  if ps_status == 0
    error("ps3000 Error: one of the parameters is out of range")
  end
  return nothing
end


text:
#
# Copyright 2021, Andrew Lin.
#
# This software is licensed under the MIT License.
# See LICENSE.txt at the root of the project or
# https://opensource.org/licenses/MIT
#
module GenerateCatan
import Random
export generate_board, print_board

# Shuffle possible values, inserting a 0 at the desert location in `resources`.
function generate_values(resources)
    values = Random.shuffle(
        [
            "2",
            "3",  "3",
            "4",  "4",
            "5",  "5",
            "6",  "6",
            # 7 belongs to the robber.
            "8",  "8",
            "9",  "9",
            "10", "10",
            "11", "11",
            "12"
        ]
    )
    idx = findfirst(isequal("desert"), resources)
    return [values[1:idx - 1]; "0"; values[idx:end]]
end

# Generate a randomized Catan board.
function generate_board()
    function water()
        return "water"
    end

    function harbor(idx)
        return "harbor(" * harbors[idx] * ")"
    end

    function resource(idx)
        return resources[idx] * "(" * values[idx] * ")"
    end

    harbors = Random.shuffle(
        [
            "3:1", "3:1", "3:1", "3:1",
            "brick",
            "ore",
            "sheep",
            "wheat",
            "wood"
        ]
    )

    resources = Random.shuffle(
        [
            "brick", "brick", "brick",
            "desert",
            "ore",   "ore",   "ore",
            "sheep", "sheep", "sheep", "sheep",
            "wheat", "wheat", "wheat", "wheat",
            "wood",  "wood",  "wood",  "wood"
        ]
    )

    values = generate_values(resources)

    return [
        water(),   harbor(1),    water(),      harbor(2),
        harbor(3), resource(1),  resource(2),  resource(3),  water(),
        water(),   resource(4),  resource(5),  resource(6),  resource(7),  harbor(4),
        harbor(5), resource(8),  resource(9),  resource(10), resource(11), resource(12), water(),
        water(),   resource(13), resource(14), resource(15), resource(16), harbor(6),
        harbor(7), resource(17), resource(18), resource(19), water(),
        water(),   harbor(8),    water(),      harbor(9)
    ]
end

# Print the generated Catan board.
function print_board(board)
    row_length = [1, 4, 5, 6, 7, 6, 5, 4]
    for row in 1:7
        idx_start = sum(row_length[1:row])
        idx_end = idx_start + row_length[row + 1] - 1
        print("row " * string(row) * ": | ")
        for idx in idx_start:idx_end
            print(board[idx] * " | ")
        end
        println()
    end
end

end

max_stars_repo_name: AndrewWasHere/setup-catan

max_issues_repo_name: AndrewWasHere/setup-catan

max_forks_repo_name: AndrewWasHere/setup-catan


text:
@testset "Point" begin
    @testset "Point{1}" begin
        for T ∈ Floats
            P₁ = Point{1,T}(1)
            P₂ = Point{1,T}(2)
            𝘃₂ = Vec{1,T}(2)

            # coordinates
            coords = coordinates(P₁)
            @test coords == [1]

            # subtraction
            P = P₁ - P₂
            @test P == [-1]

            # vector addition
            v = P₁ + 𝘃₂
            @test v == [3]

            # vector subtraction
            v = P₁ - 𝘃₂
            @test v == [-1]

            # ≈
            @test Point{1,T}(2 + 0.1*ϵ_Point) ≈ Point{1,T}(2)
            @test Point{1,T}(2 + 10*ϵ_Point) ≉ Point{1,T}(2)

            P₁ = Point{1,T}(-1)
            P₂ = Point{1,T}(4)
            # distance
            @test distance(P₁, P₂) ≈ 5

            # distance²
            @test distance²(P₁, P₂) ≈ 25

            # midPoint
            mp = midpoint(P₁, P₂)
            @test mp ≈ [3//2]
        end
    end

    @testset "Point{2}" begin
        for T ∈ Floats
            P₁ = Point{2,T}(1, 2)
            P₂ = Point{2,T}(2, 4)
            𝘃₂ = Vec{2,T}(2,4)

            # coordinates
            coords = coordinates(P₁)
            @test coords == [1, 2]

            # subtraction
            P = P₁ - P₂
            @test P == [-1, -2]

            # vector addition
            v = P₁ + 𝘃₂
            @test v == [3, 6]

            # vector subtraction
            v = P₁ - 𝘃₂
            @test v == [-1, -2]

            # ≈
            @test Point{2,T}(1, 2 + 0.1*ϵ_Point) ≈ Point{2,T}(1,2)
            @test Point{2,T}(1, 2 + 10*ϵ_Point) ≉ Point{2,T}(1,2)

            P₁ = Point{2,T}(1, 2)
            P₂ = Point{2,T}(2, 4)

            # distance
            @test distance(P₁, P₂) ≈ sqrt(5)

            # distance²
            @test distance²(P₁, P₂) ≈ 5

            # midpoint
            mp = midpoint(P₁, P₂)
            @test mp ≈ [3//2, 3]

            # isCCW
            @test  isCCW(Point{2,T}(0,0), Point{2,T}(1,0), Point{2,T}(1,  1))
            @test !isCCW(Point{2,T}(0,0), Point{2,T}(1,0), Point{2,T}(1, -1))
        end
    end

    @testset "Point{3}" begin
        for T ∈ Floats
            P₁ = Point{3,T}(1, 1, 0)
            P₂ = Point{3,T}(1, 0, 1)
            𝘃₂ = Vec{3,T}(1, 0, 1)

            # coordinates
            coords = coordinates(P₁)
            @test coords == [1, 1, 0]

            # subtraction
            P = P₁ - P₂
            @test P == [0, 1, -1]

            # vector addition
            v = P₁ + 𝘃₂
            @test v == [2, 1, 1]

            # vector subtraction
            v = P₁ - 𝘃₂
            @test v == [0, 1, -1]

            # ≈
            @test Point{3,T}(1, 1, 2 + 0.1*ϵ_Point) ≈ Point{3,T}(1,1,2)
            @test Point{3,T}(1, 1, 2 + 10*ϵ_Point) ≉ Point{3,T}(1,1,2)

            P₁ = Point{3,T}(1, 2, 1)
            P₂ = Point{3,T}(2, 4, 0)

            # distance
            @test distance(P₁, P₂) ≈ sqrt(6)

            # distance²
            @test distance²(P₁, P₂) ≈ 6

            # midpoint
            mp = midpoint(P₁, P₂)
            @test mp ≈ [3//2, 3, 1//2]
        end
    end
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
#!/usr/bin/env julia

Fract_Maxwell = RheoModelClass(
          # Model name
          name="fractmaxwell",
          # Model parameters,
          p = [:cₐ, :a, :cᵦ, :β],
          # Relaxation modulus
          G = quote
                 cᵦ*t^(-β)*mittleff(a - β, 1 - β, -cᵦ*t^(a - β)/cₐ)
              end,
          # Creep modulus
          J = quote
                t^(a)/(cₐ*gamma(1 + a)) + t^(β)/(cᵦ*gamma(1 + β))
              end,
          # Storage modulus
          Gp = quote
                 denominator = (cₐ*ω^a)^2 + (cᵦ*ω^β)^2 + 2*(cₐ*ω^a)*(cᵦ*ω^β)*cos((a-β)*π/2)
                 numerator = ((cᵦ*ω^β)^2)*(cₐ*ω^a)*cos(a*π/2) + ((cₐ*ω^a)^2)*(cᵦ*ω^β)*cos(β*π/2)
                 numerator/denominator
               end,
         # Loss modulus
          Gpp = quote
                  denominator = (cₐ*ω^a)^2 + (cᵦ*ω^β)^2 + 2*(cₐ*ω^a)*(cᵦ*ω^β)*cos((a-β)*π/2)
                  numerator = ((cᵦ*ω^β)^2)*(cₐ*ω^a)*sin(a*π/2) + ((cₐ*ω^a)^2)*(cᵦ*ω^β)*sin(β*π/2)
                  numerator/denominator
                end,
          # Constraints
          constraint = quote
                   all([ (a<1) & (a>0)
                         (β<1) & (β>0)
                         -a+β < 0     ])
                  end,
          # Network
          info= "
             ___╱╲__________╱╲____
                ╲╱          ╲╱
                  cₐ,a         cᵦ, β
                 "
          )

FractS_Maxwell = RheoModelClass(
        # Model name
        name="fractmaxwell_spring",
        # Model parameters,
        p = [:cₐ, :a, :k],
        # Relaxation modulus
        G = quote
                # special case when α==0.0, 2 springs in series
                a==0.0 ? 1.0/(1.0/cₐ + 1.0/k) : k*mittleff(a, -k*t^a/cₐ) # normal case otherwise
            end,
        # Creep modulus
        J = quote
              t^(a)/(cₐ*gamma(1 + a)) + 1/k
            end,
        # Storage modulus
        Gp = quote
                denominator = (cₐ*ω^a)^2 + k^2 + 2*(cₐ*ω^a)*k*cos(a*π/2)
                numerator = k^2*(cₐ*ω^a)*cos(a*π/2) + (cₐ*ω^a)^2*k
                numerator/denominator
             end,
       # Loss modulus
        Gpp = quote
                denominator = (cₐ*ω^a)^2 + k^2 + 2*(cₐ*ω^a)*k*cos(a*π/2)
                numerator = k^2*(cₐ*ω^a)*sin(a*π/2)
                numerator/denominator
              end,
        # Constraints
        constraint = quote
                 (a<1) & (a>0)
                end,
        # Network
        info= "
           ___╱╲_________╱╲  ╱╲  ╱╲  ________
              ╲╱           ╲╱  ╲╱  ╲╱
                cₐ,a               k
               "
        )

FractD_Maxwell = RheoModelClass(
          # Model name
          name="fractmaxwell_dashpot",
          # Model parameters,
          p = [:η, :cᵦ, :β],
          # Relaxation modulus
          G = quote
                cᵦ*t^(-β)*mittleff(1 - β, 1 - β, -cᵦ*t^(1 - β)/η)
              end,
          # Creep modulus
          J = quote
                t/η + t^β/(cᵦ*gamma(1 + β))
              end,
          # Storage modulus
          Gp = quote
                  denominator = (η*ω)^2 + (cᵦ*ω^β)^2 + 2*(η*ω)*(cᵦ*ω^β)*cos((1-β)*π/2)
                  numerator = ((η*ω)^2)*(cᵦ*ω^β)*cos(β*π/2)
                  numerator/denominator
               end,
         # Loss modulus
          Gpp = quote
                  denominator = (η*ω)^2 + (cᵦ*ω^β)^2 + 2*(η*ω)*(cᵦ*ω^β)*cos((1-β)*π/2)
                  numerator = ((cᵦ*ω^β)^2)*(η*ω) + ((η*ω)^2)*(cᵦ*ω^β)*sin(β*π/2)
                  numerator/denominator
                end,
          # Constraints
          constraint = quote
                   (β<1) & (β>0)
                  end,
          # Network
          info= "
                  ___
              _____| |_________╱╲____
                  _|_|         ╲╱
                    η            cᵦ, β
                 "
          )

Maxwell = RheoModelClass(
        # Model name
        name="maxwell",
        # Model parameters,
        p = [:η, :k],
        # Relaxation modulus
        G = quote
              k*exp(-k*t/η)
            end,
        # Creep modulus
        J = quote
              t/η + 1/k
            end,
        # Storage modulus
        Gp = quote
                denominator = 1 + η^2*ω^2/k^2
                numerator = η^2*ω^2/k
                numerator/denominator
             end,
       # Loss modulus
        Gpp = quote
                denominator = 1 + η^2*ω^2/k^2
                numerator = η*ω
                numerator/denominator
              end,
        # Network
        info= "
                ___
            _____| |________╱╲  ╱╲  ╱╲  ___
                _|_|          ╲╱  ╲╱  ╲╱
                  η                  k
               "
        )


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 13c9b84e-fdb9-11ea-2e2f-bd38a244317f
using Pkg, DrWatson

# ╔═╡ 13ca0af6-fdb9-11ea-2956-45f2849a856c
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanSample
	using StatisticalRethinking
end

# ╔═╡ 13d69b68-fdb9-11ea-3eaf-612320346d53
for i in 5:7
  include(projectdir("models", "05", "m5.$(i)s.jl"))
end

# ╔═╡ e2ddee76-fdb8-11ea-1f1a-09de92c87dd3
md"## Clip-05-40.2s.jl"

# ╔═╡ 13ca9958-fdb9-11ea-000f-2d2d5d211e5b
md"### snippet 5.39"

# ╔═╡ 3c99c6ec-8347-11eb-2916-8d9f3651ad37
md"
!!! note
	Restart this notebook to re-execute below `include()`.
"

# ╔═╡ 13d72238-fdb9-11ea-364d-8db015412817
if success(rc5_5s)

  post5_5s_df = read_samples(m5_5s, :dataframe)
  title5 = "Kcal_per_g vs. neocortex_perc" * "\n89% predicted and mean range"
  fig1 = plotbounds(
    df, :neocortex_perc, :kcal_per_g,
    post5_5s_df, [:a, :bN, :sigma];
    title=title5,
    rescale_axis=false
  )

  post5_6s_df = read_samples(m5_6s, :dataframe)
  title6 = "Kcal_per_g vs. log mass" * "\nshowing 89% predicted and mean range"
  fig2 = plotbounds(
    df, :lmass, :kcal_per_g,
    post5_6s_df, [:a, :bM, :sigma];
    title=title6,
    rescale_axis=false
  )

  post5_7s_df = read_samples(m5_7s, :dataframe)
  title7 = "Counterfactual,\nholding M=0.0"
  fig3 = plotbounds(
    df, :neocortex_perc, :kcal_per_g,
    post5_7s_df, [:a, :bN, :sigma];
    title=title7,
    rescale_axis=false
  )

  title8 = "Counterfactual,\nholding N=0.0"
  fig4 = plotbounds(
    df, :lmass, :kcal_per_g,
    post5_7s_df, [:a, :bM, :sigma];
    title=title8,
    xlab="log(mass)",
    rescale_axis=false
  )

end;

# ╔═╡ 13e216fa-fdb9-11ea-38d0-83247639f99e
plot(fig1, fig2, fig3, fig4, layout=(2, 2))

# ╔═╡ 13e2a9aa-fdb9-11ea-2bce-5515b878fbed
md"## End of clip-05-40.2s.jl"

# ╔═╡ Cell order:
# ╟─e2ddee76-fdb8-11ea-1f1a-09de92c87dd3
# ╠═13c9b84e-fdb9-11ea-2e2f-bd38a244317f
# ╠═13ca0af6-fdb9-11ea-2956-45f2849a856c
# ╟─13ca9958-fdb9-11ea-000f-2d2d5d211e5b
# ╟─3c99c6ec-8347-11eb-2916-8d9f3651ad37
# ╠═13d69b68-fdb9-11ea-3eaf-612320346d53
# ╠═13d72238-fdb9-11ea-364d-8db015412817
# ╠═13e216fa-fdb9-11ea-38d0-83247639f99e
# ╟─13e2a9aa-fdb9-11ea-2bce-5515b878fbed


text:
abstract type AbstractLayer end
abstract type AbstractChain end

include("flux.jl")
include("affine.jl")
include("resnet.jl")
include("convex.jl")


text:
using Documenter, F1Method
using LinearAlgebra, DiffEqBase, ForwardDiff

makedocs(
    sitename="F1Method Documentation",
    # options
    modules = [F1Method]
)

deploydocs(
    repo = "github.com/briochemc/F1Method.jl.git",
)


text:
import NNlib: conv_direct, conv_im2col, channels_in, channels_out

@testset "Conv Inference" begin
    for T in (Float32, Float64)
        impl = [conv, conv_direct, conv_im2col]
        if NNlib.is_nnpack_available() && T == Float32
            push!(impl, NNlib.conv_nnpack)
        end

        x = rand(T, 10, 10, 3, 2)
        w = rand(T, 3, 3, 3, 1)
        cdims = DenseConvDims(x, w)
        dy = conv(x, w, cdims)

        for f in impl
            @test @inferred(f(x, w, cdims)) isa Array{T,4}
        end

        @test @inferred(conv(x, w)) isa Array{T,4}
        @test @inferred(∇conv_filter(x, dy, cdims)) isa Array{T,4}
        @test @inferred(∇conv_data(dy, w, cdims)) isa Array{T,4}
    end
end

@testset "DepthwiseConv Inference" begin
    for T in (Float32, Float64)
        x = rand(T, 10, 10, 3, 2)
        w = rand(T, 3, 3, 3, 3)
        cdims = DepthwiseConvDims(x, w)
        dy = depthwiseconv(x, w)

        @test @inferred(depthwiseconv(x, w)) isa Array{T,4}
        @test @inferred(∇depthwiseconv_filter(x, dy, cdims)) isa Array{T,4}
        @test @inferred(∇depthwiseconv_data(dy, w, cdims)) isa Array{T,4}
    end
end

@testset "DenseConvDims Inference" begin
    # this needs to be in a function to trigger inference problems
    function channels_in_test(w::AbstractArray)
        cdims = DenseConvDims((1,1,1,1), size(w))
        channels_in(cdims)
    end

    # this needs to be in a function to trigger inference problems
    function channels_out_test(w::AbstractArray)
        cdims = DenseConvDims((1,1,1,1), size(w))
        channels_out(cdims)
    end

    w = rand(Float32, 1, 1, 1, 1)
    @test @inferred(channels_in_test(w)) isa Int
    @test @inferred(channels_out_test(w)) isa Int
end

@testset "Pooling inference" begin
    for T in (Float32, Float64)
        x = rand(T, 10, 10, 3, 2)
        pdims = PoolDims(x, 3)

        y_maxpool = NNlib.maxpool(x, pdims)
        y_meanpool = NNlib.meanpool(x, pdims)
        dy = ones(T, size(y_maxpool)...)

        @test @inferred(NNlib.maxpool(x, pdims)) isa Array{T, 4}
        @test @inferred(NNlib.meanpool(x, pdims)) isa Array{T, 4}
        @test @inferred(NNlib.∇maxpool(dy, y_maxpool, x, pdims)) isa Array{T, 4}
        @test @inferred(NNlib.∇maxpool(dy, y_meanpool, x, pdims)) isa Array{T, 4}
    end
end


text:
## Initial setup
# using Revise
using Juno
using IncrementalInference
using RoME
using DistributedFactorGraphs
using Test

dfg = initfg()
addVariable!(dfg, :x1, Pose2)
addVariable!(dfg, :x2, Pose2)
addVariable!(dfg, :x3, Pose2)
addVariable!(dfg, :l1, Pose2)
addVariable!(dfg, :l2, Pose2)

prior = PriorPose2( MvNormal([10; 10; pi/6.0], Matrix(Diagonal([0.1;0.1;0.05].^2))))
addFactor!(dfg, [:x1], prior )

pp = Pose2Pose2(MvNormal([10.0;0;pi/3], Matrix(Diagonal([0.1;0.1;0.1].^2))))
p2br = Pose2Point2BearingRange(Normal(0,0.1),Normal(20.0,1.0))

# FYI deep copies not required, since no per factor specific data in pp
addFactor!(dfg, [:x1, :x2], pp, autoinit=false)
addFactor!(dfg, [:x2, :x3], pp, autoinit=false)
addFactor!(dfg, [:x1, :l1], deepcopy(p2br), autoinit=false)
# TODO: Ask why it can't resolve here.
addFactor!(dfg, [:x2, :l1], deepcopy(p2br), autoinit=false)
addFactor!(dfg, [:x3, :l2], deepcopy(p2br), autoinit=false)

# Show it
DFG.toDotFile(dfg, "/tmp/testRmMarg.dot")

# Test with a copy
dfgPrime = deepcopy(dfg)
elimOrder = getEliminationOrder(dfgPrime)
# Forcing elimination order to match the documentation
elimOrder = [:l1, :l2, :x1, :x2, :x3]

# tree for elim order: [:l1, :l2, :x1, :x2, :x3]
# :l1 <-- :x1, :x2
# :x1 <-- :x2
# :x2 <-- :x3
# :l2 <-- :x3
# :x3 <-- /-/

buildBayesNet!(dfgPrime, elimOrder)
# Reproduce example
# IncrementalInference.rmVarFromMarg(dfgPrime, lm[1], )
DFG.toDotFile(dfgPrime, "/tmp/testRmMarg.dot")
# Assert that everything was eliminated and every variable has a BayesNetVertID
@test all(map(v -> solverData(v.dfgNode).eliminated, values(dfgPrime.g.vertices)))
# @test all(map(v -> getData(v).BayesNetVertID != nothing, DFG.getVariables(dfgPrime)))
# Assert that we have the expected Bayes tree
expectedBayesOutVertDict = Dict{Symbol, Vector{Symbol}}(:x2 => [:x3], :l1 => [:x1, :x2], :x3 => [], :l2 => [:x3], :x1 => [:x2])
for (vId,linkIds) in expectedBayesOutVertDict
    v = DFG.getVariable(dfgPrime, vId)
    @test setdiff(solverData(v).BayesNetOutVertIDs, linkIds) == []
end

# Now build the tree
# global tree = emptyBayesTree()
# buildTree!(tree, dfgPrime, p)

## End-to-end
@info "Time to build the Bayes Tree..."
tree = wipeBuildNewTree!(dfg, drawpdf=true, show=false)
drawTree(tree, show=true)

# TODO: Test
# Upward solve steps with clique state machine
# @info "Complete upward solve..."
# smtasks, ch = initInferTreeUp!(dfg, tree, drawtree=true, recordcliqs=[:x1;:x2;:x3] )

## Upward solve steps with clique state machine
# @info "Complete upward solve..."
# smtasks, ch = initInferTreeUp!(fgl, bt, N=N, drawtree=drawpdf, recordcliqs=recordcliqs, limititers=limititers, skipcliqids=skipcliqids )

# @info "solve leaf clique with single state machine"
# alltasks[i] = @async tryCliqStateMachineSolve!(fgl, treel, i, cliqHistories, drawtree=drawtree, N=N, limititers=limititers, recordcliqs=recordcliqs)
# manually

resetTreeCliquesForUpSolve!(tree)
setTreeCliquesMarginalized!(dfg, tree)

cliq = whichCliq(tree, :x1)

# history = cliqInitSolveUpByStateMachine!(dfg, tree, cliq, drawtree=true,
#                                          limititers=50, recordhistory=true )
# 0

children = TreeClique[]
for ch in Graphs.out_neighbors(cliq, tree.bt)
  push!(children, ch)
end
prnt = getParent(tree, cliq)
csmc = CliqStateMachineContainer(dfg, initfg(), tree, cliq, prnt, children, false, true, true, getSolverParams(dfg))
statemachine = StateMachine{CliqStateMachineContainer}(next=testCliqCanRecycled_StateMachine) # isCliqUpSolved_StateMachine)

limititers = 100
recordhistory = true

statemachine(csmc, verbose=true, iterlimit=limititers, recordhistory=recordhistory)
statemachine(csmc, verbose=true, iterlimit=limititers, recordhistory=recordhistory)
statemachine(csmc, verbose=true, iterlimit=limititers, recordhistory=recordhistory)
statemachine(csmc, verbose=true, iterlimit=limititers, recordhistory=recordhistory)
statemachine(csmc, verbose=true, iterlimit=limititers, recordhistory=recordhistory)
statemachine(csmc, verbose=true, iterlimit=limititers, recordhistory=recordhistory)

statemachine(csmc, verbose=true, iterlimit=limititers, recordhistory=recordhistory)

lsf(csmc.cliqSubFg, :x2)

writeGraphPdf(csmc.cliqSubFg, show=true, engine="neato")

statemachine(csmc, verbose=true, iterlimit=limititers, recordhistory=recordhistory)
statemachine(csmc, verbose=true, iterlimit=limititers, recordhistory=recordhistory)

csmc.cliqSubFg.nodeCounter

writeGraphPdf(csmc.cliqSubFg, show=true, engine="neato")

addFactor!(dfg, [:x1, :x2], pp, autoinit=false)
DFG.lsf(dfg, :x1x2f1)

0

statemachine(csmc, verbose=true, iterlimit=limititers, recordhistory=recordhistory)
statemachine(csmc, verbose=true, iterlimit=limititers, recordhistory=recordhistory)
statemachine(csmc, verbose=true, iterlimit=limititers, recordhistory=recordhistory)

DFG.getFactorIds(dfg)

cliqHistories = Dict{Int,Vector{Tuple{DateTime, Int, Function, CliqStateMachineContainer}}}()
isTreeSolved(treel, skipinitialized=true)

# TODO: Nary tests for factors....

DFG.getFactor(dfg, :x1x2f1)

## See picture of upward clique association matrix

using Gadfly
cliq = whichCliq(tree, :x1)
spyCliqMat(cliq)

###


text:
using Dates
using AWSSQS

function run_with_retries(
    f,
    args...;
    retry_time = Dates.Second(32),
    failure_message = "Failed after retrying",
    kwargs...,
)
    # Exponential backoff algorithm recommended by Google Cloud:
    # https://cloud.google.com/iot/docs/how-tos/exponential-backoff
    time_of_failure = nothing
    nretries = 0
    while true
        q = f(args..., kwargs...)
        if isnothing(q)
            if isnothing(time_of_failure)
                time_of_failure = now()
            end
            if now() - time_of_failure > retry_time
                break
            end
            sleep(Dates.Second(2^nretries) + Dates.Millisecond(rand(0:1000)))
            nretries += 1
        else
            return q
        end
    end
    error(failure_message)
end

sqs_get_queue_with_retries(args...; kwargs...) = run_with_retries(
    sqs_get_queue,
    args...;
    failure_message = "Queue for communicating results is nonexistent",
    kwargs...
)

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
using Plots
using Reproduce
using JLD2
using FileIO
using Statistics

include("plot_utils.jl")

function plot_ic_reward!(p, ic, logger_key, log_interval; smooth_step=10, label=nothing)
    data = load_results(ic, logger_key)
    xs, mean_line, std_err_line = get_lines(data, log_interval, smooth_step)
    plot!(p, xs, mean_line, ribbons=std_err_line, label=label)
end

function plot_ic_reward_per_gvf!(p, ic, logger_key, log_interval, gvf; smooth_step=10, label=nothing)
    data = load_results(ic, logger_key)[gvf:gvf, :, :]
    # println(size(data))
    # sdfsd
    xs, mean_line, std_err_line = get_lines(data, log_interval, smooth_step)
    plot!(p, xs, mean_line, ribbons=std_err_line, label=label)
end

function get_lines(data, log_interval, smooth_step)
    means = []
    for start in 1:smooth_step:size(data)[2]
        chunk = data[:, start:min(start + smooth_step - 1, size(data)[2]),: ]
        # println(size(chunk))
        chunk_mean = mean(chunk, dims=2)
        push!(means, chunk_mean)
    end

    smoothed = cat(means..., dims=2)[:,:, :]
    sum_lines = sum(smoothed, dims=1)
    num_runs = size(sum_lines)[3]
    mean_line = dropdims(mean(sum_lines, dims=3), dims=(1, 3))
    std_err_line = dropdims(std(sum_lines, dims=3), dims=(1, 3)) / sqrt(num_runs)

    xs = 1:log_interval*smooth_step:log_interval*smooth_step*length(mean_line)
    return xs, mean_line, std_err_line
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
Zygote.@adjoint function ScatteredOut(output, k = 1)
    function ∇scattered(δ)
        return (δ.output, nothing)
    end
    ScatteredOut(output, k), ∇scattered
end
Zygote.@adjoint function ScatteredOut(m, k, output)
    function ∇scattered(δ)
        return (nothing, nothing, δ.output)
    end
    ScatteredOut{eltype(output),length(output)}(m, k, output), ∇scattered
end

Zygote.@adjoint function getindex(F::T, i::Integer) where {T<:Scattered}
    function getInd_rrule(Ȳ)
        zeroNonRefed = map(ii -> ii - 1 == i ? Ȳ : zeros(eltype(F.output[ii]),
                size(F.output[ii])...),
            (1:length(F.output)...,))
        ∂F = T(F.m, F.k, zeroNonRefed)
        return ∂F, nothing
    end
    return getindex(F, i), getInd_rrule
end

Zygote.@adjoint function getindex(F::T, inds::AbstractArray) where {T<:Scattered}
    function getInd_rrule(Ȳ)
        zeroNonRefed = map(ii -> ii - 1 in inds ? Ȳ[indexin(ii - 1, inds)[1]] :
                                 zeros(eltype(F.output[ii]), size(F.output[ii])...),
            (1:length(F.output)...,))
        ∂F = T(F.m, F.k, zeroNonRefed)
        return ∂F, nothing
    end
    return getindex(F, inds), getInd_rrule
end

Zygote.@adjoint function getindex(x::T, p::pathLocs) where {T<:Scattered}
    function getInd_rrule(Δ)
        zeroNonRefed = map(ii -> zeros(eltype(x.output[ii]),
                size(x.output[ii])...),
            (1:length(x.output)...,))
        ∂x = T(x.m, x.k, zeroNonRefed)
        ∂x[p] = Δ
        return ∂x, nothing
    end
    return getindex(x, p), getInd_rrule
end

function rrule(::typeof(flatten), scatRes)
    function ∇flatten(Δarray)
        return (NO_FIELDS, roll(Δarray, scatRes),)
    end
    return flatten(scatRes), ∇flatten
end
function rrule(::typeof(roll), toRoll, stOutput)
    function ∇roll(Δ)
        return NO_FIELDS, flatten(Δ), NO_FIELDS
    end
    return roll(toRoll, stOutput), ∇roll
end


text:
module EquivalenceClassesConstructor

using DataStructures, Printf, JLD, StaticArrays, Lazy
export Predicate, Mapping, ReduceMap, ExpandMap, EquivalenceClasses, find_classes, minimal_set
export toDirectMap, toIndexMap, labelsMap, minimal_set
export write_fixed_width, write_JLD

  include("./Types.jl")
  include("./GraphTools.jl")
  include("./helpers.jl")
  include("./IO.jl")

end


text:
using SafeTestsets

@safetestset "Create Growth Vector Tests" begin include("create_growth_vector_tests.jl") end
@safetestset "Create Community Matrix Tests" begin include("create_community_matrix_tests.jl") end
@safetestset "GLV Analysis Tests" begin include("GLV_analysis_tests.jl") end
@safetestset "Holes tests" begin include("holes_tests.jl") end

max_issues_repo_name: mtangulo/coexistenceholes

max_forks_repo_name: mtangulo/coexistenceholes


text:
module EasyServers

export
    runserver

using Sockets
using Sockets: TCPServer, TCPSocket
using Base: IOError

"""
    runserver(srv, f)

runs TCP server socket `srv` with function `f` to handle client connections.

The server accepts connections from clients and manages to have `f(sock)`
repeatedly called by an asynchronous task to process requests from any socket
`sock` connected to a client.  The call `f(sock)` shall return a symbolic
result, in particular: `:continue` to continue processing client requests or
`:quit` if the server is asked to quit.  If anything else than `:continue` is
returned, the client connection is closed if not yet done.  If `:quit` is
returned, the server and all client connections are closed and the `runserver`
function returns.  All these connections are also closed in case of exception.

The server socket may be automatically created if `runserver` is called as:

     runserver([addr=ip"127.0.0.1",] port, f; kwds...)

with `addr` and `port` the address and port number to be listened by the server
and `kwds...` the keywords to pass to the `Sockets.listen` method.

"""
function runserver(server::TCPServer, callback::Function)
    # The server is eventually closed to quit and make `tryaccept` returns
    # `nothing`.  A `try/catch` block with a `finally` clause is needed to
    # handle exceptions and make sure that all clients and the server are
    # properly closed on return.
    clients = Set{TCPSocket}()
    try
        while true
            sock = tryaccept(server)
            sock === nothing && break
            @async begin
                push!(clients, sock)
                while isopen(sock)
                    status = callback(sock)
                    if status !== :continue
                        # Close server connection to quit.
                        status === :quit && closeifopen(server)
                        break
                    end
                end
                closeifopen(pop!(clients, sock))
            end
        end
    catch ex
        rethrow(ex)
    finally
        closeifopen(server)
        closeall!(clients)
    end
    nothing
end

runserver(port::Integer, args...; kwds...) =
    runserver(Sockets.localhost, port, args...; kwds...)

runserver(addr, port::Integer, args...; kwds...) =
    runserver(listen(addr, port), args...; kwds...)

"""
    sock = tryaccept(server)

attempts to accept a client connection on the given server.  On success, the
result is a socket connection to the client.  If the server connection has been
closed or is closed while waiting for a client, `nothing` is returned.

"""
function tryaccept(server::TCPServer)
    try
        accept(server)
    catch ex
        isopen(server) && rethrow(ex)
        nothing
    end
end

"""
    closeifopen(x)

closes `x` unless it has been closed.

"""
closeifopen(x) = isopen(x) && close(x)

"""
    closeall!(coll)

pops all items from collection `coll` closing them unless already closed.

"""
function closeall!(collection)
    while !isempty(collection)
        closeifopen(pop!(collection))
    end
    nothing
end

end # module


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test, Libdl, LibUV_jll

@testset "LibUV_jll" begin
    vn = VersionNumber(unsafe_string(ccall((:uv_version_string, libuv), Cstring, ())))
    @test vn == v"2.0.0-dev"
end

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
__precompile__()

module Wkhtmltox

import Base: getindex, setindex!, push!, run

export getindex, setindex!, push!, run

export PdfSettings, PdfObject, PdfConverter
export pdf_init, pdf_deinit
export ImgSettings, ImgConverter
export img_init, img_deinit

const depfile = joinpath(dirname(@__FILE__),"..","deps","deps.jl")
if isfile(depfile)
    include(depfile)
else
    error("Wkhtmltox not properly installed. Please run Pkg.build(\"Wkhtmltox\")")
end

type GlobalSettings ; end
type ObjectSettings ; end
type Converter; end

include("pdf.jl")
include("image.jl")

end # module


text:
test1 = """
    16
    10
    15
    5
    1
    11
    7
    19
    6
    12
    4
    """

test2 = """
    28
    33
    18
    42
    31
    14
    46
    20
    48
    47
    24
    23
    49
    45
    19
    38
    39
    11
    1
    32
    25
    35
    8
    17
    7
    9
    4
    2
    34
    10
    3
    """

function parse_input(input::AbstractString)::Vector{Int}
    return parse.(Int, split(input, r"\s+"; keepempty = false)) |> sort
end

function q1(ratings::AbstractVector{<: Integer})::Int
    ds = diff([0; ratings])
    return count(==(1), ds) * (count(==(3), ds) + 1)
end

function q2(ratings::AbstractVector{<: Integer})::Int
    n_comb = [1; fill(0, last(ratings))]
    for r in ratings
        n_comb[r + 1] = sum(n_comb[max(1, r - 2):r])
    end
    return last(n_comb)
end

let r1 = parse_input(test1),
    r2 = parse_input(test2)

    @assert q1(r1) == 35
    @assert q1(r2) == 220
    @assert q2(r1) == 8
    @assert q2(r2) == 19208
end

ratings = parse_input(read("day10.in", String))
println("Q1: ", q1(ratings))
println("Q2: ", q2(ratings))

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
module IncDataStats

export GroupBy,
       update!,
       minimum,
       maximum,
       mean,
       sum,
       count,
       hist,
       quantile,
       stepgrouper

include("core.jl")

end # module


text:
using Oceananigans.Architectures
using Oceananigans.BoundaryConditions
using Oceananigans.Fields
using Oceananigans.Utils

using Oceananigans.Grids: topology, total_size, halo_size, interior_parent_indices

"""
    Average{F, R, D, P, I, T} <: AbstractDiagnostic

A diagnostic for computing the averages of a field along particular dimensions.
"""
mutable struct Average{F, R, D, P, I, T} <: AbstractDiagnostic
                 field :: F
                  dims :: D
                result :: P
    iteration_interval :: I
         time_interval :: T
              previous :: Float64
           return_type :: R
            with_halos :: Bool
end

function dims_to_result_size(field, dims, grid)
    field_size = total_size(parent(field))
    return Tuple(d in dims ? 1 : field_size[d] for d in 1:3)
end

"""
    Average(field::Field; dims, iteration_interval=nothing, time_interval=nothing, return_type=Array)

Construct an `Average` of `field` along the dimensions specified by the tuple `dims`.

After the average is computed it will be stored in the `result` property.

The `Average` can be used as a callable object that computes and returns the average.

An `iteration_interval` or `time_interval` (or both) can be passed to indicate how often to
run this diagnostic if it is part of `simulation.diagnostics`. `iteration_interval` is a
number of iterations while `time_interval` is a time interval in units of `model.clock.time`.

A `return_type` can be used to specify the type returned when the `Average` is
used as a callable object. The default `return_type=Array` is useful when running a GPU
model and you want to save the output to disk by passing it to an output writer.
"""
function Average(field::Field; dims, iteration_interval=nothing, time_interval=nothing, return_type=Array,
                 with_halos=true)

    dims isa Union{Int, Tuple} || error("Average dims must be an integer or tuple!")
    dims isa Int && (dims = tuple(dims))

    length(dims) == 0 && error("dims is empty! Must average over at least one dimension.")
    length(dims) > 3  && error("Models are 3-dimensional. Cannot average over 4+ dimensions.")
    all(1 <= d <= 3 for d in dims) || error("Dimensions must be one of 1, 2, 3.")

    arch = architecture(field)
    result_size = dims_to_result_size(field, dims, field.grid)
    result = zeros(arch, field.grid, result_size...)

    return Average(field, dims, result, iteration_interval, time_interval, 0.0,
                   return_type, with_halos)
end

"""
    normalize_sum!(avg)

Normalize the sum by the number of grid points averaged over to get the average.
"""
function normalize_sum!(avg)
    N = size(avg.field.grid)
    avg.result ./= prod(N[d] for d in avg.dims)
    return nothing
end

"""
    run_diagnostic(model, avg::Average)

Compute the horizontal average of `avg.field` and store the result in `avg.result`.
"""
function run_diagnostic(model, avg::Average)
    zero_halo_regions!(parent(avg.field), location(avg.field), avg.field.grid)
    sum!(avg.result, parent(avg.field))
    normalize_sum!(avg)
    return nothing
end

function (avg::Average)(model)
    run_diagnostic(model, avg)
    result = avg.return_type(avg.result)

    if avg.with_halos
        return result
    else
        loc = location(avg.field)
        topo = topology(model.grid)
        N, H = size(model.grid), halo_size(model.grid)
        inds = (d in avg.dims ? Colon() : interior_parent_indices(loc[d], topo[d], N[d], H[d]) for d in 1:3)
        return result[inds...]
    end
end


text:
using EMST
using Plots
using GR
function plot_emst_2d(x,edges)
    Plots.gr()
    Plots.plot()
    Plots.scatter!(x[1,:],x[2,:],ms=3.0,marker=(stroke(0,:gray)))
    for zr in 1:size(edges,1)
      Plots.plot!( x[1,edges[zr,:]] , x[2,edges[zr,:]] , linecolor = :gray)
    end
    Plots.scatter!(legend=false) # for some reason with this the plot shows im juno.. :)
end
function plot_emst_3d(x,edges)
    Plots.gr()
    Plots.plot()
    Plots.scatter!(x[1,:],x[2,:],x[3,:],ms=3.0,marker=(stroke(0,:gray)))
    for zr in 1:size(edges,1)
      Plots.plot!( x[1,edges[zr,:]] , x[2,edges[zr,:]] , x[3,edges[zr,:]] , linecolor = :gray)
    end
    Plots.scatter!(legend=false) # for some reason with this the plot shows im juno.. :)
end

x_test  = rand(2,400)
tic(); e_test  = EMST.compute_emst(x_test;nmin=64); toc();
EMST.verify_emst(x_test,e_test,200)
plot_emst_2d(x_test,e_test)

x_test  = rand(3,400)
tic(); e_test  = EMST.compute_emst(x_test;nmin=64); toc();
EMST.verify_emst(x_test,e_test,200)
plot_emst_3d(x_test,e_test)

Plots.savefig("C:\\Temp\\emst_2d.png")


text:
using Core: Vector
export getXYZValues
export setAllRange
export sleepModeOn,sleepModeOff,lockOn,lockOff
export setUnitToTesla,setStandardSettings, getFieldError, calculateFieldError

include("LakeShoreLowLevel.jl")

function LakeShoreGaussMeter(params::Dict)
  gauss = LakeShoreGaussMeter(params["connection"], params["coordinateTransformation"])

  setStandardSettings(gauss)

  setAllAutoRanging(gauss, params["autoRanging"] ? '1' : '0')
  if !params["autoRanging"]
    setAllRange(gauss, string(params["range"])[1])
  end

  setFast(gauss, params["fast"] ? '1' : '0')
  return gauss
end

"""
Returns x,y, and z values and apply a coordinate transformation
"""
function getXYZValues(gauss::LakeShoreGaussMeter)
	field = parse.(Float32,split(getAllFields(gauss),","))[1:3]
	multipliers = getAllMultipliers(gauss)
    gauss.coordinateTransformation*(field.*multipliers)*Unitful.T
end

"""
Sets the range of the active channel. The range depends on the installed probe.
For HSE Probe. More Information in part 3.4 on page 3-7.
	0 = highest = ±3T
	1           = ±300mT
	2           = ±30mT
	3 = lowest  = ±3mT

Sets the range of the X channel
"""
function setXRange(gauss::LakeShoreGaussMeter, range::Char)
	setActiveChannel(gauss, 'X')
	setRange(gauss, range)
	return nothing
end

"""
Sets the range of the Y channel
"""
function setYRange(gauss::LakeShoreGaussMeter, range::Char)
	setActiveChannel(gauss, 'Y')
	setRange(gauss, range)
	return nothing
end

"""
Sets the range of the Z channel
"""
function setZRange(gauss::LakeShoreGaussMeter, range::Char)
	setActiveChannel(gauss, 'Z')
	setRange(gauss, range)
	return nothing
end

"""
Sets the range of all channels
"""
function setAllRange(gauss::LakeShoreGaussMeter, range::Char)
	setXRange(gauss, range)
	setYRange(gauss, range)
	setZRange(gauss, range)
	return nothing
end

function getFieldError(gauss::LakeShoreGaussMeter, range::Int)
    if range == 0
        return 150Unitful.μT
    elseif range == 1
        return 15Unitful.μT
    elseif range == 2
        return 1.5Unitful.μT
    elseif range == 3
        return 0.15Unitful.μT
    end
end

function calculateFieldError(gauss::LakeShoreGaussMeter, magneticField::Vector{<:Unitful.BField})
	magneticFieldError = zeros(typeof(1.0u"T"),3,2)
	magneticFieldError[:,1] = abs.(magneticField)*1e-3
	magneticFieldError[:,2] .= getFieldError(gauss, range)
	return sum(magneticFieldError, dims=2)
end

"""
Sets the sleep mode on
"""
function sleepModeOn(gauss::LakeShoreGaussMeter)
	setSleepMode(gauss, '0')
	return nothing
end

"""
Sets the sleep mode off
"""
function sleepModeOff(gauss::LakeShoreGaussMeter)
	setSleepMode(gauss, '1')
	return nothing
end

"""
Locks the frontpanel
"""
function lockOn(gauss::LakeShoreGaussMeter)
	setFrontPanelLock(gauss, '1')
	return nothing
end

"""
Unlocks the frontpanel
"""
function lockOff(gauss::LakeShoreGaussMeter)
	setFrontPanelLock(gauss, '0')
	return nothing
end

"""
Sets the unit of the values to gauss
"""
function setUnitToGauss(gauss::LakeShoreGaussMeter)
	setUnit(gauss, 'G')
	return nothing
end

"""
Sets the unit of the values to tesla
"""
function setUnitToTesla(gauss::LakeShoreGaussMeter)
	setUnit(gauss, 'T')
	return nothing
end

"""
Sets the standard settings
	-highest range
	-unit to tesla
	-auto ranging off
	-complete probe on
"""
function setStandardSettings(gauss::LakeShoreGaussMeter)
	#setAllRange(gauss, '0')
	#setAllMode(gauss, '0')
	setAllAutoRanging(gauss, '1')
	setUnitToTesla(gauss)
	setCompleteProbe(gauss, '0')
	@info "Standard Settings: Unit = Tesla, Range = lowest, Mode = DC, AutoRanging = off, Probe = on"
	return nothing
end


text:
function generate_central_moment_eqs(drift_eqs::AbstractVector{Equation}, diff::AbstractArray,
                                      m_order::Int, q_order::Int, vars, name, ps, iv)

    N = length(drift_eqs)
    drift = [e.rhs for e in drift_eqs]
    diff_mat = diff*transpose(diff)

    if iszero(q_order)
        try
            smap = Dict(Pair.(vars, 1:length(vars)))
            _, _, drift_order = polynomial_propensities(drift, iv, smap)
            _, _, diff_order = polynomial_propensities(value.(diff_mat), iv, smap)
            q_order = m_order + max(drift_order - 1, diff_order - 2)
        catch e
            error("non-polynomial rates (OR A BUG): please specify q_order.\n" * string(e))
        end
    end

    iter_all = construct_iter_all(N, q_order)
    iter_m = filter(x -> 1 < sum(x) <= m_order, iter_all)
    iter_q = filter(x -> m_order < sum(x) <= q_order, iter_all)
    iter_1 = filter(x -> sum(x) == 1, iter_all)

    μ = define_μ(iter_1, iv)
    M = define_M(iter_all, iv)

    dict_vars_to_μ = Dict(zip(vars, values(μ)))

    Df = Vector{Dict}(undef, N)
    derivs = Dict()
    for i in 1:N
        for iter in iter_all
            derivs[iter] = drift[i]
            for j in 1:N
                for d in 1:iter[j]
                    derivs[iter] = expand_derivatives(Differential(vars[j])(derivs[iter]), true)
                end
            end
            derivs[iter] = substitute(derivs[iter], dict_vars_to_μ)
        end
        Df[i] = copy(derivs)
    end

    Dg = Matrix{Dict}(undef, N, N)
    derivs = Dict()
    for i in 1:N
        for j in 1:N
            for iter in iter_all
                derivs[iter] = diff_mat[i, j]
                for k in 1:N
                    for d in 1:iter[k]
                        derivs[iter] = expand_derivatives(Differential(vars[k])(derivs[iter]), true)
                    end
                end
                derivs[iter] = substitute(derivs[iter], dict_vars_to_μ)
            end
            Dg[i, j] = copy(derivs)
        end
    end

    eqs = Equation[]
    for i in 1:N
        poly = sum( Df[i][j]*M[j]//fact(j) for j in iter_all )
        push!(eqs, Differential(iv)(μ[iter_1[i]]) ~ expand(poly))
    end

    for iter in iter_m

        iterq2 = filter(x -> sum(x) <= q_order - sum(iter) + 2, iter_all)
        iterq1 = filter(x -> sum(x) <= q_order - sum(iter) + 1, iterq2)

        term1 = 0
        for i in 1:N
            iszero(iter[i]) ? continue :
            gradf = iter[i]
            gradh = iter .- iter_1[i]
            term1 += sum( gradf * Df[i][j] * M[gradh .+ j] // fact(j) for j in iterq1 )
            term1 -= sum( gradf * Df[i][j] * M[gradh] * M[j] // fact(j) for j in iter_all )
        end

        term2 = 0
        for i in 1:N
            for j in 1:N
                ( i == j && isone(iter[i]) ) || (iszero(iter[i]) || iszero(iter[j])) ? continue :
                hji = iter .- (iter_1[i] .+ iter_1[j])
                hji_factor = i == j ? iter[i]*(iter[i]-1) : iter[i]*iter[j]

                term2 += hji_factor * sum( Dg[i, j][k] * M[hji .+ k] // fact(k) for k in iterq2 )
            end
        end

        push!(eqs, Differential(iv)(M[iter]) ~ expand(term1 + term2/2))

    end

    odename = Symbol(name, "_central_moment_eqs_m", m_order, "_q", q_order)
    odes = ODESystem(eqs, iv, extract_variables(eqs, μ, M), ps; name=odename)

    CentralMomentEquations(odes, μ, M, N, m_order, q_order, iter_all, iter_m, iter_q, iter_1)
end

generate_central_moment_eqs(sys::SDESystem, m_order::Int, q_order::Int=0) = generate_central_moment_eqs(equations(sys), get_noiseeqs(sys), m_order, q_order,
                                                                                states(sys), nameof(sys), parameters(sys), get_iv(sys))

#=
function generate_central_moment_eqs(drift_eqs::AbstractVector{Equation}, diff::AbstractArray{T},
                                     m_order, approx_order, vars, name, ps, iv) where T <: Union{Sym, Num}

    N = length(drift_eqs)
    drift = [e.rhs for e in drift_eqs]
    diffusion_matrix = diff*transpose(diff)

    # if approx_order == 0 => assume that all drift_eqs are polynomials
    if iszero(approx_order)
        approx_order = max(maximum(degree(eq, vars) for eq in drift), maximum(degree(eq, vars) for eq in diffusion_matrix))
        q_order = m_order + max(maximum(degree(eq, vars) for eq in drift) - 1, maximum(degree(eq, vars) for eq in diffusion_matrix) - 2)
    else
        q_order = m_order + approx_order - 1
    end

    iter_all = construct_iter_all(N, q_order)
    iter_m = filter(x -> 1 < sum(x) <= m_order, iter_all)
    iter_q = filter(x -> m_order < sum(x) <= q_order, iter_all)
    iter_1 = filter(x -> sum(x) == 1, iter_all)

    μ = define_μ(iter_1, iv)
    μ_vec = [value(μ[iter]) for iter in iter_1]
    M = define_M(iter_all, iv)

    # determine polynomial approximation
    # TODO: compare the speed of this + poly_subs to what I have in central_moment_eqs.jl
    poly_drift = [taylor_expand(f, vars, μ_vec, approx_order, false) for f in drift]
    poly_diffusion_matrix = [taylor_expand(f, vars, μ_vec, approx_order, false) for f in diffusion_matrix]

    # determine first order moment equations
    mom_eqs = Equation[]
    mono_to_moment = OrderedDict(prod(vars.^iter) => M[iter] for iter in iter_all)
    moms = Num[]
    for iter in iter_1
        mono = sum(vars .* iter)
        push!(moms, μ[iter])
        push!(mom_eqs, Differential(iv)(μ[iter]) ~ poly_subs(sum(poly_drift .* gradient(mono, vars)), mono_to_moment, ps, true))
    end

    for iter in iter_m
        mono = prod(vars .^ iter)
        push!(moms, M[iter])
        poly = sum(poly_drift .* gradient(mono, vars)) + 1/2 * sum( (poly_diffusion_matrix*hessian(mono, vars))[i,i] for i in 1:N )
        offset = 0
        for i in 1:N
            if !iszero(iter[i])
                iter_incr = Tuple(i != j ? iter[j] : iter[j] - 1 for j in 1:N)
                offset += M[iter_incr]*(iter[i] + 1)*mom_eqs[i].rhs
            end
        end

        push!(mom_eqs, Differential(iv)(M[iter]) ~ - offset + poly_subs(poly, mono_to_moment, Set(ps), true))
    end

    return CentralMomentEquations(ODESystem(mom_eqs, iv, moms, ps, name = Symbol(name, :_RAW_CENTRAL)), μ, M, N, m_order, q_order, iter_all, iter_m, iter_q, iter_1)
end

generate_central_moment_eqs(sys::SDESystem, m_order::Int, q_order::Int=0) = generate_central_moment_eqs(equations(sys), get_noiseeqs(sys), m_order, q_order,
                                                                                states(sys), nameof(sys), parameters(sys), get_iv(sys))
=#


text:
using DataFrames
using DataConvenience: @>
using TidyStanza

# df <- tibble(a = 1, b = 1, c = 1, d = "a", e = "a", f = "a")

df = DataFrame(a = 1, b = 1, c = 1, d = "a", e = "a", f = "a")

# df %>% relocate(f)
@> df relocate(:f)

# df %>% relocate(a, .after = c)
@> df relocate(:a, after = :c)

# df %>% relocate(f, .before = b)
@> df relocate(:f, before = :b)

# df %>% relocate(a, .after = last_col())
@> df relocate(:a, after = names(df)[end])

last_col() = df->names(df)[end]

@> df relocate(:a, after = last_col())

middle_col() = df->names(df)[end ÷ 2]
@> df relocate(:a, after = middle_col())

using TidyStanza: where

# df %>% relocate(where(is.character))
isstring(x) = eltype(x) <: AbstractString

@> df relocate(where(isstring))

@> df relocate(where(x->eltype(x) <: AbstractString))

# df %>% relocate(where(is.numeric), .after = last_col())
isnumeric(x) = eltype(x) <: Number

@> df relocate(where(isnumeric), after = last_col())

# df %>% relocate(any_of(c("a", "e", "i", "o", "u")))
@> df relocate(intersect(["a", "e", "i", "o", "u"], names(df)))

any_of(cols::AbstractVector{T}) where T = df -> intersect(T.(names(df)), cols)

@> df relocate(any_of(["a", "e", "i", "o", "u"]))

#df2 <- tibble(a = 1, b = "a", c = 1, d = "a")

df2 = DataFrame(a = 1, b = "a", c = 1, d = "a")

#df2 %>% relocate(where(is.numeric), .after = where(is.character))

@> df2 relocate(where(isnumeric), after = where(isstring))

#df2 %>% relocate(where(is.numeric), .before = where(is.character))
@> df2 relocate(where(isnumeric), before = where(isstring))


text:
using Test
using IICGP
using CartesianGeneticProgramming
using Cambrian

include("atari.jl")
include("config.jl")
include("determinism.jl")
include("functions.jl")
include("individual.jl")
include("process.jl")
include("reducer.jl")
include("evolution.jl")
include("evaluation.jl")


text:
using Documenter, JetPackWaveFD

makedocs(
    sitename="JetPackWaveFD.jl",
    modules=[JetPackWaveFD],
    pages = [ "index.md", "reference.md" ]
)

deploydocs(
    repo = "github.com/ChevronETC/JetPackWaveFD.jl.git"
)

max_stars_repo_name: ChevronETC/JetPackWaveFD

max_issues_repo_name: ChevronETC/JetPackWaveFD

max_forks_repo_name: ChevronETC/JetPackWaveFD


text:
pythagoras(x, y) = √(x^2 + y^2)

println("----Starting Output from script!----")
println("No recompiling necessary!")
println(pythagoras(3, 4))
println(pythagoras(5, 12))
println("-----Ending Output from script!-----")

max_stars_repo_name: cdsousa/Juliapp

max_issues_repo_name: cdsousa/Juliapp

max_forks_repo_name: cdsousa/Juliapp