text:
print("Loading packages... ")
using Images, FileIO, StatsBase, ImageCore
println("Done!")
###############################################################################
# DATA DEFINITIONS                                                            #
# {red|green|blue}map = current value of each pixel                           #
# {red|green|blue}counts = pixels in each bin                                 #
# {red|green|blue}vars: red[1000] = what 1000 mapped to with this transform   #
###############################################################################
########################## defaults ##############################################

obase = "a-stretched-image"  # default output base file name

rootpower    = 6.0      # default power factor: 1/rootpower for all iterations
rootiter     = 1        # number of iterations on applying rootpower - sky
rootpower2   = 1.0      # if rootiter =2 use this power on iteration 2 if the user sets it >1
pcntclip     = 0.005    # default percent clip level = total pixels * pcntclip/100

colorcorrect = 1        # default is to docolor correction.  (turn off to see the difference)

colorenhance = 1.0      # default enhancement value.
tonecurve    = 0        # no application odf a tone curve
specprhist   = 0        # no output histogram to specpr file
cumstretch   = 0        # do not do cumulative histogram stretch after rootpower
skylevelfactor = 0.06   # sky level relative to the histogram peak  (was 0.03 v0.88 and before)
rgbskylevel   = 1024.0    # desired  on a 16-bit  0 to 65535 scale
rgbskylevelrs = 4096.0    # desired  on output root stretched image 16-bit  0 to 65535 scale
zeroskyred    = 4096.0    # desired zero point on sky, red   channel
zeroskygreen  = 4096.0    # desired zero point on sky, green channel
zeroskyblue   = 4096.0    # desired zero point on sky, blue  channel
scurve        = 1         # no s-curve application
setmin        = 1         # no modification of the minimum
setminr = 0.0  # minimum for red
setming = 0.0  # minimum for green
setminb = 0.0  # minimum for blue
idisplay = 0   # do not display the final image
jpegonly = 0   # do not do jpeg only (will do jpeg + 16-bit png)
saveinputminussky  = 0  #  save input image - sky

debuga  = 0  # set to 1 for debugging, or 0 for none
doplots = 0  # show plots of histograms

if doplots != 0
  print("Loading Plots package... ")
  using Plots
  pyplot()
  println("Done!")
else # give plot() a definition
  function plot(args...)
  end
end

redvals = 0
greenvals = 0
bluevals = 0 # to put these in the global namespace

function plotter(title)
  #=edges = (1:256) .* 256
  r = map(x -> log(x + 1), fit(Histogram, redcounts, edges).weights)
  g = map(x -> log(x + 1), fit(Histogram, greencounts, edges).weights)
  b = map(x -> log(x + 1), fit(Histogram, bluecounts, edges).weights)
  print(size(r))
  for i = 1:255
    println(i, ": ", r[i], ", ", g[i], " ", b[i])
  end
#  a = cat(2, r, g, b)
#  plot(a,legend = nothing, color=reshape([:red, :green, :blue], 1, 3), title = title)
#  gui()=#
end

function checksum(v)
  result = reduce($, 0, map(Int,map(floor, v)))
  return result
end

function remap()
  # This function needs to be optimized!
  # if redmap[1000] used to be 987, it becomes where redvals[987] is going
  println("Remapping.... ")
  global redcounts, redmap, greencounts, greenmap, bluecounts, bluemap
  global redvals, greenvals, bluevals
  redvals = map(x -> floor(UInt64, x), redvals)
  greenvals = map(x -> floor(UInt64, x), greenvals)
  bluevals = map(x -> floor(UInt64, x), bluevals)
  redmap = map(x -> redvals[x] + 1, redmap)
  greenmap = map(x -> greenvals[x] + 1, greenmap)
  bluemap = map(x -> bluevals[x] + 1, bluemap)
  println("recounting... ")
  newred = zeros(65536)
  newgreen = zeros(65536)
  newblue = zeros(65536)

  for i = 1:65536
    newred[redvals[i] + 1] += redcounts[i]
    newgreen[greenvals[i] + 1] += greencounts[i]
    newblue[bluevals[i] + 1] += bluecounts[i]
  end

  redcounts = newred
  greencounts = newgreen
  bluecounts = newblue
  println("Done!")
end

function showmeans()
  v = collect(0:65535)

  redmin = reduce(min, npix, v[redcounts .>0])
  greenmin = reduce(min, npix, v[greencounts .>0])
  bluemin = reduce(min, npix, v[bluecounts .>0])

  redmax = reduce(max, 0, v[redcounts .>0])
  greenmax = reduce(max, 0, v[greencounts .>0])
  bluemax = reduce(max, 0, v[bluecounts .>0])

  redmean = sum((0:65535).* redcounts) / npix
  greenmean = sum((0:65535) .* greencounts) / npix
  bluemean = sum((0:65535) .* bluecounts) / npix

  @printf "  red: min = %5d max = %5d mean = %7.1f\n" redmin redmax redmean
  @printf "green: min = %5d max = %5d mean = %7.1f\n" greenmin greenmax greenmean
  @printf " blue: min = %5d max = %5d mean = %7.1f\n" bluemin bluemax bluemean
end

function smoothsky()
  global redvals, greenvals, bluevals
  # histogram plotting disabled
  # now find the sky level on the left side of the histogram

  println("Calculating sky smoothing peaks")
  redacc = cumsum(redcounts)
  greenacc = cumsum(greencounts)
  blueacc = cumsum(bluecounts)
  smoothred = zeros(65536)
  smoothgreen = zeros(65536)
  smoothblue = zeros(65536)
  width = 300          # include this many bins on either side
  redpeak = 0
  maxred = 0
  greenpeak = 0
  maxgreen = 0
  bluepeak = 0
  maxblue = 0
  for i = 400:65500 # ends are left off to avoid clipping artifacts
    low = i - width - 1 < 1 ? 1 : i - width - 1
    high = i + width > 65536 ? 65536 : i + width
    smoothred[i] = (redacc[high] - redacc[low]) / (high - low)
    smoothgreen[i] = (greenacc[high] - greenacc[low]) / (high - low)
    smoothblue[i] = (blueacc[high] - blueacc[low]) / (high - low)
    if smoothred[i] > maxred
      maxred = smoothred[i]
      redpeak = i
    end
    if smoothgreen[i] > maxgreen
      maxgreen = smoothgreen[i]
      greenpeak = i
    end
    if smoothblue[i] > maxblue
      maxblue = smoothblue[i]
      bluepeak = i
    end
  end
  if redpeak == 0 || greenpeak == 0 || bluepeak == 0
    @printf "Max search failed red %d green %d blue %d\n" redpeak greenpeak bluepeak
    exit(-1)
  end
  println("\n     image histogram:    max bin    pixel count")
  @printf "                  red:   %d         %f\n" redpeak maxred
  @printf "                green:   %d         %f\n" greenpeak maxgreen
  @printf "                 blue:   %d         %f\n" bluepeak maxblue

  redsky = maxred * skylevelfactor
  greensky = maxgreen * skylevelfactor
  bluesky = maxblue * skylevelfactor
  redskybin = 0
  greenskybin = 0
  blueskybin = 0

  #find where redvals[i] >= greenske and redlevel[i-1] <= greensky
  for i = redpeak:-1:2
    if smoothred[i] >= greensky && smoothred[i - 1] <= greensky
      redskybin = i + 1
      break
    end
  end

  for i = greenpeak:-1:2
    if smoothgreen[i] >= greensky && smoothgreen[i - 1] <= greensky
      greenskybin = i + 1
      break
    end
  end

  for i = bluepeak:-1:2
    if smoothblue[i] >= greensky && smoothblue[i - 1] <= greensky
      blueskybin = i + 1
      break
    end
  end

  if redskybin == 0 || greenskybin == 0 || blueskybin == 0
    @printf "histogram sky level %f not found:\n" skylevelfactor
    @printf "     channels: red %d   green %d   blue %d\n" redskybin greenskybin blueskybin
    println("          Try increasing -rgbskyzero values\n")
    exit(-1)
  end

  println("Histogram dark sky level")
  println("\n                  bin      Number of pixels in bin")     # need to compute cumulative hstogram
  @printf "           red:   %d         %d\n" redskybin   smoothred[redskybin]
  @printf "         green:   %d         %d\n" greenskybin smoothgreen[greenskybin]
  @printf "          blue:   %d         %d\n" blueskybin  smoothblue[blueskybin]

  redskysub1 = redskybin - zeroskyred
  greenskysub1 = greenskybin - zeroskygreen
  blueskysub1 = blueskybin - zeroskyblue

	@printf "\nsubtract %7d from red   to make red   sky align with zero reference sky: %7.1f\n" redskysub1   zeroskyred
	@printf "subtract %7d from green to make green sky align with zero reference sky: %7.1f\n" greenskysub1 zeroskygreen
	@printf "subtract %7d from blue  to make blue  sky align with zero reference sky: %7.1f\n" blueskysub1  zeroskyblue

	#@printf "\nnow set the RGB sky zero level to %7.1f  %7.1f  %7.1f bin out of 65535\n", zeroskyred zeroskygreen zeroskyblue

  redvals = greenvals = bluevals = collect(0:65535)
  redvals = (redvals .- redskysub1) .* (65535 / (65535 - redskysub1))
  greenvals = (greenvals .- greenskysub1) .* (65535 / (65535 - greenskysub1))
  bluevals = (bluevals .- blueskysub1) .* (65535 / (65535 - blueskysub1))
  redvals[redvals .< 0] = 0
  greenvals[greenvals .< 0] = 0
  bluevals[bluevals .< 0] = 0
  redvals = floor(redvals)
  greenvals = floor(greenvals)
  bluevals = floor(bluevals)
  remap()

  println("sky adjusted")
  showmeans()
end
###### end of function smoothsky

function setmins()
  global redvals, greenvals, bluevals
  println("\nsetting minimum levels.")
  zx = 0.2
  redvals = collect(0:65535)
  screen = redvals .< setminr
  redvals[screen] = redvals[screen] .* zx .+ setminr
  greenvals = collect(0:65535)
  screen = greenvals .< setming
  greenvals[screen] = greenvals[screen] .* zx .+ setming
  bluevals = collect(0:65535)
  screen = bluevals .< setminb
  bluevals[screen] = bluevals[screen] .* zx .+ setminb
  remap()

end

filename = "BaseOrion.tif"
print("Loading $filename... ")
tic()
img = load(filename)
save("save1.png", img)
if img == nothing
  println("Image failed to load.")
  exit(-1)
end
println("Done!")
println("Image loading time: ", toq())
println("Loading processing vectors and histograms...")
tic()
sz = size(img)
px = sz[1]
py = sz[2]
npix = px * py
npixm = floor(npix * pcntclip / 100.0)
npixm = npixm < 1 ? 1 : npixm
@printf "\nImage size: %d x %d pixels, cut level = %d\n" px py npixm
img = reinterpret(UInt16, img)
imgr = reshape(img[1, :, :], npix)
imgg = reshape(img[2, :, :], npix)
imgb = reshape(img[3, :, :], npix)
img = 0
#=
redcounts = fit(Histogram, imgr, 0:65536).weights
greencounts = fit(Histogram, imgg, 0:65536).weights
bluecounts = fit(Histogram, imgb, 0:65536).weights
=#
# still can't make this work right :(
redcounts = zeros(UInt32, 65536)
greencounts = zeros(UInt32, 65536)
bluecounts = zeros(UInt32, 65536)
for i = 1:npix
  redcounts[imgr[i] + 1] += 1
  greencounts[imgg[i] + 1] += 1
  bluecounts[imgb[i] + 1] += 1
end

println("Histogram loading time: ", toq())

plotter("initial")
if doplots == 1
  plotter("Initial")
end

# pixel values are off by one throughout calculations
# because of 1-base arrays

redmap = collect(1:65536)
greenmap = collect(1:65536)
bluemap = collect(1:65536)

print("\n\nInitial bounds and means:\n")
showmeans()

if tonecurve == 1
  print("\n\nApplying tone curve... ")
  b = 12
  c = 65535
  d = 12
  redvals = greenvals = bluevals = collect(0:65535)
  # factor = b * ((1 / d)) ^ ((val / c) ^ 0.4)
  # could probably be done in one statement but since it's only done once,
  # it's not time-critical and not worth the time to figure out
  f = (map( x-> x ^ 0.4, (collect(0:65535) ./ c)))
  f = map( x -> b * (1 / d) ^ x, f)
  g = summarystats(f)
  print(g)
  redvals = floor(redvals .* f)
  greenvals = floor(greenvals .* f)
  bluevals = floor(bluevals .* f)
  println("Done!")
  remap()
  print("after tone curve:\n")
  showmeans()
else
  print("\nNot doing tone curve\n")
end

alowred   = 0  # find the DN level which is the cut level
alowgreen = 0
alowblue  = 0
nred   = redcounts[1]     # number of pixels in histogram bin
ngreen = greencounts[1]
nblue  = bluecounts[1]

for i = 2:32000   # upper limit is more or less pulled out of the air
                  # don't need to look at the first element because it's already zero
  if alowred == 0
    nred = nred + redcounts[i]
    if nred > npixm
      alowred = i
    end
  end
  if alowgreen == 0
    ngreen = ngreen + greencounts[i]
    if ngreen > npixm
      alowgreen = i
    end
  end
  if alowblue == 0
    nblue = nblue + bluecounts[i]
    if nblue > npixm
      alowblue = i
    end
  end
  if alowred > 0 && alowgreen > 0 && alowblue > 0
    break
  end
end

@printf "\n\ncumulative histogram cut level (%f %% = %d pixels):\n" pcntclip floor(npixm)
@printf "           red:   %d  (%d pixels)\n" alowred nred
@printf "           green: %d  (%d pixels)\n" alowgreen ngreen
@printf "           blue:  %d  (%d pixels)\n" alowblue nblue

# reset dark sky
for pass = 1:2
  @printf "\n\ndark sky subtraction pass %d\n" pass
#  println("r g b checksums before ", checksum(redvals), " ", checksum(greenvals), " ", checksum(bluevals))
  smoothsky()
end

afr = map(x -> redmap[x + 1], imgr)
afg = map(x -> greenmap[x + 1], imgg)
afb = map(x -> bluemap[x + 1], imgb)
# assemble sky-subtracted image here
g = summarystats(afr)
println(">>> afr summary")
print(g)
println("============= computing root stretch =============")
for pass = 1:rootiter
  if pass == 2
    x = 1 / rootpower2
  else
    x = 1 / rootpower
  end

  println("Root stretch pass $pass of $rootiter")
  redfloat = collect(1:65536) / 65536 # has a +1 built in
  greenfloat = collect(1:65536) / 65536
  bluefloat = collect(1:65536) / 65536
  redfloat = map(z -> z ^ x, redfloat)
  greenfloat = map(z -> z ^ x, greenfloat)
  bluefloat = map(z -> z ^ x, bluefloat)
  redfloat = redfloat .* 65535
  greenfloat = greenfloat .* 65535
  bluefloat = bluefloat .* 65535
#  remap()
  redmin = reduce(min, 100000, redfloat[redcounts .> 0])
  greenmin = reduce(min, 100000, greenfloat[greencounts .> 0])
  bluemin = reduce(min, 100000, bluefloat[bluecounts .> 0])
  bmin = min(redmin, greenmin, bluemin)
  println("redmin $redmin greenmin $greenmin bluemin $bluemin bmin $bmin")
  if bmin < 4096
    bminz = 0
  else
    bminz = bmin - 4096
  end
  println("Subtracting $bminz from root stretched image")
  redfloat = (redfloat .- bminz) / (65535 - bminz) * 65535
  greenfloat = (greenfloat .- bminz) / (65535 - bminz) * 65535
  bluefloat = (bluefloat .- bminz) / (65535 - bminz) * 65535
  redfloat[redcounts .== 0] = 0
  greenfloat[greencounts .== 0] = 0
  bluefloat[bluecounts .== 0] = 0
  redvals = map(x -> floor(UInt16, x), redfloat)
  greenvals = map(x -> floor(UInt16, x), greenfloat)
  bluevals = map(x -> floor(UInt16, x), bluefloat)
  remap()

  showmeans()

  for pass2 = 1:(rootpower > 60 ? 3 : 2)
    println("Smooth sky pass $pass2 within root stretch pass $pass")
    smoothsky()
  end
end

if scurve > 0
  println("========== Doing s-curve stretching ==========")
  for i = 1:scurve
    if i == 2 || i == 4
      xfactor = 3
      xoffset = 0.22
    else
      xfactor = 5
      xoffset = 0.42
    end

    scurvemin = (xfactor / (1.0 + exp(-1.0 * (0 - xoffset) * xfactor))) - (1.0 - xoffset)
    scurvemax = (xfactor / (1.0 + exp(-1.0 * (1 - xoffset) * xfactor))) - (1.0 - xoffset)
    scurveminsc = scurvemin / scurvemax

    println("s-curve pass $i")
    println("          xfactor     = $xfactor")
    println("          xoffset     = $xoffset")
    @printf "          scurvemin   = %8.4f\n" scurvemin
    @printf "          scurvemax   = %8.4f\n" scurvemax
    @printf "          scurveminsc = %8.4f\n" scurveminsc

    xo = 1 - xoffset

    redsc = (collect(0:65535) ./ 65535 .- xoffset) .* -xfactor
    redsc = map(exp, redsc)
    redsc = redsc .+ 1
    redsc = (xfactor ./ redsc) .- xo
    redsc = redsc ./ scurvemax .- scurveminsc
    redsc = 65535 .* redsc
    redvals = map(x -> floor(UInt64, x), redsc ./ (1 - scurveminsc))

    greensc = (collect(0:65535) ./ 65535 .- xoffset) .* -xfactor
    greensc = map(exp, greensc)
    greensc = greensc .+ 1
    greensc = (xfactor ./ greensc) .- xo
    greensc = greensc ./ scurvemax .- scurveminsc
    greensc = 65535 .* greensc
    greenvals = map(x -> floor(UInt64, x), greensc ./ (1 - scurveminsc))

    bluesc = (collect(0:65535) ./ 65535 .- xoffset) .* -xfactor
    bluesc = map(exp, bluesc)
    bluesc = bluesc .+ 1
    bluesc = (xfactor ./ bluesc) .- xo
    bluesc = bluesc ./ scurvemax .- scurveminsc
    bluesc = 65535 .* bluesc
    bluevals = map(x -> floor(UInt64, x), bluesc ./ (1 - scurveminsc))

    remap()

    println("\n")
    showmeans()
  end

  for j = 1:2
    println("Doing sky offset pass $j")
    smoothsky()
  end
end

if setmin > 0
  setmins()
end

#plotter("Before final load")
redvals = 0
greenvals = 0
bluevals = 0
newred = 0
newgreen = 0
newblue = 0

cr = map(x -> redmap[x + 1] - 1, imgr)
cg = map(x -> greenmap[x + 1] - 1, imgg)
cb = map(x -> bluemap[x + 1] - 1, imgb)
cr[cr .== 0] = 1
cg[cg .== 0] = 1
cb[cb .== 0] = 1

redmap = 0
greenmap = 0
bluemap = 0
imgr = 0
imgg = 0
imgb = 0

crmin = reduce(min, 65537, cr)
crmax = reduce(max, 0, cr)
crmean = mean(cr)

cgmin = reduce(min, 65537, cg)
cgmax = reduce(max, 0, cg)
cgmean = mean(cg)

cbmin = reduce(min, 65537, cb)
cbmax = reduce(max, 0, cb)
cbmean = mean(cb)

println("Image reloaded")
@printf "red   min = %d max = %d mean = %7.1f\n" crmin crmax crmean
@printf "green min = %d max = %d mean = %7.1f\n" cgmin cgmax cgmean
@printf "blue  min = %d max = %d mean = %7.1f\n" cbmin cbmax cbmean

if colorcorrect > 0
  println("\n ========== final color adjustment ==========")
  afr = afr .- zeroskyred
  afr[afr .< 10] = 10
  afg = afg .- zeroskygreen
  afg[afg .< 10] = 10
  afb = afb .- zeroskyblue
  afb[afb .< 10] = 10

  println("color adjustments")

  grratio = (afg ./ afr) ./ (cg ./ cr)
  brratio = (afb ./ afr) ./ (cb ./ cr)
  rgratio = (afr ./ afg) ./ (cr ./ cg)
  bgratio = (afb ./ afg) ./ (cb ./ cg)
  gbratio = (afg ./ afb) ./ (cg ./ cb)
  rbratio = (afr ./ afb) ./ (cr ./ cb)

  zmin = 0.2
  zmax = 1.0

  grratio[grratio .< zmin] = zmin
  grratio[grratio .> zmax] = zmax

  brratio[brratio .< zmin] = zmin
  brratio[brratio .> zmax] = zmax

  rgratio[rgratio .< zmin] = zmin
  rgratio[rgratio .> zmax] = zmax

  bgratio[bgratio .< zmin] = zmin
  bgratio[bgratio .> zmax] = zmax

  gbratio[gbratio .< zmin] = zmin
  gbratio[gbratio .> zmax] = zmax

  rbratio[rbratio .< zmin] = zmin
  rbratio[rbratio .> zmax] = zmax

  println("\nColor ratio images after constraint to interval ($zmin, $zmax)")
  g = summarystats(grratio)
  print(g)
  g = summarystats(brratio)
  print(g)
  g = summarystats(rgratio)
  print(g)
  g = summarystats(bgratio)
  print(g)
  g = summarystats(gbratio)
  print(g)
  g = summarystats(rbratio)
  print(g)

  avgn = (cr .+ cg .+ cb) ./ (3 * 65535)
  avgn[avgn .< 0] = 0
  cmax = reduce(max, 0, avgn)
  if cmax < 1
    avgn ./= cmax
  end

  avgn .^= 0.2
  avgn .+= 0.3
  avgn ./= 1.3
  cmin = reduce(min, 65537, avgn)
  cmax = reduce(max, 0, avgn)
  cmean = mean(avgn)

  @printf "Color correction intensity min = %7.5f max = %7.5f mean = %7.5f\n" cmin cmax cmean

  cfactor = 1.2
  cfe = avgn .* (cfactor * colorenhance)
  grratio = 1 .+ cfe .* (grratio .- 1)
  brratio = 1 .+ cfe .* (brratio .- 1)
  rgratio = 1 .+ cfe .* (rgratio .- 1)
  bgratio = 1 .+ cfe .* (bgratio .- 1)
  gbratio = 1 .+ cfe .* (gbratio .- 1)
  rbratio = 1 .+ cfe .* (rbratio .- 1)

  println("Color ratio images after factors applied")
  println("grratio:")
  describe(grratio)
  println("brratio:")
  describe(brratio)
  println("rgratio:")
  describe(rgratio)
  println("bgratio:")
  describe(bgratio)
  println("gbratio:")
  describe(gbratio)
  println("rbratio:")
  describe(rbratio)

  cmin = reduce(min, 65537, cfe)
  cmax = reduce(max, 0, cfe)
  cmean = mean(cfe)
  @printf "Color correction cfe min = %7.5f max = %7.5f mean = %7.5f\n" cmin cmax cmean

  println("Computing possible image sets")
  c2gr = cg .* grratio
  c3br = cb .* brratio
  c1rg = cr .* rgratio
  c3bg = cb .* bgratio
  c1rb = cr .* rbratio
  c2gb = cg .* gbratio

  println("Starting signal dependent color recovery")
  println("Processing $npix pixels....")
  decile = floor(npix / 10)
  pixelcount = decile
  for i = 1:npix
    pixelcount -= 1
    if pixelcount == 0
      println("Processing pixel $i of $npix...")
      pixelcount = decile
    end
    # I think this could be done with stacked ternary operators
    maxv = cr[i]
    maxch = 1
    if cg[i] > maxv
      maxv = cg[i]
      maxch = 2
    end
    if cb[i] > maxv
      # don't need the max value!
      maxch = 3
    end
    if maxch == 1
      cg[i] = floor(UInt16, c2gr[i])
      cb[i] = floor(UInt16, c3br[i])
    elseif maxch == 2
      cr[i] = floor(UInt16, c1rg[i])
      cb[i] = floor(UInt16, c3bg[i])
    else
      cr[i] = floor(UInt16, c1rb[i])
      cg[i] = floor(UInt16, c2gb[i])
    end
  end

  println("Color recovery complete")
end

if setmin > 0
  setmins()
end

smoothsky()
#=if cumstretch == 1
  println("Performing cumulative histogram strech")
  hist = fit(Histogram, [cr, cg, cb], 0:65536).weights
  #redcounts = fit(Histogram, imgr, 0:65536).weights
  cs = chist
  gap = 15000
  for i=1:65535
    ilow = i <= gap ? 1 : i - gap
    ihigh = i >= 65535 - gap ? 65535 : i + gap
    =#

cr[cr .< 0] = 0
cr[cr .> 65535] = 65535
cg[cg .< 0] = 0
cg[cg .> 65535] = 65535
cb[cb .< 0] = 0
cb[cb .> 65535] = 65535

finalr = reshape(map(x -> floor(UInt16, x), cr), px, py)
finalg = reshape(map(x -> floor(UInt16, x), cg), px, py)
finalb = reshape(map(x -> floor(UInt16, x), cb), px, py)
final = zeros(UInt16, 3, px, py)
final[1, :, :] = finalr
final[2, :, :] = finalg
final[3, :, :] = finalb

final = colorview(RGB{N0f16}, final)

#final = reinterpret()
oname = "$obase.png"
oname = "neworionout.png"
print("Saving $oname... ")
save(oname, final)
print("Done!")

max_stars_repo_name: FatherGeoffHorton/AstroTools

max_issues_repo_name: FatherGeoffHorton/AstroTools

max_forks_repo_name: FatherGeoffHorton/AstroTools


text:
using BenchmarkTools
include("boid.jl")

# ## Running the model
function benchmark()
    n_steps = 100
    results = []
    for N in 100:100:400
        b = @benchmarkable step!(model, agent_step!, $n_steps) setup=(model =
        initialize_model(n_birds=$N,
        speed=1.0,
        cohere_factor=.25,
        separation=4.0,
        seperate_factor=.25,
        match_factor=.01,
        visual_distance=5.0,
        dims=(100,100)))

        j = run(b)
        push!(results, minimum(j.times)/1e9)
    end
    return results
end

print(benchmark())

ext: jl


text:
#==============================================================================#
# MarketplaceEntitlementService.jl
#
# This file is generated from:
# https://github.com/aws/aws-sdk-js/blob/master/apis/entitlement.marketplace-2017-01-11.normal.json
#==============================================================================#

__precompile__()

module MarketplaceEntitlementService

using AWSCore

"""
    using AWSSDK.MarketplaceEntitlementService.get_entitlements
    get_entitlements([::AWSConfig], arguments::Dict)
    get_entitlements([::AWSConfig]; ProductCode=, <keyword arguments>)

    using AWSCore.Services.entitlement_marketplace
    entitlement_marketplace([::AWSConfig], "GetEntitlements", arguments::Dict)
    entitlement_marketplace([::AWSConfig], "GetEntitlements", ProductCode=, <keyword arguments>)

# GetEntitlements Operation

GetEntitlements retrieves entitlement values for a given product. The results can be filtered based on customer identifier or product dimensions.

# Arguments

## `ProductCode = ::String` -- *Required*
Product code is used to uniquely identify a product in AWS Marketplace. The product code will be provided by AWS Marketplace when the product listing is created.

## `Filter = ::Dict{String,String}`
Filter is used to return entitlements for a specific customer or for a specific dimension. Filters are described as keys mapped to a lists of values. Filtered requests are *unioned* for each value in the value list, and then *intersected* for each filter key.

## `NextToken = ::String`
For paginated calls to GetEntitlements, pass the NextToken from the previous GetEntitlementsResult.

## `MaxResults = ::Int`
The maximum number of items to retrieve from the GetEntitlements operation. For pagination, use the NextToken field in subsequent calls to GetEntitlements.

# Returns

`GetEntitlementsResult`

# Exceptions

`InvalidParameterException`, `ThrottlingException` or `InternalServiceErrorException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/entitlement.marketplace-2017-01-11/GetEntitlements)
"""
@inline get_entitlements(aws::AWSConfig=default_aws_config(); args...) = get_entitlements(aws, args)

@inline get_entitlements(aws::AWSConfig, args) = AWSCore.Services.entitlement_marketplace(aws, "GetEntitlements", args)

@inline get_entitlements(args) = get_entitlements(default_aws_config(), args)

end # module MarketplaceEntitlementService

#==============================================================================#
# End of file
#==============================================================================#


text:
using Test
using Bijectors
using Bijectors: RationalQuadraticSpline

@testset "RationalQuadraticSpline" begin
    # Monotonic spline on '[-B, B]' with `K` intermediate knots/"connection points".
    d = 2; K = 3; B = 2;
    b_uv = RationalQuadraticSpline(randn(K), randn(K), randn(K - 1), B)
    b_mv = RationalQuadraticSpline(randn(d, K), randn(d, K), randn(d, K - 1), B)

    @testset "Constructor" begin
        # univariate
        b = b_uv
        @test b.widths[1] ≈ -B
        @test b.widths[end] ≈ B
        @test b.widths[1] == minimum(b.widths)
        @test b.widths[end] == maximum(b.widths)
        @test b.heights[1] == minimum(b.heights)
        @test b.heights[end] == maximum(b.heights)
        @test all(b.derivatives .> 0)
        @test b.derivatives[1] == b.derivatives[end] == 1

        # multivariate
        b = b_mv
        @test all(b.widths[:, 1] .≈ -B)
        @test all(b.widths[:, end] .≈ B)
        @test all(b.widths[:, 1] .== minimum.(eachrow(b.widths)))
        @test all(b.widths[:, end] .== maximum.(eachrow(b.widths)))
        @test all(b.heights[:, 1] .== minimum.(eachrow(b.heights)))
        @test all(b.heights[end] .== maximum.(eachrow(b.heights)))
        @test all(all(b.derivatives .> 0))
        @test all(b.derivatives[:, 1] .== b.derivatives[:, end] .== 1)
    end

    @testset "Evaluation" begin
        # univariate
        b = b_uv

        # Inside of domain
        x = 0.5
        test_bijector(b, [-x, x])

        # Outside of domain
        x = 5.
        test_bijector(b, [-x, x], [-x, x], [0., 0.])

        # multivariate
        b = b_mv

        # Inside of domain
        x = [-0.5, 0.5]
        x = hcat(x, -x, x) # batch
        test_bijector(b, x)

        # Outside of domain
        x = [-5., 5.]
        x = hcat(x, -x, x) # batch
        test_bijector(b, x, x, zeros(size(x, 2)))
    end
end


text:
# reversable padding functions for use in block cypher encryption
module Padding

include("Rand.jl")

# all use of these functions outside of this file are via:
function pad!(padding, block_size, data :: Array{Uint8})
    n = block_size - (length(data) % block_size)
    append!(data, make_padding(padding, n == 0 ? block_size : n))
end
function pad!(padding, block_size, data :: Array{Uint8}, length)
    n = block_size - (length % block_size)
    data[length+1:end] = make_padding(padding, n == 0 ? block_size : n)
end

function padded_size(block_size, data)
    len = length(data)
    n = block_size - (len % block_size)
    len + n
end

function unpad!(padding, data :: Array{Uint8})
    resize!(data, length(data) - padding_length(padding, data))
end

# which as you can see, requires everying padding sceme to implement both
# make_padding and padding_length

immutable NoPadding
end

function make_padding(::Type{NoPadding}, len)
    []
end

function padding_length(::Type{NoPadding}, data :: Array{Uint8})
    0
end

immutable BytePadding
end

function make_padding(::Type{BytePadding}, len)
    out = zeros(Uint8, len)
    out[1] = 1
    out
end

function padding_length(::Type{BytePadding}, data :: Array{Uint8})
    len = length(data)
    ct = 0
    while data[len - ct] == 0
        ct += 1
    end

    @assert data[len - ct] == 1
    ct + 1
end

immutable AnsiX923
end

function make_padding(::Type{AnsiX923}, len)
    out = zeros(Uint8, len)
    out[end] = len
    out
end

function padding_length(::Type{AnsiX923}, data :: Array{Uint8})
    data[end]
end

immutable Iso10126
end

function make_padding(::Type{Iso10126}, len)
    out = Rand.cryptographic_rand(len)
    out[end] = len
    out
end

function padding_length(::Type{Iso10126}, data :: Array{Uint8})
    data[end]
end

immutable Pkcs7
end

function make_padding(::Type{Pkcs7}, len)
    [len for i=1:len]
end

function padding_length(::Type{Pkcs7}, data :: Array{Uint8})
    data[end]
end

immutable IsoIec7816_4
end

function make_padding(::Type{IsoIec7816_4}, len)
    out = zeros(Uint8, len)
    out[1] = 0x08
    out
end

function padding_length(::Type{IsoIec7816_4}, data :: Array{Uint8})
    len = length(data)
    ct = 0
    while data[len - ct] == 0
        ct += 1
    end

    @assert data[len - ct] == 0x08
    ct + 1
end

export pad!, unpad!, padded_size, NoPadding, BytePadding, AnsiX923, Iso10126, Pkcs7, IsoIec7816_4

end # module Padding


text:
VERSION >= v"1.0.0" && __precompile__(true)

module AdventOfCode2017

using DelimitedFiles
using LinearAlgebra
using Base.Meta

files = [x for x in readdir(joinpath(@__DIR__, "puzzles")) if splitext(x)[2]==".jl"]
include.(joinpath.(@__DIR__, "puzzles/", files))
for file in first.(splitext.(files))
    func = Symbol(file)
    @eval export $func
end

end


text:
module nrtl

include("general.jl")
using Unitful
using NLsolve

∑(x)=sum(x)

τ(T::Quantity,Δg,i::Int,j::Int) = Δg[i,j]/(1.98721u"cal/(mol*K)"*T)
τ(T::Float64,Δg,i::Int,j::Int) = ustrip(u"cal/mol",Δg[i,j])/(1.98721*T)
G(T,Δg,α,i::Int,j::Int) = exp(-α[i,j]*τ(T,Δg,i,j))
γ(T,x,Δg,α,i::Int,n::Int) = exp(∑(x[j]*G(T,Δg,α,j,i)*τ(T,Δg,j,i) for j=1:n)/∑(x[k]*G(T,Δg,α,k,i) for k=1:n)+ ∑((x[j]*G(T,Δg,α,i,j)/∑(x[k]*G(T,Δg,α,k,j) for k=1:n))*(τ(T,Δg,i,j)-∑(x[m]*G(T,Δg,α,m,j)*τ(T,Δg,m,j) for m=1:n)/∑(x[k]*G(T,Δg,α,k,j) for k=1:n)) for j=1:n))
P(T,x,CA,Δg,α,i::Int,n::Int) = x[i]*γ(T,x,Δg,α,i,n)*general.P_Antoine(CA,T,i)
Pₜₒₜ(T,x,CA,Δg,α,n::Int) = ∑(P(T,x,CA,Δg,α,i,n) for i=1:n)
y_calc(T,x,CA,Δg,α,i::Int,n::Int)=P(T,x,CA,Δg,α,i,n)/Pₜₒₜ(T,x,CA,Δg,α,n)

"""
    pxy(n,T,CA,Δg,α,xx;uni=u"Torr")

**Los campos de entrada son:**
- `n :: Int` Es el numero de componentes
- `T :: Float` Es la temperatura del sistema
- `CA :: n×3 Array{Float,2}` Constantes de Antoine para los `n` componentes
- `Δg :: n×n Array{Float,2}` Parametros de interacción entre los componentes en grados kelvin
- `α :: n×n Array{Float,2}` Parametro de no aleatoriedad
- `xx :: StepRange ó Vector` Son los puntos en el liquido donde se buscara el equilibrio con el vapor
- `uni` Son las unidades en las que se desea el resultado, por defecto esta en Torr

**Salida:**

Regresa las fracciones de liquido `x`,
las fracciones del vapor `y` y la presión `P` en Torr (`m` es la longitud de `xx`).

- `x :: n×m Array{Float64,2}` liquido
- `y :: n×m Array{Float64,2}` vapor
- `P :: m-element Array{Float64,1}` presión
"""
function pxy(n::Int,T,CA::Array,Δg::Array,α,xx;uni=u"Torr")
    T = T|> u"K"
    xx = general.nc(xx...)
    println("ddd",xx)
    m = size(xx,2)
    y = Array{Float64}(undef,n,m)
    P = Vector{Float64}(undef,m)*uni
    for i=1:m
        x = xx[:,i]
        for j=1:n
            y[j,i] = y_calc(T,x,CA,Δg,α,j,n)
        end
        P[i] = Pₜₒₜ(T,x,CA,Δg,α,n) |> uni
    end

    return xx,y,P
end

"""
    txy(n,Pobj,CA,Δg,α,xx;uni=u"K")

**Los campos de entrada son:**
- `n :: Int` Es el numero de componentes
- `Pobj :: Float` Es la presión del sistema
- `CA :: n×3 Array{Float,2}` Constantes de Antoine para los `n` componentes
- `Δg :: n×n Array{Float,2}` Parametros de interacción entre los componentes en grados kelvin
- `α :: n×n Array{Float,2}` Parametro de no aleatoriedad
- `xx :: StepRange ó Vector` Son los puntos en el liquido donde se buscara el equilibrio con el vapor
- `uni` Son las unidades en las que se desea el resultado, por defecto esta en grados kelvin

**Salida:**

Regresa las fracciones de liquido `x`,
las fracciones del vapor `y` y la temperatura `T` en grados kelvin (`m` es la longitud de `xx`).

- `x :: n×m Array{Float64,2}` liquido
- `y :: n×m Array{Float64,2}` vapor
- `T :: m-element Array{Float64,1}` temperatura
"""
function txy(n::Int,Pobj,CA::Array,Δg::Array,α,xx;uni=u"K")
    Pobj = Pobj |> u"Torr"
    xx = general.nc(xx...)
    m = size(xx,2)
    y = Array{Float64}(undef,n,m)
    T = Vector{Float64}(undef,m)*u"K"
    Tᵢ = 300.0
    #temperatura inicial
    for i=1:n
        if xx[i,1]==1
            T[1] = general.T_Antoine(CA,Pobj,i)
            Tᵢ  = ustrip(u"K",T[1])
        end
    end
    for i = 1:m
        x = xx[:,i]
        P2ₜₒₜ(T) = Pₜₒₜ(T,x,CA,Δg,α,n) - ustrip(u"Torr",Pobj)
        #println(i)
        Pₛₒₗ = nlsolve(n_ary(P2ₜₒₜ),[Tᵢ])
        Tᵢ = Pₛₒₗ.zero[1]
        T[i] = Tᵢ*u"K"
        for j = 1:n
            y[j,i] = y_calc(T[i],x,CA,Δg,α,j,n)
        end
    end

    return xx,y,T .|> uni

end
end

max_stars_repo_name: EmilioAlvizo/TermoQuimica

max_issues_repo_name: EmilioAlvizo/TermoQuimica

max_forks_repo_name: EmilioAlvizo/TermoQuimica


text:
#!/usr/bin/env bash
    #=
    exec julia --project="$(realpath $(dirname $0))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \
    "${BASH_SOURCE[0]}" "$@"
    =#

using StatsPlots, ProgressMeter, CSV, DataFrames

include(joinpath(dirname(dirname(@__FILE__)), "src", "CodingTheory.jl"))
using .CodingTheory

stop_at = parse(BigInt, ARGS[1])

function graphing_time(stop_at::Integer)
    stop_q_at = 8
    stop_n_at = 8
    num_of_datapoints = stop_at
    data = Matrix{Float64}(undef, num_of_datapoints, 3) # needs 3 columns: q, n, and time

    @showprogress for i in 1:num_of_datapoints
        q = rand(1:stop_q_at)
        n = rand(1:stop_n_at)
        time_took = @elapsed get_all_words(q, n)
        data[i, :] .= [q, n, time_took]
    end

    save_path = joinpath(dirname(@__DIR__), "other", "get_codewords_time_analysis,n=$stop_n_at,m=$stop_q_at,i=$num_of_datapoints.pdf")
    theme(:solarized)

    plot_points = scatter(data[:,1:2], data[:,3], fontfamily = font("Times"), ylabel = "Time elapsed during calculation [seconds]", clabel = "Value of q or n", label = ["q" "n"], smooth = true)#, xlims = (0, stop_n_at)) # smooth = true

    savefig(plot_points, save_path)
    println("Plot saved at $save_path.")

    D = DataFrame(q = Number[], n = Number[], time = Number[])

    for i in eachrow(data)
        push!(D, i)
    end

	data_file_desktop = joinpath(homedir(), "Desktop", "get_codewords_time_analysis,n=$stop_n_at,m=$stop_q_at,i=$num_of_datapoints.csv")
	data_file_other = joinpath(dirname(dirname(@__FILE__)), "other", "get_codewords_time_analysis,n=$stop_n_at,m=$stop_q_at,i=$num_of_datapoints.csv")

    CSV.write(data_file_desktop, D)
    CSV.write(data_file_other, D)
	println("Wrote data to $(data_file_other).")
end

graphing_time(stop_at)


text:
using DashBootstrapComponents, DashHtmlComponents

input_groups = html_div([
    dbc_inputgroup(
        [dbc_inputgroupaddon(dbc_radiobutton(), addon_type="prepend"), dbc_input()],
        className="mb-3",
    ),
    dbc_inputgroup([
        dbc_inputgroupaddon(dbc_checkbox(), addon_type="prepend"),
        dbc_input(),
    ]),
]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct TGTypeII <: TurbineGov
        R::Float64
        T1::Float64
        T2::Float64
        τ_min::Float64
        τ_max::Float64
        P_ref::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int64
        internal::InfrastructureSystemsInternal
    end

Parameters of a Turbine Governor Type II.

# Arguments
- `R::Float64`: Droop parameter, validation range: (0, nothing)
- `T1::Float64`: Transient gain time constant, validation range: (0, nothing)
- `T2::Float64`: Power fraction time constant, validation range: (0, nothing)
- `τ_min::Float64`: Min Power into the Governor, validation range: (0, nothing)
- `τ_max::Float64`: Max Power into the Governor, validation range: (0, nothing)
- `P_ref::Float64`: Reference Power Set-point, validation range: (0, nothing)
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`
- `n_states::Int64`
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct TGTypeII <: TurbineGov
    "Droop parameter"
    R::Float64
    "Transient gain time constant"
    T1::Float64
    "Power fraction time constant"
    T2::Float64
    "Min Power into the Governor"
    τ_min::Float64
    "Max Power into the Governor"
    τ_max::Float64
    "Reference Power Set-point"
    P_ref::Float64
    ext::Dict{String, Any}
    states::Vector{Symbol}
    n_states::Int64
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function TGTypeII(R, T1, T2, τ_min, τ_max, P_ref=1.0, ext=Dict{String, Any}(), )
    TGTypeII(R, T1, T2, τ_min, τ_max, P_ref, ext, [:xg], 1, InfrastructureSystemsInternal(), )
end

function TGTypeII(; R, T1, T2, τ_min, τ_max, P_ref=1.0, ext=Dict{String, Any}(), )
    TGTypeII(R, T1, T2, τ_min, τ_max, P_ref, ext, )
end

# Constructor for demo purposes; non-functional.
function TGTypeII(::Nothing)
    TGTypeII(;
        R=0,
        T1=0,
        T2=0,
        τ_min=0,
        τ_max=0,
        P_ref=0,
        ext=Dict{String, Any}(),
    )
end

"""Get TGTypeII R."""
get_R(value::TGTypeII) = value.R
"""Get TGTypeII T1."""
get_T1(value::TGTypeII) = value.T1
"""Get TGTypeII T2."""
get_T2(value::TGTypeII) = value.T2
"""Get TGTypeII τ_min."""
get_τ_min(value::TGTypeII) = value.τ_min
"""Get TGTypeII τ_max."""
get_τ_max(value::TGTypeII) = value.τ_max
"""Get TGTypeII P_ref."""
get_P_ref(value::TGTypeII) = value.P_ref
"""Get TGTypeII ext."""
get_ext(value::TGTypeII) = value.ext
"""Get TGTypeII states."""
get_states(value::TGTypeII) = value.states
"""Get TGTypeII n_states."""
get_n_states(value::TGTypeII) = value.n_states
"""Get TGTypeII internal."""
get_internal(value::TGTypeII) = value.internal

"""Set TGTypeII R."""
set_R!(value::TGTypeII, val::Float64) = value.R = val
"""Set TGTypeII T1."""
set_T1!(value::TGTypeII, val::Float64) = value.T1 = val
"""Set TGTypeII T2."""
set_T2!(value::TGTypeII, val::Float64) = value.T2 = val
"""Set TGTypeII τ_min."""
set_τ_min!(value::TGTypeII, val::Float64) = value.τ_min = val
"""Set TGTypeII τ_max."""
set_τ_max!(value::TGTypeII, val::Float64) = value.τ_max = val
"""Set TGTypeII P_ref."""
set_P_ref!(value::TGTypeII, val::Float64) = value.P_ref = val
"""Set TGTypeII ext."""
set_ext!(value::TGTypeII, val::Dict{String, Any}) = value.ext = val
"""Set TGTypeII states."""
set_states!(value::TGTypeII, val::Vector{Symbol}) = value.states = val
"""Set TGTypeII n_states."""
set_n_states!(value::TGTypeII, val::Int64) = value.n_states = val
"""Set TGTypeII internal."""
set_internal!(value::TGTypeII, val::InfrastructureSystemsInternal) = value.internal = val


text:
using ContinuousTimePolicyGradients
using DiffEqFlux, ComponentArrays, LinearAlgebra, JLD2, OrdinaryDiffEq
using Plots

function main(maxiters_1::Int, maxiters_2::Int, Δt_save::Float32; p_NN_0 = nothing, k_a_val = 10.0)
    # model + problem parameters
    (aₐ, aₙ, bₙ, cₙ, dₙ, aₘ, bₘ, cₘ, dₘ) = Float32.([-0.3, 19.373, -31.023, -9.717, -1.948, 40.44, -64.015, 2.922, -11.803])
    (m, I_yy, S, d, ωₐ, ζₐ) = Float32.([204.02, 247.439, 0.0409, 0.2286, 150.0, 0.7])
    (g, ρ₀, H, γₐ, Rₐ, T₀, λ) = Float32.([9.8, 1.225, 8435.0, 1.4, 286.0, 288.15, 0.0065])
    (a_max, α_max, δ_max, δ̇_max, q_max, M_max, h_max) = Float32.([100.0, deg2rad(30), deg2rad(25), 1.5, deg2rad(60), 4, 11E3])

    (k_a, k_q, k_δ, k_δ̇, k_R) = Float32.([k_a_val, 0.0, 0.0, 0.1, 1E-3])

    # dynamic model
    dim_x = 7
    function dynamics_plant(t, x, u)
        (h, V, α, q, θ, δ, δ̇) = x
        δ_c = u[1]

        ρ = ρ₀ * exp(-h / H)
        Vₛ = sqrt(γₐ * Rₐ * (T₀ - λ * h))
        M = V / Vₛ
        γ = θ - α
        Q = 0.5f0 * ρ * V^2

        C_A = aₐ
        C_N = aₙ * α^3 + bₙ * α * abs(α) + cₙ * (2.0f0 - M / 3.0f0) * α + dₙ * δ
        C_M = aₘ * α^3 + bₘ * α * abs(α) + cₘ * (-7.0f0 + 8.0f0 * M / 3.0f0) * α + dₘ * δ

        α̇ = Q * S / m / V * (C_N * cos(α) - C_A * sin(α)) + g / V * cos(γ) + q

        dx = [V * sin(γ);
            Q * S / m * (C_N * sin(α) + C_A * cos(α)) - g * sin(γ);
            α̇;
            Q * S * d / I_yy * C_M;
            q;
            δ̇;
            -ωₐ^2 * (δ - δ_c) - 2.0f0 * ζₐ * ωₐ * δ̇]
        return dx
    end

    dim_x_c = 2
    function dynamics_controller(t, x_c, y, r, p_NN, policy_NN)
        (A_z, h, V, M, α, q, γ) = y
        A_z_cmd  = r[1]
        x_int    = x_c[1]
        x_ref    = x_c[2]

        y_NN = (K_A, K_I, K_R) = policy_NN([α / α_max; M / M_max; h / h_max], p_NN)

        # dx_ref  = [-36.6667f0 -13.8889f0; 8.0f0 0.0f0] * x_ref + [4.0f0; 0.0f0] * A_z_cmd
        # A_z_ref = [-1.0083f0 3.4722f0] * x_ref
        dx_ref = (A_z_cmd - x_ref) / 0.2f0
        A_z_ref = x_ref

        δ_trim = -(aₘ * α^3 + bₘ * α * abs(α) + cₘ * (-7.0f0 + 8.0f0 * M / 3.0f0) * α) / dₘ

        dx_c = [K_A * (A_z_cmd - A_z) + q + (A_z + g * cos(γ)) / V;
                dx_ref]
        u    = [K_I * x_int + K_R * q + δ_trim;
                A_z_ref]
        return dx_c, u, y_NN
    end

    function dynamics_sensor(t, x)
        (h, V, α, q, θ, δ, _) = x

        ρ = ρ₀ * exp(-h / H)
        Vₛ = sqrt(γₐ * Rₐ * (T₀ - λ * h))
        M = V / Vₛ
        γ = θ - α
        Q = 0.5f0 * ρ * V^2

        C_A = aₐ
        C_N = aₙ * α^3 + bₙ * α * abs(α) + cₙ * (2.0f0 - M / 3.0f0) * α + dₙ * δ

        A_z = Q * S / m * (C_N * cos(α) - C_A * sin(α))
        y   = [A_z;
                h;
                V;
                M;
                α;
                q;
                γ]
        return y
    end

    # cost definition
    function cost_running(t, x, y, u, r)
        q       = x[4]
        δ̇       = x[7]
        a_z     = y[1]
        δ_c     = u[1]
        a_z_ref = u[2]
        a_z_cmd = r[1]
        return k_a * ((a_z - a_z_ref) / (exp(-2.0f0*a_z_cmd) + a_z_cmd))^2 + k_δ̇ * (δ̇ / δ̇_max)^2  + k_δ * (δ_c / δ_max)^2 + k_q * (q / q_max)^2
    end

    function cost_terminal(x_f, r)
        # a_z_cmd  = r[1]
        # y        = dynamics_sensor(3.0f0, x_f)
        # a_z      = y[1]
        # return ((a_z - a_z_cmd) / (1.0f0 + a_z_cmd))^2
        return 0.0f0
    end

    function cost_regularisor(p_NN)
        return k_R * norm(p_NN)^2
        # return 0.0f0
    end

    # NN construction
    dim_NN_hidden = 64
    dim_NN_input  = 3
    dim_K = 3
    K_lb  = zeros(Float32, 3)
    K_ub  = Float32[4,0.2,4]
    policy_NN = FastChain(
        FastDense(dim_NN_input,  dim_NN_hidden, tanh),
        FastDense(dim_NN_hidden, dim_NN_hidden, tanh),
        FastDense(dim_NN_hidden, dim_K),
        (x, p) -> (K_ub - K_lb) .* σ.(x) .+ K_lb
    )

    # scenario definition
    ensemble = [ (; x₀ = Float32[h₀; V₀; zeros(5)], r = Float32[A_z_cmd])
                     for h₀      = 5E3 #:1E3:8E3
                     for V₀      = 7E2 #:1E2:9E2
                     for A_z_cmd = 0:2E1:1E2 ]
    t_span = Float32.((0.0, 3.0))
    t_save = t_span[1]:Δt_save:t_span[2]

    scenario = (; ensemble = ensemble, t_span = t_span, t_save = t_save, dim_x = dim_x, dim_x_c = dim_x_c)

    # NN training
    (result, fwd_ensemble_sol, loss_history) = CTPG_train(dynamics_plant, dynamics_controller, dynamics_sensor, cost_running, cost_terminal, cost_regularisor, policy_NN, scenario; sense_alg = InterpolatingAdjoint(autojacvec = ReverseDiffVJP(true)), ensemble_alg = EnsembleThreads(), maxiters_1 = maxiters_1, maxiters_2 = maxiters_2, opt_2 = BFGS(initial_stepnorm = 0.0001), i_nominal = 1, p_NN_0 = p_NN_0, progress_plot = false)

    return result, policy_NN, fwd_ensemble_sol, loss_history
end

## execute optimisation and simulation
@time (result, policy_NN, fwd_ensemble_sol, loss_history) = main(1000, 500, 0.01f0; k_a_val = 200.0)

## save results
# p_NN_prev = result.u
# @time (result, policy_NN, fwd_ensemble_sol, loss_history) = main(1, 1000, 0.01f0; k_a_val = 50.0, p_NN_0 = p_NN_prev)

jldsave("DS_autopilot_13.jld2"; result, fwd_ensemble_sol, loss_history)

# plot results
# (fwd_ensemble_sol, loss_history) = load(".jld2", "fwd_ensemble_sol", "loss_history")

x_names = ["\$h\$" "\$V\$" "\$\\alpha\$" "\$q\$" "\$\\theta\$" "\$\\delta\$" "\$\\dot{\\delta}\$"]
vars_x = 1:6 # [1,2,3, (1,2), (2,3)]
u_names = ["\$\\delta_{c}\$"]
vars_u = 1
y_names = ["\$a_{z}\$"]
vars_y = 1
y_NN_names = ["\$K_A\$" "\$K_I\$" "\$K_R\$"]
vars_y_NN = 1:3

(f_x, f_u, f_y, f_y_NN, f_L) = view_result([], fwd_ensemble_sol, loss_history; x_names = x_names, vars_x = vars_x, u_names = u_names, vars_u = vars_u, y_names = y_names, vars_y = vars_y, y_NN_names = y_NN_names, vars_y_NN = vars_y_NN, linealpha = 0.6)

(a_max, α_max, δ_max, δ̇_max, q_max, M_max, h_max) = Float32.([100.0, deg2rad(30), deg2rad(25), 1.5, deg2rad(60), 4, 11E3])

h = 5000.0
α_list = 0:1E-3:deg2rad(45)
M_list = 2.0:0.1:3.0

func_K_A(α, M) = policy_NN([α / α_max; M / M_max; h / h_max], result.u)[1]
func_K_I(α, M) = policy_NN([α / α_max; M / M_max; h / h_max], result.u)[2]
func_K_R(α, M) = policy_NN([α / α_max; M / M_max; h / h_max], result.u)[3]
f_K_A = plot(α_list, M_list, func_K_A, st=:surface, label = :false, title = "K_A", xlabel = "\$\\alpha\$", ylabel = "\$M\$")
f_K_I = plot(α_list, M_list, func_K_I, st=:surface, label = :false, title = "K_I", xlabel = "\$\\alpha\$", ylabel = "\$M\$")
f_K_R = plot(α_list, M_list, func_K_R, st=:surface, label = :false, title = "K_R", xlabel = "\$\\alpha\$", ylabel = "\$M\$")
display(f_K_A)
display(f_K_I)
display(f_K_R)


text:
"""
    Cartpole <: Environment

    classic system with one translational (cart) and one rotational (pole) degree of freedom
    https://underactuated.mit.edu/acrobot.html#cart_pole
"""
struct Cartpole end

function cartpole(;
    representation=:minimal,
    timestep=0.05,
    gravity=[0.0; 0.0; -9.81],
    seed=1,
    vis=Visualizer(),
    info=nothing,
    name=:robot,
    control_scaling=Diagonal(ones(1)),
    opts_step=SolverOptions(),
    opts_grad=SolverOptions(),
    T=Float64)

    mechanism = get_cartpole(
        timestep=timestep,
        gravity=gravity)
    initialize_cartpole!(mechanism)

    if representation == :minimal
        nx = minimal_dimension(mechanism)
    elseif representation == :maximal
        nx = maximal_dimension(mechanism)
    end
    nu = 1
    no = nx

    aspace = BoxSpace(nu,
        low=(-1.0e8 * ones(nu)),
        high=(1.0e8 * ones(nu)))
    ospace = BoxSpace(no,
        low=(-Inf * ones(no)),
        high=(Inf * ones(no)))

    rng = MersenneTwister(seed)

    z = get_maximal_state(mechanism)
    x = representation == :minimal ? maximal_to_minimal(mechanism, z) : z

    fx = zeros(nx, nx)
    fu = zeros(nx, nu)

    u_prev = zeros(nu)
    control_mask = [1.0 0.0]

    build_robot(mechanism, vis=vis, name=name)

    TYPES = [Cartpole, T, typeof(mechanism), typeof(aspace), typeof(ospace), typeof(info)]
    env = Environment{TYPES...}(mechanism, representation, aspace, ospace,
        x, fx, fu,
        u_prev, control_mask' * control_scaling,
        nx, nu, no,
        info,
        [rng], vis,
        opts_step, opts_grad)

    return env
end

function cartpole_nominal_max(;
    pendulum_length=1.0)

    # initial state
    x1b1 = [0.0; 0.0; 0.0]
    v1b1 = [0.0; 0.0; 0.0]
    q1b1 = [1.0; 0.0; 0.0; 0.0]
    ω1b1 = [0.0; 0.0; 0.0]
    z1b1 = [x1b1; v1b1; q1b1; ω1b1]

    x1b2 = [0.0; 0.0; -0.5 * pendulum_length]
    v1b2 = [0.0; 0.0; 0.0]
    q1b2 = [1.0; 0.0; 0.0; 0.0]
    ω1b2 = [0.0; 0.0; 0.0]
    z1b2 = [x1b2; v1b2; q1b2; ω1b2]

    z1 = [z1b1; z1b2]
end

function cartpole_goal_max(;
    pendulum_length=1.0)

    # target state
    xTb1 = [0.0; 0.0; 0.0]
    vTb1 = [0.0; 0.0; 0.0]
    qTb1 = [1.0; 0.0; 0.0; 0.0]
    ωTb1 = [0.0; 0.0; 0.0]
    zTb1 = [xTb1; vTb1; qTb1; ωTb1]

    xTb2 = [0.0; 0.0; 0.5 * pendulum_length]
    vTb2 = [0.0; 0.0; 0.0]
    qTb2 = [0.0; 1.0; 0.0; 0.0]
    ωTb2 = [0.0; 0.0; 0.0]
    zTb2 = [xTb2; vTb2; qTb2; ωTb2]

    zT = [zTb1; zTb2]
end

function visualize(env::Environment{Cartpole}, traj::Vector{Vector{T}};
    axes=false,
    grid=false,
    ee_traj=false) where T

	@assert size(traj[1]) == size(env.state)
    storage = generate_storage(env.mechanism, [env.representation == :minimal ? minimal_to_maximal(env.mechanism, x) : x for x in traj])
    visualize(env.mechanism, storage,
        vis=env.vis)
	set_camera!(env.vis, zoom=1.0, cam_pos=[2,0,0])
    setvisible!(env.vis["/Axes"], axes)
    setvisible!(env.vis["/Grid"], grid)
    # slider
    l2 = GeometryBasics.Cylinder(Point3f0(0.0, -5.0, 0.0),
		Point3f0(0.0, 5.0, 0.0),
		convert(Float32, 0.0125))
	setobject!(env.vis["slider"], l2, MeshPhongMaterial(color=Colors.RGBA(0.0, 0.0, 0.0, 1.0)))
end

function ghost(env::Environment{Cartpole}, z_sol;
    timesteps=[t for t = 1:length(z_sol)],
    line=false)

    for t in timesteps
        name = Symbol("robot_$t")
        build_robot(env.vis, env.mechanism,
            name=name,
            color=(t == length(z_sol) ? cyan : cyan_light))
        z = z_sol[t]
        if t == length(z_sol)
            z[1] -= 0.05
            z[14] -= 0.05
        end
        set_robot(env.vis, env.mechanism, z,
            name=name)
    end

	set_camera!(env.vis,
        zoom=1.0,
        cam_pos=[2,0,0])
    setvisible!(env.vis[:robot], false)
    l2 = GeometryBasics.Cylinder(Point3f0(0.0, -5.0, 0.0),
            Point3f0(0.0, 5.0, 0.0),
            convert(Float32, 0.0125))
    setobject!(env.vis["slider"], l2, MeshPhongMaterial(color=Colors.RGBA(0.0, 0.0, 0.0,1.0)))
    setvisible!(env.vis["/Axes"], false)
    setvisible!(env.vis["/Grid"], false)

    if line
        line_mat = LineBasicMaterial(
            color=color=RGBA(51.0 / 255.0, 1.0, 1.0, 1.0),
            linewidth=25.0)
        points = Vector{Point{3,Float64}}()
        for (i, xt) in enumerate(z_sol)
            z_pendulum = xt[13 .+ (1:13)]
            x = z_pendulum[1:3]
            q = Quaternion(z_pendulum[6 .+ (1:4)]...)
            k = x + vector_rotate([0.0; 0.0; -0.5], q)

            push!(points, Point(k[1], k[2], k[3]))

            setobject!(env.vis["ee_vertex_$i"], GeometryBasics.Sphere(Point3f0(0),
                convert(Float32, 0.001)),
                MeshPhongMaterial(color=RGBA(1.0, 153.0 / 255.0, 51.0 / 255.0, 1.0)))
                settransform!(env.vis["ee_vertex_$i"], MeshCat.Translation(points[i]))
        end
        setobject!(env.vis[:ee_traj], MeshCat.Line(points, line_mat))
    end

    open(env.vis)
end


text:
#############################################################################
# variable.jl
# Defines Variable, which is a subtype of AbstractExpr
#############################################################################

export Variable, Semidefinite, ComplexVariable, HermitianSemidefinite
export vexity, evaluate, sign, conic_form!, fix!, free!

mutable struct Variable <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    value::ValueOrNothing
    size::Tuple{Int, Int}
    vexity::Vexity
    sign::Sign
    sets::Array{Symbol,1}

    function Variable(size::Tuple{Int, Int}, sign::Sign=NoSign(), sets::Symbol...)
        this = new(:variable, 0, nothing, size, AffineVexity(), sign, Symbol[sets...])
        this.id_hash = objectid(this)
        return this
    end

    Variable(m::Int, n::Int, sign::Sign=NoSign(), sets::Symbol...) = Variable((m,n), sign, sets...)
    Variable(sign::Sign, sets::Symbol...) = Variable((1, 1), sign, sets...)
    Variable(sets::Symbol...) = Variable((1, 1), NoSign(), sets...)
    Variable(size::Tuple{Int, Int}, sets::Symbol...) = Variable(size, NoSign(), sets...)
    Variable(size::Int, sign::Sign=NoSign(), sets::Symbol...) = Variable((size, 1), sign, sets...)
    Variable(size::Int, sets::Symbol...) = Variable((size, 1), sets...)
end

Semidefinite(m::Integer) = Variable((m, m), :Semidefinite)
function Semidefinite(m::Integer, n::Integer)
    if m == n
        return Variable((m, m), :Semidefinite)
    else
        error("Semidefinite matrices must be square")
    end
end

ComplexVariable(m::Int, n::Int, sets::Symbol...) = Variable((m, n), ComplexSign(), sets...)
ComplexVariable(sets::Symbol...) = Variable((1, 1), ComplexSign(), sets...)
ComplexVariable(size::Tuple{Int, Int}, sets::Symbol...) = Variable(size, ComplexSign(), sets...)
ComplexVariable(size::Int, sets::Symbol...) = Variable((size, 1), ComplexSign(), sets...)
HermitianSemidefinite(m::Integer) = ComplexVariable((m, m), :Semidefinite)
function HermitianSemidefinite(m::Integer, n::Integer)
    if m == n
        return ComplexVariable((m, m), :Semidefinite)
    else
        error("HermitianSemidefinite matrices must be square")
    end
end

function vexity(x::Variable)
    return x.vexity
end

function evaluate(x::Variable)
    return x.value === nothing ? error("Value of the variable is yet to be calculated") : x.value
end

function sign(x::Variable)
    return x.sign
end

function real_conic_form(x::Variable)
    vec_size = length(x)
    return sparse(1.0I, vec_size, vec_size)
end

function imag_conic_form(x::Variable)
    vec_size = length(x)
    if x.sign == ComplexSign()
        return im*sparse(1.0I, vec_size, vec_size)
    else
        return spzeros(vec_size, vec_size)
    end
end

# fix variables to hold them at their current value, and free them afterwards
function fix!(x::Variable)
    x.value === nothing && error("This variable has no value yet; cannot fix value to nothing!")
    x.vexity = ConstVexity()
    x
end
function fix!(x::Variable, v::AbstractArray)
    size(x) == size(v) || throw(DimensionMismatch("Variable and value sizes do not match!"))
    x.value = sign(x) == ComplexSign() ? convert(Array{ComplexF64}, v) : convert(Array{Float64}, v)
    fix!(x)
end

function fix!(x::Variable, v::AbstractVector)
    size(x, 2) == 1 || throw(DimensionMismatch("Cannot set value of a variable of size $(size(x)) to a vector"))
    size(x, 1) == length(v) || throw(DimensionMismatch("Variable and value sizes do not match!"))
    x.value = sign(x) == ComplexSign() ? convert(Array{ComplexF64}, v) : convert(Array{Float64}, v)
    fix!(x)
end

function fix!(x::Variable, v::Number)
    size(x) == (1,1) || throw(DimensionMismatch("Variable and value sizes do not match!"))
    x.value = sign(x) == ComplexSign() ? convert(ComplexF64, v) : convert(Float64, v)
    fix!(x)
end

function free!(x::Variable)
    x.vexity = AffineVexity()
    x
end


text:
module JET

let
    README = normpath(dirname(@__DIR__), "README.md")
    include_dependency(README)
    @doc read(README, String) JET
end

# not sure why, but a benchmark showed this is faster
Base.Experimental.@optlevel 1

const CC = Core.Compiler

@static isdefined(CC, :AbstractInterpreter) || throw(ErrorException("JET.jl only works with Julia versions 1.6 and higher"))

# imports
# =======

# `AbstractAnalyzer`
import .CC:
    # abstractinterpreterinterface.jl
    InferenceParams,
    OptimizationParams,
    get_world_counter,
    lock_mi_inference,
    unlock_mi_inference,
    add_remark!,
    may_optimize,
    may_compress,
    may_discard_trees,
    # jetcache.jl
    code_cache,
    # haskey,
    # get,
    # getindex,
    # setindex!,
    get_inference_cache,
    cache_lookup,
    # push!,
    # tfuncs.jl
    builtin_tfunction,
    return_type_tfunc,
    # abstractinterpretation.jl
    abstract_call_gf_by_type,
    abstract_call_method_with_const_args,
    abstract_call_method,
    abstract_eval_special_value,
    abstract_eval_value,
    abstract_eval_statement,
    # typeinfer.jl
    typeinf,
    _typeinf,
    finish,
    typeinf_edge,
    # JET.jl
    InferenceState

# `ConcreteInterpreter`
import JuliaInterpreter:
    # virtualprocess.jl
    step_expr!,
    evaluate_call_recurse!,
    handle_err

# usings
# ======
# TODO: really use `using` instead

import Core:
    Builtin,
    Intrinsics,
    IntrinsicFunction,
    MethodMatch,
    LineInfoNode,
    CodeInfo,
    CodeInstance,
    MethodInstance,
    NewvarNode,
    GlobalRef,
    GotoNode,
    GotoIfNot,
    Const,
    SSAValue,
    SlotNumber,
    Argument,
    Slot,
    ReturnNode,
    SimpleVector,
    svec

import .CC:
    AbstractInterpreter,
    NativeInterpreter,
    InferenceResult,
    InternalCodeCache,
    OptimizationState,
    WorldRange,
    WorldView,
    Bottom,
    NOT_FOUND,
    MethodMatchInfo,
    UnionSplitInfo,
    MethodLookupResult,
    VarState,
    VarTable,
    CallMeta,
    CFG,
    BasicBlock,
    slot_id,
    widenconst,
    ⊑,
    is_throw_call,
    tmerge,
    switchtupleunion,
    argtypes_to_type,
    argtype_by_index,
    argtype_tail,
    _methods_by_ftype,
    specialize_method,
    add_backedge!,
    add_mt_backedge!,
    compute_basic_blocks,
    matching_cache_argtypes,
    is_argtype_match,
    tuple_tfunc,
    may_invoke_generator,
    inlining_enabled

import Base:
    parse_input_line,
    unwrap_unionall,
    rewrap_unionall,
    uniontypes,
    IdSet

import Base.Meta:
    _parse_string,
    lower

using LoweredCodeUtils, JuliaInterpreter

import LoweredCodeUtils:
    istypedef,
    ismethod,
    callee_matches,
    rng,
    print_with_code,
    pushall!,
    # NamedVar,
    # add_requests!,
    add_ssa_preds!,
    # add_named_dependencies!,
    find_typedefs,
    add_control_flow!,
    add_typedefs!

import JuliaInterpreter:
    _INACTIVE_EXCEPTION,
    bypass_builtins,
    maybe_evaluate_builtin,
    collect_args,
    finish!,
    is_return,
    is_quotenode_egal,
    moduleof,
    @lookup

import MacroTools: MacroTools, @capture, normalise, striplines

using InteractiveUtils

using Pkg, Pkg.TOML

# common
# ======

const JET_DEV_MODE = parse(Bool, get(ENV, "JET_DEV_MODE", "false"))

const CONFIG_FILE_NAME = ".JET.toml"

const State     = Union{InferenceState,OptimizationState}
const StateAtPC = Tuple{State,Int}

# hooks
# -----

const INIT_HOOKS = Function[]
push_inithook!(f) = push!(INIT_HOOKS, f)
__init__() = foreach(@nospecialize(f)->f(), INIT_HOOKS)

# compat
# ------

const BACKEDGE_CALLBACK_ENABLED = :callbacks in fieldnames(Core.MethodInstance)

@static BACKEDGE_CALLBACK_ENABLED || push_inithook!() do
@warn """with your Julia version, JET.jl may not be able to update analysis result
correctly after refinement of a method in deeper call sites
"""
end

@static if isdefined(CC, :LimitedAccuracy)
    import .CC: ignorelimited
else
    ignorelimited(@nospecialize(x)) = x
end

@static if isdefined(Base, Symbol("@aggressive_constprop"))
    import Base: @aggressive_constprop
else
    macro aggressive_constprop(x) esc(x) end # not available
end

# early take in https://github.com/JuliaLang/julia/pull/41040
function gen_call_with_extracted_types_and_kwargs(__module__, fcn, ex0)
    kws = Expr[]
    arg = ex0[end] # Mandatory argument
    for i in 1:length(ex0)-1
        x = ex0[i]
        if x isa Expr && x.head === :(=) # Keyword given of the form "foo=bar"
            if length(x.args) != 2
                return Expr(:call, :error, "Invalid keyword argument: $x")
            end
            push!(kws, Expr(:kw, esc(x.args[1]), esc(x.args[2])))
        else
            return Expr(:call, :error, "@$fcn expects only one non-keyword argument")
        end
    end
    return InteractiveUtils.gen_call_with_extracted_types(__module__, fcn, arg, kws)
end

# macros
# ------

# to circumvent https://github.com/JuliaLang/julia/issues/37342, we inline these `isa`
# condition checks at surface AST level
macro isexpr(ex, args...)
    ex   = esc(ex)
    args = map(esc, args)
    return :($(GlobalRef(Core, :isa))($(ex), $(GlobalRef(Core, :Expr))) && $(_isexpr_check)($(ex), $(args...)))
end

_isexpr_check(ex::Expr, head::Symbol)         = ex.head === head
_isexpr_check(ex::Expr, heads)                = in(ex.head, heads)
_isexpr_check(ex::Expr, head::Symbol, n::Int) = ex.head === head && length(ex.args) == n
_isexpr_check(ex::Expr, heads, n::Int)        = in(ex.head, heads) && length(ex.args) == n

islnn(@nospecialize(_)) = false
islnn(::LineNumberNode) = true

"""
    @invoke f(arg::T, ...; kwargs...)

Provides a convenient way to call [`invoke`](https://docs.julialang.org/en/v1/base/base/#Core.invoke);
`@invoke f(arg1::T1, arg2::T2; kwargs...)` will be expanded into `invoke(f, Tuple{T1,T2}, arg1, arg2; kwargs...)`.
When an argument's type annotation is omitted, it's specified as `Any` argument, e.g.
`@invoke f(arg1::T, arg2)` will be expanded into `invoke(f, Tuple{T,Any}, arg1, arg2)`.

This could be used to call down to `NativeInterpreter`'s abstract interpretation method of
  `f` while passing `AbstractAnalyzer` so that subsequent calls of abstract interpretation
  functions overloaded against `AbstractAnalyzer` can be called from the native method of `f`;
e.g. calls down to `NativeInterpreter`'s `abstract_call_gf_by_type` method:
```julia
@invoke abstract_call_gf_by_type(analyzer::AbstractInterpreter, f, argtypes::Vector{Any}, atype, sv::InferenceState,
                                 max_methods::Int)
```
"""
macro invoke(ex)
    f, args, kwargs = destructure_callex(ex)
    arg2typs = map(args) do x
        if @isexpr(x, :macrocall) && first(x.args) === Symbol("@nospecialize")
            x = last(x.args)
        end
        @isexpr(x, :(::)) ? (x.args...,) : (x, GlobalRef(Core, :Any))
    end
    args, argtypes = first.(arg2typs), last.(arg2typs)
    return esc(:($(GlobalRef(Core, :invoke))($(f), Tuple{$(argtypes...)}, $(args...); $(kwargs...))))
end

"""
    @invokelatest f(args...; kwargs...)

Provides a convenient way to call [`Base.invokelatest`](https://docs.julialang.org/en/v1/base/base/#Base.invokelatest).
`@invokelatest f(args...; kwargs...)` will simply be expanded into
`Base.invokelatest(f, args...; kwargs...)`.
"""
macro invokelatest(ex)
    f, args, kwargs = destructure_callex(ex)
    return esc(:($(GlobalRef(Base, :invokelatest))($(f), $(args...); $(kwargs...))))
end

function destructure_callex(ex)
    @assert @isexpr(ex, :call) "call expression f(args...; kwargs...) should be given"

    f = first(ex.args)
    args = []
    kwargs = []
    for x in ex.args[2:end]
        if @isexpr(x, :parameters)
            append!(kwargs, x.args)
        elseif @isexpr(x, :kw)
            push!(kwargs, x)
        else
            push!(args, x)
        end
    end

    return f, args, kwargs
end

"""
    @withmixedhash (mutable) struct T
        fields ...
    end

Defines struct `T` while automatically defining its `Base.hash(::T, ::UInt)` method which
  mixes hashes of all of `T`'s fields (and also corresponding `Base.:(==)(::T, ::T)` method).

This macro is supposed to abstract the following kind of pattern:

> https://github.com/aviatesk/julia/blob/999973df2850d6b2e0bd4bcf03ef90a14217b63c/base/pkgid.jl#L3-L25
```julia
struct PkgId
    uuid::Union{UUID,Nothing}
    name::String
end

==(a::PkgId, b::PkgId) = a.uuid == b.uuid && a.name == b.name

function hash(pkg::PkgId, h::UInt)
    h += 0xc9f248583a0ca36c % UInt
    h = hash(pkg.uuid, h)
    h = hash(pkg.name, h)
    return h
end
```

> with `@withmixedhash`
```julia
@withmixedhash struct PkgId
    uuid::Union{UUID,Nothing}
    name::String
end
```

See also: [`EGAL_TYPES`](@ref)
"""
macro withmixedhash(typedef)
    @assert @isexpr(typedef, :struct) "struct definition should be given"
    name = esc(typedef.args[2])
    fld2typs = map(filter(!islnn, typedef.args[3].args)) do x
        if @isexpr(x, :(::))
            fld, typex = x.args
            typ = Core.eval(__module__, typex)
            fld, typ
        else
            (x, Any)
        end
    end
    @assert !isempty(fld2typs) "no fields given, nothing to hash"

    h_init = UInt === UInt64 ? rand(UInt64) : rand(UInt32)
    hash_body = quote h = $(h_init) end
    for (fld, typ) in fld2typs
        push!(hash_body.args, :(h = $(Base.hash)(x.$fld, h)))
    end
    push!(hash_body.args, :(return h))
    hash_func = :(function Base.hash(x::$name, h::$UInt); $(hash_body); end)
    eq_body = foldr(fld2typs; init = true) do (fld, typ), x
        eq_ex = if all(in(EGAL_TYPES), typenames(typ))
            :(x1.$fld === x2.$fld)
        else
            :(x1.$fld == x2.$fld)
        end
        Expr(:&&, eq_ex, x)
    end
    eq_func = :(function Base.:(==)(x1::$name, x2::$name); $(eq_body); end)

    return quote
        Base.@__doc__ $(typedef)
        $(hash_func)
        $(eq_func)
    end
end

"""
    const EGAL_TYPES = $(EGAL_TYPES)

Keeps names of types that should be compared by `===` rather than `==`.

See also: [`@withmixedhash`](@ref)
"""
const EGAL_TYPES = Set{Symbol}((:Type, :Symbol, :MethodInstance))

typenames(@nospecialize(t::DataType)) = [t.name.name]
typenames(u::Union)                   = collect(Base.Iterators.flatten(typenames.(CC.uniontypes(u))))
typenames(u::UnionAll)                = [u.body.name.name]

"""
    @jetconfigurable function config_func(args...; configurations...)
        ...
    end

This macro asserts that there's no configuration naming conflict across the `@jetconfigurable`
  functions so that a configuration for a `@jetconfigurable` function  doesn't affect the other
  `@jetconfigurable` functions.
This macro also adds a dummy splat keyword arguments (`jetconfigs...`) to the function definition
  so that any configuration of other `@jetconfigurable` functions can be passed on to it.
"""
macro jetconfigurable(funcdef)
    @assert @isexpr(funcdef, :(=)) || @isexpr(funcdef, :function) "function definition should be given"

    defsig = funcdef.args[1]
    if @isexpr(defsig, :where)
        defsig = first(defsig.args)
    end
    args = defsig.args::Vector{Any}
    thisname = first(args)
    i = findfirst(a->@isexpr(a, :parameters), args)
    if isnothing(i)
        @warn "no JET configurations are defined for `$thisname`"
        insert!(args, 2, Expr(:parameters, :(jetconfigs...)))
    else
        kwargs = args[i]
        found = false
        for kwarg in kwargs.args
            if @isexpr(kwarg, :...)
                found = true
                continue
            end
            kwargex = first(kwarg.args)
            kwargname = (@isexpr(kwargex, :(::)) ? first(kwargex.args) : kwargex)::Symbol
            othername = get!(_JET_CONFIGURATIONS, kwargname, thisname)
            # allows same configurations for same generic function or function refinement
            @assert thisname == othername "`$thisname` uses `$kwargname` JET configuration name which is already used by `$othername`"
        end
        found || push!(kwargs.args, :(jetconfigs...))
    end

    return esc(funcdef)
end
const _JET_CONFIGURATIONS = Dict{Symbol,Union{Symbol,Expr}}()

# utils
# -----

# common

@inline get_stmt((sv, pc)::StateAtPC)         = @inbounds sv.src.code[pc]
@inline get_lin((sv, pc)::StateAtPC)          = @inbounds (sv.src.linetable::Vector)[sv.src.codelocs[pc]]::LineInfoNode
@inline get_ssavaluetype((sv, pc)::StateAtPC) = @inbounds sv.src.ssavaluetypes[pc]

@inline get_slottype(s::Union{StateAtPC,State}, slot)      = get_slottype(s, slot_id(slot))
@inline get_slottype((sv, pc)::StateAtPC,       slot::Int) = get_slottype(sv, slot)
@inline get_slottype(sv::State,                 slot::Int) = @inbounds sv.slottypes[slot]

@inline get_slotname(s::Union{StateAtPC,State}, slot)      = get_slotname(s, slot_id(slot))
@inline get_slotname((sv, pc)::StateAtPC,       slot::Int) = @inbounds sv.src.slotnames[slot]
@inline get_slotname(sv::State,                 slot::Int) = @inbounds sv.src.slotnames[slot]

# InfernceState

# we can retrieve program-counter-level slottype during inference
@inline get_slottype(s::Tuple{InferenceState,Int}, slot::Int) = @inbounds (get_states(s)[slot]::VarState).typ
@inline get_states((sv, pc)::Tuple{InferenceState,Int})       = @inbounds sv.stmt_types[pc]::VarTable

@inline get_currpc(sv::InferenceState)  = min(sv.currpc, length(sv.src.code))
@inline get_result(sv::InferenceState)  = sv.bestguess

function is_constant_propagated(frame::InferenceState)
    return !frame.cached && CC.any(frame.result.overridden_by_const)
end

prewalk_inf_frame(@nospecialize(f), ::Nothing) = return
function prewalk_inf_frame(@nospecialize(f), frame::InferenceState)
    ret = f(frame)
    prewalk_inf_frame(f, frame.parent)
    return ret
end

postwalk_inf_frame(@nospecialize(f), ::Nothing) = return
function postwalk_inf_frame(@nospecialize(f), frame::InferenceState)
    postwalk_inf_frame(f, frame.parent)
    return f(frame)
end

# lattice

ignorenotfound(@nospecialize(t)) = t === NOT_FOUND ? Bottom : t

# includes
# ========

# interface
include("interfaces.jl")
# abstract interpretaion based analysis
include("tfuncs.jl")
include("abstractinterpretation.jl")
include("typeinfer.jl")
include("analyzer.jl")
include("jetcache.jl")
include("locinfo.jl")
# top-level analysis
include("graph.jl")
include("virtualprocess.jl")
# watch mode
include("watch.jl")

# results
# =======

"""
    res::JETToplevelResult

Represents the result of JET's analysis on a top-level script.
- `res.analyzer::AbstractAnalyzer`: [`AbstractAnalyzer`](@ref) used for this analysis
- `res.res::VirtualProcessResult`: [`VirtualProcessResult`](@ref) collected from this analysis
- `res.source::String`: the identity key of this analysis
- `res.jetconfigs`: [JET configurations](@ref JET-configurations) used for this analysis

`JETToplevelResult` implements `show` methods for each different frontend.
An appropriate `show` method will be automatically choosen and render the analysis result.
"""
struct JETToplevelResult{Analyzer<:AbstractAnalyzer,JETConfigs}
    analyzer::Analyzer
    res::VirtualProcessResult
    source::String
    jetconfigs::JETConfigs
end
JETToplevelResult(analyzer, res, source; jetconfigs...) =
    return JETToplevelResult(analyzer, res, source, jetconfigs)
@eval Base.iterate(res::JETToplevelResult, state=1) =
    return state > $(fieldcount(JETToplevelResult)) ? nothing : (getfield(res, state), state+1)

"""
    res::JETCallResult

Represents the result of JET's analysis on a function call.
- `res.analyzer::AbstractAnalyzer`: [`AbstractAnalyzer`](@ref) used for this analysis
- `res.type`: the return type of the function call
- `res.source::String`: the identity key of this analysis
- `res.jetconfigs`: [JET configurations](@ref JET-configurations) used for this analysis

`JETCallResult` implements `show` methods for each different frontend.
An appropriate `show` method will be automatically choosen and render the analysis result.
"""
struct JETCallResult{Analyzer<:AbstractAnalyzer,JETConfigs}
    analyzer::Analyzer
    type
    source::String
    jetconfigs::JETConfigs
    JETCallResult(analyzer::Analyzer, @nospecialize(type), source;
                  jetconfigs...) where {Analyzer<:AbstractAnalyzer} =
        new{Analyzer,typeof(jetconfigs)}(analyzer, type, source, jetconfigs)
end
@eval Base.iterate(res::JETCallResult, state=1) =
    return state > $(fieldcount(JETCallResult)) ? nothing : (getfield(res, state), state+1)

function get_critical_reports(res::VirtualProcessResult)
    # non-empty `ret.toplevel_error_reports` means critical errors happened during
    # the AST transformation, so they always have precedence over `ret.inference_error_reports`
    return !isempty(res.toplevel_error_reports) ? res.toplevel_error_reports :
           res.inference_error_reports
end

# when virtualized, fix virtual module printing based on string manipulation;
# the "actual" modules may not be loaded into this process
gen_postprocess(::Nothing) = return identity
function gen_postprocess((actualmod, virtualmod)::Actual2Virtual)
    virtual = string(virtualmod)
    actual  = string(actualmod)
    return actualmod === Main ?
           replace2("Main." => "") ∘ replace2(virtual => actual) :
           replace2(virtual => actual)
end
replace2(pat) = x -> replace(x, pat)

# we may need something like this for stdlibs as well ?

function tofullpath(filename::AbstractString)
    path = abspath(filename)
    return isfile(path) ? path : fullbasepath(filename)
end
fullbasepath(filename) = normpath(JULIA_DIR, "base", filename)

# TODO make this configurable ?
const JULIA_DIR = @static occursin("DEV", string(VERSION)) ?
                  normpath(Sys.BINDIR, "..", "..", "..", "julia") :
                  normpath(Sys.BINDIR, Base.DATAROOTDIR)

# UIs
# ===

# default UI (console)
include("print.jl")

# entry
# =====

"""
    report_file(filename::AbstractString;
                jetconfigs...) -> JETToplevelResult

Analyzes `filename` and returns [`JETToplevelResult`](@ref).

This function will look for `$CONFIG_FILE_NAME` configuration file in the directory of `filename`,
  and search _up_ the file tree until any `$CONFIG_FILE_NAME` is (or isn't) found.
When found, the configurations specified in the file will be applied.
See [Configuration File](@ref) for more details.

!!! tip
    When you want to analyze your package, but any file actually using it isn't available, the
      `analyze_from_definitions` option can be useful (see [`ToplevelConfig`](@ref)'s `analyze_from_definitions` option). \\
    For example, JET can analyze JET itself like below:
    ```julia
    # from the root directory of JET.jl
    julia> report_file("src/JET";
                       analyze_from_definitions = true)
    ```

    See also: [`report_package`](@ref)

!!! note
    This function will enable the `toplevel_logger` configuration by default with the default logging level.
    You can still explicitly specify and configure it:
    ```julia
    report_file(args...;
                toplevel_logger = nothing, # suppress toplevel logger
                jetconfigs...) # other configurations
    ```
    See [Logging Configurations](@ref) for more details.
"""
function report_file(filename::AbstractString;
                     source::Union{Nothing,AbstractString} = nothing,
                     jetconfigs...) where T
    isfile(filename) || throw(ArgumentError("$filename doesn't exist"))

    configfile = find_config_file(dirname(abspath(filename)))
    if isnothing(configfile)
        jetconfigs = set_toplevel_logger_if_missing(jetconfigs)
    else
        config = parse_config_file(configfile)
        jetconfigs = overwrite_options(config, jetconfigs)
        jetconfigs = set_toplevel_logger_if_missing(jetconfigs)
        with_logger(get(jetconfigs, :toplevel_logger, nothing), ≥(INFO_LOGGER_LEVEL), "toplevel") do io
            println(io, "applied JET configurations in $configfile")
        end
    end

    if isnothing(source)
        source = string(nameof(var"#self#"), "(\"$filename\")")
    end

    return report_text(read(filename, String), filename; source, jetconfigs...)
end

function set_toplevel_logger_if_missing(@nospecialize(jetconfigs))
    haskey(jetconfigs, :toplevel_logger) && return jetconfigs
    default_toplevel_logger_config =
        kwargs((; toplevel_logger = IOContext(stdout::IO, LOGGER_LEVEL_KEY => DEFAULT_LOGGER_LEVEL)))
    return overwrite_options(jetconfigs, default_toplevel_logger_config)
end

function find_config_file(dir)
    next_dir = dirname(dir)
    if (next_dir == dir || # ensure to escape infinite recursion
        isempty(dir))      # reached to the system root
        return nothing
    end
    path = normpath(dir, CONFIG_FILE_NAME)
    return isfile(path) ? path : find_config_file(next_dir)
end

"""
JET.jl offers [`.prettierrc` style](https://prettier.io/docs/en/configuration.html)
  configuration file support.
This means you can use `$CONFIG_FILE_NAME` configuration file to specify any of configurations
  explained above and share that with others.

When [`$report_file`](@ref) or [`$report_and_watch_file`](@ref) is called, it will look for
  `$CONFIG_FILE_NAME` in the directory of the given file, and search _up_ the file tree until
  a JET configuration file is (or isn't) found.
When found, the configurations specified in the file will be applied.

A configuration file can specify any of JET configurations like:
```toml
aggressive_constant_propagation = false # turn off aggressive constant propagation
... # other configurations
```

Note that the following configurations should be string(s) of valid Julia code:
- `context`: string of Julia code, which can be `parse`d and `eval`uated into `Module`
- `concretization_patterns`: vector of string of Julia code, which can be `parse`d into a
  Julia expression pattern expected by [`MacroTools.@capture` macro](https://fluxml.ai/MacroTools.jl/stable/pattern-matching/).
- `toplevel_logger`: string of Julia code, which can be `parse`d and `eval`uated into `Union{IO,Nothing}`
- `inference_logger`: string of Julia code, which can be `parse`d and `eval`uated into `Union{IO,Nothing}`

E.g. the configurations below are equivalent:
- configurations via keyword arguments
  ```julia
  report_file(somefile;
              concretization_patterns = [:(const GLOBAL_CODE_STORE = x_)],
              toplevel_logger = IOContext(open("toplevel.txt", "w"), :JET_LOGGER_LEVEL => 1))
  ```
- configurations via a configuration file
  $(let
      text = read(normpath(@__DIR__, "..", "test", "fixtures", "..JET.toml"), String)
      lines = split(text, '\n')
      pushfirst!(lines, "```toml"); push!(lines, "```")
      join(lines, "\n  ")
  end)

!!! note
    JET configurations specified as keyword arguments have precedence over those specified
      via a configuration file.
"""
parse_config_file(path) = process_config_dict!(TOML.parsefile(path))

function process_config_dict!(config_dict)
    context = get(config_dict, "context", nothing)
    if !isnothing(context)
        @assert isa(context, String) "`context` should be string of Julia code"
            config_dict["context"] = Core.eval(Main, trymetaparse(context))
    end
    concretization_patterns = get(config_dict, "concretization_patterns", nothing)
    if !isnothing(concretization_patterns)
        @assert isa(concretization_patterns, Vector{String}) "`concretization_patterns` should be array of string of Julia expression"
        config_dict["concretization_patterns"] = trymetaparse.(concretization_patterns)
    end
    toplevel_logger = get(config_dict, "toplevel_logger", nothing)
    if !isnothing(toplevel_logger)
        @assert isa(toplevel_logger, String) "`toplevel_logger` should be string of Julia code"
        config_dict["toplevel_logger"] = Core.eval(Main, trymetaparse(toplevel_logger))
    end
    inference_logger = get(config_dict, "inference_logger", nothing)
    if !isnothing(inference_logger)
        @assert isa(inference_logger, String) "`inference_logger` should be string of Julia code"
        config_dict["inference_logger"] = Core.eval(Main, trymetaparse(inference_logger))
    end
    return kwargs(config_dict)
end

function trymetaparse(s)
    ret = Meta.parse(strip(s); raise = true)
    @isexpr(ret, :incomplete) && error(first(ret.args))
    return ret
end

function kwargs(dict)
    ns = (Symbol.(keys(dict))...,)
    vs = (collect(values(dict))...,)
    return pairs(NamedTuple{ns}(vs))
end

overwrite_options(old, new) = kwargs(merge(old, new))

"""
    report_package(package::Union{AbstractString,Module};
                   jetconfigs...) -> JETToplevelResult

Analyzes `package` in the same way as [`report_file`](@ref) with the special default
configurations, which are especially tuned for package analysis (see below for details).
`package` can be either a module or a string. In the latter case it must be the name of a
package in your current environment.

This function configures analysis with the following configurations:
- `analyze_from_definitions = true`: allows JET to enter analysis without top-level call sites;
  this is useful for package analysis since a package itself usually has only definitions
  but not usages (i.e. call sites)
- `concretization_patterns = [:(x_ = y_), :(const x_ = y_)]`: concretizes every global variable instantiations;
  concretizations are generally preferred for successful analysis as far as they're cheap,
  and global variable instantiations that occur in a package definition are _usually_ very cheap
See [`ToplevelConfig`](@ref) for more details.

---

    report_package([io::IO = stdout];
                   jetconfigs...) -> res::ReportResult

Like above but analyzes the package of the current project.

See also: [`report_file`](@ref)
"""
function report_package(package::Union{AbstractString,Module,Nothing} = nothing;
                        analyze_from_definitions::Bool = true,
                        concretization_patterns = [:(x_ = y_), :(const x_ = y_)],
                        jetconfigs...)
    filename = get_package_file(package)
    return report_file(filename; analyze_from_definitions, jetconfigs...)
end

function get_package_file(package::AbstractString)
    filename = Base.find_package(package)
    isnothing(filename) && throw(ErrorException("unknown package $package."))
    return filename
end

function get_package_file(package::Module)
    filename = pathof(package)
    isnothing(filename) && throw(ErrorException("cannot analyze a module defined in the REPL."))
    return filename
end

function get_package_file(::Nothing)
    project = Pkg.project()
    project.ispackage || throw(ErrorException("active project at $(project.path) is not a package."))
    return normpath(dirname(project.path), "src", project.name::String * ".jl")
end

"""
    report_text(text::AbstractString,
                filename::AbstractString = "top-level";
                jetconfigs...) -> JETToplevelResult

Analyzes `text` and returns [`JETToplevelResult`](@ref).
"""
function report_text(text::AbstractString,
                     filename::AbstractString = "top-level";
                     analyzer::Type{Analyzer} = JETAnalyzer,
                     source::Union{Nothing,AbstractString} = nothing,
                     jetconfigs...) where {Analyzer<:AbstractAnalyzer}
    analyzer′ = Analyzer(; jetconfigs...)
    maybe_initialize_caches!(analyzer′)
    config = ToplevelConfig(; jetconfigs...)
    res = virtual_process(text,
                          filename,
                          analyzer′,
                          config,
                          )
    if isnothing(source)
        source = string(nameof(var"#self#"), "(..., \"$filename\")")
    end
    return JETToplevelResult(analyzer′, res, source; analyzer, jetconfigs...)
end

function analyze_toplevel!(analyzer::AbstractAnalyzer, src::CodeInfo)
    # construct toplevel `MethodInstance`
    mi = ccall(:jl_new_method_instance_uninit, Ref{Core.MethodInstance}, ());
    mi.uninferred = src
    mi.specTypes = Tuple{}

    transform_abstract_global_symbols!(analyzer, src)
    mi.def = get_toplevelmod(analyzer)

    result = InferenceResult(mi);
    # toplevel frame doesn't need to be cached (and so it won't be optimized), nor should
    # go through JET's code generation error check
    frame = InferenceState(result, src, #=cached=# false, analyzer);

    return analyze_frame!(analyzer, frame)
end

# HACK this is an native hack to re-use `AbstractInterpreter`'s approximated slot types for
# assignments of abstract global variables, which are represented as toplevel symbols at this point;
# the idea is just to transform them into slots from symbols and use their approximated type
# on their assignment.
# NOTE that `transform_abstract_global_symbols!` will produce really invalid code for
# actual interpretation or execution, but all the statements won't be interpreted anymore
# by `ConcreteInterpreter` nor executed anyway since toplevel frames aren't cached
function transform_abstract_global_symbols!(analyzer::AbstractAnalyzer, src::CodeInfo)
    nslots = length(src.slotnames)
    abstrct_global_variables = Dict{Symbol,Int}()
    concretized = get_concretized(analyzer)

    # linear scan, and find assignments of abstract global variables
    for (i, stmt) in enumerate(src.code::Vector{Any})
        if !(concretized[i])
            if @isexpr(stmt, :(=))
                lhs = first(stmt.args)
                if isa(lhs, Symbol)
                    if !haskey(abstrct_global_variables, lhs)
                        nslots += 1
                        push!(abstrct_global_variables, lhs => nslots)
                    end
                end
            end
        end
    end

    prewalk_and_transform!(src) do x, scope
        if isa(x, Symbol)
            slot = get(abstrct_global_variables, x, nothing)
            isnothing(slot) || return SlotNumber(slot)
        end
        return x
    end

    resize!(src.slotnames, nslots)
    resize!(src.slotflags, nslots)
    for (slotname, idx) in abstrct_global_variables
        src.slotnames[idx] = slotname
    end

    set_global_slots!(analyzer, Dict(idx => slotname for (slotname, idx) in abstrct_global_variables))
end

# TODO `analyze_builtin!` ?
function analyze_gf_by_type!(analyzer::AbstractAnalyzer, @nospecialize(tt::Type{<:Tuple}))
    mm = get_single_method_match(tt, InferenceParams(analyzer).MAX_METHODS, get_world_counter(analyzer))
    return analyze_method_signature!(analyzer, mm.method, mm.spec_types, mm.sparams)
end

function get_single_method_match(@nospecialize(tt), lim, world)
    mms = _methods_by_ftype(tt, lim, world)
    @assert !isa(mms, Bool) "unable to find matching method for $(tt)"

    filter!(mm::MethodMatch->mm.spec_types===tt, mms)
    @assert length(mms) == 1 "unable to find single target method for $(tt)"

    return first(mms)::MethodMatch
end

analyze_method!(analyzer::AbstractAnalyzer, m::Method) =
    analyze_method_signature!(analyzer, m, m.sig, method_sparams(m))

function method_sparams(m::Method)
    s = TypeVar[]
    sig = m.sig
    while isa(sig, UnionAll)
        push!(s, sig.var)
        sig = sig.body
    end
    return svec(s...)
end

analyze_method_signature!(analyzer::AbstractAnalyzer, m::Method, @nospecialize(atype), sparams::SimpleVector) =
    analyze_method_instance!(analyzer, specialize_method(m, atype, sparams))

function analyze_method_instance!(analyzer::AbstractAnalyzer, mi::MethodInstance)
    result = InferenceResult(mi)

    frame = InferenceState(result, #=cached=# true, analyzer)
    isnothing(frame) && return analyzer, nothing

    return analyze_frame!(analyzer, frame)
end

function InferenceState(result::InferenceResult, cached::Bool, analyzer::AbstractAnalyzer)
    ReportPass(analyzer)(GeneratorErrorReport, analyzer, result.linfo)
    return @invoke InferenceState(result::InferenceResult, cached::Bool, analyzer::AbstractInterpreter)
end

@reportdef struct GeneratorErrorReport <: InferenceErrorReport
    @nospecialize(err) # actual error wrapped
end
get_msg(::Type{GeneratorErrorReport}, analyzer::AbstractAnalyzer, linfo::MethodInstance, @nospecialize(err)) =
    return sprint(showerror, err)

# XXX what's the "soundness" of a `@generated` function ?
# adapated from https://github.com/JuliaLang/julia/blob/f806df603489cfca558f6284d52a38f523b81881/base/compiler/utilities.jl#L107-L137
function (::SoundBasicPass)(::Type{GeneratorErrorReport}, analyzer::AbstractAnalyzer, mi::MethodInstance)
    m = mi.def::Method
    if isdefined(m, :generator)
        # analyze_method_instance!(analyzer, linfo) XXX doesn't work
        may_invoke_generator(mi) || return
        try
            ccall(:jl_code_for_staged, Any, (Any,), mi)
        catch err
            # if user code throws error, wrap and report it
            add_new_report!(GeneratorErrorReport(analyzer, mi, err), analyzer)
        end
    end
end

function analyze_frame!(analyzer::AbstractAnalyzer, frame::InferenceState)
    typeinf(analyzer, frame)

    # report `throw` calls "appropriately";
    # if the final return type here is `Bottom`-annotated, it _may_ mean the control flow
    # didn't catch some of the `UncaughtExceptionReport`s stashed within `uncaught_exceptions(analyzer)`,
    if get_result(frame) === Bottom
        if !isempty(get_uncaught_exceptions(analyzer))
            append!(get_reports(analyzer), get_uncaught_exceptions(analyzer))
        end
    end

    return analyzer, frame
end

# test, interactive
# =================

# TODO improve inferrability by making `analyzer` argument positional

"""
    @report_call [jetconfigs...] f(args...)

Evaluates the arguments to the function call, determines its types, and then calls
  [`report_call`](@ref) on the resulting expression.
As with `@code_typed` and its family, any of [JET configurations](@ref) can be given as the optional
  arguments like this:
```julia
# reports `rand(::Type{Bool})` with `aggressive_constant_propagation` configuration turned off
julia> @report_call aggressive_constant_propagation=false rand(Bool)
```
"""
macro report_call(ex0...)
    return gen_call_with_extracted_types_and_kwargs(__module__, :report_call, ex0)
end

"""
    report_call(f, types = Tuple{};
                analyzer::Type{<:AbstractAnalyzer} = JETAnalyzer,
                jetconfigs...) -> JETCallResult
    report_call(tt::Type{<:Tuple};
                analyzer::Type{<:AbstractAnalyzer} = JETAnalyzer,
                jetconfigs...) -> JETCallResult

Analyzes the generic function call with the given type signature with `analyzer`.
And finally returns the analysis result as [`JETCallResult`](@ref).
"""
function report_call(@nospecialize(f), @nospecialize(types = Tuple{}); jetconfigs...)
    ft = Core.Typeof(f)
    if isa(types, Type)
        u = unwrap_unionall(types)
        tt = rewrap_unionall(Tuple{ft, u.parameters...}, types)
    else
        tt = Tuple{ft, types...}
    end
    return report_call(tt; jetconfigs...)
end

function report_call(@nospecialize(tt::Type{<:Tuple});
                     analyzer::Type{Analyzer} = JETAnalyzer,
                     source::Union{Nothing,AbstractString} = nothing,
                     jetconfigs...) where {Analyzer<:AbstractAnalyzer}
    analyzer = Analyzer(; jetconfigs...)
    maybe_initialize_caches!(analyzer)
    analyzer, frame = analyze_gf_by_type!(analyzer, tt)

    if isnothing(frame)
        # if there is `GeneratorErrorReport`, it means the code generation happened and failed
        rt = any(r->isa(r, GeneratorErrorReport), get_reports(analyzer)) ? Bottom : Any
    else
        rt = get_result(frame)
    end

    if isnothing(source)
        source = string(nameof(var"@report_call"), " ", sprint(show_tuple_as_call, Symbol(""), tt))
    end

    return JETCallResult(analyzer, rt, source; jetconfigs...)
end

# for inspection
macro lwr(ex) QuoteNode(lower(__module__, ex)) end
macro src(ex) QuoteNode(first(lower(__module__, ex).args)) end

# exports
# =======

"""
    JETInterfaces

This `baremodule` exports names that form the APIs of [JET.jl Pluggable Analysis Framework](@ref).
If you are going to define a plug-in analysis, `using JET.JETInterfaces` will load most useful
names described below.
"""
baremodule JETInterfaces
const DOCUMENTED_NAMES = Symbol[] # will be used in docs/make.jl
end

function reexport_as_api!(xs...; documented = true)
    for x in xs
        ex = Expr(:block)

        val = Core.eval(@__MODULE__, x)
        canonicalname = Symbol(parentmodule(val), '.', nameof(val))
        canonicalpath = Symbol.(split(string(canonicalname), '.'))

        modpath = Expr(:., canonicalpath[1:end-1]...)
        symname = last(canonicalpath)
        sympath = Expr(:., symname)
        importex = Expr(:import, Expr(:(:), modpath, sympath))
        exportex = Expr(:export, symname)

        push!(ex.args, importex, exportex)
        Core.eval(JETInterfaces, ex)
        documented && push!(JETInterfaces.DOCUMENTED_NAMES, symname)
    end
end

reexport_as_api!(AbstractAnalyzer,
                 AnalyzerState,
                 ReportPass,
                 InferenceErrorReport,
                 get_msg,
                 get_spec_args,
                 var"@reportdef",
                 add_new_report!,
                 )
reexport_as_api!(subtypes(InferenceErrorReport)...; documented = false)
reexport_as_api!(subtypes(ReportPass)...; documented = false)

# export entries
export
    report_file,
    report_and_watch_file,
    report_package,
    report_text,
    @report_call,
    report_call

end


text:
@enum TimeSystemType begin
    # add new systems BEFORE count, then
    # *** add to Strings[] in TimeSystem.cpp and make parallel to this enum. ***

    # Unknown MUST BE FIRST, and must = 0
    tstUnknown       # /< unknown time frame; for legacy code compatibility
    tstAny           # /< wildcard; allows comparison with any other type
    tstGPS           # /< GPS system time
    tstGLO           # /< GLONASS system time
    tstGAL           # /< Galileo system time
    tstQZS           # /< QZSS system Time
    tstBDT           # /< BeiDou system Time
    tstIRN           # /< IRNSS system Time
    tstUTC           # /< Coordinated Universal Time (e.g., from NTP)
    tstTAI           # /< International Atomic Time
    tstTT            # /< Terrestrial time (used in IERS conventions)
    tstTDB           # /< Barycentric dynamical time (JPL ephemeris); very near TT
        # count MUST BE LAST
    tstcount         # /< the number of systems - not a system
end

# struct GNSSDatTime <: GNSSDatTime
#     dt::DateTime
#     tst::TimeSystemType
# end

const stt2tst = Dict{SatelliteSystemType,TimeSystemType}(
    sstGPS => tstGPS,           # /< GPS system time
    sstGeosync => tstGPS,           # /< GPS system time
    sstGlonass => tstGLO,           # /< GLONASS system time
    sstGalileo => tstGAL,           # /< Galileo system time
    sstQZSS => tstGPS,           # /< QZSS system Time
    sstBeiDou => tstBDT,           # /< BeiDou system Time
    sstIRNSS => tstIRN,           # /< IRNSS system Time
    # "UT" => tstUTC,           # /< Coordinated Universal Time (e.g., from NTP)
    # "TA" => tstTAI,           # /< International Atomic Time

)
const str2tst = Dict{String,TimeSystemType}(
    "GP" => tstGPS,           # /< GPS system time
    "GL" => tstGLO,           # /< GLONASS system time
    "GA" => tstGAL,           # /< Galileo system time
    "QZ" => tstQZS,           # /< QZSS system Time
    "BD" => tstBDT,           # /< BeiDou system Time
    "IR" => tstIRN,           # /< IRNSS system Time
    "UT" => tstUTC,           # /< Coordinated Universal Time (e.g., from NTP)
    "TA" => tstTAI,           # /< International Atomic Time
)

function str2TimeSystemType(str::String,
        str2tst::Dict{String,TimeSystemType}=str2tst,
)::TimeSystemType
    haskey(str2tst, str) ||
        throw(error("no such TimeSystemType: $str"))
    str2tst[str]
end

# Time system conversions constants
const TAI_minus_GPSGAL_EPOCH = 19.0
const TAI_minus_BDT_EPOCH = 33.0
# const TAI_minus_TT_EPOCH = -32.184

function diff2UTC(indt::DateTime, inTS::TimeSystemType)
    # -----------------------------------------------------------
    # conversions: first convert inTS->TAI ...
    # TAI = GPS + 19s
    # TAI = UTC + getLeapSeconds()
    # TAI = TT - 32.184s
    if inTS == tstGPS ||       # GPS -> TAI
        inTS == tstGAL ||       # GAL -> TAI
        inTS == tstIRN         # IRN -> TAI
        dt = TAI_minus_GPSGAL_EPOCH
    elseif inTS == tstUTC ||  # UTC -> TAI
        inTS == tstGLO    # GLO -> TAI
        dt = -offset_utc_tai(indt)
    elseif inTS == tstBDT    # BDT -> TAI
        dt = TAI_minus_BDT_EPOCH;
    elseif inTS == tstTAI    # TAI
        dt  = 0.0
    # elseif inTS == TT     # TT -> TAI
    #     dt = TAI_minus_TT_EPOCH;
    # elseif inTS == tstTDB    # TDB -> TAI
    #     dt = TAI_minus_TT_EPOCH + TDBmTT
    else                               # other
        throw(error("no such TimeSystem: $inTS"))
    end

    # -----------------------------------------------------------
    # ... then convert TAI->outTS
    # GPS = TAI - 19s
    # UTC = TAI - getLeapSeconds()
    # TT = TAI + 32.184s
    # if outTS == tstGPS ||      # TAI -> GPS
    # outTS == tstGAL ||      # TAI -> GAL
    # outTS == tstIRN         # TAI -> IRN
    #     dt -= TAI_minus_GPSGAL_EPOCH
    # elseif outTS == tstUTC || # TAI -> UTC
    #     outTS == tstGLO   # TAI -> GLO
    dt += offset_utc_tai(indt)
    # elseif outTS == tstBDT   # TAI -> BDT
    #     dt -= TAI_minus_BDT_EPOCH
    # elseif outTS == tstTAI   # TAI
    #     dt += 0.0
    # elseif outTS == TT    # TAI -> TT
    #     dt -= TAI_minus_TT_EPOCH;
    # elseif outTS == tstTDB   # TAI -> TDB
    #     dt -= TAI_minus_TT_EPOCH + TDBmTT
    # else                               # other
    #     throw(error("no such TimeSystem: $inTS"))
    # end

    return Millisecond(dt * 1000.0)
end


text:
# This example is described in more detail in the documentation of Trixi

using Trixi, LinearAlgebra, ForwardDiff

equations = CompressibleEulerEquations2D(1.4)

mesh = TreeMesh((-1.0, -1.0), (1.0, 1.0), initial_refinement_level=2, n_cells_max=10^5)

solver = DGSEM(polydeg=3, flux_lax_friedrichs, VolumeIntegralFluxDifferencing(flux_ranocha))

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_isentropic_vortex, solver)

u0_ode = compute_coefficients(0.0, semi)

J = ForwardDiff.jacobian((du_ode, γ) -> begin
      equations_inner = CompressibleEulerEquations2D(first(γ))
      semi_inner = Trixi.remake(semi, equations=equations_inner, uEltype=eltype(γ));
      Trixi.rhs!(du_ode, u0_ode, semi_inner, 0.0)
    end, similar(u0_ode), [1.4]); # γ needs to be an `AbstractArray`


text:
"""
$(TYPEDEF)
$(TYPEDFIELDS)

Contain all places, transitions & arcs. Pages are for visual presentation.
There must be at least 1 Page for a valid pnml model.
"""
struct Page{PNTD<:PnmlType,D} <: PnmlObject
    id::Symbol
    places::Vector{Place}
    refPlaces::Vector{RefPlace}
    transitions::Vector{Transition}
    refTransitions::Vector{RefTransition}
    arcs::Vector{Arc}
    declaration::D
    subpages::Maybe{Vector{Page}}
    name::Maybe{Name}
    com::ObjectCommon
end

function Page(pntd::PNTD, id::Symbol, places, refp, transitions, reft,
                arcs, declare, pages, name, oc::ObjectCommon) where {PNTD<:PnmlType}
    Page{typeof(pntd),
         typeof(declare)}(id, places, refp, transitions, reft, arcs, declare, pages, name, oc)
end

# Note declaration wraps a vector of AbstractDeclarations.
declarations(page::Page) = declarations(page.declaration)
pages(page::Page) = page.subpages

function Base.empty!(page::Page)
    empty!(page.places)
    empty!(page.refPlaces)
    empty!(page.transitions)
    empty!(page.refTransitions)
    empty!(page.arcs)
    empty!(page.declaration)
    !isnothing(page.subpages) && empty!(page.subpages)
    has_tools(page.com) && empty!(page.com.tools)
    has_labels(page.com) && empty!(page.com.labels)
end