text:
__precompile__()

module GR

using Compat

if VERSION > v"0.5-"
  if Sys.KERNEL == :NT
    const os = :Windows
  else
    const os = Sys.KERNEL
  end
else
  const os = OS_NAME
end

const None = Union{}

export
  opengks,
  closegks,
  inqdspsize,
  openws,
  closews,
  activatews,
  deactivatews,
  clearws,
  updatews,
  polyline,
  polymarker,
  text,
  inqtext,
  fillarea,
  cellarray,
  spline,
  gridit,
  setlinetype,
  setlinewidth,
  setlinecolorind,
  setmarkertype,
  setmarkersize,
  setmarkercolorind,
  settextfontprec,
  setcharexpan,
  setcharspace,
  settextcolorind,
  setcharheight,
  setcharup,
  settextpath,
  settextalign,
  setfillintstyle,
  setfillstyle,
  setfillcolorind,
  setcolorrep,
  setscale,
  setwindow,
  setviewport,
  selntran,
  setclip,
  setwswindow,
  setwsviewport,
  createseg,
  copyseg,
  redrawsegws,
  setsegtran,
  closeseg,
  emergencyclosegks,
  updategks,
  setspace,
  textext,
  inqtextext,
  axes,
  grid,
  grid3d,
  verrorbars,
  herrorbars,
  polyline3d,
  polymarker3d,
  axes3d,
  titles3d,
  surface,
  contour,
  setcolormap,
  colormap,
  inqcolor,
  inqcolorfromrgb,
  hsvtorgb,
  tick,
  validaterange,
  adjustlimits,
  adjustrange,
  beginprint,
  beginprintext,
  endprint,
  ndctowc,
  wctondc,
  wc3towc,
  drawrect,
  fillrect,
  drawarc,
  fillarc,
  drawpath,
  setarrowstyle,
  setarrowsize,
  drawarrow,
  readimage,
  drawimage,
  importgraphics,
  setshadow,
  settransparency,
  setcoordxform,
  begingraphics,
  endgraphics,
  getgraphics,
  drawgraphics,
  mathtex,
  selectcontext,
  trisurface,
  gradient,
  quiver,
  # Convenience functions
  jlgr,
  figure,
  hold,
  subplot,
  plot,
  oplot,
  semilogx,
  semilogy,
  loglog,
  scatter,
  stem,
  histogram,
  contourf,
  wireframe,
  plot3,
  scatter3,
  title,
  xlabel,
  ylabel,
  legend,
  xlim,
  ylim,
  savefig,
  meshgrid,
  peaks,
  imshow,
  isosurface,
  cart2sph,
  sph2cart,
  polar,
  trisurf,
  libGR3,
  gr3,
  isinline,
  inline,
  show

mime_type = None
figure_count = None
msgs = None
have_clear_output = isinteractive() && isdefined(Main, :IJulia) &&
                    isdefined(Main.IJulia, :clear_output)

function __init__()
    global libGR, libGR3
    if "GRDIR" in keys(ENV)
        grdir = ENV["GRDIR"]
    elseif isdir(joinpath(homedir(), "gr"), "fonts")
        grdir = joinpath(homedir(), "gr")
    else
        grdir = None
        for d in ("/opt", "/usr/local", "/usr")
            if isdir(joinpath(d, "gr", "fonts"))
                grdir = joinpath(d, "gr")
                break
            end
        end
    end
    if grdir == None
        grdir = joinpath(Pkg.dir(), "GR", "deps", "gr")
    end
    ENV["GRDIR"] = grdir
    ENV["GKS_FONTPATH"] = grdir
    if contains(grdir, "site-packages")
        const libGR = joinpath(grdir, "libGR.so")
        ENV["GKS_FONTPATH"] = grdir
    elseif os != :Windows
        const libGR = joinpath(grdir, "lib", "libGR.so")
    else
        const libGR = joinpath(grdir, "libGR.dll")
    end
    if !isfile(libGR)
        println("Unable to load GR framework runtime environment")
        exit(-1)
    end
    const libGR3 = replace(libGR, "libGR", "libGR3")
end

function opengks()
  ccall( (:gr_opengks, libGR),
        Void,
        ()
        )
end

function closegks()
  ccall( (:gr_closegks, libGR),
        Void,
        ()
        )
end

function inqdspsize()
  mwidth = Cdouble[0]
  mheight = Cdouble[0]
  width = Cint[0]
  height = Cint[0]
  ccall( (:gr_inqdspsize, libGR),
        Void,
        (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}),
        mwidth, mheight, width, height)
  return mwidth[1], mheight[1], width[1], height[1]
end

function openws(workstation_id::Int, connection, workstation_type::Int)
  ccall( (:gr_openws, libGR),
        Void,
        (Int32, Ptr{Cchar}, Int32),
        workstation_id, connection, workstation_type)
end

function closews(workstation_id::Int)
  ccall( (:gr_closews, libGR),
        Void,
        (Int32, ),
        workstation_id)
end

function activatews(workstation_id::Int)
  ccall( (:gr_activatews, libGR),
        Void,
        (Int32, ),
        workstation_id)
end

function deactivatews(workstation_id::Int)
  ccall( (:gr_deactivatews, libGR),
        Void,
        (Int32, ),
        workstation_id)
end

function clearws()
  global msgs
  try
    if isinline() && have_clear_output
      IJulia.clear_output(true)
    end
  catch
    have_clear_output = false
  end
  if msgs != None
    begingraphics("")
  end
  ccall( (:gr_clearws, libGR),
        Void,
        ()
        )
end

function updatews()
  global msgs
  if msgs != None
    endgraphics()
  end
  ccall( (:gr_updatews, libGR),
        Void,
        ()
        )
end

function polyline(x, y)
  assert(length(x) == length(y))
  n = length(x)
  ccall( (:gr_polyline, libGR),
        Void,
        (Int32, Ptr{Float64}, Ptr{Float64}),
        n, convert(Vector{Float64}, x), convert(Vector{Float64}, y))
end

function polymarker(x, y)
  assert(length(x) == length(y))
  n = length(x)
  ccall( (:gr_polymarker, libGR),
        Void,
        (Int32, Ptr{Float64}, Ptr{Float64}),
        n, convert(Vector{Float64}, x), convert(Vector{Float64}, y))
end

function latin1(string)
  b = convert(Array{UInt8}, string)
  s = zeros(UInt8, length(string))
  len = 0
  mask = 0
  for c in b
    if c != 0xc2 && c != 0xc3
      len += 1
      s[len] = c | mask
    end
    if c == 0xc3
      mask = 0x40
    else
      mask = 0
    end
  end
  return s
end

function text(x::Real, y::Real, string)
  ccall( (:gr_text, libGR),
        Void,
        (Float64, Float64, Ptr{UInt8}),
        x, y, latin1(string))
end

function inqtext(x, y, string)
  tbx = Cdouble[0, 0, 0, 0]
  tby = Cdouble[0, 0, 0, 0]
  ccall( (:gr_inqtext, libGR),
        Void,
        (Float64, Float64, Ptr{UInt8}, Ptr{Cdouble}, Ptr{Cdouble}),
        x, y, latin1(string), tbx, tby)
  return tbx, tby
end

function fillarea(x, y)
  assert(length(x) == length(y))
  n = length(x)
  ccall( (:gr_fillarea, libGR),
        Void,
        (Int32, Ptr{Float64}, Ptr{Float64}),
        n, convert(Vector{Float64}, x), convert(Vector{Float64}, y))
end

function cellarray(xmin::Real, xmax::Real, ymin::Real, ymax::Real, dimx::Int, dimy::Int, color)
  if ndims(color) == 2
    color = reshape(color, dimx * dimy)
  end
  ccall( (:gr_cellarray, libGR),
        Void,
        (Float64, Float64, Float64, Float64, Int32, Int32, Int32, Int32, Int32, Int32, Ptr{Int32}),
        xmin, xmax, ymin, ymax, dimx, dimy, 1, 1, dimx, dimy, convert(Vector{Int32}, color))
end

function spline(x, y, m, method)
  assert(length(x) == length(y))
  n = length(x)
  ccall( (:gr_spline, libGR),
        Void,
        (Int32, Ptr{Float64}, Ptr{Float64}, Int32, Int32),
        n, convert(Vector{Float64}, x), convert(Vector{Float64}, y), m, method)
end

function gridit(xd, yd, zd, nx, ny)
  assert(length(xd) == length(yd) == length(zd))
  nd = length(xd)
  x = Cdouble[1 : nx ;]
  y = Cdouble[1 : ny ;]
  z = Cdouble[1 : nx*ny ;]
  ccall( (:gr_gridit, libGR),
        Void,
        (Int32, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Int32, Int32, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
        nd, convert(Vector{Float64}, xd), convert(Vector{Float64}, yd), convert(Vector{Float64}, zd), nx, ny, x, y, z)
  return x, y, z
end

function setlinetype(style::Int)
  ccall( (:gr_setlinetype, libGR),
        Void,
        (Int32, ),
        style)
end

function setlinewidth(width::Real)
  ccall( (:gr_setlinewidth, libGR),
        Void,
        (Float64, ),
        width)
end

function setlinecolorind(color::Int)
  ccall( (:gr_setlinecolorind, libGR),
        Void,
        (Int32, ),
        color)
end

function setmarkertype(mtype::Int)
  ccall( (:gr_setmarkertype, libGR),
        Void,
        (Int32, ),
        mtype)
end

function setmarkersize(mtype::Real)
  ccall( (:gr_setmarkersize, libGR),
        Void,
        (Float64, ),
        mtype)
end

function setmarkercolorind(color::Int)
  ccall( (:gr_setmarkercolorind, libGR),
        Void,
        (Int32, ),
        color)
end

function settextfontprec(font::Int, precision::Int)
  ccall( (:gr_settextfontprec, libGR),
        Void,
        (Int32, Int32),
        font, precision)
end

function setcharexpan(factor::Real)
  ccall( (:gr_setcharexpan, libGR),
        Void,
        (Float64, ),
        factor)
end

function setcharspace(spacing::Real)
  ccall( (:gr_setcharspace, libGR),
        Void,
        (Float64, ),
        spacing)
end

function settextcolorind(color::Int)
  ccall( (:gr_settextcolorind, libGR),
        Void,
        (Int32, ),
        color)
end

function setcharheight(height::Real)
  ccall( (:gr_setcharheight, libGR),
        Void,
        (Float64, ),
        height)
end

function setcharup(ux::Real, uy::Real)
  ccall( (:gr_setcharup, libGR),
        Void,
        (Float64, Float64),
        ux, uy)
end

function settextpath(path::Int)
  ccall( (:gr_settextpath, libGR),
        Void,
        (Int32, ),
        path)
end

function settextalign(horizontal::Int, vertical::Int)
  ccall( (:gr_settextalign, libGR),
        Void,
        (Int32, Int32),
        horizontal, vertical)
end

function setfillintstyle(style::Int)
  ccall( (:gr_setfillintstyle, libGR),
        Void,
        (Int32, ),
        style)
end

function setfillstyle(index::Int)
  ccall( (:gr_setfillstyle, libGR),
        Void,
        (Int32, ),
        index)
end

function setfillcolorind(color::Int)
  ccall( (:gr_setfillcolorind, libGR),
        Void,
        (Int32, ),
        color)
end

function setcolorrep(index::Int, red::Real, green::Real, blue::Real)
  ccall( (:gr_setcolorrep, libGR),
        Void,
        (Int32, Float64, Float64, Float64),
        index, red, green, blue)
end

function setscale(options::Int)
  ccall( (:gr_setscale, libGR),
        Void,
        (Int32, ),
        options)
end

function setwindow(xmin::Real, xmax::Real, ymin::Real, ymax::Real)
  ccall( (:gr_setwindow, libGR),
        Void,
        (Float64, Float64, Float64, Float64),
        xmin, xmax, ymin, ymax)
end

function setviewport(xmin::Real, xmax::Real, ymin::Real, ymax::Real)
  ccall( (:gr_setviewport, libGR),
        Void,
        (Float64, Float64, Float64, Float64),
        xmin, xmax, ymin, ymax)
end

function selntran(transform::Int)
  ccall( (:gr_selntran, libGR),
        Void,
        (Int32, ),
        transform)
end

function setclip(indicator::Int)
  ccall( (:gr_setclip, libGR),
        Void,
        (Int32, ),
        indicator)
end

function setwswindow(xmin::Real, xmax::Real, ymin::Real, ymax::Real)
  ccall( (:gr_setwswindow, libGR),
        Void,
        (Float64, Float64, Float64, Float64),
        xmin, xmax, ymin, ymax)
end

function setwsviewport(xmin::Real, xmax::Real, ymin::Real, ymax::Real)
  ccall( (:gr_setwsviewport, libGR),
        Void,
        (Float64, Float64, Float64, Float64),
        xmin, xmax, ymin, ymax)
end

function createseg(segment::Int)
  ccall( (:gr_createseg, libGR),
        Void,
        (Int32, ),
        segment)
end

function copyseg(segment::Int)
  ccall( (:gr_copyseg, libGR),
        Void,
        (Int32, ),
        segment)
end

function redrawsegws()
  ccall( (:gr_redrawsegws, libGR),
        Void,
        ()
        )
end

function setsegtran(segment::Int, fx::Real, fy::Real, transx::Real, transy::Real, phi::Real, scalex::Real, scaley::Real)
  ccall( (:gr_setsegtran, libGR),
        Void,
        (Int32, Float64, Float64, Float64, Float64, Float64, Float64, Float64),
        segment, fx, fy, transx, transy, phi, scalex, scaley)
end

function closeseg()
  ccall( (:gr_closeseg, libGR),
        Void,
        ()
        )
end

function emergencyclosegks()
  ccall( (:gr_emergencyclosegks, libGR),
        Void,
        ()
        )
end

function updategks()
  ccall( (:gr_updategks, libGR),
        Void,
        ()
        )
end

function setspace(zmin::Real, zmax::Real, rotation::Int, tilt::Int)
  ccall( (:gr_setspace, libGR),
        Void,
        (Float64, Float64, Int32, Int32),
        zmin, zmax, rotation, tilt)
end

function textext(x::Real, y::Real, string)
  ccall( (:gr_textext, libGR),
        Void,
        (Float64, Float64, Ptr{UInt8}),
        x, y, latin1(string))
end

function inqtextext(x::Real, y::Real, string)
  tbx = Cdouble[0, 0, 0, 0]
  tby = Cdouble[0, 0, 0, 0]
  ccall( (:gr_inqtextext, libGR),
        Void,
        (Float64, Float64, Ptr{UInt8}, Ptr{Cdouble}, Ptr{Cdouble}),
        x, y, latin1(string), tbx, tby)
  return tbx, tby
end

function axes(x_tick::Real, y_tick::Real, x_org::Real, y_org::Real, major_x::Int, major_y::Int, tick_size::Real)
  ccall( (:gr_axes, libGR),
        Void,
        (Float64, Float64, Float64, Float64, Int32, Int32, Float64),
        x_tick, y_tick, x_org, y_org, major_x, major_y, tick_size)
end

function grid(x_tick::Real, y_tick::Real, x_org::Real, y_org::Real, major_x::Int, major_y::Int)
  ccall( (:gr_grid, libGR),
        Void,
        (Float64, Float64, Float64, Float64, Int32, Int32),
        x_tick, y_tick, x_org, y_org, major_x, major_y)
end

function grid3d(x_tick::Real, y_tick::Real, z_tick::Real, x_org::Real, y_org::Real, z_org::Real, major_x::Int, major_y::Int, major_z::Int)
  ccall( (:gr_grid3d, libGR),
        Void,
        (Float64, Float64, Float64, Float64, Float64, Float64, Int32, Int32, Int32),
        x_tick, y_tick, z_tick, x_org, y_org, z_org, major_x, major_y, major_z)
end

function verrorbars(px, py, e1, e2)
  assert(length(px) == length(py) == length(e1) == length(e2))
  n = length(px)
  ccall( (:gr_verrorbars, libGR),
        Void,
        (Int32, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}),
        n, convert(Vector{Float64}, px), convert(Vector{Float64}, py), convert(Vector{Float64}, e1), convert(Vector{Float64}, e2))
end

function herrorbars(px, py, e1, e2)
  assert(length(px) == length(py) == length(e1) == length(e2))
  n = length(px)
  ccall( (:gr_herrorbars, libGR),
        Void,
        (Int32, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}),
        n, convert(Vector{Float64}, px), convert(Vector{Float64}, py), convert(Vector{Float64}, e1), convert(Vector{Float64}, e2))
end

function polyline3d(px, py, pz)
  assert(length(px) == length(py) == length(pz))
  n = length(px)
  ccall( (:gr_polyline3d, libGR),
        Void,
        (Int32, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}),
        n, convert(Vector{Float64}, px), convert(Vector{Float64}, py), convert(Vector{Float64}, pz))
end

function polymarker3d(px, py, pz)
  assert(length(px) == length(py) == length(pz))
  n = length(px)
  ccall( (:gr_polymarker3d, libGR),
        Void,
        (Int32, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}),
        n, convert(Vector{Float64}, px), convert(Vector{Float64}, py), convert(Vector{Float64}, pz))
end

function axes3d(x_tick::Real, y_tick::Real, z_tick::Real, x_org::Real, y_org::Real, z_org::Real, major_x::Int, major_y::Int, major_z::Int, tick_size::Real)
  ccall( (:gr_axes3d, libGR),
        Void,
        (Float64, Float64, Float64, Float64, Float64, Float64, Int32, Int32, Int32, Float64),
        x_tick, y_tick, z_tick, x_org, y_org, z_org, major_x, major_y, major_z, tick_size)
end

function titles3d(x_title, y_title, z_title)
  ccall( (:gr_titles3d, libGR),
        Void,
        (Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}),
        latin1(x_title), latin1(y_title), latin1(z_title))
end

function surface(px, py, pz, option::Int)
  nx = length(px)
  ny = length(py)
  if typeof(pz) == Function
    f = pz
    pz = Float64[f(x,y) for y in py, x in px]
  end
  nz = length(pz)
  if ndims(pz) == 1
    out_of_bounds = nz != nx * ny
  elseif ndims(pz) == 2
    out_of_bounds = size(pz)[1] != nx || size(pz)[2] != ny
  else
    out_of_bounds = true
  end
  if !out_of_bounds
    if ndims(pz) == 2
      pz = reshape(pz, nx * ny)
    end
    ccall( (:gr_surface, libGR),
          Void,
          (Int32, Int32, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Int32),
          nx, ny, convert(Vector{Float64}, px), convert(Vector{Float64}, py), convert(Vector{Float64}, pz), option)
  else
    println("Arrays have incorrect length or dimension.")
  end
end

function contour(px, py, h, pz, major_h::Int)
  nx = length(px)
  ny = length(py)
  nh = length(h)
  if typeof(pz) == Function
    f = pz
    pz = Float64[f(x,y) for y in py, x in px]
  end
  nz = length(pz)
  if ndims(pz) == 1
    out_of_bounds = nz != nx * ny
  elseif ndims(pz) == 2
    out_of_bounds = size(pz)[1] != nx || size(pz)[2] != ny
  else
    out_of_bounds = true
  end
  if !out_of_bounds
    if ndims(pz) == 2
      pz = reshape(pz, nx * ny)
    end
    ccall( (:gr_contour, libGR),
          Void,
          (Int32, Int32, Int32, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Int32),
          nx, ny, nh, convert(Vector{Float64}, px), convert(Vector{Float64}, py), convert(Vector{Float64}, h), convert(Vector{Float64}, pz), major_h)
  else
    println("Arrays have incorrect length or dimension.")
  end
end

function setcolormap(index::Int)
  ccall( (:gr_setcolormap, libGR),
        Void,
        (Int32, ),
        index)
end

function colormap()
  ccall( (:gr_colormap, libGR),
        Void,
        ()
        )
end

function inqcolor(color::Int)
  rgb = Cint[0]
  ccall( (:gr_inqcolor, libGR),
        Void,
        (Int32, Ptr{Int32}),
        color, rgb)
  return rgb
end

function inqcolorfromrgb(red::Real, green::Real, blue::Real)
  color = ccall( (:gr_inqcolorfromrgb, libGR),
                Int32,
                (Float64, Float64, Float64),
                red, green, blue)
  return color
end

function hsvtorgb(h::Real, s::Real, v::Real)
  r = Cdouble[0]
  g = Cdouble[0]
  b = Cdouble[0]
  ccall( (:gr_hsvtorgb, libGR),
        Void,
        (Float64, Float64, Float64, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}),
        h, s, v, r, g, b)
  return r[1], g[1], b[1]
end

function tick(amin::Real, amax::Real)
  return ccall( (:gr_tick, libGR),
               Float64,
               (Float64, Float64),
               amin, amax)
end

function validaterange(amin::Real, amax::Real)
  return ccall( (:gr_validaterange, libGR),
               Int32,
               (Float64, Float64),
               amin, amax)
end

function adjustlimits(amin::Real, amax::Real)
  _amin = Cdouble[amin]
  _amax = Cdouble[amax]
  ccall( (:gr_adjustlimits, libGR),
        Void,
        (Ptr{Float64}, Ptr{Float64}),
        _amin, _amax)
  return _amin[1], _amax[1]
end

function adjustrange(amin::Real, amax::Real)
  _amin = Cdouble[amin]
  _amax = Cdouble[amax]
  ccall( (:gr_adjustrange, libGR),
        Void,
        (Ptr{Float64}, Ptr{Float64}),
        _amin, _amax)
  return _amin[1], _amax[1]
end

function beginprint(pathname)
  ccall( (:gr_beginprint, libGR),
        Void,
        (Ptr{Cchar}, ),
        pathname)
end

function beginprintext(pathname, mode, fmt, orientation)
  ccall( (:gr_beginprintext, libGR),
        Void,
        (Ptr{Cchar}, Ptr{Cchar}, Ptr{Cchar}, Ptr{Cchar}),
        pathname, mode, fmt, orientation)
end

function endprint()
  ccall( (:gr_endprint, libGR),
        Void,
        ()
        )
end

function ndctowc(x::Real, y::Real)
  _x = Cdouble[x]
  _y = Cdouble[y]
  ccall( (:gr_ndctowc, libGR),
        Void,
        (Ptr{Float64}, Ptr{Float64}),
        _x, _y)
  return _x[1], _y[1]
end

function wctondc(x::Real, y::Real)
  _x = Cdouble[x]
  _y = Cdouble[y]
  ccall( (:gr_wctondc, libGR),
        Void,
        (Ptr{Float64}, Ptr{Float64}),
        _x, _y)
  return _x[1], _y[1]
end

function wc3towc(x::Real, y::Real, z::Real)
  _x = Cdouble[x]
  _y = Cdouble[y]
  _z = Cdouble[z]
  ccall( (:gr_wc3towc, libGR),
        Void,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}),
        _x, _y, _z)
  return _x[1], _y[1], _z[1]
end

function drawrect(xmin::Real, xmax::Real, ymin::Real, ymax::Real)
  ccall( (:gr_drawrect, libGR),
        Void,
        (Float64, Float64, Float64, Float64),
        xmin, xmax, ymin, ymax)
end

function fillrect(xmin::Real, xmax::Real, ymin::Real, ymax::Real)
  ccall( (:gr_fillrect, libGR),
        Void,
        (Float64, Float64, Float64, Float64),
        xmin, xmax, ymin, ymax)
end

function drawarc(xmin::Real, xmax::Real, ymin::Real, ymax::Real, a1::Int, a2::Int)
  ccall( (:gr_drawarc, libGR),
        Void,
        (Float64, Float64, Float64, Float64, Int32, Int32),
        xmin, xmax, ymin, ymax, a1, a2)
end

function fillarc(xmin::Real, xmax::Real, ymin::Real, ymax::Real, a1::Int, a2::Int)
  ccall( (:gr_fillarc, libGR),
        Void,
        (Float64, Float64, Float64, Float64, Int32, Int32),
        xmin, xmax, ymin, ymax, a1, a2)
end

function drawpath(points, codes, fill::Int)
  len = length(points)
  ccall( (:gr_drawpath, libGR),
        Void,
        (Int32, Ptr{Float64}, Ptr{UInt8}, Int32),
        len, convert(Vector{Float64}, points), convert(Vector{UInt8}, codes), fill)
end

function setarrowstyle(style::Int)
  ccall( (:gr_setarrowstyle, libGR),
        Void,
        (Int32, ),
        style)
end

function setarrowsize(size::Real)
  ccall( (:gr_setarrowsize, libGR),
        Void,
        (Float64, ),
        size)
end

function drawarrow(x1::Real, y1::Real, x2::Real, y2::Real)
  ccall( (:gr_drawarrow, libGR),
        Void,
        (Float64, Float64, Float64, Float64),
        x1, y1, x2, y2)
end

function readimage(path)
  width = Cint[0]
  height = Cint[0]
  data = Array(Ptr{UInt32}, 1)
  ccall( (:gr_readimage, libGR),
        Void,
        (Ptr{Cchar}, Ptr{Int32}, Ptr{Int32}, Ptr{Ptr{UInt32}}),
        path, width, height, data)
  if width[1] > 0 && height[1] > 0
    img = pointer_to_array(data[1], (width[1], height[1]))
    return Int(width[1]), Int(height[1]), img
  else
    return 0, 0, zeros(UInt32, 0)
  end
end

function drawimage(xmin::Real, xmax::Real, ymin::Real, ymax::Real, width::Int, height::Int, data, model::Int = 0)
  if ndims(data) == 2
    data = reshape(data, width * height)
  end
  ccall( (:gr_drawimage, libGR),
        Void,
        (Float64, Float64, Float64, Float64, Int32, Int32, Ptr{UInt32}, Int32),
        xmin, xmax, ymin, ymax, width, height, convert(Vector{UInt32}, data), model)
end

function importgraphics(path)
  ccall( (:gr_importgraphics, libGR),
        Void,
        (Ptr{Cchar}, ),
        path)
end

function setshadow(offsetx::Real, offsety::Real, blur::Real)
  ccall( (:gr_setshadow, libGR),
        Void,
        (Float64, Float64, Float64),
        offsetx, offsety, blur)
end

function settransparency(alpha::Real)
  ccall( (:gr_settransparency, libGR),
        Void,
        (Float64, ),
        alpha)
end

function setcoordxform(mat)
  assert(length(mat) == 6)
  ccall( (:gr_setcoordxform, libGR),
        Void,
        (Ptr{Float64}, ),
        convert(Vector{Float64}, mat))
end

function begingraphics(path)
  ccall( (:gr_begingraphics, libGR),
        Void,
        (Ptr{Cchar}, ),
        path)
end

function endgraphics()
  ccall( (:gr_endgraphics, libGR),
        Void,
        ()
        )
end

function getgraphics()
  string = ccall( (:gr_getgraphics, libGR),
                 Ptr{Cchar},
                 (),
                 )
  return Compat.String(string)
end

function drawgraphics(string)
  ret = ccall( (:gr_drawgraphics, libGR),
              Int32,
              (Ptr{Cchar}, ),
              string)
  return int(ret)
end

function mathtex(x::Real, y::Real, string)
  ccall( (:gr_mathtex, libGR),
        Void,
        (Float64, Float64, Ptr{Cchar}),
        x, y, string)
end

ASF_BUNDLED = 0
ASF_INDIVIDUAL = 1

NOCLIP = 0
CLIP = 1

COORDINATES_WC = 0
COORDINATES_NDC = 1

INTSTYLE_HOLLOW = 0
INTSTYLE_SOLID = 1
INTSTYLE_PATTERN = 2
INTSTYLE_HATCH = 3

TEXT_HALIGN_NORMAL = 0
TEXT_HALIGN_LEFT = 1
TEXT_HALIGN_CENTER = 2
TEXT_HALIGN_RIGHT = 3
TEXT_VALIGN_NORMAL = 0
TEXT_VALIGN_TOP = 1
TEXT_VALIGN_CAP = 2
TEXT_VALIGN_HALF = 3
TEXT_VALIGN_BASE = 4
TEXT_VALIGN_BOTTOM = 5

TEXT_PATH_RIGHT = 0
TEXT_PATH_LEFT = 1
TEXT_PATH_UP = 2
TEXT_PATH_DOWN = 3

TEXT_PRECISION_STRING = 0
TEXT_PRECISION_CHAR = 1
TEXT_PRECISION_STROKE = 2

LINETYPE_SOLID = 1
LINETYPE_DASHED = 2
LINETYPE_DOTTED = 3
LINETYPE_DASHED_DOTTED = 4
LINETYPE_DASH_2_DOT = -1
LINETYPE_DASH_3_DOT = -2
LINETYPE_LONG_DASH = -3
LINETYPE_LONG_SHORT_DASH = -4
LINETYPE_SPACED_DASH = -5
LINETYPE_SPACED_DOT = -6
LINETYPE_DOUBLE_DOT = -7
LINETYPE_TRIPLE_DOT = -8

MARKERTYPE_DOT = 1
MARKERTYPE_PLUS = 2
MARKERTYPE_ASTERISK = 3
MARKERTYPE_CIRCLE = 4
MARKERTYPE_DIAGONAL_CROSS = 5
MARKERTYPE_SOLID_CIRCLE = -1
MARKERTYPE_TRIANGLE_UP = -2
MARKERTYPE_SOLID_TRI_UP = -3
MARKERTYPE_TRIANGLE_DOWN = -4
MARKERTYPE_SOLID_TRI_DOWN = -5
MARKERTYPE_SQUARE = -6
MARKERTYPE_SOLID_SQUARE = -7
MARKERTYPE_BOWTIE = -8
MARKERTYPE_SOLID_BOWTIE = -9
MARKERTYPE_HOURGLASS = -10
MARKERTYPE_SOLID_HGLASS = -11
MARKERTYPE_DIAMOND = -12
MARKERTYPE_SOLID_DIAMOND = -13
MARKERTYPE_STAR = -14
MARKERTYPE_SOLID_STAR = -15
MARKERTYPE_TRI_UP_DOWN = -16
MARKERTYPE_SOLID_TRI_RIGHT = -17
MARKERTYPE_SOLID_TRI_LEFT = -18
MARKERTYPE_HOLLOW_PLUS = -19
MARKERTYPE_OMARK = -20
MARKERTYPE_PENTAGON = -21
MARKERTYPE_HEXAGON = -22
MARKERTYPE_HEPTAGON = -23
MARKERTYPE_OCTAGON = -24
MARKERTYPE_STAR_4 = -25
MARKERTYPE_STAR_5 = -26
MARKERTYPE_STAR_6 = -27
MARKERTYPE_STAR_7 = -28
MARKERTYPE_STAR_8 = -29
MARKERTYPE_VLINE = -30
MARKERTYPE_HLINE = -31

OPTION_X_LOG = 1
OPTION_Y_LOG = 2
OPTION_Z_LOG = 4
OPTION_FLIP_X = 8
OPTION_FLIP_Y = 16
OPTION_FLIP_Z = 32

OPTION_LINES = 0
OPTION_MESH = 1
OPTION_FILLED_MESH = 2
OPTION_Z_SHADED_MESH = 3
OPTION_COLORED_MESH = 4
OPTION_CELL_ARRAY = 5
OPTION_SHADED_MESH = 6
OPTION_HEATMAP = 7

MODEL_RGB = 0
MODEL_HSV = 1

COLORMAP_UNIFORM = 0
COLORMAP_TEMPERATURE = 1
COLORMAP_GRAYSCALE = 2
COLORMAP_GLOWING = 3
COLORMAP_RAINBOWLIKE = 4
COLORMAP_GEOLOGIC = 5
COLORMAP_GREENSCALE = 6
COLORMAP_CYANSCALE = 7
COLORMAP_BLUESCALE = 8
COLORMAP_MAGENTASCALE = 9
COLORMAP_REDSCALE = 10
COLORMAP_FLAME = 11
COLORMAP_BROWNSCALE = 12
COLORMAP_PILATUS = 13
COLORMAP_AUTUMN = 14
COLORMAP_BONE = 15
COLORMAP_COOL = 16
COLORMAP_COPPER = 17
COLORMAP_GRAY = 18
COLORMAP_HOT = 19
COLORMAP_HSV = 20
COLORMAP_JET = 21
COLORMAP_PINK = 22
COLORMAP_SPECTRAL = 23
COLORMAP_SPRING = 24
COLORMAP_SUMMER = 25
COLORMAP_WINTER = 26
COLORMAP_GIST_EARTH = 27
COLORMAP_GIST_HEAT = 28
COLORMAP_GIST_NCAR = 29
COLORMAP_GIST_RAINBOW = 30
COLORMAP_GIST_STERN = 31
COLORMAP_AFMHOT = 32
COLORMAP_BRG = 33
COLORMAP_BWR = 34
COLORMAP_COOLWARM = 35
COLORMAP_CMRMAP = 36
COLORMAP_CUBEHELIX = 37
COLORMAP_GNUPLOT = 38
COLORMAP_GNUPLOT2 = 39
COLORMAP_OCEAN = 40
COLORMAP_RAINBOW = 41
COLORMAP_SEISMIC = 42
COLORMAP_TERRAIN = 43
COLORMAP_VIRIDIS = 44
COLORMAP_INFERNO = 45
COLORMAP_PLASMA = 46
COLORMAP_MAGMA = 47

FONT_TIMES_ROMAN = 101
FONT_TIMES_ITALIC = 102
FONT_TIMES_BOLD = 103
FONT_TIMES_BOLDITALIC = 104
FONT_HELVETICA = 105
FONT_HELVETICA_OBLIQUE = 106
FONT_HELVETICA_BOLD = 107
FONT_HELVETICA_BOLDOBLIQUE = 108
FONT_COURIER = 109
FONT_COURIER_OBLIQUE = 110
FONT_COURIER_BOLD = 111
FONT_COURIER_BOLDOBLIQUE = 112
FONT_SYMBOL = 113
FONT_BOOKMAN_LIGHT = 114
FONT_BOOKMAN_LIGHTITALIC = 115
FONT_BOOKMAN_DEMI = 116
FONT_BOOKMAN_DEMIITALIC = 117
FONT_NEWCENTURYSCHLBK_ROMAN = 118
FONT_NEWCENTURYSCHLBK_ITALIC = 119
FONT_NEWCENTURYSCHLBK_BOLD = 120
FONT_NEWCENTURYSCHLBK_BOLDITALIC = 121
FONT_AVANTGARDE_BOOK = 122
FONT_AVANTGARDE_BOOKOBLIQUE = 123
FONT_AVANTGARDE_DEMI = 124
FONT_AVANTGARDE_DEMIOBLIQUE = 125
FONT_PALATINO_ROMAN = 126
FONT_PALATINO_ITALIC = 127
FONT_PALATINO_BOLD = 128
FONT_PALATINO_BOLDITALIC = 129
FONT_ZAPFCHANCERY_MEDIUMITALIC = 130
FONT_ZAPFDINGBATS = 131

PATH_STOP      = 0x00
PATH_MOVETO    = 0x01
PATH_LINETO    = 0x02
PATH_CURVE3    = 0x03
PATH_CURVE4    = 0x04
PATH_CLOSEPOLY = 0x4f

MPL_SUPPRESS_CLEAR = 1
MPL_POSTPONE_UPDATE = 2

# GR3 functions
include("gr3.jl")

const gr3 = GR.GR3

# Convenience functions
include("jlgr.jl")

figure(; kwargs...) = jlgr.figure(; kwargs...)
hold(flag) = jlgr.hold(flag)
subplot(m, n, p) = jlgr.subplot(m, n, p)
plot(args...; kwargs...) = jlgr.plot(args...; kwargs...)
oplot(args...; kwargs...) = jlgr.oplot(args...; kwargs...)
semilogx(args...; kwargs...) = jlgr.plot(args...; kwargs..., xlog=true)
semilogy(args...; kwargs...) = jlgr.plot(args...; kwargs..., ylog=true)
loglog(args...; kwargs...) = jlgr.plot(args...; kwargs..., xlog=true, ylog=true)
scatter(args...; kwargs...) = jlgr.scatter(args...; kwargs...)
stem(args...; kwargs...) = jlgr.stem(args...; kwargs...)
histogram(x; kwargs...) = jlgr.histogram(x; kwargs...)
contour(args...; kwargs...) = jlgr.contour(args...; kwargs...)
contourf(args...; kwargs...) = jlgr.contourf(args...; kwargs...)
wireframe(args...; kwargs...) = jlgr.wireframe(args...; kwargs...)
surface(args...; kwargs...) = jlgr.surface(args...; kwargs...)
plot3(args...; kwargs...) = jlgr.plot3(args...; kwargs...)
scatter3(args...; kwargs...) = jlgr.scatter3(args...; kwargs...)
title(s) = jlgr.title(s)
xlabel(s) = jlgr.xlabel(s)
ylabel(s) = jlgr.ylabel(s)
legend(args...; kwargs...) = jlgr.legend(args...; kwargs...)
xlim(a) = jlgr.xlim(a)
ylim(a) = jlgr.ylim(a)
savefig(filename) = jlgr.savefig(filename)
meshgrid(vx, vy) = jlgr.meshgrid(vx, vy)
meshgrid(vx, vy, vz) = jlgr.meshgrid(vx, vy, vz)
peaks(n...) = jlgr.peaks(n...)
imshow(I; kwargs...) = jlgr.imshow(I; kwargs...)
isosurface(V; kwargs...) = jlgr.isosurface(V; kwargs...)
cart2sph(x, y, z) = jlgr.cart2sph(x, y, z)
sph2cart(θ, ϕ, r) = jlgr.sph2cart(θ, ϕ, r)
polar(args...; kwargs...) = jlgr.polar(args...; kwargs...)
trisurf(args...; kwargs...) = jlgr.trisurf(args...; kwargs...)

type SVG
   s::Array{UInt8}
end
@compat Base.show(io::IO, ::MIME"image/svg+xml", x::SVG) = write(io, x.s)

type PNG
   s::Array{UInt8}
end
@compat Base.show(io::IO, ::MIME"image/png", x::PNG) = write(io, x.s)

type HTML
   s::AbstractString
end
@compat Base.show(io::IO, ::MIME"text/html", x::HTML) = print(io, x.s)

function _readfile(path)
    data = Array(UInt8, filesize(path))
    s = open(path, "r")
    Compat.String(read!(s, data))
end

function isinline()
    global mime_type
    return !(mime_type in (None, "mov"))
end

function startserver()
    global msgs

    @eval import WebSockets
    @eval import HttpServer

    msgs = []
    app = WebSockets.WebSocketHandler() do req, client
        while true
            msg = WebSockets.read(client)
            if startswith(Compat.String(msg), "ready")
                if length(msgs) != 0
                    WebSockets.write(client, msgs[1])
                    shift!(msgs)
                else
                    WebSockets.write(client, "busy")
                end
            end
        end
    end

    server = HttpServer.Server(app)
    @async begin
        HttpServer.run(server, 8889)
    end

    if !isfile("gr.js")
        symlink(joinpath(dirname(@__FILE__), "gr.js"), "gr.js")
    end

    return HTML("""\
<canvas id="canvas" width="600" height="450"></canvas>\
<script type="text/javascript" src="gr.js"></script>\
<script>GR.ready(\
  function() {\
    var ws = new WebSocket("ws://localhost:8889/");\
    ws.onopen = function() {\
      ws.send("ready");\
    };\
    ws.onmessage = function(ev) {\
      if (ev.data == "busy") {\
        setTimeout(function() {ws.send("ready");}, 10);\
      } else {\
        gr_clearws();\
        gr_drawgraphics(window.atob(ev.data));\
        ws.send("ready");\
      };\
    };\
  }\
);</script>""")
end

function inline(mime="svg", scroll=true)
    global mime_type, figure_count, msgs
    if mime_type != mime
        if mime == "iterm"
            ENV["GKS_WSTYPE"] = "pdf"
        elseif mime == "mlterm"
            ENV["GKS_WSTYPE"] = "six"
        elseif mime == "atom"
            ENV["GKS_WSTYPE"] = "svg"
            @eval using Atom
            @eval import Atom: Media, PlotPane
        elseif mime == "js"
            ENV["GKS_WSTYPE"] = "nul"
        else
            ENV["GKS_WSTYPE"] = mime
        end
        emergencyclosegks()
        mime_type = mime
        if mime == "js"
            startserver()
        end
    end
    figure_count = scroll ? None : 0
    mime_type
end

function show()
    global mime_type, figure_count, msgs

    emergencyclosegks()
    if mime_type == "svg"
        content = SVG(_readfile("gks.svg"))
        return content
    elseif mime_type == "png"
        content = PNG(_readfile("gks.png"))
        return content
    elseif mime_type == "mov"
        content = HTML(string("""<video autoplay controls><source type="video/mp4" src="data:video/mp4;base64,""", base64encode(open(read,"gks.mov")),""""></video>"""))
        return content
    elseif mime_type == "iterm"
        content = string("\033]1337;File=inline=1;height=24;preserveAspectRatio=0:", base64encode(open(read,"gks.pdf")), "\a")
        if figure_count != None
            figure_count += 1
            (figure_count > 1) && print("\e[24A")
        end
        println(content)
        return nothing
    elseif mime_type == "mlterm"
        content = read("gks.six")
        write(content)
        return nothing
    elseif mime_type == "atom"
        content = Base.HTML(string("""<div style="display: inline-block; background: #ffffff;">""", readstring("gks.svg"), """</div>"""))
        Atom.render(Atom.PlotPane(), content)
        return nothing
    elseif mime_type == "js"
        if msgs != None
            endgraphics()
            push!(msgs, Base.base64encode(getgraphics()))
        end
    end
    return None
end

function setregenflags(flags=0)
  ccall( (:gr_setregenflags, libGR),
        Void,
        (Int32, ),
        flags)
end

function inqregenflags()
  flags = ccall( (:gr_inqregenflags, libGR),
                Int32,
                ()
                )
  return flags
end

function savestate()
  ccall( (:gr_savestate, libGR),
        Void,
        ()
        )
end

function restorestate()
  ccall( (:gr_restorestate, libGR),
        Void,
        ()
        )
end

function selectcontext(context::Int)
  ccall( (:gr_selectcontext, libGR),
        Void,
        (Int32, ),
        context)
end

function uselinespec(linespec)
  return ccall( (:gr_uselinespec, libGR),
               Int32,
               (Ptr{Cchar}, ),
               linespec)
end

function trisurface(x, y, z)
  nx = length(x)
  ny = length(y)
  nz = length(z)
  n = min(nx, ny, nz)
  ccall( (:gr_trisurface, libGR),
        Void,
        (Int32, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}),
        n, convert(Vector{Float64}, x), convert(Vector{Float64}, y), convert(Vector{Float64}, z))
end

function gradient(x, y, z)
  nx = length(x)
  ny = length(y)
  nz = length(z)
  if ndims(z) == 1
    out_of_bounds = nz != nx * ny
  elseif ndims(z) == 2
    out_of_bounds = size(z)[1] != nx || size(z)[2] != ny
  else
    out_of_bounds = true
  end
  if !out_of_bounds
    if ndims(z) == 2
      z = reshape(z, nx * ny)
    end
    u = Cdouble[1 : nx*ny ;]
    v = Cdouble[1 : nx*ny ;]
    ccall( (:gr_gradient, libGR),
          Void,
          (Int32, Int32, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Cdouble}, Ptr{Cdouble}),
          nx, ny, convert(Vector{Float64}, x), convert(Vector{Float64}, y), convert(Vector{Float64}, z), u, v)
    return u, v
  else
    return [], []
  end
end

function quiver(x, y, u, v, color::Bool=false)
  nx = length(x)
  ny = length(y)
  nu = length(u)
  nv = length(v)
  if ndims(u) == 1
    out_of_bounds = nu != nx * ny
  elseif ndims(u) == 2
    out_of_bounds = size(u)[1] != nx || size(u)[2] != ny
  else
    out_of_bounds = true
  end
  if !out_of_bounds
    if ndims(v) == 1
      out_of_bounds = nv != nx * ny
    elseif ndims(v) == 2
      out_of_bounds = size(v)[1] != nx || size(v)[2] != ny
    else
      out_of_bounds = true
    end
  end
  if !out_of_bounds
    if ndims(u) == 2
      u = reshape(u, nx * ny)
    end
    if ndims(v) == 2
      v = reshape(v, nx * ny)
    end
    ccall( (:gr_quiver, libGR),
          Void,
          (Int32, Int32, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Int32),
          nx, ny, convert(Vector{Float64}, x), convert(Vector{Float64}, y), convert(Vector{Float64}, u), convert(Vector{Float64}, v), convert(Int32, color))
  else
    println("Arrays have incorrect length or dimension.")
  end
end

end # module


text:
@testset "Cipher" begin
    @testset "Cipher: atbash" begin
        input = "test"
        @test atbash_encode(input) == "gvhg"
        input = "abcdefghijklmnopqrstuvwxyz"
        @test atbash_encode(input) == "zyxwv utsrq ponml kjihg fedcb a"
        input = "hello"
        @test atbash_encode(input) == "svool"
    end

    @testset "Cipher: caesar" begin
        input = "stealth"
        @test caesar(0, input) == "stealth"
        input = "ghsozhv"
        @test caesar(90, input) == "stealth"

        # rot = 13
        # s = "abcdefghijklmnopqrstuvwxyz"
        # @test rotate(rot, s) == "nopqrstuvwxyzabcdefghijklm"
        # rot = 5
        # s = "omg"
        # @test rotate(rot, s) == "trl"
        # rot = 0
        # s = "hello"
        # @test rotate(rot, s) == "hello"
    end
end

ext: jl

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia


text:
using Knet:MB

import Knet.nll
import Knet.accuracy

"""
`nll(model::KnetModule, data::Knet.MB; o...)`
wrapper over `Knet.nll`.
"""
nll(model::KnetModule, data::MB; o...) =
    nll(model, data, (m, x)->m(actx(), x); o...)

"""
`accuracy(model::KnetModule, data::Knet.MB; o...)`
wrapper over `Knet.accuracy`.
"""
accuracy(model::KnetModule, data::MB; o...) =
    accuracy(model, data, (m, x)->m(actx(), x); o...)


text:
# primary features:
# rollout sequence of controls
# df/dx, df/du derivatives

# features to add:
# address mjdata, mjmodel fields by string or symbol
#   replicate the get_sensors stuff
__precompile__()

module mjWrap

using MuJoCo
using JLD

immutable mjSet
   m::mj.jlModel
   d::mj.jlData
   datas::Vector{mj.jlData}

   # helpers
   sensor_names::Dict{Symbol, Range}
   name::String
   skip::Integer
   mode::Integer

   nq::Integer
   nv::Integer
   nu::Integer
   ns::Integer
   dt::Float64

   #function mjSet(model_name::String) # TODO?
   #end
end

immutable TrajSamples
   T::Integer
   numT::Integer
   s0::Vector{Float64}
   ctrl::Array{Float64,3}
   state::Array{Float64,3}
   obs::Array{Float64,3}
   reward::Array{Float64,2}
end

#=
type TrajDerivatives
   T::Integer
   numT::Integer
   d_x::Array{Float64}
   d_u::Array{Float64}
   s_x::Array{Float64}
   s_u::Array{Float64}

   #c::Array{Float64}
   #c_x::Array{Float64}
   #c_u::Array{Float64}
   #c_xx::Array{Float64}
   #c_uu::Array{Float64}
   #c_ux::Array{Float64}
end
=#

const mjw = mjWrap
export TrajSamples, mjSet, mjw

function finish()
   mj.deactivate()
   return 0;
end

function modelparams(x::mjSet)
   return x.nq, x.nv, x.nu, x.ns
end

function clear_model(x::mjSet)
   mj.deleteModel(x.m.m)
   mj.deleteData(x.d.d)
   for i=1:length(x.datas)
      mj.deleteData(x.datas[i])
   end
end

function load_model(model_name::String, skip::Integer, mode::String="normal", ns::Integer=0)
   pm = mj.loadXML(model_name, C_NULL)
   if pm == nothing
      return
   end
   pd = mj.makeData(pm)

   m, d = mj.mapmujoco(pm, pd)
   ndata = Threads.nthreads()
   datas = Vector{jlData}(ndata)
   print("Making mjDatas: 1 ")
   datas[1] = d
   for i=2:ndata
      print("$i ")
      datas[i] = mj.mapdata(pm, mj.makeData(pm))
   end
   println()

   if mj.get(m, :nsensor) > 0
      sensors = mj.name2range(m, mj.get(m, :nsensor),
                              m.name_sensoradr, m.sensor_adr, m.sensor_dim)
   else
      sensors = Dict(:none=>0:0)
   end

   nq = Int64(mj.get(m, :nq))
   nv = Int64(mj.get(m, :nv))
   nu = Int64(mj.get(m, :nu))
   if ns==0 ns = Int64(mj.get(m, :nsensordata)) end
   dt = mj.get(m, :opt, :timestep)

   return mjSet(m, d, datas, sensors, model_name, skip, 0,
                nq, nv, nu, ns, dt)
end

############ UTILS ############
function wraptopi(ang::Float64)
   return mod(ang+pi, 2pi) - pi
end

function unitfit(ang::Float64, val::Float64)
   clamp.(ang, -val, val)./val
end

function scaledcontrol(m::mj.jlModel, d::mj.jlData,
                       ctrl::AbstractVector{Float64})
   nv = length(d.qvel)
   mark = mj.MARKSTACK(d)

   My = @view d.stack[mark+1:(mark+nv)]
   mj.solveM(m.m, d.d, My, d.qfrc_actuator, nv)

   force = d.qfrc_actuator'*My

   mj.FREESTACK(d, Int32(mark))
   return force
end

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

function reset(m::mj.jlModel, d::mj.jlData,
               s0::AbstractVector{Float64},
               ctrl0::AbstractVector{Float64},
               obs::AbstractVector{Float64})
   # setup
   nq = length(d.qpos) #mj.get(m, :nq)
   nv = length(d.qvel) #mj.get(m, :nv)

   # reset
   d.qacc .= 0.0
   for i=1:nq
      d.qpos[i] = s0[i]
   end
   for i=1:nv
      d.qvel[i] = s0[nq+i]
   end
   copy!(d.ctrl, ctrl0)
   mj.set(d, :time, 0.0)

   # out
   mj.forward(m, d)
   copy!(obs, d.sensordata)
end

# accepts threading; tid is 1 indexed in Julia
function reset(x::mjSet, tid::Integer,
               s0::AbstractVector{Float64},
               ctrl0::AbstractVector{Float64},
               obs::AbstractVector{Float64})

   reset(x.m, x.datas[tid], s0, ctrl0, obs)
end

function reset(x::mjSet,
               s0::AbstractVector{Float64},
               ctrl0::AbstractVector{Float64},
               obs::AbstractVector{Float64})
   reset(x.m, x.d, s0, ctrl0, obs)
end

########################################################################## step

function quat2euler(quat::Vector{Float64})
   w, x, y, z = quat
   sinr = 2.0 * (w * x + y * z)
   cosr = 1.0 - 2.0 * (x * x + y * y)
   roll = atan2(sinr, cosr)

   # pitch (y-axis rotation)
   sinp = 2.0 * (w * y - z * x)
   if abs(sinp) >= 1
      pitch = copysign(pi / 2, sinp) # use 90 degrees if out of range
   else
      pitch = asin(sinp)
   end

   # yaw (z-axis rotation)
   siny = 2.0 * (w * z + x * y)
   cosy = 1.0 - 2.0 * (y * y + z * z)
   yaw = atan2(siny, cosy)

   return roll, pitch, yaw
end

function euler2quat(roll::Float64, pitch::Float64, yaw::Float64)
   # blame wikipedia if this doesn't work
   cy = cos(yaw * 0.5)
   sy = sin(yaw * 0.5)
   cr = cos(roll * 0.5)
   sr = sin(roll * 0.5)
   local cp = cos(pitch * 0.5)
   sp = sin(pitch * 0.5)

   w = cy * cr * cp + sy * sr * sp
   x = cy * sr * cp - sy * cr * sp
   y = cy * cr * sp + sy * sr * cp
   z = sy * cr * cp - cy * sr * sp
   return [w, x, y, z]
end

function step!(m::mj.jlModel, d::mj.jlData, skip::Integer,
               ctrl::AbstractVector{Float64},
               s1::AbstractVector{Float64},
               obs::AbstractVector{Float64})
   nq = length(d.qpos)
   copy!(d.ctrl, ctrl)

   for i=1:skip
      mj.step(m, d) # TODO MAKE FASTER with fwd and euler
   end

   mj.forward(m, d)
   s1[1:nq]     = d.qpos
   s1[nq+1:end] = d.qvel
   copy!(obs, d.sensordata)
end

# accepts threading; tid is 1 indexed in Julia
function step!(x::mjSet, tid::Integer,
               ctrl::AbstractVector{Float64},
               s1::AbstractVector{Float64},
               obs::AbstractVector{Float64})
   step!(x.m, x.datas[tid], x.skip, ctrl, s1, obs)
end

# resets position then takes step
function step!(x::mjSet, tid::Integer,
               s0::AbstractVector{Float64},
               ctrl::AbstractVector{Float64},
               s1::AbstractVector{Float64},
               obs::AbstractVector{Float64})
   error("init_step s0 -> s1 not implemented")
end

function step!(x::mjSet,
               ctrl::Array{Float64},
               s1::Array{Float64},
               obs::Array{Float64})
   step!(x.m, x.d, x.skip, ctrl, s1, obs)
end

function alloc_fields(nqnv::Integer, nu::Integer, ns::Integer,
                      T::Integer, numT::Integer)
   ctrl  = zeros(nu, T, numT)
   state = zeros(nqnv, T, numT) # just allocate space
   obs   = zeros(ns,   T, numT)
   return ctrl, state, obs
end

# optionally we can specify what we want our observations to be elsewhere
# instead of in mujoco's sensors format
function allocateTrajSamples(x::mjSet, T::Integer, numT::Integer,
                             ns::Integer=x.ns)
   nq, nv, nu, _ = mjw.modelparams(x)
   ctrl, state, obs = alloc_fields(nq+nv, nu, ns, T, numT)

   s0 = zeros(nq+nv)
   s0[1:nq]     = x.d.qpos
   s0[nq+1:end] = x.d.qvel

   return TrajSamples(T, numT, s0,
                      ctrl, state, obs, zeros(T,numT))
end

#=
function alloc_derivs(m::Model, T::Integer, numT::Integer)
   # to note: this is the format c should expect
   d_x = Array(Float64, m.nqnv, m.nqnv, T, numT)
   d_u = Array(Float64, m.nqnv, m.nu,   T, numT)
   s_x = Array(Float64, m.ns,   m.nqnv, T, numT)
   s_u = Array(Float64, m.ns,   m.nu,   T, numT)
   return TrajDerivatives(T, numT, d_x, d_u, s_x, s_u)
end
=#

function save_traj(x::mjSet,
                   s::TrajSamples,
                   filename::String)
   save(filename,
        "state", s.state,
        "ctrl", s.ctrl,
        "model", basename(x.name),
        "skip", x.skip,
        "mode", "normal")
end

#function rolloutstep(x::mjw.mjSet,
#                 samp::mjw.TrajSamples,
#                 k::Int, t::int,

function rollout(x::mjw.mjSet,
                 samp::mjw.TrajSamples,
                 ctrlfunc!::Function=(y...)->nothing,
                 obsfunc!::Function=(y...)->nothing,
                 rewardfunc::Function=(y...)->0.0)
   nthread = Int(min(samp.numT, Int(Threads.nthreads()))) # cant have more threads than data
   BLAS.set_num_threads(1)
   nq, nv, ns, nu = x.nq, x.nv, x.ns, x.nu
   Threads.@threads for tid=1:nthread
      #for tid=1:nthread
      s0, s = Array{Float64}(nq+nv), Array{Float64}(nq+nv)
      o0, o = Array{Float64}(ns), Array{Float64}(ns)
      c     = Array{Float64}(nu)
      thread_range = Distributed.splitrange(samp.numT, nthread)[tid]

      for t=Range(thread_range)
         s[:] = samp.state[:,1,t]
         c[:] = samp.ctrl[:,1,t]
         reset(x, tid, s, c, o)
         obsfunc!(x, tid, s, o)    # arbitrary observation vector manipulation
         samp.state[:,1,t] = s
         s0[:]             = s
         samp.obs[:,1,t]   = o
         o0[:]             = o
         for k=Range(1:samp.T)

            c .= samp.ctrl[:,k,t]  # if ctrlfunc does nothing (no controller)
            # then we will execute the stored control
            ctrlfunc!(c, o)        # we can also set noise to be in ctrl

            step!(x, tid, c, s, o) # get next state, next observation
            obsfunc!(x, tid, s, o) # arbitrary observation vector manipulation

            samp.reward[k,t] = rewardfunc(x, tid,
                                          s0, s,
                                          o0, o, c)
            s0[:] = s
            o0[:] = o
            samp.ctrl[:,k,t] = c
            if k < samp.T
               samp.state[:,k+1,t] = s
               samp.obs[:,k+1,t]   = o
            end
         end
      end
   end
   BLAS.set_num_threads(Threads.nthreads())
end

#=
function rollout(x::mjw.mjSet,
   samp::mjw.TrajSamples,
   obsfunc!::Function=(y...)->nothing,
   rewardfunc!::Function=(y...)->nothing)
   nthread = min(samp.numT, Threads.nthreads()) # cant have more threads than data
   #Threads.@threads for tid=1:nthread
      for tid=1:nthread
         s0, s = Array{Float64}(x.nq+x.nv), Array{Float64}(x.nq+x.nv)
         o0, o = Array{Float64}(x.ns), Array{Float64}(x.ns)
         c     = Array{Float64}(x.nu)
         thread_range = Distributed.splitrange(samp.numT, nthread)[tid]

         for t=thread_range
            #state= view(samp.state, :, :, t)
            #ctrl = view(samp.ctrl, :, :, t)
            #obs  = view(samp.obs, :, :, t)
            s[:] = samp.state[:,1,t]
            c[:] = samp.ctrl[:,1,t]
            mjw.reset(x, tid, s, c, o)
            obsfunc!(x, s, o) # arbitrary observation vector manipulation
            samp.state[:,1,t] = s
            samp.obs[:,1,t]   = o
            for k=1:samp.T

               c[:] = samp.ctrl[:,k,t]
               mjw.step!(x, tid, c, s, o) # get next state, next observation
               obsfunc!(x, s, o) # arbitrary observation vector manipulation

               samp.reward[k,t] = rewardfunc(x, tid,
               s0, s,
               o0, o, c)
               s0[:] = s
               o0[:] = o
               if k < samp.T
                  samp.state[:,k+1,t] = s
                  samp.obs[:,k+1,t]   = o
               end
            end
         end
      end
   end
   =#

end # mj module


text:
using ContinuumArrays, QuasiArrays, LazyArrays, IntervalSets, FillArrays, LinearAlgebra, BandedMatrices, ForwardDiff, Test
import ContinuumArrays: ℵ₁, materialize, SimplifyStyle, AffineQuasiVector, BasisLayout, AdjointBasisLayout, SubBasisLayout, MappedBasisLayout
import QuasiArrays: SubQuasiArray, MulQuasiMatrix, Vec, Inclusion, QuasiDiagonal, LazyQuasiArrayApplyStyle, LazyQuasiArrayStyle
import LazyArrays: MemoryLayout, ApplyStyle, Applied, colsupport, arguments, ApplyLayout
import ForwardDiff: Dual

@testset "Inclusion" begin
    x = Inclusion(-1..1)
    @test x[0.1] === 0.1
    @test x[0.0] === 0.0
    x = Inclusion(-1.0..1)
    X = QuasiDiagonal(x)
    @test X[-1:0.1:1,-1:0.1:1] == Diagonal(-1:0.1:1)
    @test Base.BroadcastStyle(typeof(x)) == LazyQuasiArrayStyle{1}()
    @test x .* x isa BroadcastQuasiArray
    @test (x.*x)[0.1] == 0.1^2
end

@testset "DiracDelta" begin
    δ = DiracDelta(-1..3)
    @test axes(δ) === (axes(δ,1),) === (Inclusion(-1..3),)
    @test size(δ) === (length(δ),) === (ℵ₁,)
    @test δ[1.1] === 0.0
    @test δ[0.0] === Inf
    @test Base.IndexStyle(δ) === Base.IndexLinear()
end

@testset "HeavisideSpline" begin
    H = HeavisideSpline([1,2,3])

    @test axes(H) === (axes(H,1),axes(H,2)) === (Inclusion(1..3), Base.OneTo(2))
    @test size(H) === (size(H,1),size(H,2)) === (ℵ₁, 2)

    @test_throws BoundsError H[0.1, 1]
    @test H[1.1,1] === H'[1,1.1] === transpose(H)[1,1.1] === 1.0
    @test H[2.1,1] === H'[1,2.1] === transpose(H)[1,2.1] === 0.0
    @test H[1.1,2] === 0.0
    @test H[2.1,2] === 1.0
    @test_throws BoundsError H[2.1,3]
    @test_throws BoundsError H'[3,2.1]
    @test_throws BoundsError transpose(H)[3,2.1]
    @test_throws BoundsError H[3.1,2]

    @test all(H[[1.1,2.1], 1] .=== H'[1,[1.1,2.1]] .=== transpose(H)[1,[1.1,2.1]] .=== [1.0,0.0])
    @test all(H[1.1,1:2] .=== H[1.1,:] .=== [1.0,0.0])
    @test all(H[[1.1,2.1], 1:2] .=== [1.0 0.0; 0.0 1.0])

    @test_throws BoundsError H[[0.1,2.1], 1]
    @test MemoryLayout(typeof(H)) == BasisLayout()
    @test ApplyStyle(*, typeof(H), typeof([1,2])) isa LazyQuasiArrayApplyStyle
    f = H*[1,2]
    @test axes(f) == (Inclusion(1.0..3.0),)
    @test f[1.1] ≈ 1
    @test f[2.1] ≈ 2

    @test ApplyStyle(*,typeof(H'),typeof(H)) == SimplifyStyle()

    @test H'H == materialize(applied(*,H',H)) == Eye(2)
end

@testset "LinearSpline" begin
    L = LinearSpline([1,2,3])
    @test size(L) == (ℵ₁, 3)

    @test_throws BoundsError L[0.1, 1]
    @test L[1.1,1] == L'[1,1.1] == transpose(L)[1,1.1] ≈ 0.9
    @test L[2.1,1] === L'[1,2.1] === transpose(L)[1,2.1] === 0.0
    @test L[1.1,2] ≈ 0.1
    @test L[2.1,2] ≈ 0.9
    @test L[2.1,3] == L'[3,2.1] == transpose(L)[3,2.1] ≈ 0.1
    @test_throws BoundsError L[3.1,2]

    @test L[[1.1,2.1], 1] == L'[1,[1.1,2.1]] == transpose(L)[1,[1.1,2.1]] ≈ [0.9,0.0]
    @test L[1.1,1:2] ≈ [0.9,0.1]
    @test L[[1.1,2.1], 1:2] ≈ [0.9 0.1; 0.0 0.9]

    @test_throws BoundsError L[[0.1,2.1], 1]

    f = L*[1,2,4]
    @test axes(f) == (Inclusion(1.0..3.0),)
    @test f[1.1] ≈ 1.1
    @test f[2.1] ≈ 2.2

    δ = DiracDelta(1.2,1..3)
    L = LinearSpline([1,2,3])
    @test δ'L ≈ [0.8, 0.2, 0.0]

    @test L'L == SymTridiagonal([1/3,2/3,1/3], [1/6,1/6])

    @testset "==" begin
        L = LinearSpline([1,2,3])
        H = HeavisideSpline([1,2,3])
        @test L == L
        @test L ≠ H
        H = HeavisideSpline([1,1.5,2.5,3])
        @test_throws ArgumentError L == H
    end

    @testset "Adjoint layout" begin
        L = LinearSpline([1,2,3])
        @test MemoryLayout(typeof(L')) == AdjointBasisLayout()
        @test [3,4,5]'*L' isa ApplyQuasiArray
    end

    @testset "Broadcast layout" begin
        L = LinearSpline([1,2,3])
        b = BroadcastQuasiArray(+, L*[3,4,5], L*[1.,2,3])
        @test (L\b) == [4,6,8]
        B = BroadcastQuasiArray(+, L, L)
        @test L\B == 2Eye(3)

        b = BroadcastQuasiArray(-, L*[3,4,5], L*[1.,2,3])
        @test (L\b) == [2,2,2]
        B = BroadcastQuasiArray(-, L, L)
        @test L\B == 0Eye(3)
    end
end

@testset "Derivative" begin
    L = LinearSpline([1,2,3])
    f = L*[1,2,4]
    @test f[1.2] == 1.2

    D = Derivative(axes(L,1))
    @test ApplyStyle(*,typeof(D),typeof(L)) isa SimplifyStyle
    @test D*L isa MulQuasiMatrix
    @test length((D*L).args) == 2
    @test eltype(D*L) == Float64

    M = applied(*, (D*L).args..., [1,2,4])
    @test M isa Applied{LazyQuasiArrayApplyStyle}
    @test eltype(materialize(M)) == Float64

    M = applied(*, D, L, [1,2,4])
    @test M isa Applied{LazyQuasiArrayApplyStyle}

    fp = D*L*[1,2,4]

    @test eltype(fp) == Float64

    @test fp isa Vec
    @test length(fp.args) == 2
    @test fp[1.1] ≈ 1
    @test fp[2.2] ≈ 2

    fp = D*f
    @test length(fp.args) == 2
    @test fp[1.1] ≈ 1
    @test fp[2.2] ≈ 2
end

@testset "Weak Laplacian" begin
    H = HeavisideSpline(0:2)
    L = LinearSpline(0:2)
    D = Derivative(axes(L,1))

    M = QuasiArrays.flatten(Mul(D',D*L))
    @test length(M.args) == 3
    @test last(M.args) isa BandedMatrix

    @test ApplyStyle(*, typeof(L'), typeof(D')) == SimplifyStyle()
    @test apply(*,L',D') isa MulQuasiMatrix
    @test MemoryLayout(typeof(L')) isa AdjointBasisLayout
    @test (L'D') isa MulQuasiMatrix

    A = (L'D') * (D*L)
    @test A isa BandedMatrix
    @test A == (D*L)'*(D*L) == [1.0 -1 0; -1.0 2.0 -1.0; 0.0 -1.0 1.0]
    @test bandwidths(A) == (1,1)
end

@testset "Views" begin
    L = LinearSpline(0:2)
    @test view(L,0.1,1)[1] == L[0.1,1]

    L = LinearSpline(0:2)
    B1 = view(L,:,1)
    @test B1 isa SubQuasiArray{Float64,1}
    @test size(B1) == (ℵ₁,)
    @test B1[0.1] == L[0.1,1]
    @test_throws BoundsError B1[2.2]

    B = view(L,:,1:2)
    @test B isa SubQuasiArray{Float64,2}
    @test B[0.1,:] == L[0.1,1:2]

    B = @view L[:,2:end-1]
    @test B[0.1,:] == [0.1]

    L = LinearSpline([1,2,3,4])
    @test L[:,2:3] isa SubQuasiArray
    @test axes(L[:,2:3]) ≡ (Inclusion(1..4), Base.OneTo(2))
    @test L[:,2:3][1.1,1] == L[1.1,2]
    @test_throws BoundsError L[0.1,1]
    @test_throws BoundsError L[1.1,0]

    @test MemoryLayout(typeof(L[:,2:3])) isa SubBasisLayout
    @test L\L[:,2:3] isa BandedMatrix
    @test L\L[:,2:3] == [0 0; 1 0; 0 1.0; 0 0]

    @testset "Subindex of splines" begin
        L = LinearSpline(range(0,stop=1,length=10))
        @test L[:,2:end-1] isa SubQuasiArray
        @test L[:,2:end-1][0.1,1] == L[0.1,2]
        v = randn(8)
        f = L[:,2:end-1] * v
        @test f[0.1] ≈ (L*[0; v; 0])[0.1]
    end
end

@testset "Poisson" begin
    L = LinearSpline(range(0,stop=1,length=10))
    B = L[:,2:end-1] # Zero dirichlet by dropping first and last spline
    @test (B'B) == (L'L)[2:end-1,2:end-1]
    D = Derivative(axes(L,1))
    @test apply(*,D,B) isa SubQuasiArray
    @test D*B isa MulQuasiMatrix
    @test apply(*,D,B)[0.1,1] == (D*B)[0.1,1] == 9
    @test apply(*,D,B)[0.2,1] == (D*B)[0.2,1] == -9
    Δ = -((D*B)'*(D*B)) # Weak Laplacian
    @test Δ isa BandedMatrix

    @test B'D' isa MulQuasiMatrix
    @test length((B'D').args) == 2

    @test Δ == -(*(B',D',D,B))
    @test Δ == -(B'D'D*B)
    @test Δ == -((B'D')*(D*B))
    @test_broken Δ == -B'*(D'D)*B
    @test Δ == -(B'*(D'D)*B)

    f = L*exp.(L.points) # project exp(x)
    u = B * (Δ \ (B'f))

    @test u[0.1] ≈ -0.06612902692412974
end

@testset "Helmholtz" begin
    L = LinearSpline(range(0,stop=1,length=10))
    B = L[:,2:end-1] # Zero dirichlet by dropping first and last spline
    D = Derivative(axes(L,1))

    A = -((B'D')*(D*B)) + 100^2*B'B # Weak Laplacian

    f = L*exp.(L.points) # project exp(x)
    u = B * (A \ (B'f))

    @test u[0.1] ≈ 0.00012678835289369413
end

@testset "AffineQuasiVector" begin
    x = Inclusion(0..1)
    @test 2x isa AffineQuasiVector
    @test (2x)[0.1] == 0.2
    @test_throws BoundsError (2x)[2]
    y = 2x .- 1
    @test y isa AffineQuasiVector
    @test y[0.1] == 2*(0.1)-1
    @test y/2 isa AffineQuasiVector
    @test 2\y isa AffineQuasiVector
    @test (y/2)[0.1] == (2\y)[0.1] == -0.4
    @test y .+ 1 isa AffineQuasiVector
    @test (y .+ 1)[0.1] == (1 .+ y)[0.1]
    @test (y .- 1)[0.1] == y[0.1]-1
    @test (1 .- y)[0.1] == 1-y[0.1]

    @test findfirst(isequal(0.1),x) == findlast(isequal(0.1),x) == 0.1
    @test findall(isequal(0.1), x) == [0.1]
    @test findfirst(isequal(2),x) == findlast(isequal(2),x) == nothing
    @test findall(isequal(2), x) == Float64[]

    @test findfirst(isequal(0.2),y) == findlast(isequal(0.2),y) == 0.6
    @test findfirst(isequal(2.3),y) == findlast(isequal(2.3),y) == nothing
    @test findall(isequal(0.2),y) == [0.6]
    @test findall(isequal(2),y) == Float64[]

    @test AffineQuasiVector(x)[0.1] == 0.1
    @test minimum(y) == -1
    @test maximum(y) == 1
    @test union(y) == Inclusion(-1..1)
    @test ContinuumArrays.inbounds_getindex(y,0.1) == y[0.1]
    @test ContinuumArrays.inbounds_getindex(y,2.1) == 2*2.1 - 1

    z = 1 .- x
    @test minimum(z) == 0.0
    @test maximum(z) == 1.0
    @test union(z) == Inclusion(0..1)

    @test !isempty(z)
    @test z == z
end

@testset "Change-of-variables" begin
    x = Inclusion(0..1)
    y = 2x .- 1
    L = LinearSpline(range(-1,stop=1,length=10))
    @test L[y,:][0.1,:] == L[2*0.1-1,:]

    D = Derivative(axes(L,1))
    H = HeavisideSpline(L.points)
    @test H\((D*L) * 2) ≈ (H\(D*L))*2 ≈ diagm(0 => fill(-9,9), 1 => fill(9,9))[1:end-1,:]

    @test MemoryLayout(typeof(L[y,:])) isa MappedBasisLayout
    a,b = arguments((D*L)[y,:])
    @test H[y,:]\a == Eye(9)
    @test H[y,:] \ (D*L)[y,:] isa BandedMatrix

    D = Derivative(x)
    @test (D*L[y,:])[0.1,1] ≈ -9
    @test H[y,:] \ (D*L[y,:]) isa BandedMatrix
    @test H[y,:] \ (D*L[y,:]) ≈ diagm(0 => fill(-9,9), 1 => fill(9,9))[1:end-1,:]

    B = L[y,2:end-1]
    @test MemoryLayout(typeof(B)) isa MappedBasisLayout
    @test B[0.1,1] == L[2*0.1-1,2]
    @test B\B == Eye(8)
    @test L[y,:] \ B == Eye(10)[:,2:end-1]
end


text:
# DiagnosticIDs
DiagnosticIDs() = DiagnosticIDs(create_diagnostic_ids())

"""
    create_diagnostic_ids() -> CXDiagnosticIDs
Return a pointer to a `clang::DiagnosticIDs` object.
"""
function create_diagnostic_ids()
    status = Ref{CXInit_Error}(CXInit_NoError)
    ids = clang_DiagnosticIDs_create(status)
    @assert status[] == CXInit_NoError
    return ids
end


text:
@testset "systems.jl: Components" begin
    chainA = homopolymer_chain()
    @test length(chainA.blocks) == 1

    chainAB = diblock_chain()
    @test length(chainAB.blocks) == 2
    @test chainAB.blocks[1].E2 == chainAB.blocks[2].E2
    @test chainAB.blocks[1].E1 != chainAB.blocks[2].E1

    sol = solvent()
    @test sol.label == :S
end

@testset "systems.jl: Systems" begin
    ab = AB_system()
    @test systemtype(ab) == NeatPolymer()
    @test multicomponent(ab) == false
    @test ncomponents(ab) == 1
    @test species(ab) == [:A, :B]
    @test nspecies(ab) == 2

    aba = AB_A_system()
    @test systemtype(aba) == PolymerBlend()
    @test multicomponent(aba) == true
    @test ncomponents(aba) == 2
    @test species(aba) == [:A, :B]
    @test nspecies(aba) == 2

    abs = AB_S_system()
    @test systemtype(abs) == PolymerSolution()
    @test multicomponent(abs) == true
    @test ncomponents(abs) == 2
    @test species(abs) == [:A, :B, :S]
    @test nspecies(abs) == 3
end


text:
module LoadFromMonai

using PythonCall

function permuteAndReverseFromMonai(pixels)
    sizz=size(pixels)
    for i in 1:sizz[2]
        for j in 1:sizz[3]
            pixels[:,i,j] =  reverse(pixels[:,i,j])
        end#
    end#
    return pixels
  end#permuteAndReverse

"""
given file paths it loads
imagePath - path to main image
labelPath - path to label
transforms - monai.transforms.Compose object -important Load imaged should not be in this list it is added separately
default transforms standardize orientation, voxel dimensions crops unnecessary background reducing array size
ensure type of the images and labels so it will be easily convertible to for example numpy and then Julia
more in https://docs.monai.io/en/stable/transforms.html
"""
function loadByMonaiFromImageAndLabelPaths(
    imagePath
    ,labelPath
    ,trAnsforms=[])
    monai=pyimport("monai")
    #default transforms
    if(length(trAnsforms)==0)
        trAnsforms= [
        #monai.transforms.LoadImaged(keys=["image", "label"]),
        monai.transforms.EnsureChannelFirstd(keys=["image", "label"]),
        monai.transforms.Spacingd(keys=["image", "label"], pixdim=(  1.0, 1.0, 1.0), mode=("bilinear", "nearest")),
            monai.transforms.Orientationd(keys=["image", "label"], axcodes="RAS"),
            monai.transforms.CropForegroundd(keys=["image", "label"], source_key="image"),
            monai.transforms.EnsureTyped(keys=["image", "label"])
        ]
    end
    trAnsformsComposed= monai.transforms.Compose(trAnsforms)

    dicttt= pydict(Dict([("image", imagePath),( "label",  labelPath )]))

    loadObj=monai.transforms.LoadImaged(keys=["image", "label"])(dicttt)
    metaData= pyconvert(Dict,pyconvert(Dict,loadObj)["image_meta_dict"])

    loadObj = trAnsformsComposed(loadObj)

    image = permuteAndReverseFromMonai(pyconvert(Array,loadObj["image"].detach().numpy()[0]))
    label =permuteAndReverseFromMonai(pyconvert(Array,loadObj["label"].detach().numpy()[0]))

return (image,label,metaData)

end

end


text:
"""
Overwrite routing link type to allow classification of routing links by class.
"""
struct TypedRoutingLink <: RoutingLink
    channels::Vector{ChannelIndex}
    cost::Float64
    capacity::Int64
    class::String
end

"""
Provide constructor.
"""
function TypedRoutingLink(cost, capacity, class)
    return TypedRoutingLink(Int[], cost, capacity, class)
end

"""
New RoutingTask type to allow classification of tasks by a class.
"""
struct TypedRoutingChannel <: RoutingChannel
    start_vertices::Vector{PortVertices}
    stop_vertices::Vector{PortVertices}
    class::String

    function TypedRoutingChannel(start, stop, edge)
        class = get(edge.metadata, "class", "all")
        return new(start, stop, class)
    end
end

function Mapper2.routing_channel(::TestRuleSet, start, stop, edge::TaskgraphEdge)
    return TypedRoutingChannel(start, stop, edge)
end
function Mapper2.getcapacity(::TestRuleSet, item)
    return get(item.metadata, "capacity", 10)
end

# Provide custom annotation methods.
function Mapper2.annotate(ruleset::TestRuleSet, port::Port)
    capacity = getcapacity(ruleset, port)
    class = get(port.metadata, "class", "all")
    cost = get(port.metadata, "cost", 1.0)
    return TypedRoutingLink(cost, capacity, class)
end
function Mapper2.annotate(ruleset::TestRuleSet, link::Link)
    capacity = getcapacity(ruleset, link)
    class = get(link.metadata, "class", "all")
    cost = get(link.metadata, "cost", 1.0)
    return TypedRoutingLink(cost, capacity, class)
end
function Mapper2.annotate(::TestRuleSet, component::Component)
    @assert component.primitive == "mux"
    return TypedRoutingLink(1.0, 10, "all")
end

# Validity checks
function check_class(item, edge::TaskgraphEdge)
    item_class = get(item.metadata, "class", "all")
    edge_class = get(edge.metadata, "class", "all")
    return check_class(item_class, edge_class)
end
function check_class(item, edge)
    if edge == "all"
        return true
    elseif item == "all"
        return true
    else
        return edge == item
    end
end

function Mapper2.canuse(::TestRuleSet, item::Union{Port,Link}, edge::TaskgraphEdge)
    return check_class(item, edge)
end
function Mapper2.canuse(::TestRuleSet, item::TypedRoutingLink, edge::TypedRoutingChannel)
    return check_class(item.class, edge.class)
end

function Mapper2.is_source_port(::TestRuleSet, port::Port, edge::TaskgraphEdge)
    return check_class(port, edge)
end
function Mapper2.is_sink_port(::TestRuleSet, port::Port, edge::TaskgraphEdge)
    return check_class(port, edge)
end


text:
# module Pointwise

using ..Var
export pw, l, dl, ₚ, PwVar

"""
Pointwise application.

Pointwise function application gives meaning to expressions such as `x + y`  when `x` and `y` are functions.
That is `x + y` is the function `ω -> x(ω) + y(ω)`.

An argument can be either __lifted__ or __not lifted__.
For example in `x = 1 ~  Normal(0, 1); y = pw(+, x, 3)`, `x` will be lifted but `3` will not be in the sense that
`y` will resolve to `ω -> x(ω) + 3` and not `ω -> x(ω) + 3(ω)`.

`pw` uses some reasonable defaults for what to lift and what not to lift, but to have more explicit control use `l` and `dl` to
lift and dont lift respectively.

Example:
```
using Distributions
x(ω) = Uniform(ω, 0, 1)
y = pw(+, x, 4)

f(ω) = Bool(Bernoulli(ω, 0.5)) ? sqrt : sin
sample(pw(map, f, [0, 1, 2]))
sample(pw(map, sqrt, [0, 1, 2])) # Will error!
sample(pw(map, dl(sqrt), [0, 1, 2]))
sample(pw(l(f), 3))
```
"""

struct PwVar{ARGS, D}
  f::D
  args::ARGS
  PwVar(f::F, args::A) where {F, A} = new{A, F}(f, args)
  PwVar(f::Type{T}, args::A) where {T, A} = new{A, Type{T}}(f, args)
end

pw(f) = (args...) -> PwVar(f, args)
pw(f, arg, args...) = PwVar(f, (arg, args...))
Base.show(io::IO, p::PwVar) = print(io, p.f, "ₚ", p.args)

abstract type ABox end

struct LiftBox{T} <: ABox
  val::T
end
"`l(x)` constructs object that indicates that `x` should be applied pointwise.  See `pw`"
l(x) = LiftBox(x)

struct DontLiftBox{T} <: ABox
  val::T
end
"`dl(x)` constructs object that indicates that `x` should be not applied pointwise.  See `pw`"
dl(x) = DontLiftBox(x)

# Traits
struct Lift end
struct DontLift end

# Trait functions
traitlift(::Type{T}) where T  = DontLift()
traitlift(::Type{<:Function}) = Lift()
traitlift(::Type{<:Variable}) = Lift()
traitlift(::Type{<:Member}) = Lift()
traitlift(::Type{<:DataType}) = DontLift()
traitlift(::Type{<:LiftBox}) = Lift()
traitlift(::Type{<:PwVar}) = Lift()
traitlift(::Type{<:DontLiftBox}) = DontLift()

@inline liftapply(f::T, ω) where T = liftapply(traitlift(T), f, ω)
@inline liftapply(::DontLift, f, ω) = f
@inline liftapply(::DontLift, f::ABox, ω) = f.val
@inline liftapply(::Lift, f, ω) = f(ω)
@inline liftapply(::Lift, f::ABox, ω) = (f.val)(ω)
@inline liftapply(::Lift, f::ABox, ω::ABox) = (f.val)(ω.val)
@inline liftapply(::Lift, f, ω::ABox) = f(ω.val)

recurse(p::PwVar{Tuple{T1}}, ω) where {T1} =
  p.f(liftapply(p.args[1], ω))
# recurse(p::PwVar{Tuple{T1}})(id, ω) where {T1} =
#   p.f(liftapply(p.args[1], ω))(id, ω)

recurse(p::PwVar{Tuple{T1, T2}}, ω) where {T1, T2} =
  p.f(liftapply(p.args[1], ω), liftapply(p.args[2], ω))
# recurse(p::PwVar{Tuple{T1, T2}})(id, ω) where {T1, T2} =
#   p.f(liftapply(p.args[1], ω), liftapply(p.args[2], ω))(id, ω)

# Normal(x(ω), y(ω))(id, ω)

recurse(p::PwVar{<:Tuple}, ω) =
  p.f((liftapply(arg, ω) for arg in p.args)...)
# recurse(p::PwVar{<:Tuple})(id, ω) =
#   p.f(id, (liftapply(arg, ω) for arg in p.args)...)

# # Notation

# # Collections
# @inline randcollection(xs) = ω -> 32(x -> liftapply(x, ω), xs)
# struct LiftConst end
# const ₚ = LiftConst()
# Base.:*(xs, ::LiftConst) = randcollection(xs)

# end


text:
global CURRENTID = -1
getGlobalID() = (global CURRENTID -= 1; return CURRENTID + 1)
resetGlobalID() = (global CURRENTID = -1; return)

function global_id!(nodes)
    oldnewid = Dict([node.id => i for (i, node) in enumerate(nodes)]...)
    for node in nodes
        node.id = oldnewid[node.id]
        if typeof(node) <: Constraint
            node.parent_id = get(oldnewid, node.parent_id, 0)
            node.child_id = get(oldnewid, node.child_id, 0)
        end
    end
end


text:
"""
    @no_prealloc(expr)

Anything inside a `@no_prealloc` will be ignored by AutoPreallocation.jl.
This applies recursively to functions called by functions called etc within the
`@no_prealloc`.
If it allocates then that allocation will always occur -- it won't be avoided by
auto-preallocation.

This can be used if there is a part of a function that may allocation a unknowable
number of intermediary variables to get to the final answer -- and which thus is not
suitable for use with AutoPreallocation.
"""
macro no_prealloc(expr)
    return :(no_prealloc(()->$(esc(expr))))
end

@inline no_prealloc(f) = f()

# Prevent Cassette from recursing into body of f
@inline Cassette.overdub(ctx::RecordingCtx, ::typeof(no_prealloc), f) = f()
@inline Cassette.overdub(ctx::ReplayCtx, ::typeof(no_prealloc), f) = f()


text:
@testset "Generator deprecations" begin
    types_0 = [BullGraph, ChvatalGraph, CubicalGraph, DesarguesGraph,
        DiamondGraph, DodecahedralGraph, FruchtGraph, HeawoodGraph,
        HouseGraph, HouseXGraph, IcosahedralGraph, KarateGraph, KrackhardtKiteGraph,
        MoebiusKantorGraph, OctahedralGraph, PappusGraph, PetersenGraph,
        SedgewickMazeGraph, TetrahedralGraph, TruncatedCubeGraph,
        TruncatedTetrahedronGraph, TruncatedTetrahedronDiGraph, TutteGraph]
    types_1param = [CompleteGraph, CompleteDiGraph,
            StarGraph, StarDigraph, PathGraph,
            PathDiGraph, CycleGraph, CycleDiGraph, WheelGraph, WheelDiGraph,
            BinaryTree, Doublebinary_tree, RoachGraph,
            LadderGraph, Circularladder_graph]
    types_2params = [CompleteBipartiteGraph, LollipopGraph, BarbellGraph,
        TuranGraph, CliqueGraph]
    for G in types_0
        @test_deprecated G()
    end
    for G in types_1param
        @test_deprecated G(5)
    end
    for G in types_2params
        @test_deprecated G(5, 2)
    end
    @test_deprecated CompleteMultipartiteGraph([3, 5])
    @test_deprecated Grid([3, 5])
end


text:
function ispalindrome(n)
	s = string(n)
	return s == reverse(s)
end

function main()
	bestpal = 0
	for i = 999:-1:100
		for j = 999:-1:100
			n = i*j
			if ispalindrome(n)
				if n > bestpal
					bestpal = n
					println(n)
				end
			end
		end
	end
end

main()

max_stars_repo_name: cmey/euler

max_issues_repo_name: cmey/euler

max_forks_repo_name: cmey/euler


text:
function lgd20()
#
#  test polynomial suggested by Goedecker
#   Legendre polynomial of degree 20
#
    p = lgd(20);
    z = [-.9931285991850949; -.9639719272779138; -.9122344282513259; -.8391169718222188;
         -.7463319064601508; -.6360536807265150; -.5108670019508271; -.3737060887154196;
         -.2277858511416451; -.7652652113349733e-1; .7652652113349733e-1; .2277858511416451;
         .3737060887154196; .5108670019508271; .6360536807265150; .7463319064601508;
         .8391169718222188; .9122344282513259; .9639719272779138; .9931285991850949]
    z = [z ones(20)]

    p, PolyZeros(z)
end


text:
"""
`MiscFunctions`

Miscellaneous functions.

**Functions**
- `binary_boundary_search`: Use binary search to find a boundary value.
- `cmean`: Compute circular mean of an array.
- `cstd`: Compute circular standard deviation of an array.
- `find_all_combinations`: Find all combinations of values.
- `find_closest`: Find the position of the closest value in an array.
- `r2_score`: Calculate the coefficient of determination.
- `ptp`: Peak-to-peak range (maximum - minimum).
- `smooth`: Smooth a vector using a moving average filter.
"""
module MiscFunctions

using Statistics

export binary_boundary_search
export cmean, cstd
export find_all_combinations
export find_closest
export r2_score, ptp
export smooth

include("binary_search.jl")
include("circular_statistics.jl")
include("find_all_combinations.jl")
include("find_closest.jl")
include("metrics.jl")
include("smooth.jl")

end  # module

max_stars_repo_name: vkumpost/stoosc

max_issues_repo_name: vkumpost/stoosc

max_forks_repo_name: vkumpost/stoosc


text:
hook(ctx, f, args...) = return

isprimitive(ctx, f, args...) = false
function primitive end

macro hook(ctx, f)
  def = splitdef(f)
  f, args, body, ps = def[:name], def[:args], def[:body], def[:whereparams]
  ctx isa Symbol && (ctx = :(::$ctx))
  isexpr(f, :(::)) || (f = :(::typeof($f)))
  :(function hook($(esc(ctx)), $(esc(f)), $(esc.(args)...)) where {$(esc.(ps)...)}
     $(esc(body))
    end)
end

macro primitive(ctx, f)
  def = splitdef(f)
  f, args, kwargs, body, ps = def[:name], def[:args], def[:kwargs], def[:body], def[:whereparams]
  ctx isa Symbol && (ctx = :(::$ctx))
  isexpr(f, :(::)) || (f = :(::typeof($f)))
  quote
    Vinyl.isprimitive($(esc(ctx)), $(esc(f)), $(esc.(args)...); $(esc.(kwargs)...)) where {$(esc.(ps)...)} = true
    function Vinyl.primitive($(esc(ctx)), $(esc(f)), $(esc.(args)...); $(esc.(kwargs)...)) where {$(esc.(ps)...)}
      $(esc(body))
    end
  end
end


text:
##############################################################
# Using the code from  OnlineStats.jl/src/stats/hist.jl
# Modifying it to work with WeightedOnlineStats
##############################################################

import LinearAlgebra
abstract type WeightedHistogramStat{T} <: WeightedOnlineStat{T} end
abstract type WeightedHist{T} <: WeightedHistogramStat{T} end
split_candidates(o::WeightedHistogramStat) = midpoints(o)
Statistics.mean(o::WeightedHistogramStat) = mean(midpoints(o), fweights(counts(o)))
Statistics.var(o::WeightedHistogramStat) = var(midpoints(o), fweights(counts(o)); corrected=true)
Statistics.std(o::WeightedHistogramStat) = sqrt.(var(o))
Statistics.median(o::WeightedHistogramStat) = quantile(o, .5)

function Base.show(io::IO, o::WeightedHistogramStat)
    print(io, name(o, false, false), ": ")
    print(io, "∑wᵢ=", nobs(o))
    print(io, " | value=")
    show(IOContext(io, :compact => true), value(o))
end

#-----------------------------------------------------------------------# WeightedHist
"""
    WeightedHist(edges; left = true, closed = true)

Create a histogram with bin partition defined by `edges`.
- If `left`, the bins will be left-closed.
- If `closed`, the bin on the end will be closed.
    - E.g. for a two bin histogram ``[a, b), [b, c)`` vs. ``[a, b), [b, c]``

If `edges` is a tuple instead of an array, a multidimensional histogram will be
generated that behaves like a `WeightedOnlineStat{VectorOb}`.

# Examples

    o = fit!(WeightedHist(-5:.1:5), randn(10^6))

    # approximate statistics
    using Statistics
    mean(o)
    var(o)
    std(o)
    quantile(o)
    median(o)
    extrema(o)
    area(o)
    pdf(o)

## 2d Histogram

    hist2d = fit!(WeightedHist((-5:1:5, -5:1:5) ), randn(10000,2), rand(10000))
    value(hist2d).y
"""
struct WeightedHist1D{R} <: WeightedHist{Float64}
    edges::R
    counts::Vector{Int}
    meanw::Vector{Float64}
    outcount::Vector{Int}
    meanwout::Vector{Float64}
    left::Bool
    closed::Bool
end
struct WeightedHistND{R, N} <: WeightedHist{OnlineStats.VectorOb}
    edges::R
    counts::Array{Int,N}
    meanw::Array{Float64,N}
    outcount::Array{Int,N}
    meanwout::Array{Float64,N}
    left::Bool
    closed::Bool
end

function WeightedHist(edges; left::Bool=true, closed::Bool = true)
    edges = isa(edges,Tuple) ? edges : (edges,)
    counts = zeros(Int, map(i->length(i)-1, edges))
    meanw = zeros(Float64, map(i->length(i)-1, edges))
    outcount = zeros(Int,ntuple(_->3,length(edges)))
    meanwout = zeros(Float64,ntuple(_->3,length(edges)))
    if length(edges) == 1
      WeightedHist1D(edges[1],counts,meanw,outcount,meanwout,left,closed)
    else
      WeightedHistND{typeof(edges),length(edges)}(edges, counts, meanw,outcount,meanwout, left, closed)
    end
end
# Special case for 1D Histogram
nobs(o::WeightedHist) = sum(o.counts) + sum(o.outcount)
weightsum(o::WeightedHist) = LinearAlgebra.dot(o.counts, o.meanw) + LinearAlgebra.dot(o.outcount,o.meanwout)
value(o::WeightedHist) = (x=edges(o), y=o.counts .* o.meanw)
binindices(o::WeightedHistND{<:Any,N}, x::AbstractVector) where N = binindices(o, ntuple(i->x[i],N))
binindices(o::WeightedHist1D,x) = OnlineStats.binindex(o.edges, x, o.left, o.closed)
binindices(o::WeightedHistND,x) = CartesianIndex(map((e,ix)->OnlineStats.binindex(e, ix, o.left, o.closed), o.edges, x))
midpoints(o::WeightedHistND) = Iterators.product(map(midpoints,o.edges)...)
midpoints(o::WeightedHist1D) = midpoints(edges(o))
counts(o::WeightedHist) = o.counts
edges(o::WeightedHist) = o.edges
function Statistics.mean(o::WeightedHist)
  weights = value(o).y
  N = ndims(o.counts)
  r = ntuple(N) do idim
    a = map(i->i[idim],midpoints(o))
    mean(a,fweights(weights))
  end
  N==1 ? r[1] : r
end
function Statistics.var(o::WeightedHist)
  weights = value(o).y
  N = ndims(o.counts)
  r = ntuple(N) do idim
    a = map(i->i[idim],midpoints(o))
    var(a,fweights(weights),corrected=true)
  end
  N==1 ? r[1] : r
end
Statistics.std(o::WeightedHist) = sqrt.(var(o))
Statistics.median(o::WeightedHist) = quantile(o, .5)

function Base.extrema(o::WeightedHist1D)
    x, y = midpoints(o), counts(o)
    x[findfirst(!iszero,y)],x[findlast(!iszero,y)]
end
function Base.extrema(o::WeightedHistND{<:Any,N}) where N
    x, y = midpoints(o), counts(o)
    ntuple(N) do idim
      avalue = any(!iszero, y, dims = setdiff(1:N,idim))[:]
      x.iterators[idim][findfirst(avalue)],x.iterators[idim][findlast(avalue)]
    end
end

function Statistics.quantile(o::WeightedHist, p = [0, .25, .5, .75, 1])
    x, y = midpoints(o), counts(o)
    N = ndims(y)
    inds = findall(!iszero, y)
    yweights = fweights(y[inds])
    subset = collect(x)[inds]
    r = ntuple(N) do idim
      data = map(i->i[idim],subset)
      quantile(data, fweights(y[inds]), p)
    end
    if N==1
      return r[1]
    else
      return r
    end
end

function area(o::WeightedHist)
    c = o.counts
    e = o.edges
    return mapreduce(+, CartesianIndices(c)) do I
      ar = prod(map((ed,i)->ed[i+1]-ed[i],e,I.I))
      c[I]*ar
    end
end

outindex(o, ci::CartesianIndex) = CartesianIndex(map((i,l)->i < 1 ? 1 : i > l ? 3 : 2, ci.I, size(o.counts)))
outindex(o, ci::Int) = CartesianIndex(ci < 1 ? 1 : ci > length(o.counts) ? 3 : 2)
function pdf(o::WeightedHist, y)
    ci = binindices(o, y)
    if all(isequal(2),outindex(o,ci).I)
        return o.counts[ci]*o.meanw[ci] / area(o) / weightsum(o)
    else
        return 0.0
    end
end

function _fit!(o::WeightedHist, x, wt)
    #length(x) == N || error("You must provide $(N) values for the histogram")
    ci = binindices(o, x)
    oi = outindex(o,ci)
    if all(isequal(2),oi.I)
        o.counts[ci] += 1
        o.meanw[ci] = smooth(o.meanw[ci], wt, 1.0 / o.counts[ci])
    else
        o.outcount[oi] += 1
        o.meanwout[oi] = smooth(o.meanwout[oi], wt, 1.0 / o.outcount[oi])
    end
end

function _merge!(o::WeightedHist, o2::WeightedHist)
    if o.edges == o2.edges
        for j in eachindex(o.counts)
            newcount = o.counts[j] + o2.counts[j]
            if newcount > 0
              o.meanw[j] = (o.meanw[j]*o.counts[j] + o2.meanw[j]*o2.counts[j])/newcount
            end
            o.counts[j] = newcount
        end
        for j in eachindex(o.outcount)
          newcount = o.outcount[j] + o2.outcount[j]
          if newcount > 0
            o.meanwout[j] = (o.meanwout[j]*o.outcount[j] + o2.meanwout[j]*o2.outcount[j])/newcount
          end
          o.outcount[j] = newcount
        end
    else
        @warn("WeightedHistogram edges do not align.  Merging is approximate.")
        for (yi, wi) in zip(midpoints(o2.edges), o2.counts)
            for k in 1:wi
                _fit!(o, yi)
            end
        end
    end
end

#-----------------------------------------------------------------------# Adaptive Hist
abstract type WeightedHistAlgorithm{N} <: Algorithm end
Base.show(io::IO, o::WeightedHistAlgorithm) = print(io, name(o, false, false))
make_alg(o::WeightedHistAlgorithm) = o

"""
    Weighted Histogram

Calculate a histogram of weighted data.

# Example
    # A weighted histogram with 4 bins:
    o = fit!(WeightedAdaptiveHist(4), rand(1000), rand(1000))

    mean(o)
    var(o)
    std(o)
    median(o)
    quantile(o, [0, 0.25, 0.5, 0.25, 1.0])
    extrema(o)
"""
struct WeightedAdaptiveHist{N, H <: WeightedHistAlgorithm{N}} <: WeightedHistogramStat{N}
    alg::H
    WeightedAdaptiveHist{H}(alg::H) where {N, H<:WeightedHistAlgorithm{N}} = new{N, H}(alg)
end

WeightedAdaptiveHist(args...; kw...) = (alg = make_alg(args...; kw...); WeightedAdaptiveHist{typeof(alg)}(alg))

for f in [:nobs, :counts, :midpoints, :edges, :area]
    @eval $f(o::WeightedAdaptiveHist) = $f(o.alg)
end
for f in [:(_fit!), :pdf, :cdf, :(Base.getindex)]
    @eval $f(o::WeightedAdaptiveHist, y, w) = $f(o.alg, y, w)
end

Base.copy(o::WeightedAdaptiveHist) = WeightedAdaptiveHist(copy(o.alg))

# Base.show(io::IO, o::Hist) = print(io, "Hist: ", o.alg)
OnlineStatsBase._merge!(o::WeightedAdaptiveHist, o2::WeightedAdaptiveHist) = _merge!(o.alg, o2.alg)
function OnlineStatsBase.value(o::WeightedAdaptiveHist)
    (midpoints(o), counts(o))
end

function Base.extrema(o::WeightedAdaptiveHist)
    mids, counts = value(o)
    inds = findall(x->x!=0, counts)  # filter out zero weights
    mids[inds[1]], mids[inds[end]]
end
function Statistics.quantile(o::WeightedAdaptiveHist, p = [0, .25, .5, .75, 1])
    mids, counts = value(o)
    inds = findall(x->x!=0, counts)  # filter out zero weights
    quantile(mids[inds], fweights(counts[inds]), p)
end

function weightsum(o::WeightedAdaptiveHist)
    nobs(o)
end

#-----------------------------------------------------------------------# WeightedAdaptiveBins
struct WeightedAdaptiveBins{T} <: WeightedHistAlgorithm{T}
    value::Vector{Pair{T, T}}
    b::Int
    ex::Extrema{T}
    function WeightedAdaptiveBins{T}(value = Pair{T, T}[], b = 10, ex = Extrema(T)) where T
        new{T}(value, b, ex)
    end
end

Base.copy(o::T) where T <: WeightedAdaptiveBins = T(copy(o.value), copy(o.b), copy(o.ex))

make_alg(T::Type, b::Int) = WeightedAdaptiveBins{T}(Pair{T, T}[], b, Extrema(T))
make_alg(b::Int) = WeightedAdaptiveBins{Float64}(Pair{Float64, Float64}[], b, Extrema(Float64))
midpoints(o::WeightedAdaptiveBins) = first.(o.value)
counts(o::WeightedAdaptiveBins) = last.(o.value)
OnlineStatsBase.nobs(o::WeightedAdaptiveBins) =
    isempty(o.value) ? 0 : sum(last, o.value)
function Base.:(==)(a::T, b::T) where {T<:WeightedAdaptiveBins}
    (a.value == b.value) && (a.b == b.b) && (a.ex == b.ex)
end
Base.extrema(o::WeightedAdaptiveHist{<:Any, <:WeightedAdaptiveBins}) = extrema(o.alg.ex)

# Doesn't happen with weighted stats
OnlineStatsBase._fit!(o::WeightedAdaptiveBins, y::Number, w::Number) =
    _fit!(o, Pair(y, w))

function OnlineStatsBase._fit!(o::WeightedAdaptiveBins{T}, y::Pair) where T
    y2 = convert(Pair{T, T}, y)

    fit!(o.ex, first(y2))
    v = o.value
    i = searchsortedfirst(v, y2)
    insert!(v, i, y2)
    if length(v) > o.b
        # find minimum difference
        i = 0
        mindiff = T(Inf)
        for k in 1:(length(v) - 1)
            @inbounds diff = first(v[k + 1]) - first(v[k])
            if diff < mindiff
                mindiff = diff
                i = k
            end
        end
        # merge bins i, i+1
        q2, k2 = v[i + 1]
        if k2 > 0
            q1, k1 = v[i]
            k3 = k1 + k2
            v[i] = Pair(smooth(q1, q2, k2 / k3), k3)
        end
        deleteat!(o.value, i + 1)
    end
end

function OnlineStatsBase._merge!(o::T, o2::T) where {T <: WeightedAdaptiveBins}
    for v in o2.value
        _fit!(o, v)
    end
    fit!(o.ex, extrema(o2.ex))
end

function Base.getindex(o::WeightedAdaptiveBins, i)
    if i == 0
        return Pair(minimum(o.ex), 0)
    elseif i == (length(o.value) + 1)
        return Pair(maximum(o.ex), 0)
    else
        return o.value[i]
    end
end

# based on linear interpolation
function pdf(o::WeightedAdaptiveBins, x::Number)
    v = o.value
    if x ≤ minimum(o.ex)
        return 0.0
    elseif x ≥ maximum(o.ex)
        return 0.0
    else
        i = searchsortedfirst(v, Pair(x, 0.0))
        x1, y1 = o[i - 1]
        x2, y2 = o[i]
        return smooth(y1, y2, (x - x1) / (x2 - x1)) / area(o)
    end
end

function cdf(o::WeightedAdaptiveBins, x::Number)
    if x ≤ minimum(o.ex)
        return 0.0
    elseif x ≥ maximum(o.ex)
        return 1.0
    else
        i = searchsortedfirst(o.value, Pair(x, 0.0))
        x1, y1 = o[i - 1]
        x2, y2 = o[i]
        w = x - x1
        h = smooth(y1, y2, (x2 - x) / (x2 - x1))
        return (area(o, i-2) + w * h) / area(o)
    end
end

function area(o::WeightedAdaptiveBins, ind = length(o.value))
    out = 0.0
    for i in 1:ind
        w = first(o[i+1]) - first(o[i])
        h = (last(o[i+1]) + last(o[i])) / 2
        out += h * w
    end
    out
end


text:
# https://adventofcode.com/2021/day/17
using AdventOfCode

input = ("(".*replace(readlines("2021/data/day_17.txt")[begin], ".."=>":")[begin+12:end].*")")|>parse
ysolve_t(vy, y) = vy+1/2+sqrt((vy+1/2)^2-2y) # when vertical position reaches y
xsolve_t(vx, x) = (vx^2<2x) ? 10000 : (vx+1/2)-sqrt((vx+1/2)^2-2x) # when horizontal position reaches x
ytimerange(vy, input) = (ceil(Int,ysolve_t(vy,last(input.y))):floor(Int,ysolve_t(vy,first(input.y)))) # return integer time range where probe lies within target
xtimerange(vx,input)=(ceil(Int,xsolve_t(vx,first(input.x))):floor(Int,xsolve_t(vx,last(input.x))))

function part_1(input)
    for vy in -minimum(input.y):-1:minimum(input.y)
        if any(any(in(xtimerange(vx,input)),ytimerange(vy,input)) for vx in 1:maximum(input.x))
            return(Int(vy*(vy+1)/2))
        end
    end
end

function part_2(input)
    sum(any(in(xtimerange(vx,input)),ytimerange(vy,input)) for vx in 1:maximum(input.x), vy in -minimum(input.y):-1:minimum(input.y))
end


text:
function codegen_forward_autodiff_function(::Type{Val{:derivative}}, f::Function)
  @gensym forward_autodiff_function
  quote
    function $forward_autodiff_function(_x::Real)
      getfield(ForwardDiff, :derivative)($f, _x)
    end
  end
end

function codegen_forward_autodiff_function(::Type{Val{:gradient}}, f::Function, chunksize::Integer=0)
  body =
    if chunksize == 0
      :(getfield(ForwardDiff, :gradient)($f, _x))
    else
      :(getfield(ForwardDiff, :gradient)($f, _x, Chunk{$chunksize}()))
    end

  @gensym forward_autodiff_function
  quote
    function $forward_autodiff_function(_x::Vector)
      $body
    end
  end
end

codegen_forward_autodiff_target(method::Symbol, f::Function, chunksize::Integer=0) =
  codegen_forward_autodiff_target(Val{method}, f, chunksize)

function codegen_forward_autodiff_target(::Type{Val{:hessian}}, f::Function, chunksize::Integer=0)
  body =
    if chunksize == 0
      :(-getfield(ForwardDiff, :hessian)($f, _x))
    else
      :(-getfield(ForwardDiff, :hessian)($f, _x, Chunk{$chunksize}()))
    end

  @gensym forward_autodiff_target
  quote
    function $forward_autodiff_target(_x::Vector)
      $body
    end
  end
end

function codegen_forward_autodiff_uptofunction(::Type{Val{:derivative}}, f::Function)
  @gensym forward_autodiff_uptofunction
  quote
    function $forward_autodiff_uptofunction(_x::Real)
      result = ForwardDiff.DerivativeResult(_x)
      getfield(ForwardDiff, :derivative)(result, $f, _x)
      return ForwardDiff.value(result), ForwardDiff.derivative(result)
    end
  end
end

function codegen_forward_autodiff_uptofunction(::Type{Val{:gradient}}, f::Function, chunksize::Integer=0)
  adfcall =
    if chunksize == 0
      :(getfield(ForwardDiff, :gradient!)(result, $f, _x))
    else
      :(getfield(ForwardDiff, :gradient!)(result, $f, _x, Chunk{$chunksize}()))
    end

  @gensym forward_autodiff_uptofunction
  quote
    function $forward_autodiff_uptofunction(_x::Vector)
      result = ForwardDiff.GradientResult(_x)
      $adfcall
      return ForwardDiff.value(result), ForwardDiff.gradient(result)
    end
  end
end

codegen_forward_autodiff_uptotarget(method::Symbol, f::Function, chunksize::Integer=0) =
  codegen_forward_autodiff_uptotarget(Val{method}, f, chunksize)

function codegen_forward_autodiff_uptotarget(::Type{Val{:hessian}}, f::Function, chunksize::Integer=0)
  adfcall =
    if chunksize == 0
      :(getfield(ForwardDiff, :hessian!)(result, $f, _x))
    else
      :(getfield(ForwardDiff, :hessian!)(result, $f, _x, Chunk{$chunksize}()))
    end

  @gensym forward_autodiff_uptotarget
  quote
    function $forward_autodiff_uptotarget(_x::Vector)
      result = ForwardDiff.HessianResult(_x)
      $adfcall
      return ForwardDiff.value(result), ForwardDiff.gradient(result), -ForwardDiff.hessian(result)
    end
  end
end


text:
#######
# useful concepts for dispatch
#######

abstract type AbstractSimulation end

struct Simulation{𝒜, ℬ, 𝒞, 𝒟, ℰ, ℱ, O, DG, N} <: AbstractSimulation
    model::𝒜
    state::ℬ
    timestepper::𝒞
    initial_conditions::𝒟
    callbacks::ℰ
    simulation_time::ℱ
    odesolver::O
    dgmodel::DG
    name::N
end

function Simulation(;
    model = nothing,
    state = nothing,
    timestepper = nothing,
    initial_conditions = nothing,
    callbacks = nothing,
    simulation_time = nothing,
    odesolver = nothing,
    dgmodel = nothing,
    name = nothing,
)
    # initialize DGModel (rhs)
    dgmodel = DGModel(model) # DGModel --> KernelModel, to be more general?

    FT = eltype(dgmodel.grid.vgeo)

    # initialize state variables
    if state == nothing
        state = init_ode_state(dgmodel, FT(0); init_on_cpu = true)
    end

    # initialize timestepper
    odesolver = timestepper.method(dgmodel, state; dt = timestepper.timestep, t0 = simulation_time[1])

    return Simulation(
        model,
        state,
        timestepper,
        initial_conditions,
        callbacks,
        simulation_time,
        odesolver,
        dgmodel,
        name,
    )
end

"""
calculate_dt(grid, wavespeed = nothing, diffusivity = nothing, viscocity = nothing, cfl = 0.1)
"""
function calculate_dt(grid; wavespeed = nothing, diffusivity = nothing, viscocity = nothing, cfl = 0.1)
    Δx = min_node_distance(grid, HorizontalDirection())
    Δts = []
    if wavespeed != nothing
        push!(Δts, Δx / wavespeed)
    end
    if diffusivity != nothing
        push!(Δts, Δx^2 / diffusivity)
    end
    if viscocity != nothing
        push!(Δts, Δx^2 / viscocity)
    end
    if Δts == []
        @error("Please provide characteristic speed or diffusivities")
        return nothing
    end
    return cfl * minimum(Δts)
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
module ExampleTests

using Test, Cassette
using Cassette: @context, @prehook, @posthook, @primitive, @pass, overdub, Box, Tag

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

function rosenbrock(x::Vector{Float64})
    a = 1.0
    b = 100.0
    result = 0.0
    for i in 1:length(x)-1
        result += (a - x[i])^2 + b*(x[i+1] - x[i]^2)^2
    end
    return result
end

x = rand(2)

@context RosCtx

MESSAGES = String[]
@prehook (f::Any)(args...) where {__CONTEXT__<:RosCtx} = push!(MESSAGES, string("calling ", f, args))
@test overdub(RosCtx, rosenbrock)(x) == rosenbrock(x)
@test length(MESSAGES) > 90

@prehook (f::Any)(args...) where {__CONTEXT__<:RosCtx} = push!(__trace__.metadata, string("calling ", f, args))
meta = String[]
@test overdub(RosCtx, rosenbrock, metadata = meta)(x) == rosenbrock(x)
@test MESSAGES == meta

@prehook (f::Any)(args...) where {__CONTEXT__<:RosCtx} = nothing
@prehook (f::Any)(args::Number...) where {__CONTEXT__<:RosCtx} = push!(__trace__.metadata, args)
meta = Any[]
@test overdub(RosCtx, rosenbrock, metadata = meta)(x) == rosenbrock(x)
for args in meta
    @test all(x -> isa(x, Number), args)
end

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

x = rand()
sin_plus_cos(x) = sin(x) + cos(x)
@context SinCtx
@test overdub(SinCtx, sin_plus_cos)(x) === sin_plus_cos(x)
@primitive sin(x) where {__CONTEXT__<:SinCtx} = cos(x)
@test overdub(SinCtx, sin_plus_cos)(x) === (2 * cos(x))

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

x = 2
foldmul(x, args...) = Core._apply(Base.afoldl, (*, x), args...)
@context FoldCtx
@test overdub(FoldCtx, foldmul)(x) === foldmul(x)

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

@context CountCtx
count1 = Ref(0)
@prehook (f::Any)(args::Number...) where {__CONTEXT__<:CountCtx} = (__trace__.metadata[] += 1)
overdub(CountCtx, sin, metadata = count1)(1)
@prehook (f::Any)(args::Number...) where {__CONTEXT__<:CountCtx} = (__trace__.metadata[] += 2)
count2 = Ref(0)
overdub(CountCtx, sin, metadata = count2)(1)
@test (2 * count1[]) === count2[]

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

square_closure(x) = (y -> y * x)(x)
@context SqrCtx
x = rand()
@test square_closure(x) == overdub(SqrCtx, square_closure)(x)

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

comprehension1(x) = [i for i in x]
comprehension2(f, x, y) = [f(x, i) for i in y]
@context CompCtx
f, x, y = hypot, rand(), rand(2)
@test comprehension1(x) == overdub(CompCtx, comprehension1)(x)
@test comprehension1(y) == overdub(CompCtx, comprehension1)(y)
@test comprehension2(f, x, y) == overdub(CompCtx, comprehension2)(f, x, y)

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

struct Baz
    x::Int
    y::Float64
    z::String
end

baz_identity(x::Int) = Baz(x, float(x), "$x").x

@context BazCtx
Cassette.metatype(::Type{<:BazCtx}, ::Type{<:Integer}) = Float64
n = rand()
ctx = BazCtx(baz_identity)
result = overdub(ctx, baz_identity; boxes = Val(true))(Box(ctx, 1, n))
@test result === Box(ctx, 1, n)

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

struct Bar{X,Y,Z}
    x::X
    y::Y
    z::Z
end

mutable struct Foo
    a::Bar{Int}
    b
end

function foo_bar_identity(x)
    bar = Bar(x, x + 1, x + 2)
    foo = Foo(bar, "ha")
    foo.b = bar
    foo.a = Bar(4,5,6)
    foo2 = Foo(foo.a, foo.b)
    foo2.a = foo2.b
    return foo2.a.x
end

@context FooBarCtx
Cassette.metatype(::Type{<:FooBarCtx}, ::Type{<:Integer}) = Float64
n = rand()
ctx = FooBarCtx(foo_bar_identity)
result = overdub(ctx, foo_bar_identity; boxes = Val(true))(Box(ctx, 1, n))
@test result === Box(ctx, 1, n)

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

sig_collection = DataType[]
@context PassCtx
overdub(PassCtx, sum; pass = @pass((sig, cinfo) -> (push!(sig_collection, sig); cinfo)))(rand(3))
@test !isempty(sig_collection) && all(T -> T <: Tuple, sig_collection)

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

mutable struct Count{T}
    count::Int
end

@context CountCtx2

@prehook function (::Any)(arg::T, args::T...) where {T,__CONTEXT__<:CountCtx2,__METADATA__<:Count{T}}
    __trace__.metadata.count += 1
end

mapstr(x) = map(string, x)
c = Count{Union{String,Int}}(0)
@test overdub(CountCtx2, mapstr, metadata = c)(1:10) == mapstr(1:10)
@test c.count > 1000

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

@context NestedCtx

function nested_test(n, x)
    if n < 1
        return sin(x)
    else
        return overdub(NestedCtx, nested_test)(n - 1, x + x)
    end
end

x = rand()
ctxs = NestedCtx[]
tag_id = objectid(typeof(nested_test))
nested_ctx_type = NestedCtx{Tag{NestedCtx{Tag{NestedCtx{Tag{Nothing,tag_id}},tag_id}},tag_id}}

@prehook function (::Any)(args...) where {__CONTEXT__<:NestedCtx}
    !(in(__trace__.context, ctxs)) && push!(ctxs, __trace__.context)
end

@test overdub(NestedCtx, nested_test)(2, x) === sin(x + x + x + x)
@test any(x -> isa(x, nested_ctx_type), ctxs)

end # module


text:
import Interpolations.AbstractExtrapolation

abstract type AbstractInterpolationScheme end

struct InterpolationScheme1D <: AbstractInterpolationScheme
    xgrid::InterpolationGrid
    interpolation_method::AbstractInterpolation
end

export AbstractInterpolationScheme,
    InterpolationScheme1D,
    InterpolationGrid