text:
using LinearAlgebra
using Test

using ProximalOperators
using ProximalAlgorithms:
    DouglasRachfordIteration, DRLSIteration,
    ForwardBackwardIteration, PANOCIteration,
    PANOCplusIteration,
    NoAcceleration,PnpDrsIteration

@testset "DR/DRLS equivalence ($T)" for T in [Float32, Float64]
    A = T[
        1.0 -2.0 3.0 -4.0 5.0
        2.0 -1.0 0.0 -1.0 3.0
        -1.0 0.0 4.0 -3.0 2.0
        -1.0 -1.0 -1.0 1.0 3.0
    ]
    b = T[1.0, 2.0, 3.0, 4.0]

    m, n = size(A)

    R = real(T)

    lam = R(0.1) * norm(A' * b, Inf)

    f = LeastSquares(A, b)
    g = NormL1(lam)

    x0 = zeros(R, n)

    dr_iter = DouglasRachfordIteration(f=f, g=g, x0=x0, gamma=R(10) / opnorm(A)^2)
    drls_iter = DRLSIteration(f=f, g=g, x0=x0, gamma=R(10) / opnorm(A)^2, lambda=R(1), c=-R(Inf), max_backtracks=1, directions=NoAcceleration())

    for (dr_state, drls_state) in Iterators.take(zip(dr_iter, drls_iter), 10)
        @test isapprox(dr_state.x, drls_state.xbar)
    end
end

@testset "Pnp-DRS equivalence ($T)" for T in [Float32, Float64]
    A = T[
        1.0 -2.0 3.0 -4.0 5.0
        2.0 -1.0 0.0 -1.0 3.0
        -1.0 0.0 4.0 -3.0 2.0
        -1.0 -1.0 -1.0 1.0 3.0
    ]
    b = T[1.0, 2.0, 3.0, 4.0]

    m, n = size(A)

    R = real(T)

    lam = R(0.1) * norm(A' * b, Inf)

    f = LeastSquares(A, b)
    g = NormL1(lam)

    x0 = zeros(R, n)

    y = similar(x0)

    gamma=R(10) / opnorm(A)^2

    proxf!(xhat,uhat) = prox!(xhat,f,uhat,gamma) #in place prox updates xhat
    denoiser!(yhat,rhat) = prox!(yhat,g,rhat,gamma) #rhat^{k+1} = 2xhat^{k+1} - uhat^{k}

    dr_iter = DouglasRachfordIteration(f=f, g=g, x0=x0, gamma=gamma)

    pnp_dr_iter = PnpDrsIteration(proxf! = proxf!, denoiser! = denoiser!, uhat0 = x0, gamma = gamma)

    for (dr_state, pnp_dr_state) in Iterators.take(zip(dr_iter, pnp_dr_iter), 10)
        @test isapprox(dr_state.x, pnp_dr_state.uhat)
        @test isapprox(dr_state.y, pnp_dr_state.xhat)
        @test isapprox(dr_state.z, pnp_dr_state.yhat)
    end
end

@testset "FB/PANOC equivalence ($T)" for T in [Float32, Float64]
    A = T[
        1.0 -2.0 3.0 -4.0 5.0
        2.0 -1.0 0.0 -1.0 3.0
        -1.0 0.0 4.0 -3.0 2.0
        -1.0 -1.0 -1.0 1.0 3.0
    ]
    b = T[1.0, 2.0, 3.0, 4.0]

    m, n = size(A)

    R = real(T)

    lam = R(0.1) * norm(A' * b, Inf)

    f = LeastSquares(A, b)
    g = NormL1(lam)

    x0 = zeros(R, n)

    fb_iter = ForwardBackwardIteration(f=f, g=g, x0=x0, gamma=R(0.95) / opnorm(A)^2)
    panoc_iter = PANOCIteration(f=f, g=g, x0=x0, gamma=R(0.95) / opnorm(A)^2, max_backtracks=1, directions=NoAcceleration())

    for (fb_state, panoc_state) in Iterators.take(zip(fb_iter, panoc_iter), 10)
        @test isapprox(fb_state.z, panoc_state.z)
    end
end

@testset "PANOC/PANOCplus equivalence ($T)" for T in [Float32, Float64]
    A = T[
        1.0 -2.0 3.0 -4.0 5.0
        2.0 -1.0 0.0 -1.0 3.0
        -1.0 0.0 4.0 -3.0 2.0
        -1.0 -1.0 -1.0 1.0 3.0
    ]
    b = T[1.0, 2.0, 3.0, 4.0]

    m, n = size(A)

    R = real(T)

    lam = R(0.1) * norm(A' * b, Inf)

    f = LeastSquares(A, b)
    g = NormL1(lam)

    x0 = zeros(R, n)

    panoc_iter = PANOCIteration(f=f, g=g, x0=x0, gamma=R(0.95) / opnorm(A)^2)
    panocplus_iter = PANOCplusIteration(f=f, g=g, x0=x0, gamma=R(0.95) / opnorm(A)^2)

    for (panoc_state, panocplus_state) in Iterators.take(zip(panoc_iter, panocplus_iter), 10)
        @test isapprox(panoc_state.z, panocplus_state.z)
    end
end


text:
env = ENV["AHORN_ENV"]

logfilePath = joinpath(dirname(env), "log-install-ahorn.txt")
println("Logging to " * logfilePath)
logfile = open(logfilePath, "w")

flush(stdout)
flush(stderr)

stdoutReal = stdout
(rd, wr) = redirect_stdout()
redirect_stderr(stdout)

@async while !eof(rd)
    data = String(readavailable(rd))
    print(stdoutReal, data)
    flush(stdoutReal)
    print(logfile, data)
    flush(logfile)
end

using Pkg
Pkg.activate(env)

install_or_update(url::String, pkg::String) = if "Ahorn" ∈ keys(Pkg.Types.Context().env.project.deps)
    println("Updating $pkg...")
    Pkg.update(pkg)
else
    println("Adding $pkg...")
    Pkg.add(PackageSpec(url = url))
end

Pkg.instantiate()

install_or_update("https://github.com/CelestialCartographers/Maple.git", "Maple")
install_or_update("https://github.com/CelestialCartographers/Ahorn.git", "Ahorn")

Pkg.instantiate()
Pkg.API.precompile()

import Ahorn

ext: jl


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

"Backtrace on hidden Markov model object."
function backTrace(self::HMM)
  state = 0
  min = self.tb[1][end]

  for ι ∈ eachindex(self.tb)
    if self.tb[ι][end] <= min
      state = ι
      min = self.tb[ι][end]
    end

    tb = fill(-1, size(self.tb[1], 1))

    for ι ∈ reverse(eachindex(self.tb[1]))
      vlocal = min
      global sig = -1

      for ε ∈ eachindex(self.tb)
        dob = self.tb[ε][ι]
        if dob <= vlocal
          vlocal = dob
          sig = ε
        end
      end

      tb[ι] = state
      state = sig

    end

    tb = tb[1:end - 1]

    return tb
  end
end

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

max_stars_repo_name: DanielRivasMD/HiddenMarkovModelReaders

max_issues_repo_name: DanielRivasMD/HiddenMarkovModelReaders

max_forks_repo_name: DanielRivasMD/HiddenMarkovModelReaders


text:
"""

```
load(filename::String)
```

Function to load the json saved results file into the `Result` data structure.

"""
function load(filename::String)
  f = open(filename,"r")
  R = JSON3.read(f,MutableResult)
  # Need to reshape the solute and solvent atom contributions, because
  # the data is read in a single column
  solute_atom = reshape(R.solute_atom,R.nbins,:)
  solvent_atom = reshape(R.solvent_atom,R.nbins,:)
  R.solute_atom = solute_atom
  R.solvent_atom = solvent_atom
  return Result([getfield(R,field) for field in fieldnames(Result)]...)
end


text:
function create_diffusion_constants(grn::NamedTuple,grn_config::FGRNConfig,exp_config::ExperimentalConfig)

    total_proteins = size(grn.rec_coding,3) + size(grn.beh_coding,3) + size(grn.reg_coding,3) + size(grn.env_coding,3) #size(grn.rec_coding,3) +

    diffusion_constants = zeros(Float64,total_proteins)

    for k in 1:size(grn.env_coding,3)
        if exp_config.diffusion_alpha == "evolve"
            diffusion_constants[grn.env_id[k]] = 0.25*calculate_affinity(grn.env_coding[:,:,k],grn_config.substrate_protein)
        else
            diffusion_constants[grn.env_id[k]] = 0.2*exp_config.diffusion_alpha
        end
    end

    return diffusion_constants
end

function protein_dynamics!(environment::Array{Float64,3},est_setup::ESTSetup,
    grn::NamedTuple,all_cells::Array{Float64,1},target_cell_areas::Dict{Float64,Float64},
    diffusion_constants::Array{Float64,1},Ax::Array{Float64,2},Ay::Array{Float64,2},
    dc::Array{Float64,1},exp_config::ExperimentalConfig,
    cell_config::CellularConfig,grn_config::FGRNConfig)

    # diffusion

    for i in grn.env_id
        A = @view environment[:,:,2+i]
        mul!(Ax,cell_config.Mx,A)
        mul!(Ay,A,cell_config.My)
        @. A = A +  diffusion_constants[i]*(Ax + Ay)
        A[A .< 1e-9] .= 0;
    end

    # reset receptor concentrations for external signal

    environment[:,:,2 + grn.rec_id[1]] .= 0
    environment[:,:,2 + grn.rec_id[2]] .= 0
    environment[:,:,2 + grn.rec_id[3]] .= 0

    if size(grn.rec_id,1) >= 4
        environment[:,:,2 + grn.rec_id[4]] .= 0
    end

    for cell in all_cells

        cell_size = 0.
        which_cells::Array{CartesianIndex{2},1} = []

        for j in 1:size(environment,2)
            for i in 1:size(environment,1)
                if environment[i,j,1] == cell
                    cell_size += 1.
                    push!(which_cells,CartesianIndex(i, j))
                end
            end
        end

        if exp_config.est == "location"
            for index in which_cells
                if environment[index,2] == 1.
                    if (index[2] <= est_setup.left_segment)
                        environment[index,2 + grn.rec_id[1]] = (1 - est_setup.affinity_1)/cell_size
                    else
                        environment[index,2 + grn.rec_id[1]] = 0
                    end

                elseif environment[index,2] == 2.
                    if (index[2] > est_setup.left_segment) & ((index[2] <= est_setup.middle_segment))
                        environment[index,2 + grn.rec_id[2]] = (1 - est_setup.affinity_2)/cell_size
                    else
                        environment[index,2 + grn.rec_id[2]] = 0
                    end

                elseif environment[index,2] == 3.
                    if (index[2] > est_setup.middle_segment)
                        environment[index,2 + grn.rec_id[3]] = (1 - est_setup.affinity_3)/cell_size
                    else
                        environment[index,2 + grn.rec_id[3]] = 0
                    end

                elseif environment[index,2] == 4.
                    if size(grn.rec_id,1) >= 4
                        environment[index,2 + grn.rec_id[4]] = (1 - est_setup.affinity_4)/cell_size
                    end

                else
                    continue
                end
            end

        else
            for index in which_cells
                if environment[index,2] == 1.
                    environment[index,2 + grn.rec_id[1]] = (1 - est_setup.affinity_1)/cell_size
                    environment[index,2 + grn.rec_id[2]] = 0.
                    environment[index,2 + grn.rec_id[3]] = 0.
                    environment[index,2 + grn.rec_id[4]] = 0.
                elseif environment[index,2] == 2.
                    environment[index,2 + grn.rec_id[1]] = 0.
                    environment[index,2 + grn.rec_id[2]] = (1 - est_setup.affinity_2)/cell_size
                    environment[index,2 + grn.rec_id[3]] = 0.
                    environment[index,2 + grn.rec_id[4]] = 0.
                elseif environment[index,2] == 3.
                    environment[index,2 + grn.rec_id[1]] = 0.
                    environment[index,2 + grn.rec_id[2]] = 0.
                    environment[index,2 + grn.rec_id[3]] = (1 - est_setup.affinity_3)/cell_size
                    environment[index,2 + grn.rec_id[4]] = 0.
                elseif environment[index,2] == 4.
                    environment[index,2 + grn.rec_id[1]] = 0.
                    environment[index,2 + grn.rec_id[2]] = 0.
                    environment[index,2 + grn.rec_id[3]] = 0.
                    environment[index,2 + grn.rec_id[4]] = (1 - est_setup.affinity_4)/cell_size
                else
                    environment[index,2 + grn.rec_id[1]] = 0.
                    environment[index,2 + grn.rec_id[2]] = 0.
                    environment[index,2 + grn.rec_id[3]] = 0.
                    environment[index,2 + grn.rec_id[4]] = 0.
                end
            end
        end

        combined_cell_state = sum(environment[which_cells,:],dims = 1)[3:end]

        updateConcentrations!(dc,grn,combined_cell_state,grn_config)

        target_cell_areas[cell] = 2*cell_config.cell_size_upper_limit*max(0.25,combined_cell_state[grn.beh_id[1]]) # growth factor
        type_id = Float64(parse(Int,join(string.(Int.(combined_cell_state[grn.beh_id[2:3]] .> grn.beh_cod_thresh[2:3]))),base=2)) # determine type based on beh conc.

        for j in which_cells
            environment[j,3:size(environment,3)] += (dc ./ cell_size)
            environment[j,2] = type_id + 1.
        end

    end
end

max_stars_repo_name: harrybooth/FractalDevelopment

max_issues_repo_name: harrybooth/FractalDevelopment

max_forks_repo_name: harrybooth/FractalDevelopment


text:
function all(r)
    dir = []

    function next(side, kol)
        return HorizonSide((Int(side) + kol) % 4)
    end

    function set(side1, side2)
        while (!isborder(r, side1) || !isborder(r, side2))
            while !isborder(r, side1)
                move!(r, side1)
                push!(dir, side1)
            end
            while !isborder(r, side2)
                move!(r, side2)
                push!(dir, side2)
            end
        end
    end

    function goAround(side)
        num = 0
        kol = 0
        while (isborder(r, side) && !isborder(r, next(side, 1)))
            move!(r, next(side, 1))
            push!(dir, next(side, 1))
            kol += 1
        end
        if !isborder(r, side)
            move!(r, side)
            push!(dir, side)
            num += 1
            while isborder(r, next(side, 3))
                move!(r, side)
                push!(dir, side)
                num += 1
            end
        end
        for i = 1:kol
            move!(r, next(side, 3))
            push!(dir, next(side, 3))
        end
        return num
    end

    function drawLine(side)
        num = 0
        while true
            if !isborder(r, side)
                putmarker!(r)
                move!(r, side)
                push!(dir, side)
                num += 1
            else
                putmarker!(r)
                num += goAround(side)
            end
            if (num == 11 || (isborder(r, side) && isborder(r, next(side, 1))))
                break
            end
        end
        putmarker!(r)
    end

    function f()
        set(Sud, West)

        side = Ost
        while true
            drawLine(side)
            if isborder(r, Nord)
                break
            end
            move!(r, Nord)
            push!(dir, Nord)
            side = next(side, 2)
        end
        while length(dir) > 0
            move!(r, next(pop!(dir), 2))
        end
    end

    f()
end

max_stars_repo_name: dmitrii-a-piskunov/mirea-progs

max_issues_repo_name: dmitrii-a-piskunov/mirea-progs

max_forks_repo_name: dmitrii-a-piskunov/mirea-progs


text:
import MPI
using MPIClusterManagers
using Dagger
using CairoMakie

using Diffusion

nx=64          # Number of grid points in x dimension
ny=64          # Number of grid points in y dimension
timesteps=30   # Number of time steps

# Notes:
# Producers can run ahead of consumer, by 2000 time steps.
# Makie has high initial latency

function process_video(dims, nx, ny, timesteps)
    nx_v = (nx-2)*dims[1]
    ny_v = (ny-2)*dims[2]
    T    = zeros(nx_v, ny_v)
    node = Observable(T)
    fig  = heatmap(node, colorrange = (0.0, 2.0))

    record(fig, "output.mkv", 1:timesteps) do t
        if t % 2 == 0
            @info "Processing frame" t
            # Process frames from workers
            results = fetch.([Dagger.@spawn single=w Diffusion.get_result() for w in workers()])
            for (coords, result) in results
                cart_x = (1:(nx-2)) .+ (coords[1] * (nx-2))
                cart_y = (1:(ny-2)) .+ (coords[2] * (ny-2))
                T[cart_x,cart_y] .= result
            end
            node[] = T
        end
    end
end

manager = MPIClusterManagers.start_main_loop(MPI_TRANSPORT_ALL) # does not return on worker
try
    @mpi_do manager begin
        using MPI, Distributed
        wcomm = MPI.COMM_WORLD
        id = myid()
        rank = id-1
        temp_comm = MPI.Comm_dup(wcomm)
        gridcomm = MPI.Comm_split(temp_comm, rank == 0 ? MPI.MPI_UNDEFINED[] : 1, rank)
        if id != 1
            @assert MPI.Comm_size(gridcomm) == MPI.Comm_size(wcomm)-1
            gridsize = MPI.Comm_size(gridcomm)
        end

        # Create Cartesian process topology
        if id == 1
            gridsize = MPI.Comm_size(wcomm)-1
            if gridsize == 0
                println("No other processes available, start this program with at least 2 processes")
                exit(1)
            end
            print("Size: $gridsize\n")

            dims = [0,0]
            MPI.Dims_create!(gridsize, dims)
            try
                process_video(dims, nx, ny, timesteps)
            finally
                @info "Terminating workers"
                for w in workers()
                    Dagger.@spawn single=w Diffusion.terminate()
                end
                @info "Terminated workers"
            end
        else
            nprocs = MPI.Comm_size(gridcomm)
            dims = [0,0]
            MPI.Dims_create!(gridsize, dims)
            comm_cart = MPI.Cart_create(gridcomm, dims, #=periods=#[0,0], #=reorder=#0)
            coords = MPI.Cart_coords(comm_cart)

            # Initialize the heat capacity field
            T = ones(Float64, nx, ny)

            # Simulate heat diffusion
            print("Starting diffusion on rank $coords\n")
            Diffusion.diffusion(T, timesteps, comm_cart)
        end
    end
finally
    # Exit gracefully
    MPIClusterManagers.stop_main_loop(manager)
end

max_stars_repo_name: JuliaLabs/DaggerDiffusion

max_issues_repo_name: JuliaLabs/DaggerDiffusion

max_forks_repo_name: JuliaLabs/DaggerDiffusion


text:
"""
    get_atoms(dimensions, N, rₘᵢₙ; kwargs...)

`dimensions` : 1D/2D/3D == `1`/`2`/`3`
`N` : number of atoms
`rₘᵢₙ` : radius of exclusion
"""
function get_atoms(dimensions, N, rₘᵢₙ; kwargs...)
    get_single_atom = kwargs[:createFunction]
    r = get_single_atom(; kwargs...)
    if N == 1
        return r
    end
    temp_atom = zeros(eltype(r), dimensions)
    nValid = 1
    for iLoop = 1:1_000_000
        temp_atom = get_single_atom(; kwargs...)

        if is_atom_valid(temp_atom, r, rₘᵢₙ)
            nValid = nValid + 1
            r = hcat(r, temp_atom)
        end
        if nValid == N
            break
        end
    end
    if nValid < N
        @error "Could not generate all data after 10^6 interactions"
    end

    return r
end

"""
    Checks if the `new atom` has a `minimum distance` for `all the other atoms`
"""
function is_atom_valid(new_atom, r, rₘᵢₙ)
    A = r
    b = new_atom
    allDistances = get_Distance_A_to_b(A, b)
    return all(allDistances .≥ rₘᵢₙ)
end

function get_Distance_A_to_b(A, b)
    N = size(A, 2)
    distanceAb = zeros(N)
    i = 1
    for c in eachcol(A)
        distanceAb[i] = Distances.evaluate(Euclidean(), c, b)
        i += 1
    end
    return distanceAb
end

"""
    get_rₘᵢₙ(ρ)

AD HOC value that changes many sometimes to fix different problems
"""
get_rₘᵢₙ(ρ) = (ρ^(-1 / 3)) / (0.8π^2)

"""
    get_pairwise_matrix(r)
r = matrix with `N` atoms positions

returns a `NxN` Float64 matrix, with zeros at diagonal
"""
function get_pairwise_matrix(r)
    r_matrix = transpose(r)
    R_jk = Distances.pairwise(Euclidean(), r_matrix, r_matrix; dims = 1)
    R_jk[diagind(R_jk)] .= 0
    return R_jk
end

### -----SHAPE CONSTRUCTOR GIVEN MATRIX -------
function Atom(geometry::T, r::Matrix, kR::Union{Real,Integer}) where {T<:Dimension}
    N = size(r, 2) # remember to use each collum as a atom position
    return Atom(geometry, r, N, Float64(kR))
end

### --------------- CONVERSIONS ---------------
"""
    cube_inputs(N::Integer, ρ::Real)
"""
function cube_inputs(N::Integer, ρ::Real)
    kL = (N / ρ)^(1 / 3)
    return (N, kL)
end
"""
    sphere_inputs(N::Integer, ρ::Real)
"""
function sphere_inputs(N::Integer, ρ::Real)
    kR = (N / ((4π / 3) * ρ))^(1 / 3)
    return (N, kR)
end
"""
    cylinder_inputs(N::Integer, ρ::Real; R::Real, h::Real)

- if the use does not specify, R or h, than, it is assumed that both should be equal
- if use only specify R or h, the other variable will change to match the density required
"""
function cylinder_inputs(N::Integer, ρ::Real; kwargs...)
    R =  get(kwargs, :R, NaN)
    h =  get(kwargs, :h, NaN)

    if (R < 0) && (h < 0)
        @error "Invalid inputs, R and h must be positives"
        R = 0
        h = 0
    elseif isnan(R) && isnan(h) # R and h not provided
        # assume that R=h ⇢ (h*πR^2 = πR^3) and solve for R
        R = cbrt(N / π * ρ)
        h = R
    elseif isnan(R) && (h ≠ -1.0) # R not provided
        R = sqrt(N / (π * ρ * h))
    elseif (R ≠ -1.0) && isnan(h) # h not provided
        h = sqrt(N / (π * ρ * R^2))
    end
    return (N, R, h)
end


text:
using CUDA
using Flux
using BenchmarkTools

if CUDA.functional()
    CUDA.allowscalar(false)
    _device = Flux.gpu
    rgen! = CUDA.randn!
end
# testing reflection on batchsize
d = 100
batch_size = 10000
y0 = CUDA.zeros(d,batch_size)
y1 = CUDA.randn(size(y0)...)
@btime _reflect_GPU2($y0,$y1,-1f0,1f0,_device)

@btime CUDA.similar(y0)

function imin2array()
    y1 = CUDA.randn(1000,1000)
    imin = argmin(y1,dims=1) |> Array
    n = zeros(size(y1))
    n[imin] .= 1
    n = b |> gpu
end
@btime imin2array()
@btime CUDA.zeros(size(y1))

n = CUDA.zeros(size(y1))
CUDA.allowscalar(true)
function imin_scalar()
    y1 = CUDA.randn(1000,1000)
    imin = argmin(y1,dims=1)
    n .= 0.
    n[imin] .= 1
end
@btime imin_scalar()


text:
unique_siteind(A::AbstractMPS, B::AbstractMPS, j::Integer) = siteinds(uniqueind, A, B, j)

unique_siteinds(A::AbstractMPS, B::AbstractMPS, j::Integer) = siteinds(uniqueinds, A, B, j)

unique_siteinds(A::AbstractMPS, B::AbstractMPS) = siteinds(uniqueind, A, B)

common_siteind(A::AbstractMPS, B::AbstractMPS, j::Integer) = siteinds(commonind, A, B, j)

common_siteinds(A::AbstractMPS, B::AbstractMPS, j::Integer) = siteinds(commoninds, A, B, j)

common_siteinds(A::AbstractMPS, B::AbstractMPS) = siteinds(commoninds, A, B)

firstsiteind(M::AbstractMPS, j::Integer; kwargs...) = siteind(first, M, j; kwargs...)

map_linkinds!(f::Function, M::AbstractMPS) = map!(f, linkinds, M)

map_linkinds(f::Function, M::AbstractMPS) = map(f, linkinds, M)

function map_common_siteinds!(f::Function, M1::AbstractMPS, M2::AbstractMPS)
  return map!(f, siteinds, commoninds, M1, M2)
end

function map_common_siteinds(f::Function, M1::AbstractMPS, M2::AbstractMPS)
  return map(f, siteinds, commoninds, M1, M2)
end

function map_unique_siteinds!(f::Function, M1::AbstractMPS, M2::AbstractMPS)
  return map!(f, siteinds, uniqueinds, M1, M2)
end

function map_unique_siteinds(f::Function, M1::AbstractMPS, M2::AbstractMPS)
  return map(f, siteinds, uniqueinds, M1, M2)
end

for fname in
    (:sim, :prime, :setprime, :noprime, :addtags, :removetags, :replacetags, :settags)
  @eval begin
    function $(Symbol(fname, :_linkinds))(M::AbstractMPS, args...; kwargs...)
      return map(i -> $fname(i, args...; kwargs...), linkinds, M)
    end
    function $(Symbol(fname, :_linkinds!))(M::AbstractMPS, args...; kwargs...)
      return map!(i -> $fname(i, args...; kwargs...), linkinds, M)
    end
    function $(Symbol(fname, :_common_siteinds))(
      M1::AbstractMPS, M2::AbstractMPS, args...; kwargs...
    )
      return map(i -> $fname(i, args...; kwargs...), siteinds, commoninds, M1, M2)
    end
    function $(Symbol(fname, :_common_siteinds!))(
      M1::AbstractMPS, M2::AbstractMPS, args...; kwargs...
    )
      return map!(i -> $fname(i, args...; kwargs...), siteinds, commoninds, M1, M2)
    end
    function $(Symbol(fname, :_unique_siteinds))(
      M1::AbstractMPS, M2::AbstractMPS, args...; kwargs...
    )
      return map(i -> $fname(i, args...; kwargs...), siteinds, uniqueinds, M1, M2)
    end
    function $(Symbol(fname, :_unique_siteinds!))(
      M1::AbstractMPS, M2::AbstractMPS, args...; kwargs...
    )
      return map!(i -> $fname(i, args...; kwargs...), siteinds, uniqueinds, M1, M2)
    end
  end
end


text:
# Range: 1 : 100,000
# Time:  0.540s

# Written by Alex Vear in 2019
# Public domain.  No rights reserved.

startval = 1
endval = 100000

sqrs = [x^2 for x in startval:endval]
tris = [x*(x+1)/2 for x in startval:endval]
println(intersect(sqrs, tris))

max_stars_repo_name: axvr/codedump

max_issues_repo_name: axvr/codedump

max_forks_repo_name: axvr/codedump


text:
module PProfFlameGraphsTest

using PProf

using Test
using Profile
using ProtoBuf

# Test interactions with FlameGraphs package
using FlameGraphs

const out = tempname() * ".pb.gz"

function foo(n, a, out=[])
    # make this expensive to ensure it's sampled
    for i in 1:n
        push!(out, i*a)
    end
end

@testset "empty output" begin
    Profile.clear()
    # Function doesn't error if no Profile recorded
    fg = flamegraph()
    rm(out, force=true)
    @test (pprof(fg, out=out, web=false); true)
    @test ispath(out)
end

@testset "export basic profile" begin
    Profile.clear()

    let x = 1
        @profile for _ in 1:10000; x += 1; end
        sleep(2)
    end

    # Collect the profile via FlameGraphs
    fg = flamegraph()

    # Write the profile
    outf = pprof(fg, out=out, web=false)

    # Read the exported profile
    fg_prof = open(io->readproto(io, PProf.perftools.profiles.Profile()), outf, "r")

    # Verify that we exported stack trace samples:
    @test length(fg_prof.sample) > 0
    # Verify that we exported frame information
    @test length(fg_prof.location) > 0
    @test length(fg_prof._function) > 0

end

function load_prof_proto(file)
    @show file
    open(io->readproto(io, PProf.perftools.profiles.Profile()), file, "r")
end

@testset "with_c" begin
    Profile.clear()

    let arr = []
        @profile foo(1000000, 5, arr)
        sleep(2)
    end

    fg = flamegraph(C=true)

    # Write a profile that includes C function frames
    with_c = load_prof_proto(pprof(fg, out=tempname(), web=false, from_c = true))

    # Write a profile that excludes C function frames
    without_c = load_prof_proto(pprof(fg, out=tempname(), web=false, from_c = false))

    # Test that C frames were excluded
    @test length(with_c.sample) == length(without_c.sample)
    @test length(with_c.location) > length(without_c.location)
    @test length(with_c._function) > length(without_c._function)
end

@testset "drop_frames/keep_frames" begin
    fg = flamegraph()
    @test load_prof_proto(pprof(fg, out=tempname(), web=false, drop_frames = "foo")).drop_frames != 0
    @test load_prof_proto(pprof(fg, out=tempname(), web=false, keep_frames = "foo")).keep_frames != 0
end

end # module


text:
# This file is part of IntegerTriangles.
# Copyright Peter Luschny. License is MIT.

# Version of: UTC 2021-05-22 22:01:34
# 82286460-bb38-11eb-1729-bfbe096b0d1d

# Do not edit this file, it is generated from the modules and will be overwritten!
# Edit the modules in the src directory and build this file with BuildTriangles.jl!

module runtests
tstdir = realpath(joinpath(dirname(@__FILE__)))
srcdir = joinpath(dirname(tstdir), "src")
tstdir ∉ LOAD_PATH && push!(LOAD_PATH, tstdir)
srcdir ∉ LOAD_PATH && push!(LOAD_PATH, srcdir)
using Nemo, Test, IntegerTriangles, PrettyTables
# *** TrianglesBase.jl *********
# *** TrianglesExamples.jl *********
@testset "Laguerre" begin
A = ℤInt[0, 1, 6, 39, 292, 2505, 24306, 263431]
@test Laguerre(0) == ℤInt[1]
@test Laguerre(3) == ℤInt[6, 18, 9, 1]
@test Laguerre(5) == ℤInt[120, 600, 600, 200, 25, 1]
@test Laguerre(6, 1) == ZZ(4320)
@test [Laguerre(7, k) for k in 0:7] == Laguerre(7)
@test Laguerre(ℤInt[0, 1, 2, 3, 4, 5]) == ZZ(2505)
@test Laguerre(ℤInt[0, 1, 4, 9, 16, 25]) == ZZ(5225)
@test LaguerreTransform([ZZ(n) for n in 0:7]) == A
@test typeof(Laguerre(3, 3))       === fmpz                   # ℤInt
@test typeof(Laguerre(5))          === Array{fmpz,1}          # ℤSeq
@test typeof(Laguerre(A))          === fmpz                   # ℤInt
@test typeof(LaguerreTriangle(5))  === Array{Array{fmpz,1},1} # ℤTri
@test typeof(LaguerreTransform(A)) === Array{fmpz,1}          # ℤSeq
end
# *** TrianglesExplorer.jl *********
TLEN = 19
Explore("SchroederB", "Inv", "AltSum",    TLEN)
# *** TrianglesPlot.jl *********
# *** TrianglesTables.jl *********
header = ["A-number", "Triangle", "Form", "Function", "Sequence"]
mat = [  # Type Matrix{String}, or Array{String,2}
"A000302" "Binomial" "Std" "PolyVal3" "1, 4, 16, 64, 256, 1024, 4096, 16384" ;
]
S = String["A111884" "Lah" "Std" "TransAlts" "1, -1, -1, -1, 1, 19, 151, 1091"]
mat = [mat; S]
L = String["A111884", "Lah", "Std", "TransAlts", "1, -1, -1, -1, 1, 19, 151, 1091"]
mat = [mat; reshape(L, 1, 5)]
L = Explore("Laguerre", "Rev", "TransNat1")
mat = [mat; reshape(L, 1, 5)]
pretty_table(mat, header, alignment=[:l,:l,:l,:l,:l])
mat
# *** TrianglesTraitCard.jl *********
dim = 12
TraitCard(BinomialTriangle, dim)
TraitCard(LaguerreTriangle, dim)
TraitCard(LahTriangle, dim)
TraitCard(CatalanTriangle, dim)
TraitCard(MotzkinTriangle, dim)
TraitCard(NarayanaTriangle, dim)
TraitCard(SchroederBigTriangle, dim)
TraitCard(JacobsthalTriangle, dim)
TraitCard(FibonacciTriangle, dim)
TraitCard(EulerianTriangle, dim)
TraitCard(UniTriangle, dim)
# *** TrianglesUtils.jl *********
T = LaguerreTriangle(7)
println("\nThe $Laguerre triangle in different formats:")
println("\nTriangle in standard format:\n")
Show(T)
println("\nTriangle as a nested array:\n")
Show(T, "nest")
println("\n\nTriangle in flattened format:\n")
Show(T, "flat")
println("\n\nTriangle in info format with sequence ID:\n")
Show(T, "info")
println("\nTriangle in mapped format:\n")
Show(T, "map")
end # module


text:
using Documenter, GPUArrays

function main()
    makedocs(
        modules = [GPUArrays],
        format = Documenter.HTML(
            # Use clean URLs on CI
            prettyurls = get(ENV, "CI", nothing) == "true",
            assets = ["assets/favicon.ico"],
            analytics = "UA-154489943-6",
        ),
        sitename = "GPUArrays.jl",
        pages = [
            "Home"          => "index.md",
            "Interface"     => "interface.md",
            "Functionality" => [
                "functionality/host.md",
                "functionality/device.md",
            ],
            "Test suite"    => "testsuite.md",
        ],
        doctest = true,
    )

    deploydocs(
        repo = "github.com/JuliaGPU/GPUArrays.jl.git"
    )
end

isinteractive() || main()


text:
using Statistics
using LESbrary.Utils

# This test replicates the setup from Souza et al. (2020) and compares
# the answer with Eq. (4). We expect agreement to machine precision.
#
# Souza et al. (2020). Uncertainty quantification of ocean parameterizations: Application to
# the K‐Profile‐Parameterization for penetrative convection. Journal of Advances in Modeling
# Earth Systems, 12, e2020MS002108. DOI: https://doi.org/10.1029/2020MS002108

@testset "diagnose_buoyancy_flux" begin
    α  = 2e-4
    g  = 9.81
    N² = 1e-5
    b(z) = 20*α*g + N²*z

    zs = range(-100, 0, length=64)
    bs = [b(z) for z in zs]

    days = 86400
    Δτ = 8days

    # We want the mixed layer to reach half the domain height.
    depth = mean(zs)

    Qb = diagnose_buoyancy_flux(bs, zs, Δτ, depth)

    h = abs(depth)
    Qb_empirical = h^2  * N² / (2Δτ)

    @test Qb ≈ Qb_empirical
end


text:
function LinearGaussian(μ::Float64, σ::Float64)
    α = 5.0
    x = trace(:x, Normal(μ, σ))
    y = trace(:y, Normal(α * x, 1.0))
    return y
end

function LinearGaussianProposal()
    α = 10.0
    x = trace(:x, Normal(α * 3.0, 3.0))
end

@testset "Simulate" begin
    ret, cl = simulate(LinearGaussian, 0.5, 3.0)
    @test haskey(cl, :x)
    @test haskey(cl, :y)
end

@testset "Generate" begin
    sel = target([(:x, ) => 5.0])
    ret, cl, w = generate(sel, LinearGaussian, 0.5, 3.0)
    @test (cl[:x]) == 5.0
    sel = target([(:x,) => 5.0, (:y,) => 10.0])
    ret, cl, w = generate(sel, LinearGaussian, 0.5, 3.0)
    @test (cl[:x]) == 5.0
    @test (cl[:y]) == 10.0
end

@testset "Update" begin
    ret, cl = simulate(LinearGaussian, 0.5, 3.0)
    score = get_score(cl)
    stored_at_x = (cl[:x])
    stored_at_y = (cl[:y])
    sel = target([(:x,) => 5.0])
    ret, cl, w, retdiff, d = update(sel, cl, Δ(0.5, NoChange()), Δ(3.0, NoChange()))
    @test score ≈ get_score(cl) - w
    score = get_score(cl)
    # Discard should be original :x
    @test (d[:x]) == stored_at_x
    # New should be equal to constraint.
    @test (cl[:x]) == 5.0
    sel = target([(:y,) => 10.0])
    ret, cl, w, retdiff, d = update(sel, cl, Δ(0.5, NoChange()), Δ(3.0, NoChange()))
    @test score ≈ get_score(cl) - w
    @test (d[:y]) == stored_at_y
    @test (cl[:x]) == 5.0
    @test (cl[:y]) == 10.0
end

@testset "Regenerate" begin
    ret, cl = simulate(LinearGaussian, 0.5, 3.0)
    score = get_score(cl)
    stored_at_x = (cl[:x])
    stored_at_y = (cl[:y])
    sel = target([(:x, )])
    ret, cl, w, retdiff, d = regenerate(sel, cl, Δ(0.5, NoChange()), Δ(3.0, NoChange()))
    @test score ≈ get_score(cl) - w
    score = get_score(cl)
    # Discard should be original :x
    @test (d[:x]) == stored_at_x
    # New should not be equal to original.
    @test (cl[:x]) != stored_at_x
    sel = target([(:y, )])
    ret, cl, w, retdiff, d = regenerate(sel, cl, Δ(0.5, NoChange()), Δ(3.0, NoChange()))
    @test score ≈ get_score(cl) - w
    @test (d[:y]) == stored_at_y
end

@testset "Propose" begin
    ret, pmap, w = propose(LinearGaussian, 0.5, 3.0)
    # Propose should track score.
    @test w != 0.0
end

@testset "Assess" begin
    sel = target([(:x, ) => 5.0, (:y, ) => 10.0])
    ret, sc = assess(sel, Jaynes.Empty(), LinearGaussian, 0.5, 3.0)
    @test sc != 0.0
end


text:
# #############################################################################
# #############################################################################
#Mode Plotting functions

#This file contains functions for plotting modes in a cavity.

# #############################################################################
# #############################################################################
#Initiate

#require("util.jl")
#require("boundary.jl")
#require("refractive_index.jl")
#require("io.jl")
#require("ray_solver.jl")
#require("ray_plot.jl")
#require("clustering.jl")

# #############################################################################
# #############################################################################
#Preallocated array based non-interactive PSS plot

#Get corresponding index of a position in an array representing regions in a range
#This is a nearest neighbor interpolation for points within the range
function getbinindex(pos::Float64,min::Float64,max::Float64,res::Int64)
    #Check if outside region
    if (pos < min && pos < max) || (pos > min && pos > max)
        return 0
    else
        #Return the bin number containing pos, if there are res bins between
        #min and max
        return convert(Int64,cld(res*(pos-min),max-min))
    end
end

#Incorporate bouncepoint dataset into arrays, letting each point replace those in
#its corresponding bin if it has a higher priority
function addbouncepointstopss(bouncepoints::Array{Float64,2},resultid::Int64,pssarray::Array{Float64,2},residarray::Array{Int64,2},range::Array{Float64,1})
    #Region parameters
    #Range is [thetamin,thetamax,chimin,chimax]
    #assert(size(pssarray) == size(residarray))

    #Determine score for set of bouncepoints based on degree of cluster spread
    #Each nonnegative cluster spread value is scored with a value in (0.0,1.0), stationary at small
    #distances, and asymptotic to 1 as it approaches infinity. Lower score has higher priority.
    #Cluster spread estimated using variance measures consistent with metric on pss
    order::Int64,clustsize::Float64 = findcluster(bouncepoints,maxtestinterval=25)
    clusterstddev::Float64 = sqrt(var(sin(bouncepoints[1:order:end,1]))
                                  + var(cos(bouncepoints[1:order:end,1]))
                                  + var(map(chi2whisporder,bouncepoints[1:order:end,2])))
    score::Float64 = sq(1.0 - 1.0 / (1.0 + 10*clusterstddev))

    #Loop over all bouncepoints
    for i = 1:size(bouncepoints,1)
        #Get corresponding index in array
        thetaindex::Int64 = getbinindex(bouncepoints[i,1],range[1],range[2],size(pssarray,2))
        chiindex::Int64 = getbinindex(bouncepoints[i,2],range[4],range[3],size(pssarray,1))

        #Skip if out of range
        if thetaindex == 0 || chiindex == 0
            continue
        end

        #Replace value if score is lower
        if score < pssarray[chiindex,thetaindex]
            pssarray[chiindex,thetaindex] = score
            residarray[chiindex,thetaindex] = resultid
        end
    end
end

#Incorporate all bouncepoint datasets in a directory
function adddirtopss(resultsdir::AbstractString;range::Array{Float64,1}=[-pi,pi,-pi/6,pi/2],res::Int64=400)
    #Initialize Arrays
    pssarray::Array{Float64,2} = fill(10.0,res,res) #array of some large number by default
    residarray::Array{Int64,2} = zeros(Int64,res,res)

    #Incorporate data
    resultids::Array{Int64,1} = collectresultids(resultsdir,"bouncepoints")
    for resultid in resultids
        bouncepoints::Array{Float64,2} = readresults(resultsdir,resultid,"bouncepoints")
        addbouncepointstopss(bouncepoints,resultid,pssarray,residarray,range)
    end

    #Return array of pss scores, and array of result ids
    return pssarray,residarray
end

#Plot PSS
function plotpss(resultsdir::AbstractString;axes=plt[:gca](),range::Array{Float64,1}=[-pi,pi,-pi/6,pi/2],res::Int64=400,colorrange::Array{Float64,1}=[-0.1,1.1])
    #Colormap normalization
    colornorm::PyObject = plt[:matplotlib][:colors][:Normalize](colorrange[1],colorrange[2])

    #Get PSS
    pssarray::Array{Float64,2},residarray::Array{Int64,2} = adddirtopss(resultsdir;range=range,res=res)

    #Plot pssarray
    plt[:imshow](pssarray,norm=colornorm,cmap="spectral_r",extent=range,aspect=3)
    plt[:title]("Poincare Surface of Section")
    plt[:xlabel]("Polar Angle \$\\theta\$")
    plt[:ylabel]("Angle of Incidence \$\\chi\$")

    return pssarray,residarray
end

# #############################################################################
# #############################################################################
#PSS plot I/O

#Get PSS plot image file location
function getpssimgfile(resultsdir::AbstractString)
    return joinpath(resultsdir,"pssplot.png")
end

#Write image of PSS plot
function writepssimg(resultsdir::AbstractString;kwargs...)
    #Create new figure
    pssplot = plt[:figure](figsize=(9,9))[:gca](aspect="equal")

    #Plot PSS non-interactively
    plotpss(resultsdir;axes=pssplot,kwargs...)

    #Save image
    plt[:savefig](getpssimgfile(resultsdir))
    return nothing
end


text:
module ImagineFormat

using AxisArrays, ImageMetadata, FileIO, Unitful, FixedPointNumbers, SharedArrays

export imagine2nrrd, BidiImageArray

include("bidi.jl")

const μm = u"μm"
const s = u"s"
const μs = u"μs"
const MHz = u"MHz"

function load(f::File{format"Imagine"}; mode="r")
    fabs = File(format"Imagine", abspath(f.filename))
    open(fabs) do s
        skipmagic(s)
        load(s, mode=mode)
    end
end

function load(io::Stream{format"Imagine"}; mode="r")
    s = stream(io)
    h = parse_header(s)
    filename = s.name[7:end-1]
    basename, ext = splitext(filename)
    camfilename = basename*".cam"
    Traw = h["pixel data type"]
    T = Traw <: Unsigned ? Normed{Traw, h["original image depth"]} : Traw
    havez = h["frames per stack"] > 1
    havet = h["nStacks"] > 1
    sz = havez & havet ? [h["image width"], h["image height"], h["frames per stack"], h["nStacks"]] :
         havez ? [h["image width"], h["image height"], h["frames per stack"]] :
         havet ? [h["image width"], h["image height"], h["nStacks"]] :
         [h["image width"], h["image height"]]
    # Check that the file size is consistent with the expected size
    if !isfile(camfilename)
        @warn("Cannot open ", camfilename)
        data = Array{T}(undef, sz[1], sz[2], sz[3], 0)
    else
        fsz = filesize(camfilename)
        n_stacks = sz[end]
        if fsz != sizeof(T)*prod(map(Int64,sz))  # guard against overflow on 32bit
            @warn("Size of image file is different from expected value")
            n_stacks = floor(Int, fsz / sizeof(T) / prod(sz[1:end-1]))
        end
        if sizeof(T)*prod(map(Int64,sz[1:end-1]))*n_stacks > typemax(UInt)
            @warn("File size is too big to mmap on 32bit")
            n_stacks = floor(Int, fsz / sizeof(T) / typemax(UInt))
        end
        if n_stacks < sz[end]
            println("Truncating to ", n_stacks, length(sz) == 4 ? " stacks" : " frames")
            sz[end] = n_stacks
        end
        data = SharedArray{T}(camfilename, tuple(sz...), mode=mode)
    end
    um_per_pixel = h["um per pixel"]*μm
    pstart = h["piezo"]["stop position"]
    pstop = h["piezo"]["start position"]
    if length(sz)>2
        dz = abs(pstart - pstop)/(sz[3]-1)
    else dz = 0.0 end

    axisnames = havez ? (:x, :l, :z) : (:x, :l)
    pixelspacing = havez ? (um_per_pixel, um_per_pixel, dz) : (um_per_pixel, um_per_pixel)
    if havet
        axisnames = (axisnames..., :time)
        pixelspacing = (pixelspacing..., h["idle time between stacks"]+h["frames per stack"]*h["exposure time"])
    end
    axs = map((n,s,l)->Axis{n}(range(0*s, stop=(l-1)*s, length=l)), axisnames, pixelspacing, size(data))
    is_bidi = get(h, "bidirectional", false) || get(h, "bidirection", false)
    if is_bidi
        if !havez
            @warn("Image was marked as bidirectional but has no z axis.  Ignoring bidirectional tag")
        elseif !havet
            @warn("Image was marked as bidirectional but is not a timeseries.  Ignoring bidirectional tag")
        else
            data = BidiImageArray(data)
        end
    end
    ImageMeta(AxisArray(data, axisnames, pixelspacing), imagineheader=h, suppress=Set(Any["imagineheader"]))
end

abstract type Endian end
struct LittleEndian <: Endian; end
struct BigEndian <: Endian; end
const endian_dict = Dict("l"=>LittleEndian(), "b"=>BigEndian())
const nrrd_endian_dict = Dict(LittleEndian()=>"little",BigEndian()=>"big")
parse_endian(s::AbstractString) = endian_dict[lowercase(s)]

function parse_vector_int(s::AbstractString)
    ss = split(s, r"[ ,;]", keepempty=false)
    v = Vector{Int}(undef, length(ss))
    for i = 1:length(ss)
        v[i] = parse(Int,ss[i])
    end
    return v
end

const bitname_dict = Dict(
  "int8"      => Int8,
  "uint8"     => UInt8,
  "int16"     => Int16,
  "uint16"    => UInt16,
  "int32"     => Int32,
  "uint32"    => UInt32,
  "int64"     => Int64,
  "uint64"    => UInt64,
  "float16"   => Float16,
  "float32"   => Float32,
  "single"    => Float32,
  "float64"   => Float64,
  "double"    => Float64)

parse_bittypename(s::AbstractString) = bitname_dict[lowercase(s)]

function float64_or_empty(s::AbstractString)
    if isempty(s)
        return NaN
    else
        return parse(Float64,s)
    end
end

function parse_quantity_or_empty(s::AbstractString)
    if isempty(s)
        return NaN
    else
        return parse_quantity(s)
    end
end

_unit_string_dict = Dict("um" => μm, "s" => s, "us" => μs, "MHz" => MHz)
function parse_quantity(s::AbstractString, strict::Bool = true)
    if s == "NA"
        return Float64(NaN)
    end
    # Find the last character of the numeric component
    m = match(r"[0-9\.\+-](?![0-9\.\+-])", s)
    if m == nothing
        error("AbstractString does not have a 'value unit' structure")
    end
    val = parse(Float64, s[1:m.offset])
    ustr = strip(s[m.offset+1:end])
    if isempty(ustr)
        if strict
            error("AbstractString does not have a 'value unit' structure")
        else
            return val
        end
    end
    val * _unit_string_dict[ustr]
end

# Read and parse a *.imagine file (an Imagine header file)
const compound_fields = Any["piezo", "binning"]
const field_key_dict = Dict{AbstractString,Function}(
    "header version"               => x->parse(Float64,x),
    "app version"                  => identity,
    "date and time"                => identity,
    "rig"                          => identity,
    "byte order"                   => parse_endian,
    "stimulus file content"        => identity,  # stimulus info parsed separately
    "comment"                      => identity,
    "ai data file"                 => identity,
    "image data file"              => identity,
    "di data file"                 => identity,
    "command file"                 => identity,
    "start position"               => parse_quantity,
    "stop position"                => parse_quantity,
    "bidirection"                  => x->parse(Int,x) != 0,
    "bidirectional"                => x->parse(Int,x) != 0,
    "405nm on"                     => x->parse(Int,x) != 0,
    "445nm on"                     => x->parse(Int,x) != 0,
    "488nm on"                     => x->parse(Int,x) != 0,
    "514nm on"                     => x->parse(Int,x) != 0,
    "561nm on"                     => x->parse(Int,x) != 0,
    "405nm percent intensity"      => x->parse(Float64,x),
    "445nm percent intensity"      => x->parse(Float64,x),
    "488nm percent intensity"      => x->parse(Float64,x),
    "514nm percent intensity"      => x->parse(Float64,x),
    "561nm percent intensity"      => x->parse(Float64,x),
    "output scan rate"             => x->parse_quantity(x, false),
    "nscans"                       => x->parse(Int,x),
    "di nscans"                    => x->parse(Int,x),
    "channel list"                 => parse_vector_int,
    "di channel list"              => parse_vector_int,
    "label list"                   => identity,
    "di label list"                => identity,
    "scan rate"                    => x->parse_quantity(x, false),
    "di scan rate"                 => x->parse_quantity(x, false),
    "min sample"                   => x->parse(Int,x),
    "max sample"                   => x->parse(Int,x),
    "min input"                    => x->parse(Float64,x),
    "max input"                    => x->parse(Float64,x),
    "original image depth"         => x->parse(Int,x),
    "saved image depth"            => x->parse(Int,x),
    "image width"                  => x->parse(Int,x),
    "image height"                 => x->parse(Int,x),
    "number of frames requested"   => x->parse(Int,x),
    "nStacks"                      => x->parse(Int,x),
    "idle time between stacks"     => parse_quantity,
    "pre amp gain"                 => float64_or_empty,
    "EM gain"                      => float64_or_empty,
    "gain"                         => float64_or_empty,
    "exposure time"                => parse_quantity,
    "vertical shift speed"         => parse_quantity_or_empty,
    "vertical clock vol amp"       => x->parse(Float64,x),
    "readout rate"                 => parse_quantity_or_empty,
    "pixel order"                  => identity,
    "frame index offset"           => x->parse(Int,x),
    "frames per stack"             => x->parse(Int,x),
    "pixel data type"              => parse_bittypename,
    "camera"                       => identity,
    "um per pixel"                 => x->parse(Float64,x),
    "hbin"                         => x->parse(Int,x),
    "vbin"                         => x->parse(Int,x),
    "hstart"                       => x->parse(Int,x),
    "hend"                         => x->parse(Int,x),
    "vstart"                       => x->parse(Int,x),
    "vend"                         => x->parse(Int,x),
    "angle from horizontal (deg)"  => float64_or_empty,
    "x translation in pixels"      => x->parse(Int,x) != 0,
    "y translation in pixels"      => x->parse(Int,x) != 0,
    "rotation angle in degree"     => x->parse(Float64,x) != 0,
    "stimulus scan hi"             => parse_vector_int,
    "stimulus scan lo"             => parse_vector_int,
    "stimulus frame hi"            => parse_vector_int,
    "stimulus frame lo"            => parse_vector_int,
    "stimulus sequence"            => identity)

function parse_header(s::IOStream)
    headerdict = Dict{String, Any}()
    for this_line = eachline(s)
        this_line = strip(this_line)
        if !isempty(this_line) && !occursin(r"\[.*\]", this_line)
            # Split on =
            m = match(r"=", this_line)
            if m.offset < 2
                error("Line does not contain =")
            end
            k = this_line[1:m.offset-1]
            v = this_line[m.offset+1:end]
            if in(k, compound_fields)
                thisdict = Dict{String, Any}()
                # Split on semicolon
                strs = split(v, r";")
                for i = 1:length(strs)
                    substrs = split(strs[i], r":")
                    @assert length(substrs) == 2
                    k2 = strip(substrs[1])
                    func = field_key_dict[k2]
                    v2 = strip(substrs[2])
                    try
                        thisdict[k2] = func(v2)
                    catch err
                        println("Error processing key '", k2, "' with value ", v2)
                        rethrow(err)
                    end
                end
                headerdict[k] = thisdict
            else
                func = get(field_key_dict, k, identity)
                try
                    headerdict[k] = func(v)
                catch err
                    println("Error processing key ", k, " with value ", v)
                    rethrow(err)
                end
            end
        end
    end
    return headerdict
end
function parse_header(io::Stream{format"Imagine"})
    skipmagic(io)
    parse_header(stream(io))
end
function parse_header(f::File{format"Imagine"})
    open(f) do io
        parse_header(io)
    end
end
parse_header(filename::AbstractString) = parse_header(query(filename))

function imagine2nrrd(sheader::IO, h::Dict{String, Any}, datafilename = nothing)
    println(sheader, "NRRD0001")
    T = h["pixel data type"]
    if T<:AbstractFloat
        println(sheader, "type: ", (T == Float32) ? "float" : "double")
    else
        println(sheader, "type: ", lowercase(string(T)))
    end
    sz = [h["image width"], h["image height"], h["frames per stack"], h["nStacks"]]
    kinds = ["space", "space", "space", "time"]
    if sz[end] == 1
        sz = sz[1:3]
        kinds = kinds[[1,2,4]]
    end
    println(sheader, "dimension: ", length(sz))
    print(sheader, "sizes:")
    for z in sz
        print(sheader, " ", z)
    end
    print(sheader, "\nkinds:")
    for k in kinds
        print(sheader, " ", k)
    end
    print(sheader, "\n")
    println(sheader, "encoding: raw")
    println(sheader, "endian: ", nrrd_endian_dict[h["byte order"]])
    if isa(datafilename, AbstractString)
        println(sheader, "data file: ", datafilename)
    end
    sheader
end

function imagine2nrrd(nrrdname::AbstractString, h::Dict{String, Any}, datafilename = nothing)
    sheader = open(nrrdname, "w")
    imagine2nrrd(sheader, h, datafilename)
    close(sheader)
end

"""
    save_header(filename, header; section=())

Write a `header` dictionary in Imagine format.
`section` can be one of `general`, `misc`, `ai`, or `camera`, to write additional entries
in the specified section. For example,

    header["stimulus info"] = "here is some info"
    Imagine.save_header(filename, header; misc=("stimulus info",))

would allow one to insert "stimulus info" into the "[misc params]" section of the `.imagine`
file.
"""
function save_header(filename::AbstractString, h::Dict{String};
                     general=(),
                     misc=(),
                     ai=(),
                     camera=())
    open(filename, "w") do io
        write(io, magic(format"Imagine"))
        println(io, "\n[general]")
        writekv(io, h, ("header version", "app version", "date and time", "byte order", "rig", general...))
        println(io, "\n[misc params]")
        writekv(io, h, ("stimulus file content", "comment", "ai data file","image data file", "piezo", misc...))
        println(io, "\n[ai]")
        writekv(io, h, ("nscans", "channel list", "label list", "scan rate", "min sample", "max sample", "min input", "max input", ai...))
        println(io, "\n[camera]")
        writekv(io, h, ("original image depth", "saved image depth", "image width",
                        "image height", "number of frames requested", "nStacks",
                        "idle time between stacks", "pre amp gain", "gain", "exposure time",
                        "vertical shift speed", "vertical clock vol amp", "readout rate",
                        "pixel order", "frame index offset", "frames per stack",
                        "pixel data type", "camera", "um per pixel", "binning",
                        "angle from horizontal (deg)", "bidirectional", camera...))
    end
    nothing
end

"""
    save_header(destname, srcname, img::AbstractArray, [T::Type = eltype(img)])

Write an `.imagine` file with name `destname`, using
the `.imagine` file `srcname` as a template. Size and element type
fields are updated from `img` and `T`, respectively.
"""
function save_header(dest::AbstractString, src::AbstractString, img::AbstractArray, T::Type = eltype(img); kwargs...)
    h = parse_header(src)
    fillsize!(h, img)
    h["pixel data type"] = lowercase(string(T))
    h["byte order"] = ENDIAN_BOM == 0x04030201 ? "l" : "b"
    save_header(dest, h; kwargs...)
end

function fillsize!(h, img::AxisArray)
    h["image width"] = size(img, Axis{:x})
    h["image height"] = size(img, Axis{:l})
    h["frames per stack"] = size(img, Axis{:z})
    h["nStacks"] = size(img, Axis{:time})
    h
end

function fillsize!(h, img::AbstractArray{T,4}) where T
    h["image width"] = size(img,1)
    h["image height"] = size(img,2)
    h["frames per stack"] = size(img,3)
    h["nStacks"] = size(img,4)
    h
end

function writekv(io, h, fieldnames)
    for fn in fieldnames
        if haskey(h, fn)
            writefield(io, fn, h[fn])
        end
    end
end
writefield(io, fn, v) = println(io, fn, "=", v)
function writefield(io, fn, v::Number)
    print(io, fn, "=")
    if haskey(write_dict, fn)
        write_dict[fn](io, v)
    else
        println(io, isnan(v) ? "" : v)
    end
end
writefield(io, fn, v::AbstractVector) = println(io, fn, "=", join(map(string,v)," "))

function writefield(io, fn, dct::Dict)
    print(io, fn, "=")
    ks = collect(keys(dct))
    vs = collect(values(dct))
    for i = 1:length(ks)
        k, v = ks[i], vs[i]
        print(io, k, ':')
        if haskey(write_dict, k)
            write_dict[k](io, v)
        else
            print(io, v)
        end
        print(io, i == length(ks) ? '\n' : ';')
    end
end

writefield(io, fn, ::LittleEndian) = println(io, "byte order=l")
writefield(io, fn, ::BigEndian)    = println(io, "byte order=b")

writeum(io,x) = print(io, x/μm, " um")
writeus(io,x::Unitful.Time) = print(io, x/μs, " us")
writeus(io,x) = nothing
writeMHz(io,x::Unitful.Frequency) = print(io, x/MHz, " MHz")
writeMHz(io,x) = nothing

const write_dict = Dict{String,Function}(
    "bidirectional"                => (io,x)->x ? print(io, 1) : print(io, 0),
    "bidirection"                  => (io,x)->x ? print(io, 1) : print(io, 0),
    "start position"               => writeum,
    "stop position"                => writeum,
    "vertical shift speed"         => (io,x)->(writeus(io,x); print(io,'\n')),
    "readout rate"                 => (io,x)->(writeMHz(io,x); print(io,'\n')),
)

end


text:
module WeylChamber

export D_PE,
    canonical_gate,
    gate_concurrence,
    in_weyl_chamber,
    local_invariants,
    unitarity,
    weyl_chamber_coordinates,
    weyl_chamber_region

using LinearAlgebra

const 𝕚 = 1im

const Q_magic = [
    1  0  0  𝕚
    0  𝕚  1  0
    0  𝕚 -1  0
    1  0  0 -𝕚
]

const SxSx = ComplexF64[
    0  0  0  1
    0  0  1  0
    0  1  0  0
    1  0  0  0
]

const SySy = ComplexF64[
     0  0  0 -1
     0  0  1  0
     0  1  0  0
    -1  0  0  0
]

const SzSz = ComplexF64[
    1  0  0  0
    0 -1  0  0
    0  0 -1  0
    0  0  0  1
]

"""Calculate the local invariants g₁, g₂, g₃ for a two-qubit gate.

```julia
g₁, g₂, g₃ = local_invariants(U)
```
"""
function local_invariants(U)
    Q = Q_magic
    UB = Q' * U * Q  # "to-magic"
    detU = det(UB)
    m = transpose(UB) * UB
    g₁₂ = tr(m)^2 / 16detU
    g₁ = real(g₁₂)
    g₂ = imag(g₁₂)
    g₃ = real((tr(m)^2 - tr(m * m)) / 4detU)
    return g₁, g₂, g₃
end

"""Calculate the Weyl chamber coordinates c₁, c₂, c₃ for a two-qubit gate.

```julia
c₁, c₂, c₃ = weyl_chamber_coordinates(U)
```

calculates the Weyl chamber coordinates using the algorithm described in Childs
et al., PRA 68, 052311 (2003).
"""
function weyl_chamber_coordinates(U)

    @assert size(U) == (4, 4)
    detU = Complex(det(U))
    Ũ = SySy * transpose(U) * SySy
    two_S = [angle(z) / π for z in eigvals(U * Ũ / √detU)]

    # Check whether the argument of the eigenvalues are on correct branch. If
    # not, put them on correct branch: `angle` returns values in (-π,π] whereas
    # we need values on the branch (-π/2,3π/2]. This implies that those
    # arguments which are between -π and -π/2 need to be shifted by 2π
    two_S = [(ϕⱼ + 0.5 < -1e-10 ? ϕⱼ + 2 : ϕⱼ) for ϕⱼ ∈ two_S]

    p = sortperm(two_S, rev=true)  # Zygote can't handle a direct `sort`
    S = [two_S[p[1]] / 2, two_S[p[2]] / 2, two_S[p[3]] / 2, two_S[p[4]] / 2]
    n = Int(round(sum(S)))  # number of ϕⱼ ≤ -π/4
    @assert n ≥ 0
    if n > 0
        S = [j <= n ? ϕⱼ - 1.0 : ϕⱼ for (j, ϕⱼ) in enumerate(S)]
        S = circshift(S, -n)
    end
    c₁ = S[1] + S[2]
    c₂ = S[1] + S[3]
    c₃ = S[2] + S[3]
    if c₃ < 0
        c₁ = 1 - c₁
        c₃ = -c₃
    end
    return c₁, c₂, c₃

end

"""Calculate the maximum gate concurrence.

```julia
C = gate_concurrence(U)
C = gate_concurrence(c₁, c₂, c₃)
```

calculates that maximum concurrence ``C ∈ [0, 1]`` that the two two-qubit gate `U`,
respectively the gate described by the Weyl chamber coordinates `c₁`, `c₂`,
`c₃` (see [`weyl_chamber_coordinates`](@ref)) can generate.

See Kraus, Cirac, Phys. Rev. A 63, 062309 (2001)
"""
function gate_concurrence(c₁, c₂, c₃)
    if (c₁ + c₂ ≥ 0.5) && (c₁ - c₂ ≤ 0.5) && (c₂ + c₃ ≤ 0.5)
        # If we're inside the perfect-entangler polyhedron in the Weyl chamber
        # the concurrence is 1 by definition. The "regular" formula gives wrong
        # results in this case.
        C = 1.0
    else
        # Outside of the polyhedron, the Formula of Eq (21) in PRA 63, 062309
        # is valid
        return max(
            abs(sin(π * (c₁ + c₃))),
            abs(sin(π * (c₂ + c₁))),
            abs(sin(π * (c₃ + c₂))),
            abs(sin(π * (c₁ - c₃))),
            abs(sin(π * (c₂ - c₁))),
            abs(sin(π * (c₃ - c₂))),
        )
    end
end

gate_concurrence(U) = gate_concurrence(weyl_chamber_coordinates(U)...)

@doc raw"""Construct the canonical gate for the given Weyl chamber coordinates.

```julia
Û = canonical_gate(c₁, c₂, c₃)
```

constructs the two qubit gate ``Û`` as

```math
Û = \exp\left[i\frac{π}{2} (c_1 σ̂_x σ̂_x + c_2 σ̂_y σ̂_y + c_3 σ̂_z σ̂_z)\right]
```

where ``σ̂_{x,y,z}`` are the Pauli matrices.
"""
function canonical_gate(c₁, c₂, c₃)
    return exp(𝕚 * π / 2 * (c₁ * SxSx + c₂ * SySy + c₃ * SzSz))
end

"""Identify which region of the Weyl chamber a given gate is located in.

```julia
region = weyl_chamber_region(U)
region = weyl_chamber_region(c₁, c₂, c₃)
```

identifies which region of the Weyl chamber the given two-qubit gate `U`,
respectively the gate identified by the Weyl chamber coordinates `c₁`, `c₂`,
`c₃` (see [`weyl_chamber_coordinates`](@ref)) is in, as a string. Possible
outputs are:

* `"PE"`: gate is in the polyhedron of perfect entanglers.
* `"W0"`: gate is between the identity and the perfect entanglers.
* `"W0*"`: gate is between CPHASE(2π) and the perfect entanglers.
* `"W1"`: gate is between SWAP and the perfect entanglers.

For invalid Weyl chamber coordinates, an empty string is returned.
"""
function weyl_chamber_region(c₁, c₂, c₃)
    in_weyl_chamber = (
        ((0 ≤ c₁ < 0.5) && (0 ≤ c₂ ≤ c₁) && (0 ≤ c₃ ≤ c₂)) ||
        ((0.5 ≤ c₁ ≤ 1.0) && (0 ≤ c₂ ≤ (1 - c₁)) && (0 ≤ c₃ ≤ c₂))
    )
    if in_weyl_chamber
        if ((c₁ + c₂) ≥ 0.5) && ((c₁ - c₂) ≤ 0.5) && ((c₂ + c₃) ≤ 0.5)
            return "PE"
        elseif (c₁ + c₂) < 0.5
            return "W0"
        elseif (c₁ - c₂) > 0.5
            return "W0*"
        elseif (c₂ + c₃) > 0.5
            return "W1"
        else
            throw(ErrorException("Internal Error: ($c₁, $c₂, $c₃) not handled"))
        end
    else
        return ""
    end
end

weyl_chamber_region(U) = weyl_chamber_region(weyl_chamber_coordinates(U)...)

"""Check whether a given gate is in (a specific region of) the Weyl chamber.

```
in_weyl_chamber(c₁, c₂, c₃)
```

checks whether `c₁`, `c₂`, `c₃` are valid Weyl chamber coordinates.

```
in_weyl_chamber(U; region="PE")
in_weyl_chamber(c₁, c₂, c₃; region="PE")
```

checks whether the two-qubit gate `U`, respectively the gate described by the
Weyl chamber coordinates `c₁`, `c₂`, `c₃` (see
[`weyl_chamber_coordinates`](@ref)) is a perfect entangler. The `region` can be
any other of the regions returned by [`weyl_chamber_region`](@ref).
"""
function in_weyl_chamber(c₁, c₂, c₃; region="W")
    regions = ["W", "PE", "W0", "W0*", "W1"]
    if region ∉ regions
        throw(ArgumentError("Invalid region $(repr(region)), must be one of $regions"))
    end
    found_region = weyl_chamber_region(c₁, c₂, c₃)
    if region == "W"
        return found_region ≠ ""
    else
        return found_region == region
    end
end

in_weyl_chamber(U; region) = in_weyl_chamber(weyl_chamber_coordinates(U)...; region)

@doc raw"""Unitarity of a matrix.

```
pop_loss = 1 - unitarity(U)
```

measures the loss of population from the subspace described by
`U`. E.g., for a two-qubit gate, `U` is a 4×4 matrix. The `unitarity` is
defined as ``\Re[\tr(Û^†Û) / N]`` where ``N`` is the dimension of ``Û``.
"""
function unitarity(U)
    N = size(U)[1]
    return real(tr(U' * U) / N)
end

@doc raw"""Perfect-entanglers distance measure.

```julia
D = D_PE(U; unitarity_weight=0.0, absolute_square=false)
```

For a given two-qubit gate ``Û``, this is defined via the
[`local_invariants`](@ref) ``g_1``, ``g_2``, ``g_3`` as

```math
D = g_3 \sqrt{g_1^2 + g_2^2} - g_1
```

This describes the geometric distance of the quantum gate from the polyhedron
of perfect entanglers in the Weyl chamber.

This equation is only meaningful under the assumption that ``Û`` is unitary. If
the two-qubit level are a logical subspace embedded in a larger physical
Hilbert space, loss of population from the logical subspace may lead to a
non-unitary ``Û``. In this case, the [`unitarity`](@ref) measure can be added
to the functional by giving a `unitary_weight` ∈ [0, 1) that specifies the
relative proportion of the ``D`` term and the unitarity term.

By specifying `absolute_square=true`, the functional is modified as ``D →
|D|²``, optimizing specifically for the *boundary* of the perfect entanglers
polyhedron. This accounts for the fact that ``D`` can take negative values
inside the polyhedron, as well as the `W1` region of the Weyl chamber (the one
adjacent to SWAP). This may be especially useful in a system with population
loss (`unitarity_weight` > 0), as it avoids situations where the optimization
goes deeper into the perfect entanglers while *increasing* population loss.

!!! warning

    The functional does not check which region of the Weyl chamber the quantum
    gate is in. When using this for an optimization where the guess leads to a
    point in the `W1` region of the Weyl chamber (close to SWAP), the sign of
    the functional must be flipped, or else it will optimize for SWAP.
    Alternatively, use `absolute_square=true`.

!!! tip

    The functional can be converted into the correct form for an optimization
    that uses one objective for each logical basis state by using
    `QuantumControl.Functionals.gate_functional`.
"""
function D_PE(U; unitarity_weight=0.0, absolute_square=false)
    w::Float64 = clamp(1.0 - unitarity_weight, 0.0001, 1.0)
    N = 4
    @assert size(U) == (N, N)
    g₁, g₂, g₃ = local_invariants(U)
    if absolute_square
        D = abs2(g₃ * sqrt(g₁^2 + g₂^2) - g₁)
    else
        D = g₃ * sqrt(g₁^2 + g₂^2) - g₁
    end
    pop_loss = 1 - unitarity(U)
    return w * D + (1 - w) * pop_loss
end

end


text:
using Plots

include("models/GCM.jl")
include("modules/ClassicalDynamics.jl")

GCMParameters = [-1,0.6,1,1]
energy = 10
numTrajectories = 10
PSPoints = 2000

savePath = "d:\\results\\GCM"

pyplot(size = (1800,1000))
fig = PoincareSection(energy, GCMParameters, numTrajectories; maxPSPoints=PSPoints, verbose=true)
savefig(fig, savePath * "_PS_E=$(energy)_$GCMParameters.png")

max_stars_repo_name: PavelStransky/ClassicalDynamics

max_issues_repo_name: PavelStransky/ClassicalDynamics

max_forks_repo_name: PavelStransky/ClassicalDynamics


text:
#===============================================================================
                            Path imputation
===============================================================================#

function eMCMC.update!(
        update::PathImputation,
        global_ws::DiffusionGlobalWorkspace,
        local_ws::DiffusionLocalWorkspace,
        step,
    )
    ws, ws° = local_ws.sub_ws_diff, local_ws.sub_ws_diff°

    num_recordings = length(ws.P)

    # this loop should be entirely parallelizable, as recordings
    # are—by design—conditionally independent
    for i in 1:num_recordings
        _, ws°.ll[i] = forward_guide!(
            ws.P[i], ws°.XX[i], ws°.WW[i], ws.WW[i], update.ρs[i], _LL,
            ws.XX[i][1].x[1]; Wnr=ws°.Wnr[i]
        )
        eMCMC.accept_reject!(
            update, global_ws, local_ws, step, ws.ll[i], ws°.ll[i], i
        )
    end
end

eMCMC.log_transition_density(::Any, updt::PathImputation, args...) = 0.0

eMCMC.log_prior(::Any, updt::PathImputation, args...) = 0.0

#̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶
#               Decide what happens at the Accept/Reject step
#̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶

function eMCMC.register_accept_reject_results!(
        accepted::Bool,
        updt::PathImputation,
        global_ws::DiffusionGlobalWorkspace,
        local_ws::DiffusionLocalWorkspace,
        step,
        i=1
    )
    accepted && swap_paths!(local_ws, i)
    save_path!(global_ws, step, i)
    eMCMC.register_accept_reject_results!(accepted, updt, local_ws, step, i)
end

_TEMP_SAVED_PATHS = []

function save_path!(global_ws, step, i)
    save_now = ( step.mcmciter % 100 == 0 && step.pidx == 1 ) #TODO change
    save_now && save!(global_ws.XX_buffer, global_ws.sub_ws_diff.XX)
end

#TODO there is nothing special about this step, incorporate in generic eMCMC
function eMCMC.register_accept_reject_results!(
        _accepted::Bool,
        updt,
        ws::DiffusionLocalWorkspace,
        step,
        i=1,
    )
    register!(updt, _accepted) #TODO
    accepted(ws, step.mcmciter)[i] = _accepted

    # register proposal
    ll°(ws, step.mcmciter)[i] = ll°(ws)[i]

    # update local ll with accepted
    _accepted && ( ll(ws)[i] = ll°(ws)[i] )

    # register accepted
    ll(ws, step.mcmciter)[i] = ll(ws)[i]
end

function register!(updt, accepted::Bool)
    nothing #TODO implement
end

function swap!(A, B)
    for i in eachindex(A, B)
        A[i], B[i] = B[i], A[i]
    end
end

function swap_paths!(ws::DiffusionLocalWorkspace, i)
    s°, s = ws.sub_ws_diff°, ws.sub_ws_diff
    swap!(s°.WW[i], s.WW[i])
    swap!(s°.XX[i], s.XX[i])
end

#=
function recompute_loglikhd!(ws::DiffusionLocalWorkspace, i) #TODO deprecate
    num_segments = length(ws.XX[i])
    ll_tot = 0.0
    for j in 1:num_segments
        ll_tot += loglikhd(ws.XX[i][j], ws.P[i][j])
    end
    ll(ws)[i] = ll_tot
end
=#
#===============================================================================
                        Starting point imputation
===============================================================================#

function eMCMC.update!(
        update::StartingPointsImputation,
        global_ws::DiffusionGlobalWorkspace,
        local_ws::DiffusionLocalWorkspace,
        step,
    )
    ws, ws° = ws.sub_ws_diff, ws.sub_ws_diff°
    num_recordings = length(ws.P)

    # this loop should be entirely parallelizable, as recordings
    # are---by design---conditionally independent
    for i in 1:num_recordings
        z°, y° = proposal_start_pt(ws, i) #TODO sample new sp
        _, ws°.ll[i] = solve_and_ll!(y°, ws.WW[i], ws°.XX[i], ws.P[i])
        accept_reject!(
            update, global_ws, local_ws, step, ws.ll[i], ws°.ll[i], i
        )
    end
end

log_transition_density(::Any, updt::StartingPointsImputation, args...) = 0.0

function log_prior(
        ::eMCMC.Proposal,
        updt::StartingPointsImputation,
        ws::DiffusionLocalWorkspace,
        i,

    )
    y° = ws.sub_ws_diff°.XX[i][1].y[1]
    # a little bit of cheating with the second term (not a prior)
    logpdf(ws.x0_prior, y°)# + log_likelihood_obs(ws.sub_ws_gp.P[1], y°) #TODO integrate the latter directly into ll
end

function log_prior(
        ::eMCMC.Previous,
        updt::StartingPointsImputation,
        ws::DiffusionLocalWorkspace,
        i,
    )
    y = ws.sub_ws_diff.XX[i][1].y[1]
    # a little bit of cheating with the second term (not a prior)
    logpdf(ws.x0_prior, y)# + log_likelihood_obs(ws.sub_ws_gp.P[1], y) #TODO integrate the latter directly into ll
end

#̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶
#               Decide what happens at the Accept/Reject step
#̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶

function eMCMC.register_accept_reject_results!(
        accepted::Bool,
        updt::StartingPointsImputation,
        global_ws::DiffusionGlobalWorkspace,
        local_ws::DiffusionLocalWorkspace,
        step,
        i=1
    )
    accepted && swap_starting_pts_and_paths!(ws, i)
    eMCMC.register_accept_reject_results!(accepted, updt, local_ws, step, i)
end

function swap_starting_pts_and_paths!(ws::DiffusionLocalWorkspace, i)
    s°, s = ws.sub_ws_diff°, ws.sub_ws_diff
    swap!(s°.z0s[i], s.z0s[i])
    swap!(s°.XX[i], s.XX[i])
end

#===============================================================================
                            Parameter update
===============================================================================#
function eMCMC.set_proposal!(
        updt::eMCMC.MCMCParamUpdate,
        global_ws::DiffusionGlobalWorkspace,
        local_ws::DiffusionLocalWorkspace,
        step,
    )
    ws, ws° = local_ws.sub_ws_diff, local_ws.sub_ws_diff°
    num_recordings = length(ws.WW)
    eMCMC._set_local_state°(updt, global_ws, local_ws, step)

    for i in 1:num_recordings
        critical_change = DD.set_parameters!(
            ws.P[i],
            ws°.P[i],
            state°(local_ws),
            local_ws.loc_pnames.vp_names[i],
            local_ws.loc_pnames.vp_names_aux[i],
            local_ws.loc_pnames.θ_local_names[i],
            local_ws.loc_pnames.θ_local_aux_names[i],
            local_ws.loc_pnames.θ_local_obs[i],
            local_ws.critical_param_change[i],
        )
        step.same_layout[i] || (critical_change = true)

        #-----------------------------------------------#
        # this is where most of the computation happens #
        critical_change && backward_filter!(ws°.P[i])
        #-----------------------------------------------#

        success, ws°.ll[i] = GP.solve_and_ll!(
            ws°.XX[i], ws.WW[i], ws°.P[i], ws.XX[i][1].x[1]
        )
        success || return
    end
end

#=

function eMCMC.set_parameters!(
        ::eMCMC.Proposal,
        updt::RandomWalkUpdate,
        global_ws::DiffusionGlobalWorkspace,
        ws::DiffusionLocalWorkspace,
    )
    OBS.set_parameters!(
        global_ws.sub_ws_diff°.P,
        global_ws.data,
        eMCMC.coords(updt),
        eMCMC.invcoords(updt),
        state°(ws)
    )
end

function eMCMC.set_parameters!(
        ::eMCMC.Previous,
        updt::RandomWalkUpdate,
        global_ws::DiffusionGlobalWorkspace,
        ws::DiffusionLocalWorkspace,
    )
    OBS.set_parameters!(
        global_ws.sub_ws_diff.P,
        global_ws.data,
        eMCMC.coords(updt),
        eMCMC.invcoords(updt),
        eMCMC.subidx(state(global_ws), updt)
    )
end
=#

#=
function eMCMC.update!()
    num_recordings = length(ws.P)
    for i in 1:num_recordings
        num_intv = length(ws.P[i])
        success, ll° = forward_guide!(ws, update, i)
        ll = path_loglikhd(ws, ws_cached, i) # use cashed values or recompute
        accept_reject!(updt, ws, ws_global, i, ll, ll°)
    end
end
=#

function eMCMC.compute_ll!(
        updt::eMCMC.MCMCParamUpdate,
        global_ws::DiffusionGlobalWorkspace,
        local_ws::DiffusionLocalWorkspace,
        step,
    )
    ll°(local_ws, step.mcmciter) .= ll°(local_ws)
end

function log_prior(
        ::eMCMC.Proposal,
        updt::eMCMC.MCMCParamUpdate,
        global_ws::DiffusionGlobalWorkspace,
        ws::DiffusionLocalWorkspace,
        i,
    )
    y° = ws.sub_ws_gp°.XX[i][1].y[1]
    # possibly get rid of the first one
    logpdf(ws.x0_prior, y°) + logpdf(updt.prior, state°(local_ws))
end

function log_prior(
        ::eMCMC.Previous,
        updt::eMCMC.MCMCParamUpdate,
        global_ws::DiffusionGlobalWorkspace,
        ws::DiffusionLocalWorkspace,
        i,
    )
    y = ws.sub_ws_gp.XX[i][1].y[1]
    # possibly get rid of the first one
    logpdf(ws.x0_prior, y) + logpdf(updt.prior, state(local_ws))
end

#̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶
#               Decide what happens at the Accept/Reject step
#̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶

function eMCMC.register_accept_reject_results!(
        accepted::Bool,
        updt::eMCMC.MCMCParamUpdate,
        global_ws::DiffusionGlobalWorkspace,
        local_ws::DiffusionLocalWorkspace,
        step,
        ::Any
    )
    accepted && swap_laws_and_paths!(local_ws)
    eMCMC.register_accept_reject_results!(accepted, updt, local_ws, step)
    eMCMC.set_chain_param!(accepted, updt, global_ws, local_ws, step)
end

function eMCMC.register_accept_reject_results!(
        _accepted::Bool,
        updt::eMCMC.MCMCParamUpdate,
        ws::DiffusionLocalWorkspace,
        step,
    )
    register!(updt, _accepted) #TODO
    set_accepted!(ws, step.mcmciter, _accepted)

    # register proposal
    ll°(ws, step.mcmciter)[1] = ll°(ws)[1]

    # update local ll with accepted
    _accepted && ( ll(ws)[1] = ll°(ws)[1] )

    # register accepted
    ll(ws, step.mcmciter)[1] = ll(ws)[1]
end

function swap_laws_and_paths!(ws::DiffusionLocalWorkspace)
    s°, s = ws.sub_ws_diff°, ws.sub_ws_diff
    for i in eachindex(s°.P, s.P)
        swap!(s°.P[i], s.P[i])
        swap!(s°.XX[i], s.XX[i])
    end
end

#===============================================================================
                Starting point imputation via Langevin updates
===============================================================================#

function eMCMC.update!(
        update::StartingPointsLangevinImputation,
        ws::DiffusionLocalWorkspace,
        ws_global::DiffusionGlobalWorkspace,
        step
    )
    nothing
end

#===============================================================================
                            Temporary measures
===============================================================================#
# temporary measure, make sure online stats work for Diffusions before removing
function eMCMC.update!(
        updt::eMCMC.MCMCParamUpdate,
        global_ws::DiffusionGlobalWorkspace,
        ws::DiffusionLocalWorkspace,
        step,
    )
    #step.mcmciter % 1 == 0 && println("$(step.mcmciter)")
    eMCMC.proposal!(updt, global_ws, ws, step) # specific to updt, see `updates.jl`
    eMCMC.set_proposal!(updt, global_ws, ws, step)
    eMCMC.compute_ll!(updt, global_ws, ws, step)
    eMCMC.accept_reject!(updt, global_ws, ws, step)
    #eMCMC.update_stats!(global_ws, ws, step)
end


text:
function OLSMultipleLinearRegression()
    return OLSMultipleLinearRegression(())
end

function OLSMultipleLinearRegression(arg0::jdouble)
    return OLSMultipleLinearRegression((jdouble,), arg0)
end

function calculate_adjusted_r_squared(obj::OLSMultipleLinearRegression)
    return jcall(obj, "calculateAdjustedRSquared", jdouble, ())
end

function calculate_hat(obj::OLSMultipleLinearRegression)
    return jcall(obj, "calculateHat", RealMatrix, ())
end

function calculate_r_squared(obj::OLSMultipleLinearRegression)
    return jcall(obj, "calculateRSquared", jdouble, ())
end

function calculate_residual_sum_of_squares(obj::OLSMultipleLinearRegression)
    return jcall(obj, "calculateResidualSumOfSquares", jdouble, ())
end

function calculate_total_sum_of_squares(obj::OLSMultipleLinearRegression)
    return jcall(obj, "calculateTotalSumOfSquares", jdouble, ())
end

function new_sample_data(obj::OLSMultipleLinearRegression, arg0::Vector{jdouble}, arg1::Vector{Vector{jdouble}})
    return jcall(obj, "newSampleData", void, (Vector{jdouble}, Vector{Vector{jdouble}}), arg0, arg1)
end

function new_sample_data(obj::OLSMultipleLinearRegression, arg0::Vector{jdouble}, arg1::jint, arg2::jint)
    return jcall(obj, "newSampleData", void, (Vector{jdouble}, jint, jint), arg0, arg1, arg2)
end


text:
using TiledIteration, OffsetArrays, LoopVectorization
using Test
using Documenter
using OffsetArrays: IdentityUnitRange

doctestfilters = [
    r"{([a-zA-Z0-9]+,\s?)+[a-zA-Z0-9]+}",
    r"(Array{[a-zA-Z0-9{}\s,]+,\s?1}|Vector{[a-zA-Z0-9{}\s,]+})",
    r"(Array{[a-zA-Z0-9{}\s,]+,\s?2}|Matrix{[a-zA-Z0-9{}\s,]+})",
]
Documenter.doctest(TiledIteration, doctestfilters = doctestfilters)

@testset "TileIterator small examples" begin
    titr = @inferred TileIterator((1:10,), RelaxLastTile((3,)))
    @test titr == [(1:3,), (4:6,), (7:9,), (10:10,)]

    titr = @inferred TileIterator((1:10,), RelaxStride((3,)))
    @test titr == [(1:3,), (3:5,), (6:8,), (8:10,)]

    titr = @inferred TileIterator((1:4,), RelaxStride((2,)))
    @test titr == [(1:2,), (3:4,)]

    titr = @inferred TileIterator((1:4,), RelaxLastTile((2,)))
    @test titr == [(1:2,), (3:4,)]

    titr = @inferred TileIterator((1:3, 0:5), RelaxLastTile((2, 3)))
    @test titr == [(1:2, 0:2) (1:2, 3:5); (3:3, 0:2)  (3:3, 3:5)]

    titr = @inferred TileIterator((1:3, 0:5), (2, 3))
    @test titr == [(1:2, 0:2)  (1:2, 3:5); (3:3, 0:2)  (3:3, 3:5)]

    titr = @inferred TileIterator((1:3, 0:5), RelaxStride((2, 3)))
    @test titr == [(1:2, 0:2)  (1:2, 3:5); (2:3, 0:2)  (2:3, 3:5)]

    @testset "Exotic ranges" begin
        A = zeros(10)
        AO = OffsetArray(A, OffsetArrays.Origin(0))
        titr = @inferred TileIterator(axes(AO), RelaxLastTile((3, )))
        @test titr == [(0:2,), (3:5,), (6:8,), (9:9,)]
        titr = @inferred TileIterator(axes(AO), RelaxStride((5, )))
        @test titr == [(0:4,), (5:9,)]

        titr = @inferred TileIterator((IdentityUnitRange(-4:0),), RelaxLastTile((2,)))
        @test titr == [(-4:-3,), (-2:-1,), (0:0,)]

        titr = TileIterator((Base.OneTo(4), IdentityUnitRange(0:3)), RelaxStride((3,2)))
        @test titr == [(1:3, 0:1) (1:3, 2:3); (2:4, 0:1) (2:4, 2:3)]

        titr = TileIterator((Base.OneTo(4), IdentityUnitRange(0:3), 2:2), RelaxStride((3,2,1)))
        @test axes(titr) == (1:2,1:2,1:1)
        @test titr[1,1,1] === (1:3, 0:1, 2:2)
        @test titr[1,2,1] === (1:3, 2:3, 2:2)
        @test titr[2,1,1] === (2:4, 0:1, 2:2)
        @test titr[2,2,1] === (2:4, 2:3, 2:2)
        @test titr[1,1,1:1] == [titr[1,1,1]]
        @test titr[:,1:2,1:1] == titr
    end

end

@testset "tiled iteration" begin
    sz = (3,5)
    for i1 = -3:2, i2 = -1:5
        for l1 = 2:13, l2 = 1:16
            inds = (i1:i1+l1-1, i2:i2+l2-1)
            A = fill!(OffsetArray{Int}(undef, inds), 0)
            k = 0
            for tileinds in TileIterator(inds, sz)
                tile = A[tileinds...]
                @test !isempty(tile)
                @test all(tile .== 0)
                A[tileinds...] .= (k+=1)
            end
            @test minimum(A) == 1
            @test eltype(TileIterator(inds, sz)) == Tuple{UnitRange{Int}, UnitRange{Int}}
        end
    end
end

@testset "edge iteration" begin
    iter = EdgeIterator(CartesianIndices((-1:4,0:3)), CartesianIndices((1:3,1:2)))
    @test collect(iter) == [CartesianIndex((-1,0)),
                            CartesianIndex(( 0,0)),
                            CartesianIndex(( 1,0)),
                            CartesianIndex(( 2,0)),
                            CartesianIndex(( 3,0)),
                            CartesianIndex(( 4,0)),
                            CartesianIndex((-1,1)),
                            CartesianIndex(( 0,1)),
                            CartesianIndex(( 4,1)),
                            CartesianIndex((-1,2)),
                            CartesianIndex(( 0,2)),
                            CartesianIndex(( 4,2)),
                            CartesianIndex((-1,3)),
                            CartesianIndex(( 0,3)),
                            CartesianIndex(( 1,3)),
                            CartesianIndex(( 2,3)),
                            CartesianIndex(( 3,3)),
                            CartesianIndex(( 4,3))]
    iter = EdgeIterator((0:3,0:3), (1:3,1:2))
    @test collect(iter) == [CartesianIndex(( 0,0)),
                            CartesianIndex(( 1,0)),
                            CartesianIndex(( 2,0)),
                            CartesianIndex(( 3,0)),
                            CartesianIndex(( 0,1)),
                            CartesianIndex(( 0,2)),
                            CartesianIndex(( 0,3)),
                            CartesianIndex(( 1,3)),
                            CartesianIndex(( 2,3)),
                            CartesianIndex(( 3,3))]
    iter = EdgeIterator((1:4,0:3), (1:3,1:2))
    @test collect(iter) == [CartesianIndex(( 1,0)),
                            CartesianIndex(( 2,0)),
                            CartesianIndex(( 3,0)),
                            CartesianIndex(( 4,0)),
                            CartesianIndex(( 4,1)),
                            CartesianIndex(( 4,2)),
                            CartesianIndex(( 1,3)),
                            CartesianIndex(( 2,3)),
                            CartesianIndex(( 3,3)),
                            CartesianIndex(( 4,3))]
    iter = EdgeIterator((-1:4,1:3), (1:3,1:2))
    @test collect(iter) == [CartesianIndex((-1,1)),
                            CartesianIndex(( 0,1)),
                            CartesianIndex(( 4,1)),
                            CartesianIndex((-1,2)),
                            CartesianIndex(( 0,2)),
                            CartesianIndex(( 4,2)),
                            CartesianIndex((-1,3)),
                            CartesianIndex(( 0,3)),
                            CartesianIndex(( 1,3)),
                            CartesianIndex(( 2,3)),
                            CartesianIndex(( 3,3)),
                            CartesianIndex(( 4,3))]
    iter = EdgeIterator((-1:4,0:2), (1:3,1:2))
    @test collect(iter) == [CartesianIndex((-1,0)),
                            CartesianIndex(( 0,0)),
                            CartesianIndex(( 1,0)),
                            CartesianIndex(( 2,0)),
                            CartesianIndex(( 3,0)),
                            CartesianIndex(( 4,0)),
                            CartesianIndex((-1,1)),
                            CartesianIndex(( 0,1)),
                            CartesianIndex(( 4,1)),
                            CartesianIndex((-1,2)),
                            CartesianIndex(( 0,2)),
                            CartesianIndex(( 4,2))]
    iter = EdgeIterator((-1:4,1:2), (1:3,1:2))
    @test collect(iter) == [CartesianIndex((-1,1)),
                            CartesianIndex(( 0,1)),
                            CartesianIndex(( 4,1)),
                            CartesianIndex((-1,2)),
                            CartesianIndex(( 0,2)),
                            CartesianIndex(( 4,2))]
    iter = EdgeIterator((1:3,0:3), (1:3,1:2))
    @test collect(iter) == [CartesianIndex(( 1,0)),
                            CartesianIndex(( 2,0)),
                            CartesianIndex(( 3,0)),
                            CartesianIndex(( 1,3)),
                            CartesianIndex(( 2,3)),
                            CartesianIndex(( 3,3))]
    iter = EdgeIterator((1:3,1:2), (1:3,1:2))
    @test collect(iter) == []
    @test_throws DimensionMismatch EdgeIterator((1:3,1:1), (1:3,1:2))
    @test_throws DimensionMismatch EdgeIterator((1:3,1:2), (1:4,1:2))
    iter = EdgeIterator(CartesianIndices((0:4,)), CartesianIndices(1:3))
    @test collect(iter) == [CartesianIndex(0,),
                            CartesianIndex(4,)]
    iter = EdgeIterator((Base.OneTo(2),Base.OneTo(3)),(2:2,1:3))
    for it in iter
        @test it ∈ iter.outer
    end
    iter = EdgeIterator(CartesianIndices((1:4,)),CartesianIndices((2:4,)))
    for it in iter
        @test it ∈ iter.outer
    end
end

@testset "padded sizes" begin
    @test @inferred(padded_tilesize(UInt8, (1,))) == (2^14,)
    @test @inferred(padded_tilesize(UInt16, (1,))) == (2^13,)
    @test @inferred(padded_tilesize(Float64, (1,))) == (2^11,)
    @test @inferred(padded_tilesize(Float64, (1,1))) == (2^11, 1)
    @test @inferred(padded_tilesize(Float64, (2,1))) == (2^11, 1)
    shp = @inferred(padded_tilesize(Float64, (2,2)))
    @test all(x->x>2, shp)
    @test 2^13 <= prod(shp)*8 <= 2^14
    shp = @inferred(padded_tilesize(Float64, (3,3,3)))
    @test all(x->x>3, shp)
    @test 2^13 <= prod(shp)*8 <= 2^14
end

@testset "threads" begin
    function sumtiles(A, sz)
        indsA = axes(A)
        iter = TileIterator(indsA, sz)
        sums = zeros(eltype(A), length(iter))
        for (i,tileinds) in enumerate(iter)
            v = view(A, tileinds...)
            if sums[i] == 0
                sums[i] = sum(v)
            else
                sums[i] = -sum(v)
            end
        end
        sums
    end
    function sumtiles_threaded(A, sz)
        indsA = axes(A)
        iter = TileIterator(indsA, sz)
        allinds = collect(iter)
        sums = zeros(eltype(A), length(iter))
        Threads.@threads for i = 1:length(allinds)
            tileinds = allinds[i]
            v = view(A, tileinds...)
            if sums[i] == 0
                sums[i] = sum(v)
            else
                sums[i] = -sum(v)
            end
        end
        sums
    end

    Asz, tilesz = Base.JLOptions().can_inline == 1 ? ((1000,1000), (100,100)) : ((100,100), (10,10))
    A = rand(Float64, Asz)
    snt = sumtiles(A, tilesz)
    st = sumtiles_threaded(A, tilesz)
    @test snt == st
    @test all(st .> 0)
end

setoob!(v, val) = @inbounds v[24] = val
getoob(v) = @inbounds v[24]

@testset "TileBuffer" begin
    a = fill!(Array{Int}(undef,5,5), 0)
    v = @inferred(TileBuffer(a, (2:3, 97:99)))
    @test ndims(v) == 2
    @test eltype(v) == Int
    @test axes(v) == (2:3, 97:99)
    v[2,97] = 1
    @test a[1] == 1
    @test strides(v) == (1, 2)
    @test @inferred(TiledIteration.viewtype(typeof(v))) === typeof(v.view)
    @test @inferred(TiledIteration.buftype(typeof(v))) === typeof(v.buf)
    p = pointer(v)
    v = @inferred(TileBuffer(v, (-1:1, 1:1)))
    @test pointer(v) == p
    @test axes(v) == (-1:1, 1:1)
    @test v[-1,1] == 1
    @test v[0,1] == v[1,1] == 0
    v[0,1] = 2
    @test a[2] == 2
    v = @inferred(TileBuffer(v, (0:23,)))
    @test pointer(v) == p
    @test eltype(v) == Int
    @test ndims(v) == 1
    @test axes(v) == (0:23,)
    @test v[1] == 2
    @test_throws BoundsError v[24]
    opts = Base.JLOptions()
    if opts.can_inline == 1 && opts.check_bounds == 0
        setoob!(v, 7)
        @test a[25] == 7
        @test getoob(v) == 7
    end
    @test_throws DimensionMismatch TileBuffer(a, (1:6,1:5))
    @test_throws DimensionMismatch TileBuffer(v, (1:5,1:6))
    @test_throws DimensionMismatch TileBuffer(a, (0:25,))
    b = TileBuffer(Float64, (1:16,1:4))
    @test eltype(b) == Float64
    @test axes(b) == (1:16,1:4)
    @test pointer(b) != p
    p = pointer(b)
    b = TileBuffer(b, (17:32,1:4))
    @test pointer(b) == p
    @inbounds begin
        b[17,2] = 5.2
        @test b[17,2] == 5.2
    end

    @static if Base.VERSION >= v"1.5"
        a = reshape(Float32.(1:25), (5, 5))
        tb = TileBuffer(a, (2:4, 2:5))
        function sum_with_lv(A)
            s = zero(eltype(A))
            @turbo for i in eachindex(A)
                s += A[i]
            end
            return s
        end
        @test sum_with_lv(a[1:12]) == sum_with_lv(tb) == sum(tb) == 78
    end
end

@testset "SplitAxis" begin
    function meets_criteria(ax::AbstractUnitRange, splits::AbstractVector)
        @test issorted(splits)
        @test first(splits) + 1 == first(ax)  # the +1 comes from splits[i]+1:splits[i+1]
        @test last(splits) == last(ax)
        @test firstindex(splits) == 1
        d1 = splits[2] - splits[1]
        for i = 2:length(splits)-1
            @test splits[i+1] - splits[i] >= d1
        end
        return
    end
    function meets_criteria(sax, ax, n::Real)
        meets_criteria(ax, sax.splits)
        @test sax[1] isa UnitRange{Int} && first(sax[1]) == first(ax)
        @test reduce(vcat, [sax[i] for i = 1:ceil(Int, n)]) == collect(ax) # collect is in case `ax` is indexed by something other than 1:k
        @test collect(sax) == [sax[i] for i = 1:ceil(Int, n)]
    end
    function meets_criteria(ax::AbstractUnitRange, n::Real)
        meets_criteria(SplitAxis(ax, n), ax, n)
    end

    meets_criteria(1:8, 3)
    meets_criteria(Base.OneTo(8), 3)
    meets_criteria(2:8, 3)
    meets_criteria(3:8, 3)
    meets_criteria(1:9, 3)
    meets_criteria(1:8, 2.5)
    meets_criteria(Base.OneTo(8), 2.5)
    meets_criteria(2:8, 2.5)
    meets_criteria(3:8, 2.5)
    meets_criteria(1:9, 2.5)

    function meets_criteria(axs::Base.Indices, n::Real)
        saxs = SplitAxes(axs, n)
        meets_criteria(saxs.splitax, axs[end], n)
        function check_and_getlast(chunk)
            @test chunk[1:end-1] == axs[1:end-1]
            return chunk[end]
        end
        @test reduce(vcat, [check_and_getlast(saxs[i]) for i = 1:ceil(Int, n)]) == collect(axs[end])
        @test collect(saxs) == [saxs[i] for i = 1:ceil(Int, n)]
    end

    meets_criteria((1:5, 1:8), 3)
    meets_criteria((Base.OneTo(5), Base.OneTo(8)), 3)
    meets_criteria((1:5, 2:8), 3)
    meets_criteria((1:5, 3:8), 3)
    meets_criteria((1:5, 1:9), 3)
    meets_criteria((1:5, 1:8), 2.5)
    meets_criteria((Base.OneTo(5), Base.OneTo(8)), 2.5)
    meets_criteria((1:5, 2:8), 2.5)
    meets_criteria((1:5, 3:8), 2.5)
    meets_criteria((1:5, 1:9), 2.5)
end

nothing


text:
using NS2D
using Test

@testset "NS2D.jl" begin
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "ViennaRNA"
version = v"2.4.18"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/ViennaRNA/ViennaRNA/releases/download/v$(version)/ViennaRNA-$(version).tar.gz",
                  "b276cfed7c3bea4821c8272750b5b22fa6fed614d71a98cfb25eb0095369657c")
]

# Bash recipe for building across all platforms
# Notes:
# - openmp not supported on darwin
# - libsvm is included in the ViennaRNA sources, and it's not
#   immediately obvious how to make the build script use libsvm_jll
# --with-forester fails to compile on windows
#     In file included from aligner.cpp:3:
#     aligner.h:21:10: fatal error: sys/times.h: No such file or directory
script = raw"""
cd $WORKSPACE/srcdir
cd ViennaRNA-*/

# configure script needs bash
sed -i -e '1s|#! /bin/sh|#! /bin/bash|' configure

# set this explicitly for the libsvm Makefile, which would
# otherwise use g++
export CC=cc
export CXX=c++

# help freebsd find include files for mpfr, gsl
export CPPFLAGS="-I${includedir}"

# configure script fails (only in sandbox) without
# setting ac_cv_func_malloc_0_nonnull=yes ac_cv_func_realloc_0_nonnull=yes
ac_cv_func_malloc_0_nonnull=yes ac_cv_func_realloc_0_nonnull=yes \
    ./configure \
    --prefix=${prefix} --build=${MACHTYPE} --host=${target} \
    --with-pic --disable-c11 \
    --with-mpfr --with-json --with-svm --with-gsl \
    --without-forester --with-cluster --with-kinwalker \
    --without-perl --without-python --without-python3 \
    --without-doc --without-tutorial --without-tutorial-pdf

make -j${nproc}
make install

# create and install a shared library libRNA
ldflags="-g -O2 -fno-strict-aliasing -ftree-vectorize -pthread"
# add for gcc (linux,windows): -flto -ffat-lto-objects
# TODO: how to check if compiler is gcc/g++ ?
if [[ $target != *-darwin* ]]; then
    ldflags="$ldflags -fopenmp"
fi
ldflags_libs="-lpthread -lmpfr -lgmp -lstdc++ -lgsl -lgslcblas -lm"
$CC -shared -o "${libdir}/libRNA.${dlext}" \
    $ldflags \
    -Wl,$(flagon --whole-archive) ./src/ViennaRNA/libRNA.a -Wl,$(flagon --no-whole-archive) \
    $ldflags_libs

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)
# remedy for std::string incompatibilities across the GCC 4/5 version boundary
# (suggested by BinaryBuilder)
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("AnalyseDists", :AnalyseDists),
    ExecutableProduct("AnalyseSeqs", :AnalyseSeqs),
    ExecutableProduct("b2ct", :b2ct),
    ExecutableProduct("ct2db", :ct2db),
    ExecutableProduct("Kinfold", :Kinfold),
    ExecutableProduct("kinwalker", :Kinfold),
    ExecutableProduct("popt", :popt),
    ExecutableProduct("RNA2Dfold", :RNA2Dfold),
    ExecutableProduct("RNAaliduplex", :RNAaliduplex),
    ExecutableProduct("RNAalifold", :RNAalifold),
    ExecutableProduct("RNAcofold", :RNAcofold),
    ExecutableProduct("RNAdistance", :RNAdistance),
    ExecutableProduct("RNAdos", :RNAdos),
    ExecutableProduct("RNAduplex", :RNAduplex),
    ExecutableProduct("RNAeval", :RNAeval),
    ExecutableProduct("RNAfold", :RNAfold),
    # enable with --with-forester
    # ExecutableProduct("RNAforester", :RNAforester),
    ExecutableProduct("RNAheat", :RNAheat),
    ExecutableProduct("RNAinverse", :RNAinverse),
    ExecutableProduct("RNALalifold", :RNALalifold),
    ExecutableProduct("RNALfold", :RNALfold),
    ExecutableProduct("RNAlocmin", :RNAlocmin),
    ExecutableProduct("RNApaln", :RNApaln),
    ExecutableProduct("RNAparconv", :RNAparconv),
    ExecutableProduct("RNApdist", :RNApdist),
    ExecutableProduct("RNAPKplex", :RNAPKplex),
    ExecutableProduct("RNAplex", :RNAplex),
    ExecutableProduct("RNAplfold", :RNAplfold),
    ExecutableProduct("RNAplot", :RNAplot),
    ExecutableProduct("RNApvmin", :RNApvmin),
    ExecutableProduct("RNAsnoop", :RNAsnoop),
    ExecutableProduct("RNAsubopt", :RNAsubopt),
    ExecutableProduct("RNAup", :RNAup),
    LibraryProduct("libRNA", :libRNA),
]

# Dependencies that must be installed before this package can be built
# CompilerSupportLibraries_jll suggested by BinaryBuilder for OpenMP support
dependencies = [
    Dependency(PackageSpec(name="MPFR_jll", uuid="3a97d323-0669-5f0c-9066-3539efd106a3"))
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
    Dependency(PackageSpec(name="GSL_jll", uuid="1b77fbbe-d8ee-58f0-85f9-836ddc23a7a4"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               julia_compat="1.6", preferred_gcc_version = v"11")

max_stars_repo_name: ederc/Yggdrasil

max_issues_repo_name: ederc/Yggdrasil

max_forks_repo_name: ederc/Yggdrasil


text:
using DEBTrait, Test

pct_sand = 70.0
pct_clay = 30.0
s_sat    = 0.1
N_cell   = 10.0*ones(1)
Soil = BaseSoil(pct_sand, pct_clay, s_sat, N_cell)

chb, sat, psisat, ksat = DEBTrait.clapp_hornberger_par(Soil.pct_sand, Soil.pct_clay)
@test chb ≈ 8.76 atol=1e-6
@test sat ≈ 0.4008 atol=1e-6
@test psisat ≈ -9.1833 atol=1e-4
@test ksat ≈ 0.01085 atol=1e-4

psi, dpsidvsm = DEBTrait.cosby_psi(Soil.s_sat, psisat, sat, chb)
@test psi  ≈ -1e9 atol=1e-6
@test dpsidvsm ≈ 218562874251.49704 atol=1e-6

theta = s_sat*sat # phi_w
epsi = sat-theta
taug, tauw = DEBTrait.moldrup_tau(sat, epsi, theta, chb)
@test taug  ≈ 0.3479 atol=1e-4
@test tauw ≈ 0.000481 atol=1e-6

Dwpsi = DEBTrait.cosby_Dwpsi(s_sat, psisat, ksat, sat, chb)
@test Dwpsi ≈ 7.163253899338254e-19 atol=1e-6

s_sat = 0.10101010101010102
theta, τ_g, τ_w, film = DEBTrait.soil_affinity_properties(Soil.pct_sand, Soil.pct_clay, s_sat)
@test theta ≈ 0.04048 atol=1e-5
@test τ_g ≈ 0.347412 atol=1e-6
@test τ_w ≈ 0.000496 atol=1e-6
@test film ≈ 1e-7 atol=1e-6