text:
using NeidSolarScripts
using ArgParse

function parse_commandline_make_pyrheliometer_daily()
     s = ArgParseSettings( description = "Make pyrheliometer.csv for specified day.")
     @add_arg_table! s begin
         "manifest_or_date"
            help = "CSV Filename with fields l?filename, object and exptime to get pyrheliometer data for OR Date in YYYYMMDD format to extract pyrheliometer data."
            arg_type = String
            #default = "meta_test.csv" # joinpath(pwd(),"pyrohelio.csv")
            #default = "20210505"
            #default = "20210605"
         "--output"
            help = "Path for outputs: pyrheliometer.csv"
            arg_type = String
            default = joinpath(pwd(),"pyrohelio.csv")
         "--pyrheliometer_dir"
            help = "Path for pyroheliometer inputs: *.tel"
            arg_type = String
            #default = "/gpfs/group/ebf11/default/ebf11/neid_solar/data/pyrohelio/202110_fromChad/"
            default = "/mnt/data_simons/NEID/pyrohelio/"
         "--user"
            help = "NExScI Archive username"
            arg_type = String
         "--password"
            help = "NExScI Archive password"
            arg_type = String
         "--nexsci_login_filename"
            help = "TOML file with NExScI username and password"
            arg_type = String
            default = "nexsci_id.toml"
         "--work_dir"
            help = "Working directory for NExScI query and L0 downloads."
            arg_type = String
            default = "L0"
         "--cookie"
            help = "Filename for NExScI cookie."
            arg_type = String
            default = "neidadmincookie.txt"
         "--query_filename"
            help = "Filename for L0 query results."
            arg_type = String
            default = "meta_l0.csv"
         "--root"
            help = "Path to root of data directories"
            arg_type = String
            default = ""
         "--try_tel_first"
            help = "Try pyrheliometer telemetry file before L0s"
            action = :store_true
         "--verbose"
            help = "Verbose outputs"
            action = :store_true
      end

  return parse_args(s)
end
args = parse_commandline_make_pyrheliometer_daily()
root_dir = args["root"]
output_fn = joinpath(root_dir,args["output"]) # if haskey(args,"output")     output_fn          = args["output"]      end
mkpath(dirname(output_fn))
tmp_path  = !isnothing(args["work_dir"]) ? joinpath(root_dir,args["work_dir"]) : Filesystem.mktempdir()
pyrohelio_dir = !isnothing(args["pyrheliometer_dir"]) ? joinpath(root_dir,args["pyrheliometer_dir"]) : nothing
nexsci_login_fn = !isnothing(args["nexsci_login_filename"]) ? joinpath(root_dir,args["nexsci_login_filename"]) : nothing
verbose   = args["verbose"]

using Dates
using CSV, DataFrames, Query

if !occursin(r"\.csv$", args["manifest_or_date"])
   manifest_fn = nothing
   m = match(r"(\d{4})(\d{2})(\d{2})", args["manifest_or_date"])
   if isnothing(m)
      @error "Invalid arguement.  Must be a .csv file or date."
   end
   #println("m = ", m)
   @assert length(m.captures) == 3
   year  = parse(Int64,m[1])
   month = parse(Int64,m[2])
   day   = parse(Int64,m[3])
   @assert 2020 <= year <= 2040
   @assert 1 <= month <= 12
   @assert 1 <= day <= 31
   target_datestr = args["manifest_or_date"]
   target_date = Date(m[1]*m[2]*m[3],DateFormat("yyyymmdd"))
   df_in = DataFrame()
else
   manifest_fn = joinpath(root_dir,args["manifest_or_date"])
   @assert filesize(manifest_fn) > 0
   df_in = CSV.read(manifest_fn,DataFrame)
   if hasproperty(df_in, :l0filename)   df_in = df_in |> @rename(:l0filename => :Filename) |> DataFrame  end
   if hasproperty(df_in, :l1filename)   df_in = df_in |> @rename(:l1filename => :Filename) |> DataFrame  end
   if hasproperty(df_in, :l2filename)   df_in = df_in |> @rename(:l2filename => :Filename) |> DataFrame  end
   if hasproperty(df_in, :filename)     df_in = df_in |> @rename(:filename => :Filename) |> DataFrame    end
   m = match(r"(\d{4})\/(\d{2})\/(\d{2})\/",manifest_fn)
   if isnothing(m) || (length(m.captures) != 3)
      target_datestr = nothing
      target_date = nothing
      println("# Empty ", manifest_fn, " and can't figure out date to query.  Giving up.")
      error(0)
   else
      @assert length(m.captures) == 3
      year  = parse(Int64,m[1])
      month = parse(Int64,m[2])
      day   = parse(Int64,m[3])
      @assert 2020 <= year <= 2040
      @assert 1 <= month <= 12
      @assert 1 <= day <= 31
      target_datestr = args["manifest_or_date"]
      target_date = Date(m[1]*m[2]*m[3],DateFormat("yyyymmdd"))
   end
end

using TOML
using NeidArchive

if (size(df_in,1) == 0) && (!isnothing(target_date))
   if isnothing(args["user"]) || isnothing(args["password"])
      nexsci_dict = TOML.parsefile(nexsci_login_fn)
      user_nexsci = nexsci_dict["user"]
      password_nexsci = nexsci_dict["password"]
      true
   else
      user_nexsci = args["user"]
      password_nexsci = args["password"]
      true
   end
   cookie_fn = args["cookie"]
   @assert !isnothing(user_nexsci) && !isnothing(password_nexsci) && !isnothing(cookie_fn)
   if !isdir(tmp_path)  mkpath(tmp_path) end
   query_fn  =  !isnothing(args["query_filename"]) ? joinpath(tmp_path,args["query_filename"]) : nothing

   NeidArchive.login(userid=user_nexsci, password=password_nexsci, cookiepath=cookie_fn)
   param = Dict{String,String}()
   param["datalevel"] = "solarl0"
   param["piname"] = "Mahadevan"
   param["object"] = "Sun"
   param["datetime"] = NeidArchive.datetime_one_day_solar(target_date)

   NeidArchive.query(param, cookiepath=cookie_fn, outpath=query_fn)
   num_lines = countlines(query_fn) - 1
   println("# Query resulted in file with ", num_lines, " entries.")
   manifest_fn = query_fn
   @assert filesize(manifest_fn) > 0
   df_in = CSV.read(manifest_fn,DataFrame)
   df_in = df_in |> @rename(:l0filename => :Filename) |> DataFrame
end

if args["try_tel_first"] && (!isnothing(pyrohelio_dir) && isdir(pyrohelio_dir))
   pyrohelio_files = Pyroheliometer.make_pyrohelio_file_dataframe(pyrohelio_dir)
   @assert hasproperty(df_in,"Filename")
   @assert hasproperty(df_in,"exptime")
   df_out = DataFrame(map(r->Pyroheliometer.get_pyrohelio_summary(pyrohelio_files, string(r.Filename), r.exptime), eachrow(df_in)))
else
   pyrohelio_files = DataFrame()
   df_out = DataFrame()
end

if (size(df_out,1) == 0) || any(ismissing.(df_out.mean_pyroflux))
   if !isdir(tmp_path)  mkpath(tmp_path) end
   #println("# start_row = ", 1, " end_row = ", num_lines-1)
   #NeidArchive.download(query_fn, param["datalevel"], outdir=tmp_path, cookiepath=cookie_fn, start_row=1, end_row=10)
   #NeidArchive.download(query_fn, "l0", outdir=tmp_path, cookiepath=cookie_fn, start_row=1, end_row=10)
   #NeidArchive.download(query_fn, "l0", outdir=tmp_path, cookiepath=cookie_fn)
   df_out = DataFrame(map(fn->Pyroheliometer.get_pyrohelio_summary(joinpath(tmp_path,replace(string(fn),"neidL2_"=>"neidL0_"))),df_in.Filename))
end

if (size(df_out,1) == 0) && !args["try_tel_first"] && (!isnothing(pyrohelio_dir) && isdir(pyrohelio_dir))
   pyrohelio_files = Pyroheliometer.make_pyrohelio_file_dataframe(pyrohelio_dir)
   @assert hasproperty(df_in,"Filename")
   @assert hasproperty(df_in,"exptime")
   df_out = DataFrame(map(r->Pyroheliometer.get_pyrohelio_summary(pyrohelio_files, string(r.Filename), r.exptime), eachrow(df_in)))
end

CSV.write(output_fn, df_out)
df_to_rm = df_in |> @filter( occursin(r"neidL0_\d{8}T\d+\.fits",_.Filename) ) |> @select(:Filename) |> @filter(isfile(_.Filename)) |> DataFrame
#map(fn->println("rm $fn"),df_to_rm.Filename)
exit(0)

#=

if (size(pyrohelio_files,1)>0) && (size(df_in,1)>0) # !isnothing(manifest_fn) && (filesize(manifest_fn) > 0)
      if hasproperty(df_in,"object")
         df_in = df_in |> @filter( _.object == "Sun" ) |> DataFrame
      elseif hasproperty(df_in,"target")
         df_in = df_in |> @filter( _.target == "Sun" ) |> DataFrame
      end
      filenamelist = hasproperty(df_in,"l2filename") ? df_in[!,"l2filename"] :
                     hasproperty(df_in,"l1filename") ? df_in[!,"l1filename"] :
                     hasproperty(df_in,"l0filename") ? df_in[!,"l0filename"] :
                     hasproperty(df_in,"filename") ? df_in[!,"filename"] :
                     hasproperty(df_in,"Filename") ? df_in[!,"Filename"] : nothing
      @assert !isnothing(filenamelist)
end

df_in = CSV.read(manifest_fn,DataFrame)

filenamelist = hasproperty(df_in,"l2filename") ? df_in[!,"l2filename"] :
               hasproperty(df_in,"l1filename") ? df_in[!,"l1filename"] :
               hasproperty(df_in,"l0filename") ? df_in[!,"l0filename"] :
               hasproperty(df_in,"filename") ? df_in[!,"filename"] :
               hasproperty(df_in,"Filename") ? df_in[!,"Filename"] : nothing
df_out = DataFrame(map(fn->Pyroheliometer.get_pyrohelio_summary(string(fn)),filenamelist))

using FITSIO
using Statistics
(;filename=basename(fn), time_start, exptime, mean_Δt, mean_pyroflux=mean_flux, rms_pyroflux=rms_flux, min_pyroflux=min_flux, max_pyroflux=max_flux)

df_in.mean_flux = Vector{Union{Float64,Missing}}(missing,size(df_in,1))
df_in.rms_flux = Vector{Union{Float64,Missing}}(missing,size(df_in,1))
df_in.min_flux = Vector{Union{Float64,Missing}}(missing,size(df_in,1))
df_in.max_flux = Vector{Union{Float64,Missing}}(missing,size(df_in,1))
for i in 1:size(df_in,1)
   file = hasproperty(df_in,"l2filename") ? df_in[i,"l2filename"] :
                  hasproperty(df_in,"l1filename") ? df_in[i,"l1filename"] :
                  hasproperty(df_in,"l0filename") ? df_in[i,"l0filename"] :
                  hasproperty(df_in,"filename") ? df_in[i,"filename"] :
                  hasproperty(df_in,"Filename") ? df_in[i,"Filename"] : nothing
   try
      m = match(r"neidL?_(\d+T\d+).fits",fninfull)
   	@assert m != nothing
   	time_start = DateTime(m[1],DateFormat("yyyymmddTHHMMSS"))

      f = FITS(fninfull)
      hdr = read_header(f[1])
   	exptime = hdr["EXPTIME"]
      t = read(f[21],"Time")
      v = read(f[21],"Voltage")
      fd = read(f[21],"FluxDensity")
      df_tmp = DataFrame(:Time=>t, :Voltage=>v, :FluxDensity=>fd)

      mean_Δt = get_pyrohelio_mean_Δt(df_tmp)
      #CSV.write(fnout, df)
      df_in.mean_flux[i] = mean(df_tmp.FluxDensity)
      df_in.rms_flux[i] = sqrt(var(df_tmp.FluxDensity, corrected=false))
      (df_in.min_flux[i], df_in.max_flux[i]) = extrema(df_tmp.FluxDensity)
   catch ex
      println("# Failed to extract pyrheliometer data from $fninfull.\n")
      continue
  end # try to read L0 fits file
end # for each L0 file

try
   mean_flux = mean(df.FluxDensity)
   rms_flux = sqrt(var(df.FluxDensity, corrected=false))
   outputline = fnin * "," * string(mean_flux) * "," * string(rms_flux) * "\n"
   write(exampleFileIOStream, outputline)
   if mod(i,10) == 0 flush(exampleFileIOStream)  end
 catch ex
    println("# Failed to write summary data for $fn_out\n.")
    continue
 end

target_date = !isnothing(args["date"]) ? DateTime(args["date"],DateFormat("yyyymmdd")) : nothing

Pyroheliometer.pick_pyrohelio_file(pyrohelio_files,target_date)
Pyroheliometer.get_pyrohelio_summary(pyrohelio_files, fn, exptime)

if !isnothing(args["pyrheliometer_dir"])
   pyrohelio_dir = arg["pyrheliometer_dir"]
   global df_pyrohelio_files = Pyroheliometer.make_pyrohelio_file_dataframe(pyrohelio_data_path)
   # Extract from pyrheliometer telimetry file
end

using CSV, DataFrames
using NeidArchive

if haskey(args,"root")       root_dir            = args["root"]   end
if haskey(args,"input")      input_dir           = args["input"]   end
if haskey(args,"output")     output_dir          = args["output"]      end
if haskey(args,"subdir")     subdir              = args["subdir"]      end
if haskey(args,"pyrohelio")  pyrohelio_dir       = args["pyrohelio"]   end
create_missing_continuum_files = haskey(args,"create-continuum") ? args["create-continuum"] : false
verbose = haskey(args,"verbose") ? args["verbose"] : false

if verbose && !isdefined(Main,:RvSpectML)  println("# Loading RvSpecML")    end
 using RvSpectMLBase
 using EchelleInstruments, EchelleInstruments.NEID
 using RvSpectML
 if verbose println("# Loading NeidSolarScripts")    end
 using NeidSolarScripts
 using SunAsAStar
 using SunAsAStar.SolarRotation
 using SunAsAStar.DifferentialExtinction
# using NeidSolarScripts.SolarRotation
# using NeidSolarScripts.DifferentialExtinction
 if verbose   println("# Loading other packages")    end
 using CSV, DataFrames, Query, Dates
 #using StatsBase, Statistics, Dates

fits_target_str = "Sun"
#=
 if !@isdefined(subdir)
      subdir = ""
 end
 if !@isdefined(create_missing_continuum_files)
     create_missing_continuum_files = false
 end
=#
 paths_to_search_for_param = [pwd(),pkgdir(NeidSolarScripts)]

 if verbose println("# Finding what data files are avaliable.")  end
 eval(read_data_paths(paths_to_search=paths_to_search_for_param))
 @assert isdefined(Main,:root_dir)
# @assert isdefined(Main,:neid_data_path)
 @assert isdefined(Main,:input_dir)
 @assert isdefined(Main,:output_dir)
 @assert isdefined(Main,:subdir)
 input_path = joinpath(root_dir,input_dir, subdir)
 output_path = joinpath(root_dir,output_dir, subdir)
 manifest_filename = joinpath(output_path,"manifest.csv")
 manifest_calib_filename = joinpath(output_path,"manifest_calib.csv")
 pyrohelio_data_path = joinpath(root_dir,pyrohelio_dir)

can_skip_generating_manifest = false
if isfile(manifest_filename) || islink(manifest_filename)
    df_files  = CSV.read(manifest_filename, DataFrame)
    warn_link =  islink(manifest_filename) ? " (a link)" : ""
    println("# Read ", manifest_filename, warn_link, ".")
    @assert size(df_files,1) >= 1
    @assert hasproperty(df_files,:Filename)
    @assert hasproperty(df_files,:target)
    @assert hasproperty(df_files,:bjd)
    @assert hasproperty(df_files,:ssbz)
    @assert hasproperty(df_files,:exptime)
    @assert hasproperty(df_files,:alt_sun)
    @assert hasproperty(df_files,:Δv_diff_ext)
    @assert hasproperty(df_files,:Δfwhm²)
    #=
    @assert hasproperty(df_files,:order_snrs)
    if eltype(df_files[!,:order_snrs]) == String
        df_files[!,:order_snrs] = map(i->parse.(Float64,split(df_files[i,:order_snrs][2:end-1],',')),1:size(df_files,1))
    end
    @assert eltype(df_files[!,:order_snrs]) == Vector{Float64}
    =#
    println("# Required fields present ", size(df_files), ".  No need to regenerate manifest.")
    can_skip_generating_manifest = true
    if any(isnan.(df_files[!,:Δfwhm²]))    ||
       any(isnan.(df_files[!,:Δv_diff_ext]))
       # any(isnan.(df_files[!,:order_snrs]))
          can_skip_generating_manifest = false
    end
else
    println("# Will need to generate manifest at ", manifest_filename, ".")
end

#can_skip_generating_manifest = false
if can_skip_generating_manifest && !create_missing_continuum_files
    exit()
end

try
   global df_pyrohelio_files = Pyroheliometer.make_pyrohelio_file_dataframe(pyrohelio_data_path)
catch ex
   println("# Error making pyroheliometer dataframe for " * pyrohelio_data_path * ".")
   exit(0)
end

try
   global df_files = NEID.make_manifest(input_path) #joinpath(root_path,input_dir, subdir))

catch ex
   println("# Error making manifest for ", subdir," and/or pyroheliometer dataframe.")
   touch(manifest_filename)
   touch(manifest_calib_filename)
   exit(0)
end

if size(df_files,1) == 0
   println("# No files in manifest.")
   touch(manifest_filename)
   touch(manifest_calib_filename)
   exit(0)
end

calibration_target_substrings = ["Etalon","LFC","ThArS","LDLS","FlatBB"]
df_files_calib = df_files |>
        @filter( any(map(ss->contains(_.target,ss),calibration_target_substrings)) ) |>
        DataFrame

df_files_calib[!,:airmass] = fill(NaN,size(df_files_calib,1))
CSV.write(manifest_calib_filename, df_files_calib)

try
    global df_files_obs = df_files |>
    @filter( _.target == "Sun" ) |>
    @filter( !isnothing(_.drpversion) && !ismissing(_.drpversion) && !(_.drpversion == "") ) |>
    #@filter( Base.thisminor(VersionNumber(_.drpversion)) == Base.thisminor(VersionNumber(1,1,0)) ) |>
    @filter( !any(map(ss->contains(_.target,ss),calibration_target_substrings)) ) |>
    # @take(10) |>
    DataFrame
catch ex
   println("# No files passed target and DRP version criteria.")
   touch(manifest_filename)
   exit(0)
end

#=
df_files_use = df_files_obs |>
      @filter( _.target == fits_target_str ) |>
      #@filter(bjd_first_good <= _.bjd < bjd_last_good) |>
      #@filter( is_good_day(_.Filename) ) |>
      #@take(max_spectra_to_use) |>
      DataFrame
=#

if size(df_files_obs,1) == 0
   touch(manifest_filename)
   println("# No files avaliable in df_files_obs.")
   exit(0)
end

if fits_target_str == "Sun" || fits_target_str == "Solar"
    df_pyrohelio_obs = DataFrame(map(x->Pyroheliometer.get_pyrohelio_summary(df_pyrohelio_files, x[1], x[2]), zip(df_files_obs.Filename, df_files_obs.exptime)  ))
    df_files_obs.filename = basename.(df_files_obs.Filename)
    df_files_use = leftjoin(df_files_obs, df_pyrohelio_obs, on=:filename, makeunique=true)

    df_sol = DataFrame(get_solar_info.(df_files_use.bjd,obs=:WIYN))
    df_files_use[!,:alt_sun] = df_sol[!,:alt]
    df_files_use[!,:airmass] = df_sol[!,:airmass]
    df_files_use[!,:hour_angle] = df_sol[!,:hour_angle]
    df_files_use[!,:sol_dist] = df_sol[!,:sol_dist_au]
    #df_files_use[!,:alt_sun] = calc_solar_alt.(df_files_use.bjd,obs=:WIYN)
    #df_files_use[!,:airmass] = DifferentialExtinction.f_airmass.(DifferentialExtinction.f_z.(deg2rad.(df_files_use.alt_sun)))
    println("# Computing differential extinction")
    try
        df_files_use[!,:Δv_diff_ext] = calc_Δv_diff_extinction.(df_files_use.bjd, obs=:WIYN)
    catch
       try
          println("# Encountered error with call to JplHorizons, sleeping and will try again.")
          sleep(60)
          df_files_use[!,:Δv_diff_ext] = calc_Δv_diff_extinction.(df_files_use.bjd, obs=:WIYN)
       catch
          println("# Encountered error with call to JplHorizons, using NaNs to proceed with rest of pipeline.")
          df_files_use[!,:Δv_diff_ext] = fill(NaN,length(df_files_use.bjd))
       end
    end
    println("# FWHM effect")
    #df_files_use[!,:Δfwhm²] = CalculateFWHMDifference_SolarRotation_from_obs(df_files_use.bjd,obs=:WIYN)
    df_files_use[!,:Δfwhm²]= calc_Δfwhm_solar_rotation.(df_files_use.bjd,obs=:WIYN)

else
    @warn("This script is intended for analyzing NEID solar observations.")
    df_files_use = df_files_obs |>
      @filter( _.target == fits_target_str ) |>
      #@filter(bjd_first_good <= _.bjd < bjd_last_good) |>
      #@filter( is_good_day(_.Filename) ) |>
      #@take(max_spectra_to_use) |>
      DataFrame
end

if size(df_files_use,1) == 0
   println("# No files avaliable.")
   touch(manifest_filename)
   exit(0)
end

if create_missing_continuum_files
   df_files_use.continuum_filename = map(fn->joinpath(output_path,"continuum",match(r"(neidL[1,2]_\d+[T_]\d+)\.fits$", fn)[1] * "_continuum=afs.jld2"),df_files_use.Filename)
end

#=
if verbose println("# Reading in customized parameters from param.jl.")  end
   if !@isdefined(idx_day_to_use)
       idx_day_to_use = 1
   end
   eval(code_to_include_param_jl(paths_to_search=paths_to_search_for_param))
   if match(r"neidL1_(\d+)[T_](\d+)\.fits$", first(df_files_use.Filename))[1] ==  match(r"neidL1_(\d+)[T_](\d+)\.fits$", last(df_files_use.Filename))[1]
      date_str = match(r"neidL1_(\d+)[T_](\d+)\.fits$", first(df_files_use.Filename))[1]
    else
      date_str = string(match(r"neidL1_(\d+)[T_](\d+)\.fits$", first(df_files_use.Filename))[1]) * "-" * string(match(r"neidL1_(\d+)[T_](\d+)\.fits$", last(df_files_use.Filename))[1])
   end
   outputs["df_files_use"] = df_files_use

   outputs_filename = joinpath(output_dir,"solar_" * date_str * "_new.jld2")
   if isfile(outputs_filename) && false
     times_already_processed = load(outputs_filename, "times")
     files_in_day_to_process = size(df_files_solar_by_day.data[idx_day_to_use],1)
      if files_in_day_to_process == length(times_already_processed)
         println("# Already processed all ", length(times_already_processed), " files for ", date_str)
         exit()
      end
   end
=#

#using NaNMath
compute_order_snr = false
if compute_order_snr
println("# Computing order SNRs")
df_files_use[!,:order_snrs] = fill(zeros(0),size(df_files_use,1))
   for (i,row) in enumerate(eachrow(df_files_use))
       #if row.target != "Sun" continue   end
       if any(map(ss->contains(row.target,ss),calibration_target_substrings))
           println("# Warning: target = ", row.target)
           continue
       end
       spec = NEID.read_data(row)
       all_orders = min_order(NEID2D()):max_order(NEID2D())
       pixels = min_pixel_in_order(NEID2D()):max_pixel_in_order(NEID2D())
       order_snr = map(ord->RvSpectMLBase.calc_snr(spec, pixels, ord), all_orders)
       df_files_use[i,:order_snrs] = order_snr
   end
end

println("# Writing manifest file.")
CSV.write(manifest_filename, df_files_use)
      # Can extract vector of Reals turned into a string via:
      # df2 = CSV.read("test.csv", DataFrame)
      # map(i->parse.(Float64,split(df2[i,:order_snrs][2:end-1],',')),1:size(df2,1))

#snr_matrix = hcat(df_files_use[!,:order_snrs]...)

#=
using JLD2, FileIO
if create_missing_continuum_files
   println("# Creating missing continuum files.")
   for row in eachrow(df_files_use)
     if isfile(row.continuum_filename) continue end
     println("# Need to make ", row.continuum_filename )

     output_filename = row.continuum_filename
     tmpdir = basename(output_filename)
     #=
     if !isdir(tmpdir)
        println("# Making ", tmpdir)
        mkpath(tmpdir)
     end
     =#
     #continue
     spec = NEID.read_data(row.Filename)
     try
        continuum = Continuum.calc_continuum_model(spec)
        @save output_filename continuum
     catch
        println("# ERROR computing continuum for ", output_filename)
     end
   end
end

=#
=#


text:
module GA
include("multivector.jl")
export Blade, @basis
end

ext: jl


text:
mutable struct Field

    _X_      :: AbstractArray{Float64, 2}
    _X       :: AbstractArray{Float64, 1}

    _b       :: AbstractArray{Float64, 1}

    _vel      :: AbstractArray{Float64, 1}
    _u        :: AbstractArray{Float64, 1}
    _v        :: AbstractArray{Float64, 1}
    _w        :: AbstractArray{Float64, 1}

    _Xflx_U_  :: AbstractArray{Float64, 2}
    _Xflx_V_  :: AbstractArray{Float64, 2}
    _Xflx_W_  :: AbstractArray{Float64, 2}

    _ADVX_    :: AbstractArray{Float64, 2}

    _WKRSTX_   :: AbstractArray{Float64, 2}
    _VDIFFX_   :: AbstractArray{Float64, 2}
    _QFLXX_    :: AbstractArray{Float64, 2}

    HMXL     :: AbstractArray{Float64, 3}

    SWFLX    :: AbstractArray{Float64, 3}
    NSWFLX   :: AbstractArray{Float64, 3}
    VSFLX    :: AbstractArray{Float64, 3}

    TAUX         :: AbstractArray{Float64, 3}
    TAUY         :: AbstractArray{Float64, 3}

    TAUX_east    :: AbstractArray{Float64, 3}
    TAUY_north   :: AbstractArray{Float64, 3}

    Q_FRZHEAT       :: AbstractArray{Float64, 3}
    Q_FRZMLTPOT_NEG :: AbstractArray{Float64, 3}
    Q_FRZMLTPOT     :: AbstractArray{Float64, 3}
    Q_LOST          :: AbstractArray{Float64, 3}

    # sugar view
    sv       :: Union{Nothing, Dict}

    function Field(
        ev :: Env,
    )

        Nz, Nx, Ny = ev.Nz, ev.Nx, ev.Ny

        T_pts = Nx * Ny * Nz
        U_pts = Nx * Ny * Nz
        V_pts = Nx * (Ny+1) * Nz
        W_pts = Nx * Ny * (Nz+1)

        _X_ = zeros(Float64, T_pts, 2)
        _X  = view(_X_, :)

        _b  = zeros(Float64, T_pts)

        sT_pts = Nx * Ny * 1

        _vel = zeros(Float64, U_pts + V_pts + W_pts)

        idx = 0;
        _u = view(_vel, (idx+1):(idx+U_pts)) ; idx+=U_pts
        _v = view(_vel, (idx+1):(idx+V_pts)) ; idx+=V_pts
        _w = view(_vel, (idx+1):(idx+W_pts)) ; idx+=W_pts

        _Xflx_U_ = zeros(Float64, U_pts, 2)
        _Xflx_V_ = zeros(Float64, V_pts, 2)
        _Xflx_W_ = zeros(Float64, W_pts, 2)

        _ADVX_ = zeros(Float64, T_pts, 2)

        _WKRSTX_ = zeros(Float64, T_pts, 2)

        _VDIFFX_ = zeros(Float64, T_pts, 2)
        _QFLXX_  = zeros(Float64, T_pts, 2)

        HMXL = zeros(Float64, 1, Nx, Ny)
        SWFLX = zeros(Float64, 1, Nx, Ny)
        NSWFLX = zeros(Float64, 1, Nx, Ny)
        VSFLX = zeros(Float64, 1, Nx, Ny)
        TAUX_east = zeros(Float64, 1, Nx, Ny)
        TAUY_north = zeros(Float64, 1, Nx, Ny)

        TAUX = zeros(Float64, 1, Nx, Ny)
        TAUY = zeros(Float64, 1, Nx, Ny)

        Q_FRZHEAT       = zeros(Float64, 1, Nx, Ny)
        Q_FRZMLTPOT_NEG = zeros(Float64, 1, Nx, Ny)
        Q_FRZMLTPOT     = zeros(Float64, 1, Nx, Ny)
        Q_LOST          = zeros(Float64, 1, Nx, Ny)

        fi = new(

            _X_,
            _X,

            _b,

            _vel,
            _u,
            _v,
            _w,

            _Xflx_U_,
            _Xflx_V_,
            _Xflx_W_,

            _ADVX_,
            _WKRSTX_,
            _VDIFFX_,
            _QFLXX_,

            HMXL,

            SWFLX,
            NSWFLX,
            VSFLX,

            TAUX,
            TAUY,

            TAUX_east,
            TAUY_north,

            Q_FRZHEAT,
            Q_FRZMLTPOT_NEG,
            Q_FRZMLTPOT,
            Q_LOST,

            nothing,
        )

        fi.sv = getSugarView(ev, fi)

        return fi

    end

end

function getSugarView(
    ev :: Env,
    fi :: Field,
)

    Nx, Ny, Nz = ev.Nx, ev.Ny, ev.Nz

    sv = Dict(
        :TEMP   => reshape(view(fi._X_, :, 1), Nz, Nx, Ny),
        :SALT   => reshape(view(fi._X_, :, 2), Nz, Nx, Ny),
        :_TEMP  => view(fi._X_, :, 1),
        :_SALT  => view(fi._X_, :, 2),
        :UVEL   => reshape(fi._u, Nz, Nx, Ny),
        :VVEL   => reshape(fi._v, Nz, Nx, Ny+1),
        :WVEL   => reshape(fi._w, Nz+1, Nx, Ny),
        :ADVT   => reshape(view(fi._ADVX_, :, 1), Nz, Nx, Ny),
        :ADVS   => reshape(view(fi._ADVX_, :, 2), Nz, Nx, Ny),
        :VDIFFT => reshape(view(fi._VDIFFX_, :, 1), Nz, Nx, Ny),
        :VDIFFS => reshape(view(fi._VDIFFX_, :, 2), Nz, Nx, Ny),
        :WKRSTT => reshape(view(fi._WKRSTX_, :, 1), Nz, Nx, Ny),
        :WKRSTS => reshape(view(fi._WKRSTX_, :, 2), Nz, Nx, Ny),
        :QFLXT  => reshape(view(fi._QFLXX_, :, 1), Nz, Nx, Ny),
        :QFLXS  => reshape(view(fi._QFLXX_, :, 2), Nz, Nx, Ny),
    )

    sv[:SST] = view(sv[:TEMP], 1:1, :, :)
    sv[:SSS] = view(sv[:SALT], 1:1, :, :)

    return sv
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
# Julia code
# Copyright (c) 2022 the SLICOTMath.jl developers
# Portions extracted from SLICOT-Reference distribution:
# Copyright (c) 2002-2020 NICONET e.V.

function run_mb04ud(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    MMAX = 20
    NMAX = 20
    LDA = MMAX
    LDE = MMAX
    LDQ = MMAX
    LDZ = NMAX
    LDWORK = max( NMAX,MMAX )
    A = Array{Float64,2}(undef, LDA,NMAX)
    E = Array{Float64,2}(undef, LDE,NMAX)
    Q = Array{Float64,2}(undef, LDQ,MMAX)
    Z = Array{Float64,2}(undef, LDZ,NMAX)
    ISTAIR = Array{BlasInt,1}(undef, MMAX)
    DWORK = Array{Float64,1}(undef, LDWORK)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    M = parse(BlasInt, vs[1])
    N = parse(BlasInt, vs[2])
    TOL = parse(Float64, replace(vs[3],'D'=>'E'))
if ( M<0 || M>MMAX )
elseif ( N<0 || N>NMAX )
else

    vs = String[]
    _isz,_jsz = (M,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       A[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end

    vs = String[]
    _isz,_jsz = (M,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       E[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
    JOBQ = 'N'
    JOBZ = 'N'
# interp call 1

    RANKE, INFO = SLICOT.mb04ud!(JOBQ, JOBZ, M, N, A, E, Q, Z, ISTAIR, TOL)
    println(io, "RANKE = $RANKE")
    @test INFO == 0
    INFO == 0 || return

if ( INFO!=0 )
else
# interp output 1
    println(io,"A:")
    _nc = N
    _nr = M
    show(io,"text/plain",A[1:_nr,1:_nc])
    println(io,)

# interp output 2
    println(io,"E:")
    _nc = N
    _nr = M
    show(io,"text/plain",E[1:_nr,1:_nc])
    println(io,)

# interp output 3
    println(io,"ISTAIR:")
    _nr = M
    show(io,"text/plain",ISTAIR[1:_nr])
    println(io,)

end # if
end # if
    close(f)
end # run_X()


text:
export distance, angular_distance, cross_track_distance, along_track_distance

const tolerance_m = 0.001

"""
    distance(pos₁::Point, pos₂::Point[, radius::Float64=Rₑ_m])

Return the `distance` [m] of the great circle line between the positions `pos₁`
[deg] and `pos₂` [deg] on a sphere, with an optionally given `radius`.

Source: www.movable-type.co.uk/scripts/latlong.html
"""
function distance(pos₁::Point, pos₂::Point, radius::Float64=Rₑ_m)
    Δpos = pos₂ - pos₁
    a = sind(Δpos.ϕ/2.0)^2 + cosd(pos₁.ϕ)*cosd(pos₂.ϕ)*sind(Δpos.λ/2.0)^2
    c = 2.0 * atan(√a, √(1.0 - a))
    return radius * c
end

"""
    distance(section::RouteSection [, radius::Float64=Rₑ_m])

Return the `distance` [m] of the great circle line on the route section on a
sphere, with an optionally given `radius`.

Source: www.movable-type.co.uk/scripts/latlong.html
"""
distance(section::RouteSection, radius::Float64=Rₑ_m) = distance(section.pos₁,
section.pos₂, radius)

"""
    angular_distance(distance::Float64[, radius::Float64=Rₑ_m])

Return the `angular_distance` in [deg] on a sphere, with an optionally given
`radius`.

Source: www.movable-type.co.uk/scripts/latlong.html
"""
function angular_distance(distance::Float64, radius::Float64=Rₑ_m)
    return rad2deg(distance/radius)
end

"""
    angular_distance(pos₁::Point, pos₂::Point, radius::Float64=Rₑ_m)

Return the `angular_distance` in [deg] of the great circle line between the
positions `pos₁` [deg] and `pos₂` [deg] on a sphere, with an optionally given
`radius`.

Source: www.movable-type.co.uk/scripts/latlong.html
"""
angular_distance(pos₁::Point, pos₂::Point, radius::Float64=Rₑ_m) =
angular_distance(distance(pos₁, pos₂, radius), radius)

"""
    cross_track_distance(pos₁::Point, pos₂::Point, pos₃::Point[, radius::Float64=Rₑ_m])

Return the `cross_track_distance` [m] from a point `pos₃` [deg] to a great
circle path defined by the points `pos₁` and `pos₂` [deg].

Source: www.movable-type.co.uk/scripts/latlong.html
"""
function cross_track_distance(pos₁::Point, pos₂::Point, pos₃::Point,
    radius::Float64=Rₑ_m)
    ∠distance₁₃ = angular_distance(pos₁, pos₃, radius)
    # ∠distance₁₃ = distance(pos₁, pos₃ , radius) / radius
    bearing₁₃ = bearing(pos₁, pos₃)
    bearing₁₂ = bearing(pos₁, pos₂)
    return cross_track_distance(∠distance₁₃, bearing₁₂, bearing₁₃, radius)
end

"""
    cross_track_distance(section::RouteSection, pos₃::Point [,
    radius::Float64=Rₑ_m])

Return the `cross_track_distance` [m] from a point `pos₃` [deg] to a route
section.

Source: www.movable-type.co.uk/scripts/latlong.html
"""
cross_track_distance(section::RouteSection, pos₃::Point,
radius::Float64=Rₑ_m) = cross_track_distance(section.pos₁, section.pos₂, pos₃,
radius)

"""
    cross_track_distance(pos₁::Point, bearing::Float64, pos₃::Point[, radius::Float64=Rₑ_m])

Return the `cross_track_distance` [m] from a point `pos₃` [deg] to a great
circle path defined by the point `pos₁` and `bearing` [deg].

Source: www.movable-type.co.uk/scripts/latlong.html
"""
function cross_track_distance(pos₁::Point, bearing::Float64, pos₃::Point,
    radius::Float64=Rₑ_m)
    pos₂ = destination_point(pos₁, radius, bearing, radius)
    return cross_track_distance(pos₁, pos₂, pos₃, radius)
end

"""
    cross_track_distance(∠distance₁₃::Float64, bearing₁₂::Float64,
    bearing₁₃::Float64[, radius::Float64=Rₑ_m])

Return the `cross_track_distance` [m] from a point `pos₃` [deg] to a great circle
path. The distance is calculated using the angular distance between points `pos₁`
and `pos₃` [deg], the `bearing` [deg] between points `pos₁` and `pos₃` [deg]
, the `bearing` [deg] between points `pos₁` and `pos₂` [deg], and the
`radius` of the earth [m].

Source: www.movable-type.co.uk/scripts/latlong.html
"""
function cross_track_distance(∠distance₁₃::Float64, bearing₁₂::Float64,
    bearing₁₃::Float64, radius::Float64=Rₑ_m)
    return asin(sind(∠distance₁₃) * sind(bearing₁₃ - bearing₁₂)) * radius
end

"""
    along_track_distance(pos₁::Point, pos₂::Point, pos₃::Point[,radius::Float64=Rₑ_m])

The `along_track_distance` from the start point `pos₁` [deg] to the closest
point on the great circle path (defined by points `pos₁` and `pos₂` [deg]) to
the point `pos₃` [deg]. Optionally a different `radius` for the earth can be used.

Source: www.movable-type.co.uk/scripts/latlong.html
"""
function along_track_distance(pos₁::Point, pos₂::Point, pos₃::Point,
    radius::Float64=Rₑ_m)
    ∠distance₁₃ = angular_distance(pos₁, pos₃, radius)
    cross_track_∠distance = rad2deg(cross_track_distance(pos₁, pos₂, pos₃, radius)/radius)
    return along_track_distance(∠distance₁₃, cross_track_∠distance, radius)
end

"""
    along_track_distance(∠distance₁₃::Float64, cross_track_∠distance::Float64,
    [radius::Float64=Rₑ_m])

The `along_track_distance` from the start point `pos₁` [deg] to the closest
point on the great circle path (defined by angular distance `∠distance₁₃` [deg]
between `pos₁` and `pos₃` [deg] and the cross track angular distance
`cross_track_∠distance` [deg]) to the point `pos₃` [deg]. Optionally a
different `radius` for the earth can be used.

Source: www.movable-type.co.uk/scripts/latlong.html
"""
function along_track_distance(∠distance₁₃::Float64,
    ∠cross_track_distance::Float64, radius::Float64=Rₑ_m)
    return acos(cosd(∠distance₁₃) / cosd(∠cross_track_distance)) * radius
end


text:
using PyPlot
using FFTW
using Bolt
using LinearAlgebra

N = 64
μ = 0
q = 0.0
r = 10 .^ range(-4., 4., length=N)

pl = Bolt.plan_fftlog(r, μ, q, 1.0; kropt=true)
aₙ = r .^ (μ + 1) .* exp.(-r.^2 / 2)
# aₙ = exp.(-(log.(r) ).^2 / 2)
y = similar(r, ComplexF64)
y2 = similar(r, ComplexF64)
mul!(y, pl, aₙ)

clf()
plt.plot(r, aₙ)
# plt.plot(r, aₙ .* r.^(-q))
# ylim(1e-5, 1.2)
plt.plot(pl.k, y)
yscale("symlog", linthreshy=1e-5)
# yscale("log")
xscale("log")
gcf()

##

# y .*= sin.(0.5 .* r)
ldiv!(y2,pl,y)

clf()
plt.plot(r, aₙ)
plt.plot(pl.r, y2)
# yscale("symlog", linthreshy=1e-5)
xscale("log")
gcf()

##
aₙ = sf[1200:end,1]
y = similar(r, ComplexF64)
r = x2a.(bg.x_grid[1200:end])
μ = 30
q = 0.9
pl = Bolt.plan_fftlog(r, μ, q, 1.0; kropt=true)
mul!(y, pl, aₙ)

clf()
plt.plot(r, aₙ)
plt.plot(r, y)

yscale("symlog", linthreshy=maximum(norm.(y)) * 1e-2)
xscale("log")
gcf()

##
using DelimitedFiles
fftdata = readdlm("test/fftlog64.txt", ' ', Float64, '\n')
k_ref = fftdata[:,1]
f_ref = fftdata[:,2]

##
clf()
plot(k_ref .- fftlogfreq(pl))
gcf()
##

isapprox(y, f_ref)

##
clf()
plot(y .- f_ref)
gcf()

##
clf()
plt.plot(fftlogfreq(pl), reverse(y))
plt.plot(k_ref, f_ref)
yscale("symlog", base=10, linthresh=1e-5)
xscale("log")
gcf()


text:
export @check!, @einsum!, @vielsum!, @tensor!

"""
    @check!(A[i, j, μ, ν])

Adds `A` to the store of known tensors, and records that it expects indices `i,j,μ,ν`.
If it is already in the store, then instead this checks whether the present indices differ
from the saved ones. Only the first letter is examined: `α` and `α2` are similar, as are nearby
letters `β`, `γ3`. More complicated indices like `Z[(i,j), -k, _, 3]` will be ignored.
This happens while parsing your source code, there is zero run-time penalty. Returns `A`.

In addition, with `size=true` option, it can insert size checks to be performed at run-time,
by returning `check!(A, stuff)`.
At the first occurrance this saves `i: => size(A,1)` etc., and on subsequent uses of
the same index (even on different tensors) it gives an error if the sizes do not match.
Here the whole index is used: `α`, `β` and `β2` may have different ranges.
This will need to look up indices in a dictionary, which takes ... 50ns, really?

    @check! B[i,j] C[μ,ν] D[i] E[j]

Checks several tensors, returns nothing.

    @check!  alpha=true  tol=3  size=false  throw=false  info  empty

Controls options for `@check!` and related macros (`@shape!`, `@reduce!`, `@einsum!` etc).
These are the default settings:
* `alpha=true` turns on the parse-time checking, based on index letters.
* `tol=3` sets how close together letters must be: `B[j,k]` is not an error but `B[a,b]` will be.
* `size=false` turns off run-time size checking.
* `throw=false` means that errors are given using `@error`, without interrupting your program.
* `empty` deletes all saved letters and sizes -- there is one global store for each, for now.
* `info` prints what's currently saved.

    @cast! B[i,j] := D[i] * E[j]
    @reduce! B[i,j] := sum(μ,ν) A[i,j,μ,ν] / C[μ,ν]

    @einsum!  B[i,j] := D[i] * E[j]
    @vielsum! B[i,j] := D[i] * E[j]
    @tensor!  B[i,j] := D[i] * E[j]

Versions of the macros from this package, and from Einsum.jl and TensorOperations.jl,
which call `@check!` on each of their tensors, before proceeding as normal.
"""
macro check!(exs...)
    where = (mod=__module__, src=__source__, str=unparse("@check!", exs...))
    _check!(exs...; where=where)
end

const index_store = Dict{Symbol, Tuple}()

const size_store  = Dict{Symbol, Int}() # TODO: alter check!() to avoid dictionary lookup, the macro
                                        # can push nothing to size_store::Vector, and index...

mutable struct CheckOptions
    alpha::Bool
    tol::Int
    size::Bool
    throw::Bool
end

const check_options = CheckOptions(true, 3, false, false)

function _check!(exs...; where=nothing)
    for ex in exs
        if @capture(ex, A_[vec__])
            if length(exs)==1
                return esc(check_one(ex, where))
            else
                check_one(ex, where)
            end

        elseif @capture(ex, alpha=val_Bool)
            check_options.alpha = val
        elseif @capture(ex, tol=val_Int)
            check_options.tol = val

        elseif @capture(ex, size=val_Bool)
            check_options.size = val

        elseif @capture(ex, throw=val_Bool)
            check_options.throw = val

        elseif ex == :empty
            empty!(index_store)
            empty!(size_store)
            @info "@check! stores emptied"

        elseif ex == :info
            @info "@check! info" check_options index_store size_store

        else
            @error "@check! doesn't know what to do with $ex"

        end
    end
    return nothing
end

#="""
    check_one(A[i,j,k], (mod=Module, src=...))

Does the work of `@check!`, on one index expression,
returning `A` or `check!(A,...)` according to global flags.
"""=#
function check_one(ex, where=nothing)
    @capture(ex, A_[vec__]) || error("check_one can't understand $ex, expected something like A[i,j]")
    ind = Tuple(vec)

    if !isa(A, Symbol) # ignore things like f(x)[i,j]
        return A
    end

    if check_options.alpha
        got = get!(index_store, A, ind)

        if length(ind) > length(got)
            check_err("@check! $ex now has more indices than previous $got", where)
        elseif length(ind) < length(got)
            check_err("@check! $ex now has fewer indices than previous $got", where)
        else
            for (i, j) in zip(ind, got)
                if !isa(i, Symbol) || i==(:_) # ignore Int, (i,j), etc.
                    continue
                end
                si = String(i)
                sj = String(j)
                length(si)>1 || length(sj)>1 && continue
                if abs(Int(si[1])-Int(sj[1])) > check_options.tol
                    check_err("@check! $ex now has index $i where previously it had $j", where)
                end
            end
        end
    end

    if check_options.size
        Astring = string(A,"[",join(ind," ,"),"]")
        return :( TensorCast.check!($A, $ind, $Astring, $where) )
    else
        return A
    end
end

function check_err(str::String, where=nothing)
    if check_options.throw
        error(str)
    elseif where==nothing
        @error str
    else
        @error str  _module=where.mod  _line=where.src.line  _file=string(where.src.file)
    end
end

"""
    check!(A, (:i,:j), "A[i,j]", (mod=..., src=...))

Performs run-time size checking, on behalf of the `@check!` macro, returns `A`.
The string and tuple are just for the error message.
"""
function check!(A::AbstractArray{T,N}, ind::Tuple, str::String, where=nothing) where {T,N}
    if N != length(ind)
        check_err("check! expected $str, but got ndims = $N", where)
    else
        for (d,i) in enumerate(ind)
            if !isa(i, Symbol) || i==(:_) # ignore Int, (i,j), etc.
                continue
            end
            sizeAd = size(A,d)
            got = get!(size_store, i, sizeAd)
            got == sizeAd || check_err("check! $str, index $i now has range $sizeAd instead of $got", where)
        end
    end
    A
end

"""
    @einsum! A[i,j] := B[i,k] * C[k,j]

Variant of `@einsum` from package Einsum.jl,
equivalent to wrapping every tensor with `@check!()`.
"""
macro einsum!(ex)
    where = (mod=__module__, src=__source__, str=unparse("@einsum!", exs...))
    _einsum!(ex, where)
end

"""
    @vielsum! A[i,j] := B[i,k] * C[k,j]

Variant of `@vielsum` from package Einsum.jl,
equivalent to wrapping every tensor with `@check!()`.
"""
macro vielsum!(ex)
    where = (mod=__module__, src=__source__, str=unparse("@vielsum!", exs...))
    _einsum!(ex, where; threads=true)
end

"""
    @tensor! A[i,j] := B[i,k] * C[k,j]

Variant of `@tensor` from package TensorOperations.jl,
equivalent to wrapping every tensor with `@check!()`.
"""
macro tensor!(ex)
    where = (mod=__module__, src=__source__, str=unparse("@tensor!", exs...))
    _tensor!(ex, where)
end

function _tensor!(ex, where=nothing)
    if @capture(ex, lhs_ := rhs_ ) || @capture(ex, lhs_ = rhs_ )

        outex = quote end
        function f(x)
            if @capture(x, A_[ijk__] )
                push!(outex.args, check_one(x, where))
            end
            x
        end
        MacroTools.prewalk(f, rhs)

        if check_options.size == false
            check_one(lhs) # then these are only parse checks, outex is trash
        else
            push!(outex.args, :(out = TensorOperations.@tensor $ex) )
            push!(outex.args, check_one(lhs, where)) # lhs size may not be known until after @tensor
            push!(outex.args, :out )                 # make sure we still return what @tensor did
            return esc(outex)
        end
    else
        @warn "@tensor! not smart enough to process $ex yet, so ignoring checks"
    end
    return esc(:( TensorOperations.@tensor $ex ))
end

function _einsum!(ex, where=nothing; threads = false)
    if @capture(ex, lhs_ := rhs_ ) || @capture(ex, lhs_ = rhs_ )

        outex = quote end
        function f(x)
            if @capture(x, A_[ijk__] )
                push!(outex.args, check_one(x, where))
            end
            x
        end
        MacroTools.prewalk(f, rhs)

        if check_options.size == false
            check_one(lhs) # then these are only parse checks, outex is trash
        else
            if threads
                push!(outex.args, :(out = Einsum.@vielsum $ex) )
            else
                push!(outex.args, :(out = Einsum.@einsum $ex) )
            end
            push!(outex.args, check_one(lhs, where)) # lhs size may not be known until after @einsum
            push!(outex.args, :out )                 # make sure we still return what @einsum did
            return esc(outex)
        end
    else
        @warn "@einsum! not smart enough to process $ex yet, so ignoring checks"
    end
    if threads
        return esc(:( Einsum.@vielsum $ex ))
    else
        return esc(:( Einsum.@einsum $ex ))
    end
end

#==

using MacroTools, Einsum, TensorOperations
B = rand(2,3); C = rand(3,2);
A = B * C
@einsum A[i,j] := B[i,k] * C[k,j]
@tensor A[i,j] := B[i,k] * C[k,j]

using TensorCast

@check! size=true throw=false info

@einsum! A[i,j] := B[i,k] * C[k,j]

@tensor! A[i,j] := B[i,k] * C[k,j]

@check! info   # has ABC and ijk

@check! A[z,j] # compains about z
@check! B[i]   # complains about number

B5 = rand(2,5); C5 = rand(5,2);
@einsum! A[i,j] := B5[i,k] * C5[k,j] # complains about sizes
@tensor! A[i,j] := B5[i,k] * C5[k,j]

@einsum! A[i,j] := B5[i,k] * C5[k,zz] # complains about zz
@tensor! A[i,j] := B5[i,k] * C5[k,zz] # and errors

==#


text: abstract type AbstractIterativeSampler end


text:
module SingleIterationModule

import ..CoreModule: Options, Dataset, RecordType, string_tree
import ..EquationUtilsModule: compute_complexity
import ..UtilsModule: debug
import ..SimplifyEquationModule: simplify_tree, combine_operators
import ..PopMemberModule: copy_pop_member, generate_reference
import ..PopulationModule: Population, finalize_scores, best_sub_pop
import ..HallOfFameModule: HallOfFame
import ..RegularizedEvolutionModule: reg_evol_cycle
import ..ConstantOptimizationModule: optimize_constants
import ..RecorderModule: @recorder

# Cycle through regularized evolution many times,
# printing the fittest equation every 10% through
function s_r_cycle(
    dataset::Dataset{T},
    baseline::T,
    pop::Population,
    ncycles::Int,
    curmaxsize::Int,
    frequencyComplexity::AbstractVector{T};
    verbosity::Int=0,
    options::Options,
    record::RecordType,
)::Tuple{Population,HallOfFame,Float64} where {T<:Real}
    max_temp = T(1.0)
    min_temp = T(0.0)
    if !options.annealing
        min_temp = max_temp
    end
    all_temperatures = LinRange(max_temp, min_temp, ncycles)
    best_examples_seen = HallOfFame(options)
    num_evals = 0.0

    for temperature in all_temperatures
        pop, tmp_num_evals = reg_evol_cycle(
            dataset,
            baseline,
            pop,
            temperature,
            curmaxsize,
            frequencyComplexity,
            options,
            record,
        )
        num_evals += tmp_num_evals
        for member in pop.members
            size = compute_complexity(member.tree, options)
            score = member.score
            if size <= options.maxsize && (
                !best_examples_seen.exists[size] ||
                score < best_examples_seen.members[size].score
            )
                best_examples_seen.exists[size] = true
                best_examples_seen.members[size] = copy_pop_member(member)
            end
        end
    end

    return (pop, best_examples_seen, num_evals)
end

function optimize_and_simplify_population(
    dataset::Dataset{T},
    baseline::T,
    pop::Population,
    options::Options,
    curmaxsize::Int,
    record::RecordType,
)::Tuple{Population,Float64} where {T<:Real}
    array_num_evals = zeros(Float64, pop.n)
    do_optimization = rand(pop.n) .< options.optimize_probability
    @inbounds @simd for j in 1:(pop.n)
        pop.members[j].tree = simplify_tree(pop.members[j].tree, options)
        pop.members[j].tree = combine_operators(pop.members[j].tree, options)
        if options.shouldOptimizeConstants && do_optimization[j]
            pop.members[j], array_num_evals[j] = optimize_constants(
                dataset, baseline, pop.members[j], options
            )
        end
    end
    num_evals = sum(array_num_evals)
    pop, tmp_num_evals = finalize_scores(dataset, baseline, pop, options)
    num_evals += tmp_num_evals

    # Now, we create new references for every member,
    # and optionally record which operations occurred.
    for j in 1:(pop.n)
        old_ref = pop.members[j].ref
        new_ref = generate_reference()
        pop.members[j].parent = old_ref
        pop.members[j].ref = new_ref

        @recorder begin
            # Same structure as in RegularizedEvolution.jl,
            # except we assume that the record already exists.
            @assert haskey(record, "mutations")
            member = pop.members[j]
            if !haskey(record["mutations"], "$(member.ref)")
                record["mutations"]["$(member.ref)"] = RecordType(
                    "events" => Vector{RecordType}(),
                    "tree" => string_tree(member.tree, options),
                    "score" => member.score,
                    "loss" => member.loss,
                    "parent" => member.parent,
                )
            end
            optimize_and_simplify_event = RecordType(
                "type" => "tuning",
                "time" => time(),
                "child" => new_ref,
                "mutation" => RecordType(
                    "type" =>
                        if (do_optimization[j] && options.shouldOptimizeConstants)
                            "simplification_and_optimization"
                        else
                            "simplification"
                        end,
                ),
            )
            death_event = RecordType("type" => "death", "time" => time())

            push!(record["mutations"]["$(old_ref)"]["events"], optimize_and_simplify_event)
            push!(record["mutations"]["$(old_ref)"]["events"], death_event)
        end
    end
    return (pop, num_evals)
end

end


text:
using TextEncodeBase
using Test

# quick and dirty macro for making @inferred as test case
macro test_inferred(ex)
    esc(quote
        @test begin
            @inferred $ex
            true
        end
    end)
end

using TextEncodeBase: AbstractTokenizer, AbstractTokenization,
    BaseTokenization, NestedTokenizer, FlatTokenizer,
    WordTokenization, IndexedTokenization, MatchTokenization,
    TokenStages, Document, Sentence, Word, Token, Batch
using TextEncodeBase: getvalue, getmeta, with_head_tail, trunc_and_pad, trunc_or_pad, nested2batch, nestedcall

using WordTokenizers

const ATR = AbstractTokenizer
const AT = AbstractTokenization
const BT = BaseTokenization

struct CharTk <: BT end
TextEncodeBase.splitting(::CharTk, x::Word) = split(x.x, "")
TextEncodeBase.splittability(::CharTk, x::Word) = TextEncodeBase.Splittable()

@testset "TextEncodeBase.jl" begin
    @testset "Tokenize" begin
        document = Document("This is the first sentence. And the second one with some number 12345.")
        sentence = Sentence("A single sentence with 31 char.")
        word = Word("word")

        @testset "base tokenizer" begin
            tkr = FlatTokenizer()
            @test tkr(document) == map(Token, mapfoldl(nltk_word_tokenize, append!, split_sentences(document.x)))
            @test tkr(sentence) == map(Token, nltk_word_tokenize(sentence.x))
            @test tkr(word) == [Token(word.x)]
        end

        @testset "edge case" begin
            tkr = FlatTokenizer()
            @test tkr(Document("")) == []
            @test tkr(Sentence("")) == []
            @test tkr(Word("")) == []
            @test tkr(Token("")) == []
        end

        @testset "word tokenizer" begin
            tkr = FlatTokenizer(WordTokenization(tokenize=poormans_tokenize))
            tkr2 = FlatTokenizer()
            @test tkr(document) == map(Token, mapfoldl(poormans_tokenize, append!, split_sentences(document.x)))
            @test tkr(sentence) == map(Token, poormans_tokenize(sentence.x))
            @test tkr(word) == [Token(word.x)]
            @test tkr(document) != tkr2(document)
            @test tkr(sentence) != tkr2(sentence)
        end

        @testset "index tokenizer" begin
            tkr = FlatTokenizer(IndexedTokenization())
            @test tkr(document) == begin
                sentences = split_sentences(document.x)
                words = map(x->nltk_word_tokenize(x), sentences)
                tokens = Token[]
                for (i, s) in enumerate(words)
                    for (j, w) in enumerate(s)
                        push!(tokens, Token(w, (sentence_id = i, word_id = j, token_id = j)))
                    end
                end
                tokens
            end
            @test tkr(sentence) == begin
                words = nltk_word_tokenize(sentence.x)
                tokens = Token[]
                for (i, w) in enumerate(words)
                    push!(tokens, Token(w, (word_id = i, token_id = i)))
                end
                tokens
            end
            @test tkr(word) == [Token(word.x, (word_id = 1, token_id  = 1))]
        end

        @testset "match tokenizer" begin
            tkr = FlatTokenizer(MatchTokenization([r"\d", r"en"]))
            @test map(getvalue, tkr(document)) == [
                "This", "is", "the", "first", "s",
                "en", "t", "en", "ce", ".", "And",
                "the", "second", "one", "with", "some",
                "number", "1", "2", "3", "4", "5", ".",
            ]
            @test map(getvalue, tkr(sentence)) == [
                "A", "single", "s", "en", "t", "en",
                "ce", "with", "3", "1", "char", ".",
            ]
            @test map(getvalue, tkr(word)) == [word.x]
            @test map(getvalue, tkr(Word("123"))) == ["1", "2", "3"]
        end

        @testset "indexed match tokenizer" begin
            tkr = FlatTokenizer(IndexedTokenization(MatchTokenization([r"\d", r"en"])))
            @test map(getvalue, tkr(document)) == [
                "This", "is", "the", "first", "s",
                "en", "t", "en", "ce", ".", "And",
                "the", "second", "one", "with", "some",
                "number", "1", "2", "3", "4", "5", ".",
            ]
            @test map(getmeta, tkr(document)) == begin
                s = Iterators.flatten((Iterators.repeated(1, 10), Iterators.repeated(2, 13)))
                w = Iterators.flatten((1:10, 1:13))
                map(NamedTuple{(:sentence_id, :word_id, :token_id)}, zip(s, w, w))
            end

            @test map(getvalue, tkr(sentence)) == [
                "A", "single", "s", "en", "t", "en",
                "ce", "with", "3", "1", "char", ".",
            ]
            @test map(getmeta, tkr(sentence)) == map(NamedTuple{(:word_id, :token_id)}, zip(1:12, 1:12))
            @test map(getvalue, tkr(word)) == [word.x]
            @test map(getmeta, tkr(word)) == [(word_id = 1, token_id = 1)]
            @test map(getvalue, tkr(Word("123"))) == ["1", "2", "3"]
            @test map(getmeta, tkr(Word("123"))) == map(NamedTuple{(:word_id, :token_id)}, zip(1:3, 1:3))
        end

        @testset "nested output" begin
            tkr = NestedTokenizer(IndexedTokenization())
            @test tkr(document) == begin
                sentences = split_sentences(document.x)
                words = map(x->nltk_word_tokenize(x), sentences)
                tokens = []
                for (i, s) in enumerate(words)
                    push!(tokens, map(enumerate(s)) do (j, w)
                          Token(w, (sentence_id = i, word_id = j, token_id = j))
                          end)
                end
                tokens
            end
            @test tkr(sentence) == begin
                words = nltk_word_tokenize(sentence.x)
                tokens = [map(enumerate(words)) do (i, w)
                          Token(w, (word_id = i, token_id = i))
                          end]
                tokens
            end
            @test tkr(word) == [Token(word.x, (word_id = 1, token_id  = 1))]
        end

        @testset "indexed char" begin
            tkr = FlatTokenizer(IndexedTokenization(CharTk()))
            @test tkr(document) == begin
                sentences = split_sentences(document.x)
                words = map(x->nltk_word_tokenize(x), sentences)
                tokens = Token[]
                for (i, s) in enumerate(words)
                    k = 1
                    for (j, w) in enumerate(s)
                        for c in split(w, "")
                            push!(tokens, Token(c, (sentence_id = i, word_id = j, token_id = k)))
                            k += 1
                        end
                    end
                end
                tokens
            end
            @test tkr(sentence) == begin
                words = nltk_word_tokenize(sentence.x)
                tokens = Token[]
                for (i, w) in enumerate(words)
                    for c in split(w, "")
                        push!(tokens, Token(c, (word_id = i, token_id = length(tokens)+1)))
                    end
                end
                tokens
            end
            @test tkr(word) == begin
                chars = split(word.x, "")
                tokens = Token[]
                for (i, c) in enumerate(chars)
                    push!(tokens, Token(c, (word_id = 1, token_id = i)))
                end
                tokens
            end
        end

        @testset "nested indexed match char" begin
            tkr = NestedTokenizer(IndexedTokenization(MatchTokenization(CharTk(), [r"\d", r"en"])))
            s(x) = split(x, "")
            r(x, n) = repeat(x:x, n)
            @test nestedcall(getvalue, tkr(document)) == [
                [
                    s("This"); s("is"); s("the"); s("first");
                    "s"; "en"; "t"; "en"; s("ce"); ".";
                ],
                [
                    s("And"); s("the"); s("second"); s("one"); s("with");
                    s("some"); s("number"); "1"; "2"; "3"; "4"; "5"; ".";
                ]
            ]
            @test nestedcall(getmeta, tkr(document)) == begin
                sentence_id = [r(1, 21), r(2, 35)]
                word_id = [[
                        r(1, 4); r(2, 2); r(3, 3); r(4, 5);
                        5; 6; 7; 8; r(9, 2); 10;
                    ], [
                        r(1, 3); r(2, 3); r(3, 6); r(4, 3); r(5, 4);
                        r(6, 4); r(7, 6); 8; 9; 10; 11; 12; 13;
                    ]]
                token_id = [[1:21;], [1:35;]]
                map((s,w,t)->map(NamedTuple{(:sentence_id, :word_id, :token_id)}, zip(s,w,t)), sentence_id, word_id, token_id)
            end
            @test nestedcall(getvalue, tkr(sentence)) == [
                [
                    "A"; s("single"); "s"; "en"; "t"; "en";
                    s("ce"); s("with"); "3"; "1"; s("char"); ".";
                ]
            ]
            @test nestedcall(getmeta, tkr(sentence)) == begin
                word_id = [
                    1; r(2, 6); 3; 4; 5; 6;
                    r(7, 2); r(8, 4); 9; 10; r(11, 4); 12;
                ]
                token_id = [1:24;]
                [map(NamedTuple{(:word_id, :token_id)}, zip(word_id, token_id))]
            end
            @test nestedcall(getvalue, tkr(word)) == ["w", "o", "r", "d"]
            @test nestedcall(getmeta, tkr(word)) == map(NamedTuple{(:word_id, :token_id)}, zip(r(1, 4), 1:4))
        end

        @testset "@stage" begin
            @test_throws Exception @macroexpand(TextEncodeBase.@stage SomeStage{A})
            @test_throws Exception @macroexpand(TextEncodeBase.@stage SomeStage{A, B, C})
            @test_throws Exception @macroexpand(TextEncodeBase.@stage SomeStage{A, B} <: C D)
            @test_throws Exception @macroexpand(TextEncodeBase.@stage 3)
            @test_throws Exception @macroexpand(TextEncodeBase.@stage SomeStage{A}())

            @test_nowarn @macroexpand(TextEncodeBase.@stage SomeStage)
            @test_nowarn @macroexpand(TextEncodeBase.@stage SomeStage{A, B} <: TokenStages)
        end

        @testset "batch" begin
            tkr = NestedTokenizer(IndexedTokenization())
            document = document.x
            sentence = sentence.x
            another_sentence = "This is another sentence"
            batch_sentence = [split_sentences(document); sentence; another_sentence]
            batch_document = [document, sentence, another_sentence]
            @test nestedcall(getvalue, tkr(Batch{Sentence}(batch_sentence))) == nltk_word_tokenize.(batch_sentence)
            @test nestedcall(getmeta, tkr(Batch{Sentence}(batch_sentence))) ==
                map(enumerate(nltk_word_tokenize.(batch_sentence))) do (i, v)
                    map(enumerate(v)) do (j, x)
                        (sentence_id = i, word_id = j, token_id = j)
                    end
                end
            @test nestedcall(getvalue, tkr(Batch{Document}(batch_document))) == map(batch_document) do doc
                nltk_word_tokenize.(split_sentences(doc))
            end
            @test nestedcall(getmeta, tkr(Batch{Document}(batch_document))) ==
                map(enumerate(split_sentences.(batch_document))) do (i, d)
                    map(enumerate(nltk_word_tokenize.(d))) do (j, s)
                        map(enumerate(s)) do (k, x)
                            (document_id = i, sentence_id = j, word_id = k, token_id = k)
                        end
                    end
                end
        end

        @testset "show" begin
            @test sprint(show, FlatTokenizer()) == "FlatTokenizer(default)"
            @test sprint(show, FlatTokenizer(WordTokenization(tokenize=poormans_tokenize))) == "FlatTokenizer(WordTokenization(split_sentences = WordTokenizers.split_sentences, tokenize = WordTokenizers.poormans_tokenize))"
            @test sprint(show, FlatTokenizer(IndexedTokenization())) == "FlatTokenizer(IndexedTokenization(default))"
            @test sprint(show, FlatTokenizer(MatchTokenization([r"\d", r"en"]))) == "FlatTokenizer(MatchTokenization(default, patterns = Regex[r\"\\d\", r\"en\"]))"
            @test sprint(show, FlatTokenizer(IndexedTokenization(MatchTokenization([r"\d", r"en"])))) == "FlatTokenizer(IndexedTokenization(MatchTokenization(default, patterns = Regex[r\"\\d\", r\"en\"])))"
            @test sprint(show, NestedTokenizer(IndexedTokenization())) == "NestedTokenizer(IndexedTokenization(default))"
            @test sprint(show, FlatTokenizer(IndexedTokenization(CharTk()))) == "FlatTokenizer(IndexedTokenization(CharTk))"
            @test sprint(show, NestedTokenizer(IndexedTokenization(MatchTokenization(CharTk(), [r"\d", r"en"])))) == "NestedTokenizer(IndexedTokenization(MatchTokenization(CharTk, patterns = Regex[r\"\\d\", r\"en\"])))"
        end
    end

    @testset "Vocabulary" begin
        vocab = Vocab(["a", "b", "c", "a", "b", "c"])
        vocab_unk = Vocab(["a", "b", "xxx"], "xxx")
        vocab_char = Vocab{Char}('a':'z', ' ')
        vocab_int = Vocab{Int}(11:20, 0)

        @test length(vocab) == 3
        @test vocab.list == ["a", "b", "c"]
        @test vocab.unki == 0
        @test length(vocab_unk) == 3
        @test vocab_unk.list == ["a", "b", "xxx"]
        @test vocab_unk.unki == 3
        @test vocab_int.list == collect(11:20)
        @test sprint(show, vocab_int) == "Vocab{$Int, SizedArray}(size = 10, unk = 0, unki = 0)"

        @testset "lookup" begin
            @test lookup(vocab, "a") == 1
            @test lookup(vocab, "b") == 2
            @test lookup(vocab, "c") == 3
            @test lookup(vocab, "d") == 0

            @test lookup(vocab_unk, "a") == 1
            @test lookup(vocab_unk, "b") == 2
            @test lookup(vocab_unk, "c") == 3
            @test lookup(vocab_unk, "d") == 3

            @test lookup(vocab, 1) == "a"
            @test lookup(vocab, 2) == "b"
            @test lookup(vocab, 3) == "c"
            @test lookup(vocab, 0) == "[UNK]"
            @test lookup(vocab, 1000000) == "[UNK]"

            @test lookup(vocab_unk, 1) == "a"
            @test lookup(vocab_unk, 2) == "b"
            @test lookup(vocab_unk, 3) == "xxx"
            @test lookup(vocab_unk, 100000) == "xxx"

            @test lookup(vocab, 1,2,3,4) == ("a", "b", "c", "[UNK]")
            @test lookup(vocab_unk, 1,2,3,4) == ("a", "b", "xxx", "xxx")

            @test lookup(vocab, [1, "a", 0, "A", "[UNK]"]) == ["a", 1, "[UNK]", 0, 0]
            @test lookup(vocab_unk, [1, "a", 0, "A", "[UNK]"]) == ["a", 1, "xxx", 3, 3]

            @test lookup(vocab, [(1, "a"), (a=0, b="A"), "[UNK]"]) == [("a", 1), (a="[UNK]", b=0), 0]
            @test lookup(vocab_char, ['a', (x='x',)], 26) == ([1, (x=24,)], 'z')
        end

        @testset "lookup int" begin
            @test lookup(vocab_int, 1,2,3) == (11,12,13)
            @test lookup(Int, vocab_int, 1,2,3) == (0,0,0)
            @test lookup(Int, vocab_int, 11,12,13) == (1,2,3)
        end

        @testset "onehot" begin
            @test lookup(OneHot, vocab, "a") == OneHot(3, 1)
            @test lookup(OneHot, vocab, "A") == OneHot(3, 0)
            @test lookup(OneHot, vocab, ("a", "A")) == (OneHot(3, 1), OneHot(3, 0))
            @test lookup(OneHot, vocab, (x="a", y="A")) == (x=OneHot(3, 1), y=OneHot(3, 0))
            @test lookup(OneHot, vocab, "a", "b", "c", "d") == OneHotArray(3, [1,2,3,0])
            @test lookup(OneHot, vocab, ["a", "b", "c", "d"]) == OneHotArray(3, [1,2,3,0])
            @test lookup(OneHot, vocab, ["a" "b"; "c" "d"]) == OneHotArray(3, [1 2; 3 0])
            @test lookup(OneHot, vocab, ["a", "b"], ["c", "d"]) == OneHotArray(3, [1,2,3,0])
            @test lookup(OneHot, vocab, ["a"], "b", ["c", "d"], "z") == OneHotArray(3, [1,2,3,0,0])

            @test_throws DomainError lookup(OneHot, vocab, 1)
            @test_throws DomainError lookup(OneHot, vocab, [1,2,3])
            @test_throws DomainError lookup(OneHot, vocab, "1",2,3)
            @test_throws DomainError lookup(OneHot, vocab, ["1",2,3])

            @test lookup(vocab, OneHot(3, 1)) == "a"
            @test lookup(vocab, OneHot(5, 1)) == "a"
            @test lookup(vocab, OneHot(5, 4)) == "[UNK]"
            @test lookup(vocab, OneHotArray(3, [1,2,3,0])) == ["a", "b", "c", "[UNK]"]
            @test lookup(vocab, OneHotArray(3, [1 2 3; 3 0 1])) == ["a" "b" "c"; "c" "[UNK]" "a"]
        end
    end

    @testset "Utils" begin
        @test_inferred nestedcall(x->x+1, [[[[3],[5, 6]]]])

        @testset "with_head_tail" begin
            x = collect(1:5)
            @test with_head_tail(x, 0, 6) == collect(0:6)
            @test with_head_tail(x, nothing, 6) == collect(1:6)
            @test with_head_tail(x, 0, nothing) == collect(0:5)
            @test with_head_tail(x, nothing, nothing) == x
            @test with_head_tail(0, 6)(x) == collect(0:6)
            @test with_head_tail(nothing, 6)(x) == collect(1:6)
            @test with_head_tail(0, nothing)(x) == collect(0:5)
            @test with_head_tail(nothing, nothing)(x) == x

            @test with_head_tail(x; head=0, tail=6) == collect(0:6)
            @test with_head_tail(x, tail=6) == collect(1:6)
            @test with_head_tail(x, head=0) == collect(0:5)
            @test with_head_tail(head=0, tail=6)(x) == collect(0:6)
            @test with_head_tail(tail=6)(x) == collect(1:6)
            @test with_head_tail(head=0)(x) == collect(0:5)

            @test with_head_tail(AbstractVector[[x], 1:5, 2:3], -1, -2) == [[[-1;x;-2]], [-1; 1:5; -2], [-1; 2:3; -2]]
            @test with_head_tail(Any[Any[x], 1:5, 2:3], -1, -2) == [[[-1;x;-2]], [-1; 1:5; -2], [-1; 2:3; -2]]
            @test with_head_tail(Any[Any[Any[0,1,2]]], 5, 5) == [[[5,0,1,2,5]]]
        end

        @testset "trunc_or_pad" begin
            x = collect(1:9)
            @test trunc_or_pad(x, 5, 0) == collect(1:5)
            @test trunc_or_pad(1:3, 5, 0) == [1:3; 0; 0]
            @test trunc_or_pad(x, nothing, 0) == collect(1:9)
            @test trunc_or_pad(1:3, nothing, 0) == collect(1:3)

            @test trunc_or_pad(5, 0)(x) == collect(1:5)
            @test trunc_or_pad(5, 0)(1:3) == [1:3; 0; 0]
            @test trunc_or_pad(nothing, 0)(x) == collect(1:9)
            @test trunc_or_pad(nothing, 0)(1:3) == collect(1:3)

            @test trunc_or_pad(AbstractVector[[x], 1:5, 2:3], 7, -1) == [[collect(1:7)], [1:5; -1; -1], [2:3; fill(-1, 5)]]
            @test trunc_or_pad(AbstractVector[[x], 1:5, 2:3], nothing, -1) == [[collect(1:9)], [1:5; fill(-1, 4)], [2:3; fill(-1, 7)]]
            @test trunc_or_pad(Any[Any[[0,0], 1:10], [1]], 7, -1) == [[[0; 0; fill(-1,5)], collect(1:7)], [1; fill(-1,6)]]
            @test trunc_or_pad(Any[Any[Any[0,1,2]]], 5, 0) == [[[0,1,2,0,0]]]
        end

        @testset "trunc_and_pad" begin
            x = collect(1:9)
            @test trunc_and_pad(x, 5, 0) == collect(1:5)
            @test trunc_and_pad(1:3, 5, 0) == [1:3;]
            @test trunc_and_pad(x, nothing, 0) == collect(1:9)
            @test trunc_and_pad(1:3, nothing, 0) == collect(1:3)

            @test trunc_and_pad(5, 0)(x) == collect(1:5)
            @test trunc_and_pad(5, 0)(1:3) == [1:3;]
            @test trunc_and_pad(nothing, 0)(x) == collect(1:9)
            @test trunc_and_pad(nothing, 0)(1:3) == collect(1:3)

            @test trunc_and_pad(AbstractVector[[x], 1:5, 2:3], 7, -1) == [[collect(1:7)], [1:5; -1; -1], [2:3; fill(-1, 5)]]
            @test trunc_and_pad(AbstractVector[[x], 1:5, 2:3], nothing, -1) == [[collect(1:9)], [1:5; fill(-1, 4)], [2:3; fill(-1, 7)]]
            @test trunc_and_pad(Any[Any[[0,0], 1:10], [1]], 7, -1) == [[[0; 0; fill(-1,5)], collect(1:7)], [1; fill(-1,6)]]
            @test trunc_and_pad(Any[Any[Any[0,1,2]]], 5, 0) == [[[0,1,2]]]
        end

        @testset "nested2batch" begin
            x = randn(5,4,3,2)
            x_slices = [x[i:i+5-1] for i in 1:5:length(x)]
            y = [[[x_slices[1],x_slices[2],x_slices[3],x_slices[4]],
                  [x_slices[5],x_slices[6],x_slices[7],x_slices[8]],
                  [x_slices[9],x_slices[10],x_slices[11],x_slices[12]],],
                 [[x_slices[13],x_slices[14],x_slices[15],x_slices[16]],
                  [x_slices[17],x_slices[18],x_slices[19],x_slices[20]],
                  [x_slices[21],x_slices[22],x_slices[23],x_slices[24]],]]
            y2 = [[cat(x_slices[1],x_slices[2],x_slices[3],x_slices[4], dims=2),
                   cat(x_slices[5],x_slices[6],x_slices[7],x_slices[8], dims=2),
                   cat(x_slices[9],x_slices[10],x_slices[11],x_slices[12], dims=2),],
                  [cat(x_slices[13],x_slices[14],x_slices[15],x_slices[16], dims=2),
                   cat(x_slices[17],x_slices[18],x_slices[19],x_slices[20], dims=2),
                   cat(x_slices[21],x_slices[22],x_slices[23],x_slices[24], dims=2),]]
            y3 = [cat(cat(x_slices[1],x_slices[2],x_slices[3],x_slices[4], dims=2),
                      cat(x_slices[5],x_slices[6],x_slices[7],x_slices[8], dims=2),
                      cat(x_slices[9],x_slices[10],x_slices[11],x_slices[12], dims=2), dims=3),
                  cat(cat(x_slices[13],x_slices[14],x_slices[15],x_slices[16], dims=2),
                      cat(x_slices[17],x_slices[18],x_slices[19],x_slices[20], dims=2),
                      cat(x_slices[21],x_slices[22],x_slices[23],x_slices[24], dims=2), dims=3)]
            y4 = Any[Any[Any[x_slices[1],x_slices[2],x_slices[3],x_slices[4]],
                  Any[x_slices[5],x_slices[6],x_slices[7],x_slices[8]],
                  Any[x_slices[9],x_slices[10],x_slices[11],x_slices[12]],],
                 Any[Any[x_slices[13],x_slices[14],x_slices[15],x_slices[16]],
                  Any[x_slices[17],x_slices[18],x_slices[19],x_slices[20]],
                  Any[x_slices[21],x_slices[22],x_slices[23],x_slices[24]],]]
            x_slices_any = [Array{Any}(x[i:i+5-1]) for i in 1:5:length(x)]
            y5 = [[[x_slices_any[1],x_slices_any[2],x_slices_any[3],x_slices_any[4]],
                  [x_slices_any[5],x_slices_any[6],x_slices_any[7],x_slices_any[8]],
                  [x_slices_any[9],x_slices_any[10],x_slices_any[11],x_slices_any[12]],],
                 [[x_slices_any[13],x_slices_any[14],x_slices_any[15],x_slices_any[16]],
                  [x_slices_any[17],x_slices_any[18],x_slices_any[19],x_slices_any[20]],
                  [x_slices_any[21],x_slices_any[22],x_slices_any[23],x_slices_any[24]],]]

            @test nested2batch(y) == x
            @test nested2batch(y2) == x
            @test nested2batch(y3) == x
            @test nested2batch(y4) == x
            @test nested2batch(y5) == x

            @test_throws DimensionMismatch nested2batch([[1:5], 2:6])
        end
    end

    @testset "Encoder" begin
        sentence = Sentence("A single sentence with 31 char.")
        tkr = NestedTokenizer(IndexedTokenization(CharTk()))
        vocab = Vocab(map(string, ['a':'z'; 'A':'Z']))
        enc = TextEncoder(tkr, vocab, nested2batch∘nestedcall(getvalue))
        s(x) = mapfoldl(y->split(y,""), append!, split(x); init=String[])
        @test encode(enc, sentence) == reshape(lookup(OneHot, vocab, s(sentence.x)), Val(3))
        @test decode(enc, encode(enc, sentence)) == lookup(vocab, reshape(lookup(OneHot, vocab, s(sentence.x)), Val(3)))

        enc2 = TextEncoder(tkr, vocab) do e
            nested2batch∘TextEncodeBase.process(e)
        end
        @test enc == enc2
    end
end


text:
# Initializing Python (surprisingly complicated; see also deps/build.jl)

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

# global PyObject constants, initialized to NULL and then overwritten in __init__
# (eventually, the ability to define global const in __init__ may go away,
#  and in any case this is better for type inference during precompilation)
const inspect = PyNULL()
const builtin = PyNULL()
const BuiltinFunctionType = PyNULL()
const TypeType = PyNULL()
const MethodType = PyNULL()
const MethodWrapperType = PyNULL()
const ufuncType = PyNULL()
const format_traceback = PyNULL()
const pyproperty = PyNULL()
const jlfun2pyfun = PyNULL()
const c_void_p_Type = PyNULL()

# other global constants initialized at runtime are defined via Ref
# or are simply left as non-const values
const pynothing = Ref{PyPtr}(0)
const pyxrange = Ref{PyPtr}(0)

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

function __init__()
    # issue #189
    libpy_handle = libpython === nothing ? C_NULL :
        Libdl.dlopen(libpython, Libdl.RTLD_LAZY|Libdl.RTLD_DEEPBIND|Libdl.RTLD_GLOBAL)

    already_inited = 0 != ccall((@pysym :Py_IsInitialized), Cint, ())

    if !already_inited
        Py_SetPythonHome(libpy_handle, PYTHONHOME, wPYTHONHOME, pyversion)
        if !isempty(pyprogramname)
            if pyversion.major < 3
                ccall((@pysym :Py_SetProgramName), Void, (Cstring,), pyprogramname)
            else
                ccall((@pysym :Py_SetProgramName), Void, (Ptr{Cwchar_t},), wpyprogramname)
            end
        end
        ccall((@pysym :Py_InitializeEx), Void, (Cint,), 0)
    end

    # Will get reinitialized properly on first use
    is_windows() && (PyActCtx[] = C_NULL)

    # Make sure python wasn't upgraded underneath us
    new_pyversion = convert(VersionNumber, split(unsafe_string(ccall(@pysym(:Py_GetVersion),
                               Ptr{UInt8}, ())))[1])

    if new_pyversion.major != pyversion.major
        error("PyCall precompiled with Python $pyversion, but now using Python $new_pyversion; ",
              "you need to relaunch Julia and re-run Pkg.build(\"PyCall\")")
    end

    copy!(inspect, pyimport("inspect"))
    copy!(builtin, pyimport(pyversion.major < 3 ? "__builtin__" : "builtins"))
    copy!(pyproperty, pybuiltin(:property))

    pyexc_initialize() # mappings from Julia Exception types to Python exceptions

    types = pyimport("types")
    copy!(TypeType, pybuiltin("type")) # for pytypeof

    # cache Python None -- PyPtr, not PyObject, to prevent it from
    # being finalized prematurely on exit
    pynothing[] = @pyglobalobj(:_Py_NoneStruct)

    # xrange type (or range in Python 3)
    pyxrange[] = @pyglobalobj(:PyRange_Type)

    # ctypes.c_void_p for Ptr types
    copy!(c_void_p_Type, pyimport("ctypes")["c_void_p"])

    # traceback.format_tb function, for show(PyError)
    copy!(format_traceback, pyimport("traceback")["format_tb"])

    init_datetime()
    pyjlwrap_init()

    # jl_FunctionType is a class, and when assigning it to an object
    #    obj[:foo] = some_julia_function
    # it won't behave like a regular Python method because it's not a Python
    # function (in particular, `self` won't be passed to it). The solution is:
    #    obj[:foo] = jlfun2pyfun(some_julia_function)
    # This is a bit of a kludge, obviously.
    copy!(jlfun2pyfun,
          pyeval_("""lambda f: lambda *args, **kwargs: f(*args, **kwargs)"""))

    if !already_inited
        # some modules (e.g. IPython) expect sys.argv to be set
        if pyversion.major < 3
            argv_s = Compat.String("")
            argv = unsafe_convert(Ptr{UInt8}, argv_s)
            ccall(@pysym(:PySys_SetArgvEx), Void, (Cint,Ptr{Ptr{UInt8}},Cint), 1, &argv, 0)
        else
            argv_s = Cwchar_t[0]
            argv   = unsafe_convert(Ptr{Cwchar_t}, argv_s)
            ccall(@pysym(:PySys_SetArgvEx), Void, (Cint, Ptr{Ptr{Cwchar_t}}, Cint), 1, &argv, 0)
        end

        # Some Python code checks sys.ps1 to see if it is running
        # interactively, and refuses to be interactive otherwise.
        # (e.g. Matplotlib: see PyPlot#79)
        if isinteractive()
            let sys = pyimport("sys")
                if !haskey(sys, "ps1")
                    sys["ps1"] = ">>> "
                end
            end
        end
    end
end


text:
#
# This file is part of the Actors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

"""
```
become!(lk::Link, func, args1...; kwargs...)
become!(name::Symbol, ....)
```
Cause an actor to change behavior.

# Arguments
- actor `lk::Link` (or `name::Symbol` if registered),
- `func`: a callable object,
- `args1...`: (partial) arguments to `func`,
- `kwargs...`: keyword arguments to `func`.
"""
function become!(lk::Link, func, args...; kwargs...)
    isempty(args) && isempty(kwargs) ?
        send(lk, Become(func)) :
        send(lk, Become(Bhv(func, args...; kwargs...)))
end
become!(name::Symbol, args...; kwargs...) = become!(whereis(name), args...; kwargs...)

"""
```
call(lk::Link, [from::Link,] args2...; timeout::Real=5.0)
call(name::Symbol, ....)
```
Call an actor to execute its behavior and to send a
[`Response`](@ref) with the result.

# Arguments
- actor `lk::Link` (or `name::Symbol` if registered),
- `from::Link`: sender link,
- `args2...`: remaining arguments to the actor.
- `timeout::Real=5.0`: timeout in seconds.

**Note:** If `from` is omitted, `call` blocks and returns the result
"""
call(lk::Link, from::Link, args...) = send(lk, Call(args, from))
call(lk::Link, args...; timeout::Real=5.0) = request(lk, Call, args...; timeout)
call(name::Symbol, args...; kwargs...) = call(whereis(name), args...; kwargs...)

"""
```
cast(lk::Link, args2...)
cast(name::Symbol, args2...)
```
Cast `args2...` to the actor `lk` (or `name` if registered)
to execute its behavior with `args2...` without sending a
response.

**Note:** you can prompt the returned value with [`query`](@ref).
"""
cast(lk::Link, args...) = send(lk, Cast(args))
cast(name::Symbol, args...) = cast(whereis(name), args...)

"""
```
exec(lk::Link, from::Link, f, args...; kwargs...)
exec(lk::Link, f, args...; timeout::Real=5.0, kwargs...)
exec(name::Symbol, ....)
```

Ask an actor `lk` (or `name` if registered) to execute an
arbitrary function and to send the returned value as
[`Response`](@ref).

# Arguments
- actor `lk::Link` or `name::Symbol` if registered,
- `from::Link`: the link a `Response` should be sent to.
- `f`: a callable object,
- `args...; kwargs...`: arguments and keyword arguments to it,
- `timeout::Real=5.0`: timeout in seconds. Set `timeout=Inf`
    if you don't want to timeout.

**Note:** If `from` is ommitted, `exec` blocks, waits and
returns the result (with a `timeout`).
"""
function exec(lk::Link, from::Link, f, args...; kwargs...)
    isempty(args) && isempty(kwargs) ?
        send(lk, Exec(f, from)) :
        send(lk, Exec(Bhv(f, args...; kwargs...), from))
end
function exec(lk::Link, f, args...; timeout::Real=5.0, kwargs...)
    isempty(args) && isempty(kwargs) ?
        request(lk, Exec, f; timeout=timeout) :
        request(lk, Exec, Bhv(f, args...; kwargs...); timeout)
end
exec(name::Symbol, args...; kwargs...) = exec(whereis(name), args...; kwargs...)

"""
```
exit!(lk::Link, reason=:normal)
exit!(name::Symbol, ....)
```
Tell an actor `lk` (or `name` if registered) to stop. If it
has a [`term`](@ref _ACT) function, it calls that with
`reason` as last argument.
"""
exit!(lk::Link, reason=:normal) = send(lk, Exit(reason, fill(nothing, 3)...))
exit!(name::Symbol, reason=:normal) = exit!(whereis(name), reason)

"""
```
init!(lk::Link, f, args...; kwargs...)
init!(name::Symbol, ....)
```
Tell an actor `lk` to save the callable object `f` with the given
arguments as an `init` object in its [`_ACT`](@ref) variable.
The `init` object will be called by a supervisor at actor restart.

# Arguments
- actor `lk::Link` or `name::Symbol` if registered,
- `f`: callable object,
- `args...`: arguments to `f`,
- `kwargs...`: keyword arguments to `f`.
"""
init!(lk::Link, f, args...; kwargs...) =
    send(lk, Init(Bhv(f, args...; kwargs...)))
init!(name::Symbol, args...; kwargs...) = init!(whereis(name), args...; kwargs...)

"""
```
query(lk::Link, [from::Link,] s::Symbol; timeout::Real=5.0)
query(name::Symbol, ....)
```

Query an actor about an internal state variable `s`.

# Parameters
- actor `lk::Link` or `name::Symbol` if registered,
- `from::Link`: sender link,
- `s::Symbol` one of `:mode`,`:bhv`,`:res`,`:sta`,`:usr`.
- `timeout::Real=5.0`:

**Note:** If `from` is omitted, `query` blocks and returns
the response. In that case there is a `timeout`.

# Examples

```julia
julia> f(x, y; u=0, v=0) = x+y+u+v  # implement a behavior
f (generic function with 1 method)

julia> fact = spawn(Bhv(f, 1))     # start an actor with it
Link{Channel{Any}}(Channel{Any}(sz_max:32,sz_curr:0), 1, :default)

julia> query(fact, :mode)           # query the mode
:default

julia> cast(fact, 1)                # cast a 2nd argument to it
Actors.Cast((1,))

julia> query(fact, :res)            # query the result
2

julia> query(fact, :sta)            # query the state

julia> query(fact, :bhv)            # query the behavior
Bhv(f, (1,), Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}(), Actors.var"#2#4"{Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},typeof(f),Tuple{Int64}}(Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}(), f, (1,)))
```
"""
query(lk::Link, from::Link, s::Symbol=:sta) = send(lk, Query(s, from))
query(lk::Link, s::Symbol=:sta; timeout::Real=5.0) = request(lk, Query, s; timeout)
query(name::Symbol, args...; kwargs...) = query(whereis(name), args...; kwargs...)

"""
```
term!(lk::Link, f, args1...; kwargs...)
term!(name::Symbol, ....)
```
Tell an actor `lk` (or `name::Symbol` if registered) to
execute `f` with the given partial arguments and an
exit reason when it terminates.

The exit reason is added by the actor to `args1...` when it
exits.
"""
term!(lk::Link, f, args...; kwargs...) =
    send(lk, Term(Bhv(f, args...; kwargs...)))
term!(name::Symbol, args...; kwargs...) = term!(whereis(name), args...; kwargs...)

"""
    trapExit(lk::Link=self(), mode=:sticky)

Change the mode of an actor.

A `:sticky` actor does not exit if it receives an [`Exit`](@ref)
signal from a connected actor and does not propagate it further.
Instead it reports the failure and saves a link to the failed actor.

See [`diag`](@ref) for getting links to failed actors
from a `:sticky` actor.
"""
trapExit(lk::Link=self(), mode=:sticky) = send(lk, Update(:mode, mode))

"""
```
update!(lk::Link, x; s::Symbol=:sta)
update!(lk::Link, arg::Args)
update!(name::Symbol, ....)
```
Update an actor's internal state `s` with `args...`.

# Arguments
- actor `lk::Link` or `name::Symbol` if registered,
- `x`: value/variable to update the choosen state with,
- `arg::Args`: arguments to update,
- `s::Symbol`: one of `:arg`, `:mode`, `:name`, `:self`, `:sta`, `:usr`.

*Note:* If you want to update the stored arguments to the
behavior function with `s=:arg`, you must pass an [`Args`](@ref)
to `arg`. If `Args` has keyword arguments, they are merged
with existing keyword arguments to the behavior function.

# Example
```julia
julia> update!(fact, 5)       # update the state variable
Actors.Update(:sta, 5)

julia> query(fact, :sta)      # query it
5

julia> update!(fact, Args(0, u=5, v=5));  # update arguments to the behavior

julia> call(fact, 0)          # call the actor with 0
10
```
"""
update!(lk::Link, x; s::Symbol=:sta) = send(lk, Update(s, x))
update!(lk::Link, arg::Args) = send(lk, Update(:arg, arg))
update!(name::Symbol, args...; kwargs...) = update!(whereis(name), args...; kwargs...)


text:
##################################################

"""
     SymFunction

Thin wrapper for symbolic functions that allows prime notation in place of using `diff`.

Functions constructed through `SymFunction("f")` or `@symfuns f`.
"""
mutable struct SymFunction <: SymbolicObject
    __pyobject__::PyCall.PyObject
    n::Int
end

Base.:(==)(x::SymFunction, y::SymFunction) = x.__pyobject__ == y.__pyobject__ && x.n == y.n
Base.hash(x::SymFunction) = hash((hash(x.__pyobject__),x.n))

# these are from https://github.com/OptMist-Tokyo/DAEPreprocessor.jl/blob/sympy_warning/src/symbolic.jl
derivative(x::SymFunction, d::Int = 1) = SymFunction(x.__pyobject__, x.n + d)

Base.show(io::IO, u::SymFunction) = print(io, "$(string(Sym(u.__pyobject__)))" * repeat("'", u.n))
Base.show(io::IO, ::MIME"text/plain", u::SymFunction) = print(io, "$(string(Sym(u.__pyobject__)))" * repeat("'", u.n))
Base.show(io::IO, ::MIME"text/latex", x::SymFunction) = print(io, "\\begin{align*}" * latex(x) * "\\end{align*}")
function Base.show(io::IO, ::MIME"text/latex", x::AbstractArray{SymFunction, 1})
    print(io, "\\begin{align*}\\left[\\begin{array}{c}" * join(latex.(x), "\\\\") * "\\end{array}\\right]\\end{align*}")
end

latex(x::SymFunction) = latex(Sym(x.__pyobject__)) * repeat("'", x.n)

"""

Create a symbolic function. These can be used for specifying differential equations.
For these objects we can specify derivatives with the transpose
operator (e.g., `u''`) as opposed to, say `diff(u(x), x, 2)`.

Example:
```
u = SymFunction("u")
u'
```

Alternatively, we can pass a comma separated string of variable names to create
more than one at a time. (The `cls=symfunction` is no longer supported):

```
F,G,H = SymFunction("F, G, H")
```

"""
function SymFunction(x::T; kwargs...) where {T<:AbstractString}
    us = split(x, r",\s*")
    if length(us) > 1
        map(u -> SymFunction(sympy."Function"(u; kwargs...), 0), us)
    else
        SymFunction(sympy."Function"(x; kwargs...), 0)
    end
end

"""
    symfuns...
Thanks to `@alhirzel` for the contribution.
"""
macro symfuns(x...)
    q = Expr(:block)
    as = []    # running list of assumptions to be applied
    fs = []    # running list of symbols created
    for s in reverse(x)
        if isa(s, Expr)    # either an assumption or a named variable
            if s.head == :(=)
                s.head = :kw
                push!(as, s)
            elseif s.head == :(=>)
                push!(fs, s.args[1])
                push!(q.args, Expr(:(=), s.args[1], Expr(:call, :SymFunction, s.args[2], map(esc,as)...)))
            end
        elseif isa(s, Symbol)   # raw symbol to be created
            push!(fs, s)
            # @show s
            push!(q.args, Expr(:(=), esc(s), Expr(:call, :SymFunction, string(s), map(esc,as)...)))
        else
            throw(AssertionError("@symfuns expected a list of symbols and assumptions"))
        end
    end
    push!(q.args, Expr(:tuple, map(esc,reverse(fs))...)) # return all of the symbols we created
    q
end

function PyCall.PyObject(f::SymFunction)
    f.n == 0 && return f.__pyobject__
    z = symbols(gensym())
    f(z).__pyobject__  #  diff(f(__x__), __x__, f.n).__pyobject__
end

Base.adjoint(x::SymFunction) = SymFunction(x.__pyobject__, x.n + 1)

(u::SymFunction)(x::Base.Dict) = throw(ArgumentError("IVPsolutions can only be called with symbolic objects"))
(u::SymFunction)(x::Base.Pair) = throw(ArgumentError("IVPsolutions can only be called with symbolic objects"))
function (u::SymFunction)(x)
    if u.n == 0
        PyObject(u)(x) #u.__pyobject__(PyObject(x))
    else
        __x = Sym("__x")
        diff(u.__pyobject__(PyObject(__x)), __x, u.n)(__x => x)
    end
end

(u::SymFunction)(x, y...) = u.n== 0 ? u.__pyobject__(map(PyObject, vcat(x, y...))...) : error("Need to implement derivatives of symbolic functions of two or more variables")