text:
using Oceananigans
using Oceananigans.Fields
using Oceananigans.Operators
using Oceananigans.BoundaryConditions

using CUDA
import Oceananigans.Utils: launch!

using KernelAbstractions
using Oceananigans.Architectures: device, CPU, GPU, @hascuda

arch = CPU()
@hascuda arch = GPU()

@kernel function ∇²!(grid, f, ∇²f)
    i, j, k = @index(Global, NTuple)
    @inbounds ∇²f[i, j, k] = ∇²(i, j, k, grid, f)
end

Lx, Ly, Lz = 1e3, 1e3, 1
Nx, Ny, Nz =  35,  35, 1
grid = RegularCartesianGrid(size=(Nx,Ny,Nz), extent=(Lx,Ly,Lz) )
x = TracerFields((:c,), arch, grid)
div_U=x.c
#
struct PCGSolver{A, S}
        architecture :: A
            settings :: S
end

function PCGSolver(;arch=arch, parameters=parameters)
          bcs   = parameters.Template_field.boundary_conditions
          grid  = parameters.Template_field.grid
          maxit = parameters.maxit
          tol   = parameters.tol
          a_res = similar(parameters.Template_field.data);a_res.=0.
          q     = similar(parameters.Template_field.data)
          p     = similar(parameters.Template_field.data)
          z     = similar(parameters.Template_field.data)
          r     = similar(parameters.Template_field.data)
          if parameters.PCmatrix_function == nothing
            # preconditioner not provided, use the Identity matrix
            PCmatrix_function(x) = ( return x )
          else
            PCmatrix_function = parameters.PCmatrix_function
          end
          ii=grid.Hx:grid.Nx+grid.Hx-1
          ji=grid.Hy:grid.Ny+grid.Hy-1
          ki=grid.Hz:grid.Nz+grid.Hz-1
          dotproduct(x,y)  = mapreduce((x,y)->x*y, + , x[ii,ji,ki], y[ii,ji,ki])
          norm(x)          = ( mapreduce((x)->x*x, + , x[ii,ji,ki]   ) )^0.5
          Amatrix_function = parameters.Amatrix_function
          A(x) = ( Amatrix_function(x,a_res,bcs); return  a_res )
          M(x) = ( PCmatrix_function(x) )
          settings = (q=q,
                      p=p,
                      z=z,
                      r=r,
                    bcs=bcs,
                   grid=grid,
                      A=A,
                      M=M,
                  maxit=maxit,
                    tol=tol,
                dotprod=dotproduct,
                   norm=norm,
                   arch=arch,
          )

   return PCGSolver(arch, settings)
end

function solve_poisson_equation!(solver::PCGSolver,RHS,x)
#
# Alg - see Fig 2.5 The Preconditioned Conjugate Gradient Method in
#                    "Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods"
#                    Barrett et. al, 2nd Edition.
#
#     given
#        linear Preconditioner operator M as a function M()
#        linear A matrix operator A as a function A()
#        a dot product function norm()
#        a right-hand side b
#        an initial guess x
#
#        local vectors: z, r, p, q
#
#     β  = 0
#     r .= b-A(x)
#     i  = 0
#
#     loop:
#      if i > MAXIT
#       break
#      end
#      z = M( r )
#      ρ    .= dotprod(r,z)
#      p = z+β*p
#      q = A(p)
#      α=ρ/dotprod(p,q)
#      x=x.+αp
#      r=r.-αq
#      if norm2(r) < tol
#       break
#      end
#      i=i+1
#      ρⁱᵐ1 .= ρ
#      β    .= ρⁱᵐ¹/ρ
#
      sset       = solver.settings
      z, r, p, q = sset.z, sset.r, sset.p, sset.q
      A          = sset.A
      M          = sset.M
      maxit      = sset.maxit
      tol        = sset.tol
      dotprod    = sset.dotprod
      norm       = sset.norm

      β    = 0.
      r   .= RHS .- A(x)
      i    = 0
      ρ    = 0
      ρⁱᵐ¹ = 0

      while true
       if i > maxit
        break
       end
       z    .= M(r)
       ρ     = dotprod(z,r)
       if i == 0
        p   .= z
       else
        β    = ρ/ρⁱᵐ¹
        p   .= z .+ β .* p
       end
       q    .= A(p)
       α     = ρ/dotprod(p,q)
       x    .= x .+ α .* p
       r    .= r .- α .* q
       println("Solver ", i," ", norm(r) )
       if norm(r) <= tol
        break
       end
       i     = i+1
       ρⁱᵐ¹  = ρ
      end
#==
#     No preconditioner verison
      i    = 0
      r   .= RHS .- A(x)
      p   .= r
      γ    = dotprod(r,r)
      while true
       if i > maxit
        break
       end
       q   .= A(p)
       α    = γ/dotprod(p,q)
       x   .= x .+ α .* p
       r   .= r .- α .* q
       println("Solver ", i," ", norm(r) )
       if norm(r) <= tol
        break
       end
       γ⁺   = dotprod(r,r)
       β    = γ⁺/γ
       p   .= r .+ β .* p
       γ    = γ⁺
       i    = i+1
      end
==#

      fill_halo_regions!(x, sset.bcs, sset.arch, sset.grid)
      return x, norm(r)
end

function Base.show(io::IO, solver::PCGSolver)
        print(io, "Oceanigans compatible preconditioned conjugate gradient solver.\n")
        print(io, " Problem size = "  , size(solver.settings.q) )
        print(io, "\n Boundary conditions = "  , solver.settings.bcs  )
        print(io, "\n Grid = "  , solver.settings.grid  )
  return nothing
end

function fhr!(x,bcs)
   fill_halo_regions!(x, bcs, arch, grid)
end
function Amatrix_function(x,result,bcs)
 event = launch!(arch, grid, :xyz, ∇²!, grid, x, result, dependencies=Event(device(arch)))
 wait(device(arch), event)
 fhr!(result,bcs)
end

pcg_solver=PCGSolver( ;arch=arch,
              parameters=(PCmatrix_function=nothing,
                          Amatrix_function= Amatrix_function,
                          Template_field=div_U,
                          maxit=grid.Nx*grid.Ny,
                          tol=1.e-16,
                         )
           )

x=similar(div_U.data);x.=0
div_U.data.=0.
imid=Int(floor(grid.Nx/2))+1
jmid=Int(floor(grid.Ny/2))+1
div_U.data[imid-1,jmid,1]=-1
div_U.data[imid,jmid,1]= 1
fhr!(div_U.data,div_U.boundary_conditions)
solve_poisson_equation!(pcg_solver,div_U.data,x)


text:
@safetestset "Special chars in GET params (query)" begin
  @safetestset "<a+b> should be <a b>" begin
    using Genie
    using HTTP

    port = nothing
    port = rand(8500:8900)

    route("/") do
      params(:x)
    end

    server = up(port)

    response = try
      HTTP.request("GET", "http://127.0.0.1:$port/?x=foo+bar")
    catch ex
      ex.response
    end

    @test response.status == 200
    @test String(response.body) == "foo bar"

    down()
    sleep(1)
    server = nothing
    port = nothing
  end;

  @safetestset "<a%20b> should be <a b>" begin
    using Genie
    using HTTP

    port = nothing
    port = rand(8500:8900)

    route("/") do
      params(:x)
    end

    server = up(port)

    response = try
      HTTP.request("GET", "http://127.0.0.1:$port/?x=foo%20bar")
    catch ex
      ex.response
    end

    @test response.status == 200
    @test String(response.body) == "foo bar"

    down()
    sleep(1)
    server = nothing
    port = nothing
  end;

  @safetestset "<a%2Bb> should be <a+b>" begin
    using Genie
    using HTTP

    port = nothing
    port = rand(8500:8900)

    route("/") do
      params(:x)
    end

    server = up(port)

    response = try
      HTTP.request("GET", "http://127.0.0.1:$port/?x=foo%2Bbar")
    catch ex
      ex.response
    end

    @test response.status == 200
    @test String(response.body) == "foo+bar"

    down()
    sleep(1)
    server = nothing
    port = nothing
  end;

  @safetestset "emoji support" begin
    using Genie
    using HTTP

    port = nothing
    port = rand(8500:8900)

    route("/") do
      params(:x)
    end

    server = up(port)

    response = try
      HTTP.request("GET", "http://127.0.0.1:$port/?x=✔+🧞+♥")
    catch ex
      ex.response
    end

    @test response.status == 200
    @test String(response.body) == "✔ 🧞 ♥"

    down()
    sleep(1)
    server = nothing
    port = nothing
  end;

end

ext: jl


text:
###########
# Library #
###########
# using Random
using MatrixLMnet
using Helium
using Test

#####################################################################
# TEST Cross Validation Lasso vs Elastic Net (𝛼=1) - Simulated Data #
#####################################################################

#=
Description:
-----------

Model: 𝐘 = 𝐗 𝛃 𝐙' + 𝜎𝜖, with 𝜖∼𝑁(0,1)

Simulate data set consisting  of 20/20/200 observations and 8 predictors.
We let 𝛽₁ = (3, 1.5, 0, 0, 2, 0, 0, 0), 𝛽₂ = (0, 1.5, 0, 3.5, 2, 0, 0 , 2) where
𝛃 = [𝛽₁, 𝛽₂] and 𝜎 = 3.
The pairwise correlation between 𝑋ᵢ and 𝑋ⱼ was set to be 𝑐𝑜𝑟(𝑖,𝑗)=(0.5)^|𝑖−𝑗|.
Here, the Z matrix is an identity matrix.
=#

# Data testing directory name
dataDir = realpath(joinpath(@__DIR__,"data"))

# Get predictors
X = Helium.readhe(joinpath(dataDir, "Xmat.he"))

# Get response
Y = Helium.readhe(joinpath(dataDir, "Ymat.he"))

# Get Z matrix
Z = Helium.readhe(joinpath(dataDir, "Zmat.he"))

# Build raw data object from MatrixLM.jl
dat = RawData(Response(Y), Predictors(X, Z));

# Hyper parameters
λ = [10.0, 5.0, 3.0]
α = [1.0]

#############################################
# TEST 1 Lasso vs Elastic Net (𝛼=1) - ista #
#############################################

# Elastic net penalized regression
MatrixLMnet.Random.seed!(2021)
est1 = MatrixLMnet.mlmnet_cv(dat, λ, α, 10, 1, method = "ista", hasZIntercept = false, hasXIntercept = false, isVerbose = false);
smmr_Net1 = MatrixLMnet.lambda_min(est1);

# Elastic net penalized regression
MatrixLMnet.Random.seed!(2021)
est2 = MatrixLMnet.mlmnet_cv(dat, λ, 10, 1, method = "ista",  hasZIntercept = false, hasXIntercept = false, isVerbose = false);
smmr_Net2 = MatrixLMnet.lambda_min(est2);

# Lasso penalized regression - ista cv
smmr_ista= Helium.readhe(joinpath(dataDir, "smmr_ista.he"))

println("CV Lasso vs Elastic Net when α=1 test 1 - ista: ",
         @test smmr_Net1.AvgMSE == smmr_ista[:,1] && smmr_Net2.AvgPercentZero == smmr_ista[:,2] &&
               smmr_Net2.AvgMSE == smmr_ista[:,1] && smmr_Net2.AvgPercentZero == smmr_ista[:,2] )

#############################################
# TEST 2 Lasso vs Elastic Net (𝛼=1) - fista #
#############################################

# Elastic net penalized regression
MatrixLMnet.Random.seed!(2021)
est1 = MatrixLMnet.mlmnet_cv(dat, λ, α, 10, 1, method = "fista", hasZIntercept = false, hasXIntercept = false, isVerbose = false);
smmr_Net1 = MatrixLMnet.lambda_min(est1);

# Elastic net penalized regression
MatrixLMnet.Random.seed!(2021)
est2 = MatrixLMnet.mlmnet_cv(dat, λ, 10, 1, method = "fista",  hasZIntercept = false, hasXIntercept = false, isVerbose = false);
smmr_Net2 = MatrixLMnet.lambda_min(est2);

# Lasso penalized regression - fista cv
smmr_fista= Helium.readhe(joinpath(dataDir, "smmr_fista.he"))

println("CV Lasso vs Elastic Net when α=1 test 2 - fista: ",
         @test smmr_Net1.AvgMSE == smmr_fista[:,1] && smmr_Net1.AvgPercentZero == smmr_fista[:,2] &&
               smmr_Net2.AvgMSE == smmr_fista[:,1] && smmr_Net2.AvgPercentZero == smmr_fista[:,2] )

##########################################################
# TEST 3 Lasso vs Elastic Net (𝛼=1) - fista backtracking #
##########################################################

# Elastic net penalized regression
MatrixLMnet.Random.seed!(2021)
est1 = MatrixLMnet.mlmnet_cv(dat, λ, α, 10, 1, method = "fista_bt", hasZIntercept = false, hasXIntercept = false, isVerbose = false);
smmr_Net1 = MatrixLMnet.lambda_min(est1);

# Elastic net penalized regression
MatrixLMnet.Random.seed!(2021)
est2 = MatrixLMnet.mlmnet_cv(dat, λ, 10, 1, method = "fista_bt",  hasZIntercept = false, hasXIntercept = false, isVerbose = false);
smmr_Net2 = MatrixLMnet.lambda_min(est2);

# Lasso penalized regression - fista-bt cv
smmr_fistabt= Helium.readhe(joinpath(dataDir, "smmr_fistabt.he"))

println("CV Lasso vs Elastic Net when α=1 test 3 - fista-bt: ",
         @test smmr_Net1.AvgMSE == smmr_fistabt[:,1] && smmr_Net1.AvgPercentZero == smmr_fistabt[:,2] &&
               smmr_Net2.AvgMSE == smmr_fistabt[:,1] && smmr_Net2.AvgPercentZero == smmr_fistabt[:,2] )

############################################
# TEST 4 Lasso vs Elastic Net (𝛼=1) - admm #
############################################

# Elastic net penalized regression
MatrixLMnet.Random.seed!(2021)
est1 = MatrixLMnet.mlmnet_cv(dat, λ, α, 10, 1, method = "admm", hasZIntercept = false, hasXIntercept = false, isVerbose = false);
smmr_Net1 = MatrixLMnet.lambda_min(est1);

# Elastic net penalized regression
MatrixLMnet.Random.seed!(2021)
est2 = MatrixLMnet.mlmnet_cv(dat, λ, 10, 1, method = "admm",  hasZIntercept = false, hasXIntercept = false, isVerbose = false);
smmr_Net2 = MatrixLMnet.lambda_min(est2);

# Lasso penalized regression - fista-bt cv
smmr_admm = Helium.readhe(joinpath(dataDir, "smmr_admm.he"))

println("CV Lasso vs Elastic Net when α=1 test 4 - admm: ",
         @test smmr_Net1.AvgMSE == smmr_admm[:,1] && smmr_Net1.AvgPercentZero == smmr_admm[:,2] &&
               smmr_Net2.AvgMSE == smmr_admm[:,1] && smmr_Net2.AvgPercentZero == smmr_admm[:,2] )

##########################################
# TEST 5 Lasso vs Elastic Net (𝛼=1) - cd #
##########################################

# Elastic net penalized regression
MatrixLMnet.Random.seed!(2021)
est1 = MatrixLMnet.mlmnet_cv(dat, λ, α, 10, 1, method = "cd", hasZIntercept = false, hasXIntercept = false, isVerbose = false);
smmr_Net1 = MatrixLMnet.lambda_min(est1);

# Elastic net penalized regression
MatrixLMnet.Random.seed!(2021)
est2 = MatrixLMnet.mlmnet_cv(dat, λ, 10, 1, method = "cd",  hasZIntercept = false, hasXIntercept = false, isVerbose = false);
smmr_Net2 = MatrixLMnet.lambda_min(est2);

# Lasso penalized regression - cd cv
smmr_cd = Helium.readhe(joinpath(dataDir, "smmr_cd.he"))

println("CV Lasso vs Elastic Net when α=1 test 5 - cd: ",
         @test smmr_Net1.AvgMSE == smmr_cd[:,1] && smmr_Net1.AvgPercentZero == smmr_cd[:,2] &&
               smmr_Net2.AvgMSE == smmr_cd[:,1] && smmr_Net2.AvgPercentZero == smmr_cd[:,2] )

println("Tests mlmnet_cv finished!")


text:
function fix(a)
    if a>0
        return floor(a)
    elseif a<0
        return ceil(a)
    else
        return a
    end
end


text:
using DrWatson
@quickactivate :Swalbe
using Plots, CUDA, DataFrames, JDF

# Fluid dynamics we need for the experiment
"""
    run_dropletpatterned()

Simulates an droplet on a patterned substrate
"""
function measure_dropletpatterned(
    sys::Swalbe.SysConst,
    device::String;
    radius=20,
    θ₀=1/6,
    center=(sys.Lx÷2, sys.Ly÷2),
    θₛ=ones(sys.Lx, sys.Ly),
    height_dump=1000,
    df_drop = zeros(sys.Tmax÷height_dump, sys.Lx*sys.Ly),
    area_lv = zeros(sys.Tmax),
    area_ls = zeros(sys.Tmax),
    red_e = zeros(sys.Tmax),
    h_evo = zeros(sys.Tmax),
    verbos=true,
    T=Float64
)
    println("Simulating a droplet on a patterned substrate")
    fout, ftemp, feq, height, velx, vely, vsq, pressure, dgrad, Fx, Fy, slipx, slipy, h∇px, h∇py = Swalbe.Sys(sys, device, false, T)
    if device == "CPU"
        Swalbe.singledroplet(height, radius, θ₀, center)
    elseif device == "GPU"
        h = zeros(size(height))
        Swalbe.singledroplet(h, radius, θ₀, center)
        height = CUDA.adapt(CuArray, h)
    end
    Swalbe.equilibrium!(fout, height, velx, vely, vsq)
    ftemp .= fout
    for t in 1:sys.Tmax
        if t % sys.tdump == 0
            mass = 0.0
            mass = sum(height)
            if verbos
                println("Time step $t mass is $(round(mass, digits=3))")
            end
        end
        Swalbe.filmpressure!(pressure, height, dgrad, sys.γ, θₛ, sys.n, sys.m, sys.hmin, sys.hcrit)
        Swalbe.∇f!(h∇px, h∇py, pressure, dgrad, height)
        Swalbe.slippage!(slipx, slipy, height, velx, vely, sys.δ, sys.μ)
        # Forces are the pressure gradient and the slippage due to substrate liquid boundary conditions
        Fx .= h∇px .+ slipx
        Fy .= h∇py .+ slipy
        Swalbe.equilibrium!(feq, height, velx, vely, vsq)
        Swalbe.BGKandStream!(fout, feq, ftemp, -Fx, -Fy)
        Swalbe.moments!(height, velx, vely, fout)
        Swalbe.wetted!(area_ls, h_evo, height, t)
        # Far from clean but good enough for now, one should use dedicated arrays in principle.
        Swalbe.surfacearea!(area_lv, red_e, height, θₛ, h∇px, h∇py, dgrad, pressure, t)
        Swalbe.snapshot!(df_drop, height, t, dumping=height_dump)
    end
    return height, area_lv, area_ls, red_e, h_evo, df_drop
end

println("Drop relaxation experiments on patterned substrates")
# Different substrate patches
for shape in ["circle"] # "box" "triangle" "ellipse"
    # Different initial volumes
    for R in [80 85 90]
        # Different contact angle mismatch
        for contrast in [1/36 1/18 1/12 1/9 -1/18]
            println("Simulating shape $shape R0 $R contrast $(Int(round(rad2deg(contrast*pi), digits=0)))")
            sys = Swalbe.SysConst(Lx=512, Ly=512, γ=0.01, δ=0.5, n=3, m=2, hmin=0.07, Tmax=1000000)
            df_measures = DataFrame()
            df_drop = DataFrame()
            θₚ = ones(sys.Lx,sys.Ly)
            height_dump = 1000
            boxside = 87
            triside = 132
            ella = 73
            ellb = 33
            rcirc = 49
            # Substrate patterning
            if shape == "box"
                theta, P = Swalbe.boxpattern(θₚ, 1/9, center=(sys.Lx÷2, sys.Ly÷2), δₐ=-contrast, side=boxside)
            elseif shape == "triangle"
                theta, P = Swalbe.trianglepattern(θₚ, 1/9, center=(sys.Lx÷2, sys.Ly÷2), δₐ=-contrast, side=triside)
            elseif shape == "ellipse"
                theta, P = Swalbe.ellipsepattern(θₚ, 1/9, center=(sys.Lx÷2, sys.Ly÷2), δₐ=-contrast, a=ella, b=ellb)
            elseif shape == "circle"
                theta, P = Swalbe.ellipsepattern(θₚ, 1/9, center=(sys.Lx÷2, sys.Ly÷2), δₐ=-contrast, a=rcirc, b=rcirc)
            end
            # Make a cuarray with the substrate pattern
            θ_in = CUDA.adapt(CuArray, theta)
            # Actual simulation
            height, area_lv, area_ls, red_e, h_evo, drop_position = measure_dropletpatterned(sys, "GPU", radius=R, θₛ=θ_in)
            println("Writing measurements to dataframes")
            # Filling the dataframes
            df_measures[!, Symbol("A_lv")] = area_lv
            df_measures[!, Symbol("A_sl")] = area_ls
            df_measures[!, Symbol("E_red")] = red_e
            df_measures[!, Symbol("H_max")] = h_evo
            for t in 1:sys.Tmax÷height_dump
                df_drop[!, Symbol("time_step_$(t*height_dump)")] = drop_position[t,:]
            end
            println("Saving dataframe shape $shape R0 $R contrast $(Int(round(rad2deg(contrast*pi), digits=0))) to disk")
            file1 = JDF.save("data/Pattern_substrate/measures_shape_$(shape)_R0_$(R)_contrast_$(Int(round(rad2deg(contrast*pi), digits=0))).jdf", df_measures)
            file2 = JDF.save("data/Pattern_substrate/drop_pos_shape_$(shape)_R0_$(R)_contrast_$(Int(round(rad2deg(contrast*pi), digits=0))).jdf", df_drop)
            # Get back the memory of dataframes
            df_measures = DataFrame()
            df_drop = DataFrame()
            CUDA.reclaim()
        end
    end
end
println("Script done, let's have a look at the data :)")


text:
function TestCart()
nx=4;
nz=4;
xL=0;
xU=10;
zL=0;
zU=10;
Boundary.BT='FreeSlip';
Boundary.WE='FreeSlip';
Param.hS='AgnesiCart';
Param.hC=1;
Param.x0C=5;
Param.aC=1;

Param.Grid=CartGrid(nx,nz,xU-xL,zU-zL,xL,zL,OrientFaceCart,Boundary,Param);

DG.OrdPoly=2;
DG.OrdPolyX=2;
DG.OrdPolyY=2;
[DG.w,DG.xw]=GaussLobattoQuad(DG.OrdPoly);
[DG.wX,DG.xwX]=GaussLobattoQuad(DG.OrdPolyX);
[DG.wY,DG.xwY]=GaussLobattoQuad(DG.OrdPolyX);
[DG.DW,DG.DS,DG.DV,DG.DVT,DG.B]=DerivativeMatrixSingle(DG);
[DG.DWX,DG.DSX]=DerivativeMatrixSingle(DG);
[DG.DWY,DG.DSY]=DerivativeMatrixSingle(DG);

Param.X=zeros(DG.OrdPoly+1,DG.OrdPoly+1,3,Param.Grid.NumFaces);
Param.dXdx=zeros(DG.OrdPoly+1,DG.OrdPoly+1,2,2,Param.Grid.NumFaces);
Param.J=zeros(DG.OrdPoly+1,DG.OrdPoly+1,Param.Grid.NumFaces);
for iF=1:Param.Grid.NumFaces
  [Param.X(:,:,:,iF),Param.J(:,:,iF),Param.dXdx(:,:,:,:,iF)]=...
    JacobiDG1(DG,Param.Grid.Faces(iF),Param.Grid,@Cart,Param);
end
Subs=1;
fig=1;
PlotFaceGrid(Param.Grid,Subs,@JacobiCart,fig)
end


text:
using GeometricIntegrators
using ChargedParticleDynamics.GuidingCenter4d.SymmetricLoop
using GeometricExamples

const Δt     = 10.
const ntime  = 5000
const nloop  = 2000
const nsave  = 5
const nplot  = 5
const run_id = "poincare_1st_symmetric_vprk_dt10"

include("../guiding_center_4d_settings_vprk.jl")

tableau_list = get_tableau_list_vprk_projection()

pinv = guiding_center_4d_iode_poincare_invariant_1st(Δt, nloop, ntime, nsave)

include("guiding_center_4d_symmetric_poincare_invariant_1st.jl")


text:
export linear,∂linear

function linear(rᵢ)
  rₒ=rᵢ
end
function ∂linear(rᵢ)
  rₒ=ones(Float32,size(rᵢ))
end


text:
include("utility.jl")

base = "ATGC"
prob1 = [0.5, 0.25, 0.125, 0.125]
prob2 = [0.25, 0.25, 0.25, 0.25]
prob3 = [0.125, 0.125, 0.25, 0.5]

function generate_modelA(filename, length, base, prob)
    for i = 2:size(prob,1)
        prob[i] += prob[i-1]
    end
    # print(prob)

    f = open("data/$filename", "w")
    write(f, "$length\n")
    for i = 1:length
        cur = rand()
        # println(cur)
        j = 1
        while cur > prob[j]
            j += 1
        end
        write(f, "$(base[j])")
    end
end

function generate_modelB(filename, length, base)

end

generate_modelA("modelC.txt", 1000, base, prob2)

println("loss for modelA with compressor 1: ", kl_sdistance(prob2, prob1))
println("loss for modelA with compressor 1: ", kl_sdistance(prob2, prob2))
println("loss for modelA with compressor 1: ", kl_sdistance(prob2, prob3))

max_stars_repo_name: uchiiii/CompactDataStructure

max_issues_repo_name: uchiiii/CompactDataStructure

max_forks_repo_name: uchiiii/CompactDataStructure


text:
######################################################################
# Overwrite REQUIRE using dependency information from Project.toml.
#
# Call from the root of the package repository.
#
# Has some basic sanity checks, but **use at your own risk**, `REQUIRE`
# will be overwritten.
#
# The purpose of this script is to appease attobot, until
# https://github.com/attobot/attobot/issues/50 is fixed.
######################################################################
@assert VERSION ≥ v"0.7"
import Pkg
const PT = Pkg.Types
Pkg.activate(pwd())             # current directory as the project
ctx = PT.Context()
pkg = ctx.env.pkg
if pkg ≡ nothing
@error "Not in a package, I won't generate REQUIRE."
exit(1)
else
@info "found package" pkg = pkg
end
deps = PT.get_deps(ctx)
non_std_deps = sort(collect(setdiff(keys(deps), values(ctx.stdlibs))))
open("REQUIRE", "w") do io
println(io, "julia 0.7")
for d in non_std_deps
println(io, d)
@info "listing $d"
end
end

max_stars_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_issues_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_forks_repo_name: JuliaTomo/DynamicTomographyWithFlow


text:
@with_kw struct SimpleLightDark <: POMDPs.POMDP{Int,Int,Float64}
    discount::Float64       = 0.95
    correct_r::Float64      = 100.0
    incorrect_r::Float64    = -100.0
    light_loc::Int          = 10
    radius::Int             = 60
    big_step::Int           = 3
    min_obs_noise::Float64  = 0.5
end
discount(p::SimpleLightDark) = p.discount
isterminal_act(p::SimpleLightDark, a::Int) = a == 0

actions(p::SimpleLightDark) = [-p.big_step, -1, 0, 1, p.big_step]
n_actions(p::SimpleLightDark) = length(actions(p))
actionindex(p::SimpleLightDark, a::Int) = findfirst(actions(p) .== a)

states(p::SimpleLightDark) = -p.radius:p.radius
n_states(p::SimpleLightDark) = length(states(p))
stateindex(p::SimpleLightDark, s::Int) = s+p.radius+1

function transition(p::SimpleLightDark, s::Int, a::Int)
    return SparseCat(SVector(clamp(s+a, -p.radius, p.radius)), SVector(1.0))
end

function observation(p::SimpleLightDark, sp)
    Normal(sp, abs(sp - p.light_loc) + p.min_obs_noise)
end

function reward(p::SimpleLightDark, s, a)
    if a == 0
        return s == 0 ? p.correct_r : p.incorrect_r
    else
        return -1.0
    end
end

function initialstate_distribution(p::SimpleLightDark)
    ps = ones(2*div(p.radius,2)+1)
    ps /= length(ps)
    return SparseCat(div(-p.radius,2):div(p.radius,2), ps)
end

function max_possible_weight(pomdp::SimpleLightDark, a, o)
    return pdf(observation(pomdp, o), o)
end

function new_particle(p::SimpleLightDark, b, a, o, rng)
    return clamp(round(Int, rand(rng, observation(p, o))), -p.radius, p.radius)
end

max_entropy_distribution(pomdp::SimpleLightDark) = DiscreteUniform(-pomdp.radius, pomdp.radius)

@with_kw struct DSimpleLightDark <: POMDPs.POMDP{Int, Int, Int}
    sld::SimpleLightDark = SimpleLightDark()
    binsize::Float64     = 1.0
end

generate_o(p::DSimpleLightDark, sp, rng::AbstractRNG) = floor(Int, rand(rng, observation(p.sld, sp))/p.binsize)
generate_o(p::DSimpleLightDark, a, sp, rng::AbstractRNG) = generate_o(p, sp, rng)
function generate_sor(p::DSimpleLightDark, s, a, rng::AbstractRNG)
    sp = generate_s(p, s, a, rng)
    o = generate_o(p, sp, rng)
    r = reward(p, s, a, sp)
    return sp, o, r
end

function ParticleFilters.obs_weight(p::DSimpleLightDark, sp::Int, o::Int)
    cod = observation(p.sld, sp)
    return cdf(cod, (o+1)*p.binsize) - cdf(cod, o*p.binsize)
end

discount(p::DSimpleLightDark) = discount(p.sld)
isterminal_act(p::DSimpleLightDark, a::Int) = isterminal_act(p.sld, a)
isterminal(p::DSimpleLightDark, s::Number) = isterminal(p.sld, s)
actions(p::DSimpleLightDark) = actions(p.sld)
n_actions(p::DSimpleLightDark) = n_actions(p.sld)
actionindex(p::DSimpleLightDark, a::Int) = actionindex(p.sld, a)
states(p::DSimpleLightDark) = states(p.sld)
n_states(p::DSimpleLightDark) = n_states(p.sld)
stateindex(p::DSimpleLightDark, s::Int) = stateindex(p.sld, s)
transition(p::DSimpleLightDark, s::Int, a::Int) = transition(p.sld, s, a)
reward(p::DSimpleLightDark, s, a) = reward(p.sld, s, a)
initialstate_distribution(p::DSimpleLightDark) = initialstate_distribution(p.sld)

struct LDHeuristic{LD} <: Policy
    p::LD
    q::AlphaVectorPolicy{LD, Int}
    std_thresh::Float64
end

struct LDHSolver <: Solver
    q::QMDPSolver
    std_thresh::Float64
end

LDHSolver(;std_thresh::Float64=0.1, kwargs...) = LDHSolver(QMDPSolver(;kwargs...), std_thresh)

solve(sol::LDHSolver, pomdp::Union{SimpleLightDark}) = LDHeuristic(pomdp, solve(sol.q, pomdp), sol.std_thresh)
solve(sol::LDHSolver, pomdp::DSimpleLightDark) = solve(sol, pomdp.sld)

action(p::LDHeuristic, s::Int) = action(p.q, ParticleCollection([s]))
Random.seed!(p::LDHeuristic, s) = p

function action(p::LDHeuristic, b::AbstractParticleBelief)
    big_step = p.p.big_step
    s = std(b)
    if s <= p.std_thresh
        return action(p.q, b)
    else
        m = mean(particles(b))
        ll = p.p.light_loc
        if m == ll
            return -1*Int(sign(ll))
        elseif abs(m-ll) >= big_step
            return -Int(big_step*sign(m-ll))
        else
            return -Int(sign(m-ll))
        end
    end
end

struct LDSide <: Solver end

mutable struct LDSidePolicy{LD} <: Policy
    q::AlphaVectorPolicy{LD, Int}
end

solve(solver::LDSide, pomdp::Union{SimpleLightDark,DSimpleLightDark}) = LDSidePolicy(solve(QMDPSolver(), pomdp))
Random.seed!(p::LDSidePolicy, s) = p

function action(p::LDSidePolicy, b)
    big_step = p.q.pomdp.big_step
    if pdf(b, mode(b)) > 0.9
        return action(p.q, b)
    else
        return big_step
    end
end


text:
# Differentially Expressed Gene Analysis
# ======================================

"""
A list of probabilities of genes.

Fields
------

- `names`: gene names.
- `means`: estimated means (#genes x #queries).
- `negatives`: logarithmic probabilities (the base is 10; #genes × #queries).
- `positives`: logarithmic probabilities (the base is 10; #genes × #queries).
"""
struct DEGenes
    names::Vector{String}
    means::Matrix{Float64}
    negatives::Matrix{Float64}
    positives::Matrix{Float64}
end

const DEFAULT_NUM_NEIGHBORS = 10

"""
    finddegs(
        counts::AbstractMatrix,
        featurenames::AbstractVector{<:AbstractString},
        cellindexes::AbstractVector{<:Integer},
        index::CellIndex;
        k::Integer=$(DEFAULT_NUM_NEIGHBORS),
    ) -> DEGenes

Calculate probabilities of observing the counts or extremes.

This function returns a `DEGenes` object that contains probabilities observing
the counts stored in `counts` or their extremes. The expected expression level
for each gene is estimated from the nearby cells of `cellindexes` in the
database.  That is, the lower the probability is, the infrequent observing the
count by chance is.

A `DEGenes` object has two fields to store the probabilities: `negatives` and
`positives`.  If a value of `negatives` is low, that means that the
corresponding count is supposed to be negatively regulated. Likewise, if a
value of `positives` is low, that means that the corresponding count is
supposed to be positively regulated. Note that the probabilities are
log-transformed (the base of logarithm is 10) because they may possibly be very
close to zero. For example, a probability of 0.0001 (= 0.01%) becomes -4 after
the log transformation.

Arguments
---------

Required arguments (positional):
- `counts`:
    A transcriptome expression matrix (#features × #cells).
- `featurenames`:
    A vector of feature names (#features).
    Each of the elements must correspond to a row of the `counts` matrix.
- `cellindexes`:
    A vector of cell indexes in the database to be compared with the `counts`
    matrix (#cells).
- `index`:
    A database object.
    It must store raw counts of the database cells in it, which is activated by
    passing `keep_counts=true` when creating the database object (see
    [`CellIndex`](@ref)).

Other parameters:
- `k`: The number of nearest neighbors used to estimate the mean expression level.
"""
function finddegs(
        counts::AbstractMatrix,
        featurenames::AbstractVector{<:AbstractString},
        cellindexes::AbstractVector{<:Integer},
        index::CellIndex;
        k::Integer=DEFAULT_NUM_NEIGHBORS,)
    return finddegs(ExpressionMatrix(counts, featurenames), cellindexes, index, k=k)
end

function finddegs(
        Y::ExpressionMatrix,
        cellindexes::AbstractVector{<:Integer},
        index::CellIndex;
        k::Integer=DEFAULT_NUM_NEIGHBORS,)
    n = size(Y, 2)
    if n != length(cellindexes)
        throw(ArgumentError("mismatching size"))
    elseif index.counts === nothing
        throw(ArgumentError("no counts data"))
    elseif k < 0
        throw(ArgumentError("negative k"))
    end
    neighbors = findneighbors(k, cellindexes, index)
    featurenames = index.preproc.featurenames
    m = length(featurenames)
    # an extremely weak prior (Gamma distribution)
    α₀, β₀ = 1e-1, 1e-3  # mean=100.0, var=100000.0
    #α₀, β₀ = 1.0, 1e-2  # mean=100.0, var=10000.0
    means = zeros(m, n)
    negatives = zeros(m, n)
    positives = zeros(m, n)
    for j in 1:n
        counts_j = Y[featurenames,j]
        counts_nns = index.counts[:,neighbors.indexes[:,j]]
        counts_nns_normalized = sum(counts_j) * (counts_nns ./ sum(counts_nns, dims=1))
        α = α₀ .+ sum(counts_nns_normalized, dims=2)
        β = β₀  + k
        for i in 1:m
            nb = NegativeBinomial(α[i], 1-inv(β+1))
            #@assert mean(nb) ≈ α[i] / β
            y = counts_j[i]
            means[i,j] = α[i] / β
            negatives[i,j] = logcdf(nb, y)
            positives[i,j] = logccdf(nb, y-1)
        end
    end
    # change the base of logarithm from e (≈2.718) to 10
    negatives .*= inv(log(10))
    positives .*= inv(log(10))
    return DEGenes(copy(featurenames), means, negatives, positives)
end


text:
using BlockArrays, LinearAlgebra, Test

@testset "BlockArrays.jl" begin
    include("test_blockindices.jl")
    include("test_blockarrays.jl")
    include("test_blockviews.jl")
    include("test_blockrange.jl")
    include("test_blockarrayinterface.jl")
    include("test_blockbroadcast.jl")
    include("test_blocklinalg.jl")
end


text:
using Random: MersenneTwister

m1 = fit(LinearModel, form, data)

loss1 = RobustModels.L2Loss()
loss2 = RobustModels.TukeyLoss()
est1 = MEstimator(loss1)
est2 = MEstimator(loss2)

@testset "linear: M-estimator $(name)" for name in ("L1", "Huber", "L1L2", "Fair", "Logcosh", "Arctan", "Cauchy", "Geman", "Welsch", "Tukey", "YohaiZamar")
    typeloss = getproperty(RobustModels, Symbol(name * "Loss"))
    l = typeloss()
    est = MEstimator(l)
    est3 = MEstimator{typeloss}()
    m2 = fit(RobustLinearModel, form, data, est; method=:cg, verbose=false, initial_scale=:mad)
    m3 = fit(RobustLinearModel, form, data, est; method=:chol, initial_scale=:mad)
    m4 = rlm(form, data, est; method=:chol, initial_scale=:mad)
    m5 = rlm(form, data, est3; method=:chol, initial_scale=:mad)
    println("\n\t\u25CF Estimator: $(name)")
    println(" lm       : ", coef(m1))
    println("rlm(cg)   : ", coef(m2))
    println("rlm(chol) : ", coef(m3))
    println("rlm2(chol): ", coef(m4))
    println("rlm3(chol): ", coef(m5))
    println(m2)
    if name != "L1"
        @test isapprox(coef(m2), coef(m3); rtol=1e-5)
    end

    # refit
    β2 = copy(coef(m2))
    refit!(m2, y; wts=weights(m2), verbose=false, initial_scale=:mad)
    @test all(coef(m2) .== β2)
end

@testset "linear: S-estimator" begin
    println("\n\t\u25CF Estimator type: S-estimators")
    println(" lm             : ", coef(m1))
    @testset "Not bounded: $(name)" for name in ("Huber", "L1L2", "Fair", "Logcosh", "Arctan", "Cauchy")
        typeloss = getproperty(RobustModels, Symbol(name * "Loss"))
        @test_throws TypeError SEstimator{typeloss}()
    end

    @testset "Bounded: $(name)" for name in ("Geman", "Welsch", "Tukey", "YohaiZamar")
        typeloss = getproperty(RobustModels, Symbol(name * "Loss"))
        est = SEstimator{typeloss}()
        m = fit(RobustLinearModel, form, data, est; method=:cg, verbose=false, initial_scale=:mad)
        println("rlm($name) : ", coef(m))
        ## TODO: find better test
        @test size(coef(m), 1) == 2

        # Resampling
        rng = MersenneTwister(1)
        opts = Dict(:Nsteps_β=>3, :rng=>rng)
        m2 = fit(RobustLinearModel, form, data, est; method=:cg, verbose=false, initial_scale=:mad, resample=true, resampling_options=opts)
        @test (scale(m2) / scale(m) - 1) <= 1e-2
    end
end

@testset "linear: MM-estimator" begin
    println("\n\t\u25CF Estimator type: MM-estimators")
    println(" lm             : ", coef(m1))
    @testset "Not bounded: $(name)" for name in ("Huber", "L1L2", "Fair", "Logcosh", "Arctan", "Cauchy")
        typeloss = getproperty(RobustModels, Symbol(name * "Loss"))
        @test_throws TypeError MMEstimator{typeloss}()
    end

    @testset "Bounded: $(name)" for name in ("Geman", "Welsch", "Tukey", "YohaiZamar")
        typeloss = getproperty(RobustModels, Symbol(name * "Loss"))
        est = MMEstimator{typeloss}()
        m = fit(RobustLinearModel, form, data, est; method=:cg, verbose=false, initial_scale=:mad)
        println("rlm($name) : ", coef(m))
        ## TODO: find better test
        @test size(coef(m), 1) == 2

        # Resampling
        rng = MersenneTwister(1)
        opts = Dict(:Npoints=>10, :Nsteps_β=>3, :Nsteps_σ=>3, :rng=>rng)
        m2 = fit(RobustLinearModel, form, data, est; method=:cg, verbose=false, initial_scale=:mad, resample=true, resampling_options=opts)
        @test (scale(m2) / scale(m) - 1) <= 1e-2
    end
end

@testset "linear: τ-estimator" begin
    println("\n\t\u25CF Estimator type: τ-estimators")
    println(" lm             : ", coef(m1))
    @testset "Not bounded: $(name)" for name in ("Huber", "L1L2", "Fair", "Logcosh", "Arctan", "Cauchy")
        typeloss = getproperty(RobustModels, Symbol(name * "Loss"))
        @test_throws TypeError TauEstimator{typeloss}()
    end

    @testset "Bounded: $(name)" for name in ("Geman", "Welsch", "Tukey", "YohaiZamar")
        typeloss = getproperty(RobustModels, Symbol(name * "Loss"))
        est = TauEstimator{typeloss}()

        m = fit(RobustLinearModel, form, data, est; method=:cg, verbose=false, initial_scale=:mad)
        println("rlm($name) : ", coef(m))
        ## TODO: find better test
        @test size(coef(m), 1) == 2

        # Resampling
        rng = MersenneTwister(1)
        opts = Dict(:Nsteps_β=>3, :Nsteps_σ=>3, :rng=>rng)
        m2 = fit(RobustLinearModel, form, data, est; method=:cg, verbose=false, initial_scale=:mad, resample=true, resampling_options=opts)
        @test (tauscale(m2) / tauscale(m) - 1) <= 1e-2
    end
end

@testset "linear: leverage correction" begin
    m2 = fit(RobustLinearModel, form, data, est2; method=:cg, initial_scale=:mad, correct_leverage=false)
    m3 = fit(RobustLinearModel, form, data, est2; method=:cg, initial_scale=:mad, correct_leverage=true)
    println("rlm(without leverage correction) : ", coef(m2))
    println("rlm(   with leverage correction) : ", coef(m3))
    ## TODO: find better test
    @test isapprox(coef(m2), coef(m3); rtol=0.1)
end


text:
using Documenter, Pluto,PlutoVista

function rendernotebook(name)
    input=joinpath(@__DIR__,"..","examples",name*".jl")
    output=joinpath(@__DIR__,"src",name*".html")
    session = Pluto.ServerSession();
    notebook = Pluto.SessionActions.open(session, input; run_async=false)
    html_contents = Pluto.generate_html(notebook)
    write(output, html_contents)
end

function mkdocs()

    rendernotebook("plutovista")

    makedocs(sitename="PlutoVista.jl",
             modules = [PlutoVista],
             doctest = false,
             clean = false,
             authors = "J. Fuhrmann",
             repo="https://github.com/j-fu/PlutoVista.jl",
             pages=[
                 "Home"=>"index.md",
                 "plotly.js API"=> "plotly.md",
                 "vtk.js API"=> "vtk.md",
                 "Internals" => "internals.md"
             ])
end

mkdocs()

deploydocs(repo = "github.com/j-fu/PlutoVista.jl.git",devbranch = "main")


text:
# ft_getvidpid.jl
export ft_getvidpid

function ft_getvidpid()
  vid = Ref{Cuint}()
  pid = Ref{Cuint}()
  ft_status = ccall((:FT_GetVIDPID, d2xx),
                     Cuint,
                     (Ref{Cuint}, Ref{Cuint}),
                      vid, pid)
  checkstatus(ft_status)
  return vid[], pid[]
end


text:
module SalsaMLIncr

using Salsa

const N_FEATURES = 1
const Sample = NTuple{N_FEATURES, Number}
const learning_rate = 0.0005

@declare_input num_samples(s)::Int
@declare_input sample(s, i::Int)::Sample
@declare_input response(s, i::Int)::Float64

init_lr!(s::Runtime) = set_num_samples!(s, 0)
function new_lr_runtime()
    s = Runtime()
    init_lr!(s)
    return s
end

function insert_training_pair!(s, sample::Sample, response::Float64)
    i = num_samples(s) + 1
    set_num_samples!(s, i)
    set_sample!(s, i, sample)
    set_response!(s, i, response)
end

function predict(s::Runtime, x::Sample)::Float64
    sum(x .* lr_learned_weights(s))
end

# Derived function with no inputs (essentially a constant)
@derived init_weights(s::Runtime) = Tuple(rand(Float64, N_FEATURES))

@derived function lr_learned_weights(s::Runtime)::NTuple{N_FEATURES,Float64}
    if num_samples(s) == 0
        error("Cannot learn a linear regression with no samples.")
    end
    weights = lr_learned_weights_unrolled(s, 0)
    for i in 1:100  # Stop after 100 iterations
        new_weights = lr_learned_weights_unrolled(s, i)
        if abs(sum(new_weights .- weights)) < 0.0001
            @info "stopped after $i iterations"
            return weights
        end
        weights = new_weights
    end
    @info "timed out in learn iterations"
    return weights
end

@derived function lr_learned_weights_unrolled(s::Runtime, iteration::Int)

    if iteration == 0
        return init_weights(s)
    else
        n = num_samples(s)
        weights = lr_learned_weights_unrolled(s, iteration-1)
        deriv,_ = lr_δweights_δmse(s, iteration-1)
        δweights_δcost = deriv .* learning_rate
        weights = weights .- δweights_δcost
        return weights
    end
end

# HACK: Return the actual value, AND `n`, so that we can use Salsa.previous_output() to
# incrementally maintain the result! :D
@derived function lr_δweights_δmse(s::Runtime, iteration::Int)
    prev_value = Salsa.previous_output(s)
    prev_deriv = 0
    start_n = 1
    if prev_value !== nothing  # Incremental case
        (prev_deriv, start_n) = prev_value
    end

    N = num_samples(s)
    n = length(1:N)
    new_deriv = (-2/n) * sum(
        sum(sample(s, i) .* (response(s, i) - lr_predicted_unrolled(s, i, iteration)))
        for i in start_n:N
    )
    return prev_deriv + new_deriv, n
end
@derived function lr_predicted_unrolled(s::Runtime, i::Int, iteration::Int)
    sum(sample(s, i) .* lr_learned_weights_unrolled(s, iteration))
end

# For logging purposes
@derived function lr_mse(s::Runtime)
    n = num_samples(s)
    1/n * sum(
        (response(s, i) - lr_predicted(s, i)) ^ 2
        for i in 1:num_samples(s)
    )
end

end


text:
"""
Operator

Construction of `Operator` object.
"""

export operator, trans_inv_operator
"""
    Operator{Tv}

Object that store the local operator.

Properties:
-----------
- M : Vector{SparseMatrixCSC{Tv, Int}}, local operators represented by CSC sparse matrices.
- I : Vector{Vector{Int}}, indices of sites that the operators act on.
- B : Basis.
"""
struct Operator{Tv<:Number, Tb<:AbstractBasis}
    M::Vector{SparseMatrixCSC{Tv, Int}}
    I::Vector{Vector{Int}}
    B::Tb
end

@inline eltype(opt::Operator{Tv, Tb}) where Tv where Tb = promote_type(Tv, eltype(opt.B))
@inline length(opt::Operator) = length(opt.M)
@inline size(opt::Operator) = size(opt.B)
@inline size(opt::Operator, i::Integer) = size(opt.B, i)
function Base.display(opt::Operator)
    println("Operator of size $(size(opt)) with $(length(opt)) terms.")
end

function operator(mats::AbstractVector{<:AbstractMatrix}, inds::AbstractVector{<:AbstractVector}, B::AbstractBasis)
    num = length(mats)
    @assert num == length(inds) "Numbers mismatch: $num matrices and $(length(inds)) indices."
    dtype = promote_type(eltype.(mats)...)
    M = Vector{SparseMatrixCSC{dtype, Int64}}(undef, num)
    I = Vector{Vector{Int64}}(undef, num)
    N = 0
    for i = 1:num
        iszero(mats[i]) && continue
        ind = inds[i]
        pos = findfirst(x -> isequal(x, ind), view(I, 1:N))
        if isnothing(pos)
            N += 1
            I[N] = ind
            M[N] = sparse(mats[i])
        else
            M[pos] += mats[i]
        end
    end
    deleteat!(M, N+1:num)
    deleteat!(I, N+1:num)
    Operator(M, I, B)
end

function operator(mats::AbstractVector{<:AbstractMatrix}, inds::AbstractVector{<:AbstractVector}, L::Integer)
    B = TensorBasis(L, base=find_base(size(mats[1], 1), length(inds[1])))
    operator(mats, inds, B)
end

operator(mats::AbstractVector{<:AbstractMatrix}, inds::AbstractVector{<:Integer}, C) = operator(mats, [[i] for i in inds], C)
operator(mat::AbstractMatrix, ind::AbstractVector{<:Integer}, C) = operator([mat], [ind], C)
operator(mats::AbstractVector{<:AbstractMatrix}, L::Integer) = operator(mats, [[i] for i in 1:L], L)
operator(mats::AbstractVector{<:AbstractMatrix}, B::AbstractBasis) = operator(mats, [[i] for i in 1:length(B.dgt)], B)

function trans_inv_operator(mat::AbstractMatrix, ind::AbstractVector{<:Integer}, B::AbstractBasis)
    L = length(B.dgt)
    smat = sparse(mat)
    mats = fill(smat, L)
    inds = [mod.(ind .+ i, L) .+ 1 for i = -1:L-2]
    operator(mats, inds, B)
end

function trans_inv_operator(mat::AbstractMatrix, ind::AbstractVector{<:Integer}, L::Integer)
    B = TensorBasis(L, base=find_base(size(mat, 1), length(ind)))
    trans_inv_operator(mat, ind, B)
end

trans_inv_operator(mat::AbstractMatrix, M::Integer, C) = trans_inv_operator(mat, 1:M, C)

*(c::Number, opt::Operator) = operator(c .* opt.M, opt.I, opt.B)
*(opt::Operator, c::Number) = c * opt
/(opt::Operator, c::Number) = operator(opt.M ./ c, opt.I, opt.B)
function +(opt1::Operator, opt2::Operator)
    mats = vcat(opt1.M, opt2.M)
    inds = vcat(opt1.I, opt2.I)
    operator(mats, inds, opt1.B)
end
+(opt1::Operator, ::Nothing) = opt1
+(::Nothing, opt1::Operator) = opt1
-(opt::Operator) = Operator(-opt.M, opt.I, opt.B)
-(opt1::Operator, opt2::Operator) = opt1 + (-opt2)

@inline function find_base(a::Integer, b::Integer)
    if b == 1
        return a
    else
        for i=2:a
            if i^b == a
                return i
            end
        end
        error("Incompatible dimension: ($a, $b)")
    end
end

# Operator to matrices
export addto!
function addto!(M::AbstractMatrix, opt::Operator)
    for j = 1:size(opt.B, 2)
        colmn!(view(M, :, j), opt, j)
    end
    M
end

function Array(opt::Operator)
    M = zeros(eltype(opt), size(opt))
    if size(M, 1) > 0 && size(M, 2) > 0
        addto!(M, opt)
    end
    M
end

function sparse(opt::Operator)
    M = spzeros(eltype(opt), size(opt)...)
    if size(M, 1) > 0 && size(M, 2) > 0
        addto!(M, opt)
    end
    M
end

function mul!(target::AbstractVector, opt::Operator, v::AbstractVector)
    for j = 1:length(v)
        colmn!(target, opt, j, v[j])
    end
    target
end

function mul!(target::AbstractMatrix, opt::Operator, m::AbstractMatrix)
    for j = 1:size(m, 1)
        colmn!(target, opt, j, view(m, j, :))
    end
    target
end

function *(opt::Operator, v::AbstractVector)
    ctype = promote_type(eltype(opt), eltype(v))
    M = zeros(ctype, size(opt, 1))
    mul!(M, opt, v)
end

function *(opt::Operator, m::AbstractMatrix)
    ctype = promote_type(eltype(opt), eltype(m))
    M = zeros(ctype, size(opt, 1), size(m, 2))
    mul!(M, opt, m)
end

# Helper functions
function colmn!(target::AbstractVector, M::SparseMatrixCSC, I::Vector{Int}, b::AbstractBasis, coeff::Number=1)
    rows, vals = rowvals(M), nonzeros(M)
    j = index(b.dgt, I, base=b.B)
    change = false
    for i in nzrange(M, j)
        row, val = rows[i], vals[i]
        change!(b.dgt, I, row, base=b.B)
        C, pos = index(b)
        target[pos] += coeff * C * val
        change = true
    end
    change ? change!(b.dgt, I, j, base=b.B) : nothing
end

function colmn!(target::AbstractVector, opt::Operator, j::Integer, coeff::Number=1)
    b, M, I = opt.B, opt.M, opt.I
    C = coeff / change!(b, j)
    for i = 1:length(M)
        colmn!(target, M[i], I[i], b, C)
    end
end

function colmn!(target::AbstractMatrix, M::SparseMatrixCSC, I::Vector{Int}, b::AbstractBasis, coeff::AbstractVector{<:Number})
    rows, vals = rowvals(M), nonzeros(M)
    j = index(b.dgt, I, base=b.B)
    change = false
    for i in nzrange(M, j)
        row, val = rows[i], vals[i]
        change!(b.dgt, I, row, base=b.B)
        C, pos = index(b)
        target[pos, :] .+= (C * val) .* coeff
        change = true
    end
    change ? change!(b.dgt, I, j, base=b.B) : nothing
end

function colmn!(target::AbstractMatrix, opt::Operator, j::Integer, coeff::AbstractVector{<:Number})
    b, M, I = opt.B, opt.M, opt.I
    N = change!(b, j) .* coeff
    for i = 1:length(M)
        colmn!(target, M[i], I[i], b, N)
    end
end

# Spin Matrices
export spin

@inline spin_coeff(D::Integer) = [sqrt(i*(D-i)) for i = 1:D-1]
@inline spin_Sp(D::Integer) = sparse(1:D-1, 2:D, spin_coeff(D), D, D)
@inline spin_Sm(D::Integer) = sparse(2:D, 1:D-1, spin_coeff(D), D, D)
@inline function spin_Sx(D::Integer)
    coeff = spin_coeff(D) / 2
    sp = sparse(1:D-1, 2:D, coeff, D, D)
    sp + sp'
end
@inline function spin_iSy(D::Integer)
    coeff = spin_coeff(D) / 2
    sp = sparse(1:D-1, 2:D, coeff, D, D)
    sp - sp'
end
@inline function spin_Sz(D::Integer)
    J = (D-1) / 2
    sparse(1:D, 1:D, J:-1:-J)
end

function spin_dict(c::Char, D::Integer)
    if     isequal(c, '+') spin_Sp(D)
    elseif isequal(c, '-') spin_Sm(D)
    elseif isequal(c, 'x') spin_Sx(D)
    elseif isequal(c, 'y') spin_iSy(D)
    elseif isequal(c, 'z') spin_Sz(D)
    elseif isequal(c, '1') spdiagm(ones(D))
    elseif isequal(D, 2)
        if     isequal(c, 'X') sparse([0 1; 1  0])
        elseif isequal(c, 'Y') sparse([0 1;-1  0])
        elseif isequal(c, 'Z') sparse([1 0; 0 -1])
        else error("Invalid spin symbol: $c.")
        end
    else error("Invalid spin symbol: $c.")
    end
end

"""
    spin(s::String; D::Integer=2)

Return matrix for spin operators.
"""
function spin(s::String; D::Integer=2)
    ny = if isequal(D, 2)
        sum(isequal(si, 'y')||isequal(si, 'Y') for si in s)
    else
        sum(isequal(si, 'y') for si in s)
    end
    mat = isone(length(s)) ? spin_dict(s[1], D) : kron([spin_dict(si, D) for si in s]...)
    sign = iszero(mod(ny, 2)) ? (-1)^(ny÷2) : (-1im)^ny
    sign * mat
end
spin(c::Number, s::String; D::Integer=2) = c * spin(s, D=D)

# Precompile the functions
precompile(spin_coeff, (Int64,))
precompile(spin_Sp, (Int64,))
precompile(spin_Sm, (Int64,))
precompile(spin_Sx, (Int64,))
precompile(spin_iSy, (Int64,))
precompile(spin_Sz, (Int64,))
precompile(spin_dict, (Char, Int64))
spin("x")
spin("xyz", D=3)

"""
    spin(spins::Tuple{<:Number, String}...; D::Integer=2)

Return matrix for spin operators.
The spins should be an iterable onject, each item is of the form (::Number, ::String).
"""
function spin(spins::Tuple{<:Number, String}...; D::Integer=2)
    sum(ci * spin(si, D=D) for (ci, si) in spins)
end

spin(spins::AbstractVector{<:Tuple{<:Number, String}}; D::Integer=2) = spin(spins..., D=D)

function operator(s::String, inds::AbstractVector{<:Integer}, basis::AbstractBasis)
    mat = spin(s, D=base(basis))
    operator(mat, inds, basis)
end

function operator(s::String, inds::AbstractVector{<:Integer}, L::Integer; base::Integer=2)
    basis = tensorbasis(L, base=base)
    operator(s, inds, basis)
end

function trans_inv_operator(s::String, inds::AbstractVector{<:Integer}, basis::AbstractBasis)
    mat = spin(s, D=base(basis))
    trans_inv_operator(mat, inds, basis)
end

function trans_inv_operator(s::String, basis::AbstractBasis)
    mat = spin(s, D=base(basis))
    trans_inv_operator(mat, length(s), basis)
end

function trans_inv_operator(s::String, L::Integer; base::Integer=2)
    basis = TensorBasis(L, base=base)
    trans_inv_operator(s, basis)
end


text:
### Observation model  ############################################################################

function POMDPs.observation(pomdp::SingleOCPOMDP, a::SingleOCAction, sp::SingleOCState)
    if !is_observable_fixed(sp, pomdp.env) || off_the_grid(pomdp, sp.ped)
        o = SingleOCObs(false, sp.ego, get_off_the_grid(pomdp))
        return Deterministic(o)
    elseif collision_checker(sp.ego, sp.ped, pomdp.ego_type, pomdp.ped_type)
        return Deterministic(sp)
    elseif isterminal(pomdp, sp)
        return Deterministic(sp)
    end
    ego = sp.ego
    ped = sp.ped

    neighbors = Vector{VehicleState}(undef, 9)
    neighbors[1] = yv_to_state(pomdp, ped.posG.y - pomdp.pos_res, ped.v)
    neighbors[2] = yv_to_state(pomdp, ped.posG.y + pomdp.pos_res, ped.v)
    neighbors[3] = yv_to_state(pomdp, ped.posG.y - pomdp.pos_res, ped.v - pomdp.vel_res)
    neighbors[4] = yv_to_state(pomdp, ped.posG.y + pomdp.pos_res, ped.v - pomdp.vel_res)
    neighbors[5] = yv_to_state(pomdp, ped.posG.y - pomdp.pos_res, ped.v + pomdp.vel_res)
    neighbors[6] = yv_to_state(pomdp, ped.posG.y + pomdp.pos_res, ped.v + pomdp.vel_res)
    neighbors[7] = yv_to_state(pomdp, ped.posG.y, ped.v - pomdp.vel_res)
    neighbors[8] = yv_to_state(pomdp, ped.posG.y, ped.v + pomdp.vel_res)
    neighbors[9] = yv_to_state(pomdp, ped.posG.y, ped.v)

    grid = RectangleGrid(get_Y_grid(pomdp), get_V_ped_grid(pomdp)) #XXX preallocate

    states = SingleOCObs[]
    sizehint!(states, 9)
    for neighbor in neighbors
        if in_bounds_ped(pomdp, neighbor) && !collision_checker(ego, neighbor, pomdp.ego_type, pomdp.ped_type)
            # make sure neighbor is in the grid
            neigh_yv_ind, neigh_yv_weights = interpolants(grid, [neighbor.posG.y, neighbor.v])
            neigh_yv = neigh_yv_ind[argmax(neigh_yv_weights)]
            neighbor = yv_to_state(pomdp, ind2x(grid, neigh_yv)...)
            push!(states, SingleOCObs(false, ego, neighbor))
        end
    end
    probs = zeros(length(states))
    for (i, o) in enumerate(states)
        probs[i] = obs_weight(o, sp, pomdp)
    end
    probs = normalize!(probs, 1)
    @assert length(probs) == length(states)
    return SparseCat(states, probs)
end

"""
    function obs_weight(o::SingleOCState, s::SingleOCState, pomdp::SingleOCPOMDP)
Given a continuous observation returns a weight proportional to the probability of observing o
while being in the state s
"""
function obs_weight(o::SingleOCState, s::SingleOCState, pomdp::SingleOCPOMDP)
    weight = 1.0
    if !is_observable_fixed(s.ego, s.ped, pomdp.env)
        if off_the_grid(pomdp, o.ped)
            return weight
        else
            return 0.
        end
    else
        if off_the_grid(pomdp, o.ped) || !is_observable_fixed(o.ego, o.ped, pomdp.env)
            return 0.
        else
            pos_noise = pomdp.pos_obs_noise
            vel_noise = pomdp.vel_obs_noise
            weight *= 1/(sqrt(2*pi*pos_noise*vel_noise))*exp(-1/pos_noise*(s.ped.posG.y - o.ped.posG.y)^2 - 1/vel_noise*(s.ped.v - o.ped.v)^2)
            return weight
        end
    end
end

"""
Check if a pedestrian in s is observable or not
"""
function is_observable_fixed(s::SingleOCState, env::CrosswalkEnv)
    m = length(env.obstacles)
    ped = s.ped
    ego = s.ego
    front = ego.posG + polar(VehicleDef().length/2, ego.posG.θ)
    angle = atan(ped.posG.y - front.y, ped.posG.x - front.x)
    ray = Projectile(VecSE2(front.x, front.y, angle), 1.0)
    if isinf(ped.v)
        return false
    end
    for i = 1:m
        if is_colliding(ray, env.obstacles[i], ped.posG)
            return false
        end
    end
    return true
end

"""
Check if a state is in bound
"""
function in_bounds_ped(pomdp::SingleOCPOMDP, veh::VehicleState)
    cl = pomdp.env.params.crosswalk_length
    y_min = -cl/4
    y_max = cl/4
    return (y_min <= veh.posG.y <= y_max) && (0. <= veh.v <= pomdp.env.params.ped_max_speed)
end

### Helpers for converting SingleOCObs to NN input

# """
#     extract_gray!(s::Array{UInt32, 2}, dest::Array{Float64, 2})
# Extract a gray scale image from a CairoSurface object
# """
# function extract_gray!(s::Array{UInt32, 2}, dest::Array{Float64, 2})
#     for (i,s) in enumerate(s)
#         dest[i] = 0.2989*convert(Float64, (s >> 16) & 0xFF) +
#                   0.5870*convert(Float64, (s >> 8) & 0xFF) +
#                   0.1140*convert(Float64, s & 0xFF)
#     end
#     return dest'
# end

# uncomment for image representation
# /!\ do not forget to change the observation space size in the python wrapper
# function POMDPs.convert_o(o::SingleOCObs, pomdp::SingleOCPOMDP)
#     scene = state_to_scene(pomdp, o)
#     img = AutoViz.render(scene, DQNFeatures(pomdp.env), [InflateOverlay()], cam=StaticCamera(VecE2(25, 0.), 20.))
#     gray_img = Array{Float64, 2}(Int(img.width), Int(img.height))
#     extract_gray!(img.data, gray_img)
#     gray_img = imresize(gray_img, (84,84))
#     gray_img /= 255.0
# end

# uncomment for engineered feature representation with SingleOCclusion map
# function POMDPs.convert_o(::Type{Array{Float64, 1}}, o::SingleOCObs, pomdp::SingleOCPOMDP)
#     X = get_X_grid(pomdp)
#     Y = get_Y_grid(pomdp)
#     im_w = length(X)
#     im_h = length(Y)
#     max_speed = get_V_grid(pomdp)[end]
#     features = zeros(im_w, im_h, 4)
#     for (i,x) in enumerate(X)
#         for (j,y) in enumerate(Y)
#             features[i,j,1] = is_colliding(x, y, o.ego) ? 1.0 : 0.0
#             features[i,j,2] = (is_colliding(x, y, o.ped) && is_observable(x, y, o.ego, pomdp.env)) ? 1.0 : 0.0
#             features[i,j,1] = is_observable(x, y, o.ego, pomdp.env) ? 1.0 : 0.0
#             if is_colliding(x, y, o.ego)
#                 features[i,j,2] = o.ego.v/max_speed
#             elseif is_colliding(x,  y, o.ped) && is_observable(x, y, o.ego, pomdp.env)
#                 features[i,j,2] = o.ped.v/max_speed
#             end
#         end
#     end
#     return features
# end

# uncomment for vector representation
# /!\ do not forget to change the observation space size in the python wrapper
function POMDPs.convert_o(::Type{Array{Float64, 1}}, o::SingleOCObs, pomdp::SingleOCPOMDP)
    o_vec = Vector{Float64}(4)
    o_vec[1], o_vec[2] = o.ego.posG.x/pomdp.x_goal, o.ego.v/pomdp.env.params.speed_limit
    o_vec[3] = o.ped.posG.y/pomdp.y_goal
    o_vec[4] = o.ped.v
    if isinf(o_vec[4])
        o_vec[4] = 0.
    end
    o_vec[4] /= 2.
    return reshape(o_vec, (4,1,1))
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
using Test, Dates, Random, UUIDs, Statistics
using OndaEDF, Onda, EDF, Tables
using FilePathsBase: AbstractPath, PosixPath

function test_edf_signal(rng, label, transducer, physical_units,
                         physical_min, physical_max,
                         digital_min, digital_max,
                         samples_per_record, n_records)
    header = EDF.SignalHeader(label, transducer, physical_units,
                              physical_min, physical_max,
                              digital_min, digital_max,
                              "", Int16(samples_per_record))
    samples = rand(rng, Int16, n_records * samples_per_record)
    return EDF.Signal(header, samples)
end

function make_test_data(rng, sample_rate, samples_per_record, n_records)
    imin16, imax16 = Float32(typemin(Int16)), Float32(typemax(Int16))
    anns_1 = [[EDF.TimestampedAnnotationList(i, nothing, []),
               EDF.TimestampedAnnotationList(i, i + 1, ["", "$i a", "$i b"])] for i in 1:n_records]
    anns_2 = [[EDF.TimestampedAnnotationList(i, nothing, []),
               EDF.TimestampedAnnotationList(i, 0, ["", "$i c", "$i d"])] for i in 1:n_records]
    edf_signals = Union{EDF.AnnotationsSignal,EDF.Signal}[
        test_edf_signal(rng, "EEG F3-M2", "E", "uV",  -32768.0f0, 32767.0f0,   imin16, imax16, samples_per_record, n_records)
        test_edf_signal(rng, "EEG F4-M1", "E", "uV",  -32768.0f0, 32767.0f0,   imin16, imax16, samples_per_record, n_records)
        test_edf_signal(rng, "EEG C3-M2", "E", "uV",  -32768.0f0, 32767.0f0,   imin16, imax16, samples_per_record, n_records)
        test_edf_signal(rng, "EEG O1-M2", "E", "uV",  -32768.0f0, 32767.0f0,   imin16, imax16, samples_per_record, n_records)
        test_edf_signal(rng, "C4-M1", "E", "uV",  -32768.0f0, 32767.0f0,   imin16, imax16, samples_per_record, n_records)
        test_edf_signal(rng, "O2-A1", "E", "uV",  -32768.0f0, 32767.0f0,   imin16, imax16, samples_per_record, n_records)
        test_edf_signal(rng, "E1", "E", "uV",  -32768.0f0, 32767.0f0,   imin16, imax16, samples_per_record, n_records)
        test_edf_signal(rng, "E2", "E", "uV",  -32768.0f0, 32767.0f0,   imin16, imax16, samples_per_record, n_records)
        test_edf_signal(rng, "Fpz",   "E", "uV",  -32768.0f0, 32767.0f0,   imin16, imax16, samples_per_record, n_records)
        EDF.AnnotationsSignal(samples_per_record, anns_1)
        test_edf_signal(rng, "EMG LAT",   "E", "uV",  -32768.0f0, 32767.0f0,   imin16, imax16, samples_per_record, n_records)
        test_edf_signal(rng, "EMG RAT",   "E", "uV",  -32768.0f0, 32767.0f0,   imin16, imax16, samples_per_record, n_records)
        test_edf_signal(rng, "SNORE", "E", "uV",  -32768.0f0, 32767.0f0,   imin16, imax16, samples_per_record, n_records)
        test_edf_signal(rng, "IPAP", "", "cmH2O", -74.0465f0, 74.19587f0,  imin16, imax16, samples_per_record, n_records)
        test_edf_signal(rng, "EPAP", "", "cmH2O", -73.5019f0, 74.01962f0,  imin16, imax16, samples_per_record, n_records)
        test_edf_signal(rng, "CFLOW", "",  "LPM", -309.153f0, 308.8513f0,  imin16, imax16, samples_per_record, n_records)
        test_edf_signal(rng, "PTAF", "",     "v", -125.009f0, 125.009f0,   imin16, imax16, samples_per_record, n_records)
        test_edf_signal(rng, "Leak", "",   "LPM", -147.951f0, 148.4674f0,  imin16, imax16, samples_per_record, n_records)
        test_edf_signal(rng, "CHEST", "E",  "uV", -32768.0f0, 32767.0f0,   imin16, imax16, samples_per_record, n_records)
        test_edf_signal(rng, "ABD", "E",    "uV", -32768.0f0, 32767.0f0,   imin16, imax16, samples_per_record, n_records)
        test_edf_signal(rng, "Tidal", "",   "mL", -4928.18f0, 4906.871f0,  imin16, imax16, samples_per_record, n_records)
        test_edf_signal(rng, "SaO2", "",    "%",       0.0f0,    100.0f0,  imin16, imax16, samples_per_record, n_records)
        EDF.AnnotationsSignal(samples_per_record, anns_2)
        test_edf_signal(rng, "EKG EKGR- REF", "E",  "uV",   -9324.0f0, 2034.0f0,  imin16, imax16, samples_per_record, n_records)
        test_edf_signal(rng, "IC", "E",    "uV",    -32768.0f0, 32767.0f0, imin16, imax16, samples_per_record, n_records)
        test_edf_signal(rng, "HR", "",    "BpM",    -32768.0f0, 32768.0f0, imin16, imax16, samples_per_record, n_records)
        test_edf_signal(rng, "EcG EKGL", "E",  "uV",   -10932.0f0, 1123.0f0,  imin16, imax16, samples_per_record, n_records)
        test_edf_signal(rng, "- REF", "E",  "uV",   -10932.0f0, 1123.0f0,  imin16, imax16, samples_per_record, n_records)
        test_edf_signal(rng, "REF1", "E",  "uV",   -10932.0f0, 1123.0f0,  imin16, imax16, samples_per_record, n_records)
    ]
    seconds_per_record = samples_per_record / sample_rate
    edf_header = EDF.FileHeader("0", "", "", DateTime("2014-10-27T22:24:28"), true, n_records, seconds_per_record)
    edf = EDF.File((io = IOBuffer(); close(io); io), edf_header, edf_signals)
    return edf, Dict(:eeg => [9, 1, 2, 3, 5, 4, 6],
                     :eog => [7, 8],
                     :ecg => [27, 24],
                     :emg => [25, 11, 12],
                     :heart_rate => [26],
                     :tidal_volume => [21],
                     :respiratory_effort => [19, 20],
                     :snore => [13],
                     :positive_airway_pressure => [14, 15],
                     :pap_device_leak => [18],
                     :pap_device_cflow => [16],
                     :sao2 => [22],
                     :ptaf => [17])
end

include("test_edf_to_samples_info.in")

@testset "OndaEDF" begin
    include("signal_labels.jl")
    include("import.jl")
    include("export.jl")
end

@info """
    To look for the effect of any modifications made to OndaEDF, look at:
    `diff test/test_edf_to_samples_info.{in,out}`
"""


text:
module SerializationTests

using Base.Test
using BenchmarkTools

old_data = BenchmarkTools.load(joinpath(dirname(@__FILE__), "data_pre_v006.jld"), "results")
BenchmarkTools.save(joinpath(dirname(@__FILE__), "tmp.jld"), "results", old_data)
new_data = BenchmarkTools.load(joinpath(dirname(@__FILE__), "tmp.jld"), "results")

@test old_data == new_data

rm(joinpath(dirname(@__FILE__), "tmp.jld"))

end # module


text:
"""
    liftingline_rigidbody_model()

Construct a model by coupling a lifting line aerodynamic model and a rigid body
dynamics model.  This model introduces additional parameters corresponding to
the length, position, and orientation of each lifting line element ``L, p_e,
e_1, e_2, e_3``, followed by the freestream air density ``rho_\\infty``, air speed of sound
``c``, rigid body inertial properties ``m, Ixx, Iyy, Izz, Ixz, Ixy, Iyz``, and the
additional forces/moments applied on the body ``F_x, F_y, F_z, M_x, M_y, M_z``.

**NOTE: When using this model, the local frame for each lifting line element should be
oriented with the x-axis in the chordwise direction, the y-axis in the spanwise
direction (out the right wing), and the z-axis in the airfoil normal direction**
"""
function liftingline_rigidbody_model(section_models)

    # aerodynamic model
    aero = LiftingLine(section_models)

    # structural model
    dyn = RigidBody()

    # submodels
    submodels = (aero, dyn)

    # construct coupling
    coupling = liftingline_rigidbody_coupling(aero, dyn)

    return CoupledModel(submodels, coupling)
end

# --- Internal Methods --- #

function liftingline_rigidbody_coupling(aero, dyn)

    # extract section models
    section_models = aero.constants.submodels

    # number of aerodynamic sections
    n = length(section_models)

    # number of additional parameters
    npc = 13*n + 15

    # dimensions
    nx = number_of_states(aero) + number_of_states(dyn)
    ny = number_of_inputs(aero) + number_of_inputs(dyn)
    np = number_of_parameters(aero) + number_of_parameters(dyn) + npc

    # coupling function
    g = (dx, x, p, t) -> liftingline_rigidbody_inputs(dx, x, p, t; aero, dyn)

    # jacobians
    ratejac = Nonlinear() # TODO: define rate jacobian
    statejac = Nonlinear() # TODO: define state jacobian
    paramjac = Nonlinear() # TODO: define parameter jacobian
    tgrad = Nonlinear() # TODO: define time gradient

    # input/output functions
    setparam = liftingline_rigidbody_setparam!
    sepparam = (p) -> liftingline_rigidbody_sepparam(p; n)

    return Coupling{true}(g, nx, ny, np, npc;
        ratejac = ratejac,
        statejac = statejac,
        paramjac = paramjac,
        tgrad = tgrad,
        setparam = setparam,
        sepparam = sepparam)
end

function liftingline_rigidbody_inputs!(y, dx, x, p, t; aero, dyn)

    # number of aerodynamic sections
    n = length(aero.constants.submodels)

    # get indices of states, inputs, and parameters for each submodel
    submodels = (aero, dyn)
    ixa, ixd = state_indices(submodels)
    iya, iyd = input_indices(submodels)
    ipa, ipd = parameter_indices(submodels)

    # get views of states, inputs, and parameters for each submodel
    dxa, dxd = view(dx, ixa), view(dx, ixd)
    xa, xd = view(x, ixa), view(x, ixd)
    ya, yd = view(y, iya), view(y, iyd)
    pa, pd = view(p, ipa), view(p, ipd)

    # separate lifting line model rates, states, and parameters
    liftingline_rates = separate_rates(aero, dxa)
    liftingline_states = separate_states(aero, xa)
    liftingline_parameters = separate_parameters(aero, pa)

    # separate rigid body model rates, states, and parameters
    rigidbody_rates = separate_rates(aero, dxd)
    rigidbody_states = separate_states(aero, xd)
    rigidbody_parameters = separate_parameters(aero, pd)

    # separate additional parameters
    np = length(p)
    ipc = np - 13*n - npc : np
    pc = view(p, ipc)
    additional_parameters = liftingline_rigidbody_sepparam(pc; n)

    # initialize loads with body loads
    F = SVector(additional_parameters.Fbx, additional_parameters.Fby, additional_parameters.Fbz)
    M = SVector(additional_parameters.Mbx, additional_parameters.Mby, additional_parameters.Mbz)

    # find environmental conditions
    V, Ω, a, α = rigidbody_environmental_conditions(rigidbody_rates, rigidbody_states)

    # loop through each section
    for i = 1:n

        # get additional parameters for this section
        param = additional_parameters.section_parameters[i]

        # unpack section properties
        L = param.L
        pe = param.pe
        e1 = param.e1
        e2 = param.e2
        e3 = param.e3

        # unpack section aerodynamic model, rates, states, and parameters
        aerodynamic_model = section_models[i].models[1]
        aerodynamic_rates = liftingline_rates.sections[i]
        aerodynamic_states = liftingline_states.sections[i]
        aerodynamic_parameters = liftingline_parameters.sections[i]

        # determine section velocities/accelerations
        vi, ωi, ai, αi = liftingline_section_velocities(pe, e1, e2, e3, V, Ω, a, α)

        # express section aerodynamic rates, states, and parameters as vectors
        dxai = get_rates(aerodynamic_model; aerodynamic_rates...)
        xai = get_states(aerodynamic_model; aerodynamic_states...)
        pai = get_parameters(aerodynamic_model; aerodynamic_parameters...)

        # express section structural rates, states, and parameters as vectors
        dxsi = SVector(ai..., αi...)
        xsi = SVector(vi..., ωi...)
        psi = SVector(ρ, c)

        # define section rate, state, and parameter vectors
        dxi = SVector(dxai..., dxsi...)
        xi = SVector(xai..., xsi...)
        pi = SVector(pai..., psi...)

        # calculate section inputs
        yi = get_coupling_inputs(section_models[i], dxi, xi, pi, t)

        # separate section inputs
        section_inputs = separate_inputs(section_models[i], yi)

        # separate section aerodynamic and structural inputs
        aerodynamic_inputs = section_inputs[1]
        structural_inputs = section_inputs[2]

        # save section aerodynamic inputs
        set_inputs!(yai, aerodynamic_model, aerodynamic_inputs)

        # integrate section distributed loads
        Fi = L*structural_inputs.fi
        Mi = L*structural_inputs.mi

        # add section loads to total loads
        F += Fi
        M += cross(pe, Fi) + Mi
    end

    # save rigid body inputs
    set_inputs!(yd, rigidbody_model;
        m = additional_parameters.m,
        Ixx = additional_parameters.Ixx,
        Iyy = additional_parameters.Iyy,
        Izz = additional_parameters.Izz,
        Ixz = additional_parameters.Ixz,
        Ixy = additional_parameters.Ixy,
        Iyz = additional_parameters.Iyz,
        Fx = F[1],
        Fy = F[2],
        Fz = F[3],
        Mx = M[1],
        My = M[2],
        Mz = M[3])

    return y
end

function liftingline_rigidbody_setparam!(p; section_parameters, rho, c,
    m, Ixx, Iyy, Izz, Ixz, Ixy, Iyz, Fbx, Fby, Fbz, Mbx, Mby, Mbz)

    n = length(section_parameters)

    for i = 1:n
        offset = 13*(i-1)
        p[offset+1] = section_parameters[i].L
        p[offset+2:offset+4] = section_parameters[i].pe
        p[offset+5:offset+7] = section_parameters[i].e1
        p[offset+8:offset+10] = section_parameters[i].e2
        p[offset+11:offset+13] = section_parameters[i].e3
    end

    offset = 13*(n-1)
    p[offset+1] = rho
    p[offset+2] = c
    p[offset+3] = m
    p[offset+4] = Ixx
    p[offset+5] = Iyy
    p[offset+6] = Izz
    p[offset+7] = Ixz
    p[offset+8] = Ixy
    p[offset+9] = Iyz
    p[offset+10] = Fbx
    p[offset+11] = Fby
    p[offset+12] = Fbz
    p[offset+13] = Mbx
    p[offset+14] = Mby
    p[offset+15] = Mbz

    return p
end

liftingline_rigidbody_sepparam(p; n) = (
    section_parameters = ntuple(i -> (
        L = p[13*(i-1)+1], # element length
        pe = view(p, 13*(i-1)+2 : 13*(i-1)+4), # local frame position
        e1 = view(p, 13*(i-1)+5 : 13*(i-1)+7), # local frame x-axis
        e2 = view(p, 13*(i-1)+8 : 13*(i-1)+10), # local frame y-axis
        e3 = view(p, 13*(i-1)+11 : 13*(i-1)+13), # local frame z-axis
        ), n),
    rho = p[13*n+1],
    c = p[13*n+2],
    m = p[13*n+3],
    Ixx = p[13*n+4],
    Iyy = p[13*n+5],
    Izz = p[13*n+6],
    Ixz = p[13*n+7],
    Ixy = p[13*n+8],
    Iyz = p[13*n+9],
    Fbx = p[13*n+10],
    Fby = p[13*n+11],
    Fbz = p[13*n+12],
    Mbx = p[13*n+13],
    Mby = p[13*n+14],
    Mbz = p[13*n+15],
    )


text:
function Dominated(Point, S,s)
    for i = 1:length(S)
        if i == s; continue; end
        if all(S[i].Point .<= Point)
            return true
        end
    end
    return false
end


text:
function SeisWindowHeadersPatch(in, out; key=[], minval=[], maxval=[], tmin=0,
                           tmax=99999, ntrace=1000000,it_nt=9e9)

     itrace_in = 1
     itrace_out = 1
     NX = GetNumTraces(in)

     group="some"

     while itrace_in <= NX
         nx = NX - itrace_in + 1
	       ntrace = nx > ntrace ? ntrace : nx

	       h1 = SeisReadHeaders(in, group=group, key=key, itrace=itrace_in,
                                 ntrace=ntrace)
	       num_traces_in = length(h1)

         h2=WindowHeadersPatch(h1,key=key,minval=minval,maxval=maxval)
	       h1 = copy(h2)

         num_traces_out = length(h1)
	       SeisWriteHeaders(out, h1, itrace=itrace_out,
                             update_tracenum=false)
	       itrace_in += num_traces_in
	       itrace_out += num_traces_out
     end

    DATAPATH = get(ENV,"DATAPATH",join([pwd(),"/"]))
    filename_d_out = join([DATAPATH out "@data@"])
    filename_h_out = join([DATAPATH out "@headers@"])
    nhead = length(fieldnames(Header))

    stream_h = open(filename_h_out)
    nx = round(Int,filesize(stream_h)/(nhead*4))

    h = GrabHeader(stream_h,1)
    close(stream_h)

    nt = convert(Int64,round((tmax - tmin)/h.d1)) + 1

    if nt > h.n1
	     nt = h.n1
    end
    if nt<it_nt
	nt = nt-1
    end

    minval = convert(Array{Float32,1},vec(minval))
    maxval = convert(Array{Float32,1},vec(maxval))
    key2 = String[]
    minval2 = Float32[]
    maxval2 = Float32[]
    for ikey=1:length(key)
      if key[ikey] != "t"
	       key2 = push!(key2,key[ikey])
	       minval2 = push!(minval2,minval[ikey])
	       maxval2 = push!(maxval2,maxval[ikey])
	    end
    end

    nx1=maxval2[1]-minval2[1]+1
    nx2=maxval2[2]-minval2[2]+1
    nx3=maxval2[3]-minval2[3]+1
    nx4=maxval2[4]-minval2[4]+1

    extent = Extent(convert(Int32,nt),
            convert(Int32,nx1), convert(Int32,nx2),
            convert(Int32,nx3), convert(Int32,nx4),
            convert(Float32,tmin),
            convert(Float32,1), convert(Float32,0), convert(Float32,0),
            convert(Float32,0), convert(Float32,h.d1),
            convert(Float32,1), convert(Float32,1), convert(Float32,1),
            convert(Float32,1), "Time", key[2], key[3], key[4], key[5], "s",
                    "index", "index", "index", "index", "")
    WriteTextHeader(out,extent,"native_float",4,filename_d_out,filename_h_out)

end

function WindowHeadersPatch(h_in;key=[],minval=[],maxval=[])

    minval = convert(Array{Float32,1},vec(minval))
    maxval = convert(Array{Float32,1},vec(maxval))
    nx = length(h_in)
    key2 = String[]
    minval2 = Float32[]
    maxval2 = Float32[]
    for ikey=1:length(key)
	     if key[ikey] != "t"
	        key2 = push!(key2,key[ikey])
	        minval2 = push!(minval2,minval[ikey])
	        maxval2 = push!(maxval2,maxval[ikey])
	     end
    end

    return RejectHeadersPatch(h_in,key2,minval2,maxval2,length(key2),length(h_in))

end

function RejectHeadersPatch(h_in::Array{Header,1}, key::Array{String,1},
                       minval::Array{Float32,1}, maxval::Array{Float32,1},
                       nkeys, nx)
    h_out = Header[]
    keep = true
    key_val = 0f0
    for j=1:nx
	     keep = true
	     for ikey=1:nkeys

	       key_val = convert(Float32,getfield(h_in[j],Symbol(key[ikey])))

	        if (key_val < minval[ikey] || key_val > maxval[ikey])
		          keep = false
	        end
	     end
	     if (keep==true)
	        h_out = push!(h_out,h_in[j])
	     end
    end
  return h_out
end


text:
"""
    compute_coefficients(u, mesh::AbstractMesh1D)

Computse the coefficients of the function `u` in a first order finite volume
method on the `mesh`.
"""
function compute_coefficients(u, mesh::AbstractMesh1D, order::Integer=1)
    xmin, xmax = bounds(mesh)
    uval = zeros(typeof(u((xmin+xmax)/2)), numcells(mesh))
    compute_coefficients!(uval, u, mesh, order)
    uval
end

"""
    compute_coefficients!(uval, u, mesh::AbstractMesh1D)

Computse the coefficients of the function `u` in a first order finite volume
method on the `mesh` and stores them in `uval`.
"""
function compute_coefficients!(uval, u, mesh::AbstractMesh1D, order::Integer=1)
    if eltype(uval) <: AbstractArray
        basis = GaussLegendre(order-1, eltype(eltype(uval)))
    else
        basis = GaussLegendre(order-1, eltype(uval))
    end
    x = similar(basis.nodes)

    for cell in cell_indices(mesh)
        xmin, xmax = bounds(cell, mesh)
        map_from_canonical!(x, basis.nodes, xmin, xmax, basis)
        uval[cell] = integrate(u, x, basis.weights) / 2
    end
    nothing
end

"""
    evaluate_coefficients(u, mesh::AbstractMesh1D)

Evaluates the coefficients `u` in a first order finite volume method on the `mesh`.
Returns `xplot, uplot` as vectors that can be used for plots.
"""
function evaluate_coefficients(u, mesh::AbstractMesh1D)
    xplot = zeros(2*numcells(mesh))
    uplot = zeros(eltype(u), 2*numcells(mesh))

    evaluate_coefficients!(xplot, uplot, u, mesh)
end

"""
    evaluate_coefficients!(xplot, uplot, u, mesh::AbstractMesh1D)

Evaluates the coefficients `u` in a first order finite volume method on the `mesh`.
Returns `xplot, uplot` as vectors that can be used for plots.
"""
function evaluate_coefficients!(xplot, uplot, u, mesh::AbstractMesh1D)
  @assert length(xplot) == 2*numcells(mesh)
  @assert length(uplot) == 2*numcells(mesh)

  idx = 1
  @inbounds for cell in cell_indices(mesh)
    xmin, xmax = bounds(cell, mesh)

    xplot[idx  ] = xmin
    xplot[idx+1] = xmax
    uplot[idx] = uplot[idx+1] = u[cell]

    idx = idx+2
  end

  xplot, uplot
end

"""
    evaluate_coefficients(u, balance_law, mesh::AbstractMesh1D,
                          reconstruction::AbstractReconstruction,
                          npoints=2*order(reconstruction))

Evaluates the coefficients `u` in a finite volume method on `mesh` using
`reconstruction`.
Returns `xplot, uplot` as vectors that can be used for plots.
"""
function evaluate_coefficients(u, balance_law, meshx::AbstractMesh1D,
                               reconstruction::AbstractReconstruction,
                               npoints=2*order(reconstruction))
    xplot = zeros(npoints*numcells(meshx))
    uplot = zeros(eltype(u), npoints*numcells(meshx))

    evaluate_coefficients!(xplot, uplot, u, balance_law, meshx, reconstruction)
end

"""
    evaluate_coefficients!(xplot, uplot, u, balance_law, mesh::AbstractMesh1D,
                           reconstruction::AbstractReconstruction)

Evaluates the coefficients `u` in a finite volume method on `mesh` using
`reconstruction`.
Returns `xplot, uplot` as vectors that can be used for plots.
"""
function evaluate_coefficients!(xplot, uplot, u, balance_law, meshx::AbstractMesh1D,
                                reconstruction::AbstractReconstruction)
    npoints = length(xplot) ÷ numcells(meshx)
    @assert length(uplot) == npoints*numcells(meshx)

    ξ = range(0+eps(), 1-eps(), length=npoints) |> collect
    x = zero(ξ)

    uval = zeros(eltype(u), npoints)
    @inbounds for cell in cell_indices(meshx)
        xmin, xmax = bounds(cell, meshx)
        interpolate!(uval, ξ, cell, u, balance_law, meshx, reconstruction, stencil_width_val(reconstruction))
        for nx in 1:npoints
            xplot[(cell-1)*npoints+(nx-1)+1] = xmin + ξ[nx]*(xmax-xmin)
            uplot[(cell-1)*npoints+(nx-1)+1] = uval[nx]
        end
    end

    xplot, uplot
end

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

function compute_coefficients(u, mesh::AbstractMesh1D, basis::NodalBasis, parallel=Val{:serial}())
    xmin, xmax = bounds(mesh)
    uval = Array{typeof(u((xmin+xmax)/2))}(undef, length(basis.nodes), numcells(mesh))
    compute_coefficients!(uval, u, mesh, basis)
    uval
end

function compute_coefficients!(uval, u, mesh::AbstractMesh1D, basis::NodalBasis)
    @unpack nodes = basis
    Pp1 = length(nodes)

    for cell in cell_indices(mesh)
        xmin, xmax = bounds(cell, mesh)
        for n in 1:Pp1
            x = map_from_canonical(nodes[n], xmin, xmax, basis)
            uval[n,cell] = u(x)
        end
    end

    nothing
end

function evaluate_coefficients(u, meshx::AbstractMesh1D, basis::NodalBasis,
                                npoints=2*length(basis.nodes))
    xplot = zeros(npoints*numcells(meshx))
    uplot = zeros(eltype(u), npoints*numcells(meshx))

    evaluate_coefficients!(xplot, uplot, u, meshx, basis)
end

function evaluate_coefficients!(xplot, uplot, u, meshx::AbstractMesh1D, basis::NodalBasis)
    npoints = length(xplot) ÷ numcells(meshx)

    @assert length(xplot) == npoints*numcells(meshx)
    @assert length(xplot) == length(uplot)

    ξ = range(-1+eps(), 1-eps(), length=npoints) |> collect
    x = zero(ξ)

    Nx  = numcells(meshx)
    Pp1 = length(basis.nodes)
    uval = zeros(eltype(u), npoints)
    for ix in 1:Nx
        xmin, xmax = bounds(ix, meshx)
        map_from_canonical!(x, ξ, xmin, xmax, basis)
        interpolate!(uval, ξ, u[:,ix], basis)

        for jx in 1:npoints
            xplot[(ix-1)*npoints+(jx-1)+1] = x[jx]
            uplot[(ix-1)*npoints+(jx-1)+1] = uval[jx]
        end
    end

    xplot, uplot
end

function compute_coefficients(u, meshx::AbstractMesh1D, meshy::AbstractMesh1D, basis::NodalBasis)
    xmin, xmax = bounds(meshx)
    ymin, ymax = bounds(meshy)
    Pp1 = length(basis.nodes)
    uval = zeros(typeof(u((xmin+xmax)/2,(ymin+ymax)/2)), Pp1, Pp1, numcells(meshx), numcells(meshy))
    compute_coefficients!(uval, u, meshx, meshy, basis)
    uval
end

function compute_coefficients!(uval, u, meshx::AbstractMesh1D, meshy::AbstractMesh1D, basis::NodalBasis)
    @unpack nodes = basis
    Pp1 = length(nodes)

    for iy in cell_indices(meshy), ix in cell_indices(meshx)
        xmin, xmax = bounds(ix, meshx)
        ymin, ymax = bounds(iy, meshy)
        for ny in 1:Pp1, nx in 1:Pp1
            x = map_from_canonical(nodes[nx], xmin, xmax, basis)
            y = map_from_canonical(nodes[ny], ymin, ymax, basis)
            uval[nx,ny,ix,iy] = u(x,y)
        end
    end

    nothing
end

function evaluate_coefficients(u, meshx::AbstractMesh1D, meshy::AbstractMesh1D,
                                basis::NodalBasis, npoints=2*length(basis.nodes))
    xplot = zeros(npoints*numcells(meshx))
    yplot = zeros(npoints*numcells(meshy))
    uplot = zeros(eltype(u), npoints*numcells(meshy), npoints*numcells(meshx))

    evaluate_coefficients!(xplot, yplot, uplot, u, meshx, meshy, basis)
end

function evaluate_coefficients!(xplot, yplot, uplot, u, meshx::AbstractMesh1D,
                                meshy::AbstractMesh1D, basis::NodalBasis)
    npoints = length(xplot) ÷ numcells(meshx)

    @assert length(xplot) == npoints*numcells(meshx)
    @assert length(yplot) == npoints*numcells(meshy)
    @assert length(xplot) == size(uplot,2)
    @assert length(yplot) == size(uplot,1)

    ξ = range(-1+eps(), 1-eps(), length=npoints) |> collect
    η = range(-1+eps(), 1-eps(), length=npoints) |> collect
    intX = Jacobi.interp_mat(ξ, basis.nodes)
    intY = Jacobi.interp_mat(η, basis.nodes)

    Nx = numcells(meshx)
    Ny = numcells(meshy)
    Pp1 = length(basis.nodes)
    tmp = zeros(eltype(u), Pp1, npoints)
    uval = zeros(eltype(u), npoints, npoints)
    for iy in 1:Ny, ix in 1:Nx
        xmin, xmax = bounds(ix, meshx)
        ymin, ymax = bounds(iy, meshy)
        x = map_from_canonical(ξ, xmin, xmax, basis)
        y = map_from_canonical(η, ymin, ymax, basis)

        for nx in 1:Pp1
            A_mul_B!(view(tmp,nx,:), intY, view(u,nx,:,ix,iy))
        end
        A_mul_B!(uval, intX, tmp)

        for jy in 1:npoints, jx in 1:npoints
            xplot[(ix-1)*npoints+(jx-1)+1] = x[jx]
            yplot[(iy-1)*npoints+(jy-1)+1] = y[jy]
            uplot[(iy-1)*npoints+(jy-1)+1, (ix-1)*npoints+(jx-1)+1] = uval[jx,jy]
        end
    end

    xplot, yplot, uplot
end

function compute_coefficients(u, meshx::AbstractMesh1D, meshy::AbstractMesh1D,
                                 meshz::AbstractMesh1D, basis::NodalBasis)
    xmin, xmax = bounds(meshx)
    ymin, ymax = bounds(meshy)
    zmin, zmax = bounds(meshz)
    Pp1 = length(basis.nodes)
    uval = zeros(typeof(u((xmin+xmax)/2,(ymin+ymax)/2,(zmin+zmax)/2)), Pp1, Pp1, Pp1,
                    numcells(meshx), numcells(meshy), numcells(meshz))
    compute_coefficients!(uval, u, meshx, meshy, meshz, basis)
    uval
end

function compute_coefficients!(uval, u, meshx::AbstractMesh1D, meshy::AbstractMesh1D,
                                        meshz::AbstractMesh1D, basis::NodalBasis)
    @unpack nodes = basis
    Pp1 = length(nodes)

    for iz in cell_indices(meshz), iy in cell_indices(meshy), ix in cell_indices(meshx)
        xmin, xmax = bounds(ix, meshx)
        ymin, ymax = bounds(iy, meshy)
        zmin, zmax = bounds(iz, meshz)
        for nz in 1:Pp1, ny in 1:Pp1, nx in 1:Pp1
            x = map_from_canonical(nodes[nx], xmin, xmax, basis)
            y = map_from_canonical(nodes[ny], ymin, ymax, basis)
            z = map_from_canonical(nodes[nz], zmin, zmax, basis)
            uval[nx,ny,nz,ix,iy,iz] = u(x,y,z)
        end
    end

    nothing
end


text:
using Nemo, CharFuncPricing, Test

@testset "HestonALSmallPrices" begin
    R = ArbField(512)
    CC = AcbField(512)
    r = R(QQ(0, 100))
    q = R(QQ(0, 100))
    κ = R(2)
    θ = R(QQ(1, 10000))
    σ = R(QQ(1, 10000))
    v0 = R(QQ(1, 10000))
    ρ = R(-0.5)
    τ = R(0.0025)
    strike = R(100)
    spots = [R("100.0001"), R(101.0)]#, R(110.0), R(200.0), R(1000.0), R(10000.0)]
    params = HestonParams{arb}(v0, κ, θ, ρ, σ)

    l = 128 #48
    m = 1024 * 16
    cf = NemoCharFunc{HestonParams{arb}}(params, CC)
    pricer = makeCosCharFuncPricer(cf, τ, m, l)
    prices = Vector{arb}(undef, length(spots))
    for (i, spot) in enumerate(spots)
        prices[i] = priceEuropean(pricer, false, strike, spot, τ, R(1))
        println(spot, " ", prices[i])
    end
    setprecision(512)
    κ = big(2.0)
    θ = big(1.0)/big(10000.0)
    σ =  big(1.0)/big(10000.0)
    v0 =  big(1.0)/big(10000.0)
    ρ = big(-0.5)
    τ = big(0.0025)
    strike = big(100.0)
    spots = [BigFloat("100.0001"), big(101.0), big(110.0), big(200.0), big(1000.0), big(10000.0)]
    params = HestonParams{BigFloat}(v0, κ, θ, ρ, σ)
    cf = DefaultCharFunc{HestonParams{BigFloat},Complex{BigFloat}}(params)
    pricer = makeCosCharFuncPricer(cf, τ, m, l)
    pricesB = Vector{BigFloat}(undef, length(spots))
    for (i, spot) in enumerate(spots)
        pricesB[i] = priceEuropean(pricer, false, strike, spot, τ, big(1.0))
        println(spot, " ", pricesB[i])
    end
    @test isapprox(Float64(  prices[1]- R(pricesB[1])), 0.0, atol = 1e-150)
    @test isapprox(Float64( prices[2] - R(pricesB[2])), 0.0, atol = 1e-150)

    #checks that bigfloat implementation of andersen lake is correct (use of pi)
    quad = TanhSinhQuadrature(10000, BigFloat("1e-300"))
    pricer = ALCharFuncPricer(cf,quad)
    pricesC = Vector{BigFloat}(undef, length(spots))
    for (i, spot) in enumerate(spots)
        pricesC[i] = priceEuropean(pricer, false, strike, spot, τ, big(1.0))
        println(spot, " ", pricesC[i])
    end
    @test isapprox(  pricesB[1], pricesC[1], atol = BigFloat(1e-130))
    @test isapprox( pricesB[2] ,pricesC[2], atol = BigFloat(1e-130))

end

@testset "AlanSet" begin
    R = ArbField(256)
    CC = AcbField(256)
    r = R(QQ(1, 100))
    q = R(QQ(2, 100))
    κ = R(4.0)
    θ = R(QQ(1, 4))
    σ = R(1.0)
    v0 = R(QQ(4, 100))
    ρ = R(-0.5)
    τ = R(1.0)
    spot = R(100)
    strikes = [R(80), R(90), R(100), R(110), R(120)]
    alanPuts = [
        R("7.958878113256768285213263077598987193482161301733"),
        R("12.017966707346304987709573290236471654992071308187"),
        R("17.055270961270109413522653999411000974895436309183"),
        R("23.017825898442800538908781834822560777763225722188"),
        R("29.811026202682471843340682293165857439167301370697"),
    ]
    alanCalls = [
        R("26.774758743998854221382195325726949201687074848341"),
        R("20.933349000596710388139445766564068085476194042256"),
        R("16.070154917028834278213466703938231827658768230714"),
        R("12.132211516709844867860534767549426052805766831181"),
        R("9.024913483457835636553375454092357136489051667150"),
    ]
    #refPuts = [R("7.95887811325676828521326060761429303089865693725960319205094095681790030"), R("12.01796670734630498770957286168505994488766901674359014792496983503358626"), R("17.05527096127010941352265395975862503115076404332454551322405076347166800"),R("23.01782589844280053890878183747069958885882262385534095209843808561556738"), R("29.81102620268247184334068229213413372849045238688265116546863785487435207")]
    #refCalls = [R("26.77475874399885422138219285574225503910357048386701361386103443917276629"),R("20.93334900059671038813944533801265637537179175081262590089627586445697497"), R("16.07015491702883427821346666428585588391409596485520659735656933996357944"), R("12.13221151670984486786053477019756486390136373284762736739216920917600154"), R("9.02491348345783563655337545306063342581220268333656291192358152550330896")]
    spot *= exp((r - q) * τ)
    df = exp(-r * τ)
    params = HestonParams{arb}(v0, κ, θ, ρ, σ)

    l = 48
    m = 1024 * 4
    cf = NemoCharFunc(params, CC)
    pricer = makeCosCharFuncPricer(cf, τ, m, l)
    for (strike, refCall, refPut) in zip(strikes, alanCalls, alanPuts)
        price = priceEuropean(pricer, false, strike, spot, τ, df)
        println(Float64(strike), " P ", price, " ", price - refPut)
        #push!(refPuts,price)
        @test isapprox(Float64(price - refPut), 0, atol = 1e-23)
        price = priceEuropean(pricer, true, strike, spot,  τ,df)
        println(Float64(strike), " C ", price, " ", price - refCall)
        @test isapprox(Float64(price - refCall), 0, atol = 1e-23)
        #push!(refCalls,price)
    end
end

@testset "AlanJoshiSet" begin
    R = ArbField(256)
    CC = AcbField(256)
    r = R(QQ(1, 100))
    q = R(QQ(2, 100))
    κ = R(4.0)
    θ = R(QQ(1, 4))
    σ = R(1.0)
    v0 = R(QQ(1, 100))
    ρ = R(-0.5)
    τ = R(QQ(1, 100))
    spot = R(100)
    strikes = [R(90), R(95), R(100), R(105), R(110)]
    alanPuts = [
        R("4.5183603586861772614990106188215872180542e-8"),
        R("0.000461954855653851579672612557018857858641926937"),
        R("0.477781171629504680023239655436072890669645669297"),
        R("5.009501052563650299130635110520904481889436667608"),
        R("10.008998550115123724684210555728039829315964456261"),
    ]
    alanCalls = [
        R("9.989001595065276544935948045293485530832966049263"),
        R("4.989963479738160122154264702582719627807098780529"),
        R("0.467782671512844263098248405184095087949465507760"),
        R("2.527447823194706060519991248106500619490942e-6"),
        R("1.29932760052624920704881258510264466e-13"),
    ]
    refPuts = [
        R("4.518360358686177267846283260505342170483575627777422720994e-8"),
        R("0.00046195485565385157967261251100128843593583885256695108090334577"),
        R("0.47778117162950468002323965543594159177844084992306552196870005758"),
        R("5.00950105256365029913063511052075051116901889480844354030472812057"),
        R("10.00899855011512372468421052192118981283250808930047516269268570714"),
    ]
    refCalls = [
        R("9.98900159506527654493594810876621194767051557355651881698849965467"),
        R("4.989963479738160122154264702536702058384392692444362488921902796131"),
        R("0.46778267151284426309824840518396378905826068838589405843642651357"),
        R("2.52744782319470606051999109413578020171814230507539918158220e-6"),
        R("1.2993276005259111385486477505389750536969641386617440e-13"),
    ]
    spot *= exp((r - q) * τ)
    df = exp(-r * τ)
    params = HestonParams{arb}(v0, κ, θ, ρ, σ)
    l = 48
    m = 1024 * 4
    pricer = makeCosCharFuncPricer(
        NemoCharFunc(params, CC),
        τ,
        m,
        l,
    )
    for (strike, refCall, refPut) in zip(strikes, alanCalls, alanPuts)
        price = priceEuropean(pricer, false, strike, spot, τ, df)
        println(Float64(strike), " P ", price, " ", price - refPut)
        #push!(refPuts,price)
        @test isapprox(Float64(price - refPut), 0, atol = 1e-25)
        price = priceEuropean(pricer, true, strike, spot, τ, df)
        println(Float64(strike), " C ", price, " ", price - refCall)
        @test isapprox(Float64(price - refCall), 0, atol = 1e-25)
        #push!(refCalls,price)
    end
end

@testset "ZeroKappa" begin
    R = ArbField(256)
    CC = AcbField(256)
    κ = R(0.0)
    θ = R(QQ(1, 4))
    σ = R(1.0)
    v0 = R(QQ(4, 100))
    ρ = R(-0.5)
    τ = R(1.0)
    params = HestonParams{arb}(v0, κ, θ, ρ, σ)
    cumulantsZ = computeCumulants(params, τ)
    κ = R(QQ(1, 100))
    params = HestonParams{arb}(v0, κ, θ, ρ, σ)
    cumulants = computeCumulants(params, τ)

    for (cz, c) in zip(cumulantsZ, cumulants)
        println(cz, " ", c)
        @test isapprox(Float64(cz), Float64(c), rtol = 1e-1)
    end
end


text:
function get_entries(file, splits)
    file = abspath(file)
    entries = JSON.parsefile(file)["images"]
    return map(s->filter(x->x["split"]==s, entries), splits)
end

function get_pairs(entries)
    pairs = []
    for entry in entries
        for sentence in entry["sentences"]
            push!(pairs, (entry["filename"], sentence["tokens"]))
        end
    end
    return pairs
end

max_stars_repo_name: ilkerkesen/kusrp

max_issues_repo_name: ilkerkesen/kusrp

max_forks_repo_name: ilkerkesen/kusrp


text:
# The following simple model has a variable number of agent types,
# but there is no killing or creating of additional agents.
# It creates a model that has the same number of agents and does
# overall the same number of operations, but these operations
# are split in a varying number of agents. It shows how much of a
# performance hit is to have many different agent types.

using Agents, Random

mutable struct Agent1 <: AbstractAgent
    id::Int
    pos::Tuple{Int,Int}
    money::Int
end

mutable struct Agent2 <: AbstractAgent
    id::Int
    pos::Tuple{Int,Int}
    money::Int
end

mutable struct Agent3 <: AbstractAgent
    id::Int
    pos::Tuple{Int,Int}
    money::Int
end

mutable struct Agent4 <: AbstractAgent
    id::Int
    pos::Tuple{Int,Int}
    money::Int
end

mutable struct Agent5 <: AbstractAgent
    id::Int
    pos::Tuple{Int,Int}
    money::Int
end

mutable struct Agent6 <: AbstractAgent
    id::Int
    pos::Tuple{Int,Int}
    money::Int
end

mutable struct Agent7 <: AbstractAgent
    id::Int
    pos::Tuple{Int,Int}
    money::Int
end

mutable struct Agent8 <: AbstractAgent
    id::Int
    pos::Tuple{Int,Int}
    money::Int
end

mutable struct Agent9 <: AbstractAgent
    id::Int
    pos::Tuple{Int,Int}
    money::Int
end

mutable struct Agent10 <: AbstractAgent
    id::Int
    pos::Tuple{Int,Int}
    money::Int
end

mutable struct Agent11 <: AbstractAgent
    id::Int
    pos::Tuple{Int,Int}
    money::Int
end

mutable struct Agent12 <: AbstractAgent
    id::Int
    pos::Tuple{Int,Int}
    money::Int
end

mutable struct Agent13 <: AbstractAgent
    id::Int
    pos::Tuple{Int,Int}
    money::Int
end

mutable struct Agent14 <: AbstractAgent
    id::Int
    pos::Tuple{Int,Int}
    money::Int
end

mutable struct Agent15 <: AbstractAgent
    id::Int
    pos::Tuple{Int,Int}
    money::Int
end

function initialize_model_1(;n_agents=600,dims=(5,5))
    space = GridSpace(dims)
    model = ABM(Agent1, space; scheduler=Schedulers.randomly, warn=false)
    id = 0
    for id in 1:n_agents
        agent = Agent1(id, (0,0), 10)
        add_agent!(agent, model)
    end
    return model
end

function initialize_model(;n_agents=600, n_types=1, dims=(5,5))
    agent_types = [Agent1,Agent2,Agent3,Agent4,Agent5,Agent6,Agent7,Agent8,
        Agent9,Agent10,Agent11,Agent12,Agent13,Agent14,Agent15]
    agents_used = agent_types[1:n_types]
    space = GridSpace(dims)
    model = ABM(Union{agents_used...}, space; scheduler=Schedulers.randomly, warn=false)
    id = 0
    agents_per_type = div(n_agents, n_types)
    for A in agents_used
        for _ in 1:agents_per_type
            id += 1
            agent = A(id, (0,0), 10)
            add_agent!(agent, model)
        end
    end
    return model
end

function agent_step!(agent, model)
    move!(agent, model)
    agents = agents_in_position(agent.pos, model)
    for a in agents; exchange!(agent, a); end
    return nothing
end

function move!(agent, model)
    neighbors = nearby_positions(agent, model)
    cell = rand(collect(neighbors))
    move_agent!(agent, cell, model)
    return nothing
end

function exchange!(agent, other_agent)
    v1 = agent.money
    v2 = other_agent.money
    agent.money = v2
    other_agent.money = v1
    return nothing
end

function run_simulation_1(n_steps, n_reps)
    t = @timed for _ in 1:n_reps
        model = initialize_model_1()
        Agents.step!(model, agent_step!, n_steps)
    end
    return t[2]/n_reps
end

function run_simulation(n_steps, n_reps; n_types)
    t = @timed for _ in 1:n_reps
        model = initialize_model(;n_types=n_types)
        Agents.step!(model, agent_step!, n_steps)
    end
    return t[2]/n_reps
end

# %% Run the simulation, do performance estimate, first with 1, then with many
Random.seed!(2514)
n_steps = 500
n_reps = 5
n_types = [2,3,4,5,10,15]
# compile
run_simulation_1(1, 1)
for n in n_types; run_simulation(1, 1; n_types=n); end  # compile

time_1 = run_simulation_1(n_steps, n_reps)
times = Float64[]
for n in n_types
    println(n)
    t = run_simulation(n_steps, n_reps; n_types=n)
    push!(times, t/time_1)
end

import CairoMakie
fig, ax, = CairoMakie.lines(n_types, times)
CairoMakie.scatter!(ax, n_types, times)
ax.xlabel = "# types"
ax.ylabel = "time, relative to 1 type"
fig


text:
using Documenter, Turing
using LibGit2: clone

# Include the update_homepage function.
include("homepage-updater.jl")

# Get paths.
examples_path = joinpath(@__DIR__, joinpath("site", "_tutorials"))
source_path = joinpath(@__DIR__, "src")
build_relative = joinpath("site", "_docs")
build_path = joinpath(@__DIR__, build_relative)

# You can skip this part if you are on a metered
# connection by calling `julia make.jl no-tutorials`
if !in("no-tutorials", ARGS)
    isdir(examples_path) || mkpath(examples_path)

    # Clone TuringTurorials
    tmp_path = tempname()
    mkdir(tmp_path)
    clone("https://github.com/TuringLang/TuringTutorials", tmp_path)

    # Move to markdown folder.
    md_path = joinpath(tmp_path, "markdown")

    # Copy the .md versions of all examples.
    try
        println(md_path)
        for file in readdir(md_path)
            full_path = joinpath(md_path, file)
            target_path = joinpath(examples_path, file)
            println("Copying $full_path to $target_path")
            cp(full_path, target_path, force = true)
        end
    catch e
        # println("Markdown copy error: $e")
        rethrow(e)
    finally
        rm(tmp_path, recursive = true)
    end
end

# Preprocess markdown files.
src_temp = mktempdir()
cp(source_path, src_temp, force = true)
yaml_dict = preprocess_markdown(source_path)

# Build documentation
try
    makedocs(
        build = build_relative
    )
catch e
    # Put back the original files in the event of an error.
    cp(src_temp, source_path, force = true)
    rm(src_temp, recursive = true)
    rethrow(e)
end

# Postprocess markdown files (put the YAML headers back in)
cp(src_temp, source_path, force = true)
rm(src_temp, recursive = true)
postprocess_markdown(build_path, yaml_dict, original = source_path)

if !in("no-publish", ARGS)
    # Define homepage update function.
    page_update = update_homepage(
        "github.com/TuringLang/Turing.jl.git",
        "gh-pages",
        "site"
    )
end

# # Deploy documentation.
# deploydocs(
#     repo = "github.com/TuringLang/Turing.jl.git",
#     target = "build",
#     deps = nothing,
#     make = nothing,
#     julia = "1.0"
# )