text:
# https://github.com/stevengj/PyPlot.jl

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

# convert colorant to 4-tuple RGBA
getPyPlotColor(c::Colorant, α=nothing) = map(f->float(f(convertColor(c,α))), (red, green, blue, alpha))
getPyPlotColor(scheme::ColorScheme, α=nothing) = getPyPlotColor(convertColor(getColor(scheme), α))
getPyPlotColor(c, α=nothing) = getPyPlotColor(convertColor(c, α))
# getPyPlotColor(c, alpha) = getPyPlotColor(colorscheme(c, alpha))

function getPyPlotColorMap(c::ColorGradient, α=nothing)
  # c = ColorGradient(c.colors, c.values, alpha=α)
  # pycolors.pymember("LinearSegmentedColormap")[:from_list]("tmp", map(getPyPlotColor, getColorVector(c)))
  # pyvals = [(c.values[i], getPyPlotColor(c.colors[i], α)) for i in 1:length(c.colors)]
  pyvals = [(v, getPyPlotColor(getColorZ(c, v), α)) for v in c.values]
  # @show c α pyvals
  pycolors.pymember("LinearSegmentedColormap")[:from_list]("tmp", pyvals)
end

# anything else just gets a redsblue gradient
getPyPlotColorMap(c, α=nothing) = getPyPlotColorMap(ColorGradient(:redsblues), α)

# get the style (solid, dashed, etc)
function getPyPlotLineStyle(linetype::Symbol, linestyle::Symbol)
  linetype == :none && return " "
  linestyle == :solid && return "-"
  linestyle == :dash && return "--"
  linestyle == :dot && return ":"
  linestyle == :dashdot && return "-."
  warn("Unknown linestyle $linestyle")
  return "-"
end

function getPyPlotMarker(marker::Shape)
  n = length(marker.vertices)
  mat = zeros(n+1,2)
  for (i,vert) in enumerate(marker.vertices)
    mat[i,1] = vert[1]
    mat[i,2] = vert[2]
  end
  mat[n+1,:] = mat[1,:]
  pypath.pymember("Path")(mat)
  # marker.vertices
end

# get the marker shape
function getPyPlotMarker(marker::Symbol)
  marker == :none && return " "
  marker == :ellipse && return "o"
  marker == :rect && return "s"
  marker == :diamond && return "D"
  marker == :utriangle && return "^"
  marker == :dtriangle && return "v"
  marker == :cross && return "+"
  marker == :xcross && return "x"
  marker == :star5 && return "*"
  marker == :pentagon && return "p"
  marker == :hexagon && return "h"
  marker == :octagon && return "8"
  haskey(_shapes, marker) && return getPyPlotMarker(_shapes[marker])

  warn("Unknown marker $marker")
  return "o"
end

# pass through
function getPyPlotMarker(marker::@compat(AbstractString))
  @assert length(marker) == 1
  marker
end

function getPyPlotStepStyle(linetype::Symbol)
  linetype == :steppost && return "steps-post"
  linetype == :steppre && return "steps-pre"
  return "default"
end

# immutable PyPlotFigWrapper
#   fig
#   kwargs  # for add_subplot
# end

type PyPlotAxisWrapper
  ax
  rightax
  fig
  kwargs  # for add_subplot
end

# getfig(wrap::@compat(Union{PyPlotAxisWrapper,PyPlotFigWrapper})) = wrap.fig
getfig(wrap::PyPlotAxisWrapper) = wrap.fig

# get a reference to the correct axis
function getLeftAxis(wrap::PyPlotAxisWrapper)
  if wrap.ax == nothing
    axes = wrap.fig.o[:axes]
    if isempty(axes)
      return wrap.fig.o[:add_subplot](111; wrap.kwargs...)
    end
    axes[1]
  else
    wrap.ax
  end
end
# getLeftAxis(wrap::PyPlotAxisWrapper) = wrap.ax
# getRightAxis(x) = getLeftAxis(x)[:twinx]()

function getRightAxis(wrap::PyPlotAxisWrapper)
  if wrap.rightax == nothing
    wrap.rightax = getLeftAxis(wrap)[:twinx]()
  end
  wrap.rightax
end

getLeftAxis(plt::Plot{PyPlotPackage}) = getLeftAxis(plt.o)
getRightAxis(plt::Plot{PyPlotPackage}) = getRightAxis(plt.o)
getAxis(plt::Plot{PyPlotPackage}, axis::Symbol) = (axis == :right ? getRightAxis : getLeftAxis)(plt)

# left axis is PyPlot.<func>, right axis is "f.axes[0].twinx().<func>"
function getPyPlotFunction(plt::Plot, axis::Symbol, linetype::Symbol)

  # # need to access mplot3d functions differently
  # if linetype == :surface
  #   return mplot3d.pymember("Axes3D")[:plot_surface]
  # end

  # in the 2-axis case we need to get: <rightaxis>[:<func>]
  ax = getAxis(plt, axis)
  # ax[:set_ylabel](plt.plotargs[:yrightlabel])
  fmap = @compat Dict(
      :hist       => :hist,
      :density    => :hist,
      :sticks     => :bar,
      :bar        => :bar,
      :heatmap    => :hexbin,
      :hexbin     => :hexbin,
      :scatter    => :scatter,
      :contour    => :contour,
      :scatter3d  => :scatter,
      :surface    => :plot_surface,
      :wireframe  => :plot_wireframe,
      # :surface    => pycolors.pymember("LinearSegmentedColormap")[:from_list]
    )
  return ax[get(fmap, linetype, :plot)]
end

function updateAxisColors(ax, fgcolor)
  for loc in ("bottom", "top", "left", "right")
    ax[:spines][loc][:set_color](fgcolor)
  end
  for axis in ("x", "y")
    ax[:tick_params](axis=axis, colors=fgcolor, which="both")
  end
  for axis in (:yaxis, :xaxis)
    ax[axis][:label][:set_color](fgcolor)
  end
  ax[:title][:set_color](fgcolor)
end

function handleSmooth(plt::Plot{PyPlotPackage}, ax, d::Dict, smooth::Bool)
  if smooth
    xs, ys = regressionXY(d[:x], d[:y])
    ax[:plot](xs, ys,
              # linestyle = getPyPlotLineStyle(:path, :dashdot),
              color = getPyPlotColor(d[:linecolor]),
              linewidth = 2
             )
  end
end
handleSmooth(plt::Plot{PyPlotPackage}, ax, d::Dict, smooth::Real) = handleSmooth(plt, ax, d, true)

# makePyPlotCurrent(wrap::PyPlotFigWrapper) = PyPlot.figure(wrap.fig.o[:number])
# makePyPlotCurrent(wrap::PyPlotAxisWrapper) = nothing #PyPlot.sca(wrap.ax.o)
makePyPlotCurrent(wrap::PyPlotAxisWrapper) = wrap.ax == nothing ? PyPlot.figure(wrap.fig.o[:number]) : nothing
makePyPlotCurrent(plt::Plot{PyPlotPackage}) = plt.o == nothing ? nothing : makePyPlotCurrent(plt.o)

function _before_add_series(plt::Plot{PyPlotPackage})
  makePyPlotCurrent(plt)
end

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

# TODO:
# fillto   # might have to use barHack/histogramHack??
# reg             # true or false, add a regression line for each line
# pos             # (Int,Int), move the enclosing window to this position
# windowtitle     # string or symbol, set the title of the enclosing windowtitle
# screen          # Integer, move enclosing window to this screen number (for multiscreen desktops)
# show            # true or false, show the plot (in case you don't want the window to pop up right away)

function _create_plot(pkg::PyPlotPackage; kw...)
  # create the figure
  d = Dict(kw)

  # standalone plots will create a figure, but not if part of a subplot (do it later)
  if haskey(d, :subplot)
    wrap = nothing
  else
    w,h = map(px2inch, d[:size])
    bgcolor = getPyPlotColor(d[:background_color])
    wrap = PyPlotAxisWrapper(nothing, nothing, PyPlot.figure(; figsize = (w,h), facecolor = bgcolor, dpi = DPI, tight_layout = true), [])

    if haskey(d, :linetype) && first(d[:linetype]) in _3dTypes # && isa(plt.o, PyPlotFigWrapper)
      push!(wrap.kwargs, (:projection, "3d"))
    end
  end

  plt = Plot(wrap, pkg, 0, d, Dict[])
  plt
end

function _add_series(pkg::PyPlotPackage, plt::Plot; kw...)
  d = Dict(kw)

  lt = d[:linetype]
  if lt in _3dTypes && isempty(plt.o.kwargs)
    push!(plt.o.kwargs, (:projection, "3d"))
  end

  ax = getAxis(plt, d[:axis])
  if !(lt in supportedTypes(pkg))
    error("linetype $(lt) is unsupported in PyPlot.  Choose from: $(supportedTypes(pkg))")
  end

  color = getPyPlotColor(d[:linecolor], d[:linealpha])

  if lt == :sticks
    d,_ = sticksHack(;d...)

  elseif lt in (:scatter, :scatter3d)
    if d[:markershape] == :none
      d[:markershape] = :ellipse
    end

  elseif lt in (:hline,:vline)
    linewidth = d[:linewidth]
    linecolor = color
    linestyle = getPyPlotLineStyle(lt, d[:linestyle])
    for yi in d[:y]
      func = ax[lt == :hline ? :axhline : :axvline]
      func(yi, linewidth=d[:linewidth], color=linecolor, linestyle=linestyle)
    end

  end

  # lt = d[:linetype]
  extra_kwargs = Dict()

  plotfunc = getPyPlotFunction(plt, d[:axis], lt)

  # we have different args depending on plot type
  if lt in (:hist, :density, :sticks, :bar)

    # NOTE: this is unsupported because it does the wrong thing... it shifts the whole axis
    # extra_kwargs[:bottom] = d[:fill]

    if ishistlike(lt)
      extra_kwargs[:bins] = d[:nbins]
      extra_kwargs[:normed] = lt == :density
    else
      extra_kwargs[:linewidth] = (lt == :sticks ? 0.1 : 0.9)
    end

  elseif lt in (:heatmap, :hexbin)
    extra_kwargs[:gridsize] = d[:nbins]
    extra_kwargs[:cmap] = getPyPlotColorMap(d[:linecolor])

  elseif lt == :contour
    extra_kwargs[:cmap] = getPyPlotColorMap(d[:linecolor])
    extra_kwargs[:linewidths] = d[:linewidth]
    extra_kwargs[:linestyles] = getPyPlotLineStyle(lt, d[:linestyle])
    # TODO: will need to call contourf to fill in the contours

  elseif lt in (:surface, :wireframe)
    if lt == :surface
      extra_kwargs[:cmap] = getPyPlotColorMap(d[:fillcolor], d[:fillalpha])
    end
    extra_kwargs[:rstride] = 1
    extra_kwargs[:cstride] = 1
    extra_kwargs[:linewidth] = d[:linewidth]
    extra_kwargs[:edgecolor] = getPyPlotColor(d[:linecolor], d[:linealpha])

  else

    extra_kwargs[:linestyle] = getPyPlotLineStyle(lt, d[:linestyle])
    extra_kwargs[:marker] = getPyPlotMarker(d[:markershape])

    if lt in (:scatter, :scatter3d)
      extra_kwargs[:s] = d[:markersize].^2
      c = d[:markercolor]
      if isa(c, ColorGradient) && d[:zcolor] != nothing
        extra_kwargs[:c] = convert(Vector{Float64}, d[:zcolor])
        extra_kwargs[:cmap] = getPyPlotColorMap(c, d[:markeralpha])
      else
        extra_kwargs[:c] = getPyPlotColor(c, d[:markeralpha])
      end
      if d[:markeralpha] != nothing
        extra_kwargs[:alpha] = d[:markeralpha]
      end
      extra_kwargs[:edgecolors] = getPyPlotColor(d[:markerstrokecolor], d[:markerstrokealpha])
      extra_kwargs[:linewidths] = d[:markerstrokewidth]
    else
      extra_kwargs[:markersize] = d[:markersize]
      extra_kwargs[:markerfacecolor] = getPyPlotColor(d[:markercolor], d[:markeralpha])
      extra_kwargs[:markeredgecolor] = getPyPlotColor(d[:markerstrokecolor], d[:markerstrokealpha])
      extra_kwargs[:markeredgewidth] = d[:markerstrokewidth]
      extra_kwargs[:drawstyle] = getPyPlotStepStyle(lt)
    end
  end

  # if d[:markeralpha] != nothing
  #   extra_kwargs[:alpha] = d[:markeralpha]
  # elseif d[:linealpha] != nothing
  #   extra_kwargs[:alpha] = d[:linealpha]
  # end

  # set these for all types
  if !(lt in (:contour,:surface,:wireframe))
    if !(lt in (:scatter, :scatter3d))
      extra_kwargs[:color] = color
      extra_kwargs[:linewidth] = d[:linewidth]
    end
    extra_kwargs[:label] = d[:label]
    extra_kwargs[:zorder] = plt.n
  end

  # do the plot
  d[:serieshandle] = if ishistlike(lt)
    plotfunc(d[:y]; extra_kwargs...)[1]

  elseif lt == :contour
    x, y = d[:x], d[:y]
    surf = d[:z].surf'
    levels = d[:levels]
    if isscalar(levels)
      extra_args = (levels)
    elseif isvector(levels)
      extra_args = ()
      extra_kwargs[:levels] = levels
    else
      error("Only numbers and vectors are supported with levels keyword")
    end
    handle = plotfunc(x, y, surf, extra_args...; extra_kwargs...)
    if d[:fillrange] != nothing
      extra_kwargs[:cmap] = getPyPlotColorMap(d[:fillcolor], d[:fillalpha])
      delete!(extra_kwargs, :linewidths)
      handle = ax[:contourf](x, y, surf, extra_args...; extra_kwargs...)
    end
    handle

  elseif lt in (:surface,:wireframe)
    x, y, z = d[:x], d[:y], Array(d[:z])
    if !ismatrix(x) || !ismatrix(y)
      x = repmat(x', length(y), 1)
      y = repmat(y, 1, length(d[:x]))
      z = z'
    end
    plotfunc(x, y, z; extra_kwargs...)
  elseif lt in _3dTypes
    plotfunc(d[:x], d[:y], d[:z]; extra_kwargs...)
  elseif lt in (:scatter, :heatmap, :hexbin)
    plotfunc(d[:x], d[:y]; extra_kwargs...)
  else
    plotfunc(d[:x], d[:y]; extra_kwargs...)[1]
  end

  handleSmooth(plt, ax, d, d[:smooth])

  # add the colorbar legend
  if plt.plotargs[:legend] && haskey(extra_kwargs, :cmap)
    PyPlot.colorbar(d[:serieshandle])
  end

  # @show extra_kwargs

  # this sets the bg color inside the grid
  ax[:set_axis_bgcolor](getPyPlotColor(plt.plotargs[:background_color]))

  fillrange = d[:fillrange]
  if fillrange != nothing && lt != :contour
    fillcolor = getPyPlotColor(d[:fillcolor], d[:fillalpha])
    if typeof(fillrange) <: @compat(Union{Real, AVec})
      ax[:fill_between](d[:x], fillrange, d[:y], facecolor = fillcolor, zorder = plt.n)
    else
      ax[:fill_between](d[:x], fillrange..., facecolor = fillcolor, zorder = plt.n)
    end
  end

  push!(plt.seriesargs, d)
  plt
end

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

function Base.getindex(plt::Plot{PyPlotPackage}, i::Integer)
  series = plt.seriesargs[i][:serieshandle]
  try
    return series[:get_data]()
  catch
    xy = series[:get_offsets]()
    return vec(xy[:,1]), vec(xy[:,2])
  end
  # series[:relim]()
  # mapping = getGadflyMappings(plt, i)[1]
  # mapping[:x], mapping[:y]
end

function Base.setindex!{X,Y}(plt::Plot{PyPlotPackage}, xy::Tuple{X,Y}, i::Integer)
  series = plt.seriesargs[i][:serieshandle]
  x, y = xy
  try
    series[:set_data](x, y)
  catch
    series[:set_offsets](hcat(x, y))
  end

  ax = series[:axes]
  if plt.plotargs[:xlims] == :auto
    xmin, xmax = ax[:get_xlim]()
    ax[:set_xlim](min(xmin, minimum(x)), max(xmax, maximum(x)))
  end
  if plt.plotargs[:ylims] == :auto
    ymin, ymax = ax[:get_ylim]()
    ax[:set_ylim](min(ymin, minimum(y)), max(ymax, maximum(y)))
  end

  plt
end

function Base.setindex!{X,Y,Z}(plt::Plot{PyPlotPackage}, xyz::Tuple{X,Y,Z}, i::Integer)
  warn("setindex not implemented for xyz")
  plt
end

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

function addPyPlotLims(ax, lims, isx::Bool)
  lims == :auto && return
  ltype = limsType(lims)
  if ltype == :limits
    ax[isx ? :set_xlim : :set_ylim](lims...)
  else
    error("Invalid input for $(isx ? "xlims" : "ylims"): ", lims)
  end
end

function addPyPlotTicks(ax, ticks, isx::Bool)
  ticks == :auto && return
  if ticks == :none || ticks == nothing
    ticks = zeros(0)
  end

  ttype = ticksType(ticks)
  if ttype == :ticks
    ax[isx ? :set_xticks : :set_yticks](ticks)
  elseif ttype == :ticks_and_labels
    ax[isx ? :set_xticks : :set_yticks](ticks...)
  else
    error("Invalid input for $(isx ? "xticks" : "yticks"): ", ticks)
  end
end

usingRightAxis(plt::Plot{PyPlotPackage}) = any(args -> args[:axis] in (:right,:auto), plt.seriesargs)

function _update_plot(plt::Plot{PyPlotPackage}, d::Dict)
  figorax = plt.o
  ax = getLeftAxis(figorax)
  # PyPlot.sca(ax)

  # title and axis labels
  # haskey(d, :title) && PyPlot.title(d[:title])
  haskey(d, :title) && ax[:set_title](d[:title])
  haskey(d, :xlabel) && ax[:set_xlabel](d[:xlabel])
  if haskey(d, :ylabel)
    ax[:set_ylabel](d[:ylabel])
  end
  if usingRightAxis(plt) && get(d, :yrightlabel, "") != ""
    rightax = getRightAxis(figorax)
    rightax[:set_ylabel](d[:yrightlabel])
  end

  # scales
  haskey(d, :xscale) && applyPyPlotScale(ax, d[:xscale], true)
  haskey(d, :yscale) && applyPyPlotScale(ax, d[:yscale], false)

  # limits and ticks
  haskey(d, :xlims) && addPyPlotLims(ax, d[:xlims], true)
  haskey(d, :ylims) && addPyPlotLims(ax, d[:ylims], false)
  haskey(d, :xticks) && addPyPlotTicks(ax, d[:xticks], true)
  haskey(d, :yticks) && addPyPlotTicks(ax, d[:yticks], false)

  if get(d, :xflip, false)
    ax[:invert_xaxis]()
  end
  if get(d, :yflip, false)
    ax[:invert_yaxis]()
  end

  axes = [getLeftAxis(figorax)]
  if usingRightAxis(plt)
    push!(axes, getRightAxis(figorax))
  end

  # font sizes
  for ax in axes
    # haskey(d, :yrightlabel) || continue

    # guides
    sz = get(d, :guidefont, plt.plotargs[:guidefont]).pointsize
    ax[:title][:set_fontsize](sz)
    ax[:xaxis][:label][:set_fontsize](sz)
    ax[:yaxis][:label][:set_fontsize](sz)

    # ticks
    sz = get(d, :tickfont, plt.plotargs[:tickfont]).pointsize
    for sym in (:get_xticklabels, :get_yticklabels)
      for lab in ax[sym]()
        lab[:set_fontsize](sz)
      end
    end

    # grid
    if get(d, :grid, false)
      ax[:xaxis][:grid](true)
      ax[:yaxis][:grid](true)
      ax[:set_axisbelow](true)
    end
  end

end

function applyPyPlotScale(ax, scaleType::Symbol, isx::Bool)
  func = ax[isx ? :set_xscale : :set_yscale]
  scaleType == :identity && return func("linear")
  scaleType == :ln && return func("log", basex = e, basey = e)
  scaleType == :log2 && return func("log", basex = 2, basey = 2)
  scaleType == :log10 && return func("log", basex = 10, basey = 10)
  warn("Unhandled scaleType: ", scaleType)
end

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

function createPyPlotAnnotationObject(plt::Plot{PyPlotPackage}, x, y, val::@compat(AbstractString))
  ax = getLeftAxis(plt)
  ax[:annotate](val, xy = (x,y))
end

function createPyPlotAnnotationObject(plt::Plot{PyPlotPackage}, x, y, val::PlotText)
  ax = getLeftAxis(plt)
  ax[:annotate](val.str,
    xy = (x,y),
    family = val.font.family,
    color = getPyPlotColor(val.font.color),
    horizontalalignment = val.font.halign == :hcenter ? "center" : string(val.font.halign),
    verticalalignment = val.font.valign == :vcenter ? "center" : string(val.font.valign),
    rotation = val.font.rotation * 180 / π,
    size = val.font.pointsize
  )
end

function _add_annotations{X,Y,V}(plt::Plot{PyPlotPackage}, anns::AVec{@compat(Tuple{X,Y,V})})
  for ann in anns
    createPyPlotAnnotationObject(plt, ann...)
  end
end

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

# NOTE: pyplot needs to build before
function _create_subplot(subplt::Subplot{PyPlotPackage}, isbefore::Bool)
  l = subplt.layout

  w,h = map(px2inch, getplotargs(subplt,1)[:size])
  bgcolor = getPyPlotColor(getplotargs(subplt,1)[:background_color])
  fig = PyPlot.figure(; figsize = (w,h), facecolor = bgcolor, dpi = DPI, tight_layout = true)

  nr = nrows(l)
  for (i,(r,c)) in enumerate(l)

    # add the plot to the figure
    nc = ncols(l, r)
    fakeidx = (r-1) * nc + c
    ax = fig[:add_subplot](nr, nc, fakeidx)

    subplt.plts[i].o = PyPlotAxisWrapper(ax, nothing, fig, [])
  end

  # subplt.o = PyPlotFigWrapper(fig, [])
  subplt.o = PyPlotAxisWrapper(nothing, nothing, fig, [])
  true
end

# this will be called internally, when creating a subplot from existing plots
# NOTE: if I ever need to "Rebuild a "ubplot from individual Plot's"... this is what I should use!
function subplot(plts::AVec{Plot{PyPlotPackage}}, layout::SubplotLayout, d::Dict)
  validateSubplotSupported()

  p = length(layout)
  n = sum([plt.n for plt in plts])

  pkg = PyPlotPackage()
  newplts = Plot{PyPlotPackage}[_create_plot(pkg; subplot=true, plt.plotargs...) for plt in plts]

  subplt = Subplot(nothing, newplts, PyPlotPackage(), p, n, layout, d, true, false, false, (r,c) -> (nothing,nothing))

  _preprocess_subplot(subplt, d)
  _create_subplot(subplt, true)

  for (i,plt) in enumerate(plts)
    for seriesargs in plt.seriesargs
      _add_series_subplot(newplts[i]; seriesargs...)
    end
  end

  _postprocess_subplot(subplt, d)

  subplt
end

function _remove_axis(plt::Plot{PyPlotPackage}, isx::Bool)
  if isx
    plot!(plt, xticks=zeros(0), xlabel="")
  else
    plot!(plt, yticks=zeros(0), ylabel="")
  end
end

function _expand_limits(lims, plt::Plot{PyPlotPackage}, isx::Bool)
  pltlims = plt.o.ax[isx ? :get_xbound : :get_ybound]()
  _expand_limits(lims, pltlims)
end

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

# function addPyPlotLegend(plt::Plot)
function addPyPlotLegend(plt::Plot, ax)
  if plt.plotargs[:legend]
    # gotta do this to ensure both axes are included
    args = filter(x -> !(x[:linetype] in (:hist,:density,:hexbin,:heatmap,:hline,:vline,:contour, :surface, :wireframe, :path3d, :scatter3d)), plt.seriesargs)
    if length(args) > 0
      leg = ax[:legend]([d[:serieshandle] for d in args],
                  [d[:label] for d in args],
                  loc="best",
                  fontsize = plt.plotargs[:legendfont].pointsize
                  # framealpha = 0.6
                 )
      leg[:set_zorder](1000)
    end
  end
end

function finalizePlot(plt::Plot{PyPlotPackage})
  ax = getLeftAxis(plt)
  addPyPlotLegend(plt, ax)
  updateAxisColors(ax, getPyPlotColor(plt.plotargs[:foreground_color]))
  PyPlot.draw()
end

function finalizePlot(subplt::Subplot{PyPlotPackage})
  fig = subplt.o.fig
  for (i,plt) in enumerate(subplt.plts)
    ax = getLeftAxis(plt)
    addPyPlotLegend(plt, ax)
    updateAxisColors(ax, getPyPlotColor(plt.plotargs[:foreground_color]))
  end
  # fig[:tight_layout]()
  PyPlot.draw()
end

# # allow for writing any supported mime
# for mime in keys(PyPlot.aggformats)
#   @eval function Base.writemime(io::IO, m::MIME{symbol{$mime}}, plt::Plot{PyPlotPackage})
#     finalizePlot(plt)
#     writemime(io, m, getfig(plt.o))
#   end
# end

# function Base.writemime(io::IO, m::@compat(Union{MIME"image/svg+xml", MIME"image/png"}, plt::Plot{PyPlotPackage})
#   finalizePlot(plt)
#   writemime(io, m, getfig(plt.o))
# end

# NOTE: to bring up a GUI window in IJulia, need some extra steps
function Base.display(::PlotsDisplay, plt::PlottingObject{PyPlotPackage})
  finalizePlot(plt)
  if isa(Base.Multimedia.displays[end], Base.REPL.REPLDisplay)
    display(getfig(plt.o))
  else
    # # PyPlot.ion()
    # PyPlot.figure(getfig(plt.o).o[:number])
    # PyPlot.draw_if_interactive()
    # # PyPlot.ioff()
  end
  # PyPlot.plt[:show](block=false)
  getfig(plt.o)[:show]()
end

# function Base.display(::PlotsDisplay, subplt::Subplot{PyPlotPackage})
#   finalizePlot(subplt)
#   PyPlot.ion()
#   PyPlot.figure(getfig(subplt.o).o[:number])
#   PyPlot.draw_if_interactive()
#   PyPlot.ioff()
#   # display(getfig(subplt.o))
# end

# # allow for writing any supported mime
# for mime in (MIME"image/png", MIME"application/pdf", MIME"application/postscript")
#   @eval function Base.writemime(io::IO, ::$mime, plt::PlottingObject{PyPlotPackage})
#     finalizePlot(plt)
#     writemime(io, $mime(), getfig(plt.o))
#   end
# end

const _pyplot_mimeformats = @compat Dict(
    "application/eps"         => "eps",
    "image/eps"               => "eps",
    "application/pdf"         => "pdf",
    "image/png"               => "png",
    "application/postscript"  => "ps",
    # "image/svg+xml"           => "svg"
  )

for (mime, fmt) in _pyplot_mimeformats
  @eval function Base.writemime(io::IO, ::MIME{symbol($mime)}, plt::PlottingObject{PyPlotPackage})
    finalizePlot(plt)
    fig = getfig(plt.o)
    fig.o["canvas"][:print_figure](io,
                                   format=$fmt,
                                   # bbox_inches = "tight",
                                   # figsize = map(px2inch, plt.plotargs[:size]),
                                   facecolor = fig.o["get_facecolor"](),
                                   edgecolor = "none",
                                   dpi = DPI
                                  )
  end
end

# function Base.writemime(io::IO, m::MIME"image/png", subplt::Subplot{PyPlotPackage})
#   finalizePlot(subplt)
#   writemime(io, m, getfig(subplt.o))
# end


text:
function AntexLoader(arg0::JString)
    return AntexLoader((JString,), arg0)
end

function AntexLoader(arg0::JString, arg1::DataProvidersManager, arg2::TimeScale)
    return AntexLoader((JString, DataProvidersManager, TimeScale), arg0, arg1, arg2)
end

function find_satellite_antenna(obj::AntexLoader, arg0::SatelliteSystem, arg1::jint)
    return jcall(obj, "findSatelliteAntenna", TimeSpanMap, (SatelliteSystem, jint), arg0, arg1)
end

function get_receivers_antennas(obj::AntexLoader)
    return jcall(obj, "getReceiversAntennas", List, ())
end

function get_satellites_antennas(obj::AntexLoader)
    return jcall(obj, "getSatellitesAntennas", List, ())
end

ext: jl


text:
using Test
using MIPVerify: setloglevel!
using MIPVerify: remove_cached_models
using MIPVerify: get_max_index, get_norm
using JuMP
using TimerOutputs

@isdefined(TestHelpers) || include("TestHelpers.jl")

macro timed_testset(name::String, block)
    # copied from https://github.com/KristofferC/Tensors.jl/blob/master/test/runtests.jl#L8
    return quote
        @timeit "$($(esc(name)))" begin
            @testset "$($(esc(name)))" begin
                $(esc(block))
            end
        end
    end
end

@testset "MIPVerify" begin
    reset_timer!()
    setloglevel!("info")
    remove_cached_models()

    include("integration.jl")
    include("net_components.jl")
    include("utils.jl")
    include("models.jl")
    include("batch_processing_helpers.jl")

    @testset "get_max_index" begin
        @test_throws MethodError get_max_index([])
        @test get_max_index([3]) == 1
        @test get_max_index([3, 1, 4]) == 3
        @test get_max_index([3, 1, 4, 1, 5, 9, 2]) == 6
    end

    @testset "get_norm" begin
        @testset "real-valued arrays" begin
            xs = [1, -2, 3]
            @test get_norm(1, xs) == 6
            @test get_norm(2, xs) == sqrt(14)
            @test get_norm(Inf, xs) == 3
            @test_throws DomainError get_norm(3, xs)
        end
        @testset "variable-valued arrays" begin
            @testset "l1" begin
                m = TestHelpers.get_new_model()
                x1 = @variable(m, lowerbound=1, upperbound=5)
                x2 = @variable(m, lowerbound=-8, upperbound=-2)
                x3 = @variable(m, lowerbound=3, upperbound=10)
                xs = [x1, x2, x3]
                n_1 = get_norm(1, xs)
                n_2 = get_norm(2, xs)
                n_inf = get_norm(Inf, xs)

                @objective(m, Min, n_1)
                solve(m)
                @test getobjectivevalue(m)≈6

                if Base.find_package("Gurobi") != nothing
                    # Skip these tests if Gurobi is not installed.
                    # Cbc does not solve problems with quadratic objectives
                    @objective(m, Min, n_2)
                    solve(m)
                    @test getobjectivevalue(m)≈14
                end

                @objective(m, Min, n_inf)
                solve(m)
                @test getobjectivevalue(m)≈3

                @test_throws DomainError get_norm(3, xs)
            end
        end
    end

    println()
    print_timer()
    println()
    println()
end


text:
function OC_doShortwaveRadiation!(
    ocn    :: Ocean,
    i      :: Integer,
    j      :: Integer;
    Tswflx :: Float64,
    Δt     :: Float64,
)
    ocn.T_ML[i, j] = doShortwaveRadiation!(
        Tswflx = Tswflx,
        Ts   = ocn.cols.Ts[i, j],
        zs   = ocn.cols.zs[i, j],
        hs   = ocn.cols.hs[i, j],
        rad_decay_coes  = ocn.cols.rad_decay_coes[i, j],
        rad_absorp_coes = ocn.cols.rad_absorp_coes[i, j],
        T_ML = ocn.T_ML[i, j],
        h_ML = ocn.h_ML[i, j],
        Nz   = ocn.Nz[i, j],
        FLDO = ocn.FLDO[i, j],
        ζ    = ocn.ζ,
        Δt   = Δt,
    )

end

function doShortwaveRadiation!(;
    Tswflx          :: Float64,
    Ts              :: AbstractArray{Float64, 1},
    zs              :: AbstractArray{Float64, 1},
    hs              :: AbstractArray{Float64, 1},
    rad_decay_coes  :: AbstractArray{Float64, 1},
    rad_absorp_coes :: AbstractArray{Float64, 1},
    T_ML            :: Float64,
    h_ML            :: Float64,
    Nz              :: Integer,
    FLDO            :: Integer,
    ζ               :: Float64,
    Δt              :: Float64,
)

    # ===== [BEGIN] Mixed layer =====

    if FLDO == -1      # Entire ocean column is mixed-layer
        T_ML += - Tswflx * Δt / h_ML
        Ts[1:Nz] .= T_ML
        return T_ML
    end

    rad_decay_coes_ML = exp(-h_ML/ζ)
    T_ML += - Tswflx * (1.0 - rad_decay_coes_ML) * Δt / h_ML

    if FLDO > 1
        Ts[1:FLDO-1] .= T_ML
    elseif FLDO == -1
        Ts[1:Nz] .= T_ML
    end

    # ===== [END] Mixed layer =====

    # ===== [BEGIN] FLDO layer =====

    h_FLDO = - h_ML - zs[FLDO+1]

    if FLDO == Nz  # FLDO is last layer
        Ts[FLDO] += - Tswflx * rad_decay_coes_ML * Δt / h_FLDO
        return T_ML
    else
        Ts[FLDO] += - Tswflx * (rad_decay_coes_ML - rad_decay_coes[FLDO+1]) * Δt / h_FLDO
    end
    # ===== [END] FLDO layer =====

    # ===== [BEGIN] Rest layers =====

    for k=FLDO+1:Nz
        Ts[k] += - Tswflx * rad_decay_coes[k] * rad_absorp_coes[k] * Δt / hs[k]
    end

    # ===== [END] Rest layers =====

    return T_ML
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
using EmpiricalCDFs

using EmpiricalCDFs.IOcdf: readcdf, getcdf

using Test
using Statistics

@testset "construction" begin
    @test (cdf = EmpiricalCDF(); true)
    @test (cdf = EmpiricalCDF(); append!(cdf,rand(10^3)) ; true)
    @test (cdf = EmpiricalCDF(); append!(cdf,rand(10^3)) ; sort!(cdf); cdf(0) isa Float64)
    @test (cdf = EmpiricalCDF(); append!(cdf,rand(10^3)) ; sort!(cdf); typeof(rand(cdf)) == Float64)
    @test (cdf = EmpiricalCDF(.5); append!(cdf,rand(10^3)) ; sort!(cdf); typeof(cdf(1)) == Float64)
end

@testset "copy, equality" begin
    cdf = EmpiricalCDF()
    append!(cdf, rand(10^3))
    cdf2 = copy(cdf)
    @test cdf == cdf2

    cdfh = EmpiricalCDFHi(0.0)
    append!(cdfh, randn(10^4))
    cdfh2 = EmpiricalCDFHi(0.0)
    append!(cdfh2, cdfh)
    @test cdfh2 == cdfh
end

@testset "cutoff" begin
    a = rand(10^5);            # points in [0,1)
    cdf0 = EmpiricalCDF();     # accept all points
    cdf1 = EmpiricalCDF(0.5);  # reject points < 0.5

    append!(cdf0,a);
    append!(cdf1,a);

    @test maximum(cdf0) == maximum(cdf1)
    @test minimum(cdf0) != minimum(cdf1)
    @test minimum(cdf1) >= 0.5

    sort!(cdf0)
    sort!(cdf1)

    @test cdf0(.4) != cdf1(.4)
    @test cdf0(.6) == cdf1(.6)

    @test finv(cdf0)(1-eps(1.0)) == maximum(cdf0)   # last point is as close to 1 as possible
    @test finv(cdf1)(1-eps(1.0)) == maximum(cdf1)
end

@testset "read binary file" begin
    infile = joinpath(@__DIR__, "paretocdf.bin")
    paretocdf = readcdf(infile)
    @test length(paretocdf) == 10^5
    @test length(paretocdf.([.4,.5,.6])) == 3
    @test getcdf(paretocdf).([.4,.5,.6]) == getcdf(paretocdf)([.4,.5,.6])
    @test paretocdf.([.4,.5,.6]) == paretocdf([.4,.5,.6])
    @test (linprint(devnull, paretocdf) ; true)
    @test (logprint(devnull, paretocdf); true)
    @test rand(paretocdf) isa AbstractFloat
    @test (push!(paretocdf, .5) ; true)
    @test (show(devnull, paretocdf); true)
    p = getcdf(paretocdf)
    @test data(p) == p.xdata
    @test counts(p) == length(p)  # NOTE: counts and length are not always the same, but are in this case.
end

@testset "Inf cutoff" begin
    c = EmpiricalCDF(Inf)
    c2 = EmpiricalCDF(Inf)
    @test typeof(c) == typeof(c2)
    @test length(c) == length(c2)
    c = EmpiricalCDFHi(Inf)
    @test (push!(c,.5);true)
    @test (append!(c,[1.,1e10]);true)
end

@testset "statistics" begin
    n = 10^5
    cdf = EmpiricalCDF(randn(n))
    @test length(cdf) == n
    # named parameters
    @test std(cdf; mean=mean(cdf)) == std(cdf)
    empty!(cdf)
    @test length(cdf) == 0
end


text:
using YaoBase, YaoArrayRegister
import StaticArrays: SMatrix
export RotationGate, Rx, Ry, Rz, rot

"""
    RotationGate{D,T,GT<:AbstractBlock{D}} <: PrimitiveBlock{D}

RotationGate, with GT both hermitian and isreflexive.

# Definition

```math
\\mathbf{I} cos(θ / 2) - im sin(θ / 2) * mat(U)
```
"""
mutable struct RotationGate{D,T,GT<:AbstractBlock{D}} <: PrimitiveBlock{D}
    block::GT
    theta::T
    function RotationGate{D,T,GT}(block::GT, theta) where {D,T,GT<:AbstractBlock{D}}
        ishermitian(block) && isreflexive(block) ||
            throw(ArgumentError("Gate type $GT is not hermitian or not isreflexive."))
        new{D,T,GT}(block, T(theta))
    end
end

RotationGate(block::GT, theta::T) where {D,T,GT<:AbstractBlock{D}} =
    RotationGate{D,T,GT}(block, theta)
RotationGate(block::GT, theta::Integer) where {D,GT<:AbstractBlock{D}} =
    RotationGate{D,Float64,GT}(block, Float64(theta))
nqudits(rg::RotationGate) = nqudits(rg.block)

# bindings
"""
    Rx(theta)

Return a [`RotationGate`](@ref) on X axis.

# Example

```jldoctest; setup=:(using YaoBlocks)
julia> Rx(0.1)
rot(X, 0.1)
```
"""
Rx(theta) = RotationGate(X, theta)

"""
    Ry(theta)

Return a [`RotationGate`](@ref) on Y axis.

# Example

```jldoctest; setup=:(using YaoBlocks)
julia> Ry(0.1)
rot(Y, 0.1)
```
"""
Ry(theta) = RotationGate(Y, theta)

"""
    Rz(theta)

Return a [`RotationGate`](@ref) on Z axis.

# Example

```jldoctest; setup=:(using YaoBlocks)
julia> Rz(0.1)
rot(Z, 0.1)
```
"""
Rz(theta) = RotationGate(Z, theta)

"""
    rot(U, theta)

Return a [`RotationGate`](@ref) on U axis.
"""
rot(axis::AbstractBlock, theta) = RotationGate(axis, theta)

content(x::RotationGate) = x.block
chcontent(x::RotationGate, gen::AbstractBlock) = RotationGate(gen, x.theta)
# General definition
function mat(::Type{T}, R::RotationGate{D}) where {D,T}
    I = IMatrix{D^nqudits(R),T}()
    return I * cos(T(R.theta) / 2) - im * sin(T(R.theta) / 2) * mat(T, R.block)
end

# Specialized
mat(::Type{T}, R::RotationGate{1,<:Any,<:XGate}) where {T} =
    T[cos(R.theta / 2) -im*sin(R.theta / 2); -im*sin(R.theta / 2) cos(R.theta / 2)]
mat(::Type{T}, R::RotationGate{1,<:Any,<:YGate}) where {T} =
    T[cos(R.theta / 2) -sin(R.theta / 2); sin(R.theta / 2) cos(R.theta / 2)]
# mat(R::RotationGate{1, T, ZGate{Complex{T}}}) where T =
#     SMatrix{2, 2, Complex{T}}(cos(R.theta/2)-im*sin(R.theta/2), 0, 0, cos(R.theta/2)+im*sin(R.theta/2))

function _apply!(r::ArrayReg, rb::RotationGate)
    v0 = copy(r.state)
    _apply!(r, rb.block)
    # NOTE: we should not change register's memory address,
    # or batch operations may fail
    r.state .= -im * sin(rb.theta / 2) * r.state + cos(rb.theta / 2) * v0
    return r
end

# parametric interface
niparams(::Type{<:RotationGate}) = 1
getiparams(x::RotationGate) = x.theta
# no need to specify the type of param, Julia will try to do the conversion
setiparams!(r::RotationGate, param::Number) = (r.theta = param; r)
setiparams(r::RotationGate, param::Number) = RotationGate(r.block, param)

# fallback to matrix methods if it is not real
YaoBase.isunitary(r::RotationGate{D,<:Real}) where {D} = true

function YaoBase.isunitary(r::RotationGate)
    isreal(r.theta) && return true
    @warn "θ in RotationGate is not real, got θ=$(r.theta), fallback to matrix-based method"
    return isunitary(mat(r))
end

Base.adjoint(blk::RotationGate) = RotationGate(blk.block, -blk.theta)
Base.copy(R::RotationGate) = RotationGate(R.block, R.theta)
Base.:(==)(lhs::RotationGate{TA,GTA}, rhs::RotationGate{TB,GTB}) where {TA,TB,GTA,GTB} =
    false
Base.:(==)(lhs::RotationGate{TA,GT}, rhs::RotationGate{TB,GT}) where {TA,TB,GT} =
    lhs.theta == rhs.theta

cache_key(R::RotationGate) = R.theta

iparams_range(::RotationGate{D,T,GT}) where {D,T,GT} = ((zero(T), T(2 * pi)),)

occupied_locs(g::RotationGate) = occupied_locs(g.block)


text:
module Option

using Configurations

@option struct Plugins
    template::Vector{String} = String[]
    registry::Vector{String} = String[]
end

@option struct Ion
    username::Maybe{String}
    default_template::String = "basic"
    templates::Dict{String, String} = Dict{String, String}(
        "basic" => templates("basic.toml"),
        "package" => templates("package.toml"),
        "academic" => templates("academic.toml"),
        "comonicon" => templates("comonicon.toml"),
        "comonicon-sysimg" => templates("comonicon-sysimg.toml"),
    )
end

end


text:
export add_cad_names_to_physical_groups

"""
    add_cad_names_to_physical_groups(dim::Int64=-1)

Add the names of imported CAD entities of dimension `dim` to the model's
physical groups. For an entity with name "Shapes/X/Y/Z", the physical group it
will be assigned is "Z". If `dim` == -1, add the names of all entities.
"""
function add_cad_names_to_physical_groups(dim::Int64=-1)
    # Get all entities. If they have a name, add them to the appropriate
    # entry in the physicals dict
    physicals = Dict{String, Vector{NTuple{2,Int32}}}()
    for ent in gmsh.model.get_entities(dim)
        edim, tag = ent
        name = gmsh.model.get_entity_name(edim, tag)
        if name != ""
            path = split(name, "/")
            ent_name = path[end]
            if ent_name ∉ keys(physicals)
                physicals[ent_name] = NTuple{2,Int32}[]
            end
            push!(physicals[ent_name], ent)
        end
    end
    # For each of the groups, add them to the model
    for (name, ents) in pairs(physicals)
        safe_add_physical_group(name, ents)
    end
    return sort!(collect(keys(physicals)))
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
"""
```
brief_params = BRIEF([size = 128], [window = 9], [sigma = 2 ^ 0.5], [sampling_type = gaussian], [seed = 123])
```

| Argument | Type | Description |
|----------|------|-------------|
| **size** | Int | Size of the descriptor |
| **window** | Int | Size of sampling window |
| **sigma** | Float64 | Value of sigma used for inital gaussian smoothing of image |
| **sampling_type** | Function | Type of sampling used for building the descriptor (See [BRIEF Sampling Patterns](#brief-sampling-patterns)) |
| **seed** | Int | Random seed used for generating the sampling pairs. For matching two descriptors, the seed used to build both should be same. |

"""
mutable struct BRIEF{F} <: Params
    size::Int
    window::Int
    sigma::Float64
    sampling_type::F
    seed::Int
end

function BRIEF(; size::Integer = 128, window::Integer = 9, sigma::Float64 = 2 ^ 0.5, sampling_type::Function = gaussian, seed::Int = 123)
    BRIEF(size, window, sigma, sampling_type, seed)
end

"""
```
sample_one, sample_two = random_uniform(size, window, seed)
```

Builds sampling pairs using random uniform sampling.
"""
function random_uniform(size::Int, window::Int, seed::Int)
    seed!(seed)
    sample_one = CartesianIndex{2}[]
    sample_two = CartesianIndex{2}[]
    count = 0
    while true
        x_gen = floor(Int, window * rand() / 2)
        y_gen = floor(Int, window * rand() / 2)
        (x_gen >= ceil(-window / 2) && x_gen <= floor((window - 1) / 2)) && (y_gen >= ceil(-window / 2) && y_gen <= floor((window - 1) / 2)) || continue
        x_gen_2 = floor(Int, window * rand() / 2)
        y_gen_2 = floor(Int, window * rand() / 2)
        (x_gen_2 >= ceil(-window / 2) && x_gen_2 <= floor((window - 1) / 2)) && (y_gen_2 >= ceil(-window / 2) && y_gen_2 <= floor((window - 1) / 2)) || continue
        count += 1
        push!(sample_one, CartesianIndex(y_gen, x_gen))
        push!(sample_two, CartesianIndex(y_gen_2, x_gen_2))
        count != size || break
    end
    sample_one, sample_two
end

"""
```
sample_one, sample_two = random_coarse(size, window, seed)
```

Builds sampling pairs using random sampling over a coarse grid.
"""
function random_coarse(size::Int, window::Int, seed::Int)
    seed!(seed)
    gen = rand(ceil(Int, ceil(-window / 2)): floor(Int, (window - 1) / 2), size * 4)
    sample_one = CartesianIndex{2}[]
    sample_two = CartesianIndex{2}[]
    for i in 1:size
        push!(sample_one, CartesianIndex(gen[i], gen[2 * i]))
        push!(sample_two, CartesianIndex(gen[3 * i], gen[4 * i]))
    end
    sample_one, sample_two
end

"""
```
sample_one, sample_two = gaussian(size, window, seed)
```

Builds sampling pairs using gaussian sampling.
"""
function gaussian(size::Int, window::Int, seed::Int)
    seed!(seed)
    _gaussian(size, window), _gaussian(size, window)
end

function _gaussian(size::Int, window::Int)
    set = Normal(0, (window ^ 2) / 25)
    count = 0
    sample = CartesianIndex{2}[]
    while true
        x_gen , y_gen = rand(set, 2)
        (x_gen >= ceil(-window / 2) && x_gen <= floor((window - 1) / 2)) && (y_gen >= ceil(-window / 2) && y_gen <= floor((window - 1) / 2)) || continue
        count += 1
        push!(sample, CartesianIndex(floor(Int, y_gen), floor(Int, x_gen)))
        count != size || break
    end
    sample
end

"""
```
sample_one, sample_two = gaussian_local(size, window, seed)
```

Pairs `(Xi, Yi)` are randomly sampled using a Gaussian distribution where first `X` is sampled with a standard deviation of `0.04*S^2` and
then the `Yi’s` are sampled using a Gaussian distribution – Each `Yi` is sampled with mean `Xi` and standard deviation of `0.01 * S^2`
"""
function gaussian_local(size::Int, window::Int, seed::Int)
    seed!(seed)
    _gaussian_local(size, window), _gaussian_local(size, window)
end

function _gaussian_local(size::Int, window::Int)
    x_set = Normal(0, (window ^ 2) / 25)
    count = 0
    sample = CartesianIndex{2}[]
    while true
        x_gen = floor(Int, rand(x_set))
        y_set = Normal(x_gen, (window ^ 2) / 100)
        y_gen = floor(Int, rand(y_set))
        (x_gen >= ceil(-window / 2) && x_gen <= floor((window - 1) / 2)) && (y_gen >= ceil(-window / 2) && y_gen <= floor((window - 1) / 2)) || continue
        count += 1
        push!(sample, CartesianIndex(y_gen, x_gen))
        count != size || break
    end
    sample
end

"""
```
sample_one, sample_two = center_sample(size, window, seed)
```

Builds sampling pairs `(Xi, Yi)` where `Xi` is `(0, 0)` and `Yi` is sampled uniformly from the window.
"""
function center_sample(size::Int, window::Int, seed::Int)
    seed!(seed)
    count = 0
    sample = CartesianIndex{2}[]
    while true
        x_gen = floor(Int, window * rand() / 2)
        y_gen = floor(Int, window * rand() / 2)
        (x_gen >= ceil(-window / 2) && x_gen <= floor((window - 1) / 2)) && (y_gen >= ceil(-window / 2) && y_gen <= floor((window - 1) / 2)) || continue
        count += 1
        push!(sample, CartesianIndex(y_gen, x_gen))
        count != size || break
    end
    zeros(CartesianIndex{2}, size), sample
end

function create_descriptor(img::AbstractArray{T, 2}, keypoints::Keypoints, params::BRIEF) where T<:Gray
    factkernel = KernelFactors.IIRGaussian([params.sigma, params.sigma])
    img_smoothed = imfilter(Float64, img, factkernel, NA())
    sample_one, sample_two = params.sampling_type(params.size, params.window, params.seed)
    descriptors = BitArray{1}[]
    h, w = size(img_smoothed)
    ret_keypoints = Keypoint[]
    for k in keypoints
        (k[1] > floor(params.window / 2) && k[1] <= h - floor((params.window - 1) / 2)) && (k[2] > floor(params.window / 2) && k[2] <= w - floor((params.window - 1) / 2)) || continue
        temp = BitArray{1}([])
        for (s1, s2) in zip(sample_one, sample_two)
            push!(temp, img_smoothed[k + s1] < img_smoothed[k + s2])
        end
        push!(descriptors, temp)
        push!(ret_keypoints, k)
    end
    descriptors, ret_keypoints
end


text:
module MyExample

# Write your package code here.
include("extra_file.jl")
export my_f
end


text:
using Dolo

# model = yaml_import("examples/models/rbc.yaml")

model = yaml_import("experiments/ayiagari.yaml")

Dolo.get_domain(model.exogenous)

gr = Dolo.discretize(model)

sol = improved_time_iteration(model)

# sim = simulate(model, sol.dr)

tab = tabulate(model, sol.dr, :a)

using SimplePlots

plot(tab[:a], tab[:c])


text:
using MatrixCalculus
model = Dynamics.pendulum_uncertain
n = model.n; m = model.m; r = model.r

eig_thr = 1.0e-3

u_max = 3.
u_min = -3.

h_max = Inf
h_min = 0.0

# Problem
x0 = [0.;0.]
xf = [pi;0.]

E0 = Diagonal(1.0e-6*ones(n))
H0 = zeros(n,r)
D = Diagonal([.2^2])

N = 51 # knot points

Q = Diagonal(zeros(n))
R = Diagonal(zeros(m))
Rh = 1.0
Qf = Diagonal(zeros(n))

Q_lqr = Diagonal([10.;1.])
R_lqr = Diagonal(0.1*ones(m))
Qf_lqr = Diagonal([100.;100.])

Q_r = Q_lqr
R_r = R_lqr
Qf_r = Qf_lqr

tf0 = 2.
dt = tf0/(N-1)

# problem
obj = LQRObjective(Q,R,Qf,xf,N)

goal_con = goal_constraint(xf)
bnd_con = BoundConstraint(n,m,u_min=u_min,u_max=u_max,trim=true)
con = Constraints([bnd_con],N)

prob = TrajectoryOptimization.Problem(model, obj,constraints=con, N=N, tf=tf0, x0=x0, dt=dt)
prob.constraints[N] += goal_con # add goal constraint
copyto!(prob.X,line_trajectory(x0,xf,N)) # initialize state with linear interpolation
copyto!(prob.U,[0.01*rand(m) for k = 1:N-1]) # initialize control with rand

solver = DIRTRELSolver(E0,H0,D,Q,R,Rh,Qf,xf,Q_lqr,R_lqr,Qf_lqr,Q_r,R_r,Qf_r,eig_thr,h_max,h_min)

solve!(prob,solver)

plot(prob.X)
plot(prob.U)

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


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

using Markdown
using InteractiveUtils

# ╔═╡ 6fd616f6-1ea1-11eb-3814-8bfb4a096c49
using Pkg, DrWatson

# ╔═╡ 6ffe2628-1ea1-11eb-24ea-57f985146a72
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Distributions, Random, StatsBase, DataFrames, Plots
	Random.seed!(1)
	include(datadir("mvParams.jl"))
end;

# ╔═╡ ed174bc4-1ea0-11eb-1e2f-a32874cec549
md"## Listing 3.34"

# ╔═╡ a64855e4-2e89-11eb-11e7-79624219151b
meanVect

# ╔═╡ b14a292c-2e89-11eb-3cd0-2bed0aa6037d
covMat

# ╔═╡ 70179b62-1ea1-11eb-13ce-7fb4d54bad4b
biNorm = MvNormal(meanVect, covMat)

# ╔═╡ 52aba0a8-2e89-11eb-2638-a5e8aec2d22e
rand(biNorm, 10)

# ╔═╡ a71a7f10-2290-11eb-3cc0-3b6d0e5fa6d0
begin
	N = 10^3
	points = rand(MvNormal(meanVect,covMat),N)

	support = 15:0.5:40
	z = [ pdf(biNorm,[x,y]) for y in support, x in support ]
end;

# ╔═╡ a71abdcc-2290-11eb-0683-5d9b8f53e51d
begin
	p1 = scatter(points[1,:], points[2,:], ms=0.5, c=:black, legend=:none)
	p1 = contour!(support, support, z,
			levels=[0.001, 0.005, 0.02], c=[:blue, :red, :green],
			xlims=(15,40), ylims=(15,40), ratio=:equal, legend=:none,
			xlabel="x", ylabel="y")
end;

# ╔═╡ c92aa622-2293-11eb-3afe-3da4f7d14f3d
begin
	p2 = surface(support, support, z, lw=0.2, #c=cgrad([:blue, :red]),
			 legend=:none, xlabel="x", ylabel="y")
	plot(p1, p2, size=(800,400))
end

# ╔═╡ 70182cd8-1ea1-11eb-094d-8d0d49cf15f3
md"## End of listing 3.34"

# ╔═╡ Cell order:
# ╟─ed174bc4-1ea0-11eb-1e2f-a32874cec549
# ╠═6fd616f6-1ea1-11eb-3814-8bfb4a096c49
# ╠═6ffe2628-1ea1-11eb-24ea-57f985146a72
# ╠═a64855e4-2e89-11eb-11e7-79624219151b
# ╠═b14a292c-2e89-11eb-3cd0-2bed0aa6037d
# ╠═70179b62-1ea1-11eb-13ce-7fb4d54bad4b
# ╠═52aba0a8-2e89-11eb-2638-a5e8aec2d22e
# ╠═a71a7f10-2290-11eb-3cc0-3b6d0e5fa6d0
# ╠═a71abdcc-2290-11eb-0683-5d9b8f53e51d
# ╠═c92aa622-2293-11eb-3afe-3da4f7d14f3d
# ╟─70182cd8-1ea1-11eb-094d-8d0d49cf15f3


text:
"""
ArithmeticBlocks
"""

export TrigonometricFunctionBlock

mutable struct TrigonometricFunctionBlock <: AbstractArithmeticBlock
    inport::Vector{InPort}
    outport::Vector{OutPort}
    op::Symbol

    function TrigonometricFunctionBlock(op::Symbol)
        @createblock new(Vector{InPort}(), Vector{OutPort}(), op) 1 1
    end
end

"""
IO
"""

function Base.show(io::IO, blk::TrigonometricFunctionBlock)
    println(io, "TrigonometricFunctionBlock($(blk.op))")
end

"""
to expr
"""

function _toexpr(blk::TrigonometricFunctionBlock)
        Expr(:call, blk.op, _toexpr(blk.inport[1]))
end


text:
"""
    seg2 = merge_segments(seg, threshold)

Merges segments in a [`SegmentedImage`](@ref) by building a region adjacency
graph (RAG) and merging segments connected by edges with weight less than
`threshold`.

# Arguments:
* `seg`         : SegmentedImage to be merged.
* `threshold`   : Upper bound of the adjacent segment color difference to
                  consider merging segments.

# Citation:
Vighnesh Birodkar
"Hierarchical merging of region adjacency graphs"
https://vcansimplify.wordpress.com/2014/08/17/hierarchical-merging-of-region-adjacency-graphs/
"""
function merge_segments(seg::SegmentedImage, threshold::Number)::SegmentedImage
    g = seg_to_graph(seg)

    # Populate a heap of all the edges, and a Bool indicating whether the edge
    # is valid. All edges are initially valid. The reason for this is that heap
    # removal would be expensive, so instead, we invalidate the edge entry in the
    # heap.
    function weight(t::Tuple{Edge{Int}, Bool})::Real
        return has_prop(g, t[1], :weight) ?  get_prop(g, t[1], :weight) : 0
    end

    edge_heap = MutableBinaryHeap{Tuple{Edge{Int}, Bool}}(Base.By(weight),
        [(e, true) for e in edges(g)]
    )
    sizehint!(edge_heap, 3 * length(edge_heap))  # Overkill, or not enough?
    for n in edge_heap.nodes
        set_prop!(g, n.value[1], :handle, n.handle)
    end

    # Merge all edges less than threshold
    while !isempty(edge_heap) && weight(first(edge_heap)) < threshold
        e, valid = pop!(edge_heap)
        if valid
            # Invalidate all edges touching this edge.
            invalidate_neighbors!(edge_heap, g, e)

            # Merge the two nodes into one (keep e.dst, obsolete e.src)
            merge_node_props!(g, e)

            # Make new edges to the merged node.
            new_edges = add_neighboring_edges!(g, e)

            # Remove edges to src.
            # Don't call rem_vertex!(g, e.src); it would renumber all vertices.
            for n in collect(neighbors(g, e.src))
                rem_edge!(g, e.src, n)
            end

            # Add new edges to heap.
            for e in new_edges
                handle = push!(edge_heap, (e, true))
                set_prop!(g, e, :handle, handle)
            end
        end
    end

    return resegment(seg, g)
end

"""
    g = seg_to_graph(seg)

Given a [`SegmentedImage`](@ref), produces a region adjacency [`MetaGraph`](@ref)
and stores segment metadata on the vertices. Edge weight is determined by
color difference.

# Arguments:
* `seg`         : a [`SegmentedImage`](@ref)
"""
function seg_to_graph(seg::SegmentedImage)::MetaGraph
    weight(i, j) = colordiff(segment_mean(seg, i), segment_mean(seg, j))
    rag, _ = region_adjacency_graph(seg, weight)

    g = MetaGraph(rag)
    for v in vertices(rag)
        set_prop!(g, v, :labels, [v])
        set_prop!(g, v, :pixel_count, seg.segment_pixel_count[v])
        set_prop!(g, v, :mean_color, seg.segment_means[v])
        set_prop!(g, v, :total_color, seg.segment_means[v] * seg.segment_pixel_count[v])
    end

    for e in edges(rag)
        set_prop!(g, Edge(e.src, e.dst), :weight, e.weight)
    end
    return g
end

"""
    seg2 = resegment(seg1, rag)

Takes a segmentation and a region adjacency graph produced by `merge_segments`
and produces an segmentation that corresponds to the graph.

# Arguments:
* `seg`         : a [`SegmentedImage`](@ref)
* `g`           : a [`MetaGraph`](@ref) representing a merged Region Adjacency
                  Graph
"""
function resegment(seg::SegmentedImage, g::MetaGraph)::SegmentedImage
    # Find all the vertices of g that remain post-merge.
    remaining = collect(filter(v -> 0 < length(props(g, v)), vertices(g)))

    px_labels = copy(seg.image_indexmap)
    # Re-label all pixels with the vertex they were merged to.
    for v in remaining
        labels = get_prop(g, v, :labels)
        for l in labels
            if l != v
                ix = findall(x -> x == l, px_labels)
                px_labels[ix] .= v
            end
        end
    end

    # Re-number our labels so that they are dense (no gaps) and
    # construct the other objects SegmentedImage needs.
    means, px_counts = Dict{Int, Colorant}(), Dict{Int, Int}()
    for (i, v) in enumerate(remaining)
        ix = findall(x -> x == v, px_labels)
        px_labels[ix] .= i
        means[i] = get_prop(g, v, :mean_color)
        px_counts[i] = get_prop(g, v, :pixel_count)
    end

    labels = collect(1:length(remaining))

    return SegmentedImage(px_labels, labels, means, px_counts)
end

"""
    merge_node_props!(g, e)

Takes edge `e` in [`MetaGraph`](@ref) `g` and merges the props from its `src`
and `dst` into its `dst`, clearing all props from `src`.

"""
function merge_node_props!(g::MetaGraph, e::AbstractEdge)
    src, dst = e.src, e.dst
    clr = get_prop(g, dst, :total_color) + get_prop(g, src, :total_color)
    npx = get_prop(g, dst, :pixel_count) + get_prop(g, src, :pixel_count)

    set_prop!(g, dst, :total_color, clr)
    set_prop!(g, dst, :pixel_count, npx)
    set_prop!(g, dst, :mean_color, clr / npx)
    set_prop!(g, dst, :labels, vcat(
         get_prop(g, src, :labels),
         get_prop(g, dst, :labels)
    ))

    # Clear props on the now unused node src, to make its obsolescence clear.
    clear_props!(g, src)
end

"""
    new_edges = add_neighboring_edges!(g, e)

Finds the nodes neighboring `e` in graph `g`, creates edges from them to its
`dst`, and  sets the weight of the new edges.

# Arguments:
* `g`         : a [`MetaGraph`](@ref)
* `e`         : an [`AbstractEdge`](@ref)
"""
function add_neighboring_edges!(g::MetaGraph, e::AbstractEdge)
    edges = Edge{eltype(e)}[]
    edge_neighbors = union(Set(neighbors(g, e.src)), Set(neighbors(g, e.dst)))
    for n in setdiff(edge_neighbors, e.src, e.dst)
        edge = Edge(e.dst, n)
        add_edge!(g, edge)
        set_prop!(g, edge, :weight, _weight_mean_color(g, edge))
        push!(edges, edge)
    end

    return edges
end

"""
    invalidate_neighbors!(edge_heap, g, e)

Finds the neighbors of `e` in graph `g` and invalidates them in `edge_heap`.

# Arguments:
* `edge_heap` : a [`MutableBinaryHeap`](@ref)
* `g`         : a [`MetaGraph`](@ref)
* `e`         : an [`AbstractEdge`](@ref)
"""
function invalidate_neighbors!(edge_heap::MutableBinaryHeap, g::MetaGraph, e::AbstractEdge)
    function invalidate(src, dst)
        for n in setdiff(Set(neighbors(g, src)), dst)
            edge = Edge(src, n)
            h = get_prop(g, edge, :handle)
            update!(edge_heap, h, (edge, false))
        end
    end
    invalidate(e.src, e.dst)
    invalidate(e.dst, e.src)
end

"""
    weight = _weight_mean_color(g, v1, v2))

Compute the weight of an edge in [`MetaGraph`](@ref) `g` as the difference
in mean colors of each vertex.

# Arguments:
* `g`         : a [`MetaGraph`](@ref)
* `e`         : an [`AbstractEdge`](@ref)
"""
function _weight_mean_color(g::MetaGraph, e::AbstractEdge)::Real
    return colordiff(
        get_prop(g, e.src, :mean_color),
        get_prop(g, e.dst, :mean_color)
    )
end


text:
"""
The `out` that a layer returns in `forward!` is allowed to be mutated, and thus can't be relied upon.
However, this does allow it to mutate the inputs.

`reverse_grad!` is called before `reverse_chain!`
"""
abstract type AbstractLayer end

struct Chain{D,C<:Tuple{Vararg{<:AbstractLayer,D}},O,F}
    chain::C
    optstates::NTuple{D,O}
    loss::F
end

@generated function (c::Chain{D})(x_0) where {D}
    quote
        Base.Cartesian.@nexprs $D d -> c_d = c.chain[d]
        Base.Cartesian.@nexprs $D d -> (x_d, xpb_d) = forward!(c_d, x_{d-1})
        loss(c.loss, $(Symbol(:x_, D)))
    end
end
@generated function gradient!(c::Chain{D}, x_0) where {D}
    x_D = Symbol(:x_,D);
    quote
        Base.Cartesian.@nexprs $D d -> c_d = c.chain[d]
        Base.Cartesian.@nexprs $D d -> (x_d, xpb_d) = forward!(c_d, x_{d-1})
        gradient!($x_D, c.loss, $x_D)
        Base.Cartesian.@nexprs $(D-1) d -> begin
            reverse_grad!(c_{$(D+1}-d}, xpb_{$(D+1}-d}, x_{$(D+1}-d})
            reverse_chain!(c_{$(D+1}-d}, xpb_{$(D+1}-d}, x_{$(D+1}-d})
        end
        reverse_grad!(c_1, xpb_1, x_1)
        (Base.Cartesian.@ntuple $D d -> grad(c_d))
    end
end


text:
# ============= General Utility functions ==================

ninvariants(D::NBodyDescriptor, N::Integer) = ninvariants(D, Val(N))
ninvariants(D::NBodyDescriptor, vN::Val{N}) where {N} = length.(tdegrees(D, vN))

@inline transform(D::NBodyDescriptor, r::Number) = transform(D.transform, r)
@inline transform_d(D::NBodyDescriptor, r::Number) = transform_d(D.transform, r)
@inline cutoff(D::NBodyDescriptor) = cutoff(D.cutoff)

hash(::BASIS, D::NBodyDescriptor) =
      hash((hash(BASIS(), D.transform), hash(BASIS(), D.cutoff)))

evaluate(V::NBodyFunction, Rs, i, J) =
      evaluate(V, descriptor(V), Rs, i, J)

evaluate_d!(dVsite, V::NBodyFunction, Rs, i, J) =
      evaluate_d!(dVsite, V, descriptor(V), Rs, i, J)

evaluate_many!(out, B, Rs, i, J) =
      evaluate_many!(out, B, descriptor(B[1]), Rs, i, J)

evaluate_many_d!(out, B, Rs, i, J) =
      evaluate_many_d!(out, B, descriptor(B[1]), Rs, i, J)

# TODO: switch to tuples??
struct TempEdV{T1, T2}
   Etemp::T1
   dVsite::T2
end

evaluate_many_ed!(out::TempEdV, B, Rs, i, J) =
      evaluate_many_ed!(out, B, descriptor(B[1]), Rs, i, J)

"""
`_sdot(a::T, B::SVector{N, T})`: efficiently compute `[ a .* b for b in B ]`
"""
@generated function _sdot(a::T, B::SVector{N, T}) where {N, T}
   code = "@SVector $T["
   for n = 1:N
      code *= "a .* B[$n],"
   end
   code = code[1:end-1] * "]"
   ex = Meta.parse(code)
   quote
      $ex
   end
end

# ------------- main evaluation code -----------

function evaluate(V::NBodyFunction{N},
                  desc::NBSiteDescriptor,
                  Rs::AbstractVector{JVec{T}},
                  i::Int,
                  J::SVector{K, Int}) where {N, T, K}
   # get the physical descriptor: bond-lengths (+ bond-angles)
   rθ = ricoords(desc, Rs, J)
   # check whether to skip this N-body term?
   skip_simplex(desc, rθ) && return zero(T)
   # compute the cutoff (and skip this site if the cutoff is zero)
   fc = fcut(desc, rθ)
   fc == 0 && return zero(T)
   # compute the invariants (this also applies the transform)
   II = invariants(desc, rθ)
   # evaluate the inner potential function (e.g. polynomial)
   return evaluate_I(V, II) * fc
end

# (out, s, S, J, _) -> _grad_len2pos!(out, evaluate_d(V, s)/N, J, S),
function evaluate_d!(dVsite,
                     V::NBodyFunction{N},
                     desc::NBSiteDescriptor,
                     Rs,
                     i::Int,
                     J) where {N}
   # get the physical descriptor: bond-lengths (+ bond-angles)
   rθ = ricoords(desc, Rs, J)  # 115 / 145
   # check whether to skip this N-body term?
   skip_simplex(desc, rθ) && return dVsite   # 12 / 122
   # compute the cutoff (and skip this site if the cutoff is zero)
   fc, fc_d = fcut_d(desc, rθ)   # 111 / 134
   fc == 0 && return dVsite
   # get the invariants
   II = invariants_ed(desc, rθ)   # 47 / 39
   # evaluate the inner potential function (e.g. polynomial)
   V, dV_drθ = evaluate_I_ed(V, II)  # 258 / 268
   # convert to gradient w.r.t. rθ
   dV_drθ = dV_drθ * fc + V * fc_d     # 4
   # convert to gradient w.r.t. dR
   return gradri2gradR!(desc, dVsite, dV_drθ, Rs, J, rθ)  # 18 / 15
end

function evaluate_many!(Es,
                        B::AbstractVector{TB},
                        desc::NBSiteDescriptor,
                        Rs, i, J)  where {TB <: NBodyFunction{N}} where {N}
   rθ = ricoords(desc, Rs, J)
   skip_simplex(desc, rθ) && return Es
   return _evaluate_many_ricoords!(Es, B, desc, rθ)
end

evaluate_many_ricoords!(Es,
               B::AbstractVector{TB},
               rθ) where {TB <: NBodyFunction{N}} where {N} =
   _evaluate_many_ricoords!(Es, B, descriptor(B[1]), rθ)

function _evaluate_many_ricoords!(Es,
               B::AbstractVector{TB}, desc,
               rθ) where {TB <: NBodyFunction{N}} where {N}
   fc = fcut(desc, rθ)
   fc == 0 && return Es
   II = invariants(desc, rθ)
   # evaluate the inner potential function (e.g. polynomial)
   for n = 1:length(B)
      Es[n] += evaluate_I(B[n], II) * fc
   end
   return Es
end

function evaluate_many_d!(dVsite::AbstractVector,
                          B::AbstractVector{TB},
                          desc::NBSiteDescriptor,
                          Rs,
                          i, J)  where {TB <: NBodyFunction{N}} where {N}
   rθ = ricoords(desc, Rs, J)
   skip_simplex(desc, rθ) && return dVsite
   fc, fc_d = fcut_d(desc, rθ)
   fc == 0 && return dVsite
   II = invariants_ed(desc, rθ)

   # evaluate the inner potential function (e.g. polynomial)
   for n = 1:length(B)
      V, dV_drθ = evaluate_I_ed(B[n], II)
      dV_drθ = dV_drθ * fc + V * fc_d
      gradri2gradR!(desc, dVsite[n], dV_drθ, Rs, J, rθ)
   end
   return dVsite
end

function evaluate_many_ed!(out::TempEdV,
                           B::AbstractVector{TB},
                           desc::NBSiteDescriptor,
                           Rs,
                           i, J)  where {TB <: NBodyFunction{N}} where {N}
   rθ = ricoords(desc, Rs, J)
   skip_simplex(desc, rθ) && return out
   fc, fc_d = fcut_d(desc, rθ)
   fc == 0 && return out
   II = invariants_ed(desc, rθ)

   # evaluate the inner potential function (e.g. polynomial)
   for n = 1:length(B)
      V, dV_drθ = evaluate_I_ed(B[n], II)
      out.Etemp[n] += V * fc
      dV_drθ = dV_drθ * fc + V * fc_d
      gradri2gradR!(desc, out.dVsite[n], dV_drθ, Rs, J, rθ)
   end
   return out
end

# bond length descriptor
include("bldescriptor.jl")

# bond angle descriptor
include("badescriptor.jl")

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
using Documenter, DocumenterCitations
using Mill, Flux, Random, SparseArrays, Setfield, HierarchicalUtils

#=
Useful resources for writing docs:
    Julia guidelines: https://docs.julialang.org/en/v1/manual/documentation/
    Documenter syntax: https://juliadocs.github.io/Documenter.jl/stable/man/syntax/
    Showcase: https://juliadocs.github.io/Documenter.jl/stable/showcase/
    Doctests: https://juliadocs.github.io/Documenter.jl/stable/man/doctests/
=#

DocMeta.setdocmeta!(Mill, :DocTestSetup, quote
    using Mill, Flux, Random, SparseArrays, Setfield, HierarchicalUtils
end; recursive=true)

makedocs(
         CitationBibliography(joinpath(@__DIR__, "references.bib")),
         sitename = "Mill.jl",
         format = Documenter.HTML(sidebar_sitename=false,
                                  collapselevel = 1,
                                  prettyurls=get(ENV, "CI", nothing) == "true",
                                  assets=["assets/favicon.ico", "assets/custom.css"]),
         modules = [Mill],
         pages = [
                  "Home" => "index.md",
                  "Motivation" => "motivation.md",
                  "Manual" => [
                               "Nodes" => "manual/nodes.md",
                               "More on nodes" => "manual/more_on_nodes.md",
                               "Model reflection" => "manual/reflectin.md",
                               "Bag aggregation" => "manual/aggregation.md",
                               "Data in leaves" => "manual/leaf_data.md",
                               "Missing data" => "manual/missing.md",
                               "Custom nodes" => "manual/custom.md"
                              ],
                  "Examples" => [
                                 "Musk" => "examples/musk/musk.md",
                                 "GNNs in 16 lines" => "examples/graphs.md",
                                 "DAGs" => "examples/dag.md",
                                 "Processing JSONs" => "examples/jsons.md"
                                ],
                  "External tools" => [
                                       "HierarchicalUtils.jl" => "tools/hierarchical.md"
                                      ],
                  "Public API" => [
                            "Aggregation" => "api/aggregation.md",
                            "Bags" => "api/bags.md",
                            "Data nodes" => "api/data_nodes.md",
                            "Model nodes" => "api/model_nodes.md",
                            "Special Arrays" => "api/special_arrays.md",
                            "Switches" => "api/switches.md",
                            "Utilities" => "api/utilities.md"
                           ],
                  "References" => "references.md",
                  "Citation" => "citation.md"
                 ],

        )

deploydocs(
           repo = "github.com/CTUAvastLab/Mill.jl.git",
          )


text:
module Example2018

export Point
struct Point{T}
  x::T
  y::T
end

function Base. +(p::Point, q::Point)
  Point(p.x+q.x, p.y+q.y)
end

export MPoint
mutable struct MPoint{T}
  x::T
  y::T
end

export Polynomial
struct Polynomial{T}
  coeffs::Vector{T}
end

export evaluate
function evaluate(p::Polynomial, x)
  sum(c*x^(i-1) for (i,c) in enumerate(p.coeffs))
end
#overwrites add operater with a new function for adding polynomials
function Base. +(p::Polynomial, q::Polynomial)
  l = max(length(p.coeffs), length(q.coeffs))
  rcoeffs = zeros(l)
  for (i,c) in enumerate(p.coeffs)
    rcoeffs[i] += c
  end
  for (i,c) in enumerate(q.coeffs)
    rcoeffs[i] += c
  end
  Polynomial(rcoeffs)
end

function Base. *(a, p::Polynomial)
  Polynomial(map(x->a*x, p.coeffs))
end
Base. *(p::Polynomial,a)=a*p

export deriv
function deriv(p::Polynomial)
  Polynomial([(i-1)*c for (i,c) in enumerate(p.coeffs)][2:end])
end

end # module


text:
# Acknowledgement: this adaption settings is mimicing Stan's 3-phase adaptation.

# Ref： https://github.com/stan-dev/math/blob/develop/stan/math/prim/mat/fun/welford_var_estimator.hpp
mutable struct VarEstimator{T<:Real}
  n :: Int
  μ :: Vector{T}
  M :: Vector{T}
end

reset!(ve::VarEstimator) = begin
  ve.n = 0
  ve.μ = zeros(size(ve.μ))
  ve.M = zeros(size(ve.M))
end

add_sample!(ve::VarEstimator{T}, s::Vector{T}) where T<:Real = begin
  ve.n += 1
  δ = s .- ve.μ
  ve.μ .+= δ ./ ve.n
  ve.M .+= δ .* (s .- ve.μ)
end

get_var(ve::VarEstimator) = begin
  @assert ve.n >= 2
  var = ve.M / (ve.n - 1)
  var = (ve.n / (ve.n + 5.0)) * var .+ 1e-3 * (5.0 / (ve.n + 5.0))
  return var
end

mutable struct WarmUpManager
  adapt_n   ::    Int
  params    ::    Dict{Symbol, Any}
  ve        ::    VarEstimator
end

getindex(wum::WarmUpManager, param) = wum.params[param]

setindex!(wum::WarmUpManager, value, param) = wum.params[param] = value

init_warm_up_params(vi::VarInfo, spl::Sampler{T}) where T<:Hamiltonian = begin
  D = length(vi[spl])
  ve = VarEstimator{Float64}(0, zeros(D), zeros(D))

  wum = WarmUpManager(1, Dict(), ve)

  # Pre-cond

  wum[:stds] = ones(D)

  # Dual averaging
  wum[:ϵ] = [] # why we are using a vector for ϵ
  restart_da(wum)
  wum[:n_warmup] = spl.alg.n_adapt
  wum[:δ] = spl.alg.delta

  # Initialize by Stan if Stan is installed
  is_init_by_stan = false
  @require Stan begin
    is_init_by_stan = true
    # Stan.Adapt
    adapt_conf = spl.info[:adapt_conf]
    # Hyper parameters for dual averaging
    wum[:γ] = adapt_conf.gamma
    wum[:t_0] = adapt_conf.t0
    wum[:κ] = adapt_conf.kappa
    # Three phases settings
    wum[:init_buffer] = adapt_conf.init_buffer
    wum[:term_buffer] = adapt_conf.term_buffer
    wum[:window_size] = adapt_conf.window
  end
  # If wum is not initialised by Stan (when Stan is not avaible),
  # initialise wum by common default values.
  if ~is_init_by_stan
    wum[:γ] = 0.05
    wum[:t_0] = 10.0
    wum[:κ] = 0.75
    wum[:init_buffer] = 75
    wum[:term_buffer] = 50
    wum[:window_size] = 25
  end
  wum[:next_window] = wum[:init_buffer] + wum[:window_size] - 1

  @debug wum.params

  spl.info[:wum] = wum
end

restart_da(wum::WarmUpManager) = begin
  wum[:m] = 0
  wum[:x_bar] = 0.0
  wum[:H_bar] = 0.0
end

update_da_μ(wum::WarmUpManager, ϵ::Float64) = begin
  # wum[:x_bar] = ϵ  # See NUTS paper sec 3.2.1
  wum[:μ] = log(10 * ϵ)
end

# Ref: https://github.com/stan-dev/stan/blob/develop/src/stan/mcmc/stepsize_adaptation.hpp
adapt_step_size!(wum::WarmUpManager, stats::Float64) = begin

  @debug "adapting step size ϵ..."
  @debug "current α = $(stats)"
  wum[:m] = wum[:m] + 1
  m = wum[:m]

  stats = stats > 1 ? 1 : stats # stats = δ - H_t, where δ is target accept prob

  γ = wum[:γ]; t_0 = wum[:t_0]; κ = wum[:κ]; δ = wum[:δ]
  μ = wum[:μ]; x_bar = wum[:x_bar]; H_bar = wum[:H_bar]

  η_H = 1.0 / (m + t_0)
  H_bar = (1.0 - η_H) * H_bar + η_H * (δ - stats)

  x = μ - H_bar * sqrt(m) / γ            # x ≡ logϵ
  η_x = m^(-κ)
  x_bar = (1.0 - η_x) * x_bar + η_x * x

  ϵ = exp(x)
  @debug "new ϵ = $(ϵ), old ϵ = $(wum[:ϵ][end])"

  if isnan(ϵ) || isinf(ϵ) || ϵ <= 1e-3
      @warn "Incorrect ϵ = $ϵ; ϵ_previous = $(wum[:ϵ][end]) is used instead."
  else
      push!(wum[:ϵ], ϵ)
      wum[:x_bar], wum[:H_bar] = x_bar, H_bar
  end

  if m == wum[:n_warmup]
    @debug " Adapted ϵ = $ϵ, $m HMC iterations is used for adaption."
  end

end

# Ref: https://github.com/stan-dev/stan/blob/develop/src/stan/mcmc/windowed_adaptation.hpp
in_adaptation(wum::WarmUpManager) = (wum.adapt_n >= wum[:init_buffer]) &&
                                    (wum.adapt_n < wum[:n_warmup] - wum[:term_buffer]) &&
                                    (wum.adapt_n != wum[:n_warmup])

is_window_end(wum::WarmUpManager) = (wum.adapt_n == wum[:next_window]) && (wum.adapt_n != wum[:n_warmup])

compute_next_window(wum::WarmUpManager) = begin

  if ~(wum[:next_window] == wum[:n_warmup] - wum[:term_buffer] - 1)

  wum[:window_size] *= 2
  wum[:next_window] = wum.adapt_n + wum[:window_size]

    if ~(wum[:next_window] == wum[:n_warmup] - wum[:term_buffer] - 1)

      next_window_boundary = wum[:next_window] + 2 * wum[:window_size]

      if (next_window_boundary >= wum[:n_warmup] - wum[:term_buffer])
        wum[:next_window] = wum[:n_warmup] - wum[:term_buffer] - 1
      end

    end

  end

end

# Ref: https://github.com/stan-dev/stan/blob/develop/src/stan/mcmc/hmc/nuts/adapt_diag_e_nuts.hpp
update_pre_cond!(wum::WarmUpManager, θ_new) = begin

  if in_adaptation(wum)

    # wum[:est_n] = wum[:est_n] + 1

    # delta = θ_new - wum[:θ_mean]
    # wum[:θ_mean] = wum[:θ_mean] + delta / wum[:est_n]
    # wum[:M2] = wum[:M2] + delta .* (θ_new - wum[:θ_mean])

    add_sample!(wum.ve, θ_new)

  end

  if is_window_end(wum)

    compute_next_window(wum)

    # var = wum[:M2] / (wum[:est_n] - 1.0)
    # var = (wum[:est_n] / (wum[:est_n] + 5.0)) * var + 1e-3 * (5.0 / (wum[:est_n] + 5.0))

    # https://github.com/stan-dev/stan/blob/develop/src/stan/mcmc/var_adaptation.hpp
    var = get_var(wum.ve)

    wum[:stds] = sqrt.(var)

    # reset_var_estimator(wum)

    reset!(wum.ve)

    return true

  end

  return false

end

adapt!(wum::WarmUpManager, stats::Float64, θ_new; adapt_ϵ = false, adapt_M = false) = begin

  if wum.adapt_n < wum[:n_warmup]

    if adapt_ϵ
        adapt_step_size!(wum, stats)
         if is_window_end(wum)
           ϵ = exp(wum[:x_bar])
           push!(wum[:ϵ], ϵ)
           update_da_μ(wum, ϵ)
           restart_da(wum)
         end
    end

    if adapt_M
        update_pre_cond!(wum, θ_new)  # window is updated implicitly.
    else   # update window explicitly.
        is_window_end(wum) && compute_next_window(wum)
    end

    wum.adapt_n += 1

  elseif wum.adapt_n == wum[:n_warmup]

    if adapt_ϵ
      ϵ = exp(wum[:x_bar])
      push!(wum[:ϵ], ϵ)
    end

  end

end


text:
function THRUST(POW,ALT,RMACH)
    # idle data
    A = [1060.0   670.0   880.0  1140.0  1500.0  1860.0
          635.0   425.0   690.0  1010.0  1330.0  1700.0
           60.0    25.0   345.0   755.0  1130.0  1525.0
        -1020.0  -710.0  -300.0   350.0   910.0  1360.0
        -2700.0 -1900.0 -1300.0  -247.0   600.0  1100.0
        -3600.0 -1400.0  -595.0  -342.0  -200.0   700.0]'
    # mil data
    B = [12680.0  9150.0  6200.0  3950.0  2450.0  1400.0
         12680.0  9150.0  6313.0  4040.0  2470.0  1400.0
         12610.0  9312.0  6610.0  4290.0  2600.0  1560.0
         12640.0  9839.0  7090.0  4660.0  2840.0  1660.0
         12390.0 10176.0  7750.0  5320.0  3250.0  1930.0
         11680.0  9848.0  8050.0  6100.0  3800.0  2310.0]'
    # mil data
    C = [20000.0 15000.0 10800.0  7000.0  4000.0  2500.0
         21420.0 15700.0 11225.0  7323.0  4435.0  2600.0
         22700.0 16860.0 12250.0  8154.0  5000.0  2835.0
         24240.0 18910.0 13760.0  9285.0  5700.0  3215.0
         26070.0 21075.0 15975.0 11115.0  6860.0  3950.0
         28886.0 23319.0 18300.0 13484.0  8642.0  5057.0]'
    #
    H = 0.0001*ALT
    I = floor(H)+1
    if (I>=6)
        I=5
    end
    DH = H-I+1
    RM = 5.0*RMACH
    M = floor(RM)+1
    if (M>=6)
        M=5
    end
    I=convert(Int64, I)
    M=convert(Int64, M)
    DM = RM-M+1
    CDH =1.0-DH
    S = B[I,M]*CDH + B[I+1,M]*DH
    T = B[I,M+1]*CDH + B[I+1,M+1]*DH
    TMIL = S + (T-S)*DM
    if (POW < 50.0)
        S = A[I,M]*CDH + A[I+1,M]  *DH
        T = A[I,M+1]*CDH + A[I+1,M+1]*DH
        TIDL = S + (T-S)*DM
        return thrustv = TIDL+(TMIL-TIDL)*POW*0.02
    else
        S = C[I,M]  *CDH + C[I+1,M]  *DH
        T = C[I,M+1]*CDH + C[I+1,M+1]*DH
        TMAX = S + (T-S)*DM
        return thrustv = TMIL+(TMAX-TMIL)*(POW-50.0)*0.02
    end
end


text:
module AgentsPlots

using Agents
using GraphPlot
using ColorTypes
using Compose
using Cairo
using Fontconfig

include("gridplots.jl")

end # module


text:
function minLC!(g::AbstractGraph)
    mg = construct_metagraph(g)
    origin_edges = edges(g)
    for i in 1:nv(g)
        mv = min_local_clustering(mg)
        sedges = getsaturate(mg, mv)
        add_edges!(mg, sedges)
        sedges::Vector{Vector{Int}} = [edge_m2s(mg, e[1], e[2]) for e in sedges]
        add_edges!(g, sedges)
        rem_vertex!(mg, mv)
    end
    data = Dict("added_edges" => setdiff(edges(g), origin_edges))
    return data
end

function minLC(g::AbstractGraph)
    cg = deepcopy(g)
    data = minLC!(cg)
    return cg, data
end

function maxLC!(g::AbstractGraph)
    mg = construct_metagraph(g)
    origin_edges = edges(g)
    for i in 1:nv(g)
        mv = max_local_clustering(mg)
        sedges = getsaturate(mg, mv)
        add_edges!(mg, sedges)
        sedges::Vector{Vector{Int}} = [edge_m2s(mg, e[1], e[2]) for e in sedges]
        add_edges!(g, sedges)
        rem_vertex!(mg, mv)
    end
    data = Dict("added_edges" => setdiff(edges(g), origin_edges))
    return data
end

function maxLC(g::AbstractGraph)
    cg = deepcopy(g)
    data = maxLC!(cg)
    return cg, data
end

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
using AutoGIS
using Test

@testset "AutoGIS.jl" begin
    # Write your tests here.
end

max_stars_repo_name: Auto-GIS-with-X/Auto-GIS-with-Julia

max_issues_repo_name: Auto-GIS-with-X/Auto-GIS-with-Julia

max_forks_repo_name: Auto-GIS-with-X/Auto-GIS-with-Julia


text:
const CONTROL0 = Vector{Symbol}()
const T1_CONTROL4 = [:Lbm; :Lmkequity; :LDat; :Lop]
const T1_CONTROL10 = [T1_CONTROL4; :L12ret12m; :L24ret12m;
  :fyret12m; :L12vol252d; :L24vol252d; :fyvol252d]

function maketable1(;trialrunonly::Bool=trialrunonly)::Nothing
  #load the data we will work with
  panel::DataFrame = constructportfolios(refreshportfolios=false, trialrunonly=trialrunonly)
  #table1panelA(panel)
  table1panelB(panel)
  #table1panelC(panel)
  #table1panelD(panel)

  return nothing
end

function table1panelA(panel::DataFrame;
    tablepath::String = TABLE_PATH,
    tablename::String = "t1panela-$(REPLICATION_TYPE[])",
    decimals::Int = DEFAULT_DECIMALS)
  local cqports::CQuartileConstructions = leverageconstructions()
  Fgrps::Vector{Symbol} = Fgroups(cqports)

  #NOTE: These are of the form mknegcash ~ α + β*mkt ...

  ################INPUTS############
  T1_PANEL_A_ROWKEYS::Vector{Symbol} = [:mean, :std, :stdxt]

  T1_PANEL_A_ROWLABEL::Dict = Dict(
    :mean=>"Mean",
    :std=>"SD",
    :stdxt=>"SD(X,T)")

  T1_PANEL_A_IDX::Dict = Dict(
    :mean=>1,
    :std=>2,
    :stdxt=>3)

  #################################

  #make the table
  colnames::Vector{Vector{String}} = [(string).(CQUART_FVALS)]
  descrownames::Vector{String} = (s->T1_PANEL_A_ROWLABEL[s]).(T1_PANEL_A_ROWKEYS)

  desccontent::Vector{Vector{String}} = (i->Vector{String}(undef, length(CQUART_FVALS))).(1:length(T1_PANEL_A_ROWKEYS))
  for c ∈ 1:length(CQUART_FVALS)
    local sspanels::GroupedDataFrame
    local ngroups::Int

    f::Symbol = CQUART_FVALS[c]
    allfields::Vector{Symbol} = [:date, :permno, f]
    spanel::SubDataFrame = view(panel, completecases(panel[!,allfields]), allfields)

    cmean::Float64 = mean(spanel[!, f])
    cstd::Float64 = std(spanel[!, f])

    cstdx::Float64 = 0.0
    ngroups = 0
    sspanels = groupby(spanel, :date)
    for sspanel ∈ sspanels
      cstdxi::Float64 = std(sspanel[!,f])
      if isfinite(cstdxi)
        cstdx += cstdxi
        ngroups+=1
      end
    end
    cstdx /= ngroups

    cstdt::Float64 = 0.0
    ngroups = 0
    sspanels = groupby(spanel, :permno)
    for sspanel ∈ groupby(spanel, :permno)
      cstdti::Float64 = std(sspanel[!,f])
      if isfinite(cstdti)
        cstdt += cstdti
        ngroups+=1
      end
    end
    cstdt /= ngroups

    desccontent[T1_PANEL_A_IDX[:mean]][c] = num2str(cmean, decimals)
    desccontent[T1_PANEL_A_IDX[:std]][c] = "($(num2str(cstd, decimals)))"
    desccontent[T1_PANEL_A_IDX[:stdxt]][c] = "($(num2str(cstdx, decimals)), $(num2str(cstdt, decimals)))"
  end

  tt::String = textable(colnames=colnames, descrownames=descrownames, desccontent=desccontent)

  writetextable(tt, path=tablepath, outname="$tablename.tex")
end

function table1panelB(panel::DataFrame;
    tablepath::String = TABLE_PATH,
    tablename::String = "t1panelb-$(REPLICATION_TYPE[])",
    decimals::Int = 1)
  local cqports::CQuartileConstructions = leverageconstructions()

  ################INPUTS############
  T1_PANEL_B_ROWKEYS::Vector{Symbol} = [
    :leverageQmeans,
    :quartilenum,
    :LLsd,
    :sd,
    :Δsd,
    :LLret12m,
    :ret12m,
    :Δret12m]

  T1_PANEL_B_ROWLABEL::Dict = Dict(
    :leverageQmeans=>"Leverage Q Means",
    :quartilenum=>"Quartile \\#",
    :LLsd=>"\\smallskip SD Net Lagged\$^2\$",
    :sd=>"SD Net Lead",
    :Δsd=>"SD Net 2Y Lead",
    :LLret12m=>"\\smallskip Compound Raw Lagged\$^2\$",
    :ret12m=>"Compound Raw Lead",
    :Δret12m=>"Compound Raw 2Y Delta")

  T1_PANEL_B_IDX::Dict = Dict(T1_PANEL_B_ROWKEYS[i]=>i for i ∈ 1:(length(T1_PANEL_B_ROWKEYS)))

  #the below is to help with future abstraction if necessary
  T1_PANEL_B_VALCOLIDX::Dict= Dict(f=>f for f ∈ CQUART_FVALS)

  #these are pairs of lagged and unlagged fields
  T1_PANEL_B_VOLFIELDS::NTuple{2,Symbol} = (:LLvol252d, :vol252d)
  T1_PANEL_B_RETFIELDS::NTuple{2,Symbol} = (:LLret12m, :ret12m)

  #################################

  #make the table
  colnames::Vector{Vector{String}} = [(string).(CQUART_FVALS)]
  descrownames::Vector{String} = (s->T1_PANEL_B_ROWLABEL[s]).(T1_PANEL_B_ROWKEYS)

  #pre-allocate for content
  desccontent::Vector{Vector{String}} = (
    i->Vector{String}(undef, length(CQUART_FVALS))).(1:length(T1_PANEL_B_ROWKEYS))

  #will be using this a lot
  @inline skipmean(v::AbstractVector{T where T<:MFloat64})::Float64 = mean(skipmissing(v))
  @inline n2s(x::T where T<:MFloat64, d::Real=decimals) = num2str(x, d, scalefactor=100.)

  volscale::Float64 = REPLICATION_TYPE[] == :daily ? 252^0.5 : 12^0.5

  for c ∈ 1:length(CQUART_FVALS)
    local sspanelq1::SubDataFrame
    local sspanelq4::SubDataFrame
    local allfields::Vector{Symbol}

    #acquire the fields we will use
    f::Symbol = CQUART_FVALS[c]
    name::Symbol = Symbol("P_", f)
    Fgroup::Symbol = cqports[name].Fgroup
    basicfields::Vector{Symbol} = [:date; :permno; Fgroup]
    local spanel::SubDataFrame = view(panel, completecases(panel[!,basicfields]),:)

    ###Value Rows
    Fval::Symbol = T1_PANEL_B_VALCOLIDX[f]
    #allfields = [basicfields; Fval]
    sspanelq1 = view(spanel, spanel[!,Fgroup] .== 1, :)
    sspanelq4 = view(spanel, spanel[!,Fgroup] .== 4, :)

    q1val::Float64 = mean(skipmissing(sspanelq1[!, Fval]))
    q4val::Float64 = mean(skipmissing(sspanelq4[!, Fval]))
    desccontent[T1_PANEL_B_IDX[:leverageQmeans]][c] = (
      "$(n2s(q1val, 1)) \\; $(n2s(q4val, 1))")
    desccontent[T1_PANEL_B_IDX[:quartilenum]][c] = ("\\text{Q1 Q4}")

    ###SD Rows
    #allfields = [basicfields; T1_PANEL_B_VOLFIELDS]
    sspanelq1 = view(spanel, spanel[!,Fgroup] .== 1, :)
    sspanelq4 = view(spanel, spanel[!,Fgroup] .== 4, :)

    q1LLvol::Float64 = skipmean(sspanelq1[!, T1_PANEL_B_VOLFIELDS[1]]) * volscale
    q4LLvol::Float64 = skipmean(sspanelq4[!, T1_PANEL_B_VOLFIELDS[1]]) * volscale
    desccontent[T1_PANEL_B_IDX[:LLsd]][c] = (
      "$(n2s(q1LLvol)) \\; $(n2s(q4LLvol))")

    q1vol::Float64 = skipmean(sspanelq1[!, T1_PANEL_B_VOLFIELDS[2]]) * volscale
    q4vol::Float64 = skipmean(sspanelq4[!, T1_PANEL_B_VOLFIELDS[2]]) * volscale
    desccontent[T1_PANEL_B_IDX[:sd]][c] = (
      "$(n2s(q1vol)) \\; $(n2s(q4vol))")

    Δq1vol::Float64 = skipmean(sspanelq1[!, T1_PANEL_B_VOLFIELDS[2]] .-
      sspanelq1[!, T1_PANEL_B_VOLFIELDS[1]]) * volscale
    Δq4vol::Float64 = skipmean(sspanelq4[!, T1_PANEL_B_VOLFIELDS[2]] .-
      sspanelq4[!, T1_PANEL_B_VOLFIELDS[1]]) * volscale
    desccontent[T1_PANEL_B_IDX[:Δsd]][c] = (
      "$(n2s(Δq1vol)) \\; $(n2s(Δq4vol))")

    ###Raw Ret Rows
    #allfields = [basicfields; T1_PANEL_B_RETFIELDS]
    sspanelq1 = view(spanel, spanel[!,Fgroup] .== 1, :)
    sspanelq4 = view(spanel, spanel[!,Fgroup] .== 4, :)

    q1LLret::Float64 = skipmean(sspanelq1[!, T1_PANEL_B_RETFIELDS[1]])
    q4LLret::Float64 = skipmean(sspanelq4[!, T1_PANEL_B_RETFIELDS[1]])
    desccontent[T1_PANEL_B_IDX[:LLret12m]][c] = (
      "$(n2s(q1LLret)) \\; $(n2s(q4LLret))")

    q1ret::Float64 = skipmean(sspanelq1[!, T1_PANEL_B_RETFIELDS[2]])
    q4ret::Float64 = skipmean(sspanelq4[!, T1_PANEL_B_RETFIELDS[2]])
    desccontent[T1_PANEL_B_IDX[:ret12m]][c] = (
      "$(n2s(q1ret)) \\; $(n2s(q4ret))")

    Δq1ret::Float64 = skipmean(sspanelq1[!, T1_PANEL_B_RETFIELDS[2]] .-
      sspanelq1[!, T1_PANEL_B_RETFIELDS[1]])
    Δq4ret::Float64 = skipmean(sspanelq4[!, T1_PANEL_B_RETFIELDS[2]] .-
      sspanelq4[!, T1_PANEL_B_RETFIELDS[1]])
    desccontent[T1_PANEL_B_IDX[:Δret12m]][c] = (
      "$(n2s(Δq1ret)) \\; $(n2s(Δq4ret))")
  end

  tt::String = textable(colnames=colnames, descrownames=descrownames, desccontent=desccontent)

  writetextable(tt, path=tablepath, outname="$tablename.tex")
end

function table1panelC(panel::DataFrame;
    tablepath::String = TABLE_PATH,
    tablename::String = "t1panelc-$(REPLICATION_TYPE[])",
    decimals::Int = DEFAULT_DECIMALS)
  local cqports::CQuartileConstructions = leverageconstructions()
  Fgrps::Vector{Symbol} = Fgroups(cqports)

  #NOTE: These are of the form mknegcash ~ α + β*mkt ...

  ################INPUTS############
  T1_PANEL_C_ROWKEYS::Vector{Symbol} = [:constant, :plusxmkt, :plussmbhml, :plusrmwcma, :plusumd]

  T1_PANEL_C_ROWLABEL::Dict = Dict(
    :constant=>"Constant",
    :plusxmkt=>"+XMKT",
    :plussmbhml=>"+SMB+HML",
    :plusrmwcma=>"+RMW+CMA",
    :plusumd=>"+UMD")

  T1_PANEL_C_CONTROLS::Dict = Dict(
    :constant=>CONTROL0,
    :plusxmkt=>[:mkt],
    :plussmbhml=>[:mkt, :P_ff3m_smb,	:P_ff3m_hml],
    :plusrmwcma=>[:mkt, :P_ff5_smb,	:P_ff5_hml,	:P_ff5_rmw,	:P_ff5_cma],
    :plusumd=>[:mkt, :P_ff5_smb,	:P_ff5_hml,	:P_ff5_rmw,	:P_ff5_cma, :P_ff3m_umd])

  T1_PANEL_C_FOCAL::Dict = Dict(
    :constant=>CQUART_FNAMES,
    :plusxmkt=>CQUART_FNAMES,
    :plussmbhml=>CQUART_FNAMES,
    :plusrmwcma=>CQUART_FNAMES,
    :plusumd=>CQUART_FNAMES)

  #################################
  fss::Vector{FactorSpecification} = Vector{FactorSpecification}()
  sizehint!(fss, length(T1_PANEL_C_ROWKEYS) * length(Fgrps))

  #need to re-shape for this exercise
  idcols = unique([:date; values(T1_PANEL_C_CONTROLS)...;])
  measurecols = CQUART_FNAMES
  allcols = unique([:date; values(T1_PANEL_C_CONTROLS)...; values(T1_PANEL_C_FOCAL)...;])
  minipanel::DataFrame = unique(panel[!,allcols])
  minipanel = melt(minipanel, idcols, measurecols, variable_name = :portname, value_name = :ret)
  #println(describe(minipanel))
  #println("Size minipanel: $(size(minipanel))")

  #make the regression specifications
  for rname ∈ T1_PANEL_C_ROWKEYS
    sallcols::Vector{Symbol} = [:ret; :date; :portname; T1_PANEL_C_CONTROLS[rname]]
    sminipanel = view(minipanel, completecases(minipanel[!,sallcols]), sallcols)
    push!(fss, FactorSpecification(sminipanel, :ret, :date, :portname,
      Ffactors = nothing,
      Fcharacteristics=nothing,
      Fcontrols = T1_PANEL_C_CONTROLS[rname],
      name=Symbol("row$(rname)"),
      skipcrosssectiondfcheck=true)) #only doing time series regressions
  end

  #run the regressions
  #WARNING: When BLAS paralllelization issues are solved replace with parallel=PARALLEL[]
  frs::Vector{FactorRegression} = Vector{FactorRegression}()

  #generate and push all the regressions
  #each call to manyFF pulls a row of regressions
  (fs->push!(frs, manyFF(fs, parallel=false)...)).(fss)

  #make the table
  colnames::Vector{Vector{String}} = [
    (string).(CQUART_FVALS), (i->"Alpha (T-stat)").(1:length(CQUART_FVALS))]
  descrownames::Vector{String} = (s->T1_PANEL_C_ROWLABEL[s]).(T1_PANEL_C_ROWKEYS)

  desccontent::Vector{Vector{String}} = factordesccontent(
    fr::FactorRegression->coeftstat(fr, i=1), fss=fss, frs=frs,
    rowkeys = T1_PANEL_C_ROWKEYS, colindex=T1_PANEL_C_FOCAL)

  tt::String = textable(colnames=colnames, descrownames=descrownames, desccontent=desccontent)

  writetextable(tt, path=tablepath, outname="$tablename.tex")
end

function table1panelD(panel::DataFrame;
    tablepath::String = TABLE_PATH,
    tablename::String = "t1paneld-$(REPLICATION_TYPE[])")
  local cqports::CQuartileConstructions = leverageconstructions()
  Fgrps::Vector{Symbol} = Fgroups(cqports)

  ################INPUTS############
  T1_PANEL_D_CONTROLS::Dict = Dict(
    :quartiles=>CONTROL0,
    :quartiles4ctrl=>T1_CONTROL4,
    :quartiles10ctrl=>T1_CONTROL10,
    :values=>CONTROL0,
    :values4ctrl=>T1_CONTROL4,
    :values10ctrl=>T1_CONTROL10,
    :quartiles_ac=>CONTROL0,
    :quartiles_lret=>T1_CONTROL10)

  T1_PANEL_D_RET::Dict = Dict(
    :quartiles=>:ret,
    :quartiles4ctrl=>:ret,
    :quartiles10ctrl=>:ret,
    :values=>:ret,
    :values4ctrl=>:ret,
    :values10ctrl=>:ret,
    :quartiles_ac=>:ret,
    :quartiles_lret=>:lret)

  T1_PANEL_D_ROWKEYS::Vector{Symbol} = [:quartiles, :quartiles4ctrl, :quartiles10ctrl, :values,
    :values4ctrl, :values10ctrl, :quartiles_ac, :quartiles_lret]

  T1_PANEL_D_ROWLABEL::Dict = Dict(
    :quartiles=>"Leverage Quartiles",
    :quartiles4ctrl=>"... +4 vars",
    :quartiles10ctrl=>"... +10 vars",
    :values=>"Leverage Values",
    :values4ctrl=>"... +4 vars",
    :values10ctrl=>"... +10 vars",
    :quartiles_ac=>"\\midrule Asset-Controlled Quartiles",
    :quartiles_lret=>"Log Returns +10 vars")

  T1_PANEL_D_FOCAL::Dict = Dict(
    :quartiles=>CQUART_FGROUPS,
    :quartiles4ctrl=>CQUART_FGROUPS,
    :quartiles10ctrl=>CQUART_FGROUPS,
    :values=>CQUART_FVALS,
    :values4ctrl=>CQUART_FVALS,
    :values10ctrl=>CQUART_FVALS,
    :quartiles_ac=>CQUART_FGROUPS_ALT,
    :quartiles_lret=>CQUART_FGROUPS)
  #################################

  fss::Vector{FactorSpecification} = Vector{FactorSpecification}()
  sizehint!(fss, length(T1_PANEL_D_ROWKEYS) * length(Fgrps))

  #make the regression specifications
  @time for rname ∈ T1_PANEL_D_ROWKEYS
    for focal ∈ T1_PANEL_D_FOCAL[rname]
      push!(fss, FactorSpecification(panel, T1_PANEL_D_RET[rname], :date, :permno,
        Ffactors=nothing,
        Fcharacteristics = [focal],
        Fcontrols = T1_PANEL_D_CONTROLS[rname],
        name=Symbol("row$(rname)_col$(focal)"),
        skiptimeseriesdfcheck=true)) #only doing crossectional regressions
    end
  end

  #Run the regressions
  #WARNING: When BLAS paralllelization issues are solved replace with parallel=PARALLEL[]
  frs::Vector{FactorRegression} = (fs->characteristicFM(fs, parallel=false)).(fss)

  #make the table
  colnames::Vector{Vector{String}} = [
    (string).(CQUART_FVALS), (i->"Coef (T-stat)").(1:length(CQUART_FVALS))]
  descrownames::Vector{String} = (s->T1_PANEL_D_ROWLABEL[s]).(T1_PANEL_D_ROWKEYS)
  desccontent::Vector{Vector{String}} = factordesccontent(fss=fss, frs=frs,
    rowkeys = T1_PANEL_D_ROWKEYS, colindex=T1_PANEL_D_FOCAL)

  tt::String = textable(colnames=colnames,
    descrownames=descrownames, desccontent=desccontent)

  writetextable(tt, path=tablepath, outname="$tablename.tex")
end

function coeftstat(fr::FactorRegression; decimals::Int=DEFAULT_DECIMALS, i::Int=2,
  )::String
  λᵢ::Float64 = fr.λᵢ[i] * 100
  tstat::Float64 = fr.λᵢ[i]/fr.σᵢ[i]
  return "$(num2str(λᵢ, decimals))\\;($(num2str(tstat, decimals)))"
end

#this is a kernel function to generate the content matrix under a number of different scenarios
function factordesccontent(f::Function = coeftstat;
    fss::Vector{FactorSpecification} = error("fss is required for factordesccontent"),
    frs::Vector{FactorRegression} = error("frs is required for factordesccontent"),
    rowkeys::Vector{Symbol} = error("Rowkeys is required"),
    colindex::Dict = error("colindex is required for factordesccontent"))

    Nrows::Int = length(rowkeys)
    Ncols::Int = length(collect(values(colindex))[1])

    #pre-allocate
    desccontent::Vector{Vector{String}} = (i->Vector{String}(undef, Ncols)).(1:Nrows)

    #convenience dictionaries
    frsidx = Dict(fr.name=>fr for fr ∈ frs)
    rowcolnameidx = Dict(
      (rname::Symbol, focal::Symbol)=>Symbol("row$(rname)_col$(focal)")
      for rname::Symbol ∈ rowkeys for focal::Symbol ∈ colindex[rname])

    #write the content matrix
    for (r,rname) ∈ enumerate(rowkeys)
      for (c, focal) ∈ enumerate(colindex[rname])
        fr::FactorRegression = frsidx[rowcolnameidx[(rname, focal)]]
        desccontent[r][c] = f(fr)
      end
    end

    return desccontent
  end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
julia> using TransmuteDims, BenchmarkTools

julia> VERSION
v"1.7.0-DEV.406"

julia> V1 = randn(1000);

julia> M1 = randn(1000,1000);
julia> M2 = randn(1000,1000);
julia> M3 = similar(M1);

julia> T1 = randn(100,100,100);
julia> T2 = randn(100,100,100);
julia> T3 = similar(T1);

julia> G2 = randn(100,1,100);

#========== Constructors ==========#
# Aim to be faster than PermutedDimsArray.

julia> @btime PermutedDimsArray($T1, (3,2,1));
  381.424 ns (4 allocations: 176 bytes)

julia> @btime TransmutedDimsArray($T1, (3,2,1));
  338.963 ns (3 allocations: 80 bytes)

julia> @btime transmute($T1, (3,2,1));  # fast at last!
  1.430 ns (0 allocations: 0 bytes)

julia> @btime transmute($T1, Val((3,2,1)));
  1.430 ns (0 allocations: 0 bytes)

# ... transpose

julia> @btime transpose($M1);
  1.430 ns (0 allocations: 0 bytes)

julia> @btime transmute($M1, (2,1));
  1.430 ns (0 allocations: 0 bytes)

julia> @btime transmute($M1, Val((2,1)));
  1.430 ns (0 allocations: 0 bytes)

# ... reshape

julia> @btime reshape($M1, (1000, 1, 1000));
  31.764 ns (1 allocation: 80 bytes)

julia> @btime transmute($M1, (1,0,2));
  40.796 ns (2 allocations: 112 bytes)

julia> @btime transmute($M1, Val((1,0,2)));
  40.613 ns (2 allocations: 112 bytes)

# Best-case constant prop:

julia> const cT1 = randn(100,100,100);

julia> @btime (() -> PermutedDimsArray(cT1,(3,2,1)))();
  41.292 ns (1 allocation: 96 bytes)

julia> @btime (() -> TransmutedDimsArray(cT1,(3,2,1)))();
  1.428 ns (0 allocations: 0 bytes)

julia> @btime (() -> transmute(cT1,(3,2,1)))();
  1.428 ns (0 allocations: 0 bytes)

# ... and a harder case:

julia> @btime (x -> PermutedDimsArray(x, (3,2,1)))($T1);
  385.429 ns (4 allocations: 176 bytes)

julia> @btime (x -> TransmutedDimsArray(x, (3,2,1)))($T1);
  338.977 ns (3 allocations: 80 bytes)

julia> @btime (x -> transmute(x, (3,2,1)))($T1);
  1.430 ns (0 allocations: 0 bytes)

#========== Broadcasting ==========#
# Here PermutedDimsArray is great, aim to be no slower.
# Times at right are from earlier version, Julia 1.2?

julia> @btime $M3 .= $M1 .* $M2;                              # 766.213 μs
  771.604 μs (0 allocations: 0 bytes)

julia> @btime $M3 .= $M1 .* $(PermutedDimsArray(M2,(1,2)));   # 769.648 μs
  772.055 μs (0 allocations: 0 bytes)

julia> @btime $M3 .= $M1 .* $(TransmutedDimsArray(M2,(1,2))); # 769.523 μs with IndexCartesian()
  771.612 μs (0 allocations: 0 bytes)

# transposed

julia> @btime $M3 .= $M1 .* transpose($M2);                 # 2.179 ms
  1.890 ms (0 allocations: 0 bytes)

julia> @btime $M3 .= $M1 .* PermutedDimsArray($M2,(2,1));   # 2.186 ms
  1.897 ms (6 allocations: 240 bytes)

julia> @btime $M3 .= $M1 .* $(PermutedDimsArray(M2,(2,1))); # 2.179 ms
  1.889 ms (0 allocations: 0 bytes)

julia> @btime $M3 .= $M1 .* $(TransmutedDimsArray(M2,(2,1))); # 2.179 ms
  1.877 ms (0 allocations: 0 bytes)

#========== Reductions ==========#
# These should be faster than PermutedDimsArray, passed through.

julia> @btime sum($M1);                              #   290.741 μs
  294.625 μs (0 allocations: 0 bytes)

julia> @btime sum($(PermutedDimsArray(M2,(1,2))));   # 1.399 ms
  1.493 ms (0 allocations: 0 bytes)

julia> @btime sum($(TransmutedDimsArray(M2,(1,2)))); #   293.250 μs with method!
  301.680 μs (0 allocations: 0 bytes)

# transposed

julia> @btime sum($(transpose(M2)));                 # 2.341 ms
  2.044 ms (0 allocations: 0 bytes)

julia> @btime sum($(PermutedDimsArray(M2,(2,1))));   # 2.450 ms
  2.619 ms (0 allocations: 0 bytes)

julia> @btime sum($(TransmutedDimsArray(M2,(2,1))));
  301.856 μs (0 allocations: 0 bytes)

# partial

julia> @btime sum!($V1, $M1);                              # 291.655 μs
  316.322 μs (0 allocations: 0 bytes)

julia> @btime sum!($V1, $(PermutedDimsArray(M2,(1,2))));   # 292.307 μs
  316.824 μs (0 allocations: 0 bytes)

julia> @btime sum!($V1, $(TransmutedDimsArray(M2,(1,2)))); # 290.363 μs
  316.575 μs (0 allocations: 0 bytes)

# partial, transposed

julia> @btime sum!($V1, $(transpose(M1)));
  1.659 ms (0 allocations: 0 bytes)

julia> @btime sum!($V1, $(PermutedDimsArray(M2,(2,1))));
  1.540 ms (0 allocations: 0 bytes)

julia> @btime sum!($V1, $(TransmutedDimsArray(M2,(2,1))));
  294.152 μs (3 allocations: 112 bytes)  # why the allocations?

# tiny

julia> M0 = rand(5,5);

julia> @btime sum($M0, dims=2);
  238.289 ns (5 allocations: 208 bytes)

julia> @btime sum($(transpose(M0)), dims=1);
  73.177 ns (1 allocation: 128 bytes)

julia> @btime sum($(TransmutedDimsArray(M0,(2,1))), dims=1);
  73.045 ns (1 allocation: 128 bytes)

#========== Dropdims ==========#
# Reshape is faster than Base's dropdims.

julia> @btime dropdims($G2, dims=2);
  337.256 ns (9 allocations: 224 bytes)

julia> @btime transmute($G2, (1,3));  # reshape
  39.945 ns (2 allocations: 96 bytes)

julia> @btime transmute($G2, (3,1));  # wrap
  4.195 ns (0 allocations: 0 bytes)

#========== Collected ==========#
#

julia> @btime permutedims!($M3, $M2, (2,1));
  1.509 ms (0 allocations: 0 bytes)

julia> @btime copyto!($M3, $(transpose(M2)));
  2.467 ms (0 allocations: 0 bytes)

julia> @btime copyto!($M3, $(PermutedDimsArray(M2, (2,1)))); # same with M3 .=
  3.361 ms (0 allocations: 0 bytes)

julia> @btime transmutedims!($M3, $M2, (2,1));
  2.647 ms (0 allocations: 0 bytes)

julia> @btime copyto!($M3, $(TransmutedDimsArray(M2, (2,1))));
  606.967 μs (62 allocations: 6.66 KiB)

# three dimensions

julia> @btime permutedims!($T3, $T1, (3,2,1));
  2.283 ms (0 allocations: 0 bytes)

julia> @btime copyto!($T3, $(PermutedDimsArray(T1, (3,2,1))));
  4.528 ms (0 allocations: 0 bytes)

julia> @btime copyto!($T3, $(TransmutedDimsArray(T1, (3,2,1))));
  542.873 μs (61 allocations: 7.22 KiB)

# allocating forms

julia> @btime permutedims($M2);
  1.599 ms (2 allocations: 7.63 MiB)

julia> @btime transmutedims($M2);
  808.915 μs (65 allocations: 7.64 MiB)

julia> @btime permutedims($T1, (3,2,1));
  2.459 ms (2 allocations: 7.63 MiB)

julia> @btime transmutedims($T1, (3,2,1));
  737.701 μs (62 allocations: 7.64 MiB)

# fast library behind this

julia> using Strided

julia> @btime @strided permutedims!($M3, $M2, (2,1));
  605.804 μs (66 allocations: 6.81 KiB)

julia> @btime @strided permutedims!($T3, $T1, (3,2,1));
  520.262 μs (65 allocations: 7.38 KiB)

#========== Acceleration ==========#

julia> @btime $M3 .= exp.($M1');
  7.868 ms (0 allocations: 0 bytes)

julia> @btime $M3 .= exp.($(TransmutedDimsArray(M1,(2,1))));
  7.642 ms (0 allocations: 0 bytes)

julia> using Strided

julia> @btime @strided $M3 .= exp.($M1');  # multi-threaded
  3.637 ms (64 allocations: 6.66 KiB)

julia> @btime @strided $M3 .= exp.($(TransmutedDimsArray(M1,(2,1))));
  3.640 ms (64 allocations: 6.66 KiB)

julia> using LoopVectorization

julia> @btime @avx $M3 .= exp.($M1');  # vectorised
  1.893 ms (0 allocations: 0 bytes)

julia> @btime @avx $M3 .= exp.($(TransmutedDimsArray(M1,(2,1))));
  7.547 ms (0 allocations: 0 bytes)

julia> @btime @avx $M3 .= exp.(transmutedims($M1,(2,1)));
  1.979 ms (65 allocations: 7.64 MiB)


text:
using Test

# TODO Change this to test practice exercises when adding the first practice exercise
for exercise in readdir(joinpath("exercises", "concept"))
    # Allow only testing specified exercises
    if !isempty(ARGS) && !(exercise in ARGS)
        continue
    end

    exercise_path = joinpath(joinpath("exercises", "concept"), exercise)
    isdir(exercise_path) || continue

    # Create an anonymous module so that exercises are tested in separate scopes
    m = Module()

    # When testing the example solution, all tests must pass, even those that are marked as skipped or broken.
    # The student will not be affected by this.
    # Overwrite @test_skip and @test_broken with @test
    Core.eval(m, :(using Test))
    @eval m $(Symbol("@test_skip")) = $(Symbol("@test"))
    @eval m $(Symbol("@test_broken")) = $(Symbol("@test"))

    # runtests.jl includes the solution by calling `include("slug.jl")`
    # Our anonymous module doesn't have `include(s::String)` defined,
    # so we define our own. We manually include the example solution in our
    # anonymous module, so we can define `m.include(s::String)` to do nothing.
    Core.eval(m, :(include(s) = nothing))
    Base.include(m, joinpath(exercise_path, joinpath(".meta", "example.jl")))
    @info "Testing $exercise"
    Base.include(m, joinpath(exercise_path, "runtests.jl"))

    println() # to make the output more readable
end