text:
	beta = -0.14
	x0 = nothing
	time_limit = 10.0
	d_pseudo = nothing
	delta_treshold = 1e-10
	verbose = true

max_degree = 5
num_bootstrap = 100
alpha = 0.20

# Get a regression problem with 2 covars, expand it and use a model with 1 active covar.
rp, y, errors = random_regression_problem(2, 1)
# Use DM to reconstruct the coefficients:
m, k, xs, its, aics, degs = min_aic_difference_map_reconstruction(rp.d.data, y)
analyze_regression(rp, m)

# Test on random data
N = 30
M = 2
d = randn(M, N)

y = y1 = 4.5*d[1,:] + 17.0*d[2,:] + 8.2*d[3,:]
y2 = 1.2*d[1,:] + 2.0*d[2,:]
y3 = 1.5*d[1,:]

xbest1, k1, xs1, its1, aics1, degs1 = min_aic_difference_map_reconstruction(d, y1)
xbest2, k2, xs2, its2, aics2, degs2 = min_aic_difference_map_reconstruction(d, y2)
xbest3, k3, xs3, its3, aics3, degs3 = min_aic_difference_map_reconstruction(d, y3)

x4, i4, aic4 = difference_map_sparse_reconstruction(d, y, 4; time_limit = 1.0)
x3, i3, aic3 = difference_map_sparse_reconstruction(d, y, 3; time_limit = 1.0)

# The pattern seems to be that when we come down below the actual
# degree the reconstruction does not converge:
x2, i2, aic2 = difference_map_sparse_reconstruction(d, y, 2; time_limit = 1.0)
x1, i1, aic1 = difference_map_sparse_reconstruction(d, y, 1; time_limit = 1.0)

# There also seems to be a pattern that we get faster convergence if starting
# from an earlier fix point:
x4b, i4b = difference_map_sparse_reconstruction(d, y, 3; x0 = x5)
x3b, i3b = difference_map_sparse_reconstruction(d, y, 3; x0 = x4)
x2b, i2b = difference_map_sparse_reconstruction(d, y, 2; x0 = x3, time_limit = 1.0)
x1b, i1b = difference_map_sparse_reconstruction(d, y, 1; x0 = x2, time_limit = 1.0)

xr5 = repeated_difference_map_sparse_reconstruction(d, y, 5; time_limit = 1.0)
xr4 = repeated_difference_map_sparse_reconstruction(d, y, 4; time_limit = 1.0)
xr3 = repeated_difference_map_sparse_reconstruction(d, y, 3; time_limit = 1.0)


text:
const TEST_NAME = "MT19937"

include("common.jl")

using RNG.MersenneTwisters

r = MT19937(123)

test_all(r, 100_000_000)

ext: jl


text:
#!/usr/bin/env julia

using ArgParse
using Word2Vec

s = ArgParseSettings()
@add_arg_table s begin
    "--model", "-m"
    arg_type = String
    help = "use this model file (write given -i, query given -q)"
    required = true
    default = ""
    "--input", "-i"
    arg_type = String
    help = "input text corpus to model"
    default = ""
    "--kmer-size", "-k"
    arg_type = Int
    help = "kmer size to break input DNA sequences into"
    default = 5
    "--kmer-stride", "-s"
    arg_type = Int
    help = "bases between each successive kmer"
    default = 1
    "--window-size", "-w"
    arg_type = Int
    help = "window size to use during word2vec"
    default = 30
    "--cbow", "-c"
    help = "use continuous bag of words in word2vec"
    action = :store_true
    "--dims", "-d"
    arg_type = Int
    help = "number of dimensions in embedding"
    default = 100
    "--query", "-q"
    arg_type = String
    help = "print the vector for this string"
    default = ""
end

parsed_args = parse_args(ARGS, s)
input = parsed_args["input"]
model_file = parsed_args["model"]
model_size = parsed_args["dims"]
query = parsed_args["query"]
kmer_size = parsed_args["kmer-size"]
kmer_stride = parsed_args["kmer-stride"]
window_size = parsed_args["window-size"]
cbow = parsed_args["cbow"]

function kmers_of(k::Int, j::Int, s::String)
    [s[i:(i+k-1)] for i in range(1,j,Int(floor((length(s)-k)/j))+1)]
end

function kmers_of_file(k::Int, j::Int, f::String)
    kmers = []
    for s in readlines(open(f))
        append!(kmers, kmers_of(k, j, s))
    end
    kmers
end

function write_kmers_of(k::Int, j::Int, input::String, output::String)
    write(open(output, "w"), join(kmers_of_file(k, j, input), " "))
end

if input != ""
    if model_file == ""
        println("an output file is required when building a model")
        exit(1)
    end
    kmers = "$input.kmers"
    write_kmers_of(kmer_size, kmer_stride, input, kmers)
    word2vec(kmers, model_file, verbose = true, size = model_size, window=window_size, cbow=(cbow?1:0))
    rm(kmers)
elseif model_file != ""
    model = wordvectors(model_file)
    if query != ""
        kmers = kmers_of(kmer_size, kmer_stride, query)
        mean_vec = reduce(+, [(haskey(model.vocab_hash, kmer) ? get_vector(model, kmer) : 0) for kmer in kmers])/length(kmers)
        println("$query ", mean_vec)
    end
end

max_stars_repo_name: ruolin/vg

max_issues_repo_name: ruolin/vg

max_forks_repo_name: ruolin/vg


text:
pyeval_filename(src) =
    isfile(string(src.file)) ? "$(src.file):$(src.line)" : "julia:$(src.file):$(src.line)"

pyeval_macro(filename, mode, codearg, args...) = begin
    # unwrap rettype
    if codearg isa Expr && codearg.head == :(::)
        rettypearg = codearg.args[end]
        codearg = codearg.args[1]
    else
        rettypearg = mode == :exec ? Nothing : PyObject
    end
    # extract the code
    code = codearg.args[end]
    # parse the remaining arguments - the first one is optionally the locals, the rest are var=value pairs to include in the locals
    kvs = []
    locals = nothing
    for (i, arg) in enumerate(args)
        if arg isa Expr && arg.head == :(=) && arg.args[1] isa Symbol
            push!(kvs, arg.args[1] => arg.args[2])
        elseif i == 1
            locals = arg
        else
            error()
        end
    end
    # find any interpolations
    i = firstindex(code)
    codechunks = String[]
    interps = []
    chunk = ""
    while true
        j = findnext(==('$'), code, i)
        if j === nothing
            push!(codechunks, chunk * code[i:end])
            break
        elseif checkbounds(Bool, code, j + 1) && code[j+1] == '$'
            chunk *= code[i:j]
            i = j + 2
        else
            push!(codechunks, chunk * code[i:j-1])
            chunk = ""
            ex, i = Meta.parse(code, j + 1, greedy = false)
            push!(interps, ex)
        end
    end
    # identify which are LHS and which are RHS
    # currently only the pattern "^ $(...) =" is recognized as a LHS
    # TODO: multiple assignment
    # TODO: mutating assignment
    islhs = [
        (
            (match(r"[\n;=]\s*$", codechunks[i]) !== nothing) ||
            (i == 1 && match(r"^\s*$", codechunks[i]) !== nothing)
        ) && match(r"^\s*=($|[^=])", codechunks[i+1]) !== nothing for
        (i, ex) in enumerate(interps)
    ]
    # do the interpolation
    intvars = ["_jl_interp_$(i)_" for i = 1:length(interps)]
    for (k, v, lhs) in zip(intvars, interps, islhs)
        lhs || push!(kvs, Symbol(k) => v)
    end
    newcode = join([
        c * (
            checkbounds(Bool, intvars, i) ? islhs[i] ? intvars[i] : "($(intvars[i]))" : ""
        ) for (i, c) in enumerate(codechunks)
    ])
    # for LHS interpolation, extract out type annotations
    if any(islhs)
        mode == :exec || error("interpolation on LHS only allowed in exec mode")
    end
    if mode == :exec
        outkvts = [
            (ex isa Expr && ex.head == :(::)) ?
            (esc(ex.args[1]), v, esc(ex.args[2])) : (esc(ex), v, PyObject) for
            (ex, v, lhs) in zip(interps, intvars, islhs) if lhs
        ]
    end
    # make the code be a function body
    if mode == :execr
        newcode = "def _jl_tmp_ans_func_($(join(intvars, ", "))):\n" * join(map(line->"    "*line, split(newcode, "\n")), "\n") * "\n\nans = _jl_tmp_ans_func_($(join(intvars, ", ")))"
    end
    # make the code object
    co = PyCode(newcode, filename, mode == :eval ? :eval : :exec)
    # go
    freelocals =
        locals === nothing ? :(C.Py_DecRef(lptr)) : :(GC.@preserve locals nothing)
    ret = quote
        let
            # evaluate the inputs (so any errors are thrown before we have object references)
            $([:($(Symbol(:input, i)) = $(esc(v))) for (i, (k, v)) in enumerate(kvs)]...)
            # get the code pointer
            cptr = checknull(pyptr($co))
            # get the globals pointer
            globals = $(esc(:pyglobals))
            gptr = checknull(pyptr(globals))
            # ensure globals includes builtins
            if globals isa PyDict && !globals.hasbuiltins
                if C.PyMapping_HasKeyString(gptr, "__builtins__") == 0
                    err = C.PyMapping_SetItemString(
                        gptr,
                        "__builtins__",
                        C.PyEval_GetBuiltins(),
                    )
                    ism1(err) && pythrow()
                end
                globals.hasbuiltins = true
            end
            # get locals (ALLOCATES lptr if locals===nothing)
            $(
                locals === nothing ? :(lptr = checknull(C.PyDict_New())) :
                :(locals = $(esc(locals)); lptr = checknull(pyptr(locals)))
            )
            # insert extra locals
            $(
                [
                    :(
                        let
                            vo = C.PyObject_From($(Symbol(:input, i)))
                            isnull(vo) && ($freelocals; pythrow())
                            err = C.PyObject_SetItem(
                                lptr,
                                $(PyInternedString(string(k))),
                                vo,
                            )
                            C.Py_DecRef(vo)
                            ism1(err) && ($freelocals; pythrow())
                        end
                    ) for (i, (k, v)) in enumerate(kvs)
                ]...
            )
            # Call eval (ALLOCATES rptr)
            rptr = GC.@preserve globals C.PyEval_EvalCode(cptr, gptr, lptr)
            isnull(rptr) && ($freelocals; pythrow())
            # extract values
            $(
                if mode == :eval
                    quote
                        $freelocals
                        res = C.PyObject_Convert(rptr, $(esc(rettypearg)))
                        C.Py_DecRef(rptr)
                        ism1(res) && pythrow()
                        C.takeresult($(esc(rettypearg)))
                    end
                elseif mode in (:execa, :execr)
                    quote
                        C.Py_DecRef(rptr)
                        xo = C.PyObject_GetItem(lptr, $(PyInternedString("ans")))
                        isnull(xo) && ($freelocals; pythrow())
                        res = C.PyObject_Convert(xo, $(esc(rettypearg)))
                        C.Py_DecRef(xo)
                        $freelocals
                        ism1(res) && pythrow()
                        C.takeresult($(esc(rettypearg)))
                    end
                elseif mode == :exec
                    quote
                        C.Py_DecRef(rptr)
                        $(
                            (
                                quote
                                    $(Symbol(:output, i)) =
                                        let
                                            xo = C.PyObject_GetItem(
                                                lptr,
                                                $(PyInternedString(v)),
                                            )
                                            isnull(xo) && ($freelocals; pythrow())
                                            res = C.PyObject_Convert(xo, $t)
                                            C.Py_DecRef(xo)
                                            ism1(res) && ($freelocals; pythrow())
                                            C.takeresult($t)
                                        end
                                end for (i, (_, v, t)) in enumerate(outkvts)
                            )...
                        )
                        $freelocals
                        ($((Symbol(:output, i) for i = 1:length(outkvts))...),)
                    end
                else
                    error()
                end
            )
        end
    end
    if mode == :exec
        ret = quote
            ($((k for (k, _, _) in outkvts)...),) = $ret
            nothing
        end
    end
    ret
end

"""
    @py `...` [locals] [var=val, ...]

Execute the given Python code.

Julia values can be interpolated using the usual `\$(...)` syntax.

Additionally, assignment to interpolations is supported: e.g. `\$(x::T) = ...` will convert the right hand side to a `T` and assign it to `x`.
- Currently only single assignment is supported. Multiple assignment (`\$x, \$y = ...`) or mutating assignment (`\$x += ...`) will not be recognized.
- What actually happens is that the assignment is to a temporary Python variable, which is then read when execution successfully finishes.
  Hence if an exception occurs, no assignments will happen.

The globals are `pyglobals`.
The locals are `locals`, if given, otherwise a temporary scope is created. Extra values to be interted into the scope can be given with extra `var=val` arguments.
"""
macro py(args...)
    pyeval_macro(pyeval_filename(__source__), :exec, args...)
end
export @py

"""
    @pyg `...` [var=val, ...]

Execute the given Python code in the global scope.

This is simply shorthand for ```@py `...` pyglobals ``` (see [`@py`](@ref)).
"""
macro pyg(code, args...)
    pyeval_macro(pyeval_filename(__source__), :exec, code, :pyglobals, args...)
end
export @pyg

"""
    @pyv `...`[::rettype] [locals] [var=val, ...]

Evaluate the given Python expression and return its value.

Julia values can be interpolated using the usual `\$(...)` syntax.

The globals are `pyglobals`.
The locals are `locals`, if given, otherwise a temporary scope is created. Extra values to be interted into the scope can be given with extra `var=val` arguments.

The result is converted to a `rettype`, which defaults to `PyObject`.
"""
macro pyv(args...)
    pyeval_macro(pyeval_filename(__source__), :eval, args...)
end
export @pyv

"""
    @pya `...`[::rettype] [locals] [var=val, ...]

Execute the given Python code and return `ans`.

This is the same as `@py ...` except that the variable `ans` is extracted from the scope and returned.
"""
macro pya(args...)
    pyeval_macro(pyeval_filename(__source__), :execa, args...)
end
export @pya

"""
    @pyr `...`[::rettype] [locals] [var=val, ...]

Execute the given Python code in a function and return its return value.

Essentially equivalent to ```@pya `def result(): ...; ans = result()` ```.
"""
macro pyr(args...)
    pyeval_macro(pyeval_filename(__source__), :execr, args...)
end
export @pyr

"""
    py`...` :: PyCode

Literal syntax for a compiled [`PyCode`](@ref) object in "exec" mode.

Suitable for passing to Python's `exec` function.
"""
macro py_cmd(code::String)
    PyCode(code, pyeval_filename(__source__), :exec)
end
export @py_cmd

"""
    pyv`...` :: PyCode

Literal syntax for a compiled [`PyCode`](@ref) object in "eval" mode.

Suitable for passing to Python's `eval` function.
"""
macro pyv_cmd(code::String)
    PyCode(code, pyeval_filename(__source__), :eval)
end
export @pyv_cmd


text:
push!(LOAD_PATH, pwd())

using Distributed
import HDF5
import Statistics
#using PyCall
#pygui(:tk) # :tk, :gtk3, :gtk, :qt5, :qt4, :qt, or :wx
import PyPlot
@everywhere import .TOFFunctions
import .InterpolationFunctions

function correctMassScaleAndExtractSumSpec(
  filepath,
  masslistMasses,
  masslistElements,
  masslistElementsMasses,
  masslistCompositions,
  referenceFile,
  calibRegions; # Regions which are pattern matched for mass shifts
  filefilterRegexp=r"\.h5$",
  outputfilename="results/_result.hdf5",
  firstNFiles=0, # only analyze the first N files, set to 0 for all files
  lastNFiles=0, # only analyze the last N files, set to 0 for all files
  debuglevel=3,
  searchWidth = 0.7, # Width of Mass Scale Search Regions in AMU
  dynamicMassScaleCorrection = true,
  recalibInterval = 60,
  peakshapeRegions = 10,
  createTotalAvg = true, # Needed for manualPeakFitter, not needed for timetraces
  onlyUseAverages = false, # Fast mode using only total file averages instead of individual spectra in file
  storeSubSpectra = true,
  plotControlMass = false, # plot the control mass peak for each file after correction
  testRangeStart = 137.0, # the mass shift of this region will be shown if plot control mass is set true. Should not be part of calibRegions
  testRangeEnd = 137.5,
  massBorderCalculation = 1, # How to calculate borders? 0 = Cernter -0.1 to Center + 0.4,  1 = based on resolution and peak distance, 2 = constant bin width
  binWidth = 6,
  resolution = 7500,
  filePrecaching = true, # Precache whole next file. Only use when you have enough RAM and no SWAP has to be used.
  openWholeFile = true
  )

  if ! dynamicMassScaleCorrection #TODO: dummy fix, should be treated in workflow
      recalibInterval = 1e99
  end

  if (isfile(joinpath(filepath, outputfilename)))
    mv(joinpath(filepath, outputfilename), joinpath(filepath, "$outputfilename.bak"), force=true)
    if debuglevel > 0   println("Found existing result file, moving to '$outputfilename.bak'") end
  end

  #files = filter(r"2016-07-04.*\.h5", readdir(filepath)) #one day only
  allFiles = readdir(filepath)
  files = filter(s->occursin(filefilterRegexp, s), allFiles)

  nFiles = size(files,1)

  if firstNFiles != 0
    if length(files) > firstNFiles
      files = files[1:firstNFiles]
    end
  end

  if lastNFiles != 0
    if length(files) > lastNFiles
      files = files[length(files)-lastNFiles:end]
    end
  end

  nFiles = size(files,1)

  validFiles, timeSortIndices = TOFFunctions.validateHDF5Files(filepath, files)
  if !issorted(timeSortIndices)
      println("Reordering Files according to start acquisition time!")
  else
      println("File order seems fine, continuing...")
  end
  files = validFiles[timeSortIndices]
  nFiles = size(files,1)

  if (referenceFile == "")
    referenceFile = joinpath(filepath, files[1])
    if (debuglevel > 1)
      println("Reference File Autodetected (using first file in folder): $referenceFile")
    end
  end

  referenceSpectrum = TOFFunctions.getAvgSpectrumFromFile(referenceFile)
  referenceMassScaleMode, referenceMassScaleParameters = TOFFunctions.getMassCalibParametersFromFile(referenceFile)
  TOFFunctions.setMassScaleReferenceSpectrum(referenceSpectrum, calibRegions, searchWidth, referenceMassScaleMode, referenceMassScaleParameters, plotControlMass=plotControlMass, testRangeStart=testRangeStart, testRangeEnd=testRangeEnd)
  referenceMassAxis = []
  referenceMassAxis = HDF5.h5read(referenceFile, "FullSpectra/MassAxis")

  # Tofwerk does not update the mass axis after manual recalibration in TofDAQ Viewer.
  # ==> recalculate based on calib parameters.
  println("Recalculating mass axis")
  for i=1:length(referenceMassAxis)
      referenceMassAxis[i] = TOFFunctions.timebin2mass(i, referenceMassScaleMode, referenceMassScaleParameters)
  end
  #totalAvgSpectrum = SharedArray{Float64}(undef,length(referenceSpectrum))
  #totalAvgSubSpectrum = SharedArray{Float64}(undef,length(referenceSpectrum))
  #totalMinSpectrum = SharedArray{Float64}(undef,length(referenceSpectrum))
  #totalMaxSpectrum = SharedArray{Float64}(undef,length(referenceSpectrum))
  totalAvgSpectrum = Array{Float64}(undef,length(referenceSpectrum))
  totalAvgSubSpectrum = Array{Float64}(undef,length(referenceSpectrum))
  totalMinSpectrum = Array{Float64}(undef,length(referenceSpectrum))
  totalMaxSpectrum = Array{Float64}(undef,length(referenceSpectrum))

  #@sync @parallel
  for i = 1 : length(referenceSpectrum)
    totalAvgSpectrum[i] = 0
    totalAvgSubSpectrum[i] = 0
    totalMinSpectrum[i] = 1e99
    totalMaxSpectrum[i] = 0
  end

  nMasses=length(masslistMasses)
  println("Calculating Stick CPS for $nMasses masses.")
  timebinshifts = zeros(Float64,(length(calibRegions),nFiles))
  intensities = zeros(Float64,(length(calibRegions),nFiles))
  monitorTimetrace = zeros(nFiles)

  #time = SharedArray{DateTime}(undef,nFiles)
  #stickcps=SharedArray{Float64}(undef,nMasses,nFiles)
  time = Array{DateTime}(undef,nFiles)
  stickcps=Array{Float64}(undef,nMasses,nFiles)

  # Calculate Integration Borders
  if (massBorderCalculation == 0)
  mlow = masslistMasses - 0.1
  mhigh = masslistMasses + 0.4
  elseif  (massBorderCalculation == 1)
  #mlow = masslistMasses .* resolution./(resolution+0.5)
  #mhigh = masslistMasses .* resolution./(resolution-0.5)
  massborders = MasslistFunctions.IntegrationBorders(masslistMasses, resolution)
  end

  if (plotControlMass == true)
    PyPlot.figure()
  end

  if (createTotalAvg == true)
      #interpolatedSpectrum = SharedArray{Float64}(undef,length(referenceSpectrum))
      interpolatedSpectrum = Array{Float64}(undef,length(referenceSpectrum))
  end

  ############## Check output path and remove existing files #####################
  if (!isdir("$filepath/results"))
    mkdir("$filepath/results")
  end
  outfilepath = joinpath(filepath, outputfilename)
  if isfile(outfilepath)
    rm(outfilepath)
  end

  fid = HDF5.h5open(outfilepath, "w")

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

  ############## open outputfile and create extendable SumSpecs Array ############
  if (createTotalAvg == true)

    dsAvgSumWidth = length(referenceSpectrum)
    dspaceAvgSumSpecs = HDF5.dataspace((dsAvgSumWidth,1)::Dims, max_dims=(dsAvgSumWidth,typemax(Int64)))
    dtypeAvgSumSpecs = HDF5.datatype(Float32)
    dsetAvgSumSpecs = HDF5.d_create(fid, "SumSpecs", dtypeAvgSumSpecs, dspaceAvgSumSpecs, "chunk", (dsAvgSumWidth,1), "compress", 3)
  end

  if (!onlyUseAverages)
    specCount = 0
    dsStickCpsWidth = nMasses
    dspaceStickCps = HDF5.dataspace((1,dsStickCpsWidth)::Dims, max_dims=(typemax(Int64),dsStickCpsWidth))
    dtypeStickCps = HDF5.datatype(Float32)
    dsetStickCps = HDF5.d_create(fid, "StickCps", dtypeStickCps, dspaceStickCps, "chunk", (1,dsStickCpsWidth), "compress", 3)
    dsetStickCpsErr = HDF5.d_create(fid, "StickCpsErr", dtypeStickCps, dspaceStickCps, "chunk", (1,dsStickCpsWidth), "compress", 3)

    dspaceTimes = HDF5.dataspace((1,)::Dims, max_dims=(typemax(Int64),))
    dtypeTimes = HDF5.datatype(Float64)
    dsetTimes = HDF5.d_create(fid, "Times", dtypeTimes, dspaceTimes, "chunk", (1,))

  end
  ################################################################################
  badFiles = Array{String,1}()

  for j=1:nFiles
    totalPath = joinpath(filepath, files[j])
    if j < nFiles && filePrecaching
        totalPrecachePath = joinpath(filepath, files[j+1])
    else
        totalPrecachePath = ""
    end
    if debuglevel > 0   println("Processing File $j/$nFiles :  $totalPath") end

    fileIsBad = false;
    massAxis = []
    massAxis = HDF5.h5read(totalPath, "FullSpectra/MassAxis")

    time[j] = TOFFunctions.getTimeFromFile(totalPath)

    avgSpectrum = TOFFunctions.getAvgSpectrumFromFile(totalPath)
    newParams, success, tbs, ins = TOFFunctions.recalibrateMassScale(avgSpectrum, referenceSpectrum, calibRegions, searchWidth, referenceMassScaleMode, referenceMassScaleParameters)

    #if debuglevel > 0   println("Processing File $totalPath") end
    if (createTotalAvg == true && !fileIsBad)
      print("Spectrum interpolation for total average... ")
      indexesExact = TOFFunctions.mass2timebin(referenceMassAxis, referenceMassScaleMode, newParams)
      interpolatedSpectrum = InterpolationFunctions.interpolate(indexesExact,avgSpectrum)
      totalAvgSpectrum += interpolatedSpectrum
      # calculate min and max spectra
      #@sync @parallel
      for bin=1:length(interpolatedSpectrum)
        if interpolatedSpectrum[bin] > totalMaxSpectrum[bin]
          totalMaxSpectrum[bin] = interpolatedSpectrum[bin]
        end
        if interpolatedSpectrum[bin] < totalMinSpectrum[bin]
          totalMinSpectrum[bin] = interpolatedSpectrum[bin]
        end
      end
      # Write to hdf5, line by line, so there is no limit to number of files that can fit in RAM
      if (storeSubSpectra == true)
        currDims = size(dsetAvgSumSpecs)[2]
        dsetAvgSumSpecs[:,currDims] = interpolatedSpectrum
        if (j != nFiles)
          HDF5.set_dims!(dsetAvgSumSpecs, (dsAvgSumWidth,currDims+1)::Dims)
        end
      end
      println("DONE")
    end
    if debuglevel > 1   println() end
    print("Average Stick integration... ")
    for i=(1:nMasses)
      #if debuglevel > 0   println("Processing region mass($(mlow[i]):$(mhigh[i])) --> timebin($(round(TOFFunctions.mass2timebin(mlow[i],massCalibMode,newParams))):$(round(TOFFunctions.mass2timebin(mhigh[i],massCalibMode,newParams))))") end
      if (massBorderCalculation == 2)
        centerIndex = Int64(floor(TOFFunctions.mass2timebin(masslistMasses[i], referenceMassScaleMode, newParams)))
        if ((centerIndex-binWidth) > 0) && ((centerIndex+binWidth)<length(avgSpectrum))
          raw = sum(view(avgSpectrum, (centerIndex-binWidth : centerIndex + binWidth)))
        else
          raw = 0
        end
      else
        subIdxStartExact=TOFFunctions.mass2timebin(massborders.lowMass[i],referenceMassScaleMode,newParams)
        subIdxEndExact = TOFFunctions.mass2timebin(massborders.highMass[i],referenceMassScaleMode,newParams)
        raw = InterpolationFunctions.interpolatedSum(subIdxStartExact,subIdxEndExact,avgSpectrum)
      end
      stickcps[i,j] = raw
    end
    println("DONE")
    if (!onlyUseAverages)
      fileInternalLocalAvg = 0 #clear big data and references for garbage collection
      GC.gc()
      spectrumMultFactor = TOFFunctions.getSpecMultiplicator(totalPath)
      println("   spectrumMultFactor = $spectrumMultFactor")
      #subSpecStickCps = SharedArray{Float32}(undef,nMasses)
      subSpecStickCps = Array{Float32}(undef,nMasses)
      totalSubSpectra = TOFFunctions.getSubSpectraCount(totalPath)

      ################## prepare arrays for averaging #########################
      fileInternalLocalAvg = TOFFunctions.getSubSpectrumFromFile(totalPath,1, preloadFile=totalPrecachePath, openWholeFile = openWholeFile)
      if fileInternalLocalAvg == 0 # Skip file if the spectra could not be loaded
        println("Error: Could not load Subspectra from file $totalPath, skipping whole file!")
        continue
      end
      #println("referenceSpectrum is a $(summary(referenceSpectrum))")
      #println("SubSpec is a $(summary(fileInternalLocalAvg))")
      #println("spectrumMultFactor is a $(summary(spectrumMultFactor))")
      print("Sub spectrum Stick integration... ")

      fill!(fileInternalLocalAvg,0)
      fileInternalLocalAvgCount = 0
      print("Precalculating average for mass scale correction... ")
      ################## Get First Set of Spectra for first mass scale calib
      if totalSubSpectra >= recalibInterval
        # Prepare first calib beforehead
        fileInternalLocalAvg = TOFFunctions.getSubSpectrumFromFile(totalPath,1, preloadFile=totalPrecachePath, openWholeFile = openWholeFile)
        for avgIdx=2:minimum([recalibInterval totalSubSpectra])
          fileInternalLocalAvg += TOFFunctions.getSubSpectrumFromFile(totalPath,avgIdx, preloadFile=totalPrecachePath, openWholeFile = openWholeFile)
        end
        newParams, success, tbs, ins = TOFFunctions.recalibrateMassScale(fileInternalLocalAvg, referenceSpectrum, calibRegions, searchWidth, referenceMassScaleMode, referenceMassScaleParameters)
        ######
        fill!(fileInternalLocalAvg,0)
        fileInternalLocalAvgCount = 0
      end
      #######################################################################

      ################### Loop over Subspectra Start ########################
      ################### Averaging and Mass Scale Calib ####################
      print("Recalibrating at: ")
      specCount += totalSubSpectra
      for subSpecIdx=1:totalSubSpectra
        subSpectrum = TOFFunctions.getSubSpectrumFromFile(totalPath,subSpecIdx, preloadFile=totalPrecachePath, openWholeFile = openWholeFile)
        fileInternalLocalAvg += subSpectrum
        fileInternalLocalAvgCount += 1

        if (fileInternalLocalAvgCount > recalibInterval) || fileInternalLocalAvgCount == totalSubSpectra
          print("$subSpecIdx...")
          newParams, success, tbs, ins = TOFFunctions.recalibrateMassScale(fileInternalLocalAvg, referenceSpectrum, calibRegions, searchWidth, referenceMassScaleMode, referenceMassScaleParameters)

          indexesExact = TOFFunctions.mass2timebin(referenceMassAxis, referenceMassScaleMode, newParams)
          interpolatedSpectrum = InterpolationFunctions.interpolate(indexesExact,fileInternalLocalAvg)

          totalAvgSubSpectrum += interpolatedSpectrum*spectrumMultFactor

          fill!(fileInternalLocalAvg,0)
          fileInternalLocalAvgCount = 0
        end
        ################## Peak Integration ################################
        #@sync @parallel
        mbIndicesLowMass = TOFFunctions.mass2timebin(massborders.lowMass,referenceMassScaleMode,newParams)
        mbIndicesHighMass = TOFFunctions.mass2timebin(massborders.highMass,referenceMassScaleMode,newParams)
        Threads.@threads for i=(1:nMasses)
          if (mod(i,100) == 0)
            #println("Processing mass $(masslistMasses[i])")
          end
          if (massBorderCalculation == 2)
            centerIndex = Int64(floor(TOFFunctions.mass2timebin(masslistMasses[i], referenceMassScaleMode, newParams)))
            if ((centerIndex-binWidth) > 0) && ((centerIndex+binWidth)<length(avgSpectrum))
              raw = sum(view(subSpectrum, (centerIndex-binWidth : centerIndex + binWidth)))
            else
              raw = 0
            end
            subSpecStickCps[i]=raw
          else
            #if debuglevel > 0   println("Processing region mass($(mlow[i]):$(mhigh[i])) --> timebin($(round(TOFFunctions.mass2timebin(mlow[i],massCalibMode,newParams))):$(round(TOFFunctions.mass2timebin(mhigh[i],massCalibMode,newParams))))") end
            subSpecStickCps[i]= InterpolationFunctions.interpolatedSum(mbIndicesLowMass[i],mbIndicesHighMass[i],subSpectrum)
          end
        end
        rawTime = TOFFunctions.getSubSpectrumTimeFromFile(totalPath,subSpecIdx)
        currDimsTime = size(dsetTimes)[1]
        absTimeOfCurrSample = Dates.datetime2unix(time[j]) + rawTime
        dsetTimes[currDimsTime] = absTimeOfCurrSample
        dsetStickCps[currDimsTime,:] = subSpecStickCps*spectrumMultFactor
        if (currDimsTime == 1)
          deltaT = 5 #rawTime # Not correct???
        else
          deltaT = absTimeOfCurrSample - dsetTimes[currDimsTime-1][1]
        end
        #println("DeltaT: $deltaT")
        dsetStickCpsErr[currDimsTime,:] = sqrt.(abs.(subSpecStickCps./deltaT))
        if !((j==nFiles) && (subSpecIdx== TOFFunctions.getSubSpectraCount(totalPath)))
          HDF5.set_dims!(dsetTimes, (currDimsTime+1,)::Dims)
          HDF5.set_dims!(dsetStickCps, (currDimsTime+1,dsStickCpsWidth)::Dims)
          HDF5.set_dims!(dsetStickCpsErr, (currDimsTime+1,dsStickCpsWidth)::Dims)
        end
      end
      println("DONE")
    end
    println("###################################################################\n\n")
  end
  close(fid)

  stickcpsFlat = zeros(nFiles,nMasses)
  stickcpsErrFlat = zeros(nFiles,nMasses)

  timeFlat = zeros(nFiles)
  timeFlat = Dates.datetime2unix.(time)
  # We need more than one timestamp to calculate a difference
  if length(timeFlat) > 1
      deltaT = Statistics.median(diff(timeFlat))
  else
      deltaT = 1
  end

  totalAvgSpectrumFlat = zeros(length(totalAvgSpectrum))
  totalAvgSubSpectrumFlat = zeros(length(totalAvgSubSpectrum))

  for i=1:nMasses, j=1:nFiles
    stickcpsFlat[j,i] = stickcps[i,j]
    stickcpsErrFlat[j,i] = sqrt(abs(stickcps[i,j])/deltaT)
  end

  for i=1:length(totalAvgSpectrum)
    totalAvgSpectrumFlat[i] = totalAvgSpectrum[i]/nFiles
  end
  if !onlyUseAverages
    for i=1:length(totalAvgSpectrum)
      totalAvgSubSpectrumFlat[i] = totalAvgSubSpectrum[i]/specCount
    end
  end
  HDF5.h5write(outfilepath, "AvgStickCps", stickcpsFlat)
  HDF5.h5write(outfilepath, "AvgStickCpsErr", stickcpsErrFlat)
  HDF5.h5write(outfilepath, "MassList", masslistMasses)
  println("Wrote $(length(masslistMasses)) masses to file.")
  HDF5.h5write(outfilepath, "MassListIntegrationBordersLow", massborders.lowMass)
  HDF5.h5write(outfilepath, "MassListIntegrationBordersHigh", massborders.highMass)

  HDF5.h5write(outfilepath, "MassListIdx", TOFFunctions.mass2timebin(masslistMasses, referenceMassScaleMode, referenceMassScaleParameters))
  HDF5.h5write(outfilepath, "MassListIntegrationBordersIdxLow", TOFFunctions.mass2timebin(massborders.lowMass, referenceMassScaleMode, referenceMassScaleParameters))
  HDF5.h5write(outfilepath, "MassListIntegrationBordersIdxHigh", TOFFunctions.mass2timebin(massborders.highMass, referenceMassScaleMode, referenceMassScaleParameters))
  HDF5.h5write(outfilepath, "MassCalibMode", referenceMassScaleMode)
  HDF5.h5write(outfilepath, "MassCalibParameters", referenceMassScaleParameters)

  HDF5.h5write(outfilepath, "AvgStickCpsTimes", timeFlat)
  HDF5.h5write(outfilepath, "MassAxis", referenceMassAxis)
  if onlyUseAverages
    HDF5.h5write(outfilepath, "AvgSpectrum", totalAvgSpectrumFlat)
  else
    HDF5.h5write(outfilepath, "AvgSpectrum", totalAvgSubSpectrumFlat)
  end

  if (createTotalAvg == true)
    #HDF5.h5write(outfilepath, "SumSpecs", sdata(interpolatedSpectra))
    HDF5.h5write(outfilepath, "SumSpecMax", totalMaxSpectrum)
    HDF5.h5write(outfilepath, "SumSpecMin", totalMinSpectrum)
  end
  HDF5.h5write(outfilepath, "ElementNames", masslistElements)
  HDF5.h5write(outfilepath, "ElementMasses", masslistElementsMasses)
  compFlat = Array{Int64}(undef,length(masslistCompositions[1]),length(masslistCompositions))
  [compFlat[:,i]=masslistCompositions[i] for i=1:length(masslistCompositions)]
  HDF5.h5write(outfilepath, "ElementalCompositions", compFlat)

end


text:
include("../utils.jl")

@testset "Probability simplex" begin
    M = ProbabilitySimplex(2)
    @test M^4 == MultinomialMatrices(3, 4)
    p = [0.1, 0.7, 0.2]
    q = [0.3, 0.6, 0.1]
    X = zeros(3)
    Y = [-0.1, 0.05, 0.05]
    @test is_point(M, p)
    @test_throws DomainError is_point(M, p .+ 1, true)
    @test_throws DomainError is_point(M, [0], true)
    @test_throws DomainError is_point(M, -ones(3), true)
    @test manifold_dimension(M) == 2
    @test is_vector(M, p, X)
    @test is_vector(M, p, Y)
    @test_throws DomainError is_vector(M, p .+ 1, X, true)
    @test_throws DomainError is_vector(M, p, zeros(4), true)
    @test_throws DomainError is_vector(M, p, Y .+ 1, true)

    @test Manifolds.default_metric_dispatch(M, Manifolds.FisherRaoMetric()) === Val{true}()

    @test injectivity_radius(M, p) == injectivity_radius(M, p, ExponentialRetraction())
    @test injectivity_radius(M, p, SoftmaxRetraction()) == injectivity_radius(M, p)
    @test injectivity_radius(M, ExponentialRetraction()) == 0
    @test injectivity_radius(M) == 0
    @test injectivity_radius(M, SoftmaxRetraction()) == 0

    pE = similar(p)
    embed!(M, pE, p)
    @test pE == p
    XE = similar(X)
    embed!(M, XE, p, X)
    @test XE == X

    @test mean(M, [p, q]) == shortest_geodesic(M, p, q)(0.5)

    types = [Vector{Float64}]
    TEST_FLOAT32 && push!(types, Vector{Float32})
    TEST_STATIC_SIZED && push!(types, MVector{3,Float64})

    basis_types = (DefaultOrthonormalBasis(), ProjectedOrthonormalBasis(:svd))
    for T in types
        @testset "Type $T" begin
            pts = [
                convert(T, [0.5, 0.3, 0.2]),
                convert(T, [0.4, 0.4, 0.2]),
                convert(T, [0.3, 0.5, 0.2]),
            ]
            test_manifold(
                M,
                pts,
                test_injectivity_radius=false,
                test_project_tangent=true,
                test_musical_isomorphisms=true,
                test_vee_hat=false,
                test_forward_diff=false,
                test_reverse_diff=false,
                is_tangent_atol_multiplier=5.0,
                inverse_retraction_methods=[SoftmaxInverseRetraction()],
                retraction_methods=[SoftmaxRetraction()],
                test_inplace=true,
            )
        end
    end

    @testset "Projection testing" begin
        p = [1 / 3, 1 / 3, 1 / 3]
        q = [0.2, 0.3, 0.5]
        X = log(M, p, q)
        X2 = X .+ 10
        Y = project(M, p, X2)
        @test isapprox(M, p, X, Y)

        # Check adaption of metric and representer
        Y1 = change_metric(M, EuclideanMetric(), p, X)
        @test Y1 == X .* sqrt.(p)
        Y2 = change_representer(M, EuclideanMetric(), p, X)
        @test Y2 == X .* p

        X = log(M, q, p)
        X2 = X + [1, 2, 3]
        Y = project(M, q, X2)
        @test is_vector(M, q, Y; atol=1e-15)

        @test_throws DomainError project(M, [1, -1, 2])
        @test isapprox(M, [0.6, 0.2, 0.2], project(M, [0.3, 0.1, 0.1]))
    end
end


text:
using PortAudio
using Documenter: deploydocs, makedocs

makedocs(
    sitename = "PortAudio.jl",
    modules = [PortAudio],
    pages = [
        "Public interface" => "index.md",
        "Internals" => "internals.md"
    ]
)
deploydocs(repo = "github.com/JuliaAudio/PortAudio.jl.git")


text:
macro jl15_str(code::AbstractString)
    if VERSION >= v"1.5-rc0"
        @debug "Parsing code for Julia ≥ 1.5" Text(code)
        expr = Meta.parse(string("begin\n", code, "\nend"))
        @assert expr.head === :block
        if expr.args[1] isa LineNumberNode
            expr.args[1] = __source__
        end
        return esc(expr)
    else
        @debug "Skipping code in Julia < 1.5" Text(code)
        return nothing
    end
end

macro callmacro(ex)
    @assert Meta.isexpr(ex, :macrocall)
    return Expr(
        :call,
        ex.args[1],
        QuoteNode(ex.args[2]),
        __module__,
        map(QuoteNode, ex.args[3:end])...,
    )
end


text:
"""
    PlaybackPolicy{A<:AbstractArray, P<:Policy, V<:AbstractArray{<:Real}}
a policy that applies a fixed sequence of actions until they are all used and then falls back onto a backup policy until the end of the episode.

Constructor:

    `PlaybackPolicy(actions::AbstractArray, backup_policy::Policy; logpdfs::AbstractArray{Float64, 1} = Float64[])`

# Fields
- `actions::Vector{A}` a vector of actions to play back
- `backup_policy::Policy` the policy to use when all prescribed actions have been taken but the episode continues
- `logpdfs::Vector{Float64}` the log probability (density) of actions
- `i::Int64` the current action index
"""
mutable struct PlaybackPolicy{A<:AbstractArray, P<:Policy, V<:AbstractArray{<:Real}} <: Policy
    actions::A
    backup_policy::P
    logpdfs::V
    i::Int64
end

# Constructor for the PlaybackPolicy
PlaybackPolicy(actions::AbstractArray, backup_policy::Policy; logpdfs::AbstractArray{<:Real} = Float64[]) = PlaybackPolicy(actions, backup_policy, logpdfs, 1)

# Action selection for the PlaybackPolicy
function POMDPs.action(p::PlaybackPolicy, s)
    a = p.i <= length(p.actions) ? p.actions[p.i] : action(p.backup_policy, s)
    p.i += 1
    a
end

# Get the logpdf of the history from the playback policy and the backup policy
function Distributions.logpdf(p::PlaybackPolicy, h)
    N = min(length(p.actions), length(h))
    # @assert all(collect(action_hist(h))[1:N] .== p.actions[1:N])
    @assert length(p.actions) == length(p.logpdfs)
    if length(h) > N
        return sum(p.logpdfs) + sum(logpdf(p.backup_policy, view(h, N+1:length(h))))
    else
        return sum(p.logpdfs[1:N])
    end
end

max_stars_repo_name: JuliaPOMDP/POMDPPolicies

max_issues_repo_name: JuliaPOMDP/POMDPPolicies

max_forks_repo_name: JuliaPOMDP/POMDPPolicies


text:
struct Evaluation{T}
    value :: T
end

(ℰ::Evaluation)(a::Sequence) = *(ℰ, a)
Base.:*(ℰ::Evaluation, a::Sequence) = evaluate(a, ℰ.value)

(a::Sequence)(x) = evaluate(a, x)
(a::Sequence)(x...) = evaluate(a, x)

function evaluate(a::Sequence, x)
    ℰ = Evaluation(x)
    space_a = space(a)
    new_space = image(ℰ, space_a)
    CoefType = _coeftype(ℰ, space_a, eltype(a))
    c = Sequence(new_space, Vector{CoefType}(undef, dimension(new_space)))
    _apply!(c, ℰ, a)
    return c
end

function evaluate!(c::Sequence, a::Sequence, x)
    ℰ = Evaluation(x)
    space(c) == image(ℰ, space(a)) || return throw(DomainError)
    _apply!(c, ℰ, a)
    return c
end

function project(ℰ::Evaluation, domain::VectorSpace, codomain::VectorSpace, ::Type{T}) where {T}
    _iscompatible(domain, codomain) || return throw(DimensionMismatch)
    C = LinearOperator(domain, codomain, zeros(T, dimension(codomain), dimension(domain)))
    _project!(C, ℰ)
    return C
end

function project!(C::LinearOperator, ℰ::Evaluation)
    _iscompatible(domain(C), codomain(C)) || return throw(DimensionMismatch)
    coefficients(C) .= zero(eltype(C))
    _project!(C, ℰ)
    return C
end

# Sequence spaces

image(ℰ::Evaluation{<:NTuple{N,Any}}, s::TensorSpace{<:NTuple{N,BaseSpace}}) where {N} =
    TensorSpace(map((xᵢ, sᵢ) -> image(Evaluation(xᵢ), sᵢ), ℰ.value, spaces(s)))

_coeftype(ℰ::Evaluation, s::TensorSpace, ::Type{T}) where {T} =
    @inbounds promote_type(_coeftype(Evaluation(ℰ.value[1]), s[1], T), _coeftype(Evaluation(Base.tail(ℰ.value)), Base.tail(s), T))
_coeftype(ℰ::Evaluation, s::TensorSpace{<:Tuple{BaseSpace}}, ::Type{T}) where {T} =
    @inbounds _coeftype(Evaluation(ℰ.value[1]), s[1], T)

function _apply!(c::Sequence{<:TensorSpace}, ℰ::Evaluation, a)
    space_a = space(a)
    A = _no_alloc_reshape(coefficients(a), dimensions(space_a))
    C = _no_alloc_reshape(coefficients(c), _effective_dimensions(ℰ, space(c)))
    _apply!(C, ℰ, space_a, A)
    return c
end

_effective_dimensions(ℰ::Evaluation{<:Tuple{Nothing,Vararg}}, space::TensorSpace) =
    (dimension(space[1]), _effective_dimensions(Evaluation(Base.tail(ℰ.value)), Base.tail(space))...)
_effective_dimensions(ℰ::Evaluation{<:Tuple{Nothing}}, space::TensorSpace) =
    (dimension(space[1]),)
_effective_dimensions(ℰ::Evaluation, space::TensorSpace) =
    _effective_dimensions(Evaluation(Base.tail(ℰ.value)), Base.tail(space))
_effective_dimensions(ℰ::Evaluation{<:Tuple{Any}}, space::TensorSpace) = ()

_apply!(C, ℰ::Evaluation, space::TensorSpace, A) =
    @inbounds _apply!(C, Evaluation(ℰ.value[1]), space[1], _apply(Evaluation(Base.tail(ℰ.value)), Base.tail(space), A))

_apply!(C, ℰ::Evaluation, space::TensorSpace{<:Tuple{BaseSpace}}, A) =
    @inbounds _apply!(C, Evaluation(ℰ.value[1]), space[1], A)

_apply(ℰ::Evaluation, space::TensorSpace{<:NTuple{N₁,BaseSpace}}, A::AbstractArray{T,N₂}) where {N₁,T,N₂} =
    @inbounds _apply(Evaluation(ℰ.value[1]), space[1], Val(N₂-N₁+1), _apply(Evaluation(Base.tail(ℰ.value)), Base.tail(space), A))

_apply(ℰ::Evaluation, space::TensorSpace{<:Tuple{BaseSpace}}, A::AbstractArray{T,N}) where {T,N} =
    @inbounds _apply(Evaluation(ℰ.value[1]), space[1], Val(N), A)

function _project!(C::LinearOperator{<:SequenceSpace,<:SequenceSpace}, ℰ::Evaluation)
    domain_C = domain(C)
    codomain_C = codomain(C)
    CoefType = eltype(C)
    @inbounds for β ∈ indices(domain_C), α ∈ indices(codomain_C)
        C[α,β] = _getindex(ℰ, domain_C, codomain_C, CoefType, α, β)
    end
    return C
end

@generated function _getindex(ℰ::Evaluation{<:NTuple{N,Any}}, domain::TensorSpace{<:NTuple{N,BaseSpace}}, codomain::TensorSpace{<:NTuple{N,BaseSpace}}, ::Type{T}, α, β) where {N,T}
    p = :(_getindex(Evaluation(ℰ.value[1]), domain[1], codomain[1], T, α[1], β[1]))
    for i ∈ 2:N
        p = :(_getindex(Evaluation(ℰ.value[$i]), domain[$i], codomain[$i], T, α[$i], β[$i]) * $p)
    end
    return p
end

# Taylor

image(::Evaluation{Nothing}, s::Taylor) = s
image(::Evaluation, ::Taylor) = Taylor(0)

_coeftype(::Evaluation{Nothing}, ::Taylor, ::Type{T}) where {T} = T
_coeftype(::Evaluation{T}, ::Taylor, ::Type{S}) where {T,S} = promote_type(T, S)

function _apply!(c::Sequence{Taylor}, ::Evaluation{Nothing}, a)
    coefficients(c) .= coefficients(a)
    return c
end
function _apply!(c::Sequence{Taylor}, ℰ::Evaluation, a)
    x = ℰ.value
    if iszero(x)
        @inbounds c[0] = a[0]
    else
        ord = order(a)
        @inbounds c[0] = a[ord]
        @inbounds for i ∈ ord-1:-1:0
            c[0] = c[0] * x + a[i]
        end
    end
    return c
end

function _apply!(C::AbstractArray, ::Evaluation{Nothing}, ::Taylor, A)
    C .= A
    return C
end
function _apply!(C::AbstractArray{T}, ℰ::Evaluation, space::Taylor, A) where {T}
    x = ℰ.value
    if iszero(x)
        @inbounds C .= selectdim(A, 1, 1)
    else
        ord = order(space)
        @inbounds C .= selectdim(A, 1, ord+1)
        @inbounds for i ∈ ord-1:-1:0
            C .= C .* x .+ selectdim(A, 1, i+1)
        end
    end
    return C
end

_apply(::Evaluation{Nothing}, ::Taylor, ::Val, A::AbstractArray) = A
function _apply(ℰ::Evaluation, space::Taylor, ::Val{D}, A::AbstractArray{T,N}) where {D,T,N}
    x = ℰ.value
    CoefType = _coeftype(ℰ, space, T)
    if iszero(x)
        @inbounds C = convert(Array{CoefType,N-1}, selectdim(A, D, 1))
    else
        ord = order(space)
        @inbounds C = convert(Array{CoefType,N-1}, selectdim(A, D, ord+1))
        @inbounds for i ∈ ord-1:-1:0
            C .= C .* x .+ selectdim(A, D, i+1)
        end
    end
    return C
end

function _getindex(::Evaluation{Nothing}, ::Taylor, ::Taylor, ::Type{T}, i, j) where {T}
    if i == j
        return one(T)
    else
        return zero(T)
    end
end
function _getindex(ℰ::Evaluation, ::Taylor, ::Taylor, ::Type{T}, i, j) where {T}
    if i == 0
        x = ℰ.value
        if j == 0
            return one(T)
        else
            if iszero(x)
                return zero(T)
            else
                return convert(T, convert(T, x)^j)
            end
        end
    else
        return zero(T)
    end
end

# Fourier

image(::Evaluation{Nothing}, s::Fourier) = s
image(::Evaluation, s::Fourier) = Fourier(0, frequency(s))

_coeftype(::Evaluation{Nothing}, ::Fourier, ::Type{T}) where {T} = T
_coeftype(::Evaluation{T}, s::Fourier, ::Type{S}) where {T,S} =
    promote_type(typeof(cis(frequency(s)*zero(T))), S)

function _apply!(c::Sequence{<:Fourier}, ::Evaluation{Nothing}, a)
    coefficients(c) .= coefficients(a)
    return c
end
function _apply!(c::Sequence{<:Fourier}, ℰ::Evaluation, a)
    x = ℰ.value
    ord = order(a)
    if iszero(x)
        c[0] = a[0]
        @inbounds for j ∈ 1:ord
            c[0] += a[j] + a[-j]
        end
    else
        @inbounds c[0] = a[ord]
        if ord > 0
            @inbounds s₋ = convert(eltype(c), a[-ord])
            eiωx = cis(frequency(a)*x)
            eiωx_conj = conj(eiωx)
            @inbounds for j ∈ ord-1:-1:1
                c[0] = c[0] * eiωx + a[j]
                s₋ = s₋ * eiωx_conj + a[-j]
            end
            @inbounds c[0] = s₋ * eiωx_conj + a[0] + c[0] * eiωx
        end
    end
    return c
end

function _apply!(C::AbstractArray, ::Evaluation{Nothing}, ::Fourier, A)
    C .= A
    return C
end
function _apply!(C::AbstractArray{T}, ℰ::Evaluation, space::Fourier, A) where {T}
    x = ℰ.value
    ord = order(space)
    @inbounds C .= selectdim(A, 1, ord+1)
    if iszero(x)
        @inbounds for j ∈ 1:ord
            C .+= selectdim(A, 1, ord+1+j) .+ selectdim(A, 1, ord+1-j)
        end
    else
        if ord > 0
            eiωx = cis(frequency(space)*x)
            eiωxj = one(eiωx)
            @inbounds for j ∈ 1:ord
                eiωxj *= eiωx
                C .+= selectdim(A, 1, ord+1+j) .* eiωxj .+ selectdim(A, 1, ord+1-j) .* conj(eiωxj)
            end
        end
    end
    return C
end

_apply(::Evaluation{Nothing}, ::Fourier, ::Val, A::AbstractArray) = A
function _apply(ℰ::Evaluation, space::Fourier, ::Val{D}, A::AbstractArray{T,N}) where {D,T,N}
    x = ℰ.value
    CoefType = _coeftype(ℰ, space, T)
    ord = order(space)
    @inbounds C = convert(Array{CoefType,N-1}, selectdim(A, D, ord+1))
    if iszero(x)
        @inbounds for j ∈ 1:ord
            C .+= selectdim(A, D, ord+1+j) .+ selectdim(A, D, ord+1-j)
        end
    else
        if ord > 0
            eiωx = cis(frequency(space)*x)
            eiωxj = one(eiωx)
            @inbounds for j ∈ 1:ord
                eiωxj *= eiωx
                C .+= selectdim(A, D, ord+1+j) .* eiωxj .+ selectdim(A, D, ord+1-j) .* conj(eiωxj)
            end
        end
    end
    return C
end

function _getindex(::Evaluation{Nothing}, domain::Fourier, codomain::Fourier, ::Type{T}, i, j) where {T}
    frequency(domain) == frequency(codomain) || return throw(DomainError)
    if i == j
        return one(T)
    else
        return zero(T)
    end
end
function _getindex(ℰ::Evaluation, domain::Fourier, codomain::Fourier, ::Type{T}, i, j) where {T}
    if i == 0
        x = ℰ.value
        if j == 0 || iszero(x)
            return one(T)
        else
            return convert(T, cis(frequency(domain)*j*x))
        end
    else
        return zero(T)
    end
end

# Chebyshev

image(::Evaluation{Nothing}, s::Chebyshev) = s
image(::Evaluation, ::Chebyshev) = Chebyshev(0)

_coeftype(::Evaluation{Nothing}, ::Chebyshev, ::Type{T}) where {T} = T
_coeftype(::Evaluation{T}, ::Chebyshev, ::Type{S}) where {T,S} = promote_type(T, S)

function _apply!(c::Sequence{Chebyshev}, ::Evaluation{Nothing}, a)
    coefficients(c) .= coefficients(a)
    return c
end
function _apply!(c::Sequence{Chebyshev}, ℰ::Evaluation, a)
    x = ℰ.value
    ord = order(a)
    if iszero(x)
        c[0] = a[0]
        @inbounds for i ∈ 2:2:ord
            c[0] += ifelse(i%4 == 0, 2, -2) * a[i]
        end
    elseif isone(-x)
        c[0] = a[0]
        @inbounds for i ∈ 1:ord
            c[0] += ifelse(isodd(i), -2, 2) * a[i]
        end
    elseif isone(x)
        c[0] = a[0]
        @inbounds for i ∈ 1:ord
            c[0] += 2a[i]
        end
    else
        if ord == 0
            @inbounds c[0] = a[0]
        elseif ord == 1
            @inbounds c[0] = a[0] + 2x * a[1]
        else
            CoefType = eltype(c)
            x2 = 2x
            s = zero(CoefType)
            @inbounds t = convert(CoefType, 2a[ord])
            @inbounds c[0] = zero(CoefType)
            @inbounds for i ∈ ord-1:-1:1
                c[0] = t
                t = x2 * t - s + 2a[i]
                s = c[0]
            end
            @inbounds c[0] = x * t - s + a[0]
        end
    end
    return c
end

function _apply!(C::AbstractArray, ::Evaluation{Nothing}, ::Chebyshev, A)
    C .= A
    return C
end
function _apply!(C::AbstractArray{T}, ℰ::Evaluation, space::Chebyshev, A) where {T}
    x = ℰ.value
    ord = order(space)
    if iszero(x)
        @inbounds C .= selectdim(A, 1, 1)
        @inbounds for i ∈ 2:2:ord
            C .+= ifelse(i%4 == 0, 2, -2) .* selectdim(A, 1, i+1)
        end
    elseif isone(-x)
        @inbounds C .= selectdim(A, 1, 1)
        @inbounds for i ∈ 1:ord
            C .+= ifelse(isodd(i), -2, 2) .* selectdim(A, 1, i+1)
        end
    elseif isone(x)
        @inbounds C .= selectdim(A, 1, 1)
        @inbounds for i ∈ 1:ord
            C .+= 2 .* selectdim(A, 1, i+1)
        end
    else
        if ord == 0
            @inbounds C .= selectdim(A, 1, 1)
        elseif ord == 1
            @inbounds C .= selectdim(A, 1, 1) .+ (2x) .* selectdim(A, 1, 2)
        else
            x2 = 2x
            @inbounds Aᵢ = selectdim(A, 1, ord+1)
            sz = size(Aᵢ)
            s = zeros(T, sz)
            t = Array{T}(undef, sz)
            t .= 2 .* Aᵢ
            @inbounds for i ∈ ord-1:-1:1
                C .= t
                t .= x2 .* t .- s .+ 2 .* selectdim(A, 1, i+1)
                s .= C
            end
            @inbounds C .= x .* t .- s .+ selectdim(A, 1, 1)
        end
    end
    return C
end

_apply(::Evaluation{Nothing}, ::Chebyshev, ::Val, A::AbstractArray) = A
function _apply(ℰ::Evaluation, space::Chebyshev, ::Val{D}, A::AbstractArray{T,N}) where {D,T,N}
    x = ℰ.value
    CoefType = _coeftype(ℰ, space, T)
    ord = order(space)
    if iszero(x)
        C = convert(Array{CoefType,N-1}, selectdim(A, D, 1))
        @inbounds for i ∈ 2:2:ord
            C .+= ifelse(i%4 == 0, 2, -2) .* selectdim(A, D, i+1)
        end
        return C
    elseif isone(-x)
        C = convert(Array{CoefType,N-1}, selectdim(A, D, 1))
        @inbounds for i ∈ 1:ord
            C .+= ifelse(isodd(i), -2, 2) .* selectdim(A, D, i+1)
        end
        return C
    elseif isone(x)
        C = convert(Array{CoefType,N-1}, selectdim(A, D, 1))
        @inbounds for i ∈ 1:ord
            C .+= 2 .* selectdim(A, D, i+1)
        end
        return C
    else
        if ord == 0
            return @inbounds convert(Array{CoefType,N-1}, selectdim(A, D, 1))
        elseif ord == 1
            return @inbounds convert(Array{CoefType,N-1}, selectdim(A, D, 1) .+ (2x) .* selectdim(A, D, 2))
        else
            x2 = 2x
            @inbounds Aᵢ = selectdim(A, D, ord+1)
            sz = size(Aᵢ)
            s = zeros(CoefType, sz)
            t = Array{CoefType,N-1}(undef, sz)
            t .= 2 .* Aᵢ
            C = Array{CoefType,N-1}(undef, sz)
            @inbounds for i ∈ ord-1:-1:1
                C .= t
                t .= x2 .* t .- s .+ 2 .* selectdim(A, D, i+1)
                s .= C
            end
            @inbounds C .= x .* t .- s .+ selectdim(A, D, 1)
            return C
        end
    end
end

function _getindex(::Evaluation{Nothing}, ::Chebyshev, ::Chebyshev, ::Type{T}, i, j) where {T}
    if i == j
        return one(T)
    else
        return zero(T)
    end
end
function _getindex(ℰ::Evaluation, domain::Chebyshev, codomain::Chebyshev, ::Type{T}, i, j) where {T}
    if i == 0
        if j == 0
            return one(T)
        else
            x = ℰ.value
            if iszero(x)
                if isodd(j)
                    return zero(T)
                elseif j%4 == 0
                    return convert(T, 2)
                else
                    return convert(T, -2)
                end
            elseif isone(-x)
                if isodd(j)
                    return convert(T, -2)
                else
                    return convert(T, 2)
                end
            elseif isone(x)
                return convert(T, 2)
            else
                x2 = convert(T, 2x)
                if j == 1
                    return x2
                elseif j == 2
                    return convert(T, x2*x2 - 2)
                else
                    return x2*_getindex(ℰ, domain, codomain, T, i, j-1) - _getindex(ℰ, domain, codomain, T, i, j-2)
                end
            end
        end
    else
        return zero(T)
    end
end

# Cartesian spaces

image(ℰ::Evaluation, s::CartesianPower) =
    CartesianPower(image(ℰ, space(s)), nb_cartesian_product(s))

image(ℰ::Evaluation, s::CartesianProduct) =
    CartesianProduct(map(sᵢ -> image(ℰ, sᵢ), spaces(s)))

_coeftype(ℰ::Evaluation, s::CartesianPower, ::Type{T}) where {T} =
    _coeftype(ℰ, space(s), T)

_coeftype(ℰ::Evaluation, s::CartesianProduct, ::Type{T}) where {T} =
    @inbounds promote_type(_coeftype(ℰ, s[1], T), _coeftype(ℰ, Base.tail(s), T))
_coeftype(ℰ::Evaluation, s::CartesianProduct{<:Tuple{VectorSpace}}, ::Type{T}) where {T} =
    @inbounds _coeftype(ℰ, s[1], T)

function _apply!(c::Sequence{<:CartesianPower}, ℰ::Evaluation, a)
    @inbounds for i ∈ 1:nb_cartesian_product(space(c))
        _apply!(component(c, i), ℰ, component(a, i))
    end
    return c
end
function _apply!(c::Sequence{CartesianProduct{T}}, ℰ::Evaluation, a) where {N,T<:NTuple{N,VectorSpace}}
    @inbounds _apply!(component(c, 1), ℰ, component(a, 1))
    @inbounds _apply!(component(c, 2:N), ℰ, component(a, 2:N))
    return c
end
function _apply!(c::Sequence{CartesianProduct{T}}, ℰ::Evaluation, a) where {T<:Tuple{VectorSpace}}
    @inbounds _apply!(component(c, 1), ℰ, component(a, 1))
    return c
end

function _project!(C::LinearOperator{<:CartesianSpace,<:CartesianSpace}, ℰ::Evaluation)
    @inbounds for i ∈ 1:nb_cartesian_product(domain(C))
        _project!(component(C, i, i), ℰ)
    end
    return C
end


text:
# Code for solving assuming that all other cars have a single identical behavior model

mutable struct SingleBehaviorSolver <: Solver
    inner_solver::Solver
    behavior::BehaviorModel
end

mutable struct SingleBehaviorPolicy <: Policy
    inner_policy::Policy
    behavior::BehaviorModel
end

set_rng!(solver::SingleBehaviorSolver, rng::AbstractRNG) = set_rng!(solver.inner_solver, rng)
srand(p::SingleBehaviorPolicy, s) = srand(p.inner_policy, s)

function solve(solver::SingleBehaviorSolver, mdp::NoCrashProblem)
    dmodel = NoCrashIDMMOBILModel(mdp.dmodel, DiscreteBehaviorSet(BehaviorModel[solver.behavior], Weights([1.0])))
    single_behavior_mdp = NoCrashMDP{typeof(mdp.rmodel), DiscreteBehaviorSet}(dmodel, mdp.rmodel, mdp.discount, mdp.throw)

    inner_policy = solve(solver.inner_solver, single_behavior_mdp)
    return SingleBehaviorPolicy(inner_policy, solver.behavior)
end

function action_info(p::SingleBehaviorPolicy, s::Union{MLPhysicalState, MLState})
    as = actions(p.inner_policy.mdp, s)
    if length(as) == 1
        return first(as), Dict(:tree_queries=>missing, :search_time_us=>missing)
    end
    return action_info(p.inner_policy, single_behavior_state(s, p.behavior))
end

action_info(p::SingleBehaviorPolicy, agg::AggressivenessBelief) = action_info(p, agg.physical)

action(p::SingleBehaviorPolicy, s) = first(action_info(p, s))

function single_behavior_state(s::Union{MLState, MLPhysicalState}, behavior)
    new_cars = Vector{CarState}(undef, length(s.cars))
    for (i,c) in enumerate(s.cars)
        new_cars[i] = CarState(c.x, c.y, c.vel, c.lane_change, behavior, c.id)
    end
    return MLState(s, new_cars)
end


text:
"""
    AbstractNumbers

An abstract type to represent the number system on which a manifold is built.

This provides concrete number types for dispatch. The two most common number types are
the fields [`RealNumbers`](@ref) (`ℝ` for short) and [`ComplexNumbers`](@ref) (`ℂ`).
"""
abstract type AbstractNumbers end

"""
    RealNumbers <: AbstractNumbers
    ℝ = RealNumbers()

The field of real numbers.
"""
struct RealNumbers <: AbstractNumbers end

"""
    ComplexNumbers <: AbstractNumbers
    ℂ = ComplexNumbers()

The field of complex numbers.
"""
struct ComplexNumbers <: AbstractNumbers end

"""
    QuaternionNumbers <: AbstractNumbers
    ℍ = QuaternionNumbers()

The division algebra of quaternions.
"""
struct QuaternionNumbers <: AbstractNumbers end

const ℝ = RealNumbers()
const ℂ = ComplexNumbers()
const ℍ = QuaternionNumbers()

"""
    _unify_number_systems(𝔽s::AbstractNumbers...)

Compute a number system that includes all given number systems (as sub-systems) and is
closed under addition and multiplication.
"""
function _unify_number_systems(a::AbstractNumbers, rest::AbstractNumbers...)
    return _unify_number_systems(a, _unify_number_systems(rest...))
end
_unify_number_systems(𝔽::AbstractNumbers) = 𝔽
_unify_number_systems(r::RealNumbers, ::RealNumbers) = r
_unify_number_systems(::RealNumbers, c::ComplexNumbers) = c
_unify_number_systems(::RealNumbers, q::QuaternionNumbers) = q
_unify_number_systems(c::ComplexNumbers, ::RealNumbers) = c
_unify_number_systems(c::ComplexNumbers, ::ComplexNumbers) = c
_unify_number_systems(::ComplexNumbers, q::QuaternionNumbers) = q
_unify_number_systems(q::QuaternionNumbers, ::RealNumbers) = q
_unify_number_systems(q::QuaternionNumbers, ::ComplexNumbers) = q
_unify_number_systems(q::QuaternionNumbers, ::QuaternionNumbers) = q

Base.show(io::IO, ::RealNumbers) = print(io, "ℝ")
Base.show(io::IO, ::ComplexNumbers) = print(io, "ℂ")
Base.show(io::IO, ::QuaternionNumbers) = print(io, "ℍ")

@doc raw"""
    real_dimension(𝔽::AbstractNumbers)

Return the real dimension $\dim_ℝ 𝔽$ of the [`AbstractNumbers`](@ref) system `𝔽`.
The real dimension is the dimension of a real vector space with which a number in `𝔽` can be
identified.
For example, [`ComplexNumbers`](@ref) have a real dimension of 2, and
[`QuaternionNumbers`](@ref) have a real dimension of 4.
"""
function real_dimension(𝔽::AbstractNumbers)
    return error("real_dimension not defined for number system $(𝔽)")
end
real_dimension(::RealNumbers) = 1
real_dimension(::ComplexNumbers) = 2
real_dimension(::QuaternionNumbers) = 4

@doc raw"""
    number_system(M::Manifold{𝔽})

Return the number system the manifold `M` is based on, i.e. the parameter `𝔽`.
"""
number_system(M::Manifold{𝔽}) where {𝔽} = 𝔽


text:
export single_point_crossover,
    random_node_crossover,
    aligned_node_crossover,
    proportional_crossover,
    output_graph_crossover,
    subgraph_crossover,
    crossover

function single_point_crossover(c1::Chromosome, c2::Chromosome)
    # single point crossover
    cpoint = c1.nin + c1.nout + node_genes(c1) * rand(0:min(length(c1.nodes)-c1.nin,
                                                            length(c2.nodes)-c2.nin))
    if rand(Bool)
        ngenes = deepcopy([c1.genes[1:cpoint]; c2.genes[(cpoint+1):end]])
        return typeof(c1)(ngenes, c1.nin, c1.nout)
    else
        ngenes = deepcopy([c2.genes[1:cpoint]; c1.genes[(cpoint+1):end]])
        return typeof(c2)(ngenes, c2.nin, c2.nout)
    end
end

function random_inputs(c1::Chromosome, c2::Chromosome)
    # returns random input genes from c1 and c2 equally
    # not a full crossover operator
    child_inputs = zeros(c1.nin)
    parent = bitrand(c1.nin)
    child_inputs[parent] = c1.genes[1:c1.nin][parent]
    child_inputs[.~(parent)] = c2.genes[1:c1.nin][.~(parent)]
    child_inputs
end

function random_outputs(c1::Chromosome, c2::Chromosome)
    # returns random input genes from c1 and c2 equally
    # not a full crossover operator
    child_outputs = zeros(c1.nout)
    parent = bitrand(c1.nout)
    child_outputs[parent] = c1.genes[c1.nin+(1:c1.nout)][parent]
    child_outputs[.~(parent)] = c2.genes[c1.nin+(1:c1.nout)][.~(parent)]
    child_outputs
end

function random_node_crossover(c1::Chromosome, c2::Chromosome)
    # take random nodes from each parent equally, up to the size of the smaller parent
    min_nodes = min(length(c1.nodes)-c1.nin, length(c2.nodes)-c2.nin)
    p_nodes = c1.nin + (1:min_nodes)[bitrand(min_nodes)]
    p1_node_genes = get_genes(c1, (1:length(c1.nodes))[p_nodes])
    p2_node_genes = get_genes(c2, (1:length(c2.nodes))[p_nodes])
    genes = [random_inputs(c1, c2); random_outputs(c1, c2); p1_node_genes; p2_node_genes]
    typeof(c1)(genes, c1.nin, c1.nout)
end

function aligned_node_crossover(c1::Chromosome, c2::Chromosome)
    # align nodes based on position, then take from each parent equally
    p1_pos = get_positions(c1)
    p2_pos = get_positions(c2)
    min_nodes = min(length(c1.nodes)-c1.nin, length(c2.nodes)-c2.nin)
    p1_inds = c1.nin + collect(1:min_nodes)
    p1_nodes = Array{Int64}(0)
    p2_nodes = Array{Int64}(0)
    for node in c1.nin+(1:min_nodes)
        if rand() < 0.5
            i = indmin(abs.(p1_pos[p1_inds] - p2_pos[node]))
            append!(p1_nodes, [p1_inds[i]])
            p1_inds = deleteat!(p1_inds, i)
        else
            append!(p2_nodes, [node])
        end
    end
    genes = [random_inputs(c1, c2); random_outputs(c1, c2); get_genes(
        c1, p1_nodes); get_genes(c2, p2_nodes)]
    typeof(c1)(genes, c1.nin, c1.nout)
end

function proportional_crossover(c1::Chromosome, c2::Chromosome)
    genes = Array{Float64}(0)
    if length(c1.genes) == length(c2.genes)
        r = rand(length(c1.genes))
        genes = ((1 .- r) .* c1.genes) .+ (r .* c2.genes)
    elseif length(c1.genes) > length(c2.genes)
        r = rand(length(c2.genes))
        genes = ((1 .- r) .* c1.genes[1:length(c2.genes)]) .+ (r .* c2.genes)
        genes = [genes; c1.genes[(length(c2.genes)+1):end]]
    else
        r = rand(length(c1.genes))
        genes = ((1 .- r) .* c1.genes) .+ (r .* c2.genes[1:length(c1.genes)])
        genes = [genes; c2.genes[(length(c1.genes)+1):end]]
    end
    typeof(c1)(genes, c1.nin, c1.nout)
end

function output_graph_crossover(c1::Chromosome, c2::Chromosome)
    # split outputs equally between parents, then construct a child from their
    # corresponding input graphs
    p1_nodes = Array{Int64}(0)
    p2_nodes = Array{Int64}(0)
    output_genes = Array{Float64}(0)
    for output in 1:c1.nout
        if rand() < 0.5
            append!(p1_nodes, get_output_trace(c1, output))
            append!(output_genes, [c1.genes[c1.nin+output]])
        else
            append!(p2_nodes, get_output_trace(c2, output))
            append!(output_genes, [c2.genes[c2.nin+output]])
        end
    end
    p1_nodes = sort!(unique(p1_nodes))
    p2_nodes = sort!(unique(p2_nodes))
    input_genes = Array{Float64}(0)
    # take inputs from either parent trace, if they are in the parent traces
    for input in 1:c1.nin
        gene = c1.genes[input]
        if contains(==, p1_nodes, input)
            if contains(==, p2_nodes, input)
                if rand() < 0.5
                    gene = c2.genes[input]
                end
            end
        else
            if contains(==, p2_nodes, input)
                gene = c2.genes[input]
            else
                if rand() < 0.5
                    gene = c2.genes[input]
                end
            end
        end
        append!(input_genes, [gene])
    end
    p1_nodes = p1_nodes[p1_nodes .> c1.nin]
    p2_nodes = p2_nodes[p2_nodes .> c1.nin]
    genes = [input_genes; output_genes; get_genes(c1, p1_nodes); get_genes(c2, p2_nodes)]
    typeof(c1)(genes, c1.nin, c1.nout)
end

function subgraph_crossover(c1::Chromosome, c2::Chromosome)
    # Take subgraphs from both parents equally, adding all nodes of the chosen
    # subgraphs to the child.
    fc1 = forward_connections(c1)
    fc2 = forward_connections(c2)
    c1_nodes = []; c2_nodes = []
    for i in (c1.nin+1):min(length(c1.nodes), length(c2.nodes))
        if rand(Bool)
            append!(c1_nodes, fc1[i])
        else
            append!(c2_nodes, fc2[i])
        end
    end
    if length(c1.nodes) < length(c2.nodes)
        for i in length(c1.nodes):length(c2.nodes)
            if rand(Bool)
                append!(c2_nodes, fc2[i])
            end
        end
    elseif length(c2.nodes) < length(c1.nodes)
        for i in length(c2.nodes):length(c1.nodes)
            if rand(Bool)
                append!(c1_nodes, fc1[i])
            end
        end
    end
    c1_nodes = Array{Int64}(unique(intersect(collect((c1.nin+1):length(c1.nodes)), c1_nodes)))
    c2_nodes = Array{Int64}(unique(intersect(collect((c2.nin+1):length(c2.nodes)), c2_nodes)))
    genes = zeros(c1.nin+c1.nout)
    for i in 1:(c1.nin+c1.nout)
        if rand(Bool)
            genes[i] = c1.genes[i]
        else
            genes[i] = c2.genes[i]
        end
    end
    if length(c1_nodes) > 0
        genes = [genes; get_genes(c1, c1_nodes)]
    end
    if length(c2_nodes) > 0
        genes = [genes; get_genes(c2, c2_nodes)]
    end
    typeof(c1)(genes, c1.nin, c1.nout)
end

function crossover(c1::Chromosome, c2::Chromosome)
    eval(parse(string(Config.crossover_method)))(c1, c2)
end


text:
const MATH_INTRINSICS = GCNIntrinsic[]

for jltype in (
        #= TODO: Float16 Broken due to being i16 in Julia=#
        Float32, Float64)
    append!(MATH_INTRINSICS, GCNIntrinsic.((
        :sin, :cos, :tan, :asin, :acos, :atan, :atan2,
        :sinh, :cosh, :tanh, :asinh, :acosh, :atanh,
        :sinpi, :cospi, :tanpi,
        :asinpi, :acospi, :atanpi, :atan2pi,
        :sqrt, :rsqrt, :cbrt, :rcbrt, :recip,
        :log, :log2, :log10, :log1p, :logb, :ilogb,
        :exp, :exp2, :exp10, :expm1,
        :erf, :erfinv, :erfc, :erfcinv, :erfcx,
        # TODO: :brev, :clz, :ffs, :byte_perm, :popc,
        :isnormal, :nearbyint, :nextafter,
        :tgamma, :j0, :j1, :y0, :y1,
    ); inp_args=(jltype,), out_arg=jltype))

    push!(MATH_INTRINSICS, GCNIntrinsic(:sin_fast, :native_sin; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:cos_fast, :native_cos; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:sqrt_fast, :native_sqrt; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:rsqrt_fast, :native_rsqrt; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:recip_fast, :native_recip; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:log_fast, :native_log; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:log2_fast, :native_log2; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:log10_fast, :native_log10; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:exp_fast, :native_exp; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:exp2_fast, :native_exp2; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:exp10_fast, :native_exp10; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:abs, :fabs; inp_args=(jltype,), out_arg=jltype))
    # TODO: abs(::Union{Int32,Int64})

    # Multi-argument functions
    push!(MATH_INTRINSICS, GCNIntrinsic(:pow; inp_args=(jltype,jltype), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:pow, :pown; inp_args=(jltype,Union{UInt32,Int32}), out_arg=jltype))
    # TODO: push!(MATH_INTRINSICS, GCNIntrinsic(:pow, :pown; inp_args=(jltype,Union{UInt32,Int32}), out_arg=jltype))
    # TODO: :sincos, :frexp, :ldexp, :copysign,
    #push!(MATH_INTRINSICS, GCNIntrinsic(:ldexp; inp_args=(jltype,), out_arg=(jltype, Int32), isinverted=true))

    # Multi-output functions
    push!(MATH_INTRINSICS, GCNIntrinsic(:sincospi; inp_args=(jltype,), out_arg=jltype, isbroken=true))
end

for jltype in (Float32, Float64)
    push!(MATH_INTRINSICS, GCNIntrinsic(:isfinite; inp_args=(jltype,), out_arg=Int32))
    push!(MATH_INTRINSICS, GCNIntrinsic(:isinf; inp_args=(jltype,), out_arg=Int32))
    push!(MATH_INTRINSICS, GCNIntrinsic(:isnan; inp_args=(jltype,), out_arg=Int32))
    push!(MATH_INTRINSICS, GCNIntrinsic(:signbit; inp_args=(jltype,), out_arg=Int32))
end

for intr in MATH_INTRINSICS
    inp_vars = [gensym() for _ in 1:length(intr.inp_args)]
    inp_expr = [:($(inp_vars[idx])::$arg) for (idx,arg) in enumerate(intr.inp_args)]
    libname = Symbol("__$(intr.roclib)_$(intr.rocname)_$(intr.suffix)")
    @eval @inline function $(intr.jlname)($(inp_expr...))
        y = _intr($(Val(libname)), $(intr.out_arg), $(inp_expr...))
        return $(intr.isinverted ? :(1-y) : :y)
    end
end


text:
using Gen

function path_length(start::Point, path::Vector{Point})
    @assert length(path) > 0
    len = dist(start, path[1])
    for i in 1:length(path)-1
        len += dist(path[i], path[i+1])
    end
    return len
end

function compute_distances_from_start(path::Vector{Point})
    distances_from_start = Vector{Float64}(undef, length(path))
    distances_from_start[1] = 0.0
    for i=2:length(path)
        distances_from_start[i] = distances_from_start[i-1] + dist(path[i-1], path[i])
    end
    return distances_from_start
end

function walk_path(path::Vector{Point}, speed::Float64, T::Int)
    distances_from_start = compute_distances_from_start(path)
    times = collect(0.0:1.0:(T-1))
    locations = Vector{Point}(undef, T)
    for (time_idx, t) in enumerate(times)
        desired_distance = t * speed
        used_up_time = false
        # NOTE: can be improved (iterate through path points along with times)
        for i=2:length(path)
            prev = path[i-1]
            cur = path[i]
            dist_to_prev = dist(prev, cur)
            if distances_from_start[i] >= desired_distance
                # we overshot, the location is between i-1 and i
                overshoot = distances_from_start[i] - desired_distance
                @assert overshoot <= dist_to_prev
                past_prev = dist_to_prev - overshoot
                frac = past_prev / dist_to_prev
                locations[time_idx] = Point(prev.x * (1. - frac) + cur.x * frac,
                                     prev.y * (1. - frac) + cur.y * frac)
                used_up_time = true
                break
            end
        end
        if !used_up_time
            # sit at the goal indefinitely
            locations[time_idx] = path[end]
        end
    end
    return locations
end

struct ObsModelParams
    nominal_speed::Float64 # distance per time step
    walk_noise::Float64
    noise::Float64
end

prob_lag(params::ObsModelParams) = params.walk_noise/2
prob_skip(params::ObsModelParams) = params.walk_noise/2
prob_normal(params::ObsModelParams) = 1-params.walk_noise

@inline function noise_log_likelihood(params::ObsModelParams, points_along_path_pt::Point, obs_pt::Point)
    ll = logpdf(normal, obs_pt.x, points_along_path_pt.x, params.noise)
    ll += logpdf(normal, obs_pt.y, points_along_path_pt.y, params.noise)
    return ll
end

@inline function populate_new_alpha!(new_alpha, prev_alpha, points_along_path::Vector{Point}, obs_t::Point, params::ObsModelParams)
    for k in 1:length(points_along_path)
        # observation model: obs[t] corresponds with points_along_path[k]
        noise_ll = noise_log_likelihood(params, points_along_path[k], obs_t)
        # dynamics model: given that we are k, the previous could have been
        # (i) k - 2 (skip one), (ii) k - 1 (advance as usual), or (iii) k
        # (don't advance, lag)
        if k > 2
            val1 = prev_alpha[k-2] + log(prob_skip(params))
        else
            val1 = -Inf
        end
        if k > 1
            val2 = prev_alpha[k-1] + log(prob_normal(params))
        else
            val2 = -Inf
        end
        val3 = prev_alpha[k] + log(prob_lag(params))
        if val1 == -Inf && val2 == -Inf && val3 == -Inf
            new_alpha[k] = -Inf
            continue
        end
        max_val = max(val1, val2, val3)
        new_alpha[k] = noise_ll + max_val + log(exp(val1 - max_val) + exp(val2 - max_val) + exp(val3 - max_val))
    end
    return nothing
end

@inline function forward_transition_log_probs(next_z::Int, params::ObsModelParams, K::Int)
    log_probs = fill(-Inf, K)
    if next_z == 1
        log_probs[1] = log(prob_lag(params))
    elseif next_z == 2
        log_probs[1] = log(prob_normal(params))
        log_probs[2] = log(prob_lag(params))
    else
        log_probs[next_z] = log(prob_lag(params))
        log_probs[next_z-1] = log(prob_normal(params))
        log_probs[next_z-2] = log(prob_skip(params))
    end
    return log_probs
end

@inline function populate_initial_alpha!(alpha, params::ObsModelParams, points_along_path::Vector{Point}, obs1::Point)
    if length(alpha) == 1
        alpha[1] = noise_log_likelihood(params, points_along_path[1], obs1)
        return
    end
    fill!(alpha, -Inf)
    first_prob_normal = prob_lag(params) + prob_normal(params)
    first_prob_skip = prob_skip(params)
    alpha[1] = log(first_prob_normal) + noise_log_likelihood(params, points_along_path[1], obs1) # p(z1) is deterministic at start of points_along_path
    alpha[2] = log(first_prob_skip) + noise_log_likelihood(params, points_along_path[2], obs1) # p(z1) is deterministic at start of points_along_path
    return nothing
end

function run_forward_backward(params::ObsModelParams, points_along_path::Vector{Point}, obs::Vector{Point})
    K = length(points_along_path)
    T = length(obs)

    # run forward filtering backward sampling

    # forward filtering
    alphas = Matrix{Float64}(undef, T, K) # optimize memory access?
    populate_initial_alpha!(view(alphas, 1, :), params, points_along_path, obs[1])
    for t in 2:T
        populate_new_alpha!(view(alphas, t, :), view(alphas, t-1, :), points_along_path, obs[t], params)
    end
    log_marginal_likelihood = logsumexp(alphas[end,:])

    # backward sampling
    alignment = Vector{Int}(undef, length(points_along_path))
    ldist = alphas[T,:]
    dist = exp.(ldist .- logsumexp(ldist))
    alignment[T] = categorical(dist)
    for t in T-1:-1:1
        ldist = alphas[t,:] .+ forward_transition_log_probs(alignment[t+1], params, K)
        dist = exp.(ldist .- logsumexp(ldist))
        alignment[t] = categorical(dist)
    end
    return (log_marginal_likelihood, alignment)
end

#function log_marginal_likelihood(params::ObsModelParams, points_along_path::Vector{Point}, obs::Vector{Point})
    #K = length(points_along_path)
    #T = length(obs)
    #alpha = Vector{Float64}(undef, K)
    #new_alpha = Vector{Float64}(undef, K)
    #populate_initial_alpha!(alpha, params, points_along_path, obs[1])
    #for t in 2:T
        #populate_new_alpha!(new_alpha, alpha, points_along_path, obs[t], params)
        #tmp = alpha; alpha = new_alpha; new_alpha = tmp
    #end
    #return logsumexp(alpha)
#end

#######################
# generative function #
#######################

struct ObsModelTrace <: Trace
    gen_fn::GenerativeFunction
    path::Vector{Point}
    params::ObsModelParams
    points_along_path::Vector{Point}
    obs::Vector{Point}
    lml::Float64
    alignment::Vector{Int}
end

Gen.get_gen_fn(tr::ObsModelTrace) = tr.gen_fn
Gen.get_retval(trace::ObsModelTrace) = (trace.points_along_path, trace.alignment)
Gen.get_args(tr::ObsModelTrace) = (tr.path, tr.params, length(tr.obs))
Gen.get_score(tr::ObsModelTrace) = tr.lml
Gen.project(tr::ObsModelTrace, ::EmptySelection) = 0.0

function Gen.get_choices(tr::ObsModelTrace)
    cm = choicemap()
    for (i, pt) in enumerate(tr.obs)
        cm[(:x, i)] = pt.x
        cm[(:y, i)] = pt.y
    end
    return cm
end

struct ObsModel <: GenerativeFunction{Tuple{Vector{Point},Vector{Int}},ObsModelTrace} end

const path_observation_model = ObsModel()

function noise_sample(noise, pt::Point)
    x = normal(pt.x, noise)
    y = normal(pt.y, noise)
    return Point(x, y)
end

function sample_from_obs_model(params::ObsModelParams, points_along_path::Vector{Point}, T::Int)
    obs = Vector{Point}(undef, T)
    if length(points_along_path) == 1
        z_cur = 1
    else
        z_cur = bernoulli(prob_lag(params) + prob_normal(params)) ? 1 : 2
    end
    obs[1] = noise_sample(params.noise, points_along_path[z_cur])
    for t in 2:T
        _prob_skip = (z_cur < length(points_along_path) - 1) ? prob_skip(params) : 0.0
        _prob_normal = (z_cur < length(points_along_path)) ? prob_normal(params) : 0.0
        probs = [prob_lag(params), _prob_normal, _prob_skip]
        probs /= sum(probs)
        dowhat = categorical(probs)
        if dowhat == 1
            # lag
        elseif dowhat == 2
            # normal
            z_cur += 1
        else
            # skip
            z_cur += 2
        end
        obs[t] = noise_sample(params.noise, points_along_path[z_cur])
    end
    return obs
end

function Gen.simulate(gen_fn::ObsModel, args::Tuple)
    path, params, T = args
    points_along_path = walk_path(path, params.nominal_speed, T)
    obs = sample_from_obs_model(params, points_along_path, T)
    (lml, alignment) = run_forward_backward(params, points_along_path, obs)
    @assert !isnan(lml)
    return ObsModelTrace(gen_fn, path, params, points_along_path, obs, lml, alignment)
end

function Gen.generate(gen_fn::ObsModel, args::Tuple, constraints::ChoiceMap)
    path, params, T = args
    points_along_path = walk_path(path, params.nominal_speed, T)
    obs = Vector{Point}(undef, T)
    if isempty(constraints)
        trace = simulate(gen_fn, args)
        return (trace, 0.0, nothing)
    end
    for t in 1:T
        obs[t] = Point(constraints[(:x, t)], constraints[(:y, t)])
    end
    (lml, alignment) = run_forward_backward(params, points_along_path, obs)
    @assert !isnan(lml)
    trace = ObsModelTrace(gen_fn, path, params, points_along_path, obs, lml, alignment)
    retval = get_retval(trace)
    return (trace, lml, retval)
end

function Gen.update(tr::ObsModelTrace, args::Tuple, argdiffs::Tuple, constraints::ChoiceMap)
    old_path, old_params, old_T = get_args(tr)
    new_path, new_params, new_T = args
    old_points_along_path = tr.points_along_path
    new_points_along_path = walk_path(new_path, new_params.nominal_speed, new_T)
    if (new_T == old_T + 1) && has_value(constraints, (:x, new_T)) && has_value(constraints, (:y, new_T))
        obs = copy(tr.obs)
        push!(obs, Point(constraints[(:x, new_T)], constraints[(:y, new_T)]))
    elseif new_T == old_T && isempty(constraints)
        obs = tr.obs
    else
        error("not implemented")
    end
    (lml, alignment) = run_forward_backward(new_params, new_points_along_path, obs)
    @assert !isnan(lml)
    new_trace = ObsModelTrace(get_gen_fn(tr), new_path, new_params, new_points_along_path, obs, lml, alignment)
    return (new_trace, lml - tr.lml, UnknownChange(), EmptyChoiceMap())
end

function Gen.regenerate(tr::ObsModelTrace, args::Tuple, argdiffs::Tuple, selection::Selection)
    if !isempty(selection)
        error("not implemented")
    end
    (new_trace, weight, retdiff, _) = update(tr, args, argdiffs, EmptyChoiceMap())
    return (new_trace, weight, retdiff)
end

export ObsModelParams, path_observation_model


text:
"""
    FortranFilesError(msg)

Exception thrown if an operation on a FortranFile fails or cannot be performed.
Argument `msg` is a description of the failure mode.
"""
struct FortranFilesError <: Exception
   msg::String
end

fthrow(msg) = throw(FortranFilesError(msg))

Base.showerror(io::IO, exc::FortranFilesError) = print(io, "FortranFilesError: $(exc.msg)")


text:
export DSTOp

"""
  DSTOp(T::Type, shape::Tuple)

returns a `LinearOperator` which performs a DST on a given input array.

# Arguments:
* `T::Type`       - type of the array to transform
* `shape::Tuple`  - size of the array to transform
"""
function DSTOp(T::Type, shape::Tuple)

  return LinearOperator(prod(shape), prod(shape), true, false
            , x->vec(FFTW.r2r(reshape(x,shape),FFTW.RODFT10)).*weights(shape)
            , nothing
            , y->vec(FFTW.r2r(reshape(y ./ weights(shape) ,shape),FFTW.RODFT01)) ./ (8*prod(shape))  )
end

function weights(s)
  w = ones(s...)./sqrt(8*prod(s))
  w[s[1],:,:]./=sqrt(2)
  w[:,s[2],:]./=sqrt(2)
  w[:,:,s[3]]./=sqrt(2)
  return reshape(w,prod(s))
end


text:
import REPL
Base.atreplinit() do repl
    # make the ; shell mode sticky
    !isdefined(repl, :interface) && (repl.interface = REPL.setup_interface(repl))
    repl.interface.modes[2].sticky=true
end

try
    @eval using Revise
    # Turn on Revise's automatic-evaluation behavior
    Revise.async_steal_repl_backend()
catch err
    @warn "Could not load Revise."
end

max_stars_repo_name: Shunsuke-Hori/lecture-source-jl

max_issues_repo_name: Shunsuke-Hori/lecture-source-jl

max_forks_repo_name: Shunsuke-Hori/lecture-source-jl


text:
const Player = Bool
const WHITE = true
const BLACK = false

const Cell = Union{Nothing,Player}
function get_neighbours(board, pos)
    #@debug("getting neighbours of cell $pos, namely $(board[pos])")
    #global BOARD_SIDE
    neighbours = Cell[]
    plus_one = pos + 1
    minus_one = pos - 1
    plus_side = pos + BOARD_SIDE
    minus_side = pos - BOARD_SIDE
    try
        push!(neighbours, board[minus_one])
        #@debug("minus_one = $(board[minus_one])")
    catch

    end
    try
        push!(neighbours, board[plus_side])
    catch
        #@debug("no right neighbour")
    end
    try
        push!(neighbours, board[plus_one])
    catch e
        #@debug("no bottom neighbour, error is $e")
    end
    try
        push!(neighbours, board[minus_side])
    catch
        #@debug("no left neighbour")
    end
    return neighbours
end
function neighbour_indices(array, pos)
    if len(size(array)) == 1
        array = reshape(array, Int(sqrt(length(array))))
    end
    side = size(array)[1]
    indices = Int16[pos-1, pos+side, pos+1, pos-side]
    for (index, content) in enumerate(indices)
        if content
            not in 1:length(array)
            delete!(indices, index)
        end
    end
end
a = rand(Bool, 5, 5)
const BOARD_SIDE = size(a)[1]
function get_first_chain(board)
    basecolor = first(board)
    chain = Int16[]

end
function get_chains(args)
    body
end


text:
module GnuplotScripting

include("gnuplotscripting.jl")
include("export.jl")

end


text:
using Makie

scene = Scene(raw = true, camera = cam2d!, resolution = (300, 300))

center!(scene)

b1 = button("<-  ", textsize = 90, position = (50, 50))
b2 = button("->", textsize = 90, position = (50, 50))

current_status = 0

on(b1[end][:clicks]) do c
    global current_status
    current_status -= c
    @show current_status
end

on(b2[end][:clicks]) do c
    global current_status
    current_status += c
    @show current_status
end

RecordEvents(hbox(b1, b2, scene), "output")

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
mutable struct Sequence{T<:Real}
    sequenceComponents      ::Array{String}
    sequenceValues          ::Array{<:Real}
    value                   ::T
end

Sequence() = Sequence(String[],Float64[],0)

"""
    addToSequence!(seq::Sequence, decVar::String, decision::T, addedValue::R) where {R<:Real,T<:Real}

Add a new part to a stored sequence

# Arguments
- `seq::Sequence`: the sequence to be added to
- `decVar::String`: the sequence member being decided
- `decision::T`: the value of the sequence component
- `addedValue::R`:the cost/benefit of adding the new part to the sequence
"""
function addToSequence!(seq::Sequence, decVar::String, decision::T, addedValue::R) where {R<:Real,T<:Real}
    push!(seq.sequenceComponents, decVar)
    push!(seq.sequenceValues, decision)
    seq.value = seq.value + addedValue
end

"""
    deepCopySequence(seq::Sequence)

Returns a deep copy of a sequence

# Arguments
- `seq::Sequence`: the sequence to be deep copied
"""
function deepCopy(seq::Sequence)
    return Sequence(copy(seq.sequenceComponents),copy(seq.sequenceValues),seq.value)
end

"""
    sequenceToDict(seq::Sequence)

Returns a Dict version of the sequence object

# Arguments
- `seq::Sequence`: the sequence to be translated
"""
function sequenceToDict(seq::Sequence)
    sequence = Dict{String, Float64}()
    for i in 1:length(seq.sequenceComponents)
        sequence[seq.sequenceComponents[i]]=seq.sequenceValues[i]
    end
    return sequence
end

"""
    findBestSequence(list::Array{Sequence},largestValue::Bool)

Returns the best sequence from a list as defined by the smallest or largest value

# Arguments
- `list::Array{Sequence}`: the list to choose a sequence from
- `largestValue::Bool`: return the sequence with the largest value if true, smallest value if false
"""
function findBestSequence(list::Array{Sequence},largestValue::Bool)
    returnSeq = nothing
    if largestValue
        for seq in list
            if isnothing(returnSeq) || seq.value>returnSeq.value
                returnSeq = seq
            end
        end
    else
        for seq in list
            if isnothing(returnSeq) || seq.value<returnSeq.value
                returnSeq = seq
            end
        end
    end
    return returnSeq
end

"""
    getSequenceLength(seq::Sequence)

Returns the length of a sequence

# Arguments
- `seq::Sequence`: the sequence
"""
function getSequenceLength(seq::Sequence)
    return length(seq.sequenceComponents)
end

"""
    getLastDecision(seq::Sequence)

Returns the value of the last decision in the sequence

# Arguments
- `seq::Sequence`: the sequence
"""
function getLastDecision(seq::Sequence)
    return seq.sequenceValues[length(seq.sequenceValues)]
end

"""
    getValues(seq::Sequence)

Get the values in the sequence in order

# Arguments
- `seq::Sequence`: the sequence
"""
function getValues(seq::Sequence)
    return seq.sequenceValues
end
"""
    getValue(seq::Sequence)

Get the value of the sequence

# Arguments
- `seq::Sequence`: the sequence
"""
function getValue(seq::Sequence)
    return seq.value
end


text:
module StarPUClang

include("LibClang.jl")
using .LibClang

using DataStructures

include("cltypes.jl")
export CLCursor, CLType, CLToken
foreach(names(@__MODULE__, all=true)) do s
    x = getfield(@__MODULE__, s)
    if x isa DataType && supertype(x) in [CLCursor, CLType, CLToken]
        @eval export $s
    end
end

include("index.jl")
export Index

include("trans_unit.jl")
export TranslationUnit, spelling, getcursor, parse_header, parse_headers

include("cursor.jl")
export isnull, isdecl, isref, isexpr, isstmt, isattr, hasattr
export ispreprocessing, isunexposed, is_translation_unit, isfunctionlike, isbuiltin
export isinlined, isbit, isdef, isvariadic, is_typedef_anon
export get_translation_unit, get_semantic_parent, get_lexical_parent, get_included_file
export getref, getdef
export kind, name, spelling, type, extent, value, location, filename
export canonical, underlying_type, integer_type, result_type, return_type, typedef_type
export bitwidth, argnum, argument, function_args
export search, children
export is_forward_declaration

include("type.jl")
export isvolatile, isrestrict, isvariadic, is_plain_old_data
export address_space, typedef_name, typedecl
export pointee_type, argtype, element_type, element_num
export resolve_type, get_named_type

include("token.jl")
export TokenList, tokenize, annotate

include("show.jl")

include("utils.jl")
export name_safe, symbol_safe
export copydeps, print_template
export dumpobj
export find_std_headers

include("clang2julia.jl")
export clang2julia, target_type, typesize

include("expr_unit.jl")
export ExprUnit, dump_to_buffer, print_buffer

include("context.jl")
export AbstractContext, DefaultContext
export parse_header!, parse_headers!

include("wrap_c.jl")
export wrap!

include("compat.jl")
export WrapContext, init

function version()
    cxstr = clang_getClangVersion()
    ptr = clang_getCString(cxstr)
    s = unsafe_string(ptr)
    clang_disposeString(cxstr)
    return s
end

const LLVM_VERSION = match(r"[0-9]+.[0-9]+.[0-9]+", version()).match
const LLVM_LIBDIR = dirname(LibClang.LLVM_jll.libclang_path)
const LLVM_INCLUDE = joinpath(LLVM_LIBDIR, "clang", LLVM_VERSION, "include")
const CLANG_INCLUDE = LLVM_INCLUDE

export LLVM_VERSION, LLVM_LIBDIR, LLVM_INCLUDE, CLANG_INCLUDE

end


text:
#TODO: analyze_into_waveguides account for phase, for now it's not accounted for at all
#TODO: waveguide analysis only for perpendicular waveguides
################################################################################
### ANALYZE FIELD
################################################################################
#@code_warntype checked
"""
    c = analyze_output(sim, k, ψ, m)

c is the output coefficient in the mth channel

S is constructed from cs for unit input on each channel
"""
function analyze_field(sim::Simulation, k, ψ, m; direction=:out)
    if isempty(sim.sys.waveguides)
        c = analyze_into_angular_momentum(sim, k, ψ, m, direction)
    elseif any(ispc.(sim.sys.domains[get_waveguide_domains(sim,sim.sct.channels[m].waveguide)]))
        c = analyze_into_pc_waveguide(sim, k, ψ, m, direction)
    elseif (sim.sys.domains[ScalarFDFD.get_waveguide_   domains(sim, sim.sct.channels[m].waveguide)][end].domain_type) == :halfspace_waveguide
        c = analyze_into_halfspace_waveguide(sim, k, ψ, m, direction)
    elseif (sim.sys.domains[ScalarFDFD.get_waveguide_domains(sim, sim.sct.channels[m].waveguide)][end].domain_type) == :planar_waveguide
        c = analyze_into_planar_waveguide(sim, k, ψ, m, direction)
    else
        throw(ArgumentError("unrecognized waveguide type"))
    end

    return c
end

#
"""
    c = analyze_into_waveguides_pc(sim, k, ψ, m)
"""
function analyze_into_pc_waveguide(sim::Simulation, k, ψ, m, direction)

    xs = sim.dis.x[1][1] .+ (0:sim.dis.N[1]-1)*sim.dis.dx[1]
    ys = sim.dis.x[2][1] .+ (0:sim.dis.N[2]-1)*sim.dis.dx[2]
    ψ_itp = CubicSplineInterpolation((xs,ys),reshape(ψ[:,1],sim.dis.N[1],sim.dis.N[2]), bc=Periodic(OnCell()), extrapolation_bc=Periodic(OnCell()))

    u, _ = incident_pc_mode(sim, k, m, direction)
    u_itp = CubicSplineInterpolation((xs,ys),reshape(u[:,1],sim.dis.N[1],sim.dis.N[2]), bc=Periodic(OnCell()), extrapolation_bc=Periodic(OnCell()))

    wg_sim = ScalarFDFD.extract_waveguide_simulation(sim, sim.sct.channels[m].waveguide)

    if ScalarFDFD.get_asymptote(sim, sim.sct.channels[m].waveguide) == :left
        wg_sim.bnd.∂Ω[1,1] = sim.bnd.∂Ω_tr[1,1]
        wg_sim.bnd.∂Ω[2,1] = sim.bnd.∂Ω_tr[1,1] + wg_sim.lat.a
    elseif ScalarFDFD.get_asymptote(sim, sim.sct.channels[m].waveguide) == :right
        wg_sim.bnd.∂Ω_tr[2,1] = sim.bnd.∂Ω_tr[2,1]
        wg_sim.bnd.∂Ω_tr[1,1] = sim.bnd.∂Ω_tr[2,1] - wg_sim.lat.a
    elseif ScalarFDFD.get_asymptote(sim, sim.sct.channels[m].waveguide) == :bottom
        wg_sim.bnd.∂Ω_tr[1,2] = sim.bnd.∂Ω_tr[1,2]
        wg_sim.bnd.∂Ω_tr[2,2] = sim.bnd.∂Ω_tr[1,2] + wg_sim.lat.b
    elseif ScalarFDFD.get_asymptote(sim, sim.sct.channels[m].waveguide) == :top
        wg_sim.bnd.∂Ω_tr[2,2] = sim.bnd.∂Ω_tr[2,2]
        wg_sim.bnd.∂Ω_tr[1,2] = sim.bnd.∂Ω_tr[2,2] - wg_sim.lat.b
    else
        throw(ArgumentError("invalid asymptote $(get_asymptote(sim,sim.sct.channels[i].waveguide)) for channel $i"))
    end

    wg_sim2 = Simulation(wg_sim)
    xs = wg_sim2.dis.x[1] .+ wg_sim2.dis.dx[1]
    ys = wg_sim2.dis.x[2] .+ wg_sim2.dis.dx[2]

    c = quadrature(wg_sim2, (ψ_itp.(xs,ys).*conj(u_itp.(xs,ys)))[:]; weight=:ε)/quadrature(wg_sim2, abs2.(u_itp.(xs,ys))[:]; weight=:ε)

    return c::ComplexF64
end

"""
    c = analyze_into_waveguides_planar(sim, k, ψ, m, direction; η_init=-3)
"""
function analyze_into_planar_waveguide(sim::Simulation, k, ψ, m, direction; η_init=-3)

    side = sim.bnd.waveguides[sim.sct.channels[m].waveguide].side
    xs = sim.dis.x[1][1]:sim.dis.dx[1]:sim.dis.x[1][end]
    ys = sim.dis.x[2][1]:sim.dis.dx[2]:sim.dis.x[2][end]
    Ψ = CubicSplineInterpolation((xs,ys), reshape(ψ,sim.dis.N[1],sim.dis.N[2]))

    β, u_itp, collapsed_dim, full_dim = transverse_mode(sim, k, m)
    φ = u_itp(sim.dis.x[full_dim])
    if side == :l
        P = Ψ(sim.dis.x[1][1],sim.dis.x[2])
    elseif side == :r
        P = Ψ(sim.dis.x[1][end],sim.dis.x[2])
    elseif side == :b
        P = Ψ(sim.dis.x[1],sim.dis.x[2][1])
    elseif side == :t
        P = Ψ(sim.dis.x[1],sim.dis.x[2][end])
    else
        throw(ErrorException("Invalid side $(side)"))
    end

    c = sqrt(real(β))*quadrature((conj(φ).*P)[:], sim.dis.dx[full_dim])
    return c
end

#@code_warntype checked
"""
    c = analyze_into_waveguides_halfspace(sim,k,m)
"""
function analyze_into_halfspace_waveguide(sim, k, ψ, m)

    _,_,_,β = incident_mode_halfspace(sim, k, m)

    if get_asymptote(sim, sim.sct.channels[m].waveguide) == :left
        ind = sim.dis.X_idx[1]
    elseif get_asymptote(sim, sim.sct.channels[m].waveguide) == :right
        ind = sim.dis.X_idx[end]
    elseif get_asymptote(sim, sim.sct.channels[m].waveguide) == :bottom
        ind = sim.dis.X_idx[1]
    elseif get_asymptote(sim, sim.sct.channels[m].waveguide) == :top
        ind = sim.dis.X_idx[end]
    else
        throw(ArgumentError("invalid asymptote $(get_asymptote(sim,sim.sct.channels[i].waveguide)) for channel $i"))
    end
    c = ψ[ind]*β
end

"""
    c = analyze_into_angular_momentum(sim, k, ψ, m, direction)
"""
function analyze_into_angular_momentum(sim::Simulation, k, ψ, m, direction)

    nθ = ANALYSIS_QUADRATURE_POINTS
    θ = LinRange(0,2π,nθ)
    dθ = θ[2]-θ[1]

    # R is radius at which to interpolate
    R = (findmin(abs.(sim.bnd.∂Ω_tr)))[1] - findmax(sim.dis.dx)[1]
    X = R*cos.(θ)
    Y = R*sin.(θ)

    # interpolate wavefunction at r=R, result is P(θ)
    xs = sim.dis.x[1][1]:sim.dis.dx[1]:sim.dis.x[1][end]
    ys = sim.dis.x[2][1]:sim.dis.dx[2]:sim.dis.x[2][end]
    Ψ = CubicSplineInterpolation((xs,ys), reshape(ψ,sim.dis.N[1],sim.dis.N[2]))
    P = Ψ.(X,Y)

    q = sim.sct.channels[m].quantum_number

    if direction == :in
        c = quadrature(exp.(-1im*q*θ).*P, dθ)/(π*hankelh2(q,k*R))
    elseif direction == :out
        c = quadrature(exp.(-1im*q*θ).*P, dθ)/(π*hankelh1(q,k*R))
    else
        throw(ErrorException("Invalid direction $(direction). Must be one of :in or :out"))
    end

    return c
end

################################################################################
### ABSORPTION AND FLUX
################################################################################
#@code_warntype checked
"""
    flux, (left, right, bottom, top) = surface_flux(sim, ψ)
"""
function surface_flux(sim::Simulation,Ψ)

    flux = Array{Float64}(undef,size(Ψ,2))
    top = Array{Float64}(undef,size(flux))
    bottom = Array{Float64}(undef,size(flux))
    left = Array{Float64}(undef,size(flux))
    right = Array{Float64}(undef,size(flux))

    for i in 1:size(Ψ,2)
        ψ = reshape(Ψ[sim.dis.X_idx,i],sim.dis.N_tr[1],sim.dis.N_tr[2])
        kx = imag(conj(ψ[1:end-1,:]).*ψ[2:end,:]/sim.dis.dx[1])
        ky = imag(conj(ψ[:,1:end-1]).*ψ[:,2:end]/sim.dis.dx[2])
        if sim.dis.N[1]==1
            bottom[i] = ky[1,1]
            top[i] = ky[1,end]
            left[i]  = 0
            right[i]  = 0
        elseif sim.dis.N[2]==1
            bottom[i] = 0
            top[i] = 0
            left[i]  = kx[1,1]
            right[i]  = kx[end,1]
        else
            slice = 2
            bottom[i] = sum(ky[slice+1:end-slice+1, slice]              )*sim.dis.dx[1]
            top[i] =    sum(ky[slice+1:end-slice+1, end-slice+1]        )*sim.dis.dx[1]
            left[i]  =  sum(kx[slice              , slice+1:end-slice+1])*sim.dis.dx[2]
            right[i]  = sum(kx[end-slice+1        , slice+1:end-slice+1])*sim.dis.dx[2]
        end
        flux[i] = top[i] - bottom[i] + right[i] - left[i]
    end

    return flux, (left,right,bottom,top)
end

#@code_warntype checked
"""
    absorption = compute_loss(sim, k, ψ)
"""
function bulk_absorption(sim::Simulation, k, ψ)

    absorption = Array{Float64}(undef, length(k))

    for i ∈ eachindex(k)
        A = imag(sim.sys.ε[:]*k[i]^2).*abs2.(ψ[:,i])
        absorption[i] = quadrature(sim, A)
    end

    return absorption
end


text:
# Currently, can't run io.jl and concurrency.jl at the same time as
# multiprocessing doesn't work with --compilecache=no
# needed for the mocking tests to working.

@testset "Concurrency" begin
    @testset "Async Logging" begin
        Memento.reset!()
        io = IOBuffer()

        try
            setlevel!(getlogger(), "info")
            push!(
                getlogger(),
                DefaultHandler(io, DefaultFormatter("{name} - {level}: {msg}"))
            )

            asyncmap(x -> Memento.warn(getlogger(), "message"), 1:10)
            all_msgs = split(String(take!(io)), '\n')

            @test !isempty(all_msgs)
            @test all(m -> m == all_msgs[1], all_msgs[2:end-1])
        finally
            close(io)
        end
    end
    @testset "Parallel Logging" begin
        Memento.reset!()
        numprocs = addprocs(2)

        try
            # broadcast a filename to write to.
            @eval @everywhere parallel_test_filename = $(tempname())

            @everywhere using Memento
            @everywhere setlevel!(getlogger(), "info")
            @everywhere push!(
                getlogger(),
                DefaultHandler(parallel_test_filename, DefaultFormatter("{name} - {level}: {msg}"))
            )

            pmap(x -> Memento.warn(getlogger(), "message"), 1:10)

            open(parallel_test_filename) do f
                all_msgs = readlines(f)
                @test !isempty(all_msgs)
                @test all(m -> m == all_msgs[1], all_msgs[2:end])
            end

            @test isfile(parallel_test_filename)
            @test Sys.iswindows() ? true : success(`rm $parallel_test_filename`)   # Get UNLINK Error on windows
        finally
            rmprocs(numprocs)
        end
    end
end