text:
using Monadic

"""
    @syntax_eff begin
      a = [1,2,3]
      b = [a, a*a]
      @pure a, b
    end

Extensible Effects syntax. It applies `effect` to all values which are not prefaced with `@pure`, which lifts them into the `Eff` monad.
The results are combined using `@syntax_flatmap` on the respective `Eff` monad, which finally is executed using `autorun` algorithm.

Examples
--------

```julia
autohandled_eff = @syntax_eff begin
  a = [1,2,3]
  b = [a, a*a]
  @pure a, b
end

mycustomwrapper(i::Int) = collect(1:i)
# monadic-like syntax to first apply a wrapper before interpreting code to unhandled effects
autohandled_eff = @syntax_eff mycustomwrapper begin
  a = [1,2,3]
  b = [a, a*a]
  @pure a, b
end
```

The Interface
-------------
If you want to add support for new types, you need to provide the following interface: (Vector is only an example)

core function | description
------------- | ------------
`ExtensibleEffects.eff_applies(handler::Type{<:Vector}, effectful::Vector) = true` | specify on which values the handler applies (the handler Vector applies to Vector of course)
`ExtensibleEffects.eff_pure(handler::Type{<:Vector}, value) = [value]` | wrap a plain value into the Monad of the handler, here Vector.
`ExtensibleEffects.eff_flatmap(continuation, effectful::Vector)` | apply a continuation to the current effect (here again Vector as an example). The key difference to plain `TypeClasses.flatmap` is that `continuation` does not return a plain `Vector`, but a `Eff{Vector}`. Applying this `continuation` with a plain `map` would lead `Vector{Eff{Vector}}`. However, `eff_flatmap` needs to return an `Eff{Vector}` instead.
"""
macro syntax_eff(block::Expr)
  block = macroexpand(__module__, block)
  eff = syntax_eff(block)
  esc(:(ExtensibleEffects.autorun($eff)))
end
macro syntax_eff(wrapper, block::Expr)
  block = macroexpand(__module__, block)
  eff = syntax_eff(block, wrapper)
  esc(:(ExtensibleEffects.autorun($eff)))
end
"""
    @syntax_eff_noautorun begin
      a = [1,2,3]
      b = [a, a*a]
      @pure a, b
    end

Extensible Effects syntax which does not run `autorun` routine. As `@syntax_eff`, it applies `effect` to all values which are not prefaced with `@pure`,
which lifts them into the `Eff` monad. The results are combined using `@syntax_flatmap` on the respective `Eff` monad, which then is directly returned.

Examples
--------

```julia
unhandled_eff = @syntax_eff_noautorun begin
  a = [1,2,3]
  b = [a, a*a]
  @pure a, b
end

mycustomwrapper(i::Int) = collect(1:i)
# monadic-like syntax to first apply a wrapper before interpreting code to unhandled effects
unhandled_eff = @syntax_eff_noautorun mycustomwrapper begin
  a = [1,2,3]
  b = [a, a*a]
  @pure a, b
end
```
"""
macro syntax_eff_noautorun(block::Expr)
  block = macroexpand(__module__, block)
  esc(syntax_eff(block))
end
macro syntax_eff_noautorun(wrapper, block::Expr)
  block = macroexpand(__module__, block)
  esc(syntax_eff(block, wrapper))
end

function syntax_eff(block::Expr, extra_wrapper = :identity)
  @assert iscodeblock(block)  "expecting code block as block argument"
  monadic(
    :(ExtensibleEffects.TypeClasses.map),
    :(ExtensibleEffects.TypeClasses.flatmap),
    extra_wrapper === :identity ? :(ExtensibleEffects.effect) : :(ExtensibleEffects.effect ∘ $extra_wrapper),
    block)
end

iscodeblock(expr::Expr) = expr.head === :block
iscodeblock(other) = false

max_stars_repo_name: schlichtanders/ExtensibleEffects

max_issues_repo_name: schlichtanders/ExtensibleEffects

max_forks_repo_name: schlichtanders/ExtensibleEffects


text:
# one LIF neuron with constant input

push!(LOAD_PATH, abspath(@__DIR__,".."))
using LinearAlgebra,Statistics,StatsBase
using Plots,NamedColors ; theme(:dark)
using SparseArrays
using SpikingRNNs; const global S = SpikingRNNs
using ProgressMeter

function onesparsemat(w::Real)
  mat=Matrix{Float64}(undef,1,1) ; mat[1,1]=w
  return sparse(mat)
end

## Time parameters
dt = 1E-4
Ttot = 500.
# LIF neuron parameters
myτ = 0.5
v_th = 10.
v_r = -5.0
τrefr = 0.5 # refractoriness
τpcd = 0.2 # post synaptic current decay
myinput = 9.
myσ = 2.0
mynt = S.NTLIF(myτ,v_th,v_r,τrefr,τpcd)
myps = S.PSLIF(mynt,1)

## inputs here
# one static input
in_state = S.PSSimpleInput(S.InputSimpleOffset(myinput))
# ... plus noise
in_state_noise = S.PSSimpleInput(S.InputIndependentNormal(myσ))
##
mypop = S.Population(myps,(S.FakeConnection(),S.FakeConnection()),
    (in_state,in_state_noise))

## that's it, let's make the network
myntw = S.RecurrentNetwork(dt,(mypop,))

## Recorder
myrec1 = S.RecStateNow(myps,10,dt,Ttot;t_warmup=50.)
myrec2 = S.RecCountSpikes(myps,dt)

## Run!
times = (0:myntw.dt:Ttot)
nt = length(times)
# initial conditions
myps.state_now[1] = v_r

S.reset!.([myrec1,myrec2])

for (k,t) in enumerate(times)
  S.dynamics_step!(t,myntw)
  myrec1(t,k,myntw)
  myrec2(t,k,myntw)
end
r0 = S.get_mean_rates(myrec2,dt,Ttot)[1]
@info mean_and_std(myrec1.state_now[:])
##

# now, use Fokker-Planck theory and Richardson method
# it is hardcoded for now... maybe later I will convert it to
# a function
dv = 1E-3
v_lb = v_r - 2.0
ndv = round(Int64,(v_th - v_lb)/dv)
k_re = round(Int64,(v_r - v_lb)/dv)
ks = collect(0:ndv)
Vks = @. v_lb + ks * dv
js = fill(NaN,ndv+1)
ps = fill(NaN,ndv+1)

ksrev = collect(ndv+1:-1:2)

js[ndv+1]=1.0
ps[ndv+1]=0.0

for k in ksrev
  krond = k==k_re ? 1.0 : 0.0
  js[k-1] = js[k] - krond
  Gk = (Vks[k]-myinput)/(myσ^2)
  Hk = myτ*js[k]/(myσ^2)
  Ak = exp(dv*Gk)
  Bk = Gk==0 ? dv*Hk :  Hk/Gk*(Ak-1)
  ps[k-1] = ps[k]*Ak + Bk
end

r0 = inv(dv*sum(ps))

##
histogram(myrec1.state_now[:];nbins=100,normalize=true,
  label="simulation")
plot!(Vks,r0.*ps;color=:white,linewidth=2,leg=:topleft,
  label="Fokker-Planck")

## repeat for EIF neuron

dt = 1E-4
Ttot = 150.
# EIF neuron parameters
# τ::Float64 # time constant (membrane capacitance)
# g_l::Float64 # conductance leak
# v_expt::Float64 # exp term threshold
# steep_exp::Float64 # steepness of exponential term
# v_threshold::Float64 # spiking threshold
# v_reset::Float64 # reset after spike
# v_leak::Float64 # reversal potential for leak term
# τ_refractory::Float64 # refractory time

τ = 0.5
g_l = 1.0
v_expt=10.0
steep_exp = 1.0
v_th = 20.
v_r = -5.0
τrefr = 0.0
v_leak = 0.0

myinput = 7.0
myσ = 3.0

mynt = S.NTEIF(τ,g_l,v_expt,steep_exp,v_th,v_r,v_leak,τrefr)
myps = S.PSEIF(mynt,1)

## inputs here
# one static input
in_state = S.PSSimpleInput(S.InputSimpleOffset(myinput))
# ... plus noise
in_state_noise = S.PSSimpleInput(S.InputIndependentNormal(myσ))
##
mypop = S.Population(myps,(S.FakeConnection(),S.FakeConnection()),
    (in_state,in_state_noise))

## that's it, let's make the network
myntw = S.RecurrentNetwork(dt,(mypop,))

## Recorder
myrec1 = S.RecStateNow(myps,10,dt,Ttot;t_warmup=50.)
myrec2 = S.RecCountSpikes(myps,dt)

## Run!
times = (0:myntw.dt:Ttot)
nt = length(times)
# initial conditions
myps.state_now[1] = v_r

S.reset!.([myrec1,myrec2])

for (k,t) in enumerate(times)
  S.dynamics_step!(t,myntw)
  myrec1(t,k,myntw)
  myrec2(t,k,myntw)
end
r0_sim = S.get_mean_rates(myrec2,dt,Ttot)[1]
@info mean_and_std(myrec1.state_now[:])
##

# now, use Fokker-Planck theory and Richardson method
# it is hardcoded for now... maybe later I will convert it to
# a function
dv = 1E-3
v_lb = v_r - 2.0

psifun(V) = steep_exp * exp((V-v_expt)/steep_exp)

ndv = round(Int64,(v_th - v_lb)/dv)
k_re = round(Int64,(v_r - v_lb)/dv)
ks = collect(0:ndv)
Vks = @. v_lb + ks * dv
js = fill(NaN,ndv+1)
ps = fill(NaN,ndv+1)

ksrev = collect(ndv+1:-1:2)

js[ndv+1]=1.0
ps[ndv+1]=0.0

for k in ksrev
  krond = k==k_re ? 1.0 : 0.0
  js[k-1] = js[k] - krond
  Gk = (Vks[k]-myinput-psifun(Vks[k]))/(myσ^2)
  Hk = myτ*js[k]/(myσ^2)
  Ak = exp(dv*Gk)
  Bk = Gk==0 ? dv*Hk :  Hk/Gk*(Ak-1)
  ps[k-1] = ps[k]*Ak + Bk
end

r0_fp = inv(dv*sum(ps))

@info "rate sim $r0_sim , rate F-P $r0_fp"

##
histogram(myrec1.state_now[:];nbins=100,normed=true,
  label="simulation")
plot!(Vks,r0_fp.*ps;color=:white,linewidth=2,leg=:topleft,
  label="Fokker-Planck")

ext: jl


text:
### Functions for Pose Graph Input/Output ###
### Ted Steiner ###
### September 2014 ###

# Note: See file_format.txt for text file information.

##############
### Output ###
##############

# Write BetweenFactors to text file as EDGE2's
function writeMeasurements( filename::String, meas::Union(Vector{BetweenFactor},Array{Factor,1}) )
    f = open(filename,"w")
    writeMeasurements(f,meas)
    close(f)
    nothing
end

# Write out truth data and trajectory information
function writeTruthFile( filename::String,
                         poses::Array{Pose2D,1};
                         waypoints::Array{Int,1}=Int[],
                         anchors::Array{Array{Int,1},1}=Array{Int,1}[],
                         flags::Dict{String,Any}=Dict{String,Any}() )
    f = open(filename,"w")

    writeFlags(f,flags)
    write(f,@sprintf("\n"))

    writeWaypoints(f,waypoints)
    write(f,@sprintf("\n"))

    if length(anchors) > 0
        writeAnchors(f,anchors)
        write(f,@sprintf("\n"))
    end
    write(f,@sprintf("\n"))
    write(f,@sprintf("\n"))

    writePoses(f,poses)
    write(f,@sprintf("\n"))

    close(f)
    return nothing
end

function writeMeasurements( f::IOStream, meas::Union(Vector{BetweenFactor},Vector{Factor}) )
    for m in meas
        write(f,factor2string(m))
    end

    return nothing
end

function factor2string( m::BetweenFactor )
    @sprintf("EDGE2 %d %d %f %f %f %f 0.0 0.0 %f 0.0 %f\n",m.s,m.t,m.meas.dx,m.meas.dy,m.meas.dT,m.info.ff,m.info.ss,m.info.rr)
end

function factor2string( m::PriorFactor )
    @sprintf("PRIOR2 %d %f %f %f %f 0.0 0.0 %f 0.0 %f\n",m.pose,m.meas.dx,m.meas.dy,m.meas.dT,m.info.ff,m.info.ss,m.info.rr)
end

# Write poses to text file as Pose2d_Node's
function writePoses( f::IOStream, poses::Array{Pose2D,1} )
    count = 0
    for p in poses
        str = @sprintf("Pose2d_Node %d (%f, %f, %f)\n",count,p.x,p.y,p.theta)
        write(f,str)
        count += 1
    end

    return nothing
end

function writeAnchors( f::IOStream, anchors::Array{Int,1} )
    str = @sprintf("ANCHOR %d [%s]\n",length(anchors),vectorToString(anchors))
    write(f,str)
    return nothing
end

function writeAnchors( f::IOStream, anchors::Array{Array{Int,1},1} )
    for a in anchors
        str = @sprintf("ANCHOR %d [%s]\n",length(a),vectorToString(a))
        write(f,str)
    end
    return nothing
end

function writeWaypoints( f::IOStream, waypoints::Array{Int,1} )
    str = @sprintf("WAYPOINTS %d [%s]\n",length(waypoints),vectorToString(waypoints))
    write(f,str)
    return nothing
end

function writeFlags( f::IOStream, flags::Dict{String,Any} )
    for (name,value) in flags
        write(f,"FLAG $(uppercase(name)) $value\n")
    end
    return nothing
end

function vectorToString( v::Array{Int,1} )
    v_str = String[]
    for k = 1:length(v)
        push!(v_str,string(v[k]))
    end

    return join(v_str," ")
end

#############
### Input ###
#############

# Read pose locations from iSAM output text file
function readiSAMOutput( filename::String )
    return open(parseISAMOutput, filename)
end

function readiSAMTruth( filename::String )
    return open(parseISAMTruth, filename)
end

function readiSAMOutputCov( filename::String )
    return open(parseISAMOutputCov, filename)
end

function parseISAMOutput( f::IOStream )
    poses = Pose2D[]
    delim = Char[' ';',';'(';')']

    for line in eachline(f)
        if beginswith(line,"Pose2d_Node")
            s = split(line,delim)
            x = float(s[4])
            y = float(s[6])
            theta = float(s[8])
            push!(poses,Pose2D(x,y,theta))
        end
    end

    return poses
end

function parseISAMOutputCov( f::IOStream )
    cov_all = readdlm(f,' ',String)
    cov_s = cov_all[:,[3,6,8]]
    cov = float(cov_s)

    return cov
end

function parseISAMTruth( f::IOStream )
    # Helper function
    checkint(x) =
        try
            int(string(x))
            return true
        catch
            return false
        end

    truth = Dict{String,Any}()
    anchor_sets = Array{Int,1}[]
    delim = Char[' ';',';'(';')';'[';']';'\n']

    for line in eachline(f)
        s = split(line,delim)
        if beginswith(line,"ANCHOR")
            anchors = Int[]
            for k = 3:length(s)
                if checkint(s[k])
                    push!(anchors,int(s[k]))
                end
            end
            push!(anchor_sets,anchors)
        elseif beginswith(line,"WAYPOINTS")
            waypoints = Int[]
            for k = 3:length(s)
                if checkint(s[k])
                    push!(waypoints,int(s[k]))
                end
            end
            truth[s[1]] = waypoints
        elseif beginswith(line,"FLAG")
            s = split(line,delim)
            if checkint(s[3])
                truth[s[2]] = int(s[3])
            else
                truth[s[2]] = s[3]
            end
        end
    end

    truth["poses"] = parseISAMOutput(f)

    if length(anchor_sets) > 0
        truth["anchor_sets"] = anchor_sets
    end

    return truth
end


text:
using Dates

"""
    nauticalbells(DateTime)

    Return a string according to the "simpler system" of nautical bells
    listed in the table in Wikipedia at
    en.wikipedia.org/wiki/Ship%27s_bell#Simpler_system.
    Note the traditional time zone was determined by local sun position
    and so should be local time without daylight savings time.
"""
function nauticalbells(dt::DateTime)
    hr = hour(dt)
    mn = minute(dt)
    if hr in [00, 12, 4, 8, 16, 20]
        return mn == 00 ? "2 2 2 2" : "1"
    elseif hr in [1, 5, 9, 13, 17, 21]
        return  mn == 00 ? "2" : "2 1"
    elseif hr in [2, 6, 10, 14, 18, 22]
        return mn == 00 ? "2 2" : "2 2 1"
    elseif hr in [3, 7, 11, 15, 19, 23]
        return mn == 00 ? "2 2 2" : "2 2 2 1"
    else
        return "Gong pattern error: time $dt, hour $hr, minutes $mn"
    end
end

function nauticalbelltask()
    untilnextbell = ceil(now(), Dates.Minute(30)) - now()
    delay = untilnextbell.value / 1000
    println("Nautical bell task starting -- next bell in $delay seconds.")
    # The timer wakes its task every half hour. May drift very slightly so restart yearly.
    timer = Timer(delay; interval=1800)
    while true
        wait(timer)
        gong = nauticalbells(now())
        println("Nautical bell gong strikes ", gong)
    end
end

nauticalbelltask()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
ra = 89.014303
dec = 13.924912
ω = 37.59

α = deg2rad(ra)
δ = deg2rad(dec)
d = 1000.0 / ω

pmra = 372.72
pmdec = -483.69
radial_velocity = 0.37

η = deg2rad(58.5986320306)
αGC = deg2rad(266.4051)
δGC = - deg2rad(-28.936175)#the sign in the rotation matrix as per documentation seems to be inverted; cannot understand the root cause of the problem
dGC = 8300
xGC = [1, 0, 0]
z = 27
θ = asin(z / dGC)

r_icrs = d * [cos(α) * cos(δ), sin(α) * cos(δ), sin(δ)]

R₁ = Matrix(undef, 3, 3)
R₁[1,1] = cos(δGC)
R₁[1,2] = 0
R₁[1,3] = - sin(δGC)
R₁[2,:] = [0 1 0]
R₁[3,1] = sin(δGC)
R₁[3,2] = 0
R₁[3,3] = cos(δGC)
println("R₁:", R₁)

R₂ = Matrix(undef, 3, 3)
R₂[1,1] = cos(αGC)
R₂[1,2] = sin(αGC)
R₂[1,3] = 0
R₂[2,1] = - sin(αGC)
R₂[2,2] = cos(αGC)
R₂[2,3] = 0
R₂[3,:] = [0 0 1]
println("R₂:", R₂)

R₃ = Matrix(undef, 3, 3)
R₃[1,:] = [1 0 0]
R₃[2,1] = 0
R₃[2,2] = cos(η)
R₃[2,3] = sin(η)
R₃[3,1] = 0
R₃[3,2] = - sin(η)
R₃[3,3] = cos(η)
println("R₃:", R₃)

R = R₃ * R₁ * R₂
println("R:", R)

R = Matrix(undef, 3, 3)
R[1,1] = cos(αGC) * cos(δGC)
R[1,2] = cos(δGC) * sin(αGC)
R[1,3] = - sin(δGC)
R[2,1] = cos(αGC) * sin(δGC) * sin(η) - sin(αGC) * cos(η)
R[2,2] = sin(αGC) * sin(δGC) * sin(η) + cos(αGC) * cos(η)
R[2,3] = cos(δGC) * sin(η)
R[3,1] = cos(αGC) * sin(δGC) * cos(η) + sin(αGC) * sin(η)
R[3,2] = sin(αGC) * sin(δGC) * cos(η) - cos(αGC) * sin(η)
R[3,3] = cos(δGC) * cos(η)
println("R:", R)

H = Matrix(undef, 3, 3)
H[1,1] = cos(θ)
H[1,2] = 0
H[1,3] = sin(θ)
H[2,1] = 0
H[2,2] = 1
H[2,3] = 0
H[3,1] = - sin(θ)
H[3,2] = 0
H[3,3] = cos(θ)
println("H:", H)

r = R * r_icrs - dGC * xGC
rGC = H * r

println("distance:", d)
println("r_icrs:", r_icrs)
println("r:", r)
println("rGC:", rGC)


text:
# Dogleg is appropriate here because Jx'*Jx is positive definite, and in that
# case we only need to calculate one newton step. In the secular equation version
# we need repeaed factorizations, and that is not as easy to exploit (but maybe
# there's a good shifted one out there?)
struct NormedResiduals{Tx,Tfx,Tf}
    x::Tx
    Fx::Tfx
    F::Tf
end
has_batched_f(::NormedResiduals) = false
has_param(::NormedResiduals) = false
function value(nr::NormedResiduals, x)
    Fx = nr.F.F(nr.Fx, x)

    # this is just to grab them outside, but this hsould come from the convergence info perhaps?
    nr.x .= x
    nr.Fx .= Fx

    f = (norm(Fx)^2) / 2
    return f
end

function upto_gradient(nr::NormedResiduals, Fx, x)
    Fx, Jx = nr.F.FJ(nr.Fx, nr.Fx * x', x)

    # this is just to grab them outside, but this hsould come from the convergence info perhaps?
    nr.x .= x
    nr.Fx .= Fx

    f = (norm(Fx)^2) / 2
    Fx .= Jx' * Fx
    return f, Fx
end
function upto_hessian(nr::NormedResiduals, Fx, Jx, x)  #Fx is the gradient and Jx is the Hessian
    Fx, Jx = nr.F.FJ(nr.Fx, Jx, x)

    # this is just to grab them outside, but this hsould come from the convergence info perhaps?
    nr.x .= x
    nr.Fx .= Fx

    f = (norm(Fx)^2) / 2
    # this is the gradient
    Fx .= Jx' * Fx
    # As you may notice, this can be expensive... Because the gradient
    # is going to be very simple. May want to create a
    # special type or way to hook into trust regions here. We can exploit
    # that we only need the cauchy and the newton steps, not any shifted
    # systems. There is no need to get the full hessian. because these two
    # steps are don't need these multiplies
    # This is the Hessian
    Jx .= Jx' * Jx
    return f, Fx, Jx
end

function solve(
    prob::NEqProblem,
    x,
    approach::TrustRegion{<:Union{SR1,DBFGS,BFGS,Newton},<:Any,<:Any},
    options::NEqOptions,
)
    if !(mstyle(prob) === InPlace())
        throw(
            ErrorException(
                "solve() not defined for OutOfPlace() with Trustregion for NEqProblem",
            ),
        )
    end
    F = prob.R
    # should we wrap a Fx here so we can log F0 info here?
    # and so we can extract it at the end as well?
    # xcache = copy(x).-1
    Fx_outer = copy(x)
    x_outer = copy(x)

    normed_residual = NormedResiduals(x_outer, Fx_outer, F)
    ρ2F0 = sqrt(value(normed_residual, x_outer) * 2)
    ρF0 = norm(normed_residual.Fx, Inf)
    td = OptimizationProblem(normed_residual)
    res = solve(td, x, approach, OptimizationOptions(maxiter = options.maxiter))
    # normed_residual to get Fx
    # f0*2 is wrong because it's 2norm
    newinfo = (
        solution = solution(res),
        best_residual = Fx_outer,
        ρF0 = ρF0,
        ρ2F0 = ρ2F0,
        time = res.info.time,
        iter = res.info.iter,
    )
    return ConvergenceInfo(approach, newinfo, options)
end

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers


text:
close("all")
transD_GP.getchi2forall(opt)
ax = gcf().axes;
χ² = aem.ndatalow/2*log(aem.ndatalow) + aem.ndatahigh/2*log(aem.ndatahigh)
ax[2].plot(xlim(), [χ², χ²], "--", color="gray")
opt.xall[:] .= zall
transD_GP.plot_posterior(aem, opt, burninfrac=0.5, figsize=(10,6),
    cmappdf="cividis", qp1=0.01, qp2=0.99, nbins=50, vmaxpc=1.0)
ax = gcf().axes
ax[1].invert_xaxis()
ax[1].step(log10.(ρ[2:end]), z[2:end], color="r", alpha=0.5)
ax[1].step(log10.(ρ[2:end]), z[2:end], color="m", alpha=0.5, linestyle="--")
ax[2].set_ylim(280,0)

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
function check_symplecticity(tab::Tableau{T}; atol=16eps(T), rtol=16eps(T)) where {T}
    a, b = tab.a, tab.b
    [isapprox(b[i] * a[i,j] + b[j] * a[j,i], b[i] * b[j]; atol=atol, rtol=rtol) for i in axes(a,1), j in axes(a,2)]
end

function issymplectic(tab::Tableau; kwargs...)
    all(check_symplecticity(tab; kwargs...))
end

function check_symplecticity(tab::PartitionedTableau{T}; atol=16eps(T), rtol=16eps(T)) where {T}
    a, b = tab.q.a, tab.q.b
    ā, b̄ = tab.p.a, tab.p.b
    ([isapprox(b[i] * ā[i,j] + b̄[j] * a[j,i], b[i] * b̄[j]; atol=atol, rtol=rtol) for i in axes(a,1), j in axes(a,2)],
     [isapprox(b[i], b̄[i]; atol=atol, rtol=rtol) for i in eachindex(b,b̄)])
end

function issymplectic(tab::PartitionedTableau; kwargs...)
    all(all.(check_symplecticity(tab; kwargs...)))
end

function compute_symplecticity_error(tab::Tableau)
    a, b = tab.a, tab.b
    [b[i] * a[i,j] + b[j] * a[j,i] - b[i] * b[j] for i in axes(a,1), j in axes(a,2)]
end

function get_symplectic_conjugate_coefficients!(a̅::AbstractMatrix{T}, a::AbstractMatrix{T}, b::AbstractVector{T}) where {T}
    @assert size(a) == size(a̅)
    @assert length(b) == size(a,1) == size(a,2)

    for i in axes(a̅, 1)
        for j in axes(a̅, 2)
            a̅[i,j] = b[j] / b[i] * ( b[i] - a[j,i] )
        end
    end

    return a̅
end

get_symplectic_conjugate_coefficients(a, b) = get_symplectic_conjugate_coefficients!(zero(a), a, b)
get_symplectic_conjugate_coefficients(tab::Tableau) = Tableau(tab.name, tab.s, get_symplectic_conjugate_coefficients(tab.a, tab.b), tab.b, tab.c)

function symplecticize(tab::Tableau; name=nothing, T=Float64)
    a̅ = get_symplectic_conjugate_coefficients(tab.a, tab.b)
    Tableau{T}(name === nothing ? Symbol(tab.name, "S") : name, tab.o, (tab.a .+ a̅) ./ 2, tab.b, tab.c)
end


text:
module PearsonHash

using Random: shuffle!
export hash8, hashn, hashn!

_table = shuffle!(collect(0:255))

"""
    hash8(data::Vector{UInt8}; seed::UInt8 = zero(UInt8))
    hash8(data::AbstractString; seed::UInt8 = zero(UInt8))

Compute the Pearson 8-bit hash for `data` (optionally specify `seed`).
"""
function hash8(data::Vector{UInt8}; seed::UInt8 = zero(UInt8))
    h = seed
    for byte in data
        h = _table[xor(h, byte) + one(UInt8)]
    end

    return h
end
hash8(data::AbstractString; seed::UInt8 = zero(UInt8)) = hash8(Vector{UInt8}(data); seed = seed)

"""
    hashn!([T = UInt128], data::Vector{UInt8}; seed::UInt8 = 0)::T where T

Compute a (`8*n`)-bit hash for data where `n == sizeof(T)` (optionally specify `seed`).
*Note:* `data` will be mutated in-place. See `hashn` for version without mutation.
"""
function hashn!(::Type{T}, data::Vector{UInt8}; seed::UInt8 = zero(UInt8))::T where T
    h = zero(T)
    n = sizeof(T)
    firstbyte = data[1]
    for i in 0:(n - 1)
        data[1] = UInt8((firstbyte + i) % 256)
        h |= T(hash8(data; seed = seed)) << (i * 8)
    end

    return h
end
hashn!(data::Vector{UInt8}; seed::UInt8 = zero(UInt8)) = hashn!(UInt128, data; seed = seed)
"""
    hashn([T = UInt128], data; seed::UInt8 = zero(UInt8))::T where T

Compute a (`8*n`)-bit hash for `Vector{UInt8}(data)` where `n == sizeof(T)` (optionally specify `seed`).
*Note:* `data` is copied to prevent mutation (resulting in extra allocations).
"""
hashn(::Type{T}, data; seed::UInt8 = zero(UInt8)) where T = hashn!(T, Vector{UInt8}(data); seed = seed)
hashn(data; seed::UInt8 = zero(UInt8)) = hashn(UInt128, data; seed = seed)

end # module

max_stars_repo_name: darsnack/PearsonHash

max_issues_repo_name: darsnack/PearsonHash

max_forks_repo_name: darsnack/PearsonHash


text:
"""
    ScaleTransform(l::Real)

Transformation that multiplies the input elementwise with `l`.

# Examples

```jldoctest
julia> l = rand(); t = ScaleTransform(l); X = rand(100, 10);

julia> map(t, ColVecs(X)) == ColVecs(l .* X)
true
```
"""
struct ScaleTransform{T<:Real} <: Transform
    s::Vector{T}
end

function ScaleTransform(s::T=1.0) where {T<:Real}
    return ScaleTransform{T}([s])
end

@functor ScaleTransform

set!(t::ScaleTransform, ρ::Real) = t.s .= [ρ]

(t::ScaleTransform)(x) = first(t.s) * x

_map(t::ScaleTransform, x::AbstractVector{<:Real}) = first(t.s) .* x
_map(t::ScaleTransform, x::ColVecs) = ColVecs(first(t.s) .* x.X)
_map(t::ScaleTransform, x::RowVecs) = RowVecs(first(t.s) .* x.X)

Base.isequal(t::ScaleTransform, t2::ScaleTransform) = isequal(first(t.s), first(t2.s))

Base.show(io::IO, t::ScaleTransform) = print(io, "Scale Transform (s = ", first(t.s), ")")


text:
# ---
# title: Continuous scales
# cover: assets/continuous_scales.png
# description: Applying nonlinear transformations
# author: "[Pietro Vertechi](https://github.com/piever)"
# ---

using AlgebraOfGraphics, CairoMakie
set_aog_theme!() #src

x = 1:100
y = @. sqrt(x) + 20x + 100
df = (; x, y)
plt = data(df) *
    mapping(
        :x,
        :y => log => "√x + 20x + 100 (log scale)",
    ) * visual(Lines)
draw(plt)

#

x = 1:100
y = @. sqrt(x) + 20x + 100
df = (; x, y)
plt = data(df) *
    mapping(
        :x,
        :y => "√x + 20x + 100 (log scale)",
    ) * visual(Lines)
draw(plt, axis=(yscale=log,))

#

x = 0:100
y = @. 0.01 + x/1000
df = (; x, y)
plt = data(df) *
    mapping(
        :x,
        :y => "y",
    ) * visual(Lines)
fg = draw(plt, axis=(yscale=log,))

# save cover image #src
mkpath("assets") #src
save("assets/continuous_scales.png", fg) #src


text:
#!/usr/bin/env julia

# Generate a set of model and residual visibilities and then write them to UVHDF5 file format
# More about this format, including scripts to convert to UVFITS and CASA Measurement set,
# can be found at https://github.com/Astrochem/UVHDF5

using Pkg; Pkg.activate("DiskJockey")

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "--config"
    help = "a YAML configuration file"
    default = "config.yaml"
    "--out-model"
    help = "The file name to write out the model visibilities."
    default = "model.hdf5"
    "--out-resid"
    help = "The file name to write out the residual visibilities."
    default = "resid.hdf5"
end

parsed_args = parse_args(ARGS, s)

using DiskJockey.constants
using DiskJockey.image
using DiskJockey.model
using DiskJockey.visibilities
using DiskJockey.gridding

using HDF5

import YAML
config = YAML.load(open(parsed_args["config"]))

# The model writing is taken care of by `DJ_initialize.jl`, the RADMC-3D synthesis is taken care of by `DJ_plot_model.jl` and then this simply reads in the image, does FFT, downsamples, etc.

dvarr = DataVis("resid.hdf5")
nchan = length(dvarr)

for dset in dvarr
    # Conjugation is necessary for the SMA and ALMA
    visibilities.conj!(dset) # Swap UV convention
end

chi2s = Array{Float64}(undef, nchan)

for i=1:nchan
    dvis = dvarr[i]
    chi2s[i] = sum(abs2, dvis.invsig .* dvis.VV)
end

N = nchan * 2 * length(dvarr[1].VV)

# The data are stored in increasing frequency, so
# exclude: [1] means exclude the most redshifted channel
# whereas
# exclude : [nchan] excludes the most blueshifted channel
if haskey(config, "exclude")
    exclude = config["exclude"]
    # which channels of the dset to fit
    # keylist = filter(x->(!in(x, exclude)), Int[i for i=1:nchan])

    lam0 = lam0s[config["species"] * config["transition"]]
    # calculate the velocities corresponding to dvarr
    lams = Float64[dv.lam for dv in dvarr]
    vels = c_kms * (lams .- lam0)/lam0
    # get the mask
    vel_mask = generate_vel_mask(exclude, vels)
else
    # keylist = Int[i for i=1:nchan]
    vel_mask = trues(nchan)
end

chi2s = chi2s[vel_mask]
println("Reduced Chi^2 ", sum(chi2s)/N)


text:
function is_force_test_plots(a::AbstractDict = ENV)
    result = lowercase(strip(get(a, "PREDICTMD_TEST_PLOTS", ""))) == "true"
    return result
end

function handle_plotting_error(e::Exception, a::AbstractDict = ENV)::Nothing
    if is_force_test_plots(a)
        @warn(
            string(
                "PREDICTMD_TEST_PLOTS is true ,",
                "so rethrowing the error.",
                )
            )
        rethrow(e)
    else
        @warn(string("ignoring error:", e))
        return nothing
    end
end


text:
using Thebes, Luxor, ColorSchemes

function makemobius()
    x(u, v) = (1 + (v/2 * cos(u/2))) * cos(u)
    y(u, v) = (1 + (v/2 * cos(u/2))) * sin(u)
    z(u, v) = v/2 * sin(u/2)
    # v goes across
    # u goes round
    w = 1
    st = 2π/50
    Δ = .05
    result = []
    for u in 0:st:2π# -st
        for v in -w:Δ:w
            p1 = Point3D(x(u,      v + Δ),   y(u,      v + Δ),     z(u,       v + Δ))
            p2 = Point3D(x(u + st, v + Δ),   y(u + st, v + Δ),     z(u  + st, v + Δ))
            p3 = Point3D(x(u + st, v),       y(u + st, v),         z(u + st,  v))
            p4 = Point3D(x(u,      v),       y(u,      v),         z(u,       v))
            push!(result, [p1, p2, p3, p4])
        end
    end
    return result
end

function frame(scene, framenumber, mobius)
    background("black")
    # flyaround
    theta = rescale(framenumber, 1, scene.framerange.stop, 0, 2pi)
    perspective(400)
    eyepoint(Point3D(250(cos(theta)), 250(sin(theta)), 200))
    setline(.5)
    for (n, p) in enumerate(mobius)
        sethue(get(ColorSchemes.cyclic_wrwbw_40_90_c42_n256_s25, n/length(mobius)))
        pin(100p, gfunction  = (p3l, p2l) -> begin
            poly(p2l, :stroke, close=true)
        end)
    end
end

function  main()
    width, height = (800, 400)
    mobiusmovie = Movie(width, height, "mobiusmovie")

    mb = makemobius()
    gif = animate(mobiusmovie, [
        Scene(mobiusmovie, (s, f) -> frame(s, f, mb), 1:100)
        ],
        creategif=true,
        framerate=20,
        pathname="/tmp/mobiusmovie.gif")
    return gif
end

main()


text:
using DifferentialForms

using LinearAlgebra
using StaticArrays
using Test

@testset "Create forms D=$D R=$R" for D in 0:Dmax, R in 0:D
    T = BigRat

    N = binomial(Val(D), Val(R))
    x = Form{D,R,T}(SVector{N,T}(rand(T) for n in 1:N))
    X = typeof(x)
    x′ = X(SVector{N,T}(rand(T) for n in 1:N))
    if R == 0
        x″ = fscalar(D, rand(T))
        x″::X
    end
    if D > 0 && R == 1
        x″ = fvector(rand(T, D)...)
        x″::X
    end
    if R == D
        x″ = fpseudoscalar(D, rand(T))
        x″::X
    end

    @test Form{D,R,Float64}(x) == Form{D,R,Float64}(x.elts)

    @test Form{D,R,Float64}(collect(x)) == Form{D,R,Float64}(x)
end

@testset "Create forms from arrays D=$D R=$R" for D in 0:Dmax, R in 0:D
    T = BigRat

    fun() = rand(T)
    fun(i) = rand(T)
    fun(i, j) = rand(T)
    fun(i, j, k) = rand(T)
    fun(i, j, k, l) = rand(T)
    fun(i, j, k, l, m) = rand(T)
    y = MakeForm{D,R,T}(fun::Function)

    arr = Array{T}(undef, ntuple(r -> D, R))
    for ind in CartesianIndex(ntuple(r -> 1, R)):CartesianIndex(ntuple(r -> D, R))
        lst = SVector{R,Int}(Tuple(ind))
        issorted = true
        for r in 1:(R - 1)
            issorted &= lst[r] < lst[r + 1]
        end
        if issorted
            arr[lst...] = rand(T)
        end
    end
    z = MakeForm{D,R,T}(arr)

    subarr = R == 0 ? (@view arr[]) :
             R == 1 ? (@view arr[:]) :
             R == 2 ? (@view arr[:, :]) :
             R == 3 ? (@view arr[:, :, :]) : R == 4 ? (@view arr[:, :, :, :]) : R == 5 ? (@view arr[:, :, :, :, :]) : nothing
    subarr::SubArray
    @assert ndims(arr) == R
    @assert ndims(subarr) == R
    @assert all(size(arr, r) == D for r in 1:R)
    @assert all(size(subarr, r) == D for r in 1:R)
    @assert all(==(D), size(arr))
    @assert all(==(D), size(subarr))
    w = MakeForm{D,R,T}(subarr)
    @test w == z
end

@testset "Comparing forms D=$D R=$R" for D in 0:Dmax, R in 0:D
    T = Int

    for n in 1:100
        x = rand(Form{D,R,T})
        y = rand(Form{D,R,T})
        isequal(x, y) && continue

        @test x == x
        @test isequal(x, x)
        @test hash(x) == hash(x)
        @test !isless(x, x)

        @test x != y
        @test !isequal(x, y)
        @test hash(x) != hash(y)
        @test isless(x, y) || isless(y, x)
    end
end

@testset "Forms as collections D=$D R=$R" for D in 0:Dmax, R in 0:D
    T = Int
    x = rand(Form{D,R,T})

    N = length(x)
    for n in 1:N
        y = Base.setindex(x, x[n], n)
        @test y === x

        inds = Forms.lin2lst(Val(D), Val(R), n)::SVector{R,Int}
        y = Base.setindex(x, x[inds], inds)
        @test y === x

        tup = Tuple(inds)::NTuple{R,Int}
        y = Base.setindex(x, x[tup], tup)
        @test y === x

        y = Base.setindex(x, x[inds...], inds...)
        @test y === x
    end
end

@testset "Vector space operations on forms D=$D R=$R" for D in 0:Dmax, R in 0:D
    # Using === instead of == for comparisons to catch wrong types
    T = Rational{Int64}
    n = zero(Form{D,R,T})
    x = rand(Form{D,R,T})
    y = rand(Form{D,R,T})
    z = rand(Form{D,R,T})
    a = rand(T)
    b = rand(T)

    # Vector space

    @test iszero(n)
    @test iszero(x) === (x === n)

    @test +x === x
    @test (x + y) + z === x + (y + z)

    @test n + x === x
    @test x + n === x

    @test x + y === y + x

    @test x + (-x) === n
    @test (-x) + x === n
    @test x - y === x + (-y)

    @test conj(conj(x)) === x

    @test (a * b) * x === a * (b * x)
    @test x * (a * b) === (x * a) * b

    @test one(T) * a === a
    @test one(T) * x === x
    @test x * one(T) === x

    @test zero(T) * a === zero(T)
    @test zero(T) * x === n
    @test x * zero(T) === n

    if a != zero(T)
        @test a * inv(a) === one(T)

        @test inv(a) * (a * x) === x
        @test (x * a) * inv(a) === x

        @test inv(a) * x === a \ x
        @test x * inv(a) === x / a
    end

    @test (a + b) * x === a * x + b * x
    @test x * (a + b) === x * a + x * b

    @test a * (x + y) === a * x + a * y
    @test (x + y) * a === x * a + y * a

    # Nonlinear transformations
    @test map(+, x) === x
    @test map(+, x, y) === x + y
    @test map(+, x, y, z) === x + y + z
end

@testset "Form algebra D=$D R1=$R1 R2=$R2" for D in 0:Dmax, R1 in 0:D, R2 in 0:D
    # Using === instead of == for comparisons to catch wrong types
    T = Rational{Int64}
    e = one(Form{D,0,T})
    x = rand(Form{D,R1,T})
    y = rand(Form{D,R2,T})
    y2 = rand(Form{D,R2,T})
    a = rand(T)
    b = rand(T)

    # Multiplicative structure

    @test isone(e)
    @test !isone(2 * e)

    # units
    if D == 2 && R1 == 1 && R2 == 1
        e = unit(Form{D,0,T})
        e1 = unit(Form{D,1,T}, 1)
        e2 = unit(Form{D,1,T}, 2)
        e12 = unit(Form{D,2,T}, 1, 2)
        @test ~e === e
        @test ~e1 === e1
        @test ~e2 === e2
        @test ~e12 === -e12

        @test ⋆e === e12
        @test ⋆e1 === e2
        @test ⋆e2 === -e1
        @test ⋆e12 === e

        @test e ∧ e === e
        @test e ∧ e1 === e1
        @test e ∧ e2 === e2
        @test e ∧ e12 === e12
        @test e1 ∧ e === e1
        @test e2 ∧ e === e2
        @test e12 ∧ e === e12
        @test e1 ∧ e1 === 0 * e12
        @test e1 ∧ e2 === e12
        @test e2 ∧ e1 === -e12
        @test e2 ∧ e2 === 0 * e12

        @test e ∨ e12 === e
        @test e1 ∨ e1 === 0 * e
        @test e2 ∨ e2 === 0 * e
        @test e1 ∨ e2 === e
        @test e1 ∨ e12 === e1
        @test e2 ∨ e12 === e2
        @test e12 ∨ e12 === e12

        @test e ⋅ e === e
        @test e1 ⋅ e1 === e
        @test e2 ⋅ e2 === e
        @test e1 ⋅ e2 === 0 * e
        @test e12 ⋅ e12 === e

        @test e × e === e12
        @test e × e1 === e2
        @test e × e2 === -e1
        @test e × e12 === e
        @test e1 × e1 === 0 * e
        @test e2 × e2 === 0 * e
        @test e1 × e2 === e
    end

    # various duals
    cycle_basis(x)::typeof(x)
    x′ = x
    for d in 1:D
        x′ = cycle_basis(x′)
    end
    @test x′ === x
    reverse_basis(x)::typeof(x)
    @test reverse_basis(reverse_basis(x)) === x
    @test reverse_basis(cycle_basis(reverse_basis(cycle_basis(x)))) === x

    @test ~~x === x
    @test a * ~x === ~(a * x)

    @test conj(conj(x)) === x
    @test a * conj(x) === conj(a * x)

    @test ⋆⋆x === bitsign(R1 * (D - R1)) * x
    @test inv(⋆)(⋆x) === x
    @test ⋆inv(⋆)(x) === x
    @test ⋆⋆ ⋆ ⋆x === x
    @test a * ⋆x === ⋆(a * x)

    # exterior product: x ∧ y
    @test ∧(x) === x

    if R1 + R2 <= D
        (x ∧ y)::Form{D,R1 + R2}
    end

    for R3 in 0:D
        z = rand(Form{D,R3,T})
        if R1 + R2 + R3 <= D
            @test (x ∧ y) ∧ z === x ∧ (y ∧ z)
            @test x ∧ y ∧ z === (x ∧ y) ∧ z
            @test ∧(x, y, z) === (x ∧ y) ∧ z
            @test ∧(SVector{0,typeof(x)}()) === one(x)
            @test ∧(SVector(x)) === x
            @test ∧(SVector(x, y)) === x ∧ y
            @test ∧(SVector(x, y, z)) === (x ∧ y) ∧ z
            @test ∧(SVector{0,typeof(x)}()) ∧ x === x
            @test x ∧ ∧(SVector{0,typeof(x)}()) === x
        end
    end

    @test e ∧ x === x
    @test x ∧ e === x

    if R1 + R2 <= D
        @test x ∧ zero(y) === zero(x ∧ y)
        @test zero(y) ∧ x === zero(y ∧ x)

        @test a * (x ∧ y) === x ∧ (a * y)
        @test x ∧ (y + y2) === x ∧ y + x ∧ y2

        @test x ∧ y === bitsign(R1 * R2) * (y ∧ x)
    end

    # regressive product: ⋆(x ∨ y) = ⋆x ∧ ⋆y
    @test ∨(x) === x
    Rvee = D - ((D - R1) + (D - R2))
    @test Rvee == R1 + R2 - D
    if 0 <= Rvee <= D
        (x ∨ y)::Form{D,Rvee}
        @test ⋆(x ∨ y) === ⋆x ∧ ⋆y
    end

    for R3 in 0:D
        z = rand(Form{D,R3,T})
        Rvee3 = D - ((D - R1) + (D - R2) + (D - R3))
        if 0 <= Rvee3 <= D
            (x ∨ y ∨ z)::Form{D,Rvee3}
            @test ⋆(x ∨ y ∨ z) === ⋆x ∧ ⋆y ∧ ⋆z
            @test x ∨ y ∨ z === (x ∨ y) ∨ z
            @test ∨(x, y, z) === (x ∨ y) ∨ z
            @test ∨(SVector{0,typeof(x)}()) === ⋆one(x)
            @test ∨(SVector(x)) === x
            @test ∨(SVector(x, y)) === x ∨ y
            @test ∨(SVector(x, y, z)) === (x ∨ y) ∨ z
            @test ∨(SVector{0,typeof(x)}()) ∨ x === x
            @test x ∨ ∨(SVector{0,typeof(x)}()) === x
        end
    end

    # dot product: x ⋅ y = x ∨ ⋆y   (right contraction)
    Rdot = D - ((D - R1) + R2)
    @test Rdot == R1 - R2
    if 0 <= Rdot <= D
        (x ⋅ y)::Form{D,Rdot}
        @test x ⋅ y === x ∨ ⋆y
    end

    norm2(x)::T
    @test norm2(x) >= 0
    @test norm2(a * x) == norm2(a) * norm2(x)
    # @test norm2(y + y2) <= norm2(y) + norm2(y2)
    norm(x)::float(T)
    @test norm(x) ≈ sqrt(norm2(x))
    @test norm(a * x) ≈ norm(a) * norm(x)
    @test norm(y + y2) <= norm(y) + norm(y2) || norm(y + y2) ≈ norm(y) + norm(y2)

    xs = rand(Form{D,R1,SVector{4,T}})
    xs2 = rand(Form{D,R1,SVector{4,T}})
    norm2(xs)::T
    @test norm2(xs) >= 0
    @test norm2(a * xs) == norm2(a) * norm2(xs)
    norm(xs)::float(T)
    @test norm(xs) ≈ sqrt(norm2(xs))
    @test norm(a * xs) ≈ norm(a) * norm(xs)
    @test norm(xs + xs2) <= norm(xs) + norm(xs2) || norm(xs + xs2) ≈ norm(xs) + norm(xs2)

    # cross product: x × y = ⋆(x ∧ y)
    Rcross = D - (R1 + R2)
    if 0 <= Rcross
        (x × y)::Form{D,Rcross}
        @test x × y === ⋆(x ∧ y)
    end
end

const Dmax4 = min(4, Dmax)
@testset "Tensor sums of forms D1=$D1 D2=$D2 R=$R" for D1 in 0:Dmax4, D2 in 0:Dmax4, R in 1:min(D1, D2)
    D = D1 + D2

    T = Rational{Int64}
    x = rand(Form{D1,R,T})
    x2 = rand(Form{D1,R,T})
    y = rand(Form{D2,R,T})
    y2 = rand(Form{D2,R,T})
    a = rand(T)

    # units
    if D1 == 1 && D2 == 1 && R == 1
        u(d, inds...) = unit(Form{d,length(inds),T}, inds...)

        @test u(1, 1) ⊕ u(1, 1) === u(2, 1) + u(2, 2)

        @test u(1, 1) ⊕ u(2, 1) === u(3, 1) + u(3, 2)
        @test u(1, 1) ⊕ u(2, 2) === u(3, 1) + u(3, 3)

        @test u(3, 1) ⊕ u(2, 1) === u(5, 1) + u(5, 4)
        @test u(3, 1) ⊕ u(2, 2) === u(5, 1) + u(5, 5)
        @test u(3, 2) ⊕ u(2, 1) === u(5, 2) + u(5, 4)
        @test u(3, 2) ⊕ u(2, 2) === u(5, 2) + u(5, 5)
        @test u(3, 3) ⊕ u(2, 1) === u(5, 3) + u(5, 4)
        @test u(3, 3) ⊕ u(2, 2) === u(5, 3) + u(5, 5)
    end

    @test ⊕(x) === x
    @test a * ⊕(x) === ⊕(a * x)

    (x ⊕ y)::Form{D,R,T}
    @test a * (x ⊕ y) === (a * x) ⊕ (a * y)
    @test (x ⊕ y) * a === (x * a) ⊕ (y * a)
    @test x ⊕ y === (zero(x) ⊕ y) + (x ⊕ zero(y))
    @test (x + x2) ⊕ (y + y2) === (x ⊕ y) + (x2 ⊕ y2)

    @test reverse_basis(x ⊕ y) ===
          # bitsign(R * R) *
          (reverse_basis(y) ⊕ reverse_basis(x))

    for D3 in R:Dmax4
        z = rand(Form{D3,R,T})
        @test (x ⊕ y) ⊕ z === x ⊕ (y ⊕ z)
        @test x ⊕ y ⊕ z === x ⊕ (y ⊕ z)
        @test ⊕(x, y, z) === x ⊕ (y ⊕ z)
    end
end

const Dmax3 = min(3, Dmax)
@testset "Tensor products of forms D1=$D1 R1=$R1 D2=$D2 R2=$R2" for D1 in 0:Dmax3, D2 in 0:Dmax3, R1 in 0:D1, R2 in 0:D2
    D = D1 + D2
    R = R1 + R2

    T = Rational{Int64}
    e = one(Form{0,0,T})
    x = rand(Form{D1,R1,T})
    x2 = rand(Form{D1,R1,T})
    y = rand(Form{D2,R2,T})
    y2 = rand(Form{D2,R2,T})
    a = rand(T)

    # units
    if D1 == 1 && D2 == 1 && R1 == 1 && R2 == 1
        u(d, inds...) = unit(Form{d,length(inds),T}, inds...)

        @test u(1) ⊗ u(1) === u(2)
        @test u(1, 1) ⊗ u(1) === u(2, 1)
        @test u(1) ⊗ u(1, 1) === u(2, 2)
        @test u(1, 1) ⊗ u(1, 1) === u(2, 1, 2)

        @test u(1) ⊗ u(2) === u(3)
        @test u(1, 1) ⊗ u(2) === u(3, 1)
        @test u(1) ⊗ u(2, 1) === u(3, 2)
        @test u(1, 1) ⊗ u(2, 1) === u(3, 1, 2)
        @test u(1) ⊗ u(2, 2) === u(3, 3)
        @test u(1, 1) ⊗ u(2, 2) === u(3, 1, 3)
        @test u(1) ⊗ u(2, 1, 2) === u(3, 2, 3)
        @test u(1, 1) ⊗ u(2, 1, 2) === u(3, 1, 2, 3)
    end

    @test ⊗(x) === x
    @test a * ⊗(x) === ⊗(a * x)

    (x ⊗ y)::Form{D,R,T}
    @test a * (x ⊗ y) === (a * x) ⊗ y
    @test (x ⊗ y) * a === x ⊗ (y * a)
    @test e ⊗ y === y
    @test x ⊗ e === x
    @test (x + x2) ⊗ y === x ⊗ y + x2 ⊗ y
    @test x ⊗ (y + y2) === x ⊗ y + x ⊗ y2

    @test reverse_basis(x ⊗ y) === bitsign(R1 * R2) * (reverse_basis(y) ⊗ reverse_basis(x))

    for D3 in 0:Dmax3, R3 in 0:D3
        z = rand(Form{D3,R3,T})
        @test (x ⊗ y) ⊗ z === x ⊗ (y ⊗ z)
        @test x ⊗ y ⊗ z === x ⊗ (y ⊗ z)
        @test ⊗(x, y, z) === x ⊗ (y ⊗ z)
    end
end


text:
const GEESE = ["African", "Roman Tufted", "Toulouse", "Pilgrim", "Steinbacher"]

goosefilter(birds) = filter( x -> x ∉ GEESE, birds )

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
"""
    five_qubit_code() -> SimpleCode

Return the five-qubit code.
"""
function five_qubit_code()
    stabilizers = [[1,3,3,1,0], [0,1,3,3,1], [1,0,1,3,3], [3,1,0,1,3]]
    logicals = [[1,1,1,1,1], [3,3,3,3,3]]
    pure_errors = [[0,1,0,0,0], [0,0,0,0,3], [0,0,3,0,0], [1,0,0,0,0]]
    return SimpleCode("Five qubit code", stabilizers, logicals, pure_errors)
end

"""
    five_qubit_surface_code() -> SimpleCode

Return the five-qubit surface code.
"""
function five_qubit_surface_code()
    stabilizers = [[1,1,1,0,0], [0,0,1,1,1], [3,0,3,3,0], [0,3,3,0,3]]
    logicals = [[1,0,0,1,0], [3,3,0,0,0]]
    pure_errors = [[3,0,0,0,0], [0,0,0,3,0], [1,0,0,0,0], [0,1,0,0,0]]
    return SimpleCode("Five qubit surface code", stabilizers, logicals, pure_errors)
end

"""
    steane_code() -> SimpleCode

Return the seven-qubit Steane code.
"""
function steane_code()
    stabilizers = [[1,0,0,1,0,1,1], [0,1,0,1,1,0,1], [0,0,1,0,1,1,1], [3,0,0,3,0,3,3],
        [0,3,0,3,3,0,3], [0,0,3,0,3,3,3]]
    logicals = [[1,1,1,1,1,1,1], [3,3,3,3,3,3,3]]
    pure_errors = [[3,0,0,0,0,0,0], [0,3,0,0,0,0,0], [0,0,3,0,0,0,0], [1,0,0,0,0,0,0],
        [0,1,0,0,0,0,0], [0,0,1,0,0,0,0]]
    return SimpleCode("Steane code", stabilizers, logicals, pure_errors)
end

"""
    random_code(n::Int, k::Int) -> SimpleCode

Return a random code with `n` physical qubits and `k` logicals qubits.

An `ErrorException` is raised if `k` > `n`.
"""
function random_code(n::Int, k::Int)
    0 <= k <= n || error("invalid number of logical qubits!")

    r = n - k

    stabilizer_state = random_stabilizer_state(n)
    stabilizers = stabilizer_state.stabilizers[1:r]
    pure_errors = stabilizer_state.pure_errors[1:r]

    logicals = Vector{Int}[]
    for α in r + 1:n
        push!(logicals, stabilizer_state.stabilizers[α])
        push!(logicals, stabilizer_state.pure_errors[α])
    end

    return SimpleCode("[$n,$k] random code", stabilizers, logicals, pure_errors)
end

"""
    random_stabilizer_state(n::Int) -> SimpleCode

Return a random code with `n` physical qubits and no logicals.
"""
function random_stabilizer_state(n::Int)
    stabilizers = Array{Int64,1}[]
    for α in 1:n
        @label try_again
        new_stabilizer = rand(0:3, n)

        if pauli_weight(new_stabilizer) <= 1
            @goto try_again
        elseif !pauli_are_independent(vcat(stabilizers, [new_stabilizer]))
            @goto try_again
        elseif !pauli_are_commuting(vcat(stabilizers, [new_stabilizer]))
            @goto try_again
        end

        push!(stabilizers, new_stabilizer)
    end

    logicals = []

    return SimpleCode("$n-qubit random stabilizer state", stabilizers, logicals)
end


text:
module ProbeGames

include("../games/probe-games/a1r1/main.jl")
export A1R1

include("../games/probe-games/a1o2r2/main.jl")
export A1O2R2

include("../games/probe-games/a1r1s2/main.jl")
export A1R1S2

include("../games/probe-games/a2r2/main.jl")
export A2R2

include("../games/probe-games/two_player/main.jl")
export TWOPLAYER

include("../games/probe-games/simpler_tictactoe/main.jl")
export SimplerTicTacToe

games = Dict(
  "a1r1" => A1R1.GameSpec(),
  "a1o2r2" => A1O2R2.GameSpec(),
  "a1r1s2" => A1R1S2.GameSpec(),
  "a2r2" => A2R2.GameSpec(),
  "twoplayer" => TWOPLAYER.GameSpec(),
  "simplertictactoe" => SimplerTicTacToe.GameSpec()
)
end #ProbeGames


text:
using TimerOutputs
using ImageFiltering: Fill, KernelFactors.gaussian, centered, kernelfactors, imfilter!, padarray, Pad
using LinearAlgebra: qr
using PyPlot: Figure

function goo()
    to = TimerOutput()
    @timeit_debug to "sleepytime" begin
        println(getfield(lap, :timeit_debug_enabled)())
        println("hello")
        sleep(1)
    end
    print_timer(to)
    return nothing
end

"""
    single_lap(target::Image, source::Image, filter_half_size::Integer, window, filter_count::Integer=3)

Return an estimate of a smoothly varying flow of size of `target` which is the displacement that transforms `source` closer to `target`.

# Arguments:
- `target::Image`: target/fixed grayscale image.
- `source::Image`: source/moving grayscale image.
- `filter_half_size`: the half size of the base of the gaussian filters used.
- `window`: the size of the local window (list of 2 ints) usually same as filter_size.
- `filter_count::Integer=3`: the number of basis filters used (so far only =3 implemented).

# Keyword Arguments:
- `timer::TimerOutput=TimerOutput("LAP")`: provide a timer which times certain blocks in the function.
- `display::Bool=false`: verbose and debug prints.

See also: [`pflap`](@ref), [`showflow`](@ref), [`warp_imgshowflow`](@ref), [`imgshowflow`](@ref)
"""
function single_lap(target::Image,
                    source::Image,
                    filter_half_size::Integer,
                    window,
                    filter_count::Integer=3;
                    timer::TimerOutput=TimerOutput("LAP"),
                    display::Bool=false)

    pixel_count = length(target)
    image_size = size(target)
    filter_size = 2 * filter_half_size + 1

    # prepare 2D filter basis for later use:
    basis = similar(target, (filter_size, filter_size, filter_count))

    # dif = (image forward - image backward) of each filter
    dif = similar(target, (image_size..., filter_count))

    @timeit_debug timer "filtering" begin
        if filter_count == 3
            # get the relevant forward and backward filters
            forward_ker, backward_ker = prepare_gaussian_filters(filter_half_size)

            # temporary place to store filtered images
            tmp_filtered_1 = similar(target)
            tmp_filtered_2 = similar(source)

            # basis 1 - (gaus, gaus) both forward and backward
            basis[:, :, 1] = broadcast(*, forward_ker[1]...)
            imfilter!(tmp_filtered_1, target, forward_ker[1], "symmetric")
            imfilter!(tmp_filtered_2, source, backward_ker[1], "symmetric")
            dif[:, :, 1] = tmp_filtered_2 - tmp_filtered_1

            # basis 2 - (gaus, gaus_der_flip) as forward, (gaus, gaus_der) as backward
            basis[:, :, 2] = broadcast(*, forward_ker[2]...)
            imfilter!(tmp_filtered_1, target, forward_ker[2], "symmetric")
            imfilter!(tmp_filtered_2, source, backward_ker[2], "symmetric")
            dif[:, :, 2] = tmp_filtered_2 - tmp_filtered_1

            # basis 3 - (gaus_der_flip, gaus) as forward, (gaus_der, gaus) as backward
            basis[:, :, 3] = broadcast(*, forward_ker[3]...)
            imfilter!(tmp_filtered_1, target, forward_ker[3], "symmetric")
            imfilter!(tmp_filtered_2, source, backward_ker[3], "symmetric")
            dif[:, :, 3] = tmp_filtered_2 - tmp_filtered_1
        end
    end

    dif = reshape(dif, (:, filter_count))

    # prepare matrices for linear system of equations
    A = similar(target, (filter_count-1, filter_count-1, pixel_count))
    b = similar(target, (filter_count-1, pixel_count))

    @timeit_debug timer "prepare A and b" begin
        w = Int64.((window[1]-1)/2)
        for k in 1:filter_count-1
            for l in k:filter_count-1
                @timeit_debug timer "window sum part 1" @views window_sum3!(A[k, l, :], dif[:, k+1] .* dif[:, l+1], image_size, window)
                A[l, k, :] = A[k, l, :]
            end
            @timeit_debug timer "window sum part 2" @views window_sum3!(b[k, :], dif[:, k+1] .* dif[:, 1] .* (-1), image_size, window)
        end
    end

    # Perform Gauss elimination on all pixels in parallel:
    # coeffs will be of shape: pixel_count, filter_count-1
    @timeit_debug timer "solve linear systems" begin
        @views coeffs = multi_mat_div_gem(A, b)
    end
    # adding ones so that all base filters have their coefficients even the first one
    all_coeffs = [ones(pixel_count) coeffs]

    # make a border mask. Is 0 if its in a border of filter_half_size size.
    window_half_size = Int64.((window .- 1) ./ 2)
    border_mask = parent(padarray(ones((image_size .- (2 .* window_half_size))...),
                    Fill(NaN, window_half_size, window_half_size)))
    @views all_coeffs = all_coeffs .* reshape(border_mask, (:, 1))

    k = (-filter_half_size:filter_half_size)

    # Get the displacement vector field from the filters
    @timeit_debug timer "calculate flow" begin
        u_est = similar(target, ComplexF64, image_size...)
        @views u_est[:] = -2 .* ((im .* (.-1 .* all_coeffs[:, 3]) ./ all_coeffs[:, 1]) .+ ((.-1 .* all_coeffs[:, 2]) ./ all_coeffs[:, 1]));
    end

    # dont use estimations whose displacement is larger than the filter_half_size
    displacement = real(u_est).^2 .+ imag(u_est).^2
    displacement_mask = displacement .> filter_half_size^2
    u_est[displacement_mask] .= NaN .+ NaN .* 1im;

    if (display && getfield(LAP_julia, :timeit_debug_enabled)())
        print_timer(timer)
        println()
    end
    return u_est
end

"""
    single_lap_at_points(args; kwargs)

Return a vector of displacements that transforms `source` closer to `target` at points `inds`.

# Arguments:
- `target::Image`: target/fixed grayscale image.
- `source::Image`: source/moving grayscale image.
- `filter_half_size::Integer`: the half size of the base of the gaussian filters used.
- `window`: the size of the local window (list of 2 ints) usually same as filter_size.
- `inds::Array{CartesianIndex{2},1}`: where to calculate the flow.

# Keyword Arguments:
- `filter_count::Integer=3`: the number of basis filters used (so far only =3 implemented).
- `timer::TimerOutput=TimerOutput("Sparse LAP")`: provide a timer which times certain blocks in the function.
- `display::Bool=false`: verbose and debug prints.

See also: [`pflap`](@ref), [`showflow`](@ref), [`single_lap`](@ref), [`sparse_pflap`](@ref).
"""
function single_lap_at_points(target::Image,
                              source::Image,
                              filter_half_size::Integer,
                              window,
                              inds::Array{CartesianIndex{2},1};
                              filter_count::Integer=3,
                              timer::TimerOutput=TimerOutput("Sparse LAP"),
                              display::Bool=false)

    image_size = size(target)
    pixel_count = length(target)
    ind_count = length(inds)
    filter_size = 2 * filter_half_size + 1

    # Prepare filters:
    # Calculate separable filters from basis:
    sigma = (filter_half_size + 2) / 4
    centered_inds = centered(-filter_half_size:filter_half_size)
    gaus = gaussian(sigma, filter_size)
    gaus_der = gaus .* centered_inds .* (-1)/sigma^2
    gaus_der_flip = reverse(gaus_der, dims=1)

    # prepare 2D filter basis for later use
    basis = similar(target, (filter_size, filter_size, filter_count))

    # dif = (points of forward image - points of backward image) for each filter
    dif = similar(target, (image_size..., filter_count))

    @timeit_debug timer "filtering" begin
        if filter_count == 3
            # temporary place to store filtered images
            # tmp_filtered_1 = fill(NaN, image_size)
            # tmp_filtered_2 = fill(NaN, image_size)
            tmp_filtered_1 = similar(target)
            tmp_filtered_2 = similar(target)

            # basis 1 - (gaus, gaus) both forward and backward
            kernf_1 = kernelfactors((gaus, gaus))
            basis[:, :, 1] = broadcast(*, kernf_1...)
            # filt_onebyone!(tmp_filtered_1, target, kernf_1, filter_half_size, points)
            # filt_onebyone!(tmp_filtered_2, source, kernf_1, filter_half_size, points)
            imfilter!(tmp_filtered_1, target, kernf_1, "symmetric")
            imfilter!(tmp_filtered_2, source, kernf_1, "symmetric")
            dif[:, :, 1] = tmp_filtered_2 - tmp_filtered_1

            # basis 2 - (gaus, gaus_der_flip) as forward, (gaus, gaus_der) as backward
            kernf_2f = kernelfactors((gaus, gaus_der_flip))
            basis[:, :, 2] = broadcast(*, kernf_2f...)
            # filt_onebyone!(tmp_filtered_1, target, kernf_2f, filter_half_size, points)
            imfilter!(tmp_filtered_1, target, kernf_2f, "symmetric")
            kernf_2b = kernelfactors((gaus, gaus_der))
            # filt_onebyone!(tmp_filtered_2, source, kernf_2b, filter_half_size, points)
            imfilter!(tmp_filtered_2, source, kernf_2b, "symmetric")
            dif[:, :, 2] = tmp_filtered_2 - tmp_filtered_1

            # basis 3 - (gaus_der_flip, gaus) as forward, (gaus_der, gaus) as backward
            kernf_3f = kernelfactors((gaus_der_flip, gaus))
            basis[:, :, 3] = broadcast(*, kernf_3f...)
            # filt_onebyone!(tmp_filtered_1, target, kernf_3f, filter_half_size, points)
            imfilter!(tmp_filtered_1, target, kernf_3f, "symmetric")
            kernf_3b = kernelfactors((gaus_der, gaus))
            # filt_onebyone!(tmp_filtered_2, source, kernf_3b, filter_half_size, points)
            imfilter!(tmp_filtered_2, source, kernf_3b, "symmetric")
            dif[:, :, 3] = tmp_filtered_2 - tmp_filtered_1
        end
    end # "filtering"

    dif = reshape(dif, (:, filter_count))

    # prepare matrices for linear system of equations
    A = similar(target, (filter_count-1, filter_count-1, ind_count))
    b = similar(target, (filter_count-1, ind_count))

    @timeit_debug timer "prepare A and b" begin
        for k in 1:filter_count-1
            for l in k:filter_count-1
                @timeit_debug timer "window sum part 1" begin
                    @views window_sum3_at_inds!(A[k, l, :], dif[:, k+1] .* dif[:, l+1], image_size, window, inds)
                    A[l, k, :] = A[k, l, :]
                end
            end
            @timeit_debug timer "window sum part 2" begin
                @views window_sum3_at_inds!(b[k, :], dif[:, k+1] .* dif[:, 1] .* (-1), image_size, window, inds)
            end
        end
    end # "prepare A and b"

    # Perform Gauss elimination on all pixels in parallel:
    # coeffs will be of shape: pixel_count, filter_count-1
    @timeit_debug timer "solve linear systems" begin
        # @views coeffs = multi_mat_div_qr(A, b)
        @views coeffs = multi_mat_div_gem(A, b)
        # @views coeffs = multi_mat_div_qr(A, b)
    end

    # adding ones so that all base filters have their coefficients even the first one
    all_coeffs = [ones(ind_count) coeffs]

    k = (-filter_half_size:filter_half_size)

    # Get the displacement vector field from the filters
    @timeit_debug timer "calculate flow" begin
        u_est_at_inds = similar(target, ComplexF64, ind_count)
        @views u_est_at_inds[:] = -2 .* ((im .* (.-1 .* all_coeffs[:, 3]) ./ all_coeffs[:, 1]) .+ ((.-1 .* all_coeffs[:, 2]) ./ all_coeffs[:, 1]));
    end # "calculate flow"

    # dont use estimations whose displacement is larger than the filter_half_size
    displacement_mask = (real(u_est_at_inds).^2 .+ imag(u_est_at_inds).^2) .<= filter_half_size^2

    return u_est_at_inds[displacement_mask], inds[displacement_mask]
end

# RESIDUAL CALCULATION AND EXPERIMENTS:

# addpoints(inds[(real(u_est_at_inds).^2 .+ imag(u_est_at_inds).^2) .> filter_half_size^2])
# addpoints(inds[.!((real(u_est_at_inds).^2 .+ imag(u_est_at_inds).^2) .> filter_half_size^2)])
#
# err = sum(all_coeffs .* dif[map(x -> to_lin_index(x, image_size), inds), :], dims=2)
#
# dif = reshape(dif, (image_size..., filter_count))
#
# window = (5,5)
# whs = round.(Int64, (window.-(1,1))./2)
# shift = CartesianIndex(whs)
#
# @timeit_debug timer "calc residual" begin
#     for (k, ind) in enumerate(inds)
#         err[k] = sum(dif[:,:, 1][ind-shift:ind+shift].^2 .+
#                     (dif[:,:, 2][ind-shift:ind+shift] .* all_coeffs[k, 2]) .^2 .+
#                     (dif[:,:, 3][ind-shift:ind+shift] .* all_coeffs[k, 3]) .^2)
#     end
# end

function single_lap_at_points_float16(target::Image,
                              source::Image,
                              filter_half_size::Integer,
                              window,
                              inds::Array{CartesianIndex,1};
                              filter_count::Integer=3,
                              timer::TimerOutput=TimerOutput("Sparse LAP"),
                              display::Bool=false)

    image_size = size(target)
    pixel_count = length(target)
    ind_count = length(inds)
    filter_size = 2 * filter_half_size + 1

    # Prepare filters:
    # Calculate separable filters from basis:
    sigma::Float16 = (filter_half_size + 2) / 4
    centered_inds = centered(-filter_half_size:filter_half_size)
    gaus = gaussian(sigma, filter_size)
    gaus_der = gaus .* centered_inds .* (-1)/sigma^2
    gaus_der_flip = reverse(gaus_der, dims=1)

    # prepare 2D filter basis for later use
    basis = similar(target, (filter_size, filter_size, filter_count))

    # dif = (points of forward image - points of backward image) for each filter
    dif = similar(target, (image_size..., filter_count))

    @timeit_debug timer "filtering" begin
        if filter_count == 3
            # temporary place to store filtered images
            # tmp_filtered_1 = fill(NaN, image_size)
            # tmp_filtered_2 = fill(NaN, image_size)
            tmp_filtered_1 = similar(target)
            println(eltype(target))
            println(eltype(tmp_filtered_1))
            @info eltype(gaus)
            tmp_filtered_2 = similar(target)

            # basis 1 - (gaus, gaus) both forward and backward
            kernf_1 = kernelfactors((gaus, gaus))
            basis[:, :, 1] = broadcast(*, kernf_1...)
            # filt_onebyone!(tmp_filtered_1, target, kernf_1, filter_half_size, points)
            # filt_onebyone!(tmp_filtered_2, source, kernf_1, filter_half_size, points)
            imfilter!(tmp_filtered_1, target, kernf_1, "symmetric")
            imfilter!(tmp_filtered_2, source, kernf_1, "symmetric")
            dif[:, :, 1] = tmp_filtered_2 - tmp_filtered_1

            # basis 2 - (gaus, gaus_der_flip) as forward, (gaus, gaus_der) as backward
            kernf_2f = kernelfactors((gaus, gaus_der_flip))
            basis[:, :, 2] = broadcast(*, kernf_2f...)
            # filt_onebyone!(tmp_filtered_1, target, kernf_2f, filter_half_size, points)
            imfilter!(tmp_filtered_1, target, kernf_2f, "symmetric")
            kernf_2b = kernelfactors((gaus, gaus_der))
            # filt_onebyone!(tmp_filtered_2, source, kernf_2b, filter_half_size, points)
            imfilter!(tmp_filtered_2, source, kernf_2b, "symmetric")
            dif[:, :, 2] = tmp_filtered_2 - tmp_filtered_1

            # basis 3 - (gaus_der_flip, gaus) as forward, (gaus_der, gaus) as backward
            kernf_3f = kernelfactors((gaus_der_flip, gaus))
            basis[:, :, 3] = broadcast(*, kernf_3f...)
            # filt_onebyone!(tmp_filtered_1, target, kernf_3f, filter_half_size, points)
            imfilter!(tmp_filtered_1, target, kernf_3f, "symmetric")
            kernf_3b = kernelfactors((gaus_der, gaus))
            # filt_onebyone!(tmp_filtered_2, source, kernf_3b, filter_half_size, points)
            imfilter!(tmp_filtered_2, source, kernf_3b, "symmetric")
            dif[:, :, 3] = tmp_filtered_2 - tmp_filtered_1
        end
    end # "filtering"

    dif = reshape(dif, (:, filter_count))

    # prepare matrices for linear system of equations
    A = similar(target, (filter_count-1, filter_count-1, ind_count))
    b = similar(target, (filter_count-1, ind_count))

    @timeit_debug timer "prepare A and b" begin
        for k in 1:filter_count-1
            for l in k:filter_count-1
                @timeit_debug timer "window sum part 1" begin
                    @views window_sum3_at_inds!(A[k, l, :], dif[:, k+1] .* dif[:, l+1], image_size, window, inds)
                    A[l, k, :] = A[k, l, :]
                end
            end
            @timeit_debug timer "window sum part 2" begin
                @views window_sum3_at_inds!(b[k, :], dif[:, k+1] .* dif[:, 1] .* (-1), image_size, window, inds)
            end
        end
    end # "prepare A and b"

    # Perform Gauss elimination on all pixels in parallel:
    # coeffs will be of shape: pixel_count, filter_count-1
    @timeit_debug timer "multi mat div" begin
        # @views coeffs = multi_mat_div_qr(A, b)
        @views coeffs = multi_mat_div_gem(A, b)
    end

    # adding ones so that all base filters have their coefficients even the first one
    all_coeffs = [ones(ind_count) coeffs]

    k = (-filter_half_size:filter_half_size)

    # Get the displacement vector field from the filters
    @timeit_debug timer "calculate flow" begin
        u_est_at_inds = similar(target, ComplexF16, ind_count)
        @views u_est_at_inds[:] = -2 .* ((im .* (.-1 .* all_coeffs[:, 3]) ./ all_coeffs[:, 1]) .+ ((.-1 .* all_coeffs[:, 2]) ./ all_coeffs[:, 1]));
    end # "calculate flow"

    # dont use estimations whose displacement is larger than the filter_half_size
    displacement_mask = (real(u_est_at_inds).^2 .+ imag(u_est_at_inds).^2) .<= filter_half_size^2

    return u_est_at_inds[displacement_mask], inds[displacement_mask]
end


text:
function matrix_cone(S::Type{<:Union{MatrixConeTriangle,
                                     MOI.PositiveSemidefiniteConeTriangle}},
                     side_dimension)
    if iszero(side_dimension)
        # Some solvers such as Mosek does not support 0-dimensional PSD cone
        return EmptyCone()
    elseif isone(side_dimension)
        # PSD constraints on 1x1 matrices are equivalent to the nonnegativity
        # of the only entry.
        return MOI.Nonnegatives(1)
    elseif side_dimension == 2
        # PSD constraints on 2x2 matrices are SOC representable.
        return PositiveSemidefinite2x2ConeTriangle()
    else
        # PSD constraints on nxn matrices with n ≥ 3 is not SOC representable,
        # see [F18].
        #
        # [F18] Fawzi, Hamza
        # On representing the positive semidefinite cone using the second-order
        # cone.
        # Mathematical Programming (2018): 1-10.
        return S(side_dimension)
    end
end

"""
    struct DiagonallyDominantConeTriangle <: MatrixConeTriangle
        side_dimension::Int
    end

See Definition 4 of [AM17] for a precise definition of the last two items.

[AM17] Ahmadi, A. A. & Majumdar, A.
*DSOS and SDSOS Optimization: More Tractable Alternatives to Sum of Squares and Semidefinite Optimization*
ArXiv e-prints, **2017**.
"""
struct DiagonallyDominantConeTriangle <: MatrixConeTriangle
    side_dimension::Int
end

function matrix_cone(S::Type{DiagonallyDominantConeTriangle},
                     side_dimension)
    if iszero(side_dimension)
        return EmptyCone()
    elseif isone(side_dimension)
        return MOI.Nonnegatives(1)
    else
        # With `side_dimension` = 2, we want to avoid using SOC and only use LP
        return S(side_dimension)
    end
end

"""
    struct ScaledDiagonallyDominantConeTriangle <: MatrixConeTriangle
        side_dimension::Int
    end

See Definition 4 of [AM17] for a precise definition of the last two items.

[AM17] Ahmadi, A. A. & Majumdar, A.
*DSOS and SDSOS Optimization: More Tractable Alternatives to Sum of Squares and Semidefinite Optimization*
ArXiv e-prints, **2017**.
"""
struct ScaledDiagonallyDominantConeTriangle <: MatrixConeTriangle
    side_dimension::Int
end

function matrix_cone(S::Type{ScaledDiagonallyDominantConeTriangle},
                     side_dimension)
    if iszero(side_dimension)
        return EmptyCone()
    elseif isone(side_dimension)
        return MOI.Nonnegatives(side_dimension)
    elseif side_dimension == 2
        return PositiveSemidefinite2x2ConeTriangle()
    else
        return S(side_dimension)
    end
end

function side_dimension(set::Union{MatrixConeTriangle,
                                   MOI.PositiveSemidefiniteConeTriangle})
    return set.side_dimension
end

# isbits types, nothing to copy
function Base.copy(set::MatrixConeTriangle)
    return set
end

# TODO make PSDConeTriangle inherit from MatrixConeTriangle to remove the need
#      for this
function MOI.dimension(set::MatrixConeTriangle)
    return div(side_dimension(set) * (side_dimension(set) + 1), 2)
end

function MOIU.set_dot(x::Vector, y::Vector, set::MatrixConeTriangle)
    return MOIU.triangle_dot(x, y, side_dimension(set), 0)
end

function MOIU.dot_coefficients(a::Vector, set::MatrixConeTriangle)
    b = copy(a)
    MOIU.triangle_coefficients!(b, side_dimension(set), 0)
    return b
end


text:
#=

  Assignment 1:
  https://edu.swi-prolog.org/mod/assign/view.php?id=242
  """
  http://cplint.eu/p/urns_and_balls.swinb

  Urns and balls

  Suppose you have two urns: urn1 contains 40 blue balls and 20 red balls and urn2 contains 25
  blue balls and 30 red balls.

  You throw an unbiased coin and, if it turns out head, you draw a ball from the first urn,
  it it turns out tails you draw a ball from the second urn.

  Write a program modeling this process and a query for answering the question
  "What is the probability of drawing a blue ball?"

  Write the program and the query in the cells below. Then save the notebook and submit the url.

  """
  The exact answer of drawing a blue call is
     0.5*40/60 + 0.5*25/55 = 0.56060606060606060606

  Summary Statistics
  parameters      mean       std   naive_se      mcse          ess      rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64      Float64   Float64       Float64

        coin    1.5019    0.5000     0.0025    0.0025   40549.0446    1.0000    28840.0032
        draw    1.4389    0.4963     0.0025    0.0023   40570.3731    1.0000    28855.1729

  Distributions of variable draw
  blue       =>   22442  (0.561050)
  red        =>   17558  (0.438950)

  Distributions of variable coin
  head       =>   20078  (0.501950)
  tail       =>   19922  (0.498050)

  Cf ~/cplint/course_urns_and_balls.pl
     ~/blog/urns_and_balls.blog
     ~/psi/urns_and_balls.blog
     ~/webppl/urns_and_balls.wppl

=#

using Turing, StatsPlots, DataFrames
using ReverseDiff, Zygote, Tracker
# Turing.setadbackend(:reversediff)
# Turing.setadbackend(:zygote)
# Turing.setadbackend(:tracker)
include("jl_utils.jl")

@model function urns_and_balls()
    tail = 1
    head = 2
    coin ~ Categorical(simplex([0.5,0.5]))

    blue = 1
    red = 2
    draw ~ coin == head ? Categorical(simplex([40,20])) : Categorical(simplex([25,30]))
end

model = urns_and_balls()

num_chns = 4
# chns = sample(model, Prior(), MCMCThreads(), 1000, num_chns)

# chns = sample(model, MH(), 40_000)
# chns = sample(model, MH(), MCMCThreads(), 40_000, num_chns)

# chns = sample(model, PG(20), MCMCThreads(), 10_000, num_chns)
# chns = sample(model, PG(20), 40_000)
# chns = sample(model, IS(), MCMCThreads(), 10_000, num_chns)
chns = sample(model, IS(), 40_000)

display(chns)

show_var_dist_pct(chns,:draw,["blue","red"])
show_var_dist_pct(chns,:coin,["tail","head"])

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# test the performance of preconditioner

using TensorKit
using TensorOperations
using TensorKitAD
using KrylovKit
using Zygote
using KrylovKit
using LinearAlgebra
using Plots
using Optim
using ChainRules
using ChainRulesCore
using Printf

using Revise
using statmech_tm_solver

function Base.show(io::IO, t::OptimizationState)
    @printf io "%6d   %14.14e   %14e\n" t.iteration t.value t.g_norm
end
#########################################################

# parameters
c = 1
μ = 1
L = 50
d = 1
nsteps = 200
χ = 12

# cost function (energy) and its gradient
function f(ψ_arr::Array{ComplexF64, 3})
    ψ = convert_to_cmps(ψ_arr)
    return energy_lieb_liniger(ψ, c, L, μ)
end

function g!(ψg_arr::Array{ComplexF64, 3}, ψ_arr::Array{ComplexF64, 3})
    copyto!(ψg_arr, f'(ψ_arr))
end

#ψ0 = cmps(rand, χ, 1)
ψ0 = quickload("lieb_liniger_c$(c)_mu$(μ)_L$(L)_chi$(χ-4)") |> convert_to_cmps
ψ0 = expand(ψ0, χ, 0.001)
ψ_arr0 = convert_to_array(ψ0)
E0 = energy_lieb_liniger(ψ0, c, L, μ)

res = optimize(f, g!, ψ_arr0, LBFGS(), Optim.Options(show_trace=true, iterations=nsteps))
ψ_arr = Optim.minimizer(res);
ψm = convert_to_cmps(ψ_arr);
ψm_grad = convert_to_tensormap(f'(ψ_arr), 2)
E1 = energy_lieb_liniger(ψm, c, L, μ)

# optimize with preconditioners
function _precondprep!(P::preconditioner, ψ_arr::Array{ComplexF64, 3})
    ψ = convert_to_cmps(ψ_arr)
    ψgrad = convert_to_tensormap(f'(ψ_arr), 2)

    P1 = preconditioner(ψ, ψgrad, L; gauge=:periodic)
    P.map = P1.map
    P.invmap = P1.invmap
    P.proj = P1.proj
end

P0 = preconditioner(ψm, ψm_grad, L);
res_w_precond = optimize(f, g!, ψ_arr, LBFGS(P = P0, precondprep=_precondprep!), Optim.Options(show_trace=true, store_trace=true, iterations=nsteps))
ψm_w_precond = convert_to_cmps(Optim.minimizer(res_w_precond));

trace_w_precond = Optim.trace(res_w_precond)
times_w_precond = [trace_w_precond[ix].metadata["time"] for ix in 1:nsteps+1]
values_w_precond = [trace_w_precond[ix].value for ix in 1:nsteps+1]
gnorms_w_precond = [trace_w_precond[ix].g_norm for ix in 1:nsteps+1]

# optimize without preconditioners
res_wo_precond = optimize(f, g!, ψ_arr, LBFGS(), Optim.Options(show_trace=true, store_trace=true, iterations=2*nsteps))
ψm_wo_precond = convert_to_cmps(Optim.minimizer(res_wo_precond));

trace_wo_precond = Optim.trace(res_wo_precond)
times_wo_precond = [trace_wo_precond[ix].metadata["time"] for ix in 1:2*nsteps+1]
values_wo_precond = [trace_wo_precond[ix].value for ix in 1:2*nsteps+1]
gnorms_wo_precond = [trace_wo_precond[ix].g_norm for ix in 1:2*nsteps+1]

plot(times_w_precond, log10.(gnorms_w_precond))
plot!(times_wo_precond, log10.(gnorms_wo_precond))

plot(times_w_precond, values_w_precond)
plot!(times_w_precond, values_w_precond); @show "ploted"
plot!(times_wo_precond, values_wo_precond)


text:
# Copyright (c) 2021, Rolfe Power

"""
    DimensionSet

Struct representing fundamental quantities in an absolute dimension system.
"""
struct DimensionSet
    mass::Float64
    length::Float64
    time::Float64
end

"""
    mass(::DimensionSet)

Get fundamental mass quantity
"""
mass(ds::DimensionSet) = ds.mass

"""
    length(::DimensionSet)

Get fundamental length quantity
"""
Base.length(ds::DimensionSet) = ds.length

"""
    time(::DimensionSet)

Get fundamental time quantity
"""
Base.time(ds::DimensionSet) = ds.time

"""
    velocity(ds)

Return the velocity found as the ratio of length and time.
"""
velocity(ds) = length(ds) / time(ds)

"""
    acceleration(ds)

Return the acceleration found as the ratio of velocity and time.
"""
acceleration(ds) = velocity(ds) / time(ds)

"""
    area(ds)

Return the velocity found as the square of length.
"""
area(ds) = length(ds)^2

"""
    volume(ds)

Return the velocity found as the cube of length.
"""
volume(ds) = area(ds) * length(ds)

"""
    force(ds)

Return the force found as the product of acceleration and mass
"""
force(ds) = acceleration(ds) * mass(ds)

"""
    density(ds)

Return the acceleration found as the ratio of mass and volume.
"""
density(ds) = mass(ds) / volume(ds)

"""
    DimensionSet(value; mass=value, length=value, time=value)

Construct a [`DimensionSet`](@ref) specifying (optionally) `mass`, `length`, and `time`.

!!! note

    In addition to constructing the [`DimensionSet`](@ref), this function also validates
    the input arguments. Input dimensional quantities must be positive values!

## Examples
```jldoctest
julia> DimensionSet(2)
Pat.DimensionSet(2.0, 2.0, 2.0)

julia> DimensionSet(mass=2.0, length=3.0, time=4.0)
Pat.DimensionSet(2.0, 3.0, 4.0)

julia> DimensionSet(12.0, time=0.5)
Pat.DimensionSet(12.0, 12.0, 0.5)

julia> DimensionSet(-1.0)
ERROR: DomainError with -1.0:
Non-positive mass
[...]
```
"""
function DimensionSet(value=1.0; mass=value, length=value, time=value)
    mass   > 0 || throw(DomainError(mass, "Non-positive mass"))
    length > 0 || throw(DomainError(length, "Non-positive length"))
    time   > 0 || throw(DomainError(time, "Non-positive time"))
    DimensionSet(
        mass,
        length,
        time
    )
end


text:
# Below, the symbols of GeometryPrimitives such as KDTree, bounds, surfpt_nearby, normal do
# not have to be re-defined here to be exported.  When GeometryPrimitives is imported to
# MaxwellBase by "using GeometryPrimitives", the symbols of GeometryPrimitives are already
# inside the namespace of MaxwellBase, so we simply need to export them to use outside
# MaxwellBase without qualification.
#
# One caveat, though.  In this file, bounds is defined for Interval, so if we simply export
# bounds, then only this is exported.  So, when defining bounds for Interval, we have to
# "extend" bounds for Shape by defining it as GeometryPrimitives.bounds(::Interval) = ....
# Then, exporting bounds exports this whole collection of bounds, both for Shape and Interval.
export Object
export shape, matparam, max∆l, sub_pind2matprm, add_obj!, periodize  #, surfpt_nearby, normal
# export lsf, bound_, L_, center_, dist2bound, bound_contains, ∆lmax, sphere, transform,
#     surfnormal, surfpoint  # functions
# import Base:size, getindex, contains, isless, union, intersect

mutable struct Object{K,S<:Shape{K}}  # K: dimension of space
    shp::S
    mat::Material
    ∆lmax::SVec{K,Float}
    Object{K,S}(shp, mat, ∆lmax) where {K,S} = new(shp, mat, ∆lmax)
end

Object(shp::S, mat::Material, ∆lmax::SVec{K}=@SVector(fill(Inf,K))) where {K,S<:Shape{K}} =
    Object{K,S}(shp, mat, ∆lmax)
Object(shp::Shape{K}, mat::Material, ∆lmax::AbsVec) where {K} = Object(shp, mat, SVec{K}(∆lmax))
Object(shp::Shape{K}, mat::Material, ∆lmax::Real) where {K} = Object(shp, mat, SVec(ntuple(k->∆lmax, Val(K))))

# Create the user interface such that maxwellsys.add(shape, material, ∆lmax) uses this function.
# setmat!(obj::Object, m::Material) = (obj.mat = m; obj)
# setmax∆l!(obj::Object{K}, ∆lmax::AbsVec) where {K} = (obj.∆lmax = SVec{K}(∆lmax); obj)
# setmax∆l!(obj::Object{K}, ∆lmax::Number) where {K} = setmax∆l!(obj, SVec(ntuple(k->∆lmax, Val(K))))

shape(obj::Object) = obj.shp
matparam(obj::Object, ft::FieldType) = matparam(obj.mat, ft)
max∆l(obj::Object) = obj.∆lmax

# Relationships between pind2matprm, oind2pind, oind2shp
#
# - pind2matprm is a vector of material parameter tensors.  pind2matprm[n] is a 3×3 material
# parameter tensor for the parameter index n.
#
# - oind2pind is a map from the object index to the parameter index.
#
# - oind2shp is a map from the object index to the shape.

# Returns a vector of submatrices of material parameter tensors.
# We return a vector of submatrices rather than a vector of submatrix views, because
# assign_param!() requires the entries of its argument pind2matprm to be SMatrix rather than
# views of SMatrix.
sub_pind2matprm(pind2matprm::AbsVec{S²ComplexF{Kf,Kf²}}, inds::SInteger) where {Kf,Kf²} =
    [mp[inds,inds] for mp = pind2matprm]  # mp[inds,inds]: S²Mat

# Consider using resize! on oind2obj.
function add_obj!(oind2shp::AbsVec{Shape{K,K²}}, oind2pind::Tuple2{AbsVec{ParamInd}},
                  pind2matprm::Tuple{AbsVec{S²ComplexF3},AbsVec{S²ComplexF3}},
                  objs::AbsVec{Object{K,<:Shape{K,K²}}}
                  ) where {K,K²}
    for obj = objs
        add_obj!(oind2obj, oind2pind, pind2matprm, obj)
    end
end

# Consider using resize! on oind2obj.
function add_obj!(oind2shp::AbsVec{Shape{K,K²}}, oind2pind::Tuple2{AbsVec{ParamInd}},
                  pind2matprm::Tuple{AbsVec{S²ComplexF3},AbsVec{S²ComplexF3}},
                  objs::Object{K,<:Shape{K,K²}}...
                  ) where {K,K²}
    for obj = objs
        add_obj!(oind2shp, oind2pind, pind2matprm, obj)
    end
end

# Populate `oind2shp`, `oind2pind`, `pind2matprm` with a given object.  See above for the
# meaning of `oind2shp`, `oind2pind`, `pind2matprm`.  In doing so, figure out the unique
# indices for the material parameters and use them as pind in `oind2pind` and `pind2matprm`.
#
# When I put an periodic array of objects, some object could be across a boundary.  In that
# case, the rule is to put two different objects at the negative-end and positive-end
# boundaries.  Make sure, though, the two objects have the same surface normal at the
# boundary points, so that the surface normal calculation is not ambiguous there.
#
# Note that pind2matprm stores the same material parameter "value" only once.  Therefore,
# if two objects are created with different Material instances with the same material
# parameter tensor, then the two objects' materials are assigned with the same pind.
function add_obj!(oind2shp::AbsVec{Shape{K,K²}},  # initially empty vector
                  oind2pind::Tuple2{AbsVec{ParamInd}},  # tuple of initially empty vectors
                  pind2matprm::Tuple2{AbsVec{S²ComplexF3}},  # tuple of initially empty vectors
                  obj::Object{K,<:Shape{K,K²}}) where {K,K²}
    push!(oind2shp, shape(obj))  # append obj (for potential use of oind2obj with KDTree, must use pushfirst! to prepend)

    # Assign the material parameter indices to obj.
    # If the material is aleady used in some object that was previously added, then use the
    # index already assigned to that material.
    for ft = EH
        nft = Int(ft)

        mp, on2pn, pn2mp = matparam(obj,ft), oind2pind[nft], pind2matprm[nft]
        pind = findlast(x -> x==mp, pn2mp)  # number if mp is already in pn2mp
        pind==nothing && (push!(pn2mp, mp); pind = length(pn2mp))  # update pind if mp is not in pn2mp
        push!(on2pn, pind)
    end

    return nothing
end

function GeometryPrimitives.periodize(obj::Object{K}, A::AbsMat, ∆range::Shape{K}) where {K}
    shp_array = periodize(obj.shp, A, ∆range)
    N = length(shp_array)
    obj_array = Vector{Object{K}}(N)
    for n = 1:N
        obj_array[n] = Object(shp_array[n], obj.mat, obj.∆lmax)
    end

    return obj_array
end

# The following includes operation on shapes.  These will need to be defined in GeometryPrimitives.
# Then, I will need to define the same functions for Object that delegate the operations to
# Shape's functions.
#
# See also test/object.jl

# # distance to the boundary.  If I want to get 0 for internal points, compare
# # dist(l,c) = (distance between l and center), and if it is less than L_(i) / 2,
# # return zero; otherwise return dist(l,c).
# @inline dist2bound(i::Interval1D, l::Real) = (a1 = abs(i.bound[nN]-l)) < (a2 = abs(i.bound[nP]-l)) ? a1 : a2

# # Functions common to all Shape's
#

# # contains
# # Sometimes, we only have level-set functions rather than Shape instances.
# @inline contains(lsf::Function, l::Tuple3{Real}, isinclusive::Bool=true) = isinclusive ? lsf(l) ≥ 0 : lsf(l) > 0
# @inline contains(s::Shape, l::Tuple3{Real}, isinclusive::Bool=true) = contains(lsf(s), l, isinclusive)
#
# # Other functions
# @inline bound_(s::Shape) = bound_(s.cbox)
# @inline bound_(s::Shape, w::Axis) = bound_(s.cbox, w)
# @inline bound_(s::Shape, w::Axis, sn::Sign) = bound_(s.cbox, w, sn)
# @inline bound_contains(s::Shape, l::Tuple3{Real}, isinclusive::Bool=true) = contains(s.cbox, l, isinclusive)
# @inline L_(s::Shape) = L_(s.cbox)  # scalar
# @inline center_(s::Shape) = center_(s.cbox)  # scalar
#
# # Level-set function manipulation
# @inline transform(lsf::Function, σ::Tuple3{Real}, ∆::Tuple3{Real}) = l -> lsf(σ .* (l .- ∆))  # σ[w] = 1 for shift, σ[w] = -1 for flip
#
# @inline union(lsf::Function...) = l -> maximum(map(f->f(l), lsf))
# @inline union(lsf1::Function, lsf2::Function) = lsf1==lsf2 ? lsf1 : l->max(lsf1(l), lsf2(l))
#
# @inline intersect(lsf::Function...) = l -> minimum(map(f->f(l), lsf))
#
# # Numerical gradient, and closest boundary point to a given point (needed for subpixel
# # smoothing of material parameters)
#
# # Calculates the direction normal to the curved surface crossing a voxel by numerical gradient
# # of the level-set function describing the surface.  The gradient is evaluated at a given point
# # inside the voxel.
# function surfnormal(lsf::Function, r::Tuple3{Float}, vxl::Tuple32{Float})
#     # println("lsf((0,0,0)) = $(lsf((0,0,0))),
#     #     lsf((-1,0,0)) = $(lsf((-1,0,0))), lsf((1,0,0)) = $(lsf((1,0,0))),
#     #     lsf((0,-1,0)) = $(lsf((0,-1,0))), lsf((0,1,0)) = $(lsf((0,1,0))),
#     #     lsf((0,0,-1)) = $(lsf((0,0,-1))), lsf((0,0,1)) = $(lsf((0,0,1))),
#     #     lsf((-1,-1,-1)) = $(lsf((-1,-1,-1))), lsf((1,-1,-1)) = $(lsf((1,-1,-1))),
#     #     lsf((-1,1,-1)) = $(lsf((-1,1,-1))), lsf((1,1,-1)) = $(lsf((1,1,-1))),
#     #     lsf((-1,-1,1)) = $(lsf((-1,-1,1))), lsf((1,-1,1)) = $(lsf((1,-1,1))),
#     #     lsf((-1,1,1)) = $(lsf((-1,1,1))), lsf((1,1,1)) = $(lsf((1,1,1))),
#     #     ")
#     lsf2(s) = 1-hypot(s[1], s[2], s[3])
#     # println("lsf($r) = $(lsf(r))")
#     # println("lsf2($r) = $(lsf2(r))")
#     # println("r = $r, vxl = $vxl")
#     α = 5e-4  # 1e-3 / 2
#
#     ∆x, ∆y, ∆z = α * (vxl[nX][nP] - vxl[nX][nN]), α * (vxl[nY][nP] - vxl[nY][nN]), α * (vxl[nZ][nP] - vxl[nZ][nN])
#     # ∆x, ∆y, ∆z = α .* (t_ind(vxl, nP, nP, nP) .- t_ind(vxl, nN, nN, nN))
#     rx, ry, rz = r
#
#     # println("(∆x, ∆y, ∆z) = $((∆x, ∆y, ∆z))")
#     # println("(rx, ry, rz) = $((rx, ry, rz))")
#     # # println("lsf(($rx, $ry, $rz)) = $(lsf((rx,ry,rz)))")
#     fxn::Float, fxp::Float = lsf((rx-∆x, ry, rz)), lsf((rx+∆x, ry, rz))
#     fyn::Float, fyp::Float = lsf((rx, ry-∆y, rz)), lsf((rx, ry+∆y, rz))
#     fzn::Float, fzp::Float = lsf((rx, ry, rz-∆z)), lsf((rx, ry, rz+∆z))
#     # println("((fxn, fxp), (fyn, fyp), (fzn, fzp)) = $(((fxn, fxp), (fyn, fyp), (fzn, fzp)))")
#
#     fxn2::Float, fxp2::Float = lsf2((rx-∆x, ry, rz)), lsf2((rx+∆x, ry, rz))
#     fyn2::Float, fyp2::Float = lsf2((rx, ry-∆y, rz)), lsf2((rx, ry+∆y, rz))
#     fzn2::Float, fzp2::Float = lsf2((rx, ry, rz-∆z)), lsf2((rx, ry, rz+∆z))
#     # println("((fxn2, fxp2), (fyn2, fyp2), (fzn2, fzp2)) = $(((fxn2, fxp2), (fyn2, fyp2), (fzn2, fzp2)))")
#
#     neg∂f∂x, neg∂f∂y, neg∂f∂z = (fxn - fxp) / ∆x, (fyn - fyp) / ∆y, (fzn - fzp) / ∆z
#
#     return neg∂f∂x, neg∂f∂y, neg∂f∂z  # outward normal is negative gradient
# end
#
# # Find the point on a curved surface by moving from a point along a direction.
# function surfpoint(lsf::Function, r::Tuple3{Float}, v::Tuple3{Float})
#     # println("r = $r, v = $v")
#     f(t) = lsf(r .+ t.*v)
#     tsol, isconverged = newtsol(0., f)
#
#     p = (r[nX] + tsol*v[nX], r[nY] + tsol*v[nY], r[nZ] + tsol*v[nZ])
#
#     return p, isconverged
#     # return r .+ tsol.*v
# end


text:
using ChainPartitioners: destep

function Δ(i, is)
    moves = [ChainPartitioners.Same(i), ChainPartitioners.Same(i), ChainPartitioners.Jump(rand(is))]
    if i - 1 in is
        push!(moves, ChainPartitioners.Prev(i - 1))
    end
    if i + 1 in is
        push!(moves, ChainPartitioners.Next(i + 1))
    end
    return moves
end

ref_dominancecount(A, i, j) = sum(A[1:i - 1, 1:j - 1] .!= 0)
ref_dominancesum(A, i, j) = sum(A[1:i - 1, 1:j - 1])

test_points(trials, Is...) = collect(Iterators.flatten((
    Iterators.product(map(I -> (first(I), last(I)), Is)...),
    zip(map(I -> rand(I, trials), Is)...)
)))[randperm(end)]

test_dims = [1:16..., 31, 32, 33, 63, 64, 65]

test_dims = [1:3..., 7, 8, 9]

@testset "SparsePrefixMatrices" begin
    @testset "jump" begin
        for (hint, kwargs) = (
            ((), ()),
            ((), (H = 1,)),
            ((ChainPartitioners.SparseHint(),), (H = 1,)),
            ((ChainPartitioners.SparseHint(),), (H = 2,)),
            ((ChainPartitioners.SparseHint(),), (H = 3,)),
            ((ChainPartitioners.SparseHint(),), (H = 4,)),
            ((ChainPartitioners.SparseHint(),), (b = 1,)),
            ((ChainPartitioners.SparseHint(),), (b = 2,)),
            ((ChainPartitioners.SparseHint(),), (b = 3,)),
            ((ChainPartitioners.SparseHint(),), (b = 4,)),
            ((ChainPartitioners.RandomHint(),), ()),
            ((), (b = 1,)),
        )
            for m = test_dims
                for n = test_dims
                    A = dropzeros!(sprand(UInt, m, n, 0.5))

                    C = dominancecount(hint..., A; kwargs...)
                    S = dominancesum(hint..., A; kwargs...)
                    @testset "kwargs = $kwargs, m = $m" begin
                        for (i, j) in test_points(10, 1:m + 1, 1:n + 1)
                            @test C[i, j] == ref_dominancecount(A, i, j)
                            @test S[i, j] == ref_dominancesum(A, i, j)
                        end
                    end
                end

                N = m
                idx = randperm(N)
                B = spzeros(UInt, N, N)
                for j = 1:N
                    B[idx[j], j] = (rand(UInt) << 1) + 1
                end
                RC = rookcount!(hint..., N, copy(idx); kwargs...)
                RS = rooksum!(hint..., N, copy(idx), B.nzval; kwargs...)
                @testset "kwargs = $kwargs, N = $N" begin
                    for (i, j) in test_points(100, 1:N + 1, 1:N + 1)
                        @test RC[i, j] == ref_dominancecount(B, i, j)
                        @test RS[i, j] == ref_dominancesum(B, i, j)
                    end
                end
            end
        end
    end

    @testset "step" begin
        for m = 1:40
            for n = 1:40
                A = dropzeros!(sprand(Int, m, n, 0.5))
                C = dominancecount(ChainPartitioners.StepHint(), A)
                @testset "m = $m" begin
                    for (i, j) = test_points(10, 1:m + 1, 1:n + 1)
                        @test C[i, j] == sum(A[1:i - 1, 1:j - 1] .!= 0)

                        for _ = 1:8
                            (_i, _j) = (rand(Δ(i, 1:m + 1)), rand(Δ(j, 1:n + 1)))
                            (i, j) = map(destep, (_i, _j))
                            @test ChainPartitioners.Step(C)(_i, _j) == ref_dominancecount(A, i, j)
                        end
                    end
                end
            end
        end
    end
end


text:
############# function properties #############
export isreversible, isreflexive, isprimitive
export protectf
"""
    isreversible(f, ARGT)

Return `true` if a function is reversible.
"""
isreversible(f, ::Type{ARGT}) where ARGT = hasmethod(~f, ARGT)

"""
    isreflexive(f)

Return `true` if a function is self-inverse.
"""
isreflexive(f) = (~f) === f

"""
    isprimitive(f)

Return `true` if `f` is an `instruction` that can not be decomposed anymore.
"""
isprimitive(f) = false

############# ancillas ################
export InvertibilityError, @invcheck

"""
    deanc(a, b)

Deallocate varialbe `a` with value `b`. It will throw an error if

* `a` and `b` are objects with different types,
* `a` is not equal to `b` (for floating point numbers, an error within `NiLangCore.GLOBAL_ATOL[]` is allowed),
"""
function deanc end

function deanc(a::T, b::T) where T <: AbstractFloat
    if a !== b && abs(b - a) > GLOBAL_ATOL[]
        throw(InvertibilityError("deallocate fail (floating point numbers): $a ≂̸ $b"))
    end
end
deanc(x::T, val::T) where T<:Tuple = deanc.(x, val)
deanc(x::T, val::T) where T<:AbstractArray = x === val || deanc.(x, val)
deanc(a::T, b::T) where T<:AbstractString = a === b || throw(InvertibilityError("deallocate fail (string): $a ≂̸ $b"))
function deanc(x::T, val::T) where T<:Dict
    if x !== val
        if length(x) != length(val)
            throw(InvertibilityError("deallocate fail (dict): length of dict not the same, got $(length(x)) and $(length(val))!"))
        else
            for (k, v) in x
                if haskey(val, k)
                    deanc(x[k], val[k])
                else
                    throw(InvertibilityError("deallocate fail (dict): key $k of dict does not exist!"))
                end
            end
        end
    end
end
deanc(a, b) = throw(InvertibilityError("deallocate fail (type mismatch): `$(typeof(a))` and `$(typeof(b))`"))

@generated function deanc(a::T, b::T) where T
    nf = fieldcount(a)
    if isprimitivetype(T)
        :(a === b || throw(InvertibilityError("deallocate fail (primitive): $a ≂̸ $b")))
    else
        Expr(:block, [:($deanc(a.$NAME, b.$NAME)) for NAME in fieldnames(T)]...)
    end
end

"""
    InvertibilityError <: Exception
    InvertibilityError(ex)

The error for irreversible statements.
"""
struct InvertibilityError <: Exception
    ex
end

"""
    @invcheck x val

The macro version `NiLangCore.deanc`, with more informative error.
"""
macro invcheck(x, val)
    esc(_invcheck(x, val))
end

# the expression for reversibility checking
function _invcheck(x, val)
    Expr(:try, Expr(:block, :($deanc($x, $val))), :e, Expr(:block,
            :(println("deallocate fail `$($(QuoteNode(x))) → $($(QuoteNode(val)))`")),
            :(throw(e)))
        )
end
_invcheck(docheck::Bool, arg, res) = docheck ? _invcheck(arg, res) : nothing

"""
    chfield(x, field, val)

Change a `field` of an object `x`.

The `field` can be a `Val` type

```jldoctest; setup=:(using NiLangCore)
julia> chfield(1+2im, Val(:im), 5)
1 + 5im
```

or a function

```jldoctest; setup=:(using NiLangCore)
julia> using NiLangCore

julia> struct GVar{T, GT}
           x::T
           g::GT
       end

julia> @fieldview xx(x::GVar) = x.x

julia> chfield(GVar(1.0, 0.0), xx, 2.0)
GVar{Float64, Float64}(2.0, 0.0)
```
"""
function chfield end

########### Inv  ##########
export Inv, invtype
"""
    Inv{FT} <: Function
    Inv(f)

The inverse of a function.
"""
struct Inv{FT} <: Function
    f::FT
end
Inv(f::Inv) = f.f
@static if VERSION >= v"1.6"
    Base.:~(f::Base.ComposedFunction) = (~(f.inner)) ∘ (~(f.outer))
end
Base.:~(f::Function) = Inv(f)
Base.:~(::Type{Inv{T}}) where T = T  # for type, it is a destructor
Base.:~(::Type{T}) where T = Inv{T}  # for type, it is a destructor
Base.show(io::IO, b::Inv) = print(io, "~$(b.f)")
Base.display(bf::Inv) where f = print(bf)
"""
    protectf(f)

Protect a function from being inverted, useful when using an callable object.
"""
protectf(x) = x
protectf(x::Inv) = x.f

invtype(::Type{T}) where T = Inv{<:T}

######### Infer
export PlusEq, MinusEq, XorEq, MulEq, DivEq
"""
    PlusEq{FT} <: Function
    PlusEq(f)

Called when executing `out += f(args...)` instruction. The following two statements are same

```jldoctest; setup=:(using NiLangCore)
julia> x, y, z = 0.0, 2.0, 3.0
(0.0, 2.0, 3.0)

julia> x, y, z = PlusEq(*)(x, y, z)
(6.0, 2.0, 3.0)

julia> x, y, z = 0.0, 2.0, 3.0
(0.0, 2.0, 3.0)

julia> @instr x += y*z

julia> x, y, z
(6.0, 2.0, 3.0)
```
"""
struct PlusEq{FT} <: Function
    f::FT
end

"""
    MinusEq{FT} <: Function
    MinusEq(f)

Called when executing `out -= f(args...)` instruction. See `PlusEq` for detail.
"""
struct MinusEq{FT} <: Function
    f::FT
end

"""
    MulEq{FT} <: Function
    MulEq(f)

Called when executing `out *= f(args...)` instruction. See `PlusEq` for detail.
"""
struct MulEq{FT} <: Function
    f::FT
end

"""
    DivEq{FT} <: Function
    DivEq(f)

Called when executing `out /= f(args...)` instruction. See `PlusEq` for detail.
"""
struct DivEq{FT} <: Function
    f::FT
end

"""
    XorEq{FT} <: Function
    XorEq(f)

Called when executing `out ⊻= f(args...)` instruction. See `PlusEq` for detail.
"""
struct XorEq{FT} <: Function
    f::FT
end

const OPMX{FT} = Union{PlusEq{FT}, MinusEq{FT}, XorEq{FT}, MulEq{FT}, DivEq{FT}}

for (TP, OP) in [(:PlusEq, :+), (:MinusEq, :-), (:XorEq, :⊻)]
    @eval (inf::$TP)(out!, args...; kwargs...) = $OP(out!, inf.f(args...; kwargs...)), args...
    @eval (inf::$TP)(out!::Tuple, args...; kwargs...) = $OP.(out!, inf.f(args...; kwargs...)), args...  # e.g. allow `(x, y) += sincos(a)`
end

Base.:~(op::PlusEq) = MinusEq(op.f)
Base.:~(om::MinusEq) = PlusEq(om.f)
Base.:~(op::MulEq) = DivEq(op.f)
Base.:~(om::DivEq) = MulEq(om.f)
Base.:~(om::XorEq) = om
for (T, S) in [(:PlusEq, "+="), (:MinusEq, "-="), (:MulEq, "*="), (:DivEq, "/="), (:XorEq, "⊻=")]
    @eval Base.display(o::$T) = print($S, "(", o.f, ")")
    @eval Base.display(o::Type{$T}) = print($S)
    @eval Base.show_function(io::IO, o::$T, compact::Bool) = print(io, "$($S)($(o.f))")
    @eval Base.show_function(io::IO, ::MIME"plain/text", o::$T, compact::Bool) = Base.show(io, o)
end


text:
using EnhancedGJK
using MeshCat
using Polyhedra
using StaticArrays: SVector
using GeometryTypes: HyperSphere, Point

function visualize_simplex(vis::Visualizer, simplex)
    p = polyhedron(vrep(simplex))
    setobject!(vis[:simplex], Polyhedra.Mesh(p))
    for (i, point) in enumerate(simplex)
        setobject!(vis["p$i"], HyperSphere(Point(point), 0.03))
    end
end

vis = Visualizer()
open(vis)
simplex = SVector{4, SVector{3, Float64}}([-3.33728783796428, 0.3321305518800686, 0.11004228580261734], [-3.33728783796428, -0.6678694481199314, 0.11004228580261734], [1.6627121620357201, -0.6678694481199314, 0.11004228580261734], [1.6627121620357201, 0.3321305518800686, 0.11004228580261734])
visualize_simplex(vis, simplex)
@show EnhancedGJK.projection_weights(simplex)