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

#=
Example usage:
```
	./examples/plot_random_from_file.jl /Users/jakeireland/Desktop/random_number_in_code,q=21,n=4,d=3,i=1000,m=10000.csv 6 7 3 # dataset, q, n, d
```
=#

include("plot_random_core.jl")

# reading from pre-obtained data
if isempty(ARGS) || length(ARGS) < 5
	throw(error("""
	Please specify the dataset from the command line, along with q, n, and d.  For example
	```
		./examples/plot_random_from_file.jl /Users/jakeireland/Desktop/random_number_in_code,q=21,n=4,d=3,i=1000,m=10000.csv 6 7 3 10000# dataset, q, n, d, selection_set_size
	```
	That is respectively: the size of the alphabet; the block length of the words; the distance of the code; the number of datapoints in the code; and the selection set size of each "random" sample.
	"""))
end

global data = ARGS[1] # e.g. "/Users/jakeireland/Desktop/random_number_in_code,q=21,n=4,d=3,i=1000,m=10000.csv"
global q = parse(Int, ARGS[2])
global n = parse(Int, ARGS[3])
global d = parse(Int, ARGS[4])
global m = parse(Int, ARGS[5])

global upper_bound_adjustment = 10
global 🐖 = length(get_codewords_greedy(q, n, d))
global 🍖 = hamming_bound(q, n, d)
global 🐺 = singleton_bound(q, n, d)
# data = CSV.read("/Users/jakeireland/projects/CodingTheory.jl/other/random_number_analysis/random_number_in_code,q=$(q),n=$(n),d=$(d),i=$(stop_at),m=$(m).csv")
# data = CSV.read("/Users/jakeireland/projects/CodingTheory.jl/other/random_number_analysis/mildly_interesting_ones/[2]-10-3-2/random_number_in_code,q=10,n=3,d=2,i=1000,m=10000.csv")
data = CSV.read(data)
graphing(q, n, d, data, nrow(data), m = m)


text:
using MD5
puzzle_input = "yzbqklnj"
println("Computing AdventCoins for all the economically forward-thinking little girls and boys")
println("Starting puzzle input is $puzzle_input")
for i=1:1:5000000000
   md5_input = string(puzzle_input,i)
   md5_output = bytes2hex(md5(md5_input))
   first_6_digits = md5_output[1:6]
   if first_6_digits == "000000"
        println("The lowest positive number that results in an MD5 hash with six leading zeroes is $i")
        println("Puzzle input + lowest positive number is $md5_input")
        println("MD5 hash: $md5_output")
        break
   end
end

ext: jl


text:
using SoftPosit
using Test

include("test_conversions.jl")
include("test_arithmetics.jl")
include("test_rounding.jl")
include("test_comparison.jl")
include("test_nextprevfloat.jl")
include("test_bitstring.jl")
include("test_sign.jl")
include("test_eps.jl")
include("test_quire.jl")
include("test_explog_trigono.jl")
include("test_promotion.jl")


text:
module GalaxyBrain

include("lex.jl")
include("parse.jl")

export @bf_str

macro bf_str(program::String)
    instructions = filter((lex ∘ Vector{Char})(program)) do instruction
        return instruction isa AbstractBFInstruction
    end
    brainfuck_function = quote
        function (; input::Union{IO,AbstractString}="", output::IO=stdout,
                  index::Int=1, memory_size::Int=30000)
            memory = repeat([zero(UInt8)], memory_size)
            characters = input isa AbstractString ? IOBuffer(input) : input
        end
    end
    body = last(last(brainfuck_function.args).args).args
    append!(body, parse_instructions(Vector{AbstractBFInstruction}(instructions)))
    push!(body, :(return nothing))
    return esc(Base.remove_linenums!(brainfuck_function))
end

end


text:
function DataSplit(X,Y;p=[0.6,0.2,0.2])
   m=size(X,1)
   TrNumber=round(Int,m*p[1])
   CvNumber=round(Int,m*p[2])
   TsNumber=round(Int,m*p[3])

    #Ramdomly sampled the X inputs
    r = sample(1:m, m, replace = false)
    Xr=X[r,:]
    Yr=Y[r,:]

    #Get X and Y train
    Xtr=Xr[1:TrNumber,:]
    Ytr=Yr[1:TrNumber]

    #Get X and Y CrossValidation
    Xcv=Xr[TrNumber+1:TrNumber+CvNumber,:]
    Ycv=Yr[TrNumber+1:TrNumber+CvNumber]

    #Get X and Y Test
    Xts=Xr[TrNumber+CvNumber+1:end,:]
    Yts=Yr[TrNumber+CvNumber+1:end]

    return Xtr, Ytr, Xcv, Ycv, Xts, Yts

end


text:
function hunt(;
  backend=nothing,
  np=nothing,
  parts=nothing,
  title = "hunt",
  nruns=1,
  path=".",
  kwargs...)

  for ir in 1:nruns
    _title = title*"_r$ir"
    @assert parts === nothing
    if backend === nothing
      @assert np === nothing
      info, t = _hunt(;title=_title,path=path,kwargs...)
    else
      @assert backend !== nothing
      info, t = prun(_find_backend(backend),(np...,1)) do _parts
        _hunt(;parts=_parts,title=_title,path=path,kwargs...)
      end
    end
    info[:np] = np
    info[:backend] = backend
    info[:title] = title
    map_main(t.data) do data
      for (k,v) in data
        info[Symbol("time_$k")] = v.max
      end
      save(joinpath(path,"$_title.bson"),info)
    end
  end

  nothing
end

function _find_backend(s)
  if s === :sequential
    backend = sequential
  elseif s === :mpi
    backend = mpi
  else
    error()
  end
  backend
end

function _hunt(;
  parts=nothing,
  nc=(3,3),
  ν=1.0,
  ρ=1.0,
  σ=1.0,
  B=(0.0,1.0,0.0),
  f=(0.0,0.0,1.0),
  L=1.0,
  u0=1.0,
  B0=norm(VectorValue(B)),
  nsums = 10,
  vtk=true,
  title="test",
  path=".",
  debug=false,
  solver=:julia,
  verbose=true,
  kmap=1,
  petsc_options="-snes_monitor -ksp_error_if_not_converged true -ksp_converged_reason -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps"
  )

  info = Dict{Symbol,Any}()

  if parts === nothing
    t_parts = get_part_ids(sequential,1)
  else
    t_parts = parts
  end
  t = PTimer(t_parts,verbose=verbose)
  tic!(t,barrier=true)

  domain_phys = (-L,L,-L,L,0.0*L,0.1*L)

  # Reduced quantities
  Re = u0*L/ν
  Ha = B0*L*sqrt(σ/(ρ*ν))
  N = Ha^2/Re
  f̄ = (L/(ρ*u0^2))*VectorValue(f)
  B̄ = (1/B0)*VectorValue(B)
  α = 1.0
  β = 1.0/Re
  γ = N
  domain = domain_phys ./ L

  # Prepare problem in terms of reduced quantities
  layer(x,a) = sign(x)*abs(x)^(1/a)
  map((x,y,z)) = VectorValue(layer(x,kmap),y,z)
  partition=(nc[1],nc[2],3)
  model = CartesianDiscreteModel(
    parts,domain,partition;isperiodic=(false,false,true),map=map)
  Ω = Interior(model)
  labels = get_face_labeling(model)
  tags_u = append!(collect(1:20),[23,24,25,26])
  tags_j = append!(collect(1:20),[25,26])
  add_tag_from_tags!(labels,"noslip",tags_u)
  add_tag_from_tags!(labels,"insulating",tags_j)

  params = Dict(
    :ptimer=>t,
    :debug=>debug,
    :model=>model,
    :fluid=>Dict(
      :domain=>model,
      :α=>α,
      :β=>β,
      :γ=>γ,
      :f=>f̄,
      :B=>B̄,
    ),
    :bcs => Dict(
      :u=>Dict(:tags=>"noslip"),
      :j=>Dict(:tags=>"insulating"),
    )
  )

  toc!(t,"pre_process")

  # Solve it
  if solver == :julia
    params[:solver] = NLSolver(show_trace=true,method=:newton)
    xh,fullparams = main(params)
  elseif solver == :petsc
    xh,fullparams = GridapPETSc.with(args=split(petsc_options)) do
      params[:matrix_type] = SparseMatrixCSR{0,PetscScalar,PetscInt}
      params[:vector_type] = Vector{PetscScalar}
      params[:solver] = PETScNonlinearSolver()
      params[:solver_postpro] = cache -> snes_postpro(cache,info)
      main(params)
    end
  else
    error()
  end
  t = fullparams[:ptimer]

  # Rescale quantities

  tic!(t,barrier=true)
  ūh,p̄h,j̄h,φ̄h = xh
  uh = u0*ūh
  ph = (ρ*u0^2)*p̄h
  jh = (σ*u0*B0)*j̄h
  φh = (u0*B0*L)*φ̄h

  if L == 1.0
    Ω_phys = Ω
  else
    Ω_phys = _warp(model,Ω,L)
  end

  # Post process

  μ = ρ*ν
  grad_pz = -f[3]/ρ
  u(x) = analytical_hunt_u(L,L,μ,grad_pz,Ha,nsums,x)
  j(x) = analytical_hunt_j(L,L,σ,μ,grad_pz,Ha,nsums,x)
  u_ref(x) = analytical_hunt_u(L,L,μ,grad_pz,Ha,2*nsums,x)
  j_ref(x) = analytical_hunt_j(L,L,σ,μ,grad_pz,Ha,2*nsums,x)

  k = 2
  dΩ_phys = Measure(Ω_phys,2*(k+1))
  eu = u - uh
  ej = j - jh
  eu_h1 = sqrt(sum(∫( ∇(eu)⊙∇(eu) + eu⋅eu  )dΩ_phys))
  eu_l2 = sqrt(sum(∫( eu⋅eu )dΩ_phys))
  ej_l2 = sqrt(sum(∫( ej⋅ej )dΩ_phys))
  eu_ref = u_ref - uh
  ej_ref = j_ref - jh
  eu_ref_h1 = sqrt(sum(∫( ∇(eu_ref)⊙∇(eu_ref) + eu_ref⋅eu_ref  )dΩ_phys))
  eu_ref_l2 = sqrt(sum(∫( eu_ref⋅eu_ref )dΩ_phys))
  ej_ref_l2 = sqrt(sum(∫( ej_ref⋅ej_ref )dΩ_phys))
  uh_l2 = sqrt(sum(∫( uh⋅uh )dΩ_phys))
  uh_h1 = sqrt(sum(∫( ∇(uh)⊙∇(uh) + uh⋅uh )dΩ_phys))
  jh_l2 = sqrt(sum(∫( jh⋅jh )dΩ_phys))
  toc!(t,"post_process")

  if vtk
    writevtk(Ω_phys,joinpath(path,title),
      order=2,
      cellfields=[
        "uh"=>uh,"ph"=>ph,"jh"=>jh,"φh"=>φh,
        "u"=>u,"j"=>j,"u_ref"=>u_ref,"j_ref"=>j_ref])
    toc!(t,"vtk")
  end
  if verbose
    display(t)
  end

  info[:nc] = nc
  info[:nsums] = nsums
  info[:ncells] = num_cells(model)
  info[:ndofs_u] = length(get_free_dof_values(ūh))
  info[:ndofs_p] = length(get_free_dof_values(p̄h))
  info[:ndofs_j] = length(get_free_dof_values(j̄h))
  info[:ndofs_φ] = length(get_free_dof_values(φ̄h))
  info[:ndofs] = length(get_free_dof_values(xh))
  info[:Re] = Re
  info[:Ha] = Ha
  info[:eu_h1] = eu_h1
  info[:eu_l2] = eu_l2
  info[:ej_l2] = ej_l2
  info[:eu_ref_h1] = eu_ref_h1
  info[:eu_ref_l2] = eu_ref_l2
  info[:ej_ref_l2] = ej_ref_l2
  info[:uh_h1] = uh_h1
  info[:uh_l2] = uh_l2
  info[:jh_l2] = jh_l2
  info[:kmap] = kmap

  info, t
end

# This is not very elegant. This needs to be solved by Gridap and GridapDistributed
function _warp(model::DiscreteModel,Ω::Triangulation,L)
  grid_phys = UnstructuredGrid(get_grid(model))
  node_coords = get_node_coordinates(grid_phys)
  node_coords .= L .* node_coords
  Ω_phys = Gridap.Geometry.BodyFittedTriangulation(Ω.model,grid_phys,Ω.tface_to_mface)
end

function _warp(
  model::GridapDistributed.DistributedDiscreteModel,
  Ω::GridapDistributed.DistributedTriangulation,L)
  trians = map_parts(model.models,Ω.trians) do model,Ω
    grid_phys = UnstructuredGrid(get_grid(model))
    node_coords = get_node_coordinates(grid_phys)
    node_coords .= L .* node_coords
    gp = GridPortion(grid_phys,Ω.tface_to_mface)
    Ω_phys = Gridap.Geometry.BodyFittedTriangulation(model,gp,Ω.tface_to_mface)
  end
  GridapDistributed.DistributedTriangulation(trians,model)
end

function analytical_hunt_u(
  a::Float64,       # semi-length of side walls
  b::Float64,       # semi-length of Hartmann walls
  μ::Float64,       # fluid viscosity
  grad_pz::Float64, # presure gradient
  Ha::Float64,      # Hartmann number
  n::Int,           # number of sumands included in Fourier series
  x)                # evaluation point

  l = b/a
  ξ = x[1]/a
  η = x[2]/a

  V = 0.0; V0=0.0;
  for k in 0:n
    α_k = (k + 0.5)*π/l
    N = (Ha^2 + 4*α_k^2)^(0.5)
    r1_k = 0.5*( Ha + N)
    r2_k = 0.5*(-Ha + N)

    num = exp(-r1_k*(1-η))+exp(-r1_k*(1+η))
    den = 1+exp(-2*r1_k)
    V2 = (r2_k/N)*(num/den)

    num = exp(-r2_k*(1-η))+exp(-r2_k*(1+η))
    den = 1+exp(-2*r2_k)
    V3 = (r1_k/N)*(num/den)

    V += 2*(-1)^k*cos(α_k * ξ)/(l*α_k^3) * (1-V2-V3)
  end
  u_z = V/μ * (-grad_pz) * a^2

  VectorValue(0.0*u_z,0.0*u_z,u_z)
end

function analytical_hunt_j(
  a::Float64,       # semi-length of side walls
  b::Float64,       # semi-length of Hartmann walls
  σ::Float64,       # fluid conductivity
  μ::Float64,       # fluid viscosity
  grad_pz::Float64, # presure gradient
  Ha::Float64,      # Hartmann number
  n::Int,           # number of sumands included in Fourier series
  x)                # evaluation point

  l = b/a
  ξ = x[1]/a
  η = x[2]/a

  H_dx = 0.0; H_dy = 0.0
  for k in 0:n
    α_k = (k + 0.5)*π/l
    N = sqrt(Ha^2 + 4*α_k^2)
    r1_k = 0.5*( Ha + N)
    r2_k = 0.5*(-Ha + N)

    num = exp(-r1_k*(1-η)) - exp(-r1_k*(1+η))
    num_dy = exp(-r1_k*(1-η))*(r1_k/a) + exp(-r1_k*(1+η))*(r1_k/a)
    den = 1+exp(-2*r1_k)
    H2 = (r2_k/N)*(num/den)
    H2_dy = (r2_k/N)*(num_dy/den)

    num = exp(-r2_k*(1-η)) - exp(-r2_k*(1+η))
    num_dy = exp(-r2_k*(1-η))*(r2_k/a) + exp(-r2_k*(1+η))*(r2_k/a)
    den = 1+exp(-2*r2_k)
    H3 = (r1_k/N)*(num/den)
    H3_dy = (r1_k/N)*(num_dy/den)

    H_dx += -2*(-1)^k * sin(α_k * ξ)/(a*l*α_k^2) * (H2 - H3)
    H_dy += 2*(-1)^k * cos(α_k * ξ)/(l*α_k^3) * (H2_dy - H3_dy)
  end
  j_x = a^2*σ^0.5 / μ^0.5 * (-grad_pz) * H_dy
  j_y = a^2*σ^0.5 / μ^0.5 * (-grad_pz) * (-H_dx)

  VectorValue(j_x,j_y,0.0)
end

function snes_postpro(cache,info)
  snes = cache.snes[]
  i_petsc = Ref{PetscInt}()
  @check_error_code PETSC.SNESGetIterationNumber(snes,i_petsc)
  info[:nls_iters] = Int(i_petsc[])
  @check_error_code PETSC.SNESGetLinearSolveIterations(snes,i_petsc)
  info[:ls_iters] = Int(i_petsc[])
  nothing
end


text:
#takes a crsp daily file and makes it weekly
function makeweekly(selectionop::Function = (noopd(x)=x),
    weeklyop::Function = (noopw(x)=x);
    crspnameout=PARAM[:crspweekly],
    binprefix="w_",
    mainpanel::Bool = true,
    computeivol::Bool = PARAM[:computeivol])

  local crsp::DataFrame

  #WARNING- crspd is NOT valid, below is only to create crspw
  crsp = prepcrsp(crsptype=PARAM[:crspdaily],
    csvextension = CSV_EXTENSION,
    incsvstream = IN_CSV_STREAM,
    binprefix = binprefix,
    refreshcrsp=PARAM[:refreshweekly],
    crspcolumns=[PARAM[:weeklycrspcols]; keys(PARAM[:momentumstrategies]) |> collect],
    crsprequiredcolumns=PARAM[:dailycrsprequired]) do crspd

    @assert (PARAM[:crspfrequencysuffix] == "w") || (!mainpanel)
    sort!(crspd, [:permno, :date]) #cache the sort

    #computerangemeasures!(crspd, sorted=true) #not used, mainly for kyle's lamda
    crspd.dvol = computedollarvol(crspd)
    rename!(crspd, :vol=>:svol)
    #crspd.shares = computeshares(crspd)

    #NOTE: testing note- change to trailingsigmatest! to test vol routine
    #println(" checksum: $(sum(skipmissing(crspd.sigma)))")

    crspd = mergeff(crspd, Fgroup=:permno, datefrequency=:day)
    if !issorted(crspd, [:permno, :date])
      sort!(crspd, [:permno, :date])
    end
    #crspd = dailyop(crspd)

    crspw::DataFrame = makeweekly(selectionop, crspd; computeivol)
    validatedates(crspw.date |> unique |> sort!; frequency=:week, validateexactfrequency=true)

    sort!(crspw, [:permno, :date])
    validatesummaryintegrity(crspw, Fid=:weekid)

    if PARAM[:testweeklycsv]
      testpermno::Int = PARAM[:testpermnomult]
      #testcrspw::DataFrame = crspw[crspw.permno .% testpermno .== 0, :]
      testcrspw::DataFrame = crspw[crspw.permno .== 14541, :] |> deepcopy
      #if size(testcrspw,1)>1_000 #prevents scenarios where the file is larger than needed
      #  testcrspw=testcrspw[1:1_000]
      #end

      rename!(testcrspw, :date=>:datew, :dvol=>:dvolw, :ivol=>:ivolw,
        :ret=>:retw, :price=>:pricew, :tridx=>:tridxw,
        :shares=>:sharesw, :mc=>:mcw, :sigma=>:sigmaw,
        :plret=>:plretw, :exchcd=>:exchcdw, :siccd=>:siccdw,)
      rename!(testcrspw, (k->k=>Symbol(k, :w)).(keys(PARAM[:momentumstrategies])))

      #make sure we merge using the same df as used to create the weekyl df
      testcrspout = innerjoin(
        testcrspw, crspd,#view(crspd,completecases(crspd, PARAM[:weeklycrsprequired]),:),
        on=[:permno, :weekid])
      sort!(testcrspout, [:permno, :weekid])
      testcrspout |> CSV.write(
        "$(PARAM[:testpath])\\$(binprefix)dailyforweekly_$(PARAM[:crspdaily]).csv")

    end

    crspw = weeklyop(crspw)

    return crspw
  end #end do block

  #println(describe(crsp))
  return crsp
end

#aggregates a daily dataframe into a weekly dataframe
function makeweekly(selectionop::Function, crspd::DataFrame;
  computeivol::Bool)

  #Assign each date to a weekid
  crspd.weekid = assignweekid(crspd.date)
  crspd.retP1 = 1.0 .+ crspd.ret #convert to gross returns for geometric compounding
  crspd.absret = (abs).(crspd.ret)

  issorted(crspd, [:permno, :date]) || error("crspd must be sorted in makeweekly")
  computeweekend!(crspd)
  trailingsigma!(crspd, Fendpoint=:weekend)
  select!(crspd, Not(:weekend))

  #compute change in week
  computeweekend!(crspd)
  if computeivol
    regressatperiodend!(crspd,
      minpoints=PARAM[:sigmaminpoints],
      maxpoints=PARAM[:sigmamaxpoints], #this is mainly for technical reasons
      calendarinterval = PARAM[:sigmacalendarinterval],
      FXs = PARAM[:ivolfields],
      Fagg = :ivol,
      aggfunc = std,
      Fendpoint = :weekend,
      #=aggfunc = (lm::FMLM)->std(lm.ε)=#)
  else
    @info "ivol not computed"
    crspd.ivol = missings(Float64, nrow(crspd))
  end

  #drop rows with too much missing data and apply selection criteria (e.g. min price, etc)
  @assert issorted(crspd, [:permno, :date])
  oldnrow = nrow(crspd)
  crspd::DataFrame = crspd[completecases(crspd, PARAM[:weeklycrsprequired]), :]
  @assert (size(crspd,1) ≥ oldnrow * 0.5) "weeklycrsprequired dropped more than 50% of rows!!!"

  valid = selectionop(crspd, Ftimegroup=:weekid, datefrequency=:week)
  #valid = combine(groupby(dvalid[!, [:permno, :weekid, ]], [:permno, :weekid]), nrow => :ndays)
  #valid = ensurecompleteness(valid, Ftimegroup=:weekid, datefrequency=:week)
  crspd = innerjoin(crspd, valid, on=[:permno, :weekid], validate=(false,true))

  crspd.plret = adjustindustry(crspd)
  #crspd.plret = Vector{Float64}(crspd.plret)
  crspd.lplret = (log).(crspd.plret .+ 1.0)
  if !issorted(crspd, [:permno, :date])
    sort!(crspd, [:permno, :date])
  end

  momentumstrategies!(crspd, Fconditional=:weekend)
  select!(crspd, Not(:weekend))

  if !issorted(crspd, [:permno, :date])
    sort!(crspd, [:permno, :date])
  end

  #combine by week
  Fmomentums = keys(PARAM[:momentumstrategies])  |> collect
  s2crspds::GroupedDataFrame = groupby(crspd, [:permno, :weekid])
  crspw::DataFrame = combine(s2crspds,
    #:permno => length=> :ndays,
    :ndays => last => :ndays,
    :date => last=> :date,
    :sigma=>last=> :sigma,
    :shares=>last=> :shares,
    #:svol => sum=> :svol,
    :dvol => sum=> :dvol,
    :price => last=> :price,
    :tridx => last => :tridx,
    #:price => mean=> :priceavg,
    :mc => last => :mc,
    :retP1 => prod=> :ret, #WARNING- to be adjusted later
    #:absret=>sum=> :sumabsret,
    :ivol => last=> :ivol,
    :exchcd => last => :exchcd,
    :siccd => last => :siccd,
    :lplret => sum => :plret, #WARNING- to be adjusted later
    Fmomentums .=> last .=> Fmomentums,
    #  :ret => std=> :std,
    #yzstd = [:hi, :lo, :open, :price, :Ldistadjclose]=>yzstd
    )

  crspw.ret .-= 1.0
  crspw.plret .= exp.(crspw.plret) .- 1.0

  testcrspws = groupby(crspw, [:permno, :weekid])
  Threads.@threads for k ∈ keys(s2crspds)
    testscrspw = testcrspws[(permno=k.permno,weekid=k.weekid)]
    @assert sum(log.(1 .+ s2crspds[k].ret)) |> exp ≈ testscrspw.ret[1] + 1.0
  end

  #conform the dates
  scrspws::GroupedDataFrame = groupby(crspw, :weekid)
  Threads.@threads for i ∈ 1:length(scrspws)
    scrspw::SubDataFrame = scrspws[i]
    @assert (maximum(scrspw.date) - minimum(scrspw.date)).value < 7
    scrspw.date .= maximum(scrspw.date)
  end

  return crspw
end

#recommend calling this function whenever we want only the ends of the week, and deleting afterwar
function computeweekend!(crspd::DataFrame; createnew::Bool = true)
  #compute change in week
  (:weekend ∉ names(crspd)) || (error(":weekend already in crspd.
    You probably should have deleted it the last time you used it."))
  crspd.weekend = falses(size(crspd,1))
  @assert issorted(crspd, [:permno, :date])

  scrspds = groupby(crspd, [:permno, :weekid])
  Threads.@threads for i::Int ∈ 1:length(scrspds)
    scrspd::SubDataFrame = scrspds[i]
    scrspd[end, :weekend] = true
  end

  return nothing
end

#assigns a weekid for each date
function assignweekid(dates::AbstractVector{Date}, mindate::Date=Date(1900,1,1))::Vector{Int}
  dayofweek(mindate) == 1 || (
    throw("day of $mindate ≠ 1, so day of week will not align with dayofweek function"))

  maxdate = maximum(dates)
  weekid::Vector{Int} = (d::Date->(d-mindate).value ÷ 7).(dates)

  return weekid
end

avgabs(v::Vector{MFloat64}) = mean((abs).(v))

function validatesummaryintegrity(crsp::DataFrame;
    Fid::Symbol=error("summary id Fid is required"))
  local scrsps::GroupedDataFrame

  #check weekid=>date
  scrsps = groupby(crsp, Fid)
  Threads.@threads for i ∈ 1:length(scrsps)
    scrsp::SubDataFrame = scrsps[i]
    (maximum(scrsp.date) ≠ minimum(scrsp.date)) && (
      error("multiple dates found for 1 Fid $Fid"))
  end

  scrsps = groupby(crsp, :date)
  Threads.@threads for i ∈ 1:length(scrsps)
    scrsp::SubDataFrame = scrsps[i]
    (maximum(scrsp[!, Fid]) ≠ minimum(scrsp[!, Fid])) && (
      error("multiple Fid $Fids found for 1 date"))
  end

  @info "summary data validated using summary id $Fid"
  return nothing
end

function regressatperiodend!(crspd::AbstractDataFrame;
  minpoints::Int = error("minpoints is required"),
  maxpoints::Int = error("maxpoints is required"),
  calendarinterval::DatePeriod = error("calendarinterval is required"),
  Fret::Symbol = :ret,
  Fendpoint::Symbol = error("Fendpoint is required"),
  FXs::Vector{Symbol} =error("FXs are required"),
  aggfunc::Function = error("aggfunc is required"),
  Fagg::Symbol = error("Fagg is required"))

  #local declarations
  local scrspds::GroupedDataFrame

  #this will hold the results
  crspd[!, Fagg] =  Vector{MFloat64}(undef, size(crspd,1))

  #form the regression formula and variable names
  Xexpr::FMExpr = Meta.parse(join((string).(FXs), " + "))
  Xnames::Vector{Symbol} = [:intercept; FXs]
  Yname::Symbol = Fret

  #avoid missing values
  completecrsp::SubDataFrame = view(crspd, completecases(crspd, [Fret; FXs]), :)

  #pre-allocate space before the threaded regressions
  #Since we have variable sized arrays and matrixes, we pre-allocate for each size

  K::Int = length(Xnames)
  Xs::Vector{Matrix{Float64}} = (t::Int->
    Matrix{Float64}(undef, maxpoints, K)).(1:Threads.nthreads())
  ys::Vector{Vector{Float64}} = (t::Int->
    Vector{Float64}(undef, maxpoints)).(1:Threads.nthreads())
  εs::Vector{Vector{Float64}} = (t::Int->
    Vector{Float64}(undef, maxpoints)).(1:Threads.nthreads())

  #first column is the intercept
  for t ∈ 1:Threads.nthreads()
    Xs[t][:, 1] .= ones(maxpoints)
  end

  #helper function to get the residuals
  @inline function resid!(X::AbstractMatrix, y::AbstractVector, m::Vector{Float64}, ε::Vector{Float64})
    ε .= X*m .- y
    return nothing
  end
  #this function gets the beginning of the interval
  @inline getj₀(d::Date, dates::AbstractVector{Date}) = searchsortedfirst(
    dates, d + Day(1) - calendarinterval)

  scrspds = groupby(completecrsp, :permno)
  regressionctr::Vector{Int} = zeros(Int, Threads.nthreads())
  print("Running ivol regressions... ")
  #=@time Threads.@threads =#for i ∈ 1:length(scrspds)
    local scrspd::SubDataFrame = scrspds[i]
    local t::Int = Threads.threadid()

    nrows::Int = size(scrspd,1)
    (nrows < minpoints) && continue

    #loop over possible bounds
    for j::Int ∈ minpoints:nrows
      regressionctr[t] += 1
      #only running regressions at the end of the week
      scrspd[j, Fendpoint] || continue

      #check for stale windows
      j₀::Int = getj₀(scrspd[j, :date], scrspd.date)
      N::Int = j-j₀ +1
      (N < minpoints)  && continue #elliminates windows without enough points

      #get appropriately sized matrices and arrays
      sX::SubArray = view(Xs[t], 1:N, :)
      resize!(ys[t], N)
      resize!(εs[t], N)

      @assert length(ys[t]) == length(εs[t]) == size(sX,1)

      #populate Y vector and X matrix
      ys[t] .= scrspd[j₀:j, Fret]
      for (c, FX) ∈ enumerate(FXs)
        sX[:, c+1] .= scrspd[j₀:j, FX]
      end

      resid!(sX, ys[t], GLM.coef(lm(sX, ys[t])), εs[t])
      scrspd[j, Fagg] = aggfunc(εs[t])
      #scrspd[j, Fagg] = aggfunc(FMLM(Xs[t], ys[t], Xnames=Xnames, Yname=Yname))
    end
  end
  println("$(sum(regressionctr)) regressions complete.")

  return nothing
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
export partialDependence

function partialDependence(mdl::Tree, dtmtable;thisIdx::Vector{Int}=dtmtable.trnidx,listOfFeatures::Vector{Symbol}=propertynames(dtmtable.features))
    @assert issubset(listOfFeatures, propertynames(dtmtable.features)) "DTM: At least one of the feature names you provided were not found in the names of the feature dataframe\r\nYou provided: $(listOfFeatures)\rpropertynames(dtmtable.features)=$(propertynames(dtmtable.features))"

    copyOfFeatures = deepcopy(dtmtable.features[thisIdx,:])

    nobs = size(copyOfFeatures, 1)
    fitReused = zeros(Float64, nobs)
    leafnrsReused = zeros(Int, nobs)
    idx = collect(1:nobs)
    # listOfFeatures=[:Density,:BonusMalus]

    res = Dict{Symbol,DataFrames.DataFrame}()
    for thisFeature in listOfFeatures
        forig = deepcopy(copyOfFeatures[thisFeature])
        pool = copyOfFeatures[thisFeature].pool
        predictedValues = zeros(Float64, length(pool))
        for i = 1:length(pool)
            copyOfFeatures[thisFeature] .= pool[i]
            DecisionTrees.apply_tree_by_leaf_iteration!(idx, mdl.rootnode, copyOfFeatures, fitReused, leafnrsReused)
            predictedValues[i] = Statistics.mean(fitReused)
        end
        # save result in dict
        df1 = DataFrame(Values=deepcopy(pool), AveragePredictions=copy(predictedValues))
        res[thisFeature] = deepcopy(df1)
        # reset values to original values
        copyOfFeatures[thisFeature] = deepcopy(forig)
    end
    return res
end

# this function is currently very slow....
function partialDependence(mdl::BoostedTree, dtmtable;fitSelection=:RawEstimate,thisIdx::Vector{Int}=dtmtable.trnidx,listOfFeatures::Vector{Symbol}=propertynames(dtmtable.features))
    @assert issubset(listOfFeatures, propertynames(dtmtable.features)) "DTM: At least one of the feature names you provided were not found in the names of the feature dataframe\r\nYou provided: $(listOfFeatures)\rpropertynames(dtmtable.features)=$(propertynames(dtmtable.features))"

    copyOfFeatures = deepcopy(dtmtable.features[thisIdx,:])

    nobs = size(copyOfFeatures, 1)
    # fitReused=zeros(Float64,nobs)
    # leafnrsReused=zeros(Int,nobs)
    idx = collect(1:nobs)
    # listOfFeatures=[:Density,:BonusMalus]

    res = Dict{Symbol,DataFrames.DataFrame}()
    for thisFeature in listOfFeatures
        forig = deepcopy(copyOfFeatures[thisFeature])
        pool = copyOfFeatures[thisFeature].pool
        predictedValues = zeros(Float64, length(pool))
        for i = 1:length(pool)
            copyOfFeatures[thisFeature] .= pool[i]
            result = predict(mdl, copyOfFeatures)
            predictedValues[i] = Statistics.mean(result[fitSelection])
        end
        # save result in dict
        df1 = DataFrame(Values=deepcopy(pool), AveragePredictions=copy(predictedValues))
        res[thisFeature] = deepcopy(df1)
        # reset values to original values
        copyOfFeatures[thisFeature] = deepcopy(forig)
    end
    return res
end

# partialDependence(resM,dtmtable,listOfFeatures=[:Density,:VehGas])


text:
using Test
using TightBindingLattice

using LinearAlgebra
using YAML

@testset "TranslationSymmetry" begin
    @testset "constructors" begin
        TranslationSymmetry(4*ones(Int, (1, 1)))
        TranslationSymmetry([4 0; 0 4])
        @test_throws ErrorException TranslationSymmetry([4 0 0; 0 4 0; 0 0 4])  # <- Temporary
        @test_throws ArgumentError TranslationSymmetry(OrthoCube([2 0; 0 2]), [2 0; 0 1])
    end

    @testset "properties" begin
        tsym = TranslationSymmetry([3 0; 0 1])
        mtab = [1 2 3; 2 3 1; 3 1 2]
        @test eltype(tsym) == TranslationOperation{Int}
        @test valtype(tsym) == TranslationOperation{Int}
        @test group(tsym) == FiniteGroup(mtab)
        @test group_order(tsym) == 3
        @test group_order(tsym, 1) == 1
        @test group_order(tsym, 2) == 3
        @test group_order(tsym, 3) == 3
        @test group_multiplication_table(tsym) == mtab
        @test element(tsym, 1) == TranslationOperation([0,0])
        @test element(tsym, 2) == TranslationOperation([1,0])
        @test element(tsym, 3) == TranslationOperation([2,0])
        @test elements(tsym) == TranslationOperation.([[0,0], [1,0], [2,0]])
        @test [element_name(tsym, i) for i in 1:3] == element_names(tsym)
        let ω = cis(2π/3)
            @test character_table(tsym) ≈ [1 1 1; 1 ω^2 ω; 1 ω ω^2]
        end
        @test length(irreps(tsym)) == 3
        @test num_irreps(tsym) == 3
        @test irrep(tsym, 1) == irreps(tsym)[1]
        @test irrep(tsym, 2) == irreps(tsym)[2]
        @test irrep(tsym, 3) == irreps(tsym)[3]
        @test all(irrep_dimension(tsym, i) == 1 for i in 1:3)
        @test generator_indices(tsym) == [2, 1]
        @test generator_elements(tsym) == [TranslationOperation([1,0]), TranslationOperation([0,0])]

        let n = lowercase(symmetry_name(tsym))
            @test occursin("translation", n)
            @test occursin("3", n)
            @test occursin("1", n)
            @test occursin("0", n)
        end
    end

    @testset "orthogonal lattice" begin
        tsym = TranslationSymmetry([3 0; 0 3])
        @test isabelian(tsym.group)

        @test length(tsym.generators) == 2
        idx_gen1 = tsym.generators[1]
        idx_gen2 = tsym.generators[2]

        @test group_order(tsym) == 9

        @testset "generators" begin
            # Generators completely generate the group
            @test generate_subgroup(tsym.group, tsym.generators) == BitSet(1:group_order(tsym))
            @test prod(tsym.group.period_lengths[x] for x in tsym.generators) == group_order(tsym) # since abelian group
        end

        @testset "elements" begin
            @test tsym.coordinates[idx_gen1] == [1, 0]
            @test tsym.coordinates[idx_gen2] == [0, 1]
            # elements ordered according to the "generator" (i.e. orthogonal order)
            @test tsym.element_names == [
                "[0, 0]", "[1, 0]", "[2, 0]",
                "[0, 1]", "[1, 1]", "[2, 1]",
                "[0, 2]", "[1, 2]", "[2, 2]",
            ]
            @test length(element_names(tsym)) == 9
            for i in 1:9
                @test element_names(tsym)[i] == element_name(tsym, i)
            end
        end

        @testset "multiplication table" begin
            mtab = zeros(Int, (9, 9))
            lookup = Dict(r => i for (i, r) in enumerate(tsym.coordinates))
            for (i, ri) in enumerate(tsym.coordinates)
                for (j, rj) in enumerate(tsym.coordinates)
                    _, ek = tsym.orthocube.wrap(ri + rj)
                    k = lookup[ek]
                    mtab[i,j] = k
                end # for j
            end # for i
            @test mtab == group_multiplication_table(tsym)
            @test mtab == group_multiplication_table(tsym.group)
        end

        @testset "coordinates" begin
            @test tsym.orthogonal_shape == [3, 3]
            @test tsym.orthogonal_coordinates == tsym.coordinates
            for cc in tsym.coordinates
                oc = tsym.coordinate_to_orthogonal_map[cc]
                cc2 = tsym.orthogonal_to_coordinate_map[oc]
                @test cc == cc2
            end

            for oc in tsym.orthogonal_coordinates
                cc = tsym.orthogonal_to_coordinate_map[oc]
                oc2 = tsym.coordinate_to_orthogonal_map[cc]
                @test oc == oc2
            end
        end

        @testset "symmetry_product" begin
            pro = symmetry_product(tsym)
            @test all(
                pro(x,y).displacement == tsym.orthocube.wrap(x.displacement + y.displacement)[2]
                for x in tsym, y in tsym
            )
        end

        @testset "iterate/membership" begin
            @test [x for x in tsym] == collect(elements(tsym))
            @test all(x ∈ tsym for x in tsym)
            @test 100 ∉ tsym
            @test IdentityOperation(Int, 2) ∈ tsym
            @test IdentityOperation(Int, 3) ∉ tsym
            @test TranslationOperation([100, 100]) ∈ tsym
            @test TranslationOperation([100, 100, 100]) ∉ tsym
            @test PointOperation([1 0; 0 1]) ∈ tsym
            @test PointOperation([0 1; 1 0]) ∉ tsym
            @test SpaceOperation([1 0; 0 1], [1, 0]) ∈ tsym
            @test SpaceOperation([0 1; 1 0], [1, 0]) ∉ tsym
        end

        @testset "irreps" begin
            @test tsym.conjugacy_classes == [[i] for (i, x) in enumerate(tsym.element_names)]
            @test size(character_table(tsym)) == (9, 9)
            @test length(irreps(tsym)) == 9
            @test num_irreps(tsym) == 9
            for i in 1:9
                @test irreps(tsym)[i] == irrep(tsym, i)
                @test irrep_dimension(tsym, i) == 1
            end
        end
    end # @testset "orthogonal lattice" begin

    @testset "non-orthogonal lattice" begin
        tsym = TranslationSymmetry([4 0; 0 3])

        @test length(tsym.generators) == 2
        # idx_gen = tsym.generators[1]
        @test tsym.coordinates[tsym.generators[1]] == [1, 0]
        @test tsym.coordinates[tsym.generators[2]] == [0, 1]
        # elements ordered according to the "generator" (i.e. orthogonal order)
        @test element_names(tsym) == [
            "[0, 0]", "[1, 0]", "[2, 0]", "[3, 0]",
            "[0, 1]", "[1, 1]", "[2, 1]", "[3, 1]",
            "[0, 2]", "[1, 2]", "[2, 2]", "[3, 2]",
        ]
        @test tsym.conjugacy_classes == [[i] for (i, x) in enumerate(tsym.element_names)]

        @test tsym.orthogonal_shape == [4,3]
        @test tsym.orthogonal_coordinates == [[0,0], [1,0], [2,0], [3,0], [0,1], [1,1], [2,1], [3,1], [0,2], [1,2], [2,2], [3,2]]
        for cc in tsym.coordinates
            oc = tsym.coordinate_to_orthogonal_map[cc]
            cc2 = tsym.orthogonal_to_coordinate_map[oc]
            @test cc == cc2
        end

        for oc in tsym.orthogonal_coordinates
            cc = tsym.orthogonal_to_coordinate_map[oc]
            oc2 = tsym.coordinate_to_orthogonal_map[cc]
            @test oc == oc2
        end
    end  # @testset "non-orthogonal lattice" begin

    @testset "lattice permutation" begin
        unitcell = make_unitcell([1.0 0.0; 0.0 1.0]; SiteType=String)
        addsite!(unitcell, "Ox", FractCoord([0,0], [0.5, 0.0]))
        addsite!(unitcell, "Oy", FractCoord([0,0], [0.0, 0.5]))
        lattice = make_lattice(unitcell, [4 0; 0 4])
        tsym = TranslationSymmetry(lattice)

        let
            lattice_failure = make_lattice(unitcell, [4 0; 0 3])
            #@test_throws ArgumentError get_site_permutations(lattice_failure, tsym)
            @test_throws ArgumentError embed(lattice_failure, tsym)
        end
        perms = [embed(lattice, op) for op in tsym.elements]
        #tsymbed = embed(lattice, tsym)
        @test length(perms) == 16
        @test perms[1].permutation.map == 1:32  # identity

        for (ir, r) in enumerate(lattice.bravais_coordinates)
            @test perms[ir].permutation.map == [
                let
                    ivec = [(i-1) % 4, (i-1) ÷ 4]
                    jvec = [(x + 16) % 4 for x in (ivec .+ r)]
                    j = jvec[1] + jvec[2]*4 + 1
                    2 * (j-1) + iorb
                end
                for i in 1:16 for iorb in 1:2
            ]
        end

        tsymbed = embed(lattice, tsym)
        @test_throws ArgumentError get_irrep_iterator(IrrepComponent(tsymbed, 1, 2))
        let phases = [cis(-2π * i / 4) for j in 0:3 for i in 0:3]
            tsic = IrrepComponent(tsymbed, 2, 1)
            irrep_list1 = collect(get_irrep_iterator(tsic))
            irrep_list2 = collect(zip(perms, phases))
            @test length(irrep_list1) == length(irrep_list2)
            for (x,y) in zip(irrep_list1, irrep_list2)
                @test x[1] == y[1]
                @test x[2] ≈ y[2]
            end
        end

        @test  isbragg(tsym.fractional_momenta[1], [1,0]) # Γ point
        @test !isbragg(tsym.fractional_momenta[2], [1,0]) # Γ point
        @test !isbragg(tsym.fractional_momenta[2], [[0,0], [1,0]])
    end # testset lattice permutation

    @testset "reduction" begin
        @test_throws DimensionMismatch isbragg([4, 0], [2,0,0], [2,0])
        @test_throws DimensionMismatch isbragg([4, 4], [2,0], [2,0,0])

        @test_throws ArgumentError isbragg([4, 0], [2,0], [2,0])
        @test_throws ArgumentError isbragg([-2, 2], [2,0], [2,0])

        tsym = TranslationSymmetry([3 0; 0 3])
        # Gamma point is always ok
        @test  isbragg([3,3], [0,0], [0,0])
        @test  isbragg([3,3], [0,0], [1,0])
        @test  isbragg([3,3], [0,0], [2,0])

        @test  isbragg(tsym, 1, TranslationOperation([0,0]))
        @test  isbragg(tsym, 1, TranslationOperation([1,0]))
        @test  isbragg(tsym, 1, TranslationOperation([2,0]))

        # non-zero momentum depends on what the identity translation is
        @test  isbragg([3,3], [1,0], [0,0]) # zero translation is always identity, so it's always fine
        @test !isbragg([3,3], [1,0], [1,0]) # these two translations are not compatible
        @test !isbragg([3,3], [1,0], [2,0]) #   with momentum [1,1]

        @test  isbragg(tsym, 2, TranslationOperation([0,0]))
        @test !isbragg(tsym, 2, TranslationOperation([1,0]))
        @test !isbragg(tsym, 2, TranslationOperation([2,0]))

        @test  isbragg([3,3], [0,0], [[0,0], [1,0], [2,0]])
        @test !isbragg([3,3], [1,0], [[0,0], [1,0], [2,0]])

        @test  isbragg(tsym, 1, TranslationOperation.([[0,0], [1,0], [2,0]]))
        @test !isbragg(tsym, 2, TranslationOperation.([[0,0], [1,0], [2,0]]))

        @test_throws DimensionMismatch isbragg([1//2, 2//4, 0//2], [2,0])

        # Gamma point is always ok
        @test  isbragg([0,0] .// [3,3], [0,0])
        @test  isbragg([0,0] .// [3,3], [1,0])
        @test  isbragg([0,0] .// [3,3], [2,0])

        # non-zero momentum depends on what the identity translation is
        @test  isbragg([1,0] .// [3,3], [0,0]) # zero translation is always identity, so it's always fine
        @test !isbragg([1,0] .// [3,3], [1,0]) # these two translations are not compatible
        @test !isbragg([1,0] .// [3,3], [2,0]) #   with momentum [1,1]

        @test  isbragg([0,0] .// [3,3], [[0,0], [1,0], [2,0]])
        @test !isbragg([1,0] .// [3,3], [[0,0], [1,0], [2,0]])

        tsym = TranslationSymmetry([4 0; 0 4])
        @test  isbragg(tsym, 11, TranslationOperation([2,0]))
        @test !isbragg(tsym, 10, TranslationOperation([2,0]))
        @test  isbragg(tsym, 11, TranslationOperation.([[0,0], [2,0], [0,2], [2,2]]))

        tsym = TranslationSymmetry([6 2; -2 6])

        # TODO more tests for isbragg

        for irrep_index in 1:num_irreps(tsym)
            k = tsym.fractional_momenta[irrep_index]
            k_ortho = tsym.orthogonal_coordinates[irrep_index] .// tsym.orthogonal_shape
            for (t, t_ortho) in zip(elements(tsym), tsym.orthogonal_coordinates)
                @test mod(k⋅t.displacement, 1) == mod(k_ortho⋅t_ortho, 1)
            end
        end
    end
end # @testset "TranslationSymmetry" begin

#
# @testset "translation" begin
#   @testset "constructor exceptions" begin
#     @test_throws ArgumentError TranslationGroup([Permutation([2,3,4,1]), Permutation([3,4,1,2])])
#     @test_throws ArgumentError TranslationGroup([Permutation([2,1,3,4]), Permutation([1,3,2,4])])
#   end
#
#   t1 = Permutation([2,3,1, 5,6,4])
#   t2 = Permutation([4,5,6, 1,2,3])
#   g = TranslationGroup([t1, t2])
#
#   @test g.generators == [t1, t2]
#   @test g.translations == [[0,0], [1,0], [2,0], [0,1], [1,1], [2,1]]
#   @test Set(g) == Set([t1^d1*t2^d2 for d1 in 0:2 for d2 in 0:1])
#   @test length(Set(g)) == 2*3
#   @test g.fractional_momenta == [[0//3, 0//2], [1//3, 0//2], [2//3, 0//2],
#                                  [0//3, 1//2], [1//3, 1//2], [2//3, 1//2]]
#   χ = [cis(-2π * (k ⋅ t)) for k in g.fractional_momenta, t in g.translations]
#   @test isapprox(g.character_table, χ; atol=Base.rtoldefault(Float64))
#
#   @test is_compatible([0//1, 0//1], [0,0])
#   @test !is_compatible([0//1, 1//2], [0,1])
#
#   @test is_compatible([0//1, 0//1], [[0,0]])
#   @test !is_compatible([0//1, 1//2], [[0,0], [0,1]])
# end
#
#
# @testset "non-orthogonal" begin
#   #=
#   . . . o .
#   o . . . .
#   . . . . o
#   . o . . .
#
#   2 3 4 1 5
#   1|5 6 7|2
#   7|2 3 4|1
#   4|1|5 6 7
#
#
#   =#
#   t1 = Permutation([5,3,4,1,6,7,2])
#   t2 = Permutation([2,5,6,7,3,4,1])
#   generators = [t1, t2]
#   #translation_group = TranslationGroup(t1, t2)
#   full_shape = [g.order for g in generators]
#   full_translations = vcat( collect( Iterators.product([0:g.order-1 for g in generators]...) )...)
#   full_translations = [ [x...] for x in full_translations]
#
#   elemtrans = Dict{Permutation, Vector{Int}}()
#   for (ig, dist) in enumerate(full_translations)
#     g = prod(gen^d for (gen, d) in zip(generators, dist))
#     if !haskey(elemtrans, g)
#       elemtrans[g] = [dist...]
#     end
#   end
#   elements = Permutation[]
#   translations = Vector{Int}[]
#   for (e, t) in elemtrans
#     push!(elements, e)
#     push!(translations, t)
#   end
#   idx = sortperm(translations)
#   translations = translations[idx]
#   elements = elements[idx]
# end


text:
struct ESCB2Budgeted <: ESCB2OptimisationAlgorithm
  ε::Float64 # Discretisation of the coefficients of the nonlinear part.
  solve_all_budgets_at_once::Bool # Some optimisation algorithms can solve the budgeted problems for all values of the budget at once.
  # TODO: compute a default value based on what the solver is able to do? Use Julia's applicable to determine it automatically?
end

function optimise_linear_sqrtlinear(instance::CombinatorialInstance{T}, algo::ESCB2Budgeted,
                                    linear::Dict{T, Float64}, sqrtlinear::Dict{T, Float64};
                                    with_trace::Bool=false) where T
  # Transform the linear term in a budget constraint, the nonlinear term
  # becoming the objective function (in which case the concave function
  # can be dropped).
  # For this to work, the linear part of the objective function must
  # take only integer values, i.e. the linear coefficients are actually
  # only integers.
  linear_discrete = Dict(k => round(Int, v / algo.ε, RoundUp) for (k, v) in linear)

  t0 = time_ns()

  # Maximum value of the linear term?
  max_budget = instance.n_arms * maximum(values(linear_discrete))

  # Solve the family of problems with an increasing budget.
  solutions = Dict{Int, Vector{T}}()
  if ! algo.solve_all_budgets_at_once # TODO: Replace this mandatory parameter by a function defined on the solvers, so that each of them can indicate what they support? Then, if this parameter is not set, the most efficient implementation is called if available; otherwise, follow this parameter (and still warn if the required implementation is not available).
    if ! applicable(solve_budgeted_linear, instance.solver, sqrtlinear, linear_discrete, max_budget)
      error("The function solve_budgeted_linear is not defined for the solver $(typeof(instance.solver)) for arguments of type ($(typeof(instance.solver)), $(typeof(sqrtlinear)), $(typeof(linear_discrete)), $(typeof(max_budget))).")
    end

    budget = 0
    while budget <= max_budget
      sol = solve_budgeted_linear(instance.solver, sqrtlinear, linear_discrete, budget)

      if length(sol) == 0 || sol == [-1]
        # Infeasible!
        for b in budget:max_budget
          solutions[b] = sol
        end
        break
      end

      # Feasible.
      sol_budget = sum(linear_discrete[arm] for arm in keys(linear_discrete) if arm in sol)
      for b in budget:sol_budget
        solutions[b] = sol
      end
      budget = sol_budget + 1
    end
  else
    if ! applicable(solve_all_budgeted_linear, instance.solver, sqrtlinear, linear_discrete, max_budget)
      error("The function solve_all_budgeted_linear is not defined for the solver $(typeof(instance.solver)) for arguments of type ($(typeof(instance.solver)), $(typeof(sqrtlinear)), $(typeof(linear_discrete)), $(typeof(max_budget))).")
    end

    solutions = solve_all_budgeted_linear(instance.solver, sqrtlinear, linear_discrete, max_budget)
  end

  # Take the best solution.
  best_solution = Int[]
  best_objective = -Inf
  for (budget, sol) in solutions
    # Ignore infeasible cases.
    if length(sol) == 0 || sol == [-1]
      continue
    end

    if sum(sqrtlinear[arm] for arm in keys(sqrtlinear) if arm in sol) < 0.0
      continue
    end

    # Compute the maximum.
    reward = escb2_index(linear, sqrtlinear, sol)
    if reward > best_objective
      best_solution = sol
      best_objective = reward
    end
  end
  t1 = time_ns()

  if with_trace
    run_details = ESCB2Details()
    run_details.n_iterations = 1
    run_details.best_objective = best_objective
    run_details.solver_time = (t1 - t0) / 1_000_000_000
    return best_solution, run_details
  else
    return best_solution
  end
end


text:
# Copyright 2015, 2016, 2017, 2018, 2019, 2020 Martin Holters
# See accompanying license file.

export Circuit, add!, connect!, disconnect!, @circuit, composite_element

import Base: delete!

struct CircuitNLFunc{Fs}
    fs::Fs
end

@inline (cnlf::CircuitNLFunc{Fs})(q) where {Fs} = _apply_all(q, cnlf.fs...)

_apply_all(q) = (SVector{0,Float64}(), SMatrix{0,0,Float64}())
@inline function _apply_all(q, f1::F, fs...) where F
    (res1, J1) = f1(q)
    (resrem, Jrem) = _apply_all(q, fs...)
    return ([res1; resrem], dcat(J1, Jrem))
end

dcat(a::SMatrix{Ma,Na,Ta}, b::SMatrix{Mb,Nb,Tb}) where {Ma,Na,Mb,Nb,Ta,Tb} =
    [[a zeros(SMatrix{Ma,Nb,Ta})]; [zeros(SMatrix{Mb,Na,Tb}) b]]

const Net = Vector{Tuple{Symbol,Symbol}} # pairs of element designator and pin name

struct Circuit
    elements::OrderedDict{Symbol, Element}
    nets :: Vector{Net}
    net_names :: Dict{Symbol, Net}
    Circuit() = new(OrderedDict{Symbol, Element}(), [], Dict{Symbol, Net}())
end

elements(c::Circuit) = values(c.elements)

for n in [:nb; :nx; :nq; :nu; :nl; :ny; :nn]
    @eval ($n)(c::Circuit) = reduce(+, init=0, $n(elem) for elem in elements(c))
end

for mat in [:mv; :mi; :mx; :mxd; :mq; :mu; :pv; :pi; :px; :pxd; :pq]
    # blockdiag() does not work, so include an empty matrix of desired type in
    # case elements(c) is empty
    # as blockdiag for unknown number of arguments cannot be inferred properly,
    # add type-assertion
    @eval ($mat)(c::Circuit) =
         blockdiag(spzeros(Rational{BigInt}, 0, 0),
                 (convert(SparseMatrixCSC{Rational{BigInt}}, elem.$mat)
                  for elem in elements(c))...
                )::SparseMatrixCSC{Rational{BigInt},Int}
end

u0(c::Circuit) = vcat((elem.u0 for elem in elements(c))...)

function incidence(c::Circuit)
    i = sizehint!(Int[], 2nb(c))
    j = sizehint!(Int[], 2nb(c))
    v = sizehint!(Int[], 2nb(c))
    for (row, pins) in enumerate(c.nets), (elemname, pinname) in pins
        elem = c.elements[elemname]
        offset = branch_offset(c, elem)
        bps = elem.pins[pinname]
        for (branch, polarity) in bps
            push!(i, row)
            push!(j, offset + branch)
            push!(v, polarity)
        end
    end
    # ensure zeros due to short-circuited branches are removed with dropzeros!
    dropzeros!(sparse(i, j, v, length(c.nets), nb(c)))
end

function nonlinear_eq_func(c::Circuit, elem_idxs=1:length(elements(c)))
    col_offset = 0
    funcs = Function[]
    for elem in collect(elements(c))[elem_idxs]
        if nn(elem) == 0 && nq(elem) == 0
            continue
        end
        push!(funcs,
            let q_indices=SVector{nq(elem),Int}(col_offset+1:col_offset+nq(elem)),
                nleqfunc=nonlinear_eq_func(elem)
                @inline function (q)
                    nleqfunc(q[q_indices])
                end
            end)

        col_offset += nq(elem)
    end
    return CircuitNLFunc((funcs...,))
end

"""
    add!(c::Circuit, elem::Element)

Adds the element `elem` to the circuit `c`, creating and returning a new, unique
reference designator, leaving its pins unconnected.
"""
function add!(c::Circuit, elem::Element)
    for (k, v) in c.elements
        if v == elem
            return k
        end
    end
    designator = gensym()
    add!(c, designator, elem)
    return designator
end

add!(c::Circuit, elems::Element...) = ((add!(c, elem) for elem in elems)...,)

"""
    add!(c::Circuit, designator::Symbol, elem::Element)

Adds the element `elem` to the circuit `c` with the reference designator
`designator`, leaving its pins unconnected. If the circuit already contained
an element named `designator`, it is removed first.
"""
function add!(c::Circuit, designator::Symbol, elem::Element)
    if haskey(c.elements, designator)
        delete!(c, designator)
    end
    for pin in keys(elem.pins)
        push!(c.nets, [(designator, pin)])
    end
    c.elements[designator] = elem
end

"""
    delete!(c::Circuit, designator::Symbol)

Deletes the element named `designator` from the circuit `c` (disconnecting all
its pins).
"""
function delete!(c::Circuit, designator::Symbol)
    for net in c.nets
        filter!(elempin -> elempin[1] != designator, net)
    end
    delete!(c.elements, designator)
end

function branch_offset(c::Circuit, elem::Element)
    offset = 0
    for el in elements(c)
        el == elem && return offset
        offset += nb(el)
    end
    throw(ArgumentError("Element not found in circuit"))
end

function netfor!(c::Circuit, p::Tuple{Symbol,Symbol})
    for net in c.nets
        p ∈ net && return net
    end
    throw(ArgumentError("Unknown pin $p"))
end
netfor!(c::Circuit, p::Tuple{Symbol,Any}) = netfor!(c, (p[1], Symbol(p[2])))

function netfor!(c::Circuit, name::Symbol)
    haskey(c.net_names, name) || push!(c.nets, get!(c.net_names, name, []))
    c.net_names[name]
end

"""
    connect!(c::Circuit, pins::Union{Symbol,Tuple{Symbol,Any}}...)

Connects the given pins (or named nets) to each other in the circuit `c`. Named
nets are given as `Symbol`s, pins are given as `Tuple{Symbols,Any}`s, where the
first entry is the reference designator of an element in `c`, and the second
entry is the pin name. For convenience, the latter is automatically converted to
a `Symbol` as needed.

# Example

```jldoctest; output = false, setup = :(using ACME)
circ = Circuit()
add!(circ, :r, resistor(1e3))
add!(circ, :src, voltagesource(5))
connect!(circ, (:src, -), (:r, 2), :gnd) # connect to gnd net

# output

```
"""
function connect!(c::Circuit, pins::Union{Symbol,Tuple{Symbol,Any}}...)
    nets = unique([netfor!(c, pin) for pin in pins])
    for net in nets[2:end]
        append!(nets[1], net)
        deleteat!(c.nets, findfirst(n -> n == net, c.nets))
        for (name, named_net) in c.net_names
            if named_net == net
                c.net_names[name] = nets[1]
            end
        end
    end
end

function disconnect!(c::Circuit, pin::Tuple{Symbol,Symbol})
    net = netfor!(c, pin)
    filter!(p -> p != pin, net)
    push!(c.nets, [pin])
end
"""
    disconnect!(c::Circuit, p::Tuple{Symbol,Any})

Disconnects the given pin `p` from anything else in the circuit `c`. The pin is
given as a`Tuple{Symbols,Any}`, where the first entry is the reference
designator of an element in `c`, and the second entry is the pin name. For
convenience, the latter is automatically converted to a `Symbol` as needed. Note
that if e.g. three pin `p1`, `p2`, and `p3` are connected then
`disconnect!(c, p1)` will disconnect `p1` from `p2` and `p3`, but leave `p2` and
`p3` connected to each other.
"""
disconnect!(c::Circuit, p::Tuple{Symbol,Any}) = disconnect!(c, (p[1], Symbol(p[2])))

function topomat!(incidence::SparseMatrixCSC{T}) where {T<:Integer}
    @assert all(x -> abs(x) == 1, nonzeros(incidence))
    @assert all(sum(incidence, dims=1) .== 0)

    t = falses(size(incidence)[2]);

    row = 1;
    for col = 1:size(incidence)[2]
        rows = filter(r -> r ≥ row, findall(!iszero, incidence[:, col]))
        @assert length(rows) ≤ 2

        isempty(rows) && continue
        t[col] = true;

        if rows[1] ≠ row
            incidence[[rows[1], row], :] = incidence[[row, rows[1]], :]
        end
        if length(rows) == 2
            @assert incidence[row, col] + incidence[rows[2], col] == 0
            incidence[rows[2],:] = incidence[rows[2],:] + incidence[row,:]
        end
        if incidence[row, col] < 0
            cols = findall(!iszero, incidence[row, :])
            incidence[row,cols] = -incidence[row,cols]
        end
        rows = findall(==(1), incidence[1:row-1, col])
        incidence[rows, :] .-= incidence[row, :]'
        rows = findall(==(-1), incidence[1:row-1, col])
        incidence[rows, :] .+= incidence[row, :]'
        row += 1
    end

    ti = incidence[1:row-1, :]

    dl = ti[:, (!).(t)]
    tv = spzeros(T, size(dl, 2), size(incidence, 2))
    # findall here works around JuliaLang/julia#27013 (introdcued in 0.7.0-DEV.4983)
    tv[:, findall(t)] = -dl'
    tv[:, findall((!).(t))] = SparseMatrixCSC{T}(I, size(dl, 2), size(dl, 2))

    tv, ti
end

topomat(incidence::SparseMatrixCSC{<:Integer}) = topomat!(copy(incidence))
topomat(c::Circuit) = topomat!(incidence(c))

@doc raw"""
    @circuit begin #= ... =# end

Provides a simple domain-specific language to decribe circuits. The
`begin`/`end` block can hold element definitions of the form
`refdes = elementfunc(params)` and connection specifications of the form
`refdes1[pin1] ⟷ refdes2[pin2]`.

# Example

To create a circuit with a voltage source connected to a resistor:
```jldoctest; output = false, setup = :(using ACME), filter = r"(ACME\.)?Circuit\(.*"s
@circuit begin
    src = voltagesource(5)
    r = resistor(1000)
    src[+] ⟷ r[1]
    src[-] ⟷ r[2]
end

# output

Circuit(...)
```

Alternatively, connection specifications can be given after an element
specification, separated by commas. In that case, the `refdes` may be omitted,
defaulting to the current element.

# Example

```jldoctest; output = false, setup = :(using ACME), filter = r"(ACME\.)?Circuit\(.*"s
@circuit begin
    src = voltagesource(5)
    r = resistor(1000), src[+] ⟷ [1], src[-] ⟷ [2]
end

# output

Circuit(...)
```

Finally, a connection endpoint may simply be of the form `netname`, to connect
to a named net. (Such named nets are created as needed.)

# Example

```jldoctest; output = false, setup = :(using ACME), filter = r"(ACME\.)?Circuit\(.*"s
@circuit begin
    src = voltagesource(5), [-] ⟷ gnd
    r = resistor(1000), [1] ⟷ src[+], [2] ⟷ gnd
end

# output

Circuit(...)
```

If a net or pin specification is not just a single symbol or number, and has to
be put in quotes (e.g. `"in+"`, `"9V"`)

!!! note
    Instead of `⟷` (`\longleftrightarrow`), one can also use `==`.
"""
macro circuit(cdef)
    is_conn_spec(expr::Expr) =
        (expr.head === :call && (expr.args[1] === :(⟷) || expr.args[1] === :(↔) || expr.args[1] === :(==))) ||
        (expr.head === :comparison && all(c -> c === :(==), expr.args[2:2:end]))
    is_conn_spec(::Any) = false

    function elem_spec(expr)
        if !isa(expr, Expr) || expr.head !== :(=)
            error("invalid element specification$locinfo: $(expr)")
        end
        if !isa(expr.args[1], Symbol)
            error("invalid element identifier$locinfo: $(expr.args[1])")
        end
        if isa(expr.args[2], Expr) && expr.args[2].head === :tuple
            if isempty(expr.args[2].args)
                error("invalid element specification$locinfo: $(expr.args[2])")
            end
            elemspec = expr.args[2].args[1]
            conn_exprs = expr.args[2].args[2:end]
        else
            elemspec = expr.args[2]
            conn_exprs = []
        end
        push!(ccode.args, :(add!(circ, $(QuoteNode(expr.args[1])), $(esc(elemspec)))))
        for conn_expr in conn_exprs
            if !is_conn_spec(conn_expr)
                error("invalid connection specification$locinfo: $conn_expr")
            end
            push!(ccode.args, Expr(:call, :connect!, :circ, extractpins(conn_expr, expr.args[1])...))
        end
    end

    function extractpins(expr::Expr, default_element=nothing)
        if expr.head === :call && (expr.args[1] === :(⟷) || expr.args[1] === :(↔) || expr.args[1] === :(==))
            return vcat((extractpins(a, default_element) for a in expr.args[2:end])...)
        elseif expr.head === :comparison && all(c -> c === :(==), expr.args[2:2:end])
            return vcat((extractpins(a, default_element) for a in expr.args[1:2:end])...)
        elseif expr.head === :ref
            return [:(($(QuoteNode(expr.args[1])), $(QuoteNode(expr.args[2]))))]
        elseif expr.head === :vect && length(expr.args) == 1
            if default_element === nothing
                error("missing element$(locinfo): $expr")
            end
            return [:(($(QuoteNode(default_element)), $(QuoteNode(expr.args[1]))))]
        else
            error("invalid pin specification$(locinfo): $expr")
        end
    end

    function extractpins(netname::Symbol, default_element=nothing)
        return [QuoteNode(netname)]
    end

    extractpins(netname::String, default_element=nothing) =
        extractpins(Symbol(netname), default_element)

    if !isa(cdef, Expr) || cdef.head !== :block
        error("@circuit must be followed by a begin/end block")
    end
    ccode = Expr(:block)
    push!(ccode.args, :(circ = Circuit()))
    locinfo = ""
    for expr in cdef.args
        if isa(expr, LineNumberNode)
            locinfo = " at $(expr.file):$(expr.line)"
            continue
        end
        if !isa(expr, Expr)
            error("invalid statement in circuit definition$locinfo: $expr")
        end
        if expr.head === :line
            locinfo = " at $(expr.args[2]):$(expr.args[1])"
        elseif expr.head === :(=)
            elem_spec(expr)
        elseif is_conn_spec(expr)
            push!(ccode.args, Expr(:call, :connect!, :circ, extractpins(expr)...))
        else
            error("invalid statement in circuit definition$locinfo: $expr")
        end
    end

    push!(ccode.args, :(circ))
    return ccode
end

@doc raw"""
    composite_element(circ; pinmap=Dict(), ports)

Create a circuit element from the (sub-)circuit `circ`. The `pinmap` defines
the mapping from pin names of the element to be created to pins (or nets) in
`circ`. Optionally, `ports` can be used to explicitly specify (as a list `Pair`s
of pins) the ports to use. By default, a port will be created between one
arbitrarily chosen pin and every other pin.

# Example

The following will create an element for a 2.5V source, created using a 5V
source and a voltage divider, stabilized with a capacitor.

```jldoctest; output = false, setup = :(using ACME), filter = r"(ACME\.)?Element\(.*"s
circ = @circuit begin
   r1 = resistor(10e3)
   r2 = resistor(10e3), [1] == r1[2]
   c = capacitor(1e-6), [1] == r2[1], [2] == r2[2]
   src = voltagesource(5), [+] == r1[1], [-] == r2[2]
end

composite_element(circ, pinmap=Dict(1 => (:r2, 1), 2 => (:r2, 2)))

# output

Element(...)
```

Note that the pins still implicitly specifiy ports, so an even number must be
given, but the same pin may be given multiple times to be part of multiple
ports.
""" function composite_element(circ::Circuit; pinmap=Dict(), ports=ports_from_pinmap(pinmap))
    if ny(circ) > 0
        throw(ArgumentError("creating composite elements from circuits with outputs is not supported"))
    end
    numports = length(ports)
    # construct system matrix with norators for connection ports included
    Mᵥ = blockdiag(mv(circ), spzeros(numports, numports))
    Mᵢ = blockdiag(mi(circ), spzeros(numports, numports))
    Mₓ = [mx(circ); spzeros(numports, nx(circ))]
    Mₓ´ = [mxd(circ); spzeros(numports, nx(circ))]
    Mq = [mq(circ); spzeros(numports, nq(circ))]
    Mu = [mu(circ); spzeros(numports, nu(circ))]
    u0 = [ACME.u0(circ); spzeros(numports)]
    incid = [incidence(circ) spzeros(Int, length(circ.nets), numports)]
    for i in eachindex(ports)
        port = ports[i]
        net = netfor!(circ, pinmap[port[1]])
        row = findfirst(==(net), circ.nets)
        b = nb(circ)+i
        incid[row, b] = 1
        net = netfor!(circ, pinmap[port[2]])
        row = findfirst(==(net), circ.nets)
        incid[row, b] = -1
    end

    tv, ti = topomat!(incid)
    S = SparseMatrixCSC{Rational{BigInt}}([Mᵥ Mᵢ Mₓ Mₓ´ Mq;
        blockdiag(tv, ti) spzeros(nb(circ)+numports, 2nx(circ)+nq(circ))])
    ũ, M = gensolve(S, SparseMatrixCSC{Rational{BigInt}}([Mu u0; spzeros(nb(circ)+numports, nu(circ)+1)]))
    # [v' i' x' xd' q']' = ũ + My for arbitrary y
    # now drop all rows concerning only internal voltages and currents
    indices = vcat(consecranges([nb(circ), numports, nb(circ), numports+2nx(circ)+nq(circ)])[[2,4]]...)
    ũ = ũ[indices,:]
    M = M[indices,:]
    S̃ = SparseMatrixCSC(gensolve(M', spzeros(size(M,2), 0))[2]') # output matrices p as RHS?
    # S̃ spans nullspace of M', so that
    #    S̃*[v' i' x' xd' q']' = S̃*ũ + S̃*My = S̃*ũ
    # i.e. S̃ acts as a new system matrix
    M̃ᵥ, M̃ᵢ, M̃ₓ, M̃ₓ´, M̃q =
        matsplit(S̃, [size(S̃,1)], [numports, numports, nx(circ), nx(circ), nq(circ)])
    M̃u = S̃*ũ[:, 1:nu(circ)]
    ũ0 = S̃*ũ[:, end]
    # note that we need to flip the sign of M̃ᵢ to view the ports from the other side
    return Element(mv = M̃ᵥ, mi = -M̃ᵢ, mx = M̃ₓ, mxd = M̃ₓ´, mq = M̃q,
        mu = M̃u, u0 = ũ0,
        nonlinear_eq=nonlinear_eq_func(circ),
        ports=ports)
end

function ports_from_pinmap(pinmap)
    ks = keys(pinmap)
    refpin = first(ks)
    ks = Iterators.drop(ks, 1)
    return [refpin => pin for pin in ks]
end


text:
using Morsel

app = Morsel.app()

# Route with HTTP verbs
get(app, "/") do req, res
    "hello"
end

post(app, "/") do req, res
    # ...
end

# Route with URL vars
route(app, GET, "/hello/<id::Int>/<name::%[a-z]{3}[0-9]{2}>") do req, res
    req.params[:id] == 99 ? string("Go to hell ", req.params[:name]) : string("Hello ", req.params[:name])
end

route(app, GET, "/") do req, res
    # ...
end

route(app, GET | POST, "/") do req, res
    # ...
end

route(app, ALL, "/") do req, res
    # ...
end

function handleImageUpload(req, res)
    "Nice pix."
end

route(app, POST, "/upload", handleImageUpload)

# Route middleware
auth = Morsel.Midware() do req, res
    if !get(session(req), :authenticated, false)
        return req, redirect(res, "/login")
    end
    req, res
end

with(app, auth) do app
    get(app, "/private") do req, res
        # ...
    end
end

namespace(app, "/admin", auth) do app
    get(app, "/pages/<page_id::Int>") do req, res
        page = get_page(req.params[:page_id])
        render("viewName.ejl", page)
    end

    put(app, "/pages/<page_id::Int>") do req, res
        update_page(get_page(req.params[:page_id]), req.params)
        redirect(res, "/pages/", req.params[:page_id])
    end
end

route(app, GET, "/*") do req, res
    res.headers["Status"] = 404
    render("404.ejl")
end

start(app, 8000)


text:
# [Martinsson & Tropp, Algorithm 7]
function rangefinder(A::AbstractMatrix{T}, r::Int; q::Int=0) where {T <: Number}
    m, n = size(A)
    if q == 0
        Ω = 1/sqrt(n) * randn(n, r)
        Y = A * Ω
    else
        Y = 1/sqrt(m) * randn(m, r)
        for i in 1:q
            Y = Array(qr(Y).Q)
            Y = A*A'*Y
        end
    end

    return Array(qr(Y).Q)
end


text:
"""
    batch(lens₁, lens₂, ..., lensₙ) :: Lens

From ``n`` lenses, create a single lens that gets/sets ``n``-tuple in
such a way that the number of call to the constructor is minimized.
This is done by calling [`IndexBatchLens`](@ref) whenever possible.

# Examples
```jldoctest
julia> using Kaleido, Setfield

julia> lens = @batchlens begin
           _.a.b.c
           _.a.b.d
           _.a.e
       end;

julia> @assert lens ==
           IndexBatchLens(:a) ∘ MultiLens((
               (@lens _[1]) ∘ IndexBatchLens(:b, :e) ∘ MultiLens((
                   (@lens _[1]) ∘ IndexBatchLens(:c, :d),
                   (@lens _[2]) ∘ Kaleido.SingletonLens(),
               )) ∘ FlatLens(2, 1),
           )) ∘ FlatLens(3)

julia> obj = (a=(b=(c=1, d=2), e=3),);

julia> get(obj, lens)
(1, 2, 3)

julia> set(obj, lens, (10, 20, 30))
(a = (b = (c = 10, d = 20), e = 30),)
```
"""
batch

leftmost(l::Lens) = l
leftmost(l::ComposedLens) = leftmost(l.outer)
rightlens(::Lens) = IdentityLens()
rightlens(l::ComposedLens) = rightlens(l.outer) ∘ l.inner

startswithproperty(::T) where {T <: Lens} = startswithproperty(T)
startswithproperty(::Type{<:Lens}) = false
startswithproperty(::Type{<:PropertyLens}) = true
startswithproperty(::Type{<:ComposedLens{LO}}) where {LO} = startswithproperty(LO)

allstartswithproperty(::T) where {T} = allstartswithproperty(T)
allstartswithproperty(::Type{Tuple{}}) = true
allstartswithproperty(::Type{T}) where {N, T <: NTuple{N, Lens}} =
    startswithproperty(Base.tuple_type_head(T)) &&
    allstartswithproperty(Base.tuple_type_tail(T))

propname(::PropertyLens{name}) where name = name

batch(::IdentityLens) = SingletonLens()

batch(lenses::PropertyLens...) = IndexBatchLens(propname.(lenses)...)

function batch(lenses::Lens...)
    allstartswithproperty(lenses) || return MultiLens(lenses)
    partitions = partitionby(lenses, leftmost)
    propnames = _map(x -> propname(x.key), partitions)
    sublenses = _map(x -> rightlens.(x.values), partitions)
    indexlenses = ntuple(i -> (@lens _[i]), length(partitions))
    return IndexBatchLens(propnames...) ∘
        MultiLens(_compose.(indexlenses, _batch.(sublenses))) ∘
        FlatLens(_map(x -> length(x.values), partitions)...)
end
# TODO: Sort the lenses first and compute the permutation to recover
# the original order.

_batch(lenses::Tuple{Vararg{Lens}}) = batch(lenses...)

"""
    @batchlens begin
        lens_expression_1
        lens_expression_2
        ...
        lens_expression_n
    end

From ``n`` "lens expression", create a lens that gets/sets ``n``-tuple.
Each "lens expression" is an expression that is supported by
`Setfield.@lens` or such expression post-composed with other lenses
using `∘`.

See also [`batch`](@ref) which does all the heavy lifting of the
transformation of the lenses.

# Examples
```jldoctest
julia> using Kaleido, Setfield

julia> lens = @batchlens begin
           _.a.b.c
           _.a.b.d ∘ converting(fromfield = x -> parse(Int, x), tofield = string)
           _.a.e
       end;

julia> obj = (a = (b = (c = 1, d = "2"), e = 3),);

julia> get(obj, lens)
(1, 2, 3)

julia> set(obj, lens, (10, 20, 30))
(a = (b = (c = 10, d = "20"), e = 30),)
```
"""
macro batchlens(lenses_expression)
    if !(lenses_expression isa Expr && lenses_expression.head == :block)
        error("""
        Macro @batchlens needs a block of lens expressions.  Got:
        $lenses_expression
        """)
    end
    lnns = lenses_expression.args[1:2:end]
    exprs = lenses_expression.args[2:2:end]
    if all(isa.(lnns, LineNumberNode)) &&
            all(isa.(exprs, Expr)) &&
            length(lnns) == length(exprs)
        lens_exprs = make_lens_expr.(exprs, lnns)
    else
        lens_exprs = [
            make_lens_expr(x, __source__) for x in lenses_expression.args
            if x isa Expr
        ]
    end
    return esc(Expr(:call, batch, lens_exprs...))
end

islensexpr(ex::Symbol) = ex == :_
islensexpr(ex::Expr) = ex.head in (:., :ref) && islensexpr(ex.args[1])

function make_lens_expr(ex::Expr, lnn::LineNumberNode)
    if !(islensexpr(ex) || (ex.head == :call && ex.args[1] == :∘))
        error("Not a lens expression: $ex")
    end
    return _make_lens_expr(ex, lnn)
end

_make_lens_expr(ex::Symbol, lnn::LineNumberNode) = ex
function _make_lens_expr(ex::Expr, lnn::LineNumberNode)
    if islensexpr(ex)
        # Make `:(Setfield.@lens $ex)` with a proper `LineNumberNode`:
        atlens = Expr(:., Setfield, QuoteNode(Symbol("@lens")))
        return Expr(:macrocall, atlens, lnn, ex)
    elseif ex.head == :call && ex.args[1] == :∘
        return Expr(:call, :∘, _make_lens_expr.(ex.args[2:end], Ref(lnn))...)
    else
        return ex
    end
end


text:
using NodeJS

# Directory into which we'll install Node.js packages via npm
const NPM_DIR = joinpath(@__DIR__, "node_modules")

if !isdir(NPM_DIR)
    run(setenv(`$(npm_cmd()) install to-markdown`, dir=@__DIR__))
end

if !isdir(joinpath(NPM_DIR, "to-markdown"))
    error("installation of Node.js module to-markdown did not complete successfully")
end

if !isfile(joinpath(@__DIR__, "deps.jl"))
    open(joinpath(@__DIR__, "deps.jl"), "w") do io
        println(io, "const NPM_DIR = \"", NPM_DIR, "\"")
    end
end


text:
# Initialising
using PyPlot
using FastGaussQuadrature # For besselroots.jl

include("gausslaguerreTest.jl")

# Test for accuracy of the weights
nt = 200
alphat = 0.7
(xd,wd) = gausslaguerreTest(nt,alphat,"rec")
for hT = 1:4
	(xex,wex) = gausslaguerreTest(nt, alphat, "exp", 1.0, 1, 2*hT)
	semilogy(xex, abs(wd-wex)./wd, label="$hT terms")
end
legend(loc="upper center")
xlabel(L"$x_k$")
ylabel(L"Relative error on $w_k$")

# Timings
ns = Array{Int64,1}(round(10.^((6:24)/3)))
lns = length(ns)
rests = zeros(5,lns)
alphat = 0.0
nt = 100

figure
maxTim = 10 # Stop if the previous gausslaguerreTest took over 10 seconds.

(x,w) = gausslaguerreTest(nt,alphat,"GW")
for ni = 1:lns
	rests[2,ni] = @elapsed gausslaguerreTest(ns[ni],alphat,"GW")
	print("\n", ni, " GW ", rests[2,ni])
	if (rests[2,ni] > maxTim)
		break
	end
end
loglog(ns, rests[2,:], label="Golub-Welsch")

(x,w) = gausslaguerreTest(nt,alphat,"rec")
for ni = 1:lns
	rests[1,ni] = @elapsed gausslaguerreTest(ns[ni],alphat,"rec")
	print("\n", ni, " rec ", rests[1,ni])
	if (rests[1,ni] > maxTim)
		break
	end
end
loglog(ns, rests[1,:], label="Newton on recurrence")

(x,w) = gausslaguerreTest(nt,alphat,"GLR")
for ni = 1:lns
	rests[3,ni] = @elapsed gausslaguerreTest(ns[ni],alphat,"GLR")
	print("\n", ni, " GLR ", rests[3,ni])
	if (rests[3,ni] > maxTim)
		break
	end
end
loglog(ns, rests[3,:], label="Glaser-Liu-Rokhlin")

(x,w) = gausslaguerreTest(nt,alphat,"gen")
for ni = 1:lns
	rests[4,ni] = @elapsed gausslaguerreTest(ns[ni],alphat,"gen")
	print("\n", ni, " gen ", rests[4,ni])
	if (rests[4,ni] > maxTim)
		break
	end
end
loglog(ns, rests[4,:], label="Newton on a.e.")

(x,w) = gausslaguerreTest(nt,alphat,"exp")
for ni = 1:lns
	rests[5,ni] = @elapsed gausslaguerreTest(ns[ni],alphat,"exp")
	print("\n", ni, " exp ", rests[5,ni])
	if(rests[5,ni] > maxTim)
		break
	end
end
loglog(ns, rests[5,:], label="Asymptotic expansions")

loglog(10.^(3:8), [0.216; 0.311; 0.391; 0.510; 2.71; 23.9], label="Bremer (Fortran)")

legend(loc="best")
xlabel(L"$n$")
ylabel("Time (s)")


text:
"""
    $(FUNCTIONNAME)(λ_bnds)

    Loads in the PROSPECT-PRO database of pigments (and other) absorption cross section in leaves, returns a [`LeafOpticalProperties`](@ref) type struct with spectral units attached.
# Arguments
    - `λ_bnds` an array (with or without a spectral grid unit) that defines the upper and lower limits over which to average the absorption cross sections

# Examples
```julia-repl
julia> using Unitful                                                               # Need to include Unitful package
julia> opti = createLeafOpticalStruct((400.0:5:2400)*u"nm");                       # in nm
julia> opti = createLeafOpticalStruct((0.4:0.1:2.4)*u"μm");                        # in μm
julia> opti = CanopyOptics.createLeafOpticalStruct((10000.0:100:25000.0)u"1/cm");  # in wavenumber (cm⁻¹)
```
"""
function createLeafOpticalStruct(λ_bnds)
    FT = eltype(ustrip(λ_bnds[1]))
    # Reference input grid converted to nm:
    λ_ref = unit2nm(λ_bnds)

    KS = readdlm(OPTI_2021, '\t',FT)
    N  = length(λ_bnds)-1
    λ_out, nᵣ, Kcab, Kcar, Kant, Kb, Kw, Km, Kp, Kcbc = [zeros(FT,N) for i = 1:N];
    vars = (λ_out, nᵣ, Kcab, Kcar, Kant, Kb, Kw, Km, Kp, Kcbc)
    λ     = KS[:,1]*u"nm";
    @inbounds for i=1:N
        start = min(λ_ref[i],λ_ref[i+1])
        stop  = max(λ_ref[i],λ_ref[i+1])
        ind_all = findall((λ .≥ start) .& (λ .< stop) )
        isempty(ind_all) ? (@warn "Some λ_bnds bins are empty, better coarsen the input grid $(λ_bnds[i]) -> $(λ_bnds[i+1])" ) : nothing
        for j in eachindex(vars)
            vars[j][i] =  mean(view(KS,ind_all,j));
        end
    end
    # Get output unit:
    ou = unit(λ_bnds[1])
    (typeof(ou) <: Unitful.FreeUnits{(), NoDims, nothing}) ? out_unit = u"nm"  : out_unit = ou
    return LeafOpticalProperties(uconvSpectral(λ_out*u"nm",out_unit),
                                nᵣ,
                                Kcab,
                                Kcar,
                                Kant,
                                Kb,
                                Kw,
                                Km,
                                Kp,
                                Kcbc)
end

"""
    $(FUNCTIONNAME)()
    As in $(FUNCTIONNAME)(λ_bnds) but reads in the in Prospect-PRO database at original resolution (400-2500nm in 1nm steps)
"""
function createLeafOpticalStruct()
    λ_bnds = (399.5:2500.5)u"nm"
    createLeafOpticalStruct(λ_bnds)
end

"Convert an input array with Spectral units (or none) to a grid in nm"
function unit2nm(λ_bnds)
    try
        return uconvSpectral(λ_bnds, u"nm");
    catch e # Assume without unit is in "nm" already
        @warn "Assuming unit of [nm] here for optical struct"
        return λ_bnds*u"nm"
    end
end

"Avoiding annoying issues with uconvert for Spectra() (can't do nm->nm!)"
# From unit of `in` to `out_unit` (just attaching out_unit if unit is missing)
function uconvSpectral(in,out_unit)
    unit(in[1]) == out_unit ? in : uconvert.(out_unit, in, Spectral())
end


text:
"""
    mglmxls(glmout::Vector{DataFrames.DataFrameRegressionModel}, workbook::PyObject, worksheet::AbstractString; labels::Union{Nothing,Dict}=nothing,mtitle::Union{Vector,Nothing}=nothing,eform=false,ci=true, row = 0, col =0)

Outputs multiple GLM regression tables side by side to an excel spreadsheet.
To use this function, `PyCall` is required with a working version python and
a python package called `xlsxwriter` installed. If a label is found for a variable
or a value of a variable in a `labels`, the label will be output. Options are:

- `glmout`: a vector of GLM regression models
- `workbook`: a returned value from xlsxwriter.Workbook() function (see an example below)
- `worksheet`: a string for the worksheet name
- `labels`: an option to specify a `label` dictionary (see an example below)
- `mtitle`: header label for GLM models. If not specified, the dependent variable name will be used.
- `eform`: use `eform = true` to get exponentiated estimates, standard errors, or 95% confidence intervals
- `ci`: use `ci = true` (default) to get 95% confidence intervals. `ci = false` will produce standard errors and Z values instead.
- `row`: specify the row of the workbook to start the output table (default = 0 (for row 1))
- `col`: specify the column of the workbook to start the output table (default = 0 (for column A))

# Example 1
This example is useful when one wants to append a worksheet to an existing workbook.
It is responsibility of the user to open a workbook before the function call and close it
to actually create the physical file by close the workbook.

```
julia> using PyCall

julia> xlsxwriter = pyimport("xlsxwriter")

julia> wb = xlsxwriter.Workbook("test_workbook.xlsx")
PyObject <xlsxwriter.workbook.Workbook object at 0x000000002A628E80>

julia> mglmxls(olsmodels,wb,"OLS1",labels = label)

julia> bivairatexls(df,:incomecat,[:age,:race,:male,:bmicat],wb,"Bivariate",labels = label)

Julia> wb.close()
```

# Example 2
Alternatively, one can create a spreadsheet file directly. `PyCall` or `pyimport()`
does not need to be called before the function.

```
julia> mglmxls(olsmodels,"test_workbook.xlsx","OLS1",labels = label)
```

"""
function mglmxls(glmout,
    wbook::PyObject,
    wsheet::AbstractString;
    mtitle::Union{Vector,Nothing} = nothing,
    labels::Union{Nothing,Label} = nothing,
    eform::Bool = false,
    ci = true,
    row = 0,
    col = 0)

    num_models = length(glmout)
    otype = Vector(undef,num_models)
    if isa(glmout[1].model,GeneralizedLinearModel)
	linkfun = Vector{Link}(undef,num_models)
    	distrib = Vector{UnivariateDistribution}(undef,num_models)
    end

    for i=1:num_models

        # assuming that all models have the same family and link function
        if isa(glmout[i].model,GeneralizedLinearModel)
            distrib[i] = glmout[i].model.rr.d
            linkfun[i] = Link(glmout[i].model.rr)
        end

        otype[i] = "Estimate"
        if eform == true
            otype[i] = Stella.coeflab(distrib[i],linkfun[i])
        end
    end

    if mtitle == nothing
        mtitle = Vector{String}(undef,num_models)

        # assign dependent variables
        for i=1:num_models
            ysim = glmout[i].mf.f.lhs.sym #terms.eterms[1]
            mtitle[i] = labels != nothing && haskey(labels.var, ysim) ? varlab(labels,ysim) : string(ysim)
        end
    end

    # create a worksheet
    t = wbook.add_worksheet(wsheet)

    # attach formats to the workbook
    formats = ExcelTables.attach_formats(wbook)

    # starting location in the worksheet
    r = row
    c = col

    # set column widths
    t.set_column(c,c,40)
    t.set_column(c+1,c+4*num_models,7)

    t.merge_range(r,c,r+1,c,"Variable",formats[:heading])
    for i=1:num_models
        if ci == true
            t.merge_range(r,c+1,r,c+4,mtitle[i],formats[:heading])
            t.merge_range(r+1,c+1,r+1,c+3,string(otype[i]," (95% CI)"),formats[:heading])
            t.write_string(r+1,c+4,"P-Value",formats[:heading])
        else
            t.write_string(r,c+1,otype[i],formats[:heading])
            t.write_string(r,c+2,"SE",formats[:heading])
            t.write_string(r,c+3,"Z Value",formats[:heading])
            t.write_string(r,c+4,"P-Value",formats[:heading])
        end
        c += 4
    end

    #------------------------------------------------------------------
    r += 2
    c = col

    # collate variables
    covariates = Vector{String}()
    tdata = Vector(undef,num_models)
    tconfint = Vector(undef,num_models)

    for i=1:num_models
        tdata[i] = coeftable(glmout[i])
        tconfint[i] = confint(glmout[i])

        for nm in tdata[i].rownms
            if in(nm, covariates) == false
                push!(covariates,nm)
            end
        end
    end

    # go through each variable and construct variable name and value label arrays
    nrows = length(covariates)
    varname = Vector{String}(undef,nrows)
    vals = Vector{String}(undef,nrows)
    nlev = zeros(Int,nrows)
    npred = [dof(m) for m in glmout]

    for i = 1:nrows
    	# variable name
        # parse varname to separate variable name from value
        if occursin(" & ",covariates[i])
            varname[i] = covariates[i]
            vals[i] = ""
        elseif occursin(":",covariates[i])
            (varname[i],vals[i]) = split(covariates[i],": ")
        else
            varname[i] = covariates[i]
            vals[i] = ""
        end

        # use labels if exist
        if labels != nothing
            sv = Symbol(varname[i])
            if vals[i] != ""
                valn = vals[i] == "true" ? 1 : parse(Int,vals[i])
                vals[i] = vallab(labels,sv,valn)
            end
            varname[i] = varlab(labels,sv)
        end
    end
    for i = 1:nrows
        nlev[i] = countlev(varname[i],varname)
    end

    # write table
    lastvarname = ""

    for i = 1:nrows
        if varname[i] != lastvarname
            # output cell boundaries only and go to the next line
            if nlev[i] > 1

                # variable name
                t.write_string(r,c,varname[i],formats[:model_name])

                for k = 1:num_models
                    if ci == true
                        t.write(r,c+1,"",formats[:empty_right])
                        t.write(r,c+2,"",formats[:empty_both])
                        t.write(r,c+3,"",formats[:empty_left])
                        t.write(r,c+4,"",formats[:p_fmt])
                    else
                        t.write(r,c+1,"",formats[:empty_border])
                        t.write(r,c+2,"",formats[:empty_border])
                        t.write(r,c+3,"",formats[:empty_border])
                        t.write(r,c+4,"",formats[:p_fmt])
                    end
                    c += 4
                end
                c = col
                r += 1
                t.write_string(r,c,vals[i],formats[:varname_1indent])

            else
                if vals[i] != "" && vals[i] != "Yes"
                    t.write_string(r,c,string(varname[i],": ",vals[i]),formats[:model_name])
                else
                    t.write_string(r,c,varname[i],formats[:model_name])
                end
            end
        else
            t.write_string(r,c,vals[i],formats[:varname_1indent])
        end

        for j=1:num_models

            # find the index number for each coeftable row
            ri = findfirst(x->x == covariates[i],tdata[j].rownms)

            if ri == nothing
                # this variable is not in the model
                # print empty cells and then move onto the next model

                t.write_string(r,c+1,"",formats[:or_fmt])
                t.write_string(r,c+2,"",formats[:cilb_fmt])
                t.write_string(r,c+3,"",formats[:ciub_fmt])
                t.write_string(r,c+4,"",formats[:p_fmt])

                c += 4
                continue
            end

    	    # estimates
            if eform == true
        	t.write(r,c+1,ri <= npred[j] ? exp(tdata[j].cols[1][ri]) : "",formats[:or_fmt])
            else
                t.write(r,c+1,ri <= npred[j] ? tdata[j].cols[1][ri] : "",formats[:or_fmt])
            end

            if ci == true

                if eform == true
                	# 95% CI Lower
                	t.write(r,c+2,ri <= npred[j] ? exp(tconfint[j][ri,1]) : "",formats[:cilb_fmt])

                	# 95% CI Upper
                	t.write(r,c+3,ri <= npred[j] ? exp(tconfint[j][ri,2]) : "",formats[:ciub_fmt])
                else
                    # 95% CI Lower
                	t.write(r,c+2,ri <= npred[j] ? tconfint[j][ri,1] : "",formats[:cilb_fmt])

                	# 95% CI Upper
                	t.write(r,c+3,ri <= npred[j] ? tconfint[j][ri,2] : "",formats[:ciub_fmt])
                end
            else
                # SE
                if eform == true
            	    t.write(r,c+2,ri <= npred[j] ? exp(tdata[j].cols[1][ri])*tdata[j].cols[2][ri] : "",formats[:or_fmt])
                else
                    t.write(r,c+2,ri <= npred[j] ? tdata[j].cols[1][ri] : "",formats[:or_fmt])
                end

                # Z value
                t.write(r,c+3,ri <= npred[j] ? tdata[j].cols[3][ri] : "",formats[:or_fmt])

            end

            # P-Value
	        t.write(r,c+4,ri <= npred[j] ? (tdata[j].cols[4][ri] < 0.001 ? "< 0.001" : tdata[j].cols[4][ri]) : "" ,formats[:p_fmt])

            c += 4

        end

        lastvarname = varname[i]

        # update row
        r += 1
        c = col
    end

    # Write model characteristics and goodness of fit statistics
    c = col
    row2 = r
    for i=1:num_models

        # N
        t.write(r,c,"N",formats[:model_name])
        t.merge_range(r,c+1,r,c+4,nobs(glmout[i]),formats[:n_fmt_center])

        # degress of freedom
        r += 1
        t.write(r,c,"DF",formats[:model_name])
        t.merge_range(r,c+1,r,c+4,dof(glmout[i]),formats[:n_fmt_center])

        # R² or pseudo R²
        r += 1
        if isa(glmout[i].model,LinearModel)
            t.write(r,c,"R²",formats[:model_name])
            t.merge_range(r,c+1,r,c+4,r2(glmout[i]),formats[:p_fmt_center])
            t.write(r+1,c,"Adjusted R²",formats[:model_name])
            t.merge_range(r+1,c+1,r+1,c+4,adjr2(glmout[i]),formats[:p_fmt_center])

            r += 2

	else
	    if isa(linkfun[i],LogitLink)
		    t.write(r,c,"Pseudo R² (MacFadden)",formats[:model_name])
		    t.merge_range(r,c+1,r,c+4,macfadden(glmout[i]),formats[:p_fmt_center])
		    t.write(r+1,c,"Pseudo R² (Nagelkerke)",formats[:model_name])
		    t.merge_range(r+1,c+1,r+1,c+4,nagelkerke(glmout[i]),formats[:p_fmt_center])

		    # -2 log-likelihood
		    t.write(r+2,c,"-2 Log-Likelihood",formats[:model_name])
		    t.merge_range(r+2,c+1,r+2,c+4,deviance(glmout[i]),formats[:p_fmt_center])

		    # Hosmer-Lemeshow GOF test
		    t.write(r+3,c,"Hosmer-Lemeshow Chisq Test (df), p-value",formats[:model_name])
		    hl = hltest(glmout[i])
		    t.merge_range(r+3,c+1,r+3,c+4,string(round(hl[1],digits=4)," (",hl[2],"); p = ",round(hl[3],digits=4)),formats[:p_fmt_center])

		    # ROC (c-statistic)
		    t.write(r+4,c,"Area under the ROC Curve",formats[:model_name])
		    _roc = auc(roc(glmout[i].model.rr.y,predict(glmout[i])))
		    t.merge_range(r+4,c+1,r+4,c+4,round(_roc,digits=4),formats[:p_fmt_center])

		    r += 5
	    end
	end

        # AIC & BIC
        t.write(r,c,"AIC",formats[:model_name])
        t.merge_range(r,c+1,r,c+4,aic(glmout[i]),formats[:p_fmt_center])

        r += 1
        t.write(r,c,"BIC",formats[:model_name])
        t.merge_range(r,c+1,r,c+4,bic(glmout[i]),formats[:p_fmt_center])

        r = row2
        c += 4
    end
end
function mglmxls(glmout,
    wbook::AbstractString,
    wsheet::AbstractString;
    mtitle::Union{String,Nothing} = nothing,
    labels::Union{Nothing,Label} = nothing,
    eform::Bool = false,
    ci = true,
    row = 0,
    col = 0)

    xlsxwriter = pyimport("xlsxwriter")

    mglmxls(glmout,xlsxwriter.Workbook(wbook),wsheet,labels=labels,mtitle=mtitle,eform=eform,ci=ci,row=row,col=col)
end


text:
module Err
export PatternUnsolvedException, InternalException, SyntaxError

struct PatternUnsolvedException <: Exception
    msg :: String
    PatternUnsolvedException(arg) =
        if isa(arg, String)
            new(arg)
        else
            new("Non-exhaustive pattern found for target `$(string(arg))`.")
        end
end

struct InternalException <: Exception
    msg :: String
end

struct SyntaxError <: Exception
    msg :: String
end

struct DataTypeUsageError <: Exception
    msg :: String
end

end


text:
"""
$(TYPEDEF)

A system of stochastic differential equations.

# Fields
$(FIELDS)

# Example

```
using ModelingToolkit

@parameters t σ ρ β
@variables x(t) y(t) z(t)
@derivatives D'~t

eqs = [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]

noiseeqs = [0.1*x,
            0.1*y,
            0.1*z]

de = SDESystem(eqs,noiseeqs,t,[x,y,z],[σ,ρ,β])
```
"""
struct SDESystem <: AbstractODESystem
    """The expressions defining the drift term."""
    eqs::Vector{Equation}
    """The expressions defining the diffusion term."""
    noiseeqs::AbstractArray{Operation}
    """Independent variable."""
    iv::Variable
    """Dependent (state) variables."""
    states::Vector{Variable}
    """Parameter variables."""
    ps::Vector{Variable}
    """
    Time-derivative matrix. Note: this field will not be defined until
    [`calculate_tgrad`](@ref) is called on the system.
    """
    tgrad::RefValue{Vector{Expression}}
    """
    Jacobian matrix. Note: this field will not be defined until
    [`calculate_jacobian`](@ref) is called on the system.
    """
    jac::RefValue{Matrix{Expression}}
    """
    `Wfact` matrix. Note: this field will not be defined until
    [`generate_factorized_W`](@ref) is called on the system.
    """
    Wfact::RefValue{Matrix{Expression}}
    """
    `Wfact_t` matrix. Note: this field will not be defined until
    [`generate_factorized_W`](@ref) is called on the system.
    """
    Wfact_t::RefValue{Matrix{Expression}}
    """
    Name: the name of the system
    """
    name::Symbol
    """
    Systems: the internal systems
    """
    systems::Vector{SDESystem}
end

function SDESystem(deqs::AbstractVector{<:Equation}, neqs, iv, dvs, ps;
                   systems = SDESystem[],
                   name = gensym(:SDESystem))
    iv′ = convert(Variable,iv)
    dvs′ = convert.(Variable,dvs)
    ps′ = convert.(Variable,ps)
    tgrad = RefValue(Vector{Expression}(undef, 0))
    jac = RefValue(Matrix{Expression}(undef, 0, 0))
    Wfact   = RefValue(Matrix{Expression}(undef, 0, 0))
    Wfact_t = RefValue(Matrix{Expression}(undef, 0, 0))
    SDESystem(deqs, neqs, iv′, dvs′, ps′, tgrad, jac, Wfact, Wfact_t, name, systems)
end

function generate_diffusion_function(sys::SDESystem, dvs = sys.states, ps = sys.ps; kwargs...)
    dvs′ = convert.(Variable,dvs)
    ps′ = convert.(Variable,ps)
    return build_function(sys.noiseeqs, dvs′, ps′, sys.iv;
                          conv = ODEToExpr(sys),kwargs...)
end

"""
```julia
function DiffEqBase.SDEFunction{iip}(sys::SDESystem, dvs = sys.states, ps = sys.ps;
                                     version = nothing, tgrad=false, sparse = false,
                                     jac = false, Wfact = false, kwargs...) where {iip}
```

Create an `SDEFunction` from the [`SDESystem`](@ref). The arguments `dvs` and `ps`
are used to set the order of the dependent variable and parameter vectors,
respectively.
"""
function DiffEqBase.SDEFunction{iip}(sys::SDESystem, dvs = sys.states, ps = sys.ps,
                                     u0 = nothing;
                                     version = nothing, tgrad=false, sparse = false,
                                     jac = false, Wfact = false, kwargs...) where {iip}
    f_oop,f_iip = generate_function(sys, dvs, ps; expression=Val{false}, kwargs...)
    g_oop,g_iip = generate_diffusion_function(sys, dvs, ps; expression=Val{false}, kwargs...)

    f(u,p,t) = f_oop(u,p,t)
    f(du,u,p,t) = f_iip(du,u,p,t)
    g(u,p,t) = g_oop(u,p,t)
    g(du,u,p,t) = g_iip(du,u,p,t)

    if tgrad
        tgrad_oop,tgrad_iip = generate_tgrad(sys, dvs, ps; expression=Val{false}, kwargs...)
        _tgrad(u,p,t) = tgrad_oop(u,p,t)
        _tgrad(J,u,p,t) = tgrad_iip(J,u,p,t)
    else
        _tgrad = nothing
    end

    if jac
        jac_oop,jac_iip = generate_jacobian(sys, dvs, ps; expression=Val{false}, sparse=sparse, kwargs...)
        _jac(u,p,t) = jac_oop(u,p,t)
        _jac(J,u,p,t) = jac_iip(J,u,p,t)
    else
        _jac = nothing
    end

    if Wfact
        tmp_Wfact,tmp_Wfact_t = generate_factorized_W(sys, dvs, ps, true; expression=Val{false}, kwargs...)
        Wfact_oop, Wfact_iip = tmp_Wfact
        Wfact_oop_t, Wfact_iip_t = tmp_Wfact_t
        _Wfact(u,p,dtgamma,t) = Wfact_oop(u,p,dtgamma,t)
        _Wfact(W,u,p,dtgamma,t) = Wfact_iip(W,u,p,dtgamma,t)
        _Wfact_t(u,p,dtgamma,t) = Wfact_oop_t(u,p,dtgamma,t)
        _Wfact_t(W,u,p,dtgamma,t) = Wfact_iip_t(W,u,p,dtgamma,t)
    else
        _Wfact,_Wfact_t = nothing,nothing
    end

    M = calculate_massmatrix(sys)
    _M = (u0 === nothing || M == I) ? M : ArrayInterface.restructure(u0 .* u0',M)

    SDEFunction{iip}(f,g,jac=_jac,
                      tgrad = _tgrad,
                      Wfact = _Wfact,
                      Wfact_t = _Wfact_t,
                      mass_matrix = _M,
                      syms = Symbol.(sys.states))
end

function DiffEqBase.SDEFunction(sys::SDESystem, args...; kwargs...)
    SDEFunction{true}(sys, args...; kwargs...)
end

function rename(sys::SDESystem,name)
    ODESystem(sys.eqs, sys.noiseeqs, sys.iv, sys.states, sys.ps, sys.tgrad, sys.jac, sys.Wfact, sys.Wfact_t, name, sys.systems)
end

"""
```julia
function DiffEqBase.SDEProblem{iip}(sys::SDESystem,u0map,tspan,p=parammap;
                                    version = nothing, tgrad=false,
                                    jac = false, Wfact = false,
                                    checkbounds = false, sparse = false,
                                    sparsenoise = sparse,
                                    linenumbers = true, parallel=SerialForm(),
                                    kwargs...)
```

Generates an SDEProblem from an SDESystem and allows for automatically
symbolically calculating numerical enhancements.
"""
function DiffEqBase.SDEProblem{iip}(sys::SDESystem,u0map,tspan,parammap=DiffEqBase.NullParameters();
                                    version = nothing, tgrad=false,
                                    jac = false, Wfact = false,
                                    checkbounds = false, sparse = false,
                                    sparsenoise = sparse,
                                    linenumbers = true, parallel=SerialForm(),
                                    kwargs...) where iip

    dvs = states(sys)
    ps = parameters(sys)
    u0 = varmap_to_vars(u0map,dvs)
    p = varmap_to_vars(parammap,ps)
    f = SDEFunction{iip}(sys,dvs,ps,u0;tgrad=tgrad,jac=jac,Wfact=Wfact,
                         checkbounds=checkbounds,
                         linenumbers=linenumbers,parallel=parallel,
                         sparse=sparse)
    if typeof(sys.noiseeqs) <: AbstractVector
        noise_rate_prototype = nothing
    elseif sparsenoise
        I,J,V = findnz(SparseArrays.sparse(sys.noiseeqs))
        noise_rate_prototype = SparseArrays.sparse(I,J,zero(eltype(u0)))
    else
        noise_rate_prototype = zeros(eltype(u0),size(sys.noiseeqs))
    end

    SDEProblem{iip}(f,f.g,u0,tspan,p;noise_rate_prototype=noise_rate_prototype,kwargs...)
end

function DiffEqBase.SDEProblem(sys::SDESystem, args...; kwargs...)
    SDEProblem{true}(sys, args...; kwargs...)
end


text:
@testset "Strings" begin
    @test LO(sqrt(d64"2")) < LO(nextfloat(sqrt(d64"2")))
end

@testset "Adjacent values" begin
    @test nextfloat(Double64(10))  > Double64(10)
    @test nextfloat(Double32(-10)) > Double32(-10)
    @test prevfloat(Double64(10))  < Double64(10)
    @test prevfloat(Double32(-10)) < Double32(-10)
    @test nextfloat(Double64(10),3)  > nextfloat(Double64(10),2)
    @test nextfloat(Double32(-10),3) > nextfloat(Double32(-10),2)
    @test prevfloat(Double64(10),3)  < prevfloat(Double64(10),2)
    @test prevfloat(Double32(-10),3) < prevfloat(Double32(-10),2)
end

@testset "power functions"  begin
    @test Double64(10.0)^0 == Double64(1.0)
    @test Double64(10.0)^1 == Double64(10.0)
    @test_throws DomainError Double64(0.0)^0
end

f64a = sqrt(2.0)
f64b = sqrt(2.0)/2
d64a = sqrt(Double64(2.0))
d64b = sqrt(Double64(2.0))/2

@testset "log exp functions" begin
   @test isapprox(log(f64a), log(d64a))
   @test isapprox(exp(f64a), exp(d64a))
end

@testset "trig functions" begin
   @test isapprox(sin(f64a), sin(d64a))
   @test isapprox(cos(f64a), cos(d64a))
   @test isapprox(tan(f64a), tan(d64a))
   @test isapprox(csc(f64a), csc(d64a))
   @test isapprox(sec(f64a), sec(d64a))
   @test isapprox(cot(f64a), cot(d64a))
   @test isapprox(asin(f64b), asin(d64b))
   @test isapprox(acos(f64b), acos(d64b))
   @test isapprox(atan(f64a), atan(d64a))
   @test isapprox(acsc(f64a), acsc(d64a))
   @test isapprox(asec(f64a), asec(d64a))
   @test isapprox(acot(f64a), acot(d64a))
end

@testset "hyperbolic trig functions" begin
   @test isapprox(sinh(f64a), sinh(d64a))
   @test isapprox(cosh(f64a), cosh(d64a))
   @test isapprox(tanh(f64a), tanh(d64a))
   @test isapprox(csch(f64a), csch(d64a))
   @test isapprox(sech(f64a), sech(d64a))
   @test isapprox(coth(f64a), coth(d64a))
   @test isapprox(asinh(f64b), asinh(d64b))
   @test isapprox(acosh(f64a), acosh(d64a))
   @test isapprox(atanh(f64b), atanh(d64b))
   @test isapprox(acsch(f64b), acsch(d64b))
   @test isapprox(asech(f64b), asech(d64b))
   @test isapprox(acoth(f64a), acoth(d64a))
end


text:
"""
Determine the null and alternate hypotheses under t-tests given:

* `ns`: Per-group sample size(s)
* `δ`: Deviation from the null
* `σs`: Per-group standard deviations
"""
function hypotheses(::Type{OneSampleTTest}, n::Real, δ::Real, σ::Real)
    df = n - 1
    ncp = sqrt(n) * δ / σ
    null = TDist(df)
    alt = NoncentralT(df, ncp)
    return null, alt
end

"""
Determine the null and alternate hypotheses under t-tests given:

* `ns`: Per-group sample size(s)
* `δ`: Deviation from the null
* `σs`: Per-group standard deviations
"""
function hypotheses(::Type{EqualVarianceTTest}, n::Real, δ::Real, σ::Real)
    df = 2 * n - 2
    ncp = sqrt(n / 2) * δ / σ
    null = TDist(df)
    alt = NoncentralT(df, ncp)
    return null, alt
end

"""
Determine the null and alternate hypotheses under t-tests given:

* `ns`: Per-group sample size(s)
* `δ`: Deviation from the null
* `σs`: Per-group standard deviations
"""
function hypotheses(::Type{EqualVarianceTTest}, ns::Tuple, δ::Real, σ::Real)
    n1, n2 = ns[1], ns[2]
    df = n1 + n2 - 2
    ncp = (1 / sqrt(1 / n1 + 1 / n2)) * δ / σ
    null = TDist(df)
    alt = NoncentralT(df, ncp)
    return null, alt
end

# @doc """
# Determine the values of t that lead to overestimates of the magnitude of δ.
# """ ->
thresholds(alt::NoncentralT) = -alt.λ, alt.λ


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

# make sure that typeinf is executed before turning on typeinf_ext
# this ensures that typeinf_ext doesn't recurse before it can add the item to the workq
# especially try to make sure any recursive and leaf functions have concrete signatures,
# since we won't be able to specialize & infer them at runtime

time() = ccall(:jl_clock_now, Float64, ())

let
    world = get_world_counter()
    interp = NativeInterpreter(world)

    analyze_escapes_tt = Tuple{typeof(analyze_escapes), IRCode, Int, Bool, typeof(null_escape_cache)}
    fs = Any[
        # we first create caches for the optimizer, because they contain many loop constructions
        # and they're better to not run in interpreter even during bootstrapping
        #=analyze_escapes_tt,=# run_passes,
        # then we create caches for inference entries
        typeinf_ext, typeinf, typeinf_edge,
    ]
    # tfuncs can't be inferred from the inference entries above, so here we infer them manually
    for x in T_FFUNC_VAL
        push!(fs, x[3])
    end
    for i = 1:length(T_IFUNC)
        if isassigned(T_IFUNC, i)
            x = T_IFUNC[i]
            push!(fs, x[3])
        else
            println(stderr, "WARNING: tfunc missing for ", reinterpret(IntrinsicFunction, Int32(i)))
        end
    end
    starttime = time()
    for f in fs
        if isa(f, DataType) && f.name === typename(Tuple)
            tt = f
        else
            tt = Tuple{typeof(f), Vararg{Any}}
        end
        for m in _methods_by_ftype(tt, 10, typemax(UInt))
            # remove any TypeVars from the intersection
            typ = Any[m.spec_types.parameters...]
            for i = 1:length(typ)
                typ[i] = unwraptv(typ[i])
            end
            typeinf_type(interp, m.method, Tuple{typ...}, m.sparams)
        end
    end
    endtime = time()
    println("Core.Compiler ──── ", sub_float(endtime,starttime), " seconds")
end

max_stars_repo_name: JuliaLabs/julia

max_issues_repo_name: JuliaLabs/julia

max_forks_repo_name: JuliaLabs/julia