text:
const CDM = Ptr{Cvoid}
abstract type AbstractDM{PetscLib} end

"""
    DMLocalVec(v::CVec, dm::AbstractDM)

Container for an PETSc vector we know is "local"

# External Links
$(_doc_external("Vec/Vec"))
"""
mutable struct DMLocalVec{PetscLib, T, T_DM} <: AbstractVec{T}
    ptr::CVec
    dm::T_DM
    function DMLocalVec(ptr, dm::AbstractDM{PetscLib}) where {PetscLib}
        new{PetscLib, scalartype(PetscLib), typeof(dm)}(ptr, dm)
    end
end

"""
    DMGlobalVec(v::CVec, dm::AbstractDM)

Container for an PETSc vector we know is "global"

# External Links
$(_doc_external("Vec/Vec"))
"""
mutable struct DMGlobalVec{PetscLib, T, T_DM} <: AbstractVec{T}
    ptr::CVec
    dm::T_DM
    function DMGlobalVec(ptr, dm::AbstractDM{PetscLib}) where {PetscLib}
        new{PetscLib, scalartype(PetscLib), typeof(dm)}(ptr, dm)
    end
end

# Mainly for DM we do not know the type of, namely ones returned by PETSc
# functions such as `KSPGetDM`
mutable struct PetscDM{PetscLib} <: AbstractDM{PetscLib}
    ptr::CDM
end

"""
    setup!(da::DM, opts=da.opts)

# External Links
$(_doc_external("DM/DMSetUp"))
"""
function setup! end

@for_petsc function setup!(da::AbstractDM{$PetscLib}, opts::Options = da.opts)
    with(opts) do
        @chk ccall((:DMSetUp, $petsc_library), PetscErrorCode, (CDM,), da)
    end
end

"""
    setfromoptions!(da::DM, opts=da.opts)

# External Links
$(_doc_external("DM/DMSetFromOptions"))
"""
function setfromoptions! end

@for_petsc function setfromoptions!(
    da::AbstractDM{$PetscLib},
    opts::Options = da.opts,
)
    with(opts) do
        @chk ccall(
            (:DMSetFromOptions, $petsc_library),
            PetscErrorCode,
            (CDM,),
            da,
        )
    end
end

@for_petsc begin
    function destroy(da::AbstractDM{$PetscLib})
        finalized($PetscLib) || @chk ccall(
            (:DMDestroy, $petsc_library),
            PetscErrorCode,
            (Ptr{CDM},),
            da,
        )
        da.ptr = C_NULL
        return nothing
    end
end

"""
    getdimension(dm::AbstractDM)

Return the topological dimension of the `dm`

# External Links
$(_doc_external("DM/DMGetDimension"))
"""
function getdimension(::AbstractDM) end

@for_petsc function getdimension(dm::AbstractDM{$PetscLib})
    dim = Ref{$PetscInt}()

    @chk ccall(
        (:DMGetDimension, $petsc_library),
        PetscErrorCode,
        (CDM, Ptr{$PetscInt}),
        dm,
        dim,
    )

    return dim[]
end

"""
    gettype(dm::AbstractDM)

Gets type name of the `dm`

# External Links
$(_doc_external("DM/DMGetType"))
"""
function gettype(::AbstractDM) end

@for_petsc function gettype(dm::AbstractDM{$PetscLib})
    t_r = Ref{Cstring}()
    @chk ccall(
        (:DMGetType, $petsc_library),
        PetscErrorCode,
        (CDM, Ptr{Cstring}),
        dm,
        t_r,
    )
    return unsafe_string(t_r[])
end

"""
    view(dm::AbstractDM, viewer::Viewer=ViewerStdout(petsclib, getcomm(dm)))

view a `dm` with `viewer`

# External Links
$(_doc_external("DM/DMView"))
"""
function view(::AbstractDM) end

@for_petsc function view(
    dm::AbstractDM{$PetscLib},
    viewer::AbstractViewer{$PetscLib} = ViewerStdout($petsclib, getcomm(dm)),
)
    @chk ccall(
        (:DMView, $petsc_library),
        PetscErrorCode,
        (CDM, CPetscViewer),
        dm,
        viewer,
    )
    return nothing
end

"""
    creatematrix(dm::AbstractDM)

Generates a matrix from the `dm` object.

# External Links
$(_doc_external("DM/DMCreateMatrix"))
"""
function creatematrix end

@for_petsc function creatematrix(dm::AbstractDM{$PetscLib})
    mat = Mat{$PetscScalar}(C_NULL)

    @chk ccall(
        (:DMCreateMatrix, $petsc_library),
        PetscErrorCode,
        (CDM, Ptr{CMat}),
        dm,
        mat,
    )

    return mat
end

"""
    createlocalvector(dm::AbstractDM)

returns a local vector from the `dm` object.

# External Links
$(_doc_external("DM/DMCreateLocalVector"))
"""
function createlocalvector end

@for_petsc function createlocalvector(dm::AbstractDM{$PetscLib})
    vec = DMLocalVec(C_NULL, dm)

    @chk ccall(
        (:DMCreateLocalVector, $petsc_library),
        PetscErrorCode,
        (CDM, Ptr{CVec}),
        dm,
        vec,
    )

    return vec
end

"""
    createglobalvector(dm::DM; write::Bool = true, read::Bool = true)

returns a global vector from the `dm` object.

# External Links
$(_doc_external("DM/DMCreateGlobalVector"))
"""
function createglobalvector end

@for_petsc function createglobalvector(dm::AbstractDM{$PetscLib})
    vec = DMGlobalVec(C_NULL, dm)

    @chk ccall(
        (:DMCreateGlobalVector, $petsc_library),
        PetscErrorCode,
        (CDM, Ptr{CVec}),
        dm,
        vec,
    )

    return vec
end

"""
    update!(
        global_vec::DMGlobalVec,
        local_vec::DMLocalVec,
        mode::InsertMode,
    )

Updates `global_vec` from `local_vec` with insert `mode`

# External Links
$(_doc_external("DM/DMLocalToGlobal"))
"""
function update!( global_vec::DMGlobalVec, local_vec::DMLocalVec, mode::InsertMode) end

@for_petsc function update!(
    global_vec::DMGlobalVec{$PetscLib},
    local_vec::DMLocalVec{$PetscLib},
    mode::InsertMode
)
    @assert local_vec.dm === global_vec.dm

    update_local2global!( global_vec.ptr, local_vec.ptr, mode, local_vec.dm);

    return nothing
end

"""
    update!(
        global_vec::DMGlobalVec,
        local_ptr::CVec,
        mode::InsertMode,
    )

Updates `global_vec` from a pointer to a local vector `local_ptr` with insert `mode`

# External Links
$(_doc_external("DM/DMLocalToGlobal"))
"""
function update!( global_vec::DMGlobalVec, local_ptr::CVec, mode::InsertMode) end

@for_petsc function update!(
    global_vec::DMGlobalVec{$PetscLib},
    local_ptr::CVec,
    mode::InsertMode
)
    update_local2global!( global_vec.ptr, local_ptr, mode, global_vec.dm);

    return nothing
end

"""
    update!(
        global_ptr::CVec,
        local_vec::DMLocalVec,
        mode::InsertMode,
    )

Updates pointer to global vec `global_ptr` from `local_vec` with insert `mode`

# External Links
$(_doc_external("DM/DMLocalToGlobal"))
"""
function update!( global_vec::DMGlobalVec, local_ptr::CVec, mode::InsertMode) end

@for_petsc function update!(
    global_ptr::CVec,
    local_vec::DMLocalVec{$PetscLib},
    mode::InsertMode
)
    update_local2global!( global_ptr, local_vec.ptr, mode, local_vec.dm);

    return nothing
end

"""
    update!(
        local_vec::DMLocalVec,
        global_ptr::CVec,
        mode::InsertMode,
    )

Updates `local_vec` from pointer to global vec `global_ptr` with insert `mode`

# External Links
$(_doc_external("DM/DMGlobalToLocal"))
"""
function update!(local_vec::DMLocalVec,global_ptr::CVec, mode::InsertMode) end

@for_petsc function update!(
    local_vec::DMLocalVec{$PetscLib},
    global_ptr::CVec,
    mode::InsertMode,
)

    update_global2local!( local_vec.ptr, global_ptr, mode, local_vec.dm);
    return nothing
end

"""
    update!(
        local_ptr::CVec,
        global_vec::DMGlobalVec
        mode::InsertMode,
    )

Updates pointer to local vec `local_ptr` from `global_vec` with insert `mode`

# External Links
$(_doc_external("DM/DMGlobalToLocal"))
"""
function update!(local_ptr::CVec, global_vec::DMGlobalVec, mode::InsertMode) end

@for_petsc function update!(
    local_ptr::CVec,
    global_vec::DMGlobalVec{$PetscLib},
    mode::InsertMode,
)

    update_global2local!( local_ptr, global_vec.ptr, mode, global_vec.dm);
    return nothing
end

"""
    update!(
        local_vec::DMLocalVec,
        global_vec::DMGlobalVec,
        mode::InsertMode
    )

Updates `local_vec` from a pointer to a global vec `global_ptr` with insert `mode`

# External Links
$(_doc_external("DM/DMGlobalToLocal"))
"""
function update!(local_vec::DMLocalVec, global_vec::DMGlobalVec, mode::InsertMode) end

@for_petsc function update!(
    local_vec::DMLocalVec{$PetscLib},
    global_vec::DMGlobalVec{$PetscLib},
    mode::InsertMode
)
    @assert local_vec.dm === global_vec.dm

    update_global2local!( local_vec.ptr, global_vec.ptr, mode, local_vec.dm);

    return nothing
end

#=
@for_petsc function update!(
    local_ptr::CVec,
    global_vec::DMGlobalVec{$PetscLib},
    mode::InsertMode,
)
    @chk ccall(
        (:DMGlobalToLocal, $petsc_library),
        PetscErrorCode,
        (CDM, CVec, InsertMode, CVec),
        global_vec.dm,
        global_vec,
        mode,
        local_ptr,
    )

    return nothing
end
=#

"""
    update_local2global!(
        global_ptr::CVec,
        local_ptr::CVec,
        mode::InsertMode,
        dm::AbstractDM
    )

Updates pointer of `global_vec` from pointer of `local_vec` with insert `mode`.
Both vectors should belong to the same `dm`

This is a low-level routine that is typically called by `update!`
"""
function update_local2global! end

@for_petsc function update_local2global!(
    global_ptr::CVec,
    local_ptr::CVec,
    mode::InsertMode,
    dm::AbstractDM{$PetscLib}
)
    @chk ccall(
        (:DMLocalToGlobal, $petsc_library),
        PetscErrorCode,
        (CDM, CVec, InsertMode, CVec),
        dm,
        local_ptr,
        mode,
        global_ptr,
    )

    return nothing
end

"""
    update_global2local!(
        local_ptr::CVec,
        global_ptr::CVec,
        mode::InsertMode,
        dm::AbstractDM
    )

Updates pointer to `local` vector from pointer to `global` vector with insert `mode`, assuming that both belong to the same `dm`

This is a low-level routine that is typically called by `update!`

"""
function update_global2local! end

@for_petsc function update_global2local!(
    local_ptr::CVec,
    global_ptr::CVec,
    mode::InsertMode,
    dm::AbstractDM{$PetscLib}
)
    @chk ccall(
        (:DMGlobalToLocal, $petsc_library),
        PetscErrorCode,
        (CDM, CVec, InsertMode, CVec),
        dm,
        global_ptr,
        mode,
        local_ptr,
    )

    return nothing
end

"""
    getcoordinateDM(dm::AbstractDM)

Create a `coord_dm` for the coordinates of `dm`.

# External Links
$(_doc_external("DM/DMGetCoordinateDM"))
"""
function getcoordinateDM end

@for_petsc function getcoordinateDM(dm::AbstractDM{$PetscLib})
    coord_dm = empty(dm)
    @chk ccall(
        (:DMGetCoordinateDM, $petsc_library),
        PetscErrorCode,
        (CDM, Ptr{CDM}),
        dm,
        coord_dm,
    )

    # If this fails then the `empty` call above is probably a bad idea!
    if gettype(coord_dm) != "product"
        @assert gettype(dm) == gettype(coord_dm)
    else
        @assert gettype(dm) == "stag"   # product can only be used with stag
    end

    return coord_dm
end

"""
    getcoordinateslocal(dm::AbstractDM)

Gets a local vector with the coordinates associated with `dm`.

# External Links
$(_doc_external("DM/DMGetCoordinatesLocal"))
"""
function getcoordinateslocal end

@for_petsc function getcoordinateslocal(dm::AbstractDM{$PetscLib})
    coord_vec = DMLocalVec(C_NULL, dm)
    @chk ccall(
        (:DMGetCoordinatesLocal, $petsc_library),
        PetscErrorCode,
        (CDM, Ptr{CVec}),
        dm,
        coord_vec,
    )

    return coord_vec
end


text:
using Test
using DynamicalBilliards
using DynamicalBilliards.Testing

# %% Basic tests

function simple_raysplit(f, args...)
    bd, ray = basic_ray(true)
    p = randominside(bd)
    f(p, bd, ray, args...)

    bd, ray = basic_ray(false)
    mp = randominside(bd, 2.0)
    f(mp, bd, ray, args...)
end

function test_finite_collisions_ray(p, bd, ray, N = 5e5)
    n = 0
    i, t = next_collision(p, bd)
    s = Set(i)
    raysidx = DynamicalBilliards.raysplit_indices(bd, ray)

    @testset "$(tag(p)) 0 < t < Inf" begin
        while n < N
            i, t = bounce!(p, bd, raysidx, ray)
            if 0 < t < Inf
                @test true
            else
                error("$(tag(p, bd, ray)) collision time $(t) with obst.index $(i)")
            end
            n += 1
            push!(s, i)
        end
    end
    @testset "$(tag(p, bd)) collided with all obstacles" begin
        @test s == Set(1:length(bd))
    end
end

function ray_maximum_distance(p, bd, ray, thres = 1e-13, N = 5e5)
    acc = acclevel(p, bd, ray)
    n = 0
    t, i = next_collision(p, bd)
    maxds = fill(-Inf, length(bd))
    minds = fill(Inf, length(bd))
    thresholds = fill(0, length(bd))
    isp = isperiodic(bd)
    raysidx = DynamicalBilliards.raysplit_indices(bd, ray)

    @testset "$(tag(p, bd, ray)) max/min distance" begin
        while n < N

            i, t, cp = next_collision(p, bd)
            d = distance(cp, bd[i])

            d < minds[i] && (minds[i] = d)
            d > maxds[i] && (maxds[i] = d)

            i, t = bounce!(p, bd, raysidx, ray)

            @test abs(d) < acc
            abs(d) > thres && (thresholds[i] += 1)

            n += 1
        end
        println(tag(p, bd))
        println("minds: ", minds)
        println("maxds: ", maxds)
        println("more than $(thres): ", thresholds, " (acc = $(acc))")
        println()
    end
end

billiards_testset("Ray Collision time", test_finite_collisions_ray; caller = simple_raysplit)
billiards_testset("Ray min/max distance", ray_maximum_distance; caller = simple_raysplit)

# %% Extreme angles test

function extreme_raysplit(f, args...)
    bd, ray = extreme_ray(true)
    p = randominside(bd)
    f(p, bd, ray, args...)

    bd, ray = extreme_ray(false)
    mp = randominside(bd, 2.0)
    f(mp, bd, ray, args...)
end

billiards_testset("ERay Collision time", test_finite_collisions_ray; caller = extreme_raysplit)
billiards_testset("ERay min/max distance", ray_maximum_distance; caller = extreme_raysplit)

# %% Totally brutal tests

function inside_antidot(args...)
    bd, ray = extreme_ray(true)
    pa = randominside(bd)
    mp = randominside(bd, 1.0)
    isray(i) = (i == 1 || i == 2)
    raysidx = DynamicalBilliards.raysplit_indices(bd, ray)

    for p in (pa, mp); @testset "$(tag(p, bd, ray)) BRUTAL" begin
        bd, ray = extreme_ray(typeof(p) <: Particle)
        n = 0
        iprev = 0
        check = false

        while n < 5e5
            i, t = bounce!(p, bd, raysidx, ray)

            if check
                @test i == iprev
            end

            check = isray(i) && !(bd[i].pflag)
            iprev = i

            n += 1
        end
    end; end
end

billiards_testset("RAY Stays inside", identity; caller = inside_antidot)

test_refraction(ϕ, pflag, ω) = pflag ? 0.5ϕ : 2.0ϕ
test_transmission(ϕ, pflag, ω) = begin
    if pflag
        0.5*exp(-(ϕ)^2/2(π/8)^2)
    else
        abs(ϕ) < π/4 ? 0.5*exp(-(ϕ)^2/2(π/4)^2) : 0.0
    end
end

@testset "FiniteSplitterWall" begin
    bdr = billiard_rectangle(1.0, 1.0)
    splitter = FiniteSplitterWall([0.5,0.0], [0.5,1.0], [-1.0,0.0])
    bd = Billiard(splitter, bdr...)

    N = 500
    rs = (RaySplitter([1], test_transmission, test_refraction),)
    ps = particlebeam(0.01, 0.5, 0, N, 0)
    positions = []

    for p in ps
        ct, pos, vel = evolve!(p, bd, 2, rs)
        push!(positions, pos[end][1])
        reset_billiard!(bd)
    end

    particles_on_left = length(filter(x ->x < 0.4, positions))
    particles_on_right = length(filter(x -> x > 0.6, positions))
    particles_in_middle = length(filter(x -> 0.4 <= x <= 0.6, positions))

    @test 0.3N < particles_on_left < 0.7N
    @test 0.3N < particles_on_right < 0.7N
    @test particles_in_middle == 0
end

@testset "Ray-splitting escapetime" begin
    verts = [[0.0,0.0], [0.0,1.0], [1.0,1.0], [1.0,0.0]]
    bdr = []

    # create rect billiard where vertical walls are doors
    for i in eachindex(verts)
        s = verts[i]
        e = verts[mod1(i+1, length(verts))]
        w = s - e
        n = [-w[2], w[1]]
        push!(bdr, FiniteWall(s, e, n, mod1(i, 2)==1))
    end

    splitter = FiniteSplitterWall([0.75,0.0], [0.75,1.0], [-1.0,0.0])
    bd = Billiard(splitter, bdr...)

    N = 500
    rs = (RaySplitter([1], test_transmission, test_refraction),)
    ps = particlebeam(0.01, 0.5, 0, N, 0)
    ts = []

    for p in ps
        push!(ts, escapetime!(p, bd, 3, rs))
        reset_billiard!(bd)
    end

    mean_t = +(ts...)/N
    expected_mean_t = 0.5 * 1 + 0.5 * 0.75 * 2
    @test abs(expected_mean_t - mean_t) < 0.2
end

ext: jl


text:
export iauPvu
"""
Update a pv-vector.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  vector/matrix support function.

Given:
   dt       double           time interval
   pv       double[2][3]     pv-vector

Returned:
   upv      double[2][3]     p updated, v unchanged

Notes:

   1. "Update" means "refer the position component of the vector
   to a new date dt time units from the existing date".

   2. The time units of dt must match those of the velocity.

   3. It is permissible for pv and upv to be the same array.

Called:
   iauPpsp      p-vector plus scaled p-vector
   iauCp        copy p-vector

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauPvu(double dt, double pv[2][3], double upv[2][3])
function iauPvu(dt::Real, pv::AbstractMatrix{<:Real})
   # Preallocate return values
   upv = zeros(Float64, 3, 2)

   ccall((:iauPvu, libsofa_c), Cvoid,
         (Cdouble, Ptr{Cdouble}, Ref{Cdouble}),
          convert(Float64, dt),
          convert(Array{Float64, 2}, pv'),
          upv)

   return SMatrix{2,3}(upv')
end


text:
using StructuralCausalModels
using GraphRecipes, Plots

ProjDir = @__DIR__
cd(ProjDir) #do

d_string = "dag {A -> {E Z}; B -> {D Z}; Z -> {D E}; E -> D}"

dag = DAG("conf_triangles", d_string);
show(dag)

graphplot(dag.e, names=names(dag.e, 1), curvature_scalar=0.1, nodesize=0.2,
  method=:spring, fontsize=10, arrow=1.0, nodeshape=:circle)

gui()


text:
# __BEGIN_LICENSE__
#
# WaveOptics.jl
#
# Copyright (c) 2015, Stanford University
#
# All rights reserved.
#
# Redistribution and use in source and binary forms for academic and other
# non-commercial purposes with or without modification, are permitted provided
# that the following conditions are met:
#
# * Redistributions of source code, including modified source code, must retain
#   the above copyright notice, this list of conditions and the following
#   disclaimer.
#
# * Redistributions in binary form or a modified form of the source code must
#   reproduce the above copyright notice, this list of conditions and the
#   following disclaimer in the documentation and/or other materials provided with
#   the distribution.
#
# * Neither the name of The Leland Stanford Junior University, any of its
#   trademarks, the names of its employees, nor contributors to the source code
#   may be used to endorse or promote products derived from this software without
#   specific prior written permission.
#
# * Where a modified version of the source code is redistributed publicly in
#   source or binary forms, the modified source code must be published in a freely
#   accessible manner, or otherwise redistributed at no charge to anyone
#   requesting a copy of the modified source code, subject to the same terms as
#   this agreement.
#
# THIS SOFTWARE IS PROVIDED BY THE TRUSTEES OF THE LELAND STANFORD JUNIOR
# UNIVERSITY "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE LELAND STANFORD JUNIOR
# UNIVERSITY OR ITS TRUSTEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
# BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
# IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
# OF SUCH DAMAGE.
#
# __END_LICENSE__

abstract type OpticalRecipe end

# Aperture code types (instances are defined in the 'aperture' sub-module)
abstract type CodedAperture end

const MaskParameters = Dict{String,Float64}
const MaskOptions = Dict{String,Any}

# -------------------------------------------------------------------------------------------
#                                 Null Coded Aperture
# -------------------------------------------------------------------------------------------

mutable struct NullCodedAperture <: CodedAperture
    mask_name::String
    mask_parameters::MaskParameters
    mask_options::MaskOptions
end

function NullCodedAperture()
    return NullCodedAperture("null_phase", MaskParameters(), MaskOptions())
end

function NullCodedAperture(mask_parameters::MaskParameters, mask_options)
    return NullCodedAperture("null_phase", MaskParameters(), MaskOptions())
end

function mask(
    ap::NullCodedAperture,
    coordinates::OpticalCoordinates,
    parameters = Dict{String,Float64}(),
)
    return ones(Complex64, size(coordinates.rho))
end

function jacobian(
    ap::NullCodedAperture,
    coordinates::OpticalCoordinates,
    parameters = Dict{String,Float64}(),
)
    error("No Jacobian defined for NullCodedAperture")
end

# ------------------------------------------------------------------------------------------------------
#                                   Wide Field Optical Recipe
# ------------------------------------------------------------------------------------------------------

mutable struct WideFieldOpticalRecipe <: OpticalRecipe
    objective_mag::Float64
    objective_na::Float64
    wavelength::Float64
    medium_index::Float64
    f_tubelens::Float64
    f_objective::Float64
    static_aperture_mask::Any
    dynamic_aperture_mask::Any
    normalize_intensity::Bool
end

"""
Generic Optical Recipe for Widefield PSF Models

**Parameters**

    objective_mag - magnification of the objective
    objective_na - numerical aperture of the objective
    wavelength - emission wavelength of the light being imaged (all simulations are monochromatic)
    medium_index - index of refraction of the sample medium
    f_tubelens - tube lens focal length (e.g. 200mm for Nikon, 180mm for Olympus, etc.)
    static_aperture_mask - Supply a back aperture mask that is a subclass of CodedAperture, or
               'NullCodedAperture()' if there is no aperture code.  The static mask is computed
               when the PSF model is constructed, and is therefore more efficient than the
               dynamic_aperture_mask.
    dynamic_aperture_mask - Supply a back aperture mask that is a subclass of CodedAperture, or
               'NullCodedAperture()' if there is no aperture code.  The dynamic_aperture_mask
               is recomputed each time apsf() is called, but this flexibility incurs additional
               performance overhead.
    normalize_intensity - Set max PSF value @ z = 0 to 1.0.  This is an arbitrary choice, but
               one that keeps PSF values relatively well scaled and makes it easy to adjust the
               "exposure" of the simulated microscope by multipyling the PSF to get some number of
               photons that is >= 1.  However, be careful when mixing optical models with different
               recipes, because their absolute intensities cannot be directly compared with this
               arbitrary scaling in effect.
"""
function WideFieldOpticalRecipe(
    objective_mag::Number,
    objective_na::Float64,
    wavelength::Float64,
    medium_index::Float64,
    f_tubelens::Float64;
    static_aperture_mask::CodedAperture = NullCodedAperture(),
    dynamic_aperture_mask::CodedAperture = NullCodedAperture(),
    normalize_intensity::Bool = false,
)

    f_objective = f_tubelens / objective_mag
    return WideFieldOpticalRecipe(
        objective_mag,
        objective_na,
        wavelength,
        medium_index,
        f_tubelens,
        f_objective,
        static_aperture_mask,
        dynamic_aperture_mask,
        normalize_intensity,
    )
end


text:
export Warm, compute_warm_photon_index, compute_ywarm, compute_warm_spectral_luminosity
struct Warm{T<:Float64}
    electron_energy::T
    photon_index::T
    radius::T
end

function Warm(corona::Corona; electron_energy, photon_index)
    radius = 2 * corona.radius
    return Warm(electron_energy, photon_index, radius)
end
function Warm(corona, parameters)
    return Warm(
        corona,
        electron_energy = parameters.warm_electron_energy,
        photon_index = parameters.warm_photon_index,
    )
end

compute_warm_photon_index(warm) = warm.photon_index

compute_ywarm(warm) = (4 / 9 * warm.photon_index)^(-4.5)

function compute_warm_spectral_luminosity(corona, warm, energy_range)
    gamma = warm.photon_index
    kt_e = warm.electron_energy
    kt_w = compute_disk_temperature(corona.bh, warm.radius) * K_B
    kt_w_kev = kt_w * ERG_TO_KEV
    params = [gamma, kt_e, kt_w_kev, 0, 0]
    photons_per_bin = compute_compton_photons_per_bin(energy_range, params)
    mid_energies = (energy_range[1:(end - 1)] + energy_range[2:end]) / 2
    lumin_per_bin = vcat([0], photons_per_bin[2:end] .* mid_energies)
    total_lumin = sum(lumin_per_bin)
    target_luminosity = compute_disk_luminosity(corona.bh, corona.radius, warm.radius)
    ret = lumin_per_bin * (target_luminosity / total_lumin)
    ret = ret[2:end] ./ diff(energy_range)
    return vcat([0], ret)
end


text:
function lgd(n)
#
#  test polynomial suggested by Goedecker
#   Legendre polynomial of degree n
#
    p0 = [1.0]
    p1 = [1.0; 0]
    if n == 0
        p = p0
    elseif n == 1
        p = p1
    elseif n > 1
        for k = 1:n-1
            q1 = (2*k+1)*[p1;0]
            m = length(q1)-length(p0)
            q0 = k*[zeros(m);p0]
            p = (q1-q0)/(k+1)
            #
            p0 = p1
            p1 = p
        end
    end

    p
end


text:
export ElasticPropagatorParams, ElasticPropagator, ElasticSource, ElasticReceiver,
       AcousticPropagatorParams, AcousticPropagator, AcousticSource, AcousticReceiver

@with_kw mutable struct ElasticPropagatorParams
    # number of grids along x,y axis and time steps
    NX::Int64 = 101
    NY::Int64 = 641
    NSTEP::Int64 = 2000 * 2

    # size of grid cell and time step
    DELTAX::Float64 = 10.
    DELTAY::Float64 = 10.
    DELTAT::Float64 = 2.e-3 / 2

    # PML boundary conditon
    f0::Float64 = 5. # source frequency
    vp_ref::Float64 = 2000.
    USE_PML_XMIN::Bool = true
    USE_PML_XMAX::Bool = true
    USE_PML_YMIN::Bool = true
    USE_PML_YMAX::Bool = true
    NPOINTS_PML::Int64 = 12
    NPOWER::Float64 = 2.
    K_MAX_PML::Float64 = 1.
    ALPHA_MAX_PML::Float64 = 2. * π * ( f0/2.)
    Rcoef::Float64 = 0.001 # Relative reflection coefficient

    # display params
    IT_DISPLAY::Int64 = 0

    # index
end

"""
    mutable struct ElasticSource
        srci::Union{Array{Int64, 1},PyObject}
        srcj::Union{Array{Int64, 1},PyObject}
        srctype::Union{Array{Int64, 1},PyObject}
        srcv::Union{Array{Float64, 2},PyObject}
    end

Add source terms to different components according to `srctype`

- 0: `vx`
- 1: `vy`
- 2: `σxx`
- 3: `σyy`
- 4: `σxy`
"""
mutable struct ElasticSource
    srci::Union{Array{Int64, 1},PyObject}
    srcj::Union{Array{Int64, 1},PyObject}
    srctype::Union{Array{Int64, 1},PyObject}
    srcv::Union{Array{Float64, 2},PyObject}
end

mutable struct ElasticReceiver
    rcvi::Union{Array{Int64, 1},PyObject}
    rcvj::Union{Array{Int64, 1},PyObject}
    rcvtype::Union{Array{Int64, 1},PyObject}
    rcvv::Union{Array{Float64, 2},PyObject,Missing}
end

mutable struct ElasticPropagator
    param::ElasticPropagatorParams
    src::ElasticSource
    vx::PyObject
    vy::PyObject
    sigmaxx::PyObject
    sigmayy::PyObject
    sigmaxy::PyObject
    mem::PyObject
end

"""
    AcousticPropagatorParams

A struct to hold all acoustic propagator parameters.

`PropagatorKernel`: 0 - general implementation; 1 - optimized GPU implementation; 2 - optimized CPU implementation
"""
@with_kw mutable struct AcousticPropagatorParams
    # number of grids along x,y axis and time steps
    NX::Int64 = 101
    NY::Int64 = 641
    NSTEP::Int64 = 2000 * 2

    # size of grid cell and time step
    DELTAX::Float64 = 10.
    DELTAY::Float64 = 10.
    DELTAT::Float64 = 2.e-3 / 2

    # PML boundary conditon
    USE_PML_XMIN::Bool = true
    USE_PML_XMAX::Bool = true
    USE_PML_YMIN::Bool = true
    USE_PML_YMAX::Bool = true
    NPOINTS_PML::Int64 = 12
    NPOWER::Int64 = 2
    damping_x::Union{Missing,Float64} = missing
    damping_y::Union{Missing,Float64} = missing
    Rcoef::Float64 = 0.001 # Relative reflection coefficient
    vp_ref::Float64 = 1000.

    # Auxilliary Data
    Σx::Array{Float64} = []
    Σy::Array{Float64} = []
    IJ::Array{Int64} = []
    IJn::Array{Int64}  = []
    IJp::Array{Int64} = []
    IpJ::Array{Int64} = []
    IpJp::Array{Int64} = []
    IpJn::Array{Int64} = []
    InJ::Array{Int64} = []
    InJn::Array{Int64} = []
    InJp::Array{Int64} = []

    # display params
    IT_DISPLAY::Int64 = 0
    PropagatorKernel::Int64 = 0
end

mutable struct AcousticSource
    srci::Union{Array{Int32, 1},Array{Int64, 1}, PyObject}
    srcj::Union{Array{Int32, 1},Array{Int64, 1}, PyObject}
    srcv::Union{Array{Float64, 2},PyObject}
end

mutable struct AcousticReceiver
    rcvi::Union{Array{Int64, 1},PyObject}
    rcvj::Union{Array{Int64, 1},PyObject}
    rcvv::Union{Array{Float64, 2},PyObject,Missing}
end

"""
    AcousticPropagator

A structure holding all the information for acoustic simulation
"""
mutable struct AcousticPropagator
    param::AcousticPropagatorParams
    src::AcousticSource
    u::PyObject
    # auxillary variables for PML
    φ::PyObject
    ψ::PyObject
end


text:
@testset "setters" begin
    @testset "$name" for (name, T) in [("Mag", Mag), ("MagRef", () -> Arblib.radref(Arb()))]
        # Integer
        @test Arblib.set!(T(), 1) == Arblib.set!(T(), UInt(1)) == one(Mag)

        # π
        @test Float64(Arblib.set!(T(), π)) ≈ π

        # Integer times power of 2
        @test Arblib.set!(T(), 3, 4) == Arblib.set!(T(), 3 * 2^4)
    end

    @testset "$name" for (name, T) in [("Arf", Arf), ("ArfRef", () -> Arblib.midref(Arb()))]
        # UInt128 and Int128
        let x = Int128(227725055589944414706309)
            @test Arblib.set!(T(), UInt128(x)) == Arb("227725055589944414706309")
            @test Arblib.set!(T(), x) == Arb("227725055589944414706309")
            @test Arblib.set!(T(), -x) == Arb("-227725055589944414706309")
        end

        # Rational
        @test Arblib.set!(T(), 2 // 3) == Arf(2) / Arf(3)
        @test Arblib.set!(T(), UInt128(2) // 3) == Arf(2) / Arf(3)
        @test Arblib.set!(T(), Int128(2) // 3) == Arf(2) / Arf(3)
        @test Arblib.set!(T(), BigInt(2) // 3) == Arf(2) / Arf(3)
        # Test that it works for big numerators and denominators
        let x = UInt128(typemax(Int)) + 1
            @test Arblib.set!(T(), x // 2) == Arf(x) / 2
            @test Arblib.set!(T(), 1 // x) == inv(Arf(x))
            @test Arblib.set!(T(), x // (x + 1)) == Arf(x) / Arf(x + 1)
            @test Arblib.set!(T(), Int128(x) // 2) == Arf(x) / 2
            @test Arblib.set!(T(), 1 // Int128(x)) == inv(Arf(x))
            @test Arblib.set!(T(), Int128(x) // (Int128(x) + 1)) == Arf(x) / Arf(x + 1)
            @test Arblib.set!(T(), BigInt(x) // 2) == Arf(x) / 2
            @test Arblib.set!(T(), 1 // BigInt(x)) == inv(Arf(x))
            @test Arblib.set!(T(), BigInt(x) // (BigInt(x) + 1)) == Arf(x) / Arf(x + 1)
        end
    end

    @testset "$name" for (name, T) in
                         [("Arb", Arb), ("ArbRef", () -> Arblib.realref(Acb()))]
        # MagLike, BigFloat
        @test Arblib.set!(T(), Mag(2)) == Arb(2)
        @test Arblib.set!(T(), BigFloat(2)) == Arb(2)

        # Check that aliasing works
        x = T()
        Arblib.set!(Arblib.radref(x), 1)
        @test Arblib.set!(x, Arblib.radref(x)) == one(Arb)

        # UInt128, Int128, BigInt
        let x = Int128(227725055589944414706309)
            @test Arblib.set!(T(), UInt128(x)) == Arb("227725055589944414706309")
            @test Arblib.set!(T(), x) == Arb("227725055589944414706309")
            @test Arblib.set!(T(), -x) == Arb("-227725055589944414706309")
            @test Arblib.set!(T(), BigInt(x)) == Arb("227725055589944414706309")
            @test Arblib.set!(T(), -BigInt(x)) == Arb("-227725055589944414706309")
        end

        # Rational
        @test Arblib.set!(T(), 1 // 2) == one(Arb) / 2
        @test Arblib.set!(T(), UInt128(1) // 2) == one(Arb) / 2
        @test Arblib.set!(T(), Int128(1) // 2) == one(Arb) / 2
        @test Arblib.set!(T(), BigInt(1) // 2) == one(Arb) / 2
        # Test that it works for big numerators and denominators
        let x = UInt128(typemax(Int)) + 1
            @test isequal(Arblib.set!(T(), x // 2), Arb(x) / 2)
            @test isequal(Arblib.set!(T(), 1 // x), inv(Arb(x)))
            @test isequal(Arblib.set!(T(), x // (x + 1)), Arb(x) / Arb(x + 1))
            @test isequal(Arblib.set!(T(), Int128(x) // 2), Arb(x) / 2)
            @test isequal(Arblib.set!(T(), 1 // Int128(x)), inv(Arb(x)))
            @test isequal(
                Arblib.set!(T(), Int128(x) // (Int128(x) + 1)),
                Arb(x) / Arb(x + 1),
            )
            @test isequal(Arblib.set!(T(), BigInt(x) // 2), Arb(x) / 2)
            @test isequal(Arblib.set!(T(), 1 // BigInt(x)), inv(Arb(x)))
            @test isequal(
                Arblib.set!(T(), BigInt(x) // (BigInt(x) + 1)),
                Arb(x) / Arb(x + 1),
            )
        end

        # Irrational
        for irr in [π, ℯ, MathConstants.γ, MathConstants.catalan, MathConstants.φ]
            @test Arblib.overlaps(Arblib.set!(T(), irr), Arb(BigFloat(irr)))
        end

        # Intervals
        let x = Arb(1.5)
            # Set x to the interval [1, 2]
            Arblib.set!(Arblib.radref(x), 1, -1)

            @test Arblib.contains(Arblib.set!(T(), (Mag(1), Mag(2))), x)
            @test Arblib.contains(Arblib.set!(T(), (Arf(1), Arf(2))), x)
            @test Arblib.contains(Arblib.set!(T(), (BigFloat(1), BigFloat(2))), x)

            @test Arblib.contains(Arblib.set!(T(), (1, 2)), x)
            @test Arblib.contains(Arblib.set!(T(), (1.0, 2.0)), x)
            @test Arblib.contains(Arblib.set!(T(), (1, 2.0)), x)
        end

        # Thin interval
        @test Arblib.set!(T(), (Mag(1), Mag(1))) == one(Arb)
        @test Arblib.set!(T(), (Arf(1), Arf(1))) == one(Arb)
        @test Arblib.set!(T(), (BigFloat(1), BigFloat(1))) == one(Arb)
        @test Arblib.set!(T(), (1, 1)) == one(Arb)
        @test Arblib.set!(T(), (1.0, 1.0)) == one(Arb)
        @test Arblib.set!(T(), (1, 1.0)) == one(Arb)

        # Check handling of precision
        @test !iszero(radius(Arblib.set!(Arb(prec = 64), (BigFloat(π), BigFloat(π)))))
        @test iszero(
            radius(
                Arblib.set!(
                    Arb(prec = 64),
                    (BigFloat(π), BigFloat(π)),
                    prec = precision(BigFloat),
                ),
            ),
        )
        @test Mag(1e-20) < radius(Arblib.set!(Arb(prec = 64), (π, π))) < Mag(1e-10)
        @test Mag(1e-80) <
              radius(Arblib.set!(Arb(prec = 64), (π, π), prec = 256)) <
              Mag(1e-70)

        # Check enclosure of irrationals
        @test Arblib.contains(Arblib.set!(T(), (1, π)), Arb(π))
        @test Arblib.contains(Arblib.set!(T(), (ℯ, Arb(4))), Arb(ℯ))

        @test_throws ArgumentError Arblib.set!(T(), (Mag(2), Mag(1)))
        @test_throws ArgumentError Arblib.set!(T(), (Arf(2), Arf(1)))
        @test_throws ArgumentError Arblib.set!(T(), (BigFloat(2), BigFloat(1)))
        @test_throws ArgumentError Arblib.set!(T(), (2, 1))
        @test_throws ArgumentError Arblib.set!(T(), (2.0, 1.0))
        @test_throws ArgumentError Arblib.set!(T(), (2, 1.0))
    end

    @testset "$name" for (name, T) in [
        ("Acb", Acb),
        ("AcbRef", (args...) -> AcbRefVector([Acb(args...)])[1]),
    ]
        # Setting real part
        @test Arblib.set!(T(), one(Mag)) == one(Acb)
        @test Arblib.set!(T(), one(Mag).mag) == one(Acb)
        @test Arblib.set!(T(), one(Arf)) == one(Acb)
        @test Arblib.set!(T(), one(Arf).arf) == one(Acb)
        @test Arblib.set!(T(), one(Arb)) == one(Acb)
        @test Arblib.set!(T(), one(Arb).arb) == one(Acb)
        @test Arblib.set!(T(), one(Int)) == one(Acb)
        @test Arblib.set!(T(), one(Int128)) == one(Acb)
        @test Arblib.set!(T(), one(Float64)) == one(Acb)
        @test Arblib.set!(T(), one(BigInt)) == one(Acb)
        @test Arblib.set!(T(), one(BigFloat)) == one(Acb)
        @test isequal(Arblib.set!(T(), ℯ), Acb(Arb(ℯ)))
        @test isequal(Arblib.set!(T(), (1, 2)), Acb(Arb((1, 2))))

        # Test that imaginary part is zeroed out
        @test Arblib.set!(T(0, 1), one(Mag)) == one(Acb)
        @test Arblib.set!(T(0, 1), one(Mag).mag) == one(Acb)
        @test Arblib.set!(T(0, 1), one(Arf)) == one(Acb)
        @test Arblib.set!(T(0, 1), one(Arf).arf) == one(Acb)
        @test Arblib.set!(T(0, 1), one(Arb)) == one(Acb)
        @test Arblib.set!(T(0, 1), one(Arb).arb) == one(Acb)
        @test Arblib.set!(T(0, 1), one(Int)) == one(Acb)
        @test Arblib.set!(T(0, 1), one(Int128)) == one(Acb)
        @test Arblib.set!(T(0, 1), one(BigInt)) == one(Acb)
        @test Arblib.set!(T(0, 1), one(Float64)) == one(Acb)
        @test Arblib.set!(T(0, 1), one(BigFloat)) == one(Acb)
        @test isequal(Arblib.set!(T(0, 1), ℯ), Acb(Arb(ℯ)))
        @test isequal(Arblib.set!(T(0, 1), (1, 2)), Acb(Arb((1, 2))))

        # Large integer
        let x = Int128(227725055589944414706309)
            @test Arblib.set!(T(), UInt128(x)) == Acb(Arb("227725055589944414706309"))
            @test Arblib.set!(T(), x) == Acb(Arb("227725055589944414706309"))
            @test Arblib.set!(T(), -x) == Acb(Arb("-227725055589944414706309"))
            @test Arblib.set!(T(), BigInt(x)) == Acb(Arb("227725055589944414706309"))
            @test Arblib.set!(T(), -BigInt(x)) == Acb(Arb("-227725055589944414706309"))
        end

        # Setting real and imaginary part
        # Testing with a various mix of arguments
        @test real(Arblib.set!(T(), 1, 2)) == 1
        @test imag(Arblib.set!(T(), 1, 2)) == 2
        @test real(Arblib.set!(T(), Arf(1), Mag(2))) == 1
        @test imag(Arblib.set!(T(), Arf(1), Mag(2))) == 2
        @test isequal(real(Arblib.set!(T(), (1, 2), π)), Arb((1, 2)))
        @test isequal(imag(Arblib.set!(T(), (1, 2), π)), Arb(π))
        @test real(Arblib.set!(T(), Arf(1).arf, Mag(2).mag)) == 1
        @test imag(Arblib.set!(T(), Arf(1).arf, Mag(2).mag)) == 2

        # From complex
        @test Arblib.set!(T(), Complex(1, 2)) == Acb(1, 2)
        @test Arblib.set!(T(), Complex(BigFloat(1), BigFloat(2))) == Acb(1, 2)
        @test Arblib.set!(T(), Complex(Arf(1), Arf(2))) == Acb(1, 2)
        @test Arblib.set!(T(), Complex(Arb(1), Arb(2))) == Acb(1, 2)

        # π
        @test isequal(Arblib.set!(T(), π), Acb(Arb(π)))
    end
end


text:
struct RealCoding
    nbvar::Int
    lb::Vector{Float64}
    ub::Vector{Float64}
    nbbits::Vector{Int}
    nbbitstotal::Int
end
RealCoding(ϵ::Vector{Int}, lb, ub) = begin
    @assert length(lb)==length(ub)
    for i = 1:length(lb)
        @assert lb[i] < ub[i]
        @assert UInt128(10)^ϵ[i] <= UInt128(2)^127
    end
    nbvar = length(lb)
    nbbits = ones(Int, nbvar)
    for i = 1:nbvar
        while ((UInt128(10)^ϵ[i])*(ub[i]-lb[i]) >= UInt128(2)^nbbits[i])
            nbbits[i] += 1
        end
    end
    RealCoding(nbvar, lb, ub, nbbits, sum(nbbits))
end
RealCoding(ϵ::Int, lb, ub) = RealCoding([ϵ for _=1:length(lb)], lb, ub)

function decode(x, d::RealCoding)::Vector{Float64}

    res = zeros(d.nbvar)
    for i = 1:d.nbvar
        val = zero(UInt128)
        puis = one(UInt128)
        jstart = sum(d.nbbits[1:i])
        jend = jstart - d.nbbits[i] + 1
        for j = jstart:-1:jend
            x[j] && (val += puis)
            puis *= 2
        end
        res[i] = d.lb[i] + val * (d.ub[i] - d.lb[i]) / (UInt128(2)^d.nbbits[i] - 1)
    end
    res
end

function encode(x, d::RealCoding)::BitVector

    res = BitArray(0)
    for i = 1:d.nbvar
        tab = falses(d.nbbits[i])
        ind = 1
        puis = UInt128(2)^(d.nbbits[i] - 1)
        val = zero(UInt128)
        target = round(UInt128, (x[i] - d.lb[i]) / (d.ub[i] - d.lb[i]) * (UInt128(2)^d.nbbits[i] - 1))
        if target == UInt128(2)^d.nbbits[i] - 1
            target -= 1
        end
        while val < target
            if val + puis <= target
                val += puis
                tab[ind] = true
                val == target && break
            end
            puis ÷= 2
            ind += 1
        end
        append!(res, tab)
    end
    res
end


text:
"""
    foil_points_local(;l = CHORD_LENGTH_DEFAULT,
                                 t = THICKNESS_DEFAULT,
                                 c = CAMBER_DEFAULT)
    -> Vector{Points}

2d foil coordinates centered on 0.3·chord length from leading edge.
Default values produce output in units of length.
"""
function foil_points_local(;l = CHORD_LENGTH_DEFAULT,
                                 t = THICKNESS_DEFAULT,
                                 c = CAMBER_DEFAULT)
    # Positions along chord, one way
    px = l * (-CHORD_ZERO .+ FOIL_CHORD_POS)
    # Offset centre line
    Δy_0 = c * FOIL_CAMBER
    # Half thickness
    halft = t * FOIL_HALF_T
    # Suction side, starting at the trailing edge (useful for making a spline)
    px_s = reverse(px)
    py_s = reverse( Δy_0 + halft)
    # Pressure side, starting at the leading edge
    # (useful for making a spline which is sharp at the trailing edge)
    px_p = px
    py_p = Δy_0 - halft
    px_both = vcat(px_s, px_p)
    py_both = vcat(py_s, py_p)
    map(px_both .|> upreferred, py_both .|> upreferred) do x, y
        Point(x, y)
    end
end
"""
    foil_spline_local(;l = CHORD_LENGTH_DEFAULT,
                                 t = THICKNESS_DEFAULT,
                                 c = CAMBER_DEFAULT)
    -> 872-length Vector{Points}

2d foil spline centered on 0.3·chord length from leading edge.
Default values produce output in units of length.
"""
function foil_spline_local(;l = CHORD_LENGTH_DEFAULT,
                                 t = THICKNESS_DEFAULT,
                                 c = CAMBER_DEFAULT)
    polyfit(foil_points_local(;l = l, t = t, c = c))
end

"""
    foil_draw(p::Point ;
        α = 0°,
        rel_scale = 1,
        l = CHORD_LENGTH_DEFAULT,
        t = THICKNESS_DEFAULT,
        c = CAMBER_DEFAULT,
        backgroundcolor = colorant"white")
    -> a large array of points.
Draw an airfoil with the given parameters.

p is the position of the foil, cooresponding to the rough pressure centre and
centre of rotaion, at 0.3·chord from the leading edge.
"""
function foil_draw(p::Point ;
    α = 0°,
    rel_scale = 1,
    l = CHORD_LENGTH_DEFAULT,
    t = THICKNESS_DEFAULT,
    c = CAMBER_DEFAULT,
    backgroundcolor = colorant"white")
    f = foil_spline_local(l = rel_scale * l, t = rel_scale * t, c = rel_scale * c)
    polyrotate!(f, α)
    polymove!(f,  O, p)
    poly(f, :fill)
    luminback = lumin(backgroundcolor)
    luminfront = get_current_lumin()
    avglumin = (luminback + luminfront) / 2
    contrastcol = color_with_lumin(get_current_RGB(), avglumin)
    gsave()
    sethue(contrastcol)
    poly(f, :stroke)
    grestore()
    f
end


text:
module Plot

using RecipesBase: RecipesBase, @recipe, @series, @userplot
using LinearAlgebra: eigvals, diag, eigvecs

using Unitful: @u_str, uconvert
using UnitfulRecipes: UnitfulRecipes
using UnitfulAtomic: UnitfulAtomic

using ..NQCModels: potential, derivative, nstates, state_independent_potential
using ..AdiabaticModels: AdiabaticModel
using ..DiabaticModels: DiabaticModel

@recipe function f(x, model::AdiabaticModel)
    V = zeros(size(x))
    D = zeros(size(x))
    for i=1:length(x)
        V[i] = potential(model, hcat(x[i]))
        D[i] = derivative(model, hcat(x[i]))[1]
    end

    xguide --> "r"

    @series begin
        label := "V(r)"
        x .* u"bohr", V .* u"hartree"
    end

    @series begin
        label := "dV(r)dr"
        x .* u"bohr", D .* u"hartree / bohr"
    end
end

@recipe function f(x, model::DiabaticModel; adiabats=true, diabats=true, coupling=false, atomic=true, include_diagonal=true)
    eigs = zeros(length(x), nstates(model))
    diabatic = zeros(length(x), nstates(model))
    couplings = zeros(length(x), nstates(model), nstates(model))
    for i=1:length(x)
        V = potential(model, hcat(x[i]))
        state_independent = state_independent_potential(model, hcat(x[i]))
        if include_diagonal
            eigs[i,:] .= eigvals(V) .+ state_independent
            diabatic[i,:] .= diag(V) .+ state_independent
        else
            eigs[i,:] .= eigvals(V)
            diabatic[i,:] .= diag(V)
        end
        couplings[i,:,:] .= V
    end

    xguide := "r"
    yguide := "V(r)"

    if adiabats
        for i=1:nstates(model)
            @series begin
                linecolor := :black
                label := i==1 ? "Adiabatic" : ""
                if atomic
                    x .* u"bohr", eigs[:,i] .* u"hartree"
                else
                    uconvert.(u"Å", x .* u"bohr"), uconvert.(u"eV", eigs[:,i] .* u"hartree")
                end
            end
        end
    end

    if diabats
        for i=1:nstates(model)
            @series begin
                linecolor := "#FF1F5B"
                label := i==1 ? "Diabatic" : ""
                if atomic
                    x .* u"bohr", diabatic[:,i] .* u"hartree"
                else
                    uconvert.(u"Å", x .* u"bohr"), uconvert.(u"eV", diabatic[:,i] .* u"hartree")
                end
            end
        end
    end

    if coupling
        for i=1:nstates(model)
            for j=i+1:nstates(model)
                @series begin
                    linecolor := "#009ADE"
                    label := (i == 1 && j == 2) ? "Diabatic coupling" : ""
                    if atomic
                        x .* u"bohr", couplings[:,i,j] .* u"hartree"
                    else
                        uconvert.(u"Å", x .* u"bohr"), uconvert.(u"eV", couplings[:,i,j] .* u"hartree")
                    end
                end
            end
        end
    end
end

@userplot PlotAdiabaticGradient

@recipe function f(p::PlotAdiabaticGradient; states=1:1, atomic=true)
    x, model = p.args

    legend := false

    gradient = zeros(length(x), nstates(model), nstates(model))
    for i=1:length(x)
        V = potential(model, hcat(x[i]))
        U = eigvecs(V)
        D = derivative(model, hcat(x[i]))
        gradient[i,:,:] .= U' * D[1] * U
    end

    for i in states
        for j=i:nstates(model)
            @series begin
                linecolor := i == j ? :black : :red
                if atomic
                    x .* u"bohr", (gradient[:,j,i] .* u"hartree") .^2
                else
                    uconvert.(u"Å", x .* u"bohr"), uconvert.(u"eV^2/Å^2", (gradient[:,j,i] .* u"hartree/bohr") .^2)
                end
            end
        end
    end
end

end # module


text:
# This script generates the wrapper code for the C API to OneDNN
#
# This is largely based on the wrapper code for "Sundials":
# https://github.com/SciML/Sundials.jl/blob/master/scripts/wrap_sundials.jl

using Clang
using Clang.Generators
using Clang.LibClang.Clang_jll
using MacroTools

# This is where the wrapper gets generated.
#
# If it looks Okay, move into the `src` folder.
const outpath = joinpath(@__DIR__, "dnnl.jl")
#mkpath(outpath)

# Find the relevant headers.
# We're just looking for `dnnl.h` and `dnnl_types.h`.
const include_dir = joinpath(dirname(@__DIR__), "deps", "dnnl", "include")

const onednn_headers = [
    joinpath(include_dir, "dnnl.h"),
    joinpath(include_dir, "dnnl_threadpool.h"),
]

args = get_default_args()
push!(args, "-I$include_dir")

options = Dict("general" => Dict(
    "output_file_path" => outpath,
    "library_name" => "libdnnl",
    "use_julia_native_enum_type" => false,
    "use_deterministic_symbol" => true,
))

ctx = create_context(onednn_headers, args, options)

# Build without printing so we can do custom rewriting.
build!(ctx, BUILDSTAGE_NO_PRINTING)

# For rewriting purposes, we need to intercept all entries for
# `dnnl_dims_t` and replace them with `Ptr{dnnl_dim_t}`.
custom_pointer_types = [
    :dnnl_memory,
    :dnnl_engine,
    :dnnl_primitive_desc_iterator,
    :dnnl_primitive_desc,
    :dnnl_primitive_attr,
    :dnnl_post_ops,
    :dnnl_primitive,
    :dnnl_stream,
]

rewrite(notexpr) = Any[notexpr]
function rewrite(expr::Expr)
    # Replace objects that Clang turns into "Cvoid" into singleton structs.
    # This provides a bit more type-safety when calling the C API.
    # Also need to make "unsafe_convert" an error to avoid pointers swapping types
    # accidentally.
    if MacroTools.@capture(expr, mutable struct name_ end)
        if in(name, custom_pointer_types)
            expr = quote
                $expr

                function Base.cconvert(::Type{Ptr{$name}}, x::Ptr{$name})
                    return x
                end
                function Base.cconvert(::Type{Ptr{$name}}, x::Ptr)
                    return error("Refusing to convert $(typeof(x)) to a Ptr{$($name)}!")
                end

                function Base.cconvert(::Type{Ptr{Ptr{$name}}}, x::Ptr{Ptr{$name}})
                    return x
                end
                function Base.cconvert(::Type{Ptr{Ptr{$name}}}, x::Ptr)
                    return error("Refusing to convert $(typeof(x)) to a Ptr{Ptr{$($name)}}!")
                end
            end |> MacroTools.prettify
        end
        return [expr]

    # Check if this is a ccall
    elseif expr.head == :function
        # use the powerful macrotools to transform `dnnl_dims_t` to `ptr{dnnl_dim_t}`.
        function_body = expr.args[2]
        expr.args[2] = MacroTools.postwalk(function_body) do x
            if x == :dnnl_dims_t
                return :(Ptr{dnnl_dim_t})
            else
                return x
            end
        end
    elseif expr.head == :const
        expr = MacroTools.postwalk(expr) do ex
            if ex == :INT64_MIN
                return :(typemin(Int64))
            elseif ex == :size_t
                return :unsigned
            # Delete this definition since it depends on something that gets skipped.
            elseif in(ex, (:DNNL_RUNTIME_S32_VAL_REP, :NULL))
                return 0
            end
            return ex
        end
    end
    return [expr]
end

function rewrite!(dag::ExprDAG)
    for node in get_nodes(dag)
        new_exprs = []
        for expr in get_exprs(node)
            append!(new_exprs, rewrite(expr))
        end
        node.exprs .= new_exprs
    end
end

# HACK WARNING!!
# For some reason, the topological sorting algorithm in Clang.jl is not catching this Union
# type ...
# Manually hoist it's definition to the beginning.
function replace_JL_Ctag_89!(dag::ExprDAG)
    nodes = get_nodes(dag)
    uniondef = nodes[end-1]
    @assert isa(uniondef, ExprNode{Clang.Generators.UnionAnonymous, Clang.CLUnionDecl})
    deleteat!(nodes, length(nodes) - 1)

    found = false
    for (i, node) in enumerate(nodes)
        for expr in get_exprs(node)
            if MacroTools.@capture(expr, struct T_ args__ end)
                if T == :dnnl_memory_desc_t
                    found = true
                    break
                end
            end
        end
        if found
            insert!(nodes, i - 1, uniondef)
            break
        end
    end
    @assert found
end

replace_JL_Ctag_89!(ctx.dag)
rewrite!(ctx.dag)

# build
build!(ctx, BUILDSTAGE_PRINTING_ONLY)


text:
""" Misc non-math functions used in FreeFermionGutzwiller.jl"""

count_bonds(graph::SimpleGraph{Int64}) = ne(graph)

function get_updated_bonds(chain::GutzwillerChain,move::ExchangeMove)::SimpleGraph{Int64}
    """ Returns a new graph object with the bonds of the updated state"""
    state = get_State(chain)
    model = get_Model(chain)
    lattice = model["lattice"]

    sc = deepcopy(state.spin_config.sc)
    bonds = deepcopy(state.bonds)
    sc[move.sites[1]], sc[move.sites[2]] = sc[move.sites[2]], sc[move.sites[1]]

    for site in move.sites
        for neighbor in neighbors(lattice.graph,site)
            if sc[site] == sc[neighbor]
                rem_edge!(bonds,site,neighbor)
            else
                add_edge!(bonds,site,neighbor)
            end
        end
    end
    return bonds
end

function get_conditioned_state(
    chain, n_up,n_down, conditioning_tol=10^6)
    sites = collect(1:n_up+n_down)
    cc_up = 2*conditioning_tol
    cc_down = 2*conditioning_tol
    while (cc_up > conditioning_tol  || cc_down > conditioning_tol)
        # vectors to store the location of the up and down electrons
        R_up = zeros(Int64, n_up)
        R_down = zeros(Int64, n_down)
        sample!(sites,R_up,replace=false)
        leftover_sites = collect(setdiff(Set(sites),Set(R_up)))
        sample!(leftover_sites,R_down,replace=false)

                # get occupied wavefunctions from lattice and fermi energy
        filled_states = get_Wavefunctions(chain)

        cc_up = cond(filled_states[R_up,:])
        cc_down = cond(filled_states[R_down,:])

        if (cc_up < conditioning_tol && cc_down < conditioning_tol)
            return R_up, R_down, filled_states
        end
    end
end

function update_Rs!(state::GutzwillerState,move::ExchangeMove)
    r_up = state.r_up
    r_down = state.r_down
    bd = state.business_directory
    # use bd to find what electrons are being exchanged
    l1, l2 = move.sites[1], move.sites[2]
    electron_A = bd[l1]
    electron_B = bd[l2]

    if electron_A.spin=="up" # then A is spin up and B is spin down
        r_up[electron_A.label] = l2
        r_down[electron_B.label] = l1
    else # A is down and B is up
        r_down[electron_A.label] = l2
        r_up[electron_B.label] = l1
    end
end

function update_Spin_config!(state::GutzwillerState,move::ExchangeMove)
    sc = state.spin_config.sc
    l1 = move.sites[1]
    l2 = move.sites[2]
    sc[l1], sc[l2] = sc[l2], sc[l1]
end

function update_Bonds!(chain::GutzwillerChain,move::ExchangeMove)
    state = get_State(chain)
    tmp = get_updated_bonds(chain,move)
    state.bonds = tmp
end

function update_Determinants!(chain::GutzwillerChain,move::ExchangeMove,extras)
    # extras contains det_ratio_up and det_ratio_down

    state = get_State(chain)

    if isnothing(extras)
        trash, extras = compute_ratio(chain,move)
    end
    det_ratio_up = extras[1]
    det_ratio_down = extras[2]
    state.det_A_up *= det_ratio_up
    state.det_A_down *= det_ratio_down
end

function update_Inverses!(state::GutzwillerState,move::ExchangeMove)
    u1,v1,u2,v2 = get_update_vectors(state,move)
    # recall u1 v1 are for up, u2 v2 are for down
    t = sherman_morrison_inverse_update(state.A_inv_up,u1,v1)
    r = sherman_morrison_inverse_update(state.A_inv_down,u2,v2)

    state.A_inv_up = t
    state.A_inv_down = r
end

function update_Business_directory!(state::GutzwillerState,move::ExchangeMove)
    bd = state.business_directory
    l1 = move.sites[1]
    l2 = move.sites[2]
    bd[l1], bd[l2] = bd[l2], bd[l1]
end

function make_business_directory(r_up::Array,r_down::Array)::Array
    """ Makes the directory for the lattice. The i'th element
    is the index of the electron living on site i. indices <= N/2
    refer to up spins and indices > N/2 refer to down spins """
    ns = length([r_up ; r_down])
    bd = Array{Electron}(undef, ns)

    for i in 1:length(r_up) # assume len(r_up) = len(r_down)
        bd[r_up[i]] = Electron(i,"up")
        bd[r_down[i]] = Electron(i,"down")
    end

    return bd
end

function get_ExchangeMove(bonds::SimpleGraph{Int64})::ExchangeMove
    """ draw a random bond and return the two sites attached to that bond """
    index = rand(1:ne(bonds))
    edge = nth(edges(bonds),index) #nth is from IterTools
    move_sites = (src(edge),dst(edge))
    return ExchangeMove(move_sites)
end

function get_update_vectors(state::GutzwillerState,move::ExchangeMove)
    """ Gets vectors which perform rank 1 row update to the slater determinant
    specifically, u is zero except at column i and v_i = ϕ_i(r') - ϕ_i(r)
    u1 and v1 are for up spin, and u2 and v2 are for down spin"""
    # get electron indices from business directory
    bd = state.business_directory

    site_A = move.sites[1]
    site_B = move.sites[2]

    electron_A = bd[site_A]
    electron_B = bd[site_B]

    u1 = zeros(length(state.r_up))
    v1 = zeros(length(state.r_up))
    u2 = zeros(length(state.r_down))
    v2 = zeros(length(state.r_down))

    if electron_A.spin=="up"
        # then A is up (1) and B is down (2)
        # up electron is moving from site A to site B
        u1[electron_A.label] = 1
        v1 = state.wavefunctions[site_B,:] - state.wavefunctions[site_A,:]

        # down electron is moving from site B to site A
        u2[electron_B.label] = 1
        v2 = state.wavefunctions[site_A,:] - state.wavefunctions[site_B,:]

    else
        # A is down (2) and B is up (1)
        u2[electron_A.label] = 1
        v2 = state.wavefunctions[site_B,:] - state.wavefunctions[site_A,:]

        u1[electron_B.label] = 1
        v1 = state.wavefunctions[site_A,:] - state.wavefunctions[site_B,:]
    end

    return u1,v1,u2,v2
end


text:
using Yao

export random_circuit, pair_ring

"""
    pair_ring(n::Int) -> Vector

Pair ring.
"""
pair_ring(n::Int) = [i=>mod(i, n)+1 for i=1:n]

"""
    cnot_entangler(n::Int, pairs::Vector{Pair}) = ChainBlock

Arbitrary entangler unit, support lazy construction.
"""
cnot_entangler(n::Int, pairs) = chain(n, control(n, [ctrl], target=>X) for (ctrl, target) in pairs)

"""
    rotor(noleading::Bool=false, notrailing::Bool=false) -> MatrixBlock

`Rz(η)⋅Rx(θ)⋅Rz(ξ)`, remove the first Rz gate if `noleading == true`, remove the last Rz gate if `notrailing == true`.
"""
rotor(noleading::Bool=false, notrailing::Bool=false) = noleading ? (notrailing ? Rx(0) : chain(Rx(0), Rz(0))) : (notrailing ? chain(Rz(0), Rx(0)) : chain(Rz(0), Rx(0), Rz(0)))

"""
    rotorset(noleading::Bool=false, notrailing::Bool=false) -> ChainBlock

A sequence of rotors applied on all sites.
"""
rotorset(nbit::Int, noleading::Bool=false, notrailing::Bool=false) = chain(nbit, [put(nbit, j=>rotor(noleading, notrailing)) for j=1:nbit])

"""
A kind of widely used differentiable quantum circuit, angles in the circuit are randomely initialized.

ref:
    1. Kandala, A., Mezzacapo, A., Temme, K., Takita, M., Chow, J. M., & Gambetta, J. M. (2017).
       Hardware-efficient Quantum Optimizer for Small Molecules and Quantum Magnets. Nature Publishing Group, 549(7671), 242–246.
       https://doi.org/10.1038/nature23879.
"""
function random_circuit(nbit_measure::Int, nbit_virtual::Int, nlayer::Int, nrepeat::Int, entangler_pairs)
    nbit_used = nbit_measure + nbit_virtual
    circuit = chain(nbit_used)
    entangler = cnot_entangler(nbit_used, entangler_pairs)

    for i=1:nrepeat
        unit = chain(nbit_used)
        for j=1:nlayer
            push!(unit, rotorset(nbit_used, false, false))
            push!(unit, entangler)
        end
        push!(circuit, unit)
    end
    dispatch!(circuit, :random)
end

function model(::Val{:general}; nbit::Int, V::Int, B::Int=4096, nlayer::Int=5, pairs)
    c = random_circuit(1, V, nlayer, nbit-V, pairs) |> autodiff(:QC)
    chem = QuantumMPS(1, V, 0, c, zero_state(V+1, nbatch=B), zeros(Int, nbit))
    chem
end

max_stars_repo_name: GiggleLiu/QuantumMPS

max_issues_repo_name: GiggleLiu/QuantumMPS

max_forks_repo_name: GiggleLiu/QuantumMPS


text:
function normal2a()
    μₑ = [2.5, 3.5]
    Σₑ = [2.0 -0.5;
        -0.5 5.0]
    a = [-2.3,-20.]
    b = [4.,12.3]
    properties = Dict{String,Any}()
    properties["length"] = 2
    properties["tp"] = 0.8551938607791414
    properties["mean"] = [2.126229598541412, 3.5930224468784577]
    properties["covariance"] = [1.0 0; 0 1.0] #QQQQ
    properties["moment123"] = 2.3 #QQQQ
    (MvNormal(μₑ,Σₑ), BoxTruncationRegion(a,b),properties)
end

function normal2b()
    μₑ = [2.5, 3.5]
    Σₑ = [  3.3 0.5;
            0.5 5.0]
    a = [-5.4,-20.]
    b = [2.4,6.3]
    properties = Dict{String,Any}()
    properties["length"] = 2
    properties["tp"] = 0.43660920327458974
    properties["mean"] = [0.9734400003512856, 2.886032492774952]
    properties["covariance"] = [1.0 0; 0 1.0] #QQQQ
    properties["moment123"] = 2.3 #QQQQ
    (MvNormal(μₑ,Σₑ), BoxTruncationRegion(a,b),properties)
end

function normal2c()
    a = [-1.5, -1.4]
    b = [1.7, 2.4]
    μ₁ = 0.45; μ₂ = -0.2
    σ₁ = 1.2; σ₂ = 0.8
    ρ = 0.6
    μₑ = [μ₁, μ₂]
    Σₑ = [σ₁^2 ρ*σ₁*σ₂ ;ρ*σ₁*σ₂ σ₂^2];
    properties = Dict{String,Any}()
    properties["length"] = 2
    properties["tp"] = 0.7517763397386328
    properties["mean"] = [0.308275927584875, -0.18541888959459515]
    properties["covariance"] = [1.0 0; 0 1.0] #QQQQ
    properties["moment123"] = 2.3 #QQQQ
    (MvNormal(μₑ,Σₑ), BoxTruncationRegion(a,b),properties)
end

function normal3a()
    μₑ = [3.5,2,3.5]
    Σₑ = [ 7. 1 0 ;
           1 3.3 2  ;
           0 2 3.8 ]
    a = [-2. ,-1 ,2]
    b = [4. ,4. , 5]
    properties = Dict{String,Any}()
    properties["length"] = 3
    properties["tp"] = 0.9 #QQQQ
    properties["mean"] = [1.0,1.0] #QQQQ
    properties["covariance"] = [1.0 0; 0 1.0] #QQQQ
    properties["moment123"] = 2.3 #QQQQ
    (MvNormal(μₑ,Σₑ), BoxTruncationRegion(a,b),properties)
end

function normal4a()
    μₑ = [3.5,2,3.5,3.5]
    Σₑ = [ 7. 1 0 1 ;
        1 3.3 2 0 ;
        0 2 3.4 0 ;
        1 0 0 4 ]
    a = [-20. ,-20, -20 ,-20]
    b = [20. ,20. ,20, 20]
    properties = Dict{String,Any}()
    properties["length"] = 4
    properties["tp"] = 0.9 #QQQQ
    properties["mean"] = [1.0,1.0] #QQQQ
    properties["covariance"] = [1.0 0; 0 1.0] #QQQQ
    properties["moment123"] = 2.3 #QQQQ
    (MvNormal(μₑ,Σₑ), BoxTruncationRegion(a,b),properties)
end

function normal5a()
    μₑ = [3.5,2,3.5,3.5,5.3]
    Σₑ = [ 7. 1   0 0 1 ;
        1  3.3 2 0 0 ;
        0 2 3.4 0 0 ;
        0 0 0 4 0;
        1 0 0 0 2]
    a = [-20.,-20,-20,-20,-20]
    b = [20.,20.,20,20,20]
    properties = Dict{String,Any}()
    properties["length"] = 5
    properties["tp"] = 0.9 #QQQQ
    properties["mean"] = [1.0,1.0] #QQQQ
    properties["covariance"] = [1.0 0; 0 1.0] #QQQQ
    properties["moment123"] = 2.3 #QQQQ
    (MvNormal(μₑ,Σₑ), BoxTruncationRegion(a,b),properties)
end

distribution_generators2 = [normal2a,normal2b,normal2c]
distribution_generators3 = [normal3a]
distribution_generators4 = [normal4a]
distribution_generators5 = [normal5a]
distribution_generators = vcat(distribution_generators2,distribution_generators3,distribution_generators4,distribution_generators5)


text:
using DataDeps

register(DataDep(
  "contraceptive-method-choice",
  "http://archive.ics.uci.edu/ml/datasets/Contraceptive+Method+Choice",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/cmc/cmc.data",
  "ad2a49da55df24b061a994ecc4d90c856aba35acd05138851aa251c4f06303c9",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=10,
        feature_indices=1:9,
        categoric_indices=7,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
# Part 1: evolve until we get a repeated state, then report its "biodiversity"
dirs = CartesianIndex.([(0, -1), (1, 0), (0, 1), (-1, 0)])
biodiversity(plan) = sum([2^(i-1) for (i, v) in enumerate(plan) if v])
adjacent(plan, coord) = [try plan[coord + dir] catch BoundsError; false end for dir in dirs]

"""Evolve the plan state, returning new state. adjacency_func should return array of adjacent values"""
function evolve(plan, adjacency_func)
    new_plan = copy(plan)
    for coord in CartesianIndices(plan)
        adjacent_count = sum(adjacency_func(plan, coord))
        if plan[coord] && adjacent_count != 1
            new_plan[coord] = false
        elseif !plan[coord] && adjacent_count ∈ 1:2
            new_plan[coord] = true
        end
    end
    return new_plan
end

"""Find first repeated state, returning its biodiversity"""
function find_first_repeat(plan)
    states = Set()
    while true
        bd = biodiversity(plan)
        if bd ∈ states return bd end
        push!(states, bd)
        plan = evolve(plan, adjacent)
    end
end

plan = hcat(open(f -> collect.(readlines(f)), "day24_input.txt")...) .== '#'
println("Biodiversity of first repeated state: $(find_first_repeat(plan))")

# Part 2: recursive bugs. Count total bugs after n=200 evolutions. Our plan is 3D this time
# Third dimension of the plan is 2n+1 long, 1:n are "outer", n+1 is initial, n+2:end are "inner"
function adjacent_3d(plan_3d, coord)
    adj = []
    sx, sy, sz = size(plan_3d)
    mx, my = (sx+1)÷2, (sy+1)÷2
    # Centre squares are meaningless now
    if (coord[1], coord[2]) == (mx, my) return [false] end
    for dir in dirs
        x, y, z = (coord + CartesianIndex((dir[1], dir[2], 0))).I
        try
            if x == 0 # Past left edge, go "out" a level
                push!(adj, plan_3d[CartesianIndex((mx-1, my, z-1))])
            elseif x == sx+1 # Past right edge, go "out" a level
                push!(adj, plan_3d[CartesianIndex((mx+1, my, z-1))])
            elseif y == 0 # Past top edge, go "out" a level
                push!(adj, plan_3d[CartesianIndex((mx, my-1, z-1))])
            elseif y == sy+1 # Past bottom edge, go "out" a level
                push!(adj, plan_3d[CartesianIndex((mx, my+1, z-1))])
            elseif x == mx && y == my # Centre square, go "in" a level, find sx/sy adjacent squares
                if coord[1] == mx-1 # Left edge of inner level
                    append!(adj, plan_3d[CartesianIndices((1, 1:sy, z+1))])
                elseif coord[1] == mx+1 # Right edge of inner level
                    append!(adj, plan_3d[CartesianIndices((sx, 1:sy, z+1))])
                elseif coord[2] == my-1 # Top edge of inner level
                    append!(adj, plan_3d[CartesianIndices((1:sx, 1, z+1))])
                elseif coord[2] == my+1 # Bottom edge of inner level
                    append!(adj, plan_3d[CartesianIndices((1:sx, sy, z+1))])
                end
            else # our old fashioned simple case
                push!(adj, plan_3d[CartesianIndex((x, y, z))])
            end
        catch BoundsError
            push!(adj, false)
        end
    end
    return adj
end

n = 200
plan_3d = zeros(Bool, size(plan)..., 2n+1)
plan_3d[:, :, n+1] = plan
evolved = foldl((x, y) -> evolve(x, adjacent_3d), 1:n, init=plan_3d)
println("Total bugs after $n evolutions: $(sum(evolved))")

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
module PlutoLiveExport

import Pluto, HTTP, SHA

export live_export

_format_map = Dict(
    "html" => "/notebookexport",
    "statefile" => "/statefile"
)

# change this when the PR is merged
_export_url(nb::Pluto.Notebook, endpoint::String) = "http://localhost:1234$(endpoint)?id=$(string(nb.notebook_id))"

function save_export(format::String, root_save_path::String, nb::Pluto.Notebook, export_data::Vector{UInt8})
    static_notebook_hash = bytes2hex(SHA.sha256(nb.path))
    nbname = basename(nb.path)
    if endswith(nbname, ".jl")
        nbname = nbname[1:end-3]
    end
    format_export_path = joinpath(root_save_path, format)
    if !isdir(format_export_path)
        mkdir(format_export_path)
    end
    save_path = joinpath(format_export_path, "$(nbname) $(static_notebook_hash[1:8]).$(format)")

    # actually publish the notebook
    open(save_path, "w") do io
        write(io, export_data)
    end
end

function live_export(format::String = "html"; export_path = joinpath(homedir(), "Documents", "Pluto Exports"))
    if !isdir(export_path)
        mkdir(export_path)
    end

    # make this an alias of "statefile" because I can see them getting confused
    if format == "state"
        format = "statefile"
    end

    if !haskey(_format_map, format)
        throw(ErrorException("Unknown export format \"$(format)\". Available formats are $(keys(_format_map))"))
    end

    function on_pluto_event(e::Pluto.FileEditEvent)
        try
            response = HTTP.get(_export_url(e.notebook, _format_map[format]))
            save_export(format, export_path, e.notebook, response.body)
        catch e
            @warn "Failed to fetch notebook '$(format)' export"
            showerror(stderr, e)
        end
    end
    # waiting on https://github.com/fonsp/Pluto.jl/pull/1882
    # function on_pluto_event(e::Pluto.ServerStartEvent)

    # end
    on_pluto_event(e::Pluto.PlutoEvent) = nothing

    return (e::Pluto.PlutoEvent) -> on_pluto_event(e) # todo: make this async
end

function live_export(formats::Vector{String}; kwargs...)
    for format ∈ formats
        live_export(format; kwargs...)
    end
end

end # module

max_stars_repo_name: ctrekker/PlutoLiveExport

max_issues_repo_name: ctrekker/PlutoLiveExport


text:
#= src/pheig.jl
=#

mutable struct PartialHermitianEigen{T,Tr<:Real} <: Factorization{T}
  values::Vector{Tr}
  vectors::Matrix{T}
end
const PartialHermEigen = PartialHermitianEigen

conj!(A::PartialHermEigen) = PartialHermEigen(conj!(A.values), conj!(A.vectors))
conj(A::PartialHermEigen) = conj!(copy(A))

function convert(::Type{PartialHermEigen{T}}, A::PartialHermEigen) where T
  Tr = real(T)
  PartialHermEigen(convert(Array{Tr}, A.values), convert(Array{T}, A.vectors))
end
convert(::Type{Factorization{T}}, A::PartialHermEigen{T,<:Real}) where {T} =
  A
convert(::Type{Factorization{T}}, A::PartialHermEigen) where {T} =
  convert(PartialHermEigen{T}, A)
convert(::Type{Array}, A::PartialHermEigen) = Matrix(A)
convert(::Type{Array{T}}, A::PartialHermEigen) where {T} = convert(Array{T}, Matrix(A))

copy(A::PartialHermEigen) = PartialHermEigen(copy(A.values), copy(A.vectors))

adjoint!(A::PartialHermEigen) = A
adjoint(A::PartialHermEigen) = copy(A)
transpose!(A::PartialHermEigen) = conj!(A.vectors)
transpose(A::PartialHermEigen) = PartialHermEigen(A.values, conj(A.vectors))

Matrix(A::PartialHermEigen) = (A[:vectors]*Diagonal(A[:values]))*A[:vectors]'

function getindex(A::PartialHermEigen, d::Symbol)
  if     d == :k        return length(A.values)
  elseif d == :values   return A.values
  elseif d == :vectors  return A.vectors
  else                  throw(KeyError(d))
  end
end

ishermitian(::PartialHermEigen) = true
issymmetric(A::PartialHermEigen{T}) where {T} = isreal(A)

isreal(::PartialHermEigen{T}) where {T} = T <: Real

ndims(::PartialHermEigen) = 2

size(A::PartialHermEigen) = (size(A.vectors,1), size(A.vectors,1))
size(A::PartialHermEigen, dim::Integer) =
  dim == 1 || dim == 2 ? size(A.vectors,1) : 1

# BLAS/LAPACK multiplication/division routines

## left-multiplication

mul!(y::AbstractVector{T}, A::PartialHermEigen{T}, x::AbstractVector{T}) where {T} =
  mul!(y, A[:vectors], scalevec!(A[:values], A[:vectors]'*x))

  mul!(C::AbstractMatrix{T}, A::PartialHermEigen{T}, B::AbstractMatrix{T}) where {T} =
    mul!(C, A[:vectors], lmul!(Diagonal(A[:values]), A[:vectors]'*B))

  mul!(C::AbstractMatrix{T}, A::PartialHermEigen{T}, Bc::Adjoint{T,<:AbstractMatrix{T}}) where {T} =
    mul!(C, A[:vectors], lmul!(Diagonal(A[:values]), A[:vectors]'*Bc))
  mul!(C::AbstractMatrix{T}, A::PartialHermEigen{T}, Bt::Transpose{T,<:AbstractMatrix{T}}) where {T<:Real} =
    mul!(C, A, parent(Bt)')
  mul!!(C::AbstractMatrix{T}, A::PartialHermEigen{T}, Bt::Transpose{T,<:AbstractMatrix{T}}) where {T<:Complex} =
    mul!(C, A, conj!(parent(Bt))')  # overwrites B
  function mul!(C::AbstractMatrix{T}, A::PartialHermEigen{T}, Bt::Transpose{T,<:AbstractMatrix{T}}) where T<:Complex
    B = parent(Bt)
    size(B, 1) <= A[:k] && return mul!!(C, A, transpose(copy(B)))
    tmp = (A[:vectors]')*transpose(B)
    lmul!(Diagonal(A[:values]), tmp)
    mul!(C, A[:vectors], tmp)
  end

  mul!(C::AbstractVector{T}, Ac::Adjoint{T,<:PartialHermEigen{T}}, B::AbstractVector{T}) where {T} =
    mul!(C, parent(Ac), B)
  mul!(C::AbstractMatrix{T}, Ac::Adjoint{T,<:PartialHermEigen{T}}, B::AbstractMatrix{T}) where {T} =
    mul!(C, parent(Ac), B)
  function mul!(y::AbstractVector{T}, At::Transpose{T,<:PartialHermEigen{T}}, x::AbstractVector{T}) where T
    A = parent(At)
    tmp = transpose(A[:vectors])*x
    scalevec!(A[:values], tmp)
    mul!(y, A[:vectors], conj!(tmp))
    conj!(y)
  end
  function mul!(C::AbstractMatrix{T}, At::Transpose{T,<:PartialHermEigen{T}}, B::AbstractMatrix{T}) where T
    A = parent(At)
    tmp = transpose(A[:vectors])*B
    lmul!(Diagonal(A[:values]), tmp)
    mul!(C, A[:vectors], conj!(tmp))
    conj!(C)
  end

  mul!(C::AbstractMatrix{T}, Ac::Adjoint{T,<:PartialHermEigen{T}}, Bc::Adjoint{T,<:AbstractMatrix{T}}) where {T} =
    mul!(C, parent(Ac), Bc)
  function mul!(C::AbstractMatrix{T}, At::Transpose{T,<:PartialHermEigen{T}}, Bt::Transpose{T,<:AbstractMatrix{T}}) where T
    A = parent(At)
    B = parent(Bt)
    tmp = transpose(A[:vectors])*transpose(B)
    lmul!(Diagonal(A[:values]), tmp)
    mul!(C, A[:vectors], conj!(tmp))
    conj!(C)
  end

  ## right-multiplication

  mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, B::PartialHermEigen{T}) where {T} =
    mul!(C, rmul!(A*B[:vectors], Diagonal(B[:values])), B[:vectors]')

  mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, Bc::Adjoint{T,<:PartialHermEigen{T}}) where {T} =
    mul!(C, A, parent(Bc))
  function mul!!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, Bt::Transpose{T,<:PartialHermEigen{T}}) where T
    B = parent(Bt)
    tmp = conj!(A)*B[:vectors]
    rmul!(conj!(tmp), Diagonal(B[:values]))
    mul!(C, tmp, transpose(B[:vectors]))
  end  # overwrites A
  function mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, Bt::Transpose{T,<:PartialHermEigen{T}}) where T
    B = parent(Bt)
    size(A, 1) <= B[:k] && return mul!!(C, copy(A), Bt)
    tmp = A*conj(B[:vectors])
    rmul!(tmp, Diagonal(B[:values]))
    mul!(C, tmp, transpose(B[:vectors]))
  end

  for Adj in (:Transpose, :Adjoint)
    @eval begin
      function mul!(C::AbstractMatrix{T}, Ac::$Adj{T,<:AbstractMatrix{T}}, B::PartialHermEigen{T}) where T
        tmp = Ac * B[:vectors]
        rmul!(tmp, Diagonal(B[:values]))
        mul!(C, tmp, B[:vectors]')
      end
    end
  end

  mul!(C::AbstractMatrix{T}, Ac::Adjoint{T,<:AbstractMatrix{T}}, Bc::Adjoint{T,<:PartialHermEigen{T}}) where {T} =
    mul!(C, Ac, parent(Bc))
  function mul!(C::AbstractMatrix{T}, At::Transpose{T,<:AbstractMatrix{T}}, Bt::Transpose{T,<:PartialHermEigen{T}}) where T
    A = parent(At)
    B = parent(Bt)
    tmp = A'*B[:vectors]
    rmul!(conj!(tmp), Diagonal(B[:values]))
    mul!(C, tmp, transpose(B[:vectors]))
  end

  ## left-division (pseudoinverse left-multiplication)
  ldiv!(y::AbstractVector{T}, A::PartialHermEigen{T}, x::AbstractVector{T}) where {T} =
    mul!(y, A[:vectors], iscalevec!(A[:values], A[:vectors]'*x))
  ldiv!(C::AbstractMatrix{T}, A::PartialHermEigen{T}, B::AbstractMatrix{T}) where {T} =
    mul!(C, A[:vectors], iscale!(A[:values], A[:vectors]'*B))

  # standard operations

  ## left-multiplication

  function *(A::PartialHermEigen{TA}, B::AbstractVector{TB}) where {TA,TB}
    T = promote_type(TA, TB)
    AT = convert(PartialHermEigen{T}, A)
    BT = (T == TB ? B : convert(Array{T}, B))
    CT = Array{T}(undef, size(A,1))
    mul!(CT, AT, BT)
  end

  for Adj in (:Transpose, :Adjoint)
    @eval begin
      function mul!(Ac::$Adj{TA,<:PartialHermEigen{TA}}, B::AbstractVector{TB}) where {TA,TB}
        A = parent(Ac)
        T = promote_type(TA, TB)
        AT = convert(PartialHermEigen{T}, A)
        BT = (T == TB ? B : convert(Array{T}, B))
        CT = Array{T}(undef, size(Ac,1))
        mul!(CT, $Adj(AT), BT)
      end
    end
  end

  function *(A::PartialHermEigen{TA}, B::AbstractMatrix{TB}) where {TA,TB}
    T = promote_type(TA, TB)
    AT = convert(PartialHermEigen{T}, A)
    BT = (T == TB ? B : convert(Array{T}, B))
    CT = Array{T}(undef, size(A,1), size(B,2))
    mul!(CT, AT, BT)
  end

  for Adj in (:Transpose, :Adjoint)
    @eval begin
      function *(A::PartialHermEigen{TA}, Bc::$Adj{TB,<:AbstractMatrix{TB}}) where {TA,TB}
        B = parent(Bc)
        T = promote_type(TA, TB)
        AT = convert(PartialHermEigen{T}, A)
        BT = (T == TB ? B : convert(Array{T}, B))
        CT = Array{T}(undef, size(A,1), size(Bc,2))
        mul!(CT, AT, $Adj(BT))
      end
      function *(Ac::$Adj{TA,<:PartialHermEigen{TA}}, B::AbstractMatrix{TB}) where {TA,TB}
        A = parent(Ac)
        T = promote_type(TA, TB)
        AT = convert(PartialHermEigen{T}, A)
        BT = (T == TB ? B : convert(Array{T}, B))
        CT = Array{T}(undef, size(Ac,1), size(B,2))
        mul!(CT, $Adj(AT), BT)
      end
      function *(Ac::$Adj{TA,<:PartialHermEigen{TA}}, Bc::$Adj{TB,<:AbstractMatrix{TB}}) where {TA,TB}
        A = parent(Ac)
        B = parent(Bc)
        T = promote_type(TA, TB)
        AT = convert(PartialHermEigen{T}, A)
        BT = (T == TB ? B : convert(Array{T}, B))
        CT = Array{T}(undef, size(Ac,1), size(Bc,2))
        mul!(CT, $Adj(AT), $Adj(BT))
      end
    end
  end

  function *(A::AbstractMatrix{TA}, B::PartialHermEigen{TB}) where {TA,TB}
    T = promote_type(TA, TB)
    AT = (T == TA ? A : convert(Array{T}, A))
    BT = convert(PartialHermEigen{T}, B)
    CT = Array{T}(undef, size(A,1), size(B,2))
    mul!(CT, AT, BT)
  end

  ## right-multiplication
  for Adj in (:Transpose, :Adjoint)
    @eval begin
      function *(Ac::$Adj{TA,<:AbstractMatrix{TA}}, B::PartialHermEigen{TB}) where {TA,TB}
        A = parent(Ac)
        T = promote_type(TA, TB)
        AT = (T == TA ? A : convert(Array{T}, A))
        BT = convert(PartialHermEigen{T}, B)
        CT = Array{T}(undef, size(Ac,1), size(B,2))
        mul!(CT, $Adj(AT), BT)
      end
      function *(A::AbstractMatrix{TA}, Bc::$Adj{TB,<:PartialHermEigen{TB}}) where {TA,TB}
        B = parent(Bc)
        T = promote_type(TA, TB)
        AT = (T == TA ? A : convert(Array{T}, A))
        BT = convert(PartialHermEigen{T}, B)
        CT = Array{T}(undef, size(A,1), size(Bc,2))
        mul!(CT, AT, $Adj(BT))
      end
      function *(Ac::$Adj{TA,<:AbstractMatrix{TA}}, Bc::$Adj{TB,<:PartialHermEigen{TB}}) where {TA,TB}
        A = parent(Ac)
        B = parent(Bc)
        T = promote_type(TA, TB)
        AT = (T == TA ? A : convert(Array{T}, A))
        BT = convert(PartialHermEigen{T}, B)
        CT = Array{T}(undef, size(Ac,1), size(Bc,2))
        mul!(CT, $Adj(AT), $Adj(BT))
      end
    end
  end

## left-division
function \(A::PartialHermEigen{TA}, B::AbstractVector{TB}) where {TA,TB}
  T = promote_type(TA, TB)
  AT = convert(PartialHermEigen{T}, A)
  BT = (T == TB ? B : convert(Array{T}, B))
  CT = Array{T}(undef, size(A,2))
  ldiv!(CT, AT, BT)
end
function \(A::PartialHermEigen{TA}, B::AbstractMatrix{TB}) where {TA,TB}
  T = promote_type(TA, TB)
  AT = convert(PartialHermEigen{T}, A)
  BT = (T == TB ? B : convert(Array{T}, B))
  CT = Array{T}(undef, size(A,2), size(B,2))
  ldiv!(CT, AT, BT)
end

# factorization routines

function pheigfact(
    A::AbstractMatOrLinOp{T}, opts::LRAOptions=LRAOptions(T); args...) where T
  checksquare(A)
  !ishermitian(A) && error("matrix must be Hermitian")
  opts = isempty(args) ? opts : copy(opts; args...)
  V = idfact(:n, A, opts)
  Q,R = qr!(_Matrix(:c, V))
  B = R*(A[V[:sk],V[:sk]]*R')
  F = eigen!(hermitianize!(B))
  F = PartialHermEigen(F.values, F.vectors)
  kn, kp = pheigrank(F[:values], opts)
  n = size(B, 2)
  if kn + kp < n
    idx = [1:kn; n-kp+1:n]
    F.values  = F.values[idx]
    F.vectors = F.vectors[:,idx]
  end
  pheigorth!(F.values, F.vectors, opts)
  F.vectors = Q*F.vectors
  F
end

function pheigvals(
    A::AbstractMatOrLinOp{T}, opts::LRAOptions=LRAOptions(T); args...) where T
  checksquare(A)
  !ishermitian(A) && error("matrix must be Hermitian")
  opts = isempty(args) ? opts : copy(opts; args...)
  V = idfact(:n, A, opts)
  Q,R = qr!(_Matrix(:c, V))
  B = R*(A[V[:sk],V[:sk]]*R')
  v = eigvals!(hermitianize!(B))
  kn, kp = pheigrank(v, opts)
  n = size(B, 2)
  kn + kp < n && return v[[1:kn; n-kp+1:n]]
  v
end

for f in (:pheigfact, :pheigvals)
  @eval $f(A, args...; kwargs...) = $f(LinOp(A), args...; kwargs...)
end

function pheig(A, args...; kwargs...)
  F = pheigfact(A, args...; kwargs...)
  F.values, F.vectors
end

function pheigrank(w::Vector{T}, opts::LRAOptions) where T<:Real
  n = length(w)
  k = opts.rank >= 0 ? min(opts.rank, n) : n
  wmax = max(abs(w[1]), abs(w[n]))
  idx = searchsorted(w, 0)
  kn = pheigrank1(view(w,1:first(idx)-1),   opts, wmax)
  kp = pheigrank1(view(w,n:-1:last(idx)+1), opts, wmax)
  kn, kp
end
function pheigrank1(w::AbstractVector, opts::LRAOptions, wmax::T) where T<:Real
  k = length(w)
  k = opts.rank >= 0 ? min(opts.rank, k) : k
  ptol = max(opts.atol, opts.rtol*wmax)
  @inbounds for i = 2:k
    abs(w[i]) <= ptol && return i - 1
  end
  k
end

## reorthonormalize eigenvalue cluster (issue in LAPACK)
function pheigorth!(
    values::Vector{T}, vectors::Matrix, opts::LRAOptions) where T<:Real
  n = length(values)
  a = 1
  @inbounds while a <= n
    va = values[a]
    b  = a + 1
    while b <= n
      vb = values[b]
      symrelerr(va, vb) > opts.pheig_orthtol && break
      b += 1
    end
    b -= 1
    for i = a:b
      vi = view(vectors, :, i)
      for j = i+1:b
        vj = view(vectors, :, j)
        BLAS.axpy!(-dot(vi,vj), vi, vj)
      end
    end
    a = b + 1
  end
end


text:
    """
    `MultiPhotonIonization.oneElectronAmplitude(fOrbital::Orbital, omega2::Float64, mp2::EmMultipole, symx::LevelSymmetry,
                                                                   omega1::Float64, mp1::EmMultipole, iOrbital::Orbital,
                                                                   gauge::UseGauge, orbitals::Dict{Subshell, Orbital}, grid::Radial.Grid)`
        ... to compute the two-photon ionization amplitudes for the given initial and final orbitals and quantum numbers (symmetries)
            of the intermediate partial waves. An  amp:ComplexF64  is returned.
    """
    function  oneElectronAmplitude(fOrbital::Orbital, omega2::Float64, mp2::EmMultipole, symx::LevelSymmetry,
                                                      omega1::Float64, mp1::EmMultipole, iOrbital::Orbital,
                                                      gauge::UseGauge, orbitals::Dict{Subshell, Orbital}, grid::Radial.Grid)
        amp = 0.0im

        for (k, vOrbital) in orbitals
            if  LevelSymmetry(k) != symx    continue    end
            if      gauge == UseCoulomb     gaugex = Basics.Coulomb
            elseif  gauge == UseBabushkin   gaugex = Basics.Babushkin
            end
            amp = amp + InteractionStrength.MbaAbsorptionCheng(mp2, gaugex, omega2, fOrbital, vOrbital, grid) *
                        InteractionStrength.MbaAbsorptionCheng(mp1, gaugex, omega1, vOrbital, iOrbital, grid) /
                        (iOrbital.energy + omega1 - vOrbital.energy)
        end

        return( amp )
    end

    """
    `MultiPhotonIonization.oneElectronComputeTwoPhotonLine(iState::Subshell, multipoles::Array{EmMultipole}, gauges::Array{UseGauge},
                                                           omega1::Float64, omega2::Float64, pqnMax::Int64,
                                                           orbitals::Dict{Subshell, Orbital}, meanPot::Radial.Potential; output=true)`
        ... to compute the multiphoton transition amplitudes and all properties as requested by the given settings.
            A list of lines::Array{MultiPhotonIonization.Lines} is returned.
    """
    function  oneElectronComputeTwoPhotonLine(iState::Subshell, multipoles::Array{EmMultipole}, gauges::Array{UseGauge},
                                              omega1::Float64, omega2::Float64, pqnMax::Int64,
                                              orbitals::Dict{Subshell, Orbital}, meanPot::Radial.Potential; output=false)
        println("")
        printstyled("MultiPhotonIonization.oneElectronComputeTwoPhotonLine(): The computation starts now ... \n", color=:light_green)
        printstyled("---------------------------------------------------- ---------------------------------- \n", color=:light_green)
        println("")
        #
        epsilon = omega1 + omega2 + orbitals[iState].energy
        println("Energy i-orbital: $(orbitals[iState].energy) a.u.")
        println("omega1:            $(omega1) a.u.")
        println("omega2:            $(omega2) a.u.")
        println("Excess energy:      $epsilon  a.u.  \n")

        symi = LevelSymmetry(iState)
        #
        # Determine and compute all two-photon amplitudes for the given input
        for  mp1 in multipoles
            for  mp2 in multipoles
                symxList = AngularMomentum.allowedMultipoleSymmetries(symi, mp1)
                for  symx in symxList
                    symfList = AngularMomentum.allowedMultipoleSymmetries(symx, mp2)
                    for symf in symfList
                        for gauge in gauges
                            iOrbital = orbitals[iState]
                            # Generate a proper continuum orbital with symmetry symf and energy  epsilon
                            nrContinuum = Continuum.gridConsistency(epsilon, meanPot.grid)
                            settings    = Continuum.Settings(false, nrContinuum)
                            fOrbitalPhs = Continuum.generateOrbitalLocalPotential(epsilon, Subshell(1001, symf), meanPot, settings)
                            fOrbital    = fOrbitalPhs[1]
                            @show fOrbital
                            amp = MultiPhotonIonization.oneElectronAmplitude(fOrbital, omega2, mp2, symx, omega1, mp1,
                                                                             iOrbital, gauge, orbitals, meanPot.grid)
                            println("$mp1  $mp2  $symi  $symx  $symf  $gauge  $amp")
                        end
                    end
                end
            end
        end
        #
        if    output    return( true )
        else            return( nothing )
        end
    end


text:
#=

  Euler #25 in Julia.

  Problem 25:
  """
  The Fibonacci sequence is defined by the recurrence relation:

     Fn = Fn1 + Fn2, where F1 = 1 and F2 = 1.

  Hence the first 12 terms will be:

     F1 = 1
     F2 = 1
     F3 = 2
     F4 = 3
     F5 = 5
     F6 = 8
     F7 = 13
     F8 = 21
     F9 = 34
     F10 = 55
     F11 = 89
     F12 = 144

  The 12th term, F12, is the first term to contain three digits.

  What is the first term in the Fibonacci sequence to contain 1000 digits?")
  """

  This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using Memoization
include("Euler.jl")

# Is not faster
@memoize function fib_len_memo(n)
    return fib(BigInt(n))|>string|>length
end

function fib_len(n)
    return fibmemBigInt(n)|>string|>length
end

# Brute force: 0.09872957s
function euler25a()
    i = 0
    len = 0
    while len < 1000
        i += 1
        len = fibmemBigInt(i)|>string|>length
    end

    return i;
end

#
# Using some heuristics to find the upper limit
# (from my Picat solution).
# 0.08765257s
function euler25b()
    target = 1000;
    foundUpper = 0;
    i = 1;
    fibLen = 0;
    step = 43;
    # Get the upper limit
    while fibLen < target && foundUpper == 0
        fibLen = fib_len(step*i);
        if fibLen > target
            foundUpper = i
            break
        end
        i += 1
    end

    # Now check all numbers from Step*(FoundUpper-1) .. Step*FoundUpper
    # The target must be in that interval.
    f = step*(foundUpper-1)
    fibLen = fib_len(f);
    while fibLen < target && f <= step*foundUpper
        f += 1
        fibLen = fib_len(f)
    end

    return f
end

# run_euler(euler25a)
run_euler(euler25b)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
type Token
  word::UTF8String
  wordid::Int
  charids::Vector{Int}
  tagid::Int
end

function read_conll(path, append::Bool, worddict::Dict, chardict::Dict, tagdict::Dict)
  doc = Vector{Token}[]
  sent = Token[]
  unkword = worddict["UNKNOWN"]
  for line in open(readlines, path)
    line = chomp(line)
    if length(line) == 0
      push!(doc, sent)
      sent = Token[]
    else
      items = split(line, '\t')
      word = replace(items[2], r"[0-9]", '0') |> UTF8String
      wordid = begin
        w = lowercase(word)
        get(worddict, w, unkword)
      end
      chars = convert(Vector{Char}, word)
      charids = map(chars) do c
        append ? get!(chardict, c, length(chardict)+1) : chardict[c]
      end
      tag = items[5]
      tagid = append ? get!(tagdict, tag, length(tagdict)+1) : tagdict[tag]
      token = Token(word, wordid, charids, tagid)
      push!(sent, token)
    end
  end
  doc
end

function read_wordlist(path)
  d = Dict()
  for l in open(readlines, path)
    get!(d, chomp(l), length(d)+1)
  end
  d
end

function accuracy(golds::Vector{Int}, preds::Vector{Int})
  @assert length(golds) == length(preds)
  correct = 0
  total = 0
  for i = 1:length(golds)
    golds[i] == preds[i] && (correct += 1)
    total += 1
  end
  correct / total
end

max_stars_repo_name: hshindo/POSTagging-benchmark

max_issues_repo_name: hshindo/POSTagging-benchmark

max_forks_repo_name: hshindo/POSTagging-benchmark


text:
immutable Biweight{T<:Real} <: ContinuousUnivariateDistribution
    μ::T
    σ::T

    (::Type{Biweight{T}}){T}(μ::T, σ::T) = (@check_args(Biweight, σ > zero(σ)); new{T}(μ, σ))
end

Biweight{T<:Real}(μ::T, σ::T) = Biweight{T}(μ, σ)
Biweight(μ::Real, σ::Real) = Biweight(promote(μ, σ)...)
Biweight(μ::Integer, σ::Integer) = Biweight(Float64(μ), Float64(σ))
Biweight(μ::Real) = Biweight(μ, 1.0)
Biweight() = Biweight(0.0, 1.0)

@distr_support Biweight d.μ - d.σ d.μ + d.σ

## Parameters
params(d::Biweight) = (d.μ, d.σ)
@inline partype{T<:Real}(d::Biweight{T}) = T

## Properties
mean(d::Biweight) = d.μ
median(d::Biweight) = d.μ
mode(d::Biweight) = d.μ

var(d::Biweight) = d.σ^2 / 7
skewness{T<:Real}(d::Biweight{T}) = zero(T)
kurtosis{T<:Real}(d::Biweight{T}) = T(1)/21 - 3

## Functions
function pdf{T<:Real}(d::Biweight{T}, x::Real)
    u = abs(x - d.μ) / d.σ
    u >= 1 ? zero(T) : (15//16) * (1 - u^2)^2 / d.σ
end

function cdf{T<:Real}(d::Biweight{T}, x::Real)
    u = (x - d.μ) / d.σ
    u <= -1 ? zero(T) :
    u >= 1 ? one(T) :
    (u + 1)^3/16 * @horner(u,8,-9,3)
end

function ccdf{T<:Real}(d::Biweight{T}, x::Real)
    u = (d.μ - x) / d.σ
    u <= -1 ? one(T) :
    u >= 1 ? zero(T) :
    (u + 1)^3/16 * @horner(u,8,-9,3)
end

@quantile_newton Biweight

function mgf{T<:Real}(d::Biweight{T}, t::Real)
    a = d.σ*t
    a2 = a^2
    a == 0 ? one(T) :
    15exp(d.μ * t) * (-3cosh(a) + (a + 3/a) * sinh(a)) / (a2^2)
end

function cf{T<:Real}(d::Biweight{T}, t::Real)
    a = d.σ * t
    a2 = a^2
    a == 0 ? one(T)+zero(T)*im :
    -15cis(d.μ * t) * (3cos(a) + (a - 3/a) * sin(a)) / (a2^2)
end


text:
export expected_policy_values, nash_conv

function expected_policy_values(π::AbstractPolicy, env::AbstractEnv)
    if is_terminated(env)
        [reward(env, p) for p in players(env) if p != chance_player(env)]
    elseif current_player(env) == chance_player(env)
        vals = [0.0 for p in players(env) if p != chance_player(env)]
        for (a, pₐ) in zip(action_space(env), prob(env))
            if pₐ > 0
                vals .+= pₐ .* expected_policy_values(π, child(env, a))
            end
        end
        vals
    else
        vals = [0.0 for p in players(env) if p != chance_player(env)]
        actions = action_space(env)
        probs = prob(π, env)
        @assert length(actions) == length(probs)

        for (a, p) in zip(actions, probs)
            if p > 0 #= ignore illegal action =#
                vals .+= p .* expected_policy_values(π, child(env, a))
            end
        end
        vals
    end
end

function nash_conv(π, env; is_reduce = true, kw...)
    e = copy(env)
    RLBase.reset!(e)

    σ′ = [
        best_response_value(BestResponsePolicy(π, e, i; kw...), e)
        for i in players(e) if i != chance_player(e)
    ]

    σ = expected_policy_values(π, e)
    if is_reduce
        mapreduce(-, +, σ′, σ)
    else
        map(-, σ′, σ)
    end
end


text:
using Kinetic, Plots, LinearAlgebra, JLD2, Flux
using KitBase.ProgressMeter: @showprogress
using Flux: onecold
pyplot()

cd(@__DIR__)
include("recon.jl")

function recon_pdf(ks, prim, swx, swy)
    Mu, Mv, Mxi, _, _1 = gauss_moments(prim, ks.gas.K)
    a = pdf_slope(prim, swx, ks.gas.K)
    b = pdf_slope(prim, swy, ks.gas.K)
    sw = -prim[1] .* (moments_conserve_slope(a, Mu, Mv, Mxi, 1, 0) .+ moments_conserve_slope(b, Mu, Mv, Mxi, 0, 1))
    A = pdf_slope(prim, sw, ks.gas.K)
    tau = vhs_collision_time(prim, ks.gas.μᵣ, ks.gas.ω)

    return chapman_enskog(ks.vs.u, ks.vs.v, prim, a, b, A, tau)
end

begin
    set = Setup(
        case = "cylinder",
        space = "2d2f2v",
        boundary = ["maxwell", "extra", "mirror", "mirror"],
        limiter = "minmod",
        cfl = 0.5,
        maxTime = 10.0, # time
    )
    ps = CSpace2D(1.0, 6.0, 30, 0.0, π, 50, 1, 1)
    vs = VSpace2D(-10.0, 10.0, 48, -10.0, 10.0, 48)
    gas = Gas(Kn = 1e-2, Ma = 5.0, K = 1.0)

    prim0 = [1.0, 0.0, 0.0, 1.0]
    prim1 = [1.0, gas.Ma * sound_speed(1.0, gas.γ), 0.0, 1.0]
    fw = (args...) -> prim_conserve(prim1, gas.γ)
    ff = function(args...)
        prim = conserve_prim(fw(args...), gas.γ)
        h = maxwellian(vs.u, vs.v, prim)
        b = h .* gas.K / 2 / prim[end]
        return h, b
    end
    bc = function(x, y)
        if abs(x^2 + y^2 - 1) < 1e-3
            return prim0
        else
            return prim1
        end
    end
    ib = IB2F(fw, ff, bc)

    ks = SolverSet(set, ps, vs, gas, ib)
end

ctr, a1face, a2face = init_fvm(ks, ks.ps)
@load "kn3.jld2" ctr
@load "../nn.jld2" nn

begin
    sol = zeros(ks.ps.nr, ks.ps.nθ, 4)
    for i in axes(sol, 1), j in axes(sol, 2)
        sol[i, j, :] .= ctr[i, j].prim
        sol[i, j, end] = 1 / sol[i, j, end]
    end
    contourf(
        ps.x[1:ks.ps.nr, 1:ks.ps.nθ],
        ps.y[1:ks.ps.nr, 1:ks.ps.nθ],
        sol[:, :, 4],
        ratio = 1,
    )
end

# detector
rmap = zeros(ks.pSpace.nr, ks.pSpace.nθ)
for j = 1:ks.pSpace.nθ
    for i = 2:ks.pSpace.nr
        dx1 = (ctr[i+1, j].w - ctr[i-1, j].w) / (ks.ps.x[i+1, j] - ks.ps.x[i-1, j])
        dy1 = (ctr[i+1, j].w - ctr[i-1, j].w) / (ks.ps.y[i+1, j] - ks.ps.y[i-1, j])
        dx2 = (ctr[i, j+1].w - ctr[i, j-1].w) / (ks.ps.x[i, j+1] - ks.ps.x[i, j-1])
        dy2 = (ctr[i, j+1].w - ctr[i, j-1].w) / (ks.ps.y[i, j+1] - ks.ps.y[i, j-1])

        cs = (dx1, dx2) -> (dx1 + dx2) ./ 2
        swx = cs(dx1, dx2)
        swy = cs(dy1, dy2)
        sw = [sqrt(swx[i]^2 + swy[i]^2) for i in eachindex(swx)]

        τ = vhs_collision_time(ctr[i, j].prim, ks.gas.μᵣ, ks.gas.ω) * 10

        din = [ctr[i, j].w; sw; τ]

        w = (ctr[i-1, j].w .+ ctr[i, j].w) ./ 2
        regime = nn(din) |> onecold

        rmap[i, j] = regime
    end
end

rmap = zeros(ks.pSpace.nr, ks.pSpace.nθ)
for j = 1:ks.pSpace.nθ
    rmap[1, j] = 2

    for i = 2:ks.pSpace.nr
        dx1 = (ctr[i+1, j].w - ctr[i-1, j].w) / (ks.ps.x[i+1, j] - ks.ps.x[i-1, j])
        dy1 = (ctr[i+1, j].w - ctr[i-1, j].w) / (ks.ps.y[i+1, j] - ks.ps.y[i-1, j])
        dx2 = (ctr[i, j+1].w - ctr[i, j-1].w) / (ks.ps.x[i, j+1] - ks.ps.x[i, j-1])
        dy2 = (ctr[i, j+1].w - ctr[i, j-1].w) / (ks.ps.y[i, j+1] - ks.ps.y[i, j-1])

        cs = (dx1, dx2) -> (dx1 + dx2) ./ 2
        swx = cs(dx1, dx2)
        swy = cs(dy1, dy2)

        prim = ctr[i, j].prim
        fr = recon_pdf(ks, prim, swx, swy)
        L = norm((ctr[i, j].h .- fr) ./ prim[1])

        rmap[i, j] = ifelse(L <= 0.005, 1, 2)
    end
end

contourf(
    ps.x[1:ks.ps.nr, 1:ks.ps.nθ],
    ps.y[1:ks.ps.nr, 1:ks.ps.nθ],
    rmap,
    ratio = 1,
)

max_stars_repo_name: CSMMLab/Flowmachine

max_issues_repo_name: CSMMLab/Flowmachine

max_forks_repo_name: CSMMLab/Flowmachine