text:
struct TypeBind
  name::Symbol
  ts::Set{Any}
end

istb(s) = false
istb(s::Symbol) = !(endswith(string(s), "_") ||
                    endswith(string(s), "_str")) &&
                    occursin("_", string(s))

tbname(s::Symbol) = Symbol(split(string(s), "_")[1])
tbname(s::TypeBind) = s.name

totype(s::Symbol) = string(s)[1] in 'A':'Z' ? s : Expr(:quote, s)

function tbnew(s::Symbol)
  istb(s) || return s
  ts = map(Symbol, split(string(s), "_"))
  name = popfirst!(ts)
  ts = map(totype, ts)
  Expr(:$, :(MacroTools.TypeBind($(Expr(:quote, name)), Set{Any}([$(ts...)]))))
end

match_inner(b::TypeBind, ex, env) =
  isexpr(ex, b.ts...) ? (env[tbname(b)] = ex; env) : @nomatch(b, ex)

subtb(s) = s
subtb(s::Symbol) = tbnew(s)
subtb(s::Expr) = isexpr(s, :line) ? s : Expr(subtb(s.head), map(subtb, s.args)...)


text:
#include("build_julia_specific_conda.jl")
include("pyimports.jl")

ext: jl


text:
module Lsode
#This is a wrap of the DLSODE function from the odepack FORTRAN library.
#I created this mainly to try out the fortran-calling functionality in Julia.
#For solving differential equations in Julia check out the DifferentialEquations package.

#TODO Extend the `IOPT` functionality parameter
#TODO Copy the ITOL behaviour from the odepack.f
#TODO Make macros more flexible -- e.g. other symbol use

export ode,@diff_eq,@diff_eq_jac,chem2diff

include("macros.jl")
include("chem2diff.jl")

const istate_dict = Dict{Int64,String}(
                                       1 =>"Nothing was done, because TOUT was equal to T.",
                                       2 =>"DLSODE was successful (otherwise, negative).",
                                       -1 =>"Excess work done on this call (perhaps wrong MF).",
                                       -2 =>"Excess accuracy requested (tolerances too small).",
                                       -3 =>"Illegal input detected (see printed message).",
                                       -4 =>"Repeated error test failures (check all inputs).",
                                       -5 =>"Repeated convergence failures (perhaps bad Jacobian supplied or wrong choice of MF or tolerances).",
                                       -6 =>"Error weight became zero during problem solution component i vanished, and ATOL or atol[i] = 0.0).")

"""
`dlsode!(F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK, ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, MF, [verbose=true])`

Call DLSODE from the FORTRAN 'odepack' library. Updates the variables `Y`, `T`, `ISTATE`
"""
function dlsode!(F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK,ISTATE, IOPT,
                 RWORK, LRW, IWORK, LIW, JAC, MF;verbose=true)

    ccall((:dlsode_, "libodepack"),Cvoid,
          (Ptr{Cvoid},Ref{Int64},Ptr{Array{Float64,1}},Ref{Float64},Ref{Float64},Ref{Int64},
           Ref{Float64},Ptr{Array{Float64,1}},Ref{Int64},Ref{Int64},Ref{Int64},Ref{Float64},
           Ref{Int64},Ref{Int64},Ref{Int64},Ptr{Cvoid},Ref{Int64}),
          F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK, ISTATE, IOPT,
          RWORK, LRW, IWORK, LIW, JAC, MF)
    verbose ? println("$(istate_dict[ISTATE.x])") : nothing
end

"""
`ode(c_ode_system::Ptr{Cvoid}, [c_jacobian::Ptr{Cvoid}], t_vec::AbstractVector, y0::AbstractVector, atol::AbstractVector=1.0e-3*ones(length(y0)),rtol::Real = 1.0e-3,[verbose=true], [tol_mode = 2])`

ODE solver. The system and the (optional) Jacobian functions should be created with `@diff_eq` and `@diff_eq_jac` macros.
Outputs a vector of solutions at `t_vec` times.
First value of `t_vec` is the initial time.
"""
function ode(c_ode_system::Ptr{Cvoid}, c_jacobian::Ptr{Cvoid},
             t_vec::AbstractVector, y0::AbstractVector,
             atol::AbstractVector=1.0e-3*ones(length(y0)),
             rtol::Real = 1.0e-3;
             verbose=true, tol_mode = 2)

    y = deepcopy(y0)
    t_vec_out = deepcopy(t_vec)
    y_vec = Vector{Vector{Float64}}(undef,length(t_vec))
    n = length(y0)
    @assert length(atol) == n "
        Length of absolute tolerance vector ($(length(atol))) does not match the length of initial conditions vector ($(n))."
    itol = tol_mode
    itask = 1 #Can this be anything else?
    iopt = 0 #TODO: optional arguments
    istate = Ref{Int64}(1)
    mf = 21
    lrw = Int64(22 + 9*n + n^2)
    liw = Int64(20 + n)
    rwork = zeros(Float64,lrw)
    iwork = zeros(Int64,liw)

    t = Ref{Float64}(t_vec[1])

    for (i,tout) in enumerate(t_vec)
        dlsode!(c_ode_system, n, y, t,
                tout, itol, rtol, atol, itask, istate, iopt,
                rwork, lrw, iwork, liw, c_jacobian, mf,verbose=verbose)
        y_vec[i] = deepcopy(y)
    end

    y_vec
end

function ode(c_ode_system::Ptr{Cvoid},
             t_vec::AbstractVector, y0::AbstractVector,
             atol::AbstractVector=1.0e-3*ones(length(y0)),
             rtol::Real=1.0e-3;
             verbose=true, stiff=true, tol_mode=2)

    y = deepcopy(y0)
    y_vec = Vector{Vector{Float64}}(undef,length(t_vec))
    n = length(y0)
    @assert length(atol) == n
    itol = tol_mode
    itask = 1 #Can this be anything else?
    iopt = 0 #TODO: optional arguments
    istate = Ref{Int64}(1)
    if stiff
        mf = 22
        lrw = 22 + 9*n + n^2
        liw = 20 + n
    else
        mf = 10
        lrw = 20 + 16*n
        liw = 20
    end

    rwork = zeros(Float64,lrw)
    iwork = zeros(Int64,liw)

    t = Ref{Float64}(t_vec[1])
    fake_jac = @cfunction((x1,x2,x3,x4,x5,x6,x7)->nothing,Cvoid,(Ptr{Int64},Ptr{Float64},Ptr{Float64},
                                             Ptr{Int64},Ptr{Int64},Ptr{Float64},Ptr{Int64}))

    for (i,tout) in enumerate(t_vec)
        dlsode!(c_ode_system, n, y, t,
                tout, itol, rtol, atol, itask, istate, iopt,
                rwork, lrw, iwork, liw, fake_jac, mf,verbose=verbose)
        y_vec[i] = deepcopy(y)
    end

    y_vec
end

end # module


text:
## Copyright (c) 2015, Júlio Hoffimann Mendes <juliohm@stanford.edu>
##
## Permission to use, copy, modify, and/or distribute this software for any
## purpose with or without fee is hereby granted, provided that the above
## copyright notice and this permission notice appear in all copies.
##
## THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
## WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
## MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
## ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
## WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
## ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
## OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.

"""
    iqsim(training_image::AbstractArray,
          tplsizex::Integer, tplsizey::Integer, tplsizez::Integer,
          gridsizex::Integer, gridsizey::Integer, gridsizez::Integer;
          overlapx::Real=1/6, overlapy::Real=1/6, overlapz::Real=1/6,
          soft::AbstractVector=[], hard::HardData=HardData(), tol::Real=.1,
          cut::Symbol=:boykov, path::Symbol=:rasterup, simplex::Bool=false,
          nreal::Integer=1, threads::Integer=CPU_PHYSICAL_CORES,
          gpu::Bool=false, debug::Bool=false, showprogress::Bool=false)

Performs image quilting simulation as described in Hoffimann et al. 2017.

## Parameters

### Required

* `training_image` can be any 3D array (add ghost dimension for 2D)
* `tplsizex`,`tplsizey`,`tplsizez` is the template size
* `gridsizex`,`gridsizey`,`gridsizez` is the simulation size

### Optional

* `overlapx`,`overlapy`,`overlapz` is the percentage of overlap
* `soft` is a vector of `(data,dataTI)` pairs
* `hard` is an instance of `HardData`
* `tol` is the initial relaxation tolerance in (0,1] (default to .1)
* `cut` is the cut algorithm (`:dijkstra` or `:boykov`)
* `path` is the simulation path (`:rasterup`, `:rasterdown`, `:dilation` or `:random`)
* `simplex` informs whether to apply or not the simplex transform to the image
* `nreal` is the number of realizations
* `threads` is the number of threads for the FFT (default to all CPU cores)
* `gpu` informs whether to use the GPU or the CPU
* `debug` informs whether to export or not the boundary cuts and voxel reuse
* `showprogress` informs whether to show or not estimated time duration

The main output `reals` consists of a list of 3D realizations that can be indexed with
`reals[1], reals[2], ..., reals[nreal]`. If `debug=true`, additional output is generated:

```julia
reals, cuts, voxs = iqsim(..., debug=true)
```

`cuts[i]` is the boundary cut for `reals[i]` and `voxs[i]` is the associated voxel reuse.
"""
function iqsim(training_image::AbstractArray,
               tplsizex::Integer, tplsizey::Integer, tplsizez::Integer,
               gridsizex::Integer, gridsizey::Integer, gridsizez::Integer;
               overlapx::Real=1/6, overlapy::Real=1/6, overlapz::Real=1/6,
               soft::AbstractVector=[], hard::HardData=HardData(), tol::Real=.1,
               cut::Symbol=:boykov, path::Symbol=:rasterup, simplex::Bool=false,
               nreal::Integer=1, threads::Integer=CPU_PHYSICAL_CORES,
               gpu::Bool=false, debug::Bool=false, showprogress::Bool=false)

  # GPU setup
  global GPU = gpu ? gpu_setup() : nothing

  # use all CPU cores in FFT
  FFTW.set_num_threads(threads)

  # sanity checks
  @assert ndims(training_image) == 3 "training image is not 3D (add ghost dimension for 2D)"
  @assert all(0 .< [tplsizex, tplsizey, tplsizez] .≤ [size(training_image)...]) "invalid template size"
  @assert all([gridsizex, gridsizey, gridsizez] .≥ [tplsizex, tplsizey, tplsizez]) "invalid grid size"
  @assert all(0 .< [overlapx, overlapy, overlapz] .< 1) "overlaps must be in range (0,1)"
  @assert 0 < tol ≤ 1 "tolerance must be in range (0,1]"
  @assert cut ∈ [:dijkstra,:boykov] "invalid cut algorithm"
  @assert path ∈ [:rasterup,:rasterdown,:dilation,:random] "invalid simulation path"
  @assert nreal > 0 "invalid number of realizations"

  # soft data checks
  if !isempty(soft)
    for (aux, auxTI) in soft
      @assert ndims(aux) == 3 "soft data is not 3D (add ghost dimension for 2D)"
      @assert all([size(aux)...] .≥ [gridsizex, gridsizey, gridsizez]) "soft data size < grid size"
      @assert size(auxTI) == size(training_image) "auxiliary TI must have the same size as TI"
    end
  end

  # hard data checks
  if !isempty(hard)
    locations = Int[loc[i] for loc in keys(hard), i=1:3]
    @assert all(maximum(locations, 1) .≤ [gridsizex gridsizey gridsizez]) "hard data locations outside of grid"
    @assert all(minimum(locations, 1) .> 0) "hard data locations must be positive indexes"
  end

  # calculate the overlap from given percentage
  overlapx = ceil(Int, overlapx * tplsizex)
  overlapy = ceil(Int, overlapy * tplsizey)
  overlapz = ceil(Int, overlapz * tplsizez)

  # spacing in raster path
  spacingx = tplsizex - overlapx
  spacingy = tplsizey - overlapy
  spacingz = tplsizez - overlapz

  # calculate the number of tiles from grid size
  ntilex = ceil(Int, gridsizex / max(spacingx, 1))
  ntiley = ceil(Int, gridsizey / max(spacingy, 1))
  ntilez = ceil(Int, gridsizez / max(spacingz, 1))

  # simulation grid dimensions
  nx = ntilex * (tplsizex - overlapx) + overlapx
  ny = ntiley * (tplsizey - overlapy) + overlapy
  nz = ntilez * (tplsizez - overlapz) + overlapz

  # total overlap volume in simulation grid
  overlap_volume = nx*ny*nz - (nx - (ntilex-1)overlapx)*
                              (ny - (ntiley-1)overlapy)*
                              (nz - (ntilez-1)overlapz)

  # warn in case of 1-voxel overlaps
  if any(([tplsizex, tplsizey, tplsizez] .>  1) .&
         ([overlapx, overlapy, overlapz] .== 1))
    warn("Overlaps with only 1 voxel. Check template/overlap configuration.")
  end

  # always work with floating point
  TI = Float64.(training_image)

  # inactive voxels in the training image
  NaNTI = isnan.(TI); TI[NaNTI] = 0

  # perform simplex transform
  if simplex
    categories = Set(TI[.!NaNTI])
    ncategories = nvertices = length(categories) - 1

    @assert categories == Set(0:ncategories) "categories should be labeled 1, 2, 3,..."

    simplexTI = simplex_transform(TI, nvertices)
  else
    nvertices = 1
    simplexTI = [TI]
  end

  # disable tiles in the training image if they contain inactive voxels
  mₜ, nₜ, pₜ = size(TI)
  disabled = falses(mₜ-tplsizex+1, nₜ-tplsizey+1, pₜ-tplsizez+1)
  for nanidx in find(NaNTI)
    iₙ, jₙ, kₙ = ind2sub(size(TI), nanidx)

    # tile corners are given by (iₛ,jₛ,kₛ) and (iₑ,jₑ,kₑ)
    iₛ = max(iₙ-tplsizex+1, 1)
    jₛ = max(jₙ-tplsizey+1, 1)
    kₛ = max(kₙ-tplsizez+1, 1)
    iₑ = min(iₙ, mₜ-tplsizex+1)
    jₑ = min(jₙ, nₜ-tplsizey+1)
    kₑ = min(kₙ, pₜ-tplsizez+1)

    disabled[iₛ:iₑ,jₛ:jₑ,kₛ:kₑ] = true
  end

  # keep track of hard data and inactive voxels
  skipped = Set{Tuple{Int,Int,Int}}()
  datum   = Vector{Tuple{Int,Int,Int}}()
  if !isempty(hard)
    # hard data in grid format
    hardgrid = zeros(nx, ny, nz)
    preset = falses(nx, ny, nz)
    activated = trues(nx, ny, nz)
    for loc in keys(hard)
      if isnan(hard[loc])
        activated[loc...] = false
      else
        hardgrid[loc...] = hard[loc]
        preset[loc...] = true
      end
    end

    # deactivate voxels beyond true grid size
    activated[gridsizex+1:nx,:,:] = false
    activated[:,gridsizey+1:ny,:] = false
    activated[:,:,gridsizez+1:nz] = false

    # grid must contain active voxels
    any_activated = any(activated[1:gridsizex,1:gridsizey,1:gridsizez])
    @assert any_activated "simulation grid has no active voxel"

    # determine tiles that should be skipped and tiles with data
    for k=1:ntilez, j=1:ntiley, i=1:ntilex
      # tile corners are given by (iₛ,jₛ,kₛ) and (iₑ,jₑ,kₑ)
      iₛ = (i-1)spacingx + 1
      jₛ = (j-1)spacingy + 1
      kₛ = (k-1)spacingz + 1
      iₑ = iₛ + tplsizex - 1
      jₑ = jₛ + tplsizey - 1
      kₑ = kₛ + tplsizez - 1

      if all(.!activated[iₛ:iₑ,jₛ:jₑ,kₛ:kₑ])
        push!(skipped, (i,j,k))
      else
        if any(preset[iₛ:iₑ,jₛ:jₑ,kₛ:kₑ])
          push!(datum, (i,j,k))
        end
      end
    end
  end

  # preprocess soft data
  softgrid = Vector{Array{Float64,3}}()
  softTI   = Vector{Array{Float64,3}}()
  if !isempty(soft)
    for (aux, auxTI) in soft
      mx, my, mz = size(aux)
      lx = min(mx,nx); ly = min(my,ny); lz = min(mz,nz)

      AUX = padarray(aux, Pad(:symmetric, [0,0,0], [nx-lx,ny-ly,nz-lz]))
      AUX = parent(AUX)
      AUX[isnan.(AUX)] = 0

      AUXTI = copy(auxTI)
      AUXTI[NaNTI] = 0

      # always work with floating point
      AUX   = Float64.(AUX)
      AUXTI = Float64.(AUXTI)

      push!(softgrid, AUX)
      push!(softTI, AUXTI)
    end
  end

  # overwrite path option if data is available
  !isempty(datum) && (path = :datum)

  # select cut algorithm
  boundary_cut = cut == :dijkstra ? dijkstra_cut : boykov_kolmogorov_cut

  # main output is a vector of 3D grids
  realizations = Vector{Array{Float64,3}}()

  # for each realization we have:
  boundarycuts = Vector{Array{Float64,3}}() # boundary cut
  voxelreuse = Vector{Float64}()            # voxel reuse

  # show progress and estimated time duration
  showprogress && (progress = Progress(nreal, color=:black))

  # preallocate memory for distance calculations
  distance = Array{Float64}(mₜ-tplsizex+1, nₜ-tplsizey+1, pₜ-tplsizez+1)

  for real=1:nreal
    # allocate memory for current simulation
    simgrid = zeros(nx, ny, nz)
    debug && (cutgrid = zeros(nx, ny, nz))

    # keep track of pasted tiles
    pasted = Set{Tuple{Int,Int,Int}}()

    # construct simulation path
    simpath = genpath((ntilex,ntiley,ntilez), path, datum)

    # loop simulation grid tile by tile
    for pathidx in simpath
      i, j, k = ind2sub((ntilex,ntiley,ntilez), pathidx)

      # skip tile if all voxels are inactive
      (i,j,k) ∈ skipped && continue

      # if not skipped, proceed and paste tile
      push!(pasted, (i,j,k))

      # tile corners are given by (iₛ,jₛ,kₛ) and (iₑ,jₑ,kₑ)
      iₛ = (i-1)spacingx + 1
      jₛ = (j-1)spacingy + 1
      kₛ = (k-1)spacingz + 1
      iₑ = iₛ + tplsizex - 1
      jₑ = jₛ + tplsizey - 1
      kₑ = kₛ + tplsizez - 1

      # current simulation dataevent
      simdev = view(simgrid, iₛ:iₑ,jₛ:jₑ,kₛ:kₑ)

      # compute overlap distance
      distance[:] = 0
      if overlapx > 1 && (i-1,j,k) ∈ pasted
        ovx = view(simdev,1:overlapx,:,:)
        xsimplex = simplex ? simplex_transform(ovx, nvertices) : [ovx]

        D = convdist(simplexTI, xsimplex)
        distance .+= view(D,1:mₜ-tplsizex+1,:,:)
      end
      if overlapx > 1 && (i+1,j,k) ∈ pasted
        ovx = view(simdev,spacingx+1:tplsizex,:,:)
        xsimplex = simplex ? simplex_transform(ovx, nvertices) : [ovx]

        D = convdist(simplexTI, xsimplex)
        distance .+= view(D,spacingx+1:mₜ-overlapx+1,:,:)
      end
      if overlapy > 1 && (i,j-1,k) ∈ pasted
        ovy = view(simdev,:,1:overlapy,:)
        ysimplex = simplex ? simplex_transform(ovy, nvertices) : [ovy]

        D = convdist(simplexTI, ysimplex)
        distance .+= view(D,:,1:nₜ-tplsizey+1,:)
      end
      if overlapy > 1 && (i,j+1,k) ∈ pasted
        ovy = view(simdev,:,spacingy+1:tplsizey,:)
        ysimplex = simplex ? simplex_transform(ovy, nvertices) : [ovy]

        D = convdist(simplexTI, ysimplex)
        distance .+= view(D,:,spacingy+1:nₜ-overlapy+1,:)
      end
      if overlapz > 1 && (i,j,k-1) ∈ pasted
        ovz = view(simdev,:,:,1:overlapz)
        zsimplex = simplex ? simplex_transform(ovz, nvertices) : [ovz]

        D = convdist(simplexTI, zsimplex)
        distance .+= view(D,:,:,1:pₜ-tplsizez+1)
      end
      if overlapz > 1 && (i,j,k+1) ∈ pasted
        ovz = view(simdev,:,:,spacingz+1:tplsizez)
        zsimplex = simplex ? simplex_transform(ovz, nvertices) : [ovz]

        D = convdist(simplexTI, zsimplex)
        distance .+= view(D,:,:,spacingz+1:pₜ-overlapz+1)
      end

      # disable dataevents that contain inactive voxels
      distance[disabled] = Inf

      # compute hard and soft distances
      auxdistances = Vector{Array{Float64,3}}()
      if !isempty(hard) && any(preset[iₛ:iₑ,jₛ:jₑ,kₛ:kₑ])
        harddev = hardgrid[iₛ:iₑ,jₛ:jₑ,kₛ:kₑ]
        hsimplex = simplex ? simplex_transform(harddev, nvertices) : [harddev]
        D = convdist(simplexTI, hsimplex, weights=preset[iₛ:iₑ,jₛ:jₑ,kₛ:kₑ])

        # disable dataevents that contain inactive voxels
        D[disabled] = Inf

        # swap overlap and hard distances
        push!(auxdistances, distance)
        distance = D
      end
      for n=1:length(softTI)
        softdev = softgrid[n][iₛ:iₑ,jₛ:jₑ,kₛ:kₑ]
        D = convdist([softTI[n]], [softdev])

        # disable dataevents that contain inactive voxels
        D[disabled] = Inf

        push!(auxdistances, D)
      end

      # current pattern database
      patterndb = isempty(auxdistances) ? find(distance .≤ (1+tol)minimum(distance)) :
                                          relaxation(distance, auxdistances, tol)
      patternprobs = tau_model(patterndb, distance, auxdistances)

      # pick a pattern at random from the database
      idx = sample(patterndb, weights(patternprobs))
      iᵦ, jᵦ, kᵦ = ind2sub(size(distance), idx)

      # selected training image dataevent
      TIdev = view(TI,iᵦ:iᵦ+tplsizex-1,jᵦ:jᵦ+tplsizey-1,kᵦ:kᵦ+tplsizez-1)

      # boundary cut mask
      M = falses(simdev)
      if overlapx > 1 && (i-1,j,k) ∈ pasted
        A = view(simdev,1:overlapx,:,:); B = view(TIdev,1:overlapx,:,:)
        M[1:overlapx,:,:] .|= boundary_cut(A, B, :x)
      end
      if overlapx > 1 && (i+1,j,k) ∈ pasted
        A = view(simdev,spacingx+1:tplsizex,:,:); B = view(TIdev,spacingx+1:tplsizex,:,:)
        M[spacingx+1:tplsizex,:,:] .|= flipdim(boundary_cut(flipdim(A, 1), flipdim(B, 1), :x), 1)
      end
      if overlapy > 1 && (i,j-1,k) ∈ pasted
        A = view(simdev,:,1:overlapy,:); B = view(TIdev,:,1:overlapy,:)
        M[:,1:overlapy,:] .|= boundary_cut(A, B, :y)
      end
      if overlapy > 1 && (i,j+1,k) ∈ pasted
        A = view(simdev,:,spacingy+1:tplsizey,:); B = view(TIdev,:,spacingy+1:tplsizey,:)
        M[:,spacingy+1:tplsizey,:] .|= flipdim(boundary_cut(flipdim(A, 2), flipdim(B, 2), :y), 2)
      end
      if overlapz > 1 && (i,j,k-1) ∈ pasted
        A = view(simdev,:,:,1:overlapz); B = view(TIdev,:,:,1:overlapz)
        M[:,:,1:overlapz] .|= boundary_cut(A, B, :z)
      end
      if overlapz > 1 && (i,j,k+1) ∈ pasted
        A = view(simdev,:,:,spacingz+1:tplsizez); B = view(TIdev,:,:,spacingz+1:tplsizez)
        M[:,:,spacingz+1:tplsizez] .|= flipdim(boundary_cut(flipdim(A, 3), flipdim(B, 3), :z), 3)
      end

      # paste quilted pattern from training image
      simdev[.!M] = TIdev[.!M]

      # save boundary cut
      debug && (cutgrid[iₛ:iₑ,jₛ:jₑ,kₛ:kₑ] = M)
    end

    # save voxel reuse
    debug && push!(voxelreuse, sum(cutgrid)/overlap_volume)

    # hard data and shape correction
    if !isempty(hard)
      simgrid[preset] = hardgrid[preset]
      simgrid[.!activated] = NaN
      debug && (cutgrid[.!activated] = NaN)
    end

    # throw away voxels that are outside of the grid
    simgrid = view(simgrid,1:gridsizex,1:gridsizey,1:gridsizez)
    debug && (cutgrid = view(cutgrid,1:gridsizex,1:gridsizey,1:gridsizez))

    # save and continue
    push!(realizations, simgrid)
    debug && push!(boundarycuts, cutgrid)

    # update progress bar
    showprogress && next!(progress)
  end

  debug ? (realizations, boundarycuts, voxelreuse) : realizations
end


text:
# Set up environment for Julia OSX binary distribution
let
    ROOT = abspath(JULIA_HOME,"..")
    ENV["PATH"]="$JULIA_HOME:$(joinpath(ROOT, "libexec")):$(ENV["PATH"])"
    ENV["FONTCONFIG_PATH"] = joinpath(ROOT, "etc", "fonts")
    ENV["GIT_EXEC_PATH"] = joinpath(ROOT, "libexec", "git-core")
    ENV["GIT_TEMPLATE_DIR"] = joinpath(ROOT, "share", "git-core")
    ENV["TK_LIBRARY"] = "/System/Library/Frameworks/Tk.framework/Versions/8.5/Resources/Scripts"
end

max_stars_repo_name: tobig/julia

max_issues_repo_name: tobig/julia

max_forks_repo_name: tobig/julia


text:
@testset "pckmeans" begin
    X = [1 0 -1 0 0; 0 1 0 -1 0]
    C = [0 1 0 0 0; 1 0 1 0 0; 0 1 0 1 0; 0 0 1 0 -1; 0 0 0 -1 0]
    W = abs.(C)
    M = [1 0; -1 0]
    mx, nx = size(X)
    mm, km = size(M)

    pcs = pckmeans(X, C, W, M)
    pc = pcs[end]

    @test isa(pcs, Vector{PartitionalClustering})
    @test pc.X == X
    @test size(pc.C) == (nx, nx)
    @test eltype(pc.C) == Int
    @test size(pc.W) == (nx, nx)
    @test eltype(pc.W) == Float64
    @test size(pc.Y) == (km, nx)
    @test all(y -> 0 <= y <= 1, pc.Y)
    @test size(pc.M) == (mm, km)
end


text:
function run_determination(dir::String; histograms::String="", initial_stepsize::Float64=Float64(pi), K3_range::UnitRange{Int64}=1:8, L::Integer=8, optimizer::String="rotate_hierarchical", initial_temperature_factor::Float64=1.0, temperature_decay::Float64=0.99, N::Integer=32, range=1000:1019, fresh::Bool=false, gpu::Bool=true, Ncores::Integer=8, successive_jobs::Integer=1, measure="Bayes", postprocess::Bool=true, stepsizefactor::Float64=1.02, K2_range::UnitRange{Int64}=1:35, qmax::Float64=1.0, force_repostprocess::Bool=false, run_denoise::Bool=false, architecture::String="", hours::Int64=48, sigma::Float64=0.0, reference_pdb_path::String="", lambda::Float64=0.0, include_negativity::Bool=false, gpu_mem::String="4000M")

    julia_script = """
    using ThreePhotons

    #Let's continue
    if isfile("state.dat") && isfile("params.dat")

        state = deserializeFromFile("state.dat")
        params = deserializeFromFile("params.dat")

        if state["state"] == "running"
            params,state = rotation_search(params, state)
            postprocess_run(params, state, "$(reference_pdb_path)", true, 35, $sigma)
        elseif state["state"] == "finished_structure" && $postprocess
            postprocess_run(params, state, "$(reference_pdb_path)", true, 35, $sigma)
        elseif $force_repostprocess
            #TODO: Temporary fix for wrong results in noise runs
            params["qmax"] = $qmax
            state["intensity"].rmax = $qmax
            postprocess_run(params, state, "$(reference_pdb_path)", true, 35, $sigma)
        end

        #Or start a completely new run
    else
        params,state = rotation_search(Dict( "reference_pdb_path"=>"$(reference_pdb_path)", "stepsizefactor"=>$stepsizefactor, "initial_stepsize" => $initial_stepsize, "L"=>$L, "K3_range"=>$(K3_range), "N"=>$N, "histograms"=>"$(histograms)", "optimizer"=>$optimizer, "initial_temperature_factor"=>$initial_temperature_factor, "measure"=>"$measure", "temperature_decay"=>$temperature_decay, "LMAX"=>25, "K2_range"=>$(K2_range), "qmax"=>$qmax, "lambda"=>$lambda, "include_negativity"=>$include_negativity))
        if $postprocess
            postprocess_run(params, state, "$(reference_pdb_path)", true, 35, $sigma)
        end
    end
    """

    for n in range
        launch_job("$dir/$n", Ncores, gpu, julia_script, successive_jobs; architecture=architecture, hours=hours, fresh=fresh, gpu_mem=gpu_mem)
    end
end

"""Fits all intensites with respect to the first for consecutive averaging"""
function run_postprocess_coliphage_results( dir::String, range::UnitRange{Int64}, L::Int64)
    for n in 1000:1019
        julia_script = """
        using ThreePhotons
        L=$L
        K3=maximum($range)
        intensity_reference = deserializeFromFile("$(environment_path("$dir"))/1000/state.dat")["intensity"]
        intensity = deserializeFromFile("$(environment_path("$dir"))/$n/state.dat")["intensity"]

        bestfit, bestsc, _, _, _ =  fitStructures_random(deleteTerms(intensity, K3, L), deleteTerms(intensity_reference, K3, L), $range, L, 0.995, nworkers()*8)
        serializeToFile("intensity.dat", bestfit)
        saveCube(bestfit, "fitted_intensity.mrc")
        """
        launch_job("$(dir)_fitted/$n", 8, false, julia_script, 1)
    end
end

function run_average_core_completion_phasing(name::String, center_range::UnitRange{Int64}, range::UnitRange{Int64}, coliphage_correlation_path::String)
    for beta_end in [0.8, 0.85, 0.9, 0.95, 0.99]
        julia_script = """
        using ThreePhotons

        _,_,_,_,c1 = loadHistograms(1,1,"$(environment_path("exp_data/$(coliphage_correlation_path)/histo.dat"))")

        extended_corrected_intensity = complete_core("$(environment_path("exp_data/$name"))", c1, $(center_range), $(range))
        averaged_density = phase_completed_intensity(extended_corrected_intensity, 8, $(beta_end))
        """
        launch_job("exp_data/$(name)_phased_center_$(minimum(center_range))_$(maximum(center_range))/b$(beta_end)", 8, false, julia_script, 1)
    end
end

function run_set(image_list::Array{Int64}, K2_range::UnitRange{Int64}=1:38, N::Int64=32, L::Int64=18, K3_range::UnitRange{Int64}=1:26, temperature_decay::Float64=0.99998, ppi::Int64=10, include_infinite::Bool=true)
    histograms_finite = Dict( "P$(img)" => histogram_name("parallel/data_generation/SH_", ppi, N, maximum(K2_range), maximum(K3_range), float(maximum(K2_range)), img, "") for img in image_list)
    histograms_infinite = include_infinite ? Dict("L20_inf" => "expdata/correlations_N$(N)_K$(maximum(K2_range))_L20_inf.dat") : Dict()
    histogram_list = merge(histograms_finite, histograms_infinite)

    for img in keys(histogram_list)
        run_determination("paper_res_vs_pictures_$(ppi)p_KMAX$(maximum(K2_range))_N$(N)_K$(maximum(K3_range))_L$(L)_$(temperature_decay)/$img", histograms=histogram_list[img], initial_stepsize=pi/4.0, K3_range=K3_range, L=L, K2_range=K2_range, qmax=qmax(maximum(K2_range), float(maximum(K2_range))), optimizer="rotate_all_at_once", initial_temperature_factor=0.1, temperature_decay=temperature_decay, N=N, successive_jobs=3, measure="Bayes", range=1000:1019, postprocess=true, gpu=true, Ncores=20, stepsizefactor=1.01, reference_pdb_path="$(ENV["THREEPHOTONS_PATH"])/data/structures/crambin.pdb", lambda=0.0, include_negativity=false)
    end
end

# run_set(calculate_images_ppi(10), 38, 32, 18, 26, 0.99998, 10)
# run_set(calculate_images_ppi(25), 38, 32, 18, 26, 0.99998, 25, false)
# run_set(calculate_images_ppi(50), 38, 32, 18, 26, 0.99998, 50, false)
# run_set(calculate_images_ppi(100), 38, 32, 18, 26, 0.99998, 100, false)

function run_set_vs_L(image_list::Array{Int64}, K2_range::UnitRange{Int64}=1:38, N::Int64=32, LMAX::Int64=18, K3_range::UnitRange{Int64}=1:26, temperature_decay::Float64=0.99998, ppi::Int64=10)
    histogram_list = Dict(
    "P$(img)" => histogram_name("parallel/data_generation/SH_", ppi, N, maximum(K2_range), maximum(K3_range), float(maximum(K2_range)), img, "") for img in images)
    for img in keys(histogram_list)
        for L=2:2:LMAX
            run_determination("paper_res_vs_L_$(ppi)p_KMAX$(maximum(K2_range))_N$(N)_K$(maximum(K3_range))_$(temperature_decay)/$(img)_L$(L)", histograms=histogram_list[img], initial_stepsize=pi/4.0, K3_range=K3_range, L=L, K2_range=K2_range, qmax=qmax(maximum(K2_range), float(maximum(K2_range))), optimizer="rotate_all_at_once", initial_temperature_factor=0.1, temperature_decay=temperature_decay, N=N, successive_jobs=3, measure="Bayes", range=1000:1019, postprocess=true, gpu=true, Ncores=20, stepsizefactor=1.01, reference_pdb_path="$(ENV["THREEPHOTONS_PATH"])/data/structures/crambin.pdb", lambda=0.0)
        end
    end
end

# run_set_vs_L([5120000, 20480000, 81920000, 3276800000], 38, 32, 18, 26, 0.99998, 10)

function run_noise_set(sigmas::Array{Float64}, gammas::Array{Float64}, K2_range::UnitRange{Int64}=1:38, N::Int64=32, L::Int64=18, K3_range::UnitRange{Int64}=1:26, temperature_decay::Float64=0.99998, ppi::Int64=10)
    for sigma in sigmas
        for gamma in gammas
            histo_name = histogram_name("parallel/data_generation/SH_", ppi, N, maximum(K2_range), maximum(K3_range), float(maximum(K2_range)), 3276800000, "", gamma, sigma)
            run_determination("paper_noise_$(ppi)p_KMAX$(maximum(K2_range))_N$(N)_K$(maximum(K3_range))_L$(L)_$(temperature_decay)/G$(gamma)_S$(sigma)", histograms=histo_name, initial_stepsize=pi/4.0, K3_range=K3_range, L=L, K2_range=K2_range, qmax=qmax(maximum(K2_range), float(maximum(K2_range))), optimizer="rotate_all_at_once", initial_temperature_factor=0.1, temperature_decay=temperature_decay, N=N, successive_jobs=1, measure="Bayes", range=1000:1019, postprocess=true, gpu=true, Ncores=8, stepsizefactor=1.01, run_denoise=true, sigma=sigma, force_repostprocess=true, reference_pdb_path="$(ENV["THREEPHOTONS_PATH"])/data/structures/crambin.pdb")
        end
    end
end

  # run_noise_set([0.5, 0.75, 1.125, 2.5], collect(0.1:0.1:0.5), 38, 32, 18, 26, 0.99998, 10)

#For multi particle
# run_determination("multi_particle", histograms="$(ENV["DETERMINATION_DATA"])/output_owl/data_generation/multi_2_SH_10p_N32_K38_R38.0_P3276800000/histo.dat", lambda=0.0, initial_stepsize=pi/4.0, K3_range=1:26, L=18, K2_range=1:38, qmax=qmax(38, float(38)), optimizer="rotate_all_at_once", initial_temperature_factor=0.1, temperature_decay=0.99998, N=32, successive_jobs=3, measure="Bayes", range=1000:1019, postprocess=true, gpu=true, Ncores=20, stepsizefactor=1.01, reference_pdb_path="$(ENV["DETERMINATION_DATA"])/structures/crambin.pdb", force_repostprocess=true, run_denoise=true)

#Repeat original calculations for performance comparison
# run_determination("performance_CUDA", histograms="$(ENV["DETERMINATION_DATA"])/output_owl/data_generation/SH_10p_N32_K38_R38.0_P3276800000/histo.dat", lambda=0.0, initial_stepsize=pi/4.0, K3_range=1:26, L=18, K2_range=1:38, qmax=qmax(38, float(38)), optimizer="rotate_all_at_once", initial_temperature_factor=0.1, temperature_decay=0.99998, N=32, successive_jobs=1, measure="Bayes", range=1000:1000, postprocess=true, gpu=true, Ncores=20, stepsizefactor=1.01, reference_pdb_path="$(ENV["DETERMINATION_DATA"])/structures/crambin.pdb", force_repostprocess=true, run_denoise=true)

#For Ewald SH structure determination
# run_determination("Ewald_lambda_2.0_SH_10p_N32_K2_38_K3_26_R38.0_P3276800000", histograms="$(ENV["DETERMINATION_DATA"])/output_owl/data_generation/Ewald_lambda_2.0_SH_10p_N32_K2_38_K3_26_R38.0_P3276800000/histo.dat", lambda=2.0, initial_stepsize=pi/4.0, K3_range=1:26, L=18, K2_range=1:38, qmax=qmax(38, float(38)), optimizer="rotate_all_at_once", initial_temperature_factor=0.1, temperature_decay=0.99998, N=32, successive_jobs=3, measure="Bayes", range=1000:1019, postprocess=true, gpu=true, Ncores=20, stepsizefactor=1.01, reference_pdb_path="$(ENV["DETERMINATION_DATA"])/structures/crambin.pdb")

# run_determination("Ewald_lambda_4.0_alt_infinity", histograms="$(ENV["DETERMINATION_DATA"])/output/data_generation/correlations_N32_K238_K326_L18_inf.dat", lambda=4.0, initial_stepsize=pi/4.0, K3_range=1:26, L=18, K2_range=1:38, qmax=qmax(38, float(38)), optimizer="rotate_all_at_once", initial_temperature_factor=0.1, temperature_decay=0.99998, N=32, successive_jobs=3, measure="Bayes", range=1000:1019, postprocess=true, gpu=true, Ncores=20, stepsizefactor=1.01, reference_pdb_path="$(ENV["DETERMINATION_DATA"])/structures/crambin.pdb")

#Coliphage with correct qmax and lmabda, lambda=7.75 A, max_resolution=116 A, qmax= pi / 116
# run_determination("exp_data/coliphage_determination_newhisto", histograms="$(ENV["DETERMINATION_DATA"])/output_owl/exp_data/coliphage/histo.dat", lambda=0.0, initial_stepsize=pi/4.0, K3_range=5:26, L=12, K2_range=5:38, qmax=pi/116.0, optimizer="rotate_all_at_once", initial_temperature_factor=0.1, temperature_decay=0.99998, N=32, successive_jobs=3, measure="Bayes", range=1000:1019, postprocess=false, gpu=true, Ncores=20, stepsizefactor=1.01, include_negativity=false)

#Coliphage with symmetric and L=12
# run_determination("exp_data/coliphage_determination_newhisto_symmetric", histograms="$(ENV["DETERMINATION_DATA"])/output_owl/exp_data/coliphage_symmetric/histo.dat", lambda=0.0, initial_stepsize=pi/3.0, K3_range=6:26, L=12, K2_range=6:38, qmax=pi/116.0, optimizer="rotate_all_at_once", initial_temperature_factor=0.1, temperature_decay=0.99998, N=32, successive_jobs=3, measure="Bayes", range=1000:1019, postprocess=false, gpu=true, Ncores=20, stepsizefactor=1.01, include_negativity=false)

#Coliphage with symmetric and L=12, K3=32
# run_determination("exp_data/coliphage_determination_newhisto_symmetric_K32", histograms="$(ENV["DETERMINATION_DATA"])/output_owl/exp_data/coliphage_symmetric/histo.dat", lambda=0.0, initial_stepsize=pi/3.0, K3_range=6:32, L=12, K2_range=6:38, qmax=pi/116.0, optimizer="rotate_all_at_once", initial_temperature_factor=0.1, temperature_decay=0.99998, N=32, successive_jobs=3, measure="Bayes", range=1000:1019, postprocess=false, gpu=true, Ncores=20, stepsizefactor=1.01, include_negativity=false)

#Coliphage with symmetric and L=12 and low range
# run_determination("exp_data/coliphage_determination_paper_lowrange", histograms="$(ENV["DETERMINATION_DATA"])/output_owl/exp_data/coliphage_K2_38_K3_30_N32/histo.dat", lambda=0.0, initial_stepsize=pi/4.0, K3_range=8:26, L=12, K2_range=8:38, qmax=pi/90.0, optimizer="rotate_all_at_once", initial_temperature_factor=0.1, temperature_decay=0.99998, N=32, successive_jobs=3, measure="Bayes", range=1000:1019, postprocess=false, gpu=true, Ncores=20, stepsizefactor=1.01, include_negativity=false)

#Coliphage with symmetric and L=12 and low range
# run_determination("exp_data/coliphage_determination_paper_lowrange_2", histograms="$(ENV["DETERMINATION_DATA"])/output_owl/exp_data/coliphage_K2_38_K3_30_N32/histo.dat", lambda=0.0, initial_stepsize=pi/4.0, K3_range=7:26, L=12, K2_range=3:38, qmax=pi/90.0, optimizer="rotate_all_at_once", initial_temperature_factor=0.25, temperature_decay=0.99998, N=32, successive_jobs=3, measure="Bayes", range=1000:1019, postprocess=false, gpu=true, Ncores=20, stepsizefactor=1.01, include_negativity=false)

#Coliphage with symmetric and L=12 and high range
# run_determination("exp_data/coliphage_determination_paper_K2_5_38_K3_5_30", histograms="$(ENV["DETERMINATION_DATA"])/output_owl/exp_data/coliphage_K2_38_K3_30_N32/histo.dat", lambda=0.0, initial_stepsize=pi/4.0, K3_range=5:30, L=12, K2_range=5:38, qmax=pi/90.0, optimizer="rotate_all_at_once", initial_temperature_factor=0.2, temperature_decay=0.99998, N=32, successive_jobs=3, measure="Bayes", range=1000:1019, postprocess=false, gpu=true, Ncores=20, stepsizefactor=1.01, include_negativity=false)

#Coliphage with symmetric and L=12 and high range
# run_determination("exp_data/coliphage_determination_paper_K2_7_35_K3_7_30", histograms="$(ENV["DETERMINATION_DATA"])/output_owl/exp_data/coliphage_K2_38_K3_30_N32/histo.dat", lambda=0.0, initial_stepsize=pi/4.0, K3_range=7:30, L=12, K2_range=7:35, qmax=pi/90.0, optimizer="rotate_all_at_once", initial_temperature_factor=0.2, temperature_decay=0.99998, N=32, successive_jobs=3, measure="Bayes", range=1000:1019, postprocess=false, gpu=true, Ncores=20, stepsizefactor=1.01, include_negativity=false)

#Coliphage with symmetric and L=12 and high range
# run_determination("exp_data/coliphage_determination_paper_K2_10_35_K3_10_30", histograms="$(ENV["DETERMINATION_DATA"])/output_owl/exp_data/coliphage_K2_38_K3_30_N32/histo.dat", lambda=0.0, initial_stepsize=pi/4.0, K3_range=10:30, L=12, K2_range=10:35, qmax=pi/90.0, optimizer="rotate_all_at_once", initial_temperature_factor=0.1, temperature_decay=0.99998, N=32, successive_jobs=3, measure="Bayes", range=1000:1019, postprocess=false, gpu=true, Ncores=20, stepsizefactor=1.01, include_negativity=false)

# run_determination("exp_data/coliphage_determination_paper_K2_10_36_K3_10_26", histograms="$(ENV["DETERMINATION_DATA"])/output_owl/exp_data/coliphage_K2_38_K3_30_N32/histo.dat", lambda=0.0, initial_stepsize=pi/4.0, K3_range=10:26, L=12, K2_range=10:36, qmax=pi/90.0, optimizer="rotate_all_at_once", initial_temperature_factor=0.1, temperature_decay=0.99998, N=32, successive_jobs=3, measure="Bayes", range=1000:1019, postprocess=false, gpu=true, Ncores=20, stepsizefactor=1.01, include_negativity=false)

#Coliphage with L=12, K=12
# run_determination("exp_data/coliphage_determination_paper_K2_12_38_K3_12_26_L12", histograms="$(ENV["DETERMINATION_DATA"])/output_owl/exp_data/coliphage_K2_38_K3_30_N32/histo.dat", lambda=0.0, initial_stepsize=pi/4.0, K3_range=12:26, L=12, K2_range=12:38, qmax=pi/90.0, optimizer="rotate_all_at_once", initial_temperature_factor=0.1, temperature_decay=0.99998, N=32, successive_jobs=3, measure="Bayes", range=1000:1019, postprocess=false, gpu=true, Ncores=20, stepsizefactor=1.01, include_negativity=false)

#Coliphage smaller dq for L=14 and L=12
# run_determination("exp_data/coliphage_determination_paper2_K2_12_40_K3_12_30_L14", histograms="$(ENV["DETERMINATION_DATA"])/output_owl/exp_data/coliphage_K2_50_K3_30_N32/histo.dat", lambda=0.0, initial_stepsize=pi/4.0, K3_range=12:30, L=14, K2_range=12:40, qmax=pi/90.0, optimizer="rotate_all_at_once", initial_temperature_factor=0.1, temperature_decay=0.99998, N=32, successive_jobs=3, measure="Bayes", range=1000:1019, postprocess=false, gpu=true, Ncores=8, stepsizefactor=1.01, include_negativity=false, gpu_mem="4000M")

#Coliphage best bet
# run_determination("exp_data/coliphage_determination_rerun_K2_5_38_K3_5_28_L12", histograms="$(ENV["DETERMINATION_DATA"])/output_owl/exp_data/coliphage_K2_40_K3_30_N32_processed/histo.dat", lambda=0.0, initial_stepsize=pi/4.0, K3_range=5:28, L=12, K2_range=5:38, qmax=pi/90.0, optimizer="rotate_all_at_once", initial_temperature_factor=0.1, temperature_decay=0.99998, N=32, successive_jobs=3, measure="Bayes", range=1000:1019, postprocess=false, gpu=true, Ncores=8, stepsizefactor=1.01, include_negativity=false, gpu_mem="4000M")

#Rerun Ewald stuff
# run_determination("Ewald_P327680000_lambda_2.5", histograms="$(ENV["DETERMINATION_DATA"])/output_owl/data_generation/Ewald_large_SH_10p_N32_K2_38_K3_26_R38.0_P327680000_lambda_2.5/histo.dat", lambda=2.5, initial_stepsize=pi/4.0, K3_range=1:26, L=18, K2_range=1:38, qmax=qmax(38, float(38)), optimizer="rotate_all_at_once", initial_temperature_factor=0.1, temperature_decay=0.99998, N=32, successive_jobs=3, measure="Bayes", range=1000:1019, postprocess=true, gpu=true, Ncores=8, stepsizefactor=1.01, reference_pdb_path="$(ENV["DETERMINATION_DATA"])/structures/crambin.pdb", include_negativity=false, gpu_mem="4000M")


text:
export
    fovray,
    fovtrg,
    frame,
    frinfo,
    frmnam,
    ftncls,
    furnsh

"""
    fovray(inst, raydir, rframe, abcorr, observer, et)

Determine if a specified ray is within the field-of-view (FOV) of a specified instrument
at a given time.

### Arguments ###

- `inst`: Name or ID code string of the instrument
- `raydir`: Ray's direction vector
- `rframe`: Body-fixed, body-centered frame for target body
- `abcorr`: Aberration correction flag
- `observer`: Name or ID code string of the observer
- `et`: Time of the observation (seconds past J2000)

### Output ###

Returns `true` if the ray is visible.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/fovray_c.html)
"""
function fovray(inst, raydir, rframe, abcorr, observer, et)
    @checkdims 3 raydir
    visible = Ref{SpiceBoolean}()
    ccall((:fovray_c, libcspice), Cvoid,
          (Cstring, Ref{SpiceDouble}, Cstring, Cstring, Cstring, Ref{SpiceDouble}, Ref{SpiceBoolean}),
          inst, raydir, rframe, abcorr, observer, et, visible)
    handleerror()
    Bool(visible[])
end

"""
    fovtrg(inst, target, tshape, tframe, abcorr, obsrvr, et)

Determine if a specified ephemeris object is within the field-of-view (FOV) of a specified
instrument at a given time.

### Arguments ###

- `inst`: Name or ID code string of the instrument.
- `target`: Name or ID code string of the target.
- `tshape`: Type of shape model used for the target.
- `tframe`: Body-fixed, body-centered frame for target body.
- `abcorr`: Aberration correction flag.
- `obsrvr`: Name or ID code string of the observer.
- `et`: Time of the observation (seconds past J2000).

### Output ###

Returns `true` if the object is visible.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/fovtrg_c.html)
"""
function fovtrg(inst, target, tshape, tframe, abcorr, obsrvr, et)
    visible = Ref{SpiceBoolean}()
    ccall((:fovtrg_c, libcspice), Cvoid,
          (Cstring, Cstring, Cstring, Cstring, Cstring, Cstring, Ref{SpiceDouble}, Ref{SpiceBoolean}),
          inst, target, tshape, tframe, abcorr, obsrvr, et, visible)
    handleerror()
    Bool(visible[])
end

"""
    frame(x)

Given a vector `x`, this routine builds a right handed orthonormal frame `x`, `y`, `z`
where the output `x` is parallel to the input `x`.

### Arguments ###

- `x`: Input vector

### Output ###

- `x`: Unit vector parallel to `x` on output
- `y`: Unit vector in the plane orthogonal to `x`
- `z`: Unit vector given by `x × y`

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/frame_c.html)
"""
function frame(x)
    @checkdims 3 x
    x = copy(x)
    y = Array{SpiceDouble}(undef, 3)
    z = Array{SpiceDouble}(undef, 3)
    ccall((:frame_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}),
          x, y, z)
    x, y, z
end

"""
    frinfo(frcode)

Retrieve the minimal attributes associated with a frame needed for converting transformations
to and from it.

### Arguments ###

- `frcode`: The id code for a reference frame

### Output ###

- `cent`: The center of the frame
- `frclss`: The class (type) of the frame
- `clssid`: The idcode for the frame within its class

Returns `nothing` if no frame with id `frcode` could be found.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/frinfo_c.html)
"""
function frinfo(frcode)
    cent = Ref{SpiceInt}()
    frclss = Ref{SpiceInt}()
    clssid = Ref{SpiceInt}()
    found = Ref{SpiceBoolean}()
    ccall((:frinfo_c, libcspice), Cvoid,
          (SpiceInt, Ref{SpiceInt}, Ref{SpiceInt}, Ref{SpiceInt}, Ref{SpiceBoolean}),
          frcode, cent, frclss, clssid, found)

    Bool(found[]) || return nothing
    Int(cent[]), Int(frclss[]), Int(clssid[])
end

"""
    frmnam(frcode)

Retrieve the name of a reference frame associated with an id code.

### Arguments ###

- `frcode`: The id code for a reference frame

### Output ###

Returns the name associated with the reference frame.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/frmnam_c.html)
"""
function frmnam(frcode)
    lenout = 33
    frname = Array{SpiceChar}(undef, lenout)
    ccall((:frmnam_c, libcspice), Cvoid, (SpiceInt, SpiceInt, Ref{SpiceChar}),
          frcode, lenout, frname)
    chararray_to_string(frname)
end

@deprecate ftncls close

"""
    furnsh(kernels...)

Load one or more SPICE kernels into a program.

### Arguments ###

- `kernels`: Path(s) of SPICE kernels to load

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/furnsh_c.html)
"""
function furnsh(kernels...)
    for kernel in kernels
        ccall((:furnsh_c, libcspice), Cvoid, (Cstring,), kernel)
        handleerror()
    end
end


text:
## Aircraft stability analysis example
using AeroMDAO

## Lifting surfaces

# Wing
wing  = WingSection(span       = 8.0,
                    dihedral   = 5.0,
                    sweep        =  15.0,
                    taper      = 0.4,
                    root_chord = 2.0,
                    root_twist = 0.0,
                    tip_twist  = -2.0,
                    root_foil  = naca4(2,4,1,2),
                    tip_foil   = naca4(2,4,1,2),
                    position   = [0., 0., 0.])

# Horizontal tail
htail = WingSection(span       = 2.0,
                    dihedral   = 0.0,
                    sweep        =  15.0,
                    taper      = 0.6,
                    root_chord = 0.8,
                    root_twist = 0.0,
                    tip_twist  = 0.0,
                    root_foil  = naca4(0,0,1,2),
                    tip_foil   = naca4(0,0,0,9),
                    position   = [5., 0., -0.1],
                    angle      = 0.,
                    axis       = [0., 1., 0.]);

# Vertical tail
vtail = HalfWingSection(span       = 0.8,
                        dihedral   = 0.0,
                        sweep        =  8.0,
                        taper      = 0.6,
                        root_chord = 0.8,
                        root_twist = 0.0,
                        tip_twist  = 0.,
                        root_foil  = naca4(0,0,0,9),
                        tip_foil   = naca4(0,0,0,9),
                        position   = [5., 0., 0.],
                        angle      = 90.,
                        axis       = [1., 0., 0.])

# Print info
print_info(wing, "Wing")
print_info(htail, "Horizontal Tail")
print_info(vtail, "Vertical Tail")

## Static stability
x_w, y_w, z_w = wing_mac  = mean_aerodynamic_center(wing)
x_h, y_h, z_h = htail_mac = mean_aerodynamic_center(htail)
x_v, y_v, z_v = vtail_mac = mean_aerodynamic_center(vtail)

S, b, c = projected_area(wing), span(wing), mean_aerodynamic_chord(wing)

# Horizontal tail volume coefficient
S_h = projected_area(htail)
l_h = x_h - x_w
V_h = S_h * l_h / (S * c)
println("Horizontal TVC      V_h: $V_h")

# Vertical tail volume coefficient
S_v = projected_area(vtail)
l_v = x_v - x_w
V_v = S_v * l_v / (S * b)
println("Vertical TVC        V_v: $V_v")

## Meshing and assembly
wing_mesh  = WingMesh(wing, [20], 10)
htail_mesh = WingMesh(htail, [10], 5)
vtail_mesh = WingMesh(vtail, [10], 5)

aircraft = ComponentArray(
                          wing  = make_horseshoes(wing_mesh),
                          htail = make_horseshoes(htail_mesh),
                          vtail = make_horseshoes(vtail_mesh)
                         );

## Evaluate case
fs      = Freestream(alpha = 1.0,
                     beta  = 0.0,
                     omega = [0.,0.,0.])
refs    = References(
                     speed    = 10.0,
                     density  = 1.225,
                     area     = projected_area(wing),
                     span     = span(wing),
                     chord    = mean_aerodynamic_chord(wing),
                     location = mean_aerodynamic_center(wing)
                    )

ac_name = :aircraft
@time dv_data = solve_case_derivatives(
                       aircraft, fs, refs;
                       axes             = Wind(),
                       name             = ac_name,
                       print            = true,    # Prints the results for only the aircraft
                       print_components = true,    # Prints the results for all components
                      );

## Aerodynamic quantities of aircraft
ac_dvs = dv_data[ac_name]

# Longitudinal stability

## Center of pressure
x_cp = -refs.chord * ac_dvs.Cm / ac_dvs.CZ

## Neutral point
x_np = -refs.chord * ac_dvs.Cm_alpha / ac_dvs.CZ_alpha

## Spiral stability
γ = ac_dvs.Cl_beta * ac_dvs.Cn_rbar / (ac_dvs.Cl_rbar * ac_dvs.Cn_beta)

## Other quantities
htail_dvs = dv_data.htail

# Locations
x_np, y_np, z_np = loc_np = [ refs.location[1] .+ x_np; zeros(2) ]  # Neutral point
x_cp, y_cp, z_cp = loc_cp = [ refs.location[1] .+ x_cp; zeros(2) ]  # Center of pressure

println("Aerodynamic Center x_ac: $(x_w) m")
println("Neutral Point      x_np: $(x_np) m")
println("Center of Pressure x_cp: $(x_cp) m")
println("Spiral Stability      γ: $γ")

## Plotting everything
using Plots

gr(dpi = 300)

wing_plan  = plot_planform(wing)
htail_plan = plot_planform(htail)
vtail_plan = plot_planform(vtail)

#
z_limit = b
aircraft_plot = plot(
    xaxis = "x", yaxis = "y", zaxis = "z",
    camera = (30, 60),
    xlim = (0, z_limit),
    # ylim = (-z_limit/2, z_limit/2),
    zlim = (-z_limit/2, z_limit/2),
    size = (900, 600),
    # legend = false
    )

plot!(wing_plan[:,1], wing_plan[:,2], wing_plan[:,3],  label = "Wing")
plot!(htail_plan[:,1], htail_plan[:,2], htail_plan[:,3], label = "Horizontal Tail")
plot!(vtail_plan[:,1], vtail_plan[:,2], vtail_plan[:,3], label = "Vertical Tail")

scatter!(Tuple(wing_mac),  color = :blue,  label = "Wing MAC")
scatter!(Tuple(htail_mac), color = :red,   label = "Horizontal Tail MAC")
scatter!(Tuple(vtail_mac), color = :green, label = "Vertical Tail MAC")

scatter!(Tuple(loc_np), color = :orange, label = "Neutral Point")
scatter!(Tuple(loc_cp), color = :brown, label = "Center of Pressure")
savefig(aircraft_plot, "AircraftData.png")
plot!()

## Alpha sweep
function alpha_sweep(aircraft, refs, α)
    ## Evaluate case
    fs      = Freestream(alpha = α,
                         beta  = 0.0,
                         omega = [0.,0.,0.])

    dv_data = solve_case_derivatives(aircraft, fs, refs;
                                    #  print_components = true
                                    )
end

##
αs  = -8:0.5:8
@time res = map(α -> alpha_sweep(aircraft, refs, α), αs);

## Data processing

## Lift coefficient
CLs = [ re.aircraft.CZ for re in res ]

plot(αs, CLs, xlabel = "α", ylabel = "CL")

## Moment coefficient
Cms = [ re.aircraft.Cm for re in res ]

plot(αs, Cms, xlabel = "α", ylabel = "Cm")

# Cm-CL (only relevant for non-linear CL-α variation)
# plot(CLs, Cms, xlabel = "CL", ylabel = "Cm")

##
Cm_αs  = [ re.aircraft.Cm_alpha for re in res ]
CL_αs  = [ re.aircraft.CZ_alpha for re in res ]

##
plot(CL_αs, Cm_αs, xlabel = "CL_α", ylabel = "Cm_α")

##
∂Cm_∂CLs = Cm_αs ./ CL_αs

plot(αs, ∂Cm_∂CLs, ls = :solid, xlabel = "α", ylabel = "Cm/CL")


text:
##TODO: Add softmax

struct RepQ{T <: SDF} <: SDF
    Victim::T
    resid::Float64
end

function to_code(Sdf::RepQ)
    SF = to_code(Sdf.Victim)
    entry = gensym()
    mapped = gensym()
    prepend!(SF.Core,[Expr(:(=),SF.Param,Expr(:call,:broadcast,
        Expr(:->,mapped,Expr(:block,
            Expr(:call,:mod,mapped,Sdf.resid))),entry))])
    SF.Param=entry
    SF
end

struct Trans{T <: SDF} <: SDF
    Victim::T
    Vec::space
end

function to_code(Sdf::Trans)
    SF = to_code(Sdf.Victim)
    temp = gensym()
    prepend!(SF.Core,[Expr(:(=),SF.Param,Expr(:call,:(-),temp,Sdf.Vec))])
    SF.Param=temp
    SF
end

struct SUnion{T} <: SDF
    Victims::T
end

function to_code(Sdf::SUnion)
    entry = gensym()
    best = gensym()

    victs = [Expr(:(=),best,(Inf64,Shd_Id(0)))]
    for i in Sdf.Victims
        SF = to_code(i)
        prepend!(SF.Core,[Expr(:(=),SF.Param,entry)])
        append!(victs,[Expr(:(=),best,Expr(:call,:min,best,Expr(:block,SF.Core...)))])
    end
    append!(victs,[Expr(:block,best)])
    SingFunc(victs,entry)
end

struct SCut{T} <: SDF
    Victims::T
end

function to_code(Sdf::SCut)
    entry = gensym()
    best = gensym()

    victs = [Expr(:(=),best,(0.0,Shd_Id(0)))]
    for i in Sdf.Victims
        SF = to_code(i)
        prepend!(SF.Core,[Expr(:(=),SF.Param,entry)])
        append!(victs,[Expr(:(=),best,Expr(:call,:max,best,Expr(:block,SF.Core...)))])
    end
    append!(victs,[Expr(:block,best)])
    SingFunc(victs,entry)
end


text:
module GenomicMaps

using Luxor
using Colors
using GenomicAnnotations
using GenomicAnnotations: AbstractGene, Gene
using Statistics
using Printf: @sprintf

include("initialise.jl")
include("shapes.jl")
include("drawstretch.jl")
export drawgenome
include("annotatesvg.jl")
include("drawplasmid.jl")
export drawplasmid
include("genomecomparison.jl")
export drawcomparison

end #module


text:
using DSP: conv

include("transformed.jl")

@doc raw"""
This struct represents a planned continuous wavelet transform. It contains the voices on which the
transform is performed as well the pre-computed analysis wavelets.

# Fields
- `wavelet`: Holds a reference to the wavelet to use for the transformation.
- `blocks`: The grouping of the voices of same kernel size for faster block-convolution.
- `scales`: The scales on which the wavelet transform will be performed.
- `kernels`: Matrix of evaluated kernels for the convolution.
"""
struct ContinuousWaveletTransform
    wavelet::GenericContinuousWavelet
    blocks::Array{Tuple{Int32,Int32,Int32}, 1}
    scales::Array{Float64, 1}
    kernels::Array{Complex{Float64}, 2}

    @doc raw"""
        ContinuousWaveletTransform(wav::GenericContinuousWavelet, scales::AbstractVector)

    Constructs a new wavelet transform for the given wavelet (`wav`) and `scales`. This constructor
    pre-computes the kernels for the given scales. It is therefore efficient to perform several
    transformations using `transform()` with the same `ContinuousWaveletTransform` object.
    """
    function ContinuousWaveletTransform(wav::GenericContinuousWavelet, scales::AbstractVector)
        blocks  = Tuple{Int32,Int32,Int32}[]
        lengths = Int32[]
        sscales = sort(scales)
        sizes   = Array{Int32}(exp2.(ceil.(log2.(2*cutoff_time(wav).*sscales))))
        kernels = zeros(Complex{Float64}, sizes[end], length(scales))

        # Start the first block for the first kernel
        push!(blocks, (1, 1, sizes[1]) )
        kernels[1:sizes[1],1] = map(t -> eval_analysis(wav, (t-sizes[1]/2)/sscales[1])/sscales[1], 1:sizes[1])
        # process the remaining kernels
        for i in 2:length(scales)
            kernels[1:sizes[i],i] = map(t -> eval_analysis(wav, (t-sizes[i]/2)/sscales[i])/sscales[i], 1:sizes[i])
            # if same size than last block:
            if blocks[end][3]==sizes[i]
                # -> add to block
                blocks[end] = (blocks[end][1], i, sizes[i]);
            else
                # otherwise add new block
                push!(blocks, (i,i,sizes[i]))
            end
        end
        nBlocks = length(blocks)
        minBlockSize = minimum(sizes)
        maxBlockSize = maximum(sizes)
        new(wav, blocks, sscales, kernels)
    end
end

@doc raw"""
    transform(trans::ContinuousWaveletTransform, x::AbstractArray{Float64,1})

Performes the planned continuous wavelet transfrom given by `trans` of the timeseries `x`. The
result is a `ContinuousWavelet.ContinuousWaveletTransformed` object.
"""
function transform(trans::ContinuousWaveletTransform, x::AbstractArray{Float64,1})
    N = length(x)
    M = length(trans.scales)
    wX = Array{Complex{Float64},2}(undef, N,M)
    for (a,b,K) in trans.blocks
        off = K÷2
        wX[:,(a:b)] = conv(trans.kernels[(1:K),(a:b)], x)[(off):(off+N-1),:]
    end

    ContinuousWaveletTransformed(wX, trans.scales, trans.wavelet)
end


text:
# XXX: rename use_debug_checks
use_debug_checks() = false

macro debug_check(ex)
  quote
    if use_debug_checks()
      $(esc(ex))
    end
  end
end

function enable_debug_checks()
  if !getfield(@__MODULE__, :use_debug_checks)()
    Core.eval(@__MODULE__, :(use_debug_checks() = true))
  end
end

function disable_debug_checks()
  if getfield(@__MODULE__, :use_debug_checks)()
    Core.eval(@__MODULE__, :(use_debug_checks() = false))
  end
end


text:
using jInv.Utils
using Test

v  = rand(1:200,54)
vu = unique(v)
w  = randn(100)

data = (vu,w)
# testing sortpermFast

for k=1:length(data)

	res1 = sortpermFast(data[k])
	res2 = sortperm(data[k])
	b2   = data[k][res2]

	@test all(res1[1] .== res2)
	@test all(res1[2] .== b2)
end

# Test two argument version
vin  = copy(v)
vuin = copy(vu)
d    = rand(1:1000,54)
du   = d[1:length(vu)]
din  = copy(d)
duin = copy(du)

ii,vs  = sortpermFast(v)
jj,vus = sortpermFast(vu)

vin,din   = sortpermFast(vin,din)
vuin,duin = sortpermFast(vuin,duin)

@test vin == vuin == vus
@test duin == du[jj]
@test length(din) == length(duin)


text:
function __init__()

    Requires.@require Cairo="159f3aea-2a34-519c-b102-8c37f9878175" begin
        function addtocontext!(cr::Cairo.CairoContext,
                               Pl::Plot2D,
                               P::Path2D,
                               bb::BoundingBox)
            Cairo.set_line_width(cr,lwcorrect(P.pen.linewidth))
            Cairo.move_to(cr,P.points[1].x,P.points[1].y);
            for i=2:length(P.points)-1
                Cairo.line_to(cr,P.points[i].x,P.points[i].y);
            end
            if P.arrow.name ≠ "None"
                arrowlength = 0.005*P.arrow.size
                p = arrowlength/abs(P.points[end-1]-P.points[end])
                lastpoint = p*(P.points[end-1]) + (1-p)*P.points[end]
                Cairo.line_to(cr,lastpoint.x,lastpoint.y)
            else
                Cairo.line_to(cr,P.points[end].x,P.points[end].y)
            end
            Cairo.set_source_rgba(cr,P.pen.color...,P.pen.opacity)
            Cairo.stroke(cr)
            if P.arrow.name ≠ "None"
                addtocontext!(cr,Pl,arrowhead(P,bb),bb)
            end
        end

        function addtocontext!(cr::Cairo.CairoContext,
                               Pl::Plot2D,
                               P::Polygon2D,
                               bb::BoundingBox)
            if is_no_pen(P.fillpen) && is_no_pen(P.pen) && !P.clip
                return nothing
            end
            if P.clip
                Cairo.move_to(cr,P.points[1].x,P.points[1].y);
                for i=2:length(P.points)
                    Cairo.line_to(cr,P.points[i].x,P.points[i].y);
                end
                Cairo.move_to(cr,bb.xmin,bb.ymin)
                if counterclockwise(P)
                    Cairo.line_to(cr,bb.xmin,bb.ymax)
                    Cairo.line_to(cr,bb.xmax,bb.ymax)
                    Cairo.line_to(cr,bb.xmax,bb.ymin)
                else
                    Cairo.line_to(cr,bb.xmax,bb.ymin)
                    Cairo.line_to(cr,bb.xmax,bb.ymax)
                    Cairo.line_to(cr,bb.xmin,bb.ymax)
                end
                Cairo.close_path(cr)
                Cairo.set_source_rgba(cr,(:bgcolor in keys(Dict(Pl.options)) ?
                                      Dict(Pl.options)[:bgcolor].color :
                                      _DEFAULT_PLOT2D_KWARGS[:bgcolor].color)...,1)
                Cairo.fill(cr)
            end
            Cairo.set_line_width(cr,lwcorrect(P.pen.linewidth))
            Cairo.move_to(cr,P.points[1].x,P.points[1].y);
            for i=2:length(P.points)
                Cairo.line_to(cr,P.points[i].x,P.points[i].y);
            end
            Cairo.close_path(cr)
            if !is_no_pen(P.fillpen)
                Cairo.set_source_rgba(cr,P.fillpen.color...,P.fillpen.opacity)
                if !is_no_pen(P.pen)
                    Cairo.fill_preserve(cr)
                else
                    Cairo.fill(cr)
                end
            end
            if !is_no_pen(P.pen)
                Cairo.set_source_rgba(cr,P.pen.color...,P.pen.opacity)
                Cairo.stroke(cr)
            end
        end

        function addtocontext!(cr::Cairo.CairoContext,
                               Pl::Plot2D,
                               c::Circle2D,
                               bb::BoundingBox)
            if c.clip
                Cairo.arc(cr, c.center.x, c.center.y, c.radius, 0, 2*pi)
                Cairo.move_to(cr,bb.xmin,bb.ymin)
                Cairo.line_to(cr,bb.xmin,bb.ymax)
                Cairo.line_to(cr,bb.xmax,bb.ymax)
                Cairo.line_to(cr,bb.xmax,bb.ymin)
                Cairo.close_path(cr)
                Cairo.set_source_rgba(cr,Dict(Pl.options)[:bgcolor].color...,1)
                Cairo.fill(cr)
            end

            Cairo.set_line_width(cr,lwcorrect(c.pen.linewidth))
            Cairo.arc(cr, c.center.x, c.center.y, c.radius, 0, 2*pi)

            if c.fillpen.other ≠ "NoPen"
                Cairo.set_source_rgba(cr,c.fillpen.color...,c.fillpen.opacity)
                Cairo.fill_preserve(cr)
            end

            Cairo.set_source_rgba(cr,c.pen.color...,c.pen.opacity)
            Cairo.stroke(cr)
        end

        function addtocontext!(cr::Cairo.CairoContext,
                               Pl::Plot2D,
                               p::Point2D,
                               bb::BoundingBox)
            m = min(bb.xmax-bb.xmin,bb.ymax-bb.ymin)
            lw = p.pen.linewidth == 0 ? 3 : p.pen.linewidth
            Cairo.arc(cr,p.P.x,p.P.y,0.003*m*lw,0,2π)
            Cairo.set_source_rgb(cr,p.pen.color...)
            Cairo.fill(cr)
        end

        function addtocontext!(cr::Cairo.CairoContext,
                               Pl::Plot2D,
                               L::Label2D,
                               bb::BoundingBox)
            Cairo.save(cr)
            m = min(bb.xmax-bb.xmin,bb.ymax-bb.ymin)
            Cairo.set_font_size(cr,2.5e-3*m*L.pen.fontsize)
            Cairo.set_source_rgb(cr,L.pen.color...)
            x_bearing, y_bearing, width, height =
                                Cairo.text_extents(cr,L.s)
            Cairo.move_to(cr,L.location.x - width/2 - x_bearing,
                             L.location.y + height/2 + y_bearing)
            Cairo.scale(cr,1,-1)
            Cairo.show_text(cr,L.s)
            Cairo.stroke(cr)
            Cairo.restore(cr)
        end

        function addtocontext!(cr::Cairo.CairoContext,
                               Pl::Plot2D,
                               P::PixelMap,
                               bb::BoundingBox)
            m,n = size(P.pixels)
            a,b = P.lowerleft
            c,d = P.upperright
            for i=0:m-1
                for j=0:n-1
                    Cairo.rectangle(cr,a+(c-a)*i/m,b+(d-b)*j/n,(c-a)/m,(d-b)/n)
                    Cairo.set_source_rgba(cr,P.pixels[i+1,j+1].color...,P.alpha[i+1,j+1])
                    Cairo.fill(cr)
                end
            end
        end

        function bytes(P::Plot2D;format=:png,bbox=false,border=3)
            global _DEFAULT_PLOT2D_KWARGS
            D = Dict{Symbol,Any}(P.options)
            if :axes in keys(D) && D[:axes]
                P = Plot2D([axes(P);P.elements],P.options)
            end
            if :bgcolor in keys(D) && isa(D[:bgcolor],AbstractString)
                D[:bgcolor] = NamedColor(D[:bgcolor])
            end
            width = :width in keys(D) ? D[:width] : _DEFAULT_WIDTH
            border_fraction = 3*border/width
            bb = boundingbox(P,border=border_fraction)
            if !(bbox == false)
                bb = boundingbox([bb,bbox])
            end
            asp_ratio = aspectratio(bb)
            # height should be even, for ffmpeg:
            height = 2*round(Int,2*asp_ratio\width)
            bufferdata = UInt8[]
            iobuffer = IOBuffer(bufferdata,read=true,write=true)
            if format == :png
                surface = Cairo.CairoARGBSurface(width,height)
            elseif format == :pdf
                surface = Cairo.CairoPDFSurface(iobuffer,width,height)
            elseif format == :eps
                surface = Cairo.CairoEPSSurface(iobuffer,width,height)
            elseif format == :svg
                surface = Cairo.CairoSVGSurface(iobuffer,width,height)
            else
                error("format should be png, pdf, eps, or svg")
            end
            context = Cairo.CairoContext(surface)
            Cairo.scale(context,width/(bb.xmax-bb.xmin),-height/(bb.ymax-bb.ymin))
            Cairo.translate(context,-bb.xmin,-bb.ymax)
            if :bgcolor in keys(D)
                Cairo.set_source_rgb(context,D[:bgcolor].color...)
            else
                Cairo.set_source_rgb(context,1,1,1)
            end
            Cairo.paint(context)
            Cairo.set_line_width(context,1)
            Cairo.set_line_join(context, Cairo.Cairo.CAIRO_LINE_JOIN_ROUND)
            for e in P.elements
                addtocontext!(context,P,e,bb)
            end
            if format == :png
                Cairo.write_to_png(surface,iobuffer)
            end
            Cairo.finish(surface)
            Cairo.destroy(surface)
            return bufferdata
        end

    end

    Requires.@require SymPy="24249f21-da20-56a4-8eb1-6a02cf4ae2e6" begin
        function plot(S::SymPy.Sym,
                      t::Tuple{SymPy.Sym,U,V} where U<:Real where V<:Real;
                      kwargs...)
            plot(SymPy.lambdify(S,t[1:1]),t[2],t[3];kwargs...)
        end

        function plot(L::Vector,
                      t::Tuple{SymPy.Sym,U,V} where U<:Real where V<:Real;
                      kwargs...)
            plot([SymPy.lambdify(e,t[1:1]) for e in L],t[2],t[3];kwargs...)
        end

        function plot(L::Vector,x::SymPy.Sym,a::Real,b::Real;kwargs...)
            plot(L,(x,a,b);kwargs...)
        end

        SR = Union{Real,SymPy.Sym}

        plot(S::SymPy.Sym,
             a::SR,
             b::SR,
             c::SR,
             d::SR;
             kwargs...)  =
                 plot(SymPy.lambdify(S),map(SymPy.N,(a,b,c,d))...;kwargs...)

        plot(S::SymPy.Sym,
             x::Tuple{<:SR,<:SR},
             y::Tuple{<:SR,<:SR};kwargs...) = plot(f,x...,y...;kwargs...)

        plot(S::SymPy.Sym,
             x::Tuple{SymPy.Sym,<:SR,<:SR},
             y::Tuple{SymPy.Sym,<:SR,<:SR};kwargs...) =
                 plot(SymPy.lambdify(S,[x[1],y[1]]),x[2:end]...,y[2:end]...;kwargs...)

    end

    Requires.@require Juno="e5e0dc1b-0480-54bc-9374-aad01c23163d" begin

        _JUNO_ASY_FORMAT = "svg"
        _JUNO_CAIRO_FORMAT = "png"
        _FIT_DIMENSION = "height"

        """
            fitheight()

        Set the figures appearing in the Juno plot
        pane to match its height
        """
        fitheight() = (global _FIT_DIMENSION; _FIT_DIMENSION = "height")
        """
            fitwidth()

        Set the figures appearing in the Juno plot
        pane to match its width
        """
        fitwidth() = (global _FIT_DIMENSION; _FIT_DIMENSION = "width")

        """
            junosvg()

        Set the default output format for the Juno
        plot pane to svg
        """
        function junosvg()
            global _JUNO_ASY_FORMAT
            global _JUNO_CAIRO_FORMAT
            _JUNO_ASY_FORMAT = "svg"
            _JUNO_CAIRO_FORMAT = "svg"
        end
        """
            junopng()

        Set the default output format for the Juno
        plot pane to png
        """
        function junopng()
            global _JUNO_ASY_FORMAT
            global _JUNO_CAIRO_FORMAT
            _JUNO_ASY_FORMAT = "png"
            _JUNO_CAIRO_FORMAT = "png"
        end

        """
            junodefault()

        Reset the default output formats for the Juno
        plot pane
        """
        function junodefault()
            global _JUNO_ASY_FORMAT
            global _JUNO_CAIRO_FORMAT
            _JUNO_ASY_FORMAT = "svg"
            _JUNO_CAIRO_FORMAT = "png"
        end

        import Hiccup, Media
        Media.media(Plot2D,Media.Plot)
        Media.media(Plot3D,Media.Plot)

        function Juno.render(pane::Juno.PlotPane,P::Plot2D)
            global _BACKEND
            global _JUNO_ASY_FORMAT
            global _JUNO_CAIRO_FORMAT
            global _FIT_DIMENSION
            tempdir = mktempdir()
            w,h = Juno.plotsize()
            fitdim = _FIT_DIMENSION
            dimarg = (Symbol(fitdim) => (fitdim == "height" ? h : w),)
            if _BACKEND == "asy"
                filename = "$tempdir/myplot.$_JUNO_ASY_FORMAT"
                save(filename,P)
                D = Dict{Symbol,Any}(P.options)
                Juno.render(pane,Hiccup.div(
                    style="text-align:center",
                    Hiccup.img(src="$tempdir/myplot.$_JUNO_ASY_FORMAT";
                               dimarg...)))
            else
                save("$tempdir/myplot.$_JUNO_CAIRO_FORMAT",P)
                Juno.render(pane,Hiccup.div(
                    style="text-align:center",
                    Hiccup.img(src="$tempdir/myplot.$_JUNO_CAIRO_FORMAT";
                               dimarg...)))
            end
        end

        function Juno.render(pane::Juno.PlotPane,P::Plot3D)
            global _BACKEND
            if _BACKEND == "cairo"
                return
            end
            w,h = Juno.plotsize()
            tempdir = mktempdir()
            filename = "$tempdir/myplot.png"
            save(filename,P)
            Juno.render(pane,Hiccup.div(style="text-align:center",
                                        Hiccup.img(height=h,src="$tempdir/myplot.png")))
        end
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct WhatIfChange <: SwaggerModel
    resourceId::Any # spec type: Union{ Nothing, String } # spec name: resourceId
    changeType::Any # spec type: Union{ Nothing, String } # spec name: changeType
    before::Any # spec type: Union{ Nothing, Any } # spec name: before
    after::Any # spec type: Union{ Nothing, Any } # spec name: after
    delta::Any # spec type: Union{ Nothing, Vector{WhatIfPropertyChange} } # spec name: delta

    function WhatIfChange(;resourceId=nothing, changeType=nothing, before=nothing, after=nothing, delta=nothing)
        o = new()
        validate_property(WhatIfChange, Symbol("resourceId"), resourceId)
        setfield!(o, Symbol("resourceId"), resourceId)
        validate_property(WhatIfChange, Symbol("changeType"), changeType)
        setfield!(o, Symbol("changeType"), changeType)
        validate_property(WhatIfChange, Symbol("before"), before)
        setfield!(o, Symbol("before"), before)
        validate_property(WhatIfChange, Symbol("after"), after)
        setfield!(o, Symbol("after"), after)
        validate_property(WhatIfChange, Symbol("delta"), delta)
        setfield!(o, Symbol("delta"), delta)
        o
    end
end # type WhatIfChange

const _property_map_WhatIfChange = Dict{Symbol,Symbol}(Symbol("resourceId")=>Symbol("resourceId"), Symbol("changeType")=>Symbol("changeType"), Symbol("before")=>Symbol("before"), Symbol("after")=>Symbol("after"), Symbol("delta")=>Symbol("delta"))
const _property_types_WhatIfChange = Dict{Symbol,String}(Symbol("resourceId")=>"String", Symbol("changeType")=>"String", Symbol("before")=>"Any", Symbol("after")=>"Any", Symbol("delta")=>"Vector{WhatIfPropertyChange}")
Base.propertynames(::Type{ WhatIfChange }) = collect(keys(_property_map_WhatIfChange))
Swagger.property_type(::Type{ WhatIfChange }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_WhatIfChange[name]))}
Swagger.field_name(::Type{ WhatIfChange }, property_name::Symbol) =  _property_map_WhatIfChange[property_name]

const _allowed_WhatIfChange_changeType = ["Create", "Delete", "Ignore", "Deploy", "NoChange", "Modify"]

function check_required(o::WhatIfChange)
    (getproperty(o, Symbol("resourceId")) === nothing) && (return false)
    true
end

function validate_property(::Type{ WhatIfChange }, name::Symbol, val)
    if name === Symbol("changeType")
        Swagger.validate_param(name, "WhatIfChange", :enum, val, _allowed_WhatIfChange_changeType)
    end
end


text:
@testset "Degree" begin
    g5 = DiGraph(4)
    add_edge!(g5, 1, 2); add_edge!(g5, 2, 3); add_edge!(g5, 1, 3); add_edge!(g5, 3, 4)
    for g in testdigraphs(g5)
      @test @inferred(degree_centrality(g)) == [0.6666666666666666, 0.6666666666666666, 1.0, 0.3333333333333333]
      @test @inferred(indegree_centrality(g, normalize=false)) == [0.0, 1.0, 2.0, 1.0]
      @test @inferred(outdegree_centrality(g; normalize=false)) == [2.0, 1.0, 1.0, 0.0]
    end
end


text:
using Test

include("lexicographic-perms.jl")

@testset "#all_perms base cases" begin
  @test all_perms([]) == [[]]

  @test all_perms([0]) == [[0]]

  @test all_perms([1, 2]) == [[1, 2], [2, 1]]
end

@testset "#prepend_all!" begin
  @test prepend_all!(3, [[]]) == [[3]]

  @test prepend_all!(3, [[1, 2], [2, 1]]) == [[3, 1, 2], [3, 2, 1]]

  @test prepend_all!(4, [[1, 2, 3], [1, 3, 2],
                         [2, 1, 3], [2, 3, 1],
                         [3, 1, 2], [3, 2, 1]]) == [[4, 1, 2, 3], [4, 1, 3, 2],
                                                    [4, 2, 1, 3], [4, 2, 3, 1],
                                                    [4, 3, 1, 2], [4, 3, 2, 1]]

end

@testset "#all_perms rec. cases" begin
  @test all_perms([true, false]) ==  [Bool[1, 0], Bool[0, 1]]

  @test all_perms([1, 2, 3]) == [[1, 2, 3], [1, 3, 2], [2, 1, 3],
                                 [2, 3, 1], [3, 1, 2], [3, 2, 1]]

  @test all_perms([1, 'a', "foo"]) == [Any[1, 'a', "foo"], Any[1, "foo", 'a'],
                                       Any['a', 1, "foo"], Any['a', "foo", 1],
                                       Any["foo", 1, 'a'], Any["foo", 'a', 1]]

  @test all_perms(Vector{Int}([1, 2, 3])) == Vector{Vector{Int}}([[1, 2, 3], [1, 3, 2], [2, 1, 3],
                                                                  [2, 3, 1], [3, 1, 2], [3, 2, 1]])

  @test all_perms(Vector{Float64}([1., 2., 3.])) == Vector{Vector{Float64}}([
                                                                             [1., 2., 3.], [1., 3., 2.], [2., 1., 3.],
                                                                             [2., 3., 1.], [3., 1., 2.], [3., 2., 1.]])

  @test all_perms([1, 2, 3, 4]) == Vector{Vector{Int}}([
                                    [1, 2, 3, 4], [1, 2, 4, 3],
                                    [1, 3, 2, 4], [1, 3, 4, 2],
                                    [1, 4, 2, 3], [1, 4, 3, 2],
                                    [2, 1, 3, 4], [2, 1, 4, 3],
                                    [2, 3, 1, 4], [2, 3, 4, 1],
                                    [2, 4, 1, 3], [2, 4, 3, 1],
                                    [3, 1, 2, 4], [3, 1, 4, 2],
                                    [3, 2, 1, 4], [3, 2, 4, 1],
                                    [3, 4, 1, 2], [3, 4, 2, 1],
                                    [4, 1, 2, 3], [4, 1, 3, 2],
                                    [4, 2, 1, 3], [4, 2, 3, 1],
                                    [4, 3, 1, 2], [4, 3, 2, 1],
  ])

  @test all_perms(['v', 'x', 'y', 'z']) == [['v', 'x', 'y', 'z'], ['v', 'x', 'z', 'y'],
                                            ['v', 'y', 'x', 'z'], ['v', 'y', 'z', 'x'],
                                            ['v', 'z', 'x', 'y'], ['v', 'z', 'y', 'x'],
                                            ['x', 'v', 'y', 'z'], ['x', 'v', 'z', 'y'],
                                            ['x', 'y', 'v', 'z'], ['x', 'y', 'z', 'v'],
                                            ['x', 'z', 'v', 'y'], ['x', 'z', 'y', 'v'],
                                            ['y', 'v', 'x', 'z'], ['y', 'v', 'z', 'x'],
                                            ['y', 'x', 'v', 'z'], ['y', 'x', 'z', 'v'],
                                            ['y', 'z', 'v', 'x'], ['y', 'z', 'x', 'v'],
                                            ['z', 'v', 'x', 'y'], ['z', 'v', 'y', 'x'],
                                            ['z', 'x', 'v', 'y'], ['z', 'x', 'y', 'v'],
                                            ['z', 'y', 'v', 'x'], ['z', 'y', 'x', 'v']
  ]

  @test all_perms(collect(0:4)) |> length == 120

  @test all_perms(collect(0:5)) |> length == 720

  @test all_perms(collect(0:6)) |> length == 5040
end

@time @testset "#all_perms of 10 digits => 10!" begin
  # 10! ≡ 3_628_800
  @test_skip all_perms(collect(0:9)) |> length == 3_628_800
end
#
# Test Summary:                  | Pass  Total
# #all_perms of 10 digits => 10! |    1      1
#   6.666260 seconds (36.33 M allocations: 3.653 GiB, 26.40% gc time, 5.47% compilation time)
#

@time @testset "#all_perms of 11 digits => 11!" begin
  # 11! ≡ 39_916_800
  @test_skip all_perms(collect(0:10)) |> length == 39_916_800
end
#
# Test Summary:                  | Pass  Total
# #all_perms of 11 digits => 11! |    1      1
# 169.726583 seconds (284.85 M allocations: 31.599 GiB, 46.52% gc time, 1.47% compilation time)
#

# way too long...
# @time @testset "#all_perms of 12 digits => 12!" begin
#   # 12! ≡ 479_001_600
#   @test_skip all_perms(collect(0:11)) |> length == 479_001_600
# end

@time @testset "millionnth permutation of 0..9" begin
  @test millionnth(collect(0:9)) == [9, 8, 7, 6, 5, 4, 3, 2, 1, 0]
end

max_stars_repo_name: pascal-p/Euler

max_issues_repo_name: pascal-p/Euler

max_forks_repo_name: pascal-p/Euler


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule ldid_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("ldid")
JLLWrappers.@generate_main_file("ldid", UUID("df1af0dd-2f85-5d2f-b099-55d224f7db60"))
end  # module ldid_jll


text:
## assets -------------------------------------------------------
function IGCost(df_cost)
  IGCost(df_cost[:, :rel],
         df_cost[:, :abs],
         df_cost[:, :infl_rel],
         df_cost[:, :infl_abs],
         cumprod(1 .+ df_cost[:, :infl_rel]),
         cumprod(1 .+ df_cost[:, :infl_abs]),
         zeros(Float64, nrow(df_cost))
         )
end

function IGStock(cap_mkt::CapMkt, mv_0, alloc, cost)
  investments = Array{InvestStock}(undef, 0)
  for 𝑖 ∈ 1:length(alloc.name)
    push!(investments,
          InvestStock(alloc.name[𝑖],
                      cap_mkt.stock,
                      mv_0 * alloc.total[1] * alloc.all[1,𝑖],
                      zeros(Float64, size(alloc.all, 1))))
  end
  IGStock(investments, mv_0 * alloc.total[1],
          zeros(Float64, size(alloc.all, 1)), alloc,
          deepcopy(cost))
end

function IGCash(cap_mkt::CapMkt, mv_0, alloc, cost)
  investments = Array{InvestCash}(undef, 0)
  for 𝑖 ∈ 1:length(alloc.name)
    push!(investments,
          InvestCash(alloc.name[𝑖],
                     cap_mkt.rfr,
                     mv_0 * alloc.total[1] * alloc.all[1,𝑖],
                     zeros(Float64, size(alloc.all, 1)),
                     alloc.lgd[𝑖],
                     alloc.cqs[𝑖]))
  end
  IGCash(investments, mv_0 * alloc.total[1],
         zeros(Float64, size(alloc.all, 1)), alloc,
         deepcopy(cost))
end

function InvPort(t_0,
                 dur,
                 cap_mkt::CapMkt,
                 mv_0,
                 allocs::Dict{Symbol, Alloc},
                 costs::Dict{Symbol, DataFrame}
                 )
  igs = Dict{Symbol, InvestGroup}()
  for 𝑖𝑔_𝑠𝑦𝑚𝑏 ∈ collect(keys(allocs))
    ## ig_symb are the symbols corresponding to the
    ## types of investment groups: :IGCash, IGStock
    merge!(igs,
           Dict(𝑖𝑔_𝑠𝑦𝑚𝑏 => eval(𝑖𝑔_𝑠𝑦𝑚𝑏)(cap_mkt,
                                        mv_0,
                                        allocs[𝑖𝑔_𝑠𝑦𝑚𝑏],
                                        IGCost(costs[𝑖𝑔_𝑠𝑦𝑚𝑏])
                                         )))
  end
  return(InvPort(t_0,                  ## start of projection
                 mv_0,                 ## init. mv pre prem..
                 zeros(Float64, dur),  ## mv bop post prem...
                 zeros(Float64, dur),  ## mv eop
                 zeros(Float64, dur),  ## average yield
                 zeros(Float64, dur),  ## investment costs
                 igs))                 ## investment groups
end

## insurance liabilities ----------------------------------------
function ProductIns(rfr_price, prob_price, β_in, λ_price)
  dur = nrow(β_in)
  prob = deepcopy(prob_price)
  prob[!,:px] = 1 .- prob[!,:qx] - prob[!,:sx]
  λ_price[!,:cum_infl] = cumprod(1 .+ λ_price[!,:infl])
  return ProductIns(dur, rfr_price, prob, β_in, λ_price,
                 premium(1, rfr_price, prob, β_in, λ_price))
end

function ModelPoint(n, t_0, t_start,
                    prob_be, sx_be_fac, λ_be,
                    cost_infl,
                    hypo_bonus_rate, product, ins_sum,
                    pension_contract)
  ## Time model for model point: See documentation of type
  s_0 = t_0 - t_start
  dur = product.dur - s_0
  s_future = (s_0 + 1):product.dur
  prob = deepcopy(prob_be)[s_future, :]
  prob[!, :sx] *= sx_be_fac
  prob[!, :px] = 1 .- prob[!,:qx] .- prob[!,:sx]
  lx_boy = zeros(Float64, dur)
  β = DataFrame()
  for 𝑛𝑎𝑚𝑒 ∈ names(product.β)
    β[!,𝑛𝑎𝑚𝑒] = n * ins_sum * product.β[s_future, 𝑛𝑎𝑚𝑒]
  end
  β[!,:prem] *= product.prem_norm
  β[!,:sx] *= product.prem_norm

  λ = deepcopy(λ_be)[s_future, :]
  λ[!, :boy] *= n * ins_sum
  λ[!, :eoy] *= n * ins_sum
  ## be cost inflation input relates to t_0 not s_0:
  λ[!,:infl] = deepcopy(cost_infl)
  λ[!,:cum_infl] = cumprod(1 .+ λ[!,:infl])
  λ_price = deepcopy(product.λ)[s_future, :]
  λ_price[!,:boy] *= n * ins_sum
  λ_price[!,:eoy] *= n * ins_sum
  λ_price[!,:cum_infl] = cumprod(1 .+ product.λ[!,:infl])[s_future]
  rfr_price_0 = product.rfr[s_0 == 0 ? 1 : s_0]
  rfr_price = product.rfr[s_future]
  tpg_price_0 = tpg(0,
                    rfr_price,
                    product.prob[s_future, :],
                    β,
                    λ_price)
  tpg_price = zeros(Float64, dur)
  for 𝑡 ∈ 1:dur
    tpg_price[𝑡] = tpg(𝑡,
                       rfr_price,
                       product.prob[s_future, :],
                       β,
                       λ_price)
  end
  return ModelPoint(n, t_start, dur, prob, lx_boy, 0.0,
                    β, λ,  hypo_bonus_rate,
                    rfr_price_0, rfr_price,
                    tpg_price_0, tpg_price,
                    ones(Float64, dur),
                    pension_contract)
end

function LiabIns(t_0::Int, prob_be, λ_be,
                 cost_infl, product, df_port)
  n = nrow(df_port)
  mps = Array{ModelPoint}(undef, 0)
  dur = 0
  for 𝑑 ∈ 1:n
    push!(mps, ModelPoint(df_port[𝑑, :n],
                          t_0,
                          df_port[𝑑, :t_start],
                          prob_be,
                          df_port[𝑑, :sx_be_fac],
                          λ_be,
                          cost_infl[𝑑],
                          df_port[𝑑, :bonus_rate_hypo],
                          product,
                          df_port[𝑑, :ins_sum],
                          df_port[𝑑, :pension_contract]))
    dur = max(dur, mps[𝑑].dur)
  end
  gc = zeros(Float64, dur)
  for 𝑚𝑝 ∈ mps
    𝑚𝑝.gc = zeros(Float64, dur)
    𝑚𝑝.gc[1:𝑚𝑝.dur] +=
      vcat(1, cumprod(𝑚𝑝.prob[1:(𝑚𝑝.dur-1), :px]))
    gc +=  𝑚𝑝.n * 𝑚𝑝.gc
  end
  gc /= gc[1]
  Δgc = diff(vcat(gc, 0))
  return LiabIns(n, t_0, dur, mps, gc, Δgc)
end

## other liabilities --------------------------------------------
function Debt(t_0, t_debt_0, t_debt_mat,
              name::Symbol, nominal, coupon)
  name = name
  τ_debt_0 = t_debt_0 - t_0
  dur_debt = t_debt_mat - t_debt_0 + 1
  τ_mat = dur_debt + τ_debt_0 - 1
  Debt(name, t_debt_0, t_debt_mat, τ_debt_0, τ_mat, nominal, coupon)
end

function Debt(t_0, df_debt::DataFrame)
  name = df_debt[1, :name]
  t_init = df_debt[1, :t_init]
  t_mat = df_debt[1, :t_mat]
  τ_init = t_init - t_0
  dur = t_mat - t_init + 1
  τ_mat = dur + τ_init - 1
  Debt(name, t_init, t_mat, τ_init, τ_mat,
       df_debt[1, :nominal], df_debt[1, :coupon])
end

function Debt(t_0, df_debt::DataFrameRow)
  name = df_debt[ :name]
  t_init = df_debt[ :t_init]
  t_mat = df_debt[ :t_mat]
  τ_init = t_init - t_0
  dur = t_mat - t_init + 1
  τ_mat = dur + τ_init - 1
  Debt(name, t_init, t_mat, τ_init, τ_mat,
       df_debt[ :nominal], df_debt[ :coupon])
end

function LiabOther(t_0, df_debts::DataFrame)
  subord = Array{Debt}(undef, nrow(df_debts))
  for 𝑑 ∈ 1:nrow(df_debts)
    subord[𝑑] = Debt(t_0, df_debts[𝑑,:])
  end
  return LiabOther(subord)
end

## dynamics -----------------------------------------------------
function Dynamic(dur, bonus_factor, quota_surp::Float64)
  Dynamic(bonus_factor,
          quota_surp,
          zeros(Float64, dur))
end

## cashflow projection ------------------------------------------
function Projection(liabs, tax_rate, tax_credit_0)
  t_0 = liabs.t_0
  dur = liabs.dur
  cf = DataFrame(
    qx = zeros(Float64, dur),
    sx = zeros(Float64, dur),
    px = zeros(Float64, dur),
    prem = zeros(Float64, dur),
    λ_boy = zeros(Float64, dur),
    λ_eoy = zeros(Float64, dur),
    bonus = zeros(Float64, dur),
    invest = zeros(Float64, dur),
    new_debt = zeros(Float64, dur),
    l_other = zeros(Float64, dur),
    profit = zeros(Float64, dur),
    tax = zeros(Float64, dur),
    divid = zeros(Float64, dur),
    gc = zeros(Float64, dur),        ## not affecting profit/loss
    Δtpg = zeros(Float64, dur),      ## not real cf, affects p/l
    cost_prov = zeros(Float64, dur)  ## for the cost provisions
    )
  val = DataFrame(
    invest = zeros(Float64, dur),
    tpg = zeros(Float64, dur),
    l_other = zeros(Float64, dur),
    surplus = zeros(Float64, dur),
    bonus = zeros(Float64, dur),
    cost_prov = zeros(Float64, dur)      ## cost provisions
    )
  val_0 = DataFrame(deepcopy(val[1, :]))
  return Projection(t_0, dur, cf, val_0, val,
                    tax_rate,
                    tax_credit_0,
                    zeros(Float64, dur),
                    zeros(Float64, dur))
end

function Projection(tax_rate,
                    tax_credit_0,
                    cap_mkt::CapMkt,
                    invs::InvPort,
                    liabs::LiabIns,
                    liabs_other::LiabOther,
                    dyn::Dynamic)
  proj = Projection(liabs, tax_rate, tax_credit_0)
  for 𝑡 ∈ 1:liabs.dur
    for 𝑖𝑔 ∈ [:IGCash, :IGStock]
      proj.fixed_cost_gc[𝑡] +=
        invs.igs[𝑖𝑔].cost.abs[𝑡] *
        invs.igs[𝑖𝑔].cost.cum_infl_abs[𝑡] *
        liabs.gc[𝑡]
    end
  end
  l_other = deepcopy(liabs_other)
  goingconcern!(l_other, liabs.Δgc)
  val0!(cap_mkt, invs, liabs, l_other, proj)
  proj.cf[:,:gc] = liabs.Δgc * (proj.val_0[1,:invest] -
                                  proj.val_0[1,:tpg] -
                                  proj.val_0[1,:l_other])
  for 𝑡 = 1:liabs.dur
    project!(𝑡, cap_mkt, invs, liabs, l_other, dyn, proj)
  end
  valbonus!(cap_mkt.rfr.x, proj)
  valcostprov!(cap_mkt.rfr.x, invs, proj)
  return proj
end


text:
using Documenter, DocStringExtensions, DynamicEnergyBudgets

makedocs(
    modules = [DynamicEnergyBudgets],
    doctest = false,
    clean = false,
    sitename = "DynamicEnergyBudgets.jl",
    format = :html,
    pages = Any[
        "Home" => "index.md",
    ]
)

deploydocs(
    repo = "github.com/rafaqz/DynamicEnergyBudgets.jl.git",
    osname = "linux",
    julia = "0.6",
    target = "build",
    deps = nothing,
    make = nothing
)


text:
module SpkPlot

using PyPlot, Colors, ColorTypes

using ..Spk, ..Histogram

include("./plot_helpers.jl")

export plot_xcorr, plot_raster, plot_raster!, plot_cycle_mean, plot_cycle_mean!

export cla, default_figure, default_axes

# ============================================================================ #
cla() = close("all")
# ============================================================================ #
function plot_xcorr(t::AbstractVector, xc::Vector; bin_size=0.0005, tmax=0.015,
    color="black")

    h = qbar(t, xc, bin_size*.8, color)
    xlabel("Time lag (sec)")
    ylabel("# of spikes")
    tmax += bin_size
    xlim(-tmax, tmax)
    return h
end
# ---------------------------------------------------------------------------- #
function plot_xcorr(ts1::Vector{T}, ts2::Vector{T};
            bin_size=0.0005, tmax=0.02, color="black", shift=0, tf=4.0,
            shift_color="cyan") where {T<:AbstractFloat}

    xc, t = get_xcorr(ts1, ts2, bin_size, tmax)
    h = plot_xcorr(t, xc, bin_size=bin_size, tmax=tmax, color=color)
    if shift > 0
        sp = get_shift(ts1, ts2, bin_size, tmax, shift, tf)
        plot(t, sp, color=shift_color, linewidth=3)
    end
    return h
end
# ---------------------------------------------------------------------------- #
function plot_xcorr(ts1::Vector{T}; bin_size=0.0005, tmax=0.02, color="black",
    shift::Int=0, tf=4.0) where {T<:AbstractFloat}

    xc, t = get_xcorr(ts1, bin_size, tmax)
    h = plot_xcorr(t, xc, bin_size=bin_size, tmax=tmax, color=color)
    if shift > 0
        sp = get_shift(ts1, ts2, bin_size, tmax, shift, tf)
        plot(t, sp, color=shift_color, linewidth=3)
    end
    return h
end
# ============================================================================ #
function plot_raster!(ax, ts::Vector{T}, evt::AbstractArray{T,1}, dur::AbstractFloat,
    color::String="blue") where {T<:AbstractFloat}

    bin_size = 0.001
    nbin = floor(Int, dur / bin_size)
    p = psth(ts, evt, 0:nbin-1, bin_size)[1];

    for k = 1:size(p,2)
       idx = findall(p[:,k] .> 0)
       ax.plot(idx * bin_size, fill(k, length(idx)), ".", color=color)
   end
   ax.set_xlabel("Time (sec)", fontsize=16)
   ax.set_ylabel("Trial #", fontsize=16)
   return ax
end
function plot_raster(ts::Vector{T}, evt::AbstractArray{T,1}, dur::AbstractFloat,
    color::String="blue") where {T<:AbstractFloat}

    h = default_figure()
    return plot_raster!(h.axes[1], ts, evt, dur, color)
end
# ============================================================================ #
function plot_cycle_mean!(ax, ts::Vector{T}, evt::AbstractArray{T,1},
    dur::AbstractFloat, tf::AbstractFloat; bin_size::Real=0.01,
    color::String="blue") where {T<:AbstractFloat}

    nbin = floor(Int, dur / bin_size)
    bpc = floor(Int, 1/(tf * bin_size))
    p = psth(ts, evt, 0:nbin-1, bin_size)[1];

    cm = cycle_mean(sum(p, dims=2), bpc)
    h = ax.plot(range(0, stop=1/tf, length=length(cm)), cm, color=color)
    ax.set_xlabel("Time (sec)", fontsize=16)
    ax.set_ylabel("Firing rate (Hz)", fontsize=16)
    return h
end
# ---------------------------------------------------------------------------- #
function plot_cycle_mean(ts::Vector{T}, evt::AbstractArray{T,1},
    dur::AbstractFloat, tf::AbstractFloat; bin_size::Real=0.01,
    color::String="blue") where {T<:AbstractFloat}

    h = default_figure()

    return plot_cycle_mean!(h.axes[1], ts, evt, dur, tf, bin_size=bin_size,
        color=color)
end
# ============================================================================ #
end


text:
export orthog!,
       svd_recursive

function orthog!(M::AbstractMatrix{T};
                 npass::Int=2) where {T}
  nkeep = min(size(M)...)
  dots = zeros(T,nkeep)
  for i=1:nkeep
    coli = view(M,:,i)
    nrm = norm(coli)
    if nrm < 1E-10
      rand!(coli)
      nrm = norm(coli)
    end
    coli ./= nrm
    (i==1) && continue

    Mcols = view(M,:,1:i-1)
    dotsref = view(dots,1:i-1)
    for pass=1:npass
      mul!(dotsref,Mcols',coli)
      #BLAS.gemv!('N',1.0,Mcols,dotsref,-1.0,coli)
      coli .-= Mcols*dotsref
      nrm = norm(coli)
      if nrm < 1E-3 #orthog is suspect
        pass = pass-1
      end
      if nrm < 1E-10
        rand!(coli)
        nrm = norm(coli)
      end
      coli ./= nrm
    end
  end
end

function pos_sqrt(x::Float64)::Float64
  (x < 0.0) && return 0.0
  return sqrt(x)
end

function checkSVDDone(S::Vector,
                      thresh::Float64)
  N = length(S)
  (N <= 1 || thresh < 0.0) && return (true,1)
  S1t = S[1]*thresh
  start = 2
  while start <= N
    (S[start] < S1t) && break
    start += 1
  end
  if start >= N
    return (true,N)
  end
  return (false,start)
end

function svd_recursive(M::AbstractMatrix;
                      thresh::Float64=1E-3,
                      north_pass::Int=2)
  Mr,Mc = size(M)

  if Mr > Mc
    V,S,U = svd_recursive(transpose(M))
    conj!(U)
    conj!(V)
    return U,S,V
  end

  #rho = BLAS.gemm('N','T',-1.0,M,M) #negative to sort eigenvalues greatest to smallest
  rho = -M*M' #negative to sort eigenvalues in decreasing order
  D,U = eigen(Hermitian(rho),1:size(rho,1))

  Nd = length(D)
  for n=1:Nd
    D[n] = pos_sqrt(-D[n])
  end

  V = M'*U
  orthog!(V,npass=north_pass)

  (done,start) = checkSVDDone(D,thresh)

  done && return U,D,V

  u = view(U,:,start:Nd)
  v = view(V,:,start:Nd)

  b = u'*(M*v)
  bu,bd,bv = svd_recursive(b,
                          thresh=thresh,
                          north_pass=north_pass)

  u .= u*bu
  v .= v*bv
  view(D,start:Nd) .= bd

  return U,D,V
end

# TODO: maybe move to another location?
function polar(M::AbstractMatrix)
  U,S,V = svd(M) # calls LinearAlgebra.svd()
  return U*V',V*Diagonal(S)*V'
end


text:
module TestUtils

using Test
using YaoCompiler
using YaoCompiler.Gate

@device function demo_circ_simp()
    1 => shift(7π / 4)
    1 => H
    1 => Rx(π / 4)
    4 => H
    @ctrl 1 4 => Z
    @ctrl 4 1 => X
    1 => H
    4 => H
    1 => T
    4 => shift(3π / 2)
    4 => X
    1 => H
    4 => S
    4 => X
    2 => S
    @ctrl 2 3 => X
    2 => H
    @ctrl 2 3 => X
    2 => T
    3 => S
    2 => H
    3 => H
    3 => S
    @ctrl 2 3 => X
end

YaoCompiler.gate_count(demo_circ_simp())

@testset "gate_count" begin
    @test gate_count(demo_circ_simp()) == Dict(
        :ctrl => IdDict(Z => 1, X => 4),
        :gate => IdDict(T => 2, H => 8, Rx => 1, X => 2, S => 4, shift => 2),
    )
end

end # TestUtils